Scipy 稀疏反转或 spsolve 导致 UMFPACK_ERROR_OUT_OF_MEMORY

Scipy sparse invert or spsolve lead to UMFPACK_ERROR_OUT_OF_MEMORY

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

有没有人遇到过这种情况?

稀疏数组只能将矩阵的 non-zero 个条目放入内存。现在假设你做一个反转。这意味着矩阵的几乎所有条目都变成 non-zero。稀疏矩阵是内存优化的。

有一些操作可以应用于稀疏矩阵而不会丢失 "spare" 属性:

  • 加法,只要加一个常数就可以保持稀疏矩阵的稀疏性。

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

让我们从基础开始:永远不要反转稀疏矩阵,这是完全没有意义的。请参阅 Tim Davis 的 discussion on MATLAB central and particularly this comment

简而言之:没有用于对矩阵进行数值求逆的算法。每当您尝试以数值方式计算 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 计算的,而是通过专门的算法计算的。

针对您的具体问题,对于大型稀疏方程组,没有 black-box 求解器。只有充分了解矩阵问题的结构和性质,才能选择正确的 class 求解器。矩阵的属性反过来又是您要解决的问题的结果。例如。当您通过 FEM 离散化椭圆 PDE 系统时,您最终会得到代数方程的对称正稀疏系统。了解问题的属性后,您就可以选择正确的解决策略。

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

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