从LM算法求解最小二乘到DeepLM(1)
最小二乘问题
最小二乘是一种寻找最佳拟合的数学方法，例如有一些数据点散布在坐标系中，希望寻找一条线能够穿过这些点，并且使得所有点到这条线的垂直距离最短。
具体来说，最小二乘法的核心思想分为以下几个步骤：
- 计算每个数据点到这条线的垂直距离
- 把这些距离平方后（防止正负抵消）相加
- 找到一条使距离平方和最小的线
[fig 散点和线]
求解最小二乘问题
从最基础的例子开始，假设要拟合一条直线，也就是线性最小二乘问题 $$ y=ax+b $$ $n$个数据点的误差函数为 $$ L(a,b)=\sum_1^n(y_i-(ax_i+b))^2 $$
由于这个误差函数是一个凸二次函数，可以通过求$L$对$a$和$b$的偏导并令其等于0来求最小误差：
$$ \partial L/ \partial a=-2\sum_1^n(y_i-ax_i-b)x_i=0 \\ \partial L/ \partial b=-2\sum_1^n(y_i-ax_i-b)=0 $$求解以上方程组即可得到最优的参数a和b，这个方程组称为正规方程
在实际计算中，将多个点带入公式（1）， $$ x_1a+b=y_1\\ x_2a+b=y_2\\ \cdots\\ x_na+b=y_n $$ 可以表示为矩阵形式 $$ \begin{bmatrix} x_1&1\\ x_2&1\\ \cdots&\cdots\\ x_n&1 \end{bmatrix} \begin{bmatrix} a\\ b \end{bmatrix}= \begin{bmatrix} y_1\\ y_2\\ \cdots\\ y_n \end{bmatrix} $$
简写为 $$ \mathbf{X}\bf{\beta}=\mathbf{y} $$ 这个方程组一般没有精确解（有精确解的情况下，所有的点都刚好在直线上），所以要通过最小二乘法求最优的a和b，根据公式（2），误差函数为： $$ L(\beta)=(\mathbf{X}\bf{\beta}-\mathbf{y})^\top(\mathbf{X}\bf{\beta}-\mathbf{y}) $$ 通过求导得到正规方程 $$ \mathbf{X}^\top\mathbf{X}\bf{\beta}=\mathbf{X}^\top\mathbf{y} $$ 所以最终的解为 $$ \bf{\beta}=(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y} $$ 更加复杂的情况下，需要拟合的点可能不是一条直线，而是曲线，也就是非线性最小二乘问题 $$ y=f(x,\beta) $$ 其中$f$为非线性函数
非线性问题没有办法用正规方程求解，因为非线性问题的正规函数也是非线性的，例如，要求解指数函数的最小二乘问题 $$ y=\beta_1exp(\beta_2x) $$ 目标函数为： $$ L(\beta_1,\beta_2)=\sum(y_i-\beta_1exp(\beta_2x_i))^2 $$ 求偏导得 $$ \partial L/\partial\beta_1=-2\sum(y_i-\beta_1exp(\beta_2x_i))exp(\beta_2x_i)=0\\ \partial L/\partial\beta_2=-2\sum(y_i-\beta_1exp(\beta_2x_i))\beta_1x_iexp(\beta_2x_i)=0 $$ 得到的非线性的正规函数没有办法直接求解，因此需要换一种方法，通过泰勒展开用线性近似代替原来的非线性函数
假设当前是第k次迭代，表示为上标$(k)$
公式（12）所表示的指数函数最小二乘目标函数，在点$\beta_1^{(k)},\beta_2^{(k)}$处的一阶泰勒展开为 $$ f(x;\beta)\approx f(x;\beta^{(k)})+(\partial f/\partial \beta_1)(\beta_1-\beta_1^{(k)})+(\partial f/\partial\beta_2)(\beta_2-\beta_2^{(k)}) $$ 其中$f(x;\beta^{(k)})=\beta_1^{(k)}exp(\beta_2^{(k)}x)$，$\partial f/\partial\beta_1=exp(\beta_2^{(k)}x)$，$\partial f/\partial\beta_2=\beta_1^{(k)}xexp(\beta_2^{(k)}x)$，代入得 $$ f(x_i;\beta) \approx \beta_1^{(k)}exp(\beta_2^{(k)}x_i) + exp(\beta_2^{(k)}x_i)(\beta_1-\beta_1^{(k)}) + \beta_1^{(k)}x_iexp(\beta_2^{(k)}x_i)(\beta_2-\beta_2^{(k)}) $$ 代入残差表达式$r_i=y_i-\beta_1exp(\beta_2x_i)$，得 $$ r_i = y_i - [\beta_1^{(k)}exp(\beta_2^{(k)}x_i) + exp(\beta_2^{(k)}x_i)(\beta_1-\beta_1^{(k)}) + \beta_1^{(k)}x_iexp(\beta_2^{(k)}x_i)(\beta_2-\beta_2^{(k)})] $$ 将参数的差记为增量的形式 $$ \Delta\beta_1 = \beta_1-\beta_1^{(k)}\\ \Delta\beta_2 = \beta_2-\beta_2^{(k)} $$ 整理后得到 $$ r_i = y_i - \beta_1^{(k)}exp(\beta_2^{(k)}x_i) - exp(\beta_2^{(k)}x_i)\Delta\beta_1 - \beta_1^{(k)}x_iexp(\beta_2^{(k)}x_i)\Delta\beta_2 $$
$$ r_i = r_i^{(k)} - [exp(\beta_2^{(k)}x_i)\Delta\beta_1 + \beta_1^{(k)}x_iexp(\beta_2^{(k)}x_i)\Delta\beta_2] $$记Jacobian矩阵的第i行为： $$ J_i = [-exp(\beta_2^{(k)}x_i), -\beta_1^{(k)}x_iexp(\beta_2^{(k)}x_i)] $$ 参数的增量为： $$ \Delta\beta = [\Delta\beta_1, \Delta\beta_2]^T $$ 则可以写成矩阵形式： $$ \begin{bmatrix} r_1 \ r_2 \ \vdots \ r_n \end{bmatrix} \approx \begin{bmatrix} r_1^{(k)} \ r_2^{(k)} \ \vdots \ r_n^{(k)}
\end{bmatrix}
\begin{bmatrix} exp(\beta_2^{(k)}x_1) & \beta_1^{(k)}x_1exp(\beta_2^{(k)}x_1) \\ exp(\beta_2^{(k)}x_2) & \beta_1^{(k)}x_2exp(\beta_2^{(k)}x_2) \\ \vdots & \vdots \\ exp(\beta_2^{(k)}x_n) & \beta_1^{(k)}x_nexp(\beta_2^{(k)}x_n) \end{bmatrix}\begin{bmatrix} \Delta\beta_1 \\ \Delta\beta_2 \end{bmatrix}$$ 简写为： $$r \approx r^{(k)} - J\Delta\beta $$ 此时通过泰勒展开的线性化，原本的最小二乘问题变为： $$ \min_{\Delta\beta} |r^{(k)} - J\Delta\beta|^2 $$ 由非线性最小二乘问题，转换成了线性最小二乘问题，即可通过求解正规方程的方式求解，正规方程为： $$ (J^TJ)\Delta\beta = J^Tr^{(k)} $$ 其中$r_i^{(k)} = y_i - \beta_1^{(k)}exp(\beta_2^{(k)}x_i)$，解得： $$ \Delta\beta = (J^TJ)^{-1}J^Tr^{(k)} $$ 由于泰勒展开只是在当前点$\beta^{(k)}$附近的局部近似，只在$\beta^{(k)}$附近具有好的近似效果，因此，只能在当前点的附近使用，如果直接更新参数 $$ \beta^{(k+1)} = \beta^{(k)} + \Delta\beta $$ 则违反了只能在当前点**附近**使用的规则，因此，需要使用一个步长因子来控制更新的方式，熟悉梯度下降的同学应该对步长这个概念不陌生，添加步长后的更新方式： $$ \beta^{(k+1)} = \beta^{(k)} + \alpha\Delta\beta $$ 其中$\alpha$为步长因子，$0<\alpha<1$
在更新参数$\beta$后，再重新从公式（25）开始计算新的参数更新量，经过多次迭代后得到一个较为准确的解
Levenberg-Marquardt方法
Levenberg-Marquardt(LM)方法中，修改了公式（25）中的正规方程，修改后的方程为： $$ (J^TJ + \lambda I)\Delta\beta = J^Tr^{(k)} $$ 多了一个阻尼参数$\lambda$，$I$为单位矩阵
当$\lambda$很大时，$\lambda I$项占主导地位，因此公式（29）可以近似为： $$ \lambda I\Delta\beta \approx J^Tr^{(k)} $$ 此时的解为： $$ \Delta\beta \approx \frac{1}{\lambda}J^Tr^{(k)} $$ 目标函数（12）关于参数$\beta$的梯度为 $$ \nabla L = \frac{\partial L}{\partial \beta} = -2\sum_i r_i\frac{\partial f}{\partial \beta} = -2J^Tr $$ 而公式（31）可以写为： $$ \Delta\beta \approx -\frac{1}{2\lambda}(-2J^Tr^{(k)}) = -\frac{1}{2\lambda}\nabla L $$ 这时的形式符合梯度下降： $$ \beta^{(k+1)} = \beta^{(k)} - \alpha\nabla L $$ 可以将公式（33）看做步长为$1/2\lambda$的梯度下降，所以在$\lambda$很大时，LM方法类似于梯度下降
当$\lambda$很小时，$J^TJ$占主导地位 $$ J^TJ\Delta\beta = J^Tr^{(k)} $$
$$ \Delta\beta = (J^TJ)^{-1}J^Tr^{(k)} $$此时的更新方程实际上就是Gauss-Newton方法的更新方程，在这里不赘述
总结来说就是，当$\lambda$很大时，步长与$\lambda$成反比，下降方向通过一阶导数确定，接近最速下降方向 $$ \|\Delta\beta\| \approx \frac{\|J^Tr\|}{\lambda} $$ 当$\lambda$很小时，步长由二阶导数的近似$(J^TJ)^{-1}$确定 $$ \|\Delta\beta\| \approx \|(J^TJ)^{-1}J^Tr\| $$ 在实际的计算中，初始计算时通常使用较大的$\lambda$，随着迭代的进行，逐渐减小$\lambda$，可以让计算过程始终保持稳定