上页下页5.7 矩阵的正交三角化及应用本节介绍初等反射阵及平面旋转阵，矩阵正交约 化，它们在矩阵计算中起着重要作用.上页下页5.7.1 初等反射阵定义9 设向量wRn且wTw=1，称矩阵H(w)=I-2wwT为初等反射阵(或称为豪斯霍尔德(Householder)变换). 如果记w =(w1,w2,wn)，则上页下页定理25 设有初等反射阵H(w)=I-2wwT, 其中 wTw=1, 则(1) H是对称矩阵, 即HT=H.(2) H是正交矩阵, 即H-1=H.(3) 设A为对称矩阵, 那么A1=H-1AH=HAH亦是对称矩阵.证明 只证H的正交性, 其它显然.设向量u0, 则显然 是一个初等反射阵.上页下页下面考察初等反射阵的几何意义. 考虑以w为法 向量且过原点O的超平面S: wTx=0. 设任意向量vRn, 则v=x+y, 其中xS, yS. 于是Hx=(I-2wwT)x=x-2wwTx=x.对于yS, 易知Hy=-y,从而对任意向量vRn, 总有Hv=x-y=v,其中v为v关于平面的镜面反射(见图5-1).wvySxv图 5-1上页下页定理26 设x, y为两个不相等的n维向量, |x|2=|y|2, 则存在一个初等反射阵H, 使Hx=y. 证明 令 , 则得到一个初等反射阵而且上页下页由|x|2=|y|2, 有yTy=xTx, 而数xTy=yTx, 从而所以得 Hx=x-(x-y)=y .容易说明, w是使Hx=y成立的唯一长度等于1的向量(不计符号).上页下页定理27(约化定理) 设x=(x1,x2,xn)T0, 则存在 初等反射阵H, 使Hx=-e1, 其中证明 见书p217.算法6 见书p218.上页下页5.7.2 平面旋转阵设x, yR2, 则变换是平面上向量的一个旋转变换，其中为正交矩阵.上页下页Rn中变换：y=Px，称为Rn中平面xi, xj的旋转变换(或称为吉文斯 (Givens)变换)，P=P(i,j,)=P(i,j)称为平面旋转矩阵.其中x=(x1,x2,xn)T, y=(y1,y2,yn)T, 使而上页下页显然，P(i, j,)具有性质：(1) P与单位阵I只是在(i, i), (i, j), (j, i) , (j, j)位置元 素不一样，其它相同.(2) P为正交矩阵(P-1=PT).(3) P(i, j)A(左乘)只需计算第i行与第j行元素，即 对A=(aij)mn有其中，c=cos，s=sin.(4) AP(i, j)(右乘)只需计算第i列与第j列元素，即上页下页利用平面旋转变换，可使向量x中的指定元素变为零.定理28(约化定理) 设x=(x1,xi , xj , xn)T, 其 中xi, xj不全为零，则可选择平面旋转阵P(i, j,) ，使其中证明 取 由利用矩阵乘 法，显然有于是，由c, s的取法得上页下页5.7.3 矩阵的QR分解下面讨论用正交矩阵来约化矩阵，可得到下述结果.设有设ARmn且为非零矩阵，则存在初等反射矩阵 H1, H2, , Hs使上页下页(1) 第1步约化：如果a1=0，取H1=I，即这一步不 需要约化，不妨设a10，于是可选取初等反射阵使于是其中上页下页(2) 第k步约化：设已完成对A上述第1步第k-1 步的约化，再进行第k步约化. 即存在初等反射阵H1, H2, , Hk-1使其中上页下页这里, Rk为k-1阶上三角阵,不妨设ck0, 否则这一步不需要约化(如果A列满秩 , 则ck0). 于是, 可选取初等反射阵使令上页下页第k步约化为令s=min(m-1, n), 继续上述过程, 最后有上页下页总结上述讨论给出下述结果.定理29(矩阵的正交约化定理) 设ARmn且A0, s=min(m-1, n), 则存在初等反射阵H1, H2, , Hs使且计算量约为n2m-n3/3(当mn)次乘法运算.上页下页(2) 设ARnn为非奇异矩阵, 则A有分解定理30(矩阵的QR分解) 其中R为n阶非奇异上三角阵.(1) 设ARmn且A的秩为n(mn), 则存在初等反射 阵H1, H2, , Hn使A=QR, 其中Q为正交矩阵, R为上三角阵. 且当R具有正对角元 素时, 分解唯一.上页下页证明 (1)由定理29可得.(2) 由设及定理29存在初等反射阵H1, H2, , Hn-1使记QT=Hn-1H2H1, 则上式为QTA=R, 即 A=QR, 其中Q为正交矩阵, R为上三角阵. 上页下页唯一性, 设有A=Q1R1 =Q2R2, 其中Q1, Q2为正交矩阵, R1, R2为非奇异上三角阵, 且R1, R2具有正对角元素,则由假设, 及对称正定矩阵ATA的Cholesky分解的唯一性, 则R1=R2. 从而可得Q1=Q2.上页下页下面考虑平面旋转变换来约化矩阵.定理31(用吉文斯变换计算矩阵的QR分解) 设 ARnn为非奇异矩阵, 则证明见书p224-自看.(1) 存在正交矩阵P1, P2, , Pn-1使(2) A有QR分解: A=QR. 其中Q为正交阵, R为非奇异上三角阵. 且当R对角元素 都为正时, 分解是唯一的.上页下页第8章 矩阵特征问题的计算 8.1 引言 8.2 幂法及反幂法 8.3 豪斯霍尔德方法 8.4 QR方法上页下页8.1 引 言工程技术中有多种振动问题，如桥梁或建筑物的振动，机械零件、飞机机翼的振动，及一些稳定性分析和相关分析在数学上都可转化为求矩阵特征值与特征向量的问题.下面先复习一些矩阵的特征值和特征向量的基础 知识.上页下页定义1 已知n阶矩阵A=(aij)，则称为A的特征多项式. 一般有n个根(实的或复的，复根按重数计算)称为A的特征值. 用(A)表示A的所有特征值的集合.A的特征方程上页下页 设为A的特征值，相应的齐次方程组注：当A为实矩阵时， ()=0为实系数n次代数方程，其复根是共轭成对出现.的非零解x称为矩阵A的对应于的特征向量.例1 求A的特征值及特征向量，其中上页下页解 矩阵A的特征方程为求得矩阵A的特征值为：对应于各特征值矩阵A的特征向量分别为：上页下页定理1 设为ARnn的特征值, 且Ax=x (x0)，则有 -p为A-pI的特征值，即(A-pI)x=(-p)x ； c为的cA特征值(c0为常数)；下面叙述有关特征值的一些结论： k为Ak的特征值，即Akx=kx ； 设A为非奇异矩阵，那么0 , 且-1为A-1的特 征值，即A-1x=-1x .上页下页定理2 设i(i=1,2,n)为n阶矩阵A=(aij)的特征值，则有 称为A的迹；定理3 设ARnn，则有定理4 设A 为分块上三角矩阵，即其中每个对角块Aii均为方阵，则上页下页定理5 设A与B为相似矩阵（即存在非奇异矩阵P使B=P-1AP），则定理5说明，一个矩阵A经过相似变换，其特征值不变.一个亏损矩阵是一个没有足够特征向量的矩阵，亏损矩阵在理论上和计算上都存在困难. A与B有相同的特征值； 如果y是B的特征向量，则Py是A的特征向量.定义2 如果实矩阵A有一个重数为k的特征值, 且对应于的A的线性无关的特征向量个数|2|n|，则对任何非零向量v0(a10)，幂法的算式成立.又设A有n个线性无关的特征向量，1对应的r个线性 无关的特征向量为x1,x2,xr，则由(2.2)式有如果A的主特征值为实的重根, 即1=2=r, 且|r|r+1|n|，上页下页为A的特征向量，这说明当A的主特征值是实的重根时，定理5的结论还是正确的.应用幂法计算A的主特征值1及其对应的特征向量时，如果|1|1(或|1|2|n|，则对任意非零初始 向量v0=u0(a10)，有幂法计算公式为则有上页下页例1 用幂法计算矩阵的主特征值与其对应的特征向量.解 取 v0=u0=(0,0,1)T , 则上页下页直到k=8 时的计算结果见下表k 1 2, 4, 1, 4 0.5, 1, 0.25 2 4.5, 9, 7.75 90.5, 1, 0.8611 3 5.7222, 11.4444, 8.36111.44440.5, 1, 0.7360 4 5.4621, 10.9223, 8.2306 10.92230.5, 1, 0.7536 5 5.5075, 11.0142, 8.2576 11.01420.5, 1, 0.7494 6 5.4987, 10.9974, 8.2494 10.99740.5, 1, 0.7501 7 5.5002, 11.0005, 8.2501 11.00050.5, 1, 0.7500 8 5.5000, 11.0000, 8.2500 11.00000.5, 1, 0.7500从而见书p303-例3.上页下页8.2.2 幂法的加速方法1、原点平移法由前面讨论知道，应用幂法计算A的主特征值的 收敛速度主要由比值 r=|2/1|来决定，但当r 接近于1 时，收敛可能很慢. 这时，一个补救办法是采用加速 收敛的方法.其中p为参数，设A的特征值为i，则对矩阵B的特征 值为i-p ，而且A, B的特征向量相同.引进矩阵 B=A-pI .上页下页如果要计算A的主特征值1, 只要选择合适的数p ，使1-p为矩阵B=A-pI 的主特征值，且 那么，对矩阵B=A-pI应用幂法求其主特征值1-p, 收 敛速度将会加快. 这种通过求B=A-pI的主特征值和特 征向量，而得到A的主特征值和特征向量的方法叫原 点平移法. 对于A的特征值的某种分布，它是十分有 效的.上页下页例4 设AR44有特征值比值r=|2/1|0.9. 做变换 B=A-12I (p=12), 则B的特征值为应用幂法计算B的主特征值1的收敛速度的比值为虽然常常能够选择有利的p值, 使幂法得到加速, 但设计一个自动选择适当参数p的过程是困难的.上页下页下面考虑当A的特征值是实数时，怎样选择p使采 用幂法计算1得到加速.且使收敛速度的比值设A的特征值都是实数，且满足则对实数p，使矩阵A-pI的主特征值为1-p或n-p时 ，当我们计算1及x1时，首先应选取p使上页下页显然，当2-p=-(n-p )时，即 P=(2+n)/2P* 时为最小值，这时收敛速度的比值为当希望计算n时，应选取p=(1+n-1)/2P* 使得应用幂法计算n得到加速.当A的特征值都是实数，满足且2, n能初步估计出来，我们就能确定P*的近似值.上页下页例2 用原点平移加速法求例1中矩阵A的主特征值 与其对应的特征向量.对B应用幂法，仍 取 v0=(0,0,1)T , 则解 取p=-2.5, 做平移变换B=A-pI，则上页下页迭代5步的计算结果见下表k1 2, 4, 3.54 0.5, 1, 0.875 2 7, 14, 10.562