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

Question

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

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

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

Answer 1

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

Это не так, потому что вы все еще можете использовать mldivide, чтобы получить обратное. Обратите внимание, что A<sup>-1</sup> = A<sup>-1</sup> * I. В MATLAB это эквивалентно

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

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

---

Прежде всего, их определитель равен 1, поэтому они определенно обратимы

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

---

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

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

---

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

MATLAB использует исключение Гаусса для вычисления обратной матрицы общего порядка (полный ранг, не разреженный, без каких-либо специальных свойств), используя mldivide, и это Θ(n<sup>3</sup>), где n - размер матрицы. Итак, в вашем случае, n=5000 и есть 1.25 x 10<sup>11</sup> операций с плавающей запятой. Таким образом, на разумной машине с вычислительной мощностью около 10 Гфлопс вам понадобится не менее 12,5 секунд для вычисления обратного, и выхода из этого нет, если вы не используете «особые свойства» (если они могут использоваться )

Answer 2

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

Вот еще несколько идей из математического переполнения:

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

Answer 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 s на диагонали.Следовательно, I-A является верхним треугольником со всеми 0 s по диагонали.Следовательно, все собственные значения I-A равны 0, поэтому его характеристический многочлен равен x^dim(A), а минимальный многочлен равен x^{k+1} для некоторого k < dim(A).Поскольку матрица удовлетворяет своему минимальному (и характеристическому) полиному, это означает, что (I-A)^{k+1} = 0.Следовательно, вышеприведенный ряд равен конечному , причем самый большой ненулевой член равен (I-A)^k.Так что она сходится.

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

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