283 lines
10 KiB
Markdown
283 lines
10 KiB
Markdown
# Численная устойчивость
|
||
|
||
## Математическое число и число в памяти
|
||
|
||
Компьютер хранит большинство вещественных чисел в формате с плавающей точкой:
|
||
|
||
\[
|
||
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`.
|
||
- Математически эквивалентные формулы могут быть вычислительно неравноценны.
|
||
|
||
## Резюме
|
||
|
||
Численная надёжность требует различать чувствительность задачи и качество
|
||
алгоритма. Стабильные сдвиги, логарифмическая область, факторизации, разумное
|
||
масштабирование и проверки конечности важнее механического увеличения точности.
|