Разреженные матрицы: структуры данных и алгоритмы упаковки

В большинстве реальных инженерных задач, таких как расчет упругости методом конечных элементов (МКЭ), анализ электрических цепей или обработка социальных графов, матрицы систем линейных уравнений оказываются разреженными. Это означает, что подавляющее большинство их элементов (иногда до 99.9%) строго равны нулю. Хранить такие нули в памяти компьютера в виде классического двумерного массива — колоссальная и бессмысленная трата ресурсов. Для матриц размерностью в миллионы строк требуются сотни гигабайт оперативной памяти. Решением этой проблемы стала разработка специализированных форматов упаковки разреженных матриц и адаптированных алгоритмов линейной алгебры, которые выполняют вычисления, игнорируя нулевые элементы и избегая явления заполнения (fill-in).

Форматы COO, CSR и CSC: принципы сжатия

Самым примитивным способом хранения разреженных данных является координатный формат (COO — Coordinate list). В нем матрица хранится в виде трех одномерных массивов: массива значений ненулевых элементов, массива их номеров строк и массива номеров столбцов. Хотя этот формат удобен для первичного конструирования матрицы, он неэффективен для арифметических операций. Стандартом в вычислительной линейной алгебре стали форматы сжатого хранения строками (CSR — Compressed Sparse Row) и столбцами (CSC — Compressed Sparse Column). В формате CSR вместо хранения номера строки для каждого элемента хранится массив указателей на начало каждой строки (row pointers). Это позволяет радикально сократить объем памяти (убрав дублирование индексов строк) и обеспечивает феноменально быстрое умножение разреженной матрицы на плотный вектор.

Умножение SpMV: ядро высокопроизводительных вычислений

Операция умножения разреженной матрицы на плотный вектор (Sparse Matrix-Vector multiplication, SpMV) является вычислительным сердцем практически всех итерационных методов линейной алгебры, таких как метод сопряженных градиентов или GMRES. В формате CSR алгоритм SpMV перебирает только ненулевые элементы, выполняя ровно столько умножений, сколько реальных чисел хранится в памяти. Однако для современных процессоров и видеокарт (GPU) классический SpMV создает проблемы с доступом к памяти: из-за нерегулярной структуры матрицы процессору приходится прыгать по случайным адресам вектора-сомножителя, что приводит к частым промахам кэша (cache misses). Для преодоления этого ограничения инженеры разрабатывают блочные форматы хранения (BSR), которые группируют ненулевые элементы в крошечные плотные блоки, оптимизируя векторизацию (SIMD) и доступ к памяти.

Проблема заполнения (Fill-in) при прямых методах

Когда к разреженной матрице применяется классическое LU-разложение или разложение Холецкого, происходит катастрофическое явление: на месте нулей образуются ненулевые элементы. Это называется заполнением (fill-in). Матрица коэффициентов, которая в исходном виде занимала 10 Мегабайт, в процессе решения может раздуться до 10 Гигабайт, исчерпав всю память системы. Величина заполнения критически зависит от порядка, в котором мы исключаем неизвестные (от нумерации строк и столбцов). Таким образом, перед решением системы матрицу необходимо переупорядочить.

Эвристики переупорядочивания: Cuthill-McKee и Minimum Degree

Задача поиска идеальной перестановки строк и столбцов, минимизирующей заполнение матрицы, является NP-трудной, поэтому на практике применяются жадные эвристические алгоритмы. Одним из самых известных является алгоритм Катхилла-Макки (Cuthill-McKee) и его обратная версия (RCM). Этот алгоритм работает как обход графа в ширину: он собирает ненулевые элементы как можно ближе к главной диагонали, уменьшая «ширину ленты» матрицы, что идеально подходит для профильных решателей. Другой подход — алгоритм минимальной степени (Minimum Degree Ordering). На каждом шаге гипотетического исключения Гаусса он выбирает ту строку (вершину графа), у которой меньше всего ненулевых соседей. Этот метод не сужает ленту, но в целом порождает гораздо меньше мусорных элементов, являясь золотым стандартом в библиотеках прямых разреженных решателей (таких как UMFPACK или PARDISO).

Оценить