使用 SciPy 的带状稀疏矩阵的矩阵求逆

Question

我正在尝试以最有效的方式求解带状稀疏矩阵的逆矩阵，以便将其合并到我的实时系统中。我正在生成代表卷积运算的稀疏带状矩阵。目前，我正在使用 scipy.sparse.linalg 库中的 spsolve。我发现使用 scipy.linalg 库中的 solve_banded 是一种更好的方法。但是，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

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

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

使用 SciPy 的带状稀疏矩阵的矩阵求逆

Matrix Inversion of Banded Sparse Matrix using SciPy

python

matrix

scipy

sparse-matrix

matrix-inverse