- 前言
- 一、奇异值的性质汇总
- 1. 奇异值服从的等式关系
- 2. 奇异值服从的不等式关系
- 二、秩亏缺最小二乘解
- 1. 低秩矩阵逼近
- 2. 有效秩
- 3. S V D mathrm{SVD} SVD 的子集选择方法
- 三、奇异值分解的 Q R mathrm{QR} QR 分解算法
- 四、奇异值分解的精确计算
本文学习过程来源是《矩阵分析与应用-张贤达》一书. 可以通过 z-lib 下载.
一、奇异值的性质汇总 1. 奇异值服从的等式关系-
矩阵 A m × n A_{m times n} Am×n 和其复共轭转置矩阵 A H A^{mathrm{H}} AH 具有相同的奇异值.
-
矩阵 A m × n A_{m times n} Am×n 的非零奇异值是 A A H AA^{mathrm{H}} AAH 或者 A H A A^{mathrm{H}A} AHA 的非零特征值的正平方根.
-
σ > 0 sigma > 0 σ>0 是矩阵 A m × n A_{m times n} Am×n 的单奇异值, 当且仅当 σ 2 sigma^2 σ2 是 A A H AA^{mathrm{H}} AAH 或 A H A A^{mathrm{H}}A AHA 的单特征值.
-
若 p = min { m , n } p = min{m,n} p=min{m,n}, 且 σ 1 , σ 2 , ⋯ , σ p sigma_1,sigma_2,cdots,sigma_p σ1,σ2,⋯,σp 是矩阵 A m × n A_{m times n} Am×n 的奇异值, 则
t r ( A H A ) = ∑ i = 1 p σ i 2 {mathrm{tr}}(A^{mathrm{H}}A) = sum_{i=1}^{p} sigma_i^2 tr(AHA)=i=1∑pσi2
-
矩阵行列式的绝对值等于矩阵奇异值之乘积, 即 ∣ d e t ( A ) ∣ = σ 1 σ 2 ⋯ σ n |{mathrm{det}}(A)|=sigma_1sigma_2cdotssigma_n ∣det(A)∣=σ1σ2⋯σn.
-
矩阵 A A A 的谱范数等于 A A A 的最大奇异值, 即 ∥ A ∥ s p e c = σ m a x lVert A rVert_{mathrm{spec}}=sigma_{mathrm{max}} ∥A∥spec=σmax
-
若 m ≥ n m ge n m≥n, 对于矩阵 A m × n A_{m times n} Am×n, 有
σ m i n ( A ) = min { ( x H A H A x x H x ) 1 / 2 : x ≠ 0 } = min { ( x H A H A x ) 1 / 2 : x H x = 1 , x ∈ C n } begin{aligned} sigma_{mathrm{min}}(A) &= min left { left ( frac{x^{mathrm{H}}A^{mathrm{H}}Ax}{x^{mathrm{H}}x} right )^{1/2} : x neq 0 right } \ &= min left { (x^{mathrm{H}}A^{mathrm{H}}Ax)^{1/2}:x^{mathrm{H}}x=1,x in C^n right } end{aligned} σmin(A)=min{(xHxxHAHAx)1/2:x=0}=min{(xHAHAx)1/2:xHx=1,x∈Cn}
- 若 m ≥ n m ge n m≥n, 对于矩阵 A m × n A_{m times n} Am×n, 有
σ m a x ( A ) = max { ( x H A H A x x H x ) 1 / 2 : x ≠ 0 } = max { ( x H A H A x ) 1 / 2 : x H x = 1 , x ∈ C n } begin{aligned} sigma_{mathrm{max}}(A) &= max left { left ( frac{x^{mathrm{H}}A^{mathrm{H}}Ax}{x^{mathrm{H}}x} right )^{1/2} : x neq 0 right } \ &= max left { (x^{mathrm{H}}A^{mathrm{H}}Ax)^{1/2}:x^{mathrm{H}}x=1,x in C^n right } end{aligned} σmax(A)=max{(xHxxHAHAx)1/2:x=0}=max{(xHAHAx)1/2:xHx=1,x∈Cn}
- 若 m × m m times m m×m 矩阵 A A A 非奇异, 则
1 σ min ( A ) = max { ( x H ( A − 1 ) H A − 1 x x H x ) 1 / 2 : x ≠ 0 , x ∈ C n } frac{1}{sigma_{min}(A)} = max left { left ( frac{x^{mathrm{H}}(A^{-1})^{mathrm{H}}A^{-1}x}{x^{mathrm{H}}x} right )^{1/2}:x neq 0, x in C^n right } σmin(A)1=max{(xHxxH(A−1)HA−1x)1/2:x=0,x∈Cn}
- 若 A = U [ Σ 1 O O O ] V H A = U begin{bmatrix} Sigma_1&O \O&Oend{bmatrix}V^{mathrm{H}} A=U[Σ1OOO]VH 是 m × n m times n m×n 矩阵 A A A 的奇异值分解, 则 A A A 的 M o o r e − P e n r o s e mathrm{Moore-Penrose} Moore−Penrose 逆矩阵
A † = V [ Σ 1 − 1 O O O ] U H A^{dagger} = V begin{bmatrix} Sigma_1^{-1} & O \ O & O end{bmatrix} U^{mathrm{H}} A†=V[Σ1−1OOO]UH
- 若 σ 1 , σ 2 , ⋯ , σ p sigma_1,sigma_2,cdots,sigma_p σ1,σ2,⋯,σp 是 m × n m times n m×n 矩阵 A A A 的非奇异值 (其中, p = min { m , n } p=min{m,n} p=min{m,n}), 则矩阵 [ O A A H O ] begin{bmatrix} O&A \A^{mathrm{H}}&Oend{bmatrix} [OAHAO] 具有 2 p 2p 2p 个非奇异值 σ 1 , ⋯ , σ p , − σ 1 , ⋯ , − σ p sigma_1,cdots,sigma_p,-sigma_1,cdots,-sigma_p σ1,⋯,σp,−σ1,⋯,−σp 和 ∣ m − n ∣ |m - n| ∣m−n∣ 个零奇异值.
-
若 A A A 和 B B B 是 m × n m times n m×n 矩阵, 对于 1 ≤ i , j ≤ p , i + j ≤ p + 1 ( p = min { m , n } ) 1 le i,j le p, i+j le p+1(p =min{m,n}) 1≤i,j≤p,i+j≤p+1(p=min{m,n}), 有
σ i + j − 1 ( A + B ) ≤ σ i ( A ) + σ j ( B ) sigma_{i+j-1}(A+B) le sigma_i(A) + sigma_j(B) σi+j−1(A+B)≤σi(A)+σj(B)
特别地, 当 j = 1 j=1 j=1 时, σ i ( A + B ) ≤ σ i ( A ) + σ i ( B ) , i = 1 , 2 , ⋯ , p sigma_i(A+B) le sigma_i(A) + sigma_i(B), i = 1,2,cdots,p σi(A+B)≤σi(A)+σi(B),i=1,2,⋯,p 成立.
-
对矩阵 A m × n , B m × n A_{m times n}, B_{m times n} Am×n,Bm×n, 有
σ max ( A + B ) ≤ σ max ( A ) + σ max ( B ) sigma_{max}(A+B) le sigma_{max}(A) + sigma_{max}(B) σmax(A+B)≤σmax(A)+σmax(B)
- 若 A A A 和 B B B 是 m × n m times n m×n 矩阵, 则
∑ j = 1 p [ σ j ( A + B ) − σ j ( A ) ] 2 ≤ ∥ B ∥ F 2 , p = min { m , n } sum_{j=1}^{p}[sigma_j(A+B) - sigma_j(A)]^2 le lVert B rVert_{mathrm{F}}^2, p = min{m,n} j=1∑p[σj(A+B)−σj(A)]2≤∥B∥F2, p=min{m,n}
- 若 A m × m = [ a 1 , a 2 , ⋯ , a m ] A_{m times m}=[a_1,a_2,cdots,a_m] Am×m=[a1,a2,⋯,am] 的奇异值 σ 1 ( A ) ≥ σ 2 ( A ) ≥ ⋯ ≥ σ m ( A ) sigma_1(A) ge sigma_2(A) ge cdots ge sigma_m(A) σ1(A)≥σ2(A)≥⋯≥σm(A), 则
∑ j = 1 k [ σ m − k + j ( A ) ] 2 ≤ ∑ j = 1 k a j H a j ≤ ∑ j = 1 k [ σ j ( A ) ] 2 , k = 1 , 2 , ⋯ , m sum_{j=1}^{k}[sigma_{m-k+j}(A)]^2 le sum_{j=1}^{k}a_j^{mathrm{H}}a_j le sum_{j=1}^{k}[sigma_j(A)]^2, quad k=1,2,cdots,m j=1∑k[σm−k+j(A)]2≤j=1∑kajHaj≤j=1∑k[σj(A)]2,k=1,2,⋯,m
- 若 p = min { m , n } p=min{m,n} p=min{m,n}, 且 A m × n A_{m times n} Am×n 和 B m × n B_{m times n} Bm×n 的奇异值排列为 σ 1 ( A ) ≥ σ 2 ( A ) ≥ ⋯ ≥ σ p ( A ) , σ 1 ( B ) ≥ σ 2 ( B ) ≥ ⋯ ≥ σ p ( B ) sigma_1(A) ge sigma_2(A) ge cdots ge sigma_p(A), sigma_1(B) ge sigma_2(B) ge cdots ge sigma_p(B) σ1(A)≥σ2(A)≥⋯≥σp(A),σ1(B)≥σ2(B)≥⋯≥σp(B) 和 σ 1 ( A + B ) ≥ σ 2 ( A + B ) ≥ ⋯ ≥ σ p ( A + B ) sigma_1(A+B) ge sigma_2(A+B) ge cdots ge sigma_p(A+B) σ1(A+B)≥σ2(A+B)≥⋯≥σp(A+B), 则
σ i + j − 1 ( A B H ) ≤ σ i ( A ) σ j ( B ) , 1 ≤ i , j ≤ p , i + j ≤ p + 1 sigma_{i+j-1}(AB^{mathrm{H}}) le sigma_i(A)sigma_j(B), 1 le i,j le p,i+j le p+1 σi+j−1(ABH)≤σi(A)σj(B),1≤i,j≤p,i+j≤p+1
-
设 m × ( n − 1 ) m times (n-1) m×(n−1) 矩阵 B B B 是删去 m × n m times n m×n 矩阵 A A A 任意一列得到的矩阵, 并且它们的奇异值都按照非降顺序排列, 则
σ 1 ( A ) ≥ σ 1 ( B ) ≥ σ 2 ( A ) ≥ σ 2 ( B ) ≥ ⋯ ≥ σ h ( A ) ≥ σ h ( B ) ≥ 0 sigma_1(A) ge sigma_1(B) ge sigma_2(A) ge sigma_2(B) ge cdots ge sigma_h(A) ge sigma_h(B) ge 0 σ1(A)≥σ1(B)≥σ2(A)≥σ2(B)≥⋯≥σh(A)≥σh(B)≥0
式中, h = min { m , n − 1 } h = min{m, n-1} h=min{m,n−1}
-
设 ( m − 1 ) × n (m-1) times n (m−1)×n 矩阵 B B B 是删去 m × n m times n m×n 矩阵 A A A 任意一行得到的矩阵, 并且它们的奇异值都按照非降顺序排列, 则
σ 1 ( A ) ≥ σ 1 ( B ) ≥ σ 2 ( A ) ≥ σ 2 ( B ) ≥ ⋯ ≥ σ h ( A ) ≥ σ h ( B ) ≥ 0 sigma_1(A) ge sigma_1(B) ge sigma_2(A) ge sigma_2(B) ge cdots ge sigma_h(A) ge sigma_h(B) ge 0 σ1(A)≥σ1(B)≥σ2(A)≥σ2(B)≥⋯≥σh(A)≥σh(B)≥0
式中, h = min { m , n − 1 } h = min{m, n-1} h=min{m,n−1}
-
矩阵 A m × n A_{m times n} Am×n 的最大奇异值满足不等式
σ max ( A ) ≥ [ 1 n t r ( A H A ) ] 1 / 2 sigma_{max}(A) ge [frac{1}{n}mathrm{tr}(A^{mathrm{H}}A)]^{1/2} σmax(A)≥[n1tr(AHA)]1/2
二、秩亏缺最小二乘解 1. 低秩矩阵逼近在奇异值分析的应用中, 常常需要用一个低秩的矩阵逼近一个含噪声或扰动的矩阵. 低秩矩阵好理解, 但是这和现实世界中又有什么关系呢?
举个非常简单的例子, 当拍摄一张大海的照片, 看起来就好像一个像素经过不断复杂形成的. 想象最小复制单元是一行或一列, 那么类推到矩阵中是不是就是说明组成图像的所有向量其实是线性相关的, 然后矩阵的秩就为 1, 这就叫做低秩矩阵. 但是当图中出现轮船, 海鸟, 岛屿等事物, 图片的复杂程度就会上升, 也就不能够叫做低秩矩阵.
在信号处理中, 用低秩矩阵逼近信号矩阵, 就是为了去除噪声或扰动的信号.
下面定理给出了逼近质量的评价.
定理 1 : 令 A ∈ R m × n A in R^{m times n} A∈Rm×n 的奇异值分解由 A = ∑ i = 1 p σ i u i v i T A = sum_{i=1}^{p}sigma_iu_iv_i^{mathrm{T}} A=∑i=1pσiuiviT 给出, 其中, p = r a n k ( A ) p=mathrm{rank}(A) p=rank(A)
若
k
<
p
k < p
k
min
r
a
n
k
(
B
)
=
k
∥
A
−
B
∥
s
p
e
c
=
∥
A
−
A
k
∥
s
p
e
c
=
σ
k
+
1
(1)
min_{mathrm{rank}(B)=k}lVert A-B rVert_{mathrm{spec}} = lVert A-A_k rVert_{mathrm{spec}} = sigma_{k+1} tag{1}
rank(B)=kmin∥A−B∥spec=∥A−Ak∥spec=σk+1(1)
min
r
a
n
k
(
B
)
=
k
∥
A
−
B
∥
F
=
∥
A
−
A
k
∥
F
=
∑
i
=
k
+
1
q
σ
i
2
(2)
min_{mathrm{rank}(B)=k}lVert A-B rVert_{mathrm{F}} = lVert A-A_k rVert_{mathrm{F}} = sqrt{sum_{i=k+1}^{q}sigma_i^2} tag{2}
rank(B)=kmin∥A−B∥F=∥A−Ak∥F=i=k+1∑qσi2
(2) 式子中,
q
=
min
{
m
,
n
}
q=min{m,n}
q=min{m,n}. 在需要计算秩
r
r
r 的估计值
r
^
hat{r}
r^ 的方法, 在信号处理和系统理论中, 常将该估计值称为 “有效秩” 有效秩确定有以下两种常用方法 归一化奇异值方法 计算归一化奇异值
σ
i
ˉ
=
σ
i
^
σ
1
^
bar{sigma_i} = frac{hat{sigma_i}}{hat{sigma_1}}
σiˉ=σ1^σi^ 选择满足准则
σ
i
ˉ
≥
ϵ
bar{sigma_i} ge epsilon
σiˉ≥ϵ 的最大整数作为有效秩的估计值
r
^
hat{r}
r^. 范数比方法 令
m
×
n
m times n
m×n 矩阵
A
k
A_k
Ak 是原
m
×
n
m times n
m×n 矩阵
A
A
A 的秩
k
k
k 近似, 定义该近似矩阵与原矩阵的
F
r
o
b
e
n
i
u
s
mathrm{Frobenius}
Frobenius 范数比为
v
(
k
)
=
∥
A
k
∥
F
∥
A
∥
F
=
σ
1
2
+
σ
2
2
+
⋯
+
σ
k
2
σ
1
2
+
σ
2
2
+
⋯
+
σ
h
2
,
h
=
min
{
m
,
n
}
v(k) = frac{lVert A_k rVert_{mathrm{F}}}{lVert A rVert_{mathrm{F}}} = frac{sqrt{sigma_1^2 + sigma_2^2 + cdots + sigma_k^2}}{sigma_1^2 + sigma_2^2 + cdots + sigma_h^2}, quad h = min{m,n}
v(k)=∥A∥F∥Ak∥F=σ12+σ22+⋯+σh2σ12+σ22+⋯+σk2
,h=min{m,n} 并选择满足
v
(
k
)
≥
α
v(k) ge alpha
v(k)≥α 的最大整数作为有效秩估计值
r
^
hat{r}
r^. 其中
α
alpha
α 是接近于 1 的阈值. 算法 1 : 子集选择算法 Step 1 : 计算
A
A
A 的
S
V
D
mathrm{SVD}
SVD, 并确定
A
A
A 的有效秩
r
^
hat{r}
r^ Step 2 : 计算置换矩阵
P
P
P, 使得在
A
P
=
[
B
1
,
B
2
]
AP=[B_1,B_2]
AP=[B1,B2] 中的矩阵
B
1
∈
C
m
×
r
^
B_1 in C^{m times hat{r}}
B1∈Cm×r^ 的列是 “足够线性无关的” Step 3 : 将
A
x
=
0
Ax=0
Ax=0 的
L
S
mathrm{LS}
LS 问题变换为求
A
P
[
z
0
]
=
b
APbegin{bmatrix}z\0end{bmatrix}=b
AP[z0]=b 的
L
S
mathrm{LS}
LS 解
z
∈
C
r
^
z in C^{hat{r}}
z∈Cr^, 即求
∥
B
1
z
=
b
∥
2
lVert B_1z = b rVert_2
∥B1z=b∥2 的极小化变量
z
z
z. 算法 2 : 低秩
L
S
mathrm{LS}
LS 算法, 给定
A
∈
R
m
×
n
,
b
∈
R
m
A in R^{m times n}, b in R^m
A∈Rm×n,b∈Rm, 计算
x
=
a
r
g
min
z
∥
A
P
[
z
0
]
−
b
∥
2
x = mathrm{arg}min_z Bigg lVert APbegin{bmatrix}z\0end{bmatrix}-b Bigg rVert_2
x=argzmin∥∥∥∥∥AP[z0]−b∥∥∥∥∥2 Step 1 : 计算
S
V
D
mathrm{SVD}
SVD, 确定有效秩
r
^
hat{r}
r^, 根据
r
^
hat{r}
r^ 将
V
V
V 分块为
V
=
[
V
11
V
12
V
21
V
22
]
V = begin{bmatrix} V_{11}& V_{12}\ V_{21}& V_{22} end{bmatrix}
V=[V11V21V12V22] 其中,
V
11
∈
R
r
^
×
r
^
V_{11} in R^{hat{r} times hat{r}}
V11∈Rr^×r^, 然后存储
V
11
V_{11}
V11 和
V
21
V_{21}
V21 Step 2 : 利用列主元
Q
R
mathrm{QR}
QR 算法计算
Q
T
[
V
11
T
,
V
21
T
]
P
=
[
R
11
,
R
12
]
Q^mathrm{T}[V_{11}^mathrm{T},V_{21}^mathrm{T}]P=[R_{11},R_{12}]
QT[V11T,V21T]P=[R11,R12], 然后计算
A
P
=
[
B
1
,
B
2
]
AP=[B_1,B_2]
AP=[B1,B2], 其中
B
1
∈
R
m
×
r
^
B_1 in R^{m times hat{r}}
B1∈Rm×r^ Step 3 : 计算
z
=
(
B
1
T
B
1
)
−
1
B
1
T
b
z = (B_1^mathrm{T}B_1)^{-1}B_1^mathrm{T}b
z=(B1TB1)−1B1Tb 附 列主元
Q
R
mathrm{QR}
QR 算法 给定矩阵
A
∈
R
m
×
n
A in R^{m times n}
A∈Rm×n, 其中,
m
≥
n
m ge n
m≥n. 下面的算法计算
r
=
r
a
n
k
(
A
)
r = mathrm{rank}(A)
r=rank(A) 和列主元
Q
R
mathrm{QR}
QR 分解
Q
T
A
P
=
[
R
11
R
12
O
O
]
Q^mathrm{T}AP=begin{bmatrix} R_{11}&R_{12} \ O&O end{bmatrix}
QTAP=[R11OR12O] 其中,
R
11
∈
R
r
×
r
R_{11} in R^{r times r}
R11∈Rr×r 是上三角的非奇异矩阵, 作为输出结果,
A
A
A 的上三角部分存放
R
R
R 的上三角部分, 置换矩阵
P
P
P 用整数向量
p
i
v
piv
piv 编码 (若
p
i
v
piv
piv 的第
j
j
j 个元素等于整数
m
m
m, 则
P
P
P 的第
j
j
j 列仅第
m
m
m 个元素为 1, 而其他元素皆为零) 例: 假定
[
3
4
−
1
7
4
−
3
2
5
3
−
1
4
5
]
,
b
=
[
1
1
1
1
]
begin{bmatrix} 3& 4& -1\ 7& 4& -3\ 2& 5& 3\ -1& 4& 5 end{bmatrix}, quad b = begin{bmatrix} 1\ 1\ 1\ 1 end{bmatrix}
⎣⎢⎢⎡372−14454−1−335⎦⎥⎥⎤,b=⎣⎢⎢⎡1111⎦⎥⎥⎤ 得到
r
a
n
k
(
A
)
=
2
mathrm{rank}(A)=2
rank(A)=2 和
x
L
S
=
[
0.0815
0.1545
0.0730
]
x_{mathrm{LS}} = begin{bmatrix} 0.0815\ 0.1545\ 0.0730 end{bmatrix}
xLS=⎣⎡0.08150.15450.0730⎦⎤ 应用算法 低秩
L
S
mathrm{LS}
LS 方法得到
P
=
[
0
1
0
1
0
0
0
0
1
]
,
x
=
[
0.0845
0.2275
0.0000
]
P = begin{bmatrix} 0& 1& 0\ 1& 0& 0\ 0& 0& 1 end{bmatrix}, quad x = begin{bmatrix} 0.0845\ 0.2275\ 0.0000 end{bmatrix}
P=⎣⎡010100001⎦⎤,x=⎣⎡0.08450.22750.0000⎦⎤ 通常是两个阶段: 第一阶段为矩阵的二重对角化, 通过
H
o
u
s
e
h
o
l
d
e
r
mathrm{Householder}
Householder 变换将矩阵
A
m
×
n
A_{m times n}
Am×n 变换为二重对角矩阵 (除对角线及其上面一条对角线的元素外, 其他元素全为零). 第二阶段利用
Q
R
mathrm{QR}
QR 分解, 保存二重对角矩阵的形式不变, 利用正交变换使上一条对角线的元素逐渐减小, 使矩阵接近对角矩阵. 此时, 需要具体构造正交矩阵
V
V
V, 但
U
U
U 可以不予保留. 和之前在线性代数中学习的人工进行
S
V
D
SVD
SVD 分解不同, 这里的算法更适合于计算机. 有点难度, 需要慢慢消化. (… 未完待续) 初始矩阵
A
A
A 的元素通过观察或计算得到, 存在一定的误差. 为了使得奇异值和奇异向量的计算结果是精确的, 提出了右边
J
a
c
o
b
i
mathrm{Jacobi}
Jacobi 旋转. (… 未完待续)
3.
S
V
D
mathrm{SVD}
SVD 的子集选择方法
