非线性最小二乘

目录

• • 1. Gauss-Newton 迭代
    • 1.1 算法原理
    • 1.2 迭代步骤
    • 1.3 改进
    • 1.4 如何使用
  • 2. Levenberg–Marquardt 算法
    • 2.1 算法原理
    • 2.2 迭代步骤
    • 2.3 如何使用

涉及 Gauss-Newton 迭代 与 LM 非线性最小二乘求解算法

作者

赵曦

日期

2024/05/24

版本

1.0

上一篇教程：最小二乘法 ↑
下一篇教程：非线性方程（组）数值解与迭代法 ↓

---

\[ \def\red#1{\color{red}{#1}} \def\teal#1{\color{teal}{#1}} \def\green#1{\color{green}{#1}} \def\transparent#1{\color{transparent}{#1}} \def\orange#1{\color{orange}{#1}} \def\Var{\mathrm{Var}} \def\Cov{\mathrm{Cov}} \def\tr{\mathrm{tr}} \def\fml#1{\text{(#1)}} \def\ptl#1#2{\frac{\partial#1}{\partial#2}} \]

1. Gauss-Newton 迭代

1.1 算法原理

数据处理中，最常见的一种函数形式是

\[f(\pmb x)=\frac12\sum_{i=1}^m\varphi_i^2(\pmb x)\tag{1-1}\]

如果 \(f(x)\) 的极小点 \(x^*\) 满足 \(f(x^*)< e\)（ \(e\) 是预先给定的精度），那么可以认为 \(x^*\) 是 方程组

\[\varphi_i(x_1,x_2,\cdots,x_n)=0,\quad i=1,2,\cdots,m(m\le n)\]

的解。

对式 \(\fml{1-1}\)，可以用一阶导数运算来代替牛顿法中的二阶导数矩阵的求逆运算。因为

\[\begin{align}\nabla f(\pmb x)&=\begin{bmatrix} \sum\limits_{i=1}^m\ptl{\varphi_i(\pmb x)}{x_1}\varphi_i(\pmb x)\\ \sum\limits_{i=1}^m\ptl{\varphi_i(\pmb x)}{x_2}\varphi_i(\pmb x)\\ \vdots\\\sum\limits_{i=1}^m\ptl{\varphi_i(\pmb x)}{x_n}\varphi_i(\pmb x)\end{bmatrix}=\begin{bmatrix} \ptl{\varphi_1(\pmb x)}{x_1}&\ptl{\varphi_2(\pmb x)}{x_1}&\cdots&\ptl{\varphi_m(\pmb x)}{x_1}\\ \ptl{\varphi_1(\pmb x)}{x_2}&\ptl{\varphi_2(\pmb x)}{x_2}&\cdots&\ptl{\varphi_m(\pmb x)}{x_2}\\ \vdots&\vdots&\ddots&\vdots\\ \ptl{\varphi_1(\pmb x)}{x_n}&\ptl{\varphi_2(\pmb x)}{x_n}&\cdots&\ptl{\varphi_m(\pmb x)}{x_n} \end{bmatrix}\begin{bmatrix} \varphi_1(\pmb x)\\\varphi_2(\pmb x)\\\vdots\\\varphi_m(\pmb x)\end{bmatrix}\\ &=\pmb J^T(\pmb x)\pmb\varphi(\pmb x)\tag{1-2}\end{align}\]

式中， \(\pmb\varphi(\pmb x)\) 为 \(m\) 维的函数向量，即 \(\pmb\varphi(\pmb x)=[\varphi_1(x),\varphi_2(x),\cdots,\varphi_m(x)]^T\)； \(\pmb J(\pmb x)\) 为函数 \(\pmb\varphi(\pmb x)\ (i=1,2,\cdots,m)\) 一阶偏导数组成的矩阵，即

\[\pmb J(\pmb x)=\begin{bmatrix} \ptl{\varphi_1(\pmb x)}{x_1}&\ptl{\varphi_1(\pmb x)}{x_2}&\cdots&\ptl{\varphi_1(\pmb x)}{x_n}\\ \ptl{\varphi_2(\pmb x)}{x_1}&\ptl{\varphi_2(\pmb x)}{x_2}&\cdots&\ptl{\varphi_2(\pmb x)}{x_n}\\ \vdots&\vdots&\ddots&\vdots\\ \ptl{\varphi_m(\pmb x)}{x_1}&\ptl{\varphi_m(\pmb x)}{x_2}&\cdots&\ptl{\varphi_m(\pmb x)}{x_n} \end{bmatrix}\tag{1-3}\]

特别当 \(\pmb\varphi(\pmb x)\) 是 \(\pmb x\) 的线性函数时，式 \(\fml{1-2}\) 中的 \(\pmb J(\pmb x)\) 是常系数矩阵。这时式 \(\fml{1-1}\) 的海赛矩阵可以写成

\[\pmb H=\pmb J^T\pmb J\tag{1-4}\]

这样，关于 \(\pmb\varphi(\pmb x)\) 是 \(\pmb x\) 的线性函数时的最小二乘法的迭代公式可以写为

\[\pmb x_{k+1}=\pmb x_k-\pmb H^{-1}(\pmb x_k)\pmb J^T(\pmb x_k)\pmb\varphi(\pmb x_k)\tag{1-5}\]

当 \(\pmb\varphi(\pmb x)\) 不是 \(\pmb x\) 的线性函数时，也可以近似将式 \(\fml{1-4}\) 视为函数 \(f(\pmb x)\) 的海赛矩阵，即

\[\pmb H\approx\left[\pmb J(\pmb x_k)\right]^T\pmb J(\pmb x_k)\tag{1-6}\]

所以，关于 \(\pmb\varphi(\pmb x)\) 不是 \(\pmb x\) 的线性函数时的最小二乘法的迭代公式也可以写为 \(\fml{1-5}\) 的形式。

1.2 迭代步骤

1. 选择一初始点 \(\pmb x_0=(x_{1,0},x_{2,0},\cdots,x_{n,0})^T\)；
2. 算出

   \[\Delta\pmb x_0=-\pmb H_0^{-1}\pmb J_0^T\pmb\varphi(\pmb x_k)\tag{1-6}\]

3. 令 \(\pmb x_1\) 为函数 \(f(\pmb x)\) 的极小点的第 1 次近似，则有

   \[\pmb x_1=\pmb x_0+\Delta\pmb x_0\tag{1-7}\]

4. 以 \(\pmb x_1\) 代替前面的 \(\pmb x_0\)， \(\Delta\pmb x_1\) 代替 \(\Delta\pmb x_0\)，重复上述计算过程，直到

   \[\|\pmb\varphi(\pmb x_k)\|<\epsilon'\tag{1-8a}\]

   或

   \[\|\nabla f(\pmb x_k)\|<\epsilon''\tag{1-8b}\]

   为止。 \(\epsilon'\) 和 \(\epsilon''\) 是预先给定的精度。

1.3 改进

上述高斯－牛顿最小二乘法利用了目标函数在极小值处近似为自变量各分量的平方和的特点，用 \(\pmb J^T\pmb J\) 近似代替牛顿法中 \(f(\pmb x)\) 的二阶导数矩阵，大大节省了计算量。但是它对初始点 \(\pmb x_0\) 的要求比较严格，如果初始点 \(\pmb x_0\) 与极小点 \(\pmb x^*\) 相距很远，这个算法往往失败。原因是

1. 上述算法基于线性逼近，但在 \(\pmb x_0\) 远离极小点时，这种线性逼近无效；
2. \(\pmb J_0^T\pmb J_0\) 的最大特征值与最小特征值的比很大，致使解 \(\Delta\pmb x_0\) 变得无意义。

为此采取下述改进的办法。在求出 \(\pmb x_k\) 的校正量 \(\Delta\pmb x_k\) 后，不把 \(\pmb x_k+\Delta\pmb x_k\) 作为第 \(k+1\) 次近似点，而是将 \(\Delta\pmb x_k\) 作为下一步的一维方向搜索。求 \(\alpha_k\)，使

\[f(\pmb x_k+\alpha_k\pmb s_k)=\min_{\alpha>0}f(\pmb x_k+\alpha\pmb s_k)\]

然后令

\[\pmb x_{k+1}=\pmb x_k+\alpha_k\pmb s_k\]

以 \(\pmb x_{k+1}\) 代替 \(\pmb x_k\) 重复上述计算过程，直到 \(\|\pmb\varphi(\pmb x_k)\|<\epsilon'\) 或 \(\|\nabla f(\pmb x_k)\|<\epsilon"\) 为止。

1.4 如何使用

RMVL 提供了改进的 Gauss-Newton 迭代算法，可参考 rm::lsqnonlin 函数。例如，我们需要拟合一个正弦函数

\[y=A\sin(\omega t+\varphi_0)+b\]

其中， \(A,\omega,\varphi_0,b\) 是待拟合的参数，不妨统一写为 \(\pmb x=(A,\omega,\varphi_0,b)\)，也就是说我们需要拟合的函数是

\[\green y=x_1\sin(x_2\green t+x_3)+x_4\]

其中 \(t\) 和 \(y\) 是可以观测到的数据，我们需要通过观测的数据来拟合 \(\pmb x\) 的值。比方说，下面的 obtain 函数就可以观测每一帧的数据。

C++

double obtain();

Python

def obtain(): ...

例如经过了 20 帧的数据采集，我们得到了一个长度为 20 的队列，即

C++

std::deque<double> datas;
 
/* code */
datas.push_front(obtain());
if (datas.size() == 20)
    datas.pop_back();
/* code */

Python

datas = []
 
# code
datas.insert(0, obtain())
if len(datas) == 20:
    datas.pop()
# code

准备好数据后，可以使用下面的代码来拟合正弦函数。

C++

rm::FuncNds lsq_sine(datas.size());
for (std::size_t i = 0; i < datas.size(); ++i)
    lsq_sine.push_back([=](const std::vector<double> &x) {
        return x[0] * std::sin(x[1] * i + x[2]) + x[3] - datas[i];
    });
 
// 拟合正弦函数，初始值为 (1, 0.02, 0, 1.09)
auto x = rm::lsqnonlin(lsq_sine, {1, 0.02, 0, 1.09}); // 默认采用 Gauss-Newton 算法

Python

lsq_sine = []
 
for i, data in enumerate(datas):
    lsq_sine.append(
        lambda x, i=i, data=data: x[0] * np.sin(x[1] * i + x[2]) + x[3] - data
    )
 
# 拟合正弦函数，初始值为 (1, 0.02, 0, 1.09)
x = rm.lsqnonlin(lsq_sine, [1, 0.02, 0, 1.09]) # 默认采用 Gauss-Newton 算法

2. Levenberg–Marquardt 算法

2.1 算法原理

普通的 Gauss-Newton 迭代，在初始值附近做了一阶线性化处理，而当初始值与极小值相差较远时，曲线的非线性特性会导致迭代失败。为了解决这个问题，Levenberg–Marquardt 算法在 Gauss-Newton 迭代的基础上，引入了一个参数 \(\lambda\)，使得迭代公式变为

\[\pmb x_{k+1}=\pmb x_k-\left(\pmb J^T\pmb J+\red{\lambda\pmb I}\right)^{-1}\pmb J^T\pmb\varphi(\pmb x_k)\tag{2-1}\]

\(I\) 是单位矩阵， \(\lambda\) 是一个非负数。如果 \(\lambda\) 取值较大时， \(\lambda I\) 占主要地位，此时的 LM 算法更接近一阶梯度下降法，说明此时距离最终解还比较远，用一阶近似更合适。反之，如果 \(\lambda\) 取值较小时， \(\pmb H=\pmb J^T\pmb J\) 占主要地位，说明此时距离最终解距离较近，用二阶近似模型比较合适，可以避免梯度下降的震荡，容易快速收敛到极值点。因此参数 \(\lambda\) 不仅影响到迭代的方向还影响到迭代步长的大小。

令初值为 \(\pmb x_0\)，可以设置 \(\lambda\) 的初值 \(\lambda_0\) 为

\[\begin{align}\pmb A_0&=\pmb J^T(\pmb x_0)\pmb J(\pmb x_0)\\ \lambda_0&=\tau\max_i\left\{a_{ii}^{(0)}\right\}\end{align}\tag{2-2}\]

其中， \(\tau\) 可以自己指定， \(a_{ii}\) 表示矩阵 \(\pmb A\) 对角线元素。此外， \(\lambda\) 需要在迭代过程中不断调整，以保证迭代的收敛性。一般会判断近似的模型与实际函数之间的差异，可以使用下面的公式来判断

\[\rho_k=\frac{f(\pmb x_k+\Delta\pmb x_k)-f(\pmb x_k)}{\pmb J(\pmb x_k)\Delta\pmb x_k}\tag{2-3}\]

2.2 迭代步骤

1. 选择一初始点 \(\pmb x_0=(x_{1,0},x_{2,0},\cdots,x_{n,0})^T\)，按照式 \(\fml{2-2}\) 计算 \(\lambda_0\)；
2. 对于第 \(k\) 次迭代，根据式 \(\fml{2-1}\) 计算 \(\Delta\pmb x_k\)，并计算 \(\rho_k\)；
3. 如果
   • \(\rho_k\le0.25\)，应减小 \(\lambda_k\)，即

     \[\lambda_{k+1}=\frac{\lambda_k}2\]

   • \(0.25<\rho_k\le0.75\)，保持 \(\lambda_k\) 不变；
   • \(\rho_k>0.75\)，增大 \(\lambda_k\)，即

     \[\lambda_{k+1}=2\lambda_k\]

   • 如果 \(\rho_k\le0\)，这时不应该更新 \(\pmb x_k\)，即

     \[\pmb x_{k+1}=\pmb x_k\]

     并且和上面 \(\rho_k\le0.25\) 的情况一样，减小 \(\lambda_k\)，反之，在 \(\rho_k>0\) 的情况下，更新 \(\pmb x_k\)，即

     \[\pmb x_{k+1}=\pmb x_k\]

2.3 如何使用

还是上面的例子，我们可以使用下面的代码来拟合正弦函数。

C++

// 拟合正弦函数，初始值为 (1, 0.02, 0, 1.09)
 
rm::OptimalOptions options;
options.lsq_mode = rm::LsqMode::LM; // 使用 LM 算法
options.max_iter = 2000;            // 最大迭代次数可以设置高一点，以保证收敛
auto x = rm::lsqnonlin(lsq_sine, {1, 0.02, 0, 1.09}, options);

Python

# 拟合正弦函数，初始值为 (1, 0.02, 0, 1.09)
 
options = rm.OptimalOptions()
options.lsq_mode = rm.LsqMode.LM # 使用 LM 算法
options.max_iter = 2000          # 最大迭代次数可以设置高一点，以保证收敛
x = rm.lsqnonlin(lsq_sine, [1, 0.02, 0, 1.09], options)