向量化（手动）正向替换

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 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)

这是有关循环矩阵函数的更多信息。

现在，将所有内容整合在一起：

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应使用的下三角部分。