[机器学习] [三维重建]最小二乘法

一、线性最小二乘法

假设平面中存在一组数据点，需要找到一个数学模型（在这里是一条直线）去拟合这组数据点，这就是拟合。

而实际上，样本点的维数不会仅仅是2维的，而可以是任意的$p$维。

为了不失一般性，假设有一组样本：
D = { ( x 1 , y 1 ) , ( x 2 , y 2 ) , ( x 3 , y 3 ) , . . . , ( x n , y n ) } D=\{(x_1,y_1),(x_2,y_2),(x_3,y_3),...,(x_n,y_n)\} D={(x1​,y1​),(x2​,y2​),(x3​,y3​),...,(xn​,yn​)}
$x_i$ 是一个 $p$ 维的向量，用于表达第 $i$ 个样本被观察的 $p$ 个特征，而 $y_i$ 表示第 $i$ 个样本的取值，是一个1维的数值。线性最小二乘法拟合出来的直线也就是样本的所有 $p$ 维特征值到最终样本取值 $y_i$的一个近似的线性映射关系。

这里需要明确一下适用范围，只有当所要拟合的数学模型（$y_i$ 与 $x_i$之间可能的关系） 是线性函数的关系，才能适用线性最小二乘法。

线性函数：

从数学角度来说，对于一个函数 $y=f(x)$ 如果它是线性函数，需要满足两个法则：

（1）比例性：对于任意常数 $a$，都有 $ay=f(ax)$ 成立，即将 $x$ 扩大 $a$ 倍，$y$ 也将随之扩大 $a$ 倍。

（2）叠加性：若$y_1=f(x_1)，y_2=f(x_2)$则$y_1+y_2=f(x_1+x_2)$

但这么看来，只有$y=wx$属于线性函数，而$y=wx+b$ 不属于数学角度上的线性函数，但这里暂且把深度学习领域与数学领域分开，在深度学习领域 $y=wx+b$也属于线性函数。

回到最小二乘法，对于每个样本 $x_i$ 都可以写成一个 $p$ 维的向量：
$$x_i=\left[\begin{array}{c}{x}_{i1}\\ x_{i2}\\ x_{i3}\\ ...\\ x_{ip}\end{array}\right]$$
因此最小二乘法的目标就是估计出一个系数向量：
$$w=\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\\ ...\\ w_p\end{array}\right]$$
和一个偏置常数 $b$
来建立如下线性关系：
$$w^{T}x+b\Rightarrow y$$
为简化描述，另$b=w_0{x}_{i0}$，其中$x_{i0}$恒为1，则上述两个向量更新为：
$$x_i=\left[\begin{array}{c}1\\ x_{i1}\\ x_{i2}\\ x_{i3}\\ ...\\ x_{ip}\end{array}\right]$$
$$w=\left[\begin{array}{c}{w}_0\\ w_1\\ w_2\\ w_3\\ ...\\ w_p\end{array}\right]$$
映射关系简化成：
$$w^{T}x\Rightarrow y$$
最小二乘法就是估计参数向量 $w$ 的一种方法，在最小二乘法中，定义了如下的一个目标函数：
$$L(w)=\sum _{i=1}^N|w^T{x}_i-y_i{|}^2$$

可以看出，目标函数$L(w)$是以参数向量 $w$ 为自变量的线性函数。

其中，$w^T{x}_i$意味着拟合得到的数学模型针对输入值$x_i$得到的输出值，而$y_i$意味着真实情况下，$x_i$对应的真实值，也就是实际的观测值。找到一个系数向量$w$，也就是找到一个数学模型，使得$L(w)$的取值最小，也就是模型的输出值和真实值之间的误差的平方和最小。

具体计算步骤：
$$L(w)=\sum _{i=1}^N|w^T{x}_i-y_i{|}^2$$
$$=\left[\begin{array}{cccc}{w}^T{x}_1-y_1& w^T{x}_2-y_2& ...& w^T{x}_N-y_N\end{array}\right]\left[\begin{array}{c}{w}^T{x}_1-y_1\\ w^T{x}_2-y_2\\ ...\\ w^T{x}_N-y_N\end{array}\right]$$
对于第一块向量：
$$\left[\begin{array}{cccc}{w}^T{x}_1-y_1& w^T{x}_2-y_2& ...& w^T{x}_N-y_N\end{array}\right]$$
$$=w^T\left[\begin{array}{cccc}{x}_1& x_2& ...& x_N\end{array}\right]-\left[\begin{array}{cccc}{y}_1& y_2& ...& y_N\end{array}\right]$$
$$=w^T{X}^T-Y^T$$
$$=\left[\begin{array}{cccccc}{w}_0& w_1& w_2& w_3& ...& w_p\end{array}\right]\left[\begin{array}{ccccc}1& 1& 1& ...& 1\\ x_{11}& x_{21}& x_{31}& ...& x_{N1}\\ x_{12}& x_{22}& x_{32}& ...& x_{N2}\\ ...& & & & \\ x_{1p}& x_{2p}& x_{3p}& ...& x_{Np}\end{array}\right]–\left[\begin{array}{ccccc}{y}_1& y_2& y_3& ...& y_N\end{array}\right]$$
其中，$X$为N个样本集合：
$$X=\left[\begin{array}{c}{x}_1^T\\ x_2^T\\ x_3^T\\ ...\\ x_N^T\end{array}\right]=\left[\begin{array}{cccccc}1& x_{11}& x_{12}& x_{13}& ...& x_{1p}\\ 1& x_{21}& x_{22}& x_{23}& ...& x_{2p}\\ 1& x_{31}& x_{32}& x_{33}& ...& x_{3p}\\ ...& & & & & \\ 1& x_{N1}& x_{N2}& x_{N3}& ...& x_{Np}\end{array}\right]$$
$Y$为这N个样本对应的观测值的集合：
$$Y=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ ...\\ y_N\end{array}\right]$$
而后面一块向量是前面一块的转置，因此：
$$\left[\begin{array}{c}{w}^T{x}_1-y_1\\ w^T{x}_2-y_2\\ ...\\ w^T{x}_N-y_N\end{array}\right]={\left[\begin{array}{cccc}{w}^T{x}_1-y_1& w^T{x}_2-y_2& ...& w^T{x}_N-y_N\end{array}\right]}^T=(w^T{X}^T-Y^T{)}^T=Xw-Y$$
最终得到：
$$L(w)=(w^T{X}^T-Y^T)(Xw-Y)$$
化简得：
$$L(w)=w^T{X}^{T}Xw-w^T{X}^{T}Y-Y^{T}Xw+Y^{T}Y$$
由于L(w)的最终结果是一个标量（实数），所以式子中的每一项必然也是一个实数，式子内$w^T{X}^{T}Y$和$Y^{T}Xw$是互为转置的关系，又由于两者都是实数，所以可知两者相等，于是式子进一步化简：
$$L(w)=w^T{X}^{T}Xw-2w^T{X}^{T}Y+Y^{T}Y$$
L(w)对w求偏导$\frac{\partial }{\partial w}L(w)$，当$\frac{\partial }{\partial w}L(w)=0$时，L(w)取得极值，可能是全局最大值、全局最小值或局部极大值、极小值，如果存在多个值，只要逐个比较它们的函数值的大小即可，最终可以寻找到全局最小值（在这里可能不需要，这个目标函数似乎是个凸函数，极值便是全局最小值）：
$$\frac{\partial }{\partial w}L(w)=\frac{\partial }{\partial w}(w^T{X}^{T}Xw-2w^T{X}^{T}Y+Y^{T}Y)=2X^{T}Xw-2X^{T}Y=0$$
得：
$$2X^{T}Xw-2X^{T}Y=0$$
$$X^{T}Xw=X^{T}Y$$
最终得到用于估计参数向量w的表达式：
$$w=(X^{T}X{)}^{-1}{X}^{T}Y$$
注意到，最后一步中，对$X^{T}X$有一个取逆的步骤，因此，存在潜在假设：

$X^{T}X$是非奇异矩阵，即$X^{T}X$可逆，即$X^{T}X$的行列式不为零

但$X^{T}X$在以下情况会是不可逆的：
（1）当样本数量 $N$ 小于参数向量 $w$ 的维度的时候

（2）在所有特征中，存在一个特征与另一个特征线性相关或一个特征与若干个特征线性相关的时候

具体来说，矩阵 $X$ 是 $N\times p$ 维的矩阵，若存在线性相关的特征，即 $p$ 个特征之间不是线性无关，即列不满秩：$R(X)<p，R(X^T)<p,R(X^{T}X)<p$，所以$X^{T}X$不可逆
$$X=\left[\begin{array}{c}{x}_1^T\\ x_2^T\\ x_3^T\\ ...\\ x_N^T\end{array}\right]=\left[\begin{array}{cccccc}1& x_{11}& x_{12}& x_{13}& ...& x_{1p}\\ 1& x_{21}& x_{22}& x_{23}& ...& x_{2p}\\ 1& x_{31}& x_{32}& x_{33}& ...& x_{3p}\\ ...& & & & & \\ 1& x_{N1}& x_{N2}& x_{N3}& ...& x_{Np}\end{array}\right]$$

因此，当出现不可逆的情况时，有如下解决方法：
（1）筛选出线性无关的特征，保证特征之间是线性无关的
（2）增加样本数量，保证样本数量大于特征维度（参数向量w的维度）（样本数据的维度）
（3）采用正则化的方法

二、非线性最小二乘法

1.一阶二阶泰勒展开

当所要拟合的数学模型是非线性的情况时，便不再适用非线性最小二乘法。

假设参数向量 $w$ （维度仍然是 $p$ 维）是某个标量非线性数学模型 $f(x)$ 的参数：
$$f(x):x_i\in R^p\Rightarrow y_i\in R$$
与线性最小二乘法类似，需要一个描述模型预测值与真实观测值之间的差距的函数：
$$f(x_i)-y_i$$
需要注意的是，$x_i$ 是样本数据，$y_i$ 是样本取值，观测结果。于是，在导入数据后，$x_i$ 和 $y_i$ 是已知的，而数学模型 $f(x)$ 内包含的参数 $w$ 才是未知的，需要求解的。

因此，上述函数看作是自变量为 $w$ 的函数：
$$r_i(w)$$
目标函数：
$$L(w)=\sum _{i=1}^N{r}_i^2(w)$$

注：每个$r_i(w)$是不一样的，在这个以 $w$ 为自变量的函数中，它的参数是已知数据 $x_i$ 和 $y_i$ ，每次收集的数据都是不一样的，因此每个$r_i(w)$是不同的函数。

从线性最小二乘法的步骤可以知道，通过令目标函数的导数为零
$$\frac{dL(w)}{dw}=0$$
求解上述方程，便能得到使目标函数最小的最优值 $w$

这个方程是否容易求解，取决于$L(w)$的导函数的形式，如果$L(w)$为简单的线性函数，那么这个问题就是线性最小二乘问题。
但是前面已说，$L(w)$并不是线性函数，这使得该方程可能不容易求解。

求解这个方程需要我们知道关于目标函数的全局性质，这在非线性情况下通常是不太可能的，对于此类非线性问题，可以使用迭代的方式，从一个初始值出发，不断更新当前的需要优化的变量，使目标函数逐步下降，具体步骤如下：

1.给定某个初始值$w_0$
2.对于第 $k$ 次迭代，寻找一个增量$\Delta w_k$，使得$L(w_k+\Delta w_k)<L(w_k)$，即目标函数下降
3.若$\Delta w_k$足够小，则停止迭代
4.否则，令$w_{k+1}=w_k+\Delta w_k$，返回第二步

这使得求解导函数为零的问题（全局问题）变成了不断寻找下降增量$\Delta w_k$的问题（局部问题），从而能只关心$L(w)$在迭代值处的局部性质，因此能引入微分的思想：
$$非线性函数在某个值附近的一个小邻域内，可以视作线性的$$

具体来说，考虑第 $k$ 次迭代，假设在 $w_k$ 处想寻找到一个增量 $\Delta w_k$，使得目标函数 $L(w)$ 下降，通过将目标函数在 $w_k$ 附近作泰勒展开，近似成线性：
$$L(w_k+\Delta w_k)\approx L(w_k)+J(w_k{)}_{1\times p}^T(\Delta w_k{)}_{p\times 1}+\frac{1}{2}\Delta w_k^{T}H(w_k)\Delta w_k\text{\hspace{1em}(1)}$$
其中，$J(w_k)$（特殊情况下的雅可比矩阵、梯度向量、梯度）是 $L(w)$ 关于 $w$ 的一阶导数，$H(w_k)$ （海森矩阵）则是二阶导数，它们需要在 $w_k$ 处取值。属于多元微积分课程的知识。贴一下相关概念的表述。

关于式（1）到底是怎么来的：

目标函数：
$$L(w)=\sum _{i=1}^N{r}_i^2(w)=r_1^2(w)+r_2^2(w)+...+r_N^2(w)$$
引入函数向量：
$$r(w)={\left[\begin{array}{c}{r}_1(w)\\ r_2(w)\\ ...\\ r_N(w)\end{array}\right]}_{N\times 1}$$
$$L(w)=r(w{)}^{T}r(w)$$
目标函数显然是个关于自变量 $w$ 的多元函数（$w$是个p维参数向量），它在 $w_k$ 的泰勒展开式：
$$L(w_k+\Delta w_k)\approx L(w_k)+\nabla L(w_k{)}^T\Delta w_k+\frac{1}{2}\Delta w_k^{T}H(w_k)\Delta w_k\text{\hspace{1em}(2)}$$
$$\nabla L(w_k)={\left[\begin{array}{c}\frac{\partial L}{w_1}\\ \\ \frac{\partial L}{w_2}\\ \\ ...\\ \\ \frac{\partial L}{w_p}\end{array}\right]}_{p\times 1}$$
其中：
$$\frac{\partial L}{w_i}=\frac{\partial L}{r_1}\frac{\partial r_1}{w_i}+\frac{\partial L}{r_2}\frac{\partial r_2}{w_i}+...+\frac{\partial L}{r_N}\frac{\partial r_N}{w_i}$$
$$=\left[\begin{array}{cccc}\frac{\partial r_1}{w_i}& \frac{\partial r_2}{w_i}& ...& \frac{\partial r_N}{w_i}\end{array}\right]\left[\begin{array}{c}\frac{\partial L}{r_1}\\ \frac{\partial L}{r_2}\\ ...\\ \frac{\partial L}{r_N}\end{array}\right]$$
则：
$$\nabla L(w_k)=\left[\begin{array}{cccc}\frac{\partial r_1}{w_1}& \frac{\partial r_2}{w_1}& ...& \frac{\partial r_N}{w_1}\\ \frac{\partial r_1}{w_2}& \frac{\partial r_2}{w_2}& ...& \frac{\partial r_N}{w_2}\\ ...& ...& ...& ...\\ \frac{\partial r_1}{w_p}& \frac{\partial r_2}{w_p}& ...& \frac{\partial r_N}{w_p}\end{array}\right]\left[\begin{array}{c}\frac{\partial L}{r_1}\\ \frac{\partial L}{r_2}\\ ...\\ \frac{\partial L}{r_N}\end{array}\right]=\left[\begin{array}{cccc}\frac{\partial r_1}{w_1}& \frac{\partial r_2}{w_1}& ...& \frac{\partial r_N}{w_1}\\ \frac{\partial r_1}{w_2}& \frac{\partial r_2}{w_2}& ...& \frac{\partial r_N}{w_2}\\ ...& ...& ...& ...\\ \frac{\partial r_1}{w_p}& \frac{\partial r_2}{w_p}& ...& \frac{\partial r_N}{w_p}\end{array}\right]\left[\begin{array}{c}2r_1\\ 2r_2\\ ...\\ 2r_N\end{array}\right]$$
$$=2J(w{)}^{T}r(w)$$
2的影响不大，有时候目标函数前会特意设置个$\frac{1}{2}$，目的就是抵消此处的2

最终结果：
$$\nabla L(w_k)=J(w{)}^{T}r(w)\text{\hspace{1em}(3)}$$
其中，$J(w)$就是作为中间变量被引入的函数向量 $r(w)$
$$r(w)=\left[\begin{array}{c}{r}_1(w)\\ r_2(w)\\ ...\\ r_N(w)\end{array}\right]$$
关于自变量参数向量 $w$ 的雅可比矩阵。
$r(w)$ 就是引入的函数向量。

将式（3）代入式（2），最终得到目标函数的泰勒展开式：
$$L(w_k+\Delta w_k)\approx L(w_k)+r(w{)}_{1\times N}^{T}J(w{)}_{N\times p}(\Delta w_k{)}_{p\times 1}+\frac{1}{2}\Delta w_k^{T}H(w_k)\Delta w_k\text{\hspace{1em}(4)}$$
与式（1）对比发现，式（1）中的一次项内并没有 $r(w)$ ，这是因为推导得到式（4）时引用了函数变量 $r(w)$ ，从而使得符合雅可比矩阵的适用情况，当没有引入函数变量 $r(w)$ 的时候，$J(w)$（或者说$\nabla L(w_k)$） 实际上就是目标函数$L(w)$ 关于自变量参数向量 $w$ 的梯度向量：
$$\nabla L(w_k)=J(w)=\left[\begin{array}{c}\frac{\partial L}{w_1}\\ \\ \frac{\partial L}{w_2}\\ \\ ...\\ \\ \frac{\partial L}{w_p}\end{array}\right]$$
简单来说就是，梯度向量是雅可比矩阵的特殊情况。所以式（1）是没有引入函数向量$r(w)$时推导得到的式子，而式（4）更具一般性，是引入函数向量 $r(2)$ 后推导得到的式子。两个式子中，$J(w)$ 的意义不同。

注意：

（1）留意文中向量的横竖，全文中向量都是竖着的列向量，导致某些公式的转置符号可能与网上别的文章的转置符号略有区别。

（2）留意文中矩阵的行列维数，多元函数泰勒展开式的矩阵形式，实际上是把累加的部分写成了矩阵形式

2.$\Delta w_k$的求解

处理式（1）：
$$L(w_k+\Delta w_k)\approx L(w_k)+J(w_k{)}^T\Delta w_k+\frac{1}{2}\Delta w_k^{T}H(w_k)\Delta w_k\text{\hspace{1em}(1)}$$
在使用泰勒展开后，可以将非线性的目标函数在某个邻域内采用线性多项式逼近，对这个线性多项式的不同处理（如保留一阶项或保留二阶项），取得增量$\Delta w_k$的方向，再指定一个步长 $\lambda$ ，从而保证目标函数下降。求解$\Delta w_k$的不同处理方式，就是非线性最小二乘法的不同方法。

（1）最速下降法

已知梯度方向是函数值增长最快的方向，那么梯度方向的反方向自然是函数值下降最快的方向，如果取增量为反向的梯度，就一定可以保证函数下降（注：梯度即多元函数的一阶导数，即这里的雅可比矩阵$J(w_k)$ ）：
$$\Delta w_k=-J(w_k)$$
由于讨论都是在第 $k$ 次迭代进行，并不涉及其他迭代信息，为了简化符号，省略下标 $k$：
$$\Delta w=-J(w)$$

（2）牛顿法

最速下降法只保留了泰勒展开的一阶梯度信息，为了更加逼近非线性函数，也可以保留至泰勒展开的二阶梯度信息：
$$L(w+\Delta w)\approx L(w)+J(w{)}^T\Delta w+\frac{1}{2}\Delta w^{T}H(w)\Delta w$$
以处理线性多项式的方式处理上式，对 $\Delta w$ 求一阶导数，并令结果为0（对线性多项式来说，一阶导数为0处便是极值所在处，也就是最小值可能所在处，如果有多个结果，对多个结果进行比较，最小的那个就是最小值所在处）：
$$J+H\Delta w=0$$
解得增量$\Delta w$:
$$\Delta w=H^{-1}(-J)$$

上述的一阶和二阶梯度法都比较直观，利用泰勒展开在迭代点附近用线性的一次或二次多项式近似了目标函数，然后用该线性多项式最小值来猜测目标函数的最小值，这种方法只要目标函数在局部微小区域内是线性的，就能成立。由于微分思想保证非线性函数的局部微小区域确实可以采用线性去逼近，所以这种算法确实是可行的。

但上述两个方法存在一定的缺点，最速下降法过于贪心，在到达最小值的过程中常常走出锯齿路线，增加了迭代次数。

而牛顿法需要计算目标函数的 $H$ 矩阵，这通常是比较麻烦的，实际上人们更倾向于避免 $H$ 矩阵的计算，所以提出了几种更实用的方法（高斯牛顿法和列文伯格—马夸尔特法）。

（3）高斯牛顿法

最小二乘法的目标函数（为方便计算，消去求导后可能出现的2，所以前面添加了个系数$\frac{1}{2}$）：
$$L(w)=\frac{1}{2}\sum _{i=1}^N{r}_i^2(w)$$
记得之前引入的
$$r(w)=\left[\begin{array}{c}{r}_1(w)\\ r_2(w)\\ ...\\ r_N(w)\end{array}\right]$$后续会用到。

之前的方法都是针对目标函数$L(w)$进行泰勒展开，而这个目标函数是多个非线性函数 $r_i(w)$ 的平方的累加和，高斯牛顿的思想是先针对 $r_i(w)$
进行一阶的泰勒展开：
$$r_i(w+\Delta w)\approx r_i(w)+J_i(w{)}^T\Delta w$$
$J_i(w)$ 是 非线性函数$r_i(w)$ 关于 自变量参数向量 $w$ 的雅可比矩阵，将经过一阶泰勒展开后的 $r_i(w)$ 代入原目标函数，得到的目标函数如下：
$$L(w+\Delta w)=\sum _{i=1}^N\frac{1}{2}(r_i(w)+J_i(w{)}^T\Delta w{)}^2\text{\hspace{1em}(5)}$$
问题转变成了：在此次迭代中（别忘了迭代仍然存在，只是没写下标 $k$而已），寻找增量 $\Delta w$ ，使式（5）在此次迭代中达到最小值。

将式（5）的平方展开，并对 $\Delta w$ 求一阶导数，并令导数等于 0 ，便能取到使该式子达到最小值时的$\Delta w$ 的取值。

先将目标函数的平方展开：
$$L(w+\Delta w)=\sum _{i=1}^N\frac{1}{2}[r_i(w)+J_i(w{)}^T\Delta w){]}^T[r_i(w)+J_i(w{)}^T\Delta w)]$$
$$=\sum _{i=1}^N\frac{1}{2}[r_i^2(w)+2r_i(w)J_i(w{)}^T\Delta w+\Delta w^T{J}_i(w)J_i(w{)}^T\Delta w]$$
求上式对 $\Delta w$ 的导数，并令其为零：
$$\sum _{i=1}^N{J}_i(w)r_i(w)+\sum _{i=1}^N{J}_i(w)J_i(w{)}^T\Delta w=0$$
将累加形式写成矩阵形式：

第一个累加项：
$$J_i(w)=\left[\begin{array}{c}\frac{\partial r_i}{\partial w1}\\ \\ \frac{\partial r_i}{\partial w2}\\ \\ ...\\ \\ \frac{\partial r_i}{\partial wp}\end{array}\right]$$
$r_i(w)$ 是第 $i$ 个样本数据得到的输出值与真实值之间的误差 $error$，是个标量，是一个数。
$$J_i(w)r_i(w)=\left[\begin{array}{c}{r}_i(w)\frac{\partial r_i}{\partial w1}\\ \\ r_i(w)\frac{\partial r_i}{\partial w2}\\ \\ ...\\ \\ r_i(w)\frac{\partial r_i}{\partial wp}\end{array}\right]$$
$$\sum _{i=1}^N{J}_i(w)r_i(w)=\left[\begin{array}{c}{\sum }_{i=1}^N{r}_i(w)\frac{\partial r_i}{\partial w1}\\ \\ {\sum }_{i=1}^N{r}_i(w)\frac{\partial r_i}{\partial w2}\\ \\ ...\\ \\ {\sum }_{i=1}^N{r}_i(w)\frac{\partial r_i}{\partial wp}\end{array}\right]$$
$$=\left[\begin{array}{cccc}\frac{\partial r_1}{\partial w1}& \frac{\partial r_2}{\partial w1}& ...& \frac{\partial r_N}{\partial w1}\\ & & & \\ \frac{\partial r_1}{\partial w2}& \frac{\partial r_2}{\partial w2}& ...& \frac{\partial r_N}{\partial w2}\\ & & & \\ ...& & & \\ & & & \\ \frac{\partial r_1}{\partial wp}& \frac{\partial r_2}{\partial wp}& ...& \frac{\partial r_N}{\partial wp}\end{array}\right]\left[\begin{array}{c}{r}_1(w)\\ \\ r_2(w)\\ \\ ...\\ \\ r_N(w)\end{array}\right]$$
$$=J(w{)}^{T}r(w)$$
$r(w)$ 就是上面提到的引入的函数向量，$J(w)$ 就是这个函数向量对自变量 $p$ 维参数向量 $w$ 的雅克比矩阵。

第二个累加项的系数：
$$J_i(w)J_i(w{)}^T=\left[\begin{array}{c}\frac{\partial r_i}{\partial w1}\\ \\ \frac{\partial r_i}{\partial w2}\\ \\ ...\\ \\ \frac{\partial r_i}{\partial wp}\end{array}\right]\left[\begin{array}{cccc}\frac{\partial r_i}{\partial w1}& \frac{\partial r_i}{\partial w2}& ...& \frac{\partial r_i}{\partial wp}\end{array}\right]$$
$$=\left[\begin{array}{cccc}\frac{\partial r_i}{\partial w1}\frac{\partial r_i}{\partial w1}& \frac{\partial r_i}{\partial w1}\frac{\partial r_i}{\partial w2}& ...& \frac{\partial r_i}{\partial w1}\frac{\partial r_i}{\partial wp}\\ & & & \\ \frac{\partial r_i}{\partial w2}\frac{\partial r_i}{\partial w1}& \frac{\partial r_i}{\partial w2}\frac{\partial r_i}{\partial w2}& ...& \frac{\partial r_i}{\partial w2}\frac{\partial r_i}{\partial wp}\\ ...& & & \\ \frac{\partial r_i}{\partial wp}\frac{\partial r_i}{\partial w1}& \frac{\partial r_i}{\partial wp}\frac{\partial r_i}{\partial w2}& ...& \frac{\partial r_i}{\partial wp}\frac{\partial r_i}{\partial wp}\end{array}\right]$$
$$\sum _{i=1}^N{J}_i(w)J_i(w{)}^T==\left[\begin{array}{cccc}{\sum }_{i=1}^N\frac{\partial r_i}{\partial w1}\frac{\partial r_i}{\partial w1}& {\sum }_{i=1}^N\frac{\partial r_i}{\partial w1}\frac{\partial r_i}{\partial w2}& ...& {\sum }_{i=1}^N\frac{\partial r_i}{\partial w1}\frac{\partial r_i}{\partial wp}\\ & & & \\ {\sum }_{i=1}^N\frac{\partial r_i}{\partial w2}\frac{\partial r_i}{\partial w1}& {\sum }_{i=1}^N\frac{\partial r_i}{\partial w2}\frac{\partial r_i}{\partial w2}& ...& {\sum }_{i=1}^N\frac{\partial r_i}{\partial w2}\frac{\partial r_i}{\partial wp}\\ ...& & & \\ {\sum }_{i=1}^N\frac{\partial r_i}{\partial wp}\frac{\partial r_i}{\partial w1}& {\sum }_{i=1}^N\frac{\partial r_i}{\partial wp}\frac{\partial r_i}{\partial w2}& ...& {\sum }_{i=1}^N\frac{\partial r_i}{\partial wp}\frac{\partial r_i}{\partial wp}\end{array}\right]$$
$$=\left[\begin{array}{cccc}\frac{\partial r_1}{\partial w1}& \frac{\partial r_2}{\partial w1}& ...& \frac{\partial r_N}{\partial w1}\\ & & & \\ \frac{\partial r_1}{\partial w2}& \frac{\partial r_2}{\partial w2}& ...& \frac{\partial r_N}{\partial w2}\\ & & & \\ ...& & & \\ & & & \\ \frac{\partial r_1}{\partial wp}& \frac{\partial r_2}{\partial wp}& ...& \frac{\partial r_N}{\partial wp}\end{array}\right]\left[\begin{array}{cccc}\frac{\partial r_1}{\partial w1}& \frac{\partial r_1}{\partial w2}& ...& \frac{\partial r_1}{\partial wp}\\ & & & \\ \frac{\partial r_2}{\partial w1}& \frac{\partial r_2}{\partial w2}& ...& \frac{\partial r_2}{\partial wp}\\ & & & \\ ...& & & \\ & & & \\ \frac{\partial r_N}{\partial w1}& \frac{\partial r_N}{\partial w2}& ...& \frac{\partial r_N}{\partial wp}\end{array}\right]$$
$$=J(w{)}^{T}J(w)$$
这段累加很繁琐，注意各个矩阵的元素以及各个矩阵的维度，对具体代码实现来说很重要。

最终可以得到如下方程组：
$$J(w{)}^{T}r(w)+J(w{)}^{T}J(w)\Delta w=0$$
其中，$J(w)$ 强调过很多次了，是引入的函数向量 $r(w)$ 针对自变量参数向量 $w$ 的雅克比矩阵，$r(w)$ 是引入的函数变量。

转换得到如下方程组：
$$J(w{)}^{T}J(w)\Delta w=-J(w{)}^{T}r(w)\text{\hspace{1em}(6)}$$
这个方程组是关于增量 $\Delta w$ 的线性方程组，称作增量方程，或高斯牛顿方程，或正规方程。把左边的系数定义为 $H$ ，右边定义为 $g$，得到如下线性方程组形式：
$$H\Delta w=g$$
求解上述增量方程可以解出增量 $\Delta w$
$$\Delta w=-(J^{T}J{)}^{-1}{J}^{T}r(w)$$

这里把左侧系数记作 $H$ 是有特殊意义的，对比牛顿法可知，高斯牛顿法用 $J(w{)}^{T}J(w)$ 作为牛顿法中二阶海森矩阵的近似，从而省略了计算海森矩阵的步骤。求解线性方程组式（6）是整个算法的核心所在。

注意到求解过程中，有对矩阵 $J^{T}J$ 取逆的操作，这就要求矩阵 $J^{T}J$ 是可逆矩阵。但在实际情况中，可能出现 $J^{T}J$ 为奇异矩阵（不满秩，即不可逆）的情况，直观来说，这时原函数在这个点的局部近似并不像是一个二次函数。列文伯格——马夸尔特方法在一定程度上解决了这个问题（$J^{T}J$ 并不总是可逆），一般认为它比高斯牛顿法更健壮，但收敛速度可能会比高斯牛顿法慢，被称为阻尼牛顿法。

（4）列文伯格——马夸尔特法（LM法）

高斯牛顿法中采用的近似二阶泰勒泰勒展开的方式只在展开点附近一定范围内有较好的近似效果，当增量$\Delta w$ 过大时，就有可能超出这个范围，因此自然想到给增量 $\Delta w$ 添加一个范围，称为信赖区域。这个范围定义了在什么情况下二阶近似是有效的，在这个信赖区域内，我们认为近似是有效的，超出这个区域，近似就可能出现问题。

通过近似模型和实际函数之间的差异来确定信赖区域的范围，定义一个指标 $\rho$ 来刻画近似的好坏：
$$\rho =实际函数下降的值/近似模型下降的值$$
当 $\rho$ 接近于1，则认为近似是良好的；当 $\rho$ 太小，说明实际函数下降的值远小于近似模型下降的值，意味着近似比较差，需要缩小近似范围；当 $\rho$ 比较大，说明实际函数下降的值比近似模型下降的值要大，意味着可以放大近似范围。

这部分理论先不写了，SLAM十四讲中跳得太快，理解能力有限。

直接上结论：
列文伯格——马夸尔特法在高斯牛顿法的基础上，加上了一个不等式约束：
$$||D\Delta w_k|{|}^2\le \mu$$
$\mu$ 是信赖区域的半径，$D$ 为系数矩阵，列文伯格把 $D$ 取为单位阵 $I$ ，相当于把 $\Delta w_k$ 约束在一个高维球内，而马夸尔特特出可以将 $D$ 取为非负数的对角阵——实际上通常用 $J^{T}J$ 的对角元素的平方根，使得$\Delta w_k$ 被约束在一个高维椭圆内，椭圆短轴对应梯度小的维度。

将该不等式约束用拉格朗日乘数法的方式放入目标函数内，构成拉格朗日函数，令该拉格朗日函数关于 $\Delta w$ 的导数为零，最终的形式核心仍然是求解一个关于增量的线性方程组：
$$(H+\lambda D^{T}D)\Delta w_k=g$$

可以看到，相比于高斯牛顿法，增量方程中仅多了一项 $\lambda D^{T}D$。$\lambda$ 为拉格朗日乘子。

一方面，当参数 $\lambda$ 较小时，$H$ 占主要地位，说明二次近似模型在该范围内的近似效果比较好，此时列文伯格——马夸尔特法更接近于高斯牛顿法。

另一方面，当参数 $\lambda$ 较大时，$\lambda D^{T}D$ 占主要地位，说明二次近似模型在该范围内的近似效果并不好，若取简化形式 $D=I$ ，$(H+\lambda I)\Delta w_k=g$，此时列文伯格——马夸尔特法更接近一阶梯度下降法（即最速下降法，别忘了，右边 $g=-J(w{)}^{T}r(w)$ ，里面是有一个一阶导数的）。

列文伯格——马夸尔特法在一定程度上避免了 $H$ 不可逆的问题，提供更稳定、更准确的增量$\Delta w_k$ 的计算方式。由于细节部分尚未完全吃透，这里先直接贴一下SLAM十四讲里的算法流程：

这里的式（6.34）是上文中计算 $\rho$ 的式子。

参考：
1.什么是「最小二乘法」？
https://www.zhihu.com/question/304164814/answer/549972357?utm_source=wechat_session
2.[数值计算] 数据拟合——非线性最小二乘法
https://zhuanlan.zhihu.com/p/83320557
3.梯度的方向为什么是函数值增加最快的方向？
https://zhuanlan.zhihu.com/p/38525412
4.泰勒展开式&牛顿法
https://blog.csdn.net/weixin_42468475/article/details/114642193
5.到底什么是线性函数，什么是非线性函数
https://blog.csdn.net/qq_38973721/article/details/108709827
6.线性最小二乘法与非线性最小二乘法
https://blog.csdn.net/qq_40061206/article/details/118703747
7.多元函数的泰勒(Taylor)展开式
https://blog.csdn.net/red_stone1/article/details/70260070
8.雅可比矩阵与海森矩阵
https://blog.csdn.net/hero_cjn/article/details/87888801
9.《视觉SLAM十四讲》——高翔、张涛
10. 0 范数、1 范数、2 范数有什么区别？
https://www.zhihu.com/question/20473040
11.熟悉陌生的2-范数（向量的模）
https://blog.csdn.net/geter_CS/article/details/85100045
12.《The matrix cookbook》——Kaare Brandt Petersen、Michael Syskind Pedersen
13.【简要复习】牛顿法、高斯—牛顿法（Gauss–Newton）和其他拟牛顿法
https://zhuanlan.zhihu.com/p/162901829
14.《最优化理论与方法》——王景恒