稀疏的稀疏反转或溶解导致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的服务器上，仍然无济于事.
以前有人遇到过吗?
首先让我说，应该在 http://scicomp上更好地询问此问题.stackexchange.com ，那里有大量的计算科学和数值线性代数专家.
让我们从基础开始:从不反转稀疏矩阵，这完全没有意义.请参阅有关MATLAB Central的讨论，尤其是此<蒂姆·戴维斯(Tim Davis)的href ="http://blogs.mathworks.com/loren/2007/05/16/purpose-of-inv/#comment-16231" rel ="nofollow">评论. >
简而言之:没有用于对矩阵进行数值求逆的算法.每当您尝试对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系统时，最终得到对称正稀疏代数方程组.一旦知道了问题的性质，就可以选择正确的解决方案.
在您的情况下，您尝试使用通用直接求解器，而不对方程重新排序.众所周知，这将生成填充，破坏spsolve函数第一阶段中iBs矩阵的稀疏性(应该将其分解).请注意，全双精度150000 x 150000矩阵需要大约167 GB的内存.为了减少因式分解过程中的填充，有很多技术可以对方程进行重新排序，但是您没有提供足够的信息来给您一个明智的提示.
很抱歉，但是您应该考虑在 http://scicomp.stackexchange.com 上重新提出问题，明确指出哪个是您要解决的问题，以便提供有关矩阵结构和属性的线索.
I am trying to invert a large (150000,150000) sparse matrix as follows:
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)
leads to the following error message:
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
I rejigged the program to simply solve a system of linear differential equations:
import numpy as np
N=Bs.shape[0]
I=np.ones(N)
M=splu.spsolve(Bs,I)
and I encounter the same error again
I was using this code on a machine with 16 GB of RAM and then moved it onto a server with 32 GB of RAM, still to no avail.
has anyone encountered this before?
First let me say that this question should be better asked on http://scicomp.stackexchange.com where there is a great community of experts in computational science and numerical linear algebra.
Let's start from the basics: never invert a sparse matrix, it's completely meaningless. See this discussion on MATLAB central and particularly this comment by Tim Davis.
Briefly: there are no algorithms for numerically inverting a matrix. Whenever you try to numerically compute the inverse of a NxN matrix, you solve in fact N linear systems with N rhs vectors corresponding to the columns of the identity matrix.
In other words, when you compute
from scipy.sparse import eye
from scipy.sparse.linalg import (inv, spsolve)
N = Bs.shape[0]
iBs = inv(Bs)
iBs = spsolve(Bs, eye(N))
the last two statements (inv(eye) and spsolve(Bs, eye(N))) are equivalent. Please note that the identity matrix (eye(N)) is not a ones vector (np.ones(N)) as you question falsely assumes.
The point here is that matrix inverses are seldom useful in numerical linear algebra: the solution of Ax = b is not computed as inv(A)*b, but by a specialised algorithm.
Going to your specific problem, for big sparse system of equations there are no black-box solvers. You can chose the correct class of solvers only if you have a good understanding of the structure and properties of your matrix problem. The properties of your matrices in turn are a consequence of the problem you are trying to solve. E.g. when you discretise by the FEM a system of elliptic PDE, you end up with a symmetric positive sparse system of algebraic equations. Once you know the properties of your problem, you can choose the correct solving strategy.
In your case, you are trying to use a generic direct solver, without reordering the equations. It is well known that this will generate fill-ins which destroy the sparsity of the iBs matrix in the first phase of the spsolve function (which should be a factorisation.) Please note that a full double precision 150000 x 150000 matrix requires about 167 GB of memory. There are a lot of techniques for reordering equations in order to reduce the fill-in during factorisation, but you don't provide enough info for giving you a sensible hint.
I'm sorry, but you should considering reformulating your question on http://scicomp.stackexchange.com clearly stating which is the problem you are trying to solve, in order to give a clue on the matrix structure and properties.
这篇关于稀疏的稀疏反转或溶解导致UMFPACK_ERROR_OUT_OF_MEMORY的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
