- 线性代数一直迷的不行……so
- 这个博客主要参考MIT18.06和MIT18.065两门课的内容,构建并梳理各个知识点的联系。一定要学过一遍线性代数再看,这里着重记录“感觉”,信息密度大,多停多想,不要贪快。遇到不知道且没有解释的的名词可以查下百科。对于没有证明(或许提供了思路)的定理可以当作练习
- MIT18.06讲了线性代数,18.065讲了线性代数在数据科学中的诸多用法。本文计算机视觉和机器学习有关的内容还参考了很多Prince的Computer vision: models, learning, and inference的附录C,鼓励大家看原书,内容非常清楚。本文同时穿插了一些信号、图像方面的知识点
- 这里主要介绍为应用服务的理论工具。部分应用需要少量机器学习基础。如果对某部分应用不感兴趣,可以连带相关章节略去不看
- 作者水平非常有限,如有错误或任何批评,恳请读者一定要指出来
目录
线性变换
(逻辑上应该先讲线性空间,但是会不明意义,所以先从几何角度建立线性变换的感觉,强烈建议去看3B1B 线性代数的本质理解线性变换,无前置基础)
这里给出
T
T
T满足线性变换的条件
T
(
v
+
w
)
=
T
(
v
)
+
T
(
w
)
T
(
c
v
)
=
c
T
(
v
)
\begin{aligned} T(\bm v+\bm w)&=T(\bm v)+T(\bm w) \\ T(c\bm v)&=cT(\bm v) \end{aligned}
T(v+w)T(cv)=T(v)+T(w)=cT(v)
- 只需要确定 T T T对所有基向量的影响,就能完全掌握这个变换
- 两个定理(对几何直觉构建帮助不大,可以先略去)
- 对于同一个元素,在基
α
\alpha
α下坐标为
x
x
x,在基
β
\beta
β下坐标为
x
′
x'
x′,如果
β
=
α
P
\beta = \alpha P
β=αP,那么
x
=
P
x
′
x=Px'
x=Px′.
证明:因为 α x = β x ′ \alpha x = \beta x' αx=βx′. P P P也叫过渡矩阵或基变换矩阵 - 如果
β
=
α
P
\beta = \alpha P
β=αP,线性变换
T
T
T在两组基下的矩阵依次为
A
,
B
A,B
A,B,即
T
(
α
)
=
α
A
,
T
(
β
)
=
β
B
T(\alpha)=\alpha A, T(\beta)=\beta B
T(α)=αA,T(β)=βB,则
B
=
P
−
1
A
P
B=P^{-1}AP
B=P−1AP.
证明: β B = T ( β ) = T ( α P ) = T ( α ) P = α A P = β P − 1 α P \beta B = T(\beta)=T(\alpha P) = T(\alpha)P=\alpha AP=\beta P^{-1}\alpha P βB=T(β)=T(αP)=T(α)P=αAP=βP−1αP,因为 β \beta β满秩,所以 B = P − 1 A P B=P^{-1}AP B=P−1AP
即 A A A和 B B B是两组基下的同一个变换,它们一定相似,且 P P P就是相似变换矩阵
- 对于同一个元素,在基
α
\alpha
α下坐标为
x
x
x,在基
β
\beta
β下坐标为
x
′
x'
x′,如果
β
=
α
P
\beta = \alpha P
β=αP,那么
x
=
P
x
′
x=Px'
x=Px′.
- n n n维线性空间和 R n \mathbb R^n Rn都同构,即维数相等的线性空间都同构,线性空间的结构完全由维数决定
线性空间
数域 F F F中对加法和数乘封闭的非空集合
- Ax=b的解空间不构成线性子空间,注意 b b b的存在,这种情况不满足加法和乘法的封闭法则
子空间的加和交
对于两个子空间
S
S
S和
U
U
U
d
i
m
(
S
)
+
d
i
m
(
U
)
=
d
i
m
(
S
+
U
)
+
d
i
m
(
S
∩
U
)
dim(S) + dim(U) = dim(S+U) + dim(S\cap U)
dim(S)+dim(U)=dim(S+U)+dim(S∩U)
矩阵空间
以矩阵为元素的线性空间
4个基本子空间
- 列空间Column Space C ( A ) C(A) C(A): 所有 A x Ax Ax
- 行空间Raw Space C ( A T ) C(A^T) C(AT)
- 零空间Null space of a matrix N ( A ) N(A) N(A): A x = 0 Ax=0 Ax=0 for all x x x
- 左零空间Left Null space N ( A T ) N(A^T) N(AT): x T A = 0 x^TA=0 xTA=0 for all x x x
关系如图所示(图片摘自网络),这张图非常重要,理解了这张图就理解了很多抽象的概念。
注意
A
x
=
0
Ax=0
Ax=0就已经揭示了
C
(
A
T
)
C(A^T)
C(AT)和
N
(
A
)
N(A)
N(A)是正交的。
A=CR分解,行秩等于列秩
对于 A ∈ R m × n A\in \mathbb R^{m\times n} A∈Rm×n,如果列秩为 r r r, r r r维的列空间中的无关列向量拼成矩阵 C = [ c v 1 , ⋯ , c v r ] C=[c_{v_1}, \cdots, c_{v_r}] C=[cv1,⋯,cvr], A A A的每一列可以由 C C C中的列线性表达,则有 A = C R A=CR A=CR,其中 C ∈ R m × r , R ∈ R r × n C\in \mathbb R^{m\times r},\ R \in \mathbb R^{r\times n} C∈Rm×r, R∈Rr×n. 注意该式也可看作 A A A的每一行由 R R R中的 r r r行进行线性组合得到,所以行空间维度不超过 r r r,也即行秩不超过列秩;类似方法可证列秩不超过行秩。即证明 行 秩 = 列 秩 = r 行秩=列秩=r 行秩=列秩=r。所以 A = C R A=CR A=CR分解中 C C C是列无关的, R R R也是行无关的
- 从 C C C和 R R R两个角度看 A A A的构成,也说明了秩的乘法性质: r a n k [ A B ] ⩽ min { r a n k [ A ] , r a n k [ B ] } rank[AB]\leqslant \min\{rank[A], rank[B] \} rank[AB]⩽min{rank[A],rank[B]}
- A x = b Ax=b Ax=b有解即 b b b在 A A A的列空间中
几个性质
- 当 B B B 可逆时, N ( B A ) = N ( A ) N(BA)=N(A) N(BA)=N(A)
- 行空间和列空间维度都为 r r r
- 行空间与零空间正交,所以零空间 N ( A ) N(A) N(A)维度为 n − r n-r n−r,也即线性方程组 A x = 0 Ax=0 Ax=0解的构成维度;同理左零空间 N ( A T ) N(A^T) N(AT)与列空间正交,维度为 m − r m-r m−r
子空间投影
- 对于向量
a
a
a和
b
b
b,
b
b
b在
a
a
a上投影向量p为
a a T b a T a = a a T a T a b = P b (1) a\frac{a^Tb}{a^T a}=\frac{aa^T}{a^Ta}b=Pb \tag{1} aaTaaTb=aTaaaTb=Pb(1)
其中 P P P看作是投影矩阵,可以自行证明。
最小二乘问题
- 求解 A x = b Ax=b Ax=b时,如果无解, b b b不在 C ( A ) C(A) C(A)当中,那么改求 A x ^ = p A \hat x=p Ax^=p,其中 p p p是 b b b在 C ( A ) C(A) C(A)上的投影。注意到 b − A x ^ b - A \hat x b−Ax^垂直于 C ( A ) C(A) C(A),在Left Null Space。所以 A T ( b − A x ^ ) = 0 A^T(b-A\hat x)=0 AT(b−Ax^)=0,即 A T A x ^ = A T b A^TA\hat x=A^T b ATAx^=ATb
- 所以 x ^ = ( A T A ) − 1 A T b \hat x = (A^TA)^{-1}A^T b x^=(ATA)−1ATb, p = A x ^ = A ( A T A ) − 1 A T b = P b p=A\hat x=A (A^TA)^{-1}A^T b=Pb p=Ax^=A(ATA)−1ATb=Pb,注意该式和式(1)的相似性. P = A ( A T A ) − 1 A T P=A(A^TA)^{-1}A^T P=A(ATA)−1AT称之为投影矩阵. x ^ \hat x x^也称之为 A A A的左逆
- 注意 P P P的两个性质,这两个性质也是判定投影矩阵的充分条件: P T = P , P 2 = P P^T=P, P^2=P PT=P,P2=P,(注意 P 2 = P P^2=P P2=P不足以判定投影矩阵,例如矩阵 ( 1 , 1 ; 0 , 0 ) (1,1;0,0) (1,1;0,0)). 因为 P 2 = P P^2=P P2=P,特征值只能为0或1. 此外, I − P I-P I−P也是投影矩阵
- 因为 b = p + e b=p+e b=p+e. ( e e e是投影误差),而 p = P b p=Pb p=Pb,所以 e = ( I − P ) b e=(I-P)b e=(I−P)b. 且有 ( I − P ) 2 = I − P , ( I − P ) T = I − P (I-P)^2=I-P, (I-P)^T=I-P (I−P)2=I−P,(I−P)T=I−P. 注意和 P P P的类似性
- A T A A^TA ATA 可逆当且仅当 A A A的所有列独立. 一个可靠的证明方法是两者的Null Space相同. 这里似乎 r ( A T A ) = r ( A ) r(A^TA)=r(A) r(ATA)=r(A)是成立的
分析角度看最小二乘问题
A
x
=
b
Ax=b
Ax=b
以最小二乘形式看待
x
^
=
arg min
x
[
(
A
x
−
b
)
T
(
A
x
−
b
)
]
=
arg min
x
[
x
T
A
T
A
x
−
2
x
T
A
T
b
]
\begin{aligned} \hat x &= \argmin_x [(Ax-b)^T(Ax-b)] \\ &= \argmin_x [x^TA^TAx - 2x^TA^Tb] \end{aligned}
x^=xargmin[(Ax−b)T(Ax−b)]=xargmin[xTATAx−2xTATb]
标量对
x
x
x求导得
x
=
(
A
T
A
)
−
1
A
T
b
x=(A^TA)^{-1}A^Tb
x=(ATA)−1ATb
结果和投影方法是一样的。
其中
A
A
A的行数如果比
x
x
x少,则奇异
应用:线性回归
这一节摘自PRML P143,为了与PRML配图一致,符号标注不符合本文的惯用标注。
几何解释
记训练集标注为 t = ( t 1 , . . . , t N ) T \bf t = (t_1, ..., t_N)^T t=(t1,...,tN)T,并构成标注空间 R N \mathbb R^N RN, S \mathcal{S} S是能在训练集的标注空间中用广义线性回归张成的超平面,也即训练集的列空间
- 这里线性回归的基可以是带核
φ
(
X
)
\varphi (X)
φ(X)的,实际上带核的仍然是张成超平面,而不是曲面,超平面的第
i
i
i个基由
φ
i
(
X
)
\varphi_i(X)
φi(X)决定,
φ
i
\varphi_i
φi表示第
i
i
i个特征,
X
X
X表示所以的N个数据。
这样线性回归是求了标注空间中训练集所在位置在超平面上的投影,垂直距离即为误差向量,范数为最小二乘的结果。(图片摘自Bishop的Pattern Recognition and Machine Learning)
- 注意误差向量垂直于 S \mathcal{S} S,所以如果截距项全1向量在 S S S当中,那么必有误差向量自身之和为0. 也即线性回归拟合均值一定等于标注均值
多重共线性缺陷
之前只知道多重共线性不好,到底哪里不好一直说不清楚。这里把它讲清楚。
多重共线性的灾难在于参数值爆炸。
我们记(可选经过核变换后)训练集为
Φ
∈
R
N
×
M
\Phi \in \mathbb{R}^{N \times M}
Φ∈RN×M,其中
M
M
M是特征维度。
r
a
n
k
[
Φ
]
<
M
rank[\Phi]<M
rank[Φ]<M时,即产生了多重共线性问题,也即特征之间线性相关。因为
r
a
n
k
[
Φ
]
=
r
a
n
k
[
Φ
T
Φ
]
=
r
a
n
k
[
Φ
Φ
T
]
rank[\Phi] = rank[\Phi^T \Phi] = rank[\Phi \Phi^T]
rank[Φ]=rank[ΦTΦ]=rank[ΦΦT],(注:方法为证明
Φ
x
=
0
\Phi x =0
Φx=0与
Φ
T
Φ
x
=
0
\Phi^T \Phi x= 0
ΦTΦx=0同解)。所以如何判断
r
a
n
k
[
Φ
]
rank[\Phi]
rank[Φ]与
M
M
M的关系,只需要计算
Φ
T
Φ
\Phi^T \Phi
ΦTΦ是否奇异。
- 此外,如果 Φ T Φ \Phi^T \Phi ΦTΦ接近奇异,即行列式很小,那么线性回归的参数闭式解 ( Φ T Φ ) − 1 Φ T t (\Phi^T \Phi)^{-1} \Phi^T \bf t (ΦTΦ)−1ΦTt会非常大。关于行列式和逆的定义参考下一节。
- 从几何角度解释,即两个基向量方向非常近,那么为了表达出与这两个基向量几乎垂直的方向上的位置,这两个向量需要不断抵消, 系数会增长非常快!
这里如果无法求逆,则可以求伪逆,后文伪逆部分会继续讨论最小二乘求逆问题。
行列式和逆
最基本的性质
- d e t [ I ] = 1 det[I]=1 det[I]=1
- 交换两行,行列式取反
- a) 某一行乘
k
k
k倍,行列式乘
k
k
k倍.
b) d e t [ ⋮ a i 1 + a i 2 ⋮ ] = d e t [ ⋮ a i 1 ⋮ ] + d e t [ ⋮ a i 2 ⋮ ] det \begin{bmatrix}\vdots \\ a_{i_1}+a_{i_2} \\ \vdots\end{bmatrix}=det \begin{bmatrix}\vdots \\ a_{i_1} \\ \vdots\end{bmatrix} +det \begin{bmatrix}\vdots \\ a_{i_2} \\ \vdots\end{bmatrix} det⎣⎢⎢⎡⋮ai1+ai2⋮⎦⎥⎥⎤=det⎣⎢⎢⎡⋮ai1⋮⎦⎥⎥⎤+det⎣⎢⎢⎡⋮ai2⋮⎦⎥⎥⎤
用这三条性质能推出剩下一大堆东西,包括基本算法和其他性质,都能推出来,列几个常用的性质
- d e t [ A T ] = d e t [ A ] det[A^T]=det[A] det[AT]=det[A]
- d e t [ A B ] = d e t [ A ] d e t [ B ] det[AB]=det[A]det[B] det[AB]=det[A]det[B],所以 d e t [ A − 1 ] = 1 / d e t [ A ] det[A^{-1}]=1/det[A] det[A−1]=1/det[A]
- 利用性质2和3可以证明
d
e
t
[
⋮
a
i
1
⋮
a
i
2
⋮
]
=
d
e
t
[
⋮
a
i
1
⋮
a
i
1
+
a
i
2
⋮
]
det \begin{bmatrix}\vdots \\ a_{i_1} \\ \vdots \\ a_{i_2} \\\vdots\end{bmatrix}=det \begin{bmatrix}\vdots \\ a_{i_1} \\ \vdots \\ a_{i_1}+a_{i_2} \\\vdots\end{bmatrix}
det⎣⎢⎢⎢⎢⎢⎢⎢⎡⋮ai1⋮ai2⋮⎦⎥⎥⎥⎥⎥⎥⎥⎤=det⎣⎢⎢⎢⎢⎢⎢⎢⎡⋮ai1⋮ai1+ai2⋮⎦⎥⎥⎥⎥⎥⎥⎥⎤
由此进一步可证不满秩矩阵的行列式为0。
行列式表示矩阵组成的体积
可以证明体积也满足上述三个重要性质。于是得证。(3b不是很好证)
建议参考3B1B——线性代数的本质
- 三阶行列式即混合积: a × ( b ⋅ c ) = d e t [ [ a , b , c ] ] a\times (b\cdot c)=det[[a,b,c]] a×(b⋅c)=det[[a,b,c]]
行列式算法
对于方阵
A
∈
R
n
×
n
A \in \mathbb{R}^{n\times n}
A∈Rn×n,一种写法:
d
e
t
[
A
]
=
∑
n
!
terms
±
a
1
α
a
2
β
⋯
a
n
ω
det[A]=\sum_{n! \text{ terms}} \pm a_{1\alpha}a_{2\beta}\cdots a_{n\omega}
det[A]=n! terms∑±a1αa2β⋯anω
其中
(
α
,
β
,
γ
,
⋯
,
ω
)
=
perm of
(
1
,
2
,
⋯
,
n
)
(\alpha, \beta,\gamma, \cdots,\omega) = \text{perm of } (1,2,\cdots,n)
(α,β,γ,⋯,ω)=perm of (1,2,⋯,n)
另一种写法:
d
e
t
[
A
]
=
a
11
C
11
+
a
12
C
12
+
⋯
+
a
1
n
C
1
n
det[A]=a_{11}C_{11} + a_{12} C_{12} + \cdots + a_{1n} C_{1n}
det[A]=a11C11+a12C12+⋯+a1nC1n
这里只以第一行为例,也可以取其他行。
C
i
j
C_{ij}
Cij是第
i
i
i行第
j
j
j列的代数余子式cofactor,注意有正有负。
逆矩阵
A
C
T
=
d
e
t
[
A
]
I
AC^T=det[A]I
ACT=det[A]I
注意
C
C
C是代数余子式形成的矩阵,
C
T
C^T
CT称之为伴随矩阵Adjugate Matrix. 正确性可以自行验证.
- 伴随矩阵和逆矩阵只差一个系数。然而,伴随矩阵对不可逆的矩阵也有定义,并且不需要用到除法(摘自百度百科)
克拉默法则
A
x
=
b
Ax=b
Ax=b
有唯一解时
x
=
A
−
1
b
=
1
d
e
t
[
A
]
C
T
b
x=A^{-1}b=\frac{1}{det[A]}C^Tb
x=A−1b=det[A]1CTb
注意
C
T
b
C^Tb
CTb中代数余子式的含义,可以得到
x
i
=
d
e
t
[
B
i
]
d
e
t
[
A
]
B
i
=
[
A
⋅
1
,
⋯
,
A
⋅
(
i
−
1
)
,
b
,
A
⋅
(
i
+
1
)
,
⋯
,
A
⋅
n
]
x_i=\frac{det[B_i]}{det[A]} \\ B_i = \begin{bmatrix} A_{\cdot 1}, \cdots, A_{\cdot (i-1)}, b, A_{\cdot (i+1)}, \cdots, A_{\cdot n} \end{bmatrix}
xi=det[A]det[Bi]Bi=[A⋅1,⋯,A⋅(i−1),b,A⋅(i+1),⋯,A⋅n]
几何解释参考3B1B——线性代数的本质
特殊等式
- d e t [ I n + u v T ] = d e t [ 1 + u T v ] det[I_n+uv^T]=det[1+u^Tv] det[In+uvT]=det[1+uTv]。它的扩展形式为 d e t [ I n + A B T ] = d e t [ I m + A T B ] det[I_n+AB^T]=det[I_m+A^T B] det[In+ABT]=det[Im+ATB],其中 A , B ∈ R n × m A,B\in \mathbb R^{n\times m} A,B∈Rn×m. 只需证明左右特征值连乘相等即可。注意 A T B A^TB ATB与 A B T AB^T ABT特征值相等(Sylvester定理)
正交矩阵
Q
Q
Q是方阵,且每一列的范数为1,不同列正交。故有
Q
T
Q
=
I
,
Q
−
1
=
Q
T
Q^TQ=I,Q^{-1}=Q^T
QTQ=I,Q−1=QT
例子:哈达玛矩阵,参考博客
- 实际上每一行范数也为1,也与其他行正交
- 复空间的正交矩阵叫做酉矩阵unitary,即满足 A ˉ T A = I \bar A^TA=I AˉTA=I
- 正交矩阵特征值绝对值都为1,因为 ( Q x ) ‾ T Q x = λ ˉ λ x ˉ T x \overline {(Qx)}^TQx=\bar \lambda\lambda \bar x^Tx (Qx)TQx=λˉλxˉTx, x ˉ T x > 0 \bar x^T x>0 xˉTx>0,所以 ∥ λ ∥ = 1 \|\lambda\|=1 ∥λ∥=1. 注意旋转矩阵这一类,考虑其旋转和特征向量的物理意义,除了 I I I外,无法找出实特征向量和特征值
- 注意:如果 Q T Q = I Q^TQ=I QTQ=I,当 Q Q Q不为方阵的时候,不一定能推出 Q Q T = I QQ^T=I QQT=I
旋转矩阵与正交变换
正交矩阵行列式绝对值为1,所以线性变换体积不变。左乘的线性变换效果是可看作是一种反演(改变手性)+旋转,这个在SVD中会有感受。行列式为1时,为旋转矩阵,为-1则包含了反演
正交矩阵的线性变换叫正交变换,正交变换有以下好性质:
- 距离不变
- 夹角不变
- 上述两条也能推出内积不变,即 x 1 T x 2 = ( A x 1 ) T ( A x 2 ) x_1^Tx_2=(Ax_1)^T(Ax_2) x1Tx2=(Ax1)T(Ax2)
反射矩阵
即正交对称阵 r T = r r^T=r rT=r.
- 注意这意味着特征值都为实数±1,且特征向量都存在。这是一个不带旋转的变换,如果带旋转的话不可能特征值都为±1
A=QR与Gram-schmitt正交化
A=QR
Q
Q
Q是正交矩阵,其中
R
R
R是上三角矩阵.
这个上三角矩阵成立的原因可以思考一下Gram-schmitt正交化的过程,对
Q
Q
Q的列向量
q
1
,
q
2
,
⋯
,
q
n
q_1,q_2,\cdots, q_n
q1,q2,⋯,qn一个一个进行正交。每一次正交只依赖于前面的列向量
数值计算:Gram-schmitt中的过小分母问题
- 在Gram-schmitt时,会有分母 q i T q i q_i^T q_i qiTqi的出现,如果两个向量离得很近,这个值有很小,会在数值计算中出现问题。一种解决方法是:假定上一个被正交化的是 q i q_i qi,在计算下一个要被正交化的向量 q j q_j qj时,从没有正交化的向量当中把 a ⋅ − q i q i T q i T q i a ⋅ a_\cdot - \frac{q_i q_i^T}{q_i^T q_i}a_\cdot a⋅−qiTqiqiqiTa⋅都计算一遍,选范数最大的。注意这一步计算不引入额外开销,因为这个计算是必须要经历的。
应用:信号处理中的变换
信号处理中有一大类变换是线性变换,其中的部分变换是正交变换,可以看作是对原始信号 x x x施加了一组新的标准正交基,例如 y = A ˉ T x y=\bar A^Tx y=AˉTx,新的基为正交矩阵 A A A中的各列,而且 A A A和 A T A^T AT互为正反变换。实际例子包括:
- 傅里叶变换,记
ω
=
e
2
π
/
n
\omega=e^{2\pi/n}
ω=e2π/n,则
A = 1 n [ 1 1 1 ⋯ 1 1 ω ω 2 ⋯ ω n − 1 1 ω 2 ω 4 ⋯ ω 2 ( n − 1 ) ⋮ 1 ω n − 1 ω 2 ( n − 1 ) ⋯ ω ( n − 1 ) ( n − 1 ) ] A=\frac{1}{\sqrt n}\begin{bmatrix} 1&1&1&\cdots&1 \\ 1 & \omega & \omega^2 & \cdots& \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \cdots &\omega^{2(n-1)} \\ &&&\vdots& \\ 1&\omega^{n-1}&\omega^{2(n-1)}&\cdots& \omega^{(n-1)(n-1)} \end{bmatrix} A=n1⎣⎢⎢⎢⎢⎢⎡11111ωω2ωn−11ω2ω4ω2(n−1)⋯⋯⋯⋮⋯1ωn−1ω2(n−1)ω(n−1)(n−1)⎦⎥⎥⎥⎥⎥⎤
注意 A A A是复数矩阵,虽然对称,但不是共轭对称,所以不是埃尔米特矩阵,但它是酉矩阵,反变换为 A ˉ T \bar A^T AˉT - 小波变换。参考我的另一篇博客:小波变换——公式整理和简单介绍
特征值与特征向量
对于方阵
A
∈
R
n
×
n
A \in \mathbb{R}^{n\times n}
A∈Rn×n
A
v
=
λ
v
Av=\lambda v
Av=λv
v
v
v为特征向量,
λ
\lambda
λ为对应特征值
- ∑ λ = t r [ A ] \sum \lambda=tr[A] ∑λ=tr[A], t r tr tr表示矩阵的迹trace:对角线元素之和
- ∏ λ = d e t [ A ] \prod \lambda = det[A] ∏λ=det[A]
- 对称矩阵特征值一定为实数,反对称矩阵则为零或纯虚数;其他矩阵介于对称与反对称之间
- A A A和 A T A^T AT的特征值相同,但特征向量不一定,因为Left Null Space not equals to Null Space.
- Sylvester定理:
A
B
AB
AB和
B
A
BA
BA的非零特征值相等
(摘自程云鹏《矩阵论》,笔记来自zhenliang师兄。感谢!)
特征值分解
A
=
V
Λ
V
−
1
A=V\Lambda V^{-1}
A=VΛV−1
其中,
V
V
V是特征向量列拼接组成的矩阵,
Λ
\Lambda
Λ是对应特征值。这里需要假设
n
n
n个线性无关的特征向量存在,即
A
A
A可逆。
- A k = V Λ k V − 1 A^k=V \Lambda ^k V^{-1} Ak=VΛkV−1,所以 A k A^k Ak特征向量不变,特征值乘 k k k次。
- A k u = ( V Λ k ) ( V − 1 u ) = ( λ 1 v 1 , ⋯ , λ n v n ) ( c 1 , ⋯ , c n ) T = ∑ i = 1 n c i λ i k v i A^k u=(V\Lambda^k)( V^{-1}u)=(\lambda_1 v_1,\cdots,\lambda_n v_n)(c_1,\cdots,c_n)^T=\sum_{i=1}^{n}c_i \lambda_i^k v^i Aku=(VΛk)(V−1u)=(λ1v1,⋯,λnvn)(c1,⋯,cn)T=∑i=1nciλikvi
- 若 A A A的所有特征值都不同,那么一定有 n n n个线性无关特征向量(可对角化);但特征值相同时(代数重数大于1),需要再研究,线性无关的特征向量可能不够,即几何重数低于代数重数
- 特征向量组唯一当且仅当所有的特征值都是唯一的
- 给定特征值和特征向量唯一确定一个矩阵,可以从几何角度感知,参考3B1B线性代数本质
实对称矩阵的特征值分解
谱定理Spectral Theorem
如果
S
S
S是实对称矩阵,
S
S
S一定可以分解成
n
n
n个线性无关的正交实特征向量和实特征值(不一定都不同)
S
=
Q
Λ
Q
T
S = Q\Lambda Q^T
S=QΛQT
正规矩阵Normal Matrix
实际上当
A
A
A满足
A
A
ˉ
T
=
A
ˉ
T
A
A\bar A^T=\bar A^TA
AAˉT=AˉTA时,特征向量都正交,这样的矩阵叫正规矩阵normal matrix,其满足
A
=
Q
Λ
Q
T
A=Q\Lambda Q^T
A=QΛQT
正规矩阵的典型特例:
- 实对称矩阵, λ ∈ R \lambda \in \mathbb R λ∈R
- 正交矩阵, ∣ λ ∣ = 1 |\lambda|=1 ∣λ∣=1
- A = A ˉ T A=\bar A^T A=AˉT,其中 A ˉ \bar A Aˉ是 A A A的共轭矩阵。满足上述条件的矩阵也叫埃尔米特矩阵Hermitian Matrix或自共轭矩阵
- A A A是反对称矩阵,特征向量也正交,注意特征值一定为0或纯虚数
主轴定理
Q
Q
Q是
A
A
A的特征向量组成的正交矩阵,作用在单位圆上后,可以将
A
A
A看作是沿
v
i
v_i
vi方向扩展
λ
i
\lambda_i
λi倍。如图所示(摘自GoodFellow《深度学习》)
解释:
A
u
=
Q
Λ
Q
T
u
Au=Q\Lambda Q^T u
Au=QΛQTu,从右往左看,这里正交矩阵
Q
T
Q^T
QT对单位圆
u
u
u的旋转作用看不出来,
Λ
\Lambda
Λ把坐标轴拉长,
Q
Q
Q再把拉长后的椭圆旋转
- Λ \Lambda Λ通常按降序排列
- 矩阵奇异当且仅当有0特征值
- 所有特征值为正,则正定;不为负,则半正定;不为正,则半负定;为负,则负定。感受一下正负定对线性变换的影响
实对称矩阵的谱分解
实对称矩阵的特征分解:
S
=
(
Q
Λ
)
Q
T
=
∑
i
λ
i
q
i
q
i
T
S=(Q\Lambda) Q^T = \sum_i \lambda_i q_i q_i^T
S=(QΛ)QT=i∑λiqiqiT
注意谱是对称的,并且全是秩一矩阵
且有
S
q
i
=
(
∑
i
λ
i
q
i
q
i
T
)
q
i
=
λ
q
i
Sq_i =( \sum_i \lambda_i q_i q_i^T) q_i= \lambda q_i
Sqi=(i∑λiqiqiT)qi=λqi
应用:瑞利商Rayleigh Quotient
R
(
S
,
x
)
=
x
T
S
x
x
T
x
R(S,x)=\frac{x^TSx}{x^Tx}
R(S,x)=xTxxTSx
其中
S
S
S是实对称矩阵。
-
易证 λ min ⩽ R ( S , x ) ⩽ λ max \lambda_{\min}\leqslant R(S,x)\leqslant \lambda_{\max} λmin⩽R(S,x)⩽λmax, x x x分别指向最大和最小特征值对应的特征向量时,取最大和最小值。
-
当 x x x指向其他特征向量方向时, x x x位于鞍点saddle point,可以把 x + Δ x x + \Delta x x+Δx用特征向量表示,然后发现函数增减和 Δ x \Delta x Δx方向有关。不过鞍点可以是maxmin问题的解,例如当 S S S是对角阵时,对于 i i i维子空间 V i V^{i} Vi, max V i min x ∈ V i R ( S , x ) = λ i \max_{V^i}\min_{x \in V^i}R(S,x)=\lambda_i Vimaxx∈ViminR(S,x)=λi
其中 λ i \lambda_i λi是第 i i i大的特征值。 -
广义瑞利商:
R ( S , B , x ) = x T S x x T B x = x ~ T ( B − 1 / 2 S B − 1 / 2 ) x ~ x ~ T x ~ R(S, B,x)=\frac{x^TSx}{x^TBx}=\frac{\tilde x^T(B^{-1/2}SB^{-1/2}) \tilde x}{\tilde x^T \tilde x} R(S,B,x)=xTBxxTSx=x~Tx~x~T(B−1/2SB−1/2)x~
仍然能化成一般形式 -
很多聚类或降维任务,包括线性判别分析、主成分分析等,经过推导,往往能推导出瑞利商的形式
应用:马尔可夫矩阵Markov Matrix(随机矩阵Stochastic Matrix)与转移过程
马尔可夫矩阵 A A A:
- 每列的和为1
- 所有元素大于0
物理意义是列到行概率转移。
- 上述定义1保证了1是特征值,原因在于 A − I A-I A−I奇异. 并且其对应的特征向量每个值都不小于0(总之不能异号),根据马尔可夫矩阵的物理意义可以理解。特征值为1的特征向量也是稳态 A ∞ u 0 A^{\infty}u_0 A∞u0的组成部分
- 两条定义共同保证了特征值不大于1. 注意马尔可夫矩阵在无限次放后的稳态性, A n = V Λ k V − 1 A^n=V \Lambda ^k V^{-1} An=VΛkV−1,特征值必须收敛
数值计算:快速特征值计算方法
记 A 0 = A A_0=A A0=A,做QR分解得到 A 0 = Q 0 R 0 A_0=Q_0 R_0 A0=Q0R0,取 A 1 = R 0 Q 0 A_1=R_0 Q_0 A1=R0Q0. 则 A 1 = R 0 A 0 R − 1 A_1 = R_0 A_0 R^{-1} A1=R0A0R−1,和 A 0 A_0 A0相似,所以特征值相同,同样方法构造 A 2 , A 3 , ⋯ A_2, A_3, \cdots A2,A3,⋯,下三角部分会越来越小(证明略,可能比较复杂),从而对角线越来越接近特征值。
改进
引入平移矩阵,得到 A 0 − s I A_0-sI A0−sI,进行QR分解,取 A 1 = R 0 Q 0 + s I = R 0 ( A 0 − s I ) R 0 − 1 + s I = R 0 A 0 R 0 − 1 A_1=R_0 Q_0 +sI=R_0 (A_0 - sI )R_0^{-1}+ sI=R_0 A_0R_0^{-1} A1=R0Q0+sI=R0(A0−sI)R0−1+sI=R0A0R0−1,会发现这样这样的平移不改变特征值。这样做的好处是得到的特征值收敛更快。
这一节内容参考了博客:MIT 18.065—机器学习中的矩阵方法12 计算特征值和奇异值
这样的计算方法在SVD分解中也有应用。MATLAB对这一类问题就用的是这种方法。这里Mark一下,不再展开。
相似矩阵
A
∼
B
⇔
A
=
M
B
M
−
1
A \sim B \Leftrightarrow A=MBM^{-1}
A∼B⇔A=MBM−1
注意相似的传递性,这有点像是抽象代数里的群作用在集合上形成的“轨道”,相似的矩阵在一个轨道里
- 相似矩阵特征值相同。对于 A x = λ x Ax=\lambda x Ax=λx,有 M − 1 A M M − 1 x = λ M − 1 x M^{-1}AMM^{-1}x=\lambda M^{-1}x M−1AMM−1x=λM−1x,进而 B ( M − 1 x ) = λ M − 1 x B(M^{-1}x)=\lambda M^{-1}x B(M−1x)=λM−1x,这说明 A A A的特征值 λ \lambda λ也是 B B B的特征值。特征向量不同。
- 上述命题反过来不成立,如果特征值有重根,特征向量个数有可能不同。当特征向量不够的时候,甚至无法对角化。例如
[ 4 1 0 4 ] \begin{bmatrix} 4 & 1\\ 0 & 4 \end{bmatrix} [4014]
这也叫Jordan标准型Jordan form - 每个矩阵都像相似于一个Jordan矩阵,Jordan矩阵由Jordan块组成。而Jordan块数则等于不同的特征向量数。这里仅mark一下,不再展开,等用到时再详细了解吧。
- A B AB AB和 B A BA BA特征值相同,证明方法是 B ( A B ) B − 1 = B A B(AB)B^{-1}=BA B(AB)B−1=BA,两者相似,这里假定 B B B可逆. SVD和这里有一定的相通之处,即 A T A , A A T A^TA, AA^T ATA,AAT有相同非零特征值。
正定矩阵
一种特殊的实对称矩阵。几个等价定义:
- ∀ x ≠ 0 , x T S x > 0 \forall x\neq 0, x^TSx>0 ∀x=0,xTSx>0,即二次型quadratic form(相对于线型linear form)大于0
- 特征值全为正,即 S = Q Λ Q T S=Q\Lambda Q^T S=QΛQT中 Λ \Lambda Λ对角线元素全正。所以迹恒正
- 主元pivots全为正。主元实际上和LU分解,和二次型配方都是一一对应的
- 顺序主子式leading determinants全为正。注意主元和顺序主子式的关系,前 k k k主元的乘积是第 k k k顺序主子式,可以用高斯消元法解释,高斯消元法不会影响行列式。如果有一个主元为负,那么该主子式为负
-
S
=
A
T
A
S=A^TA
S=ATA,
A
A
A列无关.
如果 A A A不满秩,亦能半正定positive semi definite.考虑定义1和 N ( A ) N(A) N(A)易证. 所以协方差总是半正定,就像是方差总是正数
性质:
- 如果 A , B A,B A,B都正定,那么 A + B A+B A+B正定,用定义1秒证
- 正定矩阵可逆,且逆矩阵正定
说明:
- 二次型 x T S x x^TSx xTSx也叫做能量Energy,许多机器学习和深度学习的问题就是在找参数 x x x来最小化能量
- 二阶偏导数的Hessian矩阵半正定时,原函数是凸函数convex function
应用:图的拉普拉斯矩阵Laplacian Matrix
对于
v
v
v个点的图简无向单simple graph,拉普拉斯矩阵
L
∈
R
v
×
v
L \in \mathbb R^{v\times v}
L∈Rv×v存储了图的信息
L
=
D
−
A
L=D-A
L=D−A
其中
D
∈
R
v
×
v
D \in \mathbb R^{v\times v}
D∈Rv×v是对角矩阵,表示每个顶点的度;
A
A
A是图的邻接矩阵,仅由1或0组成。
-
L
L
L一定是半正定矩阵。
证明方法为 L = M T M L=M^TM L=MTM,其中 M ∈ R e × v M\in \mathbb R^{e\times v} M∈Re×v是关联矩阵Incidence matrix, e e e是图的边数, M i ⋅ M_{i\cdot} Mi⋅是第 i i i条边(连接 p , q p,q p,q两点,且 p > q p > q p>q), M i p = 1 , M i q = − 1 M_{ip}=1, M_{iq}=-1 Mip=1,Miq=−1,行中其他元素为0.
另外注意 r a n k [ L ] = v − 1 rank[L]=v-1 rank[L]=v−1,所以只能半正定. L L L有特征值0,对应特征向量为 1 n \textbf 1_{n} 1n - 拉普拉斯矩阵还有一些性质,和图的部分数量有关,另外和谱聚类等方法有联系,可以自行查阅相关资料,如ESL第14章。
扩展:隔行扫描定理Cauchy Interlacing Theorem
该定理揭示了秩1正定矩阵所带来的变化会如何影响特征值。
对于对称矩阵
S
S
S,特征值为
λ
1
⩾
λ
2
⩾
⋯
\lambda_1\geqslant \lambda_2\geqslant \cdots
λ1⩾λ2⩾⋯,对于
S
+
θ
u
u
T
S+\theta uu^T
S+θuuT得到特征值
μ
1
⩾
μ
2
⩾
⋯
\mu_1\geqslant \mu_2\geqslant \cdots
μ1⩾μ2⩾⋯,其中
θ
\theta
θ是一个系数,
u
u
u是单位向量。
u
u
T
uu^T
uuT一定正定,并可以看作是positive change,会导致特征值上升。但是特征值的上升有上界,实际上有
λ
1
+
θ
⩾
μ
1
⩾
λ
1
⩾
μ
2
⩾
λ
2
⩾
⋯
(3)
\lambda_1 + \theta \geqslant \mu_1\geqslant \lambda_1\geqslant \mu_2\geqslant \lambda2 \geqslant \cdots \tag{3}
λ1+θ⩾μ1⩾λ1⩾μ2⩾λ2⩾⋯(3)
如何证明?
Weyl不等式
(该定理不再证明)
对于对称矩阵
S
,
T
S,T
S,T,降特征值降序排列,有
λ
i
+
j
−
1
S
+
T
⩽
λ
i
S
+
λ
j
T
\lambda^{S+T}_{i+j-1}\leqslant \lambda^S_{i} + \lambda^T_j
λi+j−1S+T⩽λiS+λjT.
- 当
j
=
1
j=1
j=1时,
λ
i
S
+
T
⩽
λ
i
S
+
λ
max
T
\lambda^{S+T}_{i}\leqslant \lambda^S_{i} + \lambda^T_{\max}
λiS+T⩽λiS+λmaxT.
可以发现式(3)中 λ 1 + θ ⩾ μ 1 \lambda_1 + \theta \geqslant \mu_1 λ1+θ⩾μ1
该定理给出了对称矩阵带来变化的特征值变化上界 - 当
j
=
2
j=2
j=2时,
λ
i
+
1
S
+
T
⩽
λ
i
S
+
λ
2
T
\lambda^{S+T}_{i+1} \leqslant \lambda^S_i + \lambda^T_{2}
λi+1S+T⩽λiS+λ2T
式(3)中秩1矩阵的 λ 2 T = 0 \lambda^T_2=0 λ2T=0,所以序关系得证
奇异值分解Singular Value Decomposition(SVD)
奇异值分解非常重要,在后文应用非常多。
对任意
m
×
n
m\times n
m×n矩阵
A
A
A,总能
A
=
U
Σ
V
T
A=U\Sigma V^T
A=UΣVT
其中,
U
U
U是
m
m
m阶正交矩阵,
Σ
\Sigma
Σ是
m
×
n
m\times n
m×n对角矩阵且
σ
1
⩾
σ
2
⩾
⋯
⩾
0
\sigma_1 \geqslant \sigma_2 \geqslant \cdots \geqslant 0
σ1⩾σ2⩾⋯⩾0,
V
V
V是
n
n
n阶正交矩阵。奇异值分解并不唯一,往往让
Σ
\Sigma
Σ对角元素沿着从大到小的顺序排,在某些情况下唯一。参考:可逆矩阵的奇异值分解唯一吗?存在奇异值分解唯一的矩阵吗?
- 性质: A T A = V ( Σ T Σ ) V T A^TA=V(\Sigma^T\Sigma)V^T ATA=V(ΣTΣ)VT, A A T = U ( Σ Σ T ) U T AA^T=U(\Sigma\Sigma^T)U^T AAT=U(ΣΣT)UT
- 如果 A A A是正定,则SVD即是正交特征分解
代数解释
如果
r
a
n
k
[
A
]
=
r
rank[A]=r
rank[A]=r,奇异值分解让
A
A
A把行空间
C
(
A
T
)
C(A^T)
C(AT)的一组
r
r
r个单位正交基同构映射到列空间
C
(
A
)
C(A)
C(A)中去得到另一组
r
r
r个单位正交基,注意
C
(
A
)
C(A)
C(A)和
C
(
A
T
)
C(A^T)
C(AT)维度相同,即
A
[
v
1
,
v
2
,
⋯
,
v
r
]
=
[
u
1
,
u
2
,
⋯
,
u
r
]
[
σ
1
⋱
σ
r
0
⋱
0
]
(2)
A[ v_1, v_2, \cdots , v_r] =[u_1,u_2,\cdots,u_r]\begin{bmatrix} \sigma_1 \\ &\ddots \\ &&\sigma_r \\ &&&0 \\ &&&&\ddots \\ &&&&&0 \end{bmatrix} \tag{2}
A[v1,v2,⋯,vr]=[u1,u2,⋯,ur]⎣⎢⎢⎢⎢⎢⎢⎡σ1⋱σr0⋱0⎦⎥⎥⎥⎥⎥⎥⎤(2)
也可写成
A
V
=
U
Σ
AV=U\Sigma
AV=UΣ
其中
V
V
V中多余的列如
v
r
+
1
,
⋯
,
v
n
v_{r+1}, \cdots, v_n
vr+1,⋯,vn从
A
A
A的零空间
N
(
A
)
N(A)
N(A)中拿即可,其维度正好是
n
−
r
n-r
n−r,而且所有列都正交;
U
U
U中多余的列类似,从左零空间
N
(
A
T
)
N(A^T)
N(AT)中拿. 对于
i
>
r
i>r
i>r,
Σ
i
i
=
0
\Sigma_{ii} = 0
Σii=0,
r
r
r列之后等式两边全是0
-
如何找出 V V V的前 r r r列?
从 A T A A^TA ATA的非零特征值对应的特征向量中拿, A T A = V ( Σ T Σ ) V T A^TA=V(\Sigma^T \Sigma)V^T ATA=V(ΣTΣ)VT, V V V也即实对称矩阵 A T A A^TA ATA的特征向量,从中也能看出 Σ \Sigma Σ;- σ i 2 v i = A T ( A v i ) \sigma_i^2 v_i=A^T(Av_i) σi2vi=AT(Avi),所以 v i v_i vi一定在行空间 C ( A T ) C(A^T) C(AT)中
-
如何找出 U U U的前 r r r列?
令 u i = A v i σ i u_i=\frac{Av_i}{\sigma_i} ui=σiAvi- 可以看出 u i u_i ui一定在 A A A的列空间里
- 可以证明这
r
r
r列是单位正交的,不妨记前
r
r
r列的各个矩阵为
V
r
∈
R
n
×
r
,
U
r
∈
R
m
×
r
,
Σ
r
∈
R
r
×
r
V_r \in \mathbb R^{n\times r}, \ U_r \in \mathbb R^{m\times r},\ \Sigma_r \in \mathbb R^{r\times r}
Vr∈Rn×r, Ur∈Rm×r, Σr∈Rr×r,则
U
r
=
A
V
r
Σ
r
−
1
U
r
T
U
r
=
(
A
V
r
Σ
r
−
1
)
T
A
V
r
Σ
r
−
1
=
Σ
r
−
1
V
r
T
(
V
Σ
T
Σ
V
T
)
V
r
Σ
r
−
1
=
(
Σ
r
−
1
V
r
T
V
Σ
T
)
(
Σ
r
−
1
V
r
T
V
Σ
T
)
T
U_r=AV_r \Sigma_r^{-1} \\ U_r^T U_r = (AV_r \Sigma_r^{-1})^TAV_r \Sigma_r^{-1}=\Sigma_r^{-1} V_r^T(V\Sigma^T \Sigma V^T) V_r \Sigma_r^{-1}=(\Sigma_r^{-1}V_r^TV\Sigma^T)( \Sigma_r^{-1} V_r^TV\Sigma^T)^T
Ur=AVrΣr−1UrTUr=(AVrΣr−1)TAVrΣr−1=Σr−1VrT(VΣTΣVT)VrΣr−1=(Σr−1VrTVΣT)(Σr−1VrTVΣT)T.
其中 V r T V = I r × n V_r^TV=I_{r\times n} VrTV=Ir×n, Σ r − 1 I r × n = Σ r − 1 \Sigma_r^{-1}I_{r\times n} = \Sigma_r^{-1} Σr−1Ir×n=Σr−1, Σ r − 1 Σ T = I r × r \Sigma_r^{-1}\Sigma^T=I_{r\times r} Σr−1ΣT=Ir×r,所以 U r T U r = I r × r U_r^TU_r=I_{r\times r} UrTUr=Ir×r. 即证明了这 r r r列一定单位正交,式(2)是一定可以构造出来的
-
注意 V V V和 U U U一起构成了四个基本子空间,SVD实现了四个基本子空间的分解,
-
U : [ C ( A ) , N ( A T ) ] U: [C(A), N(A^T)] U:[C(A),N(AT)]: , V : [ C ( A T ) , N ( A ) ] V: [C(A^T), N(A)] V:[C(AT),N(A)]
-
从矩阵分解的角度揭示了矩阵从四个记本子空间中如何生成的,这也是SVD各种秒用的原理之一
几何解释
从线性变换的角度看奇异值分解揭示了矩阵的组成方法, 如下图所示(摘自Prince的Computer vision: models, learning, and inference):
- V T V^T VT 旋转(+反演,反射与否可用行列式判断,下同)
- L L L 伸缩每个维度(假定 L L L对角线都是正值,也即把反演放到旋转中去)
- U U U 旋转(+反演)
- 只改变某部分,对最终矩阵
A
3
A_3
A3 的影响
奇异值 σ \sigma σ
- 个数: m i n ( m , n ) min(m, n) min(m,n)
- 非0奇异值的数量是矩阵的秩
- 最大奇异值与最小奇异值之比反映了可逆性的度量,称为条件数。这个概念和后文的数值秩有联系
- 根据行列式的定义,知,奇异值的乘积为行列式(因为只有 L L L拉伸了矩阵, U U U和 V V V不改变面积)
逆
A
−
1
=
V
L
−
1
U
T
A^{-1}=VL^{-1} U^T
A−1=VL−1UT
从几何角度和数值角度都很直观
极分解Polar Decomposition生成SVD
A
=
Q
S
A=QS
A=QS
其中
Q
Q
Q正交(酉矩阵),
S
S
S半正定矩阵(半正定埃尔米特矩阵)。
所以
A
=
Q
V
Λ
V
ˉ
T
A=QV\Lambda \bar V^T
A=QVΛVˉT
Q
V
=
U
QV=U
QV=U正交,
V
V
V正交,这也即SVD分解
- 关于求 S S S, A ˉ T A = S ˉ T Q ˉ T Q S = S 2 \bar A^TA=\bar S^T \bar Q^T QS=S^2 AˉTA=SˉTQˉTQS=S2,用矩阵开方的知识,正规矩阵都可以开方, S = A ˉ T A = V Λ 2 V ˉ T = V Λ V ˉ T S=\sqrt {\bar A^TA} = \sqrt {V\Lambda ^2 \bar V^T}=V\Lambda \bar V^T S=AˉTA=VΛ2VˉT=VΛVˉT
矩阵范数
这里插入矩阵范数,以便应用部分使用
范数定义
满足这三条的叫做范数
- ∥ A ∥ ⩾ 0 \|A\| \geqslant 0 ∥A∥⩾0
- ∥ c A ∥ = ∣ c ∣ ∥ A ∥ \|cA\|=|c| \|A\| ∥cA∥=∣c∣∥A∥
- ∥ A + B ∥ ⩽ ∥ A ∥ + ∥ B ∥ \|A+B\| \leqslant \|A\|+\|B\| ∥A+B∥⩽∥A∥+∥B∥
如果定义
∥
A
∥
=
(
∑
i
,
j
a
i
,
j
p
)
1
/
p
\|A\|=(\sum_{i,j} a_{i,j}^p)^{1/p}
∥A∥=(∑i,jai,jp)1/p,当
p
<
1
p<1
p<1时,会违背第三条性质,所以不是范数;
p
⩾
1
p\geqslant 1
p⩾1是可以的。
接下来给出3个重要范数,这些范数都可以用奇异值表示。
2范数(谱范数)
∥
A
∥
2
=
max
x
∥
A
x
∥
2
∥
x
∥
2
=
σ
1
\|A\|_2=\max_x\frac{\|Ax\|_2}{\|x\|_2}=\sigma_1
∥A∥2=xmax∥x∥2∥Ax∥2=σ1
即最大的奇异值,不妨假定
∥
x
∥
2
=
1
\|x\|_2=1
∥x∥2=1,证明方法为记
x
=
∑
i
λ
i
v
i
x=\sum_i {\lambda_i}v_i
x=∑iλivi,则
∑
i
λ
i
2
=
1
\sum_i \lambda_i^2=1
∑iλi2=1,根据SVD分解性质有
A
x
=
∑
i
λ
i
σ
i
u
i
Ax=\sum_i \lambda_i\sigma_i u_i
Ax=∑iλiσiui,所以
∥
A
x
∥
2
=
∑
i
λ
i
2
σ
i
2
⩽
σ
1
\|Ax\|_2=\sqrt{\sum_{i} \lambda_i^2\sigma_i^2}\leqslant \sigma_1
∥Ax∥2=∑iλi2σi2⩽σ1
F范数Frobenius Norm
∥ A ∥ F \|A\|_F ∥A∥F是所有元素平方和再开方, ∥ A ∥ F = t r [ A T A ] ∣ = t r [ V Σ T Σ V T ] \|A\|_F=\sqrt {tr[A^TA]|}=\sqrt {tr[V\Sigma^T \Sigma V^T]} ∥A∥F=tr[ATA]∣=tr[VΣTΣVT],根据迹的轮换对称性 t r [ A B ] = t r [ B A ] tr[AB]=tr[BA] tr[AB]=tr[BA],所以 ∥ A ∥ F = t r [ V T V Σ T Σ ] = ∑ i = 1 r σ 2 \|A\|_F=\sqrt {tr[V^TV\Sigma^T \Sigma ]}=\sqrt {\sum_{i=1}^r \sigma^2} ∥A∥F=tr[VTVΣTΣ]=i=1∑rσ2
核范数Nuclear Norm
∥ A ∥ N = ∑ i = 1 r σ i \|A\|_{N}=\sum_{i=1}^r \sigma_i ∥A∥N=i=1∑rσi
这三种范数在正交变换下不变
Q A = ( Q U ) Σ V T QA=(QU)\Sigma V^T QA=(QU)ΣVT,这是 Q A QA QA的SVD分解,如果范数只由 Σ \Sigma Σ决定,那么 Q Q Q不影响范数
向量特例
这三个范数相等,所以经常省略下标,简写为 ∥ x ∥ \|x\| ∥x∥
- 基追踪问题Basis pursuit
m i n ∥ x ∥ p s . t . A x = b p = 1 o r 2 min\|x\|_p\qquad s.t. \quad Ax=b \\ p=1\ or\ 2 min∥x∥ps.t.Ax=bp=1 or 2
这个问题似乎不太好解
应用:SVD矩阵低秩压缩
这里矩阵可以是一幅图像或DataFrame型数据。如果
r
a
n
k
[
A
]
=
r
rank[A]=r
rank[A]=r,那么
A
=
U
Σ
V
T
=
∑
i
=
1
r
σ
i
u
i
v
i
T
A=U\Sigma V^T = \sum_{i=1}^r \sigma_i u_i v_i^T
A=UΣVT=i=1∑rσiuiviT
有点像谱分解,上述
σ
\sigma
σ按绝对值从大到小排列,越大的
σ
\sigma
σ占据了数据的主要成分。
- 如果 r ≪ min ( m , n ) r \ll \min(m,n) r≪min(m,n),则该矩阵可以很好的无损压缩。
- 注意 ∥ u i v i T ∥ F = t r [ v i u i T u i v i T ] = t r [ v i v i T ] = t r [ v i T v i ] = 1 \left\|u_i v_i ^T \right\|_F=tr[v_i u_i^T u_i v_i^T]=tr[v_i v_i^T]=tr[v_i^T v_i]=1 ∥∥uiviT∥∥F=tr[viuiTuiviT]=tr[viviT]=tr[viTvi]=1. 所以每一个谱的数据量主要由 σ \sigma σ衡量。可以只保留最大的 k k k个 σ \sigma σ,对 A A A进行有损压缩得到 A k = U k Σ k V k T A_k=U_k \Sigma_k V_k^T Ak=UkΣkVkT
Eckart-Young定理
如果
A
k
=
U
k
Σ
k
V
k
T
A_k=U_k \Sigma_k V_k^T
Ak=UkΣkVkT,对任意秩
k
k
k矩阵,
∥
A
−
B
∥
⩾
∥
A
−
A
k
∥
\|A-B \| \geqslant\|A-A_k\|
∥A−B∥⩾∥A−Ak∥
A
k
A_k
Ak是秩
k
k
k矩阵中最接近
A
A
A的。该定理不再证明。它说明了SVD矩阵低秩压缩和PCA的最优性。
数值秩Numerical Rank
对于奇异值很小的部分,数值秩采取软处理。定义容忍度
ϵ
∈
(
0
,
1
)
\epsilon \in (0,1)
ϵ∈(0,1),
r
a
n
k
ϵ
[
A
]
=
r
rank_\epsilon [A]=r
rankϵ[A]=r含义是:
σ
k
+
1
⩽
ϵ
σ
1
(
x
)
σ
k
>
ϵ
σ
1
(
x
)
\begin{aligned} \sigma_{k+1} \leqslant \epsilon \sigma_1(x) \\ \sigma_{k} > \epsilon \sigma_1(x) \end{aligned}
σk+1⩽ϵσ1(x)σk>ϵσ1(x)
- r a n k 0 [ x ] = r a n k [ x ] rank_0[x]=rank[x] rank0[x]=rank[x]
- σ k + 1 ( x ) = ∥ A − A k ∥ 2 \sigma_{k+1}(x)=\|A-A_k\|_2 σk+1(x)=∥A−Ak∥2,根据Eckart-Young定理, σ k + 1 \sigma_{k+1} σk+1度量了最优拟合差异的2范数
数值秩类似于模糊数学的概念,如果一个矩阵满秩,但是很多奇异值很小,导致矩阵不容易求逆。在数值秩看来这就是一个低秩矩阵。典型例子:
- 希尔伯特矩阵 H H H,其中 H i j = 1 / ( i + j − 1 ) H_{ij}=1/(i+j-1) Hij=1/(i+j−1),该矩阵满秩,但数值秩很低。对于1000维希尔伯特矩阵, r a n k [ H 1000 ] = 1000 , r a n k 1 0 − 15 [ H 1000 ] = 28 rank[H_{1000}]=1000, rank_{10^{-15}}[H_{1000}]=28 rank[H1000]=1000,rank10−15[H1000]=28
- 范德蒙矩阵。注意范德蒙矩阵也是多项式线性回归中扩展后的数据矩阵
其他判断低数值秩矩阵的方法包括Sylvester等式,不再展开细讲。
应用:主成分分析Principal Components Analysis(PCA)
假设存在数据集
X
=
[
x
1
,
⋯
,
x
n
]
∈
R
m
×
n
X=[x_1,\cdots,x_n] \in \mathbb R^{m\times n}
X=[x1,⋯,xn]∈Rm×n,
m
m
m是特征数,
n
n
n是样本数,希望通过线性变换
W
∈
R
k
×
m
W \in \mathbb R^{k\times m}
W∈Rk×m把特征降低到
k
k
k维,左乘得到降维结果
W
X
WX
WX,并希望降维再升维度的数据能通过线性变换尽可能代表原始数据,距离最小。
因为降维后的矩阵秩不超过
k
k
k,如果再升维回去,秩仍然不超过
k
k
k,根据Eckart-Young定理,
X
k
X_k
Xk是最接近
X
X
X的矩阵。
其中
X
k
=
U
k
Σ
k
V
k
T
X_k=U_k \Sigma_k V_k^T
Xk=UkΣkVkT通过SVD分解得到。
什么样的线性变换能在降维和升维后得到
X
k
X_k
Xk?
- 降维:实际上可以取 U k T U_k^T UkT左乘, U k T X U_k^T X UkTX把数据集降维到 R k × n \mathbb R^{k\times n} Rk×n,即 u i T x u_i^Tx uiTx是降维后第 i i i维的结果
- 升维:取矩阵 U k U_k Uk作为升维矩阵, U k ( U k T X ) = U k ( U k T U ) Σ V T = U k I k × m Σ V T = U k [ Σ k , 0 k × ( n − k ) ] V T = U k Σ k V k T = X k U_k(U_k^TX)=U_k(U_k^TU)\Sigma V^T=U_kI_{k\times m}\Sigma V^T=U_k [\Sigma_{k}, \textbf{0}_{k\times (n-k)}]V^T=U_k\Sigma_kV_k^T=X_k Uk(UkTX)=Uk(UkTU)ΣVT=UkIk×mΣVT=Uk[Σk,0k×(n−k)]VT=UkΣkVkT=Xk.
所以这样的变换是成立的,即对
X
X
T
XX^T
XXT进行特征分解,取前
k
k
k大特征值对应的特征向量
u
1
,
⋯
,
u
r
u_1, \cdots, u_r
u1,⋯,ur作为降维的投影向量。
如果变换之前,先对
X
X
X减去各个特征均值,得到
X
^
=
X
−
μ
m
1
1
×
n
\hat X=X-\mu_{m} \textbf 1_{1\times n}
X^=X−μm11×n,即
X
^
1
n
=
0
\hat X \textbf 1_{n}=0
X^1n=0,再进行降维,这就是PCA,其中
1
1
×
n
\textbf 1_{1\times n}
11×n和
1
n
\textbf 1_{n}
1n表示全1行向量和列向量。注意
X
^
X
^
T
/
(
m
−
1
)
\hat X \hat X^T/(m-1)
X^X^T/(m−1)是特征间的协方差矩阵。
应用:多维缩放Multiple Dimensional Scaling
对于 n n n个坐标未知的样本点,给定两两之间的距离平方,构成距离矩阵 D ∈ R n × n D \in \mathbb R^{n\times n} D∈Rn×n,其中 D i j = ∥ x i − x j ∥ 2 D_{ij}=\|x_i-x_j\|^2 Dij=∥xi−xj∥2. 求出这 n n n个样本点的坐标。往往我们要求坐标空间的维度 m < n m<n m<n.
- 这个问题在降维中有实际意义。参考低维嵌入和流形学习。考虑高维空间中的点,这些点可能分布在低维流形manifold上,如图所示(摘自周志华《机器学习》),通过计算相邻样本点的距离,构成图,计算任意两点间的测地线距离(最短路径长度)当作这两点的距离,就获得了距离矩阵,但是不知道低维样本点坐标。
已知
D
i
j
=
x
i
T
x
i
+
x
j
T
x
j
−
2
x
i
T
x
j
D_{ij}=x_i^Tx_i +x_j^Tx_j-2x_i^Tx_j
Dij=xiTxi+xjTxj−2xiTxj,思路是先求出
X
T
X
X^TX
XTX. 假设
∑
i
=
1
n
x
i
=
0
\sum_{i=1}^n x_i=0
∑i=1nxi=0,则
∑
i
=
1
n
∑
j
=
1
n
D
i
j
=
2
n
⋅
t
r
[
X
T
X
]
\sum_{i=1}^n\sum_{j=1}^n D_{ij}=2n\cdot tr[X^TX]
i=1∑nj=1∑nDij=2n⋅tr[XTX]
可以求出
t
r
[
X
T
X
]
tr[X^TX]
tr[XTX]. 列求和则得到
∑
j
=
1
n
D
i
j
=
n
x
i
T
x
i
+
t
r
[
X
T
X
]
\sum_{j=1}^n D_{ij}=nx_i^Tx_i + tr[X^TX]
j=1∑nDij=nxiTxi+tr[XTX]
又可以求出任何
x
i
T
x
i
x_i^Tx_i
xiTxi,从而任何
x
i
T
x
j
x_i^Tx_j
xiTxj都可以求出。写成矩阵结果为
X
i
T
X
j
=
−
1
2
(
D
−
1
n
D
1
n
×
n
−
1
n
1
n
×
n
D
+
1
n
2
t
r
[
D
T
D
]
)
X_i^TX_j=-\frac{1}{2}(D-\frac{1}{n}D \textbf 1_{n\times n}-\frac{1}{n}\textbf 1_{n\times n}D +\frac{1}{n^2}tr[D^TD])
XiTXj=−21(D−n1D1n×n−n11n×nD+n21tr[DTD])
其中
1
n
×
n
\textbf 1_{n\times n}
1n×n表示全1矩阵。
如果
D
D
D是满足三角不等式的合理距离矩阵,那么
X
T
X
X^TX
XTX是半正定对称阵。
从而特征值分解(也即SVD分解)得到
X
T
X
=
V
Λ
V
T
X^TX=V\Lambda V^T
XTX=VΛVT,如果我们希望降到
m
m
m维,根据前文Eckart-Young定理,取特征值最大的前
m
m
m个特征值和特征向量,得到
X
^
=
Λ
m
1
/
2
V
m
T
\hat X=\Lambda^{1/2}_mV^T_m
X^=Λm1/2VmT.
(注意找低维近似的时候,没有办法做矩阵开方)
应用:主方向/最小值方向问题
主方向问题:
arg max
b
∥
A
b
∥
s
.
t
.
∥
b
∥
=
1
\argmax_b \|Ab\| \qquad s.t. \quad \|b\|=1
bargmax∥Ab∥s.t.∥b∥=1
最小值方向问题则是求
arg
min
∥
A
b
∥
\arg\min \|Ab\|
argmin∥Ab∥。
- 作为工具,该问题和最小二乘的“地位”类似。可以用于求解几何模型中针孔摄像机的外在参数,双视图匹配等问题。下一节正交Procrustes问题也在这些类似问题上求参数时会用到。详见Computer vision: models, learning, and inference的14-16章。
- 根据限制,
b
b
b必须在单位圆上,问题变为找一个方向,对应变换后的椭圆长轴。
根据上述SVD分解,待求方向经过 V T V^T VT变换后,为坐标系主轴方向 X = ( 1 , 0 , . . . , 0 ) X=(1, 0, ..., 0) X=(1,0,...,0)
把它变回去,该方向为 ( V T ) − 1 X = V X (V^T)^{-1}X=VX (VT)−1X=VX,即 V V V的第一列。 - 同理,最小值方向为 V V V的最后一列
应用:正交Procrustes问题
Procrustes是希腊语Προκρούστης.
寻找正交矩阵
Q
Q
Q,使
Q
^
=
arg min
Q
∥
Q
A
−
B
∥
F
\hat Q =\argmin_Q \left\|Q A-B \right\|_F
Q^=Qargmin∥QA−B∥F
首先
Q
^
=
arg min
Q
t
r
[
(
Q
A
−
B
)
T
(
Q
A
−
B
)
]
=
arg max
Q
t
r
[
Q
T
B
A
T
]
\begin{aligned} \hat Q &=\argmin_Q tr\left[(Q A-B)^T(Q A-B)\right] \\ &=\argmax_Q tr \left[Q^TBA^T\right] \end{aligned}
Q^=Qargmintr[(QA−B)T(QA−B)]=Qargmaxtr[QTBAT]
计算SVD分解
B
A
T
=
U
Σ
V
T
BA^T=U\Sigma V^T
BAT=UΣVT,则
Q
^
=
arg max
Q
t
r
[
Q
T
U
Σ
V
T
]
=
arg max
Q
t
r
[
V
T
Q
T
U
Σ
]
\begin{aligned} \hat Q &=\argmax_Q tr\left[Q^T U\Sigma V^T\right]\\ &= \argmax_Q tr\left[V^TQ^T U\Sigma \right] \end{aligned}
Q^=Qargmaxtr[QTUΣVT]=Qargmaxtr[VTQTUΣ]
注意到
t
r
[
V
T
Q
T
U
Σ
]
=
t
r
[
Z
Σ
]
=
∑
i
=
1
I
z
i
i
l
i
i
tr \left[V^TQ^T U\Sigma \right] = tr \left[Z\Sigma \right] = \sum_{i=1}^I z_{ii}l_{ii}
tr[VTQTUΣ]=tr[ZΣ]=i=1∑Iziilii
其中
Z
=
V
T
Q
T
U
Z=V^TQ^TU
Z=VTQTU,
Z
Z
Z是三个正交矩阵乘积,所以也正交,对角线上每个数值小于等于1,所以选择
Z
=
I
Z=I
Z=I来最大化上述目标,全局解为
Q
^
=
U
V
T
\hat Q = UV^T
Q^=UVT
- 本问题的一个特殊形式是给定方阵
B
B
B,寻找最接近的正交矩阵
Q
Q
Q,即最优化
Q ^ = arg min Q ∥ Q − B ∥ F \hat Q=\argmin_Q \|Q - B\|_F Q^=Qargmin∥Q−B∥F
用上述方法,得到解为 Q ^ = U V T \hat Q=UV^T Q^=UVT,其中 U Σ V T = B U\Sigma V^T=B UΣVT=B
应用:压缩感知Compressed Sensing
压缩感知的目的是通过部分信息恢复全部信息。思想是把全信息空间转换成稀疏表达空间(例如空域转频域、把读者喜好的书籍条目转为读者喜好的书籍分类等)。这种表达理应是稀疏的,所以当只观察到部分数据时,可以以此去寻找稀疏表达空间中尽量稀疏的解。
压缩感知可以用于矩阵补全,进一步可用于填补数据缺失值。对于矩阵
A
A
A,部分数据缺失,压缩感知的处理方式是填补缺失值并使
r
a
n
k
[
A
]
rank[A]
rank[A]最小。这是NP难问题,注意到
r
a
n
k
[
A
]
rank[A]
rank[A]在集合
{
X
∈
R
m
×
n
:
∥
X
∥
F
2
⩽
1
}
\{X\in \mathbb R^{m\times n}:\|X\|^2_F \leqslant 1\}
{X∈Rm×n:∥X∥F2⩽1}的凸包是
X
X
X的核范数,所以可通过最小化核范数的形式近似求解原问题。可用半正定规划求解。秩和核范数的关系有点像是0范数和1范数的关系。
更多讨论参考周志华《机器学习》11.6节。
求逆技巧
左右逆
对于 A ∈ R m × n A \in \mathbb R ^{m\times n} A∈Rm×n
- 左逆: ( A T A ) − 1 A T ∈ R n × m (A^TA)^{-1}A^T \in \mathbb R^{n\times m} (ATA)−1AT∈Rn×m,这需要 A T A A^TA ATA可逆,也即 r a n k [ A ] = n rank[A]=n rank[A]=n
- 右逆: A T ( A A T ) − 1 ∈ R n × m A^T(AA^T)^{-1} \in \mathbb R^{n\times m} AT(AAT)−1∈Rn×m,需要 r a n k [ A T ] = r a n k [ A ] = m rank[A^T]=rank[A]=m rank[AT]=rank[A]=m
伪逆
考虑行空间和列空间,从行空间中拿元素 x x x, A x Ax Ax一定在列空间中。而且 x x x和 A x Ax Ax能构成一一对应,两个空间维度相同。注意该定理和SVD代数解释中的原理联系。证明:
- 对于 x 1 , x 2 ∈ C ( A T ) x_1,x_2 \in C(A^T) x1,x2∈C(AT), A x 1 ≠ A x 2 Ax_1\neq Ax_2 Ax1=Ax2,否则 A ( x 1 − x 2 ) = 0 A(x_1-x_2)=0 A(x1−x2)=0, x 1 − x 2 x_1-x_2 x1−x2不可能在 N ( A ) N(A) N(A)当中,所以单射!
- 另外只取
x
∈
C
(
A
T
)
,
A
x
x\in C(A^T), Ax
x∈C(AT),Ax就能张成列空间,证明:
对于 x ′ ∈ R n x'\in \mathbb R^{n} x′∈Rn, A x ′ Ax' Ax′是表示各种列向量的组合情况,必然能张成列空间。另一方面,必有分解 x ′ = x + x N x'=x+x_N x′=x+xN,其中 x ∈ C ( A T ) , x N ∈ N ( A ) x\in C(A^T),x_N \in N(A) x∈C(AT),xN∈N(A), A x ′ = A x + A x N = A x Ax'=Ax+Ax_N=Ax Ax′=Ax+AxN=Ax,也即只需要 A x Ax Ax就能张成列空间。满射! - 所以双射,即一一对应成立。
此时限制在这两空间上的逆就是伪逆,即对于 A x = y Ax=y Ax=y,伪逆 A + A^+ A+使 x = A + y = ( A + A ) x x=A^{+}y=(A^+A)x x=A+y=(A+A)x
应用:无需 A T A A^TA ATA可逆的最小二乘问题
先来介绍下SVD伪逆
SVD伪逆
伪逆有很多求法,一种求法是SVD分解:
A
=
U
Σ
V
T
A=U\Sigma V^T
A=UΣVT,
U
,
V
U,V
U,V可逆,
Σ
+
\Sigma^+
Σ+对角线元素为
1
/
σ
1
,
1
/
σ
2
,
⋯
,
1
/
σ
r
,
0
,
⋯
,
0
1/\sigma_1, 1/\sigma_2, \cdots, 1/\sigma_r, 0, \cdots,0
1/σ1,1/σ2,⋯,1/σr,0,⋯,0,其余元素为0,注意
Σ
∈
R
m
×
n
,
Σ
+
∈
R
n
×
m
\Sigma \in \mathbb R^{m\times n},\Sigma^+ \in \mathbb R^{n\times m}
Σ∈Rm×n,Σ+∈Rn×m,
A
+
=
V
Σ
+
U
T
A^+=V\Sigma^+ U^T
A+=VΣ+UT
-
A
V
=
U
Σ
=
[
u
1
σ
1
,
⋯
,
u
r
σ
r
,
0
,
⋯
,
0
]
AV=U\Sigma=[u_1\sigma_1, \cdots, u_r\sigma_r, 0,\cdots,0]
AV=UΣ=[u1σ1,⋯,urσr,0,⋯,0],
A
A
A把
v
r
+
1
,
⋯
,
v
n
v_{r+1}, \cdots, v_n
vr+1,⋯,vn变换成0,这是必然的,因为
v
r
+
1
,
⋯
,
v
n
v_{r+1},\cdots,v_n
vr+1,⋯,vn是在
N
(
A
)
N(A)
N(A)中。有趣的是
A
+
U
=
V
Σ
+
=
[
v
1
/
σ
1
,
⋯
,
u
r
/
σ
r
,
0
,
⋯
,
0
]
A^+U=V\Sigma^+=[v_1/\sigma_1, \cdots, u_r/\sigma_r, 0, \cdots, 0]
A+U=VΣ+=[v1/σ1,⋯,ur/σr,0,⋯,0],这说明
A
+
A^+
A+把
u
r
+
1
,
⋯
,
u
m
u_{r+1},\cdots,u_m
ur+1,⋯,um变换为0,进一步说明:
- 性质1: A + A^+ A+把 N ( A T ) N(A^T) N(AT)变换成0
对于
A
x
=
b
Ax=b
Ax=b,有最优解
x
^
=
(
A
T
A
)
−
1
A
T
b
\hat x = (A^TA)^{-1}A^Tb
x^=(ATA)−1ATb,这需要
A
T
A
A^TA
ATA可逆。用伪逆可以消除该限制,而且思考过程非常漂亮!
重新审视最小二乘问题,仍用投影的思路,但是换一个思考方式。记
b
=
b
C
+
b
N
b=b_C+b_N
b=bC+bN,其中
b
C
∈
C
(
A
)
,
b
N
∈
N
(
A
T
)
b_C \in C(A), b_N \in N(A^T)
bC∈C(A),bN∈N(AT),根据投影,
A
x
=
b
Ax=b
Ax=b与
A
x
=
b
C
Ax=b_C
Ax=bC同解。
如果取
x
=
A
+
b
C
x=A^+b_C
x=A+bC,则
A
x
=
A
A
+
b
C
Ax=AA^+b_C
Ax=AA+bC,因为
b
C
∈
C
(
A
)
b_C \in C(A)
bC∈C(A),所以
A
A
+
b
C
=
b
C
AA^+b_C=b_C
AA+bC=bC,说明
x
=
A
+
b
C
x=A^+b_C
x=A+bC是
A
x
=
b
C
Ax=b_C
Ax=bC的解,也即
A
x
=
b
Ax=b
Ax=b的解。
根据性质1,
A
+
b
C
=
A
+
b
A^+b_C=A^+b
A+bC=A+b,所以
x
=
A
+
b
=
V
Σ
+
U
b
x=A^+b=V\Sigma^+Ub
x=A+b=VΣ+Ub是
A
x
=
b
Ax=b
Ax=b的一个最小二乘可行解
- 讨论:
- 如果
A
T
A
A^TA
ATA可逆,
r
a
n
k
[
A
]
=
n
rank[A]=n
rank[A]=n,意味着
m
⩾
n
m\geqslant n
m⩾n,则
x
=
A
+
b
=
(
A
T
A
)
−
1
A
T
b
x=A^+b=(A^TA)^{-1}A^Tb
x=A+b=(ATA)−1ATb
此时 A + = ( A T A ) − 1 A T A^+=(A^TA)^{-1}A^T A+=(ATA)−1AT;进一步如果 A A A可逆,则 A + = A − 1 A^+=A^{-1} A+=A−1 - 如果 A T A A^TA ATA不可逆,则 x = A + b x=A^+b x=A+b是可行解中的一个,是唯一的 C ( A T ) C(A^T) C(AT)中的解
- 如果
A
T
A
A^TA
ATA可逆,
r
a
n
k
[
A
]
=
n
rank[A]=n
rank[A]=n,意味着
m
⩾
n
m\geqslant n
m⩾n,则
x
=
A
+
b
=
(
A
T
A
)
−
1
A
T
b
x=A^+b=(A^TA)^{-1}A^Tb
x=A+b=(ATA)−1ATb
注意当 σ \sigma σ从0扰动到非0时,伪逆变化非常大。一种处理方法是把 σ \sigma σ小于某个阈值的都当作0. 另一种方法是引入L2正则化
L2正则化
对于
A
x
=
b
Ax=b
Ax=b问题,其中
A
∈
R
m
×
n
A\in \mathbb R^{m\times n}
A∈Rm×n,如果
r
a
n
k
[
A
]
<
n
rank[A]< n
rank[A]<n或者
A
A
A接近不列满秩,
A
T
A
A^TA
ATA行列式很小,求逆后结果爆炸,则可以改求
min
x
∥
A
x
−
b
∥
2
2
+
δ
2
∥
x
∥
2
2
(
δ
>
0
)
\min_x \|Ax-b\|^2_2+\delta^2\|x\|_2^2 \qquad (\delta > 0)
xmin∥Ax−b∥22+δ2∥x∥22(δ>0)
上式对
x
x
x进行正则化regularization,这是统计学家的观点。但从线性代数的角度可以写成成矩阵形式,即最小二乘
A
∗
x
=
[
A
δ
I
n
×
n
]
x
=
[
b
0
n
×
1
]
=
b
∗
A^*x=\begin{bmatrix} A \\ \delta I_{n\times n} \end{bmatrix} x= \begin{bmatrix} b \\ {0}_{n\times 1} \end{bmatrix}=b^*
A∗x=[AδIn×n]x=[b0n×1]=b∗
这里
A
∗
T
A
∗
A^{*T}A^*
A∗TA∗的可逆性很好
- r a n k [ A ∗ ] = n rank[A^*]=n rank[A∗]=n,从而 A ∗ T A ∗ A^{*T}A^* A∗TA∗可逆
- 另一种更推荐的解释方法: A ∗ T A ∗ = A T A + δ 2 I n × n A^{*T}A^*=A^TA+\delta^2 I_{n\times n} A∗TA∗=ATA+δ2In×n,这里 A T A A^TA ATA半正定, δ \delta δ越大, A ∗ T A ∗ A^{*T}A^* A∗TA∗越正定,越远离不满秩。
从而 x = ( A ∗ T A ∗ ) − 1 A ∗ T b ∗ = ( A ∗ T A ∗ ) − 1 A T b = ( A T A + δ 2 I n × n ) − 1 A T b x=(A^{*T}A^*)^{-1}A^{*T}b^*=(A^{*T}A^*)^{-1}A^{T}b=(A^TA+\delta^2 I_{n\times n})^{-1}A^{T}b x=(A∗TA∗)−1A∗Tb∗=(A∗TA∗)−1ATb=(ATA+δ2In×n)−1ATb
逆关系
Sherman-Morrison恒等式
对于
U
,
V
∈
R
n
×
k
U,V \in \mathbb R^{n\times k}
U,V∈Rn×k,
r
a
n
k
[
U
]
=
r
a
n
k
[
V
]
=
k
rank[U]=rank[V]=k
rank[U]=rank[V]=k
(
I
n
−
U
V
T
)
−
1
=
I
n
+
U
(
I
k
−
V
T
U
)
−
1
V
T
(I_n-UV^T)^{-1}=I_n+U(I_k -V^TU)^{-1}V^T
(In−UVT)−1=In+U(Ik−VTU)−1VT
只需证明
I
n
=
I
n
−
U
V
T
+
(
I
n
−
U
V
T
)
U
(
I
k
−
V
T
U
)
−
1
V
T
=
I
n
−
U
V
T
+
U
(
I
k
−
V
T
U
)
−
1
V
T
−
U
V
T
U
(
I
k
−
V
T
U
)
−
1
V
T
=
I
n
−
U
V
T
+
U
(
I
k
−
V
T
U
)
(
I
k
−
V
T
U
)
−
1
V
T
=
I
n
\begin{aligned} I_n&=I_n-UV^T+(I_n-UV^T)U(I_k -V^TU)^{-1}V^T \\ &=I_n-UV^T+U(I_k -V^TU)^{-1}V^T-UV^TU(I_k -V^TU)^{-1}V^T \\ &=I_n-UV^T+U(I_k-V^TU)(I_k -V^TU)^{-1}V^T =I_n\\ \end{aligned}
In=In−UVT+(In−UVT)U(Ik−VTU)−1VT=In−UVT+U(Ik−VTU)−1VT−UVTU(Ik−VTU)−1VT=In−UVT+U(Ik−VTU)(Ik−VTU)−1VT=In
- 秩1矩阵
u
,
v
u,v
u,v
( I − u v T ) − 1 = I + u v T 1 − v T u (I-uv^T)^{-1}=I+ \frac{uv^T}{1-v^Tu} (I−uvT)−1=I+1−vTuuvT
Sherman–Morrison–Woodbury恒等式
该式也叫矩阵求逆引理、Woodbury矩阵恒等式。
令Sherman-Morrison恒等式中
U
=
A
X
,
V
T
=
C
−
1
Y
U=AX, V^T=C^{-1}Y
U=AX,VT=C−1Y,其中
A
∈
R
n
×
n
,
X
∈
R
n
×
k
,
Y
∈
R
k
×
n
,
C
∈
R
k
×
k
A \in \mathbb R ^{n\times n}, X \in \mathbb R ^{n\times k}, Y \in \mathbb R^{k\times n}, C \in \mathbb R ^{k\times k}
A∈Rn×n,X∈Rn×k,Y∈Rk×n,C∈Rk×k,推导几步,易得
(
A
−
1
+
X
C
−
1
Y
)
−
1
=
A
−
A
X
(
C
+
Y
A
X
)
−
1
Y
A
(A^{-1}+XC^{-1}Y)^{-1}=A-AX(C+YAX)^{-1}YA
(A−1+XC−1Y)−1=A−AX(C+YAX)−1YA
- 该式的一个用处是降低求逆维度,可以把左侧对 n n n维求逆改为右侧对 k k k维求逆,如果 n > k n>k n>k,则有助于求逆。
如果
C
=
I
k
C=I_k
C=Ik,易得
(
A
−
1
+
X
Y
)
−
1
=
A
−
A
X
(
I
k
+
Y
A
X
)
−
1
Y
A
(A^{-1}+XY)^{-1}=A-AX(I_k+YAX)^{-1}YA
(A−1+XY)−1=A−AX(Ik+YAX)−1YA
- 该公式揭示了如果对矩阵 A − 1 A^{-1} A−1有秩 k k k的修改 X Y XY XY,例如 X Y XY XY是修改的 C R CR CR分解,那么其逆矩阵也会有秩 k k k的修改过程。
如果
Y
=
X
T
Y=X^T
Y=XT,易得
(
A
−
1
+
X
C
−
1
X
T
)
−
1
=
A
−
A
X
(
X
T
A
X
+
C
)
−
1
X
T
A
(A^{-1} + X C^{-1} X^T)^{-1}=A-AX(X^TAX+C)^{-1}X^TA
(A−1+XC−1XT)−1=A−AX(XTAX+C)−1XTA
应用:因子分析Factor Analysis的概率密度函数
对于
n
n
n维高斯分布,如果维度过大(例如把整个图像像素建模成高斯分布),协方差参数太多,达到
O
(
n
2
)
\mathcal O(n^2)
O(n2)量级,可以用因子分析模型进行有损压缩,把参数降低到线性。该模型形式化为
P
(
x
∣
h
)
=
N
(
x
∣
μ
+
Φ
h
,
Σ
)
P
(
h
)
=
N
(
h
∣
0
,
I
)
\begin{aligned} P( x| h) &= \mathcal N( x| \mu + \Phi h, \Sigma) \\ P( h) &= \mathcal N( h| 0, I) \end{aligned}
P(x∣h)P(h)=N(x∣μ+Φh,Σ)=N(h∣0,I)
其中观测数据
x
x
x的维度为
n
n
n,隐变量
h
h
h的维度为
k
(
n
>
k
)
k\quad (n>k)
k(n>k),参数包括:
Φ
∈
R
n
×
k
\Phi \in \mathbb R^{n\times k}
Φ∈Rn×k,
Σ
\Sigma
Σ 是对角阵.
- x = μ + Φ h + ϵ x=\mu+\Phi h + \epsilon x=μ+Φh+ϵ, x x x可由隐变量 h h h解释
- 当 Σ \Sigma Σ是单位阵的常数倍时,该模型对球形协方差建模,称为概率主成分分析probabilistic principal component analysis或概率PCA.
- 因子分析模型和PCA类似,可以用于人脸形状建模等问题。
可以求出 P ( x ) = N ( x ∣ μ , Φ Φ T + Σ ) P( x) = \mathcal N( x| \mu , \Phi\Phi^T +\Sigma) P(x)=N(x∣μ,ΦΦT+Σ) ,仍然是高斯分布,且协方差参数进行了压缩。如何求解呢?
求解方法1:高斯线性模型
注意到这是高斯线性模型(linear Gaussian model)!直接套结论。结果为
Φ
Φ
T
+
Σ
\Phi\Phi^T +\Sigma
ΦΦT+Σ. 可以参考PRML 2.3.3节,用了高斯分布逆辨识,比较tricky.
这个过程中需要用舒尔补恒等式Schur complement identity对分块矩阵求逆,而舒尔补恒等式则和Sherman-Morrison-Woodbury恒等式有很大联系。
求解方法2:配平方法:边缘分布推导
直接暴力求解,需要耐心
P
(
x
)
=
∫
N
(
x
∣
μ
+
Φ
h
,
Σ
)
N
(
h
∣
0
,
I
)
d
h
=
C
1
∫
e
x
p
{
−
1
2
(
x
−
μ
−
Φ
h
)
T
Σ
−
1
(
x
−
μ
−
Φ
h
)
}
e
x
p
{
−
1
2
h
T
h
}
d
h
=
C
2
∫
e
x
p
{
−
1
2
[
(
x
−
μ
)
T
Σ
−
1
(
x
−
μ
)
−
2
(
x
−
μ
)
T
Σ
−
1
Φ
h
+
h
T
(
Φ
T
Σ
−
1
Φ
+
I
)
h
]
}
d
h
=
C
3
exp
{
−
1
2
[
(
x
−
μ
)
T
[
Σ
−
1
−
Σ
−
1
Φ
(
Φ
T
Σ
−
1
Φ
+
I
)
−
1
Φ
T
Σ
−
1
]
(
x
−
μ
)
]
}
×
∫
e
x
p
{
−
1
2
[
h
−
(
Φ
T
Σ
−
1
Φ
+
I
)
−
1
Φ
T
Σ
−
1
(
x
−
μ
)
]
T
(
Φ
T
Σ
−
1
Φ
+
I
)
[
h
−
(
Φ
T
Σ
−
1
Φ
+
I
)
−
1
Φ
T
Σ
−
1
(
x
−
μ
)
]
}
d
h
\begin{aligned} P( x) &= \int \mathcal N( x| \mu + \Phi h, \Sigma) \mathcal N ( h|0, I) dh\\ =& C_1 \int exp \left\{-\frac{1}{2}( x- \mu - \Phi h)^T \Sigma^{-1}( x- \mu - \Phi h) \right\} exp \left\{-\frac{1}{2} h^T h\right\} d h \\ =& C_2 \int exp\left\{-\frac{1}{2} [( x- \mu)^T \Sigma^{-1}( x- \mu ) -2( x - \mu)^T \Sigma^{-1} \Phi h + h^T ( \Phi^T \Sigma^{-1} \Phi + I) h]\right\} d h \\ =& C_3 \exp \left \{ -\frac{1}{2}[( x- \mu)^T[ \Sigma^{-1} - \Sigma^{-1} \Phi ( \Phi^T \Sigma^{-1} \Phi + I)^{-1} \Phi^T \Sigma^{-1}]( x- \mu )] \right \} \\& \times \int exp \left \{-\frac{1}{2} [ h - ( \Phi^T \Sigma ^{-1} \Phi + I)^{-1} \Phi^T \Sigma^{-1}( x - \mu)]^T ( \Phi^T \Sigma^{-1} \Phi + I)[ h - ( \Phi^T \Sigma ^{-1} \Phi + I)^{-1} \Phi^T \Sigma^{-1}( x - \mu)] \right \} d h \end{aligned}
P(x)====∫N(x∣μ+Φh,Σ)N(h∣0,I)dhC1∫exp{−21(x−μ−Φh)TΣ−1(x−μ−Φh)}exp{−21hTh}dhC2∫exp{−21[(x−μ)TΣ−1(x−μ)−2(x−μ)TΣ−1Φh+hT(ΦTΣ−1Φ+I)h]}dhC3exp{−21[(x−μ)T[Σ−1−Σ−1Φ(ΦTΣ−1Φ+I)−1ΦTΣ−1](x−μ)]}×∫exp{−21[h−(ΦTΣ−1Φ+I)−1ΦTΣ−1(x−μ)]T(ΦTΣ−1Φ+I)[h−(ΦTΣ−1Φ+I)−1ΦTΣ−1(x−μ)]}dh
其中
C
1
,
C
2
,
C
3
C_1,C_2,C_3
C1,C2,C3表示
x
,
h
x,h
x,h无关的系数。注意到把
x
x
x有关而
h
h
h无关的项提到积分外,积分内是关于
h
h
h的一个高斯分布。积分后得到
C
4
exp
{
−
1
2
[
(
x
−
μ
)
T
[
Σ
−
1
−
Σ
−
1
Φ
(
Φ
T
Σ
−
1
Φ
+
I
)
−
1
Φ
T
Σ
−
1
]
(
x
−
μ
)
]
}
\begin{aligned} C_4 \exp \left \{ -\frac{1}{2}[( x- \mu)^T[ \Sigma^{-1} - \Sigma^{-1} \Phi ( \Phi^T \Sigma^{-1} \Phi + I)^{-1} \Phi^T \Sigma^{-1}]( x- \mu )] \right \} \end{aligned}
C4exp{−21[(x−μ)T[Σ−1−Σ−1Φ(ΦTΣ−1Φ+I)−1ΦTΣ−1](x−μ)]}
这是一个高斯分布,均值为
μ
\mu
μ,方差为
[
Σ
−
1
−
Σ
−
1
Φ
(
Φ
T
Σ
−
1
Φ
+
I
)
−
1
Φ
T
Σ
−
1
]
−
1
[ \Sigma^{-1} - \Sigma^{-1} \Phi ( \Phi^T \Sigma^{-1} \Phi + I)^{-1} \Phi^T \Sigma^{-1}]^{-1}
[Σ−1−Σ−1Φ(ΦTΣ−1Φ+I)−1ΦTΣ−1]−1
这里套用Sherman–Morrison–Woodbury恒等式,上式化为
Φ
Φ
T
+
Σ
\Phi\Phi^T +\Sigma
ΦΦT+Σ
应用:卡尔曼滤波
卡尔曼滤波在跟踪问题中作为时序模型常见。卡尔曼滤波中含有两个高斯线性模型,和因子分析模型类似,Sherman–Morrison–Woodbury恒等式可以化简卡尔曼滤波的表达式。更多内容可以参考Prince的Computer vision: models, learning, and inference第19章或Bishop的Pattern Recognition and Machine Learning第13.3节
分解方法Factorizations汇总
A = C R 列/行向量无关 A = L U elimination,下上三角分解 A = Q R Gram-Schmitt, Q 正交 A = V Λ V − 1 一 般 特 征 分 解 ( Λ 为 J o r d a n 标 准 型 ) S = Q Λ Q T S 对称, 特征分解, Q 正交 A = U Σ V T SVD A = Q S Polar Decomposition \begin{aligned} A=&CR &\qquad \text{ 列/行向量无关}\\ A=& LU&\qquad \text{ elimination,下上三角分解}\\ A=&QR &\qquad \text{Gram-Schmitt, $Q$正交} \\ A=&V\Lambda V^{-1} & 一般特征分解(\Lambda为Jordan标准型)\\ S=&Q\Lambda Q^T & \qquad \text{$S$对称, 特征分解,$Q$正交} \\ A= &U\Sigma V^T &\qquad \text{SVD} \\ A = &QS &\qquad \text{Polar Decomposition} \end{aligned} A=A=A=A=S=A=A=CRLUQRVΛV−1QΛQTUΣVTQS 列/行向量无关 elimination,下上三角分解Gram-Schmitt, Q正交一般特征分解(Λ为Jordan标准型)S对称, 特征分解,Q正交SVDPolar Decomposition
注意每种分解前后的自由度相等,可以自行验证。
对于矩阵
A
∈
R
m
×
n
,
r
a
n
k
[
A
]
=
r
A \in \mathbb R^ {m \times n}, rank[A]=r
A∈Rm×n,rank[A]=r,其自由度为
(
m
+
n
−
r
)
r
(m+n-r)r
(m+n−r)r,可以用SVD分解或者从几何角度找子空间中的向量组来证明。
A=LU的解释与高斯消元法
L
U
=
[
l
1
,
l
2
,
⋯
]
[
u
1
u
2
⋮
]
=
l
1
u
1
+
l
2
u
2
+
⋯
LU=[l_1, l_2, \cdots]\begin{bmatrix} u_1\\ u_2 \\ \vdots \end{bmatrix} \\ =l_1u_1+l_2u_2+\cdots
LU=[l1,l2,⋯]⎣⎢⎡u1u2⋮⎦⎥⎤=l1u1+l2u2+⋯
L
L
L是单位下三角矩阵(对角线为1)。
l
i
u
i
l_iu_i
liui维度与
A
A
A同,可以看作是
l
1
u
1
l_1u_1
l1u1剥离了
A
A
A的第一行第一列,……这样,
l
i
u
i
l_i u_i
liui只需面对已经消除了
i
−
1
i-1
i−1行和
i
−
1
i-1
i−1列的矩阵,所以
l
i
l_i
li的前
i
−
1
i-1
i−1个元素都是0,
u
i
u_i
ui同理。这样即实现了
L
U
LU
LU分解.
- 当A的所有顺序主子式都不为0时,主元不为0,矩阵
A
A
A可以进行LU分解,且是唯一分解(摘自百度百科)
上图也可以从高斯消元的角度看待, U U U是消元结果, L L L是初等矩阵的拼接
矩阵微分
一般情况
矩阵微分不难理解,但是很多资料疏于整理,导致这套技术不统一也难成体系,学习者不易掌握一套通用方法。
矩阵求导有很多种形式:
- 标量对矩阵(向量)求(偏)导:这是机器学习最常用的形式,即Loss对参数求导
- 矩阵(向量)对标量求导
- 向量对向量求导
推荐两个文章:
- 矩阵求导术(上):这篇博客讲授Jacobian辨识,并将其作为一种较为通用的标量对矩阵求导的工具。可以处理点乘、矩阵乘法等多种运算。
- Matrix Calculus:师兄的文章,除了标量对矩阵求导外,对各类矩阵求导形式进行了精炼的推导和总结。适合自学和速查
杂知识:
- 标量对向量求导后,再求二阶导的Jacobian矩阵即为Hessian矩阵
- 对于, f ( A ) = ln d e t [ A ] f(A)=\ln det[A] f(A)=lndet[A], ∇ f = ( A − 1 ) T \nabla f=(A^{-1})^T ∇f=(A−1)T,可以自行验证
应用:神经网络反向传播向量化
利用该方法可以推导出向量化的反向传播,见神经网络反向传播向量化(CS231n A1 Q4)
应用:动量Momentum优化二次型
对于函数
f
(
x
)
=
1
2
x
T
S
x
f(x)=\frac{1}{2}x^TSx
f(x)=21xTSx,
S
S
S是正定矩阵。
采用梯度下降的线搜索会以zigzag的形式逼近极小值。(如果每次线搜索都搜到最好的位置,那么zigzag每两个相邻迭代的方向是垂直的)
更高效的方法是采用momentum优化器,迭代方式为
x
k
+
1
=
x
k
−
α
z
k
z
k
+
1
=
∇
f
k
+
1
+
β
z
k
x_{k+1}=x_k - \alpha z_k \\ z_{k+1}=\nabla f_{k+1} + \beta z_{k}
xk+1=xk−αzkzk+1=∇fk+1+βzk
用向量化写法可形式化为
[
I
0
−
S
I
]
[
x
z
]
k
+
1
=
[
I
−
α
I
0
β
I
]
[
x
z
]
k
\begin{bmatrix} I & 0 \\ -S & I \end{bmatrix} \begin{bmatrix} x\\ z \end{bmatrix}_{k+1} = \begin{bmatrix} I & -\alpha I \\ 0 & \beta I \end{bmatrix} \begin{bmatrix} x\\ z \end{bmatrix}_{k}
[I−S0I][xz]k+1=[I0−αIβI][xz]k
记
S
S
S特征值
λ
\lambda
λ,对应特征向量
q
q
q。注意
S
S
S是正定矩阵,只考虑一个特征分量,如果
x
k
=
c
k
q
,
z
k
=
d
k
q
x_k=c_kq,z_k=d_k q
xk=ckq,zk=dkq,那么
S
x
k
=
c
k
λ
q
Sx_k=c_k\lambda q
Sxk=ckλq,易证
[
1
0
−
λ
1
]
[
c
k
+
1
d
k
+
1
]
=
[
1
−
α
0
β
]
[
c
k
d
k
]
[
c
k
+
1
d
k
+
1
]
=
[
1
−
α
λ
β
−
λ
α
]
[
c
k
d
k
]
=
R
[
c
k
d
k
]
\begin{aligned} \begin{bmatrix} 1 & 0 \\ -\lambda & 1 \end{bmatrix} \begin{bmatrix} c_{k+1}\\ d_{k+1} \end{bmatrix}& = \begin{bmatrix} 1 & -\alpha \\ 0 & \beta \end{bmatrix} \begin{bmatrix} c_k\\ d_k \end{bmatrix} \\ \begin{bmatrix} c_{k+1}\\ d_{k+1} \end{bmatrix} &= \begin{bmatrix} 1 & -\alpha \\ \lambda & \beta -\lambda \alpha \end{bmatrix} \begin{bmatrix} c_k\\ d_k \end{bmatrix}\\ &=R\begin{bmatrix} c_k\\ d_k \end{bmatrix} \end{aligned}
[1−λ01][ck+1dk+1][ck+1dk+1]=[10−αβ][ckdk]=[1λ−αβ−λα][ckdk]=R[ckdk]
S
S
S是正定的,记最小和最大特征值分别为
m
,
M
(
>
0
)
m, M(>0)
m,M(>0),
M
m
\frac{M}{m}
mM即条件数,为1时,
S
S
S是单位阵;当条件数很大时,则不好用梯度下降优化。
这里希望找到最好的
α
,
β
\alpha, \beta
α,β,使对任意
λ
\lambda
λ(满足
m
⩽
λ
⩽
M
m\leqslant \lambda \leqslant M
m⩽λ⩽M),
c
k
,
d
k
c_k, d_k
ck,dk减少最快,即
R
R
R的所有特征值尽可能小。
可以求得最优解为
α
o
p
t
=
(
2
M
+
m
)
2
β
o
p
t
=
(
M
−
m
M
+
m
)
2
\alpha_{opt}=\left(\frac{2}{\sqrt M + \sqrt m}\right)^2 \\ \beta_{opt}=\left(\frac{\sqrt M - \sqrt m}{\sqrt M + \sqrt m}\right)^2
αopt=(M+m2)2βopt=(M+mM−m)2
此时
R的特征值
<
(
1
−
b
1
+
b
)
2
\text{R的特征值} <\left (\frac{1-\sqrt b}{1+\sqrt b} \right)^2
R的特征值<(1+b1−b)2
其中
b
=
m
M
b=\frac{m}{M}
b=Mm。
(该式正确性存疑,望有凸优化的大佬告知)
逆矩阵
对于矩阵函数
A
(
t
)
A(t)
A(t),易得
(
A
+
d
A
)
−
1
−
A
−
1
=
(
A
+
d
A
)
−
1
[
A
−
(
A
+
d
A
)
]
A
−
1
=
−
A
−
1
d
A
A
−
1
\begin{aligned} (A+dA)^{-1}-A^{-1} &=(A+dA)^{-1}[A-(A+dA)]A^{-1} \\ &=-A^{-1}dAA^{-1} \end{aligned}
(A+dA)−1−A−1=(A+dA)−1[A−(A+dA)]A−1=−A−1dAA−1
所以
d
A
−
1
d
t
=
−
A
−
1
d
A
d
t
A
−
1
\frac{dA^{-1}}{dt}=-A^{-1}\frac{dA}{dt}A^{-1}
dtdA−1=−A−1dtdAA−1
特征值
对于矩阵
A
(
t
)
A(t)
A(t),转置后特征值不变,但特征向量则有可能会变,所以对于同一特征值,有列特征向量和行特征向量:
A
(
t
)
x
(
t
)
=
λ
(
t
)
x
(
t
)
y
(
t
)
T
A
(
t
)
=
λ
(
t
)
y
(
t
)
T
A(t)x(t)=\lambda(t)x(t) \\ y(t)^TA(t)=\lambda(t) y(t)^T
A(t)x(t)=λ(t)x(t)y(t)TA(t)=λ(t)y(t)T
对相同特征值对应的行列特征向量,引入限制
y
(
t
)
T
x
(
t
)
=
1
y(t)^Tx(t)=1
y(t)Tx(t)=1
则
y
(
t
)
T
A
(
t
)
x
(
t
)
=
λ
(
t
)
y(t)^TA(t)x(t) =\lambda(t)
y(t)TA(t)x(t)=λ(t)
所以
d
λ
d
t
=
d
y
T
d
t
A
x
+
y
T
d
A
d
t
x
+
y
T
A
d
x
d
t
=
y
T
d
A
d
t
x
+
d
y
T
d
t
λ
x
+
λ
y
T
d
x
d
t
=
y
T
d
A
d
t
x
+
λ
d
(
y
T
x
)
d
t
=
y
(
t
)
T
d
A
(
t
)
d
t
x
(
t
)
\begin{aligned} \frac{d\lambda}{d t} &= \frac{dy^T}{dt}Ax+y^T\frac{dA}{dt}x+y^TA\frac{dx}{dt} \\ &=y^T\frac{dA}{dt}x + \frac{dy^T}{dt}\lambda x + \lambda y^T \frac{dx}{dt} \\ &= y^T\frac{dA}{dt}x + \lambda \frac{d(y^Tx)}{dt} \\ &= y(t)^T\frac{dA(t)}{dt}x(t) \end{aligned}
dtdλ=dtdyTAx+yTdtdAx+yTAdtdx=yTdtdAx+dtdyTλx+λyTdtdx=yTdtdAx+λdtd(yTx)=y(t)TdtdA(t)x(t)
该式揭示了特征值的变化如何对应于矩阵的变化
奇异值
对于
A
A
A的同一奇异值
σ
\sigma
σ对应的左右奇异向量
u
,
v
u,v
u,v,有
u
(
t
)
T
A
(
t
)
v
(
t
)
=
u
(
t
)
T
[
σ
(
t
)
u
(
t
)
]
=
σ
(
t
)
u(t)^TA(t)v(t)=u(t)^T[\sigma(t)u(t)]=\sigma(t)
u(t)TA(t)v(t)=u(t)T[σ(t)u(t)]=σ(t),所以
d
σ
d
t
=
d
u
T
d
t
A
v
+
u
T
d
A
d
t
v
+
u
T
A
d
v
d
t
=
σ
d
u
T
d
t
u
+
u
T
d
A
d
t
v
+
σ
v
T
d
v
d
t
\begin{aligned} \frac{d\sigma}{dt} &= \frac{du^T}{dt}Av + u^T\frac{dA}{dt}v + u^TA\frac{dv}{dt} \\ &= \sigma\frac{du^T}{dt}u + u^T\frac{dA}{dt}v + \sigma v^T\frac{dv}{dt} \end{aligned}
dtdσ=dtduTAv+uTdtdAv+uTAdtdv=σdtduTu+uTdtdAv+σvTdtdv
注意到
d
(
u
T
u
)
d
t
=
0
=
d
u
T
d
t
u
+
u
T
d
u
d
t
\frac{d(u^Tu)}{dt}=0=\frac{du^T}{dt} u+u^T\frac{du}{dt}
dtd(uTu)=0=dtduTu+uTdtdu,所以
d
u
T
d
t
u
=
0
\frac{du^T}{dt}u=0
dtduTu=0,同理
d
v
T
d
t
=
0
\frac{dv^T}{dt}=0
dtdvT=0,所以
d
σ
d
t
=
u
T
(
t
)
d
A
(
t
)
d
t
v
(
t
)
\begin{aligned} \frac{d\sigma}{dt}& =u^T(t)\frac{dA(t)}{dt}v(t) \end{aligned}
dtdσ=uT(t)dtdA(t)v(t)
其中
u
,
v
u,v
u,v是对应
σ
\sigma
σ的左右奇异向量。
行列式
d
∣
A
∣
=
∣
A
∣
t
r
(
A
−
1
d
A
)
d|A|=|A|tr(A^{-1}dA)
d∣A∣=∣A∣tr(A−1dA)
从而
d
ln
∣
A
∣
=
t
r
(
A
−
1
d
A
)
d\ln |A|=tr(A^{-1}dA)
dln∣A∣=tr(A−1dA)
特殊矩阵
循环矩阵Circulant Matrix
行或列逐级轮换的矩阵
C
=
[
c
0
c
n
−
1
⋯
c
2
c
1
c
1
c
0
c
n
−
1
c
2
⋮
c
1
c
0
⋱
⋮
c
n
−
2
⋱
⋱
c
n
−
1
c
n
−
1
c
n
−
2
⋯
c
1
c
0
]
C=\begin{bmatrix} c_0 & c_{n-1} & \cdots & c_2 &c_1 \\ c_1 & c_0 & c_{n-1} & &c_2 \\ \vdots & c_1 & c_0 & \ddots&\vdots \\ c_{n-2} & &\ddots &\ddots & c_{n-1} \\ c_{n-1}&c_{n-2}&\cdots& c_1 & c_0 \end{bmatrix}
C=⎣⎢⎢⎢⎢⎢⎡c0c1⋮cn−2cn−1cn−1c0c1cn−2⋯cn−1c0⋱⋯c2⋱⋱c1c1c2⋮cn−1c0⎦⎥⎥⎥⎥⎥⎤
C
C
C是一个正规矩阵。
一种特殊的循环矩阵是
P
4
=
[
1
1
1
1
]
P_4=\begin{bmatrix} &1&&\\ &&1&\\ &&&1\\ 1&&& \end{bmatrix}
P4=⎣⎢⎢⎡1111⎦⎥⎥⎤
这是一个正交矩阵。该矩阵作用的效果是轮换
P
4
[
x
0
x
1
x
2
x
3
]
=
[
x
1
x
2
x
3
x
0
]
P_4\begin{bmatrix} x_0 \\x_1\\x_2\\x_3 \end{bmatrix}=\begin{bmatrix} x_1\\x_2\\x_3\\x_0 \end{bmatrix}
P4⎣⎢⎢⎡x0x1x2x3⎦⎥⎥⎤=⎣⎢⎢⎡x1x2x3x0⎦⎥⎥⎤
-
如果 C 1 , C 2 C_1,C_2 C1,C2是两个循环矩阵,那么 C 1 C 2 C_1C_2 C1C2仍然是循环矩阵。
注意 P 2 P^2 P2相当于 P P P对自己的列轮换,所以
P 4 2 = [ 1 1 1 1 ] P^2_4=\begin{bmatrix} &&1&\\ &&&1\\ 1&&&\\ &1&& \end{bmatrix} P42=⎣⎢⎢⎡1111⎦⎥⎥⎤ -
任意循环矩阵 C ∈ R n × n C\in \mathbb R^{n\times n} C∈Rn×n,可以分解为 C = ∑ i = 0 n − 1 c i P n i C=\sum_{i=0}^{n-1}c_i P^i_n C=∑i=0n−1ciPni
-
循环矩阵是一种特殊的常对角矩阵Toeplitz matrix
应用:离散傅里叶变换DFT
正交矩阵一节提到了傅里叶变换,这里我们从循环矩阵的角度来考察。
思考
P
n
P_n
Pn的特征值和特征向量,因为
P
n
P_n
Pn的效果是轮换,所以可以想到特征值等分复平面单位圆,一个特征向量的每个元素依特征值等差取在复平面单位圆上。形式化:记
ω
=
e
2
π
/
n
\omega = e^{2\pi/n}
ω=e2π/n,则特征值为
0
,
ω
,
ω
2
,
⋯
,
ω
n
−
1
0, \omega, \omega^2, \cdots, \omega^{n-1}
0,ω,ω2,⋯,ωn−1,对应特征向量的矩阵为
V
=
1
n
[
1
1
1
⋯
1
1
ω
ω
2
⋯
ω
n
−
1
1
ω
2
ω
4
⋯
ω
2
(
n
−
1
)
⋮
1
ω
n
−
1
ω
2
(
n
−
1
)
⋯
ω
(
n
−
1
)
(
n
−
1
)
]
V=\frac{1}{\sqrt n}\begin{bmatrix} 1&1&1&\cdots&1 \\ 1 & \omega & \omega^2 & \cdots& \omega^{n-1} \\ 1 & \omega^2 & \omega^4 & \cdots &\omega^{2(n-1)} \\ &&&\vdots& \\ 1&\omega^{n-1}&\omega^{2(n-1)}&\cdots& \omega^{(n-1)(n-1)} \end{bmatrix}
V=n1⎣⎢⎢⎢⎢⎢⎡11111ωω2ωn−11ω2ω4ω2(n−1)⋯⋯⋯⋮⋯1ωn−1ω2(n−1)ω(n−1)(n−1)⎦⎥⎥⎥⎥⎥⎤
这恰是离散傅里叶变换的基。
另外,因为任意循环矩阵
C
=
∑
i
=
0
n
−
1
c
i
P
n
i
C=\sum_{i=0}^{n-1}c_i P^i_n
C=∑i=0n−1ciPni,所以特征向量和
P
P
P都是相同的。可以证明特征值则满足
d
i
a
g
[
Λ
]
=
V
[
c
0
c
1
⋮
c
n
−
1
]
diag[\Lambda] = V\begin{bmatrix}c_0 \\ c_1\\\vdots\\ c_{n-1} \end{bmatrix}
diag[Λ]=V⎣⎢⎢⎢⎡c0c1⋮cn−1⎦⎥⎥⎥⎤
克罗内克积Kronecker Product(直积)
A
⊗
B
=
[
a
11
B
⋯
a
1
n
B
⋮
⋱
⋮
a
m
1
B
⋯
a
m
n
B
]
A\otimes B=\begin{bmatrix}a_{11}B &\cdots& a_{1n}B \\ \vdots & \ddots & \vdots \\ a_{m1}B & \cdots & a_{mn}B \end{bmatrix}
A⊗B=⎣⎢⎡a11B⋮am1B⋯⋱⋯a1nB⋮amnB⎦⎥⎤
其中
A
∈
R
m
×
n
,
B
∈
R
p
×
q
,
A
⊗
B
∈
R
m
p
×
n
q
A \in \mathbb R^{m\times n}, B\in \mathbb R^{p\times q}, A\otimes B \in \mathbb R^{mp \times nq}
A∈Rm×n,B∈Rp×q,A⊗B∈Rmp×nq
有性质
(
A
⊗
B
)
T
=
A
T
⊗
B
T
(A\otimes B)^T=A^T\otimes B^T
(A⊗B)T=AT⊗BT
克罗内克积在求二阶导中常见,例如矩阵求导术(下)
参考文献:
[1] MIT18.06 Linear Algebra
[2] MIT18.065 Matrix Methods in Data Analysis, Signal Processing, and Machine Learning
[3] 3Blue1Brown 线性代数的本质
[4] Prince S J D. Computer vision: models, learning, and inference. Cambridge University Press, 2012.
[5] Goodfellow I, Bengio Y, Courville A. Deep learning. MIT press, 2016.
[6] Gonzales, Rafael C., and Richard E. Woods. Digital image processing. Fourth Edition 2018.
[7] Bishop C . Pattern Recognition and Machine LearningStat Sci. 2006.
[8] 周志华. 机器学习. 清华大学出版社. 2016.
[9] 同济大学数学系. 线性代数第五版. 高等教育出版社. 2007.
[10] 程云鹏等. 矩阵论. 西北工业大学出版社. 2006.
[11] 居余马等. 线性代数第2版. 清华大学出版社. 2002.
[12] 维基百科,如Woodbury matrix identity等