Можно ли вычислить инверсию разреженной матрицы в Python так же быстро, как в Matlab?

Question

Для Matlab требуется 0,02 секунды для вычисления инверсии диагональной матрицы с использованием разреженной команды.

P = diag(1:10000); 
P = sparse(P); 
tic; 
A = inv(P); 
toc 

Однако для кода Python требуется навсегда - несколько минут.

import numpy as np 
import time 

startTime = time.time() 
P = np.diag(range(1,10000)) 
A = np.linalg.inv(P) 
runningTime = (time.time()-startTime)/60 
print "The script was running for %f minutes" % runningTime 

Я попытался использовать модуль Scipy.sparse, но это не помогло. Время работы снизилось, но только до 40 секунд.

import numpy as np 
import time 
import scipy.sparse as sps 
import scipy.sparse.linalg as spsl 

startTime = time.time() 
P = np.diag(range(1,10000)) 
P_sps = sps.coo_matrix(P) 
A = spsl.inv(P_sps) 
runningTime = (time.time()-startTime)/60 
print "The script was running for %f minutes" % runningTime 

Возможно ли запустить код так быстро, как он работает в Matlab?

Answer 1

Вот ответ. Когда вы запускаете inv в matlab для разреженной матрицы, matlab проверяет различные свойства матрицы для оптимизации вычисления. Для редкой диагональной матрицы, вы можете запустить followin код, чтобы увидеть, что делает MATLAB

n = 10000; 
a = diag(1:n); 
a = sparse(a); 
I = speye(n,n); 
spparms('spumoni',1); 
ainv = inv(a); 
spparms('spumoni',0); 

Matlab напечатает следующее:

sp\: bandwidth = 0+1+0. 
sp\: is A diagonal? yes. 
sp\: do a diagonal solve. 

Так MATLAB является инвертирование только по диагонали.

Как Scipy инвертирует матрицу? Здесь мы имеем code:

... 
from scipy.sparse.linalg import spsolve 
... 

def inv(A): 
    """ 
    Some comments... 
    """ 
    I = speye(A.shape[0], A.shape[1], dtype=A.dtype, format=A.format) 
    Ainv = spsolve(A, I) 
    return Ainv 

и spsolve

# Cover the case where b is also a matrix 
    Afactsolve = factorized(A) 
    tempj = empty(M, dtype=int) 
    x = A.__class__(b.shape) 
    for j in range(b.shape[1]): 
     xj = Afactsolve(squeeze(b[:, j].toarray())) 
     w = where(xj != 0.0)[0] 
     tempj.fill(j) 
     x = x + A.__class__((xj[w], (w, tempj[:len(w)])), 
          shape=b.shape, dtype=A.dtype) 

т.е., scipy factorize A, а затем решить набор линейных систем, где правые части являются координатными векторами (образующими единичную матрицу). Заказывая все решения в матрице, получим обратную исходную матрицу.

Если Matlab используется диагональной структурой матрицы, но scipy не является (конечно, scipy также использует структуру матрицы, но менее эффективным способом, по крайней мере, для примера), Matlab должен быть намного быстрее.

EDIT Чтобы быть уверенным, что propossed @ P.Escondido, мы будем пробовать незначительные изменения в матрице А, чтобы проследить процедуру MatLab, когда матрица не является диагональной:

n = 10000; a = diag(1:n); a = sparse(a); ainv = sparse(n,n); 
spparms('spumoni',1); 
a(100,10) = 500; a(10,1000) = 200; 
ainv = inv(a); 
spparms('spumoni',0); 

Он печатает из следующих действий:

sp\: bandwidth = 90+1+990. 
sp\: is A diagonal? no. 
sp\: is band density (0.00) > bandden (0.50) to try banded solver? no. 
sp\: is A triangular? no. 
sp\: is A morally triangular? yes. 
sp\: permute and solve. 
sp\: sprealloc in sptsolve: 10000 10000 10000 15001 

Answer 2

Вы Witholding важной информации из вашего программного обеспечения: тот факт, что матрица является диагональной делает его очень легко инвертировать: вы просто инвертировать каждый элемент диагонали:

P = np.diag(range(1,10000)) 
A = np.diag(1.0/np.arange(1,10000)) 

Конечно, это справедливо только для диагональных матриц ...

Answer 3

Как насчет splu(), это быстрее, но нужен плотный массив и вернуться плотный массив:

Создать случайную матрицу:

import numpy as np 
import time 
import scipy.sparse as sps 
import scipy.sparse.linalg as spsl 
from numpy.random import randint 
N = 1000 
i = np.arange(N) 
j = np.arange(N) 
v = np.ones(N) 

i2 = randint(0, N, N) 
j2 = randint(0, N, N) 
v2 = np.random.rand(N) 

i = np.concatenate((i, i2)) 
j = np.concatenate((j, j2)) 
v = np.concatenate((v, v2)) 

A = sps.coo_matrix((v, (i, j))) 
A = A.tocsc() 

%time B = spsl.inv(A) 

вычислить обратную матрицу с помощью splu():

%%time 
lu = spsl.splu(A) 
eye = np.eye(N) 
B2 = lu.solve(eye) 

проверка результат:

np.allclose(B.todense(), B2.T) 

Вот выход% Время:

inv: 2.39 s 
splv: 193 ms 

Answer 4

Если вы пытаетесь с этим результатом будет лучше:

import numpy as np 
import time 
import scipy.sparse as sps 
import scipy.sparse.linalg as spsl 

P = np.diag(range(1,10000)) 
P_sps = sps.coo_matrix(P) 
startTime = time.time() 
A = spsl.inv(P_sps) 
runningTime = (time.time()-startTime)/60 
print "The script was running for %f minutes" % runningTime 

Теперь вы можете сравнить с вашим скриптом matlab.