Is it possible to compute an inverse of sparse matrix in Python as fast as in Matlab?

Question

It takes 0.02 seconds for Matlab to compute the inverse of a diagonal matrix using the sparse command.

P = diag(1:10000);
P = sparse(P);
tic;
A = inv(P);
toc

However, for the Python code it takes forever - several minutes.

import numpy as np
import time

startTime = time.time()
P = np.diag(range(1,10000))
A = np.linalg.inv(P)
runningTime = (time.time()-startTime)/60
print "The script was running for %f minutes" % runningTime

I tried to use Scipy.sparse module but it did not help. The running time dropped, but only to 40 seconds.

import numpy as np
import time
import scipy.sparse as sps
import scipy.sparse.linalg as spsl

startTime = time.time()
P = np.diag(range(1,10000))
P_sps = sps.coo_matrix(P)
A = spsl.inv(P_sps)
runningTime = (time.time()-startTime)/60
print "The script was running for %f minutes" % runningTime

Is it possible to run the code as fast as it runs in Matlab?

Answer 1

Here is the answer. When you run inv in matlab for a sparse matrix, matlab check different properties of the matrix to optimize the calculation. For a sparse diagonal matrix, you can run the followin code to see what is matlab doing

n = 10000;
a = diag(1:n);
a = sparse(a);
I = speye(n,n);
spparms('spumoni',1);
ainv = inv(a);
spparms('spumoni',0);

Matlab will print the following:

sp\: bandwidth = 0+1+0.
sp\: is A diagonal? yes.
sp\: do a diagonal solve.

So matlab is inverting only the diagonal.

How does Scipy invert the matrix?? Here we have the code:

...
from scipy.sparse.linalg import spsolve
...

def inv(A):
    """
    Some comments...
    """
    I = speye(A.shape[0], A.shape[1], dtype=A.dtype, format=A.format)
    Ainv = spsolve(A, I)
    return Ainv

and spsolve

    # Cover the case where b is also a matrix
    Afactsolve = factorized(A)
    tempj = empty(M, dtype=int)
    x = A.__class__(b.shape)
    for j in range(b.shape[1]):
        xj = Afactsolve(squeeze(b[:, j].toarray()))
        w = where(xj != 0.0)[0]
        tempj.fill(j)
        x = x + A.__class__((xj[w], (w, tempj[:len(w)])),
                            shape=b.shape, dtype=A.dtype)

i.e., scipy factorize A and then solve a set of linear systems where the right hand sides are the coordinate vectors (forming the identity matrix). Ordering all the solutions in a matrix we obtain the inverse of the initial matrix.

If matlab is taken advantage of the diagonal structure of the matrix, but scipy is not (of course scipy is also using the structure of the matrix, but in a less efficient way, at least for the example), matlab should be much faster.

EDIT To be sure, as @P.Escondido propossed, we will try a minor modification in matrix A, to trace the matlab procedure when the matrix is not diagonal:

n = 10000; a = diag(1:n); a = sparse(a); ainv = sparse(n,n);
spparms('spumoni',1);
a(100,10) = 500; a(10,1000) = 200; 
ainv = inv(a);
spparms('spumoni',0);

It prints out the following:

sp\: bandwidth = 90+1+990.
sp\: is A diagonal? no.
sp\: is band density (0.00) > bandden (0.50) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? yes.
sp\: permute and solve.
sp\: sprealloc in sptsolve: 10000 10000 10000 15001

Answer 2

How about splu(), it's faster but need a dense array and return dense array:

Create a random matrix:

import numpy as np
import time
import scipy.sparse as sps
import scipy.sparse.linalg as spsl
from numpy.random import randint
N = 1000
i = np.arange(N)
j = np.arange(N)
v = np.ones(N)

i2 = randint(0, N, N)
j2 = randint(0, N, N)
v2 = np.random.rand(N)

i = np.concatenate((i, i2))
j = np.concatenate((j, j2))
v = np.concatenate((v, v2))

A = sps.coo_matrix((v, (i, j)))
A = A.tocsc()

%time B = spsl.inv(A)

calculate inverse matrix by splu():

%%time
lu = spsl.splu(A)
eye = np.eye(N)
B2 = lu.solve(eye)

check the result:

np.allclose(B.todense(), B2.T)

Here is the %time output:

inv: 2.39 s
splv: 193 ms