ошибки округления инверсии матрицы numpy

Когда вы вычисляете обратное, ваши массивы преобразуются в float64, чей машинный эпсилон равен 1e-15. Эпсилон — это относительный шаг квантования числа с плавающей запятой.

Если есть сомнения, мы можем запросить у numpy информацию о типе данных с плавающей запятой, используя функцию finfo. В таком случае

np.finfo('float64')
finfo(resolution=1e-15, 
      min=-1.7976931348623157e+308, max=1.7976931348623157e+308, 
      dtype=float64)

Итак, технически ваше значение меньше eps является очень точным представлением 0 для типа float64.

Если вас беспокоит только представление, вы можете указать numpy не печатать маленькие числа с плавающей запятой (1 eps или меньше от 0) с помощью:

np.set_printoptions(suppress=True)

После этого ваш оператор печати возвращает:

[[ 0.  -2.5]
 [-2.  -6.5]]

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

• Почему числа с плавающей запятой неточны?

или в сети:

• Проблемы точности с плавающей запятой
• Что каждый программист должен знать об арифметике с плавающей запятой

Это не полный ответ, но он может указать вам правильное направление. Что вам действительно нужно, так это массивы numpy, которые используют десятичные числа для математики. Вы можете разумно подумать, чтобы попробовать:

import numpy as np
from decimal import Decimal
A = np.array([[2,2],[4,-1]],np.int)
for i, a in np.ndenumerate(A):
    A[i] = Decimal(a)
    print type(A[i])

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

Одной из возможностей было бы установить тип данных, таким образом:

def decimal_array(arr):
    X = np.array(arr, dtype = Decimal)
    for i, x in np.ndenumerate(X): X[i] = Decimal(x)
    return X

A = decimal_array([[2,2],[4,-1]])
B = decimal_array([[1,3],[-1,-1]])

A = A.transpose()
B = B.transpose()
Binv = np.linalg.inv(B) #calculate the inverse

Но теперь, если вы

print Binv.dtype

вы увидите, что инверсия преобразовала его обратно в плавающий. Причина в том, что linalg.inv (как и многие другие функции) ищет "common_type" B, который является скаляром, к которому, по ее мнению, он может принудить элементы вашего массива.

Хотя, возможно, это не безнадежно. Я посмотрел, сможете ли вы решить эту проблему, создав собственный тип dtype, но оказалось, что скаляры (целые числа, числа с плавающей запятой и т. д.) вообще не являются типами dtype. Вместо этого вы, вероятно, захотите зарегистрировать новый скаляр — это Decimal — как указано в файле статья о скалярах. Вы увидите ссылку на Numpy C-API (не бойтесь). Найдите на странице слова «регистрация» и «скаляр», чтобы начать.