Скупой редкий инвертированный или spsolve привести к 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 ГБ ОЗУ, а затем перенес его на сервер с 32 ГБ ОЗУ, но все еще безрезультатно.

кто-нибудь сталкивался с этим раньше?

Источник

2 ответа

Решение

Прежде всего позвольте мне сказать, что этот вопрос лучше задать на http://scicomp.stackexchange.com/ где есть большое сообщество экспертов в области вычислительной науки и числовой линейной алгебры.

Давайте начнем с основ: никогда не инвертируйте разреженную матрицу, это совершенно бессмысленно. Смотрите эту дискуссию о MATLAB Central и особенно этот комментарий Тима Дэвиса.

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

В вашем случае вы пытаетесь использовать общий прямой решатель, не переупорядочивая уравнения. Хорошо известно, что это приведет к заполнению, которое разрушит разреженность iBs матрица на первом этапе spsolve функция (которая должна быть факторизацией). Обратите внимание, что для матрицы с двойной точностью 150000 x 150000 требуется около 167 ГБ памяти. Существует много способов переупорядочения уравнений, чтобы уменьшить заполнение во время факторизации, но вы не предоставляете достаточно информации, чтобы дать вам разумную подсказку.

Извините, но вы должны рассмотреть вопрос о переформулировании своего вопроса на http://scicomp.stackexchange.com/ четко указав, какую проблему вы пытаетесь решить, чтобы дать представление о структуре и свойствах матрицы.

Источник

Разреженный массив только помещает ненулевые записи вашей матрицы в память. Теперь предположим, что вы делаете инверсию. Это означает, что почти все элементы матрицы становятся ненулевыми. Разреженные матрицы оптимизированы по памяти.

Есть несколько операций, которые вы можете применить к разреженным матрицам, не теряя свойство "spare":

  • Кроме того, простое добавление константы может сохранить разреженную матрицу разреженной.
Источник
Другие вопросы по тегам