线性方程组的数值解法市公开课获奖课件

1、数值计算办法第二章 线性方程组数值解法第1页第1页第二章 线性方程组数值解法2.0 引 言2.1 Gauss消去法2.2 矩阵三角分解2.3 QR分解和奇异值分解第2页第2页 在自然科学和工程技术中诸多问题处理经常归结为解线性代数方程组。 比如：电学中网络问题 用最小二乘法求试验数据曲线拟合问题 解非线性方程组问题 用差分法或者有限元办法解常微分方程 偏微分方程边值问题等都造成求解线性代数方程组。 2.0 引 言第3页第3页这些方程组系数矩阵大体分为两种一个是低阶稠密矩阵(比如，阶数大约为150)另一个是大型稀疏矩阵(即矩阵阶数高且零元素较多) 2.0 引 言第4页第4页设有线性方程组Ax =

2、 b，其中 为非奇异阵， ， 关于线性方程组数值解法普通有两类：直接法与迭代法。2.0 引 言第5页第5页1. 直接法 就是通过有限步算术运算，可求得方程组准确解办法(若计算过程中没有舍入误差)。 但实际计算中由于舍入误差存在和影响，这种办法也只能求得线性方程组近似解。 本章将阐述这类算法中最基本高斯消去法及其一些变形。 这类办法是解低阶稠密矩阵方程组有效办法，近十几年来直接法在求解含有较大型稀疏矩阵方程组方面取得了较大进展。2.0 引 言第6页第6页2. 迭代法 就是用某种极限过程去逐步迫近线性方程组准确解办法。 迭代法含有需要计算机存贮单元较少、程序设计简朴、原始系数矩阵在计算过程中始终不

3、变等长处，但存在收敛性及收敛速度问题。 迭代法是解大型稀疏矩阵方程组(尤其是由微分方程离散后得到大型方程组)主要办法。 第6章简介迭代法解线性方程组。2.0 引 言第7页第7页 高斯(Gauss)消去法是解线性方程组最惯用办法之一 它基本思想是通过逐步消元（行初等变换），把方程组化为系数矩阵为三角形矩阵同解方程组，然后用回代法解此三角形方程组（简朴形式）得原方程组解。比如： 直接法基础2.1 Gauss消去法第8页第8页 下面讨论普通解n阶方程组高斯消去法。1. 消去过程 将原方程组记为 A(1)x =b(1)其中A(1)=(aij(1)nn=(aij)nn ，b(1)=b(1) 第一次消元。

4、2.1 Gauss消去法第9页第9页1. 消去过程(1) 第一次消元。其中2.1 Gauss消去法注：若a11(1) = 0，则在第1列中至少有一个元素不为0，可互换行后再消元第10页第10页(2) 第k次消元。2.1 Gauss消去法第11页第11页1. 消去过程(2) 第k次消元。注：为减少计算量，令 ，则2.1 Gauss消去法第12页第12页1. 消去过程(3) 当k = n 1时得完毕第n 1次消元后得到与原方程组等价三角形方程组A(n)x = b(n)注：当det(A) 0时，显然有aii(i) 0，(i = 1，n)称aii(i)为主元素。2.1 Gauss消去法第13页第13页

5、2. 回代过程求解三角形方程组A(n)x = b(n)，得到求解公式注：求解过程称为回代过程。2.1 Gauss消去法第14页第14页3. 算法设计令 i = k+1,k+2,n2.1 Gauss消去法输入A，n，epsFor (k = 1 to n - 1)for (i = k+1 to n)A(i，k) = A(i，k)/ A(k，k)for (j = k+1 to n+1)A(i,j) = A(i,j) - A(i,k)* A(k,j)A(n，n+1) = A(n，n+1) / A(n，n)for (k = n 1 to 1 step -1)w = 0for (j = k+1 to n)

6、w = w + A(k，j)* A(j，n+1)A(k，n+1) = A(k，n+1) wA(k，n+1) = A(k，n+1)/A(k，k)第15页第15页3. 算法设计function x = gauss(a)m=size(a);n=m(1);for k=1:n-1 for i=k+1:n a(i,:)=a(i,:)-a(k,:)*a(i,k)/a(k,k); endendfor j=n:-1:1 a(j,n+1)=(a(j,n+1)-a(j,j+1:n)*a(j+1:n,n+1)/a(j,j)endx=a(1:n,n+1);2.1 Gauss消去法输入A，n，epsFor (k = 1

7、to n - 1)for (i = k+1 to n)A(i，k) = A(i，k)/ A(k，k)for (j = k+1 to n+1)A(i,j) = A(i,j) - A(i,k)* A(k,j)A(n，n+1) = A(n，n+1) / A(n，n)for (k = n 1 to 1 step -1)w = 0for (j = k+1 to n)w = w + A(k，j)* A(j，n+1)A(k，n+1) = A(k，n+1) wA(k，n+1) = A(k，n+1)/A(k，k)第16页第16页3. 算法设计 i = k+1,k+2,n2.1 Gauss消去法第17页第17页3

8、. 算法设计 在Excel中1A1B1C1D1E12A2B2C2D2E23A3B3C3D3E34A4B4C4D4E35=A1=B1=C1=D1=E16=B2-B1*A2/A1=C2-C1*A2/A1=D2-D1*A2/A1=E2-E1*A2/A17=B3-B1*A3/A1=C3-C1*A3/A1=D3-D1*A3/A1=E3-E1*A3/A18=B4-B1*A4/A1=C4-C1*A4/A1=D4-D1*A4/A1=E4-E1*A4/A19=A5=B5=C5=D5=E510=B6=C6=D6=E611=C7-C6*B7/B6=D7-D6*B7/B6=E7-E6*B7/B612=C8-C6*B8

9、/B6=D8-D6*B8/B6=E8-E6*B8/B613=A9=B9=C9=D9=E9=(E13-D13*F16-C13*F15-B13*F14)/A1314=B10=C10=D10=E10=(E14-D14*F16-C14*F15)/B1415=C11=D11=E11=(E15-D15*F16)/C1516=D12-D11*C12/C11=E12-E11*C12/C11=E16/D162.1 Gauss消去法第18页第18页1A1B1C1D1E12A2B2C2D2E23A3B3C3D3E34A4B4C4D4E35=A1=B1=C1=D1=E16=B2-B$1*$A2/$A$1=C2-C$1

10、*$A2/$A$1=D2-D$1*$A2/$A$1=E2-E$1*$A2/$A$17=B3-B$1*$A3/$A$1=C3-C$1*$A3/$A$1=D3-D$1*$A3/$A$1=E3-E$1*$A3/$A$18=B4-B$1*$A4/$A$2=C4-C$1*$A4/$A$1=D4-D$1*$A4/$A$1=E4-E$1*$A4/$A$19=A5=B5=C5=D5=E510=B6=C6=D6=E611=C7-C$6*$B7/$B$6=D7-D$6*$B7/$B$6=E7-E$6*$B7/$B$612=C8-C$6*$B8/$B$6=D8-D$6*$B8/$B$6=E8-E$6*$B8/$B$

11、613=A9=B9=C9=D9=E9=(E13-D13*F16-C13*F15-B13*F14)/A1314=B10=C10=D10=E10=(E14-D14*F16-C14*F15)/B1415=C11=D11=E11=(E15-D15*F16)/C1516=D12-D$11*$C12/$C$11=E12-E$11*$C12/$C$11=E16/D162.1 Gauss消去法3. 算法设计 在Excel中第19页第19页4. Gauss消去法计算量以乘除法次数为主(1) 消元过程：第k步时(n k) + (n k) (n k + 1) = (n k) (n k + 2)共有2.1 Gauss

12、消去法第20页第20页4. Gauss消去法计算量(1) 消元过程：(2) 回代过程：求xi中, 乘ni次,除1次,共ni+1次(i = 1,n1)共有2.1 Gauss消去法第21页第21页4. Gauss消去法计算量(1) 消元过程：(2) 回代过程：(3) 总次数为注：当n = 20时约为2670次，比克莱姆法则9.71020次大大减少。2.1 Gauss消去法第22页第22页5. 阐明(1) 若消元过程中出现akk(k) = 0，则无法继续(2) 若akk(k) 0，但较小，则小主元做除数将产生大误差(3) 每次消元都选择绝对值最大者作主元，称为高斯主元消去法2.1 Gauss消去法第

13、23页第23页5. 阐明(4) 通常第k步选择第k列主对角元下列元素绝对值最大者作主元（该行与第k行对调），称为列主元消去法。2.1 Gauss消去法第24页第24页5. 阐明【例1】用舍入三位有效数字求解线性方程组（准确解是x1 = 10.0，x2 = 1.00）解：1) 不选主元Gauss消去法计算结果： x1 = -10.0，x2 = 1.01，此解无效； 2) 按列选主元Gauss消去法计算结果：x1 = 10.0，x2 = 1.00.2.1 Gauss消去法第25页第25页5. 阐明【例2】求解线性方程组 和（准确解是x1 = 1.0，x2 = 1.0）2.1 Gauss消去法第26

14、页第26页5. 阐明(5) 行标度化 当线性方程组系数矩阵元素大小相差很大时，则也许引进因丢失有效数字而产生误差，并且舍入误差影响严重，即使用Gauss主元消去法得到解也不可靠 为避免这一问题，可将系数矩阵行元素按百分比增减以使元素改变减小2.1 Gauss消去法第27页第27页5. 阐明(5) 行标度化 如用每行元素最大模除该行各元素，使它们模都小于1，这叫行标度化，其目的是要找到真正主元 消元过程仍是对原方程组进行，只但是在Gauss列主元消去法算法中，按列选主元ck时，应修改为其中2.1 Gauss消去法第28页第28页6. 列主元消去法算法输入A，n，epsFor (k = 1 to

15、n - 1)选主元A(I0，k) 拟定行号p = A(I0，k)若| p | = epsText = 1，退出循环若I0 kT换行（I0 k）消元若| A(n，n) | | A(I0，k)|TI0 = ifor (i = k+1 to n)A(i，k) = A(i，k)/ A(k，k)for (j = k+1 to n+1)A(i,j) = A(i,j) - A(i,k)* A(k,j)A(n，n+1) = A(n，n+1) / A(n，n)for (k = n 1 to 1 step -1)w = 0for (j = k+1 to n)w = w + A(k，j)* A(j，n+1)A(k，

16、n+1) = A(k，n+1) wA(k，n+1) = A(k，n+1)/A(k，k)for (j = k to n)n=A(k，j), A(k，j)=A(I0，j), A(I0，j)=n2.1 Gauss消去法第29页第29页6. 列主元消去法算法function x=gaussa(a)m=size(a); n=m(1); x=zeros(n,1);for k=1:n-1 c,i=max(abs(a(k:n,k); q=i+k-1; if q=k d=a(q,:);a(q,:)=a(k,:);a(k,:)=d end for i=k+1:n a(i,:)=a(i,:)-a(k,:)*a(i,

17、k)/a(k,k) endendfor j=n:-1:1 x(j)=(a(j,n+1)-a(j,j+1:n)*x(j+1:n)/a(j,j)end2.1 Gauss消去法输入A，n，epsFor (k = 1 to n - 1)选主元A(I0，k) 拟定行号p = A(I0，k)若| p | = epsText = 1，退出循环若I0 kT换行（I0 k）消元若| A(n，n) | 行标k时,lkm=0 j=k, k+1, ., nk=1, 2, ., n2.2 矩阵三角分解第35页第35页2.2.1 LU分解和LDU分解2. LU分解设A各阶顺序主子式不为0，且A = LU，即(2)主对角线

18、下边：当行标m 列标k时, umk= 0 i=k+1,k+2,.,nk=1,2,.,n1欲求lik和ukj2.2 矩阵三角分解第36页第36页2.2.1 LU分解和LDU分解2. LU分解(1)主对角线(含)上边:当列标m行标k时,lkm=0 j=k, k+1, .,nk=1, 2,.,n(2)主对角线下边：当行标m 列标k时，umk=0 i=k+1, k+2,.,nk=1, 2, ., n1设k = 1, a1j= l11u1j u1j = a1j j = 1, 2, ., n ai1= li1u11 li1 = ai1/u11 i = 2, 3, ., n欲求lik和ukj2.2 矩阵三角

19、分解第37页第37页2.2.1 LU分解和LDU分解2. LU分解普通地，最后：u11u12.u1n第1步l21u22.u2n第2步l31l32.ln1ln2unn第n步2.2 矩阵三角分解第38页第38页2.2.1 LU分解和LDU分解3. 求解方程组下三角 Ly = b：上三角 Ux = y：2.2 矩阵三角分解第39页第39页2.2.1 LU分解和LDU分解3. 求解方程组【例3】用杜丽特尔法解方程组解：2.2 矩阵三角分解第40页第40页2.2.1 LU分解和LDU分解3. 求解方程组【例3】用杜丽特尔法解方程组解：2.2 矩阵三角分解第41页第41页2.2.1 LU分解和LDU分解3

20、. 求解方程组【例3】用杜丽特尔法解方程组解：2.2 矩阵三角分解第42页第42页2.2.1 LU分解和LDU分解3. 求解方程组【例3】用杜丽特尔法解方程组解：2.2 矩阵三角分解第43页第43页2.2.1 LU分解和LDU分解3. 求解方程组【例3】用杜丽特尔法解方程组解：2.2 矩阵三角分解第44页第44页2.2.1 LU分解和LDU分解3. 求解方程组【例3】用杜丽特尔法解方程组解：2.2 矩阵三角分解第45页第45页2.2.1 LU分解和LDU分解4. 紧凑格式2.2 矩阵三角分解第46页第46页2.2.1 LU分解和LDU分解5. 代码function x=lua(a,b)m=si

21、ze(a);n=m(1);x=zeros(n,1);y=zeros(n,1);for k=1:n a(k,k:n) = a(k,k:n)-a(k,1:k-1)*a(1:k-1,k:n) a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*a(1:k-1,k)/a(k,k)endU = triu(a,0)L =eye(n)+ tril(a,-1)for i=1:n y(i) = b(i)-L(i,1:i-1)*y(1:i-1)endfor i=n:-1:1 x(i)=(y(i)-U(i,i+1:n)*x(i+1:n)/U(i,j)end2.2 矩阵三角分解第47页第47页

22、2.2.1 LU分解和LDU分解5. 代码紧凑格式：function x=luj(a)m=size(a);n=m(1);x=zeros(n,1)for k=1:n a(k,k:n+1) = a(k,k:n+1)-a(k,1:k-1)*a(1:k-1,k:n+1) a(k+1:n,k)=(a(k+1:n,k)-a(k+1:n,1:k-1)*a(1:k-1,k)/a(k,k)endU = triu(a,0)for j=n:-1:1 x(j)=(a(j,n+1)-U(j,j+1:n)*x(j+1:n)/U(j,j)end2.2 矩阵三角分解第48页第48页5791026681091871089225

23、7659=A1=B1=C1=D1=E1=(E5-(B5*F6+C5*F7+D5*F8)/A5=A2/A$5=B2-$A6*B5=C2-$A6*C5=D2-$A6*D5=E2-$A6*E5=(E6-(C6*F7+D6*F8)/B6=A3/A$5=(B3-A7*B$5)/B$6=C3-($A7*C5+$B7*C6)=D3-($A7*D5+$B7*D6)=E3-($A7*E5+$B7*E6)=(E7-D7*F8)/C7=A4/A$5=(B4-A8*B$5)/B$6=(C4-A8*C$5-B8*C$6)/C$7=D4-($A8*D5+$B8*D6+$C8*D7)=E4-($A8*E5+$B8*E6+$

24、C8*E7)=E8/D82.2 矩阵三角分解2.2.1 LU分解和LDU分解6. Excel实现第49页第49页2.2 矩阵三角分解2.2.1 LU分解和LDU分解6. Excel实现第50页第50页2.2.1 LU分解和LDU分解7. LDU分解设A = LU，令D = diag(u11, u22, ., unn)，则U1 = D-1U仍是一个单位上三角阵故矩阵三角分解除了杜丽特尔分解，尚有下列：(1) 克劳特分解：A = LU，L为下三角阵，U为单位上三角阵(2) LDU分解：A = LDU，L为单位下三角阵，D为对角阵，U为单位上三角阵上述分解均含有唯一性！2.2 矩阵三角分解第51页第

25、51页2.2.2 乔列斯基分解1. 原理定理：若A对称，则A可唯一分解为A = LDLT其中L为单位下三角，D为元素全是正数对角阵证实：设A = LDU，L为单位下三角阵，D为对角阵，U为单位上三角阵。由对称性LDU = UTDLT，由唯一性知L = UT，因此A = LDLT2.2 矩阵三角分解第52页第52页2.2.2 乔列斯基分解1. 原理定理：若A对称，则A可唯一分解为A = LDLT其中L为单位下三角，D为元素全是正数对角阵证实：设A = LDLT，D = diag(d1, d2, ., dn)对于ei = (0, ., 0, 1, 0, ., 0)T, 存在xi 0, 使得LTxi

26、 = ei另外, xiTAxi = xiT(LDLT)xi = (LTxi)TD(LTxi) = eiTDei = di由于A正定，有di = xiTAxi 0 2.2 矩阵三角分解第53页第53页2.2.2 乔列斯基分解1. 原理定理：若A对称，则A可唯一分解为A = LDLT其中L为单位下三角，D为元素全是正数对角阵乔列斯基分解: 若A正定，则A可唯一分解为A = GGT其中G为下三角阵记则有A = LDLT = LD1/2D1/2LT = (LD1/2)(LD1/2)T = GGT2.2 矩阵三角分解第54页第54页2.2.2 乔列斯基分解2. 算法设A = GGT则有：从而2.2 矩阵

27、三角分解第55页第55页2.2.2 乔列斯基分解2. 算法设A = GGT则有：注：由上式看出 ，因此 即算法是稳定2.2 矩阵三角分解第56页第56页2.2.2 乔列斯基分解2. 算法Ax = b GGTx = b Gy = b，GTx = y求解公式：上述办法称为“平方根法”注：不能选主元作行互换破坏对称性2.2 矩阵三角分解第57页第57页2.2.2 乔列斯基分解【例4】用平方根法解方程组解:2.2 矩阵三角分解第58页第58页2.2.2 乔列斯基分解3. 代码function x=choleskey(a,b)m=size(a);n=m(1);x=zeros(n,1);y=zeros(n

28、,1);G=zeros(m)for j=1:n G(j,j)=sqrt(a(j,j)-G(j,1:j-1)*G(j,1:j-1) G(j+1:n,j)=(a(j+1:n,j)-G(j+1:n,1:j-1)*G(j,1:j-1)/G(j,j)endfor i=1:n y(i)=(b(i)-G(i,1:i-1)*y(1:i-1)/G(i,i)endfor i=n:-1:1 x(i)=(y(i)-G(i+1:n,i)*x(i+1:n)/G(i,i)end2.2 矩阵三角分解第59页第59页2.2.2 乔列斯基分解4. Excel中计算2.2 矩阵三角分解第60页第60页2.2.3 追赶法在实际问题中

29、，经常碰到下列形式方程组这种方程组系数矩阵称为三对角矩阵。下列针对该方程组特点提供一个简便有效算法 追赶法2.2 矩阵三角分解第61页第61页2.2.3 追赶法作克劳特分解，A=LU，易知L，U具下列形式：比较第i行元素得2.2 矩阵三角分解第62页第62页2.2.3 追赶法作克劳特分解，A=LU，易知L，U具下列形式：计算过程：l1u1l2u2.ln2.2 矩阵三角分解第63页第63页2.2.3 追赶法有了ALU分解后，线性方程组Ax = d等价于两个简朴方程组：Ly = d，Ux = y Ly = f 即2.2 矩阵三角分解第64页第64页2.2.3 追赶法有了ALU分解后，线性方程组Ax

30、 = d等价于两个简朴方程组：Ly = d，Ux = y Ux = y即2.2 矩阵三角分解第65页第65页2.2.3 追赶法分解公式是计算y公式是计算过程：l1u1y1l2u2y2.lnyn (追)2.2 矩阵三角分解第66页第66页2.2.3 追赶法分解公式是计算y公式是计算过程：2.2 矩阵三角分解第67页第67页2.2.3 追赶法分解公式是计算x公式是(赶)2.2 矩阵三角分解第68页第68页2.2.3 追赶法2.2 矩阵三角分解第69页第69页2.2.3 追赶法【例5】用追赶法解解：追2.2 矩阵三角分解第70页第70页2.2.3 追赶法【例5】用追赶法解解：赶2.2 矩阵三角分解第

31、71页第71页2.2.3 追赶法2.2 矩阵三角分解输入a, b, c, dl1 = b1, y1 = d1/l1i = 1 to n 1ui = ci/lili+1 = bi+1 ai+1uiyi+1 = (di+1 ai+1yi)/li+1xn = yni = n 1 to 1xi = yi uixi+1第72页第72页2.2.3 追赶法function x=chase(a,b,c,d)m=size(b);n=m(2);L=linspace(0,0,n);U=L;x=L;y=L;L(1)=b(1);y(1)=d(1)/L(1);for k=1:n-1 U(k)=c(k)/L(k) L(k+

32、1)=b(k+1)-a(k)*U(k) y(k+1)=(d(k+1)-a(k)*y(k)/L(k+1)endx(n)=y(n)for i=n-1:-1:1 x(i)=y(i)-U(i)*x(i+1)end2.2 矩阵三角分解输入a, b, c, dl1 = b1, y1 = d1/l1i = 1 to n 1ui = ci/lili+1 = bi+1 ai+1uiyi+1 = (di+1 ai+1yi)/li+1xn = yni = n 1 to 1xi = yi uixi+1第73页第73页2.2.3 追赶法2.2 矩阵三角分解输入a, b, c, dd1 = d1/b1i = 1 to n

33、 1ci = ci/bibi+1 = bi+1 ai+1cidi+1 = (di+1 ai+1di)/bi+1i = n 1 to 1di = di cidi+1输入a, b, c, dl1 = b1, y1 = d1/l1i = 1 to n 1ui = ci/lili+1 = bi+1 ai+1uiyi+1 = (di+1 ai+1yi)/li+1xn = yni = n 1 to 1xi = yi uixi+1第74页第74页2.2.3 追赶法function d=chase(a,b,c,d)m=size(b);n=m(2);d(1)=d(1)/b(1);for k=1:n-1 c(k)

34、=c(k)/b(k); b(k+1)=b(k+1)-a(k)*c(k); d(k+1)=(d(k+1)-a(k)*d(k)/b(k+1);endfor i=n-1:-1:1 d(i)=d(i)-c(i)*d(i+1);end2.2 矩阵三角分解输入a, b, c, dd1 = d1/b1i = 1 to n 1ci = ci/bibi+1 = bi+1 ai+1cidi+1 = (di+1 ai+1di)/bi+1i = n 1 to 1di = di cidi+1第75页第75页2.2.3 追赶法2.2 矩阵三角分解第76页第76页定理：若A为对角占优(diagonally dominant

35、)三对角阵，且满足|b1| |c1| 0，|bi| |ai|+|ci| ，aici 0，i =2，3，n-1|bn| |an| 0.则追赶法可解以A为系数矩阵方程组注：1. 假如A是严格对角占优阵，则不要求三对角线上所有元素非零。2.2 矩阵三角分解第77页第77页注：1. 假如A是严格对角占优阵，则不要求三对角线上所有元素非零。2. 依据不等式可知：分解过程中，矩阵元素不会过度增大，算法确保稳定。3. 运算量为 O(6n)。2.2 矩阵三角分解第78页第78页2.3 QR分解和奇异值分解 矩阵分解是将矩阵分解为数个含有特殊性质矩阵因子乘积. 除了三角分解以外，尚有本节要简介QR分解（QR F

36、actorization）和奇异值分解法（Singular Value Decompostion）.第79页第79页2.3 QR分解和奇异值分解2.3.1 正交矩阵【定义2.3.1】若矩阵QRnn, 且满足QQT = QTQ = I则称矩阵Q为正交矩阵正交矩阵Q有下列性质： Q-1 = QT; det(Q) = 1; Qx长度与x长度相等.下面简介几类特殊正交矩阵.第80页第80页2.3 QR分解和奇异值分解2.3.1 正交矩阵1. 置换矩阵将单矩阵任意两行(列)互换得到矩阵，称为置换矩阵. 譬如, 将单位矩阵第i行和第j列互换, 得到置换矩阵Pij：任意个置换矩阵乘积仍然是置换矩阵.第81页

37、第81页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）对于某个角度，记s = sin, c = cos那么， 是一个正交矩阵. 记w = (x, y)T为二维平面中一个向量，用极坐标表示为w = (rsin, rcos)T. 那么即Gw表示将向量w顺时针旋转角所得到向量，如图2-1所表示.第82页第82页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）推广到n n情形，形如矩阵称为Givens矩阵或Givens变换，或称(平面)旋转矩阵（旋转变换），其中为旋转角度. 显然，G(i, j, )也是正交矩阵.第83页第83页2.

38、3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）若xRnn，y = G(i, j, )x，则y分量为假如要使yj = 0只要选择 满足即可第84页第84页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵。解：首先消去A(2, 1)，结构Givens变换G(1, 2, ), 其中第85页第85页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵

39、。解：从而G(1, 2, ) = 第86页第86页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵。解：A1 = G(1, 2, )A = 第87页第87页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵。解：另一方面消去A1(3, 2), 结构Givens变换G(2, 3, )，其中第88页第88页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Giv

40、ens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵。解：从而得A2 = G(2, 3, )A1 = 第89页第89页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵。解：最后消去A2(4, 3)。结构Givens变换G(3, 4, )，其中第90页第90页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵。解：从而，上三角矩阵R为

41、 A3 = G(3, 4, )A2 = 第91页第91页2.3 QR分解和奇异值分解2.3.1 正交矩阵2. 旋转矩阵（Givens变换）【例6】用Givens变换将海森伯格(Hdssenberg)型矩阵化为上三角矩阵。第92页第92页2.3 QR分解和奇异值分解2.3.1 正交矩阵3. 反射矩阵（Householder变换）设wRn, 且|w|2 = 1, 则P =I 2wwT称为Householder 变换或者Householder矩阵Householder矩阵有下列性质：(1) PT = P，即P是对称矩阵；(2) PPT =P2 =I 2wwT 2wwT + 4w(wTw)wT = I

42、, 即P是正交阵(3) 如图2-2，设w是R3上一个单位向量，并设S为过原点且与w垂直平面，则一切vRn能够分解成v = v1 + v2, 其中v1 S, v2S.第93页第93页2.3 QR分解和奇异值分解2.3.1 正交矩阵3. 反射矩阵（Householder变换）Householder矩阵有下列性质：(1) PT = P，即P是对称矩阵；(2) PPT =P2 =I 2wwT 2wwT + 4w(wTw)wT = I, 即P是正交阵(3) 如图2-2，设w是R3上一个单位向量，并设S为过原点且与w垂直平面，则一切vRn能够分解成v = v1 + v2, 其中v1 S, v2S.不难验证

43、Pv1 = v2, Pv2 = v2, 因此Pv = v1 v2这样，v经变换后象Pv是v关于S对称向量。第94页第94页2.3 QR分解和奇异值分解2.3.1 正交矩阵3. 反射矩阵（Householder变换）Householder矩阵有下列性质：(1) PT = P，即P是对称矩阵；(2) PPT =P2 =I 2wwT 2wwT + 4w(wTw)wT = I, 即P是正交阵(3) 如图2-2，设w是R3上一个单位向量，并设S为过原点且与w垂直平面，则一切vRn能够分解成v = v1 + v2, 其中v1 S, v2S.因此，Householder 变换又称镜面反射变换Househol

44、der矩阵也称初等反射矩阵。第95页第95页2.3 QR分解和奇异值分解2.3.1 正交矩阵3. 反射矩阵（Householder变换）一个主要应用是对x0, 求Householder矩阵P, 使得Px = ke1.由正交矩阵性质可知|Px|2 = |ke1|2 = |x|2, 即k =|x|2. 由上面所讨论P结构，有(令)u = x ke1, w = u / |u|2设x = (x1, , xn)T, 为了使x ke1计算时不损失有效数位取k = sign(x 1)|x|2, 则u = ( x 1+sign(x 1)|x|2, x 2, , xn)T从而P = I uuT。其中 = (|u

45、|22)1 = (|x|2(| x|2 + | x1|)-1e1= (1, 0, , 0)T第96页第96页2.3 QR分解和奇异值分解2.3.1 正交矩阵3. 反射矩阵（Householder变换）【例7】已知x = (3, 5, 1, 1)T, 求Householder矩阵P, 使得Px = 6e1，其中|x|2 = 6.解: 取k = 6, u = x ke1 = (9, 5, 1, 1)T, |u|2 = 108, = 1/54则第97页第97页2.3 QR分解和奇异值分解2.3.2 QR分解本节给出正交三角分解（又称QR分解）存在性定理和唯一性定理定理2.3.4 设ARnn, 则存在

46、正交阵P，使得PA = R, 其中R为上三角阵.第98页第98页2.3 QR分解和奇异值分解2.3.2 QR分解定理2.3.4 设ARnn, 则存在正交阵P，使得PA = R, 其中R为上三角阵.证实：首先考虑A第一列a1=(a11, a21, , an1)T, 可找到Householder矩阵P1, 使得P1a1元素除了第1个以外都为0.同理，找到P2使得P2P1A第二列对角元下列元素为0，而第一列对角元下列元素与P1A同样是0. 依次这样下去，能够得到Pn-1Pn-2P1A = R,其中R为上三角形矩阵，P = Pn-1Pn-2P1为正交阵. 给出结构性证实第99页第99页2.3 QR分解

47、和奇异值分解2.3.2 QR分解定理2.3.4 设ARnn, 则存在正交阵P，使得PA = R, 其中R为上三角阵.定理2.3.5 设ARnn，且A非奇异，则存在正交阵Q与上三角阵R，使得A有下列分解A = QR且当R对角元均为正时，分解是唯一.该定理确保了A可分解为A = QR，若A非奇异，则R也非奇异. 假如不要求R对角元为正，则分解不是唯一.第100页第100页2.3 QR分解和奇异值分解2.3.2 QR分解【例8】用Householder变换作矩阵AQR分解解：找Householder矩阵P1R33, 使则有第101页第101页2.3 QR分解和奇异值分解2.3.2 QR分解【例8】用

48、Householder变换作矩阵AQR分解解：再找 R22, 使 (1.44949, 3.44949)T = (*，0)T，得且第102页第102页2.3 QR分解和奇异值分解2.3.2 QR分解【例8】用Householder变换作矩阵AQR分解解：这是一个下三角矩阵, 但对角元皆为负数.只要令D = - I, R = - P2P1A就是对角元为正上三角矩阵, 使得A = QR, 其中第103页第103页2.3 QR分解和奇异值分解2.3.2 QR分解【例8】用Householder变换作矩阵AQR分解第104页第104页2.3 QR分解和奇异值分解2.3.2 QR分解 QR分解是计算特性值有力工具，也是用于其它矩阵计算问题，包括解方程组Ax = b. 这只要令y = QTb, 再解上三角形组Rx = y. 这个计算过程是稳定，也不必选主元，但是计算量比高斯消去法快要大一倍.第105页第105页2.3 QR分解和奇异值分解2.3.2 QR分解Matlab以qr函数来执行QR分解法，其语法为Q, R = qr(A)其中Q正交矩阵