One: matrix LU decomposition

The LU decomposition of matrix is to decompose A non-singular matrix AAA into the form A=LUA=LUA=LU, where LLL is A lower triangular matrix whose main diagonal is 111. UUU is an upper triangular matrix.

For example, A=[124372233]A= begin{bmatrix} 1&2&4 \\ 3&7&2 \\ 2&3&3 \\ end{bmatrix}A=⎣⎢, 132273423 \ ⎥⎤.


[ 1 0 0 3 1 0 2 – 1 1 ] ⋅ [ 1 2 4 0 1 – 10 0 0 – 15 ] \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}

Now the main question is how to calculate the matrices LLL and UUU from matrix AAA? We will discuss this in detail below.

1.1 LU decomposition principle

Firstly, we start with the matrix UUU, because it is an upper triangular matrix, so it is easy to think of the Gaussian elimination method, and then we get UUU by eliminating the elements in the lower left corner of the main diagonal of matrix AAA with 000.

Then calculate the matrix LLL, here is a trick, you can think of it as: because of LLL, so the lower left part of UUU can be eliminated to 000, so we record the multiples multiplied by each row of matrix AAA when the lower left part of UUU is eliminated to 000, the multiples subtracted is the value of the lower left element of LLL!

1.2 Example of LU decomposition calculation


A = [ 1 2 4 3 7 2 2 3 3 ] – ( 2 ) – 3 x ( 1 ) [ 1 2 4 0 1 – 10 2 3 3 ] – ( 3 ) – 2 x ( 1 ) [ 1 2 4 0 1 – 10 0 – 1 – 5 ] – ( 3 ) + 1 x ( 2 ) [ 1 2 4 0 1 – 10 0 0 – 15 ] = U 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

In the process of operation, the multiple of the corresponding element subtracted from the lower left (red number above) is the element in the lower left corner of matrix LLL, which can be obtained:


L = [ 1 0 0 3 1 0 2 – 1 1 ] L= \begin{bmatrix} 1 & 0 & 0 \\ \color{red}{3} & 1 & 0 \\ \color{red}{2} & \color{red}{-1} & 1 \\ \end{bmatrix}

1.3 Calculation formula summary

General formulas are important because they make programming a lot easier. We can sort out the following pseudocodes according to the above derivation process:


f o r   i = 1 : n f o r   j = i : n At this time i Is the row subscript, j For the column subscript U i j = A i j – ∑ k = 1 i – 1 L i k U k j f o r   x = i + 1 : n At this time x Is the row subscript, i For the column subscript L x i = ( A x i – ∑ k = 1 i – 1 L x k U k i ) / U i i for \text{ } i = 1 : n \hspace{6cm} \\ for \text{ } j = i : N \ quad I row subscript, at this time for column subscript j \ \ \ qquad U_ {ij} = A_ {ij} – \ sum_ {k = 1} ^ {1} I – L_ U_ {ik} {kj} \ img tags like hspace {1 cm} \ \ \ qquad for \ text {} x = I + 1: N \ quad for row subscript x at this time, I for column subscript \ \ \ qquad L_ {xi} = (A_ {xi} – \ sum_ {k = 1} ^ {1} I – L_ U_ {xk} {ki})/U_ {2} \ img tags like hspace \ \ {0 cm}

Where n is the row or column length of square matrix, it can be seen that the first row of matrix UUU is calculated, then the first column of matrix LLL is calculated, then the second row of matrix UUU is calculated, then the second column of matrix LLL is calculated, and so on.

Two: matrix LU decomposition MATLAB implementation

clc,clear all,close all
LU factorization of % matrix

%% implements itself
A = [1 2 4;3 7 2;2 3 3] 
[n,n] = size(A);
L = eye(n,n); % L is initialized to the identity matrix
U = zeros(n,n); % U is initialized to the zero matrix

for i = 1 : n % according to the calculation formula
    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
%% built-in function implementation

[L1,U1] = lu(A)
Copy the code

Three: matrix LU decomposition C++ implementation

#include <iostream>
#include <vector>
using namespace std;

int main(a)
{
    vector<vector<double>> a = { {1.2.4}, {3.7.2}, {2.3.3}};int n = a.size(a); vector<vector<double>> u(n, vector<double>(n));
    vector<vector<double>> l(n, vector<double>(n));

    for (int i = 0; i < n; i++) // Initialize the matrices L and 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; // Compute the matrix 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]; // Compute the matrix L
            sum = 0;
        }
    }

    cout << "A:" << endl; // Output matrix 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; // Output matrix 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; // Output matrix U
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
                printf("%.3f ", u[i][j]);
        }
        cout << endl;
    }
    return 0;
}
Copy the code