數(shù)值線性代數(shù) 課程設(shè)計報告

1、數(shù)值線性代數(shù)課程設(shè)計報告 課 程： 數(shù)值線性代數(shù) 題 目：正則化方法VS正交變換法 年 級： 2011級 專 業(yè)： 信息與計算科學 學 號： 11213013 姓 名： 劉 蘭 指導(dǎo)教師： 賈志剛 江蘇師范大學數(shù)學與統(tǒng)計學院摘要 正則化方法是求解最小二乘問題最古老的方法，然而在使用的過程中要注意最小二乘問題在化為正則化方程后的條件數(shù)是原問題的平方，使得求解過程增加了對舍入誤差的敏感性問題。而正交變換法可以說是求解最小二乘問題更實用的方法，本文將給出一個具體的算例，通過對用Cholesky分解、列主元Gauss消去法求解正則化方程組的方法和正交變換法求解結(jié)果的分析，直觀地給出這兩種方法的比較。關(guān)

2、鍵詞：最小二乘問題；正則化方程組；Cholesky分解法；列主元Gauss消去法；正交變換法.一、最小二乘問題定義3.1.1 給定矩陣及向量，確定使得        (3.1.3)該問題稱為最小二乘問題，簡記為LS(Least Squares)問題，其中稱為殘向量。在所討論的問題中，若線性地依賴于，則稱其為線性最小二乘問題；若非線性依賴于，則稱其為非線性最小二乘問題.二、正則化方法記最小二乘問題的解集為LS，即 LS=:x是LS問題（3.1.3）的解定理3.1.4 xLS當且僅當 (3.1.5)方程組(3.1.5)常常被稱為

3、最小二乘問題的正則化方程組，它是一個含有n個方程的線性方程組。在A的列向量線性無關(guān)的條件下，AA對稱正定，故可用平方根法求解線性方程組(3.1.5)，這樣我們就得到了求解最小二乘問題最古老的算法-正則化方法，其基本步驟如下：(1) 計算C=，d=；(2) 用平方根法計算C的Cholesky分解：C=LL；(3) 求解三角方程組Ly=d和Lx=y.三、正交變換法設(shè),.由于2范數(shù)具有正交不變性，故對任意的正交矩陣QR，有 這樣，最小二乘問題 就等價于原最小二乘問題(3.1.3).因此，我們就可望通過適當選取正交矩陣Q，使院問題轉(zhuǎn)化為較容易求解的最小二乘問題(3.1.3).這就是正交變換法的基本思想

4、.現(xiàn)在考慮如何求正交矩陣Q. 定理3.3.1(QR分解定理) 設(shè)（mn）,則A有QR分解： （3.3.2）其中是正交矩陣，是具有非負對角元的上三角陣；而且當m = n且A非奇異時，上述的分解還是唯一的。我們假定已知的QR分解（3.3.2）?，F(xiàn)將分塊為并且令那么由此即知，是最小二乘問題(3.1.3)的解當且僅當是的解。這樣一來，最小二乘問題(3.1.3)的解就可很容易從上三角方程組求得。綜上可得正交變換法的基本步驟為：（1） 計算A的QR分解（3.3.2）；（2） 計算；（3） 求解上三角方程組.四、具體算例令 其中k是一個大于1的參數(shù)，而 顯然，矩陣A的條件數(shù)為| A |2| A -1|2 =

5、 k,而最小二乘問題 r(x) = |Ax-b| 2= min的唯一解為 其中u1,u2和v1,v2分別是矩陣U和V的列向量，且直接計算有r (x*) = |Ax* - b|2 = 61/2注意：x* 和 r (x*)與k無關(guān). 對k = 105 , 107, 1010分別用如下的三種方法求解這一最小二乘問題： 用Cholesky分解法求解正則化方程組（簡記為NC）； 用列主元Gauss消去法求解正則化方程組（簡記為NG）； 正交變換法（簡記為QR）.五、結(jié)果分析計算結(jié)果列在了下表中，其中x表示計算解，x*表示精確解.從表中可以看出，隨著系數(shù)矩陣條件數(shù)的增大，正則化方法得到的計算解誤差越來越大

6、，而正交變換法得到的計算解就沒有受到任何影響，可見正則化方法和正交變換法之間有著明顯的差異。此外，就這一例子來看，用Cholesky分解法求解正則化方程組和用列主元Gauss消去法求解正則化方程組所得結(jié)果幾乎一樣。六、本次設(shè)計的總結(jié)短短的幾頁紙，卻花了很久的時間才做好。無論是程序的修改運行還是寫報告時各種公式符號的輸入，都琢磨了很長時間。也因為是第一次做矩陣計算的課程設(shè)計，很多東西不熟悉，而且自己目前所學到的知識也有限，所以是下了功夫才能出來的成果，雖然頁面不太美觀，但因為幾乎都是自己一個字一個符號打上去或從別的文檔上一個公式一個公式找了粘貼來的，還是很有成就感?？梢哉f通過這次課程設(shè)計真的收獲

7、了很多。這次的課程設(shè)計不僅檢驗并鞏固了我所學習的知識，也培養(yǎng)了我如何去把握一件事情，如何去完成一件事情的能力。其實每次寫電子文檔最讓人糾結(jié)的就是數(shù)學上的各種符號和公式的輸入了，更別說這次都是和矩陣有關(guān)的符號了，嘗試了各種方法，雖然最終的頁面效果還是差強人意。但通過這次課程設(shè)計，我在很多方面都有所提高。因為課程設(shè)計綜合運用了本專業(yè)所學課程的理論知識，可以幫助我一邊學習新的知識一邊回顧并加強已學的知識，讓我對矩陣計算的理論知識有了更直觀的理解，當然更重要的是培養(yǎng)和提高了我們對解決矩陣計算問題的興趣，我想學好矩陣計算的關(guān)鍵還在于當應(yīng)用它求解實際問題時，能最大程度的節(jié)省開銷以及減小誤差。所以課程設(shè)計對

8、我們來說是一個很好的鍛煉方式。參考文獻：1 張德豐 .MATLAB數(shù)值分析與應(yīng)用M.北京：國防工業(yè)出版社.2007.1 2 李慶揚，王能超，易大義.數(shù)值分析M.武漢：華中科技大學出版社.2006.7附錄：%前代法function y=qiandai(L,d,n)y(1)=d(1)/L(1,1);for k=2:n sum=L(k,1:k-1)*y(1:k-1)' y(k)=(d(k)-sum)/L(k,k);end%回帶法function x=huidai(U,y,n)x(n)=y(n)/U(n,n);for k=n:-1:1 sum=U(k,k+1:n)*x(k+1:n)'

9、x(k)=(y(k)-sum)/U(k,k); end% Cholesky分解function L,K=Cholesky(A,n)for k=1:n A(k,k)=sqrt(A(k,k); A(k+1:n,k)=A(k+1:n,k)/A(k,k); for j=k+1:n A(j:n,j)=A(j:n,j)-A(j:n,k)*A(j,k); endendL=tril(A);K=L'調(diào)用：U=1/sqrt(3),1/sqrt(2);-1/sqrt(3),0;1/sqrt(3),-1/sqrt(2);V=sqrt(3)/2,1/2;-1/2,sqrt(3)/2;S=105*sqrt(2),

10、0;0,sqrt(2); %k=105時A=U*S*V' %生成算例中的矩陣Ab=3;2;-1;C=A'*A; %正則化方法d=A'*b;L,K=cholesky(C,2)y=inv(L)*dx=inv(K)*yr=norm(A*x-b,2)x0=1;sqrt(3);r1=norm(x-x0,2)結(jié)果：A = 1.0e+004 * 7.0711 -4.0824 -7.0711 4.0825 7.0710 -4.0826C = 1.0e+010 * 1.5000 -0.8660 -0.8660 0.5000d = 2.0000 3.4641B = -1.1547r = 2

11、.4495r1 = 4.7677e-007%k=107時L = 1.0e+007 * 1.2247 0 -0.7071 0.0000K = 1.0e+007 * 1.2247 -0.7071 0 0.0000y = 0.0000 2.8174x = 0.9922 1.7186r = 2.4496r1 = 0.0155%k=1010時L = 1.0e+010 * 1.2247 0 -0.7071 0 + 0.0000iK = 1.0e+010 * 1.2247 -0.7071 0 0 - 0.0000iy = 0.0000 0 - 0.0361ix = 1.0e-003 * 0.1628 0.2

12、819r = 3.7420r1 = 2.0003%列主元Gauss消去法 function L,U1=lGauss(A,n)B=zeros(1,n);for k=1:n-1 p=k; for i=k+1:n if abs(A(i,k)>abs(A(k,k); p=i; end end B(1,:)=A(p,:); A(p,:)=A(k,:); A(k,:)=B(1,:); u(k)=p; if A(k,k)=0 A(k+1:n,k)=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); else f

13、printf('A is singular.'); break; endendL=tril(A,-1)+eye(n);U1=triu(A);調(diào)用：U=1/sqrt(3),1/sqrt(2);-1/sqrt(3),0;1/sqrt(3),-1/sqrt(2);V=sqrt(3)/2,1/2;-1/2,sqrt(3)/2;S=105*sqrt(2),0;0,sqrt(2);%k=105 %k=105時A=U*S*V' %生成算例中的矩陣Ab=3;2;-1;C=A'*Ad=A'*b;L,U1=lGauss(C,2)y=inv(L)*dx=inv(U1)*yr=

14、norm(A*x-b,2)x0=1;sqrt(3)r1=norm(x-x0,2)結(jié)果：L = 1.0000 0 -0.5774 1.0000U1 = 1.0e+010 * 1.5000 -0.8660 0 0.0000r = 2.4495r1 = 2.3849e-007%k=107時L = 1.0000 0 -0.5774 1.0000U1 = 1.0e+014 * 1.5000 -0.8660 0 0.0000r = 2.4495r1 = 0.0039%k=1010時L = 1.0000 0 -0.5774 1.0000U1 = 1.0e+020 * 1.5000 -0.8660 0 -0.

15、0000r = 3.7420r1 = 2.0003%Givens變換function G,y=givens(x,i,k)% 求解標準正交的Given變換矩陣G，使用Gx=y，其中y(k)=0，y(i)=sqrt(x(i)2+x(k)2)% x：需要進行Givens變換的列向量% i：變?yōu)閟qrt(x(i)2+x(j)2)的元素下標% k：變?yōu)?的元素的下標% G：Givens變換矩陣% y：Givens變換結(jié)果r=sqrt(x(i)2+x(k)2);c=x(i)/r;s=x(k)/r;G=eye(length(x);G(i,i)=c;G(i,k)=s;G(k,i)=-s;G(k,k)=c;y=

16、x;y(i)=r;y(k)=0;%基于Givens變換的QR分解function Q,R=QRgivens(A)n=size(A,1);R=A;Q=eye(n);for i=1:n-1 for k=2:n-i+1 x=R(i:n,i); rt=givens(x,1,k); %總是用向量的第一個元素去消下面其他的元素 r=blkdiag(eye(i-1),rt); %生成由一個以上的矩陣組成的對角陣 Q=Q*r' R=r*R; endend調(diào)用：U=1/sqrt(3),1/sqrt(2);-1/sqrt(3),0;1/sqrt(3),-1/sqrt(2);V=sqrt(3)/2,1/2;

17、-1/2,sqrt(3)/2;S=105*sqrt(2),0;0,sqrt(2); %k=105A=U*S*V' %生成算例中的矩陣Ab=3;2;-1;Q,R=QRgivens(A)Q1=Q(1,1),Q(1,2);Q(2,1),Q(2,2);Q(3,1),Q(3,2) %將Q分塊取n階的Q1c1=Q1'*bR0=R(1,1),R(1,2);R(2,1),R(2,2);x=huidai(R0,c1,2)' %最好不直接用逆求，而用回帶法避免R是奇異矩陣r=norm(A*x-b,2)x0=1;sqrt(3);r1=norm(x-x0,2)結(jié)果：Q = 0.5774 0.7

18、071 0.4082 -0.5774 0.0000 0.8165 0.5773 -0.7071 0.4082R = 1.0e+005 * 1.2247 -0.7071 -0.0000 0.0000 -0.0000 -0.0000Q1 = 0.5774 0.7071 -0.5774 0.0000 0.5773 -0.7071c1 = 0.0000 2.8284x = 1.0000 1.7321r = 2.4495r1 = 1.9202e-012%k=107時Q = 0.5774 0.7071 0.4082 -0.5774 0.0000 0.8165 0.5774 -0.7071 0.4082R = 1.0e+007 * 1.2247 -0.7071 0 0.0000 0 0Q1 = 0.5774 0.7071 -0.5774 0.0000