Алгоритмы на графах для разреженных матриц: Катхилл-Макки и минимальная степень

Скрытая угроза прямых решателей: проблема заполнения

Мы неоднократно обращали внимание на то, что огромные системы уравнений (СЛАУ), возникающие в инженерных расчетах, являются сверхразреженными: в них присутствуют миллионы нулей и лишь малая доля значащих чисел. Для их хранения созданы эффективные форматы вроде CSR. Однако если мы попытаемся решить такую систему точным, прямым методом (LU-разложение или разложение Холецкого для симметричных матриц), мы столкнемся с катастрофическим явлением, которое в вычислительной алгебре называется «заполнением» (fill-in).

В процессе прямого исключения переменных (вспомните метод Гаусса), вычитая одну строку из другой, алгоритм неизбежно складывает ненулевые элементы. Если в исходной строке стоял спасительный ноль, то в процессе прямого хода на его месте может появиться значащее число. К концу алгоритма матрица, которая изначально была почти пустой, рискует превратиться в гигантскую, абсолютно плотную структуру. Миллионы нулей исчезнут, алгоритм сожрет всю оперативную память кластера (Out of Memory) и процесс остановится. Как же предотвратить это катастрофическое разрастание данных?

Матрица как топологический граф

Гениальное решение пришло на стыке алгебры и дискретной математики. Любую симметричную матрицу можно представить в виде неориентированного графа. Каждая строка (уравнение) — это узел графа (вершина). Если матричный элемент на пересечении строки A и столбца B не равен нулю, значит между узлом A и узлом B в графе существует физическое ребро связи. Заполнение матрицы в процессе LU-разложения в терминах теории графов означает добавление новых, фиктивных ребер между узлами.

Порядок, в котором мы перечисляем уравнения (нумеруем переменные x1, x2, x3), не меняет физического решения системы, но он критически, на порядки, влияет на уровень заполнения матрицы нулями в памяти! Следовательно, перед тем как запустить вычислительно тяжелое LU-разложение, нам нужно запустить быстрый алгоритм переупорядочивания (эвристику на графе), который найдет оптимальную нумерацию строк и столбцов матрицы для минимизации будущего заполнения.

Ширина ленты и алгоритмы Reverse Cuthill-McKee (RCM)

Первым классом таких алгоритмов стали методы уменьшения ширины ленты (Bandwidth reduction). Исторически самым известным из них является алгоритм Катхилла-Макки (Cuthill-McKee), а точнее его улучшенная версия — обратный алгоритм Катхилла-Макки (RCM). Алгоритм начинает поиск с периферийного узла графа (узла с минимальным количеством связей) и послойно, методом поиска в ширину (BFS), перенумеровывает все связанные с ним узлы. В результате все ненулевые элементы матрицы «сбегаются» и плотно прижимаются к главной диагонали, образуя узкую ленту. Вне этой ленты остаются чистые нули, которые гарантированно никогда не заполнятся при разложении.

Для более сложных и неструктурированных 3D-сеток ленточные алгоритмы уступают место семейству методов минимальной степени (Minimum Degree Algorithm — AMD). Алгоритм AMD работает на опережение: на каждом шаге он симулирует процесс исключения Гаусса и всегда выбирает для исключения ту переменную, узел которой в графе в данный момент имеет наименьшее число соседей (минимальную степень). Это локально жадная эвристика, но она минимизирует создание новых ребер на каждом микрошаге. Сегодня мощнейшие программные комплексы на графах (такие как METIS или Scotch) выполняют глубокое иерархическое разбиение графа матрицы на кластеры, что является обязательным «нулевым» шагом перед любым суперкомпьютерным матричным расчетом.