矩阵的LU分解——MATLAB实现

---

前言

　　矩阵$A$的$LU$分解是将矩阵A表示成一个下三角矩阵$L$和上三角矩阵$U$的乘积，即$A=LU$，这样的结构便于科学计算。

---

一、矩阵的LU分解(三角分解)定义

　　设矩阵$A$是$m\times n$阶矩阵；$L$是$m\times m$阶下三角矩阵，主对角元素为1；$U$是$m\times n$阶上三角矩阵。称$LU$为矩阵$A$的一个三角分解，有：
$$A_{m\times n}=\left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}\\ a_{21}& a_{22}& a_{23}& \cdots & a_{2n}\\ a_{31}& a_{32}& a_{33}& \cdots & a_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& a_{m3}& \cdots & a_{mn}\end{array}\right]\text{\hspace{1em}(1-1)}$$
$$L_{m\times m}=\left[\begin{array}{ccccc}1& 0& 0& \cdots & 0\\ l_{21}& 1& 0& \cdots & 0\\ l_{31}& l_{32}& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ l_{m1}& l_{m2}& l_{m3}& \cdots & 1\end{array}\right]\text{\hspace{1em}(1-2)}$$
$$U_{m\times n}=\left[\begin{array}{ccccc}{u}_{11}& u_{12}& u_{13}& \cdots & u_{1n}\\ 0& u_{22}& u_{23}& \cdots & u_{2n}\\ 0& 0& u_{33}& \cdots & u_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & u_{mn}\end{array}\right]\text{\hspace{1em}(1-3)}$$
$$A_{m\times n}=L_{m\times m}{U}_{m\times n}\begin{array}{r}\Rightarrow \left[\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \cdots & a_{1n}\\ a_{21}& a_{22}& a_{23}& \cdots & a_{2n}\\ a_{31}& a_{32}& a_{33}& \cdots & a_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& a_{m3}& \cdots & a_{mn}\end{array}\right]=\left[\begin{array}{ccccc}{u}_{11}& u_{12}& u_{13}& \cdots & u_{1n}\\ l_{21}{u}_{11}& l_{21}{u}_{12}+u_{22}& l_{21}{u}_{13}+u_{23}& \cdots & l_{21}{u}_{1n}+u_{2n}\\ l_{31}{u}_{11}& l_{31}{u}_{12}+l_{32}{u}_{22}& l_{31}{u}_{13}+l_{32}{u}_{23}+u_{33}& \cdots & l_{31}{u}_{1n}+l_{32}{u}_{2n}+u_{3n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ l_{m1}{u}_{11}& l_{m1}{u}_{12}+l_{m2}{u}_{22}& l_{m1}{u}_{13}+l_{m2}{u}_{23}+l_{m3}{u}_{33}& \cdots & l_{m1}{u}_{1n}+l_{m2}{u}_{2n}+\cdots +l_{m,m-1}{u}_{m-1,n}+u_{mn}\end{array}\right]\end{array}$$
　　由上式，比较矩阵$A$和$LU$中的各项，可以得出通项：
$$a_{ij}=\{\begin{array}{ll}{u}_{ij}& j=1,i=1,2,\cdots ,m\\ {\sum }_{k=1}^{i-1}{l}_{ik}{u}_{kj}+u_{ij}& i\le j,i=2,3,\cdots m,j=i,i+1,\cdots ,n\\ {\sum }_{k=1}^{j-1}{l}_{ik}{u}_{kj}& i>j,j=2,3,\cdots ,n,i=j+1,\cdots ,m\end{array}\text{\hspace{1em}(1-5)}$$
　　可以求得上三角矩阵$U_{m\times n}$的表达式为：
$$u_{ij}=\{\begin{array}{ll}{a}_{ij}-{\sum }_{k=1}^{i-1}{l}_{ik}{u}_{kj}& i\le j\\ 0& i>j\end{array}\text{\hspace{1em}(1-6)}$$
　　同理，可以求得下三角矩阵$L_{m\times m}$的表达式为：
$$l_{ij}=\{\begin{array}{ll}0& i<j\\ 1& i=j\\ (a_{ij}-{\sum }_{k=1}^{j-1}{l}_{ik}{u}_{kj})/u_{jj}& i>j\end{array}\text{\hspace{1em}(1-7)}$$
　　从$L$的表达式(1-7)可以看出，当$u_{jj}=0$时，递推表达式产生了问题。返回表达式(1-4),这时候不难看出，当$u_{jj}=0$时，$l_{ij}$可以为任意值,等式恒成立。为了简化计算，当$u_{jj}=0$时，取$l_{ij}=0(i>j）$。
　　当矩阵$A$的各阶主子式$det(A_k)=det(A(1:k,1:k)),k=1,2,\cdots ,m$不为0时，$LU$分解式唯一；当矩阵$A$的$k$阶主子式为$0$时，$LU$分解式有无数个。通常在实际应用中，我们只需要一组$LU$分解式，为了简化计算，本文中，当$u_{jj}=0$时，取$l_{ij}=0(i>j)$。
　　当然了，当$u_{jj}=0$时，我们也可以通过对矩阵$A$进行初等行交换，使其各阶主子式不为$0$。上述过程相当于对矩阵左乘一个置换矩阵，来找到一组$LU$分解式，即$PA=LU$，$P$为一置换矩阵，此为MATLAB中的做法（函数lu）。
　　$LU$分解的更多信息参考：https://en.wikipedia.org/wiki/LU_decomposition#Closed_formula

---

二、LU分解实现

1.matlab代码

function [L,U] = lu_decompose(A)
%   lu decompose
%   L:下三角矩阵
%	U:上三角矩阵
%	A:输入矩阵
[m,n]=size(A);

L=eye(m);
L(:,1)=A(:,1)/A(1,1);%L第一列赋值

U=zeros(m,n);
U(1,:)=A(1,:);%U第一行赋值

for i=2:m
    for j=2:n
        if i<=j
            U(i,j)=A(i,j)-sum(L(i,1:i-1).*U(1:i-1,j)');%递推表达式(1-6)
        else
            if U(j,j)==0
                L(i,j)=0;
            else
                L(i,j)=(A(i,j)-sum(L(i,1:j-1).*U(1:j-1,j)'))/U(j,j);%递推表达式(1-7)
            end
        end
    end
end

end

2.数据验证

>> A=[1,2,3;4,5,6;7,8,9];
>> [L,U,P]=lu(A)
L =
       1              0              0       
       1/7            1              0       
       4/7            1/2            1       
U =
       7              8              9       
       0              6/7           12/7     
       0              0              0       
P =
       0              0              1       
       1              0              0       
       0              1              0  
>> [L1,U1]=lu_decompose(A)
L1 =
     1     0     0
     4     1     0
     7     2     1
U1 =
     1     2     3
     0    -3    -6
     0     0     0
>> [L1,U1]=lu_decompose(P*A)
L1 =
       1              0              0       
       1/7            1              0       
       4/7            1/2            1       
U1 =
       7              8              9       
       0              6/7           12/7     
       0              0              0   

---

总结

　　从矩阵的乘积的性质出发的得到上三角矩阵和下三角矩阵的表达式，并考虑上三角矩阵主对角元素为0的情况下上三角矩阵和下三角矩阵的表达式。

---