Почему я получаю ошибку памяти при решении разреженной системы линейных уравнений? - PullRequest
Новая
       18

Почему я получаю ошибку памяти при решении разреженной системы линейных уравнений?

1 голос
/ 20 мая 2019

При попытке решить большую разреженную систему линейных уравнений снизу я просто получаю MemoryError:.Как я могу решить эту проблему?

Кроме того, этот код основан на реализации в Matlab, которая должна работать нормально.В оригинальной версии M - это трехмерная матрица, я не знаю, могут ли возникнуть проблемы с памятью из-за моей модификации, преобразовывающей M в 2D 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:

1 Ответ

2 голосов
/ 20 мая 2019

В трассировке должно быть указано, находится ли проблема в 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. Вы должны использовать какое-то стороннее расширение или обходной путь, так как MATLAB sparse ограничен 2d (по крайней мере, так было, когда я использовал его много лет назад для работы в FEM).

Формат scipy.sparse csc аналогичен внутреннему разрежению MATLAB. Фактически это то, что вы получите, если передадите матрицу через save и scipy.io.loadmat. csr аналогично, но с ориентацией строки.

Когда я создавал матрицы жесткости FEM в MATLAB, я использовал эквивалент входных данных scipy coo. То есть создаем 3 массива data, row и col. Когда coo преобразуется в csr, добавляются дублирующие элементы, аккуратно обрабатывая перекрытие подматриц элементов FEM. (это поведение одинаково в scipy и MATLAB).

Повторное добавление матриц lil, как вы делаете, должно работать (если индексация верна), но я ожидаю, что это будет существенно медленнее.

...