枕骨稀疏倒转或散斑导致UMFPACK_ERROR_OUT_OF_MEMORY

Question

我试图将一个大型(150000,150000)稀疏矩阵反演如下：

import scipy as sp
import scipy.sparse.linalg as splu

#Bs is a large sparse matrix with shape=(150000,150000)

#calculating the sparse inverse
iBs=splu.inv(Bs)

导致以下错误消息：

Traceback (most recent call last):
    iBs=splu.inv(Bs)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/linsolve.py", line 134, in spsolve
autoTranspose=True)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/umfpack/umfpack.py", line 603, in linsolve
self.numeric(mtx)
  File "/usr/lib/python2.7/dist-packages/scipy/sparse/linalg/dsolve/umfpack/umfpack.py", line 450, in numeric
umfStatus[status]))
RuntimeError: <function umfpack_di_numeric at 0x7f2c76b1d320> failed with UMFPACK_ERROR_out_of_memory

我对程序进行了重新调整，以简单地求解一个线性微分方程组：

import numpy as np

N=Bs.shape[0]

I=np.ones(N)

M=splu.spsolve(Bs,I)

我再次遇到同样的错误

我在一台拥有16 GB RAM的机器上使用了这段代码，然后将它移动到32 GB RAM的服务器上，仍然没有效果。

以前有人遇到过这种情况吗？

Answer 1

首先，让我说，这个问题应该更好地问在http://scicomp.stackexchange.com上，那里有一个伟大的计算科学和数值线性代数的专家群体。

让我们从基础开始：从不反转稀疏矩阵，它是完全没有意义的。请参阅MATLAB上的这个讨论，特别是Tim的这个评论。

简单地说:没有任何算法来数值反演一个矩阵。当你试图数值计算一个NxN矩阵的逆时，你实际上解了N个线性系统，它的N个rhs向量对应于恒等矩阵的列。

换句话说，当你计算

from scipy.sparse import eye
from scipy.sparse.linalg import (inv, spsolve)

N = Bs.shape[0]
iBs = inv(Bs)
iBs = spsolve(Bs, eye(N))

最后两个语句(inv(eye)和spsolve(Bs, eye(N)))是等价的。请注意，恒等矩阵(eye(N))是，而不是，这是一个向量(np.ones(N))，因为您的疑问是错误的假设。

这里的要点是矩阵逆在数值线性代数中很少有用: Ax =b的解不是作为inv(A)*b计算的，而是由一个专门的算法计算的。

关于你的具体问题，对于大型稀疏方程组，没有黑匣子求解者。只有很好地理解矩阵问题的结构和性质，才能选择正确的求解类。矩阵的性质反过来又是你想要解决的问题的结果。例如，当你用有限元法离散一个椭圆PDE系统时，你会得到一个对称的正稀疏代数方程组。一旦你知道了问题的性质，你就可以选择正确的解决策略。

在您的情况下，您尝试使用通用直接求解器，而不重新排序方程。众所周知，这将生成填充，在iBs函数的第一阶段(应该是因式分解)破坏spsolve矩阵的稀疏性。请注意，一个完整的双精度150000 x 150000矩阵需要大约167 GB的内存。有许多重新排序方程的技术，以减少在分解过程中的填充，但您没有提供足够的信息来给您一个明智的提示。

对不起，您应该考虑在http://scicomp.stackexchange.com上重新定义您的问题，明确说明您要解决的问题是什么，以便给出矩阵结构和性质方面的线索。

Answer 2

稀疏数组只适合内存中矩阵的非零项。现在假设你做了一个反演。这意味着矩阵的几乎所有条目都变成非零。稀疏矩阵是内存优化的。

有些操作可以应用于稀疏矩阵，而不会丢失“备用”属性：

• 另外，只需增加一个常数就可以保持稀疏矩阵的稀疏性。