矩阵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;
}