Eigen教程-sparse

转载 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(...); //只读