Обратная матрица

Обра́тная ма́трица — такая матрица [math]\displaystyle{ A^{-1} }[/math], при умножении которой на исходную матрицу [math]\displaystyle{ A }[/math] получается единичная матрица [math]\displaystyle{ E }[/math]:

[math]\displaystyle{ AA^{-1} = A^{-1}A = E. }[/math] Если у матрицы есть обратная матрица, то она является единственной.

Обратную матрицу можно определить как:

[math]\displaystyle{ A^{-1} = \frac{\mbox{adj} A}{|A|}, }[/math]

где [math]\displaystyle{ \mbox{adj} A }[/math] — соответствующая присоединённая матрица,

[math]\displaystyle{ |A| }[/math] — определитель матрицы [math]\displaystyle{ A }[/math].

Из этого определения следует критерий обратимости: матрица обратима тогда и только тогда, когда она невырождена, то есть её определитель не равен нулю. Для неквадратных матриц и вырожденных матриц обратных матриц не существует. Однако возможно обобщить это понятие и ввести псевдообратные матрицы, похожие на обратные по многим свойствам.

Большинство квадратных матриц являются обратимыми, но нахождение обратной матрицы является сложной задачей для больших матриц. Множество всех обратимых матриц размера [math]\displaystyle{ n\times n }[/math] называется полная линейная группа порядка [math]\displaystyle{ n }[/math].

Свойства обратной матрицы

Пусть квадратные матрицы [math]\displaystyle{ A,~ B }[/math] — невырожденные. Тогда:

• Если матрица [math]\displaystyle{ A }[/math] обратимая, то и обратная ей матрица [math]\displaystyle{ A^{-1} }[/math] является обратимой. Матрица обратная матрице [math]\displaystyle{ A^{-1} }[/math] это матрица [math]\displaystyle{ A }[/math], то есть [math]\displaystyle{ \bigl(A^{-1}\bigr)^{-1} = A }[/math]
• [math]\displaystyle{ \det A^{-1} = \frac{1}{\det A} }[/math], где [math]\displaystyle{ \det }[/math] обозначает определитель.
• Если матрицы [math]\displaystyle{ A,~ B }[/math] являются обратимыми матрицами, то и их произведения [math]\displaystyle{ AB }[/math] и [math]\displaystyle{ BA }[/math] являются обратимыми матрицами (произведение матриц некоммутативно). Матрица обратная матрице [math]\displaystyle{ AB }[/math] это матрица [math]\displaystyle{ (AB)^{-1} }[/math], матрица обратная [math]\displaystyle{ BA }[/math] это матрица [math]\displaystyle{ (BA)^{-1} }[/math]. [math]\displaystyle{ (AB)^{-1} = B^{-1}A^{-1}\ }[/math], а [math]\displaystyle{ (BA)^{-1} = A^{-1}B^{-1}.\quad }[/math] Если [math]\displaystyle{ A_1, \ldots, A_m }[/math] обратимые матрицы размера [math]\displaystyle{ n\times n }[/math], то произведение [math]\displaystyle{ A_1 \cdots A_m }[/math] является обратимой матрицей, а обратной ей матрицей является [math]\displaystyle{ A^{-1}_m \cdots A^{-1}_1 }[/math].
• [math]\displaystyle{ (A^T)^{-1} = (A^{-1})^T }[/math], где [math]\displaystyle{ T }[/math] обозначает операцию транспонирования.
• [math]\displaystyle{ (kA)^{-1} = k^{-1}A^{-1} }[/math] для любого коэффициента [math]\displaystyle{ k\not=0 }[/math].
• [math]\displaystyle{ E^{-1} = E }[/math].
• Пусть необходимо решить систему линейных уравнений [math]\displaystyle{ Ax=b }[/math], где [math]\displaystyle{ x }[/math] — искомый вектор, [math]\displaystyle{ b }[/math] — ненулевой вектор. Если [math]\displaystyle{ A^{-1} }[/math] существует, то [math]\displaystyle{ x=A^{-1} b }[/math]. В противном случае либо размерность пространства решений больше нуля, либо их нет вовсе.
• [math]\displaystyle{ \left({A} - {B}\right)^{-1} = {A}^{-1} + {A}^{-1}{B}\left({A} - {B}\right)^{-1}, }[/math]
• [math]\displaystyle{ \left({A} + {B}\right)^{-1} = \left({I} + {A}^{-1}{B}\right)^{-1}{A}^{-1}, }[/math]
• [math]\displaystyle{ \left({A} - {B}\right)^{-1} = \sum_{k=0}^{\infty} \left({A}^{-1}{B}\right)^k{A}^{-1}. }[/math]

Способы нахождения обратной матрицы

Если матрица обратима, то для нахождения обратной матрицы можно воспользоваться одним из следующих способов:

Точные (прямые) методы

Метод Жордана—Гаусса

Возьмём две матрицы: саму [math]\displaystyle{ A }[/math] и единичную матрицу [math]\displaystyle{ E }[/math]. Приведём матрицу [math]\displaystyle{ A }[/math] к единичной методом Гаусса—Жордана, применяя преобразования по строкам (можно также применять преобразования и по столбцам). После применения каждой операции к первой матрице применим ту же операцию ко второй. Когда приведение первой матрицы к единичному виду будет завершено, вторая матрица окажется равной [math]\displaystyle{ A^{-1} }[/math].

При использовании метода Гаусса первая матрица будет умножаться слева на одну из элементарных матриц [math]\displaystyle{ \Lambda_i }[/math] (трансвекцию или диагональную матрицу с единицами на главной диагонали, кроме одной позиции):

[math]\displaystyle{ \Lambda_1 \cdot \dots \cdot \Lambda_n \cdot A = \Lambda A = E \Rightarrow \Lambda = A^{-1}. }[/math]

[math]\displaystyle{ \Lambda_m = \begin{bmatrix} 1 & \dots & 0 & -a_{1m} /a_{mm} & 0 &\dots & 0 \\ & & &\dots & & &\\ 0 & \dots & 1 & -a_{m-1m} /a_{mm} & 0 &\dots &0 \\ 0 & \dots & 0 & 1/a_{mm} & 0 &\dots & 0 \\ 0 & \dots & 0 & -a_{m+1m} /a_{mm} & 1 &\dots &0 \\ & & &\dots & & &\\ 0 & \dots & 0 & -a_{nm}/a_{mm} & 0 &\dots & 1 \end{bmatrix}. }[/math]

Вторая матрица после применения всех операций станет равна [math]\displaystyle{ \Lambda }[/math], то есть будет искомой. Сложность алгоритма — [math]\displaystyle{ O(n^3) }[/math].

С помощью матрицы алгебраических дополнений

Матрица, обратная матрице [math]\displaystyle{ A }[/math], представима в виде:

[math]\displaystyle{ {A}^{-1}= {{\mbox{adj} (A)}\over{\det(A)}}, }[/math]

где [math]\displaystyle{ \mbox{adj}(A) }[/math] — присоединенная матрица (матрица, составленная из алгебраических дополнений для соответствующих элементов транспонированной матрицы).

Сложность алгоритма зависит от сложности [math]\displaystyle{ O_{\det} }[/math] алгоритма расчета определителя и равна [math]\displaystyle{ O(n^2) \cdot O_{\det} }[/math].

Использование LU- или LUP-разложения

Матричное уравнение [math]\displaystyle{ AX=I_n }[/math] для обратной матрицы [math]\displaystyle{ X }[/math] можно рассматривать как совокупность [math]\displaystyle{ n }[/math] систем вида [math]\displaystyle{ Ax=b }[/math]. Обозначим [math]\displaystyle{ i }[/math]-й столбец матрицы [math]\displaystyle{ X }[/math] через [math]\displaystyle{ X_i }[/math]; тогда [math]\displaystyle{ AX_i=e_i }[/math], [math]\displaystyle{ i=1,\ldots,n }[/math], поскольку [math]\displaystyle{ i }[/math]-м столбцом матрицы [math]\displaystyle{ I_n }[/math] является единичный вектор [math]\displaystyle{ e_i }[/math]. Иными словами, нахождение обратной матрицы сводится к решению [math]\displaystyle{ n }[/math] уравнений с одной матрицей и разными правыми частями. Решение этих уравнений может быть упрощено с помощью LU- или LUP-разложения матрицы [math]\displaystyle{ A }[/math]. После выполнения LUP-разложения за время [math]\displaystyle{ O(n^3) }[/math] на решение каждого из [math]\displaystyle{ n }[/math] уравнений нужно время [math]\displaystyle{ O(n^2) }[/math], так что и этот алгоритм требует времени [math]\displaystyle{ O(n^3) }[/math].

Матрицу, обратную к заданной невырожденной матрице [math]\displaystyle{ A }[/math], можно также вычислить непосредственно с помощью матриц, полученных в результате разложения.

Результатом LUP-разложения матрицы [math]\displaystyle{ A }[/math] является равенство [math]\displaystyle{ PA=LU }[/math]. Пусть [math]\displaystyle{ PA=B }[/math], [math]\displaystyle{ B^{-1}=D }[/math]. Тогда из свойств обратной матрицы можно записать: [math]\displaystyle{ D=U^{-1}L^{-1} }[/math]. Если умножить это равенство на [math]\displaystyle{ U }[/math] и [math]\displaystyle{ L }[/math] то можно получить два равенства вида [math]\displaystyle{ UD=L^{-1} }[/math] и [math]\displaystyle{ DL=U^{-1} }[/math]. Первое из этих равенств представляет собой систему из [math]\displaystyle{ n^2 }[/math] линейных уравнений, для [math]\displaystyle{ n(n+1)/2 }[/math] из которых известны правые части (из свойств треугольных матриц). Второе также представляет систему из [math]\displaystyle{ n^2 }[/math] линейных уравнений, для [math]\displaystyle{ n(n-1)/2 }[/math] из которых известны правые части (также из свойств треугольных матриц). Вместе они представляют собой систему из [math]\displaystyle{ n^2 }[/math] равенств. С их помощью можно рекуррентно определить все [math]\displaystyle{ n^2 }[/math] элементов матрицы [math]\displaystyle{ D }[/math]. Тогда из равенства [math]\displaystyle{ (PA)^{-1} = A^{-1} P^{-1} = B^{-1} = D }[/math] получаем равенство [math]\displaystyle{ A^{-1} = DP }[/math].

В случае использования LU-разложения ([math]\displaystyle{ A=LU }[/math]) не требуется перестановки столбцов матрицы [math]\displaystyle{ D }[/math], но решение может разойтись даже если матрица [math]\displaystyle{ A }[/math] невырождена.

Сложность обоих алгоритмов — [math]\displaystyle{ O(n^3) }[/math].

Итерационные методы

Матрицу [math]\displaystyle{ A^{-1} }[/math] можно вычислить с произвольной точностью в результате выполнения следующего итерационного процесса, называющегося методом Шульца и являющегося обобщением классического метода Ньютона:

[math]\displaystyle{ X_{k+1} = 2X_k - X_k A X_k. }[/math]

Последовательность матриц [math]\displaystyle{ X_k }[/math] сходится к [math]\displaystyle{ A^{-1} }[/math] при [math]\displaystyle{ k \to \infty }[/math]. Существует также так называемый обобщённый метод Шульца, который описывается следующими рекуррентными соотношениями:

[math]\displaystyle{ \begin{cases} \Psi_k=E-AX_k, \\ X_{k+1}=X_k \sum\limits_{i=0}^n \Psi^i_k. \end{cases} }[/math]

Выбор начального приближения

Проблема выбора начального приближения [math]\displaystyle{ X_0 }[/math] в рассматриваемых здесь процессах итерационного обращения матриц не позволяет относиться к ним как к самостоятельным универсальным методам, конкурирующими с прямыми методами обращения, основанными, например, на [math]\displaystyle{ LU }[/math]-разложении матриц. Имеются некоторые рекомендации по выбору [math]\displaystyle{ X_0 }[/math], обеспечивающие выполнение условия [math]\displaystyle{ \rho(\Psi_0)\lt 1 }[/math] (спектральный радиус матрицы меньше единицы), являющегося необходимым и достаточным для сходимости итерационного процесса. Однако при этом, во-первых, требуется знать оценку сверху спектра обращаемой матрицы [math]\displaystyle{ A }[/math] либо матрицы [math]\displaystyle{ AA^T }[/math] (а именно, если [math]\displaystyle{ A }[/math] — симметричная положительно определённая матрица и [math]\displaystyle{ \rho(A)\leqslant\beta }[/math], то можно взять [math]\displaystyle{ X_0={\alpha}E }[/math], где [math]\displaystyle{ \alpha\in\left(0,2/\beta\right) }[/math]; если же [math]\displaystyle{ A }[/math] — произвольная невырожденная матрица и [math]\displaystyle{ \rho(AA^T)\leqslant\beta }[/math], то полагают [math]\displaystyle{ X_0={\alpha}A^T }[/math], где также [math]\displaystyle{ \alpha\in\left(0,2/\beta\right) }[/math]; можно, конечно, упростить ситуацию и, воспользовавшись тем, что [math]\displaystyle{ \rho(AA^T)\leqslant\mathcal{k}AA^T\mathcal{k} }[/math], положить [math]\displaystyle{ X_0=A^T/\|AA^T\| }[/math]). Во-вторых, при таком задании начальной матрицы нет гарантии, что [math]\displaystyle{ \|\Psi_0\| }[/math] будет малой (возможно, даже окажется [math]\displaystyle{ \|\Psi_0\|\gt 1 }[/math]), и высокий порядок скорости сходимости обнаружится далеко не сразу.

Для метода Ньютона в качестве начального приближения можно выбрать [math]\displaystyle{ X_0 = A^H/ \left( ||A||_1 ||A||_\infty \right) }[/math], где верхний индекс [math]\displaystyle{ H }[/math] обозначает эрмитово сопряжение, [math]\displaystyle{ || \cdot ||_1 }[/math] и [math]\displaystyle{ || \cdot ||_\infty }[/math] — соответствующие матричные нормы. Такое [math]\displaystyle{ X_0 }[/math] вычисляется всего за [math]\displaystyle{ O(n^2) }[/math] операций, где [math]\displaystyle{ n }[/math] — порядок матрицы, и обеспечивает сходимость алгоритма.

Примеры

Матрица 2 × 2

[math]\displaystyle{ \mathbf{A}^{-1} = \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix}^{-1} = \frac{1}{\det\mathbf{A}} \begin{bmatrix} d & -b \\ -c & a \\ \end{bmatrix} = \frac{1}{ad - bc} \begin{bmatrix} d & -b \\ -c & a \\ \end{bmatrix}. }[/math]

Обращение матрицы 2 × 2 возможно только при условии, что [math]\displaystyle{ ad - bc = \det A \neq 0 }[/math].

Ссылки

• Реализация с полным выбором ведущего элемента на C++