Есть ли быстрый способ инвертировать матрицу в Matlab?

У меня есть много больших (около 5000 x 5000) матриц, которые мне нужно инвертировать в Matlab. Мне действительно нужен обратный, поэтому я не могу использовать mldivide вместо этого, что намного быстрее для решения Ax = b всего за один b.

Мои матрицы исходят из проблемы, что означает, что у них есть хорошие свойства. Во-первых, их определитель равен 1, поэтому они определенно обратимы. Тем не менее, они не диагонализуемы, или я бы постарался их диагонализировать, инвертировать их, а затем вернуть обратно. Их записи - все действительные числа (фактически рациональные).

Я использую Matlab для получения этих матриц и для этого, что мне нужно делать со своими обратными, поэтому я предпочел бы способ ускорить Matlab. Но если есть другой язык, который я могу использовать, это будет быстрее, тогда, пожалуйста, дайте мне знать. Я не знаю много других языков (немного, но с C и немного, но из Java), поэтому, если это действительно сложно на каком-то другом языке, то я, возможно, не смогу его использовать. Пожалуйста, продолжайте предлагать, однако, в случае.

Ответ 1

Мне действительно нужен обратный, поэтому я не могу использовать mldivide,...

Это не так, потому что вы можете использовать mldivide для получения обратного. Заметим, что A^-1 = A^-1 * I. В MATLAB это эквивалентно

invA = A\speye(size(A));

На моей машине это занимает около 10,5 секунд для матрицы 5000x5000. Обратите внимание, что MATLAB имеет функцию inv для вычисления инверсии матрицы. Хотя это займет примерно столько же времени, оно менее эффективно с точки зрения числовой точности (больше информации в ссылке).

Во-первых, их определитель равен 1, поэтому они определенно обратимы

Вместо det(A)=1 это номер условия вашей матрицы, который определяет, насколько точным или стабильным будет обратный. Обратите внимание, что det(A)=∏_i=1:n λ_i. Поэтому установка λ₁=M, λ_n=1/M и λ_i≠1,n=1 даст вам det(A)=1. Однако, как M → ∞, cond(A) = M² → ∞ и λ_n → 0, что означает, что ваша матрица приближается к сингулярности и будут большие числовые ошибки при вычислении обратного.

Мои матрицы исходят из проблемы, что означает, что у них есть хорошие свойства.

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

Я бы предпочел способ ускорить работу Matlab

MATLAB использует исключение Гаусса для вычисления инверсии общей матрицы (полный ранг, не разреженный, без каких-либо специальных свойств) с использованием mldivide, и это Θ(n³), где n - размер матрицы. Итак, в вашем случае n=5000 и есть операции с плавающей запятой 1.25 x 10¹¹. Таким образом, на разумной машине с 10 Gflops вычислительной мощности вам потребуется не менее 12,5 секунд для вычисления обратного, и нет выхода из этого, если вы не используете "специальные свойства" (если они пригодны для использования )

Ответ 2

Инвертирование произвольной матрицы 5000 x 5000 не является вычислительной легкостью, независимо от того, какой язык вы используете. Я бы рекомендовал посмотреть на приближения. Если ваши матрицы имеют низкий ранг, вы можете попробовать приближение низкого ранга M = USV '

Вот еще несколько идей из math-overflow:

https://mathoverflow.net/search?q=matrix+inversion+approximation

Ответ 3

Предположим сначала, что собственные значения 1. Пусть A - йорданова каноническая форма вашей матрицы. Затем вы можете вычислить A^{-1}, используя только умножение и добавление матрицы на

A^{-1} = I + (I-A) + (I-A)^2 + ... + (I-A)^k

где k < dim(A). Почему это работает? Потому что генерация функций является удивительной. Напомним расширение

(1-x)^{-1} = 1/(1-x) = 1 + x + x^2 + ...

Это означает, что мы можем инвертировать (1-x) с помощью бесконечной суммы. Вы хотите инвертировать матрицу A, поэтому вы хотите взять

A = I - X

Решение для X дает X = I-A. Поэтому подстановкой имеем

A^{-1} = (I - (I-A))^{-1} = 1 + (I-A) + (I-A)^2 + ...

Здесь я только что использовал идентификационную матрицу I вместо числа 1. Теперь у нас есть проблема сходимости для решения, но на самом деле это не проблема. По предположению, что A находится в форме Жордана и имеет все собственные значения, равные 1, мы знаем, что A является верхней треугольной со всеми 1 по диагонали. Поэтому I-A является верхней треугольной со всеми 0 на диагонали. Поэтому все собственные значения I-A 0, поэтому его характеристический многочлен равен x^dim(A), а его минимальный многочлен x^{k+1} для некоторого k < dim(A). Так как матрица удовлетворяет ее минимальному (и характеристическому) многочлену, это означает, что (I-A)^{k+1} = 0. Поэтому приведенная выше серия конечна, причем наибольший ненулевой член равен (I-A)^k. Таким образом, он сходится.

Теперь, в общем случае, поместите вашу матрицу в форму Jordan, так что у вас есть треугольная матрица блока, например:

A 0 0
0 B 0
0 0 C

Где каждый блок имеет одно значение по диагонали. Если это значение A для A, используйте вышеуказанный трюк, чтобы инвертировать 1/a * A, а затем снова умножить A. Поскольку полная матрица является треугольной, обратный будет

A^{-1} 0      0
0      B^{-1} 0
0      0      C^{-1}

Нет ничего особенного в наличии трех блоков, поэтому это работает независимо от того, сколько у вас есть.

Обратите внимание, что этот трюк работает всякий раз, когда у вас есть матрица в форме Джордана. Вычисление инверсии в этом случае будет очень быстрым в Matlab, потому что оно включает только умножение матрицы, и вы даже можете использовать трюки, чтобы ускорить это, поскольку вам нужны только мощности одной матрицы. Это может вам не помочь, если это действительно дорого, чтобы получить матрицу в форме Джорда.