1、武漢大學數學與統計學院數值分析實驗報告實驗名稱關于正定矩陣cholesky分解的研究實驗時間2006年 12月 16 日姓名雷錦江 班級05級數類2班學號200531000174成績一、實驗目的，內容 二、相關背景知識介紹 三、代碼四、數值結果 五、計算結果的分析 六、計算中出現的問題，解決方法及體會一、 實驗目的，內容通過對正定矩陣cholesky分解問題的深入研究，掌握對矩陣進行cholesky分解的不同方法。同時對C語言的編程方法進行練習，在此基礎上，學會使用數學軟件matlab對相關問題進行處理。二、 相關背景知識介紹定理(Cholesky分解)如果A是對稱正定的，則存在惟一的一個對角

2、元全部大于零的下三角 陣，滿足。證明:由對稱矩陣的分解定理,存在單位下三角陣L和對角陣D=diag使得A=LD.由于大于零，則矩陣是對角元大于零的實下三角陣.它同時滿足.惟一性由分解的惟一性可推得。分解被稱為Cholesky分解,G被稱為Cholesky三角陣.如果我們計算Cholesky分解,然后解三角形方程組Gy=b和x=y,則b=Gy=G(x)=Gx=Ax.在上述定理中,Cholesky分解的證明是構造性的.然而,可以通過利用方程來得到計算Cholesky三角陣的更有效的方法.例：矩陣是正定的.常規Cholesky分解對于n階正定矩陣A，可導出按列的次序來計算其Cholesky分解因子G

3、的元素的計算公式。對第k=1,2,，n列，i=k+1,k+2,n 基于Gaxpy的Cholesky分解我們首先導出一個含有大量Gaxpy運算的Cholesky分解的實現方法.比較等式的第j列可得A(:,j)=.也就是說,G(j,j)(:,j)=A(:,j)-v.(*)如果G的第j-1列已知,則可計算出v.由(*)中各元素間的相等關系推出G(j:n,j)=v(j:n)/.這是數乘的Gaxpy運算。于是得到基于運算Gaxpy的Cholesky分解計算方法:For j =1:nv(j:n)=A(j:n,j)For k=1:j-1 v(j:n)=v(j:n)-G(j,k)G(j:n,k)EndG(j:

4、n,j)=v(j:n)/ End可以在計算過程中用G覆蓋A的下三角部分.算法(Cholesky分解:基于Gaxpy運算)給出對稱正定陣,本算法計算出一個下三角陣滿足.對所有ij,G(I,j)覆蓋A(i,j).For j=1:nIf j1A(j:n,j)=A(j:n,j)-A(j:n,1:j-1)EndA(j:n,j)=A(j:n,j)/End本算法需flops.基于外積的Cholesky分解另一種基于外積(秩為1)的Cholesky分解可通過對矩陣做如下劃分得到: (*)這里,且因是A正定的,故.而是的主子陣,其中故它也是正定的.如果存在Cholesky分解= ,則由(4.2.6)有,其中.因

5、此可通過反復利用(*)來得到最終的Cholesky分解,其方向和kji形式的高斯消去法一樣.算法(Cholesky分解:基于外積運算)給定一對稱正定陣.本算法計算滿足的下三角陣,對所有的ij,G(i,j)覆蓋A(i,j)For k=1:nA(k,k)=A(k+1:n,k)=A(k+1:n,k)/A(k,k)For j=k+1:nA(j:n,j)=A(j:n,j)-A(j:n,k)A(j,k)End End本算法需flops.注意,其中的j循環計算外積的下三角部分:A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k).回想在Gaxpy運算和外積運算的比較,容易得知Ga

6、xpy運算中讀寫向量的次數要比外積運算中少一半.基于分塊點積的Cholesky分解假定對稱正定,將A=()和其Cholesky因子看作含有方塊對角塊的N*N分塊陣.取出等式中關于第(I,j)塊(ij)的等式有=定義 S=-則有,當i=j時,當ij時, .通過合適的排序,這些方程可用來求得所有的算法(Cholesky:基于分塊點積運算)給定一個對稱正定陣,本算法可求得一個下三角陣滿足,A的下三角部分被G覆蓋,A被看作是含方對角塊的N*N分塊陣.For j=1:NFor i=j:NS=-If i=j計算Cholesky分解S=Else從=S解出End用覆蓋End End整個算法需flops,與前述

7、其他形式的Cholesky算法相同.假定A被適當分塊,則本算法中含有大量的矩陣乘法運算.例如,如果n=rN,且每個都是R*R的,則3級運算量約占1-(1/.).由于沒有給出基于gaxpy運算的Cholesky算法得出。將基于gaxpy的Cholesky分解執行r步后，我們已知中的和。再接下來我們不對A，而是對已明顯形成具有可利用的對稱結構的約化后的矩陣進行另外的r步基于gaxpy的Cholesky計算。按此方法進行下去，我們得到基于分塊的Cholesky計算。按此方法進行下去，我們得到基于分塊的Cholesky分解算法，去第k不是對n-(k-1)r階的矩陣進行r次基于gaxpy的Cholesk

8、y分解，接著是階為n-kr的三級算法。如果，則三級運算量約占1-3/（2N）。Cholesky分解的穩定性在精確運算下,我們知對稱正定陣存在Cholesky分解.反之,如果一個Cholesky分解過程能夠順利進行完且得到嚴格大于零的平方根,那么A是正定的.因此,判斷是否正定,我們只需用上述任一方法來試著計算其Cholesky分解.有舍入誤差是情形更為有意思. Cholesky算法的數值穩定性大致可從不等式:導出.該不等式說明Cholesky三角陣因子有很好的界.由也可推出相同的結論.Wilkinson(1968)在其經典的論文中對Cholesky分解的舍入誤差作了深入研究.利用該論文的結果可證

9、明,如果是通過上述任一Cholesky分解求得的Ax=b的計算解,則滿足擾動方程組(A+E)=b,其中=,是由決定的小常數. Wilkinson進一步指出,如果,其中是另一個小常數,則分解能夠執行到底,而不會出現對負數開平方根.三、代碼常規Cholesky分解#define N 3#includemath.hmain()float aNN,gNN,*pa,temp1,temp2; int i,j,k; pa=a0; printf(Please input the maxim:n); for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for

10、(i=1;iN;i+) for(j=0;ji;j+) aij=aji; for(k=0;kN;k+) temp1=0; for(j=0;jk;j+) temp1=temp1+gkj*gkj; gkk=sqrt(akk-temp1); for(i=k+1;iN;i+) temp2=0; for(j=0;jk;j+) temp2=temp2+gij*gkj; gik=(aik-temp2)/gkk; for(i=0;iN-1;i+) for(j=i+1;jN;j+) gij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f ,gij);

11、 getch(); 基于Gaxpy的Cholesky分解#define N 5#includemath.hmain()float aNN,*pa,temp; int i,j,k; pa=a0; printf(Please input the maxim:n); for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for(i=1;iN;i+) for(j=0;ji;j+) aij=aji; for(j=0;jN;j+) temp=0; for(k=0;kj;k+)temp=temp+ajk*ajk; ajj=sqrt(ajj-temp); f

12、or(i=j+1;iN;i+)temp=0; for(k=0;kj;k+) temp=temp+ajk*aik; aij=(aij-temp)/ajj; for(i=0;iN-1;i+) for(j=i+1;jN;j+) aij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f ,aij); getch(); 基于外積的Cholesky分解#define N 4#includemath.hmain()float aNN,*pa; int i,j,p,q; pa=a0; printf(Please input the maxim:n)

13、; for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for(i=1;iN;i+) for(j=0;ji;j+) aij=aji; for(j=0;jN;j+) for(p=j+1;pN;p+)for(q=j+1;qN;q+) apq=apq-apj*aqj/ajj; ajj=sqrt(ajj); for(i=j+1;iN;i+)aij=aij/ajj; for(i=0;iN-1;i+) for(j=i+1;jN;j+) aij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f

14、 ,aij); getch(); 基于分塊點積的Cholesky分解#define N 5#includemath.hmain()float aNN,*pa,temp; int i,j,k,p,q; pa=a0; printf(Please input the maxim:n); for(i=0;iN;i+) pa=pa+i; for(j=i;jN;j+)scanf(%f,pa+); for(i=1;iN;i+) for(j=0;ji;j+) aij=aji; p=N/3; for(q=0;qp;q+) for(j=3*q;j3*q+3;j+)temp=0; for(k=0;kj;k+) te

15、mp=temp+ajk*ajk; ajj=sqrt(ajj-temp); for(i=j+1;iN;i+) temp=0; for(k=0;kj;k+)temp=temp+ajk*aik; aij=(aij-temp)/ajj; if(N-p*3)=2) for(j=N-2;jN;j+) temp=0; for(k=0;kj;k+)temp=temp+ajk*ajk; ajj=sqrt(ajj-temp); for(i=j+1;iN;i+)temp=0; for(k=0;kj;k+) temp=temp+ajk*aik; aij=(aij-temp)/ajj; if(N-p*3)=1) tem

16、p=0; for(k=0;kN-1;k+) temp=temp+aN-1k*aN-1k; aN-1N-1=sqrt(aN-1N-1-temp); for(i=0;iN-1;i+) for(j=i+1;jN;j+) aij=0; for(i=0;iN;i+) printf(n); for(j=0;jN;j+)printf(%6.2f ,aij); getch();四數值結果常規Cholesky分解Please input the maxim:4 -2 4 10 -2 8 2.00 0.00 0.00 -1.00 3.00 0.00 2.00 0.00 2.00基于Gaxpy的Cholesky分解

17、Please input the maxim:1 1 1 1 1 2 2 2 2 3 3 3 4 4 1.00 0.00 0.00 0.00 0.00 1.00 1.00 0.00 0.00 0.00 1.00 1.00 1.00 0.00 0.001.00 1.00 1.00 1.00 0.001.00 1.00 1.00 1.00 1.00基于外積的Cholesky分解Please input the maxim:1 1 1 1 2 2 2 3 3 4 1.00 0.00 0.00 0.00 1.00 1.00 0.00 0.00 1.00 1.00 1.00 0.00 1.00 1.00

18、 1.00 1.00 基于分塊點積的Cholesky分解Please input the maxim9 6 3 12 3 29 12 13 12 41 30 47 37 40 64300 0.00 0.00 0.00 0.00 2.00 5.00 0.00 0.00 0.00 1.00 2.00 6.00 0.00 0.004.00 1.00 4.00 2.00 0.001.00 2.00 7.00 3.00 1.00五計算結果分析由于在編程中才有的是浮點數，所以由于系統的限制，該程序的精確度最高能達到6-7位有效數字，而且由于水平有限，該程序只是針對正定矩陣進行Cholesky分解的。對于標

19、準的正定矩陣，且在有效數字位數之內，該程序是有效的。在程序的輸出中，為了使輸出更為明了清晰，我選擇了保留兩位小數，再不影響精度的情況下，這種選取是可取的。六計算中出現的問題，解決方法及體會這篇論文從開始策劃到最終完成，共花了我們一個多星期的時間，從最開始的擔心無法完成課題的研究，到后來的查資料，編程，到最后的完成論文，我們體會到了其間的種種苦樂，當最終論文完成時，我們感到就如同完成了一件花了很多時間和精力完成的工藝品，也許有很多地方略顯粗糙，但這畢竟是我們努力的結果，我們對于這樣的結果是滿意的，已經沒有遺憾了。在進行Cholesky分解的研究過程中，我認為最關鍵的是對于算法的研究和程序的編寫，

20、由于相關資料略顯簡略，我花了一定的時間將算法弄透，接下來就進入了程序的編寫過程。在編程中，我確實遇到了一些困難，其中多次出現了參數設定的錯誤，而且該錯誤不能被低階的矩陣所反映，所以我曾一度以為程序編寫完成，而實際上程序中還存在著問題，當我對程序作最后測試的時候，才發現了問題，并花了相當長的時間完成了程序的矯正。從中我也得到了頗多的收獲。同時我想指出的是，在考慮基于分塊矩陣的Cholesky分解時，我曾一度想讓計算機能夠根據矩陣的階數自動將其分解，其具體方法是設定一個素數庫，然后在素數庫中按從小到大的順序，依次除以矩陣階數，直到出現能整除的素數為止，其數學基礎是應用了代數基本原則。但后來我又考慮

21、到如果矩陣的階數是一個很大的素數，那么這種方法將轉變為基于Gaxpy的Cholesky分解，而且由于素數庫的限制，該算法不能將某些矩陣分解。于是我最終采取了前面提到的關于基于分塊矩陣的Cholesky分解。最后，要感謝向老師對我們這個課題研究的關心，他在百忙之中抽空對我們的研究提出了寶貴的意見，在這里向他表示誠摯的感謝！七參考書目：1.MATRIX COMPUTATIONS（Third Edition） Gene H. Golub & Charles F. Van Loan2.數值計算方法 鄭慧嬈 陳紹林 莫忠息 黃象鼎（與邢周華同學合作完成）教 師 評 語指導教師： 年 月 日武漢大學數學與統計學院數值分析原始書記記錄實驗名稱：關于正定矩陣cholesky分解的研究 實驗時間