视觉SLAM理论到实践系列文章
下面是《视觉SLAM十四讲》学习笔记的系列记录的总链接，本人发表这个系列的文章链接均收录于此
视觉SLAM理论到实践系列文章链接
下面是专栏地址：
视觉SLAM理论到实践专栏
前言
高翔博士的《视觉SLAM14讲》学习笔记的系列记录
视觉SLAM理论到实践系列（五）——非线性优化（1）
前期准备
优化问题
组成
- 优化对象
- 目标函数/损失函数/评价函数
- 约束条件
分类
凸优化
条件
- 目标函数是凸的
- 不等式约束是凸的
- 所在的空间是凸的
- 等式约束是仿射形式(affine)
非凸优化
常见形式
（1）Linear Programming
（2）Quadratic Programming
（3）Quadratically Constrained QP
（4）Second-order Conic Programming
（5）Semi-Definite Programming
只要不是凸的，那就是非凸
判断问题是凸的还是非凸的，同样也是门学问
常见的处理思路，是借助松弛的思想，将非凸问题近似为凸问题：
此时需要清楚原非凸问题中的哪个或哪些函数导致了非凸性，之后考虑的是如何用凸优化模型来逼近原问题
说白了，就是在保留原问题部分性质的条件下，使用简单的项替代目标函数中难以处理的项，进而使得问题更易求解
数据来源
（1）批量式batch
一次性给定一批观测数据
如何根据这些数据估计出最准确的结果
（2）增量式incremental
数据是随着时间逐渐采集的
如何根据后面的数据调整之前的估计结果
在本讲，我们只考虑批量数据，在有静态的现成数据和明晰的优化目标情况下，如何求解优化问题
优化对象
合理选择
SLAM
- 位置
- 姿态
- 传感器参数
机器人抽象量（上位机）
- 时间多项式轨迹
- 推力大小方向
- 移动速度
- 转向角
执行器底层量（下位机）
- 电机转速
- 舵机转角
注意：对于相同目标，选择不同的状态量，问题的收敛的速度、构造的难度是有区别的
维度削减
微分平坦 differential flatness
描述飞机状态12个参量
经过微分平坦之后，归纳为4个量(flat output) x , y , z , y a w x ,y ,z ,yaw x,y,z,yaw ，可以理解为这四个量和它们自身的微分张成了无人机的状态空间
目标函数
目标函数的选择与要处理的问题密切相关
机器人定位与建图
状态估计角度
- 最大后验估计
- 最大似然估计
位姿优化角度
最小化重投影误差
从特征角度
点云配准 alignment
前端路径点搜索
djistra / A* / PRM / RRT* /…
生成一些路径点
后端优化生成轨迹
-
pos
如 x = a 1 ∗ t 5 + a 2 ∗ t 4 x=a_1*t^5+a_2*t^4 x=a1∗t5+a2∗t4
-
vel（速度），可由 pos 求导得到
-
acc （加速度），可由 vel 求导得到
微分平坦下与姿态关联
-
jerk，可由 acc 求导得到
-
姿态微分，表示无人机转动快慢
min_jerk => 让无人机尽可能少的抖动
-
微分平坦下与油门推力关联
-
-
snap，可由 jerk 求导得到
-
油门量微分，表示无人机油门调整速度
min_snap => 让无人机尽量节省能量
-
轨迹跟踪
最优控制 Optimal Control
Linear Quadratic Regulator (LQR)
模型预测控制 MPC
约束条件
无约束优化往往比有约束优化更容易求解，我们很想把有约束优化化成无约束优化问题
将约束加在目标函数中
- 拉格朗日乘数
- 罚函数
转换对象空间
-
三角换元
使用三角函数+角度量，来化简一些问题，并且很容易地用三角函数化简公式(倍角公式、正余弦和角公式，诱导公式，辅助角公式…)
-
流形空间(manifold)
优化位姿的时候，我们用李代数来代替李群，就是这样一种思路
问题求解
解析解
求导，找极值点
数值解
求梯度，找增量
常见优化库
（1）（常用）Ceres-Solver：http://ceres-solver.org/
代码开源，可用于求解具有边界约束的非线性最小二乘问题和一般无约束优化问题，google开发
（2）（常用）g2o：https://openslam-org.github.io/g2o.html
代码开源，用于解决非线性最小二乘问题，特定是将优化中涉及的对象用图(graph)的形式来表示
（3）（常用）gtsam：https://gtsam.org/
使用因子图和贝叶斯网络作为底层计算范式而不是稀疏矩阵来优化最可能的配置，Lio-sam
（4）（常用）OSQP：https://osqp.org/
多位大牛联合开发，用来求解QP问题的求解器，用的人很多，也很容易上手
作业会布置一个题目
（5）（常用）Casadi：https://web.casadi.org/
用于数值优化和最优控制(即涉及微分方程的优化)的开源软件工具；使用算法微分有效地生成导数信息，建立、求解和执行常微分方程(ODE)或微分代数方程(DAE)系统的正向和伴随灵敏度分析，以及制定和解决非线性程序(NLP)问题和最优控制问题(OCP)；可用于c++， Python和MATLAB，在性能上几乎没有差异；C++文档有限，python文档很全，最好从python入手；
作业会布置一个题目
（6）MOSEK：https://www.mosek.com/
适合的问题多，学生能免费用，无源代码，库仅支持x86
（7）OOQP：https://pages.cs.wisc.edu/~swright/ooqp/
快速鲁棒的QP求解器，代码开源
（8）GLPK：https://www.gnu.org/software/glpk/
快速鲁棒LP求解器，代码开源
（9）CVX：http://cvxr.com/cvx/
MATLAB封装包，内嵌多种优化库
（10）Acado：https://acado.github.io/index.html
用于自动控制和动态优化的软件环境和算法集；它为使用各种算法进行直接最优控制提供了一个通用框架，包括模型预测控制、状态和参数估计以及鲁棒优化。ACADO工具包实现为自包含的c++代码，并附带用户友好的MATLAB界面
学习建议
优化只是一种解决问题的手段
核心在于如何将问题构造成优化问题
不要拘泥于某种优化工具，去尝试各种优化工具，不同的优化库特点不同，擅长的领域也各有不同
学习优化的最好方式，是看书和思考，并时常和他人交流，视频只能是辅助，只看视频会越来浮躁，直到完全学不下去
学习资料
-
最优化：建模、算法与理论
中文教材，介绍近年来经典算法流程和理论，可用来入门
-
Numerical Optimization
理论方面很全，还关注了优化模型在计算机有限精度下优化结果的稳定性和鲁棒性，具有很好的工程实践指导
-
Lectures on Convex Optimization
思路清晰，讲凸优化的理论问题，光滑、非光滑、结构优化等
-
Lectures on modern convex optimization
关注锥优化、多项式求解复杂度问题
补充：向量、矩阵求导知识
1）通项求导：https://www.bilibili.com/video/BV1xk4y1B7RQ/?spm_id_from=333.788.video.desc.click&vd_source=01f587c37f2deec29c814f0c36677696
2）凑微分求导：https://www.bilibili.com/video/BV1vV4y1p7Nn/?spm_id_from=333.337.search-card.all.click&vd_source=0da0b7e545e1a65e82836ac4eff73077
3）规律+查表求导：https://en.wikipedia.org/wiki/Matrix_calculus
4）Dog-leg方法：https://blog.csdn.net/qq_35590091/article/details/94628887
上图中向量对矩阵求导，以及矩阵对向量或矩阵求导，结果为张量，不好显示
矩阵或者向量求导，只是批量求导数的一种方式，本质还是矩阵或者向量中的标量之间在求导数，只不过是借助矩阵或者向量的形式，同时对多个因变量关于自变量求导
因此，所谓矩阵或者向量求导，只不过是对求导结果按照约定好的规则组织起来罢了
两种布局
分母Denominator（常用）
分子Numerator
注意： 分子布局导数 = 分母布局导 数 T 分子布局导数 = 分母布局导数 ^T 分子布局导数=分母布局导数T
求导思路
（1）将原函数 写成通项或者分块
详见：https://www.bilibili.com/video/BV1xk4y1B7RQ/?p=5&vd_source=0da0b7e545e1a65e82836ac4eff73077
（2）借助微分定义凑导数
详见：https://www.bilibili.com/video/BV1vV4y1p7Nn/?spm_id_from=333.337.search-card.all.click&vd_source=0da0b7e545e1a65e82836ac4eff73077
微分求导法
（3）使用规律+查表 => 导数
详见：https://en.wikipedia.org/wiki/Matrix_calculus
求导规律
1）常量可提到导数外
2）“分子” 线性组合可拆开
3）转置符号可以提出来
导数表
1）scalar by scalar
2）（常用）vector by vector
默认列向量关于行向量求导
也可列列，行行，行列
可以用通用的矩阵方法来记
3）（常用）scalar by vector
最小二乘的二阶项常用
Identities: scalar-by-vector ∂ y ∂ x = ∇ x y \frac{\partial y}{\partial \mathbf{x}}=\nabla_{\mathbf{x}}y ∂x∂y=∇xy
4）vector by scalar
NOTE: The formulas involving the vector-by-vector derivatives ∂ g ( u ) ∂ u \frac{\partial \mathbf{g(u)}}{\partial \mathbf{u}} ∂u∂g(u) and ∂ f ( g ) ∂ g \frac{\partial \mathbf{f(g)}}{\partial \mathbf{g}} ∂g∂f(g) (whose outputs are matrices) assume the matrices are laid out consistent with the vector layout, i.e. numerator-layout matrix when numerator-layout vector and vice versa; otherwise, transpose the vector-by-vector derivatives.
注：涉及向量对向量导数 ∂ g ( u ) ∂ u \frac{\partial \mathbf{g(u)}}{\partial \mathbf{u}} ∂u∂g(u) 和 ∂ f ( g ) ∂ g \frac{\partial \mathbf{f(g)}}{\partial \mathbf{g}} ∂g∂f(g)（其输出为矩阵）的公式假设矩阵的布局与向量布局一致，即当分子布局向量时为分子布局矩阵，反之亦然；否则，用向量导数转置向量。
5）scalar by matrix
6）matrix by scalar
重点导数
vector by vector
即
∂
a
⃗
⋅
x
⃗
∂
x
⃗
=
a
⃗
\frac{\partial \vec{\mathbf{a}} \cdot\vec{\mathbf{x}}}{\partial \vec{\mathbf{x}}}=\vec{\mathbf{a}}
∂x∂a⋅x=a
- 通项+布局
- 凑微分定义
- 规律+查表
scalar by vector
即
∂
x
⃗
T
A
x
⃗
∂
x
⃗
=
(
A
+
A
T
)
x
⃗
\frac{\partial {\vec{\mathbf{x}}}^\text{T} A\vec{\mathbf{x}}}{\partial \vec{\mathbf{x}}}=(A+A^T)\vec{\mathbf{x}}
∂x∂xTAx=(A+AT)x
- 通项+布局
- 凑微分定义
- 规律+查表
注意
不管元素的排列方式如何，求导还是按照元素求导，位元素结果是确定的，只不过元素排列有差别
分子分母形式都有人用，还有人混用，不必过于纠结，关注式子传达的含义更重要（十四讲附录B）
自己写的时候按照一种规范就行，建议分母构型
矩阵和向量是不满足交换律的，所以不能随意变动矩阵或向量之间的相对位置，常量也不行（见 Matrix Calculus）
求导法则+查表 能基本解决日常求导问题
非线性最小二乘问题
上面介绍了优化与优化库，介绍了矩阵向量求导方法，下面来说一种具体优化问题的解法：无约束最小二乘是优化问题；
关于非线性，我们只需要将其在局部线性化就可以了，具体的求解方式，我们选择迭代求解（这是一种较为通用的方法，因为复杂问题通常是难以求出解析解的）
这里作者通过讲机器人状态估计，来引出非线性最小二乘问题；
非线性最小二乘的求解是比较独立的模块，是解决非线性最小二乘问题的常用方法，并不是状态估计独有的；
所以，先介绍如何求解一个非线性最小二乘问题，再和大家一起看看，我们是如何建立优化问题，来估计机器人的位姿的。
问题描述
min x F ( x ) = 1 2 ∥ f ( x ) ∥ 2 2 \min\limits_{x}F(\boldsymbol{x})=\frac{1}{2}\|f(\boldsymbol{x})\|_2^2 xminF(x)=21∥f(x)∥22
其中，自变量 x ∈ R n x\in \mathbb{R}^n x∈Rn， f f f 是任意标量非线性函数 f ( x ) : R n ↦ R f(x):\mathbb{R}^n\mapsto \mathbb{R} f(x):Rn↦R
注意这里的系数 1 2 \frac{1}{2} 21 是无关紧要的，有些文献上带有这个系数，有些文献则不带，它不会影响之后的结论。
显然，如果
f
f
f 是个数学形式上很简单的函数，那么该问题可以用解析形式来求。令目标函数的导数为零，然后求解
x
x
x 的最优值，就和求二元函数的极值一样：
d
F
d
x
=
0
.
\frac{\mathrm{d}F}{\mathrm{d}\boldsymbol{x}}=\boldsymbol{0}.
dxdF=0.
解此方程，就得到了导数为零处的极值。
这个方程是否容易求解呢？这取决于 f f f 导函数的形式。如果 f f f 为简单的线性函数，那么这个问题就是简单的线性最小二乘问题，但是有些导函数可能形式复杂，使得该方程可能不容易求解。求解这个方程需要我们知道关于目标函数的全局性质，而通常这是不大可能的。
解法
解析解
f ( x ) f(x) f(x) 的形式复杂，不一定能求出解析解
数值解
对于不方便直接求解的最小二乘问题，我们可以用迭代的方式，从一个初始值出发不断地更新当前的优化变量，使目标函数下降。具体步骤可列写如下：
- 给定某个初始值 x 0 x_0 x0
- 对于第 k k k 次迭代，寻找一个增量 Δ x k \Delta x_k Δxk ，使得 ∥ f ( x k + Δ x k ) ∥ 2 2 \|f({x}_k+\Delta{x}_k)\|_2^2 ∥f(xk+Δxk)∥22 达到极小值。
- 若 Δ x k \Delta x_k Δxk 足够小，则停止。
- 否则，令 x k + 1 = x k + Δ x k x_{k+1}=x_k+\Delta x_k xk+1=xk+Δxk，返回第2步。
这让求解导函数为零的问题变成了一个不断寻找下降增量 Δ x k \Delta x_k Δxk 的问题，我们将看到，由于可以对 f f f 进行线性化，增量的计算将简单很多。当函数下降直到增量非常小的时候，就认为算法收敛，目标函数达到了一个极小值。
在这个过程中，问题在于如何找到每次迭代点的增量，而这是一个局部的问题，我们只需要关心 f f f 在迭代值处的局部性质而非全局性质。
问题关键：求 Δ x k \Delta x_k Δxk
- 各个维度加还是减
- 加多少 or 减多少
直接法
泰勒展开：
一元：
f
(
x
)
=
f
(
x
k
)
+
f
′
(
x
k
)
(
x
−
x
k
)
+
1
2
f
′
′
(
x
)
(
x
−
x
k
)
2
+
⋯
f(x)=f(x_{k})+f^{\prime}(x_{k})(x-x_{k})+\frac{1}{2}f^{\prime\prime}(x)(x-x_{k})^{2}+\cdots
f(x)=f(xk)+f′(xk)(x−xk)+21f′′(x)(x−xk)2+⋯
多元：
f
(
X
)
=
f
(
X
k
)
+
∑
∂
f
(
X
k
)
∂
x
i
(
x
i
−
x
i
k
)
+
1
2
∑
∂
2
f
(
X
k
)
∂
x
j
∂
x
j
(
x
i
−
x
i
k
)
(
x
j
−
x
j
k
)
+
⋯
f(X)=f(X_{k})+\sum\frac{\partial f(X_{k})}{\partial x_{i}}(x_{i}-x_{i}^{k})+\frac{1}{2}\sum\frac{\partial^{2}f(X_{k})}{\partial x_{j}\partial x_{j}}(x_{i}-x_{i}^{k})(x_{j}-x_{j}^{k})+\cdots
f(X)=f(Xk)+∑∂xi∂f(Xk)(xi−xik)+21∑∂xj∂xj∂2f(Xk)(xi−xik)(xj−xjk)+⋯
其中一阶导的项为：
[
∂
f
∂
x
1
∂
f
∂
x
2
⋯
∂
f
∂
x
n
]
⌊
x
1
−
x
1
k
x
2
−
x
2
k
⋮
x
n
−
x
n
k
⌋
=
[
∇
f
(
X
k
)
]
T
⋅
(
X
−
X
k
)
=
∇
T
f
(
X
k
)
⋅
(
X
−
X
k
)
\begin{aligned} &\begin{bmatrix}\frac{\partial f}{\partial x_{1}}&\frac{\partial f}{\partial x_{2}}&\cdots&\frac{\partial f}{\partial x_{n}}\end{bmatrix}\left\lfloor\begin{array}{c}x_{1}-x_{1}^{k}\\x_{2}-x_{2}^{k}\\\vdots\\x_{n}-x_{n}^{k}\end{array}\right\rfloor \\ &=[\nabla f(X_{k})]^{T}\cdot(X-X_{k}) \\ &=\nabla^{T}f(X_{k})\cdot(X-X_{k}) \end{aligned}
[∂x1∂f∂x2∂f⋯∂xn∂f]
x1−x1kx2−x2k⋮xn−xnk
=[∇f(Xk)]T⋅(X−Xk)=∇Tf(Xk)⋅(X−Xk)
其中
∇
f
(
X
k
)
\nabla f(X_{k})
∇f(Xk) 为梯度，又叫雅可比矩阵，又写作
J
J
J，是一个列向量
二阶导的项为：
1
2
[
x
1
−
x
1
k
x
2
−
x
2
k
⋯
x
n
−
x
n
k
]
[
∂
2
f
∂
x
1
2
∂
2
f
∂
x
1
∂
x
2
.
.
.
∂
2
f
∂
x
1
∂
x
n
∂
2
f
∂
x
2
∂
x
1
∂
2
f
∂
x
2
2
.
.
.
∂
2
f
∂
x
2
∂
x
n
⋮
⋮
⋱
⋮
∂
2
f
∂
x
n
∂
x
1
∂
2
f
∂
x
n
∂
x
2
.
.
.
∂
2
f
∂
x
n
2
]
[
x
1
−
x
1
k
x
2
−
x
2
k
⋮
x
n
−
x
n
k
]
=
1
2
(
X
−
X
k
)
T
∇
2
f
(
X
k
)
(
X
−
X
k
)
\frac{1}{2}\Big[\begin{matrix}x_1-x_1^k&x_2-x_2^k&\cdots&x_n-x_n^k\\\end{matrix}\Big] \begin{bmatrix}\frac{\partial^2f}{\partial x_1^2}&\frac{\partial^2f}{\partial x_1\partial x_2}&...&\frac{\partial^2f}{\partial x_1\partial x_n}\\\frac{\partial^2f}{\partial x_2\partial x_1}&\frac{\partial^2f}{\partial x_2^2}&...&\frac{\partial^2f}{\partial x_2\partial x_n}\\\vdots&\vdots&\ddots&\vdots\\\frac{\partial^2f}{\partial x_n\partial x_1}&\frac{\partial^2f}{\partial x_n\partial x_2}&...&\frac{\partial^2f}{\partial x_n^2}\end{bmatrix} \begin{bmatrix}x_1-x_1^k\\x_2-x_2^k\\\vdots\\x_n-x_n^k\end{bmatrix} \\\\ =\frac{1}{2}(X-X_{k})^{T}\nabla^{2}f(X_{k})(X-X_{k})
21[x1−x1kx2−x2k⋯xn−xnk]
∂x12∂2f∂x2∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f∂x22∂2f⋮∂xn∂x2∂2f......⋱...∂x1∂xn∂2f∂x2∂xn∂2f⋮∂xn2∂2f
x1−x1kx2−x2k⋮xn−xnk
=21(X−Xk)T∇2f(Xk)(X−Xk)
其中
∇
2
f
(
X
k
)
\nabla^{2}f(X_{k})
∇2f(Xk) 为海森矩阵，又写作
H
H
H
一阶（梯度下降）
现在考虑第 k k k 次迭代，假设我们在 x k x_k xk 处，想要寻到增量 Δ x k \Delta x_k Δxk ，那么最直观的方式是将目标函数在 x k x_k xk 附近进行泰勒展开
在
x
k
x_k
xk 处展开
F
(
x
k
+
Δ
x
k
)
≈
F
(
x
k
)
+
J
(
x
k
)
T
Δ
x
k
+
1
2
Δ
x
k
T
H
(
x
k
)
Δ
x
k
.
F(x_{k}+\Delta x_{k})\approx F(x_{k})+J\left(x_{k}\right)^{\text{T}}\Delta x_{k}+\frac{1}{2}\Delta x_{k}^{\text{T}}H(\boldsymbol{x}_{k})\Delta x_{k}.
F(xk+Δxk)≈F(xk)+J(xk)TΔxk+21ΔxkTH(xk)Δxk.
其中
J
(
x
k
)
J(x_k)
J(xk) 是
F
(
x
)
F(x)
F(x) 关于
x
x
x 的一阶导数（也叫梯度，雅克比 Jacobian矩阵），
H
H
H 是二阶导数（也叫海森矩阵 Hessian ），
x
T
A
x
x^TAx
xTAx 是二次型
注：我们把 J ( x ) J(x) J(x) 写成列向量，它可以和 Δ x \Delta x Δx 进行内积，得到一个标量。
-
其中 J ( x k ) = [ ∂ F ( x ) ∂ x 1 , . . . , ∂ F ( x ) ∂ x n ] T \mathrm{J\left(x_k\right)=[\frac{\partial F\left(x\right)}{\partial x_1},...,\frac{\partial F\left(x\right)}{\partial x_n}]^T} J(xk)=[∂x1∂F(x),...,∂xn∂F(x)]T ，是 F ( x ) F(x) F(x) 关于 x x x 的一阶导数在 x k x_k xk 处取值， J J J 称为梯度或雅克比矩阵，需要注意的是 $ J({x_k})$ 是一个 n n n 维列向量， Δ x k \Delta x_k Δxk 也是一个 n n n 维列向量，因此泰勒展开中是 J ( x k ) J({x_k}) J(xk) 的转置乘以 $ \Delta x_k$ ，从而得到一个标量。
-
H ( x k ) H(x_k) H(xk) 是 $ F(x)$ 关于 x x x 的二阶导数在 x k x_k xk 处取值， H H H 称为海塞矩阵，是一个 n × n n\times n n×n 的方阵，且是一个对称阵（当二次连续时）。
其形式为：
H = [ ∂ 2 F ∂ x 1 2 ∂ 2 F ∂ x 1 ∂ x 2 . . . ∂ 2 F ∂ x 1 ∂ x n ∂ 2 F ∂ x 2 ∂ x 1 ∂ 2 F ∂ x 2 2 . . . ∂ 2 F ∂ x 2 ∂ x n ⋮ ⋮ ⋱ ⋮ ∂ 2 F ∂ x n ∂ x 1 ∂ 2 F ∂ x n ∂ x 2 . . . ∂ 2 F ∂ x n 2 ] H= \begin{bmatrix}\frac{\partial^2F}{\partial x_1^2}&\frac{\partial^2F}{\partial x_1\partial x_2}&...&\frac{\partial^2F}{\partial x_1\partial x_n}\\\frac{\partial^2F}{\partial x_2\partial x_1}&\frac{\partial^2F}{\partial x_2^2}&...&\frac{\partial^2F}{\partial x_2\partial x_n}\\\vdots&\vdots&\ddots&\vdots\\\frac{\partial^2F}{\partial x_n\partial x_1}&\frac{\partial^2F}{\partial x_n\partial x_2}&...&\frac{\partial^2F}{\partial x_n^2}\end{bmatrix} H= ∂x12∂2F∂x2∂x1∂2F⋮∂xn∂x1∂2F∂x1∂x2∂2F∂x22∂2F⋮∂xn∂x2∂2F......⋱...∂x1∂xn∂2F∂x2∂xn∂2F⋮∂xn2∂2F
我们可以选择保留泰勒展开的一阶或二阶项，那么对应的求解方法则称为一阶梯度或二阶梯度法
如果保留一阶梯度，那么取增量为反向的梯度，即可保证函数下降：
Δ
x
∗
=
−
J
(
x
k
)
=
−
∇
F
(
x
)
\Delta x^* =-J\left(x_{k}\right)=-\nabla F(x)
Δx∗=−J(xk)=−∇F(x)
当然这只是个方向，通常我们还要再指定一个步长
λ
\lambda
λ ，步长可以根据一定的条件来计算。
这种方法被称为最速下降法，它的直观意义非常简单，只要我们沿着反向梯度方向前进，在一阶（线性）的近似下，目标函数必定会下降。
注意：
（1）这里的梯度是目标函数的梯度
F
(
x
)
=
1
2
∥
f
(
x
)
∥
2
2
∂
F
(
x
)
∂
x
=
1
2
∂
f
(
x
)
∂
x
∂
(
f
(
x
)
)
2
∂
f
(
x
)
=
J
(
x
)
f
(
x
)
F(\boldsymbol{x})=\frac{1}{2}\|f(\boldsymbol{x})\|_2^2 \\\\ \frac{\partial F(x)}{\partial x}=\frac{1}{2}\frac{\partial f(x)}{\partial x}\frac{\partial (f(x))^2}{\partial f(x)} =J(x)f(x)
F(x)=21∥f(x)∥22∂x∂F(x)=21∂x∂f(x)∂f(x)∂(f(x))2=J(x)f(x)
这里
J
(
x
)
J(x)
J(x) 为
∂
f
(
x
)
∂
x
\frac{\partial f(x)}{\partial x}
∂x∂f(x)
（2）导数(梯度)是增量的方向
（3）梯度取反就是梯度下降的方向
（4）通常不会直接让
J
J
J 代表步长，会加因子or饱和函数（深度学习中把这个因子叫做学习率），即
Δ
x
∗
=
−
α
J
(
x
k
)
=
−
α
∇
F
(
x
)
\Delta x^* =-\alpha J\left(x_{k}\right)=-\alpha \nabla F(x)
Δx∗=−αJ(xk)=−α∇F(x)
因子可以根据一些算法来计算，也可以设置为一个比较小的值
二阶（牛顿法）
仅使用一阶梯度作为下降方向有时并不能很好地拟合函数的局部趋势
对二阶展开，将步长
Δ
x
k
\Delta x_k
Δxk 作为
x
=
x
k
x=x_k
x=xk 优化变量
Δ
x
∗
=
arg
min
(
F
(
x
)
+
J
(
x
)
T
Δ
x
+
1
2
Δ
x
T
H
Δ
x
)
\Delta\boldsymbol{x}^*=\arg\min\bigg(F\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^{T}\Delta\boldsymbol{x}+\frac{1}{2}\Delta\boldsymbol{x}^{T}H\Delta\boldsymbol{x}\bigg)
Δx∗=argmin(F(x)+J(x)TΔx+21ΔxTHΔx)
右侧只含
Δ
x
\Delta x
Δx 的零次、一次和二次项。
求右侧等式关于
Δ
x
\Delta x
Δx 的导数
∂
F
(
x
)
+
J
(
x
)
T
Δ
x
+
1
2
Δ
x
T
H
Δ
x
∂
Δ
x
=
0
+
∂
J
(
x
)
⋅
Δ
x
∂
Δ
x
+
1
2
Δ
x
T
H
Δ
x
∂
Δ
x
=
J
(
x
)
+
1
2
(
H
+
H
T
)
Δ
x
=
J
(
x
)
+
H
Δ
x
\begin{aligned} & \frac{\partial F\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^{T}\Delta\boldsymbol{x}+\frac{1}{2}\Delta\boldsymbol{x}^{T}H\Delta\boldsymbol{x}}{\partial \Delta x} \\ &=0+\frac{\partial J(x)\cdot \Delta x}{\partial \Delta x}+\frac{1}{2}\frac{\Delta\boldsymbol{x}^{T}H\Delta\boldsymbol{x}}{\partial \Delta x} \\ &=J(x)+\frac{1}{2}(H+H^T)\Delta x \\ &=J(x)+H\Delta x \end{aligned}
∂Δx∂F(x)+J(x)TΔx+21ΔxTHΔx=0+∂Δx∂J(x)⋅Δx+21∂ΔxΔxTHΔx=J(x)+21(H+HT)Δx=J(x)+HΔx
令它为零，得到
J
+
H
Δ
x
=
0
⇒
H
Δ
x
=
−
J
J+H\Delta x=0\Rightarrow H\Delta x=-J
J+HΔx=0⇒HΔx=−J
求解这个线性方程，就得到了增量，该方法又称为牛顿法
这个式子很神奇，这是在 F ( x ) F(x) F(x) 的基础上计算最合适的 Δ x \Delta x Δx，使得 F ( x + Δ x ) ~F(x+\Delta x) F(x+Δx) 达到最小
x x x 在这里是一个参数，真正的变量是 Δ x \Delta x Δx，当 x x x 确定下来后， Δ x \Delta x Δx 也就直接确定下来了
如果能求出 Δ x \Delta x Δx 的解析解，那么在任意位置 x x x 都能准确知道 Δ x \Delta x Δx，使得 F ( x + Δ x ) ~F(x+\Delta x) F(x+Δx) 达到最小的最佳 Δ x \Delta x Δx
退一步讲，限于 F ( x ) F(x) F(x) 的形式，我们无法求出解析解，将 argmin 作为目标函数，以 Δ x \Delta x Δx 为优化变量，也会让我们找到更加符合函数下降趋势的 Δ x \Delta x Δx
就这样，边寻找当前 x k x_k xk 下的最佳 Δ x \Delta x Δx，边更新下一时刻的 x k + 1 x_{k+1} xk+1，虽然二阶近似也不会非常准确，但是只要步长不过于大，迭代次数多，就能找到近乎准确的数值极小值
总结
一阶和二阶梯度法都十分直观，只要把函数在迭代点附近进行泰勒展开，并针对更新量做最小化即可。
事实上，我们用一个一次或二次的函数近似了原函数，然后用近似函数的最小值来猜测原函数的极小值。只要原目标函数局部看起来像一次或二次函数，这类算法就是成立的（这也是现实中的情形）。
缺点：
（1）最速下降法过于贪心，容易走出锯齿路线，反而增加了迭代次数。
（2）牛顿法则需要计算目标函数的 H H H 矩阵，这在问题规模较大时非常困难，我们通常倾向于避免 H H H 的计算。海森矩阵 H H H 计算不易，求逆也不易
对于一般的问题，一些拟牛顿法可以得到较好的结果，而对于最小二乘问题，还有几类更实用的方法：高斯牛顿法和列文伯格—马夸尔特方法。
改进二阶
高斯-牛顿法 G-N
参考：高斯牛顿法详解
高斯牛顿法是最优化算法中最简单的方法之一，它的思想是将 f ( x ) f(x) f(x) 进行一阶的泰勒展开
注意这里不是目标函数
F
(
x
)
F(x)
F(x) 而是
f
(
x
)
f(x)
f(x) ，否则就变成牛顿法了！
f
(
x
+
Δ
x
)
≈
f
(
x
)
+
J
(
x
)
T
Δ
x
.
f\left(\boldsymbol{x}+\Delta\boldsymbol{x}\right)\approx f\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^\text{T}\Delta\boldsymbol{x}.
f(x+Δx)≈f(x)+J(x)TΔx.
这里
J
(
x
)
T
J(x)^T
J(x)T 为
f
(
x
)
f(x)
f(x) 关于
x
x
x 的导数，为
n
×
1
n\times 1
n×1 的列向量。根据前面的框架，当前的目标是寻找增量
Δ
x
\Delta x
Δx ，使得
∥
f
(
x
+
Δ
x
)
∥
2
\|f(x+\Delta x)\|^2
∥f(x+Δx)∥2 达到最小。
为了求
Δ
x
\Delta x
Δx，我们需要解一个线性的最小二乘问题：
Δ
x
∗
=
arg
min
Δ
x
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
2
.
\Delta x^*=\arg\min\limits_{\Delta x}\frac{1}{2}\left\|f\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^{\mathsf{T}}\Delta\boldsymbol{x}\right\|^2.
Δx∗=argΔxmin21
f(x)+J(x)TΔx
2.
这个方程与之前的有什么不一样呢？对比之前的式子如下：
Δ
x
∗
=
arg
min
(
F
(
x
)
+
J
(
x
)
T
Δ
x
+
1
2
Δ
x
T
H
Δ
x
)
\Delta\boldsymbol{x}^*=\arg\min\bigg(F\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^{T}\Delta\boldsymbol{x}+\frac{1}{2}\Delta\boldsymbol{x}^{T}H\Delta\boldsymbol{x}\bigg)
Δx∗=argmin(F(x)+J(x)TΔx+21ΔxTHΔx)
可以看到高斯-牛顿法不带
H
H
H 矩阵，用一阶泰勒展开的平方来得到二阶项
根据极值条件，将上述目标函数对 Δ x \Delta x Δx 求导，并令导数为零。
为此，先展开目标函数的平方项
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
2
=
1
2
(
f
(
x
)
+
Δ
x
T
J
(
x
)
)
T
(
f
(
x
)
+
J
(
x
)
T
Δ
x
)
=
1
2
(
f
(
x
)
T
+
J
(
x
)
T
Δ
x
)
(
f
(
x
)
+
J
(
x
)
T
Δ
x
)
=
1
2
(
∥
f
(
x
)
∥
2
2
+
2
f
(
x
)
J
(
x
)
T
Δ
x
+
Δ
x
T
J
(
x
)
J
(
x
)
T
Δ
x
)
.
\begin{aligned} &\frac{1}{2}\|f\left(x\right)+J\left(\boldsymbol{x}\right)^{\mathsf{T}}\Delta x\|^{2} \\ &=\frac{1}{2}\Bigl(f\left(\boldsymbol{x}\right)+\Delta\boldsymbol{x}^{\mathsf{T}}J\left(\boldsymbol{x}\right)\Bigr)^{\mathsf{T}}\Bigl(f\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^{\mathsf{T}}\Delta\boldsymbol{x}\Bigr) \\ &=\frac{1}{2}\Bigl(f\left(\boldsymbol{x}\right)^{\mathsf{T}}+J\left(\boldsymbol{x}\right)^{\mathsf{T}}\Delta\boldsymbol{x}\Bigr)\Bigl(f\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^{\mathsf{T}}\Delta\boldsymbol{x}\Bigr) \\ &=\frac{1}{2}\left(\|f(\boldsymbol{x})\|_{2}^{2}+2f\left(\boldsymbol{x}\right)J(\boldsymbol{x})^{\mathsf{T}}\Delta\boldsymbol{x}+\Delta x^{\mathsf{T}}J(\boldsymbol{x})J(\boldsymbol{x})^{\mathsf{T}}\Delta\boldsymbol{x}\right). \end{aligned}
21∥f(x)+J(x)TΔx∥2=21(f(x)+ΔxTJ(x))T(f(x)+J(x)TΔx)=21(f(x)T+J(x)TΔx)(f(x)+J(x)TΔx)=21(∥f(x)∥22+2f(x)J(x)TΔx+ΔxTJ(x)J(x)TΔx).
注意：
（1）这里的计算方法是把 f ( x ) + J ( x ) T Δ x f\left(x\right)+J\left(\boldsymbol{x}\right)^{T}\Delta x f(x)+J(x)TΔx 当成向量来计算，实际上其为标量，直接做乘法即可
（2）二阶的牛顿法如下
Δ x ∗ = arg min ( F ( x ) + J ( x ) T Δ x + 1 2 Δ x T H Δ x ) \Delta\boldsymbol{x}^*=\arg\min\bigg(F\left(\boldsymbol{x}\right)+J\left(\boldsymbol{x}\right)^{T}\Delta\boldsymbol{x}+\frac{1}{2}\Delta\boldsymbol{x}^{T}H\Delta\boldsymbol{x}\bigg) Δx∗=argmin(F(x)+J(x)TΔx+21ΔxTHΔx)
对比可得，高斯-牛顿法是使用 J ( x ) J ( x ) T J(\boldsymbol{x})J(\boldsymbol{x})^{\mathsf{T}} J(x)J(x)T 来近似代替 H H H
求上式关于
Δ
x
\Delta x
Δx 的导数，并令其为零：
J
(
x
)
f
(
x
)
+
J
(
x
)
J
T
(
x
)
Δ
x
=
0
.
J(\boldsymbol{x})f\left(\boldsymbol{x}\right)+J(\boldsymbol{x})J^{T}\left(\boldsymbol{x}\right)\Delta\boldsymbol{x}=\boldsymbol{0}.
J(x)f(x)+J(x)JT(x)Δx=0.
可以得到如下方程组：
J
(
x
)
J
T
(
x
)
⏟
H
(
x
)
Δ
x
=
−
J
(
x
)
f
(
x
)
⏟
g
(
x
)
\begin{gathered} \underbrace{J(x)J^{\mathrm{T}}\left(\boldsymbol{x}\right)}_{H(x) }\Delta x=\underbrace{-J(\boldsymbol{x})f\left(\boldsymbol{x}\right)}_{g(x)} \\ \end{gathered}
H(x)
J(x)JT(x)Δx=g(x)
−J(x)f(x)
这个方程是关于变量
Δ
x
\Delta x
Δx 的线性方程组，我们称它为增量方程，也可以称为高斯牛顿方程(Gauss–Newton equation)或者正规方程(Normal equation)。
我们把左边的系数定义为
H
H
H ，右边定义为
g
g
g ，那么上式变为
H
Δ
x
=
g
H\Delta x=g
HΔx=g
这里把左侧记作
H
H
H 是有意义的。
此时我们可以得到
Δ
x
∗
\Delta\boldsymbol{x}^*
Δx∗
Δ
x
∗
=
−
(
J
(
x
)
J
T
(
x
)
)
−
1
J
(
x
)
f
(
x
)
\Delta\boldsymbol{x}^*=-\left(J(\boldsymbol{x})J^{T}\left(\boldsymbol{x}\right)\right)^{-1}J(\boldsymbol{x})f\left(\boldsymbol{x}\right)
Δx∗=−(J(x)JT(x))−1J(x)f(x)
注意：
（1）这里 J ( x ) J T ( x ) J(\boldsymbol{x})J^{T}\left(\boldsymbol{x}\right) J(x)JT(x) 是一个半正定矩阵，即 λ i ≥ 0 \lambda_i\geq0 λi≥0 ，即不一定有逆矩阵，当矩阵 J ( x ) J T ( x ) J(\boldsymbol{x})J^{T}\left(\boldsymbol{x}\right) J(x)JT(x) 出现一些病态时， Δ x ∗ \Delta\boldsymbol{x}^* Δx∗ 会出现不稳定，可能会导致运算的结果不收敛
（2）高斯牛顿是对 f ( x ) f(x) f(x) 展开，牛顿法是对 F ( x ) F(x) F(x) 展开，所以其实是使用 f ( x ) f(x) f(x) 的雅可比去近似 F ( x ) F(x) F(x) 的海森矩阵，也就是要注意两个展开式里的 J J J 不是同一个，一个是小 f f f 的 J J J ，一个是大 F F F 里的 J J J
对比牛顿法可见，高斯牛顿法用 J J T JJ^T JJT 作为牛顿法中二阶Hessian 矩阵的近似，从而省略了计算 H H H 的过程。
求解增量方程是整个优化问题的核心所在。
如果我们能够顺利解出该方程，那么高斯牛顿法的算法步骤可以写成
- 给定初始值 x 0 x_0 x0。
- 对于第 k k k 次迭代，求出当前的雅可比矩阵 J ( x k ) J(x_k) J(xk) 和误差 f ( x k ) f(x_k) f(xk) 。
- 求解增量方程： H Δ x k = g H\Delta x_k=g HΔxk=g 。
- 若 Δ x k \Delta x_k Δxk 足够小，则停止。否则，令 Δ x k + 1 = x k + Δ x k \Delta x_{k+1}=x_k+\Delta x_k Δxk+1=xk+Δxk ，返回第2步。
缺点
（1） J J T JJ^T JJT 是半正定的，不一定可逆
直观上看，原函数在局部不近似二次函数
这会使得迭代增量稳定性变差，算法不以收敛
（2）步长 Δ x \Delta x Δx 的尺度无法很好的把握，容易造成不收敛
为了求解增量方程，我们需要求解H-1,这需要H矩阵可逆，但实际数据中计算得到的 J J T JJ^T JJT 却只有半正定性。
也就是说，在使用高斯牛顿法时，可能出现
J
J
T
JJ^T
JJT 为奇异矩阵或者病态(ill-condition)的情况，此时增量的稳定性较差，导致算法不收敛。直观地说，原函数在这个点
的局部近似不像一个二次函数。
更严重的是，就算我们假设 H H H 非奇异也非病态，如果我们求出来的步长 Δ x \Delta x Δx 太大，也会导致我们采用的局部近似式不够准确，这样一来我们甚至无法保证它的迭代收敛，哪怕是让目标函数变得更大都是有可能的。
列文伯格-马夸尔特 L-M
列文伯格一马夸尔特方法在一定程度上修正了这些问题。一般认为它比高斯牛顿法更为健壮，但它的收敛速度可能比高斯牛顿法更慢，被称为阻尼牛顿法（Damped Newton Method)。
高斯牛顿法中采用的近似二阶泰勒展开只能在展开点附近有较好的近似效果，所以我们很自然地想到应该给 Δ x \Delta x Δx 添加一个范围，称为信赖区域(Trust Region)。
这个范围定义了在什么情况下二阶近似是有效的，这类方法也称为信赖区域方法（Trust Region Method)。
在信赖区域里，我们认为近似是有效的；出了这个区域，近似可能会出问题。
上图实线表示实际的函数，虚线表示近似的函数，可以看到信赖域范围内的近似是有效的
一个比较好的方法是根据我们的近似模型跟实际函
数之间的差异来确定：如果差异小，说明近似效果好，我们扩大近似的范围；反之，如果差异大，就缩小近似的范围。
我们定义一个指标
ρ
\rho
ρ 来刻画近似的好坏程度：
ρ
=
f
(
x
+
Δ
x
)
−
f
(
x
)
J
(
x
)
T
Δ
x
.
\rho=\frac{f\left(\boldsymbol{x}+\Delta\boldsymbol{x}\right)-f\left(\boldsymbol{x}\right)}{J\left(\boldsymbol{x}\right)^{\text{T}}\Delta\boldsymbol{x}}.
ρ=J(x)TΔxf(x+Δx)−f(x).
ρ
\rho
ρ 的分子是实际函数下降的值，分母是近似模型下降的值。
如果 ρ \rho ρ 接近于1，则近似是好的。
如果 ρ \rho ρ 太小，说明实际减小的值远少于近似减小的值，则认为近似比较差，需要缩小近似范围。
反之，如果 ρ \rho ρ 比较大，则说明实际下降的比预计的更大，我们可以放大近似范围。
简单来说：
分子 f ( x + Δ x ) f\left(\boldsymbol{x}+\Delta\boldsymbol{x}\right) f(x+Δx) 表示实际函数差，分母 J ( x ) T Δ x J\left(\boldsymbol{x}\right)^{\text{T}}\Delta\boldsymbol{x} J(x)TΔx 表示近似函数差
（1） ρ < 1 \rho<1 ρ<1 ，近似 比 实际大，激进，应缩小信任域
（2） ρ ≈ 1 \rho\approx1 ρ≈1 ，近似 贴近 实际，信任域维持不变
（3） ρ > 1 \rho>1 ρ>1 ，近似 比 实际大，保守，应扩大信任域
列-马 方法步骤：
-
给定初始值 x 0 x_0 x0 ，以及初始优化半径 μ \mu μ。
-
对于第 k k k 次迭代，在高斯牛顿法的基础上加上信赖区域，求解：
min Δ x k 1 2 ∥ f ( x k ) + J ( x k ) T Δ x k ∥ 2 , s.t. ∥ D Δ x k ∥ 2 ⩽ μ , \min\limits_{\Delta x_k}\frac{1}{2}\left\|f\left(x_k\right)+J\left(\boldsymbol{x}_k\right)^{\text{T}}\Delta x_k\right\|^2,\quad\text{s.t.}\quad\left\|D\Delta x_k\right\|^2\leqslant\mu, Δxkmin21 f(xk)+J(xk)TΔxk 2,s.t.∥DΔxk∥2⩽μ,
其中， μ \mu μ 是信赖区域的半径， D D D 为系数矩阵，将在后文说明 -
按式 ρ = f ( x + Δ x ) − f ( x ) J ( x ) T Δ x \rho=\frac{f\left(\boldsymbol{x}+\Delta\boldsymbol{x}\right)-f\left(\boldsymbol{x}\right)}{J\left(\boldsymbol{x}\right)^{\text{T}}\Delta\boldsymbol{x}} ρ=J(x)TΔxf(x+Δx)−f(x) 计算 ρ \rho ρ。
-
若 ρ > 3 4 \rho>\frac{3}{4} ρ>43 ，则设置 μ = 2 μ \mu=2\mu μ=2μ
-
.若 ρ > 1 4 \rho>\frac{1}{4} ρ>41 ，则设置 μ = 0.5 μ \mu=0.5\mu μ=0.5μ。
-
如果 ρ \rho ρ 大于某阈值，则认为近似可行。令 x k + 1 = x k + Δ x k x_{k+1}=x_k+\Delta x_k xk+1=xk+Δxk
-
判断算法是否收敛。如不收敛则返回第2步，否则结束。
这里近似范围扩大的倍数和阈值都是经验值，可以替换成别的数值。
在式 min Δ x k 1 2 ∥ f ( x k ) + J ( x k ) T Δ x k ∥ 2 , s.t. ∥ D Δ x k ∥ 2 ⩽ μ \min\limits_{\Delta x_k}\frac{1}{2}\left\|f\left(x_k\right)+J\left(\boldsymbol{x}_k\right)^{\text{T}}\Delta x_k\right\|^2,\quad\text{s.t.}\quad\left\|D\Delta x_k\right\|^2\leqslant\mu Δxkmin21 f(xk)+J(xk)TΔxk 2,s.t.∥DΔxk∥2⩽μ 中，我们把增量限定于一个半径为 μ \mu μ 的球中，认为只在这个球内才是有效的。
带上 D D D 之后，这个球可以看成一个椭球。
在列文伯格提出的优化方法中，把 D D D 取成单位阵 I I I ，相当于直接把 Δ x k \Delta x_k Δxk 约束在一个球中。随后，马夸尔特提出将 D D D 取成非负数对角阵——实际中通常用 J T J J^TJ JTJ 的对角元素平方根，使得在梯度小的维度上约束范围更大一些。
无论如何，在列文伯格一马夸尔特优化中，我们都需要解式 min Δ x k 1 2 ∥ f ( x k ) + J ( x k ) T Δ x k ∥ 2 , s.t. ∥ D Δ x k ∥ 2 ⩽ μ \min\limits_{\Delta x_k}\frac{1}{2}\left\|f\left(x_k\right)+J\left(\boldsymbol{x}_k\right)^{\text{T}}\Delta x_k\right\|^2,\quad\text{s.t.}\quad\left\|D\Delta x_k\right\|^2\leqslant\mu Δxkmin21 f(xk)+J(xk)TΔxk 2,s.t.∥DΔxk∥2⩽μ 那样一个子问题来获得梯度。
这个子问题是带不等式约束的优化问题，我们用拉格朗日乘子把约束项放到目标函数中，构成拉格朗日函数：
L
(
Δ
x
k
,
λ
)
=
1
2
∥
f
(
x
k
)
+
J
(
x
k
)
T
Δ
x
k
∥
2
+
λ
2
(
∥
D
Δ
x
k
∥
2
−
μ
)
.
\mathcal{L}(\Delta x_k,\lambda)=\frac{1}{2}\left\|f\left(\boldsymbol{x}_k\right)+\boldsymbol{J}\left(\boldsymbol{x}_k\right)^{\mathsf{T}}\Delta\boldsymbol{x}_k\right\|^2+\frac{\lambda}{2}\left(\left\|\boldsymbol{D}\Delta\boldsymbol{x}_k\right\|^2-\mu\right).
L(Δxk,λ)=21
f(xk)+J(xk)TΔxk
2+2λ(∥DΔxk∥2−μ).
这里
λ
\lambda
λ 为拉格朗日乘子。
注意：
（1）这里比较正规的做法是拉格朗日乘数法+KKT条件
（2）原书中注释道：严谨的读者可能不满意此处的叙述。信赖域原问题的约束条件除了拉格朗日函数求导为零，KKT条件还会有一些别的约束： λ > 0 \lambda>0 λ>0，且 λ ( ∥ D Δ x ∥ 2 − μ ) = 0 \lambda(\|D\Delta x\|^2-\mu)=0 λ(∥DΔx∥2−μ)=0。但是在L-M迭代中，我们不妨把它看成在原问题的目标函数上，以 λ \lambda λ 为权重的惩罚项(Argumented Lagrangian)。在每一步迭代后，若发现信赖域条件不满足，或者目标函数增加，就增加 λ \lambda λ 的权重，直到最终满足信赖域条件。所以，理论上对L-M算法存在不同的解释，但实际中我们只关心它是否顺利工作。
上式对
Δ
x
k
\Delta x_k
Δxk 求导可得
∂
L
∂
Δ
x
k
=
J
(
x
)
f
(
x
)
+
J
(
x
)
J
T
(
x
)
Δ
x
+
λ
2
∂
(
D
Δ
x
k
)
2
∂
Δ
x
k
=
J
(
x
)
f
(
x
)
+
J
(
x
)
J
T
(
x
)
Δ
x
+
λ
2
∂
Δ
x
k
T
D
T
D
Δ
x
k
∂
Δ
x
k
=
J
(
x
)
f
(
x
)
+
J
(
x
)
J
T
(
x
)
Δ
x
+
λ
2
(
D
T
D
+
D
D
T
)
Δ
x
k
=
J
(
x
)
f
(
x
)
⏟
−
g
+
J
(
x
)
J
T
(
x
)
⏟
H
Δ
x
+
λ
D
D
T
Δ
x
k
=
(
H
+
λ
D
T
D
)
Δ
x
k
−
g
\begin{aligned} \frac{\partial L}{\partial \Delta x_k} &=J({x})f\left({x}\right)+J({x})J^{T}\left({x}\right)\Delta{x}+\frac{\lambda}{2}\frac{\partial (D\Delta x_k)^2}{\partial \Delta x_k} \\ &=J({x})f\left({x}\right)+J({x})J^{T}\left({x}\right)\Delta{x}+\frac{\lambda}{2}\frac{\partial\Delta x_k^TD^T D\Delta x_k}{\partial \Delta x_k} \\ &=J({x})f\left({x}\right)+J({x})J^{T}\left({x}\right)\Delta{x}+\frac{\lambda}{2}(D^TD+DD^T)\Delta x_k \\ &=\underbrace{J({x})f\left({x}\right)}_{-g}+\underbrace{J({x})J^{T}\left({x}\right)}_{H}\Delta{x}+\lambda DD^T\Delta x_k \\ &=\left(H+\lambda D^{\mathrm{T}}D\right)\Delta x_{k}-g \end{aligned}
∂Δxk∂L=J(x)f(x)+J(x)JT(x)Δx+2λ∂Δxk∂(DΔxk)2=J(x)f(x)+J(x)JT(x)Δx+2λ∂Δxk∂ΔxkTDTDΔxk=J(x)f(x)+J(x)JT(x)Δx+2λ(DTD+DDT)Δxk=−g
J(x)f(x)+H
J(x)JT(x)Δx+λDDTΔxk=(H+λDTD)Δxk−g
其中
H
=
J
(
x
)
J
T
(
x
)
,
g
=
−
J
(
x
)
f
(
x
)
H=J({x})J^{T}\left({x}\right),g=-J({x})f\left({x}\right)
H=J(x)JT(x),g=−J(x)f(x)
类似于高斯牛顿法中的做法，令该拉格朗日函数关于
Δ
x
\Delta x
Δx 的导数为零，它的核心仍是计算增量的线性方程：
(
H
+
λ
D
T
D
)
Δ
x
k
−
g
=
0
(
H
+
λ
D
T
D
)
Δ
x
k
=
g
Δ
x
k
=
(
H
+
λ
D
T
D
)
−
1
g
Δ
x
k
=
−
(
J
(
x
)
J
T
(
x
)
+
λ
D
T
D
)
−
1
J
(
x
)
f
(
x
)
\begin{aligned} \left(H+\lambda D^{\mathrm{T}}D\right)\Delta x_{k}-g&=0 \\\\ \left(H+\lambda D^{\mathrm{T}}D\right)\Delta x_{k}&=g \\\\ \Delta x_{k}&=\left(H+\lambda D^{\mathrm{T}}D\right)^{-1}g \\\\ \Delta x_{k}&=-\left(J({x})J^{T}\left({x}\right)+\lambda D^{\mathrm{T}}D\right)^{-1}J({x})f\left({x}\right) \end{aligned}
(H+λDTD)Δxk−g(H+λDTD)ΔxkΔxkΔxk=0=g=(H+λDTD)−1g=−(J(x)JT(x)+λDTD)−1J(x)f(x)
可以看到，相比于高斯牛顿法，增量方程多了一项
λ
D
T
D
\lambda D^TD
λDTD。
如果考虑它的简化形式，即
D
=
I
D=I
D=I，那么相当于求解：
(
H
+
λ
I
)
Δ
x
k
=
g
\left(H+\lambda I\right)\Delta x_k={g}
(H+λI)Δxk=g
注意
对于
(
H
+
λ
I
)
Δ
x
k
=
g
\left(H+\lambda I\right)\Delta x_k={g}
(H+λI)Δxk=g
其中
H
=
J
(
x
)
J
T
(
x
)
,
g
=
−
J
(
x
)
f
(
x
)
H=J({x})J^{T}\left({x}\right),g=-J({x})f\left({x}\right)
H=J(x)JT(x),g=−J(x)f(x)
D D D 表示信任域的形状，单位阵表示球型域
（1）与高斯牛顿法相比，左边 H H H 上加了个 λ I \lambda I λI
（2）当 λ \lambda λ 比较小，退化为高斯牛顿方法
（3）当 λ \lambda λ 比较大，退化为梯度下降方法
L-M方法可以看成 梯度下降 与 高斯牛顿法 之间的融合，通过 λ \lambda λ 来切换两种方法
G-N等方法求方向，改变迭代步长
L-M方法是在给定信任域下，同时确定长度和方向
Dog-Leg
也是一种操作信任域方法
“人为地”定义了利用信赖域来选择增量大小
补充
初值问题
初值选择对于非凸优化的问题的求解结果十分重要【ICP，PnP，连续运动假设提供先验】
求逆运算
求逆是维度三次方复杂度操作，可以通过矩阵分解的方式简化【矩阵分析 矩阵论】
如使用 SVD分解，QR分解，Cholesky分解
1104
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
