Что означают кратные машинного эпсилона? - PullRequest
2 голосов
/ 25 сентября 2019

Я пытаюсь исследовать источник невязки задачи матричного уравнения ( Ax = b ).Чтобы проверить мой ответ, я вычитаю Ax-b , ожидая 0 .Вместо «чистых» нулей я получаю значения в том же порядке, что и машинный эпсилон, что нормально.Проблема в том, что эти остатки кажутся кратными друг другу, поэтому я не уверен, как их интерпретировать.

Я нашел некоторые детали здесь: Проблема вычисления машинного эпсилона , это не такНе проясните, почему кратные значения эпсилона возникают вместо одного или другого.

Я проверил в своей системе, используя np.finfo(float).eps, что привело к 2.220446049250313e-16.Один из остатков, которые я получаю в решении x , идентичен этому значению, однако другой, кажется, составляет половину эпсилона.

Вот код, который я использовал:

# Arbitrary Matrix A and Vector b
A = np.array([[2,-1,0],[1,-2,1],[0,-1,2]])
b = np.array([[1],[0],[1]])

# Solve for Vector x
x = np.linalg.solve(A,b)

# Calculate difference, expected to be column of zeros
diff = A.dot(x) - b
print(diff)

Вот вывод:

Output: 
[[ 0.00000000e+00]
 [-1.11022302e-16]  #-------> Is this machine epsilon...
 [-2.22044605e-16]] #-------> ...or this?

Каково объяснение / толкование этого?Я понимаю, что значения, меньшие, чем эпсилон, все еще могут быть представлены, но в этом случае, почему не оба остатка -1.11022302e-16?

Заранее спасибо!

1 Ответ

3 голосов
/ 25 сентября 2019

Так называемый машинный эпсилон - это просто единица наименьшей точности (ULP) в 1. То есть, это значение позиции младшего значащего бита в представлении 1. Когда в значении есть 53 бита,1 представлена ​​двоичной цифрой 1.000… 000 2 , где после двоичной точки находится 52 нуля.Таким образом, значение позиции самой младшей цифры составляет 2 −52 , а 2 −52 - это ULP, равное 1.

В общем, пусть ULP (x ) обозначает единицу наименьшей точности для x .Обычно формат с плавающей запятой представляет число в виде (−1) s f b e , где b - фиксированная база (два для двоичных форматов, десять для десятичного, 16 для шестнадцатеричного), s - знаковый бит (0 для+, 1 для -), e - показатель степени, а f - значение, с p цифрами, где p - фиксированное значениесумма за формат.Для двоичного кода IEEE-754 p равно 53 для 53 бит.ULP - это значение положения наименьшей точности в значениии, масштабируемое показателем степени, поэтому, если некоторое число x представлено в формате с плавающей запятой со знаковым битом s ,Значение f и показатель степени e , его ULP составляет b 1− p b е .(Я предположил, что формат значения - это одна базовая цифра - b перед точкой радиуса и p -1 цифр после точки радиуса, и поэтому его самая низкая цифра имеетзначение позиции b 1− p . Такие значения находятся в интервале [1, b ).Иногда значения по-разному масштабируются, и показатель степени корректируется для компенсации.Например, в доказательствах может быть полезно, чтобы значения и целые числа были.)

В двоичных форматах ULP (2) = 2 • ULP (1), ULP (½) = ½ • ULP (1)), ULP (¼) = ¼ • ULP (1) и т. Д.

Предположим, что вы вычислили два значения, которые находятся в интервале [1, 2), и это будет, если вычислить с действительным числомарифметические, быть равными, но они были вычислены с арифметикой с плавающей точкой и, как оказалось, немного отличаются.Из-за формата представления они могут отличаться только на кратные ULP (1).Когда вы вычитаете такие числа, вы часто получаете 0, ULP (1), 2 • ULP (1) или другое кратное ULP (1), в зависимости от обстоятельств.Когда два числа, которые были бы идентичны, если их вычислить с помощью арифметики с действительными числами, вычисляются с использованием арифметики с плавающей запятой, они могут испытывать различные ошибки округления в различных частях вычисления.

Если вы вычисляете два значения, которые находятся винтервал [½, 1), они могут отличаться только на кратные ULP (½).

Вот почему вы видите различные кратные или двоичные дроби ULP (1).Это просто артефакт квантования формата с плавающей запятой.

...