奇异值分解(Singular Value Decomposition,SVD) 是线性代数中重要的矩阵分解,是特征分解在任意矩阵上的推广,在立体视觉、三维重建领域应用非常广泛。由于可以用于求解线性方程的最小二乘解,所以在求解本质矩阵、单应性矩阵、点云刚性变换矩阵时,都能用到SVD。本篇即给大家简单介绍下奇异值分解,并通过公式推导来说明其在线性最小二乘解
A
x
=
b
Ax=b
Ax=b上的应用。
奇异值分解(SVD Decomposition)定义
对于任意矩阵
A
∈
R
m
×
n
A\in R^{m\times n}
A∈Rm×n,存在正交矩阵
U
∈
R
m
×
m
U\in R^{m\times m}
U∈Rm×m,
V
∈
R
n
×
n
V\in R^{n\times n}
V∈Rn×n,以及对角矩阵
Σ
∈
R
m
×
n
\Sigma \in R^{m\times n}
Σ∈Rm×n:
其中,对角线元素满足:
σ
1
≥
⋯
≥
σ
r
≥
σ
r
+
1
=
⋯
=
σ
min
(
m
,
n
)
=
0
\sigma_1\geq\cdots\geq\sigma_r\geq\sigma_{r+1}=\cdots=\sigma_{\min(m,n)}=0
σ1≥⋯≥σr≥σr+1=⋯=σmin(m,n)=0
使得,
A
=
U
Σ
V
T
A=U\Sigma V^T
A=UΣVT。
此分解叫做矩阵 A A A的奇异值分解(Singular Value Decomposition),是一个非常重要的矩阵分解,对角矩阵 Σ \Sigma Σ的对角线元素 σ i \sigma_i σi叫做矩阵 A A A的奇异值。矩阵 U U U的列向量成为左奇异向量,矩阵 V V V的列向量成为右奇异向量。
使用正交矩阵
V
V
V,可得到以下公式:
A
V
=
U
Σ
AV=U\Sigma
AV=UΣ
这可解释为,存在一组特殊的正交向量集(例如 V V V的列向量集),通过矩阵 A A A映射到另一组正交向量集(例如 U U U的列向量集)。
SVD的一些特性
给定矩阵
A
A
A的一组SVD分解
A
=
U
Σ
V
T
A=U\Sigma V^T
A=UΣVT
其中, σ 1 ≥ ⋯ ≥ σ r ≥ σ r + 1 = ⋯ = σ min ( m , n ) = 0 \sigma_1\geq\cdots\geq\sigma_r\geq\sigma_{r+1}=\cdots=\sigma_{\min(m,n)}=0 σ1≥⋯≥σr≥σr+1=⋯=σmin(m,n)=0
存在以下推论( R ( A ) R(A) R(A)和 N ( A ) N(A) N(A)分别为矩阵 A A A的值域空间和零空间):
- r a n k ( A ) = r rank(A)=r rank(A)=r
- R ( A ) = R ( [ u 1 , . . . u r ] ) R(A)=R([u_1,...u_r]) R(A)=R([u1,...ur])
- N ( A ) = R ( [ u r + 1 , . . . , u n ] ) N(A)=R([u_{r+1},...,u_n]) N(A)=R([ur+1,...,un])
- R ( A T ) = R ( [ v 1 , . . . v r ] ) R(A^T)=R([v_1,...v_r]) R(AT)=R([v1,...vr])
- N ( A T ) = R ( [ v r + 1 , . . . , v m ] ) N(A^T)=R([v_{r+1},...,v_m]) N(AT)=R([vr+1,...,vm])
如果我们引入
U
r
=
[
u
1
,
.
.
.
,
u
r
]
,
Σ
=
d
i
a
g
(
σ
1
,
.
.
.
,
σ
r
)
,
V
r
=
[
v
1
,
.
.
.
,
v
r
]
U_r=[u_1,...,u_r],\Sigma=diag(\sigma_1,...,\sigma_r),V_r=[v_1,...,v_r]
Ur=[u1,...,ur],Σ=diag(σ1,...,σr),Vr=[v1,...,vr]
则有
A
=
U
r
Σ
r
V
r
T
=
∑
i
=
1
r
σ
i
u
i
v
i
T
A=U_r\Sigma_rV_r^T=\sum_{i=1}^r\sigma_iu_iv_i^T
A=UrΣrVrT=i=1∑rσiuiviT
这称为矩阵
A
A
A的二进制分解(dyadic decomposition),即将秩为
r
r
r的矩阵
A
A
A分解为
r
r
r个秩为1的矩阵之和。
同时,可以得到
A
T
A
=
V
Σ
T
Σ
V
T
a
n
d
A
A
T
=
U
Σ
Σ
T
U
T
A^TA=V\Sigma^T\Sigma V^T and AA^T=U\Sigma\Sigma^TU^T
ATA=VΣTΣVTandAAT=UΣΣTUT
可知,奇异值的平方
σ
i
2
,
i
=
1
,
.
.
.
,
p
\sigma_i^2,i=1,...,p
σi2,i=1,...,p是对称矩阵
A
T
A
A^TA
ATA和
A
A
T
AA^T
AAT的特征值,
v
i
v_i
vi和
u
i
u_i
ui分别是对应的特征向量。 证明如下图:
该结论很有用,我们在解
A
x
=
0
Ax=0
Ax=0时,需要求
A
T
A
A^TA
ATA的最小特征值对应的特征向量,即为SVD分解后矩阵
V
V
V的最后一列。另一种应用场景是协方差矩阵的PCA分析时,可直接以SVD分解后的矩阵
V
V
V各列为各方向矢量。
通过该结论,我们很容易的可以计算矩阵
A
A
A的2-范数和Frobenius范数:
∣
∣
A
∣
∣
2
=
max
x
≠
0
∣
∣
A
x
∣
∣
2
∣
∣
x
∣
∣
2
=
max
x
≠
0
x
T
A
T
A
x
x
T
x
=
max
x
≠
0
x
T
λ
A
T
A
x
x
T
x
=
max
x
≠
0
λ
A
T
A
=
σ
1
∣
∣
A
∣
∣
F
=
∑
i
=
1
m
∑
j
=
1
n
a
i
j
2
=
t
r
a
c
e
(
A
T
A
)
=
σ
1
2
+
⋯
+
σ
p
2
,
p
=
min
(
m
,
n
)
\begin{aligned} ||A||_2&=\sqrt{\max_{x\neq0}\frac{||Ax||_2}{||x||_2}}=\sqrt{\max_{x\neq0}\frac{x^TA^TAx}{x^Tx}}=\sqrt{\max_{x\neq0}\frac{x^T\lambda_{A^TA}x}{x^Tx}}=\sqrt{\max_{x\neq0}\lambda_{A^TA}}=\sigma_1\\ ||A||_F&=\sqrt{\displaystyle\sum_{i=1}^m\sum_{j=1}^na_{ij}^2}=\sqrt{trace(A^TA)}=\sqrt{\sigma_1^2+\cdots+\sigma_p^2},p=\min(m,n) \end{aligned}
∣∣A∣∣2∣∣A∣∣F=x=0max∣∣x∣∣2∣∣Ax∣∣2=x=0maxxTxxTATAx=x=0maxxTxxTλATAx=x=0maxλATA=σ1=i=1∑mj=1∑naij2=trace(ATA)=σ12+⋯+σp2,p=min(m,n)
SVD求解线性最小二乘问题Ax=b
我们再来看SVD的经典应用:解线性方程的最小二乘解。
考虑一个线性最小二乘问题
A
x
=
b
Ax=b
Ax=b:
min
x
∣
∣
A
x
−
b
∣
∣
2
2
\min_x||Ax-b||_2^2
xmin∣∣Ax−b∣∣22
给定矩阵 A ∈ R m × n A\in R^{m\times n} A∈Rm×n的一组SVD分解: A = U Σ V T A=U\Sigma V^T A=UΣVT
根据矩阵
U
U
U和
V
V
V的正交性,有
∣
∣
A
x
−
b
∣
∣
2
2
=
∣
∣
U
T
(
A
x
−
b
)
∣
∣
2
2
=
∣
∣
U
T
A
x
−
U
T
b
∣
∣
2
2
=
∣
∣
Σ
V
T
x
−
U
T
b
∣
∣
2
2
||Ax-b||_2^2=||U^T(Ax-b)||_2^2=||U^TAx-U^Tb||_2^2=||\Sigma V^Tx-U^Tb||_2^2
∣∣Ax−b∣∣22=∣∣UT(Ax−b)∣∣22=∣∣UTAx−UTb∣∣22=∣∣ΣVTx−UTb∣∣22
另
z
=
V
T
x
z=V^Tx
z=VTx,则有
∣
∣
A
x
−
b
∣
∣
2
2
=
∣
∣
Σ
V
T
x
−
U
T
b
∣
∣
2
2
=
∑
i
=
1
r
(
σ
i
z
i
−
u
i
T
b
)
2
+
∑
i
=
r
+
1
m
(
u
i
T
b
)
2
||Ax-b||_2^2=||\Sigma V^Tx-U^Tb||_2^2=\displaystyle\sum_{i=1}^r(\sigma_iz_i-u_i^Tb)^2+\sum_{i=r+1}^m(u_i^Tb)^2
∣∣Ax−b∣∣22=∣∣ΣVTx−UTb∣∣22=i=1∑r(σizi−uiTb)2+i=r+1∑m(uiTb)2
因此,
min
x
∣
∣
A
x
−
b
∣
∣
2
2
=
min
x
(
∑
i
=
1
r
(
σ
i
z
i
−
u
i
T
b
)
2
+
∑
i
=
r
+
1
m
(
u
i
T
b
)
2
)
\min_x||Ax-b||_2^2=\min_x(\displaystyle\sum_{i=1}^r(\sigma_iz_i-u_i^Tb)^2+\sum_{i=r+1}^m(u_i^Tb)^2)
xmin∣∣Ax−b∣∣22=xmin(i=1∑r(σizi−uiTb)2+i=r+1∑m(uiTb)2)
显然,当
σ
i
z
i
=
u
i
T
b
\sigma_iz_i=u_i^Tb
σizi=uiTb时,取最小,则最小二乘解为:
z
i
=
u
i
T
b
σ
i
,
i
=
1
,
.
.
.
,
r
z
i
=
a
r
b
i
t
r
a
r
y
,
i
=
r
+
1
,
.
.
.
,
n
\begin{aligned} z_i&=\frac{u_i^Tb}{\sigma_i},i=1,...,r\\ z_i&=arbitrary,i=r+1,...,n \end{aligned}
zizi=σiuiTb,i=1,...,r=arbitrary,i=r+1,...,n
最小值为
min
x
∣
∣
A
x
−
b
∣
∣
2
2
=
∑
i
=
r
+
1
m
(
u
i
T
b
)
2
\min_x||Ax-b||_2^2=\sum_{i=r+1}^m(u_i^Tb)^2
xmin∣∣Ax−b∣∣22=i=r+1∑m(uiTb)2
由
z
=
V
T
x
z=V^Tx
z=VTx,可得
x
=
V
z
x=Vz
x=Vz,因此,通过SVD来求线性最小二乘解的公式为
x
=
V
z
z
i
=
u
i
T
b
σ
i
,
i
=
1
,
.
.
.
,
r
z
i
=
a
r
b
i
t
r
a
r
y
,
i
=
r
+
1
,
.
.
.
,
n
\begin{aligned} x&=Vz\\ z_i&=\frac{u_i^Tb}{\sigma_i},i=1,...,r\\ z_i&=arbitrary,i=r+1,...,n \end{aligned}
xzizi=Vz=σiuiTb,i=1,...,r=arbitrary,i=r+1,...,n
在上式中我们发现,当
r
=
n
r=n
r=n时,有唯一最小二乘解,当
r
<
n
r<n
r<n时,有无数最小二乘解,此时我们可求最小范数解:
x
†
=
V
z
†
z
i
†
=
u
i
T
b
σ
i
,
i
=
1
,
.
.
.
,
r
z
i
†
=
0
,
i
=
r
+
1
,
.
.
.
,
n
\begin{aligned} x_\dagger&=Vz_\dagger\\ z_i^\dagger&=\frac{u_i^Tb}{\sigma_i},i=1,...,r\\ z_i^\dagger&=0,i=r+1,...,n \end{aligned}
x†zi†zi†=Vz†=σiuiTb,i=1,...,r=0,i=r+1,...,n
我们知道最小二乘解也可通过 ( A T A ) − 1 ( A T b ) (A^TA)^{-1}(A^Tb) (ATA)−1(ATb)来求解,而通过SVD求解的优点是不需要复杂的求逆运算,且可处理 A T A A^TA ATA为奇异阵不可逆的情况。
博主简介:
Ethan Li 李迎松(知乎:李迎松)
武汉大学 摄影测量与遥感专业博士
主方向立体匹配、三维重建
2019年获测绘科技进步一等奖(省部级)
爱三维,爱分享,爱开源
GitHub: https://github.com/ethan-li-coding
邮箱:ethan.li.whu@gmail.com
个人微信:
欢迎交流!
关注博主不迷路,感谢!
博客主页:https://ethanli.blog.csdn.net