在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] 

我知道第一溶液是不effecient，第二個是良好的，如果數量的元件在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.46813352e-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倒數的對稱矩陣實際上來說總是不可避免地會產生比所有上述其他方法更快，這只是因爲

這需要兩次檢查：在所述矩陣厄米特（正方形）？是行列式非零？這加上numpy.dot肯定會擊敗任何其他方法。

另一方面，致電scipy的標準求解器將需要更多的檢查，決定適當的求解器，部分和/或全部樞軸等;對Cholesky求解器的調用將涉及更多的檢查，如對稱性和矩陣的正定性（本徵值爲正值），最後兩個調用一個用於Cholesky分解，一個用於三角求解器將經過許多層CPython直到它達到底層的Fortran/Clapack，從而完成對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 

順便說一句，你的矩陣是dense，所以你不應該使用for循環或list理解，因爲scipy的密集線性求解器可以處理右側的矩陣。