LU分解、矩阵求逆与解线性方程组(matlab代码)

一、前言

  参考知乎文章:LU分解、LUP分解、Cholesky分解中的LUP分解，LUP分解相比LU分解更加稳定，在矩阵求逆、求解解线性方程组有广泛运用。笔者修复了文章所述LUP算法一个BUG，当待分解的矩阵某一列均为0时，例如矩阵[0 1;0 1]，无法进行LUP分解。增加基于LUP分解求解线性方程组代码、在矩阵方程AX=B中，当B为单位阵时，利用LUP分解可以快速稳定地求解矩阵的逆。

二、LUP分解MATLAB代码

function [L, U, P] = lup(A)
if size(A, 1) ~= size(A, 2)
    error('矩阵行数与列数不等')
end

U = A;
L = eye(size(A));
P = eye(size(A));
n = size(A, 2);
for j = 1 : n - 1
    
    %Select i(>=j) that maximizes |U(i,j)|
    index = -1;
    maxValue = 0.0;
    for i = j : n
        temp = abs(U(i,j));
        if temp > maxValue
            maxValue = temp;
            index = i;
        end
    end
    
    if index == -1
        continue;
    end
    
    %Interchange rows of U: U(j, j : n) <-> U(i, j : n)
    for k = j : n
        temp = U(j, k);
        U(j, k) = U(index, k);
        U(index, k) = temp;
    end
    
    %Interchange rows of L: L(j, 1 : j - 1) <-> L(i, 1 : j - 1)
    for k = 1 : j - 1
        temp = L(j, k);
        L(j, k) = L(index, k);
        L(index, k) = temp;
    end
    
    %Interchange rows of P: P(j, 1 : n) <-> P(i, 1 : n)
    for k = 1 : n
        temp = P(j, k);
        P(j, k) = P(index, k);
        P(index, k) = temp;
    end
    
    for i = j + 1 : n
        L(i, j) = U(i, j) / U(j, j);
        for k = j : n
            U(i, k) = U(i, k) - L(i, j) * U(j, k);
        end
    end
end
end

三、LUP分解求解方程组MATLAB代码

function X = solve_matrix_equation_by_lup(A, B)
if size(A, 1) ~= size(A, 2)
    error('矩阵行数与列数不等')
end

U = A;
L = eye(size(A));
n = size(A, 2);
m = size(B, 2);
for j = 1 : n - 1
    
    %Select i(>=j) that maximizes |U(i,j)|
    index = -1;
    maxValue = 0.0;
    for i = j : n
        temp = abs(U(i,j));
        if temp > maxValue
            maxValue = temp;
            index = i;
        end
    end
    
    if index == -1
        continue;
    end
    
    %Interchange rows of U: U(j, j : n) <-> U(i, j : n)
    for k = j : n
        temp = U(j, k);
        U(j, k) = U(index, k);
        U(index, k) = temp;
    end
    
    %Interchange rows of L: L(j, 1 : j - 1) <-> L(i, 1 : j - 1)
    for k = 1 : j - 1
        temp = L(j, k);
        L(j, k) = L(index, k);
        L(index, k) = temp;
    end
    
    %Interchange rows of P: P(j, 1 : n) <-> P(i, 1 : n), C = P * B，等价于对B交换行
    for k = 1 : m
        temp = B(j, k);
        B(j, k) = B(index, k);
        B(index, k) = temp;
    end
    
    for i = j + 1 : n
        L(i, j) = U(i, j) / U(j, j);
        for k = j : n
            U(i, k) = U(i, k) - L(i, j) * U(j, k);
        end
    end
end

for i = 1 : n
    if abs(U(i, i)) < 1.0e-20
        disp('方程无解')
        return;
    end
end

%L * y = C
y = zeros(n, m);
for j = 1 : m
    for i = 1 : n
        sum = 0.0;
        for k = 1 : i - 1
            sum = sum + L(i, k) * y(k, j);
        end
        y(i, j) = B(i, j) - sum;
    end
end

%U * x = y
X = zeros(n, m);
for j = 1 : m
    for i = n :-1 : 1
        sum = 0.0;
        for k = i + 1 : n
            sum = sum + U(i, k) * X(k, j);
        end
        X(i, j) = (y(i, j) - sum) / U(i, i);
    end
end

end

四、LUP分解、矩阵求逆、求解线性方程组测试例子

clc;
clear;
A = [1 4 5;2 3 2;3 4 2];
B = [1 0 0;0 1 0;0 0 1];

[l,u,p] = lu(A)
[L,U,P] = lup(A)

invA = inv(A)
x = A\B
X = solve_matrix_equation_by_lup(A, B)