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

Question

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

此外，此代码基于 Matlab 中的实现，应该可以正常运行。原版中，M是一个三维矩阵，不知道是否是因为我修改将M转换为二维而出现内存问题scipy.sparse.lil_matrix(而不是coo)我可以迭代地填充M。

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 稀疏矩阵。您必须使用某种第 3 方扩展或解决方法，因为 MATLAB 稀疏仅限于 2d(至少是几年前我在 FEM 工作中使用它时)。

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

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

像您一样重复添加 lil 矩阵应该可以工作(如果索引正确)，但我预计它会慢得多。