Прелюдия


Численное решение линейных систем уравнений является незаменимым шагом во многих сферах прикладной математики, инженерии и IT-индустрии, будь то работа с графикой, расчёт аэродинамики какого-нибудь самолёта или оптимизация логистики. Модная нынче «машинка» без этого тоже не особо бы продвинулась. Причём решение линейных систем, как пр��вило, сжирает наибольший процент из всех вычислительных затрат. Понятно, что эта критическая составляющая требует максимальной оптимизации по скорости.

Часто работают с т.н. разреженными матрицами — теми, у которых нулей на порядки больше остальных значений. Такое, например, неизбежно, если имеешь дело с уравнениями в частных производных или с какими-нибудь другими процессами, в которых возникающие элементы в определяющих линейных соотношениях связаны лишь с «соседями». Вот возможный пример разреженной матрицы для известного в классической физике одномерного уравнения Пуассона $-\phi^{''} = f$ на равномерной сетке (да, пока в ней нулей не так много, но при измельчении сетки их будет будь здоров!):

$\begin{pmatrix}2 & -1 & 0 & 0 & 0\\ -1 & 2 & -1 & 0 & 0\\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2\end{pmatrix}$


Противоположные им матрицы — те, в которых на количество нулей не обращают внимания и учитывают все компоненты без исключения, — называют плотными.

Разреженные матрицы часто, из разнообразных соображений, представляют в сжато-столбцовом формате — Compressed Sparse Column (CSC). При таком описании используются два целочисленных массива и один с плавающей точкой. Пусть у матрицы всего $nnz_A$ ненулей и $N$ столбцов. Элементы матрицы перечисляются по столбцам слева направо, без исключений. Первый массив $i_A$ длины $nnz_A$ содержит номера строк ненулевых компонент матрицы. Второй массив $j_A$ длины $N + 1$ содержит… нет, не номера столбцов, потому что тогда матрица была бы записана в координатном формате (coordinate), или троичном (triplet). А содержит второй массив порядковые номера тех компонент матрицы, с которых начинаются столбцы, включая дополнительный фиктивный столбец в конце. Наконец, третий массив $v_A$ длины $nnz_A$ уже содержит сами компоненты, в порядке, соответствующем первому массиву. Вот, например, в предположении, что нумерация строк и столбцов матриц начинается с нуля, для конкретной матрицы

$A = \begin{pmatrix}0 & 1 & 0 & 4 & -1\\ 7 & 0 & 2.1 & 0 & 3\\ 0 & 0 & 0 & 10 & 0\end{pmatrix}$


массивы будут $i_A = \{1, 0, 1, 0, 2, 0, 1\}$, $j_A = \{0, 1, 2, 3, 5, 7\}$, $v_A = \{7, 1, 2.1, 4, 10, -1, 3\}$.

Методы решения линейных систем делятся на два больших класса — прямые и итерационные. Прямые основаны на возможности представить матрицу в виде произведения двух более простых матриц, чтобы затем расщепить решение на два простых этапа. Итерационные используют уникальные свойства линейных пространств и работают на древнем как мир способе последовательного приближения неизвестных к искомому решению, причём в самом процессе сходимости матрицы используются, как правило, только для умножения их на векторы. Итерационные методы гораздо дешевле прямых, но вяло работают на плохо обусловленных матрицах. Когда важна железобетонная надёжность — используют прямые. Вот их я тут и хочу немного затронуть.

Допустим, для квадратной матрицы $A$ нам доступно разложение вида $A = LU$, где $L$ и $U$ — соответственно нижнетреугольная и верхнетреугольная матрицы соотвественно. Первая означает, что выше диагонали у неё одни нули, вторая — что ниже диагонали. Как именно мы получили это разложение — нас сейчас не интересует. Вот простой пример разложения:

$\begin{pmatrix}1 &-1 &-1 \\ 2 &- 1 &-0.5 \\ 4 &-2 &-1.5\end{pmatrix} = \begin{pmatrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1\end{pmatrix} \begin{pmatrix}1 &-1 &-1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5\end{pmatrix}$


Как в этом случае решать систему уравнений $Ax = f$, например, с правой частью $f = \begin{pmatrix}4 \\2 \\3\end{pmatrix}$? Первый этап — прямой ход (forward solve = forward substitution). Обозначаем $y := Ux$ и работаем с системой $Ly = f$. Т.к. $L$ — нижнетреугольная, последовательно в цикле находим все компоненты $y$ сверху вниз:

$\begin{pmatrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1\end{pmatrix} \begin{pmatrix}y_1 \\y_2 \\ y_3\end{pmatrix} = \begin{pmatrix}4 \\ 2 \\ 3\end{pmatrix}\implies y = \begin{pmatrix}4 \\ -6 \\ -1\end{pmatrix}$



Центральная идея состоит в том, что, по нахождении $i$-ой компоненты вектора $y$, она умножается на столбец с тем же номер матрицы $L$, который затем вычитается из правой части. Сама матрица как будто схлопывается слева направо, уменьшаясь в размере по мере нахождения всё больше и больше компонент вектора $y$. Этот процесс называется «уничтожением столбцов» (column еlimination).

Второй этап — обратный ход (backward solve = backward substitution). Имея найденный вектор $y$, решаем $Ux = y$. Здесь мы уже идём по компонентам снизу вверх, но идея остаётся прежней: $i$-ый столбец умножается на только что найденную компоненту $x_i$ и переносится вправо, а матрица схлопывается справа налево. Весь процесс проиллюстрирован для упомянутой в примере матрицы на картинке ниже:

$\small\begin{pmatrix}1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1\end{pmatrix} \begin{pmatrix}y_1 \\y_2 \\ y_3\end{pmatrix} = \begin{pmatrix}4 \\ 2 \\ 3\end{pmatrix}\implies \begin{pmatrix} 1 & 0 \\ 2 & 1\end{pmatrix} \begin{pmatrix}y_2 \\ y_3\end{pmatrix} = \begin{pmatrix} -6 \\ -13\end{pmatrix}\implies\begin{pmatrix}1\end{pmatrix} \begin{pmatrix} y_3\end{pmatrix} = \begin{pmatrix}-1\end{pmatrix}$


$\small\begin{pmatrix}1 &-1 &-1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5\end{pmatrix}\begin{pmatrix}x_1 \\ x_2 \\ x_3\end{pmatrix} = \begin{pmatrix}4 \\ -6 \\ -1\end{pmatrix}\Rightarrow\begin{pmatrix}1 &-1 \\ 0 & 1 \end{pmatrix}\begin{pmatrix}x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix}6 \\ -9\end{pmatrix}\implies\begin{pmatrix}1 \end{pmatrix}\begin{pmatrix}x_1\end{pmatrix} = \begin{pmatrix}-3\end{pmatrix}$


и наше решение будет $x = \begin{pmatrix}-3\\-9\\2\end{pmatrix}$.

Если матрица плотная, то есть полностью представлена в виде какого-то массива, одномерного или двумерного, и доступ к конкретному элементу в ней совершается за время $O(1)$, то подобная процедура решения при уже имеющемся разложении не представляет сложности и кодится легко, поэтому даже не будем тратить на неё время. Что делать, если матрица разреженная? Тут в принципе тоже не сложно. Вот код на С++ для прямого хода, в котором решение $x$ записывается за место правой части, без проверок корректности вводимых данных (массивы CSC соответствуют нижнетреугольной матрице):

Алгоритм 1:

void forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x)
{
    size_t j, p, n = x.size();
    for (j = 0; j < n; j++) // цикл по столбцам
    {
        x[j] /= vA[jA[j]];
        for (p = jA[j]+1; p < jA[j+1]; p++)
            x[iA[p]] -= vA[p] * x[j] ;
    }
}

Всё дальнейшее обсуждение будет касаться только прямого хода для решения нижнетреугольной системы $Lx = f$.

Завязка


А что если правая часть, т.е. вектор справа от знака равенства $Lx = f$, сам имеет большое количество нулей? Тогда имеет смысл пропустить вычисления, связанные с нулевыми позициями. Изменения в коде в этом случае минимальны:

Алгоритм 2:

void fake_sparse_forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x)
{
    size_t j, p, n = x.size();
    for (j = 0; j < n; j++) // цикл по столбцам
    {
        if(x[j])
        {
            x[j] /= vA[jA [j]];
            for (p = jA[j]+1; p < jA[j+1]; p++)
                x[iA[p]] -= vA[p] * x[j];
        }
    }
}

Единственное, что мы добавили — это оператор if, целью которого является сокращение количества арифметических действий до их фактического числа. Если, например, вся правая часть состоит из нулей, то и считать ничего не надо: решение и будет правой частью.

Всё выглядит здорово и конечно же будет работать, но тут, после долгой прелюдии, наконец-то видна проблема — асимптотически низкая производительность данного решателя в случае больших систем. А всё из-за самого факта наличия цикла for. В чём же проблема? Даже если условие внутри if оказывается истинным крайне редко, от самого цикла никуда не деться, и это порождает сложность алгоритма $O(N)$, где $N$ — размер матрицы. Как бы циклы ни оптимизировались современными компиляторами, эта сложность даст о себе знать при больших $N$. Особенно обидно, когда весь вектор $f$ сплошь состоит из нулей, ведь, как мы говорили, и делать-то в этом случае тогда ничего не надо! Ни одной арифметической операции! Какие к чёрту $O(N)$ действий?

Ну хорошо, допустим. Даже если так, почему нельзя просто стерпеть прогон for в холостую, ведь реальных вычислений с вещественными числами, т.е. тех, что попадают под if, всё равно будет мало? А дело в том, что данная процедура прямого хода с разреженной правой частью сама часто используется во внешних циклах и лежит в основе разложения Холец��ого $A = LL^T$ и левостороннего (left-looking) LU-разложения. Да-да, одних из тех разложений, без умения делать которые все эти прямые и обратные ходы в решении линейных систем теряют всякий практический смысл.

Теорема. Если матрица симметричная положительно-определённая (SPD), то её можно представить в виде $A = LL^T$ единственным образом, где $L$ — нижнетреугольная матрица с положительными элементами на диагонали.

Для сильно разреженных SPD матриц используют верхнестороннее (up-looking) разложение Холецкого. Схематически представляя разложение в матрично-блочном виде

$\begin{pmatrix} \mathbf{L}_{11} & \mathbf{0} \\ \mathbf{l}_{12}^T & l_{22} \end{pmatrix}\begin{pmatrix} \mathbf{L}_{11}^T & \mathbf{l}_{12} \\ \mathbf{0} & l_{22} \end{pmatrix} = \begin{pmatrix}\mathbf{A}_{11} & \mathbf{a}_{12} \\ \mathbf{a}_{12}^T & a_{22}\end{pmatrix},$


весь процесс факторизации можно логически разбить всего на три шага.

Алгоритм 3:
  1. верхний метод разложения Холецкого меньшей размерности $\mathbf{L}_{11}\mathbf{L}_{11}^T = \mathbf{A}_{11}$ (рекурсия!)
  2. нижнетреугольная система с разреженной правой частью $\mathbf{L}_{11}\mathbf{l}_{12} = \mathbf{a}_{12}$ (вот оно!)
  3. вычисление $l_{22} = \sqrt{a_{22} - \mathbf{l}_{12}^T\mathbf{l}_{12}}$ (тривиальная операция)

На практике это реализуется так, что в одном большом цикле совершаются шаги 3 и 2, причём именно в этом порядке. Таким образом, осуществляется прогон матрицы по диагонали сверху вниз, увеличивающий матрицу $L$ cтрока за строкой на каждой итерации цикла.

Если в подобном аглоритме сложность прямого хода в шаге 2 будет $O(N)$, где $N$ — размер нижнетреугольной матрицы $\mathbf{L}_{11}$ на произвольной итерации большого цикла, то cложность всего разложения будет как минимум $O(N^2)$! Этого бы ох как не хотелось!

State of the art


Многие алгоритмы так или иначе основаны на мимикрировании человеческих действий при решении проблем. Если дать человеку нижнетреугольную линейную систему с правой частью, в которой только 1-2 ненуля, то он сначала, конечно, пробежит вектор правой части глазами сверху вниз ( тот проклятый цикл сложности $O(N)$! ), чтобы найти эти ненули. Потом он будет испо��ьзовать только их, не тратя время на нулевые компоненты, ведь на решение последние не повлияют: делить нули на диагональные компоненты матрицы нет смысла, равно как и переносить домноженный на ноль столбец вправо. Это и есть показанный выше алгоритм 2. Чудес не бывает. Но что если человеку сразу дать номера ненулевых компонент из каких-то других источников? Например, если правая часть — это столбец какой-то другой матрицы, как это обстоит дело с разложением Холецкого, то там у нас есть мгновенный доступ к его ненулевым компонентам, если их запрашивать последовательно:


// если столбец с индексом j, то пробегаем только ненулевые элементы:
for (size_t p = jA[j]; p < jA[j+1]; p++)
// получаем ненулевой элемент vA[p] матрицы на пересечении строки iA[p] и столбца j

Сложность такого доступа — $O(nnz_j)$, где $nnz_j$ — количество ненулевых компонент в столбце j. Благодарим бога за формат CSC! Как видим, он не только для экономии памяти используется.

В таком случае, мы можем уловить самую суть происходящего при решении методом прямого хода $Lx = f$ для разреженных матрицы $L$ и правой части $f$. Задерживаем дыхание: мы берём ненулевую компоненту $f_i$ в правой части, находим соответствующюю ей переменную $x_i$ делением на $L_{ii}$, а потом, домножив на эту найденную переменную весь i-ый столбец, вводим дополнительные ненули в правой части вычитанием этого столбца из правой части! Данный процесс отлично описывается на языке… графов. Причём ориентированных и нецикличных.

Определим для нижнетреугольной матрицы, у которой отсутствуют нули на диагонали, граф связности. Будем считать, что нумерация строк и столбцов начинается с нуля.
Определение. Графом связности для нижнетреугольной матрицы $L$ размера $N$, у которой отсутствуют нули на диагонали, назовём совокупность множества узлов $V = \{0, ..., N-1\}$ и ориентированных рёбер $E = \{(j, i) | L_{ij} \ne 0\}$.

Вот, например, как выглядит граф связности для нижнетреугольной матрицы

в которой цифры просто соответствуют порядковому номеру диагонального элемента, а точки обозначают ненулевые компоненты ниже диаг��нали:



Определение. Достижимостью ориентированного графа $G$ на множестве индексов $W \subset V$ назовём такое множество индексов $R_G(W) \subset V$, что в любой индекс $z\in R_G(W)$ можно попасть проходом по графу с некоторого индекса $w(z) \in W$.

Пример: для графа с картинки $R_G(\{ 0, 4\}) = \{0, 4, 5, 6\}$. Понятно, что всегда выполняется $W \subset R_G(W)$.

Если представить каждый узел графа как номер столбца породившей его матрицы, то соседние узлы, в которые направлены его рёбра, соответствуют номерам строк ненулевых компонент в этом столбце.

Пусть $nnz(x) = \{j | x_j \ne 0\}$ обозначает множество индексов, соответствующих ненулевым позициям в векторе x.

Теорема Гильберта (нет, не того, чьим именем названы пространства) Множество $nnz(x)$, где $x$ есть вектор решения разреженной нижнетреугольной системы $Lx = f$ с разреженной правой частью $f$, совпадает с $R_G(nnz(f))$.

Добавочка от себя: в теореме мы не учитываем маловероятную возможноть того, что ненулевые числа в правой части, при уничтожении стобцов, могут сократиться в ноль, например, 3 — 3 = 0. Этот эффект называется numerical cancellation. Учитывать такие спонтанно возникающие нули бессмысленная трата времени, и они воспринимаются как и все остальные числа в ненулевых позициях.

Эффективный метод проведения прямого хода с заданной разреженной правой частью, в предположении, что мы имеем прямой доступ к ненулевым её компонентам по индексам, заключается в следующем.

  1. Пробегаем граф методом "поиск в глубину" (depth first search), последовательно стартуя от индекса $i\in nnz(f)$ каждой ненулевой компоненты правой части. Запись найденных узлов в массив $R_G(nnz(f))$ при этом совершаем в том порядке, в котором мы возвращаемся по графу назад! По аналогии с армией захватчиков: оккупируем страну без жестокости, но когда нас стали гнать обратно, мы, разозлённые, уничтожаем всё на своём пути.
    Стоит отметить, что абсолютно не обязательно, чтобы список индексов $nnz(f)$ был отсортирован по возврастанию при подаче его на вход алгоритму «поиск в глубину». Стартовать можно в любом порядке на множестве $nnz(f)$. Различный порядок следования принадлежащих множеству $nnz(f)$ индексов не влияет на итоговое решение, как мы увидим на примере.

    Этот шаг не требует вообще никаких знаний о вещественночисленном массиве $v_A$! Добро пожаловать в мир символьного анализа при работе с прямыми разреженными решателями!
  2. Переходим уже к самому решению, имея в собственном распоряжении массив $R_G(nnz(f))$ с предыдущего шага. Столбцы уничтожаем в порядке, обратном порядку записи этого массива. Вуаля!


Пример


Рассмотрим пример, на котором оба шага детально продемонстрированы. Пусть у нас есть матрица размера 12 x 12 следующего вида:

Соответствующий ей граф связности имеет вид:

Пусть в правой части ненули находятся только на позициях 3 и 5, то есть $nnz(f) = \{3, 5\}$. Пробежим граф, стартуя с этих индексов в написанном порядке. Тогда метод «поиск в глубину» будет выглядеть следующим образом. Сначала мы посетим индексы в порядке $3\to 8\to 11$, не забывая помечать индексы как посещённые. На рисунке ниже они закрашены синим:

При возвращении назад, заносим индексы в наш массив индексов ненулевых компонент решения $nnz(x) := \{\color{blue}{11}, \color{blue}8, \color{blue}3\}$. Далее пробуем пробежать $5\to8\to...$, но натыкаемся на уже помеченный узел 8, поэтому этот маршрут не трогаем и переходим к ветви $5\to9\to...$. Итогом этого пробега будет $5\to9\to10$. Узел 11 посетить не можем, так как он уже меченный. Итак, возвращаемся назад и дополняем массив $nnz(x)$ новым уловом, отмеченным зелёным цветом: $nnz(x) := \{\color{blue}{11}, \color{blue}8, \color{blue}3, \color{green}{10}, \color{green}9, \color{green}5\}$. А вот рисунок:

Мутно-зелёные узлы 8 и 11 — те, которые мы хотели посетить во время второго пробега, но не смогли, т.к. уже посетили во время первого.

Таким образом, массив $nnz(x) = R_G(nnz(f))$ сформирован. Переходим ко второму шагу. А шаг простой: пробегаем массив $nnz(x)$ в обратном порядке (справа налево), находя соответствующие компоненты вектора решения $x$ делением на диагональные компоненты матрицы $L$ и перенося столбцы вправо. Остальные компоненты $x$ как были нулями, так ими и остались. Cхематически это изображено ниже, где нижняя стрелка указывает на порядок уничтожения столбцов:

Обратим внимание: в таком порядке уничтожения столбцов номер 3 встретится позже номеров 5, 9 и 10! Столбцы уничтожаются не в отсортированном по возрастанию порядке, что было бы ошибкой для плотных, т.е. неразреженных матриц. Но для разреженных матриц подобное в порядке вещей. Как видно из ненулевой структуры матрицы в данном примере, использование столбцов 5, 9 и 10 до столбца 3 не исказит в ответ�� компоненты $x_5$, $x_9$, $x_{10}$ нисколько, у них c $x_3$ просто «нет пересечений». Наш метод это учёл! Аналогично использование столбца 8 после столбца 9 не испортит компоненту $x_9$. Прекрасный алгоритм, не правда ли?

Если же мы будем делать обход графа сначала с индекса 5, а потом с индекса 3, то наш массив получится $nnz(x) := \{\color{blue}{11}, \color{blue}8, \color{blue}{10}, \color{blue}{9}, \color{blue}5, \color{green}3\}$. Уничтожение столбцов в порядке, обратном этому массиву, даст абсолютно такое же решение как и в первом случае!

Сложность всех наших операций масштабируется по количеству фактических арифметических операций и количеству ненулевых компонент в правой части, но никак не по размеру матрицы! Мы добились своего.

Критика


ОДНАКО! Критически настроенный читатель заметит, что само предположение в начале, будто у нас есть «прямой доступ к ненулевым компонентам правой части по индексам», уже означает, что когда-то ранее мы пробежали правую часть сверху вниз, чтобы найти эти самые индексы и организовать их в массиве $nnz(f)$, то есть уже затратили $O(N)$ действий! Более того, сам пробег графа требует, чтобы мы предварительно выделили память максимально возможной длины (куда-то же нужно записывать найденные поиском в глубину индексы!), чтобы не мучиться с дальнейшей переаллокацией по мере увеличения массива $nnz(x)$, а это тоже требует $O(N)$ операций! Зачем тогда, мол, весь сыр-бор?

А ведь действительно, для разового решения разреженной нижнетреугольной системы с разреженной же правой частью, изначально заданной в виде плотного вектора, нет смысла тратить время разработчика на все упомянутые выше алгоритмы. Они могут оказаться даже медленнее метода в лоб, представленного алгоритмом 2 выше. Но, как уже было сказано ранее, этот аппарат незаменим при факторизациях Холецкого, поэтому помидорами в меня кидаться не стоит. Действительно, перед запуском алгоритма 3, вся необходимая память максимальной длины выделяется сразу и требует $O(N)$ времени. В дальнейшем цикле по столбцам $A$ все данные лишь перезаписываются в массиве фиксированной длины $N$, причём только в тех позициях, в которых эта пе��езаписть актуальна, благодаря прямому доступу к ненулевым элементам. И именно за счёт этого и возникает эффективность!

Реализация на C++


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

/*
j - стартовый индекс
iA, jA - целочисленные массивы нижнетреугольной матрицы, представленной в формате CSC
top - значение текущей глубины массива nnz(x); в самом начале передаём 0
result - массив длины N, на выходе содержит nnz(x) с индекса 0 до top-1 включительно
marked - массив меток длины N; на вход подаём заполненным нулями
*/
void DepthFirstSearch(size_t j, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t& top, std::vector<size_t>& result, std::vector<int>& marked)
{
    size_t p;
    marked[j] = 1; // помечаем узел j как посещённый
    for (p = jA[j]; p < jA[j+1]; p++) // для каждого ненулевого элемента в столбце j
    {
        if (!marked[iA[p]]) // если iA[p] не помечен
        {
            DepthFirstSearch(iA[p], iA, jA, top, result, marked); // Поиск в глубину на индексе iA[p]
        }
    }
    result[top++] = j; // записываем j в массив nnz(x)
}

Если в самом первом вызове DepthFirstSearch передать переменную top равную нулю, то после завершения всей рекурсивной процедуры переменная top будет равняться количеству найденных индексов в массиве result. Домашнее задание: напишите ещё одну функцию, которая принимает массив индексов ненулевых компонент в правой части и последовательно передаёт их функции DepthFirstSearch. Без этого алогоритм не полный. Обратите внимание: массив вещественных чисел $vA$ мы вообще не передаём в функцию, т.к. он не нужен в процессе символьного анализа.

Несмотря на простоту, дефект реализации очевиден: для больших систем переполнение стека не за горами. Ну тогда есть вариант, основанный на цикле вместо рекурсии. Он посложнее для понимания, но уже годен на все случаи жизни:

/*
j - стартовый индекс
N - размер матрицы
iA, jA - целочисленные массивы нижнетреугольной матрицы, представленной в формате CSC
top - значение текущей глубины массива nnz(x)
result - массив длины N, на выходе содержит часть nnz(x) с индексов top до ret-1 включительно, где ret - возвращаемое значение функции
marked - массив меток длины N
work_stack - вспомогательный рабочий массив длины N
*/
size_t DepthFirstSearch(size_t j, size_t N, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t top, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack)
{
    size_t p, head = N - 1;
    int done;
    result[N - 1] = j; // инициализируем рекурсионный стек
    while (head < N)
    {
        j = result[head]; // получаем j с верхушки рекурсионного стека
        if (!marked[j])
        {
            marked[j] = 1; // помечаем узел j как посещённый
            work_stack[head] = jA[j];
        }
        done = 1; // покончили с узлом j в случае отсутствия непосещённых соседей
        for (p = work_stack[head] + 1; p < jA[j+1]; p++) // исследуем всех соседей j
        {
            // работаем с cоседом с номером iA[p]
            if (marked[iA[p]]) continue; // узел iA[p] посетили раньше, поэтому пропускаем
            work_stack[head] = p; // ставим на паузу поиск по узлу j
            result[--head] = iA[p]; // запускаем поиск в глубину на узле iA[p]
            done = 0; // с узлом j ещё не покончили
            break;
        }
        if (done) // поиск в глубину на узле j закончен
        {
            head++; // удаляем j из рекурсионного стека
            result[top++] = j; // помещаем j в выходной массив
        }
    }
    return (top);
}

А вот так уже выглядит сам генератор ненулевой структуры вектора решения $nnz(x)$:
/*
iA, jA - целочисленные массивы матрицы, представленной в формате CSC
iF - массив индексов ненулевых компонент вектора правой части f
nnzf - количество ненулей в правой части f
result - массив длины N, на выходе содержит nnz(x) с индексов 0 до ret-1 включительно, где ret - возвращаемое значение функции
marked - массив меток длины N, передаём заполненным нулями
work_stack - вспомогательный рабочий массив длины N
*/
size_t NnzX(const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<size_t>& iF,  size_t nnzf, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack)
{
    size_t p, N, top;
    N = jA.size() - 1;
    top = 0;
    for (p = 0; p < nnzf; p++)
        if (!marked[iF[p]]) // начинаем поиск в глубину на непомеченном узле
            top = DepthFirstSearch(iF[p], N, iA, jA, vA, top, result, marked, work_stack);
    for (p = 0; p < top; p++)
        marked[result[p]] = 0; // очищаем метки
    return (top);
}

Наконец, объединяя всё во едино, пишем сам нижнетреугольный решатель для случая разреженной правой части:

/*
iA, jA, vA - массивы матрицы, представленной в формате CSC
iF - массив индексов ненулевых компонент вектора правой части f
nnzf - количество ненулей в правой части f
vF - массив ненулевых компонент вектора правой части f
result - массив длины N, на выходе содержит nnz(x) с индексов 0 до ret-1 включительно, где ret - возвращаемое значение функции
marked - массив меток длины N, передаём заполненным нулями
work_stack - вспомогательный рабочий массив длины N
x - вектор решения длины N, которое получим на выходе; на входе заполнен нулями
*/
size_t lower_triangular_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, const std::vector<size_t>& iF,  size_t nnzf, const std::vector<double>& vF, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack, std::vector<double>& x)
{
    size_t top, p, j;
    ptrdiff_t px;
    top = NnzX(iA, jA, iF, nnzf, result, marked, work_stack);
    for (p = 0; p < nnzf; p++) x[iF[p]] = vF[p]; // заполняем плотный вектор
    for (px = top; px > -1; px--) // прогон в обратном порядке
    {
        j = result[px]; // x(j) будет ненулём
        x [j] /= vA[jA[j]]; // мгновенное нахождение x(j)
        for (p = jA[j]+1; p < jA[j+1]; p++)
        {
            x[iA[p]] -= vA[p]*x[j]; // уничтожение j-ого столбца
        }
    }
    return (top) ;
}

Видим, что цикл наш пробегает только по индексам массива $nnz(x)$, а не по всему набору $0, 1, ..., N - 1$. Готово!

Cуществует реализация, которая не использует массив меток marked с целью экономии памяти. Вместо этого, используется дополнительное множество индексов $V_1$, не пересекающееся с множеством $V = \{0, ..., N - 1\}$, в которое происходит взаимно-однозначное о��ображение простой алгебраической операцией в качестве процедуры отмечивания узла. Однако, в нашу эпоху дешевизны памяти, экономить её на одном массиве длины $N$ кажется абсолютно излишним.

В качестве заключения


Процесс решения разреженной системы линейных уравнений прямым методом, как правило, разбивается на три этапа:

  1. Символьный анализ
  2. Численная факторизация на основе данных сивольного анализа
  3. Решение полученных треугольных систем с плотной правой частью

Второй шаг — численная факторизация — наиболее ресурсозатратная часть и пожирает большую часть (>90%) расчётного времени. Цель первого шага — понизить дороговизну второго. Пример символьного анализа и был представлен в данном посте. Однако именно первый шаг требует самого длительного времени разработки и максимальных умственных затрат со стороны разработчика. Хороший символьный анализ требует знания теории графов и деревьев и владения «алгоритмическим чутьём». Второй шаг несоизмеримо проще в реализации.
Ну а третий шаг и по реализации, и по расчётному времени в большинстве случаев самый непритязательный.

Хорошее введение в прямые методы для разреженных систем есть в книге профессора факультета компьютерных наук университета Texas A&M Тима Девиса под названием "Direct Methods for Sparse Linear Systems". Там и описан показанный выше алгоритм. С самим Девисом я не знаком, хотя в своё время сам работал в том же университете (правда, на другом факультете). Если не ошибаюсь, Девис сам участвовал в разработке решателей, используемых в Матлабе(!). Более того, он причастен даже к генераторам картинок в Google Street View (метод наименьших квадратов). Кстати, на том же факультете числится никто иной как сам Бьёрн Страуструп — создатель языка C++.

Возможно, когда-нибудь напишу ещё что-нибудь на тему разреженных матриц или численных методов вообще. Всем добра!