一：矩阵LU分解
矩阵的LU分解目的是将一个非奇异矩阵(A)分解成(A=LU)的形式，其中(L)是一个主对角线为(1)的下三角矩阵；(U)是一个上三角矩阵。
比如(A= egin{bmatrix} 1 & 2 & 4 \ 3 & 7 & 2 \ 2 & 3 & 3 \ end{bmatrix})，我们最终要分解成如下形式：
[A=Lcdot U =
egin{bmatrix}
1 & 0 & 0 \
3 & 1 & 0 \
2 & -1 & 1 \
end{bmatrix}
cdot
egin{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=egin{bmatrix}
1 & 2 & 4 \
3 & 7 & 2 \
2 & 3 & 3 \
end{bmatrix}
overset{(2)- color{red}{3} imes (1)}{underset{}{ o}}
egin{bmatrix}
1 & 2 & 4 \
0 & 1 & -10 \
2 & 3 & 3 \
end{bmatrix}
overset{(3)- color{red}{2} imes (1)}{underset{}{ o}}
egin{bmatrix}
1 & 2 & 4 \
0 & 1 & -10 \
0 & -1 & -5 \
end{bmatrix}
overset{(3)+ color{red}{1} imes (2)}{underset{}{ o}}
egin{bmatrix}
1 & 2 & 4 \
0 & 1 & -10 \
0 & 0 & -15 \
end{bmatrix}
=U
]
在运算过程中左下相应元素减去的倍数(上面红色的数字)便是矩阵(L)左下角的元素，可以得到：
[L=
egin{bmatrix}
1 & 0 & 0 \
color{red}{3} & 1 & 0 \
color{red}{2} & color{red}{-1} & 1 \
end{bmatrix}]
1.3 计算公式总结
通用计算公式是很重要的，因为有了公式之后，编程起来就方便很多了。我们可以根据上面的推导过程整理出如下伪代码：
[for ext{ } i = 1 : n hspace{6cm} \
for ext{ } j = i : n quad此时i为行下标，j为列下标\
qquad U_{ij}=A_{ij}-sum_{k=1}^{i-1} L_{ik}U_{kj} hspace{1cm}\
qquad for ext{ } 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;
}