Solver time of linear equations in SciPy

Question

This is my python code for 4 functions solving the linear system equations:

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]

I know that first solution is not effecient, the second one is good if the number elements in b is small and solutions 3 and 4 are good if b is big.

But my tests reveal the opposite

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]

Output

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

Why is this the case?

Answer 1

You are getting those results because the dimensions of your matrices ( A & B) are way too small. With such small matrices the overhead of calling the functions and doing the required checks is essentially higher than the cost of actual computation. Practically speaking for a symmetric matrix of 2x2 inverse is always inevitably going to be faster than all the other aforementioned methods, simply because

This requires two checks: is the matrix Hermitian (square)? is the determinant non-zero?. This coupled with numpy.dot will certainly beat any other approach.

On the other hand, a call to scipy's standard solver will require more checks, decision on the appropriate solver, partial and/or full pivoting and so on; a call to Cholesky's solver would involve even more checks such as symmetricity and positive-definiteness of the matrix (eigen values being positive) and finally two calls one for Cholesky decomposition and one for triangular solver will go through many layers of CPython until it reaches the underlying Fortran/C lapack, whence a complete overkill for a 2x2 matrix. This is exactly what your results imply.

So increase the size of your matrix and you will get to see the actual merit of each approach

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

and here are the results on my machine for 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

Increasing size of A matrix to 1000 (i.e. n=1000), I get

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

By the way, your matrices are dense so you should not use for loops or list comprehensions, as the scipy's dense linear solver can handle matrices on the right hand side.