ai_learning/common/foundations/numerical_stability.md

Численная устойчивость

Математическое число и число в памяти

Компьютер хранит большинство вещественных чисел в формате с плавающей точкой:

x=(-1)^s\,m\,2^e,

где знак s, мантисса m и показатель e имеют конечное число битов. Поэтому представимо лишь конечное множество чисел. Результат операции обычно округляется:

\operatorname{fl}(x\circ y)=(x\circ y)(1+\delta),
\qquad |\delta|\lesssim\varepsilon_{\mathrm{mach}},

если нет переполнения или потери нормализации. Для float64 \varepsilon_{\mathrm{mach}}\approx2.22\cdot10^{-16}, для float32 — примерно 1.19\cdot10^{-7}.

import numpy as np

for dtype in (np.float32, np.float64):
    info = np.finfo(dtype)
    print(dtype.__name__, info.eps, info.max, info.tiny)

Ошибка задачи и ошибка алгоритма

Обусловленность — свойство математической задачи: насколько выход чувствителен к малому изменению входа. Устойчивость — свойство алгоритма: насколько вычисленный ответ похож на точный ответ для слегка возмущённых данных.

Относительное число обусловленности скалярной функции:

\kappa_f(x)=\left|\frac{x f'(x)}{f(x)}\right|.

Из df\approx f'(x)dx следует

\frac{|df|}{|f(x)|}
\approx\kappa_f(x)\frac{|dx|}{|x|}.

Если \kappa велико, никакой алгоритм не восстановит много точных цифр из неточного входа. Для линейной системы \mathbf{A}\mathbf{x}=\mathbf{b}:

\kappa_2(\mathbf{A})
=\|\mathbf{A}\|_2\|\mathbf{A}^{-1}\|_2
=\frac{\sigma_{\max}}{\sigma_{\min}}.

Катастрофическая потеря значимости

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

\sqrt{x+1}-\sqrt{x}

при большом x — разность близких чисел. Домножим на сопряжённое:

\sqrt{x+1}-\sqrt{x}
=\frac{(\sqrt{x+1}-\sqrt{x})(\sqrt{x+1}+\sqrt{x})}
{\sqrt{x+1}+\sqrt{x}}
=\frac1{\sqrt{x+1}+\sqrt{x}}.

Вторая форма численно надёжнее.

import math

x = 1e16
unstable = math.sqrt(x + 1.0) - math.sqrt(x)
stable = 1.0 / (math.sqrt(x + 1.0) + math.sqrt(x))
print(unstable, stable)

Переполнение, исчезновение и логарифмы

e^x быстро переполняется, а e^{-x} исчезает до нуля. Для

\operatorname{logsumexp}(\mathbf{x})
=\log\sum_i e^{x_i}

выберем m=\max_i x_i:

\log\sum_i e^{x_i}
=\log\left(e^m\sum_i e^{x_i-m}\right)
=m+\log\sum_i e^{x_i-m}.

Теперь все показатели x_i-m\leq0, поэтому переполнения нет.

import numpy as np

x = np.array([1000.0, 1001.0, 999.0])
m = x.max()
value = m + np.log(np.exp(x - m).sum())
print(value)

Тот же сдвиг стабилизирует softmax:

\operatorname{softmax}(\mathbf{x})_i
=\frac{e^{x_i-m}}{\sum_j e^{x_j-m}}.

Сдвиг не меняет результат, потому что общий множитель e^{-m} сокращается.

Вероятности лучше хранить в логарифмах

Произведение многих вероятностей стремится к нулю:

p=\prod_{i=1}^{n}p_i.

В логарифмической области:

\log p=\sum_{i=1}^{n}\log p_i.

Сумма устойчивее произведения и превращает максимизацию правдоподобия в максимизацию логарифма. Для выражений около нуля используйте специальные функции:

\log(1+x)\quad\longrightarrow\quad\texttt{log1p(x)},

e^x-1\quad\longrightarrow\quad\texttt{expm1(x)}.

x = 1e-16
print(np.log(1 + x), np.log1p(x))
print(np.exp(x) - 1, np.expm1(x))

Суммирование и порядок операций

Сложение с плавающей точкой не ассоциативно:

(a+b)+c\ne a+(b+c).

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

def kahan_sum(values):
    total = 0.0
    correction = 0.0
    for value in values:
        adjusted = value - correction
        updated = total + adjusted
        correction = (updated - total) - adjusted
        total = updated
    return total

values = [1e16, 1.0, 1.0, -1e16]
print(sum(values), kahan_sum(values))

На практике math.fsum и многие библиотечные редукции используют улучшенные стратегии, но точная реализация зависит от библиотеки и устройства.

Не вычисляйте обратную матрицу без необходимости

Для решения \mathbf{A}\mathbf{x}=\mathbf{b} запись \mathbf{x}=\mathbf{A}^{-1}\mathbf{b} математически верна, но явное построение обратной матрицы делает больше операций и обычно увеличивает ошибку.

A = np.array([[3.0, 1.0], [1.0, 2.0]])
b = np.array([9.0, 8.0])
x = np.linalg.solve(A, b)
print(x)

Для задачи наименьших квадратов не следует бездумно использовать нормальные уравнения:

\mathbf{X}^\top\mathbf{X}\mathbf{w}
=\mathbf{X}^\top\mathbf{y}.

Их число обусловленности примерно возводится в квадрат:

\kappa_2(\mathbf{X}^\top\mathbf{X})
=\kappa_2(\mathbf{X})^2.

np.linalg.lstsq использует более устойчивые факторизации.

Масштаб признаков

Если один столбец \mathbf{X} имеет порядок 10^{-6}, а другой 10^6, поверхность потерь сильно вытянута. Это ухудшает обусловленность гессиана и замедляет градиентный спуск. Стандартизация

z_j=\frac{x_j-\mu_j}{\sigma_j}

не гарантирует идеальную обусловленность, но часто существенно её улучшает. Параметры \mu_j,\sigma_j оцениваются только по обучающей выборке.

Допуски вместо точного равенства

Вещественные результаты сравнивают с абсолютным и относительным допуском:

|a-b|\leq \text{atol}+\text{rtol}|b|.

print(0.1 + 0.2 == 0.3)
print(np.isclose(0.1 + 0.2, 0.3, rtol=1e-12, atol=1e-15))

Абсолютный допуск важен около нуля, относительный — для больших величин. Допуск должен следовать из масштаба задачи, а не скрывать ошибку.

Диагностика

Проверяйте данные и промежуточные результаты:

def describe_numeric(name, array):
    a = np.asarray(array)
    print(
        name,
        "shape=", a.shape,
        "dtype=", a.dtype,
        "finite=", np.isfinite(a).all(),
        "min=", np.nanmin(a),
        "max=", np.nanmax(a),
    )

Полезные признаки проблемы: nan, inf, резкий рост нормы параметров, нулевые вероятности, огромный разброс масштабов и большое np.linalg.cond(A).

Типичные ошибки

• Добавление \varepsilon повсюду может изменить модель и замаскировать причину; сначала выведите устойчивую форму.
• float64 уменьшает округление, но не исправляет плохо обусловленную задачу.
• Нормализация не должна использовать статистики тестовой выборки.
• Проверка только конечного ответа пропускает место появления первого nan.
• Математически эквивалентные формулы могут быть вычислительно неравноценны.

Резюме

Численная надёжность требует различать чувствительность задачи и качество алгоритма. Стабильные сдвиги, логарифмическая область, факторизации, разумное масштабирование и проверки конечности важнее механического увеличения точности.