是否有可能像在Matlab中一样快地计算Python中的稀疏矩阵的逆？

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

这是答案.当您在matlab中为稀疏矩阵运行inv时,matlab会检查矩阵的不同属性以优化计算.对于稀疏对角矩阵,您可以运行以下代码来查看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如何反转矩阵？这里我们有代码:

...
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应该快得多.

编辑 可以肯定的是,正如@ P.Escondido所提议的那样,我们将尝试在矩阵A中进行微小的修改,以便在矩阵不是对角线时跟踪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