Численная линейная алгебра в ML

Netflix Prize (2006-2009): задача рекомендательной системы на матрице 480K × 17K со 100 миллионами оценок, по сути, задача матричного дополнения через SVD. Победитель использовал рандомизированный SVD и iterative методы. Google PageRank: матрица переходов 50 миллиардов × 50 миллиардов, решается итерационным методом степенной итерации (power method). Вся современная ML-инфраструктура стоит на численной линейной алгебре.

• **LLM-inference**: матричные умножения занимают >90% времени, cuBLAS и рандомизированные методы делают это возможным на GPU
• **Gaussian Processes**: Cholesky разложение ковариационной матрицы, узкое место GP при n>10⁴, рандомизированные методы открывают путь к n=10⁶
• **Sparse solvers в FEM**: ANSYS, Abaqus, OpenFOAM используют AMG и GMRES для систем 10⁸ уравнений

Предварительные знания

• Systems of Nonlinear Equations
• Direct Methods: LU, Cholesky, QR

Матричные разложения в ML: SVD, QR, Cholesky

В GPT-4 один прямой проход содержит около 1,8 трлн операций с плавающей точкой - все они сводятся к умножению матриц, основе численной линейной алгебры. ML-модели, это по существу операции с матрицами. PCA, рекомендательные системы, kernel methods, Gaussian processes, все они сводятся к нескольким фундаментальным разложениям. Понимание этих разложений отличает практика от пользователя sklearn.

Наивное решение задачи наименьших квадратов min||Ax-b||² через нормальные уравнения AᵀAx = Aᵀb: - AᵀA имеет число обусловленности κ(AᵀA) = κ(A)² - Если κ(A) = 10⁶, то κ(AᵀA) = 10¹², потеря 12 знаков точности! QR-разложение решает ту же задачу, но работает с Q и R по отдельности, κ(R) = κ(A). Всегда используйте scipy.linalg.lstsq или np.linalg.lstsq вместо явного AᵀA!

Итерационные методы для больших систем: CG, GMRES

Прямые методы (LU, Cholesky) требуют O(n³) операций и O(n²) памяти. Для матриц 10⁶×10⁶ это неприемлемо. Итерационные методы работают только с матрично-векторными произведениями Av, которые часто выполняются за O(nnz) для разреженных матриц.

**CG (Conjugate Gradient)**: только для SPD матриц. Минимизирует ошибку в A-норме. Каждая итерация: одно Av, несколько скалярных произведений, O(n) памяти. **MINRES**: для симметричных (не обязательно положительно определённых). Минимизирует невязку ||r||₂. **GMRES**: для произвольных невырожденных матриц. Минимизирует ||r||₂ в пространстве Крылова. Память растёт с числом итераций → используют restarted GMRES(k). **BiCGSTAB**: компромисс между GMRES и CG, постоянная память, работает для несимметричных. Часто используется в CFD для конвективных уравнений.

Предобусловливатели: ускорение в 10-100×

Сходимость CG и GMRES зависит от числа обусловленности κ(A): число итераций ~ √κ(A) для CG. Для плохо обусловленных задач (например, уравнения Пуассона на мелких сетках) κ ~ 1/h² → тысячи итераций. **Предобусловливатель** M преобразует систему Ax=b в M⁻¹Ax = M⁻¹b с лучшим κ.

**Jacobi**: M = diag(A). Простейший, O(n) памяти. Слабо помогает. **ILU(0)**: неполное LU без заполнения. O(nnz) памяти. Хорош для общих несимметричных задач. **IC(0)** (Incomplete Cholesky): для SPD. Аналог ILU(0) для симметричных задач. **Algebraic Multigrid (AMG)**: иерархическое огрубление матрицы. Почти O(n) для задач, связанных с эллиптическими PDE. Реализации: PyAMG, hypre. **Multigrid (геометрический)**: иерархия сеток для PDE. Сложнее в реализации, но часто O(n log n). Наиболее эффективен для структурированных сеток.

Рандомизированная линейная алгебра

Полное SVD матрицы m×n требует O(mn·min(m,n)) операций. Для матриц 10⁶×10³ это дорого. Рандомизированные методы используют случайные проекции для нахождения приближённого low-rank разложения за O(mn·k), при условии, что матрица имеет быстро убывающие сингулярные числа.

**Random sketching**: заменить большую матрицу A (m×n) на маленькую SA (k×n) через случайную матрицу S (k×m), k ≪ m. Для задачи наименьших квадратов: вместо min||Ax-b||₂ решать min||SAx-Sb||₂. Типы случайных матриц: - Gaussian: S_{ij} ~ N(0,1/k), простейший вариант - Sparse Johnson-Lindenstrauss: большинство элементов 0 → быстрое умножение - SRHT (Subsampled Randomized Hadamard Transform): FFT-структура, O(mn log k) **Теорема Джонсона-Линденштрауса**: случайная проекция с k = O(log(1/δ)/ε²) строками сохраняет все попарные расстояния с точностью ε с вероятностью 1-δ. Фундаментальный результат рандомизированных алгоритмов.

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

• **SVD/QR/Cholesky**: три ключевых разложения ML, SVD для PCA и низкоранговых аппроксимаций, QR для устойчивых least squares, Cholesky для GP и байесовских методов
• **CG/GMRES**: итерационные методы работают только с Av, критично для разреженных матриц, где хранение невозможно; CG для SPD, GMRES для несимметричных
• **Предобусловливатели**: M ≈ A⁻¹ ускоряет CG/GMRES в 10-100×; AMG достигает почти O(n) для задач с эллиптическими PDE
• **Randomized SVD**: O(mnk) вместо O(mn²) за счёт случайной проекции, эффективно когда σ быстро убывают; sklearn TruncatedSVD использует этот подход

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

Численная линейная алгебра, фундамент и для оптимизации, и для ML:

• Собственные значения и итерационные методы — QR-алгоритм, степенная итерация, Arnoldi, это численная линейная алгебра в действии для спектральных задач
• Численная устойчивость — κ(A) и κ(AᵀA) = κ(A)², понимание обусловленности критично для выбора алгоритма
• Параллельные алгоритмы — ScaLAPACK, cuBLAS, SUMMA, параллельные реализации тех самых разложений для GPU и distributed memory

Вопросы для размышления

• Когда имеет смысл использовать iterative solver вместо прямого (LU/Cholesky)?
• Если Cholesky разложение падает с ошибкой «матрица не положительно определена», что это означает в контексте GP или байесовской регрессии?
• Рандомизированный SVD возвращает немного разные результаты при каждом запуске. Как это влияет на воспроизводимость ML-экспериментов?

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

• la-01-vectors-intro