为什么旋转代码会给我的逻辑单元分解带来失败?

Question

在发言之前，我是和一个翻译一起写这篇文章的，而不是一个说英语的人。因此，即使我的英语很糟糕，我也希望你能理解。

在学习了如何编写非旋转的逻辑单元分解代码之后，我决定尝试编写一个使用旋转的逻辑单元分解代码作为作业。由于我的非旋转代码运行得很好，所以我认为制作旋转的lu分解代码会很容易。

但结果很糟，真的很糟…

% non-pivoting code
function x = mylu_np(A)
[m,~] = size(A);
L = eye(m); U = A(:,:);
for i = 1:m-1
for j = i+1:m
L(j,i) = U(j,i)/U(i,i);
U(j,i) = 0;
U(j,i+1:m) = U(j,i+1:m) - L(j,i)*U(i,i+1:m);
end
end
x = [L U];
end

我的代码如下请告诉我这个代码哪里出错了。我不知道为什么我只放了枢轴部分，却发生了如此严重的故障。

function x = mylu_pp(A)
[m,~] = size(A);
P = eye(m); L = eye(m); U = A(:,:);
for i = 1:m-1
for j = i+1:m
if abs(U(i,i)) < abs(U(j,i))
U([i j],:) = U([j i],:);
P([i j],:) = P([j i],:);
end
end
for j = i+1:m
L(j,i) = U(j,i)/U(i,i);
U(j,i) = 0;
U(j,i+1:m) = U(j,i+1:m) - L(j,i)*U(i,i+1:m);
end
end
x = [P L U];
end

我测试了4*4矩阵A，像这样

A = [0.0001 -5.0300 5.8090 7.8320;
2.2660 1.9950 1.2120 8.0080;
8.8500 5.6810 4.5520 1.3020;
6.7750 -2.2530 2.9080 3.9700];
W = mylu_pp(A)
P = W(1:4,1:4); L = W(1:4,5:8); U = W(1:4,9:12);
P*A - L*U

，答案是这样的。非常糟糕。

ans = [0.0000 0.0000 0.0000 0.0000;
6.7749 4.3489 3.4847 0.9967;
-2.2659 -7.0250 -1.6521 2.1754
-4.5090 2.6761 -1.8326 -3.1721]

告诉我哪个部分有问题。欢迎各种帮助。我将提前感谢你。

Answer 1

问题与矩阵L有关。

当交换矩阵U中的i和j行时，我们也必须交换矩阵L的i和j行中的元素，但只交换之前更新的元素(不是整个行):

if i > 1
L([i j], 1:i-1) = L([j i], 1:i-1);
end

---

我思考的方式相当于重命名方程系统中的变量。
我不是数学家，我的解释不会很好……

假设我们有方程组:A*v = b.
使用LU分解，我们将方程解为:L*U*v = b.

假设:

[x
v =  y
z
w]

在矩阵U中交换两行后，我们必须保持与原方程组等价。
保持等价的方法是重命名v中的元素。

当第2行和第4行交换时，我们可以交换y和w的名称，得到新的v:

[x
v0 =  w
z
y]

---

等价于名称交换:

• 交换矩阵U中的两行。
• 交换矩阵P中的两行。
  P*A将A中的行重新排序为新的"交换顺序"。
• 我们还必须交换矩阵L行中的元素，但只交换之前更新的元素(不是整行)。
  交换L的第2行和第4行相关元素的名称，与交换y和w的名称是等价的。

---

代码示例:

A = [0.0001 -5.0300 5.8090 7.8320;
2.2660  1.9950 1.2120 8.0080;
8.8500  5.6810 4.5520 1.3020;
6.7750 -2.2530 2.9080 3.9700];
[P, L, U] = mylu_pp(A);
P*A - L*U
function [P, L, U] = mylu_pp(A)
m = size(A, 1);
P = eye(m); L = eye(m); U = A;
for i = 1:m-1
for j = i+1:m
if abs(U(i,i)) < abs(U(j,i))
U([i j],:) = U([j i],:);
P([i j],:) = P([j i],:);

if i > 1
% Swap the relevant elements "1:i-1" in the row of matrix L (only the elements that were updated until now).
L([i j], 1:i-1) = L([j i], 1:i-1);
end
end
end
for j = i+1:m
L(j,i) = U(j,i)/U(i,i);
U(j,i) = 0;
U(j,i+1:m) = U(j,i+1:m) - L(j,i)*U(i,i+1:m);
end
end
%x = [P L U];
end

---

输出P*A - L*U:

ans =
0     0     0     0
0     0     0     0
0     0     0     0
0     0     0     0