Eigen教程-sparse

本文详细介绍了Eigen库中稀疏矩阵的基本操作,包括初始化、变形、赋值、插入等,并覆盖了算术运算、矩阵属性获取、转置、排列、逐元素操作等方面的知识。
转载 http://blog.youkuaiyun.com/xuezhisdc/article/details/54633274
  • 本文对稀疏矩阵SparseMatrix的主要操作进行了总结。首先,建议先阅读《Eigen教程2 - 稀疏矩阵操作》。
  • 关于稀疏矩阵,最重要的一点是:稀疏矩阵的存储方式,是按列优先储存,还是按行优先存储。
  • 绝大多数的稀疏矩阵的算术操作都会断言(判断)操作数的存储方式是否相同。

稀疏矩阵初始化

构造函数

  • 默认列优先。
SparseMatrix<double> sm1(1000,1000);  //创建一个1000x1000的双精度的稀疏矩阵
SparseMatrix<std::complex<double>,RowMajor> sm2; //声明一个复数,行优先的稀疏矩阵

变形和保留

  • 注意:调用reserve()时,没必要给定最终矩阵的准确的非零元素个数。但是,给定准确的非零元素的个数,可以避免多次重新分配内存。
sm1.resize(m,n);      // 将sm1变形成mxn的矩阵
sm1.reserve(nnz);     // 为nnz个非零元素分配空间  

赋值

  • 复制构造函数可以将一种存储顺序变成另一种顺序。
SparseMatrix<double,Colmajor> sm1;
//使用sm1初始化sm2
SparseMatrix<double,Rowmajor> sm2(sm1), sm3;        
// 赋值(自动修改存储顺序)
sm3 = sm1; 

逐元素插入

  • insert()假定该元素不存在;如果元素已经存在,使用coeffRef()。
// 插入一个新元素 
 sm1.insert(i, j) = v_ij;  
// 更新value v_ij
 sm1.coeffRef(i,j) = v_ij;
 sm1.coeffRef(i,j) += v_ij;
 sm1.coeffRef(i,j) -= v_ij;

批插入

std::vector< Eigen::Triplet<double> > tripletList; //三元组列表
tripletList.reserve(estimation_of_entries); //保留非零元素空间
//先填充tripletList,然后将其存储到稀疏矩阵中
sm1.setFromTriplets(TripletList.begin(), TripletList.end());

常量或随机插入

  • 去除所有的非零元素
sm1.setZero()

矩阵属性

  • 下面这些函数可以获取矩阵的信息:
sm1.rows();         // 行数
sm1.cols();         // 列数 
sm1.nonZeros();     // 非零元素的个数   
sm1.outerSize();    // 列优先的稀疏矩阵的行数,行优先的稀疏矩阵的列数
sm1.innerSize();    // 列优先的稀疏矩阵的列数,行优先的稀疏矩阵的行数
sm1.norm();          // 矩阵的欧式范数
sm1.squaredNorm();  //矩阵的平方范数(平方和)
sm1.blueNorm();   // 
sm1.isVector();     // 检查是否是稀疏向量或稀疏矩阵
sm1.isCompressed(); // 检查是否是压缩格式

算术操作

  • 在具有相同存储顺序和维度已知的稀疏矩阵上执行算术操作非常容易。
  • 注意:在不同存储顺序的矩阵间也可以执行算术运算。
  • sm表示稀疏矩阵,dm表示密集矩阵,dv表示密集向量。

加减操作

  • sm1和sm2应具有相同的存储顺序。
sm3 = sm1 + sm2; 
sm3 = sm1 - sm2;
sm2 += sm1; 
sm2 -= sm1; 

标量乘法

  • 只要维度和存储顺序允许,可以组合多种运算。
sm3 = sm1 * s1;   sm3 *= s1; 
sm3 = s1 * sm1 + s2 * sm2; sm3 /= s1;

稀疏乘法

sm3 = sm1 * sm2;
dm2 = sm1 * dm1;
dv2 = sm1 * dv1;

转置 和 伴随矩阵

  • 注意:转置运算改变存储顺序,不支持原地转置运算transposeInPlace()。
sm2 = sm1.transpose();
sm2 = sm1.adjoint();

排列

perm.indices();      // Reference to the vector of indices
sm1.twistedBy(perm); // Permute rows and columns
sm2 = sm1 * perm;    // Permute the columns
sm2 = perm * sm1;    // Permute the columns

逐元素操作

  • sm1和sm2的存储顺序应相同。
sm1.cwiseProduct(sm2);
sm1.cwiseQuotient(sm2);
sm1.cwiseMin(sm2);
sm1.cwiseMax(sm2);
sm1.cwiseAbs();
sm1.cwiseSqrt();

其它支持的操作

子矩阵

  • 相对于密集矩阵,这里的所有函数都是只读的。
sm1.block(startRow, startCol, rows, cols); 
sm1.block(startRow, startCol); 
sm1.topLeftCorner(rows, cols); 
sm1.topRightCorner(rows, cols);
sm1.bottomLeftCorner( rows, cols);
sm1.bottomRightCorner( rows, cols);

Range

  • 在行优先稀疏矩阵中,内向量是指一个行;在列优先稀疏矩阵中,内向量是指一个列。
  • 对于可以读写(RW)的子矩阵,可以在不同存储顺序的矩阵上运算。
sm1.innerVector(outer);           // RW
sm1.innerVectors(start, size);    // RW
sm1.leftCols(size);               // RW
sm2.rightCols(size);              // RO 只读,因为sm2是行优先的
sm1.middleRows(start, numRows);   // RO 只读,因为sm1是列优先的
sm1.middleCols(start, numCols);   // RW
sm1.col(j);                       // RW

三角和自伴随视图

  • 允许三角视图和块视图之间的组合使用。
sm2 = sm1.triangularview<Lower>();
sm2 = sm1.selfadjointview<Lower>();

三角求解

  • 对于一般的稀疏求解,使用《Eigen教程4 - 求解稀疏线性方程组》中合适的模块。
dv2 = sm1.triangularView<Upper>().solve(dv1);
dv2 = sm1.topLeftCorner(size, size).triangularView<Lower>().solve(dv1);

低级API

  • 如果矩阵是非压缩模式,应首先调用makeCompressed()函数。
  • 注意:下面这些函数主要用于和外部库进行交互。
  • 更好的访问矩阵的元素值的方法是使用InnerIterator类,详情查看《Eigen教程2 - 稀疏矩阵操作》
sm1.valuePtr();      // 指向数值的指针
sm1.innerIndextr();  // 指向索引的指针
sm1.outerIndexPtr(); // 指向每一个内向量开头的指针

映射外部缓存

int outerIndexPtr[cols+1];
int innerIndices[nnz];
double values[nnz];
Map<SparseMatrix<double> > sm1(rows,cols,nnz,outerIndexPtr, innerIndices,values); //RW 可读写
Map<const SparseMatrix<double> > sm2(...); //只读
### 使用 Eigen 库中的 SparseMatrix 类 #### 创建和初始化稀疏矩阵对象 为了创建一个 `SparseMatrix` 对象,可以指定其维度以及动态分配的空间大小。通常建议预先估计非零元素的数量以优化性能。 ```cpp #include <Eigen/Sparse> // 定义一个双精度浮点数类型的稀疏矩阵, 行优先存储方式 Eigen::SparseMatrix<double> m(rows, cols); m.reserve(Eigen::VectorXi::Constant(cols, 8)); // 预留每列大约有8个非零项 ``` #### 插入非零元素到稀疏矩阵中 通过调用 `.insert()` 方法可以在特定位置插入新的非零条目;如果该位置已经存在,则会更新现有值。注意,在完成所有插入之前不应访问矩阵数据结构。 ```cpp for(int k = 0; k < num_entries; ++k){ int i = row_indices[k]; int j = col_indices[k]; double v = values[k]; // 如果(i,j)处还没有设置过数值则插入新值v m.insert(i, j) = v; } ``` #### 执行基本线性代数运算 一旦构建好稀疏矩阵之后就可以像对待常规密集型矩阵那样执行各种算术计算了。例如求解线性方程组 Ax=b: ```cpp // 构建系数矩阵A (假设已知) Eigen::SparseMatrix<double> A; // 设置右侧向量b Eigen::VectorXd b(n); // 解决Ax=b问题得到x Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver(A); Eigen::VectorXd x = solver.solve(b); ``` 对于稠密矩阵,类 `Map<SparseMatrixType>` 可用于将外部缓冲区视为 EigenSparseMatrix 对象[^1]。这允许更灵活的数据交换机制而无需复制整个数组内容。 在许多应用场景下(比如有限元分析),经常会遇到规模巨大但大部分区域都接近于零的大尺寸矩阵。这时采用专门设计用来保存少量有效成分的形式——即所谓的“稀疏矩阵”,能够显著降低资源占用并加快算法效率[^3]。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值