如何在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。
发布于 2016-12-14 20:02:41
默认情况下,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
。
发布于 2019-04-17 15:12:27
你可以使用这个技巧(尽管如前所述,你可能会失去数值稳定性):
[L, U] = lu(sparse(A), 0)
发布于 2019-05-31 20:39:36
您可能想要考虑进行LDU分解,而不是未旋转的LU。看,没有旋转的LU在数字上是不稳定的-即使对于满秩和可逆的矩阵也是如此。上面提供的简单算法说明了原因-有除以涉及的矩阵的每个对角元素。因此,如果对角线上的任何地方都有零,分解就会失败,即使矩阵可能仍然是非奇异的。
Wikipedia在这里谈了一点LDU分解:
https://en.wikipedia.org/wiki/LU_decomposition#LDU_decomposition
而不引用算法。它引用了以下教科书作为存在的证明:
罗杰·霍恩;约翰逊·查尔斯·R. (1985年),矩阵分析,剑桥大学出版社,ISBN 978-0-521-38632-6。请参见第3.5节。
LDU被保证存在(至少对于可逆矩阵),它在数值上是稳定的,并且它也是唯一的(假设L和U都被约束为在对角线上具有单位元素)。
然后,如果由于任何原因"D“挡住了你的路,你可以将对角矩阵D吸收到L (L:=L_D)或U (U:=D_U)中,或者在L和U之间对称地分割它(比如L:=L*sqrt(D)和U:=sqrt(D)*U),或者你想怎么做都行。有无数种方法可以将LDU拆分成LU,这就是为什么LU分解不是唯一的。
https://stackoverflow.com/questions/41150997
复制