python - 向量化(手动)正向替换

我在手写代码（我知道NumPy可以用np.linalg.solve或类似的代码为我解决这个问题）来解决线性系统。我要编写的功能之一是正向替换-即为Ly=b求解y，其中L是单位下三角矩阵，而b是列向量。

我想出了以下解决方案

def solve_forward(L, b):
    y = b.copy()
    for r in range(1, L.shape[0]):
        y[r] -= L[r, :r] @ y[:r]
    return y

我想知道进行某种减法累加来删除循环是否可行，或者这看起来是否像“矢量化”一样。

最佳答案

方程Ly=b的解可以通过对矩阵L求逆并在两边进行左乘运算得到，从而得到：y=L'b，其中L'是L的逆矩阵。可以用例如np.linalg.inv。

在不使用numpy的情况下进行矩阵求逆将是乏味的。但是，我怀疑您可能会做得很好，因为您有一个较低的三角形单位矩阵。

下三角单元矩阵L的逆（对角线也有1）可以表示为

def Linv(i,j):
    if i==j:
        return 1
    elif j==i-1:
        return -1

可能有更好的方法来计算它，但这是一种方法：

import numpy as np
from scipy.linalg import circulant
L=np.tril(np.ones((4,4)))
dim=L.shape[0]
Linv=[np.concatenate([np.array([1,-1]), np.zeros(dim-2)])]
Linv=np.tril(circulant(Linv))
print(Linv)

这是有关circulant matrix function的更多信息。

现在，将它们放在一起：

import numpy as np
from scipy.linalg import circulant

def L_inv(l_dim):
    Linv=[np.concatenate([np.array([1,-1]), np.zeros(l_dim-2)])]
    Linv=tril(circulant(Linv))

def solve_forward(L, b):
    y = L_inv(L.shape[0]) @ b
    return y

哪个应该按预期工作。

编辑：在这种情况下，以前的toeplitz将不起作用。用更合适的circulant进行切换。

编辑2：仅应使用circulant的下部三角形部分。

关于python - 向量化(手动)正向替换，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/50977894/