Matrix Inversion of Banded Sparse Matrix using SciPy

Question

I am trying to solve the inverse of a banded sparse matrix in the most efficient way so that I can incorporate this in my real-time system. I am generating sparse-banded matrices which represent a convolution operation. Currently, I am using spsolve from scipy.sparse.linalg library. I found that there is a better way by using solve_banded from the scipy.linalg library. However, solve_banded requires (l,u) which is the number of non-zero lower and upper diagonals and ab which (l + u + 1, M) array like banded matrix. I am not sure how to convert my code so that I can use solve_banded. Any help with this regard is highly appreciated.

import numpy as np
from scipy import linalg
import math
import time

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


def ABC(deg, fc, N):
    r"""Generate sparse-banded matrices
    """

    omc = 2*math.pi*fc
    t = ((1-math.cos(omc))/(1+math.cos(omc)))**deg

    p = 1
    for k in np.arange(deg):
        p = np.convolve(p,np.array([-1,1]),'full')
    P = spdiags(np.kron(p,np.ones((N,1))).T, np.arange(deg+1), N-deg, N)
    B = P.T.dot(P)

    q = np.sqrt(t)
    for k in np.arange(deg):
        q = np.convolve(q,np.array([1,1]),'full')
    Q = spdiags(np.kron(q,np.ones((N,1))).T, np.arange(deg+1), N-deg, N)

    C = Q.T.dot(Q)
    A = B + C

    return A,B,C

if __name__ == '__main__':

    mu = 0.1
    deg = 3
    wc = 0.1

    for i in np.arange(1,7,1):

        # some dense random vector
        x = np.random.rand(10**i,1)
        # generate sparse banded matrices
        A,_,C = ABC(deg, wc, 10**i)
        # another banded matrix
        G = mu*A.dot(A.T) + C.dot(C.T)

        # SCIPY SPSOLVE
        st = time.time()
        y = spsolve(G,x)
        et = time.time()
        print("SCIPY SPSOLVE: N = ", 10**i, "Time taken: ", et-st)

Results

SCIPY SPSOLVE: N =  10 Time taken:  0.0
SCIPY SPSOLVE: N =  100 Time taken:  0.0
SCIPY SPSOLVE: N =  1000 Time taken:  0.015689611434936523
SCIPY SPSOLVE: N =  10000 Time taken:  0.020943641662597656
SCIPY SPSOLVE: N =  100000 Time taken:  0.16722917556762695
SCIPY SPSOLVE: N =  1000000 Time taken:  1.7254831790924072

Answer 1

Solved it using solveh_banded from the scipy library. Very fast matrix inversion technique for extremely large sparse-banded matrices when the matrix is symmetric and positive definite banded matrix.

from scipy.linalg import solveh_banded

def sp_inv(A, x):

    A = A.toarray()
    N = np.shape(A)[0]
    D = np.count_nonzero(A[0,:])
    ab = np.zeros((D,N))
    for i in np.arange(1,D):
        ab[i,:] = np.concatenate((np.diag(A,k=i),np.zeros(i,)),axis=None)
    ab[0,:] = np.diag(A,k=0)
    y = solveh_banded(ab,x,lower=True)
    return y