Когда следует использовать log1p и expm1?

30

У меня есть простой вопрос, который действительно сложен для Google (кроме канонического « Что должен знать каждый учёный-компьютерщик» об арифметической работе с плавающей точкой ).

Когда следует использовать такие функции, как log1pили, expm1а не logи exp? Когда их не следует использовать? Чем отличаются реализации этих функций с точки зрения их использования?

Тим
источник
2
Добро пожаловать в Scicomp.SE! Это очень разумный вопрос, но было бы легче ответить, если бы вы немного объяснили, на что log1p вы ссылаетесь (особенно на то, как это реализовано, поэтому нам не нужно догадываться).
Кристиан Клэйсон
4
Для вещественных аргументов, log1p и expm1 должны использоваться, когда мало, например, когда в точности с плавающей запятой. См., Например, docs.scipy.org/doc/numpy/reference/generated/numpy.expm1.html и docs.scipy.org/doc/numpy/reference/generated/numpy.log1p.html . ( х ) х 1 + х = 1(x)(x)x1+x=1
GoHokies
@ChristianClason спасибо, я имею в виду в основном C ++ std или R, но, как вы спросите, я начинаю думать, что изучение различий в реализациях также было бы очень интересно.
Тим
Вот пример: scicomp.stackexchange.com/questions/8371/…
Хуан М. Белло-Ривас
1
@ user2186862 «когда мало» правильно, но не только «когда с точностью с плавающей запятой» (что происходит для в обычной арифметике двойной точности). Страницы документации, на которые вы ссылаетесь, показывают, что они уже полезны, например, для . 1 + x = 1 x 10 - 16 x 10 - 10x1+x=1x1016x1010
Федерико Полони,

Ответы:

25

Мы все знаем, что подразумевает, что для , мы имеем . Это означает, что если нам нужно вычислить в плавающей точке , дляможет произойти катастрофическая отмена.

exp(x)=n=0xnn!=1+x+12x2+
|x|1exp(x)1+xexp(x)1|x|1

Это может быть легко продемонстрировано в python:

>>> from math import (exp, expm1)

>>> x = 1e-8
>>> exp(x) - 1
9.99999993922529e-09
>>> expm1(x)
1.0000000050000001e-08

>>> x = 1e-22
>>> exp(x) - 1
0.0
>>> expm1(x)
1e-22

Точные значения:

exp(108)1=0.000000010000000050000000166666667083333334166666668exp(1022)1=0.000000000000000000000100000000000000000000005000000

В общем, «точная» реализация expи expm1должна быть правильной, чтобы не более 1ULP (то есть одна единица последнего места). Однако, поскольку достижение этой точности приводит к «медленному» коду, иногда доступна быстрая и менее точная реализация. Например, в CUDA у нас есть expfи expm1f, где fозначает быстрый. Согласно руководству по программированию CUDA C, приложение. Дexpf имеет погрешность 2ULP.

Если вас не волнуют ошибки порядка нескольких ULPS, обычно разные реализации экспоненциальной функции эквивалентны, но помните, что ошибки могут быть где-то спрятаны ... (Помните ошибку Pentium FDIV ?)

Таким образом, довольно ясно, что expm1следует использовать для вычисления для малого . Использование его для общего не вредно, так как ожидается, что он будет точным во всем диапазоне:exp(x)1xxexpm1

>>> exp(200)-1 == exp(200) == expm1(200)
True

(В приведенном выше примере значительно ниже 1ULP, равного , поэтому все три выражения возвращают одно и то же число с плавающей запятой.)1exp(200)

Аналогичное обсуждение справедливо для обратных функций logи, log1pпоскольку для .log(1+x)x|x|1

Стефано М
источник
1
Этот ответ уже содержался в комментариях к вопросу ОП. Однако я счел полезным дать более длинный (хотя и базовый) отчет только для ясности, в надежде, что он будет полезен для некоторых читателей.
Стефано М
Хорошо, но тогда можно просто сделать вывод: «Я всегда могу использовать expm1 вместо exp» ...
Тим
1
@tim ваш вывод неверен: вы всегда можете использовать expm1(x)вместо exp(x)-1. Конечно exp(x) == exp(x) - 1не держится в общем.
Стефано М
ОК, это понятно. И есть ли четкие критерии для ? x1
Тим
1
@Tim нет четкого порога, и ответ зависит от точности реализации с плавающей запятой (и решаемой проблемы). Хотя expm1(x)должно быть с точностью до 1ULP во всем диапазоне , постепенно теряется точность от нескольких ULP при до полного сбоя при , где - машинно-эпсилон. 0x1exp(x) - 1x1x<ϵϵ
Стефано М
1

Чтобы расширить разницу между logи log1pэто может помочь вспомнить график, если логарифм:

Логарифм

Если ваши данные содержат нули, то вы, вероятно, не хотите использовать, logпоскольку они не определены в ноль. И когда приближается к , значение приближается к . Таким образом, если ваши значения близки к , то значение потенциально может быть большим отрицательным числом. Например, и и т. Д. Это может быть полезно, но также может искажать ваши данные в сторону больших отрицательных чисел, особенно если ваш набор данных также содержит числа, которые намного больше нуля.x0ln(x)x0ln(x)ln(1e)=1ln(1e10)=10

С другой стороны, когда приближается к , значение приближается к с положительного направления. Например, и . Таким образом, выдает только положительные значения и устраняет «опасность» больших отрицательных чисел. Обычно это обеспечивает более однородное распределение, когда набор данных содержит числа, близкие к нулю.x0ln(x+1)0ln(1+1e)0.31ln(1+1e10)0.000045log1p

Короче говоря, если набор данных больше , то обычно все в порядке. Но если набор данных имеет числа от до , то обычно лучше.10 1log01log1p

sfmiller940
источник