[Eigen中文文档] 稀疏矩阵操作

这篇具有很好参考价值的文章主要介绍了[Eigen中文文档] 稀疏矩阵操作。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

文档总目录

英文原文(Sparse matrix manipulations)

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

模块 头文件 内容
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 指的是内部向量,它是列主矩阵的列,或行主矩阵的行。Outer 指的是另一个方向。

如下例子中更好的解释了这个存储方案:

有矩阵:
[ 0 3 0 0 0 22 0 0 0 17 7 5 0 1 0 0 0 0 0 0 0 0 14 0 8 ] \begin{bmatrix} 0&3&0&0&0 \\ 22&0&0&0&17 \\ 7&5&0&1&0 \\ 0&0&0&0&0 \\ 0&0&14&0&8 \\ \end{bmatrix} 0227003050000001400100017008
其可能的稀疏列主元表示之一为:

Values:			22	7	_	3	5	14	_	_	1	_	17	8
InnerIndices:	1	2	_	0	2	4	_	_	2	_	1	4
OuterStarts:	0	3	5	8	10	12
InnerNNZs:		2	2	1	1	2	

目前,保证给定内部向量的元素始终按递增的内部索引排序。_ 表示可用于快速插入新元素的可用空间。假设不需要重新分配,随机元素的插入复杂度为 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 对外部库的大多数包装器都需要压缩矩阵作为输入。

Eigen 运算的结果总是产生压缩的稀疏矩阵。另一方面,将新元素插入 SparseMatrix 会将其转换为未压缩模式。

这是先前以压缩模式表示的矩阵:

Values:			22	7	3	5	14	1	17	8
InnerIndices:	1	2	0	2	4	2	1	4
OuterStarts:	0	2	4	5	6	8

SparseVectorSparseMatrix 的特例,其中仅存储 ValuesInnerIndices 数组。SparseVector 没有压缩/未压缩模式的概念。

第一个示例

在描述每个单独的类之前,让我们从以下典型示例开始:使用有限差分格式和 Dirichlet 边界条件在规则二维网格上求解拉普拉斯方程 Δ u = 0 Δu=0 Δu=0 。这样的问题可以在数学上表示为以下形式的线性问题, A x = b Ax = b Ax=b,其中 x x x 是大小为 m m m 的未知向量(在我们的例子中,是像素的值), b b b 是由边界条件产生的右侧矢量,并且 A A A 是一个大小为 m ∗ m m*m mm 由拉普拉斯算子离散化产生的仅包含少数非零元素的矩阵。

示例代码如下:

#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> 的三元组列表。三元组是将非零条目表示为三元组的简单对象:行索引、列索引、值。

在 main 函数中,我们声明了一个三元组系数列表(作为标准矢量)和由 buildProblem 函数填充的右侧矢量 b。然后将非零系数的原始和平面列表转换为真正的 SparseMatrix 对象 A。请注意,列表的元素不必排序,可能的重复系数将被合并。

最后一步是有效地解决组装的问题。由于生成的矩阵 A 在结构上是对称的,可以通过 SimplicialLDLT 类执行直接的 Cholesky 分解,该类的行为类似于密集对象的 LDLT 对应操作。

生成的向量 x 包含作为一维数组的像素值,该数组保存到如下图所示的 jpeg 文件中。

[Eigen中文文档] 稀疏矩阵操作

buildProblemsaveAsBitmap 函数超出了本教程的范围。出于好奇和可重复性的目的,可在这里下载。

SparseMatrix 类

矩阵和向量属性

SparseMatrixSparseVector 类采用三个模板参数:标准类型(例如 double) 存储顺序(ColMajorRowMajor,默认为 ColMajor) 内部索引类型(默认为 int)。

对于稠密矩阵对象,构造函数需要传入对象的大小。这里有些例子:

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

在本教程的其余部分,matvec分别表示任何稀疏矩阵和稀疏向量对象。

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

// 标准维度
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表示一个类型;更多详细信息参见 The template and typename keywords in C++

填充稀疏矩阵

由于 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!

std::vector中的三元组元素可能以任意顺序存储,并且可能包含重复元素,这些元素将被setFromTriplets()函数进行求和处理。有关更多详细信息,请参阅SparseMatrix :: setFromTriplets()函数和类Triplet。

然而,在某些情况下,可以通过直接将非零元素插入目标矩阵来实现稍高的性能和较低的内存消耗。这种方法的典型场景如下所示:

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
  • 这里的关键是第2行,我们为每列预留了6个非零项的空间。在许多情况下,每列或每行的非零元素数量可以很容易地预先知道。如果每个内部向量的非零元素数量差异很大,则可以通过提供具有operator[](int j)方法(返回第 j 个内部向量的保留大小)的向量对象(例如通过VectorXistd::vector<int>实现)来为每个内部向量指定保留大小。如果只能得到每个内部向量的非零元素数量的粗略估计,则强烈建议将其高估而不是低估。如果省略此行,则将为新元素的第一次插入每个内部向量保留2个元素的空间。

  • 第4行执行了一个排序插入。在这个例子中,理想情况是第 j 列不是满的,并且包含内部索引比 i 小的非零元素。在这种情况下,此操作可以简化为 O(1) 操作。

  • 调用insert(i,j)时,元素ij不能已经存在,否则使用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;

sm1sm2sm3必须全是行优先或者全是列优先。另一方面,目标矩阵sm4没有限制。例如,这意味着对于 A T + A A^T+A AT+A,矩阵 A T A^T AT 必须被计算成一个具有兼容存储顺序的临时矩阵:

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

二项式系数的操作符也可以混合稀疏和稠密表达式:

m2 = 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 sm1 are stored
dm2 = sm1.selfadjointView<Upper>() * dm1;   // if only the upper part of sm1 is stored
dm2 = sm1.selfadjointView<Lower>() * dm1;   // if only the lower part of sm1 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。有关详细介绍,请参见块操作。但是,出于性能考虑,向子稀疏矩阵的写入要受到更多限制,目前仅限于连续的列(或行)集,这些列(或行)是列主(或行主)稀疏矩阵可写的。此外,这些信息必须在编译时知道,不能使用诸如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) = ...;

此外,稀疏矩阵提供了SparseMatrixBase::innerVector()SparseMatrixBase::innerVectors()方法,它们是针对列优先存储的col/middleCols方法和针对行优先存储的row/middleRows方法的别名。

三角形视图和自共轭视图

与密集矩阵一样,triangularView()函数可用于处理矩阵的三角部分,并对具有密集右侧的矩阵进行三角求解:

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

selfadjointView()函数允许进行各种操作。文章来源地址https://www.toymoban.com/news/detail-497914.html

  • 优化稀疏-稠密矩阵乘积
dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of sm1 are stored
dm2 = sm1.selfadjointView<Upper>() * dm1;   // if only the upper part of sm1 is stored
dm2 = sm1.selfadjointView<Lower>() * dm1;   // if only the lower part of sm1 is stored
  • 复制三角形部分
// makes a full selfadjoint matrix from the upper triangular part
sm2 = sm1.selfadjointView<Upper>();  
// copies the upper triangular part to the lower triangular part
sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>();      
  • 对称排列的应用
PermutationMatrix<Dynamic,Dynamic> P = ...;
// compute P S P' from the upper triangular part of A, and make it a full matrix
sm2 = A.selfadjointView<Upper>().twistedBy(P);     
// compute P S P' from the lower triangular part of A, and then only compute the lower part
sm2.selfadjointView<Lower>() = A.selfadjointView<Lower>().twistedBy(P);       

到了这里,关于[Eigen中文文档] 稀疏矩阵操作的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • [Eigen中文文档] 在 CMake 项目中使用 Eigen

    文档总目录 英文原文(Using Eigen in CMake Projects) Eigen提供了CMake支持,使得该库可以轻松地在CMake项目中使用。 注意:启用这个功能需要CMake 3.0(或更高版本)。 Eigen提供了一个CMake示例,名为 Eigen3::Eigen ,可以使用 find_package CMake 命令导入,并通过调用 target_link_libraries 来使用,

    2024年02月17日
    浏览(43)
  • [Eigen中文文档] 线性代数与分解

    文档总目录 英文原文(Linear algebra and decomposition) 本节说明如何求解线性系统,计算各种分解,如 LU 、 QR 、 SVD 、 特征分解 …… 求解基本线性系统 问题 :有一个方程组,写成矩阵方程如下: A x = b Ax = b A x = b 其中 A A A 和 b b b 是矩阵(作为一种特殊情况, b b b 也可以是一个

    2024年02月07日
    浏览(43)
  • Eigen 矩阵Matrix及其简单操作

    在Eigen,所有的矩阵和向量都是Matrix模板类的对象,Vector只是一种特殊的矩阵(一行或者一列)。 Matrix有6个模板参数,主要使用前三个参数,剩下的有默认值。 Scalar是表示元素的类型,RowsAtCompileTime为矩阵的行,ColsAtCompileTime为矩阵的列。 库中提供了一些类型便于使用,比如

    2024年02月12日
    浏览(35)
  • 三元组操作(相加)——稀疏矩阵(c语言)

     运行环境:TDM-GCC 三元组用来存储 稀疏矩阵 比较节省空间,因为稀疏矩阵大部分都是零元素,而三元组只记录非零元素。 这里两个相加的矩阵有着同样的 i行 j列。 运行结果:         行数:3 列数:4  元素个数:4         --------------------         第0 行  第1 列  1    

    2024年02月05日
    浏览(48)
  • SLAM——Eigen函数库之矩阵块运算,高阶操作middleCols与segment用法

    Eigen/四元数/欧拉角/旋转矩阵 相关系列文章 Eigen/Matlab 使用小结 SLAM——之Eigen入门(矩阵运算及几何模块) SLAM——之Eigen函数库,一个相对复杂的EIgen使用实例 SLAM——Eigen函数库:矩阵块运算,block操作 SLAM——Eigen函数库之 Eigen::Ref 使用实例 欧拉角和旋转矩阵相互转换 四元数

    2024年02月10日
    浏览(43)
  • 用三元组表实现稀疏矩阵的基本操作

    目录 问题描述 数据结构 算法设计 算法流程图  源代码  运行结果      ​    编写程序用三元组表实现稀疏矩阵的按列转置操作。 本设计使用三元组表实现。 程序中设计了三个函数: 1.函数InitSPNode()用来建立一个稀疏矩阵的三元组表。     首先输入行数、列数和非零元的

    2024年02月03日
    浏览(38)
  • 【数据结构】数组和字符串(十):稀疏矩阵的链接存储:十字链表的矩阵操作(加法、乘法、转置)

    【数据结构】数组和字符串(一):矩阵的数组表示   矩阵是以按行优先次序将所有矩阵元素存放在一个一维数组中。但是对于特殊矩阵,如对称矩阵、三角矩阵、对角矩阵和稀疏矩阵等, 如果用这种方式存储,会出现大量存储空间存放重复信息或零元素的情况,这样会造

    2024年02月08日
    浏览(60)
  • Android修行手册-POI操作中文API文档

    Unity3D特效百例 案例项目实战源码 Android-Unity实战问题汇总 游戏脚本-辅助自动化 Android控件全解手册 再战Android系列 Scratch编程案例 软考全系列 Unity3D学习专栏 蓝桥系列 ChatGPT和AIGC 专注于 Android/Unity 和各种游戏开发技巧,以及 各种资源分享 (网站、工具、素材、源码、游戏等

    2024年01月16日
    浏览(44)
  • 矩阵分解及其Eigen实现

    主要是用来记录自己的学习过程,内容也主要来自于网上的各种资料,然后自己总结而来,参考的资料都以注明,感谢这些作者的分享。如果内容有误,请大家指点。 定义        将矩阵等价为两个矩阵 L L L 和 U U U 的乘积 ,其中 L L L 和 U U U 分别是单位下三角矩阵和上三角

    2024年02月03日
    浏览(47)
  • Eigen-Matrix矩阵

    在Eigen中,所有矩阵和向量都是矩阵模板类的对象。向量只是矩阵的一种特殊情况,要么有一行,要么有一列。矩阵就是一个二维数表,可以有多行多列。 Matrix类有六个模板参数,但现在只需要了解前三个参数就足够了。剩下的三个参数都有默认值,我们暂时不碰它们,我们

    2024年03月09日
    浏览(66)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包