Построить обратную матрицу методом Гаусса

Question

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

Спасибо: ранее здесь подсказали с первой частью - массив дает треугольную матрицу. Однако пытаясь выполнить обратный ход, ловлю ошибку, свидетельствующую о проблемах с размерностью: shapes (1,6) and (1,6) not aligned: 6 (dim 1) != 1 (dim 0)

import numpy as np
def inverse_matrix(matrix_origin):
    """
    Функция получает на вход матрицу, затем добавляет к ней единичную матрицу, 
    проводит элементарные преобразования по строкам с первоначальной, добиваясь получения слева единичной матрицы.
    В этом случае справа окажется матрица, которая является обратной к заданнй первоначально 
    """
    # Склеиваем 2 матрицы: слева - первоначальная, справа - единичная
    m = np.hstack((matrix_origin, 
                np.matrix(np.diag([1.0 for i in range(matrix_origin.shape[0])]))))
#     n = matrix_origin.shape[1]
    n = m.shape[1]
    
    for nrow, row in enumerate(m):
        # nrow равен номеру строки
        # row содержит саму строку матрицы
        divider = row[nrow] # диагональный элемент
        # делим на диагональный элемент:
        row /= divider
        # теперь вычитаем приведённую строку из всех нижележащих строк:
        for lower_row in m[nrow+1:]:
            factor = lower_row[nrow] # элемент строки в колонке nrow
            lower_row -= factor*row # вычитаем, чтобы получить ноль в колонке nrow

# обратный ход:
    for k in range(n - 1, 0, -1):
        for row_ in range(k - 1, -1, -1):
            if m[row_, k]:
                # 1) Все элементы выше главной диагонали делаем равными нулю
                m[row_, :] -= m[k, :] * m[row_, k]
    return np.hsplit(m, n // 2)[1]
    
            
matrix = np.array([[3.8, 6.7, -1.2], 
                   [6.4, 1.3, -2.7], 
                   [2.4, -4.5, 3.5]])

inverse_matrix(np.copy(matrix))

Answer 1

У вас несколько ошибок в коде.

Вы конструируете правую часть как np.matrix, поэтому после np.hstack вы получаете тоже объект типа np.matrix. Это очень специальный вид массивов. В частности, итератор m возвращает не одномерные строки чисел, а двумерные массивы формы (1,6). В результате row[nrow] оказывается не числом, а одномерным массивом.

В numpy есть специальная функция для построения единичной матрицы: np.eye(n) строит единичную матрицу размером n x n. Поэтому вам лучше конструировать m как m = np.hstack((matrix_origin, np.eye(len(matrix_origin)))
Вы инициализировали n как число колонок в m: m = m.shape[1]. Но тогда вот эта строчка неверна: for row_ in range(k - 1, -1, -1) - у вас нет строк с номерами 5,4,3. Если n определяет число строк, то нужно присваивать вот так: n = m.shape[0]
Но в таком случае сломается ваш код выделения правой части np.hsplit(m, n // 2)[1]. Я предлагаю использовать индексирование вместо hsplit: m[:, n:].copy() есть двумерный массив, каждая строка которого есть строка из m начиная с элемента с номером n, то есть элементы №№ 3,4 и 5. Как раз правая часть матрицы m. Метод copy() вызывается для того, чтобы не держать указатель на m, в противном случае m будет висеть в памяти до тех пор, пока "жив" указатель на результат функции inverse_matrix

После исправления этих ошибок получится что-то вроде

def inverse_matrix(matrix_origin):
    """
    Функция получает на вход матрицу, затем добавляет к ней единичную матрицу, 
    проводит элементарные преобразования по строкам с первоначальной, добиваясь получения слева единичной матрицы.
    В этом случае справа окажется матрица, которая является обратной к заданнй первоначально 
    """
    # Склеиваем 2 матрицы: слева - первоначальная, справа - единичная
    n = matrix_origin.shape[0]
    m = np.hstack((matrix_origin, np.eye(n)))
    
    for nrow, row in enumerate(m):
        # nrow равен номеру строки
        # row содержит саму строку матрицы
        divider = row[nrow] # диагональный элемент
        # делим на диагональный элемент:
        row /= divider
        # теперь вычитаем приведённую строку из всех нижележащих строк:
        for lower_row in m[nrow+1:]:
            factor = lower_row[nrow] # элемент строки в колонке nrow
            lower_row -= factor*row # вычитаем, чтобы получить ноль в колонке nrow
    # обратный ход:
    for k in range(n - 1, 0, -1):
        for row_ in range(k - 1, -1, -1):
            if m[row_, k]:
                # 1) Все элементы выше главной диагонали делаем равными нулю
                m[row_, :] -= m[k, :] * m[row_, k]
    return m[:,n:].copy()

Результат инвертирования вашей матрицы

array([[ 0.04128819,  0.09805945,  0.08980182],
       [ 0.15689513, -0.08790039, -0.01401625],
       [ 0.1734104 , -0.18025555,  0.206115  ]])

Я добавил jupyter notebook со своим вариантом, он быстрее на 10-15% за счёт выбора операций индексирования.

БЛОГ НА HUSL

Построить обратную матрицу методом Гаусса

Ответы (1 шт):