特征值相关概念
基础概念
矩阵A\mathbf AA的特征值和特征向量:
Ax=λx, x≠0\mathbf {Ax}=\lambda\mathbf x,\ \mathbf x\neq\mathbf 0
Ax=λx, x=0
特征值的求法:
λ: det(λI−A)=0\lambda:\ \det(\lambda\mathbf I-\mathbf A)=0
λ: det(λI−A)=0
特征向量的求法:
x: (λI−A)x=0\mathbf x:\ (\lambda\mathbf I-\mathbf A)\mathbf x=\mathbf 0
x: (λI−A)x=0
特征向量构成特征子空间
矩阵的特征值谱(即特征值组成的集合)记作λ(A)\lambda(\mathbf A)λ(A)
基本结论
定理 矩阵奇异当且仅当矩阵特征值含有000
定理 λ(A)=λ(AT)\lambda(\mathbf A)=\lambda(\mathbf A^T)λ(A)=λ(AT)
证明 转置不改变行列式。det(λI−A)=det(λI−AT)\det(\lambda\mathbf I-\mathbf A)=\det(\lambda\mathbf I-\mathbf A^T)det(λI−A)=det(λI−AT)
定理 对角阵或上(下)三角阵的特征值为对角元
定理 分块对角阵或分块上(下)三角阵的特征值谱是各个对角块特征值谱的并集
定义 对矩阵A,B\mathbf A,\mathbf BA,B,若存在可逆阵X\mathbf XX使得A=XBX−1\mathbf A=\mathbf{XB}\mathbf X^{-1}A=XBX−1,则称A,B\mathbf A,\mathbf BA,B相似
定理 相似矩阵特征值谱相同
证明 对A\mathbf AA的特征对(λ,v)(\lambda,\mathbf v)(λ,v)
Bv=X−1AXv=X−1(λXv)=λv\mathbf{Bv}=\mathbf{X}^{-1}\mathbf{AXv}=\mathbf{X}^{-1}(\lambda\mathbf{Xv})=\lambda\mathbf v
Bv=X−1AXv=X−1(λXv)=λv
定理(矩阵运算结果的特征值)
设λj\lambda_jλj为A\mathbf AA的特征值,则
λ(cA)={cλj}\lambda(c\mathbf A)=\{c\lambda_j\}λ(cA)={cλj}
λ(A+cI)={λj+c}\lambda(\mathbf A+c\mathbf I)=\{\lambda_j+c\}λ(A+cI)={λj+c}
λ(P(A))={P(λj)}\lambda(P(\mathbf A))=\{P(\lambda_j)\}λ(P(A))={P(λj)},其中PPP为多项式
λ(A−1)={λj−1}\lambda(\mathbf A^{-1})=\{\lambda_j^{-1}\}λ(A−1)={λj−1}
特征值的重数、对角化、亏损阵
定义(特征值的重数)
对A\mathbf AA的特征值λ0\lambda_0λ0,此特征值在特征方程det(λI−A)=0\det(\lambda\mathbf I-\mathbf A)=0det(λI−A)=0中的重数定义为代数重数,特征子空间的维数dimN(λ0I−A)\dim\mathcal N(\lambda_0\mathbf I-\mathbf A)dimN(λ0I−A)定义为几何重数
定理(代数重数和几何重数的关系) 特征值的代数重数大于或等于其几何重数
证明
设矩阵A\mathbf AA特征值λj\lambda_jλj的代数重数为njn_jnj,几何重数为kjk_jkj,则可以在其特征子空间N(λjI−A)\mathcal N(\lambda_j\mathbf I-\mathbf A)N(λjI−A)中找到一组基:x1,⋯ ,xkj\mathbf x_1,\cdots,\mathbf x_{k_j}x1,⋯,xkj
将此基扩充成Rn\mathbb R^nRn的一组基:x1,⋯ ,xn\mathbf x_1,\cdots,\mathbf x_nx1,⋯,xn
令X=[x1⋯xn]\mathbf X=\begin{bmatrix}\mathbf x_1&\cdots&\mathbf x_n\end{bmatrix}X=[x1⋯xn],则有
AX=A[x1⋯xn]=[λjx1⋯λjxkjAxkj+1⋯Axn]=X[λjIkjBC]\mathbf {AX}=\mathbf A\begin{bmatrix}\mathbf x_1&\cdots&\mathbf x_n\end{bmatrix}=\begin{bmatrix}\lambda_j\mathbf x_1&\cdots&\lambda_j\mathbf x_{k_j}&\mathbf{Ax}_{k_j+1}&\cdots&\mathbf{Ax}_n\end{bmatrix}=\mathbf X\begin{bmatrix}\lambda_j\mathbf I_{k_j}&\mathbf B\\&\mathbf C\end{bmatrix}
AX=A[x1⋯xn]=[λjx1⋯λjxkjAxkj+1⋯Axn]=X[λjIkjBC]
于是
det(λI−A)=det(X−1(λI−A)X)=det[(λ−λj)IkjBλIn−kj−C]=(λ−λj)kjdet(λIn−kj−C)\det(\lambda\mathbf I-\mathbf A)=\det(\mathbf X^{-1}(\lambda\mathbf I-\mathbf A)\mathbf X)=\det\begin{bmatrix}(\lambda-\lambda_j)\mathbf I_{k_j}&\mathbf B\\&\lambda\mathbf I_{n-k_j}-\mathbf C\end{bmatrix}=(\lambda-\lambda_j)^{k_j}\det(\lambda\mathbf I_{n-k_j}-\mathbf C)
det(λI−A)=det(X−1(λI−A)X)=det[(λ−λj)IkjBλIn−kj−C]=(λ−λj)kjdet(λIn−kj−C)
即(λ−λj)kj(\lambda-\lambda_j)^{k_j}(λ−λj)kj整除A\mathbf AA的特征多项式PA(λ)P_{\mathbf A}(\lambda)PA(λ),故nj≥kjn_j\geq k_jnj≥kj
定理 不同特征值的特征向量线性无关
证明 利用特征值的定义显然
定义 如果所有特征值的几何重数都等于代数重数,则矩阵称为非亏损阵,否则为亏损阵
非亏损阵的特征向量可以组成全空间的一组基
定理 非亏损阵和可对角化矩阵是同样的概念,即,矩阵A\mathbf AA非亏损等价于存在X\mathbf XX使X−1AX=Λ\mathbf X^{-1}\mathbf {AX}=\mathbf{\Lambda}X−1AX=Λ
证明 取一组特征向量,得到的矩阵就是X\mathbf XX
实对称矩阵的对角化
定理 实对称矩阵的特征值是实数
证明 对A\mathbf AA的特征对(λ,x)(\lambda, \mathbf x)(λ,x),有
Ax=λx\mathbf {Ax}=\lambda\mathbf x
Ax=λx
两边同取共轭,有
Axˉ=λˉxˉ\mathbf A\bar{\mathbf x}=\bar\lambda \bar{\mathbf x}
Axˉ=λˉxˉ
对第一个式子左乘xˉT\bar{\mathbf x}^TxˉT,有
xˉTAx=λxˉTx\bar{\mathbf x}^T\mathbf {Ax}=\lambda\bar{\mathbf x}^T\mathbf x
xˉTAx=λxˉTx
对第二个式子左乘xT\mathbf x^TxT,有
xTAxˉ=λˉxTxˉ\mathbf x^T\mathbf A\bar{\mathbf x}=\bar\lambda \mathbf x^T\bar{\mathbf x}
xTAxˉ=λˉxTxˉ
由于A\mathbf AA对称
xˉTAx=xTAxˉ\bar{\mathbf x}^T\mathbf {Ax}=\mathbf x^T\mathbf A\bar{\mathbf x}
xˉTAx=xTAxˉ
于是
λxˉTx=λˉxˉTx\lambda\bar{\mathbf x}^T\mathbf x=\bar\lambda\bar{\mathbf x}^T\mathbf x
λxˉTx=λˉxˉTx
又由于x≠0\mathbf x\neq\mathbf 0x=0
xˉTx>0\bar{\mathbf x}^T\mathbf x>0
xˉTx>0
于是
λ=λˉ\lambda=\bar\lambda
λ=λˉ
定理(实对称矩阵的谱分解/正交对角化) 设A\mathbf AA实对称,则存在正交矩阵Q\mathbf QQ,使A=QΛQT\mathbf A=\mathbf{Q\Lambda Q}^TA=QΛQT
证明 对阶数nnn归纳,n=1n=1n=1时显然,假设n−1n-1n−1阶成立
对nnn阶的A\mathbf AA,可以找到一组实特征对(λ1,q1)(\lambda_1,\mathbf q_1)(λ1,q1),不妨设∣∣q1∣∣=1||\mathbf q_1||=1∣∣q1∣∣=1,将q1\mathbf q_1q1扩充成一组实标准正交基q1,⋯ ,qn\mathbf q_1,\cdots,\mathbf q_nq1,⋯,qn,设
Q1=[q1⋯qn]\mathbf Q_1=\begin{bmatrix}\mathbf q_1&\cdots&\mathbf q_n\end{bmatrix}
Q1=[q1⋯qn]
则
AQ1=[Aq1⋯Aqn]=[λ1q1Aq2⋯Aqn]=Q1[λ1bTC]\mathbf{AQ}_1=\begin{bmatrix}\mathbf A\mathbf q_1&\cdots&\mathbf A\mathbf q_n\end{bmatrix}=\begin{bmatrix}\lambda_1\mathbf q_1&\mathbf A\mathbf q_2&\cdots&\mathbf A\mathbf q_n\end{bmatrix}=\mathbf Q_1\begin{bmatrix}\lambda_1&\mathbf b^T\\&\mathbf C\end{bmatrix}
AQ1=[Aq1⋯Aqn]=[λ1q1Aq2⋯Aqn]=Q1[λ1bTC]
于是
Q1TAQ1=[λ1bTC]\mathbf Q_1^T\mathbf{AQ}_1=\begin{bmatrix}\lambda_1&\mathbf b^T\\&\mathbf C\end{bmatrix}
Q1TAQ1=[λ1bTC]
左端是实对称矩阵,于是
[λ1bTC]=[λ1C]\begin{bmatrix}\lambda_1&\mathbf b^T\\&\mathbf C\end{bmatrix}=\begin{bmatrix}\lambda_1&\\&\mathbf C\end{bmatrix}
[λ1bTC]=[λ1C]
且C\mathbf CC也是实对称矩阵,根据归纳假设知存在
C=Q~2Λ2Q~2T\mathbf C=\tilde{\mathbf Q}_2\mathbf{\Lambda}_2\tilde{\mathbf Q}_2^T
C=Q~2Λ2Q~2T
令
Q2=[1Q~2],Λ=[λ1Λ2]\mathbf Q_2=\begin{bmatrix}1&\\&\tilde{\mathbf Q}_2\end{bmatrix},\quad\mathbf{\Lambda}=\begin{bmatrix}\lambda_1&\\&\mathbf\Lambda_2\end{bmatrix}
Q2=[1Q~2],Λ=[λ1Λ2]
则
Q1TAQ1=[λ1C]=Q2ΛQ2T\mathbf Q_1^T\mathbf{AQ}_1=\begin{bmatrix}\lambda_1&\\&\mathbf C\end{bmatrix}=\mathbf Q_2\mathbf\Lambda\mathbf Q_2^T
Q1TAQ1=[λ1C]=Q2ΛQ2T
于是
A=Q1Q2ΛQ2TQ1T\mathbf A=\mathbf Q_1\mathbf Q_2\mathbf{\Lambda Q}_2^T\mathbf Q_1^T
A=Q1Q2ΛQ2TQ1T
令
Q=Q1Q2\mathbf Q=\mathbf Q_1\mathbf Q_2
Q=Q1Q2
证毕
定理 实对称矩阵属于不同特征值的特征向量正交
证明 可以直接根据正交对角化得到此结论,也可以直接证明:
设A\mathbf AA有两个特征对(λ1,x1),(λ2,x2)(\lambda_1,\mathbf x_1),(\lambda_2,\mathbf x_2)(λ1,x1),(λ2,x2),满足λ1≠λ2\lambda_1\neq\lambda_2λ1=λ2,则
λ1x1Tx2=x2TAx1=x1TAx2=λ2x1Tx2\lambda_1\mathbf x_1^T\mathbf x_2=\mathbf x_2^T\mathbf{Ax}_1=\mathbf x_1^T\mathbf{Ax}_2=\lambda_2\mathbf x_1^T\mathbf x_2
λ1x1Tx2=x2TAx1=x1TAx2=λ2x1Tx2
于是
x1Tx2=0\mathbf x_1^T\mathbf x_2=0
x1Tx2=0
Jordan分解
定理(Jordan分解) 任意方阵A\mathbf AA可以进行如下分解:
A=XJX−1\mathbf A=\mathbf{XJX}^{-1}
A=XJX−1
其中
J=[J1⋱Jp]\mathbf J=\begin{bmatrix}\mathbf J_1&&\\&\ddots&\\&&\mathbf J_p\end{bmatrix}
J=J1⋱Jp
每个Js\mathbf J_sJs称为Jordan块,形如
Js=[λs1λs⋱⋱1λs]\mathbf J_s=\begin{bmatrix}\lambda_s&1\\&\lambda_s&\ddots\\&&\ddots&1\\&&&\lambda_s\end{bmatrix}
Js=λs1λs⋱⋱1λs
每个特征值可能有多个Jordan块,但其Jordan块总数为其几何重数,总阶数为其代数重数
Jordan分解是谱分解A=XΛX−1\mathbf A=\mathbf{X\Lambda X}^{-1}A=XΛX−1对亏损阵的推广,得到的J\mathbf JJ称为Jordan标准形
Jordan分解的理论基础和求法与数值分析课程关系不大, 在这里不展开,详见Jordan分解
特征值的分布
特征值的分布对于迭代法解线性方程组的收敛性有很强的影响,例如
ρ(B)=maxj∣λj(B)∣\rho(\mathbf B)=\max_j|\lambda_j(\mathbf B)|
ρ(B)=jmax∣λj(B)∣
cond(A)2=λmax(ATA)λmin(ATA){\rm cond}(\mathbf A)_2=\sqrt{\frac{\lambda_{\rm max}(\mathbf A^T\mathbf A)}{\lambda_{\rm min}(\mathbf A^T\mathbf A)}}
cond(A)2=λmin(ATA)λmax(ATA)
定义(Gerschgorin圆盘) 设A∈Cn×n\mathbf A\in\mathbb C^{n\times n}A∈Cn×n,rk=∑j=1,j≠kn∣akj∣r_k=\sum\limits_{j=1,j\neq k}^n|a_{kj}|rk=j=1,j=k∑n∣akj∣,在复平面上以akka_{kk}akk为圆心,rkr_krk为半径的圆称为A\mathbf AA的Gerschgorin圆盘
定理(圆盘定理) A\mathbf AA的特征值在nnn个圆盘的并集上。圆盘并集的每个连通部分内包含的特征值数量等于其包含的圆盘数量(计算重数)
证明
对特征对(λ,x)(\lambda,\mathbf x)(λ,x),设xkx_kxk是x\mathbf xx的最大分量,由于
Ax=λx\mathbf {Ax}=\lambda\mathbf x
Ax=λx
知
λxk=∑jakjxj\lambda x_k=\sum\limits_j a_{kj}x_j
λxk=j∑akjxj
于是
(λ−akk)xk=∑j≠kakjxj(\lambda-a_{kk})x_k=\sum\limits_{j\neq k}a_{kj}x_j
(λ−akk)xk=j=k∑akjxj
于是
∣λ−akk∣≤∑j≠k∣akj∣|\lambda-a_{kk}|\leq\sum\limits_{j\neq k}|a_{kj}|
∣λ−akk∣≤j=k∑∣akj∣
证毕
注 圆盘定理的一个直接推论是按行严格对角占优矩阵可逆,因为这样的矩阵的所有圆盘都不包含原点;且由于转置不改变特征值,这里的按行也可以是按列
推论 对对角元均正的实矩阵:
严格对角占优的情况下,特征值实部均正
严格对角占优且对称的情况下,对称正定
圆盘定理可以用来估计矩阵的特征值分布,可以使用的技巧有:
注意到对A\mathbf AA和AT\mathbf A^TAT使用圆盘定理可以得到不同的圆盘(行圆盘和列圆盘不同),可以选择二者中范围较小的来使用
实矩阵如果有虚特征值,一定共轭成对出现,因此如果实矩阵一个圆盘内只有一个特征值,此特征值一定是实的
幂法和反幂法
幂法与规格化幂法
定义 矩阵模最大的特征值称为主特征值,其特征向量称为主特征向量
注 主特征值不一定唯一
算法(幂法)
随机取非零向量v0\mathbf v_0v0,作迭代法
vk+1=Avk\mathbf v_{k+1}=\mathbf {Av}_k
vk+1=Avk
可以发现相邻的vk\mathbf v_kvk逐渐呈倍数关系,其倍数为一个主特征值,而vk\mathbf v_kvk则趋近于主特征值的特征向量
定理 若矩阵有唯一的主特征值,则幂法收敛于主特征值的特征向量
证明 将v0\mathbf v_0v0对谱分解或Jordan分解中矩阵X\mathbf XX构成的基(即特征向量基或Jordan向量基)分解即可
注 幂法中,vk≈λ1kx1\mathbf v_k\approx \lambda_1^k\mathbf x_1vk≈λ1kx1,此时迭代次数较多时很容易上溢或下溢,因此实际计算中对每个迭代向量规格化(即使模最大的分量模为111,与标准化类似但计算量比标准化小,会得到长度略大于111的一个向量),以防止溢出
注 决定幂法收敛速度的主要因素是模次大特征值(“次主特征值”)与主特征值的模之比
当幂法加入规格化后,规格化向量经一步迭代后得到的最大模长元素即可作为主特征值的近似值(因为此元素在迭代前为111),减少了一次除法,这也是规格化相比标准化的优势
加速幂法的收敛
原点位移技术
考虑到B=A−sI\mathbf B=\mathbf A-s\mathbf IB=A−sI的特征值是λ(A)−s\lambda(\mathbf A)-sλ(A)−s,如果
∣λ2−sλ1−s∣<∣λ2λ1∣\left|\frac{\lambda_2-s}{\lambda_1-s}\right|<\left|\frac{\lambda_2}{\lambda_1}\right|
λ1−sλ2−s<λ1λ2
则对B\mathbf BB应用幂法可以加速收敛
Rayleigh商加速
定义 实矩阵A\mathbf AA对实向量x\mathbf xx的Rayleigh商定义为
R(x)=xTAxxTxR(\mathbf x)=\frac{\mathbf x^T\mathbf{Ax}}{\mathbf x^T\mathbf x}
R(x)=xTxxTAx
定理 A\mathbf AA对任何向量的瑞利商不大于其最大特征值,不小于其最小特征值
定理 幂法中的规格化向量uk\mathbf u_kuk满足
R(uk)=λ1+O((λ2λ1)2k)R(\mathbf u_k)=\lambda_1+O\left(\left(\frac{\lambda_2}{\lambda_1}\right)^{2k}\right)
R(uk)=λ1+O((λ1λ2)2k)
注 对于幂法中取uk\mathbf u_kuk最大模分量的操作,其得到的结果为λ1+O((λ2λ1)k)\lambda_1+O\left(\left(\frac{\lambda_2}{\lambda_1}\right)^{k}\right)λ1+O((λ1λ2)k),收敛不如瑞利商法快
反幂法
由于A−1\mathbf A^{-1}A−1的特征值是A\mathbf AA的倒数,对A−1\mathbf A^{-1}A−1应用幂法即可得到模最小特征值
此时每步迭代可以通过求解线性方程组来计算A−1\mathbf A^{-1}A−1乘向量(避免求逆),但计算量仍然比正幂法大
适用条件是最小模特征值唯一
反幂法可以与原点位移技术结合来求任何一个已知近似范围的特征值,例如利用圆盘定理得到某个λj≈p\lambda_j\approx pλj≈p,则可以对B=A−pI\mathbf B=\mathbf A-p\mathbf IB=A−pI应用反幂法求得λj\lambda_jλj
矩阵的QR分解
Householder变换
定义 对单位向量w\mathbf ww,可以基于w\mathbf ww构造一个Householder矩阵(初等反射矩阵) :
H(w)=I−2wwT\mathbf H(\mathbf w)=\mathbf I-2\mathbf{ww}^T
H(w)=I−2wwT
说明 对任意v\mathbf vv构造分解
v=v∥+v⊥, v∥=kw, v⊥⊥w\mathbf v=\mathbf v^{\parallel}+\mathbf v^{\perp},\ \mathbf v^\parallel=k\mathbf w,\ \mathbf v^\perp\perp\mathbf w
v=v∥+v⊥, v∥=kw, v⊥⊥w
则
Hv∥=v∥−2wwTv∥=kw−2wwT(kw)=kw−2kw=−kw=−v∥\mathbf{Hv}^\parallel=\mathbf v^\parallel-2\mathbf{ww}^T\mathbf v^\parallel=k\mathbf w-2\mathbf {ww}^T(k\mathbf w)=k\mathbf{w}-2k\mathbf w=-k\mathbf w=-\mathbf v^\parallel
Hv∥=v∥−2wwTv∥=kw−2wwT(kw)=kw−2kw=−kw=−v∥
Hv⊥=v⊥−2wwTv⊥=v⊥\mathbf{Hv}^\perp=\mathbf v^{\perp}-2\mathbf{ww}^T\mathbf v^\perp=\mathbf v^\perp
Hv⊥=v⊥−2wwTv⊥=v⊥
即H=I−2wwT\mathbf H=\mathbf I-2\mathbf{ww}^TH=I−2wwT满足以w\mathbf ww为镜面法向的反射矩阵的语义,而满足上述语义的变换应当是唯一的
命题(Householder矩阵的性质)
H(w)=H(−w)\mathbf H(\mathbf w)=\mathbf H(-\mathbf w)H(w)=H(−w)
H\mathbf HH是对称的、正交的
定理 对任意两个等长的向量x,y\mathbf x,\mathbf yx,y,存在Householder矩阵H\mathbf HH,使Hx=y\mathbf {Hx}=\mathbf yHx=y
证明 只需用w=1∣∣x−y∣∣2(x−y)\mathbf w=\frac{1}{||\mathbf x-\mathbf y||_2}(\mathbf x-\mathbf y)w=∣∣x−y∣∣21(x−y)构造Householder矩阵即可
这说明,我们可以使用Householder矩阵来消元:令y=∣∣x∣∣2e1\mathbf y=||\mathbf x||_2\mathbf e_1y=∣∣x∣∣2e1,则由此构造出的H\mathbf HH可以作为一个正交矩阵消去某个列向量x\mathbf xx
实际使用中,一般令y=−σe1=−sign(x1)∣∣x∣∣2e1\mathbf y=-\sigma\mathbf e_1=-{\rm sign}(x_1)||\mathbf x||_2\mathbf e_1y=−σe1=−sign(x1)∣∣x∣∣2e1,这样可以达到数值上更稳定的结果
使用Householder变换实现QR分解
定义(矩阵的QR分解) 对矩阵A∈Rm×n\mathbf A\in\mathbb R^{m\times n}A∈Rm×n,则如下形式的分解为其QR分解:
A=QR, Q∈Rm×m, R∈Rm×n\mathbf A=\mathbf{QR},\ \mathbf Q\in\mathbb R^{m\times m},\ \mathbf R\in\mathbb R^{m\times n}
A=QR, Q∈Rm×m, R∈Rm×n
其中Q\mathbf QQ正交,R\mathbf RR形如下面两种形状:
m>nm>nm>n时,R=[R0O]\mathbf R=\begin{bmatrix}\mathbf R_0\\\mathbf O\end{bmatrix}R=[R0O],其中R0\mathbf R_0R0为上三角方阵
m m=nm=nm=n时,R\mathbf RR为上三角方阵 在上面的三种情形中,m>nm>nm>n时有信息冗余,此时不必要要求Q\mathbf QQ是完整的正交矩阵,于是有: 定义(经济的QR分解) 对矩阵A∈Rm×n\mathbf A\in\mathbb R^{m\times n}A∈Rm×n,其中m>nm>nm>n,则如下形式的分解是其经济的QR分解: A=QR, Q∈Rm×n, R∈Rn×n\mathbf A=\mathbf{QR},\ \mathbf Q\in\mathbb R^{m\times n},\ \mathbf R\in\mathbb R^{n\times n} A=QR, Q∈Rm×n, R∈Rn×n 其中Q\mathbf QQ列正交(即其列向量组成标准正交向量组),R\mathbf RR为上三角方阵 定理 任意实矩阵存在QR分解;非奇异方阵存在唯一使R\mathbf RR对角元均正的QR分解 Householder矩阵构造QR分解的方式 显然,对矩阵的第一列,直接令y1=−σ1e1=−sign(a11)∣∣a1∣∣2e1\mathbf y_1=-\sigma_1\mathbf e_1=-{\rm sign}(a_{11})||\mathbf a_1||_2\mathbf e_1y1=−σ1e1=−sign(a11)∣∣a1∣∣2e1,构造Householder矩阵H1\mathbf H_1H1使H1a1=y1\mathbf H_1\mathbf a_1=\mathbf y_1H1a1=y1即可 此时,原矩阵将变成 H1A=[−σ1r1T0A2]\mathbf H_1\mathbf A=\begin{bmatrix}-\sigma_1&\mathbf r_1^T\\\mathbf 0&\mathbf A_2\end{bmatrix} H1A=[−σ10r1TA2] 此时可以用n−1n-1n−1阶的Householder矩阵H2′\mathbf H_2'H2′消去A2\mathbf A_2A2的第一列,这等价于用 H2=[1H2′]\mathbf H_2=\begin{bmatrix}1&\\&\mathbf H_2'\end{bmatrix} H2=[1H2′] 显然,H2\mathbf H_2H2仍然是正交矩阵 因此,我们类似构造即可得到 Hn−1Hn−2⋯H2H1A=R\mathbf H_{n-1}\mathbf H_{n-2}\cdots\mathbf H_2\mathbf H_1\mathbf A=\mathbf R Hn−1Hn−2⋯H2H1A=R 于是 Q=H1−1H2−1⋯Hn−1−1\mathbf Q=\mathbf H_1^{-1}\mathbf H_2^{-1}\cdots\mathbf H_{n-1}^{-1} Q=H1−1H2−1⋯Hn−1−1 是正交矩阵 QR分解的实际计算过程 考虑Householder矩阵的构造: H1=I−2wwT\mathbf H_1=\mathbf I-2\mathbf {ww}^T H1=I−2wwT 其中 w=a1+σe1∣∣a1+σe1∣∣2\mathbf w=\frac{\mathbf a_1+\sigma\mathbf e_1}{||\mathbf a_1+\sigma\mathbf e_1||_2} w=∣∣a1+σe1∣∣2a1+σe1 如果我们令 v1=a1+σe1\mathbf v_1=\mathbf a_1+\sigma\mathbf e_1 v1=a1+σe1 则 w=v1∣∣v1∣∣2\mathbf w=\frac{\mathbf v_1}{||\mathbf v_1||_2} w=∣∣v1∣∣2v1 于是 H1=I−2v1v1Tv1Tv1\mathbf H_1=\mathbf I-2\frac{\mathbf v_1\mathbf v_1^T}{\mathbf v_1^T\mathbf v_1} H1=I−2v1Tv1v1v1T 这就是不要求法向量单位的情况下的HouseHolder矩阵 在此情况下作矩阵向量乘则转化为 H1aj=aj−2v1Tajv1Tv1v1\mathbf H_1\mathbf a_j=\mathbf a_j-2\frac{\mathbf v_1^T\mathbf a_j}{\mathbf v_1^T\mathbf v_1}\mathbf v_1 H1aj=aj−2v1Tv1v1Tajv1 这就是我们熟悉的Gauss-Schmidt正交化过程 此时,v1Tv1\mathbf v_1^T\mathbf v_1v1Tv1实际可以只计算一次,每次矩阵乘向量的实际计算复杂度变为了向量内积和向量线性运算 这样进行分解的总复杂度为O(mn2)O(mn^2)O(mn2) Givens旋转变换进行QR分解 正交矩阵有两种形态:反射和旋转,前面研究的是利用反射变换进行QR分解的思路,下面我们讨论旋转变换 定义(二阶Givens矩阵) 二阶Givens矩阵形态如下: G=[cosθsinθ−sinθcosθ]\mathbf G=\begin{bmatrix}\cos\theta&\sin\theta\\-\sin\theta&\cos\theta\end{bmatrix} G=[cosθ−sinθsinθcosθ] 定义(n阶Givens矩阵) nnn阶Givens矩阵为二阶矩阵嵌入nnn阶单位阵的结果,例子如下: G=[1cosθsinθ1−sinθcosθ1]\mathbf G=\begin{bmatrix}1\\&\cos\theta&&\sin\theta\\&&1\\&-\sin\theta&&\cos\theta\\&&&&1\end{bmatrix} G=1cosθ−sinθ1sinθcosθ1 简化的表示为 G=[1cs1−sc1], c2+s2=1\mathbf G=\begin{bmatrix}1\\&c&&s\\&&1\\&-s&&c\\&&&&1\end{bmatrix},\ c^2+s^2=1 G=1c−s1sc1, c2+s2=1 我们可以通过选择ccc和sss的值将一个二维向量化为基向量的形状: G[x1x2]=[α0]\mathbf G\begin{bmatrix}x_1\\x_2\end{bmatrix}=\begin{bmatrix}\alpha\\0\end{bmatrix} G[x1x2]=[α0] 只需令 c=x1α, s=x2αc=\frac{x_1}\alpha,\ s=\frac{x_2}\alpha c=αx1, s=αx2 注意,此时必须满足 α2=x12+x22\alpha^2=x_1^2+x_2^2 α2=x12+x22 与Householder变换不同,Givens变换只影响向量的两行,因此,如果要将一个nnn维稠密向量进行消去,需要n−1n-1n−1个Givens变换,分别在第111维和第k=2,⋯ ,nk=2,\cdots,nk=2,⋯,n维构成的平面上旋转;而Householder变换只需要1个 因此,虽然单个Givens变换的计算量更小,但其对稠密矩阵的计算量很大。Givens变换进行QR分解更适合稀疏矩阵 求矩阵的所有特征值 其思路是利用正交相似变换(不改变特征值)将矩阵转化为三角阵或分块三角阵(容易求出特征值) 思路:收缩技术 假设我们已知矩阵A\mathbf AA的一个特征向量x\mathbf xx(这可以通过幂法得到) 那么我们可以构造一个正交变换H\mathbf HH使得Hx=σe1\mathbf {Hx}=\sigma\mathbf e_1Hx=σe1,则 HAHTe1=HA(1σx)=1σH(λx)=λe1\mathbf {HAH}^T\mathbf e_1=\mathbf{HA}(\frac1\sigma\mathbf x)=\frac1\sigma\mathbf H(\lambda\mathbf x)=\lambda\mathbf e_1 HAHTe1=HA(σ1x)=σ1H(λx)=λe1 于是e1\mathbf e_1e1是正交相似于A\mathbf AA的矩阵HAHT\mathbf {HAH}^THAHT的同一特征值的特征向量 于是我们可以得到 HAHT=[λr1TA2]\mathbf {HAH}^T=\begin{bmatrix}\lambda&\mathbf r_1^T\\&\mathbf A_2\end{bmatrix} HAHT=[λr1TA2] 此时问题就变小了 然而,利用收缩技术和幂法来求解,效率不高,而且误差会累积 QR迭代算法 定理(实Schur分解) 对于任何A∈Rn×n\mathbf A\in\mathbb R^{n\times n}A∈Rn×n,存在正交阵Q\mathbf QQ使得 QTAQ=S\mathbf Q^T\mathbf {AQ}=\mathbf S QTAQ=S 其中S\mathbf SS是拟上三角阵,即对角块不超过2阶的分块上三角阵 注意:在复数域下,Schur分解可以得到真上三角阵,但因为实矩阵的特征值可能有复数,因此无法得到全实的真Schur分解,但由于实矩阵的复特征值一定共轭成对出现,因此可以将成对出现的复特征值包裹在一个2阶实对角块中 QR迭代算法 对Ak\mathbf A_kAk做QR分解:Ak=QkRk\mathbf A_k=\mathbf Q_k\mathbf R_kAk=QkRk 作迭代:Ak+1=RkQk\mathbf A_{k+1}=\mathbf R_k\mathbf Q_kAk+1=RkQk 这即Ak+1=QkTAQk\mathbf A_{k+1}=\mathbf Q_k^T\mathbf {AQ}_kAk+1=QkTAQk,构成了一组正交相似矩阵序列 定理(QR迭代算法的收敛定理) 若矩阵A\mathbf AA的所有等模特征值组或者是实重特征值,或者是复共轭特征值,则QR迭代所得的矩阵序列基本收敛于拟上三角阵 每一步QR迭代的计算量很大,且QR迭代可能不收敛或慢收敛 为此,提出以下算法: 实用的QR算法改进技术 降低QR分解的计算量:Hessenberg化 上Hessenberg型: [×××××××××××××××××××]\begin{bmatrix}\times&\times&\times&\times&\times\\\times&\times&\times&\times&\times\\ &\times&\times&\times&\times\\&&\times&\times&\times\\&&&\times&\times \end{bmatrix} ××××××××××××××××××× 由于正交相似变换不改变矩阵特征值,因此我们使用正交相似变换来进行上述转化。首先将矩阵分块为 A=[a11r1Tc1A22(1)]\mathbf A=\begin{bmatrix}a_{11}&\mathbf r_1^T\\\mathbf c_1&\mathbf A_{22}^{(1)}\end{bmatrix} A=[a11c1r1TA22(1)] 可以针对左下角块构造一个Householder变换使得 H1′c1=σ1e1\mathbf H_1'\mathbf c_1=\sigma_1\mathbf e_1 H1′c1=σ1e1 令 H1=[10T0H1′]\mathbf H_1=\begin{bmatrix}1&\mathbf 0^T\\\mathbf 0&\mathbf H_1'\end{bmatrix} H1=[100TH1′] 于是 H1A=[a11r1Tσ1e1H1′A22(1)]\mathbf H_1\mathbf A=\begin{bmatrix}a_{11}&\mathbf r_1^T\\\sigma_1\mathbf e_1&\mathbf H_1'\mathbf A_{22}^{(1)}\end{bmatrix} H1A=[a11σ1e1r1TH1′A22(1)] H1AH1=[a11r1TH1′σ1e1H1′A22(1)H1′]\mathbf H_1\mathbf A\mathbf H_1=\begin{bmatrix}a_{11}&\mathbf r_1^T\mathbf H_1'\\\sigma_1\mathbf e_1&\mathbf H_1'\mathbf A_{22}^{(1)}\mathbf H_1'\end{bmatrix} H1AH1=[a11σ1e1r1TH1′H1′A22(1)H1′] 由于H1\mathbf H_1H1是对称矩阵,上述结果就是对A\mathbf AA进行一步正交相似变换得到的结果 我们可以对其右下角块迭代上述正交相似变换,从而将A\mathbf AA正交相似于Hessenberg矩阵 对Hessenberg矩阵进行QR分解的计算量小于一般的稠密矩阵 改善收敛性:原点位移技术 常规QR迭代的单步迭代操作为: {QkRk=AkAk+1=RkQk\begin{cases}\mathbf Q_k\mathbf R_k=\mathbf A_k\\ \mathbf A_{k+1}=\mathbf R_k\mathbf Q_k \end{cases} {QkRk=AkAk+1=RkQk 此时可以在其中加入位移因子,使得每步迭代前后仍然正交相似: {QkRk=Ak−skIAk+1=RkQk+skI\begin{cases}\mathbf Q_k\mathbf R_k=\mathbf A_k-s_k\mathbf I\\ \mathbf A_{k+1}=\mathbf R_k\mathbf Q_k+s_k\mathbf I \end{cases} {QkRk=Ak−skIAk+1=RkQk+skI 对于位移因子sks_ksk的选取,可以使用Rayleigh策略:取sks_ksk为Ak\mathbf A_kAk的(n,n)(n,n)(n,n)元 矩阵的奇异值分解 基本概念 奇异值和奇异向量是针对非方阵的,特征值和特征向量的扩展定义 定义(奇异值、奇异向量) 对矩阵A∈Rm×n\mathbf A\in\mathbb R^{m\times n}A∈Rm×n,如果 Av=σuATu=σvσ≥0, u≠0∈Rm, v≠0∈Rn\mathbf {Av}=\sigma\mathbf u\\ \mathbf A^T\mathbf u=\sigma\mathbf v\\ \sigma\geq0,\ \mathbf u\neq\mathbf 0\in\mathbb R^{m},\ \mathbf v\neq\mathbf 0\in\mathbb R^{n} Av=σuATu=σvσ≥0, u=0∈Rm, v=0∈Rn 则称σ\sigmaσ为A\mathbf AA的奇异值,u,v\mathbf u,\mathbf vu,v分别为对应的左奇异向量和右奇异向量 定理(奇异值分解, SVD) 对任何矩阵A∈Rm×n, m≥n\mathbf A\in\mathbb R^{m\times n},\ m\geq nA∈Rm×n, m≥n,可以找到正交矩阵U,V\mathbf U,\mathbf VU,V和类对角矩阵Σ∈Rm×n\mathbf \Sigma\in\mathbb R^{m\times n}Σ∈Rm×n,使得 A=UΣVT\mathbf A=\mathbf {U\Sigma V}^T A=UΣVT 其中 Σ=[σ1σ2⋱σn0⋯], σ1≥σ2≥⋯≥σn≥0\mathbf \Sigma=\begin{bmatrix}\sigma_1\\&\sigma_2\\&&\ddots\\&&&\sigma_n\\0&\cdots\end{bmatrix},\ \sigma_1\geq\sigma_2\geq\cdots\geq\sigma_n\geq0 Σ=σ10σ2⋯⋱σn, σ1≥σ2≥⋯≥σn≥0 讨论 上述分解等价于 AV=UΣ\mathbf {AV}=\mathbf{U\Sigma} AV=UΣ 即 [Av1⋯Avn]=[σ1u1⋯σnun]\begin{bmatrix}\mathbf{Av}_1&\cdots&\mathbf{Av}_n\end{bmatrix}=\begin{bmatrix}\sigma_1\mathbf u_1&\cdots&\sigma_n\mathbf u_n\end{bmatrix} [Av1⋯Avn]=[σ1u1⋯σnun] 这就是奇异向量的第一部分定义,而第二部分定义则来自于 UTA=ΣVT\mathbf U^T\mathbf A=\mathbf \Sigma\mathbf V^T UTA=ΣVT 于是,奇异值分解存在等价于说,矩阵属于不同奇异值的奇异向量正交 讨论 这里的m≥nm\geq nm≥n不是必须的,只是为了叙述方便,如果m SVD的证明 对ATA\mathbf A^T\mathbf AATA作谱分解: ATA=VΛVT\mathbf A^T\mathbf A=\mathbf{V\Lambda V}^T ATA=VΛVT 令 Λ=ΣTΣ\mathbf \Lambda=\mathbf \Sigma^T\mathbf \Sigma Λ=ΣTΣ 于是 ATA=VΣTΣVT\mathbf A^T\mathbf A=\mathbf {V\Sigma}^T\mathbf{\Sigma V}^T ATA=VΣTΣVT 假设Σ\mathbf\SigmaΣ的前rrr个对角元非零,则取出这部分的对角矩阵作为Σr\mathbf\Sigma_rΣr,而V\mathbf VV的前rrr列作为Vr\mathbf V_rVr 于是令 Ur=AVrΣr−1\mathbf U_r=\mathbf{AV}_r\mathbf\Sigma_r^{-1} Ur=AVrΣr−1 (这里A,Vr,Σr\mathbf A,\mathbf V_r,\mathbf\Sigma_rA,Vr,Σr大小分别为m×nm\times nm×n,n×rn\times rn×r,r×rr\times rr×r,得到m×rm\times rm×r的Ur\mathbf U_rUr) 由于 UrTUr=Σr−1VrTATAVrΣr−1=Σr−1[σ12⋱σr2]Σr−1=Ir\mathbf U_r^T\mathbf U_r=\mathbf\Sigma_r^{-1}\mathbf{V}_r^T\mathbf{A}^T\mathbf{AV}_r\mathbf\Sigma_r^{-1}=\mathbf\Sigma_r^{-1}\begin{bmatrix}\sigma_1^2\\&\ddots\\&&\sigma_r^2\end{bmatrix}\mathbf\Sigma_r^{-1}=\mathbf I_r UrTUr=Σr−1VrTATAVrΣr−1=Σr−1σ12⋱σr2Σr−1=Ir 可知Ur\mathbf U_rUr列正交,将其扩充至m×mm\times mm×m即可得到U\mathbf UU 简化的SVD 上面的证明中,我们选择了非零奇异值及其对应的奇异向量,作为如下形式的分解: A=UrΣrVrT\mathbf A=\mathbf U_r\mathbf \Sigma_r\mathbf V_r^T A=UrΣrVrT 其中零奇异值对应的部分对乘积结果毫无影响,因此这就是SVD的简化形式 同时,我们还可以将其展开: A=[u1⋯ur][σ1⋱σr][v1T⋮vrT]=σ1u1v1T+⋯+σrurvrT\mathbf A=\begin{bmatrix}\mathbf u_1&\cdots&\mathbf u_r\end{bmatrix} \begin{bmatrix}\sigma_1\\&\ddots\\&&\sigma_r\end{bmatrix} \begin{bmatrix}\mathbf v_1^T\\\vdots\\\mathbf v_r^T\end{bmatrix}=\sigma_1\mathbf u_1\mathbf v_1^T+\cdots+\sigma_r\mathbf u_r\mathbf v_r^T A=[u1⋯ur]σ1⋱σrv1T⋮vrT=σ1u1v1T+⋯+σrurvrT 这是矩阵的一种秩一分解,从此也可以看出,非零奇异值的个数rrr对应矩阵的秩 SVD的几何意义 我们回顾一下谱分解的几何意义: A=XΛX−1\mathbf A=\mathbf{X\Lambda X}^{-1} A=XΛX−1 假设我们有一个向量 v=[v1⋮vn]\mathbf v=\begin{bmatrix}v_1\\\vdots\\v_n\end{bmatrix} v=v1⋮vn 现在我们找到了一组基 X=[x1⋯xn]\mathbf X=\begin{bmatrix}\mathbf x_1&\cdots&\mathbf x_n\end{bmatrix} X=[x1⋯xn] 于是,可以将v\mathbf vv用这组基表示: v=k1x1+k2x2+⋯+knxn=X[k1⋮kn]=[v1⋮vn]\mathbf v=k_1\mathbf x_1+k_2\mathbf x_2+\cdots+k_n\mathbf x_n=\mathbf X\begin{bmatrix}k_1\\\vdots\\k_n\end{bmatrix}=\begin{bmatrix}v_1\\\vdots\\v_n\end{bmatrix} v=k1x1+k2x2+⋯+knxn=Xk1⋮kn=v1⋮vn 于是 [k1⋮kn]=X−1[v1⋮vn]\begin{bmatrix}k_1\\\vdots\\k_n\end{bmatrix}=\mathbf X^{-1}\begin{bmatrix}v_1\\\vdots\\v_n\end{bmatrix} k1⋮kn=X−1v1⋮vn 即A=XΛX−1\mathbf A=\mathbf{X\Lambda X}^{-1}A=XΛX−1中的第一步X−1\mathbf X^{-1}X−1是将输入向量从标准基上的坐标转化为X\mathbf XX定义的基上的坐标 而根据特征向量的意义,A\mathbf AA对特征向量的操作仅仅是伸缩变换,因此A\mathbf AA在X\mathbf XX定义的基上的形式就是对角矩阵Λ\mathbf\LambdaΛ 而最后一步X\mathbf XX,则是第一步的逆过程,将输出向量从X\mathbf XX基上的坐标转换回原始坐标 于是,谱分解的几何意义是将一个映射A\mathbf AA拆分成了以下三步: 第一步:X−1\mathbf X^{-1}X−1,将输入向量从标准基换到特征基 第二步:Λ\mathbf\LambdaΛ,将特征基下的输入向量变换得到输出向量 第三步:X\mathbf XX,将输出向量从特征基换回标准基 对于实对称矩阵而言,我们有更好的谱分解形式: A=QΛQT\mathbf A=\mathbf{Q\Lambda Q}^T A=QΛQT 此形式与一般谱分解的区别是,将特征基取成了一组标准正交基,这是实对称矩阵才具有的独特性质 接下来我们回到SVD。SVD相比谱分解而言,不要求矩阵是方阵,这意味着,矩阵的定义域即输入空间,和陪域即输出空间,从本质上不再是同一个空间,于是,我们也不再能够找到一组同时对输入空间和输出空间有效的特征基 此时,我们就需要对输入空间和输出空间分别取一组基(可以分别称为右奇异基和左奇异基),SVD就意味着将原始矩阵的映射进行如下拆分: 第一步:VT\mathbf V^TVT,将输入向量从标准基换到右奇异基 第二步:Σ\mathbf\SigmaΣ,将右奇异基下的输入向量变换,得到左奇异基下的输出向量 第三步:U\mathbf UU,将输出向量从左奇异基换到标准基 这样看起来,这个右奇异基和左奇异基的选取似乎就非常任意了,但注意,Σ\mathbf\SigmaΣ是拟对角阵,这意味着,从输入空间(右奇异基)到输出空间(左奇异基)的映射,其实是沿着坐标轴的伸缩映射,但由于输入输出空间维数可能不相等,这个过程中可能会舍弃一些维度,或增加一些值为000的维度 接下来,我们将U\mathbf UU和V\mathbf VV分别拆分成前rrr列和后若干列: U=[UrUm−r]\mathbf U=\begin{bmatrix}\mathbf U_r&\mathbf U_{m-r}\end{bmatrix} U=[UrUm−r] V=[VrVn−r]\mathbf V=\begin{bmatrix}\mathbf V_r&\mathbf V_{n-r}\end{bmatrix} V=[VrVn−r] 根据Σ\mathbf\SigmaΣ的非零值只有前rrr个对角元,我们可以知道这四个子矩阵的列空间与映射A\mathbf AA的一些关系: R(Vn−r)\mathcal R(\mathbf V_{n-r})R(Vn−r)经A\mathbf AA映射后为0\mathbf 00 所有向量经A\mathbf AA映射后必然落入R(Ur)\mathcal R(\mathbf U_r)R(Ur) 于是,我们可以得到SVD与四个基本子空间的关系: 列空间:R(Ur)\mathcal R(\mathbf U_r)R(Ur),列空间的基是U\mathbf UU的前rrr列 零空间:R(Vn−r)\mathcal R(\mathbf V_{n-r})R(Vn−r),零空间的基是V\mathbf VV的后n−rn-rn−r列 行空间:零空间的正交补,即V\mathbf VV的前rrr列 左零空间:列空间的正交补,即U\mathbf UU的后m−rm-rm−r列 此时,我们就可以从简化SVD的角度再次讨论几何意义: 第一步:VrT\mathbf V_r^TVrT,将输入向量中属于零空间的部分舍弃,属于零空间正交补(即行空间)的部分映射到一个rrr维空间中,得到的结果实际上就是输入向量在行空间的右奇异基下的坐标表示 第二步:Σr\mathbf \Sigma_rΣr,对行空间下用右奇异基表示的输入向量进行对角矩阵变换,得到列空间下用左奇异基表示的输出向量 第三步:Ur\mathbf U_rUr,将列空间下左奇异基表示的输出向量转化为陪域下标准基表示的输出向量,此时不属于列空间的那部分固定为000 从这个角度,每个向量在简化SVD的视角其实经历了以下四个空间: 此时我们也可以考虑矩阵的广义逆: A∗=VrΣr−1UrT\mathbf A^*=\mathbf {V}_r\mathbf{\Sigma}_r^{-1}\mathbf U_r^T A∗=VrΣr−1UrT 显然,它将上述空间变换的过程取了逆序,同时将行空间到列空间的对角变换也取了逆。广义逆在行空间到列空间的作用上与原矩阵是互逆的,而对零空间和左零空间中的内容会舍弃 SVD与矩阵范数 矩阵的2-范数为矩阵的最大奇异值 矩阵的Frobinius范数为所有奇异值平方和的算术平方根 低秩逼近 假设我们已经有了上述矩阵秩一分解(回忆:这其实是简化SVD的另一种表示形式): A=σ1u1v1T+⋯+σrurvrT\mathbf A=\sigma_1\mathbf u_1\mathbf v_1^T+\cdots+\sigma_r\mathbf u_r\mathbf v_r^T A=σ1u1v1T+⋯+σrurvrT 对于r0 Ar0=σ1u1v1T+⋯+σr0ur0vr0T\mathbf A_{r_0}=\sigma_1\mathbf u_1\mathbf v_1^T+\cdots+\sigma_{r_0}\mathbf u_{r_0}\mathbf v_{r_0}^T Ar0=σ1u1v1T+⋯+σr0ur0vr0T 由于所有奇异值是按从大到小排序的,上面取到的就是最佳的低秩逼近结果,可以证明,它与原矩阵的误差的范数(在2-范数或F-范数下)是最小的