Матричная инверсия полосчатой ​​разреженной матрицы с использованием SciPy - PullRequest
Новая
       34

Матричная инверсия полосчатой ​​разреженной матрицы с использованием SciPy

0 голосов
/ 09 февраля 2019

Я пытаюсь решить инверсию полосчатой ​​разреженной матрицы наиболее эффективным способом, чтобы я мог включить ее в свою систему реального времени.Я генерирую разреженные матрицы, которые представляют операцию свертки.В настоящее время я использую 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

1 Ответ

0 голосов
/ 13 февраля 2019

Решил, используя 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
...