Научные вычисления
Системы линейных уравнений
Цели урока
- Знать LU-разложение и его преимущество при множестве правых частей
- Понимать когда Cholesky быстрее LU и почему он стабилен без pivoting
- Объяснить разницу между Jacobi, Gauss-Seidel и CG
- Понять форматы разреженных матриц и экономию памяти
Предварительные знания
2003 год. Boeing симулирует краш-тест нового крыла в виртуальной аэродинамической трубе. FEM-модель - 2 миллиона уравнений. Инженеры решают 400 таких систем за день, меняя параметры нагрузки. Без LU-разложения и разреженных матриц на это ушли бы годы. С ними - часы.
- **FEM в инженерии** - расчёт напряжений в конструкции сводится к решению Ax=b с разреженной матрицей жёсткости на миллионы уравнений
- **Компьютерная графика** - симуляция ткани, жидкости и дыма решает разреженные системы на каждом кадре
- **PageRank Google** - ранжирование миллиардов страниц через итеративное решение системы с разреженной матрицей связей
От Чолеского до разреженных матриц
1910 - Cholesky создаёт метод для геодезических расчётов, погибает на войне в 1918. 1952 - Hestenes и Stiefel независимо изобретают метод сопряжённых градиентов (CG). 1970-е - появляются форматы CSR/CSC для хранения разреженных матриц, открывая путь к FEM-расчётам промышленного масштаба. Одна цепочка идей от геодезиста 1910 года до авиасимуляторов 2020-х.
LU-разложение
Авиасимулятор Boeing считает нагрузку на крыло. FEM-модель - матрица жёсткости 50 тысяч на 50 тысяч элементов, одно решение Ax = b на каждую конфигурацию крена. Задача одна - разложить A один раз, потом решать с разными b за секунды. Это и есть суть **LU-разложения**: A = LU, где L нижнетреугольная, U верхнетреугольная.
**PA = LU** - версия с перестановками строк (pivoting). P - матрица перестановок, гарантирующая численную стабильность. Без pivoting деление на маленький диагональный элемент катастрофически усиливает ошибки округления.
Сложность LU-разложения - $O(2n^3/3)$. Главное преимущество: разложив A один раз, можно решить Ax = b для множества разных b за $O(n^2)$ каждый. Именно так PyTorch.solve() работает при обучении с L-BFGS - одна матрица Гессе, сотни правых частей.
| Операция | Сложность | Когда нужно |
|---|---|---|
| LU-разложение | O(2n^3/3) | Один раз для матрицы A |
| Forward substitution | O(n^2) | Для каждого нового b |
| Back substitution | O(n^2) | Для каждого нового b |
| Полное solve | O(2n^3/3) | Если b всего одно |
В чём главное преимущество LU-разложения перед прямым solve для нескольких правых частей?
Разложение Cholesky
Ковариационные матрицы в статистике, матрицы жёсткости в FEM, kernel matrices в machine learning - все они симметричные и положительно определённые (SPD). Для таких матриц существует разложение **вдвое быстрее** LU с вдвое меньшей памятью: **Cholesky** A = LL^T.
**Положительно определённая** матрица - для любого ненулевого вектора x выполняется $x^T A x > 0$. Все собственные значения строго положительны. Gaussian Process в ML - это ковариационная матрица, и именно Cholesky даёт сэмплирование из GP за $O(n^3/3)$.
| Свойство | LU | Cholesky |
|---|---|---|
| Требование к A | Квадратная, невырожденная | Симметричная, положительно определённая |
| Формула | PA = LU | A = LL^T |
| Сложность | O(2n^3/3) | O(n^3/3) |
| Числовая стабильность | Нужен pivoting | Стабильна без pivoting |
| Память | Две матрицы L, U | Одна матрица L |
Cholesky вдвое быстрее LU и требует вдвое меньше памяти. Pivoting не нужен - разложение гарантированно стабильно для SPD матриц. В PyTorch: `torch.linalg.cholesky` - именно так считают posterior в Gaussian Process и сэмплируют из multivariate normal.
Andre-Louis Cholesky (1875-1918)
Французский офицер-геодезист Cholesky разработал метод в 1910 году для решения систем уравнений в геодезических измерениях. Погиб в Первой мировой войне в 1918 году. Метод опубликован посмертно коллегой в 1924 году - через 6 лет после гибели автора.
Для какой матрицы Cholesky-разложение гарантированно работает?
Итеративные методы
Матрица $10^6 \times 10^6$. LU требует $10^{18}$ операций - при 10 TFLOPS это 100 000 секунд. Нереально. Но в реальных задачах (FEM, CFD, графы) матрица **разреженная**: 99.99% элементов - нули. **Итеративные методы** начинают с приближения и улучшают его, умножая только на ненулевые элементы.
**Jacobi** - обновляет все компоненты параллельно, используя значения предыдущей итерации. **Gauss-Seidel** - обновляет последовательно, используя уже обновлённые компоненты. Gauss-Seidel обычно сходится быстрее, но не параллелизуется - Jacobi легко запускается на GPU.
**Conjugate Gradient (CG)** - самый важный итеративный метод для SPD матриц. Сходится за n шагов (в точной арифметике), на практике - за $O(\sqrt{\kappa(A)})$ шагов с предобуславливателем, где $\kappa$ - число обусловленности. Именно CG используется в scipy.sparse.linalg для FEM-задач.
| Метод | Тип матрицы | Скорость сходимости | Параллелизм |
|---|---|---|---|
| Jacobi | Диагонально доминантная | Медленная | Высокий (GPU) |
| Gauss-Seidel | Диагонально доминантная | Быстрее Jacobi в ~2x | Низкий |
| CG | SPD | O(sqrt(cond(A))) итераций | Средний |
| GMRES | Любая невырожденная | Зависит от спектра | Средний |
Почему для матрицы $10^6 \times 10^6$ используют итеративные методы, а не LU?
Разреженные матрицы
Матрица $10^6 \times 10^6$ в плотном формате - 8 TB в float64. Но трёхдиагональная матрица того же размера содержит лишь 3 миллиона ненулевых элементов. Отношение: 3 из $10^{12}$. Хранить нужно только ненулевые - это **разреженный** формат, и он экономит сотни тысяч раз.
**CSR** (Compressed Sparse Row) - хранит значения, индексы столбцов и указатели строк. Быстрый для умножения Ax и строчного доступа. **CSC** (Compressed Sparse Column) - аналог по столбцам, быстрый для решения треугольных систем. PageRank Google работает именно с CSR: матрица переходов 10 млрд страниц хранится в разреженном формате.
| Формат | Хранение | Быстрая операция | Типичное применение |
|---|---|---|---|
| CSR | По строкам | A @ x, срезы строк | Итеративные солверы |
| CSC | По столбцам | Срезы столбцов, solve | Прямые солверы (UMFPACK) |
| COO | Тройки (i, j, val) | Быстрое построение | Сборка матрицы из элементов |
| LIL | Списки списков | Быстрая вставка | Инкрементальное построение |
Практика: создавать матрицу в COO или LIL (удобно добавлять элементы), затем конвертировать в CSR для вычислений. Никогда не вызывать `.toarray()` для больших разреженных матриц - это попытка создать плотную матрицу размером в терабайты.
numpy.linalg.solve подходит для любых систем линейных уравнений
np.linalg.solve - прямой метод O(n^3), работающий с плотными матрицами. Для систем размером 10^6 и более нужны итеративные методы (CG, GMRES) в сочетании с разреженными форматами (scipy.sparse).
Матрица 10^6 x 10^6 в плотном формате занимает 8 TB памяти, а LU-разложение требует 10^18 операций. Итеративные методы для разреженных матриц работают с O(nnz) памяти и сходятся за сотни итераций.
Трёхдиагональная матрица $10^6 \times 10^6$ содержит ~3 миллиона ненулевых элементов. Сколько памяти сэкономит CSR по сравнению с плотным хранением?
Ключевые идеи
- **LU-разложение** PA=LU - универсальный прямой метод O(n^3); выгоден при множестве правых частей b
- **Cholesky** A=LL^T - вдвое быстрее LU для SPD матриц, не нужен pivoting
- **Итеративные методы** (Jacobi, Gauss-Seidel, CG) - единственный вариант для n > 10^5
- **Разреженные матрицы** (CSR/CSC) хранят только ненулевые элементы, экономя в сотни тысяч раз
Связанные темы
Системы линейных уравнений - ядро численных методов:
- NumPy и линейная алгебра — np.linalg.solve - прямой метод; здесь разбираем что происходит под капотом
- Введение в научные вычисления — Системы уравнений возникают при дискретизации моделей (FEM, разностные схемы)
- Линейное программирование — Simplex-метод выполняет LU-разложение на каждом pivot-шаге
Вопросы для размышления
- Почему Cholesky не нуждается в pivoting, а LU - нуждается? Что в структуре SPD матриц гарантирует стабильность?
- Если проектировать солвер для FEM-сетки 1000x1000x1000 (10^9 узлов), какую комбинацию методов выбрать?
- Почему спектральный радиус итерационной матрицы должен быть меньше 1 для сходимости?
Связанные уроки
- sci-02 — NumPy линейная алгебра - базис, без него нет смысла в LU
- sci-01 — Дискретизация моделей порождает системы Ax=b
- sci-04 — Собственные значения строятся поверх решения линейных систем
- nm-01 — Численные методы идут глубже: обусловленность, нормы, погрешности
- ml-13 — SVM через QP - квадратичное программирование опирается на LU/CG
- opt-03 — LP-солверы используют LU для построения симплекс-шагов
- dm-01 — Графовые матрицы - разреженные структуры того же типа
- nm-03