矩阵之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:[e-mail] allclcA=[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;   end   settle_num_row=row-i;   for j=1:settle_num_row      L(i+j,i)=A(i+j,i)/A(i,i);      A(i+j,:)=A(i+j,:)-L(i+j,i).*A(i,:);   endendL=L+eye(row);for i=1:row    for j=i:colum        U(i,j)=A(i,j);    endendfor i=1:row    P_matrix(i,P(i))=1;endPb=P_matrix*b;for i=1:row    if i==1       y(i)=Pb(i);    else        temp=0;        for k=1:i-1           temp=temp+L(i,k)*y(k);        end        y(i)=Pb(i)-temp;    endendy=y';for i=row:-1:1    if i==row        x(i)=y(i)/U(i,i);    else        temp=0;        for k=i:row            temp=temp+U(i,k)*x(k);        end         x(i)=(y(i)-temp)/U(i,i);    endendx=x';clear  j k index max_data settle_num_row tempfprintf('x=\n');for i=1:row    fprintf('%d\n',x(i));endclear i

最终结果如下：

x=
12
6
-13
-15