Что означают кратные машинному эпсилону?

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

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

Я проверил свою систему, используя 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?

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


python numpy precision floating-point epsilon
person Matija Milenovic    schedule 25.09.2019    source источник

Ответы (1)


arrow_upward
3
arrow_downward

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

В общем, пусть ULP(x) обозначает единицу наименьшей точности для x. Как правило, формат с плавающей запятой представляет число как (−1)s • f • be, где b — фиксированная база (два для двоичных форматов, десять для десятичных, 16 для шестнадцатеричных), s — бит знака (0 для +, 1 для —), e — показатель степени, а f — мантиссы с p цифрами, где p — фиксированная сумма для формата. Для IEEE-754 binary32 p равно 53 для 53 бит. ULP — это значение позиции наименьшей точности в мантиссе, масштабированное по показателю степени, поэтому, если некоторое число x представлено в формате с плавающей запятой со знаковым битом s, мантиссы f и показатель степени e, его ULP равен b1−p • бд. (Я предположил, что формат мантиссы представляет собой одну цифру с основанием b перед точкой счисления и p−1 цифр после точки счисления, и поэтому ее младшая цифра имеет значение позиции b1−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). Это просто артефакт квантования формата с плавающей запятой.

person Eric Postpischil    schedule 25.09.2019
comment
Спасибо за ответ Эрик. Не могли бы вы объяснить, почему ULP(x) связан с интервалом, которому принадлежат значения? Если интервал равен [0,1), означает ли это, что следует ожидать 0 или ULP(1)? - person Matija Milenovic; 25.09.2019
comment
@MatijaMilenovic: ULP пропорционален масштабированию из-за показателя степени. Я обновил ответ по этому поводу. - person Eric Postpischil; 25.09.2019