14

我正在尝试反转一个大的(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 的服务器上，但仍然无济于事。

有没有人遇到过这个？

4

2 回答 2

4

首先让我说，这个问题最好在http://scicomp.stackexchange.com上提出，那里有一个庞大的计算科学和数值线性代数专家社区。

让我们从基础开始：永远不要反转稀疏矩阵，它完全没有意义。请参阅有关 MATLAB 中心的讨论，尤其是 Tim Davis 的此评论。

简而言之：没有用于数值反转矩阵的算法。每当您尝试对 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 计算的，而是通过专门的算法计算的。

针对您的具体问题，对于大型稀疏方程组，没有黑盒求解器。只有充分了解矩阵问题的结构和属性，才能选择正确的求解器类别。反过来，矩阵的属性是您试图解决的问题的结果。例如，当您通过 FEM 离散椭圆 PDE 系统时，您最终会得到一个对称的正稀疏代数方程系统。一旦了解了问题的性质，就可以选择正确的解决策略。

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

抱歉，您应该考虑在http://scicomp.stackexchange.com上重新表述您的问题，明确说明您要解决的问题，以便为矩阵结构和属性提供线索。

于 2016-01-23T23:26:13.757 回答
0

稀疏数组仅将矩阵的非零条目放入内存。现在假设你做一个反转。这意味着矩阵的几乎所有条目都变为非零。稀疏矩阵是内存优化的。

您可以在不丢失“备用”属性的情况下对稀疏矩阵应用一些操作：

  • 另外，只要加上一个常数就可以保持稀疏矩阵的稀疏。
于 2016-01-21T12:35:30.613 回答