Методы Ланцоша и Арнольди: поиск собственных значений гигантских матриц

Охота за экстремальными резонансами

В предыдущих материалах мы исследовали проблему поиска всех собственных значений матрицы с помощью QR-алгоритма. QR-разложение блестяще работает для плотных матриц умеренной размерности (до нескольких тысяч строк). Однако в реальных физических задачах инженеры сталкиваются со сверхбольшими разреженными матрицами, содержащими миллионы уравнений. Например, при анализе устойчивости плазмы в термоядерном реакторе, при расчете резонансных частот гигантского моста или в алгоритмах поиска кластеров в графах социальных сетей. Применение прямого QR-алгоритма к матрице миллион на миллион физически невозможно — оно потребует экзабайты оперативной памяти и столетия вычислений.

К счастью, в таких задачах инженеру почти никогда не нужны абсолютно все миллионы частот. Обычно критически важно знать лишь несколько самых малых собственных значений (они определяют фундаментальную частоту колебаний моста) или несколько самых больших (они определяют алгоритмическую устойчивость системы). Для поиска этих экстремальных собственных значений у сверхбольших разреженных матриц были созданы мощнейшие итерационные алгоритмы: метод Ланцоша (для симметричных матриц) и метод Арнольди (для несимметричных).

Подпространство Крылова и ортогональный базис

Оба метода базируются на элегантной концепции подпространств Крылова, которую мы упоминали в контексте решения СЛАУ. Алгоритм выбирает случайный начальный вектор q1 и начинает многократно умножать его на гигантскую матрицу A. В результате формируется последовательность векторов: q1, A*q1, A^2*q1, A^3*q1. Эти векторы несут в себе колоссальную математическую информацию о доминирующих собственных векторах матрицы (так как при многократном умножении вектор вытягивается вдоль главного собственного направления).

Проблема в том, что эти векторы очень быстро становятся почти коллинеарными (параллельными), что уничтожает их вычислительную полезность (матрица их скалярных произведений становится вырожденной). Гениальность метода Арнольди (Уолтер Арнольди, 1951) заключается в применении процесса ортогонализации Грама-Шмидта прямо на лету! Каждый новый сгенерированный вектор тут же жестко проецируется и очищается от компонентов всех предыдущих векторов. В результате мы получаем идеальный ортонормированный базис подпространства Крылова. В этом новом, маленьком базисе исходная гигантская матрица проецируется (сжимается) в крошечную матрицу Хессенберга размера m x m (где m — число итераций, например, 50 или 100).

Метод Ланцоша и проблема потери ортогональности

Если исходная матрица системы была симметричной (как в задачах сопротивления материалов), процесс Арнольди превращается в Метод Ланцоша (Корнелий Ланцош, 1950). Из-за свойств симметрии матрица Хессенберга магическим образом схлопывается в симметричную трехдиагональную матрицу. Это означает, что для генерации каждого нового ортогонального вектора алгоритму Ланцоша не нужно вспоминать и вычитать все предыдущие 100 векторов — ему достаточно использовать короткую рекуррентную формулу, опирающуюся только на два последних вектора! Это делает метод Ланцоша феноменально быстрым и экономным по памяти.

Получив эту крошечную матрицу проекции, мы просто решаем ее стандартным QR-алгоритмом за миллисекунды. Ее собственные значения (называемые числами Ритца) с фантастической скоростью сходятся к истинным экстремальным собственным значениям гигантской исходной матрицы. Однако у метода Ланцоша есть вычислительный недуг: из-за конечной разрядности компьютеров (ошибок округления float) генерируемые векторы постепенно теряют свою взаимную перпендикулярность. Это порождает появление фиктивных (призрачных) собственных значений, которые являются математической галлюцинацией. Для борьбы с этим современная вычислительная алгебра использует методы полной или выборочной переортогонализации, которые встроены в ядро таких сверхмощных индустриальных пакетов, как ARPACK (используемых функциями `eigs` в MATLAB и SciPy).