Итерационные методы подпространств Крылова для гигантских матриц
Как решать системы с миллиардами уравнений?
При дискретизации трехмерных дифференциальных уравнений в частных производных (например, в аэродинамике или гидрогеологии) методом конечных элементов возникают колоссальные системы линейных алгебраических уравнений (СЛАУ). Размерность таких матриц может достигать сотен миллионов и даже миллиардов строк. Использование прямого метода Гаусса для таких гигантов невозможно: даже суперкомпьютеру не хватит оперативной памяти для хранения матрицы, не говоря уже о годах, необходимых для расчетов.
К счастью, такие матрицы являются крайне разреженными. Более 99% их элементов — это нули. Для их хранения используются специальные форматы (например, CSR — Compressed Sparse Row), в которых хранятся только ненулевые элементы и их координаты. Однако прямые методы в процессе исключения неизбежно превращают нули в ненулевые элементы (проблема заполнения). Поэтому для решения таких огромных разреженных СЛАУ применяются исключительно современные итерационные методы, вершиной эволюции которых являются алгоритмы на основе подпространств Крылова.
Подпространства Крылова: выжимание максимума из умножения матрицы на вектор
В 1931 году русский математик Алексей Крылов опубликовал работу, которая заложила основы нового класса методов. Идея состоит в следующем: поскольку мы можем быстро умножать нашу гигантскую разреженную матрицу (назовем ее A) на вектор, почему бы не использовать операцию умножения для построения базиса, в котором мы будем искать решение?
Если мы возьмем начальный вектор невязки (r) и начнем последовательно умножать его на матрицу A, мы получим последовательность векторов: r, Ar, A^2r, A^3r и так далее. Линейная оболочка (пространство, натянутое на эти векторы) называется подпространством Крылова. Гениальность подхода заключается в том, что вместо решения системы в исходном пространстве с миллионами измерений, мы ищем приближенное решение (проекцию) в подпространстве Крылова, размерность которого мала (например, 50 или 100). На каждом шаге (итерации Арнольди или Ланцоша) к подпространству добавляется новый ортогональный вектор, и приближение становится все точнее.
GMRES и Метод сопряженных градиентов (CG)
В зависимости от свойств матрицы A применяются различные алгоритмы подпространств Крылова. Если исходная матрица симметричная и положительно определенная (что типично для задач теплопроводности и теории упругости), используется легендарный Метод сопряженных градиентов (CG). Он использует короткие рекуррентные формулы: для вычисления нового вектора нужно хранить в памяти только данные с предыдущего шага. Это делает метод фантастически быстрым и экономным по памяти.
Однако если матрица несимметричная (что часто бывает в задачах конвекции-диффузии и гидродинамики), метод CG не работает. Для таких матриц в 1986 году Саадом и Шульцем был разработан метод GMRES (Generalized Minimal Residual method). Метод GMRES минимизирует норму невязки на всем подпространстве Крылова. В отличие от CG, он вынужден хранить в памяти все базисные векторы, поэтому на практике применяется с «рестартами» (GMRES(m)): после достижения m итераций накопленные векторы сбрасываются, и процесс начинается заново с текущего наилучшего приближения. Сегодня алгоритмы подпространств Крылова, оснащенные мощными предобуславливателями (preconditioners) типа ILU или Algebraic Multigrid (AMG), являются абсолютным индустриальным стандартом при суперкомпьютерном моделировании.