Предобуславливание (Preconditioning): ключ к решению гигантских СЛАУ
Узкое место итерационных методов Крылова
Ранее мы обсуждали, что гигантские разреженные системы линейных алгебраических уравнений (СЛАУ), возникающие из метода конечных элементов или разностей, невозможно решить прямым методом Гаусса из-за проблемы заполнения нулей. Индустрия использует итерационные методы подпространств Крылова, такие как Метод сопряженных градиентов (CG) для симметричных матриц и GMRES для несимметричных.
Однако эффективность этих методов напрямую зависит от одного фундаментального математического свойства матрицы — числа обусловленности (Condition Number). Число обусловленности характеризует «растянутость» спектра матрицы (отношение максимального собственного значения к минимальному). В реальных физических задачах (особенно при сильных перепадах свойств материалов, например, металл и воздух, или при сильном измельчении сетки) число обусловленности стремится к бесконечности. Матрица становится «плохо обусловленной». В такой ситуации метод сопряженных градиентов начинает топтаться на месте: невязка колеблется, и алгоритму требуются сотни тысяч итераций для сходимости, что делает расчет недопустимо долгим.
Математическая трансформация системы
Для преодоления этого вычислительного тупика применяется предобуславливание (Preconditioning). Идея заключается в преобразовании исходной СЛАУ (Ax = b) в другую, математически эквивалентную систему, которая имеет точно такое же решение, но гораздо лучшее число обусловленности. Для этого вводится матрица предобуславливателя (обозначим ее M).
Мы умножаем левую и правую части системы на обратную матрицу предобуславливателя: M^(-1)Ax = M^(-1)b. Теперь алгоритм Крылова решает новую систему с матрицей M^(-1)A. Идеальным предобуславливателем была бы сама матрица A (тогда M^(-1)A превратилась бы в единичную матрицу, и алгоритм сошелся бы за 1 итерацию). Но обращение A — это именно то, чего мы избегаем. Следовательно, матрица M должна удовлетворять двум противоречивым требованиям: она должна быть «максимально похожа» на матрицу A (чтобы сгруппировать собственные значения), и при этом решение системы My = z должно вычисляться алгоритмически очень быстро и дешево (так как это придется делать на каждой итерации).
Популярные алгоритмы: ILU и диагональное масштабирование
Самым простым является диагональный предобуславливатель (метод Якоби), где M — это просто матрица, состоящая только из диагональных элементов исходной матрицы A. Обратить такую матрицу тривиально (просто поделить на диагональные числа), но ее эффективность для сложных задач крайне мала.
Золотым стандартом в индустрии стало Неполное LU-разложение (Incomplete LU factorization, ILU). Алгоритм пытается сделать стандартное разложение матрицы A на нижнюю и верхнюю треугольные матрицы, но принудительно отбрасывает (заменяет нулями) те новые элементы (fill-ins), которые возникают в процессе исключения. В результате мы получаем приближенные матрицы L и U, которые сохраняют идеальную разреженность исходной системы. Решение системы с такими треугольными матрицами выполняется мгновенно прямой и обратной подстановкой. Баланс между точностью ILU (задаваемой порогом отбрасывания) и скоростью одной итерации является главным искусством настройки современных коммерческих решателей. Для самых тяжелых задач применяются многосеточные предобуславливатели (AMG), которые обеспечивают теоретическую линейную масштабируемость при решении миллиардных систем.