QR分解.

【引理】设 x⃗=(x1,x2,⋯,xi,⋯,xj,⋯,xn)T\vec x=(x_1,x_2,\cdots,x_i,\cdots,x_j,\cdots,x_n)^Tx=(x1,x2,⋯,xi,⋯,xj,⋯,xn)T，其中 xi⋅xj≠0x_i·x_j\neq0xi⋅xj=0，则存在一个Givens变换矩阵 PijP_{ij}Pij 使得Pijx⃗=(y1,y2,⋯,yi,⋯,yj,⋯,yn)TP_{ij}\vec x=(y_1,y_2,\cdots,y_i,\cdots,y_j,\cdots,y_n)^TPijx=(y1,y2,⋯,yi,⋯,yj,⋯,yn)T其中yi=xi2+xj2y_i=\sqrt{x_i^2+x_j^2}yi=xi2+xj2yj=0y_j=0yj=0
上述Givens矩阵 PijP_{ij}Pij 中P[i,i]=cos⁡θ,P[i,j]=sin⁡θP[i,i]=\cos\theta~,~P[i,j]=\sin\thetaP[i,i]=cosθ , P[i,j]=sinθP[j,i]=−sin⁡θ,P[j,j]=cos⁡θP[j,i]=-\sin\theta~,~P_[j,j]=\cos\thetaP[j,i]=−sinθ , P[j,j]=cosθcos⁡θ=xixi2+xj2,sin⁡θ=xjxi2+xj2\cos\theta=\frac{x_i}{\sqrt{x_i^2+x_j^2}}~,~\sin\theta=\frac{x_j}{\sqrt{x_i^2+x_j^2}}cosθ=xi2+xj2xi , sinθ=xi2+xj2xj
实际操作时，可以通过对 x⃗\vec xx 的规范化操作来避免数值溢出。

【定理】对于 nnn 阶非奇异矩阵 AAA 而言，可以分解为正交阵 QQQ 和上三角阵 RRR 的乘积，即 A=QR.A=QR.A=QR. 该分解是唯一的，当矩阵 RRR 中的元素满足∀i∈[1,n],R[i,i]>0\forall~i\in[1,n],R[i,i]>0∀ i∈[1,n],R[i,i]>0

【定理】对于 nnn 阶非奇异矩阵 AAA 而言，存在一系列平面旋转变换矩阵 Pk∣k=1,2,⋯,n−1P_k|k=1,2,\cdots,n-1Pk∣k=1,2,⋯,n−1，使得 Pn−1⋯P2P1A=RP_{n-1}\cdots P_2P_1A=RPn−1⋯P2P1A=R其中 RRR 是上三角阵且 rii>0,i=1,2,⋯,n.r_{ii}>0,i=1,2,\cdots,n.rii>0,i=1,2,⋯,n.
【证明】由于 AAA 非奇异，所以不可能存在全零列，因此对于第一列，可以做出一系列平面旋转变换 P1=∏P1kP_1=\prod P_{1k}P1=∏P1k使得第一列约化为 (r11,0,⋯,0)T(r_{11},0,\cdots,0)^T(r11,0,⋯,0)T，同理对于后面的第 mmm列，可以做一系列平面旋转变换 PkP_kPk 使得该列约化为 (r1m,r2m,⋯,rmm,0,⋯,0)T(r_{1m},r_{2m},\cdots,r_{mm},0,\cdots,0)^T(r1m,r2m,⋯,rmm,0,⋯,0)T，最终得到上三角阵，并引入对角阵 D={±1,±1,⋯,±1}D=\{\pm1,\pm1,\cdots,\pm1\}D={±1,±1,⋯,±1} 使得主对角元全正。
显然如果令 QT=∏i=1n−1PiQ^T=\prod^{n-1}_{i=1}P_iQT=i=1∏n−1Pi那么QTA=RQ^TA=RQTA=R因为 QQQ 是正交阵，所以 QT=Q−1Q^T=Q^{-1}QT=Q−1，所以有A=(QT)−1R=QRA=\Big(Q^T\Big)^{-1}R=QRA=(QT)−1R=QR
具体地如何进行QR分解，得到正交阵 QQQ 和上三角阵 RRR，主流的有三种方案：Schmidt正交化、Givens变换和Householder变换，在第四部分对前两者有简介。

QR迭代算法.

设 AAA 是 nnn 阶非奇异阵，记 A1=AA_1=AA1=A，对其进行QR分解得到 A1=Q1R1A_1=Q_1R_1A1=Q1R1按照如下公式进行迭代：Ak+1=RkQk=QkTAkQk=Qk+1Rk+1(*)A_{k+1}=R_kQ_k=Q_k^TA_kQ_k=Q_{k+1}R_{k+1}\tag{*}Ak+1=RkQk=QkTAkQk=Qk+1Rk+1(*)
上述利用矩阵QR分解和递推式 (∗)(*)(∗) 构造矩阵序列 {Ak}\{A_k\}{Ak} 的过程即为QR算法。

Ak+1=QkTAkQkA_{k+1}=Q^T_{k}A_kQ_kAk+1=QkTAkQk
Ak+1=(Q1Q2⋯Qk)TA1(Q1Q2⋯Qk)A_{k+1}=(Q_1Q_2\cdots Q_k)^TA_1(Q_1Q_2\cdots Q_k)Ak+1=(Q1Q2⋯Qk)TA1(Q1Q2⋯Qk)
若记 Qk~=∏i=1kQi,Rk~=∏i=k1Ri\tilde{Q_k}=\prod^k_{i=1}Q_i~,~\tilde{R_k}=\prod^1_{i=k}R_iQk~=i=1∏kQi , Rk~=i=k∏1Ri那么矩阵 AkA^kAk 的QR分解式为 Ak=Qk~Rk~A^k=\tilde{Q_k}\tilde{R_k}Ak=Qk~Rk~
【证明】数学归纳法
①当 k=1k=1k=1 时显然有 A=Q1R1=Q1~R1~A=Q_1R_1=\tilde{Q_1}\tilde{R_1}A=Q1R1=Q1~R1~
②设 Ak−1A^{k-1}Ak−1 分解式为 Ak−1=Qk−1~Rk−1~A^{k-1}=\tilde{Q_{k-1}}\tilde{R_{k-1}}Ak−1=Qk−1~Rk−1~
③那么 Qk~Rk~=Q1Q2⋯Qk−1(QkRk)Rk−1⋯R2R1=Q1Q2⋯Qk−1AkRk−1⋯R2R1=Qk−1~AkRk−1~=Qk−1~(Qk−1~)TAQk−1~Rk−1~=AAk−1=Ak\begin{aligned}\tilde{Q_k}\tilde{R_k} &=Q_1Q_2\cdots Q_{k-1}(Q_kR_k)R_{k-1}\cdots R_2R_1\\ &=Q_1Q_2\cdots Q_{k-1}A_kR_{k-1}\cdots R_2R_1\\ &=\tilde{Q_{k-1}}A_k\tilde{R_{k-1}}\\ &=\tilde{Q_{k-1}}\Big(\tilde{Q_{k-1}}\Big)^TA\tilde{Q_{k-1}}\tilde{R_{k-1}}\\ &=AA^{k-1}\\ &=A^k \end{aligned}Qk~Rk~=Q1Q2⋯Qk−1(QkRk)Rk−1⋯R2R1=Q1Q2⋯Qk−1AkRk−1⋯R2R1=Qk−1~AkRk−1~=Qk−1~(Qk−1~)TAQk−1~Rk−1~=AAk−1=Ak

收敛性.

对于 nnn 阶非奇异阵 AAA，若 AAA 的特征值满足 ①∣λ1∣>∣λ2∣>⋯>∣λn∣①~|\lambda_1|>|\lambda_2|>\cdots>|\lambda_n|① ∣λ1∣>∣λ2∣>⋯>∣λn∣；②A②~A② A 有标准型 A=XDX−1A=XDX^{-1}A=XDX−1，其中 D=diagλ1,λ2,⋯,λnD=diag{\lambda_1,\lambda_2,\cdots,\lambda_n}D=diagλ1,λ2,⋯,λn，且 X−1X^{-1}X−1 有三角分解 X−1=LUX^{-1}=LUX−1=LU，则由QR算法产生的矩阵序列 {Ak}\{A_k\}{Ak} 有下列极限：lim⁡k→∞Ak=R\lim_{k\rightarrow\infin}A_k=Rk→∞limAk=R 其中 RRR 为上三角阵，并且 rii=λi.r_{ii}=\lambda_i.rii=λi.
【证明】对于非奇异矩阵 AAA，存在可逆矩阵 XXX 使得 X−1AX=DX^{-1}AX=DX−1AX=D，即 A=XDX−1A=XDX^{-1}A=XDX−1，因此有Ak=XDkX−1(1)A^k=XD^kX^{-1}\tag{1}Ak=XDkX−1(1)条件中认为 X−1X^{-1}X−1 存在LU分解 X−1=LUX^{-1}=LUX−1=LU，因此 (1)(1)(1) 式可写为 Ak=XDkLU=X(DkLD−k)DkU(2)A^k=XD^kLU=X(D^kLD^{-k})D^{k}U\tag{2}Ak=XDkLU=X(DkLD−k)DkU(2)将 (2)(2)(2) 中 DkLD−kD^kLD^{-k}DkLD−k 做加性分解，得到DkLD−k=E+ΛkD^kLD^{-k}=E+\Lambda_kDkLD−k=E+Λk其中 Λk\Lambda_kΛk 是一个下三角阵，且其主对角元 Λ[i,i]=0\Lambda[i,i]=0Λ[i,i]=0，非零元素表达式为Λ[i,j]=L[i,j]⋅(λiλj)k,i>j\Lambda[i,j]=L[i,j]·\Big(\frac{\lambda_i}{\lambda_j}\Big)^k~,~i>jΛ[i,j]=L[i,j]⋅(λjλi)k , i>j由于特征值按绝对值降序排列，因此有下面的极限存在lim⁡k→∞Λ=0\lim_{k\rightarrow\infin}\Lambda=0k→∞limΛ=0这里的 000 代表零矩阵。由于矩阵 XXX 非奇异，对其进行QR分解，有X=QRX=QRX=QR至此 (2)(2)(2) 式可以写为Ak=QR(E+Λk)DkU=Q(E+RΛkR−1)RDkU(3)A^k=QR(E+\Lambda_k)D^{k}U=Q(E+R\Lambda_kR^{-1})RD^kU\tag{3}Ak=QR(E+Λk)DkU=Q(E+RΛkR−1)RDkU(3)当 k→∞k\rightarrow\infink→∞ 时，矩阵 R(E+Λk)R(E+\Lambda_k)R(E+Λk) 可逆，从而 R(E+Λ)R−1=E+RΛkR−1R(E+\Lambda)R^{-1}=E+R\Lambda_kR^{-1}R(E+Λ)R−1=E+RΛkR−1 可逆，对其进行QR分解，有E+RΛkR−1=Qk′Rk′E+R\Lambda_kR^{-1}=Q_k'R_k'E+RΛkR−1=Qk′Rk′从而 (3)(3)(3) 式写为Ak=(QQk′)(Rk′RDkU)(4)A^k=(QQ_k')(R_k'RD^kU)\tag{4}Ak=(QQk′)(Rk′RDkU)(4)显然 (4)(4)(4) 式已经是 AAA 的QR分解，为确保唯一性，引入对角阵D∗=diag{±1,±1,⋯,±1}D^*=diag\{\pm1,\pm1,\cdots,\pm1\}D∗=diag{±1,±1,⋯,±1}将 (4)(4)(4) 写为Ak=(QQk′D∗)(D∗Rk′RDkU)(5)A^k=(QQ_k'D^*)(D^*R_k'RD^kU)\tag{5}Ak=(QQk′D∗)(D∗Rk′RDkU)(5)可得Qk~=QQk′D∗\tilde{Q_k}=QQ_k'D^*Qk~=QQk′D∗Rk~=D∗Rk′RDkU\tilde{R_k}=D^*R_k'RD^kURk~=D∗Rk′RDkU因此 Ak+1=(Q1Q2⋯Qk)TA1(Q1Q2⋯Qk)=(Qk~)TAQk~=D∗(Qk′)TQTAQQk′D∗(6)A_{k+1}=(Q_1Q_2\cdots Q_k)^TA_1(Q_1Q_2\cdots Q_k)=\Big(\tilde{Q_k}\Big)^TA\tilde{Q_k}=D^*\Big(Q_k'\Big)^TQ^TAQQ_k'D^*\tag{6}Ak+1=(Q1Q2⋯Qk)TA1(Q1Q2⋯Qk)=(Qk~)TAQk~=D∗(Qk′)TQTAQQk′D∗(6)因为 A=XDX−1,X=QRA=XDX^{-1},X=QRA=XDX−1,X=QR，所以 QTAQ=RDR−1Q^TAQ=RDR^{-1}QTAQ=RDR−1，因此 (6)(6)(6) 式可以写为Ak+1=D∗(Qk′)T(RDR−1)Qk′D∗(7)A_{k+1}=D^*\Big(Q_k'\Big)^T(RDR^{-1})Q_k'D^*\tag{7}Ak+1=D∗(Qk′)T(RDR−1)Qk′D∗(7)记 R0=RDR−1,gk=Qk′D∗R_0=RDR^{-1},g_k=Q_k'D^*R0=RDR−1,gk=Qk′D∗R0[i,i]=λiR_0[i,i]=\lambda_iR0[i,i]=λi所以 (7)(7)(7) 式可以写为Ak+1=gkTR0gk(8)A_{k+1}=g_k^TR_0g_k\tag{8}Ak+1=gkTR0gk(8)由于(E+Λk→E)⇒(R(E+Λ)R−1=E+RΛkR−1→E)\Big(E+\Lambda_k\rightarrow E\Big)\Rightarrow \Big(R(E+\Lambda)R^{-1}=E+R\Lambda_kR^{-1}\rightarrow E\Big)(E+Λk→E)⇒(R(E+Λ)R−1=E+RΛkR−1→E)所以 Qk′→EQ'_k\rightarrow EQk′→E，所以gk→D∗=diag{±1,±1,⋯,±1}g_k\rightarrow D^*=diag\{\pm1,\pm1,\cdots,\pm1\}gk→D∗=diag{±1,±1,⋯,±1}所以 k→∞k\rightarrow\infink→∞ 时，lim⁡Ak+1[i,i]=λi,i=1,2,⋯,n\lim A_{k+1}[i,i]=\lambda_i~,~i=1,2,\cdots,nlimAk+1[i,i]=λi , i=1,2,⋯,n其收敛速度由 ∣∣Λk∣∣∞||\Lambda_k||_{\infin}∣∣Λk∣∣∞ 决定，有上限∣∣Λk∣∣∞≤C⋅max⁡j∈[1,n−1]∣λj+1λj∣k||\Lambda_k||_{\infin}≤C·\max_{j\in[1,n-1]}\Big|\frac{\lambda_{j+1}}{\lambda_j}\Big|^k∣∣Λk∣∣∞≤C⋅j∈[1,n−1]max∣∣∣λjλj+1∣∣∣k
上述收敛速度定义式 rn=∣λn/λn−1∣r_n=|\lambda_n/\lambda_{n-1}|rn=∣λn/λn−1∣ 与幂法中的收敛速度定义式相似度极高，当 rn≈1r_n\approx1rn≈1 时收敛极慢，可以效仿幂法中采取原点位移策略，即对于矩阵 A−s⋅EA-s·EA−s⋅E 使用QR分解，那么收敛速度即为 rn∗=∣λn−sλn−1−s∣r_n^*=\Big|\frac{\lambda_n-s}{\lambda_{n-1}-s}\Big|rn∗=∣∣∣λn−1−sλn−s∣∣∣上述算法称为带原点位移的QR迭代算法。

【Reference】UTEP-Math 5330.

逻辑极度清晰，可惜我弄丢了上述Lecture Note的网址，关于Householder变换可以参考《Householder变换》。
另外注意到Note中讨论了 a11=a112+a212+⋯+an12=0a_{11}=\sqrt{a_{11}^2+a_{21}^2+\cdots+a_{n1}^2}=0a11=a112+a212+⋯+an12=0 的情况，不难发现上述等式成立意味着 ai1=0,i=1,2,⋯,na_{i1}=0,i=1,2,\cdots,nai1=0,i=1,2,⋯,n，那么 AAA 是奇异矩阵，无法进行QR分解。
第一个重点是 Algorithm3.1 的伪代码，一般情况下我们假定 AAA 是 nnn 阶非奇异阵，即 m=nm=nm=n，伪代码中 β\betaβ 是当次迭代列向量的模值，γ\gammaγ 是向量 vvv 的分母，后续计算出向量 vvv 以及对应的豪氏矩阵 Hv=E−2vvTH_v=E-2vv^THv=E−2vvT与矩阵 AAA 进行左乘 A1=HvA0A_1=H_vA_0A1=HvA0，变量 QQQ 初值是 EmE_mEm 单位阵，迭代过程中对 HvH_vHv 左乘，所以最终停止时 Q=∏i=1n−1HiQ=\prod_{i=1}^{n-1}H_iQ=i=1∏n−1Hi即为QR分解中 A=QRA=QRA=QR 的【QQQ】，而此时保存 AAA 的单元被【RRR】所覆盖。
11-13行代码作用是在矩阵本身已是上三角化时直接进行下一次迭代。这篇Note的一个特点是，开篇已经假定 AAA 是非奇异阵，但后续的分析和伪码又是针对 A∈Rm×nA\in R^{~m\times n}A∈R m×n 考虑的，个人猜测是出于算法健壮性。
Note的最后一部分给出了基于Householder变化进行QR分解复杂度分析：
上半部分关于 AAA 阵迭代无需显式计算豪氏矩阵以及 QQQ 阵迭代的浮点运算 flops 挺靠谱的，后半部分关于整个算法的复杂度分析持保留态度，尤其是最后给出的 23r3\frac23r^332r3，这里不妨认为 r=min⁡m=n(m,n)=nr=\min_{m=n}(m,n)=nr=m=nmin(m,n)=n那么复杂度为 23r3\frac23r^332r3，但这个已经出现 4p24p^24p2 项，求和后应为 43r3\frac43r^334r3，个人倾向于这个答案，大致的 O(n3)O(n^3)O(n3) 是无疑的。
给出结果为 43n3\frac43n^334n3 的文献如下，来自UTexas(这俩名字应该是一个学校吧)：

Summary.

需要明确的是，第四部分内容是对于QR分解以及基于 Householder变换实现QR分解的，并不是用于求解矩阵全部特征值的QR迭代法。
【存疑】QR迭代法思路之一就是在迭代公式中使用Householder变换实现每一次迭代的QR分解，使用Householder变换进行QR分解的时间复杂度为 23n3\frac23n^332n3，大约是第一部分构造性证明中使用Givens变换的 23\frac2332 倍。
另一种QR迭代法的思路是首先使用Householder变换，将一般实矩阵(实对称矩阵)约化为上海森伯格矩阵(对称三对角矩阵)，而后再进行QR迭代，这样做的优势在于可以使用对于稀疏矩阵较为高效的Givens变换阵(常数时间即可计算出)实现QR分解，而无需使用需要 n2n^2n2 才能计算出的豪斯荷尔德矩阵。

【亟待解决】直接使用Householder变换进行的QR分解时间复杂度，有文献表明是 43n3\frac43n^334n3，也有文献认为是 23n3\frac23n^332n3，个人倾向于前一种说法，理由在上面。
另外的两种QR分解方法，格拉姆-施密特正交化GS(以及修正格拉姆-施密特正交化MGS)和吉文斯变换，其复杂度也为 O(n3)O(n^3)O(n3)，有说法认为前者的常数因子为 111，后者为 43.\frac43.34.
使用豪斯荷尔德变换将实矩阵约化为上Hessenberg矩阵需要的乘法次数：53n3.\frac53n^3.35n3.
那么下图中的说法就不能成立：

【NA】基于QR分解的特征值迭代法相关推荐

基于QR分解迭代求解方阵特征值和特征向量
基于QR分解迭代求解方阵特征值和特征向量一.特征值与特征向量求解的难点线性代数的知识告诉我们如果要求一个方阵的特征值,只需要求解如下的特征方程的根即可: f ( λ ) = ( λ − λ 1 ) ...
Julia 矩阵QR分解和特征值
Julia 矩阵QR分解和特征值前言 1. 施密特正交 (1) 利用施密特正交求出正交矩阵Q (2) 求出上三角矩阵R (3) 改进的消减QR分解 2. 完全QR分解 3. 矩阵QR分解的作用 (1 ...
双步位移求解特征值matlab,数值分析——带双步位移的QR分解求特征值算法
C语言实现数值分析中带双步位移的QR分解求特征值算法. 数值分析(B) 大作业(二) 1.算法设计: ①矩阵的拟上三角化: 对实矩阵A进行相似变换化为拟上三角矩阵A(n 1),其变换矩阵采用 ...
人工智能里的数学修炼 | 矩阵的花样分解：特征值分解（EVD）、相似对角化、QR分解、Schur分解、奇异值分解（SVD）的概念纠缠与详解
前言在高等代数里,矩阵分解是一个十分基础与重要的内容,任何一个学校对于理工科的研究生教育都会开设相应的课程,如:矩阵分析.矩阵论.线性系统等.看了不少社区的问答.笔记和博客,在它们的基础上加入一些自 ...
数值分析中的QR分解及其代码实现
QR分解若 A∈Cm×k A ∈ C m × k A\in C^{m×k} 是一个列满秩的矩阵,rank(A) = k,则矩阵A 可以分解为 A=QR A = Q R A = QR , Q∈Cm×k ...
5.12 QR分解的阻尼倒数法和正则化方法区别
5.12 QR分解的阻尼倒数法和正则化方法区别基于QR分解的阻尼倒数法目的是改善矩阵 AAA 的病态,正则化方法目的也是改善矩阵 AAA 的病态.区别是阻尼倒数法能获得稀疏解,即最优解分量趋近 00 ...
斯密特正交化进行QR分解
矩阵类: #include<iostream> #include <stdlib.h> #include <math.h> #define iter 60 //迭代 ...
基于Householder变换的QR分解
01.function [Q,R]=qrhs(A)02.% 基于Householder变换,将<strong><span style="color:#993399;&quo ...
基于Givens变换的QR分解
01.function [Q,R]=qrgv(A)02.% 基于Givens变换,将方阵A分解为A=QR,其中Q为正交矩阵,R为上三角阵03.%04.% 参数说明05.% A:需要进行QR分解的方阵0 ...

【NA】基于QR分解的特征值迭代法

文章目录

QR分解.

QR迭代算法.

收敛性.

【Reference】UTEP-Math 5330.

Summary.

【NA】基于QR分解的特征值迭代法相关推荐

最新文章

热门文章