Обработка и решение разреженных задач включает различные модули, которые кратко описаны ниже:

#include <Eigen/SparseCore>

#include <Eigen/SparseCholesky>

#include<Eigen/SparseLU> 

#include<Eigen/SparseQR>

#include <Eigen/IterativeLinearSolvers>

#include <Eigen/Sparse>

Формат разреженной матрицы

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

Класс SparseMatrix

Класс SparseMatrix — основное представление разреженных матриц в модуле разреженных матриц Eigen; он обеспечивает высокую производительность и низкое потребление памяти. Он реализует более универсальный вариант широко используемой схемы хранения по столбцам (или строкам). Он состоит из четырёх компактных массивов:

• Values: хранит значения коэффициентов, отличных от нуля.
• InnerIndices: хранит индексы строк (соответственно, столбцов) ненулевых элементов.
• OuterStarts: хранит для каждого столбца (соответственно, строки) индекс первого ненулевого элемента в предыдущих двух массивах.
• InnerNNZs: хранит количество ненулевых элементов в каждом столбце (соответственно, строке). Слово inner относится к внутреннему вектору, который является столбцом для матрицы, хранящейся по столбцам, или строкой для матрицы, хранящейся по строкам. Слово outer относится к другому направлению.

Эта схема хранения лучше объясняется на примере. Следующая матрица

и одно из возможных разреженных представлений, по столбцам:

В настоящее время элементы данного внутреннего вектора гарантированно всегда упорядочены по возрастанию внутренних индексов. "_" указывает доступное свободное пространство для быстрого вставки новых элементов. Предполагая, что перераспределение не требуется, вставка произвольного элемента, следовательно, занимает O(nnz_j), где nnz_j — количество ненулевых элементов соответствующего внутреннего вектора. С другой стороны, вставка элементов с возрастающими внутренними индексами в данный внутренний вектор намного эффективнее, поскольку это требует только увеличения соответствующей InnerNNZs записи, что является операцией O(1).

Случай, когда свободного места нет, является особым случаем и называется режимом сжатия. Он соответствует широко используемым схемам хранения по столбцам (или строкам) (CCS или CRS). Любую SparseMatrix можно привести к этому виду, вызвав функцию SparseMatrix::makeCompressed(). В этом случае можно отметить, что массив InnerNNZs избыточен с массивом OuterStarts, так как выполняется равенство: InnerNNZs[j] = OuterStarts[j+1]-OuterStarts[j]. Поэтому на практике вызов SparseMatrix::makeCompressed() освобождает этот буфер.

Стоит отметить, что большинство наших обёрток для внешних библиотек требуют сжатые матрицы в качестве входных данных.

Результаты операций Eigen всегда создают сжатые разреженные матрицы. С другой стороны, вставка нового элемента в SparseMatrix преобразует её в режим нескомпрессированный.

Вот предыдущая матрица, представленная в сжатом виде:

SparseVector — это частный случай SparseMatrix, где хранятся только массивы Values и InnerIndices. Для SparseVector нет понятия сжатого/несжатого режима.

Первый пример

Прежде чем описывать каждый отдельный класс, давайте начнём с типичного примера: решение уравнения Лапласа \( \Delta u = 0 \) на регулярной 2D-сетке с помощью схемы конечных разностей и граничных условий Дирихле. Такая задача математически может быть выражена как линейная задача вида \( Ax=b \), где \( x \) — вектор m неизвестных (в нашем случае, значения пикселей), \( b \) — вектор правой части, полученный из граничных условий, и \( A \) — \( m \times m \) матрица, содержащая лишь несколько ненулевых элементов, полученных из дискретизации оператора Лапласа.

#include <Eigen/Sparse> #include <vector> #include <iostream> typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double typedef Eigen::Triplet<double> T; void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n); void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename); int main(int argc, char** argv) { if(argc!=2) { std::cerr << "Error: expected one and only one argument.\n"; return -1; } int n = 300; // size of the image int m = n*n; // number of unknowns (=number of pixels) // Assembly: std::vector<T> coefficients; // list of non-zeros coefficients Eigen::VectorXd b(m); // the right hand side-vector resulting from the constraints buildProblem(coefficients, b, n); SpMat A(m,m); A.setFromTriplets(coefficients.begin(), coefficients.end()); // Solving: Eigen::SimplicialCholesky<SpMat> chol(A); // performs a Cholesky factorization of A Eigen::VectorXd x = chol.solve(b); // use the factorization to solve for the given right hand side // Export the result to a file: saveAsBitmap(x, n, argv[1]); return 0; }

В этом примере мы начинаем с определения типа разреженной матрицы по столбцам с типом данных double SparseMatrix<double>, и списком троек того же типа данных Triplet<double>. Тройка — это простая структура, представляющая ненулевой элемент как тройку: индекс row, индекс column, value.

В основной функции мы объявляем список coefficients троек (как std::vector) и вектор правой части \( b \), которые заполняются функцией buildProblem. Затем неотсортированный и плоский список ненулевых элементов преобразуется в настоящий объект SparseMatrix A. Обратите внимание, что элементы списка не обязательно должны быть отсортированы, и возможные повторяющиеся элементы будут суммированы.

Последний шаг заключается в фактическом решении составленной задачи. Поскольку полученная матрица A является симметричной по построению, мы можем выполнить прямое разложение Холецкого с помощью класса SimplicialLDLT, который ведёт себя как его аналог LDLT для плотных объектов.

Полученный вектор x содержит значения пикселей как одномерный массив, который сохраняется в файл jpeg, показанный справа от кода выше.

Описание функций buildProblem и save выходит за рамки этого учебника. Они приведены здесь для любознательных и целей воспроизводимости.

Класс SparseMatrix

Свойства матрицы и вектора
Классы SparseMatrix и SparseVector принимают три шаблонных аргумента: тип скаляра (например, double), порядок хранения (ColMajor или RowMajor, по умолчанию ColMajor) и тип индекса внутреннего вектора (по умолчанию int).

Как и для плотных объектов Matrix, конструкторы принимают размер объекта. Вот несколько примеров:

SparseMatrix<std::complex<float> > mat(1000,2000);         // declares a 1000x2000 column-major compressed sparse matrix of complex<float>
SparseMatrix<double,RowMajor> mat(1000,2000);              // declares a 1000x2000 row-major compressed sparse matrix of double
SparseVector<std::complex<float> > vec(1000);              // declares a column sparse vector of complex<float> of size 1000
SparseVector<double,RowMajor> vec(1000);                   // declares a row sparse vector of double of size 1000

В остальной части учебника mat и vec представляют собой любые разреженные матрицы и разреженные векторы, соответственно.

Размеры матрицы можно запросить с помощью следующих функций:

mat.rows()
mat.cols()

vec.size() 

mat.innerSize()
mat.outerSize()

mat.nonZeros() 

vec.nonZeros() 

Итерация по ненулевым коэффициентам
Случайный доступ к элементам разреженного объекта можно получить через функцию coeffRef(i,j). Однако эта функция включает достаточно дорогостоящий бинарный поиск. В большинстве случаев необходимо только итерировать по элементам, отличным от нуля. Это достигается стандартным циклом по внешнему измерению, а затем итерацией по ненулевым элементам текущего внутреннего вектора с помощью InnerIterator. Таким образом, ненулевые элементы должны посещаться в том же порядке, что и порядок хранения. Вот пример:

SparseMatrix<double> mat(rows,cols);
for (int k=0; k<mat.outerSize(); ++k)
  for (SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
  {
    it.value();
    it.row();   // row index
    it.col();   // col index (here it is equal to k)
    it.index(); // inner index, here it is equal to it.row()
  }

SparseVector<double> vec(size);
for (SparseVector<double>::InnerIterator it(vec); it; ++it)
{
  it.value(); // == vec[ it.index() ]
  it.index();
}

Для записываемого выражения значение, на которое ссылаются, можно изменить с помощью функции valueRef(). Если тип разреженной матрицы или вектора зависит от параметра шаблона, то необходимо использовать ключевое слово typename, чтобы указать, что InnerIterator обозначает тип; см. Ключевые слова шаблона и typename в C++ для получения дополнительных сведений.

Заполнение разреженной матрицы

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

Простейший способ создания разреженной матрицы с гарантией хорошей производительности — сначала создать список так называемых троек, а затем преобразовать его в SparseMatrix.

Вот типичный пример использования:

typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;
tripletList.reserve(estimation_of_entries);
for(...)
{
  // ...
  tripletList.push_back(T(i,j,v_ij));
}
SparseMatrixType mat(rows,cols);
mat.setFromTriplets(tripletList.begin(), tripletList.end());
// mat is ready to go!

Список троек может содержать элементы в произвольном порядке и даже дублированные элементы, которые будут суммированы функцией setFromTriplets(). См. функцию SparseMatrix::setFromTriplets() и класс Triplet для получения дополнительных сведений.

В некоторых случаях немного лучшая производительность и меньшее потребление памяти могут быть достигнуты путем непосредственного вставки ненулевых элементов в целевую матрицу. Типичный сценарий такого подхода показан ниже:

1: SparseMatrix<double> mat(rows,cols);         // default is column major
2: mat.reserve(VectorXi::Constant(cols,6));
3: for each i,j such that v_ij != 0
4:   mat.insert(i,j) = v_ij;                    // alternative: mat.coeffRef(i,j) += v_ij;
5: mat.makeCompressed();                        // optional

• Ключевым элементом здесь является строка 2, где мы резервируем место для 6 ненулевых элементов на столбец. Во многих случаях количество ненулевых элементов на столбец или строку можно легко узнать заранее. Если оно сильно меняется для каждого внутреннего вектора, то можно указать размер резерва для каждого внутреннего вектора, предоставив объект вектора с оператором [](int j), возвращающим размер резерва для j-th внутреннего вектора (например, через VectorXi или std::vector<int>). Если можно получить только приблизительную оценку количества ненулевых элементов на внутренний вектор, то очень рекомендуется переоценить ее, а не наоборот. Если эта строка опущена, то при первой вставке нового элемента будет зарезервировано место для 2 элементов на внутренний вектор.
• Строка 4 выполняет сортированную вставку. В этом примере идеальный случай, когда столбец j-th не заполнен и содержит ненулевые элементы, внутренние индексы которых меньше i. В этом случае данная операция сводится к тривиальной операции O(1).
• При вызове insert(i,j) элемент i ,j не должен уже существовать, в противном случае используйте метод coeffRef(i,j), который позволит, например, накапливать значения. Этот метод сначала выполняет бинарный поиск и, наконец, вызывает insert(i,j), если элемент уже не существует. Он более гибкий, но также и более дорогостоящий, чем insert().
• Строка 5 подавляет оставшееся пустое пространство и преобразует матрицу в хранилище сжатых столбцов.

Поддерживаемые операторы и функции

Из-за своего специального формата хранения разреженные матрицы не могут предложить такой же уровень гибкости, как плотные матрицы. В модуле разреженных матриц Eigen мы выбрали, чтобы предоставить только подмножество API плотных матриц, которое может быть эффективно реализовано. В дальнейшем sm обозначает разреженную матрицу, sv — разреженный вектор, dm — плотную матрицу, а dv — плотный вектор.

Основные операции

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

sm1.real()   sm1.imag()   -sm1                    0.5*sm1
sm1+sm2      sm1-sm2      sm1.cwiseProduct(sm2)

Однако, существует сильное ограничение: порядки хранения должны совпадать. Например, в следующем примере:

sm4 = sm1 + sm2 + sm3;

sm1, sm2 и sm3 должны быть все строчно-главными или все столбцово-главными. С другой стороны, нет ограничений на целевую матрицу sm4. Например, это означает, что для вычисления \( A^T + A \) матрица \( A^T \) должна быть вычислена во временной матрице с совместимым порядком хранения:

SparseMatrix<double> A, B;
B = SparseMatrix<double>(A.transpose()) + A;

Бинарные операторы по коэффициентам могут также объединять разреженные и плотные выражения:

sm2 = sm1.cwiseProduct(dm1);
dm2 = sm1 + dm1;
dm2 = dm1 - sm1;

С точки зрения производительности, сложение/вычитание разреженных и плотных матриц лучше выполнять в два этапа. Например, вместо dm2 = sm1 + dm1, лучше написать:

dm2 = dm1;
dm2 += sm1;

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

Разреженные выражения также поддерживают транспонирование:

sm1 = sm2.transpose();
sm1 = sm2.adjoint();

Однако метод transposeInPlace() отсутствует.

Матричные произведения

Eigen поддерживает различные типы произведений разреженных матриц, которые обобщены ниже:

• разреженная-плотная:

  dv2 = sm1 * dv1;
  dm2 = dm1 * sm1.adjoint();
  dm2 = 2. * sm1 * dm1;
  

• симметричная разреженная-плотная. Произведение симметричной разреженной матрицы на плотную матрицу (или вектор) также можно оптимизировать, указав симметрию с помощью selfadjointView():

  dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
  dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
  dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored

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

  sm3 = sm1 * sm2;
  sm3 = 4 * sm1.adjoint() * sm2;
  

  Второй алгоритм удаляет явные нули или значения, меньшие заданного порога на лету. Он включен и контролируется функциями prune():

  sm3 = (sm1 * sm2).pruned();                  // removes numerical zeros
  sm3 = (sm1 * sm2).pruned(ref);               // removes elements much smaller than ref
  sm3 = (sm1 * sm2).pruned(ref,epsilon);       // removes elements smaller than ref*epsilon
  

• перестановки. Наконец, перестановки также могут быть применены к разреженным матрицам:

  PermutationMatrix<Dynamic,Dynamic> P = ...;
  sm2 = P * sm1;
  sm2 = sm1 * P.inverse();
  sm2 = sm1.transpose() * P;
  

Блочные операции

Что касается чтения, разреженные матрицы предоставляют тот же API, что и для плотных матриц, для доступа к подматрицам, таким как блоки, столбцы и строки. См. Блочные операции для подробного введения. Однако по причинам производительности запись в подматрицу разреженной матрицы намного ограниченнее, и в настоящее время только непрерывные наборы столбцов (соответственно, строк) столбцово-главной (соответственно, строчно-главной) SparseMatrix могут быть записываемыми. Более того, эта информация должна быть известна на этапе компиляции, что исключает такие методы, как block(...) и corner*(...). Доступный API для записи в SparseMatrix приведен ниже:

SparseMatrix<double,ColMajor> sm1;
sm1.col(j) = ...;
sm1.leftCols(ncols) = ...;
sm1.middleCols(j,ncols) = ...;
sm1.rightCols(ncols) = ...;
 
SparseMatrix<double,RowMajor> sm2;
sm2.row(i) = ...;
sm2.topRows(nrows) = ...;
sm2.middleRows(i,nrows) = ...;
sm2.bottomRows(nrows) = ...;

Кроме того, разреженные матрицы предоставляют методы SparseMatrixBase::innerVector() и SparseMatrixBase::innerVectors(), которые являются псевдонимами методов col/middleCols для столбцово-главного хранения и методов row/middleRows для строчно-главного хранения.

Треугольные и самосопряжённые представления

Так же, как и с плотными матрицами, функция triangularView() может использоваться для обращения к треугольной части матрицы и выполнения треугольных решений с плотным правым членом:

dm2 = sm1.triangularView<Lower>(dm1);
dv2 = sm1.transpose().triangularView<Upper>(dv1);

Функция selfadjointView() позволяет выполнять различные операции:

• оптимизированные разреженные-плотные матричные произведения:

  dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
  dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
  dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored
  

• копирование треугольных частей:

  sm2 = sm1.selfadjointView<Upper>();                               // makes a full selfadjoint matrix from the upper triangular part
  sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>();      // copies the upper triangular part to the lower triangular part
  

• применение симметричных перестановок:

  PermutationMatrix<Dynamic,Dynamic> P = ...;
  sm2 = A.selfadjointView<Upper>().twistedBy(P);                                // compute P S P' from the upper triangular part of A, and make it a full matrix
  sm2.selfadjointView<Lower>() = A.selfadjointView<Lower>().twistedBy(P);       // compute P S P' from the lower triangular part of A, and then only compute the lower part
  

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