Scipy разреженный инверс или 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 ГБ ОЗУ, но все же безрезультатно.

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

Ответ 1

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

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

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

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

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

Ответ 2

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

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

Добавление, просто добавление константы может уменьшить разреженную матрицу.