Численные методы: точность, устойчивость, арифметика с плавающей точкой

Вся математика предыдущих статей трека жила в идеальном мире. Вектор в линейной алгебре состоял из настоящих вещественных чисел. Определитель в статье про матрицы вычислялся точно. Вероятность в теории вероятностей была числом из отрезка [0, 1], а не из 2⁶⁴ возможных битовых комбинаций.

Компьютер в этот мир не пускают. У него конечная память, а вещественных чисел континуум — несчётно много. Значит, любая программа, работающая с вещественными числами, работает не с ними, а с их конечной моделью, и каждая арифметическая операция чуть-чуть врёт.

Численные методы — это дисциплина о том, как врать контролируемо. Не «избежать ошибок» (нельзя), а знать их величину, не дать им накапливаться и уметь отличить «мой алгоритм плох» от «эта задача в принципе не решается точнее».

Цена невежества здесь измеряется не в неудобстве, а в деньгах и жизнях: ракета Ariane 5 (1996) взорвалась через 37 секунд после старта из-за неперехваченного переполнения при конвертации 64-битного float в 16-битное знаковое целое; батарея Patriot (1991) промахнулась по «Скаду» в Дахране из-за накопленного за 100 часов дрейфа в 0.34 секунды — представление числа 0.1 в 24-битном регистре давало ошибку примерно 9.5·10⁻⁸ на такт (отчёт GAO IMTEC-92-26).

Часть 0. Три сцены из обычной жизни разработчика

Сцена 1. Тест на равенство. assert 0.1 + 0.2 == 0.3 падает. Джуниор пишет issue «баг в Python». Это не баг: ни 0.1, ни 0.2, ни 0.3 не представимы в двоичной плавающей точке, и сумма двух ближайших приближений оказывается ближайшим приближением к 0.30000000000000004.

Сцена 2. Сумма денег не сходится. Отчёт по 2 млн транзакций расходится с бухгалтерией на 3 копейки. Причина не в логике, а в порядке суммирования float64: сложение вещественных чисел в плавающей точке не ассоциативно, и параллельный reduce по 16 потокам даёт другой результат, чем последовательный цикл.

Сцена 3. Модель обучается-обучается и выдаёт NaN. На 12-й тысяче шагов loss становится nan. Где-то посчитался log(0), или exp(1000) дал inf, или градиент в float16 ушёл в субнормальную зону и обнулился.

Все три — один и тот же корень: вы считали, что работаете с R, а работали с конечным множеством F ⊂ Q. Разберёмся, что это за множество.

Численные методы отвечают за две нижние ветви — и учат оценивать вклад верхних.

Часть 1. Что такое число с плавающей точкой

Интуиция: научная нотация в двоичной системе

Астроном не пишет массу Солнца как 1989100000000000000000000000000 кг, он пишет 1.9891 · 10³⁰. Пять значащих цифр плюс порядок. Точность фиксирована относительно величины: и для массы Солнца, и для массы электрона у нас пять значащих цифр.

Float — ровно это, только основание 2, и всё упаковано в фиксированное число бит.

Строгое определение

Система с плавающей точкой F(β, p, e_min, e_max) — множество чисел вида

    x = ± ( d₀ . d₁ d₂ … d_{p−1} )_β · β^e,

где β — основание (у IEEE 754: β = 2),
      p — точность, число значащих разрядов,
      0 ≤ dᵢ < β,  e_min ≤ e ≤ e_max.

Число нормализовано, если d₀ ≠ 0. В двоичной системе это значит d₀ = 1
всегда — и поэтому его НЕ ХРАНЯТ («скрытый бит»), выигрывая один разряд.

Плюс специальные значения: ±0, ±∞, NaN и субнормальные числа (d₀ = 0, e = e_min).

Стандарт IEEE 754-2019 фиксирует конкретные форматы:

Обратите внимание на bfloat16: у него та же экспонента, что у float32, но втрое меньше мантиссы. Это сознательный размен, сделанный Google для TPU: в глубоком обучении важнее не потерять динамический диапазон (градиенты гуляют на 10 порядков), чем сохранить значащие цифры. Конвертация float32 → bfloat16 — это просто отбрасывание младших 16 бит.

Ключевое следствие: сетка неравномерна

Между 1 и 2 умещается ровно 2⁵² чисел float64. Между 2 и 4 — тоже 2⁵², но на вдвое более широком интервале, значит расстояние между соседями вдвое больше. Между 2⁵² и 2⁵³ соседи отстоят на 1 — целые числа. А выше 2⁵³ целые числа уже не все представимы:

>>> 2**53
9007199254740992
>>> float(2**53) == float(2**53 + 1)   # два разных целых — один float
True

Отсюда практическое правило: не храните идентификаторы в float64, если они могут превысить 2⁵³ ≈ 9·10¹⁵. Классическая боль — JSON: в JavaScript единственный числовой тип это float64, и снежинка Twitter/X (64-битный ID) при парсинге теряет младшие биты. Поэтому такие API отдают ID строкой.

ulp (unit in the last place) — расстояние до соседнего представимого числа. Оно зависит от величины. А вот единичное округление u = eps/2 = 2⁻⁵³ ≈ 1.11·10⁻¹⁶ — постоянно и характеризует относительную точность.

Классы значений и как между ними попасть

Два практических следствия из этой диаграммы.

Субнормальные числа тормозят. Аудиодвижки и физические симуляции, где сигнал экспоненциально затухает к нулю, регулярно проваливаются в субнормальную зону и получают внезапную просадку производительности в десятки раз (микрокод вместо аппаратного пути). Лечение — режим FTZ (flush-to-zero) через _MM_SET_FLUSH_ZERO_MODE или флаг компилятора, либо добавление к сигналу крошечного «dither».

NaN распространяется, но не сигнализирует. По умолчанию IEEE 754 использует quiet NaN: операция не бросает исключение, а тихо возвращает NaN, который просачивается через весь пайплайн. Отладка «где родился первый NaN» — отдельное искусство: в PyTorch есть torch.autograd.set_detect_anomaly(True), в numpy — np.seterr(all='raise'), в C — feenableexcept(FE_INVALID).

Модель арифметики с округлением

Всё дальнейшее строится на одном свойстве, которое гарантирует IEEE 754 для + − × ÷ √:

Каждая базовая операция даёт ПРАВИЛЬНО ОКРУГЛЁННЫЙ точный результат:

    fl(a ⊙ b) = (a ⊙ b)(1 + δ),   |δ| ≤ u,   ⊙ ∈ {+, −, ×, ÷}

где u — единичное округление (2⁻⁵³ для float64).

Это НЕ «результат примерно правильный», а точная гарантия:
машина вычисляет a⊙b в бесконечной точности и округляет один раз.

Из этой модели выводится всё остальное — оценки накопления погрешности, границы для скалярного произведения, теоремы об устойчивости LU-разложения. Именно поэтому IEEE 754 — одно из важнейших инженерных достижений XX века: до него каждая архитектура округляла по-своему, и переносимый численный анализ был невозможен. Историю борьбы за стандарт стоит прочитать у самого Уильяма Кэхена, получившего за него премию Тьюринга.

Чего модель НЕ гарантирует: трансцендентные функции (sin, exp, log, pow). Стандарт лишь рекомендует правильное округление; на практике libm даёт 0.5–1 ulp, но разные реализации (glibc, musl, Intel SVML, CUDA) дают разные последние биты. Отсюда невоспроизводимость результатов между CPU и GPU.

Что ломается: законы алгебры

Множество F с операциями ⊕, ⊗ — не поле и даже не полугруппа в привычном смысле (см. абстрактную алгебру):

# Сложение НЕ ассоциативно
>>> (0.1 + 0.2) + 0.3
0.6000000000000001
>>> 0.1 + (0.2 + 0.3)
0.6

# Умножение НЕ дистрибутивно относительно сложения
>>> a, b, c = 1e16, 1.0, -1e16
>>> (a + b) + c        # 1e16 + 1 округлилось обратно в 1e16
0.0
>>> a + (b + c)
1.0

# x + 1 == x при достаточно большом x
>>> 1e17 + 1 == 1e17
True

# Но кое-что сохраняется:
# ассоциативность и коммутативность умножения — нет,
# зато коммутативность сложения ⊕ — ДА (округление симметрично).

Не-ассоциативность — не академический курьёз. Она означает, что параллельная редукция принципиально даёт другой ответ, чем последовательная, и что менять порядок суммирования — это менять результат. Все споры про «детерминизм обучения нейросетей» упираются сюда: атомарные atomicAdd на GPU складывают в недетерминированном порядке.

Часть 2. Погрешности: откуда, сколько, и когда это катастрофа

Абсолютная, относительная, значащие цифры

Пусть x — точное значение, x̂ — вычисленное.

    абсолютная ошибка:   E_abs = |x̂ − x|
    относительная:       E_rel = |x̂ − x| / |x|   (при x ≠ 0)

x̂ имеет примерно k верных значащих десятичных цифр,
если E_rel ≲ 5 · 10^(−k).

Относительная ошибка — единственная осмысленная мера, когда масштаб задачи неизвестен. Абсолютная ошибка 1 мм катастрофична для литографии чипа и смехотворна для геодезии.

Катастрофическое сокращение — враг номер один

Вычитание двух близких чисел не добавляет ошибки (оно даже точно при a/2 ≤ b ≤ 2a — лемма Стерблца), но обнажает уже имеющуюся: старшие верные разряды уничтожаются, младшие мусорные всплывают наверх.

Разберём руками, в десятичной системе с 6 значащими цифрами:

a = 1.23456789…  →  fl(a) = 1.23457   (ошибка ~4e-6, относительная ~3e-6)
b = 1.23453210…  →  fl(b) = 1.23453   (ошибка ~2e-6, относительная ~2e-6)

точная разность:      a − b = 0.00003579…
вычисленная:      fl(a) ⊖ fl(b) = 0.00004000

Абсолютная ошибка выросла всего до 4e-6 — вычитание её не увеличило.
Но ОТНОСИТЕЛЬНАЯ ошибка: 4e-6 / 3.6e-5 ≈ 12%.
Было 6 верных цифр — стала МЕНЕЕ ОДНОЙ.

Мораль: сокращение не создаёт ошибку, оно её повышает в ранге. Ошибка, которая была шумом в 16-м разряде, становится главным членом.

Классический разбор: корни квадратного уравнения

Формула из школы x = (−b ± √(b²−4ac)) / 2a численно неустойчива, когда b² ≫ 4ac: тогда √(b²−4ac) ≈ |b|, и один из корней считается как разность почти равных чисел.

import math

def roots_naive(a, b, c):
    """Школьная формула. Работает, пока b² не доминирует."""
    d = math.sqrt(b*b - 4*a*c)
    return ((-b + d) / (2*a), (-b - d) / (2*a))

def roots_stable(a, b, c):
    """Устойчивый вариант: считаем «безопасный» корень, второй — через теорему Виета.

    Идея: x1·x2 = c/a. Один корень получаем сложением одинаковых знаков
    (сокращения нет), второй — делением, а не вычитанием.
    """
    d = math.sqrt(b*b - 4*a*c)
    q = -0.5 * (b + math.copysign(d, b))   # знаки b и d совпадают → сложение
    return (q / a, c / q)

print(roots_naive(1.0, 1e8, 1.0))   # (-7.450580596923828e-09, -100000000.0)
print(roots_stable(1.0, 1e8, 1.0))  # (-100000000.0, -1e-08)
# Точный малый корень: -1.0000000000000001e-08
# Наивная формула ошиблась на 25% — при том что все входные данные точны!

Наивный код потерял ноль верных цифр в большом корне и все, кроме одной, в малом. Никакого предупреждения, никакого исключения — просто неверное число. Этот приём (заменить вычитание на деление через Виета) — из Numerical Recipes, §5.6, и он же реализован в LAPACK.

Второй классический разбор: дисперсия одним проходом

«Оптимизация», которую пишут все и которая ломается на реальных данных:

def var_naive(xs):
    """Однопроходная формула E[x²] − (E[x])². Численно ужасна."""
    n = len(xs)
    s1 = sum(xs)
    s2 = sum(x*x for x in xs)
    return (s2 - s1*s1/n) / (n - 1)

def var_two_pass(xs):
    """Два прохода: сначала среднее, потом отклонения. Устойчиво, но два прохода."""
    n = len(xs)
    mean = sum(xs) / n
    return sum((x - mean)**2 for x in xs) / (n - 1)

def var_welford(xs):
    """Алгоритм Уэлфорда: один проход, устойчив, работает потоково.

    Инвариант после k элементов: mean = среднее, M2 = сумма квадратов отклонений.
    Сложность: O(n) по времени, O(1) по памяти.
    """
    n = 0
    mean = 0.0
    M2 = 0.0
    for x in xs:
        n += 1
        delta = x - mean
        mean += delta / n
        M2 += delta * (x - mean)   # намеренно: старая delta × новая
    return M2 / (n - 1)

data = [1e9 + i for i in range(1, 6)]   # 1000000001 … 1000000005
print(var_naive(data))      # 0.0     ← ПОЛНАЯ ЧУШЬ
print(var_two_pass(data))   # 2.5
print(var_welford(data))    # 2.5     ← правильно, и в один проход

s2 ≈ 5·10¹⁸, s1²/n ≈ 5·10¹⁸, их разность — 12 в 19-м разряде, которого в float64 просто нет. Формула математически верна и вычислительно мертва.

Уэлфорд (1962) вычитает среднее на лету, работая с малыми числами. Это тот самый метод, который стоит внутри pandas.Series.var, torch.nn.BatchNorm и потоковых агрегаторов метрик. Оригинал: B. P. Welford, Note on a Method for Calculating Corrected Sums of Squares and Products, Technometrics 4(3), 1962.

Компенсированное суммирование

Сложение миллиона чисел наивным циклом даёт погрешность, растущую как O(n·u) в худшем случае и O(√n·u) при случайных знаках. Компенсация Кэхена держит её на уровне O(u), почти независимо от n:

def kahan_sum(xs):
    """Компенсированное суммирование Кэхена.

    c хранит «потерянный хвост» — то, что не поместилось в мантиссу s.
    На следующем шаге хвост возвращается в игру.
    Стоимость: 4 операции вместо 1, погрешность ~2u независимо от n.
    """
    s = 0.0
    c = 0.0
    for x in xs:
        y = x - c              # вернуть накопленный долг
        t = s + y              # здесь теряются младшие биты y
        c = (t - s) - y        # ровно то, что потерялось (знак инвертирован)
        s = t
    return s

def neumaier_sum(xs):
    """Улучшение Ноймайера: корректно работает, когда очередной x
    БОЛЬШЕ текущей суммы (случай, на котором классический Кэхен ломается)."""
    s = 0.0
    c = 0.0
    for x in xs:
        t = s + x
        if abs(s) >= abs(x):
            c += (s - t) + x
        else:
            c += (x - t) + s
        s = t
    return s + c

xs = [1e16, 1.0, -1e16]
acc = 0.0
for x in xs:
    acc += x
print(acc)                # 0.0  ← наивно
print(kahan_sum(xs))      # 0.0  ← Кэхен ТОЖЕ ошибается на этом входе
print(neumaier_sum(xs))   # 1.0  ← верно
import math; print(math.fsum(xs))   # 1.0 ← точное (Шевчук), гарантированно

Три уровня качества, три цены:

Кстати, sum() в CPython 3.12+ сам использует компенсацию Ноймайера для float — поэтому пример выше в интерактивной сессии даст 1.0, если написать sum(xs), но 0.0 в явном цикле. Это ровно та деталь, которая ломает интуицию людей, читающих «одинаковый» код.

numpy.sum использует попарное суммирование (pairwise) — рекурсивное деление пополам. Это даёт O(log n) рост погрешности бесплатно, без лишних flop. Детали — в документации numpy по np.sum.

Часть 3. Обусловленность задачи против устойчивости алгоритма

Это центральное понятие всей дисциплины, и его чаще всего путают.

Определения

Пусть задача — вычислить f(x). Алгоритм вернул ŷ.

ПРЯМАЯ ошибка (forward error):     ‖ŷ − f(x)‖ / ‖f(x)‖
    «насколько ответ отличается от правильного»

ОБРАТНАЯ ошибка (backward error):  min ‖Δx‖ / ‖x‖  по всем Δx с f(x+Δx) = ŷ
    «для каких слегка изменённых данных мой ответ был бы ТОЧНЫМ»

ЧИСЛО ОБУСЛОВЛЕННОСТИ задачи:
    cond(f, x) = lim_{ε→0} sup_{‖Δx‖≤ε‖x‖} (‖f(x+Δx) − f(x)‖/‖f(x)‖) / (‖Δx‖/‖x‖)

Для гладкой скалярной f:   cond(f, x) = |x · f′(x) / f(x)|

Главное неравенство:
    прямая ошибка  ≲  cond(f, x) · обратная ошибка

Алгоритм называется ОБРАТНО УСТОЙЧИВЫМ, если обратная ошибка ~ O(u)
для любых входных данных.

Интуиция: разделение ответственности

Представьте, что вы переводчик. Обратная ошибка — насколько вы переврали исходный текст. Число обусловленности — насколько эта область чувствительна к переводу (юридический договор против прогноза погоды). Прямая ошибка — насколько плохо всё кончилось.

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

Так же и в вычислениях: устойчивость — свойство алгоритма, обусловленность — свойство задачи. Обратно устойчивый алгоритм на плохо обусловленной задаче ОБЯЗАН давать большую прямую ошибку, и это не его вина. Требовать от него точности — значит требовать невозможного.

Число обусловленности на пальцах: скалярный случай

# cond(f, x) = |x f'(x) / f(x)|

# 1) Вычитание близких чисел: f(x) = x − 1 около x = 1
#    f'(x) = 1, cond = |x / (x−1)| → ∞ при x → 1.
#    ЗАДАЧА плохо обусловлена. Никакой алгоритм не спасёт.

# 2) f(x) = sqrt(x):  cond = |x · (1/(2√x)) / √x| = 1/2.
#    Идеально обусловлена: относительная ошибка входа даже уменьшается.

# 3) f(x) = exp(x):   cond = |x|.
#    При x = 1000 относительная ошибка входа усиливается в 1000 раз.
#    Вот почему softmax по «сырым» логитам — плохая идея.

# 4) f(x) = log(x) около x = 1: cond = |1/log(x)| → ∞.
#    Отсюда специальная функция log1p(x) = log(1+x), точная у нуля:
import math
x = 1e-16
print(math.log(1 + x))   # 0.0        — катастрофа: 1+x округлилось до 1
print(math.log1p(x))     # 1e-16      — верно
# Симметрично: expm1(x) = exp(x) − 1
print(math.exp(1e-16) - 1)   # 0.0
print(math.expm1(1e-16))     # 1e-16

log1p, expm1, hypot, fma, copysign — это не украшения стандартной библиотеки, а готовые устойчивые переформулировки известных ловушек. Их наличие в C99/math — прямое следствие численного анализа. Если вы видите в коде math.sqrt(x*x + y*y) — это баг: при x = 1e200 будет inf, тогда как math.hypot(x, y) вернёт правильный ответ.

Число обусловленности матрицы

Для системы Ax = b (см. статью о матрицах):

cond(A) = ‖A‖ · ‖A⁻¹‖ ,   в спектральной норме:  cond₂(A) = σ_max / σ_min

Оценка чувствительности решения:
    ‖Δx‖/‖x‖  ≤  cond(A) · ( ‖ΔA‖/‖A‖ + ‖Δb‖/‖b‖ )

Правило большого пальца:
    если cond(A) ≈ 10^k, то в решении теряется примерно k
    верных десятичных цифр из 16 доступных во float64.

Отношение σ_max/σ_min — прямая связь с SVD: обусловленность это «во сколько раз самое сильное направление растяжения превосходит самое слабое».

Проверим на матрице Гильберта H[i][j] = 1/(i+j+1) — эталонно плохо обусловленной:

import numpy as np

def hilbert(n):
    return np.array([[1.0 / (i + j + 1) for j in range(n)] for i in range(n)])

for n in (5, 10, 15):
    H = hilbert(n)
    x_true = np.ones(n)          # известный точный ответ
    b = H @ x_true               # правая часть
    x_hat = np.linalg.solve(H, b)
    print(f"n={n:2d}  cond={np.linalg.cond(H):.3e}  "
          f"max|x̂−x|={np.max(np.abs(x_hat - x_true)):.3e}")

# n= 5  cond=4.766e+05  max|x̂−x|=2.091e-12
# n=10  cond=1.602e+13  max|x̂−x|=2.428e-05
# n=15  cond=3.605e+17  max|x̂−x|=9.227e+02      ← ответ 923 вместо 1

Смотрите на предсказательную силу правила: при cond ≈ 1e13 из 16 цифр остаётся 3 — и правда ошибка 1e-5 при значениях порядка 1. При cond ≈ 3.6e17 > 1/u система численно вырождена, и LAPACK честно возвращает мусор. LAPACK при этом не ошибся — он обратно устойчив, его обратная ошибка порядка u. Задача такая.

Мини-эксперимент: чувствительность видна глазами

import numpy as np
A = np.array([[1.0, 1.0],
              [1.0, 1.0000001]])           # почти вырожденная
print(f"cond = {np.linalg.cond(A):.3e}")   # 4.000e+07

print(np.linalg.solve(A, np.array([2.0, 2.0000001])))   # [1. 1.]
print(np.linalg.solve(A, np.array([2.0, 2.0000002])))   # [0. 2.]
# Изменили b в 8-м знаке — решение изменилось ПОЛНОСТЬЮ.

Это не численная проблема. Это геометрия: две прямые почти параллельны, и точка их пересечения безумно чувствительна к наклону. Никакая арифметика не сделает такую задачу устойчивой — только смена постановки (регуляризация, наименьшие квадраты, дополнительные данные).

Зачем нужен пивотинг: устойчивость руками

Гауссово исключение без выбора главного элемента не является обратно устойчивым. Демонстрация в четыре строки:

import numpy as np

def lu_solve_nopivot(A, b):
    """Гаусс без пивотинга. Сложность O(n³) по времени, O(n²) по памяти."""
    A = A.astype(float).copy()
    b = b.astype(float).copy()
    n = len(b)
    for k in range(n - 1):
        for i in range(k + 1, n):
            m = A[i, k] / A[k, k]        # если A[k,k] крошечный — m гигантский
            A[i, k:] -= m * A[k, k:]
            b[i] -= m * b[k]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (b[i] - A[i, i+1:] @ x[i+1:]) / A[i, i]
    return x

A = np.array([[1e-20, 1.0],
              [1.0,   1.0]])
b = np.array([1.0, 2.0])

print(lu_solve_nopivot(A, b))   # [0. 1.]  ← неверно
print(np.linalg.solve(A, b))    # [1. 1.]  ← верно (LAPACK делает пивотинг)

Что произошло: множитель m = 1/1e-20 = 1e20, второе уравнение превращается в 1e20·x₂ ≈ 1e20, исходная информация x₁ + x₂ = 2 полностью затирается округлением. Матрица при этом прекрасно обусловлена (cond ≈ 2.6) — виноват целиком алгоритм.

Частичный выбор главного элемента (переставить строки так, чтобы |A[k,k]| был максимален в столбце) держит все множители |m| ≤ 1 и делает метод устойчивым на практике. Строгая теория — Wilkinson, Rounding Errors in Algebraic Processes (1963), и глава 9 в Higham.

Часть 4. Итерационные методы: сходимость и когда остановиться

Нахождение корня: бисекция против Ньютона

Бисекция. Дано: f непрерывна, f(a)·f(b) < 0.
Инвариант: корень всегда внутри [a, b]. Длина интервала делится пополам.
Сходимость ЛИНЕЙНАЯ: error_{k+1} ≈ 0.5 · error_k.
За 52 шага интервал шириной 1 сжимается до ulp. Гарантирована ВСЕГДА.

Метод Ньютона. x_{k+1} = x_k − f(x_k)/f′(x_k).
Сходимость КВАДРАТИЧНАЯ вблизи простого корня: error_{k+1} ≈ C · error_k².
Число верных цифр удваивается за шаг. Но глобальной гарантии НЕТ.

ПСЕВДОКОД: гибрид (так устроены brentq, TOMS748, реальные solver-ы)

вход: f, [a, b] с разными знаками, tol
пока (b − a) > tol:
    предложить шаг Ньютона (или секущей / обратной квадратичной интерполяции)
    если предложенная точка ВНЕ [a, b] или сходимость медленная:
        взять середину интервала           # откат к гарантии
    вычислить f в новой точке, сузить скобку
вернуть середину

Итог: скорость Ньютона + надёжность бисекции.

import math

def newton_sqrt2():
    """Ньютон для x² − 2 = 0. Смотрим на удвоение верных цифр."""
    f  = lambda x: x*x - 2
    fp = lambda x: 2*x
    x = 1.0
    for k in range(6):
        print(f"k={k}  x={x!r:22}  err={abs(x - math.sqrt(2)):.3e}")
        x = x - f(x) / fp(x)

newton_sqrt2()
# k=0  x=1.0                    err=4.142e-01
# k=1  x=1.5                    err=8.579e-02
# k=2  x=1.4166666666666667     err=2.453e-03
# k=3  x=1.4142156862745099     err=2.124e-06
# k=4  x=1.4142135623746899     err=1.595e-12
# k=5  x=1.4142135623730951     err=0.000e+00

Ошибки: 4e-1 → 9e-2 → 2e-3 → 2e-6 → 2e-12 → 0. Показатель ошибки примерно удваивается — это и есть квадратичная сходимость. Пять итераций до машинной точности против пятидесяти двух у бисекции.

Но: Ньютон теряет квадратичность на кратных корнях, потому что там f′ тоже обращается в ноль:

g  = lambda x: (x - 1)**2      # корень кратности 2
gp = lambda x: 2*(x - 1)
x = 2.0
for _ in range(5):
    x = x - g(x)/gp(x)
    print(x)     # 1.5, 1.25, 1.125, 1.0625, 1.03125 — просто деление пополам

Сходимость упала до линейной с коэффициентом 1/2. Лечение — модифицированный Ньютон x − m·f/f′ при известной кратности m, или работа с f/f′ вместо f.

Как правильно останавливать итерации

Самая частая ошибка в самописных солверах — критерий abs(x_new - x_old) < 1e-10. Он ломается в обе стороны: для x ~ 1e6 это недостижимо (ulp там уже 1e-10), для x ~ 1e-12 — срабатывает сразу.

def converged(x_new, x_old, rtol=1e-12, atol=1e-300):
    """Смешанный критерий: относительный + абсолютный (страховка около нуля).
    Ровно так устроены np.allclose, scipy.optimize и большинство ODE-солверов."""
    return abs(x_new - x_old) <= atol + rtol * abs(x_new)

Правила, выстраданные практикой:

1. Никогда не сравнивайте float через == — кроме проверки на точное «то же самое значение» (например, сторожевые константы) и на NaN через x != x.
2. Всегда ставьте лимит итераций. Итерация, которая «не сошлась», должна вернуть ошибку, а не крутиться вечно и не молча вернуть последнее значение.
3. Не требуйте точности лучше u·cond — вы просите невозможного, и цикл не завершится.
4. Для тестов сравнивайте в ULP (math.ulp, numpy.testing.assert_allclose(..., rtol=...)), а не в абсолютных величинах.

Часть 5. Дискретизация: производные, интегралы, ОДУ

Численная производная и оптимальный шаг

Здесь два вида погрешности воюют друг с другом, и это лучший в дисциплине пример компромисса.

Правая разность:      D₊(h) = (f(x+h) − f(x)) / h
    ошибка усечения (Тейлор):  ~ (h/2)·|f″|          — падает с h
    ошибка округления:         ~ 2u·|f| / h          — РАСТЁТ при h → 0
    сумма минимальна при  h* ≈ 2√(u · |f|/|f″|) ≈ √u ≈ 1.5e-8
    достижимая точность:  ~ √u ≈ 1e-8   (половина цифр потеряна!)

Центральная разность: D₀(h) = (f(x+h) − f(x−h)) / (2h)
    ошибка усечения:  ~ (h²/6)·|f‴|
    оптимум:  h* ≈ ∛u ≈ 6e-6,  точность ~ u^(2/3) ≈ 4e-11

import math
f, x0, exact = math.exp, 1.0, math.e
print(f"{'h':>8} {'|forward−e|':>14} {'|central−e|':>14}")
for h in (1e-2, 1e-5, 1e-8, 1e-11, 1e-14):
    fwd = (f(x0 + h) - f(x0)) / h
    ctr = (f(x0 + h) - f(x0 - h)) / (2*h)
    print(f"{h:8.0e} {abs(fwd-exact):14.3e} {abs(ctr-exact):14.3e}")

#        h    |forward−e|    |central−e|
#    1e-02      1.364e-02      4.530e-05
#    1e-05      1.359e-05      5.859e-11   ← центральная у своего оптимума
#    1e-08      6.603e-09      6.603e-09   ← правая у своего оптимума
#    1e-11      3.263e-05      1.043e-05   ← округление берёт верх
#    1e-14      9.338e-03      9.338e-03   ← полный шум

График ошибки — буква V: слева усечение, справа округление. Уменьшать шаг «для точности» после дна V — вредно. Это ловушка, в которую попадают все, кто пишет численный градиент для проверки backprop; правильный ответ — использовать центральную разность с h ≈ 1e-5 для float64, либо перейти на комплексно-шаговое дифференцирование (Im f(x + ih)/h), у которого вообще нет вычитания и потому нет ошибки округления.

Обыкновенные дифференциальные уравнения и жёсткость

Задача Коши:  y′ = F(t, y),  y(0) = y₀.

Явный Эйлер:    y_{k+1} = y_k + h·F(t_k, y_k)      — просто, дёшево, УСЛОВНО устойчив
Неявный Эйлер:  y_{k+1} = y_k + h·F(t_{k+1}, y_{k+1})  — на каждом шаге решать уравнение,
                                                          зато БЕЗУСЛОВНО устойчив

Тестовое уравнение y′ = λy, λ < 0 (затухание):
    явный:   y_{k+1} = (1 + hλ)·y_k   →  устойчив ⟺ |1 + hλ| < 1 ⟺ h < 2/|λ|
    неявный: y_{k+1} = y_k / (1 − hλ) →  устойчив при ЛЮБОМ h > 0

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

import math

LAM, T, Y0 = -1000.0, 0.05, 1.0
exact = Y0 * math.exp(LAM * T)          # 1.93e-22

def euler_explicit(h):
    y = Y0
    for _ in range(int(T / h)):
        y = y + h * LAM * y
    return y

def euler_implicit(h):
    """Для линейного F решаем шаг аналитически: y_{k+1}(1 − hλ) = y_k."""
    y = Y0
    for _ in range(int(T / h)):
        y = y / (1 - h * LAM)
    return y

# граница устойчивости явной схемы: h < 2/1000 = 0.002
for h in (0.0021, 0.002, 0.0005):
    print(f"h={h}  explicit={euler_explicit(h):.4e}")
# h=0.0021  explicit=-8.9543e+00   ← РАСХОДИТСЯ, знак скачет, модуль растёт
# h=0.002   explicit=-1.0000e+00   ← ровно на границе: колеблется, не затухает
# h=0.0005  explicit=7.8886e-31    ← устойчиво

print(f"implicit h=0.0021: {euler_implicit(0.0021):.4e}")   # 4.9967e-12 — устойчиво
print(f"exact:             {exact:.4e}")                    # 1.9287e-22

Обратите внимание: при h = 0.0021 явная схема выдаёт −8.95 вместо 1.9e-22. Не «неточно» — качественно неверно: решение должно монотонно затухать, а получается растущая осцилляция. Неявная схема при том же шаге даёт 5e-12 — не идеально точно (шаг-то грубый), но качественно правильно: затухание есть.

Это фундаментальный водораздел: точность и устойчивость — разные вещи. Схема может быть точной и неустойчивой (тогда она бесполезна), устойчивой и неточной (тогда она хотя бы не врёт про характер решения). Теорема Лакса для линейных задач связывает их: согласованность + устойчивость ⟺ сходимость.

В проде это значит: если ваша симуляция «взрывается», первым делом смотрите не на точность метода, а на область устойчивости и спектр якобиана (см. собственные значения). Готовые инструменты — scipy.integrate.solve_ivp с методами Radau или BDF для жёстких задач, RK45/DOP853 для нежёстких (документация SciPy).

Часть 6. Где это всплывает в реальном коде

Деньги: почему float для валюты — почти всегда баг

0.1 не представим в двоичной системе, а вся бухгалтерия человечества десятичная. Накопление ошибок при миллионах операций даёт расхождения, которые аудит не прощает.

from decimal import Decimal, getcontext, ROUND_HALF_UP

# ПЛОХО
total = 0.0
for _ in range(1000):
    total += 0.1
print(total)                    # 99.9999999999986  ← не 100.0

# ХОРОШО: целые копейки/сатоши/микроценты
total_cents = sum(10 for _ in range(1000))
print(total_cents / 100)        # 100.0

# ИЛИ: десятичная арифметика с явным правилом округления
getcontext().rounding = ROUND_HALF_UP
total = sum(Decimal("0.1") for _ in range(1000))
print(total)                    # 100.0

Тот же принцип в БД: в PostgreSQL для денег берут NUMERIC(19,4), а не DOUBLE PRECISION. NUMERIC — десятичная арифметика произвольной точности, медленнее в разы, но точная в десятичном смысле (документация PostgreSQL по числовым типам). Ловушка: SELECT 0.1::float8 + 0.2::float8 = 0.3::float8 вернёт false, и индекс по float-колонке с равенством работать «как ожидается» не будет.

Ещё тоньше: GROUP BY + SUM(float8) в распределённой БД (ClickHouse, BigQuery, Spark) даёт разные результаты между запусками, потому что порядок слияния партиций недетерминирован. Если ваш дашборд «прыгает» в последнем знаке — вот причина.

Машинное обучение

log-sum-exp. Прямое вычисление log(Σ exp(zᵢ)) переполняется уже при z ≈ 710 для float64 и при z ≈ 89 для float32:

import math

def logsumexp_naive(zs):
    return math.log(sum(math.exp(z) for z in zs))

def logsumexp_stable(zs):
    """Тождество: log Σ exp(zᵢ) = m + log Σ exp(zᵢ − m), где m = max zᵢ.
    После вычитания максимума все экспоненты ≤ 1 — переполнение невозможно,
    а слагаемое exp(0) = 1 гарантирует, что сумма не занулится."""
    m = max(zs)
    return m + math.log(sum(math.exp(z - m) for z in zs))

zs = [1000.0, 1001.0, 1002.0]
try:
    logsumexp_naive(zs)
except OverflowError as e:
    print("naive:", e)                  # naive: math range error
print("stable:", logsumexp_stable(zs))  # stable: 1002.4076059644444

Ровно этот приём стоит внутри torch.logsumexp, scipy.special.logsumexp и любой корректной реализации softmax и cross-entropy. Именно поэтому в PyTorch есть CrossEntropyLoss, принимающий логиты, а не вероятности: композиция log_softmax + nll устойчива, а log(softmax(x)) — нет.

Смешанная точность. Обучение в float16 даёт 2× память и до 8× скорость на тензорных ядрах, но float16 имеет минимальное нормальное значение 6e-5 — градиенты проваливаются в субнормальную зону и обнуляются. Решение — loss scaling: умножить loss на S ≈ 2¹⁵ перед backward, поделить градиенты на S перед шагом оптимизатора. Аккумуляторы и веса при этом держат в float32 (Micikevicius et al., Mixed Precision Training, arXiv:1710.03740).

Не сравнивайте эмбеддинги на равенство. Косинусная близость двух «одинаковых» векторов, посчитанных на CPU и на GPU, будет 0.9999999 — и хеш-дедупликация по float-значению не сработает.

Компиляторы и оптимизации

-ffast-math (GCC/Clang) разрешает компилятору притворяться, что float — это R: переставлять слагаемые, заменять x/y на x·(1/y), считать, что NaN и −0.0 не бывает. Ускорение может быть двукратным, а результат — неверным, причём невоспроизводимо (зависит от уровня оптимизации и от того, влез ли цикл в вектор).

Конкретный побочный эффект: -ffast-math включает -ffinite-math-only, после чего проверка if (isnan(x)) может быть удалена компилятором как заведомо ложная. Ваша защита от NaN просто исчезает из бинарника.

Отдельная история — FMA (fused multiply-add): a*b + c с одним округлением вместо двух. Точнее, быстрее — но меняет результат. Один и тот же исходник, скомпилированный с -march=native на машине с FMA и без, даст разные последние биты. Если у вас есть регрессионные тесты с точным сравнением чисел — они начнут падать при смене железа в CI.

Историческая аномалия: x87 считал во внутренних 80-битных регистрах, и double x = a*b; мог дать разный результат в зависимости от того, попала ли x в регистр или в память. Это отравляло жизнь два десятилетия, пока SSE2 не сделал 64-битные вычисления честными. Подробности — в классической статье Гольдберга What Every Computer Scientist Should Know About Floating-Point Arithmetic, которую стоит прочитать целиком хотя бы раз в жизни.

Геометрия и вычислительная надёжность

Предикаты вроде «лежит ли точка слева от прямой» вычисляются через знак определителя. При почти вырожденной конфигурации знак получается случайным — и алгоритм построения выпуклой оболочки или триангуляции Делоне падает с нарушением инварианта или зацикливается. Это не редкость, а норма: любой геометрический код на float обязан либо использовать адаптивные точные предикаты (Shewchuk, Robust Predicates), либо считать в целых числах на сетке.

Криптография: почему там нет float

В криптографии (см. связь с абстрактной алгеброй) вся арифметика — модульная целочисленная, точная по построению. Float там запрещён двумя причинами: он неточен, и он не постоянен по времени — субнормальные операции длятся дольше, что даёт побочный канал. Известны атаки по времени на реализации, где float просочился в обработку секретных данных.

Часть 7. Типичные заблуждения

«Float64 достаточно точен, чтобы не думать». 16 цифр — много, пока задача не имеет cond = 1e14. Матрица Гильберта 15×15 съедает всё за один вызов solve.

«Проблема в том, что float неточен». Нет: базовые операции IEEE 754 идеально точны — они дают правильно округлённый результат. Проблема в накоплении и в сокращении, то есть в структуре ваших формул.

«Чем меньше шаг, тем точнее». Только до дна V-образной кривой. Дальше округление растёт быстрее, чем падает усечение.

«Округлю всё через Decimal и буду спокоен». Decimal точен для десятичных дробей, но Decimal(1)/Decimal(3) всё равно округляется, и все проблемы с обусловленностью и сокращением остаются. Меняется основание, не природа.

«Сравню с эпсилоном — abs(a-b) < 1e-9». Абсолютный порог работает только около единицы. Для a ≈ 1e12 он недостижим, для a ≈ 1e-15 он объявляет равными всё подряд. Нужен смешанный критерий.

«NaN != NaN — это баг языка». Это требование IEEE 754. NaN означает «не число», и никакие два «не числа» не равны. Отсюда sorted() со смешанными NaN даёт мусор — компаратор нарушает транзитивность, и порядок становится неопределённым.

«Если тесты зелёные, численный код верный». Численный баг обычно проявляется на конкретной конфигурации данных (близкие числа, большой разброс масштабов, почти вырожденная матрица). Тестируйте на плохих входах: 1e-300, 1e300, inf, nan, -0.0, почти равные числа, матрицы с известной обусловленностью.

«Ассоциативность сложения можно предполагать». Именно на этом предположении компилятор с -ffast-math векторизует ваш цикл и меняет ответ.

Часть 8. Практический чек-лист

Прежде чем считать численный код готовым:

1. Знаю ли я cond своей задачи? Хотя бы порядок. Для линейных систем — np.linalg.cond. Для скалярной функции — |x f′(x)/f(x)|. Без этого нельзя сказать, хорош ли ответ.
2. Есть ли в коде вычитание близких чисел? Ищите a - b, где a ≈ b. Есть ли эквивалентная формула без него (Виета, log1p, expm1, hypot, Уэлфорд)?
3. Отмасштабированы ли данные? Признаки в диапазонах [0,1] и [0,1e9] в одной матрице — это cond ≈ 1e9 из воздуха. Центрирование и нормализация — не только про ML, это ещё и численная гигиена.
4. Длинные аккумуляторы компенсированы? Больше 1e6 слагаемых — попарное как минимум, Ноймайер если важно.
5. Есть ли эталон? Пересчитайте на маленьком примере в mpmath с 50 знаками или в Decimal и сравните. Это единственный честный способ узнать реальную ошибку.
6. Проверены ли граничные значения? inf, nan, -0.0, субнормальные, 1e-308, 1e308.
7. Итерации ограничены и имеют смешанный критерий останова?
8. Тесты сравнивают через allclose с осмысленными rtol/atol, а не через ==?
9. Задокументирована ли ожидаемая точность? «Функция возвращает результат с относительной ошибкой не хуже 1e-12 при cond(A) < 1e6» — вот полезный докстринг численной функции.
10. Разобрался ли я, откуда взялось расхождение CPU/GPU/версий BLAS, прежде чем объявить его «шумом»?

Мини-итог

• Компьютер работает не с R, а с конечным множеством F — неравномерной сеткой с постоянной относительной точностью u ≈ 1.1e-16 для float64.
• IEEE 754 гарантирует правильное округление базовых операций: fl(a⊙b) = (a⊙b)(1+δ), |δ| ≤ u. Это фундамент всех оценок.
• В F не работают ассоциативность и дистрибутивность. Порядок операций — часть спецификации, а не деталь реализации.
• Главный источник катастроф — сокращение близких чисел: оно не создаёт ошибку, а повышает её в ранге до главного члена.
• Обусловленность — свойство задачи, устойчивость — свойство алгоритма. прямая ошибка ≲ cond · обратная ошибка. Хороший алгоритм отвечает только за второй множитель.
• Дискретизация всегда даёт компромисс усечение/округление с V-образной кривой ошибки. Есть оптимальный шаг, и он не «как можно меньше».
• Устойчивость схемы для ОДУ важнее её точности: неточный, но устойчивый ответ качественно верен, точный и неустойчивый — бесполезен.
• В проде это проявляется как: деньги в NUMERIC, softmax через log-sum-exp, loss scaling в fp16, отказ от -ffast-math там, где важны NaN, точные предикаты в геометрии.

Источники

• Nicholas J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002 — библия дисциплины, страница книги на SIAM.
• David Goldberg, What Every Computer Scientist Should Know About Floating-Point Arithmetic, ACM Computing Surveys, 1991 — полный текст.
• IEEE 754-2019 Standard for Floating-Point Arithmetic.
• Lloyd N. Trefethen, David Bau, Numerical Linear Algebra, SIAM, 1997 — лучшее введение в обусловленность и обратную устойчивость.
• W. Kahan, Lecture Notes on the Status of IEEE 754 — история стандарта от автора.
• Press, Teukolsky, Vetterling, Flannery, Numerical Recipes, 3rd ed. — numerical.recipes (алгоритмы отличные, лицензия недружелюбная — читайте, не копируйте).
• Jonathan Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates.
• Micikevicius et al., Mixed Precision Training, arXiv:1710.03740.
• Документация Python math — fsum, ulp, nextafter, isclose, log1p, expm1, hypot.
• Floating Point Guide — короткий FAQ для тех, кому надо объяснить коллеге за 5 минут.
• Отчёт GAO IMTEC-92-26 о Patriot в Дахране и отчёт комиссии по Ariane 5 Flight 501.

Что дальше

Мы разобрались, как считать надёжно, когда противник — конечная арифметика. Дальше противник становится разумным: он выбирает стратегию, реагирует на вашу и максимизирует свою выгоду. Как рассуждать о таких системах строго — в следующей статье трека.

Теория игр: равновесия, механизмы, приложения в распределённых системах