Разреженные матрицы: форматы хранения и SpMV

1. Матрицы, состоящие из пустоты

Матрица называется разреженной (sparse), когда ненулевых элементов в ней настолько мало, что это выгодно использовать, — на практике доли процента. Это не экзотика, а нормальное состояние больших матриц реального мира:

• Уравнения физики на сетках. В уравнении теплопроводности температура узла связана только с четырьмя соседями: в строке матрицы 576×576 — максимум 5 ненулевых (0.84%, посчитано для картинки ниже), и доля падает с ростом сетки.
• Графы. Матрица смежности соцсети: у пользователя сотни друзей из миллиарда строк — миллиардные доли заполнения.
• Рекомендации. Таблица «пользователь × товар»: каждый купил ничтожную долю ассортимента.

У каждой задачи — свой «портрет»: сетки дают аккуратные диагонали, графы — россыпь. Портрет определяет и выбор формата хранения, и поведение алгоритмов. Вот, для сравнения, матрица настоящего инженерного расчёта:

Теперь арифметика масштаба — числа измерены, скрипт в папке сode. Матрица 10 000 × 10 000 double при 0.1% ненулевых (сто тысяч чисел из ста миллионов клеток):

Выигрыш по времени — прямое следствие выигрыша по памяти: умножение на вектор упирается в пропускную способность памяти (урок C06), а разреженный формат просто не читает нули. Причём 800 МБ — это ещё скромно: матрица смежности соцсети в плотном виде не поместилась бы ни в одну машину мира.

2. Форматы: COO и CSR на пальцах

Идея всех форматов одна: хранить список ненулевых, а не таблицу всех клеток. Разберём два главных на одной матрице:

        столбцы:  0  1  2  3
              ┌               ┐
    строка 0  │  5  0  0  1   │
    строка 1  │  0  8  0  0   │
    строка 2  │  0  0  3  0   │
    строка 3  │  0  6  0  4   │
              └               ┘        ненулевых: 6 из 16

COO (coordinate) — самый наивный: три массива длины nnz («number of nonzeros»), по тройке (строка, столбец, значение) на элемент:

row = [0, 0, 1, 2, 3, 3]
col = [0, 3, 1, 2, 1, 3]
val = [5, 1, 8, 3, 6, 4]

COO удобен для сборки матрицы (дописывай тройки в любом порядке), но плох для вычислений: элементы одной строки нужно ещё найти. CSR (compressed sparse row) чинит это, сортируя элементы по строкам и сжимая массив row: вместо номера строки у каждого элемента храним, где в data начинается каждая строка:

data    = [5, 1, 8, 3, 6, 4]      # значения, строка за строкой
indices = [0, 3, 1, 2, 1, 3]      # столбец каждого значения
indptr  = [0, 2, 3, 4, 6]         # границы строк: строка r живёт
                                  # в data[indptr[r] : indptr[r+1]]

Проверьте себя по indptr: строка 0 — элементы data[0:2] = {5, 1} (столбцы 0 и 3); строка 1 — data[2:3] = {8}; строка 2 — data[3:4] = {3}; строка 3 — data[4:6] = {6, 4}. Длина indptr — всегда «строк + 1», а последний элемент равен nnz. Эти массивы сверены со scipy.sparse.csr_matrix — она строит в точности их.

3. SpMV: умножение на вектор и его дисбаланс

SpMV (sparse matrix–vector multiplication, y = A·x) — рабочая лошадь разреженного мира: итерационные решатели, PageRank, графовые нейросети внутри крутят именно её. По CSR она пишется в четыре строки — и заметьте, насколько это близко к скалярному произведению из C05, только «своя» часть данных у каждой строки своя:

for (int row = 0; row < n; row++) {
    float sum = 0.0f;
    for (int k = indptr[row]; k < indptr[row + 1]; k++)
        sum += data[k] * x[indices[k]];       // только ненулевые!
    y[row] = sum;
}

Строки независимы — это map по строкам (C04), и первое CUDA-ядро очевидно: строка на тред:

__global__ void spmv_csr(int n, const int* indptr, const int* indices,
                         const float* data, const float* x, float* y)
{
    int row = blockIdx.x * blockDim.x + threadIdx.x;
    if (row < n) {
        float sum = 0.0f;
        for (int k = indptr[row]; k < indptr[row + 1]; k++)
            sum += data[k] * x[indices[k]];
        y[row] = sum;
    }
}

Ядро верное, но у него врождённая болезнь — дисбаланс нагрузки. Длина внутреннего цикла у каждого треда своя: строка-«звезда» соцсети с миллионом связей и строка-одиночка с двумя попадают в один варп, и 31 тред ждёт, пока тридцать второй дожуёт миллион слагаемых, — SIMT из C02 во всей красе, тот же эффект, что у границы множества Мандельброта в C08, только злее: там разброс был в сотни раз, здесь бывает в миллионы. Второй удар — доступ x[indices[k]]: индексы скачут по вектору произвольно, кэш почти не помогает; вот почему SpMV даже у профессионалов выжимает лишь малую долю пиковой пропускной способности.

Стандартное лечение дисбаланса — варп на строку: все 32 треда варпа совместно жуют одну строку (каждый берёт элементы с шагом 32), а потом складывают частичные суммы редукцией — тем самым деревом из C05, только внутри варпа. Длинная строка перестаёт держать соседей: варп занят ею целиком, а не одним тредом. Для матриц с ровными короткими строками (сетки!) хорош формат ELL с ядром «тред на строку»; универсальные библиотеки комбинируют подходы автоматически.

Что подводит к главному практическому совету: в продакшене SpMV не пишут руками — берут cuSPARSE, библиотеку NVIDIA из состава CUDA Toolkit (там же живут cuBLAS для плотных матриц из C07 и Thrust для сортировок из C08). Она выбирает алгоритм под портрет матрицы и обгоняет самописные ядра почти всегда. Своё ядро выше — учебное: оно нужно, чтобы понимать, почему у cuSPARSE такие ручки настройки и откуда берутся её ограничения.

4. Python: scipy.sparse за пять минут

В Python разреженный мир давно собран в scipy.sparse — и все числа раздела 1 измерены именно им (скрипт c09_sparse_bench.py в папке сode):

import numpy as np
import scipy.sparse as sp

# случайная разреженная матрица: 10^8 клеток, 10^5 ненулевых
A = sp.random(10_000, 10_000, density=0.001, format="csr")

x = np.random.rand(10_000)
y = A @ x                          # SpMV: 0.076 мс против 16.4 мс у плотной

# те же три массива, что в разделе 2:
A.data, A.indices, A.indptr

Рабочие правила те же, что в C++, только именованные по-питоновски: собирать — в COO/LIL, считать — в CSR/CSC, конвертация между форматами дешёвая (.tocsr(), .tocoo()). Плотную матрицу из разреженной делает .toarray() — и это обычно последняя строка перед MemoryError: вызывайте её только на маленьких кусках. Связка с Numba из линии P прямая: @njit-функции охотно принимают тройку data/indices/indptr как обычные numpy-массивы — последовательный SpMV из раздела 3, переписанный на Python, в нашей проверке совпал со scipy до последнего знака.

5. Тренажёр: CSR на лету

Контрольные вопросы

• Восстановите по CSR-массивам data=[5,1,8,3,6,4], indices=[0,3,1,2,1,3], indptr=[0,2,3,4,6] содержимое строки 3.
  Строка 3 живёт в data[indptr[3]:indptr[4]] = data[4:6] = {6, 4}, столбцы — indices[4:6] = {1, 3}. То есть в строке 3 стоят 6 в столбце 1 и 4 в столбце 3, остальное — нули.
• Почему CSR-умножение быстрее плотного в сотни раз, если процессор тот же?
  Умножение на вектор упирается в память, а не в арифметику (C06). Плотное умножение читает все 10⁸ клеток, включая 99.9% нулей; CSR читает только 10⁵ ненулевых с их индексами. Меньше байтов — пропорционально меньше времени: ×215 в нашем измерении.
• Почему матрицу собирают в COO, а считают в CSR?
  В COO тройки дописываются в конец в любом порядке — сборка дешёвая. В CSR данные сжаты и отсортированы по строкам: вычисления по строкам быстрые, но вставка элемента в середину требует сдвига всех трёх массивов. Отсюда протокол: накопить в COO/LIL, один раз конвертировать в CSR.
• Откуда в ядре «строка на тред» дисбаланс и как его лечит схема «варп на строку»?
  Длина цикла треда равна числу ненулевых его строки; в варпе SIMT все ждут самого длинного — строка-«звезда» останавливает 31 соседа. В схеме «варп на строку» все 32 треда обрабатывают одну строку с шагом 32 и складывают частичные суммы редукцией: длинная строка занимает варп целиком, а не один тред.
• Почему для SpMV «10% пиковой пропускной способности» — не приговор и с чем сравнивать?
  Доступ x[indices[k]] нерегулярен: обращения к памяти не сливаются, кэш-линии подтягиваются ради одного числа — физика задачи, а не ошибка кода. Честная линейка — cuSPARSE на той же матрице: если своё учебное ядро в разы медленнее библиотеки, вот тогда есть что чинить.
• Когда формат ELL выгоден, а когда губителен?
  ELL хранит прямоугольник «строки × максимум ненулевых в строке» — доступ регулярный, GPU доволен. Выгоден при почти одинаковой длине строк (сеточные матрицы). Губителен при рваных строках: одна строка-звезда с миллионом элементов раздует паддинг всех остальных строк до миллиона — памяти уйдёт больше, чем у плотной.

Источники

1. Bell N., Garland M. Efficient Sparse Matrix-Vector Multiplication on CUDA : NVIDIA Technical Report NVR-2008-004. — NVIDIA Corporation, 2008. — URL: https://research.nvidia.com/publication/2008-12_efficient-sparse-matrix-vector-multiplication-cuda (дата обращения: 09.07.2026).
2. cuSPARSE Library // NVIDIA Docs : [сайт]. — URL: https://docs.nvidia.com/cuda/cusparse/ (дата обращения: 09.07.2026).
3. Sparse arrays (scipy.sparse) // SciPy Documentation : [сайт]. — URL: https://docs.scipy.org/doc/scipy/reference/sparse.html (дата обращения: 09.07.2026).
4. Sparse matrix // Wikipedia : [сайт]. — URL: https://en.wikipedia.org/wiki/Sparse_matrix (дата обращения: 09.07.2026).