在 MATLAB 中执行 LU 分解而不旋转

Question

如何在MATLAB中实现函数lu(A)，让L*U直接A，同时得到真正的L矩阵？

当我使用 [L,U] = lu(A) 时，MATLAB 没有给我正确的 L 矩阵。当我使用[L,U,P] = lu(A)时，我需要实现P*A = L*U，但我只想乘L*U得到A.

Answer 1

MATLAB 的 lu 默认情况下始终执行旋转。例如，如果您在尝试执行传统的 LU 分解算法时有一个等于 0 的对角系数，则它将不起作用，因为在执行高斯消去法以创建上三角矩阵时需要对角系数 U所以你会得到除以零的错误。需要旋转以确保分解稳定。

但是，如果你能保证你的矩阵的对角线系数不为零，这很简单，但你必须自己写这个。您所要做的就是对矩阵进行高斯消去，将矩阵降为降阶梯形。结果降阶梯形矩阵为U，而高斯消元中去除L的下三角部分所需的系数将放在下三角一半，使U.

假设您的矩阵存储在 A 中，这样的方法可能有效。请记住，我在这里假设一个方阵。非旋转 LU 分解算法的实现放在名为 lu_nopivot:

的 MATLAB 函数文件中

function [L, U] = lu_nopivot(A)

n = size(A, 1); % Obtain number of rows (should equal number of columns)
L = eye(n); % Start L off as identity and populate the lower triangular half slowly
for k = 1 : n
    % For each row k, access columns from k+1 to the end and divide by
    % the diagonal coefficient at A(k ,k)
    L(k + 1 : n, k) = A(k + 1 : n, k) / A(k, k);

    % For each row k+1 to the end, perform Gaussian elimination
    % In the end, A will contain U
    for l = k + 1 : n
        A(l, :) = A(l, :) - L(l, k) * A(k, :);
    end
end
U = A;

end

作为一个运行示例，假设我们有以下 3 x 3 矩阵：

>> rng(123)
>> A = randi(10, 3, 3)

A =

     7     6    10
     3     8     7
     3     5     5

运行算法给我们：

>> [L,U] = lu_nopivot(A)

L =

    1.0000         0         0
    0.4286    1.0000         0
    0.4286    0.4474    1.0000   

U =

    7.0000    6.0000   10.0000
         0    5.4286    2.7143
         0         0   -0.5000

L 和 U 相乘得到：

>> L*U

ans =

     7     6    10
     3     8     7
     3     5     5

...这是原始矩阵A.

Answer 2

你可以使用这个 hack（尽管如前所述，你可能会失去数值稳定性）：

[L, U] = lu(sparse(A), 0)

Answer 3

您可能要考虑进行 LDU 分解而不是非透视 LU。看，没有旋转的 LU 在数值上是不稳定的——即使对于满秩和可逆的矩阵也是如此。上面提供的简单算法说明了原因 - 涉及矩阵的每个对角线元素。因此，如果对角线上的任何位置都为零，则分解失败，即使矩阵可能仍然是 non-singular.

维基百科在这里略微讨论了 LDU 分解：

https://en.wikipedia.org/wiki/LU_decomposition#LDU_decomposition

没有引用算法。它引用了以下教科书作为存在证明：

霍恩，罗杰 A.； Johnson, Charles R. (1985)，矩阵分析，剑桥大学出版社，ISBN 978-0-521-38632-6。请参阅第 3.5 节。

LDU是保证存在的（至少对于可逆矩阵而言），它是数值稳定的，也是唯一的（前提是L和U都被约束在对角线上有单位元素）。

然后，如果由于任何原因 "D" 妨碍您，您可以将对角矩阵 D 吸收到 L (L:=LD) 或 U (U:=D) U)，或者在 L 和 U 之间对称分割（例如 L:=L*sqrt(D) 和 U:=sqrt(D)*U），或者你想怎么做。 LDU 分解成 LU 的方法有无数种，这就是 LU 分解不唯一的原因。

在 MATLAB 中执行 LU 分解而不旋转

Perform LU decomposition without pivoting in MATLAB

matlab

matrix

linear-algebra