Eigen关于稀疏矩阵

处理和解决稀疏问题涉及各种模块，总结如下：

模块	头文件	内容
SparseCore	#include <Eigen/SparseCore>	SparseMatrix 和 SparseVector 类、矩阵组装、 基本稀疏线性代数（包括稀疏三角求解器）
SparseCholesky	#include <Eigen/SparseCholesky>	直接用稀疏 LLT 和 LDLT Cholesky 分解来解决稀疏自伴随正定问题
SparseLU	#include<Eigen/SparseLU>	稀疏 LU 分解来解决一般平方稀疏系统
SparseQR	#include<Eigen/SparseQR	用于解决稀疏线性最小二乘问题的稀疏 QR 分解
IterativeLinearSolvers	#include <Eigen/IterativeLinearSolvers>	求解大型一般线性平方问题的迭代求解器（包括自伴随正定问题）
Sparse	#include <Eigen/Sparse>	包含以上所有模块

矩阵稀疏格式

在许多应用中（例如，有限元方法），通常处理非常大的矩阵，其中只有少数系数不为零。 在这种情况下，可以通过使用仅存储非零系数的专用表示来减少内存消耗并提高性能。 这样的矩阵称为稀疏矩阵。

SparseMatrix

SparseMatrix 类是 Eigen 稀疏模块的主要稀疏矩阵表示； 它提供高性能和低内存使用率。 它实现了广泛使用的压缩列（或行）存储方案的更通用的变体。 它由四个紧凑的数组组成：

• Values:存储非零系数值。
• InnerIndices:存储非零的行（或列）索引。
• OuterStarts:为每一列（或行）存储前两个数组中第一个非零的索引。
• InnerNNZs:存储每列（相应行）的非零数。 inner指的是一个内部向量，它是列主矩阵的列，或行主矩阵的行。 外部这个词指的是另一个方向。

目前，给定内部向量的元素保证总是通过增加内部索引进行排序。 “_”表示可以快速插入新元素的可用空间。 假设不需要重新分配，随机元素的插入因此在 O(nnz_j) 中，其中 nnz_j 是相应内部向量的非零数。 另一方面，在给定的内部向量中插入具有增加内部索引的元素效率更高，因为这仅需要增加相应的 InnerNNZs 条目，即 O(1) 操作。

没有可用空间的情况是一种特殊情况，称为压缩模式。 它对应于广泛使用的压缩列（或行）存储方案（CCS 或 CRS）。 通过调用 SparseMatrix::makeCompressed() 函数，可以将任何 SparseMatrix 转换为这种形式。 在这种情况下，可以注意到 InnerNNZs 数组与 OuterStarts 是冗余的，因为我们等式：InnerNNZs[j] = OuterStarts[j+1]-OuterStarts[j]。 因此，实际上调用 SparseMatrix::makeCompressed() 会释放此缓冲区。

例子

在描述每个单独的类之前，让我们从以下典型示例开始：使用有限差分格式和狄利克雷边界条件在规则 2D 网格上求解拉普拉斯方程 Δu=0。 这样的问题可以数学表达为 Ax=b 形式的线性问题，其中 x 是 m 个未知数的向量（在我们的例子中，像素的值），b 是由边界条件产生的右侧向量，并且 A 是一个 m×m 矩阵，仅包含少数非零元素，这些元素是由 Laplacian 算子离散化产生的。

稀疏矩阵类

SparseMatrix 和 SparseVector 类采用三个模板参数：标量类型（例如，double）、存储顺序（ColMajor 或 RowMajor，默认为 ColMajor）、内部索引类型（默认为 int）。

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

可以使用以下函数查询矩阵的维度：

Standard dimensions	mat.rows()mat.cols()	vec.size()
Sizes along the inner/outer dimensions	mat.innerSize()mat.outerSize()
Number of non zero coefficients	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()

}

填充稀疏矩阵

由于 SparseMatrix 的特殊存储方案，在添加新的非零条目时必须特别小心。 例如，单个纯随机插入到 SparseMatrix 的成本是 O(nnz)，其中 nnz 是当前非零系数的数量。

创建稀疏矩阵同时保证良好性能的最简单方法是首先构建所谓的三元组列表，然后将其转换为 SparseMatrix。

一下是一个简单的例子

typedef Eigen::Triplet<double>T;

std::vector.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!

The std::vector of triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets(). See the SparseMatrix::setFromTriplets() function and class Triplet for more details.

但是，在某些情况下，可以通过将非零值直接插入目标矩阵来实现稍微更高的性能和更低的内存消耗。 这种方法的典型场景如下图所示：

SparseMatrix<double> mat(rows,cols); // default is column major

mat.reserve(VectorXi::Constant(cols,6));

for each i,j such that v_ij != 0

mat.insert(i,j) = v_ij; // alternative: mat.coeffRef(i,j) += v_ij;

mat.makeCompressed(); // optional

第 4 行执行排序插入。 在此示例中，理想情况是第 j 列未满且包含内部索引小于 i 的非零值。 在这种情况下，此操作归结为简单的 O(1) 操作。

调用 insert(i,j) 时，元素 i ,j 必须不存在，否则使用 coeffRef(i,j) 方法，该方法将允许例如累加值。 此方法首先执行二分查找，如果元素不存在，则最后调用 insert(i,j)。 它比 insert() 更灵活，但成本也更高。

第 5 行抑制了剩余的空白空间，并将矩阵转换为压缩的列存储。

支持的运算符和函数

由于其特殊的存储格式，稀疏矩阵无法提供与密集矩阵相同级别的灵活性。 在 Eigen 的稀疏模块中，下面sm表示稀疏矩阵，sv表示稀疏向量，dm表示稠密矩阵，dv表示稠密向量。

Eigen 支持各种稀疏矩阵乘积，总结如下：

• sparse-dense:

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

Block operations

关于读取访问，稀疏矩阵公开与密集矩阵相同的 API，以访问子矩阵，例如块、列和行。 详细介绍见块操作。 然而，出于性能原因，写入子稀疏矩阵的限制要大得多，目前只有列主（分别为行主）稀疏矩阵的连续列集（分别为行）是可写的。 此外，必须在编译时知道这些信息，省去诸如 block(...) 和 corner*(...) 之类的方法。 对 SparseMatrix 进行写访问的可用 API 总结如下：

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) = ...;