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

Question

Matlab 使用 sparse 命令计算对角矩阵的逆矩阵需要 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 如何反转矩阵？ 这里我们有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 分解 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

Answer 2

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)

这是 %time 输出：

inv: 2.39 s
splv: 193 ms

Answer 3

您正在从您的软件中提取关键信息：矩阵是对角线这一事实使其非常容易反转：您只需反转其对角线的每个元素：

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

当然，这仅对对角矩阵有效...

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 脚本进行比较。