向量化（手动）正向替换

Question

我正在手写代码（我知道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 

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

Answer 1

方程的解Ly=b可以通过将矩阵求逆L并在两边都左乘而得出y=L'b，从而得到：，其中L'是的逆矩阵L。可以使用例如来对这个矩阵求逆np.linalg.inv。在不使用numpy的情况下进行矩阵求逆将是乏味的。但是，我怀疑您可能会做得很好，因为您有一个较低的三角形单位矩阵。下三角单元矩阵的逆L（1对角线也带有s）可以表示为def Linv(i,j):    if i==j:        return 1    elif j==i-1:        return -1可能有更好的方法来计算它，但这是一种方法：import numpy as npfrom scipy.linalg import circulantL=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)这是有关循环矩阵函数的更多信息。现在，将所有内容整合在一起：import numpy as npfrom scipy.linalg import circulantdef 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应使用的下三角部分。