特征值相关概念

基础概念

矩阵A\mathbf A的特征值和特征向量：

Ax=λx, x≠0\mathbf {Ax}=\lambda\mathbf x,\ \mathbf x\neq\mathbf 0

特征值的求法：

λ: det⁡(λI−A)=0\lambda:\ \det(\lambda\mathbf I-\mathbf A)=0

特征向量的求法：

x: (λI−A)x=0\mathbf x:\ (\lambda\mathbf I-\mathbf A)\mathbf x=\mathbf 0

特征向量构成特征子空间

矩阵的特征值谱（即特征值组成的集合）记作λ(A)\lambda(\mathbf A)

基本结论

定理 矩阵奇异当且仅当矩阵特征值含有00

定理 λ(A)=λ(AT)\lambda(\mathbf A)=\lambda(\mathbf A^T)

证明 转置不改变行列式。det⁡(λI−A)=det⁡(λI−AT)\det(\lambda\mathbf I-\mathbf A)=\det(\lambda\mathbf I-\mathbf A^T)

定理 对角阵或上（下）三角阵的特征值为对角元

定理 分块对角阵或分块上（下）三角阵的特征值谱是各个对角块特征值谱的并集

定义 对矩阵A,B\mathbf A,\mathbf B，若存在可逆阵X\mathbf X使得A=XBX−1\mathbf A=\mathbf{XB}\mathbf X^{-1}，则称A,B\mathbf A,\mathbf B相似

定理 相似矩阵特征值谱相同

证明 对A\mathbf A的特征对(λ,v)(\lambda,\mathbf v)

Bv=X−1AXv=X−1(λXv)=λv\mathbf{Bv}=\mathbf{X}^{-1}\mathbf{AXv}=\mathbf{X}^{-1}(\lambda\mathbf{Xv})=\lambda\mathbf v

定理（矩阵运算结果的特征值）

设λj\lambda_j为A\mathbf A的特征值，则

  • λ(cA)={cλj}\lambda(c\mathbf A)=\{c\lambda_j\}
  • λ(A+cI)={λj+c}\lambda(\mathbf A+c\mathbf I)=\{\lambda_j+c\}
  • λ(P(A))={P(λj)}\lambda(P(\mathbf A))=\{P(\lambda_j)\}，其中PP为多项式
  • λ(A−1)={λj−1}\lambda(\mathbf A^{-1})=\{\lambda_j^{-1}\}

特征值的重数、对角化、亏损阵

定义（特征值的重数）

对A\mathbf A的特征值λ0\lambda_0，此特征值在特征方程det⁡(λI−A)=0\det(\lambda\mathbf I-\mathbf A)=0中的重数定义为代数重数，特征子空间的维数dim⁡N(λ0I−A)\dim\mathcal N(\lambda_0\mathbf I-\mathbf A)定义为几何重数

定理（代数重数和几何重数的关系） 特征值的代数重数大于或等于其几何重数

证明

设矩阵A\mathbf A特征值λj\lambda_j的代数重数为njn_j，几何重数为kjk_j，则可以在其特征子空间N(λjI−A)\mathcal N(\lambda_j\mathbf I-\mathbf A)中找到一组基：x1,⋯ ,xkj\mathbf x_1,\cdots,\mathbf x_{k_j}

将此基扩充成Rn\mathbb R^n的一组基：x1,⋯ ,xn\mathbf x_1,\cdots,\mathbf x_n

令X=[x1⋯xn]\mathbf X=\begin{bmatrix}\mathbf x_1&\cdots&\mathbf x_n\end{bmatrix}，则有

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}

于是

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)

即(λ−λj)kj(\lambda-\lambda_j)^{k_j}整除A\mathbf A的特征多项式PA(λ)P_{\mathbf A}(\lambda)，故nj≥kjn_j\geq k_j

定理 不同特征值的特征向量线性无关

证明 利用特征值的定义显然

定义 如果所有特征值的几何重数都等于代数重数，则矩阵称为非亏损阵，否则为亏损阵

非亏损阵的特征向量可以组成全空间的一组基

定理 非亏损阵和可对角化矩阵是同样的概念，即，矩阵A\mathbf A非亏损等价于存在X\mathbf X使X−1AX=Λ\mathbf X^{-1}\mathbf {AX}=\mathbf{\Lambda}

证明 取一组特征向量，得到的矩阵就是X\mathbf X

实对称矩阵的对角化

定理 实对称矩阵的特征值是实数

证明 对A\mathbf A的特征对(λ,x)(\lambda, \mathbf x)，有

Ax=λx\mathbf {Ax}=\lambda\mathbf x

两边同取共轭，有

Axˉ=λˉxˉ\mathbf A\bar{\mathbf x}=\bar\lambda \bar{\mathbf x}

对第一个式子左乘xˉT\bar{\mathbf x}^T，有

xˉTAx=λxˉTx\bar{\mathbf x}^T\mathbf {Ax}=\lambda\bar{\mathbf x}^T\mathbf x

对第二个式子左乘xT\mathbf x^T，有

xTAxˉ=λˉxTxˉ\mathbf x^T\mathbf A\bar{\mathbf x}=\bar\lambda \mathbf x^T\bar{\mathbf x}

由于A\mathbf A对称

xˉTAx=xTAxˉ\bar{\mathbf x}^T\mathbf {Ax}=\mathbf x^T\mathbf A\bar{\mathbf x}

于是

λxˉTx=λˉxˉTx\lambda\bar{\mathbf x}^T\mathbf x=\bar\lambda\bar{\mathbf x}^T\mathbf x

又由于x≠0\mathbf x\neq\mathbf 0

xˉTx>0\bar{\mathbf x}^T\mathbf x>0

于是

λ=λˉ\lambda=\bar\lambda

定理（实对称矩阵的谱分解/正交对角化） 设A\mathbf A实对称，则存在正交矩阵Q\mathbf Q，使A=QΛQT\mathbf A=\mathbf{Q\Lambda Q}^T

证明 对阶数nn归纳，n=1n=1时显然，假设n−1n-1阶成立

对nn阶的A\mathbf A，可以找到一组实特征对(λ1,q1)(\lambda_1,\mathbf q_1)，不妨设∣∣q1∣∣=1||\mathbf q_1||=1，将q1\mathbf q_1扩充成一组实标准正交基q1,⋯ ,qn\mathbf q_1,\cdots,\mathbf q_n，设

Q1=[q1⋯qn]\mathbf Q_1=\begin{bmatrix}\mathbf q_1&\cdots&\mathbf q_n\end{bmatrix}

则

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}

于是

Q1TAQ1=[λ1bTC]\mathbf Q_1^T\mathbf{AQ}_1=\begin{bmatrix}\lambda_1&\mathbf b^T\\&\mathbf C\end{bmatrix}

左端是实对称矩阵，于是

[λ1bTC]=[λ1C]\begin{bmatrix}\lambda_1&\mathbf b^T\\&\mathbf C\end{bmatrix}=\begin{bmatrix}\lambda_1&\\&\mathbf C\end{bmatrix}

且C\mathbf C也是实对称矩阵，根据归纳假设知存在

C=Q~2Λ2Q~2T\mathbf C=\tilde{\mathbf Q}_2\mathbf{\Lambda}_2\tilde{\mathbf Q}_2^T

令

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}

则

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

于是

A=Q1Q2ΛQ2TQ1T\mathbf A=\mathbf Q_1\mathbf Q_2\mathbf{\Lambda Q}_2^T\mathbf Q_1^T

令

Q=Q1Q2\mathbf Q=\mathbf Q_1\mathbf Q_2

证毕

定理 实对称矩阵属于不同特征值的特征向量正交

证明 可以直接根据正交对角化得到此结论，也可以直接证明：

设A\mathbf A有两个特征对(λ1,x1),(λ2,x2)(\lambda_1,\mathbf x_1),(\lambda_2,\mathbf x_2)，满足λ1≠λ2\lambda_1\neq\lambda_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

于是

x1Tx2=0\mathbf x_1^T\mathbf x_2=0

Jordan分解

定理（Jordan分解） 任意方阵A\mathbf A可以进行如下分解：

A=XJX−1\mathbf A=\mathbf{XJX}^{-1}

其中

J=[J1⋱Jp]\mathbf J=\begin{bmatrix}\mathbf J_1&&\\&\ddots&\\&&\mathbf J_p\end{bmatrix}

每个Js\mathbf J_s称为Jordan块，形如

Js=[λs1λs⋱⋱1λs]\mathbf J_s=\begin{bmatrix}\lambda_s&1\\&\lambda_s&\ddots\\&&\ddots&1\\&&&\lambda_s\end{bmatrix}

每个特征值可能有多个Jordan块，但其Jordan块总数为其几何重数，总阶数为其代数重数

Jordan分解是谱分解A=XΛX−1\mathbf A=\mathbf{X\Lambda X}^{-1}对亏损阵的推广，得到的J\mathbf J称为Jordan标准形

Jordan分解的理论基础和求法与数值分析课程关系不大， 在这里不展开，详见Jordan分解

特征值的分布

特征值的分布对于迭代法解线性方程组的收敛性有很强的影响，例如

ρ(B)=max⁡j∣λj(B)∣\rho(\mathbf B)=\max_j|\lambda_j(\mathbf 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)}}

定义（Gerschgorin圆盘） 设A∈Cn×n\mathbf A\in\mathbb C^{n\times n}，rk=∑j=1,j≠kn∣akj∣r_k=\sum\limits_{j=1,j\neq k}^n|a_{kj}|，在复平面上以akka_{kk}为圆心，rkr_k为半径的圆称为A\mathbf A的Gerschgorin圆盘

定理（圆盘定理） A\mathbf A的特征值在nn个圆盘的并集上。圆盘并集的每个连通部分内包含的特征值数量等于其包含的圆盘数量（计算重数）

证明

对特征对(λ,x)(\lambda,\mathbf x)，设xkx_k是x\mathbf x的最大分量，由于

Ax=λx\mathbf {Ax}=\lambda\mathbf x

知

λxk=∑jakjxj\lambda x_k=\sum\limits_j a_{kj}x_j

于是

(λ−akk)xk=∑j≠kakjxj(\lambda-a_{kk})x_k=\sum\limits_{j\neq k}a_{kj}x_j

于是

∣λ−akk∣≤∑j≠k∣akj∣|\lambda-a_{kk}|\leq\sum\limits_{j\neq k}|a_{kj}|

证毕

注 圆盘定理的一个直接推论是按行严格对角占优矩阵可逆，因为这样的矩阵的所有圆盘都不包含原点；且由于转置不改变特征值，这里的按行也可以是按列

推论 对对角元均正的实矩阵：

  • 严格对角占优的情况下，特征值实部均正
  • 严格对角占优且对称的情况下，对称正定

圆盘定理可以用来估计矩阵的特征值分布，可以使用的技巧有：

  • 注意到对A\mathbf A和AT\mathbf A^T使用圆盘定理可以得到不同的圆盘（行圆盘和列圆盘不同），可以选择二者中范围较小的来使用
  • 实矩阵如果有虚特征值，一定共轭成对出现，因此如果实矩阵一个圆盘内只有一个特征值，此特征值一定是实的

幂法和反幂法

幂法与规格化幂法

定义 矩阵模最大的特征值称为主特征值，其特征向量称为主特征向量

注 主特征值不一定唯一

算法（幂法）

随机取非零向量v0\mathbf v_0，作迭代法

vk+1=Avk\mathbf v_{k+1}=\mathbf {Av}_k

可以发现相邻的vk\mathbf v_k逐渐呈倍数关系，其倍数为一个主特征值，而vk\mathbf v_k则趋近于主特征值的特征向量

定理 若矩阵有唯一的主特征值，则幂法收敛于主特征值的特征向量

证明 将v0\mathbf v_0对谱分解或Jordan分解中矩阵X\mathbf X构成的基（即特征向量基或Jordan向量基）分解即可

注 幂法中，vk≈λ1kx1\mathbf v_k\approx \lambda_1^k\mathbf x_1，此时迭代次数较多时很容易上溢或下溢，因此实际计算中对每个迭代向量规格化（即使模最大的分量模为11，与标准化类似但计算量比标准化小，会得到长度略大于11的一个向量），以防止溢出

注 决定幂法收敛速度的主要因素是模次大特征值（“次主特征值”）与主特征值的模之比

当幂法加入规格化后，规格化向量经一步迭代后得到的最大模长元素即可作为主特征值的近似值（因为此元素在迭代前为11），减少了一次除法，这也是规格化相比标准化的优势

加速幂法的收敛

原点位移技术

考虑到B=A−sI\mathbf B=\mathbf A-s\mathbf I的特征值是λ(A)−s\lambda(\mathbf A)-s，如果

∣λ2−sλ1−s∣<∣λ2λ1∣\left|\frac{\lambda_2-s}{\lambda_1-s}\right|<\left|\frac{\lambda_2}{\lambda_1}\right|

则对B\mathbf B应用幂法可以加速收敛

Rayleigh商加速

定义 实矩阵A\mathbf A对实向量x\mathbf x的Rayleigh商定义为

R(x)=xTAxxTxR(\mathbf x)=\frac{\mathbf x^T\mathbf{Ax}}{\mathbf x^T\mathbf x}

定理 A\mathbf A对任何向量的瑞利商不大于其最大特征值，不小于其最小特征值

定理 幂法中的规格化向量uk\mathbf u_k满足

R(uk)=λ1+O((λ2λ1)2k)R(\mathbf u_k)=\lambda_1+O\left(\left(\frac{\lambda_2}{\lambda_1}\right)^{2k}\right)

注 对于幂法中取uk\mathbf u_k最大模分量的操作，其得到的结果为λ1+O((λ2λ1)k)\lambda_1+O\left(\left(\frac{\lambda_2}{\lambda_1}\right)^{k}\right)，收敛不如瑞利商法快

反幂法

由于A−1\mathbf A^{-1}的特征值是A\mathbf A的倒数，对A−1\mathbf A^{-1}应用幂法即可得到模最小特征值

此时每步迭代可以通过求解线性方程组来计算A−1\mathbf A^{-1}乘向量（避免求逆），但计算量仍然比正幂法大

适用条件是最小模特征值唯一

反幂法可以与原点位移技术结合来求任何一个已知近似范围的特征值，例如利用圆盘定理得到某个λj≈p\lambda_j\approx p，则可以对B=A−pI\mathbf B=\mathbf A-p\mathbf I应用反幂法求得λj\lambda_j

矩阵的QR分解

Householder变换

定义 对单位向量w\mathbf w，可以基于w\mathbf w构造一个Householder矩阵（初等反射矩阵） ：

H(w)=I−2wwT\mathbf H(\mathbf w)=\mathbf I-2\mathbf{ww}^T

说明 对任意v\mathbf v构造分解

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

则

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⊥=v⊥\mathbf{Hv}^\perp=\mathbf v^{\perp}-2\mathbf{ww}^T\mathbf v^\perp=\mathbf v^\perp

即H=I−2wwT\mathbf H=\mathbf I-2\mathbf{ww}^T满足以w\mathbf w为镜面法向的反射矩阵的语义，而满足上述语义的变换应当是唯一的

命题（Householder矩阵的性质）

  • H(w)=H(−w)\mathbf H(\mathbf w)=\mathbf H(-\mathbf w)
  • H\mathbf H是对称的、正交的

定理 对任意两个等长的向量x,y\mathbf x,\mathbf y，存在Householder矩阵H\mathbf H，使Hx=y\mathbf {Hx}=\mathbf y

证明 只需用w=1∣∣x−y∣∣2(x−y)\mathbf w=\frac{1}{||\mathbf x-\mathbf y||_2}(\mathbf x-\mathbf y)构造Householder矩阵即可

这说明，我们可以使用Householder矩阵来消元：令y=∣∣x∣∣2e1\mathbf y=||\mathbf x||_2\mathbf e_1，则由此构造出的H\mathbf H可以作为一个正交矩阵消去某个列向量x\mathbf x

实际使用中，一般令y=−σe1=−sign(x1)∣∣x∣∣2e1\mathbf y=-\sigma\mathbf e_1=-{\rm sign}(x_1)||\mathbf x||_2\mathbf e_1，这样可以达到数值上更稳定的结果

使用Householder变换实现QR分解

定义（矩阵的QR分解） 对矩阵A∈Rm×n\mathbf A\in\mathbb R^{m\times 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}

其中Q\mathbf Q正交，R\mathbf R形如下面两种形状：

  • m>nm>n时，R=[R0O]\mathbf R=\begin{bmatrix}\mathbf R_0\\\mathbf O\end{bmatrix}，其中R0\mathbf R_0为上三角方阵
  • m<nm<n时，R=[R1R2]\mathbf R=\begin{bmatrix}\mathbf R_1&\mathbf R_2\end{bmatrix}，其中R1\mathbf R_1为上三角方阵
  • m=nm=n时，R\mathbf R为上三角方阵

在上面的三种情形中，m>nm>n时有信息冗余，此时不必要要求Q\mathbf Q是完整的正交矩阵，于是有：

定义（经济的QR分解） 对矩阵A∈Rm×n\mathbf A\in\mathbb R^{m\times n}，其中m>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}

其中Q\mathbf Q列正交（即其列向量组成标准正交向量组），R\mathbf R为上三角方阵

定理 任意实矩阵存在QR分解；非奇异方阵存在唯一使R\mathbf R对角元均正的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_1，构造Householder矩阵H1\mathbf H_1使H1a1=y1\mathbf H_1\mathbf a_1=\mathbf y_1即可

此时，原矩阵将变成

H1A=[−σ1r1T0A2]\mathbf H_1\mathbf A=\begin{bmatrix}-\sigma_1&\mathbf r_1^T\\\mathbf 0&\mathbf A_2\end{bmatrix}

此时可以用n−1n-1阶的Householder矩阵H2′\mathbf H_2'消去A2\mathbf A_2的第一列，这等价于用

H2=[1H2′]\mathbf H_2=\begin{bmatrix}1&\\&\mathbf H_2'\end{bmatrix}

显然，H2\mathbf H_2仍然是正交矩阵

因此，我们类似构造即可得到

Hn−1Hn−2⋯H2H1A=R\mathbf H_{n-1}\mathbf H_{n-2}\cdots\mathbf H_2\mathbf H_1\mathbf A=\mathbf R

于是

Q=H1−1H2−1⋯Hn−1−1\mathbf Q=\mathbf H_1^{-1}\mathbf H_2^{-1}\cdots\mathbf H_{n-1}^{-1}

是正交矩阵

QR分解的实际计算过程

考虑Householder矩阵的构造：

H1=I−2wwT\mathbf H_1=\mathbf I-2\mathbf {ww}^T

其中

w=a1+σe1∣∣a1+σe1∣∣2\mathbf w=\frac{\mathbf a_1+\sigma\mathbf e_1}{||\mathbf a_1+\sigma\mathbf e_1||_2}

如果我们令

v1=a1+σe1\mathbf v_1=\mathbf a_1+\sigma\mathbf e_1

则

w=v1∣∣v1∣∣2\mathbf w=\frac{\mathbf v_1}{||\mathbf v_1||_2}

于是

H1=I−2v1v1Tv1Tv1\mathbf H_1=\mathbf I-2\frac{\mathbf v_1\mathbf v_1^T}{\mathbf v_1^T\mathbf v_1}

这就是不要求法向量单位的情况下的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

这就是我们熟悉的Gauss-Schmidt正交化过程

此时，v1Tv1\mathbf v_1^T\mathbf v_1实际可以只计算一次，每次矩阵乘向量的实际计算复杂度变为了向量内积和向量线性运算

这样进行分解的总复杂度为O(mn2)O(mn^2)

Givens旋转变换进行QR分解

正交矩阵有两种形态：反射和旋转，前面研究的是利用反射变换进行QR分解的思路，下面我们讨论旋转变换

定义（二阶Givens矩阵） 二阶Givens矩阵形态如下：

G=[cos⁡θsin⁡θ−sin⁡θcos⁡θ]\mathbf G=\begin{bmatrix}\cos\theta&\sin\theta\\-\sin\theta&\cos\theta\end{bmatrix}

定义（n阶Givens矩阵） nn阶Givens矩阵为二阶矩阵嵌入nn阶单位阵的结果，例子如下：

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=[1cs1−sc1], c2+s2=1\mathbf G=\begin{bmatrix}1\\&c&&s\\&&1\\&-s&&c\\&&&&1\end{bmatrix},\ c^2+s^2=1

我们可以通过选择cc和ss的值将一个二维向量化为基向量的形状：

G[x1x2]=[α0]\mathbf G\begin{bmatrix}x_1\\x_2\end{bmatrix}=\begin{bmatrix}\alpha\\0\end{bmatrix}

只需令

c=x1α, s=x2αc=\frac{x_1}\alpha,\ s=\frac{x_2}\alpha

注意，此时必须满足

α2=x12+x22\alpha^2=x_1^2+x_2^2

与Householder变换不同，Givens变换只影响向量的两行，因此，如果要将一个nn维稠密向量进行消去，需要n−1n-1个Givens变换，分别在第11维和第k=2,⋯ ,nk=2,\cdots,n维构成的平面上旋转；而Householder变换只需要1个

因此，虽然单个Givens变换的计算量更小，但其对稠密矩阵的计算量很大。Givens变换进行QR分解更适合稀疏矩阵

求矩阵的所有特征值

其思路是利用正交相似变换（不改变特征值）将矩阵转化为三角阵或分块三角阵（容易求出特征值）

思路：收缩技术

假设我们已知矩阵A\mathbf A的一个特征向量x\mathbf x（这可以通过幂法得到）

那么我们可以构造一个正交变换H\mathbf H使得Hx=σe1\mathbf {Hx}=\sigma\mathbf e_1，则

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

于是e1\mathbf e_1是正交相似于A\mathbf A的矩阵HAHT\mathbf {HAH}^T的同一特征值的特征向量

于是我们可以得到

HAHT=[λr1TA2]\mathbf {HAH}^T=\begin{bmatrix}\lambda&\mathbf r_1^T\\&\mathbf A_2\end{bmatrix}

此时问题就变小了

然而，利用收缩技术和幂法来求解，效率不高，而且误差会累积

QR迭代算法

定理（实Schur分解） 对于任何A∈Rn×n\mathbf A\in\mathbb R^{n\times n}，存在正交阵Q\mathbf Q使得

QTAQ=S\mathbf Q^T\mathbf {AQ}=\mathbf S

其中S\mathbf S是拟上三角阵，即对角块不超过2阶的分块上三角阵

注意：在复数域下，Schur分解可以得到真上三角阵，但因为实矩阵的特征值可能有复数，因此无法得到全实的真Schur分解，但由于实矩阵的复特征值一定共轭成对出现，因此可以将成对出现的复特征值包裹在一个2阶实对角块中

QR迭代算法

  • 对Ak\mathbf A_k做QR分解：Ak=QkRk\mathbf A_k=\mathbf Q_k\mathbf R_k
  • 作迭代：Ak+1=RkQk\mathbf A_{k+1}=\mathbf R_k\mathbf Q_k

这即Ak+1=QkTAQk\mathbf A_{k+1}=\mathbf Q_k^T\mathbf {AQ}_k，构成了一组正交相似矩阵序列

定理（QR迭代算法的收敛定理） 若矩阵A\mathbf A的所有等模特征值组或者是实重特征值，或者是复共轭特征值，则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}

可以针对左下角块构造一个Householder变换使得

H1′c1=σ1e1\mathbf H_1'\mathbf c_1=\sigma_1\mathbf e_1

令

H1=[10T0H1′]\mathbf H_1=\begin{bmatrix}1&\mathbf 0^T\\\mathbf 0&\mathbf H_1'\end{bmatrix}

于是

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}

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}

由于H1\mathbf H_1是对称矩阵，上述结果就是对A\mathbf A进行一步正交相似变换得到的结果

我们可以对其右下角块迭代上述正交相似变换，从而将A\mathbf A正交相似于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=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}

对于位移因子sks_k的选取，可以使用Rayleigh策略：取sks_k为Ak\mathbf A_k的(n,n)(n,n)元

矩阵的奇异值分解

基本概念

奇异值和奇异向量是针对非方阵的，特征值和特征向量的扩展定义

定义（奇异值、奇异向量） 对矩阵A∈Rm×n\mathbf A\in\mathbb R^{m\times 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}

则称σ\sigma为A\mathbf A的奇异值，u,v\mathbf u,\mathbf v分别为对应的左奇异向量和右奇异向量

定理（奇异值分解, SVD） 对任何矩阵A∈Rm×n, m≥n\mathbf A\in\mathbb R^{m\times n},\ m\geq n，可以找到正交矩阵U,V\mathbf U,\mathbf V和类对角矩阵Σ∈Rm×n\mathbf \Sigma\in\mathbb R^{m\times n}，使得

A=UΣVT\mathbf A=\mathbf {U\Sigma V}^T

其中

Σ=[σ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

讨论 上述分解等价于

AV=UΣ\mathbf {AV}=\mathbf{U\Sigma}

即

[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}

这就是奇异向量的第一部分定义，而第二部分定义则来自于

UTA=ΣVT\mathbf U^T\mathbf A=\mathbf \Sigma\mathbf V^T

于是，奇异值分解存在等价于说，矩阵属于不同奇异值的奇异向量正交

讨论 这里的m≥nm\geq n不是必须的，只是为了叙述方便，如果m<nm<n仍然可以得到相同形式的分解，不过Σ\Sigma的零部分在右侧

SVD的证明

对ATA\mathbf A^T\mathbf A作谱分解：

ATA=VΛVT\mathbf A^T\mathbf A=\mathbf{V\Lambda V}^T

令

Λ=ΣTΣ\mathbf \Lambda=\mathbf \Sigma^T\mathbf \Sigma

于是

ATA=VΣTΣVT\mathbf A^T\mathbf A=\mathbf {V\Sigma}^T\mathbf{\Sigma V}^T

假设Σ\mathbf\Sigma的前rr个对角元非零，则取出这部分的对角矩阵作为Σr\mathbf\Sigma_r，而V\mathbf V的前rr列作为Vr\mathbf V_r

于是令

Ur=AVrΣr−1\mathbf U_r=\mathbf{AV}_r\mathbf\Sigma_r^{-1}

（这里A,Vr,Σr\mathbf A,\mathbf V_r,\mathbf\Sigma_r大小分别为m×nm\times n，n×rn\times r，r×rr\times r，得到m×rm\times r的Ur\mathbf U_r）

由于

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

可知Ur\mathbf U_r列正交，将其扩充至m×mm\times m即可得到U\mathbf U

简化的SVD

上面的证明中，我们选择了非零奇异值及其对应的奇异向量，作为如下形式的分解：

A=UrΣrVrT\mathbf A=\mathbf U_r\mathbf \Sigma_r\mathbf V_r^T

其中零奇异值对应的部分对乘积结果毫无影响，因此这就是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

这是矩阵的一种秩一分解，从此也可以看出，非零奇异值的个数rr对应矩阵的秩

SVD的几何意义

我们回顾一下谱分解的几何意义：

A=XΛX−1\mathbf A=\mathbf{X\Lambda X}^{-1}

假设我们有一个向量

v=[v1⋮vn]\mathbf v=\begin{bmatrix}v_1\\\vdots\\v_n\end{bmatrix}

现在我们找到了一组基

X=[x1⋯xn]\mathbf X=\begin{bmatrix}\mathbf x_1&\cdots&\mathbf x_n\end{bmatrix}

于是，可以将v\mathbf v用这组基表示：

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}

于是

[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}

即A=XΛX−1\mathbf A=\mathbf{X\Lambda X}^{-1}中的第一步X−1\mathbf X^{-1}是将输入向量从标准基上的坐标转化为X\mathbf X定义的基上的坐标

而根据特征向量的意义，A\mathbf A对特征向量的操作仅仅是伸缩变换，因此A\mathbf A在X\mathbf X定义的基上的形式就是对角矩阵Λ\mathbf\Lambda

而最后一步X\mathbf X，则是第一步的逆过程，将输出向量从X\mathbf X基上的坐标转换回原始坐标

于是，谱分解的几何意义是将一个映射A\mathbf A拆分成了以下三步：

  • 第一步：X−1\mathbf X^{-1}，将输入向量从标准基换到特征基
  • 第二步：Λ\mathbf\Lambda，将特征基下的输入向量变换得到输出向量
  • 第三步：X\mathbf X，将输出向量从特征基换回标准基

对于实对称矩阵而言，我们有更好的谱分解形式：

A=QΛQT\mathbf A=\mathbf{Q\Lambda Q}^T

此形式与一般谱分解的区别是，将特征基取成了一组标准正交基，这是实对称矩阵才具有的独特性质

接下来我们回到SVD。SVD相比谱分解而言，不要求矩阵是方阵，这意味着，矩阵的定义域即输入空间，和陪域即输出空间，从本质上不再是同一个空间，于是，我们也不再能够找到一组同时对输入空间和输出空间有效的特征基

此时，我们就需要对输入空间和输出空间分别取一组基（可以分别称为右奇异基和左奇异基），SVD就意味着将原始矩阵的映射进行如下拆分：

  • 第一步：VT\mathbf V^T，将输入向量从标准基换到右奇异基
  • 第二步：Σ\mathbf\Sigma，将右奇异基下的输入向量变换，得到左奇异基下的输出向量
  • 第三步：U\mathbf U，将输出向量从左奇异基换到标准基

这样看起来，这个右奇异基和左奇异基的选取似乎就非常任意了，但注意，Σ\mathbf\Sigma是拟对角阵，这意味着，从输入空间（右奇异基）到输出空间（左奇异基）的映射，其实是沿着坐标轴的伸缩映射，但由于输入输出空间维数可能不相等，这个过程中可能会舍弃一些维度，或增加一些值为00的维度

接下来，我们将U\mathbf U和V\mathbf V分别拆分成前rr列和后若干列：

U=[UrUm−r]\mathbf U=\begin{bmatrix}\mathbf U_r&\mathbf U_{m-r}\end{bmatrix}

V=[VrVn−r]\mathbf V=\begin{bmatrix}\mathbf V_r&\mathbf V_{n-r}\end{bmatrix}

根据Σ\mathbf\Sigma的非零值只有前rr个对角元，我们可以知道这四个子矩阵的列空间与映射A\mathbf A的一些关系：

  • R(Vn−r)\mathcal R(\mathbf V_{n-r})经A\mathbf A映射后为0\mathbf 0
  • 所有向量经A\mathbf A映射后必然落入R(Ur)\mathcal R(\mathbf U_r)

于是，我们可以得到SVD与四个基本子空间的关系：

  • 列空间：R(Ur)\mathcal R(\mathbf U_r)，列空间的基是U\mathbf U的前rr列
  • 零空间：R(Vn−r)\mathcal R(\mathbf V_{n-r})，零空间的基是V\mathbf V的后n−rn-r列
  • 行空间：零空间的正交补，即V\mathbf V的前rr列
  • 左零空间：列空间的正交补，即U\mathbf U的后m−rm-r列

此时，我们就可以从简化SVD的角度再次讨论几何意义：

  • 第一步：VrT\mathbf V_r^T，将输入向量中属于零空间的部分舍弃，属于零空间正交补（即行空间）的部分映射到一个rr维空间中，得到的结果实际上就是输入向量在行空间的右奇异基下的坐标表示
  • 第二步：Σr\mathbf \Sigma_r，对行空间下用右奇异基表示的输入向量进行对角矩阵变换，得到列空间下用左奇异基表示的输出向量
  • 第三步：Ur\mathbf U_r，将列空间下左奇异基表示的输出向量转化为陪域下标准基表示的输出向量，此时不属于列空间的那部分固定为00

从这个角度，每个向量在简化SVD的视角其实经历了以下四个空间：

image-20250425162003808

此时我们也可以考虑矩阵的广义逆：

A∗=VrΣr−1UrT\mathbf A^*=\mathbf {V}_r\mathbf{\Sigma}_r^{-1}\mathbf U_r^T

显然，它将上述空间变换的过程取了逆序，同时将行空间到列空间的对角变换也取了逆。广义逆在行空间到列空间的作用上与原矩阵是互逆的，而对零空间和左零空间中的内容会舍弃

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

对于r0<rr_0<r，如果我们希望找到一个矩阵Ar0\mathbf A_{r_0}，使得它在所有秩为r0r_0的矩阵中“最接近”A\mathbf A，我们就可以选择A\mathbf A秩一分解中的前r0r_0项：

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

由于所有奇异值是按从大到小排序的，上面取到的就是最佳的低秩逼近结果，可以证明，它与原矩阵的误差的范数（在2-范数或F-范数下）是最小的