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

At a first glance, the ab matrix that solve_banded wants looks a lot like the .data attribute of a scipy.dia_matrix format. Both define the diagonals as a padded rectangular array (the alternative being the ragged list that sparse.diags accepts). — hpaulj
Let me look into it. Also, it appears that spsolve does not give the right solution for the inverse problem. I checked the solution with np.linalg.inv and the error seems to be huge. — Maxtron

Maxtron Maxtron · Accepted Answer · 2019-02-13T05:17:02

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

Matrix Inversion of Banded Sparse Matrix using SciPy

1 Answers