使用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

- Maxtron

---

乍一看，solve_banded 函数所需的 ab 矩阵看起来很像 scipy.dia_matrix 格式的 .data 属性。两者都将对角线定义为填充的矩形数组（另一种选择是 sparse.diags 接受的不规则列表）。 - hpaulj

让我来看看。此外，似乎spsolve对于逆问题并没有给出正确的解决方案。我用np.linalg.inv检查了解决方案，误差似乎很大。 - Maxtron

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

- Maxtron