python - SciPy 中线性方程的求解器时间

Question

这是我的 4 个函数求解线性系统方程的 python 代码:

def inverse_solution(A, B):
    inv_A = scipy.linalg.inv(A)
    return [numpy.dot(inv_A, b) for b in B]

def scipy_standart_solution(A, B):
    return [scipy.linalg.solve(A, b) for b in B]

def cholesky_solution(A, B):
    K = scipy.linalg.cholesky(A, lower = True)
    t_K = K.T
    return [scipy.linalg.solve_triangular(t_K, scipy.linalg.solve_triangular(K, b, lower = True)) for b in B]

def scipy_cholesky_solution(A, B):
    K = scipy.linalg.cho_factor(A)
    return [scipy.linalg.cho_solve(K, b) for b in B]

我知道第一个解决方案效率不高，如果 b 中的数字元素很小，第二个解决方案很好，如果 b 很大，解决方案 3 和 4 很好.

但我的测试结果恰恰相反

A = numpy.array([[1,0.000000001],[0.000000001,1]])
for length in range(5):
    B = numpy.random.rand(length + 1,2)
    r_1 = %timeit -o -q inverse_solution(A,B)
    r_2 = %timeit -o -q scipy_standart_solution(A,B)
    r_3 = %timeit -o -q cholesky_solution(A, B)
    r_4 = %timeit -o -q scipy_cholesky_solution(A,B)
    print[r_1.best, r_2.best, r_3.best, r_4.best, length + 1]
print("////////////////////////////")
A = scipy.linalg.hilbert(12)
for length in range(5):
    B = numpy.random.rand(length + 1,12)
    r_1 = %timeit -o -q inverse_solution(A,B)
    r_2 = %timeit -o -q scipy_standart_solution(A,B)
    r_3 = %timeit -o -q cholesky_solution(A, B)
    r_4 = %timeit -o -q scipy_cholesky_solution(A,B)
    print[r_1.best, r_2.best, r_3.best, r_4.best, length + 1]

输出

[3.4187317043965046e-05, 4.0246286353324035e-05, 0.00010478259103165283, 5.702318102410473e-05, 1]
[5.8342068180991904e-05, 7.186096089739067e-05, 0.00015128208746822017, 8.083242991273707e-05, 2]
[3.369307648390851e-05, 9.87348262759614e-05, 0.00020628010956514232, 0.00012435874243853036, 3]
[3.79271575715201e-05, 0.00013554678863379266, 0.00027863672228316006, 0.0001421170191610699, 4]
[3.5385220680527143e-05, 0.00017145862333669016, 0.0003393085250830197, 0.00017440956920201814, 5]
////////////////////////////
[4.571321046813352e-05, 4.6984071999949605e-05, 9.794712904695047e-05, 5.9280641995266595e-05, 1]
[4.815794144123942e-05, 9.101890875345049e-05, 0.00017026901620170064, 9.290563584772826e-05, 2]
[4.604331714660219e-05, 0.00013565361678288213, 0.0002540085146054736, 0.00012585519183521684, 3]
[4.8241120284303915e-05, 0.0001758718917520369, 0.00031790739992344183, 0.00016162940724917405, 4]
[4.9397840771318616e-05, 0.00021475323253511647, 0.00037772389328304714, 0.00020302321951302815, 5]

为什么会这样？

Answer 1

您得到这些结果是因为矩阵(A 和 B)的维度太小了。对于如此小的矩阵，调用函数和执行所需检查的开销基本上高于实际计算的成本。实际上，对于 2x2 的对称矩阵逆总是不可避免地比所有其他上述方法更快，仅仅是因为

这需要两个检查:矩阵是 Hermitian(正方形)吗？ 行列式是否非零？这与 numpy.dot 相结合肯定会击败任何其他方法。

另一方面，调用 scipy 的标准求解器将需要更多检查、决定合适的求解器、部分和/或全部 旋转 等等;对 Cholesky 求解器的调用将涉及更多检查，例如矩阵的对称性 和正定性(特征值为正)，最后两次调用，一个用于 Cholesky 分解，另一个调用for triangular solver 将通过许多层 CPython 直到它到达底层 Fortran/C lapack，从那里完全矫枉过正 2x2 矩阵。这正是您的结果所暗示的。

因此增加矩阵的大小，您将看到每种方法的实际优点

n = 100
A = np.random.rand(n,n)
A = 0.5*(A+A.T)           # symmetrize
np.fill_diagonal(A,5)     # make sure matrix is positive-definite

B = numpy.random.rand(5,n)

%timeit inverse_solution(A,B)
%timeit scipy_standard_solution(A,B)
%timeit cholesky_solution(A, B)
%timeit scipy_cholesky_solution(A,B)
%timeit scipy.linalg.solve(A,B.T).T

下面是我机器上 n=100 的结果:

1000 loops, best of 3: 998 µs per loop # inverse
100 loops, best of 3: 2.21 ms per loop # standard solver with loops
1000 loops, best of 3: 972 µs per loop # decompose and solve lower tri
1000 loops, best of 3: 505 µs per loop # Cholesky 
1000 loops, best of 3: 530 µs per loop # standard solver w/o loops

将 A 矩阵的大小增加到 1000(即 n=1000)，我得到了

1 loops, best of 3: 736 ms per loop   # inverse
1 loops, best of 3: 1.24 s per loop   # standard solver with loops
1 loops, best of 3: 209 ms per loop   # decompose and solve lower tri
10 loops, best of 3: 171 ms per loop  # Cholesky
1 loops, best of 3: 254 ms per loop   # standard solver w/o loops

顺便说一句，你的矩阵是密集的，所以你不应该使用for 循环或list 理解，因为scipy 的密集线性求解器可以处理右侧的矩阵。