数值计算之 最小二乘法(3)最小二乘的矩阵解法
- 前言
- 回顾最小二乘的线性解
- 列满秩矩阵的最小二乘解法
- Cholesky分解求线性最小二乘解
- QR分解求线性最小二乘解
- 亏秩矩阵的最小二乘解法
- SVD分解求亏秩最小二乘解
- 补充1:超定齐次方程组的线性最小二乘解法
- 补充2:欠定行满秩方程组的线性最小二乘解法
- 后记
前言
之前将最小二乘法与线性方程组求解关联,得到了线性最小二乘的矩阵形式,以及线性最小二乘的几何意义。
本篇将介绍线性最小二乘的矩阵解法。
回顾最小二乘的线性解
对于线性超定方程组
A
x
=
b
,
A
∈
R
m
×
n
,
m
>
n
Ax=b,A\in R^{m\times n},m>n
Ax=b,A∈Rm×n,m>n,其最小二乘解可表示为:
arg min
x
∣
∣
A
x
−
b
∣
∣
2
2
=
arg min
x
(
A
x
−
b
)
T
(
A
x
−
b
)
→
A
T
A
x
=
A
T
b
i
f
r
a
n
k
(
A
)
=
n
,
x
=
(
A
T
A
)
−
1
A
T
b
\argmin_x ||Ax-b||^2_2=\argmin_x (Ax-b)^T(Ax-b) \\ \to A^TAx=A^Tb \\ \quad \\ if \quad rank(A)=n, x=(A^TA)^{-1}A^Tb
xargmin∣∣Ax−b∣∣22=xargmin(Ax−b)T(Ax−b)→ATAx=ATbifrank(A)=n,x=(ATA)−1ATb
以上解法具有两个问题:① A T A A^TA ATA求逆运算的效率;② r a n k ( A ) < n rank(A)n,rank(A)=n Ax=b,A∈Rm×n,m>n,rank(A)=n,有 A T A A^TA ATA对称且可逆。回顾矩阵分解的知识,可知 A T A A^TA ATA矩阵能够进行Cholesky分解和QR分解。
Cholesky分解求线性最小二乘解
通过将
A
T
A
A^TA
ATA分解为下三角矩阵
L
L
L的乘积
L
L
T
LL^T
LLT:
A
T
A
=
L
L
T
A^TA=LL^T
ATA=LLT
则可以将复杂的线性方程组求解:
A
T
A
x
=
A
T
b
A^TAx=A^Tb
ATAx=ATb
转化为两个简单的三角线性方程组求解:
L
L
T
x
=
A
T
b
→
L
T
x
=
y
,
L
y
=
A
T
b
LL^Tx=A^Tb \to L^Tx=y, Ly=A^Tb
LLTx=ATb→LTx=y,Ly=ATb
Cholesky分解速度很快,但精度一般,稳定性差。适合在限定时间内的大规模超定线性方程计算求解。
QR分解求线性最小二乘解
对于列满秩矩阵
A
A
A而言,可以唯一分解为正交矩阵与对角元为正的上三角矩阵的乘积:
A
=
Q
R
=
Q
[
R
0
]
A=QR=Q\begin{bmatrix} R \\ 0 \end{bmatrix}
A=QR=Q[R0]
现在考虑最小二乘法的QR分解法:
arg min
x
∣
∣
A
x
−
b
∣
∣
2
2
=
arg min
x
∣
∣
Q
[
R
0
]
x
−
b
∣
∣
=
arg min
x
∣
∣
[
R
0
]
x
−
Q
T
b
∣
∣
2
2
=
arg min
x
∣
∣
[
R
x
0
]
−
[
β
1
β
2
]
∣
∣
2
2
=
arg min
x
∣
∣
R
x
−
β
1
∣
∣
2
2
+
∣
∣
β
2
∣
∣
2
2
⟺
arg min
x
∣
∣
R
x
−
β
1
∣
∣
2
2
→
R
x
=
β
1
,
x
=
R
−
1
β
1
[
β
1
β
2
]
=
Q
T
b
\argmin_x ||Ax-b||^2_2=\argmin_x ||Q\begin{bmatrix} R \\ 0 \end{bmatrix}x-b|| \\ = \argmin_x ||\begin{bmatrix} R \\ 0 \end{bmatrix}x-Q^Tb||^2_2 \\ = \argmin_x ||\begin{bmatrix} Rx \\ 0 \end{bmatrix} - \begin{bmatrix} \beta_1 \\ \beta_2 \end{bmatrix}||^2_2 \\ = \argmin_x ||Rx-\beta_1||^2_2 + ||\beta_2||^2_2 \\ \iff \argmin_x ||Rx-\beta_1||^2_2 \\ \to Rx=\beta_1,x=R^{-1}\beta_1 \\ \quad \\ \begin{bmatrix} \beta_1 \\ \beta_2 \end{bmatrix} = Q^Tb
xargmin∣∣Ax−b∣∣22=xargmin∣∣Q[R0]x−b∣∣=xargmin∣∣[R0]x−QTb∣∣22=xargmin∣∣[Rx0]−[β1β2]∣∣22=xargmin∣∣Rx−β1∣∣22+∣∣β2∣∣22⟺xargmin∣∣Rx−β1∣∣22→Rx=β1,x=R−1β1[β1β2]=QTb
QR分解的速度较快,精度一般。不过由于存在高精度的QR分解方式,因此适合需要高精度解的小规模超定方程组计算。
亏秩矩阵的最小二乘解法
对于线性方程组 A x = b , A ∈ R m × n , m > n , r a n k ( A ) < n Ax=b,A\in R^{m\times n},m>n,rank(A)n,rank(A) n , r a n k ( A ) = r < n A_{m\times n},m>n,rank(A)=rn,rank(A)=r n A_{m\times n}x=0,m>n Am×nx=0,m>n而言,其最小二乘解就是 A A A的SVD分解后的 V V V的最后一个列向量。
证明:
A
=
U
[
Σ
0
]
V
T
A
x
=
U
[
Σ
0
]
V
T
x
=
0
y
=
V
T
x
,
Σ
y
=
0
也
就
是
说
,
求
A
x
=
0
,
转
换
为
求
Σ
y
=
0
的
最
小
二
乘
解
arg min
y
∣
∣
Σ
y
∣
∣
2
2
=
arg min
y
y
T
Σ
T
Σ
y
=
arg min
y
∑
i
=
1
n
σ
i
2
y
i
2
s
.
t
.
∣
∣
y
∣
∣
>
0
Σ
对
角
线
元
素
由
大
到
小
排
列
,
则
满
足
y
=
[
0
,
0
,
…
,
0
,
y
n
]
T
,
y
n
≠
0
时
,
获
得
Σ
y
的
最
小
二
乘
解
V
T
x
=
y
,
x
=
V
y
=
y
n
v
n
i
f
∣
∣
y
∣
∣
2
=
1
t
h
e
n
x
=
v
n
A=U\begin{bmatrix} \Sigma \\ 0 \end{bmatrix} V^T \\ \quad \\ Ax= U\begin{bmatrix} \Sigma \\ 0 \end{bmatrix} V^Tx=0 \\ \quad \\ y=V^Tx,\Sigma y=0 \\ \quad \\ 也就是说,求Ax=0,转换为求\Sigma y=0的最小二乘解 \\ \quad \\ \argmin_y ||\Sigma y||^2_2 \\ =\argmin_y y^T\Sigma^T\Sigma y \\ =\argmin_y \sum_{i=1}^n \sigma_i^2 y_i^2 \\ s.t. \quad ||y||>0 \\ \quad \\ \Sigma 对角线元素由大到小排列,则满足 \\ \quad \\ y=[0,0,\dots,0,y_n]^T,y_n\ne 0 \\ \quad \\ 时,获得\Sigma y的最小二乘解 \\ \quad \\ V^Tx=y,x=Vy=y_nv_n \\ \quad \\ if \quad ||y||_2=1 \\ then \quad x=v_n \\
A=U[Σ0]VTAx=U[Σ0]VTx=0y=VTx,Σy=0也就是说,求Ax=0,转换为求Σy=0的最小二乘解yargmin∣∣Σy∣∣22=yargminyTΣTΣy=yargmini=1∑nσi2yi2s.t.∣∣y∣∣>0Σ对角线元素由大到小排列,则满足y=[0,0,…,0,yn]T,yn=0时,获得Σy的最小二乘解VTx=y,x=Vy=ynvnif∣∣y∣∣2=1thenx=vn
进一步,超定齐次线性方程组
A
x
=
0
Ax=0
Ax=0的最小二乘解还等于
A
T
A
A^TA
ATA的最小特征值对应的特征向量:
证明:
A
=
U
[
Σ
0
]
V
T
A
T
A
=
V
[
Σ
0
]
[
Σ
0
]
V
T
=
V
Σ
2
V
T
A
T
A
x
=
V
Σ
2
V
T
x
=
0
y
=
V
T
x
,
Λ
=
Σ
2
,
Λ
y
=
0
A=U\begin{bmatrix} \Sigma \\ 0 \end{bmatrix} V^T \\ \quad \\ A^TA = V \begin{bmatrix} \Sigma & 0 \end{bmatrix}\begin{bmatrix} \Sigma \\ 0 \end{bmatrix}V^T=V\Sigma^2V^T \\ \quad \\ A^TAx = V\Sigma^2V^Tx=0 \\ y=V^Tx,\Lambda=\Sigma^2,\Lambda y=0 \\ \quad \\
A=U[Σ0]VTATA=V[Σ0][Σ0]VT=VΣ2VTATAx=VΣ2VTx=0y=VTx,Λ=Σ2,Λy=0
后面的证明就与上面一模一样了。
补充2:欠定行满秩方程组的线性最小二乘解法
还需要补充的是欠定方程组的最小二乘解法。虽然用的少,但却实实在在的碰到了这个问题:对极几何中本质矩阵的求解。
首先需要确定,线性齐次欠定方程组必然存在无穷组解析解。
对于非齐次欠定方程组 A m × n x = b , m < n , r a n k ( A ) = r A_{m\times n}x=b,m
