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

Question

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

此外，此代码基于 Matlab 中的实现，应该可以正常运行。在原始版本中，M 是一个三维矩阵，我不知道是否会因为我的修改将M 转换为二维scipy.sparse.lil_matrix（而不是coo）而出现内存问题。我可以反复填写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)

根据这些非零值的结构，@ 乘积的非零值可能比 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 稀疏矩阵。您必须使用某种 3rd 方扩展或解决方法，因为 MATLAB sparse 仅限于 2d（至少在几年前我将它用于 FEM 工作时）。

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

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

按照您的做法重复添加 lil 矩阵应该可以工作（如果索引正确），但我预计它会慢很多。