一：矩阵LU分解
矩阵的LU分解目的是将一个非奇异矩阵\(A\)分解成\(A=LU\)的形式，其中\(L\)是一个主对角线为\(1\)的下三角矩阵；\(U\)是一个上三角矩阵。
比如\(A= \begin{bmatrix}
1 & 2 & 4 \\
3 & 7 & 2 \\
2 & 3 & 3 \\
\end{bmatrix}\)，我们最终要分解成如下形式：
\[A=L\cdot U =
\begin{bmatrix}
1 & 0 & 0 \\
3 & 1 & 0 \\
2 & -1 & 1 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
0 & 0 & -15 \\
\end{bmatrix}
\]
\begin{bmatrix}
1 & 0 & 0 \\
3 & 1 & 0 \\
2 & -1 & 1 \\
\end{bmatrix}
\cdot
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
0 & 0 & -15 \\
\end{bmatrix}
\]
现在主要的问题是如何由矩阵\(A\)计算得到矩阵\(L\)和\(U\)呢？我们将在下面详细讨论。
1.1 LU分解原理
首先从矩阵\(U\)入手，因为它是一个上三角矩阵，所以很容易想到高斯消元法，依次把矩阵\(A\)主对角线左下角的元素消为\(0\)就得到\(U\)了。
然后计算矩阵\(L\)，这里有个技巧，可以这样想，正是因为有了\(L\)，所以\(U\)的左下部分才能被消为\(0\)，所以我们记录一下把\(U\)的左下部分消为\(0\)时矩阵\(A\)每行所乘的倍数，这个减去的倍数便是\(L\)左下元素的值！
1.2 LU分解计算举例
\[A=\begin{bmatrix}
1 & 2 & 4 \\
3 & 7 & 2 \\
2 & 3 & 3 \\
\end{bmatrix}
\overset{(2)- \color{red}{3} \times (1)}{\underset{}{\to}}
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
2 & 3 & 3 \\
\end{bmatrix}
\overset{(3)- \color{red}{2} \times (1)}{\underset{}{\to}}
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
0 & -1 & -5 \\
\end{bmatrix}
\overset{(3)+ \color{red}{1} \times (2)}{\underset{}{\to}}
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
0 & 0 & -15 \\
\end{bmatrix}
=U
\]
1 & 2 & 4 \\
3 & 7 & 2 \\
2 & 3 & 3 \\
\end{bmatrix}
\overset{(2)- \color{red}{3} \times (1)}{\underset{}{\to}}
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
2 & 3 & 3 \\
\end{bmatrix}
\overset{(3)- \color{red}{2} \times (1)}{\underset{}{\to}}
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
0 & -1 & -5 \\
\end{bmatrix}
\overset{(3)+ \color{red}{1} \times (2)}{\underset{}{\to}}
\begin{bmatrix}
1 & 2 & 4 \\
0 & 1 & -10 \\
0 & 0 & -15 \\
\end{bmatrix}
=U
\]
在运算过程中左下相应元素减去的倍数(上面红色的数字)便是矩阵\(L\)左下角的元素，可以得到：
\[L=
\begin{bmatrix}
1 & 0 & 0 \\
\color{red}{3} & 1 & 0 \\
\color{red}{2} & \color{red}{-1} & 1 \\
\end{bmatrix}\]
\begin{bmatrix}
1 & 0 & 0 \\
\color{red}{3} & 1 & 0 \\
\color{red}{2} & \color{red}{-1} & 1 \\
\end{bmatrix}\]
1.3 计算公式总结
通用计算公式是很重要的，因为有了公式之后，编程起来就方便很多了。我们可以根据上面的推导过程整理出如下伪代码：
\[for \text{ } i = 1 : n \hspace{6cm} \\
for \text{ } j = i : n \quad此时i为行下标，j为列下标\\
\qquad U_{ij}=A_{ij}-\sum_{k=1}^{i-1} L_{ik}U_{kj} \hspace{1cm}\\
\qquad for \text{ } x = i+1 : n \quad 此时x为行下标，i为列下标\\
\qquad L_{xi}=(A_{xi}-\sum_{k=1}^{i-1} L_{xk}U_{ki}) /U_{ii} \hspace{0cm}\\
\]
for \text{ } j = i : n \quad此时i为行下标，j为列下标\\
\qquad U_{ij}=A_{ij}-\sum_{k=1}^{i-1} L_{ik}U_{kj} \hspace{1cm}\\
\qquad for \text{ } x = i+1 : n \quad 此时x为行下标，i为列下标\\
\qquad L_{xi}=(A_{xi}-\sum_{k=1}^{i-1} L_{xk}U_{ki}) /U_{ii} \hspace{0cm}\\
\]
其中\(n\)为方阵的行或列长度，可以看出先计算矩阵\(U\)的第一行，再计算矩阵\(L\)的第一列，再计算矩阵\(U\)的第二行，再计算矩阵\(L\)的第二列，依此类推。
二：矩阵LU分解MATLAB实现
clc,clear all,close all
% 矩阵的LU分解
%% 自己实现
A = [1 2 4;3 7 2;2 3 3]
[n,n] = size(A);
L = eye(n,n); % L初始化为单位矩阵
U = zeros(n,n); % U初始化为零矩阵
for i = 1 : n % 根据计算公式实现
for j = i : n
U(i,j) = A(i,j) - sum(L(i,1 : i - 1) .* U(1 : i - 1,j)');
end
for x = i + 1 : n
L(x,i) = (A(x,i) - sum(L(x,1 : i - 1) .* U(1 : i - 1,i)')) ./ U(i,i);
end
end
L
U
%% 内置函数实现
[L1,U1] = lu(A)
三：矩阵LU分解C++实现
#include <iostream>
#include <vector>
using namespace std;
int main()
{
vector<vector<double>> a = { {1,2,4},{3,7,2},{2,3,3} };
int n = a.size();
vector<vector<double>> u(n, vector<double>(n));
vector<vector<double>> l(n, vector<double>(n));
for (int i = 0; i < n; i++) //初始化矩阵L和矩阵U
for (int j = 0; j < n; j++)
{
u[i][j] = 0;
if (i == j) l[i][j] = 1;
}
for (int i = 0; i < n; i++)
{
double sum = 0;
for (int j = i; j < n; j++)
{
for (int k = 0; k <= i - 1; k++)
sum += l[i][k] * u[k][j];
u[i][j] = a[i][j] - sum; //计算矩阵U
sum = 0;
}
for (int x = i + 1; x < n; x++)
{
for (int k = 0; k <= i - 1; k++)
sum += l[x][k] * u[k][i];
l[x][i] = (a[x][i] - sum) / u[i][i]; //计算矩阵L
sum = 0;
}
}
cout << "A:" << endl; //输出矩阵A
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
printf("%.3f ", a[i][j]);
}
cout << endl;
}
cout << "L:" << endl; //输出矩阵L
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
printf("%.3f ", l[i][j]);
}
cout << endl;
}
cout << "U:" << endl; //输出矩阵U
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
printf("%.3f ", u[i][j]);
}
cout << endl;
}
return 0;
}