2.6 QR分解和QR算法
为一个给定的矩阵 A∈Mn 提供一种 Schur 酉上三角化 (2.3.1) 的特殊计算方法, 以及 (在某些假设下) 计算特征值的一个通用的数值方法, 称为 QR 算法, QR 算法基于一般矩阵 A∈Mn,m 的所谓 QR 分解.
2.6.1 定理 (QR 分解) 如果 A∈Mn,m 且 n⩾m , 那么存在有标准正交列的矩阵 Q∈Mn,m 和上三角矩阵 R∈Mm , 使得 A−QR . 如果 n=m , 那么 Q 是酉矩阵; 此外, 如果 A 是非奇异矩阵, 则可以选取 R 为具有正对角元的上三角矩阵, 并且在这种情形, 因子 Q 和 R 都是唯一的. 如果 A∈Mn,m(R) , 那么 Q 和 R 都可以取实矩阵.
证明:如果 A∈Mn,m ,且 rankA=m ,则 A 的各列构成 Cn 的一个无关组,把 Gram-Schmidt 过程(0.6.4)应用于 A 的各列,用矩阵记号描述所得结果正是 A 的 QR 分解。Gram
Schmidt算法的自然推广使同样的矩阵记号描述能够应用于 A 的各列可能是相关的一般情形. 设 A=[a1⋯am] 写成了以它的列 ai∈Cn 为子块的分块形式. 如果 a1=0 , 令 q1=0 ; 否则令 q1=a1/(a1∗a1)1/2 . 对每个 k=2,3,⋯,m , 如同在通常的 Gram-Schmidt 过程中那样, 算出
yk=ak−i=1∑k−1(qi∗ak)qi. 如果 yk−0 (它成立,当且仅当 ak 是 a1,a2,…,ak−1 的线性组合),则令 qk=0 ;否则令 qk=yk/(yk∗yk)1/2 ,这时,向量 q1,…,qm 是正交组,其中的每个元素或者是单位向量,或者是零向量.每个向量 qj 是 a1,…,aj 的线性组合,反过来,上述做法保证每个列 aj 是 q1,…,qi 的线性组合.因此存在纯量 rkj 使得
aj=k=1∑jrkjqk,j=1,2,…,m.(2.6.2) 如果对所有 k>j ,令 rkj=0 ,而对于所有使 qi=0 的每个 i ,对所有 j=1,2,…,m ,令 rij=0 ,那么,经过上述过程,上三角矩阵 R=[rij]∈Mm 和向量 q1,q2,…,qm 由 a1,a2,…,am 所确定。矩阵 Q≡[q1,…,qm]∈Mm,m 有正交列(其中有些可能是零),因而(2.6.2)表明 A=QR .
如果 rankA=m ,则 Q 有标准正交列,因而所要求的分解形式已经得到。特别地,如果 m=n ,且 A 是非奇异矩阵,那么,根据(2.1.4e), Q 必须是酉矩阵,而非奇异矩阵 R=Q∗A 的诸对角元必须是非零的。在这种情形,因为要求 R 是上三角矩阵,故 q1 是 a1 的纯量倍数,并且对于 i=2,⋯,n,qi 位于这样一个二维空间内,它在由 a1,⋯,ai 生成的空间中是 a1,⋯,ai ,生成空间的正交补。因此,除了差一个绝对值为 1 的比例因子以外,每个 qi 是唯一确定的。于是用 R′≡diag(∣r11∣/r11,⋯,∣rmn∣/rmn)R 代替 R 以及用 Q′≡Qdiag(r11/∣r11∣,⋯,rmn/∣rmn∣) 代替 Q ,便给出唯一的分解 A=Q′R′ ,这正是定理的叙述中所断言的。
如果 A 的各列不是无关的,取由 Q 的非零列组成的(标准正交)组,然后把它扩充为 C 的一个标准正交基;用上述方法得来的新向量记作 z1,z2,⋯,zp 。现在,用 z1 代替 Q 的第1个零列,用 z2 代替第2个零列,如此做下去,直到所有零列都被这种方式所取代;用 Q′ 表示所得到的矩阵。因为 Q′ 的新列与 R 的零行相对应,于是 Q′ 有标准正交列,且 QR=Q′R ,因此 A=Q′R 是所要求的分解式。
如果 A 是实矩阵,注意到所有必要的运算都可经实的运算来完成,因而所得到的因子是实矩阵. □
练习 如果 A∈Mn,m ,且 n⩽m ,证明, A 可以分解成 A=LP ,其中, L∈Mn 是下三角矩阵,而 P∈Mn,m 具有标准正交行;至于其余的论断,这个分解有与(2.6.1)相平行的叙述。
练习 证明,任一个形如 B=A∗A 的矩阵 B∈Mn (其中 A∈Mn )都可以写成 B=LL∗ ,其中 L∈Mn 是具有非负对角元的下三角矩阵,并证明,如果 A 是非奇异矩阵,那么这个分解是唯一的。称这个分解为 B 的 Cholesky 分解;每个正定矩阵都可以按这种方式分解(见第7章)。提示:写出 A=QR 。
QR分解有重要的数值意义[例如,(2.6.3)],然而作为一种理论方法,它同样是有价值的。例如,已知矩阵 A∈Mn 经酉相似的上三角化,可以直接从它的通常的相似上三角化得到。
112
113
如果 S−1AS=T 是上三角矩阵,且如(2.6.1)中, S−QR ,那么 R−1Q∗AQR=T 且 Q∗AQ=RTR−1 ,而它作为上三角矩阵的乘积仍是上三角矩阵。用同样的方式可以推出,有关同时三角化的定理,例如(2.4.15),实际上是同时两三角化定理。这就是说,如果 Mn 中的某个矩阵族可经相似同时三角化,那么它同样可经西等价同时三角化。
下面,不加证明地叙述计算特征值的QR算法,并简要地说明它的某些特性。
2.6.3 QR 算法 设 Aij∈Mn 已知,写出 Aij−QijRi ,其中 Qij 和 Rij 如(2.6.1)所示,然后定义 Aij=RijQij ,再写出 Aij−QijRi ,其中 Qij 是酉矩阵,而 Rij 是上三角矩阵,这样继续做下去,一般地,如果因子 Aii=QiRi ,就定义 Ak−1=RkQk .
练习 证明,用QR算法得到的每个 Ak 西等价于 A1,k=1,2,⋯
在某些情形(例如,如果 Aii 的所有特征值有不同的绝对值),当 k→∞ 时,QR迭代值 Ak 将收敛于一个上三的矩阵,因为这个上三角矩阵西等价于 Aii ,所以 Aii 的各特征值便被计算出来。
如果 A0 是实矩阵, QR 迭代值 Ak 可以通过实的运算来计算。如果 A0 有非实特征值,那么就不能指望 QR 迭代值会收敛于一个上一角矩阵。这是因为这个上三角极限必须是实矩阵。不过,在某些情形,可以选取迭代值 Ak ,使它们收敛于一个具有 1×1 和 2×2 主对角子块的实分块上三角矩阵。而出现这种情形的充分条件是,除了各个成双成对的非实复共轭特征值有相同的模以外, A0 的所有实特征值都有不同的模。因为分块三角矩阵的各特征值是诸对角子块特征值的集合之并,所以 A0 的特征值是作为 QR 迭代值 Ak 的分块一角极限的 1×1 对角元以及该极限的 2×2 对角子块的特征值(的复共轭对)而出现的,它们可以很容易地用实数运算和二次公式计算出来。
2.6.4 例 下述例子说明, QR 算法不一定总收敛于一个三角矩阵。设
A=[0−110]. 于是, σ(A)={i} ,特征值没有不同的模。如果 Aij=A ,这个算法的一个可能的序列是
A1=[0110][−1001],A1−[−1001][0110]−[01−10]−[0110][1001],A0−[100…1][0110]=A0. 另一个是
A0=A0[1001]=A1. 在这两种情形都出现了循环,因而序列 {Ak} 不收敛于上三角矩阵。但是,存在 {Ak} 的一个选择,它收敛于一个分块上三角矩阵。
习题
设 x1,⋯,xm 是 Cn 中的给定无关向量组,且设 X≡[x1x2⋯xm]∈Mn,m ,假定对向量 x1,⋯,xm 实施(0.6.4)中所描述的Gram-Schmidt过程得到标准正交组 z1,⋯,zm 且设 Z≡∣z1⋯zm∣∈Mn,m . (a) 对 k=1,2,⋯,m ,设 Zk=[z1z2⋯zkxk−1xk+2⋯xm] ,其中 zk 是经Gram-Schmidt过程的第 k 步而得到的单位向量,因而 Zm=Z 。证明, Z1=XΔ1 , Z2=Z1Δ2 ,…, Zm=Zm−1Δm ,其中,每个 Δi 都是非奇异上三角矩阵,而 Δi=I−一个除第 i 列以外其余各列全为零的上三角矩阵. (b) 设 Ti=Δ1Δ2⋯Δk , k=1,2,⋯,m 。说明 Tk 是上三角矩阵且 Zk=XTk , k=1,2,⋯,m 。设 T=Tm ,因而 Z=XT 。(c)矩阵 T 与定理(2.6.1)证明中的上三角矩阵 R 之间的相互关系是什么?(d)证明,对 j−k+1,k+2,⋯,m , Zi 和 Ti 的前 k 列不会变动,因此,(gram-Schmidt过程的第 k 步就产生最终得到的矩阵 Z 和 T 的第 k 列.
设 r1,⋯,rm 是 Cn 中给定的线性无关向量组,又设 X≡∣x1⋯xn∣∈Mn,m 。考虑下面的算法:
令 Z=X ,且记 Z=[z1⋯zm] ,一开始先让每列 zi=xi,i=1,⋯,m .
Ⅱ. 对 k=1,2,…,m ,实施以下代换:
(1)首先,列 zk 用 (zk,zk)1,2 代替;
(ii)对 j=k+1,k+2,…,m ,每列 zi 用 zi−(zi,zk)zk 代替.
我们用 (x,y)=y′x 表示 Cn 上的普通内积。(a)证明这个过程的最后结果是矩阵 Z 具有标准正交列,它与习题1中用Gram-Schmidt过程得到的矩阵 Z 是相同的。(b)若 Zk 表示该算法在第 k 步中止时矩阵 Z 的存储信息, k−1,2,…,m ,证明 Z1−XΔ1,Z2−Z1Δ2,…,Zm−Zm+1Δm ;这里,每个 Δi 是一个非奇异上三角矩阵,其形式是 Δi−I−1 个除第 i 行以外其余各行全为零的上三角矩阵。(c)设 Tk=Δ1Δ2…Δk,k−1,2,…,m 。注意到 Tk 是上三角矩阵,而 Zk=XTk,k=1,2,…,m 。证明每个 Tk 的前 k 列与习题1中矩阵 Ti 的前 k 列是相同的,即使 Δk 和 Zk 可能各不相同。设 T−Tm 。(d)证明,对 j−k+1,k⋅2,…,m,Zj 和 Tj 的前 k 列不会变动,因此该算法的第 k 步就产生最终得到的矩阵 Z 和 T 的第 k 列。这个算法称为修改的Gram-Schmidt过程;它通过重新调整计算得到与通常的Gram-Schmidt过程相同的结果。虽然修改的算法与通常的Gram-Schmidt算法在数学上是等价的,但是对于数值计算来说,修改的Gram-Schmidt过程更受欢迎,因为它需要的存储较少,在遇到 X 的某些列接近平行这样一些难题时,它有助于得到一个能算出的 Z 。这个 Z 的诸列要比用通常的Gram-Schmidt算法得到的 Z 更接近于正交。在遇到难题时可以方便地通过列旋转策略使它的性能进一步改进:在实施步骤Ⅱ(I)之前,首先选作 zk 的是一个余下的列 zj,j⩾k ,其长度的平方 zj′zj 最大。在数值计算中每一步(不要求反演)实际上得到 Δi1 ,然后累积这些因子的乘积可以计算 X 的QR分解中的三角因子。
说明如何用一系列由 Householder 变换组成的乘积来实现 QR 分解。证明上述过程必须用 n 个 Householder 变换来完成,且 Q 是这些变换的乘积。一般认为,这个方法在计算上优于(2.6.1)证明中所采用的 Gram-Schmidt 算法。
如果应用于 A1∈Mn 的QR算法收敛于一个上三角矩阵,可以怎样计算 A0 的特征向量
呢?提示:需要解一个右边也是零的(奇异)三角形方程组
设将 QR 算法应用于某个矩阵 A∈Mn ,且设 {Ak} 是 QR 迭代值的序列。如果序列收敛,且 limk→∞Ak=B ,用选择原理(2.1.8)仔细说明为什么 B 酉等价于 A 。这为什么是重要的?
进一步阅读 关于Gram-Schmidt、修改的Gram-Schmidt以及几个其他的正交化过程的另外的参考资料和详细的叙述,可参看[GVI]的pp.146-169.关于QR算法和证明方法,以及另外的参考资料可见由D.Watkins写的综述,“Understanding the QR Algorithm,"SIAM Rev.24(1982);也可参看[Ste],427-440.