Численные методы в ML

torch.linalg.solve использует LAPACK routines, написанные в 1970-х. Нейросеть GPT-4 работает на фундаменте 50-летнего кода для численной линейной алгебры. Понимание этого фундамента - разница между дебаггингом NaN через Stack Overflow и исправлением за 10 минут.

• **Mixed precision (bf16):** BERT, GPT, LLaMA тренируются в bf16 на A100/H100 - 2x ускорение при сохранении стабильности благодаря диапазону fp32. NVIDIA Transformer Engine автоматически управляет этим.
• **LoRA:** fine-tuning LLaMA-65B через rank-16 приближение весов - обучается 0.1% параметров с качеством full fine-tuning. Math base - алгоритм Халько-Мартинссона-Троппа.
• **K-FAC в AlphaFold 2:** оптимизатор второго порядка для структуры белков использовал кронекерову факторизацию матрицы Фишера - 10x меньше шагов чем Adam.

Численная стабильность в глубоком обучении

**torch.linalg.solve использует LAPACK routines, написанные в 1970-х.** Каждая нейросеть стоит на фундаменте 50-летнего кода для численной линейной алгебры - и этот код написан людьми, которые знали о плавающей точке всё, что нужно.

Форматы fp16, bf16 и fp32 - это не просто «меньше бит». Это разные компромиссы между диапазоном и точностью мантиссы. fp16 максимум 65504 - выше идёт Inf. bf16 имеет тот же диапазон что fp32 (±3.4×10³⁸), но меньшую точность мантиссы (7 бит вместо 23).

Типичный loss в трансформерах - порядка 1e-4. В fp16 minimum representable positive normal - 6.1e-5. Если loss ≈ 1e-4, он попадает в denormal range - и градиенты буквально обнуляются. **GradScaler** умножает loss на S=2¹⁵=32768 перед backward, делит градиенты обратно перед шагом оптимизатора.

Softmax(x) при x=[1000, 999, 0]: exp(1000)=Inf в fp32. Решение - вычесть max(x) перед softmax. PyTorch делает это автоматически в F.softmax и F.cross_entropy. В кастомных CUDA kernels это легко пропустить - и получить NaN на всём batch.

L-BFGS, K-FAC и пространства Крылова

Adam и SGD используют только первые производные. L-BFGS и K-FAC учитывают **кривизну** - вторые производные. На задачах с маленькими датасетами или physics-informed neural networks это даёт сходимость за десятки шагов вместо тысяч.

L-BFGS хранит m=10-20 последних пар (sₖ, yₖ) вместо полной матрицы гессиана N×N. Стоимость: O(m·n) памяти против O(n²). Пара: sₖ = xₖ₊₁ - xₖ (шаг), yₖ = ∇f(xₖ₊₁) - ∇f(xₖ) (изменение градиента). Гессиан восстанавливается рекуррентно из этих пар.

**K-FAC (Kronecker-Factored Approximate Curvature):** используется в AlphaFold 2. Матрица Фишера F ≈ A ⊗ G - кронекерово произведение. Обращение: (A ⊗ G)⁻¹ = A⁻¹ ⊗ G⁻¹. Стоимость в 100x дешевле полной матрицы Фишера. Сходимость в 10x меньше шагов чем Adam.

Пространство Крылова Kₖ(A, b) = span{b, Ab, A²b, ..., Aᵏ⁻¹b} - основа методов CG и GMRES для Ax=b. L-BFGS неявно строит аппроксимацию в пространстве Крылова через пары (s, y). Отсюда общая природа всех этих методов.

Рандомизированный SVD и LoRA

Полный SVD матрицы 12288×12288 (размер весов GPT-3) стоит O(mn·min(m,n)) ≈ O(10¹²) операций. Алгоритм Халько-Мартинссона-Троппа (2011) вычисляет rank-k приближение за O(mn·log k) - в 100x быстрее при k ≪ min(m,n). Именно на этом алгоритме построен **LoRA** - стандарт fine-tuning LLM.

Идея: матрица A проецируется на случайное подпространство Ω ∈ ℝⁿˣ(k+p). Если A имеет быстро убывающие сингулярные значения, случайное подпространство поймёт главные направления.

**Алгоритм Халько (рандомизированный SVD):** 1. Ω ← случайная матрица ℝⁿˣ(k+p), p=10 (oversampling) 2. Y = A·Ω (sketch, O(mn·k)) 3. Y = Q·R (QR-разложение Y, Q - ортонормированный базис) 4. B = Qᵀ·A (малая матрица (k+p)×n) 5. B = U_hat·S·Vᵀ (SVD малой матрицы) 6. U = Q·U_hat Ошибка: ||A - Aₖ||₂ ≤ C·σₖ₊₁(A). Power iterations (Y ← (AAᵀ)^q·Y) усиливают зазор между сингулярными значениями.

Автодифференцирование: forward vs reverse mode

PyTorch и JAX вычисляют точные производные через цепное правило на вычислительном графе - не аппроксимацию. Числовое дифференцирование (конечные разности) используется только для верификации. Различие между forward и reverse mode AD - это O(n) против O(1) по числу параметров.

У нейросети n=10⁹ параметров и один скаляр loss. Forward mode потребует 10⁹ прямых проходов - по одному на каждый параметр. Reverse mode даёт ∂L/∂θᵢ для всех n параметров за один обратный проход. Стоимость backprop - около 3-5x стоимости forward pass.

JAX предоставляет grad(grad(f)) - вторые производные за O(n), не O(n²). Это работает через dual numbers в forward mode, вложенный в reverse mode. Используется для meta-learning (MAML) и Physics-Informed Neural Networks.

Ключевые идеи

• **fp16 vs bf16:** bf16 имеет диапазон fp32 - нет переполнения больших loss/gradient values, GradScaler не нужен; fp16 max=65504, нужен GradScaler с scale=2¹⁵
• **Gradient clipping:** если ||g||₂ > max_norm, g ← g·(max_norm/||g||₂). Защита от gradient explosion, стандарт для трансформеров (max_norm=1.0)
• **L-BFGS:** аппроксимация гессиана через m пар (s,y); O(m·n) памяти; требует closure и Wolfe line search; не работает с stochastic mini-batch
• **Рандомизированный SVD (Halko):** rank-k за O(mn·log k); power iterations усиливают зазор сингулярных значений; основа LoRA
• **Reverse mode AD:** один обратный проход → ∂L/∂θ для всех n параметров; точность machine epsilon; gradcheck для верификации реализаций

Связанные темы

Численные методы пронизывают весь ML стек:

• Конечные элементы — Physics-informed neural networks (PINN) используют те же weak formulations - FEM как loss в нейросети
• Численные методы на собеседовании — Типичные вопросы Meta/Google: чем отличается L-BFGS от Adam, как дебаггировать NaN в тренировке

Связанные уроки

• la-01-vectors-intro