Как использовать scipy.linalg.solve_banded, чтобы решить Ax = b, когда A является полосатой матрицей? - PullRequest
Новая
       11

Как использовать scipy.linalg.solve_banded, чтобы решить Ax = b, когда A является полосатой матрицей?

0 голосов
/ 01 октября 2018

Для решения системы линейных уравнений 

Ax=b

, в которой A - полосатая матрица, подобная

image

, имеющей 5 ненулевых диагоналей.

Но в отличие от полосовой матрицы, A имеет три ненулевые диагонали, смещения которых равны 0, -1 и 1, и две ненулевые диагонали со смещениями -m и m.

Iпопытался решить его напрямую с помощью 

diagonals = [Ap, As, An, Aw, Ae]
A = diags(diagonals, [0, -1, 1, -m, m]).toarray()
x = linalg.solve(A, b)

Этот метод создал целое число A. Но запасные части A, так что этот метод потратил слишком много памяти, чтобы сохранить нулевые элементы.

Поэтому я попытался использовать solve_banded 

A = np.zeros((2 * m + 1, len(initial)))
A[0] = Ae
A[m - 1] = An
A[m] = Ap
A[m + 1] = As
A[2 * m] = Aw
x = linalg.solve_banded((m, m), A, b)

Этот метод стоит меньше памяти, чем предыдущие, но он все равно потратил некоторое количество на (2m-4) нулевых векторов.Есть ли более разумные методы, которые используют только пять ненулевых векторов?

1 Ответ

0 голосов
/ 26 апреля 2019

Эй, я могу частично ответить на этот вопрос.Чтобы уменьшить объем памяти, вы можете сохранить форму разреженной матрицы, не превращая ее в большую матрицу, выполнив команду to.array(): 

A = diags(diagonals, [0, -1, 1, -m, m])

Теперь разреженная матрица имеет собственный метод spsolve и, следовательно, scipy.sparse.linalg.spsolve(A,b) будет работать.

Чтобы еще больше повысить производительность, вы можете sparse.linalg.splu(A).solve(b), который преобразует A с помощью LU-декомпозиции, а затем объект SuperLU снова имеет метод solve.(Как я уже проверял, этот метод немного быстрее, чем непосредственно spsolve, который также использует LU-подобную декомпозицию. См., Например, здесь .)

Единственная проблема здесь - в серединеLU-декомпозиции использование памяти также будет огромным из-за декомпозиции.

Мне также было интересно, есть ли способ синтезировать метод solve_banded с разреженной матрицей, так как не существует встроенного метода.

...