Неожиданное поведение log1p numpy

Документы numpy говорят:

Для входных данных с действительным знаком log1p является точным также для x, настолько малого, что 1 + x == 1 с точностью с плавающей запятой.

Кажется, это косвенный способ сказать, что они «отказались» от сложных входных данных. Действительно, в простом Python (без cmath.log1p) в Windows:

>>> import cmath
>>> p = 1e-14 * (1 + 1j)
>>> p
(1e-14+1e-14j)
>>> cmath.log(1 + p)
(9.992007221626409e-15+9.9999999999999e-15j)

что очень близко к тому, что вы получили. Итак, похоже, что numpy.log1p(c) для сложного c добавляет 1,0 к c перед вызовом numpy.log().

Результат «должен быть» близок к p - p**2/2 для такого маленького p (только первые два члена разложения Тейлора log(1+p) около p=0):

>>> p - p**2/2
(1e-14+9.9999999999999e-15j)

который по сути предоставляется модулем расширения mpmath:

>>> import mpmath
>>> p = mpmath.mpc(1, 1) * 1e-14
>>> p
mpc(real='1.0e-14', imag='1.0e-14')
>>> mpmath.log1p(p)
mpc(real='1.0e-14', imag='9.9999999999999006e-15')

Так что это похоже на (недостаточно документированное) ограничение текущей реализации numpy в log1p().

Поэтому, если вам нужно использовать для этого numpy, вам придется предоставить собственную альтернативу.

Что numpy делает

Под обложкой numpy.log1p(c) преобразует c+1 в полярную форму, а затем использует:

log(r * e**(j*t)) = log(r) + j*t

чтобы построить результат. Итак, если c = a + b*j, результат имеет

• реальная часть: log(hypot(a + 1, b))
• часть изображения: atan2(b, a + 1)

Если |a| настолько мало, что a + 1 теряет часть (или всю!) информации из-за округления, не повезло.

Предупреждение о точности

Рассмотрим комплекс c = -5,884298465609164e-91 + 1,0848313692199447e-45j.

Это было сгенерировано «случайным образом» во время тестирования. К сожалению, ни один из предлагаемых подходов не справится с этой задачей, за исключением настройки mpmath на использование точности в несколько сотен бит.

Почему: пусть R — действительная часть, а J — мнимая часть. Настоящая часть log1p(c) — это log(hypot(R+1, J)). Вот! ;-)

>>> import mpmath
>>> from mpmath import log, hypot, mpf, workprec. mpc
>>> R = -5.884298465609164e-91
>>> J = 1.0848313692199447e-45
>>> hypot(R+1, J)
mpf('1.0')
>>> with workprec(400):
...     hypot(R+1, J)
mpf('1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000005884295498218099788929146684839')

Поскольку он продолжает вести журнал, они дают очень разные результаты.

Это потому, что hypot(1+R. J), в свою очередь, sqrt((1 + R)^^2 + J**2) = sqrt(1 + R**2 + 2*R + J**2), и на самом деле нам не помогает знать, что log(sqrt(1+f)) примерно равно f/2, когда f очень мало.

Забудьте R**2, потому что в данном случае он очень маленький. Нас ранит 2*R + J**2:

>>> 2*R
-1.1768596931218329e-90
>>> J**2
1.17685909964362e-90
>>> 2*R + J**2
-5.934782129408062e-97

2*R и J**2 настолько близки к отрицаниям друг друга, что почти семь старших десятичных цифр сокращаются, когда мы их складываем. С большей точностью эти биты заменяются ведущими битами той части J**2, которая была потеряна при округлении:

>>> with workprec(200):
...     2*mpf(R) + mpf(J)**2
mpf('-5.9347821292576139709922469679882199756565314433969258094957309e-97')

Это соответствует только 53-битному результату в первых 10 цифрах. Без дополнительной точности конечные биты J**2 просто теряются навсегда.

>>> c = mpc(R, J)
>>> with workprec(400): # as good as it gets
...     mpmath.log1p(c)
mpc(real='-2.9673910646288069854961234839941099878282657216984629047478654331954101753008213969808976137411482956635728432029405164533e-97', imag='1.0848313692199446892733858946092178721423949541818749887436188379566165369930349882511849209687548440293467928718096798443e-45')

>>> mpmath.log1p(c) # last 4+ digits are trash
mpc(real='-2.9673910646284523e-97', imag='1.0848313692199447e-45')

>>> mpmath.log(c + 1) # total trash
mpc(real='5.8842954982180997e-91', imag='1.0848313692199447e-45')

>>> c2 = c + 1 # this trick is a little worse than plain log1p
>>> c * mpmath.log(c2) / (c2 - 1)
mpc(real='-2.9673910647040311e-97', imag='1.0848313692199447e-45')

>>> xxtest(c) #  # 4-term Taylor series the same
mpc(real='-2.9673910647040311e-97', imag='1.0848313692199447e-45')

В этом труднее разобраться, чем в одномерном (реальнозначном log1p) случае: он «глубокий», а не «поверхностный», потому что катастрофическое сокращение происходит не в том, что мы вычисляем, а под прикрытием в недрах реализация. Недостаточно использовать приемы, которые сохраняют полную точность R, несмотря на добавление 1 к R, потому что в данном случае проблема не в этом: вместо этого мы не сохраняем достаточную точность J**2, чтобы компенсировать возможность того, что его ведущее значение биты могут быть потеряны из-за отмены при добавлении к нему 2*R. Но даже для того, чтобы добраться до этих конечных битов J**2, требуется какой-то способ получить эффект от использования дополнительной точности.

Последующие действия

Мне еще предстоит увидеть надежную log1p(complex) реализацию, не позволяющую работать mpmath с абсурдно (в контексте) высокой точностью.

Частично это похоже на ошибку в mpmath.log1p(). Я обнаружил случай, когда каждая часть реальной части результата была неверной. Итак, открыл отчет об ошибке. Похоже, это тот случай, когда часть реализации временно немного повышает точность, чтобы компенсировать ошибки округления, не осознавая, что это увеличение убедит более позднюю часть реализации выбрать неправильную ветвь. В частности, этот случай был бы хорош, если бы он использовал разложение в ряд Тейлора с двумя членами, но он ошибочно полагает, что точность достаточно велика, поэтому вместо этого можно использовать log(c+1) с удвоенной точностью. Но это не так. Даже утроенной точности будет недостаточно.

Таким образом, пока это не будет исправлено, возможно, что точность mpmath придется увеличить до нескольких тысяч бит, чтобы охватить все случаи.

Кстати, я написал свой собственный lop1p(complex) на Python. Это единственный надежный вариант, который я знаю ;-) Раньше я писал коммерческие математические библиотеки в рамках своей повседневной работы, поэтому у меня высокие стандарты «надежности». В частности, меня не особо волнует относительная погрешность норм. Вместо этого я считаю, что существует бесконечно точный правильный результат, а строгая мера ошибки - это то, на сколько ULP (единиц на последнем месте) действительная и мнимая части удалены от этого идеала. Это неумолимо.

Даже в «обычных» случаях другие реализации, которые я видел, обычно дают результаты, которые отличаются на 50 ULP от идеальных (суммирование абсолютных ошибок ULP в реальной и мнимой частях).

Реализация, которую я сейчас имею, ни разу не достигла общей ошибки 3 ULP в более чем 200 миллионах рандомизированных тестовых случаев. Это худший результат на сегодняшний день:

new max 2.894273405102144 at (0.0019532098478106533-5.356595846417641e-17j)
   (0.0019513048136734744  -5.3461536863894766e-17j) got
   (0.0019513048136734748  -5.346153686389476e-17j)  ideal

Так что это территория «в худшем случае небольшая разница в последних отображаемых цифрах». Достижение этого было полностью вопросом контроля распространения ошибок в реальной части и требовало эмуляции повышенной точности. Я считаю, что для достижения большего результата потребуется заменить платформенные функции atan2(), log() и log1p() с действительными значениями на всегда правильно округленные.

Я бы предложил использовать здесь ряд Тейлора из трех или четырех членов. Вы используете Numba, поэтому получить что-то одновременно точное и производительное довольно просто.

Вот как это можно реализовать:

import numba as nb


@nb.vectorize(fastmath=True)
def log1p_acc(x):
    if np.abs(x) >= 1e-3:
        # np.log1p is accurate for large values of x
        return np.log1p(x)
    else:
        terms = 4
        x_pow = x
        sign = 1
        sum_ = 0
        for i in range(1, terms + 1):
            # Note: use * (1 / i) to allow optimizer to avoid fdiv
            sum_ += sign * x_pow * (1 / i)
            sign *= -1
            x_pow *= x
        return sum_

Чтобы изменить количество терминов, вы можете изменить переменную terms.

Точность

Давайте измерим, насколько точна эта функция.

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

Я выбрал реализацию mpmath.log1p в качестве эталонной реализации и установил точность 1000 цифр.¹

Далее нам нужно задаться вопросом, в какой области нам важна точность. Я предположил, что вас интересуют только небольшие числа, поэтому я сгенерировал выборки из логарифмически-равномерного распределения, охватывающего диапазон от 1e-30 до 1. По сути, это предполагает, что вас волнует диапазон от 1e-30 до 1e-29 настолько же, насколько вас волнует примерно в диапазоне от 0,1 до 1. Вот пример того, как выглядит логарифмически равномерное распределение в виде гистограммы с линейным и логарифмическим масштабом. https://en.wikipedia.org/wiki/Reciprocal_distribution#/media/File:Reciprocal_Histogram.svg

Используя это логарифмически-равномерное распределение, я отбирал образцы из одного из четырех случаев одинаковое количество раз:

1. Реальные и сложные – это независимые логарифмические однородные выборки.
2. Реальный и комплексный – это одна и та же логарифмическая однородная выборка.
3. Действительное число логарифмически однородно, а комплексное равно нулю.
4. Действительное значение равно нулю, а комплексное является логарифмически равномерным.

Это формирует мой тестовый набор для оценки производительности.

В-третьих, я выбрал метрику ошибок. По некоторым показателям ваш результат с np.log1p не так уж и плох: он отличается от истинного примерно на 10^-17, что примерно соответствует машинному эпсилону. В относительном выражении это плохо – скидка около 0,07%.

По этой причине, я думаю, вас больше интересует относительная ошибка. Я измерил это по формуле error = abs(true - pred)/abs(true), где abs() — комплексное абсолютное значение.

Затем я сравнил свою функцию, log1p SciPy и log1p NumPy. На следующем графике показана относительная ошибка в зависимости от абсолютного значения входных данных для трех различных способов ее расчета.

В этом сюжете есть несколько примечательных моментов:

• Для некоторых входных данных относительная ошибка NumPy составляет почти 100%. Например, код np.log1p((1e-16+1e-30j)) возвращает значение с нулем в качестве реального члена.
• SciPy намного точнее, чем NumPy — точность большинства входных данных составляет примерно 1e-16.
• От x=1e-30 до 1e-4 ошибка пользовательского метода составляет менее 1e-16.
• При x=1e-4 ошибка пользовательского метода превышает SciPy. Если бы он не переключился обратно на NumPy при x=1e-3, ошибка будет продолжать расти, особенно по мере выхода за пределы интервала сходимости.

Вот код, используемый для создания этого графика:

import numba as nb
import numpy as np
import scipy.special as sc
import scipy
import matplotlib.pyplot as plt
import mpmath
mpmath.mp.dps = 1000

plt.rcParams["figure.figsize"] = (12,8)



def generate_log_distribution(N):
    """Generate small numbers spanning many orders of magnitude"""
#     return 10 ** np.random.uniform(-3, -30, size=N)
    return scipy.stats.loguniform.rvs(1e-30, 1e-0, size=N)

test_set_size = 100000 // 4
test_set1 = generate_log_distribution(test_set_size) + 1j * generate_log_distribution(test_set_size)
identical_real_complex = generate_log_distribution(test_set_size)
test_set2 = identical_real_complex + 1j * identical_real_complex
test_set3 = generate_log_distribution(test_set_size) + 0j
test_set4 = 1j * generate_log_distribution(test_set_size) + 0
test_set = np.concatenate([test_set1, test_set2, test_set3, test_set4])
test_set_ref = [mpmath.log1p(c) for c in test_set]


@nb.vectorize(fastmath=True)
def log1p_acc(x):
    if np.abs(x) >= 1e-3:
        # np.log1p is accurate for large values of x
        return np.log1p(x)
    else:
        terms = 4
        x_pow = x
        sign = 1
        sum_ = 0
        for i in range(1, terms + 1):
            # Note: use (1 / i) to allow optimizer to avoid fdiv
            sum_ += sign * x_pow * (1 / i)
            sign *= -1
            x_pow *= x
        return sum_


functions = [
    ("NumPy log1p", np.log1p),
    ("SciPy log1p", sc.log1p),
    ("Custom log1p", log1p_acc),
]
for func_name, func in functions:
    x = []
    y = []
    for i in range(len(test_set)):
        number = test_set[i]
        logged = func(number)
        # Do error calculation in arbitrary precision
        logged = mpmath.mpc(logged)
        out = logged - test_set_ref[i]
        out_abs = mpmath.fabs(out)
        number_abs = mpmath.fabs(number)
        relative_err = out_abs / mpmath.fabs(test_set_ref[i])
        x.append(number_abs)
        y.append(relative_err)
    print(f"Report for {func_name}")
    print(f"Average relative error: {float(np.mean(y))}")
#     print("Worst", np.max(y), "for", test_set[np.argmax(y)])
    print()
    plt.scatter(x, y, label=func_name, s=1, alpha=0.1)
    plt.xscale('log')
    plt.yscale('log')
    plt.xlabel('Norm of x')
    plt.ylabel('Relative error in log1p(x)')
    leg = plt.legend()
    for lh in leg.legend_handles:
        lh.set_alpha(1)
        lh.set_sizes([20])

¹ Обратите внимание, что установка dps на 1000 означает только то, что промежуточные вычисления выполняются с точностью до 1000 знаков. Это не обязательно означает, что окончательный результат имеет точность до тысячи цифр. Это может быть точно только вполовину меньше. Я не проверял, насколько численно стабилен mpmath.log1p().

Производительность

Я сравнил производительность этой функции с версиями SciPy и NumPy, взяв log1p всего тестового набора. Это примерно в 2 раза быстрее, чем NumPy, и в 3 раза быстрее, чем SciPy. На оценку log1p требуется около 14 наносекунд.

Код:

print("NumPy log1p")
%timeit np.log1p(test_set)
print("SciPy log1p")
%timeit sc.log1p(test_set)
print("Custom log1p")
%timeit log1p_acc(test_set)

Выход:

NumPy log1p
2.98 ms ± 63.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
SciPy log1p
3.69 ms ± 5.58 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Custom log1p
1.42 ms ± 2.74 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Чтобы сделать это быстро, я применил следующие идеи:

• Я начал с серии Тейлор и превратил ее в цикл. Вместо возведения в степень для вычисления каждого члена я умножил предыдущий член на x.
• Деление с плавающей запятой обычно происходит немного медленнее, чем умножение с плавающей запятой. По этой причине я переписал sign * x_pow / i как sign * x_pow * (1 / i)
• Поскольку terms является константой, Numba достаточно умен, чтобы полностью развернуть этот цикл, что делает это так же быстро, как явное запись ряда Тейлора.
• Я использовал fastmath, который позволяет Numba изменить порядок выполнения математических вычислений.
• Я обнаружил, что использование переменной sign происходит немного быстрее, чем выбор -= или += с индексом цикла.