用SciPy进行带状稀疏矩阵的矩阵求逆

Question

我试图以最有效的方式解决带状稀疏矩阵的逆，以便我可以将其合并到我的实时系统中。我正在生成代表卷积运算的稀疏带状矩阵。目前，我正在使用spsolve图书馆的scipy.sparse.linalg。我发现使用solve_banded库中的scipy.linalg有更好的方法。然而，solve_banded需要(l,u)，这是非零的上下对角线和ab的数量，(l + u + 1, M)阵列像带状矩阵。我不知道如何转换我的代码，以便我可以使用solve_banded。对此方面的任何帮助都非常感谢。

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)

结果

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

使用solveh_banded库中的scipy解决了这个问题。当矩阵是对称的和正定带状矩阵时，非常快速的矩阵求逆技术用于极大的稀疏带状矩阵。

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