使用相对稀疏的矩阵（0.4％的数据）计算最小二乘法的最快方法

Question

上下文：我通过最小二乘法求解方程

A_des.x = b

，找到

x

。我必须执行此操作数千次。

A_des

每次迭代之间都会发生变化，但保持稀疏（始终约为矩阵中实际数据的 4%）。

问题：

A_des

很大（

(28106, 1185)

），求解线性系统大约需要1秒。我需要执行这个操作300000次。

问题： 您是否有什么可以加快此过程的线索？

我尝试过的：我在 GPU 和 CPU 上编写了我的函数。我使用

pytorch

进行 GPU 计算，但加速效果不一致。我无法并行化我的函数，因为求解器函数通常已经是线程化的。我尝试在CPU上使用多个求解器函数：

np.linalg.solve

、

np.linalg.lstsq

、

scipy.linalg.solve

、

scipy.optimize.least_squares

、

scipy.sparse.linalg.lsqr

（预先将

A_des

与

scipy.sparse.coo_array(A_des)

转换）。

我发现的最佳方法：我发现的最好的妥协是使用

np.linalg.solve

，我只需要手动进行一些矩阵操作。例如，使用最小二乘函数求解系统将为：

x = np.linalg.lstsq(A_des, b)

，而使用

solve

则为：

x = np.linalg.solve(A_des.T@A_des, A_des.T@b)

。我发现可以更快地解决系统问题的唯一其他选项是：

x = scipy.sparse.linalg.lsqr(A_des, b, atol = 1e-3, btol = 1e-3)

，但结果并不好。

我的功能：

def Inverter(GPU, A_des, b):

    if GPU:
        # Migrate b to torch
        b= torch.from_numpy(b).to(device).double()
        # Migrate design matrix to GPU
        A_des = torch.from_numpy(A_des).to(device)
        x= torch.linalg.solve(A_des.T@A_des,A_des.T@(b)).cpu()
    else:
        #x= scipy.sparse.coo_array(A_des)
        #x= scipy.sparse.linalg.lsqr(A_des, b)[0]
        #x = scipy.linalg.solve(A_des.T@A_des,A_des.T@(b))
        x= np.linalg.solve(A_des.T@A_des,A_des.T@b)
        

    return x

Answer 1

您可以尝试使用 Numba 库在 GPU 上并行化函数

https://numba.readthedocs.io/en/stable/user/jit.html