为什么在求解稀疏线性方程组时会出现内存错误？

Question

当尝试从下方求解大型稀疏线性方程组时，我只得到了 MemoryError:。我该如何解决这个问题？

此外，此代码基于 Matlab 中的实现，应该运行没问题。原版中，M是一个三维矩阵，不知道是不是因为我把M修改成了二维的scipy.sparse.lil_matrix（而不是coo) s.t。我可以迭代地填写 M in.

def dectrans(features, faces, template):
    """
    Decode transformation matrix.

    features : shape (3*N) numpy.ndarray
    faces : (Nx3) array
    template : (Nx3) array
    """
    ftrs = features

    weight = 1

    fixvertex = 1
    fixto = np.zeros((3))

    # M shape originally (len(template), len(template), 10 * len(template))
    M = scipy.sparse.lil_matrix((len(template)+fixvertex, len(template) * 10 * len(template)))
    dx = scipy.sparse.lil_matrix((len(template)+fixvertex,3))

    # build laplacian system
    for i in range(len(faces)):
        v = faces[i,:]
        ...
        M[v][:,v] = M[v][:,v] + WIJ # WIJ some 3x3 matrix
        dx[v,:] = dx[v,:] + WIJ.dot(x.T)

    weight = np.ones((fixvertex)) * weight

    for i in range(fixvertex):
        M[len(template)+i, fixvertex-1] = weight[i]

    dx[len(template):len(template),:] = fixto.dot(np.tile(weight, (3))) 

    M = np.real(M)
    dx = np.real(dx)
    Mt = M.T
    model = scipy.sparse.linalg.spsolve(Mt @ M, Mt.dot(dx)) # here I get the error

    return model

这是我得到的错误的回溯：

MemoryError                               Traceback (most recent call last)
<ipython-input-10-9aa6e73eb179> in <module>
     20         rr = encrelrot(v_positions, faces, r_v_positions, f_neighbors)
     21 
---> 22         modelout = dectrans(decrelrot(rr, f_neighbors), faces, r_v_positions)

<ipython-input-8-cdb51dd3cadf> in dectrans(features, faces, template)
    616     print("Size dx", dx.nnz)
    617     #M = M.tocsr()
--> 618     model = scipy.sparse.linalg.spsolve(Mt @ M, Mt.dot(dx))
    619 
    620     return model

~/anaconda3/lib/python3.6/site-packages/scipy/sparse/base.py in __matmul__(self, other)
    560             raise ValueError("Scalar operands are not allowed, "
    561                              "use '*' instead")
--> 562         return self.__mul__(other)
    563 
    564     def __rmatmul__(self, other):

~/anaconda3/lib/python3.6/site-packages/scipy/sparse/base.py in __mul__(self, other)
    480             if self.shape[1] != other.shape[0]:
    481                 raise ValueError('dimension mismatch')
--> 482             return self._mul_sparse_matrix(other)
    483 
    484         # If it's a list or whatever, treat it like a matrix

~/anaconda3/lib/python3.6/site-packages/scipy/sparse/compressed.py in _mul_sparse_matrix(self, other)
    494                                      other.indptr, other.indices),
    495                                     maxval=M*N)
--> 496         indptr = np.empty(major_axis + 1, dtype=idx_dtype)
    497 
    498         fn = getattr(_sparsetools, self.format + '_matmat_pass1')

MemoryError:

Answer 1

回溯应该显示问题是在 spsolve 中还是在创建一个或两个参数时，Mt@M 或 Mt.dot(dx)。

具有 M 和 dx 形状 ((6891, 474721000)、(6891, 3)

 Mt@M
 (474721000,6891) + (6891, 474721000) => (474721000, 474721000)
 Mt.dot(dx)   # why not Mt@dx?
 (474721000,6891) + (6891, 3) => (474721000, 3)

根据这些非零值的结构，@product 的非零值可能比 M 多得多，这可能会产生内存错误。一个或另一个的回溯可能有助于我们对此进行诊断。

更常见的内存错误是由于试图从稀疏数组创建密集数组而导致的，但这里似乎确实是这种情况。但同样，追溯可以帮助排除这种情况。

如果增量填充矩阵值，建议使用

lil 格式。 csr 用于矩阵乘积，但如果需要，sparse 很容易将 lil 转换为 csr。所以这应该不是问题。

===

创建一个包含 1 个非零元素的稀疏矩阵：

In [263]: M=sparse.lil_matrix((1000,100000))                                 
In [264]: M                                                                  
Out[264]: 
<1000x100000 sparse matrix of type '<class 'numpy.float64'>'
    with 0 stored elements in LInked List format>
In [265]: M[0,0]=1                                                           
In [266]: M                                                                  
Out[266]: 
<1000x100000 sparse matrix of type '<class 'numpy.float64'>'
    with 1 stored elements in LInked List format>

这个 @ 没有产生内存错误，并且结果只有 1 个非零项，正如预期的那样。但是运行这个有一个明显的延迟，表明它正在做一些大的计算：

In [267]: M.T@M                                                              
Out[267]: 
<100000x100000 sparse matrix of type '<class 'numpy.float64'>'
    with 1 stored elements in Compressed Sparse Row format>

在 csr 等效项上执行相同的 @ 没有时间延迟：

In [268]: M1=M.tocsr()                                                       
In [269]: M1.T@M1                                                            
Out[269]: 
<100000x100000 sparse matrix of type '<class 'numpy.float64'>'
    with 1 stored elements in Compressed Sparse Column format>

===

您在 MATLAB 中提到了 3d 稀疏矩阵。您必须使用某种第 3 方扩展或解决方法，因为 MATLAB 稀疏仅限于 2d（至少我几年前使用它进行 FEM 工作时是这样）。

scipy.sparsecsc格式类似于MATLAB的内部稀疏。事实上，如果您通过 save 和 scipy.io.loadmat 传输矩阵，您就会得到这样的结果。 csr 类似，但具有行方向。

当我在 MATLAB 中创建 FEM 刚度矩阵时，我使用了 scipy coo 输入的等价物。即创建data、row、col3个数组。当一个 coo 转换为 csr 时，添加了重复元素，巧妙地处理了 FEM 元素的子矩阵重叠。（这种行为在 scipy 和 MATLAB 中是相同的）。

重复添加 lil 矩阵应该可行（如果索引正确），但我预计它会慢得多。

为什么在求解稀疏线性方程组时会出现内存错误？

Why do I get a Memory Error when solving a sparse system of linear equations?

python

matlab

numpy

scipy

sparse-matrix