矩阵LU分解的MATLAB与C++实现

一：矩阵LU分解

矩阵的LU分解目的是将一个非奇异矩阵\(A\)分解成\(A=LU\)的形式，其中\(L\)是一个主对角线为\(1\)的下三角矩阵；\(U\)是一个上三角矩阵。

比如\(A= \begin{bmatrix}
1 & 2 & 4 \\
3 & 7 & 2 \\
2 & 3 & 3 \\
\end{bmatrix}\)，我们最终要分解成如下形式：

现在主要的问题是如何由矩阵\(A\)计算得到矩阵\(L\)和\(U\)呢？我们将在下面详细讨论。

1.1 LU分解原理

首先从矩阵\(U\)入手，因为它是一个上三角矩阵，所以很容易想到高斯消元法，依次把矩阵\(A\)主对角线左下角的元素消为\(0\)就得到\(U\)了。

然后计算矩阵\(L\)，这里有个技巧，可以这样想，正是因为有了\(L\)，所以\(U\)的左下部分才能被消为\(0\)，所以我们记录一下把\(U\)的左下部分消为\(0\)时矩阵\(A\)每行所乘的倍数，这个减去的倍数便是\(L\)左下元素的值！

1.2 LU分解计算举例

在运算过程中左下相应元素减去的倍数(上面红色的数字)便是矩阵\(L\)左下角的元素，可以得到：

1.3 计算公式总结

通用计算公式是很重要的，因为有了公式之后，编程起来就方便很多了。我们可以根据上面的推导过程整理出如下伪代码：

其中\(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;
}