1 📖 概念入门:最小二乘估计(LSE)
最小二乘估计,简写为LSE(Leart Squares Estimate)
1-1 🔖 最小二乘模型的引出
线性模型如下:
Y
=
A
X
+
ε
\mathbf{Y}=\mathbf{A} \boldsymbol{X}+\boldsymbol{\varepsilon}
Y=AX+ε
其中:
Y = ( y 1 , y 2 , … , y n ) T \mathbf{Y}=\left(y_{1}, y_{2}, \ldots, y_{n}\right)^{T} Y=(y1,y2,…,yn)T :为观测向量
A \bf{A} A:为一个 n × n {n}\times{n} n×n的矩阵
ε = ( ε 1 , ε 2 , … , ε n ) T \varepsilon=\left(\varepsilon_{1}, \varepsilon_{2}, \ldots, \varepsilon_{n}\right)^{T} ε=(ε1,ε2,…,εn)T :为随机误差向量
注意:
上面式常常写做
f ( X ) = A X − b f(X)=\mathbf{A} {X}-b f(X)=AX−b
- b b b:这里相当于上面的 Y \bf{Y} Y
如果
(
Y
−
A
X
L
S
E
)
T
(
Y
−
A
X
L
S
E
)
=
min
X
(
Y
−
A
X
)
T
(
Y
−
A
X
)
\left(\mathbf{Y}-\mathbf{A} \boldsymbol{X}_{L S E}\right)^{T}\left(\mathbf{Y}-\mathbf{A} \boldsymbol{X}_{L S E}\right)=\min _{\boldsymbol{X}}(\mathbf{Y}-\mathbf{A} \boldsymbol{X})^{T}(\mathbf{Y}-\mathbf{A} \boldsymbol{X})
(Y−AXLSE)T(Y−AXLSE)=Xmin(Y−AX)T(Y−AX)
称
X
L
S
E
\boldsymbol{X}_{L S E}
XLSE 是
X
\boldsymbol{X}
X的线性估计,LSE的意思是最小二乘估计(Leart Squares Estimate)
1-2 🔖 最小二乘模型求解(假设A满秩)
下面开始求解
β
L
S
E
\boldsymbol{\beta}_{L S E}
βLSE,首先设:
Q
(
X
)
=
(
Y
−
A
X
)
T
(
(
Y
−
A
X
)
Q(\boldsymbol{X})=(\mathbf{Y}-\mathbf{A} \boldsymbol{X})^{T}((\mathbf{Y}-\mathbf{A} \boldsymbol{X})
Q(X)=(Y−AX)T((Y−AX)
求LSE等价于求
Q
(
X
)
Q(\boldsymbol{X})
Q(X)的最小值,令
∂
Q
(
X
)
∂
X
=
−
2
A
T
Y
+
2
A
T
A
X
=
0
\frac{\partial Q(\boldsymbol{X})}{\partial \boldsymbol{X}}=-2 \mathbf{A}^{T} \mathbf{Y}+2 \mathbf{A}^{T} \mathbf{A} \boldsymbol{X}=0
∂X∂Q(X)=−2ATY+2ATAX=0
可得:
A
T
A
X
=
A
T
Y
\mathbf{A}^{T} \mathbf{A} \boldsymbol{X}=\mathbf{A}^{T} \mathbf{Y}
ATAX=ATY
这个方程组叫做正规方程(normal equations)组。
设
A
\bf{A}
A是列满秩矩阵,则正规方程的解唯一:
X
L
S
E
=
(
A
T
A
)
−
1
A
T
Y
\boldsymbol{X}_{L S E}=\left(\mathbf{A}^{T} \mathbf{A}\right)^{-1} \mathbf{A}^{T} \mathbf{Y}
XLSE=(ATA)−1ATY
显然,
X
L
S
E
\boldsymbol{X}_{L S E}
XLSE是
Y
\bf{Y}
Y的线性函数。
1-3 🔖 高斯马尔科夫定理介绍
首先插入一下什么是线性估计和最佳线性无偏估计:
线性估计:若估计量是观测值的线性函数,则称它为线性估计。
最佳线性无偏估计(BLUE, Best Linear unbiased estimator):设 θ ^ \hat{\theta} θ^是参数 θ {\theta} θ 的线性无偏估计。如果对 θ {\theta} θ的任意一个线性无偏估计 θ ∗ \theta^{*} θ∗,有 var ( θ ∗ ) ≥ var ( θ ^ ) \operatorname{var}\left(\theta^{*}\right) \geq \operatorname{var}(\hat{\theta}) var(θ∗)≥var(θ^),则称 θ ^ \hat{\theta} θ^是参数 θ {\theta} θ的 BLUE 。
在统计学中,高斯-马尔可夫定理是指:
在线性回归模型中,如果误差满足零均值、同方差且互不相关,则回归系数的最佳线性无偏估计(BLUE, Best Linear Unbiased Estimator)就是普通最小二乘法估计。
高斯-马尔可夫定理的假设条件是:
-
零均值: E ( ε ) = 0 \mathrm{E}(\varepsilon)=0 E(ε)=0
-
同方差且不相关: var ( ε ) = E ( ε ε T ) = σ 2 I n \operatorname{var}(\varepsilon)=\mathrm{E}\left(\varepsilon \varepsilon^{T}\right)=\sigma^{2} \mathbf{I}_{\mathrm{n}} var(ε)=E(εεT)=σ2In (其中 I n \bf{I}_{n} In为n阶单位矩阵(Identity Matrix)。)
2 📖 雅克比矩阵、海森矩阵与非线性最小二乘间的关系
2-1 🔖 雅克比矩阵
雅可比矩阵:
[
∂
y
1
∂
x
1
⋯
∂
y
1
∂
x
n
⋮
⋱
⋮
∂
y
m
∂
x
1
⋯
∂
y
m
∂
x
n
]
\left[\begin{array}{ccc}\frac{\partial y_{1}}{\partial x_{1}} & \cdots & \frac{\partial y_{1}}{\partial x_{n}} \\ \vdots & \ddots & \vdots \\ \frac{\partial y_{m}}{\partial x_{1}} & \cdots & \frac{\partial y_{m}}{\partial x_{n}}\end{array}\right]
⎣⎢⎡∂x1∂y1⋮∂x1∂ym⋯⋱⋯∂xn∂y1⋮∂xn∂ym⎦⎥⎤
此矩阵表示为:
J
(
x
1
,
…
,
x
n
)
J\left(x_{1}, \ldots, x_{n}\right)
J(x1,…,xn)
或者
∂
(
y
1
,
…
,
y
m
)
∂
(
x
1
,
…
,
x
n
)
\frac{\partial\left(y_{1}, \ldots, y_{m}\right)}{\partial\left(x_{1}, \ldots, x_{n}\right)}
∂(x1,…,xn)∂(y1,…,ym)
2-2 🔖 海森矩阵
在数学中,海森矩阵(Hessian matrix 或 Hessian)是一个自变量为向量的实值函数的二阶偏导数组成的方块矩阵,此函数如下:
f
(
x
1
,
x
2
,
…
,
x
n
)
f\left(x_{1}, x_{2}, \ldots, x_{n}\right)
f(x1,x2,…,xn)
如果 f 所有的二阶导数都存在,那么 f 的海森矩阵即:
H
(
f
)
i
j
(
x
)
=
D
i
D
j
f
(
x
)
H(f)_{i j}(x)=D_{i} D_{j} f(x)
H(f)ij(x)=DiDjf(x)
其中
x
=
(
x
1
,
x
2
,
…
,
x
n
)
T
x=\left(x_{1}, x_{2}, \ldots, x_{n}\right)^{T}
x=(x1,x2,…,xn)T,即
H
(
f
)
=
[
∂
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(f)=\left[\begin{array}{cccc} \frac{\partial^{2} f}{\partial x_{1}^{2}} & \frac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{1} \partial x_{n}} \\ \frac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{2}^{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{2} \partial x_{n}} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} f}{\partial x_{n} \partial x_{1}} & \frac{\partial^{2} f}{\partial x_{n} \partial x_{2}} & \cdots & \frac{\partial^{2} f}{\partial x_{n}^{2}} \end{array}\right]
H(f)=⎣⎢⎢⎢⎢⎢⎡∂x12∂2f∂x2∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f∂x22∂2f⋮∂xn∂x2∂2f⋯⋯⋱⋯∂x1∂xn∂2f∂x2∂xn∂2f⋮∂xn2∂2f⎦⎥⎥⎥⎥⎥⎤
2-3 🔖 解非线性最小二乘
一、最小二乘问题
min
x
∈
R
′
f
(
x
)
=
min
x
∈
R
′
∑
i
=
1
m
f
i
2
(
x
)
\min _{x \in R^{\prime}} f(x)=\min _{x \in R^{\prime}} \sum_{i=1}^{m} f_{i}^{2}(x)
x∈R′minf(x)=x∈R′mini=1∑mfi2(x)
这里
f
i
(
x
)
(
i
=
1
,
2
,
⋯
,
m
)
f_{i}(x)(i=1,2, \cdots, m)
fi(x)(i=1,2,⋯,m) 可理解为误差, 上式优化问题就是要使误差的平方和最小(即最小二乘问题)。
- 当 f i ( x ) ( i = 1 , 2 , ⋯ , m ) f_{i}(x)(i=1,2, \cdots, m) fi(x)(i=1,2,⋯,m) 都是线性函数时,上式问题是线性最小二乘问题;
- 当 f i ( x ) ( i = 1 , 2 , ⋯ , m ) f_{i}(x)(i=1,2, \cdots, m) fi(x)(i=1,2,⋯,m)有非线性函数时,上式就是非线性最小二乘问题。
二、线性最小二乘问题
f
i
(
x
)
=
p
i
T
x
−
b
i
(
i
=
1
,
2
,
⋯
,
m
)
A
=
[
p
1
T
⋮
p
m
T
]
,
b
=
[
b
1
⋮
b
m
]
\begin{array}{l} f_{i}(x)=p_{i}^{T} x-b_{i} \quad(i=1,2, \cdots, m) \\ A=\left[\begin{array}{c} p_{1}^{T} \\ \vdots \\ p_{m}^{T} \end{array}\right], b=\left[\begin{array}{c} b_{1} \\ \vdots \\ b_{m} \end{array}\right] \end{array}
fi(x)=piTx−bi(i=1,2,⋯,m)A=⎣⎢⎡p1T⋮pmT⎦⎥⎤,b=⎣⎢⎡b1⋮bm⎦⎥⎤
则
f
(
x
)
=
(
A
x
−
b
)
T
(
A
x
−
b
)
=
x
T
A
T
A
x
−
2
b
T
A
x
+
b
T
b
f(x)=(A x-b)^{T}(A x-b)=x^{T} A^{T} A x-2 b^{T} A x+b^{T} b
f(x)=(Ax−b)T(Ax−b)=xTATAx−2bTAx+bTb
令
∇
f
(
x
)
=
2
A
T
A
x
−
2
A
T
b
=
0
\nabla f(x)=2 A^{T} A x-2 A^{T} b=0
∇f(x)=2ATAx−2ATb=0
得
A
T
A
x
=
A
T
b
A^{T} A x=A^{T} b \quad
ATAx=ATb
上面式子称为法方程组
若
A
\mathrm{A}
A 列满秩, 则
A
T
A
\mathrm{A}^{\mathrm{T}} \mathrm{A}
ATA 为正定对称矩阵, 从而可逆,得到, 线性最小二成问题的全局最优解
x
∗
=
(
A
T
A
)
−
1
A
T
b
x^{*}=\left(A^{T} A\right)^{-1} A^{T} b
x∗=(ATA)−1ATb
三、非线性最小二乘问题
由于的
f
i
(
x
)
\mathrm{f}_{\mathbf{i}}(\mathbf{x})
fi(x) 非线性性, 此时按梯度等于0得到的是一个非线性方程组, 求解困难。常用的基本思想是用一些列最小二乘问题求解该非线性是小二乘问题:
设 x ( k ) \mathrm{x}^{(\mathrm{k})} x(k) 是解的第k次近似,在 x ( k ) \mathbf{x}^{(\mathrm{k})} x(k) 处将函数 f i ( x ) \mathbf{f}_{\mathbf{i}}(\mathbf{x}) fi(x) 线性化, 把问题化为线性最小二乘问题, 求出第 k + 1 \mathbf{k}+1 k+1 次近似解 x ( k + 1 ) \mathrm{x}^{(k+1)} x(k+1); 再从 x ( k + 1 ) \mathrm{x}^{(k+1)} x(k+1)出发, 重复此过程, 直到达到迭代终止难则。
2-4 🔖 牛顿法
注意:这里雅克比 J = A J=A J=A,因为这里的雅克比 J J J是相对于 f = A X f=AX f=AX来说的,并不是针对最小二乘来说的
为了能够方便理解整个流程, 我们还是拿一元函数 f ( x ) f(x) f(x) 来举例。
我们首先对其在
x
=
x
0
x=x_{0}
x=x0 处进行二阶泰勒展开:
f
(
x
)
=
f
(
x
0
)
+
f
i
(
x
0
)
(
x
−
x
0
)
+
f
′
′
(
x
0
)
2
(
x
−
x
0
)
2
+
o
(
x
−
x
0
)
2
f(x)=f\left(x_{0}\right)+f^{i}\left(x_{0}\right)\left(x-x_{0}\right)+\frac{f^{\prime \prime}\left(x_{0}\right)}{2}\left(x-x_{0}\right)^{2}+o\left(x-x_{0}\right)^{2}
f(x)=f(x0)+fi(x0)(x−x0)+2f′′(x0)(x−x0)2+o(x−x0)2
其中, 由于泰勒展开的特性,后面
o
(
x
−
x
0
)
o\left(x-x_{0}\right)
o(x−x0) 部分不予考虑,我们只考虑前面展开部分的极值问题:
g
(
x
)
=
f
(
x
0
)
+
f
′
(
x
0
)
(
x
−
x
0
)
+
f
′
′
(
x
0
)
2
(
x
−
x
0
)
2
g(x)=f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)+\frac{f^{\prime \prime}\left(x_{0}\right)}{2}\left(x-x_{0}\right)^{2}
g(x)=f(x0)+f′(x0)(x−x0)+2f′′(x0)(x−x0)2
上面的式子是一个一元二次函数,那么其极值就是一阶导数为0的时候,我们可以先微分:
g
′
(
x
)
=
f
′
(
x
0
)
+
f
′
′
(
x
0
)
(
x
−
x
0
)
g^{\prime}(x)=f^{\prime}\left(x_{0}\right)+f^{\prime \prime}\left(x_{0}\right)\left(x-x_{0}\right)
g′(x)=f′(x0)+f′′(x0)(x−x0)
令
g
′
(
x
1
)
=
0
g^{\prime}\left(x_{1}\right)=0
g′(x1)=0, 则此时的
x
1
x_{1}
x1 就是极值(为了方便说明,暂不考虑鞍点的情况)。故
f
′
(
x
0
)
+
f
′
′
(
x
0
)
(
x
1
−
x
0
)
=
0
x
1
=
x
0
−
f
′
(
x
0
)
f
′
′
(
x
0
)
f^{\prime}\left(x_{0}\right)+f^{\prime \prime}\left(x_{0}\right)\left(x_{1}-x_{0}\right)=0 \\ x_{1}=x_{0}-\frac{f^{\prime}\left(x_{0}\right)}{f^{\prime \prime}\left(x_{0}\right)}
f′(x0)+f′′(x0)(x1−x0)=0x1=x0−f′′(x0)f′(x0)
以此类推,可以得到迭代公式:
x
n
+
1
=
x
n
−
f
′
(
x
n
)
f
′
′
(
x
n
)
x_{n+1}=x_{n}-\frac{f^{\prime}\left(x_{n}\right)}{f^{\prime \prime}\left(x_{n}\right)}
xn+1=xn−f′′(xn)f′(xn)
根据这个方法就可以不断的迭代下去直到收敘,最终找到极值了。
如果是多元的情况, 则一阶导数
f
′
(
x
)
f^{\prime}(x)
f′(x) 被叫做梯度, 也称之为雅可比矩阵
J
J
J (这里不太严谨。 严格来说, 矩阵的梯度为一阶导的转置, 函数的梯度为一阶导,这里并没有进行详细的区分),二 阶导数矩阵
f
′
′
(
x
)
f^{\prime \prime}(x)
f′′(x), 也被叫做海塞矩阵
H
H
H 。 如果是收敘的话,
Δ
x
=
x
n
+
1
−
x
n
≈
0
\Delta x=x_{n+1}-x_{n} \approx 0
Δx=xn+1−xn≈0, 则式子 可以转化为:
Δ
x
=
−
f
t
(
x
n
)
f
′
′
(
x
n
)
=
−
J
H
\Delta x=-\frac{f^{t}\left(x_{n}\right)}{f^{\prime \prime}\left(x_{n}\right)}=-\frac{J}{H}
Δx=−f′′(xn)ft(xn)=−HJ
也就是说:
H
Δ
x
=
−
J
H \Delta x=-J
HΔx=−J
这样,就可以求出能够取得函数
f
(
x
)
f(x)
f(x) 的极值点, 继而算出函数
f
(
x
)
f(x)
f(x) 的极值。
由于牛顿法需要算二阶导数,如果高阶的话,需要算海塞矩阵,这里是有三个缺陷:
- 要求给定的方程需要二阶可导
- 非凸函数的海森矩阵不一定有逆
- 数据较大的时候,海塞矩阵的计算量偏大
因此,需要思考别的方法来进行最小二乘问题的优化和求解。
2-5 🔖 高斯牛顿法
https://zhuanlan.zhihu.com/p/113946848
注意:这里雅克比 J = A J=A J=A,因为这里的雅克比 J J J是相对于 f = A X f=AX f=AX来说的,并不是针对最小二乘来说的
如果代入到最小二乘问题中,牛顿法和梯度下降法都是针对目标函数
F
(
x
k
)
F\left(x_{k}\right)
F(xk)来进行求解的,这样,就不可避免的需要求得海塞矩阵
H
H
H,所以,为了避免这个问题,我们选取了误差函数
f
(
x
)
f(x)
f(x) 来进行优化求解:
min
F
(
x
)
=
1
2
∥
f
(
x
)
∥
2
2
\min F(x)=\frac{1}{2}\|f(x)\|_{2}^{2}
minF(x)=21∥f(x)∥22
那么,我们从上面的迭代步骤2中可以看到:
1、给定某个初始值 x 0 x_{0} x0;
2、对于第 k k k 次迭代,寻找一个增量 Δ x k \Delta x_{k} Δxk, 使得误差 ∥ f ( x k + Δ x k ) ∥ 2 2 \left\|f\left(x_{k}+\Delta x_{k}\right)\right\|_{2}^{2} ∥f(xk+Δxk)∥22 达到极小值
3、若 Δ x k \Delta x_{k} Δxk 足够小,则停止迭代 ;
4、否则, 令 x k + 1 = x k + Δ x k x_{k+1}=x_{k}+\Delta x_{k} xk+1=xk+Δxk, 返回第2步。
那么,我们对
f
(
x
+
Δ
x
)
f(x+\Delta x)
f(x+Δx) 进行一阶泰勒展开。
f
(
x
+
Δ
x
)
≈
f
(
x
)
+
J
(
x
)
T
Δ
x
+
o
(
Δ
x
)
f(x+\Delta x) \approx f(x)+J(x)^{T} \Delta x+o(\Delta x)
f(x+Δx)≈f(x)+J(x)TΔx+o(Δx)
我们需要求
Δ
x
\Delta x
Δx 使得上面的式子
∥
f
(
x
+
Δ
x
)
∥
2
2
\|f(x+\Delta x)\|_{2}^{2}
∥f(x+Δx)∥22 有最小值, 所以,我们可以得到最小二乘问题 为:
Δ
x
∗
=
arg
min
1
2
∥
f
(
x
+
Δ
x
)
∥
2
2
≈
arg
min
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
2
2
\Delta x^{*}=\arg \min \frac{1}{2}\|f(x+\Delta x)\|_{2}^{2} \approx \arg \min \frac{1}{2}\left\|f(x)+J(x)^{T} \Delta x\right\|_{2}^{2}
Δx∗=argmin21∥f(x+Δx)∥22≈argmin21∥∥f(x)+J(x)TΔx∥∥22
为了求极值,对其求导:
m
(
x
)
=
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
2
=
1
2
(
f
(
x
)
+
J
(
x
)
T
Δ
x
)
T
(
f
(
x
)
+
J
(
x
)
T
Δ
x
)
=
1
2
(
∥
f
(
x
)
∥
2
+
2
f
(
x
)
J
(
x
)
T
Δ
x
+
Δ
x
T
J
(
x
)
J
(
x
)
T
Δ
x
m(x)=\frac{1}{2}\left\|f(x)+J(x)^{T} \Delta x\right\|^{2}=\frac{1}{2}\left(f(x)+J(x)^{T} \Delta x\right)^{T}\left(f(x)+J(x)^{T} \Delta x\right) \\ =\frac{1}{2}\left(\|f(x)\|^{2}+2 f(x) J(x)^{T} \Delta x+\Delta x^{T} J(x) J(x)^{T} \Delta x\right.
m(x)=21∥∥f(x)+J(x)TΔx∥∥2=21(f(x)+J(x)TΔx)T(f(x)+J(x)TΔx)=21(∥f(x)∥2+2f(x)J(x)TΔx+ΔxTJ(x)J(x)TΔx
故,对其求导可以得到:
m
′
(
x
)
=
J
(
x
)
f
(
x
)
+
J
(
x
)
J
(
x
)
T
Δ
x
m^{\prime}(x)=J(x) f(x)+J(x) J(x)^{T} \Delta x
m′(x)=J(x)f(x)+J(x)J(x)TΔx
则,此时可以转化为线性求解问题:
m
′
(
x
)
=
0
→
J
(
x
)
J
(
x
)
T
Δ
x
=
−
J
(
x
)
f
(
x
)
m^{\prime}(x)=0 \quad \rightarrow \quad J(x) J(x)^{T} \Delta x=-J(x) f(x)
m′(x)=0→J(x)J(x)TΔx=−J(x)f(x)
令
J
(
x
)
J
(
x
)
T
J(x) J(x)^{T}
J(x)J(x)T 定义为
H
(
x
)
H(x)
H(x), 令
−
J
(
x
)
f
(
x
)
-J(x) f(x)
−J(x)f(x) 定义为
g
(
x
)
g(x)
g(x), 则此时变为了:
H
Δ
x
=
g
s.t
H
=
J
J
T
&
g
=
−
J
f
H \Delta x=g \quad \text { s.t } \quad H=J J^{T} \quad \& \quad g=-J f
HΔx=g s.t H=JJT&g=−Jf
这样,就可以优化求解了。上面的最小二乘的优化步骤就可以变为 :
1、给定某个初始值 x 0 x_{0} x0;
2、对于第k次迭代,求出当前的雅可比矩阵 J ( x k ) J\left(x_{k}\right) J(xk) 和误差 f ( x k ) f\left(x_{k}\right) f(xk);
3、求解增量方程 : H Δ x k = g H \Delta x_{k}=g HΔxk=g;
4、若 Δ x k \Delta x_{k} Δxk 足够小,则停止迭代 ;
5、否则,令 x k + 1 = x k + Δ x k x_{k+1}=x_{k}+\Delta x_{k} xk+1=xk+Δxk, 返回第2步。
相比较于传统的最小二乘求解方法,只更改了两个步骤。
该方法的优点和缺点如下:
优点:
- 避免了求海塞矩阵,大大减少了计算量。
缺点:- 为了求解 H − H^{-} H−, 需要 H H H 矩阵可逆, 但是实际上 J J T J J^{T} JJT 只有半正定性,所以,当为奇异矩阵的时
候,稳定性较差,算法不收敘。- 如果求出来的步长 Δ x k \Delta x_{k} Δxk 太大,会导致其局部近似不精确,严重的时候,可能无法保证迭代收 敘。
- 容易和梯度下降法一样,陷入钻齿状,导致迭代次数较长。
不过,为了能够更好的进行最小二乘问题的求解,我们可以使用列文伯格-马夸特法 ( L M ) (\mathrm{LM}) (LM) 来进行求 解。
2-6 🔖 LM法
为什么选择LM法?
原因一:如果 J ( x ) J ( x ) T Δ x = − J ( x ) f ( x ) J(x) J(x)^{T} \Delta x=-J(x) f(x) J(x)J(x)TΔx=−J(x)f(x)中的矩阵 J k ⊤ J k J_{k}^{\top} J_{k} Jk⊤Jk是病态或奇异时,线性方程组的求解就会遇到困难, 此时基本的对策是正则化的方法, 即
( J k T J k + α I ) d = − J k T f ( k ) \left(J_{k}^{T} J_{k}+\alpha I\right) d=-J_{k}^{T} f^{(k)} (JkTJk+αI)d=−JkTf(k)
原因二:在高斯牛顿法的缺点中,可以看到,有一点使容易进入锯齿状,导致迭代的次数较长。所以,为了避免其步长过大导致的问题,该方法提出了信赖区域,设定一个区域。使得步长能够受到控制
在更新迭代的过程中,为了判定近似值的好坏,我们设定了一个评判指标:
ρ
=
f
(
x
+
Δ
x
)
−
f
(
x
)
J
(
x
)
T
Δ
x
\rho=\frac{f(x+\Delta x)-f(x)}{J(x)^{T} \Delta x}
ρ=J(x)TΔxf(x+Δx)−f(x)
这个指标就是我们的近似指标,可以看到其分为以下几种情况:
- ρ \rho ρ 接近1, 近似是好的, 不需要更改;
- ρ \rho ρ 太小, 则实际减少的值小于近似减少的值,近似较大,需要缩小近似的范围;
- ρ \rho ρ 太大, 则实际减少的值大于近似减少的值,近似较小,需要扩大近似的范围。
这样的话,就可以动态调整步长了。
通过近似指标,我们可以设定信赖区域的大小 0 _{0} 0 当没有接近我们设定的间值,则不断调整动态区 域,直到找到好的近似结果。
当找到符合要求的近似结果后,就可以进行后续正常的迭代更新了。
因此,使用该信赖区域后,可以更新算法流程:
1、给定某个初始值 x 0 x_{0} x0;
2、对于第 k k k 次迭代,在高斯牛顿法的基础上加入信赖区域:
min 1 2 ∥ f ( x k ) + J ( x k ) T Δ x k ∥ 2 , s.t ∥ D Δ x k ∥ 2 ≤ μ \min \frac{1}{2}\left\|f\left(x_{k}\right)+J\left(x_{k}\right)^{T} \Delta x_{k}\right\|^{2}, \quad \text { s.t } \quad\left\|D \Delta x_{k}\right\|^{2} \leq \mu min21∥∥∥f(xk)+J(xk)TΔxk∥∥∥2, s.t ∥DΔxk∥2≤μ
其中, μ \mu μ 是信赖半径, D D D 为系数矩阵。
3、计算近似指标 ρ \rho ρ :
ρ = f ( x + Δ x ) − f ( x ) J ( x ) T Δ x \rho=\frac{f(x+\Delta x)-f(x)}{J(x)^{T} \Delta x} ρ=J(x)TΔxf(x+Δx)−f(x)
4、根据经验值,设定:
- 若 ρ > 3 4 \rho>\frac{3}{4} ρ>43, 则设置 μ = 2 μ \mu=2 \mu μ=2μ, 跳转第6步 ; ; ;
- 若 ρ < 1 4 \rho<\frac{1}{4} ρ<41, 则设置 μ = 0.5 μ \mu=0.5 \mu μ=0.5μ, 跳转第6步
- 若 ρ \rho ρ 大于设定的间值,则跳转至第5步,求解 Δ x k \Delta x_{k} Δxk, 令 x k + 1 = x k + Δ x k ∘ x_{k+1}=x_{k}+\Delta x_{k \circ} xk+1=xk+Δxk∘
5、求解增量方程 : ( H + λ I ) Δ x k = g (H+\lambda I) \Delta x_{k}=g (H+λI)Δxk=g;
6、若 Δ x k \Delta x_{k} Δxk 足够小,则停止迭代,否则,返回第2步。
至于增量方程的获取,可以通过拉格朗日函数来求解:
min
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
s.t
∣
∣
D
Δ
x
<
μ
∥
2
\min \frac{1}{2}\left\|f(x)+J(x)^{T} \Delta x\right\| \quad \text { s.t } \quad|| D \Delta x<\mu \|_{2}
min21∥∥f(x)+J(x)TΔx∥∥ s.t ∣∣DΔx<μ∥2
构建拉格朗日函数,
λ
\lambda
λ 是系数因子:
L
(
Δ
x
,
λ
)
=
1
2
∥
f
(
x
)
+
J
(
x
)
T
Δ
x
∥
2
+
λ
2
(
∥
D
Δ
x
∥
2
−
μ
)
L(\Delta x, \lambda)=\frac{1}{2}\left\|f(x)+J(x)^{T} \Delta x\right\|^{2}+\frac{\lambda}{2}\left(\|D \Delta x\|^{2}-\mu\right)
L(Δx,λ)=21∥∥f(x)+J(x)TΔx∥∥2+2λ(∥DΔx∥2−μ)
这样的话,化简后求导就可以得到:
J
(
x
)
f
(
x
)
+
J
(
x
)
J
T
(
x
)
Δ
x
+
λ
D
T
D
Δ
x
=
0
J(x) f(x)+J(x) J^{T}(x) \Delta x+\lambda D^{T} D \Delta x=0
J(x)f(x)+J(x)JT(x)Δx+λDTDΔx=0
我们化简后得到:
(
J
J
T
+
λ
D
T
D
)
Δ
x
=
−
J
f
\left(J J^{T}+\lambda D^{T} D\right) \Delta x=-J f
(JJT+λDTD)Δx=−Jf
在本文中, 我们令
H
=
J
J
T
,
g
=
−
J
f
H=J J^{T}, g=-Jf
H=JJT,g=−Jf 在实际使用中, 通常用
I
I
I 来代替
D
T
D
D^{T} D
DTD 。 所以,公式 就变为:
(
H
+
λ
I
)
Δ
x
k
=
g
(H+\lambda I) \Delta x_{k}=g
(H+λI)Δxk=g
这样,就可以得到对应的增量方程了。
代入算法流程中,最终就可以优化得到最小二乘问题的极小值了。
2.7 🔖 几种方法的联系与优缺点
- 最速下降法:直观上将本算法十分简单,直接按照梯度的反方向下降即可;缺点是过于贪心,容易呈锯齿状下降,从而增加迭代次数。
- 牛顿法:相对而言也非常直观,同时由于引入了二阶导数,可以处理一阶导为0的情况;但缺点是二阶导数具有非常大的计算量。
- 高斯牛顿法:在牛顿法的基础上进行了一定程度的简化,使用 J T J J^{T}J JTJ 代替海塞矩阵,避免了二阶导数的计算;缺点在于 J T J J^{T}J JTJ很容易病态,导致无法得到正确的结果。
- LM法:通过引入阻尼项使得 J T J J^{T}J JTJ不那么容易病态,并且可以通过调整阻尼完成在梯度法和牛顿法之间切换;缺点不太清楚。
- 牛顿法: H Δ x = − J \quad H \Delta x=-J \quad HΔx=−J s.t H = f ′ ′ ( x k ) , J = f ′ ( x k ) \quad H=f^{\prime \prime}\left(x_{k}\right), J=f^{\prime}\left(x_{k}\right) H=f′′(xk),J=f′(xk)
- 梯度下降法: Δ x = − γ J \quad \Delta x=-\gamma J \quad Δx=−γJ s.t J = f t ( x k ) \quad J=f^{t}\left(x_{k}\right) J=ft(xk)
- 高斯牛顿法: H Δ x = g H \Delta x=g \quad HΔx=g s.t H = J J T , g = − J f \quad H=J J^{T}, \quad g=-J f H=JJT,g=−Jf
- 列文伯格-马夸特法: ( H + λ I ) Δ x k = g (H+\lambda I) \Delta x_{k}=g \quad (H+λI)Δxk=g s.t H = J J T , g = − J f \quad H=J J^{T}, g=-J f H=JJT,g=−Jf
其实,这四种方法在最小二乘的问题求解中,也是有着联系的。
我们设定最小二乘问题为:
min
F
(
x
)
=
1
2
∥
f
(
x
)
∥
2
2
\min F(x)=\frac{1}{2}\|f(x)\|_{2}^{2}
minF(x)=21∥f(x)∥22
根据针对求解的是目标函数还是误差函数,可以将问题进行分类:
针对目标函数 F ( x ) F(x) F(x)优化
对于目标函数
F
(
x
)
F(x)
F(x) 进行一阶泰勒展开:
F
(
x
k
+
Δ
x
k
)
≈
F
(
x
k
)
+
J
(
x
k
)
T
Δ
x
k
F\left(x_{k}+\Delta x_{k}\right) \approx F\left(x_{k}\right)+J\left(x_{k}\right)^{T} \Delta x_{k}
F(xk+Δxk)≈F(xk)+J(xk)TΔxk
此时,此时变成求最小值的问题, 则:
Δ
x
k
∗
=
arg
min
(
F
(
x
k
)
+
J
(
x
k
)
T
Δ
x
k
)
\Delta x_{k}^{*}=\arg \min \left(F\left(x_{k}\right)+J\left(x_{k}\right)^{T} \Delta x_{k}\right)
Δxk∗=argmin(F(xk)+J(xk)TΔxk)
故,对其求最小值,可以进行求一阶导数为0:
Δ
x
k
=
−
J
(
x
k
)
T
\Delta x_{k}=-J\left(x_{k}\right)^{T}
Δxk=−J(xk)T
可以看到,如果增加一个步长 λ \lambda λ, 此时的方法就是梯度下降法。
如果对目标函数
F
(
x
)
F(x)
F(x) 其进行二阶泰勒展开:
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\left(x_{k}+\Delta x_{k}\right) \approx F\left(x_{k}\right)+J\left(x_{k}\right)^{T} \Delta x_{k}+\frac{1}{2} \Delta x_{k}^{T} H\left(x_{k}\right) \Delta x_{k}
F(xk+Δxk)≈F(xk)+J(xk)TΔxk+21ΔxkTH(xk)Δxk
则,此时的增量方程为最小二乘问题:
Δ
x
k
∗
=
arg
min
(
F
(
x
k
)
+
J
(
x
k
)
T
Δ
x
k
+
1
2
Δ
x
k
T
H
(
x
k
)
Δ
x
k
)
\Delta x_{k}^{*}=\arg \min \left(F\left(x_{k}\right)+J\left(x_{k}\right)^{T} \Delta x_{k}+\frac{1}{2} \Delta x_{k}^{T} H\left(x_{k}\right) \Delta x_{k}\right)
Δxk∗=argmin(F(xk)+J(xk)TΔxk+21ΔxkTH(xk)Δxk)
则,为了求其最小值, 对
Δ
x
k
∗
\Delta x_{k}^{*}
Δxk∗ 进行求导:
J
(
x
k
)
+
H
(
x
k
)
Δ
x
k
=
0
→
H
∗
Δ
x
=
−
J
J\left(x_{k}\right)+H\left(x_{k}\right) \Delta x_{k}=0 \rightarrow H * \Delta x=-J
J(xk)+H(xk)Δxk=0→H∗Δx=−J
则此时的方法为牛顿法。
所以,如果从目标函数 F ( x ) F(x) F(x) 入手的话,梯度下降法和牛顿法就是其一阶泰勒展开和二阶泰勒展开 的求解方法。
针对误差函数 f ( x ) f(x) f(x)优化
对误差函数的优化,可以得到两种方法:高斯牛顿法和列文伯格-马夸特法。
其中,列文伯格-马夸特法是在高斯牛顿法的基础上得来的。他们之间的联系,主要取决于参数
λ
\lambda
λ
- λ \lambda λ 较大的时候, 公式变为: λ I Δ x k = g \lambda I \Delta x_{k}=g λIΔxk=g, 此时为梯度下降法
- λ \lambda λ 较小的时候, 公式变为: H Δ x k = g H \Delta x_{k}=g HΔxk=g, 此时为高斯牛顿法。
所以, L M \bf{L M} LM 法可以在一定程度上避免线性方程组的系数矩阵的非奇异和病态的问题,从而得到更 精准,更稳定的增量 Δ x k \Delta x_{k} Δxk 。
3 📖 矩阵求解
参考链接
https://blog.csdn.net/wangshuailpp/article/details/80209863
https://www.cnblogs.com/AndyJee/p/3846455.html
3-1 🔖 内容一:直接给出 A X = B \bf{AX=B} AX=B解的情况。
(1)R(A)< r(A|B),方程组无解
(2)r(A)=r(A|B)=n,方程组有唯一解
(3)r(A)=r(A|B) < n,方程组有无穷解
(4)r(A)>r(A|B),这种情况不存在
其中r()代表矩阵的秩,A|B是增广矩阵,n是X未知数个数。
- 当A满秩时,为了提高效率通常采用QR分解、LTLD分解、Cholesky分解和SVD分解(奇异分解)等。
- 当A亏秩时,只能使用SVD分解方法,其他方法将失效。
QR分解?LTLD分解、Cholesky分解、SVD分解?
3-2 🔖内容二:线性最小二乘问题
min ∥ A x − b ∥ 2 2 A ∈ R m ∗ n x ∈ R n b ∈ R m \min \|A x-b\|_{2}^{2} \\ A \in R^{m^{*} n} \quad x \in R^{n} \quad b \in R^{m} min∥Ax−b∥22A∈Rm∗nx∈Rnb∈Rm
m个方程求解n个未知数,有三种情况:
- m = n \mathrm{m}=\mathrm{n} m=n 且A为非奇异, 则有唯一解, x = A − 1 b x=A^{-1} b x=A−1b
- m > n m>n m>n, 约束的个数大于未知数的个数, 称为超定问题 (overdetermined)
- m < n \mathrm{m}<\mathrm{n} m<n, 负定/欠定问题 (underdetermined)
通常我们遇到的都是超定问题, 此时A = = = b的解是不存在的,从而转向解最小二乘问题: f ( x ) = ∥ A x − b ∥ f(x)=\|A x-b\| f(x)=∥Ax−b∥ ,其中 f ( x ) \mathrm{f}(\mathrm{x}) f(x) 为凸函数。令一阶导数为 0, 得到: A T A x − A T b = 0 A^{T} A x-A^{T} b=0 ATAx−ATb=0,称之为正规方程
一般解: x = ( A T A ) − 1 A T b x=\left(A^{T} A\right)^{-1} A^{T} b x=(ATA)−1ATb可见一般在视觉SLAM中后端优化就是约束项大于未知数的个数(超定问题), 采用最小二乘问题求 解
4 📖 最小二乘参考链接
https://zhuanlan.zhihu.com/p/113946848
https://zhuanlan.zhihu.com/p/114820251
https://blog.csdn.net/wangshuailpp/article/details/80209863
https://www.cnblogs.com/AndyJee/p/3846455.html
https://blog.csdn.net/wangshuailpp/article/details/80209863
https://www.cnblogs.com/AndyJee/p/3846455.html