Матрица, умноженная на обратную, не дает идентичности

Question

Я пытаюсь найти обратную матрицу данной матрицы, используя функцию np.linalg.inv(). inv дает матрицу, которая выглядит хорошо, но затем при попытке умножить исходную матрицу на обратную, результат не является тождественным, как предполагается определением обратной матрицы.

from numpy.linalg import inv

M = np.random.random((4, 4))
Mi = inv(M)
I = M @ Mi # using matrix multiplication operator
I.astype(int) #as the output looks like 2.77555756e-17

>>> array([[0, 0, 0, 0],
           [0, 0, 0, 0],
           [0, 0, 0, 0],
           [0, 0, 0, 1]])

Что явно не тождество (дает несколько иной ответ при многократном запуске)

Answer 1

Какова ваша цель здесь?

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

Если это так, вы должны сделать следующее:

import numpy as np

matrix1 = np.random.random((4,4))
matrix1_inv = np.linalg.inv(matrix1)

# get what may be identity matrix
ident_pred = matrix1 @ matrix1_inv

# get what is the identity matrix
ident_true = np.identity(matrix1.shape[0])

# check that they are the same
print(np.allclose(ident_pred, ident_true))

Answer 2

Проблема в том, что функция astype не округляется, она просто усекается. Итак, причина, по которой вы не видите единичную матрицу, состоит в том, что другие значения, которые должны быть 1, были где-то около 0.99999. Вы можете использовать это:

import numpy as np
a = np.random.random((4,4))
b = np.linalg.inv(a)
# 0.00001 is the tolerance about which you which to consider values to be == 1
c = np.array(a@b + 0.00001, dtype=int)
print(c)

Если вы хотите просто округлить (высокий допуск == 0,5), используйте вместо этого:

import numpy as np
a = np.random.random((4,4))
b = np.linalg.inv(a)
c = np.array(np.round(a@b), dtype=int)
print(c)

Кроме того, рекомендуется использовать полную функцию np.linalg.inv.

Answer 3

Вопрос в том, матрица близка к np.eye(4) или нет.

Вот как вы должны проверить:

I = M@Mi
EPS = 1e-8
r = np.all(np.abs(I - np.eye(4)) < EPS)

r будет указывать, близки ли две матрицы (I и тождество) к 1e-8.

Answer 4

Когда вы печатаете I, это выглядит так:

array([[ 1.00000000e+00, -5.55111512e-17, -1.66533454e-16, 0.00000000e+00],
       [ 6.38378239e-16,  1.00000000e+00, -5.55111512e-17, 0.00000000e+00],
       [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00, 0.00000000e+00],
       [-5.55111512e-17, -1.11022302e-16, -1.24900090e-16, 1.00000000e+00]])

Однако, записи 1.00 не являются точными. Когда вы печатаете 1 - I, вы можете увидеть это:

array([[-2.22044605e-16,  1.00000000e+00,  1.00000000e+00, 1.00000000e+00],
       [ 1.00000000e+00,  2.22044605e-16,  1.00000000e+00, 1.00000000e+00],
       [ 1.00000000e+00,  1.00000000e+00,  0.00000000e+00, 1.00000000e+00],
       [ 1.00000000e+00,  1.00000000e+00,  1.00000000e+00, 0.00000000e+00]])

Положительные диагональные записи представляют значения в I, которые немного меньше единицы. Когда вы делаете целочисленное усечение (что и делает astype(int)), вы устанавливаете эти элементы на ноль. Вместо этого округлите значения до ближайшего целого числа вместо усечения:

np.around(I).astype(int)

Однако вы не всегда будете иметь целочисленные входные данные, подобные этому, в этом случае округление будет вводить в заблуждение. Numpy предоставляет функцию allclose для сравнения значений в пределах допуска:

np.allclose(I, np.identity(I.shape[0]), rtol=0, atol=1e-15)

Вы также можете выполнить поэлементную проверку, используя isclose:

np.isclose(I, np.identity(I.shape[0]), rtol=0, atol=1e-15)

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

Answer 5

Попробуйте сначала округлить его перед преобразованием в int.

np.around(I).astype(int)

Создание случайной матрицы:

>>> M = np.random.random((4,4))
>>> M
array([[0.51351957, 0.57864882, 0.0489495 , 0.85066216],
       [0.60052988, 0.93844708, 0.74651889, 0.17237584],
       [0.26191596, 0.46451226, 0.46514401, 0.81917544],
       [0.19247662, 0.82801899, 0.83839146, 0.08531949]])

В обратном порядке:

>>> from numpy.linalg import inv
>>> Mi = inv(M)
>>> Mi
array([[-1.3515514 ,  3.53647196,  1.0391335 , -3.64654487],
       [ 2.76188122, -2.23981308, -2.74634579,  3.35680468],
       [-2.44320291,  1.47102487,  2.36135635, -1.28451339],
       [ 0.2533113 , -0.69591469,  1.10498293, -0.00818495]])

Теперь умножение M и Mi должно привести к идентичности.

>>> M @ Mi
array([[ 1.00000000e+00, -4.44089210e-16, -1.11022302e-16, -6.93889390e-18],
       [-4.16333634e-17,  1.00000000e+00, -8.32667268e-17, -8.60856525e-17],
       [ 5.55111512e-17, -2.22044605e-16,  1.00000000e+00, -1.57859836e-16],
       [ 6.24500451e-17, -8.32667268e-17, -2.35922393e-16, 1.00000000e+00]])

Но это явно не идентичность. Но если вы посмотрите внимательно, диагональные значения довольно близки к 1, а все остальные значения - это действительно маленькие числа (почти нули), например -16 или -17 в показателе степени.

Эта ошибка заключается в том, что значения с плавающей точкой никогда не являются точными значениями, в них всегда есть какая-то ошибка. Взгляните на статью 15. Арифметика с плавающей точкой: проблемы и ограничения и Не нарушена ли математика с плавающей точкой? .

Теперь, если мы просто конвертируем его в int, велика вероятность, что он все равно не будет идентичностью. Поскольку значение действительно близко к 1, оно может быть немного меньше 1, что приводит к 0 при приведении к int.

>>> (M @ Mi).astype(int)
array([[1, 0, 0, 0],
       [0, 0, 0, 0],
       [0, 0, 1, 0],
       [0, 0, 0, 0]])

Но если вы округлите его до преобразования в int, вы получите личность.

>>> np.around(M @ Mi).astype(int)
array([[1, 0, 0, 0],
       [0, 1, 0, 0],
       [0, 0, 1, 0],
       [0, 0, 0, 1]])