前向替换填充数组

Question

所以我试图创建一个脚本，将采取一个数组，将其转换成包装形式，然后运行对它的Cholesky分解，最后使用向前和向后替换来解决它(Ax=b)。

我遇到的麻烦是为打包数组创建前向替换(列方向)算法。我已经正确地完成了反向替换，但是我似乎无法正确地转换算法。

这是我目前的反向替换算法(这是完美的)

x = b;
p = n*(n+1)/2;
for j = n:-1:1
  x(j) = x(j) / u(p);
  p = p - 1;
  for i = j-1:-1:1
    x(i) = x(i) - u(p) * x(j);
    p = p - 1;
  end
end

其中n为原A矩阵的大小(即平方和SPD)。A矩阵到u的映射如下:

A(i,j) = u(i+j*(j-1)/2)

我目前的前向代换算法迭代如下:

x = b;
p = 1;
for j = 1:n
    x(j) = x(j) / u(p);
    p = p + 1;
    for i = j+1:n
        x(i) = x(i) - u(p) * x(j);
        p = p + 1;
    end
end

我似乎不知道我做错了什么。无论我尝试什么，最终只给我NaN或Inf作为x的答案。有没有比我聪明的人能算出算法应该是什么样的?

Answer 1

据我所知，你正试图解决导致Ax = b解决方案的两个问题。

其中A = LL'也就是A等于下三角矩阵(L)乘以它的转置矩阵(L')

首先用正向替换求解Ly = b。

最后通过反代入求解L'x = y。

对于您的反向替换，您的u(p)定义如下(对于n =3):

     | u(1) u(2) u(4) |
L' = |  0   u(3) u(5) |
     |  0    0   u(6) |

我认为你遇到的问题与前向替换来自L将被定义为:

    | u(1)  0    0   |
L = | u(2) u(3)  0   |
    | u(4) u(5) u(6) |

你的值p需要有以下的进位来完成你想要的前向替换:

{1, 2, 4, 3, 5, 6}

我已经修改了你当前的代码来实现这个模式。

x = b;
p = 1;
for j = 1:n
    d = j-1;
    p = j + (d^2+d)/2;
    x(j) = x(j) / u(p);
    d = d + 1;
    for i = j+1:n
        p = j + (d^2+d)/2;
        x(i) = x(i) - u(p) * x(j);
        d = d + 1;
    end
end

我能够对此进行测试，并且它产生了预期的结果。我希望这对你有帮助。