Title: 非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (II)
姊妹博文
非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (I)
非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (II) ⟵ \longleftarrow ⟵ 本篇
非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (III)
非线性最小二乘问题的数值方法 —— 从牛顿迭代法到高斯-牛顿法 (IV)
↑ \uparrow ↑ 理论部分
↓ \downarrow ↓ 实例部分
目录
I.引言
II.牛顿迭代法解非线性方程
1.一元非线性方程形式的牛顿迭代法
2.高维非线性方程组形式的牛顿迭代法
3.牛顿迭代法的雅可比矩阵
III.牛顿迭代法解非线性最小二乘问题
1.从方程问题到最小二乘问题的演化
2.最小二乘问题 Jacobian 矩阵的推导
3.最小二乘问题 Jacobian 矩阵的性质
A.维度
B.整体 Hessian 矩阵
C.局部 Hessian 矩阵
D.Jacobian 矩阵/Hessian 矩阵
E.Hessian 矩阵的对称性
IV.高斯-牛顿法解非线性最小二乘问题
1.高斯-牛顿法的获得
2.高斯-牛顿法的优势
3.高斯-牛顿法的解读 —— 优化观点
V.最小二乘法与高斯的贡献
VI.总结
参考文献
II. 牛顿迭代法解非线性方程
1. 一元非线性方程形式的牛顿迭代法
第一次学牛顿迭代法是在大学的 “数值计算方法” 课上, 用来求非线性函数的根/零点, 也就是对应的非线性方程的解[1].
该方法的原理如 Fig. 1 所示. 牛顿迭代法的最终目的是求函数 y = f ( x ) y=f(x) y=f(x) 的零点 (假设零点存在的情况下), 即满足 f ( x ∗ ) = 0 f(x^{\ast})=0 f(x∗)=0 的 x ∗ x^{\ast} x∗.
初始:
因为不知道 x ∗ x^{\ast} x∗ 的确切值, 先猜测一个初始值 x [ 0 ] x_{[0]} x[0], 显然 ∣ f ( x [ 0 ] ) ∣ > ε \left|f(x_{[0]})\right|>\varepsilon f(x[0]) >ε (其中 ε \varepsilon ε 为一个判断函数值接近于 0 的门限值).
迭代1:
作函数
f
(
x
)
f(x)
f(x) 在
x
[
0
]
x_{[0]}
x[0] 处的切线, 切线的斜率为
f
′
(
x
[
0
]
)
f^{'}(x_{[0]})
f′(x[0]), 则切线方程可以写成
y
−
f
(
x
[
0
]
)
x
−
x
[
0
]
=
f
′
(
x
[
0
]
)
(II-1-1)
\frac{y-f(x_{[0]})}{x-x_{[0]}} = f^{'}(x_{[0]}) \tag{II-1-1}
x−x[0]y−f(x[0])=f′(x[0])(II-1-1)
当
y
=
0
y=0
y=0 时, 该切线与 x 轴相交于
x
[
1
]
x_{[1]}
x[1],
x
[
1
]
=
x
[
0
]
−
f
(
x
[
0
]
)
f
′
(
x
[
0
]
)
(II-1-2)
x_{[1]} = x_{[0]} - \frac{f(x_{[0]})}{f^{'}(x_{[0]})} \tag{II-1-2}
x[1]=x[0]−f′(x[0])f(x[0])(II-1-2)
每获得一个 x 轴上的交点, 都判断该交点是否足够接近零点.
如果如下判别式成立, 则
x
[
1
]
x_{[1]}
x[1] 就近似看做是零点值, 完成工作.
∣
f
(
x
[
1
]
)
∣
<
ε
(II-1-3)
\left| f(x_{[1]}) \right |< \varepsilon \tag{II-1-3}
f(x[1])
<ε(II-1-3)
如判别式不成立则继续迭代计算.
迭代2:
再以新得到的
x
[
1
]
x_{[1]}
x[1] 作为初始点, 作函数
f
(
x
)
f(x)
f(x) 在
x
[
1
]
x_{[1]}
x[1] 处的切线, 并求该切线与 x 轴的交点
x
[
2
]
x_{[2]}
x[2]
x
[
2
]
=
x
[
1
]
−
f
(
x
[
1
]
)
f
′
(
x
[
1
]
)
(II-1-4)
x_{[2]} = x_{[1]} - \frac{f(x_{[1]})}{f^{'}(x_{[1]})} \tag{II-1-4}
x[2]=x[1]−f′(x[1])f(x[1])(II-1-4)
同样需要验证判别式 (类似式 (II-1-3)), 以确定是否获得目标值.
迭代3:
如此循环迭代,
x
[
n
]
=
x
[
n
−
1
]
−
f
(
x
[
n
−
1
]
)
f
(
x
[
n
−
1
]
′
)
(II-1-5)
x_{[n]} = x_{[n-1]} - \frac{f(x_{[n-1]})}{f(x^{'}_{[n-1]})} \tag{II-1-5}
x[n]=x[n−1]−f(x[n−1]′)f(x[n−1])(II-1-5)
得到趋于
x
∗
x^{\ast}
x∗ 的序列,
x
[
0
]
→
x
[
1
]
→
x
[
2
]
→
⋯
→
x
∗
x_{[0]} \rightarrow x_{[1]} \rightarrow x_{[2]} \rightarrow \cdots \rightarrow x^{\ast}
x[0]→x[1]→x[2]→⋯→x∗.
结束:
只要中间迭代步满足
∣
f
(
x
[
n
]
)
∣
<
ε
(II-1-6)
\left| f(x_{[n]}) \right| < \varepsilon \tag{II-1-6}
f(x[n])
<ε(II-1-6)
就将
x
[
n
]
x_{[n]}
x[n] 作为
x
∗
x^{\ast}
x∗ 的近似值, 并结束求解过程.
2. 高维非线性方程组形式的牛顿迭代法
函数切线方程式 (II-1-1) 可以写成等价形式
f
(
x
)
≈
f
(
x
[
i
]
)
+
f
′
(
x
[
i
]
)
(
x
−
x
[
i
]
)
(II-2-1)
f(x) \approx f(x_{[i]})+ f^{'}(x_{[i]}) ({x-x_{[i]}}) \tag{II-2-1}
f(x)≈f(x[i])+f′(x[i])(x−x[i])(II-2-1)
这是一元函数在
x
[
i
]
x_{[i]}
x[i] 处的一阶泰勒近似 (A First-Order Taylor Approximation).
一元方程推广到非线性方程组为
f
(
x
)
=
0
(II-2-2)
\mathbf{f}(\mathbf{x}) = \mathbf{0} \tag{II-2-2}
f(x)=0(II-2-2)
其中
f
=
[
f
1
(
x
1
,
x
2
,
⋯
,
x
n
)
f
2
(
x
1
,
x
2
,
⋯
,
x
n
)
⋮
f
m
(
x
1
,
x
2
,
⋯
,
x
n
)
]
\mathbf{f}=\begin{bmatrix}f_1(x_1,x_2,\cdots,x_n)\\ f_2(x_1,x_2,\cdots,x_n) \\ \vdots \\ f_m(x_1,x_2,\cdots,x_n)\end{bmatrix}
f=
f1(x1,x2,⋯,xn)f2(x1,x2,⋯,xn)⋮fm(x1,x2,⋯,xn)
,
x
=
[
x
1
x
2
⋮
x
n
]
\mathbf{x} = \begin{bmatrix}x_1 \\ x_2\\ \vdots \\x_n \end{bmatrix}
x=
x1x2⋮xn
.
相应地, 推广到高维的一阶泰勒近似为
f
(
x
)
≈
f
(
x
[
i
]
)
+
∂
f
(
x
)
∂
x
∣
x
[
i
]
(
x
−
x
[
i
]
)
(II-2-3)
\mathbf{f}(\mathbf{x}) \approx \mathbf{f}(\mathbf{x}_{[i]})+ \left.\frac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}\right|_{\mathbf{x}_{[i]}} (\mathbf{x} - \mathbf{x}_{[i]}) \tag{II-2-3}
f(x)≈f(x[i])+∂x∂f(x)
x[i](x−x[i])(II-2-3)
其中
x
[
i
]
\mathbf{x}_{[i]}
x[i] 为
x
\mathbf{x}
x 的第 i 步迭代值, 写成分量形式为
x
[
i
]
=
[
x
1
i
,
x
2
i
,
⋯
,
x
n
i
]
T
\mathbf{x}_{[i]} = \begin{bmatrix} x_1^i, x_2^i, \cdots, x_n^i\end{bmatrix}^{\rm\small T}
x[i]=[x1i,x2i,⋯,xni]T.
取
f
(
x
)
=
0
\mathbf{f}(\mathbf{x}) = \mathbf{0}
f(x)=0 时, 得到迭代值
x
[
i
+
1
]
=
x
[
i
]
−
(
∂
f
(
x
)
∂
x
∣
x
[
i
]
)
+
f
(
x
[
i
]
)
(II-2-4)
\mathbf{x}_{[i+1]} = \mathbf{x}_{[i]} - \left(\left.\frac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}\right|_{\mathbf{x}_{[i]}} \right)^{+}\,\mathbf{f}(\mathbf{x}_{[i]}) \tag{II-2-4}
x[i+1]=x[i]−(∂x∂f(x)
x[i])+f(x[i])(II-2-4)
其中
(
⋅
)
+
(\cdot)^{+}
(⋅)+ 为矩阵伪逆 (The Matrix PseudoInverse).
以上即为牛顿迭代法的高维推广.
需要说明, 不管是一维情况还是高维情况下的牛顿迭代法, 都会根据实际情况而有各种变形实现以及分类讨论.
我们此处只是给出了原理性示意, 仅适用于最理想情况.
3. 牛顿迭代法的雅可比矩阵
定义
f
(
x
)
\mathbf{f}(\mathbf{x})
f(x) 在
x
[
i
]
\mathbf{x}_{[i]}
x[i] 处的雅可比矩阵 (The Jacobian Matrix)
J
(
x
[
i
]
)
≜
∂
f
(
x
)
∂
x
∣
x
[
i
]
=
[
∂
f
1
∂
x
1
∂
f
1
∂
x
2
⋯
∂
f
1
∂
x
n
∂
f
2
∂
x
1
∂
f
2
∂
x
2
⋯
∂
f
2
∂
x
n
⋮
⋮
⋱
⋮
∂
f
m
∂
x
1
∂
f
m
∂
x
2
⋯
∂
f
m
∂
x
n
]
x
=
x
[
i
]
(II-3-1)
\mathbf{J}(\mathbf{x}_{[i]}) \triangleq \left. \frac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}\right|_{\mathbf{x}_{[i]}} = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} &\frac{\partial f_1}{\partial x_2} &\cdots & \frac{\partial f_1}{\partial x_n}\\ \frac{\partial f_2}{\partial x_1} &\frac{\partial f_2}{\partial x_2} &\cdots & \frac{\partial f_2}{\partial x_n}\\ \vdots &\vdots &\ddots &\vdots \\ \frac{\partial f_m}{\partial x_1} &\frac{\partial f_m}{\partial x_2} &\cdots & \frac{\partial f_m}{\partial x_n}\\ \end{bmatrix}_{\mathbf{x}=\mathbf{x}_{[i]}} \tag{II-3-1}
J(x[i])≜∂x∂f(x)
x[i]=
∂x1∂f1∂x1∂f2⋮∂x1∂fm∂x2∂f1∂x2∂f2⋮∂x2∂fm⋯⋯⋱⋯∂xn∂f1∂xn∂f2⋮∂xn∂fm
x=x[i](II-3-1)
(关于向量求导可参见机器人中常用矩阵等式-I (汇总))
则式 (II-2-4) 可以写为
x
[
i
+
1
]
=
x
[
i
]
−
J
(
x
[
i
]
)
+
f
(
x
[
i
]
)
(II-3-2)
\mathbf{x}_{[i+1]} = \mathbf{x}_{[i]} - \mathbf{J}(\mathbf{x}_{[i]})^{+}\,\mathbf{f}(\mathbf{x}_{[i]}) \tag{II-3-2}
x[i+1]=x[i]−J(x[i])+f(x[i])(II-3-2)
我们知道矩阵伪逆存在着左逆和右逆 (参考四足机器人中不同优先级任务的执行——Null-Space Projection方法), 考虑到本博客围绕着最优化问题, 数据维度
m
m
m
>
>
> 参数维度
n
n
n, 所以此处的伪逆为左逆, 即
J
+
=
(
J
T
J
)
−
1
J
T
(II-3-3)
\mathbf{J}^{+} = \left( \mathbf{J}^{\small\rm T} \mathbf{J} \right)^{\rm\small {-1}} \mathbf{J}^{\small\rm T} \tag{II-3-3}
J+=(JTJ)−1JT(II-3-3)
故式 (II-3-2) 可进一步写为
x
[
i
+
1
]
=
x
[
i
]
−
[
J
(
x
[
i
]
)
T
J
(
x
[
i
]
)
]
−
1
J
(
x
[
i
]
)
T
f
(
x
[
i
]
)
(II-3-4)
\mathbf{x}_{[i+1]} = \mathbf{x}_{[i]} - \left[ \mathbf{J}(\mathbf{x}_{[i]})^{\small\rm T} \mathbf{J}(\mathbf{x}_{[i]}) \right]^{\rm\small {-1}} \mathbf{J}(\mathbf{x}_{[i]})^{\small\rm T} \mathbf{f}(\mathbf{x}_{[i]}) \tag{II-3-4}
x[i+1]=x[i]−[J(x[i])TJ(x[i])]−1J(x[i])Tf(x[i])(II-3-4)