最小二乘問題詳解4:非線性最小二乘
1. 引言
在論述最小二乘問題的時候,很多文章都喜歡用擬合直線來舉例,但是在現實中像擬合直線這樣的線性最小二乘問題往往不是常態,現實世界中更多是像投影成像這種非線性最小二乘問題。在本文中,我們就講解一下非線性最小二乘問題。
不過,在繼續閱讀本文之前,一定要先理解之前的3篇文章,因為線性最小二乘是求解非線性最小二乘問題的基礎:
2. 定義
具體來說,非線性最小二乘的目標就是找到一組參數 \(\theta\),使得非線性模型 \(f(x; \theta)\) 最好地擬合觀測數據。非線性最小二乘模型通常使用如下公式來表達:
其中:
- \(\mathbf{x} \in \mathbb{R}^d\):輸入變量(如坐標、時間等)
- \(\theta \in \mathbb{R}^p\):待估計參數向量
- \(f: \mathbb{R}^d \times \mathbb{R}^p \to \mathbb{R}^m\):非線性向量值函數
- \(\varepsilon\):觀測噪聲
- \(\mathbf{y} \in \mathbb{R}^m\):觀測輸出
我們有 \(n\) 組數據:
同樣的,非線性最小二乘的目標是最小化殘差平方和:
其中:
- \(\mathbf{r}_i(\theta) = \mathbf{y}_i - f(\mathbf{x}_i; \theta)\):第 \(i\) 個樣本的殘差向量
- \(\mathbf{r}(\theta) = [\mathbf{r}_1^T, \mathbf{r}_2^T, \dots, \mathbf{r}_n^T]^T\):所有殘差拼成的長向量,維度 \(N = n \cdot m\)
那么非線性最小二乘能不能像線性那樣直接求解呢?顯然是不行的。因為在線性情況下,殘差是 \(\mathbf{r} = A\theta - b\),目標函數是 \(\|A\theta - b\|^2\),這是一個二次函數,所以有閉式解 \((A^T A)^{-1} A^T b\)。但在非線性情況下,\(\| \mathbf{r}(\theta) \|^2\) 是一個非凸、非二次函數,無法直接求逆,也就無法計算出閉式解。更直觀的說,通過非線性函數\(f(\mathbf{x}; \theta)\)是無法寫成類似設計矩陣\(A\theta=b\)這樣的線性方程組的。
一種求解的思路是局部線性化(Local Linearization)。雖然 \(f(\mathbf{x}; \theta)\) 是非線性的,但在當前估計 \(\theta_k\) 附近,我們可以用一階泰勒展開近似模型輸出:
其中:
- \(J_i(\theta_k) = \frac{\partial f(\mathbf{x}_i; \theta)}{\partial \theta^T} \big|_{\theta = \theta_k} \in \mathbb{R}^{m \times p}\):第 \(i\) 個樣本的模型輸出對參數的雅可比塊。
- \(J(\theta_k)\) 是將所有 \(J_i\) 垂直堆疊得到的 \(N \times p\) 矩陣(\(N = n \cdot m\),\(p\) 為參數維度)。
代入殘差定義 \(\mathbf{r}_i(\theta) = \mathbf{y}_i - f(\mathbf{x}_i; \theta)\),得到殘差的線性近似:
對所有樣本拼接成向量形式:
令 \(\Delta\theta = \theta - \theta_k\),則:
這是一個關于 \(\Delta\theta\) 的線性模型,為非線性最小二乘求解奠定了基礎。
3. 雅可比矩陣
在進行正式求解之前,我們先理解一下雅可比矩陣(Jacobian Matrix),因為它在非線性最小二乘中起著核心作用。
從泰勒展開的角度看,一元函數的一階近似使用導數,而多元向量函數的一階近似則使用雅可比矩陣——它是多元函數所有一階偏導數組成的矩陣。在非線性最小二乘中,我們關注的是模型輸出 \(f(\mathbf{x}_i; \theta)\) 對參數 \(\theta\) 的變化率。
對于第 \(i\) 個樣本,其模型輸出對參數的雅可比塊為:
- \(m\):模型輸出維度(例如,二維圖像坐標 \((u,v)\) 則 \(m=2\))
- \(p\):待估參數維度
- \(J_i(\theta)\) 是一個 \(m \times p\) 矩陣,表示模型輸出的每個分量對每個參數的偏導數
將所有 \(n\) 個樣本的雅可比塊垂直堆疊,得到整體雅可比矩陣:
4. Gauss-Newton
求解非線性最小二乘問題最基礎最好理解的就是Gauss-Newton方法,它結合了牛頓法的迭代優化框架(就是高中數學中迭代逼近求解平方根的過程)和高斯的線性化思想,所以將其稱為Gauss-Newton方法。它的算法流程如下:
-
初始化:選一個初始猜測 \(\theta_0\)
-
迭代:對 \(k = 0, 1, 2, \dots\)
a. 計算當前殘差:\(\mathbf{r}_k = \mathbf{y} - f(\mathbf{x}; \theta_k)\)
b. 計算雅可比矩陣 \(J_k = J(\theta_k)\)
c. 求解線性最小二乘子問題:\[\min_{\Delta \theta} \| \mathbf{r}_k - J_k \Delta \theta \|^2 \]解為:
\[\Delta \theta = (J_k^T J_k)^{-1} J_k^T \mathbf{r}_k \]d. 更新參數:
\[\theta_{k+1} = \theta_k + \Delta \theta \] -
終止:當 \(\|\Delta \theta\|\) 很小或目標函數變化很小時停止
Gauss-Newton方法的基本思路是每一步都在當前點附近做線性近似;然后走一步,走的方向是讓殘差平方和下降最快的方向之一;最終收斂到一個局部最小值。可能這個算法流程可能還是有點抽象,我們將其中一次迭代的過程論述的更清楚一些:
-
在 \(\theta_k\) 處做一階泰勒展開。對殘差函數 \(\mathbf{r}(\theta)\) 在 \(\theta_k\) 附近線性化:
\[\mathbf{r}(\theta) \approx \mathbf{r}(\theta_k) - J(\theta_k) (\theta - \theta_k) \]令 \(\Delta \theta = \theta - \theta_k\),則:
\[\mathbf{r}(\theta) \approx \mathbf{r}_k - J_k \Delta \theta \]其中:
- \(\mathbf{r}_k = \mathbf{r}(\theta_k)\)
- \(J_k = J(\theta_k)\)
-
構造局部二次近似目標函數。代入原目標:
\[S(\theta) = \| \mathbf{r}(\theta) \|^2 \approx \| \mathbf{r}_k - J_k \Delta \theta \|^2 \]展開:
\[S(\theta) \approx \mathbf{r}_k^T \mathbf{r}_k - 2 \mathbf{r}_k^T J_k \Delta \theta + \Delta \theta^T J_k^T J_k \Delta \theta \]這是一個關于 \(\Delta \theta\) 的二次函數。
-
最小化這個二次函數。對 \(\Delta \theta\) 求導并令導數為零:
\[\frac{\partial S}{\partial (\Delta \theta)} = -2 J_k^T \mathbf{r}_k + 2 J_k^T J_k \Delta \theta = 0 \]得到正規方程(Normal Equation):
\[\boxed{ J_k^T J_k \Delta \theta = J_k^T \mathbf{r}_k } \] -
求解 Gauss-Newton 步長。如果 \(J_k^T J_k\) 可逆,則解為:
\[\boxed{ \Delta \theta = (J_k^T J_k)^{-1} J_k^T \mathbf{r}_k } \] -
更新參數
\[\theta_{k+1} = \theta_k + \Delta \theta \]
看到這個過程中的正規方程了嗎?這就是我們說的非線性最小二乘求解的基礎是線性最小二乘的原因了,非線性最小二乘問題的每次迭代過程就是一個線性最小二乘子問題。
非線性最小二乘與線性最小二乘求解過程的對比如下:
| 特性 | 線性最小二乘 | 非線性最小二乘(Gauss-Newton) |
|---|---|---|
| 模型 | \(f(x; \theta) = A \theta\) | \(f(x; \theta)\) 任意非線性 |
| 殘差 | \(\mathbf{r} = A\theta - b\) | \(\mathbf{r}(\theta) = \mathbf{y} - f(\mathbf{x}; \theta)\) |
| 雅可比 | \(J = A\)(常數) | \(J(\theta) = \frac{\partial f}{\partial \theta^T}\)(依賴 \(\theta\)) |
| 解 | \(\theta^* = (A^T A)^{-1} A^T b\) | \(\theta_{k+1} = \theta_k + (J_k^T J_k)^{-1} J_k^T \mathbf{r}_k\) |
| 是否迭代 | 否 | 是 |

浙公網安備 33010602011771號