Обратная матрица методом Гаусса

Для получения матрицы $A^{-1}$ будем исходить из того, что она является решением уравнения $AX = E$, где $E$ - единичная матрица. Представим искомую матрицу $X=(x_{ij})^n_{i, j=1}$ как набор векторов-столбцов img а единичную матрицу $E$ как набор единичных векторов img Матричное уравнение $AX = E$ в соответствии с правилами умножения матриц возможно заменить не связанной между собой системой уравнений img Каждое из этих уравнений может быть решено методом Гаусса. Заметим то обстоятельство, что все СЛАУ имеют одну и ту же матрицу коэффициентов, поэтому

Суть метода Гаусса-Жордана заключается в том, что если с единичной матрицей $E$ провести элементарные преобразования, которыми невырожденная квадратная матрица $A$ приводится к $E$, то получится обратная матрица $A^{-1}$. Для отображения всех элементарных преобразований, совершающихся над матрицей $A$, на единичную матрицу $E$, удобно "склеить" две матрицы в одну. Если матрицы $A$ и $E$ имели размер $n\times n$, то склеенная матрица будет иметь размер $n\times 2*n$. Таким образом, все преобразования над матрицей $A$ (перестановка строк, умножение строки на константу) будут автоматически совершаться и над единичной матрицей). Обозначим склеенную матрицу как $M = (m_{ij})$ Приведение матрицы $A$ к единичной совершается в два этапа: на первом этапе будем двигаться "сверху-вниз", получая в итоге треугольную матрицу, с единицами на главной диагонали и нулями ниже. На втором этапе будем двигаться "снизу-вверх", преобразуя элементы, лежащие выше главной диагонали к нулю, тем самым добиваясь единичной матрицы в результате.

Прямой ход

Пусть матрица $A$ имеет размер $n\times n$. Тогда на первом этапе будет совершено $n$ шагов. На каждом шаге необходимо совершить три действия.

1. Поменять местами строки $i$ и $i + k$, в случае если $m_{ii} = 0$. Этот шаг необходим для исключения ситуации нахождения на главной диагонали 0.
2. Преобразовать элемент, стоящий на главной диагонали к 1. Для этого необохдимо домножить всю строку на $\frac{1}{m_{ii}}$.
3. Обнулить все нижележащие элементы столбца. Для этого вычтем из каждой нижележащей строки c индексом $i + k, k=i+1..n$ результат умножения элемента $m_{i+k, i}$ c текущей строкой. Напомним, что в результате шага 2 первый элемент $m_{kk}$ 1, поэтому первым элементом вычитаемой строки будет являться значение первого элемента строки $i+k$, что и даст 0 в элементе $m_{i+k, i}$

Обратный ход

Теперь необходимо преобразовать матрицу $A$ так, чтобы все элементы n-ого столбца выше $m_{nn}$ стали равны нулю. Для этого прибавляем к n-1 строке соответсвующие элементы n-ой строки, умноженные на $-m_{n-1, n}$. К n-2 строке прибавляем соответсвующие элементы (n-1)-ой строки, умноженные на $-m_{n-2, n}$ и т.д. Аналогичные действия необходимо совершить над оставшимися $n-3, n-4, ..., n - (n-2)$ строками.

Реализация

Приведем реализацию функции inverse, принимающей в качестве аргумента исходную матрицу коэффициентов, и возвращающую матрицу, обратную к исходной. Склеим матрицу коэффициентов с единичной матрицей с помощью функции hstack

m = np.hstack((matrix_origin, 
                np.matrix(np.diag([1.0 for i in range(matrix_origin.shape[0])]))))

Прямой ход

# forward trace
    for k in range(n):
        # 1) Swap k-row with one of the underlying if m[k, k] = 0
        swap_row = pick_nonzero_row(m, k)
        if swap_row != k:
            m[k, :], m[swap_row, :] = m[swap_row, :], np.copy(m[k, :])
        # 2) Make diagonal element equals to 1
        if m[k, k] != 1:
            m[k, :] *= 1 / m[k, k]
        # 3) Make all underlying elements in column equal to zero
        for row in range(k + 1, n):
            m[row, :] -= m[k, :] * m[row, k]

Для проверки первгого условия реализуем вспомогательную функцию pick_non_zero_row, возвращающую индекс первой строки, в которой элемент интересующего нас столбца не равен нулю

def pick_nonzero_row(m, k):
    while k < m.shape[0] and not m[k, k]:
        k += 1
    return k

Обратный ход

# backward trace
    for k in range(n - 1, 0, -1):
        for row in range(k - 1, -1, -1):
            if m[row, k]:
                # 1) Make all overlying elements equal to zero in the former identity matrix
                m[row, :] -= m[k, :] * m[row, k]

Возвратим преобразованную единичную матрицу, т.е. вторую часть "склееного" массива, воспользовавшись функцией hsplit

return np.hsplit(m, n // 2)[1]