SciPy разреженные инвертировать или spsolve приводят к UMFPACK_ERROR_OUT_OF_MEMORY

Question

Я пытаюсь перевернуть большую (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) 

a d Я снова столкнулся с той же ошибкой

Я использовал этот код на машине с 16 ГБ ОЗУ, а затем переместил ее на сервер с 32 ГБ ОЗУ, но все же безрезультатно.

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

Answer 1

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

Начнем с основ: никогда инвертировать разреженную матрицу, это совершенно бессмысленно. См. Это discussion на центральном центральном устройстве 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)) равна не a один вектор (np.ones(N)), как вы сомневаетесь в ложных предположениях.

Дело в том, что матричные инверсии редко используются в числовой линейной алгебре: решение Ax = b не вычисляется как inv (A) * b, а специализированным алгоритмом.

Переходя к вашей конкретной проблеме, для большой разреженной системы уравнений нет black-box решателей. Вы можете выбрать правильный класс решателей, только если у вас есть хорошее представление о структуре и свойствах вашей матрицы. Свойства ваших матриц в свою очередь являются следствием проблемы, которую вы пытаетесь решить. Например.когда вы дискретизируете FEM системой эллиптического PDE, вы получаете симметричную положительную разреженную систему алгебраических уравнений. Как только вы узнаете о свойствах своей проблемы, вы можете выбрать правильную стратегию решения.

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

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

Answer 2

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

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

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