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 中一样快地运行代码?
最佳答案
答案在这里。当您在 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
# 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
关于python - 是否有可能像在 Matlab 中一样快地在 Python 中计算稀疏矩阵的逆矩阵?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/21693613/