最近由于数值分析实验课要求，需要通过matlab实现矩阵的LU分解。但是看了很多网友写的程序，基本上都是通过循环嵌套循环来实现矩阵的LU分解。略感琐碎，因此最近两天便一直在思考能否利用矩阵的乘v法，来简化循环，让代码更加美观，便于理解。话不多说，直接上代码。

代码($A$为待分解矩阵，要求$A$的顺序主子式不为0)1

2

3

4

5

6

7

8

9

10n = size(A,1);

L = eye(n);

U = zeros(n);

U(1,:) = A(1,:);

L(:,1) = A(:,1)/U(1,1);

for i=2:n % 每次循环计算U的一行 L的一列

U(i,i:n) = A(i,i:n) - L(i,:)*U(:,i:n);

L(i:n,i) = (A(i:n,i) - L(i:n,:)*U(:,i))/U(i,i);

L = L - diag(diag(L)-1);

end

分解思路

首先我们回顾一下矩阵的LU分解，矩阵的LU分解就是将一个矩阵$A$,分解为两个矩阵的乘积，记为$LU$,其中$L$为下三角形矩阵，$U$为上三角矩阵。以4阶矩阵为例，数学表达式如下：

$$left[

begin{matrix}

a_{11} & a_{12} & a_{13} & a_{14}\

a_{21} & a_{22} & a_{23} & a_{24} \

a_{31} & a_{32} & a_{33} & a_{34} \

a_{41} & a_{42} & a_{43} & a_{44}

end{matrix}

right] =

left[

begin{matrix}

1 & 0 & 0 & 0 \

l_{21} & 1 & 0 & 0 \

1_{31} & l_{32} & 1 & 0 \

1_{41} & l_{42} & l_{43} & 1 \

end{matrix}

right]

left[

begin{matrix}

u_{11} & u_{12} & u_{13} & u_{14} \

0 & u_{22} & u_{23} & u_{24} \

0 & 0 & u_{33} & u_{34} \

0 & 0 & 0 & u_{44} \

end{matrix}

right] $$

其中$A$矩阵已知，需要求矩阵$L$和矩阵$U$。

求解思路是通过矩阵的乘法：需要求$u_{11}$就去寻找对应位置$a_{11}$是如何被计算的，需要求$l_{21}$就去寻找$a_{21}$是如何被计算的。

例如：我们需要求$u_{13}$，由矩阵的乘法可以知道：$a_{13}$等于$L$的第一行乘以$U$的第三列，那么自然地：$a_{13} = u_{13}$。那么我们就求得：$u_{13} = a_{13}$

代码解释

1-3行代码：初始化(自己去Matlab中运行以下就知道是什么意思了)

4-5行代码：计算$U$的第一行与$L$的第一列

6-10行代码：通过循环计算$U$的第$k$行与$L$的第$k$列

重点解析后四行代码

众所周知(任何一本书上都会详细的讲到)，计算矩阵$U$的第$k$行与矩阵$L$的第$k$列的公式为：(推导思路和上述分解思路一致，可以自行验证)

$$u_{kj}=a_{kj}-sumlimits_{q=1}^{k-1}l_{kq}u_{qj}，j=k,…,n$$

$$l_{ik}=frac{a_{ik}-sumlimits_{q=1}^{k-1}l_{iq}u_{qk}}{u_{kk}}，i=k+1,k+2,…n$$

粗略分析，要计算$U$的第$k$行，那么就需要算$k-n+1$个元素，每个元素的计算中又有一个求和，再加上不仅要算$U$的第$k$行，还有$U$的第$k+1$到$n$行。因此没有分析错的话应该是三重循环。

下面来解释我通过利用矩阵的乘法来简化循环的思路：

以$n=4，i=2$为例,此时的矩阵形态为：

$$A=left[

begin{matrix}

a_{11} & a_{12} & a_{13} & a_{14}\

a_{21} & a_{22} & a_{23} & a_{24} \

a_{31} & a_{32} & a_{33} & a_{34} \

a_{41} & a_{42} & a_{43} & a_{44}

end{matrix}

right]

$$

$$

L =

left[

begin{matrix}

1 & 0 & 0 & 0 \

l_{21} & 1 & 0 & 0 \

1_{31} & 0 & 1 & 0 \

1_{41} & 0 & 0 & 1 \

end{matrix}

right],U=

left[

begin{matrix}

u_{11} & u_{12} & u_{13} & u_{14} \

0 & 0 & 0 & 0 \

0 & 0 & 0 & 0 \

0 & 0 & 0 & 0 \

end{matrix}

right] $$

需要计算矩阵$U$的第二行与矩阵$L$的第二列，可以观察到图中分块矩阵存在运算关系，因此可以利用此关系，直接运算得到$U$的第二行。

从而避免了对其单个元素计算，提升了代码的可读性。