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

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

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

\[
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}\).

```python
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}}.
\]

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

```python
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\), поэтому переполнения нет.

```python
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)}.
\]

```python
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).
\]

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

```python
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}\) математически верна, но явное
построение обратной матрицы делает больше операций и обычно увеличивает ошибку.

```python
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|.
\]

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

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

## Диагностика

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

```python
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`.
- Математически эквивалентные формулы могут быть вычислительно неравноценны.

## Резюме

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