最小二乘问题详解4：非线性最小二乘

1. 引言

在论述最小二乘问题的时候，很多文章都喜欢用拟合直线来举例，但是在现实中像拟合直线这样的线性最小二乘问题往往不是常态，现实世界中更多是像投影成像这种非线性最小二乘问题。在本文中，我们就讲解一下非线性最小二乘问题。

不过，在继续阅读本文之前，一定要先理解之前的3篇文章，因为线性最小二乘是求解非线性最小二乘问题的基础：

1. 《最小二乘问题详解1：线性最小二乘》
2. 《最小二乘问题详解2：线性最小二乘求解》
3. 《最小二乘问题详解3：线性最小二乘实例》

2. 定义

具体来说，非线性最小二乘的目标就是找到一组参数 \(\theta\)，使得非线性模型 \(f(x; \theta)\) 最好地拟合观测数据。非线性最小二乘模型通常使用如下公式来表达：

\[\mathbf{y} = f(\mathbf{x}; \theta) + \varepsilon \]

其中：

• \(\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{x}_i, \mathbf{y}_i) \}_{i=1}^n \]

同样的，非线性最小二乘的目标是最小化残差平方和：

\[\boxed{ \min_{\theta} S(\theta) = \min_{\theta} \sum_{i=1}^n \| \mathbf{r}_i(\theta) \|^2 = \min_{\theta} \| \mathbf{r}(\theta) \|^2 } \]

其中：

• \(\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\) 附近，我们可以用一阶泰勒展开近似模型输出：

\[f(\mathbf{x}_i; \theta) \approx f(\mathbf{x}_i; \theta_k) + J_i(\theta_k) (\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)\)，得到残差的线性近似：

\[\mathbf{r}_i(\theta) \approx \mathbf{r}_i(\theta_k) - J_i(\theta_k) (\theta - \theta_k) \]

对所有样本拼接成向量形式：

\[\boxed{ \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 \]

这是一个关于 \(\Delta\theta\) 的线性模型，为非线性最小二乘求解奠定了基础。

3. 雅可比矩阵

在进行正式求解之前，我们先理解一下雅可比矩阵（Jacobian Matrix），因为它在非线性最小二乘中起着核心作用。

从泰勒展开的角度看，一元函数的一阶近似使用导数，而多元向量函数的一阶近似则使用雅可比矩阵——它是多元函数所有一阶偏导数组成的矩阵。在非线性最小二乘中，我们关注的是模型输出 \(f(\mathbf{x}_i; \theta)\) 对参数 \(\theta\) 的变化率。

对于第 \(i\) 个样本，其模型输出对参数的雅可比块为：

\[J_i(\theta) = \frac{\partial f(\mathbf{x}_i; \theta)}{\partial \theta^T} \in \mathbb{R}^{m \times p} \]

• \(m\)：模型输出维度（例如，二维图像坐标 \((u,v)\) 则 \(m=2\)）
• \(p\)：待估参数维度
• \(J_i(\theta)\) 是一个 \(m \times p\) 矩阵，表示模型输出的每个分量对每个参数的偏导数

将所有 \(n\) 个样本的雅可比块垂直堆叠，得到整体雅可比矩阵：

\[\boxed{ J(\theta) = \begin{bmatrix} J_1(\theta) \\ J_2(\theta) \\ \vdots \\ J_n(\theta) \end{bmatrix}= \begin{bmatrix} \frac{\partial f(\mathbf{x}_1; \theta)}{\partial \theta^T} \\ \frac{\partial f(\mathbf{x}_2; \theta)}{\partial \theta^T} \\ \vdots \\ \frac{\partial f(\mathbf{x}_n; \theta)}{\partial \theta^T} \end{bmatrix} \in \mathbb{R}^{N \times p} } \quad \text{其中 } N = n \cdot m \]

4. Gauss-Newton

求解非线性最小二乘问题最基础最好理解的就是Gauss-Newton方法，它结合了牛顿法的迭代优化框架（就是高中数学中迭代逼近求解平方根的过程）和高斯的线性化思想，所以将其称为Gauss-Newton方法。它的算法流程如下：

1. 初始化：选一个初始猜测 \(\theta_0\)

2. 迭代：对 \(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 \]

3. 终止：当 \(\|\Delta \theta\|\) 很小或目标函数变化很小时停止

Gauss-Newton方法的基本思路是每一步都在当前点附近做线性近似；然后走一步，走的方向是让残差平方和下降最快的方向之一；最终收敛到一个局部最小值。可能这个算法流程可能还是有点抽象，我们将其中一次迭代的过程论述的更清楚一些：

1. 在 \(\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)\)

2. 构造局部二次近似目标函数。代入原目标：

   \[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\) 的二次函数。

3. 最小化这个二次函数。对 \(\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 } \]

4. 求解 Gauss-Newton 步长。如果 \(J_k^T J_k\) 可逆，则解为：

   \[\boxed{ \Delta \theta = (J_k^T J_k)^{-1} J_k^T \mathbf{r}_k } \]

5. 更新参数

   \[\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\)
是否迭代	否	是