本篇文章参考了李航老师的《统计学习方法》第二版。,同时参考了机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用
矩阵的奇异值分解是将一个矩阵分解成多个小矩阵的乘积,在PCA(principal component analysis,主成分分析)和LSA(latent semantic analysis,潜在语义分析)中都有应用,在此把矩阵的分解过程记录下来,证明过程先略过。
矩阵的SVD分解
首先看一下定义,给定非0的m
×
\times
×n实矩阵A,它可以表示成以下三个实矩阵乘积形式的运算。
A
=
U
D
V
T
A = UDV^T
A=UDVT
其中
U
U
U是
m
m
m阶正交矩阵,
V
V
V是
n
n
n阶正交矩阵。
D
D
D是由降序排列的非负的对角线元素组成的
m
×
n
m\times n
m×n对角矩阵。满足如下关系
U
U
T
=
I
UU^T=I
UUT=I
V V T = I VV^T=I VVT=I
D = d i a g ( σ 1 , σ 2 , . . . , σ p ) D = diag(\sigma_1, \sigma_2,...,\sigma_p) D=diag(σ1,σ2,...,σp)
σ 1 ≥ σ 2 ≥ . . . ≥ σ p ≥ 0 \sigma_1 \geq \sigma_2 \geq...\geq \sigma_p \geq0 σ1≥σ2≥...≥σp≥0
p = m i n ( m , n ) p = min(m,n) p=min(m,n)
计算过程
(1)首先求出
A
A
T
AA^T
AAT的特征值和特征向量
计算对称矩阵
W
=
A
T
A
W=A^TA
W=ATA。
求解特征方程
(
W
−
λ
I
)
x
=
0
(W-\lambda I)x = 0
(W−λI)x=0,得到特征值
λ
i
\lambda_i
λi并将特征值由大到小排列得到
λ
1
≥
λ
2
≥
.
.
.
≥
λ
n
≥
0
\lambda_1 \geq \lambda_2 \geq ... \geq \lambda_n \geq 0
λ1≥λ2≥...≥λn≥0
将特征值代入特征方程求得对应的特征向量。
(2)求
n
n
n阶正交矩阵
V
V
V
将特征向量单位化,得到单位特征向量
v
1
,
v
2
,
.
.
.
,
v
n
v_1, v_2,...,v_n
v1,v2,...,vn,构成
n
n
n阶正交矩阵
V
V
V
V
=
[
v
1
v
2
.
.
.
v
n
]
V = [v_1 \ v_2 \ ... \ v_n]
V=[v1 v2 ... vn]
(3)求
m
×
n
m\times n
m×n对角矩阵
D
D
D
计算
A
A
A的奇异值
σ
i
=
λ
i
,
i
=
1
,
2
,
.
.
.
,
n
\sigma_i = \sqrt{\lambda_i}, i=1,2,...,n
σi=λi,i=1,2,...,n
构造
m
×
n
m \times n
m×n矩形对角矩阵
D
D
D,主对角线元素是奇异值,其余元素是0。
D
=
d
i
a
g
(
σ
1
,
σ
2
,
.
.
.
,
σ
n
)
D = diag(\sigma_1, \sigma_2, ..., \sigma_n)
D=diag(σ1,σ2,...,σn)
(4)求
m
m
m阶正交矩阵
U
U
U
对
A
A
A的前
r
r
r个正奇异值,令
u
j
=
1
σ
j
A
v
j
u_j = \frac{1}{\sigma_j}Av_j
uj=σj1Avj
得到
U 1 = [ u 1 u 2 . . . u r ] U_1 = [u_1 \ u_2 \ ... \ u_r] U1=[u1 u2 ... ur]
求 A T A^T AT的零空间的一组标准正交基 { u r + 1 , u r + 2 , . . . , u m } \{u_{r+1}, u_{r+2}, ... ,u_m\} {ur+1,ur+2,...,um},令
U 2 = [ u r + 1 u r + 2 . . . u m ] U_2 = [u_{r+1} \ u_{r+2} \ ... \ u_m] U2=[ur+1 ur+2 ... um]
并令
U
=
[
U
1
U
2
]
U = [U_1 \ U_2]
U=[U1 U2]
(5)得到
A
A
A的奇异值分解
A = U D V T A = UDV^T A=UDVT
举例计算
上面的基本步骤,下面来看一个具体的例子。
求矩阵
A
A
A的奇异值分解。
A
=
[
1
1
2
2
0
0
]
A=\begin{bmatrix} 1 &1\\ 2 &2\\ 0 &0 \end{bmatrix}
A=⎣⎡120120⎦⎤
(1)求矩阵
A
T
A
A^TA
ATA的特征值和特征向量
A T A = [ 1 2 0 1 2 0 ] [ 1 1 2 2 0 0 ] = [ 5 5 5 5 ] A^TA=\begin{bmatrix} 1 & 2 & 0\\ 1 & 2 & 0\\ \end{bmatrix} \begin{bmatrix} 1 &1\\ 2 &2\\ 0 &0 \end{bmatrix} = \begin{bmatrix} 5 & 5\\ 5 & 5\\ \end{bmatrix} ATA=[112200]⎣⎡120120⎦⎤=[5555]
特征值 λ \lambda λ和特征向量 x x x满足特征方程
( A T A − λ I ) x = 0 (A^TA-\lambda I)x = 0 (ATA−λI)x=0
得到齐次线性方程组
{ ( 5 − λ ) x 1 + 5 x 2 = 0 5 x 1 + ( 5 − λ ) x 2 = 0 \left \{\begin{array}{lr} (5-\lambda)x_1+5x_2 = 0 \\ 5x_1+(5-\lambda)x_2 = 0 \end{array} \right. {(5−λ)x1+5x2=05x1+(5−λ)x2=0
该方程组有非零解的充要条件是
∣
5
−
λ
5
5
5
−
λ
∣
=
0
\left| \begin{array}{cccc} 5-\lambda & 5 \\ 5 & 5-\lambda \end{array} \right| = 0
∣∣∣∣5−λ555−λ∣∣∣∣=0
即
λ
2
−
10
λ
=
0
\lambda^2-10\lambda=0
λ2−10λ=0,解此方程,会得到矩阵
A
T
A
A^TA
ATA的特征值
λ
1
=
10
\lambda_1=10
λ1=10和
λ
2
=
0
\lambda_2=0
λ2=0。
将特征值
λ
1
=
10
\lambda_1=10
λ1=10代入线性方程组,得到对应的单位特征向量
v 1 = [ 1 2 1 2 ] v_1 = \begin{bmatrix} \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} v1=[2121]
同样得到特征值 λ 2 = 0 \lambda_2=0 λ2=0对应的单位特征向量
v 2 = [ 1 2 − 1 2 ] v_2 = \begin{bmatrix} \frac{1}{\sqrt{2}} \\ -\frac{1}{\sqrt{2}} \end{bmatrix} v2=[21−21]
(2)求正交矩阵 V V V
V = [ v 1 v 2 ] = [ 1 2 1 2 1 2 − 1 2 ] V = [v_1 \ v_2] = \begin{bmatrix} \frac{1}{\sqrt 2} & \frac{1}{\sqrt 2} \\ \frac{1}{\sqrt 2} & -\frac{1}{\sqrt 2}\end{bmatrix} V=[v1 v2]=[212121−21]
(3)求对角矩阵
奇异值为
σ
1
=
λ
1
=
10
\sigma_1 = \sqrt{\lambda_1} = \sqrt{10}
σ1=λ1=10和
σ
2
=
0
\sigma_2=0
σ2=0,构造对角矩阵
D
=
[
10
0
0
0
0
0
]
D = \begin{bmatrix} \sqrt{10} & 0 \\ 0 & 0 \\ 0 & 0 \\ \end{bmatrix}
D=⎣⎡1000000⎦⎤
注意要到
D
D
D中加上行向量0,使得三个矩阵可以进行相乘。
(4)求正交矩阵
U
U
U
基于
A
A
A的正奇异值
σ
1
\sigma_1
σ1计算得到列向量
u
1
u_1
u1
u
1
=
1
σ
1
A
v
1
=
1
10
[
1
1
2
2
0
0
]
[
1
2
1
2
]
=
[
1
5
2
5
0
]
u_1 = \frac{1}{\sigma_1}Av_1 = \frac{1}{\sqrt{10}} \begin{bmatrix} 1 & 1\\ 2 & 2\\ 0 & 0 \\ \end{bmatrix} \begin{bmatrix} \frac{1}{\sqrt2} \\ \frac{1}{\sqrt2}\end{bmatrix} = \begin{bmatrix} \frac{1}{\sqrt5} \\ \frac{2}{\sqrt{5}} \\ 0\end{bmatrix}
u1=σ11Av1=101⎣⎡120120⎦⎤[2121]=⎣⎡51520⎦⎤
列向量 u 2 , u 3 u_2, u_3 u2,u3是 A T A^T AT的零空间的一组标准正交基。为此,求解以下线性方程组
A T x = [ 1 2 0 1 2 0 ] [ x 1 x 2 x 3 ] = [ 0 0 ] A^Tx = \begin{bmatrix} 1 & 2 & 0 \\ 1 & 2 & 0\\ \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} ATx=[112200]⎣⎡x1x2x3⎦⎤=[00]
即
x
1
+
2
x
2
+
0
x
3
=
0
x_1+2x_2+0x_3=0
x1+2x2+0x3=0,分别取
(
x
2
,
x
3
)
(x_2, x_3)
(x2,x3)为(1,0)和(0,1),得到一组正交基
(
−
2
,
1
,
0
)
T
,
(
0
,
0
,
1
)
T
(-2, 1, 0)^T, (0,0,1)^T
(−2,1,0)T,(0,0,1)T。化成标准正交基之后
u
2
=
(
−
2
5
,
1
5
,
0
)
,
u
3
=
(
0
,
0
,
1
)
T
u_2 = (-\frac{2}{\sqrt5},\frac{1}{\sqrt5}, 0), u_3 = (0,0, 1)^T
u2=(−52,51,0),u3=(0,0,1)T
构造正交矩阵 U U U
U
=
[
1
5
−
2
5
0
2
5
1
5
0
0
0
1
]
U = \begin{bmatrix} \frac{1}{\sqrt5} & -\frac{2}{\sqrt5} & 0 \\ \frac{2}{\sqrt5} & \frac{1}{\sqrt5} & 0 \\ 0 & 0 & 1 \end{bmatrix}
U=⎣⎡51520−52510001⎦⎤
(5)矩阵
A
A
A的奇异值分解
A
=
U
D
V
T
=
[
1
5
−
2
5
0
2
5
1
5
0
0
0
1
]
[
10
0
0
0
0
0
]
[
1
2
1
2
1
2
−
1
2
]
A = UDV^T = \begin{bmatrix} \frac{1}{\sqrt5} & -\frac{2}{\sqrt5} & 0 \\ \frac{2}{\sqrt5} & \frac{1}{\sqrt5} & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \sqrt{10} & 0 \\ 0 & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} \frac{1}{\sqrt2} & \frac{1}{\sqrt2} \\ \frac{1}{\sqrt2} & - \frac{1}{\sqrt2}\end{bmatrix}
A=UDVT=⎣⎡51520−52510001⎦⎤⎣⎡1000000⎦⎤[212121−21]
再来看一段简单的程序,是不是和得到的结果一致。
import numpy as np
A = [[1,1],[2,2],[0,0]]
u,d,v = np.linalg.svd(A)
结果可以看作是一致的,因为奇异值分解的值并不唯一。
SVD的意义
根据机器学习中的数学(5)-强大的矩阵奇异值分解(SVD)及其应用这篇文章的观点,在很多情况下,前10%甚至前1%的奇异值的和就占了全部奇异值之和的99%以上,也就是说,我们可以用前r大的奇异值来近似描述矩阵,采用截断奇异值分解的话,三个小矩阵的维度可以大大的降低。来看个例子。
A
=
[
1
0
0
0
0
0
0
4
0
3
0
0
0
0
0
0
2
0
0
0
]
A = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 2 & 0 & 0 & 0 \\ \end{bmatrix}
A=⎣⎢⎢⎢⎢⎡10002003000000004000⎦⎥⎥⎥⎥⎤
矩阵的秩为3,奇异值应该有3个,如果我们只取前两个奇异值,计算可以得到前两个奇异值为
σ
1
=
4
,
σ
2
=
3
\sigma_1=4, \sigma_2 =3
σ1=4,σ2=3,为了和奇异值相匹配,
U
U
U和
V
V
V也需要做截断,并且计算出其左奇异向量
U
U
U和右奇异向量
V
V
V
U
=
[
0
0
1
0
0
1
0
0
0
0
]
,
D
=
[
4
0
0
3
]
,
V
=
[
0
0
0
1
0
0
1
0
]
U = \begin{bmatrix} 0 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 0 \\ 0 & 0 \end{bmatrix},D = \begin{bmatrix} 4 & 0 \\ 0 & 3 \end{bmatrix},V = \begin{bmatrix} 0 & 0 \\ 0 & 1 \\ 0 & 0 \\ 1 & 0 \end{bmatrix}
U=⎣⎢⎢⎢⎢⎡0100000100⎦⎥⎥⎥⎥⎤,D=[4003],V=⎣⎢⎢⎡00010100⎦⎥⎥⎤
利用截断的奇异值和奇异向量计算出
A
A
A的近似值
A
~
=
U
D
V
T
=
[
0
0
0
0
0
0
0
4
0
3
0
0
0
0
0
0
0
0
0
0
]
\widetilde{A} = UDV^T = \begin{bmatrix}0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 \\ 0 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{bmatrix}
A
=UDVT=⎣⎢⎢⎢⎢⎡00000003000000004000⎦⎥⎥⎥⎥⎤
这可以看成是对原矩阵的近似,不过是在弗罗贝尼乌斯范数(Frobenius)下的近似,即
∣
∣
A
∣
∣
F
=
(
∑
i
=
1
m
∑
j
=
1
n
(
a
i
j
2
)
)
1
2
≈
∣
∣
A
~
∣
∣
F
||A||_F = (\sum_{i=1}^m\sum_{j=1}^n(a_{ij}^2))^{\frac{1}{2}} \approx ||\widetilde A||_F
∣∣A∣∣F=(i=1∑mj=1∑n(aij2))21≈∣∣A
∣∣F
在PCA中的应用
PCA的主要思想是将高维度空间的特征转化成低维度的新特征,在减轻后续算法计算量的同时仍能使算法保持较高的精度,主要用于数据的预处理过程。在此省略掉其计算过程,看一个比较直观的例子,只描述和奇异值分解相关的,本示例来自于李航老师的《统计学习方法》第二版。
例:假设有n个同学参加四门课程的考试,将学生们的考试成绩看作随机变量的取值,对考试成绩进行标准化处理,得到样本的相关矩阵
R
R
R,如下表。
课程 | 语文 | 外语 | 数学 | 物理 |
---|---|---|---|---|
语文 | 1 | 0.44 | 0.29 | 0.33 |
外语 | 0.44 | 1 | 0.35 | 0.32 |
数学 | 0.29 | 0.35 | 1 | 0.60 |
物理 | 0.33 | 0.32 | 0.60 | 1 |
设变量
x
1
,
x
2
,
x
3
,
x
4
x_1,x_2,x_3,x_4
x1,x2,x3,x4分别表示语文、外语、数学、物理的成绩。对样本相关矩阵进行特征值分解,得到相关矩阵的特征值,并按大小排序,
λ
1
=
2.17
,
λ
2
=
0.87
,
λ
3
=
0.57
,
λ
4
=
0.39
\lambda_1 = 2.17, \lambda_2 = 0.87, \lambda_3=0.57, \lambda_4= 0.39
λ1=2.17,λ2=0.87,λ3=0.57,λ4=0.39
这些特征值就是各主成分的方差贡献率,假设要求主成分的累计方差贡献率大于75%,那么只需取前两个主成分即可,即k=2,求出对应于特征值 λ 1 , λ 2 \lambda_1,\lambda_2 λ1,λ2的单位特征向量,列于下表。
项目 | x 1 x_1 x1 | x 2 x_2 x2 | x 3 x_3 x3 | x 4 x_4 x4 | 方差贡献率 |
---|---|---|---|---|---|
y 1 y_1 y1 | 0.460 | 0.476 | 0.523 | 0.537 | 0.543 |
y 2 y_2 y2 | 0.574 | 0.486 | -0.476 | -0.456 | 0.218 |
按照上表可以得到第一、第二主成分:
y
1
=
0.460
x
1
+
0.476
x
2
+
0.523
x
3
+
0.537
x
4
y_1 = 0.460x_1+0.476x_2+0.523x_3+0.537x_4
y1=0.460x1+0.476x2+0.523x3+0.537x4
y 2 = 0.574 x 1 + 0.486 x 2 − 0.476 x 3 − 0.456 x 4 y_2 = 0.574x_1+0.486x_2-0.476x_3-0.456x_4 y2=0.574x1+0.486x2−0.476x3−0.456x4
第一主成分的实际意义在于各门成绩提高都可以使 y 1 y_1 y1成绩提高,也就是第一主成分反映了学生的整体成绩,第二主成分表明文科成绩提高可以使 y 2 y_2 y2提高,理科成绩提高可以使 y 2 y_2 y2降低,第二主成分反映了学生的文科成绩和理科成绩的关系。这就可以将4维的特征降为2维。