矩阵之LU分解
最近學習了矩陣的LU分解,LU分解是一種矩陣求費其次方程的算法,在matlab中通過矩陣求逆的方式(inv語句)就可以實現,求逆過程中,一般是通過高斯消去法進行求逆,算法復雜度n^3/2,而在其他語言中,一般求解廢棄次方程組是需要自己進行編寫的(或者從網上download一個程序~,因此編一個有效的解法可以大大降低程序的運行時間。
LU分解在本質上是高斯消元法的一種表達形式。實質上是將A通過初等行變換變成一個上三角矩陣,其變換矩陣就是一個單位下三角矩陣。這正是所謂的杜爾里特算法(Doolittle algorithm):從下至上地對矩陣A做初等行變換,將對角線左下方的元素變成零,然后再證明這些行變換的效果等同于左乘一系列單位下三角矩陣,這一系列單位下三角矩陣的乘積的逆就是L矩陣,它也是一個單位下三角矩陣。一般來說,如果LU分解只是為了單純求一個非齊次方程組,則沒有任何優勢可言,算法復雜度還要高于普通的高斯消去法,但是如果想要求解具有一些結果擾動的方程,即AX=b,b有很多情況,但這些情況只是細微的不同,此時,LU分解則在算法復雜度上具有一定的優勢。該算法主要參照了Carl D.Meyer的<Matrix Analysis and Applied Linear Algebra>一書。matlab編寫。
代碼如下:
%Programed by Lu Qi,University of Chinese Academy of Sciences %my email:qqlu1992@gmail.com clear all clc A=[1 2 -3 4;4 8 12 -8;2 3 2 1;-3 -1 1 -4]; b=[3 60 1 5]'; [row colum]=size(A); L=zeros(row); U=zeros(row); P_matrix=zeros(row); P=[1:colum]'; for i=1:row[max_data index]=max(abs(A(i:row,i))); %判斷某列最大數的位置max_data=A(index+i-1,i);if index~=i %進行主元的交換行temp=A(i,:);A(i,:)=A(index+i-1,:);A(index+i-1,:)=temp;temp=P(i);P(i)=P(index+i-1);P(index+i-1)=temp;temp=L(i,:);L(i,:)=L(index+i-1,:);L(index+i-1,:)=temp;endsettle_num_row=row-i;for j=1:settle_num_rowL(i+j,i)=A(i+j,i)/A(i,i);A(i+j,:)=A(i+j,:)-L(i+j,i).*A(i,:);end end L=L+eye(row); for i=1:rowfor j=i:columU(i,j)=A(i,j);end endfor i=1:rowP_matrix(i,P(i))=1; endPb=P_matrix*b; for i=1:rowif i==1y(i)=Pb(i);elsetemp=0;for k=1:i-1temp=temp+L(i,k)*y(k);endy(i)=Pb(i)-temp;end end y=y'; for i=row:-1:1if i==rowx(i)=y(i)/U(i,i);elsetemp=0;for k=i:rowtemp=temp+U(i,k)*x(k);endx(i)=(y(i)-temp)/U(i,i);end end x=x'; clear j k index max_data settle_num_row temp fprintf('x=\n'); for i=1:rowfprintf('%d\n',x(i)); end clear i最終結果如下:x=
12
6
-13
-15
總結
- 上一篇: 聚类之isodata算法
- 下一篇: 聚类之hierachical clust