Сложность обращения матрицы в NumPy

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

Мой текущий выбор - numpy.linalg.inv . Из моих чисел я вижу, что он масштабируется как где n - количество строк, поэтому метод, похоже, исключает Гауссу. $O(n^3)$

Согласно Википедии , существуют более быстрые алгоритмы. Кто-нибудь знает, есть ли библиотека, которая реализует это?

Интересно, а почему бы не использовать эти более быстрые алгоритмы?

numpy dense-matrix inverse physicsGuy
источник

Вы должны выполнить свои матрицы раньше. Посмотри на Сципи. Разреженный за вашу помощь. Он содержит много инструментов, которые вам нужны.

Tobal

@Tobal не уверен, что я следую ... как бы вы "выполнить" матрицу? и как именно scipy.sparseпоможет?

GoHokies

@GoHokies Scipy является дополнением к NumPy. Плотные / разреженные матрицы должны быть реализованы задолго до того, как вы выполните некоторые вычисления, это улучшает ваши вычисления. Пожалуйста, прочтите этот docs.scipy.org/doc/scipy/reference/sparse.html, он лучше, чем я, объясняет мой плохой английский.

Tobal

@Tobal Вопрос, в частности, относится к плотным матрицам, поэтому я не понимаю, как scipy.sparseздесь уместно?

Кристиан Клэйсон

@Tobal - я думаю, что до сих пор не понимаю. Что именно вы имеете в виду под «преформом ваших матриц» и «матрицы должны быть реализованы задолго до того, как вы выполните некоторые вычисления»? Что касается вашего последнего комментария, вы наверняка согласитесь, что методы, которые можно использовать для разреженных и плотных матриц, очень разные.

Вольфганг Бангерт

Ответы:

(Это слишком долго для комментариев ...)

$n$

$\mathcal{O}$ $n$ $C_1 n^3$ $C_2 n^{2.x}$ $n$
Сложность предполагает, что каждая (арифметическая) операция занимает одно и то же время, но на практике это далеко не так: умножение набора чисел на одно и то же происходит намного быстрее, чем умножение одинакового количества разных чисел. Это связано с тем, что основным узким местом в современных вычислениях является получение данных в кэш, а не фактические арифметические операции с этими данными. Таким образом, алгоритм, который можно переставить так, чтобы он имел первую ситуацию (называемую кэшированием ), будет намного быстрее, чем алгоритм, в котором это невозможно. (Это относится к алгоритму Штрассена, например.)

Кроме того, численная стабильность, по крайней мере, так же важна, как и производительность; и здесь, опять же, стандартный подход обычно побеждает.

$\mathcal{O}(n^3)$ $\mathcal{O}(n^{2.x})$

A^{- 1} b

$A^{-1}b$

A x = b

$Ax=b$ numpy.linalg.solve

x

$x$

A

$A$

A^{- 1}

$A^{-1}$

A

$A$

Кристиан Клэйсон
источник

Отличный ответ, спасибо, сэр, в частности, за то, что вы указали на дьявола в деталях (константы в больших цифрах O), который делает большую разницу между теоретической скоростью и практической скоростью.

gaborous

Я думаю, что «обратная редко нужна» часть должна быть подчеркнута больше. Если цель состоит в том, чтобы решить систему дифференциальных уравнений, маловероятно, что требуется полное обратное.

Джаред Гогуэн

@o_o Ну, это был мой первый оригинальный комментарий (который я удалил после объединения всех в один ответ). Но я подумал, что в интересах сайта (и более поздних читателей) ответ должен ответить на реальный вопрос в вопросе (который является разумным и тематическим), даже если за этим стоит проблема XY. Кроме того , я не хочу показаться слишком увещевать ...

Christian Клейсон

n

$n$

A

$A$

Вы, вероятно, должны заметить, что глубоко внутри скрытого исходного кода (см. Https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src ) процедура inv пытается вызвать функцию dgetrf из вашего системного пакета LAPACK, который затем выполняет LU-декомпозицию вашей исходной матрицы. Это морально эквивалентно исключению Гаусса, но может быть настроено на несколько меньшую сложность с помощью более быстрых алгоритмов умножения матриц в высокопроизводительной BLAS.

Если вы пойдете по этому пути, вы должны быть предупреждены, что заставить всю цепочку библиотек использовать новую библиотеку, а не системную, которая шла с вашим дистрибутивом, довольно сложно. Одной из альтернатив в современных компьютерных системах является рассмотрение параллельных методов с использованием таких пакетов, как scaLAPACK или (в мире python) petsc4py. Однако их, как правило, удобнее использовать в качестве итерационных решателей для систем линейной алгебры, чем применять к прямым методам, а PETSc в конкретных целевых разреженных системах больше, чем плотных.

origimbo
источник