Eigen教程3 - 稀疏矩阵操作

稀疏矩阵操作

操作和求解稀疏问题需要的模块：

• SparseCore
  
  • SparseMatrix 和 SparseVector 类，基本线性代数（包括三角求解器）
• SparseCholesky
  
  • 稀疏LLT和LDLTCholesky分解，解决稀疏正定问题。
• SparseLU
  
  • 稀疏LU分解
• SparseQR
  
  • 稀疏QR分解
• IterativeLinearSolvers
• Sparse
  
  • 包含上述所有的模块。

稀疏矩阵的格式

• SparseMatrix是Eigen的稀疏模块的最主要的稀疏矩阵。它实现了Compressed Column (or Row) Storage方案。

• SparseMatrix包含了4个精简的数组：

  • Values: 存储非零元素
  • InnerIndices: 存储非零元素的行/列索引.
  • OuterStarts: 存储每一列/行第一个非零元素在上面两个数组中的索引。
  • InnerNNZs: 存储每一列/行中非零元素的数量。Inner在列优先矩阵中表示一个列，在行优先矩阵中表示一个行。Outer表示另一个方向。

• 如果中间没有留未占用的空间，就是压缩模式。 对应于Compressed Column (or Row) Storage Schemes (CCS or CRS)。

• 任何SparseMatrix都可以通过SparseMatrix::makeCompressed() 方法变成稀疏模式。此时，InnerNNZs相对于OuterStarts而言，就是多余的，因为InnerNNZs[j] = OuterStarts[j+1]-OuterStarts[j]。因此SparseMatrix::makeCompressed()会释放InnerNNZ存储空间。
• Eigen的操作总是返回压缩模式的稀疏矩阵。当向一个SparseMatrix插入新元素时，会变成非压缩模式。
• SparseVector是SparseMatrix的一个特例，只有 Values 和 InnerIndices数组。没有压缩和非压缩的区别。

例 1

• 求解线性方程组：$Ax=b$。
• A是一个mxm的稀疏矩阵。

// 参考链接：http://eigen.tuxfamily.org/dox/group__TutorialSparse.html#include <Eigen/Sparse>#include <vector>//#include <QImage> //Qt#define M_PI       3.14159265358979323846typedef Eigen::SparseMatrix<double> SpMat; // 声明一个列优先的双精度稀疏矩阵类型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){    assert(argc==2);    int n = 300;  // 图像的尺寸    int m = n*n;  // 未知元素的数量（等于像素数）    // Assembly:    std::vector<T> coefficients;            // list of non-zeros coefficients    Eigen::VectorXd b(m);                   // 等号右边的向量b，根据约束条件生成    buildProblem(coefficients, b, n);    SpMat A(m,m); // 等号左边的矩阵A    A.setFromTriplets(coefficients.begin(), coefficients.end());    // 求解    Eigen::SimplicialCholesky<SpMat> chol(A);  // 执行A的 Cholesky分解    Eigen::VectorXd x = chol.solve(b);         // 使用A的Cholesky分解来求解等号右边的向量b    // Export the result to a file:    //saveAsBitmap(x, n, argv[1]);    return 0;}void insertCoefficient(int id, int i, int j, double w, std::vector<T>& coeffs,    Eigen::VectorXd& b, const Eigen::VectorXd& boundary){    int n = int(boundary.size());    int id1 = i+j*n;    if(i==-1 || i==n) b(id) -= w * boundary(j); // constrained coefficient    else  if(j==-1 || j==n) b(id) -= w * boundary(i); // constrained coefficient    else  coeffs.push_back(T(id,id1,w));              // unknown coefficient}void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n){    b.setZero();    Eigen::ArrayXd boundary = Eigen::ArrayXd::LinSpaced(n, 0,M_PI).sin().pow(2);    for(int j=0; j<n; ++j)    {        for(int i=0; i<n; ++i)        {            int id = i+j*n;            insertCoefficient(id, i-1,j, -1, coefficients, b, boundary);            insertCoefficient(id, i+1,j, -1, coefficients, b, boundary);            insertCoefficient(id, i,j-1, -1, coefficients, b, boundary);            insertCoefficient(id, i,j+1, -1, coefficients, b, boundary);            insertCoefficient(id, i,j,    4, coefficients, b, boundary);        }    }}/*void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename){    Eigen::Array<unsigned char,Eigen::Dynamic,Eigen::Dynamic> bits = (x*255).cast<unsigned char>();    QImage img(bits.data(), n,n,QImage::Format_Indexed8);    img.setColorCount(256);    for(int i=0;i<256;i++) img.setColor(i,qRgb(i,i,i));    img.save(filename);}*/

• 小结：上述代码中的函数saveAsBitmap()需要Qt，因此没有调试成功。
• Triplet是一个三元组，存储一个非零元素：（行，列，值）。

SparseMatrix类

矩阵和向量的属性

• SparseMatrix类和SparseVector类有3个模板参数：元素类型，存储顺序，（ColMajor(默认)或RowMajor）和内索引类型（默认int）。
• 对于密集Matrix类，构造函数的参数是对象的大小。

SparseMatrix<std::complex<float> > mat(1000,2000);         //声明一个 1000x2000 列优先的压缩的稀疏矩阵（元素类型：complex<float>）SparseMatrix<double,RowMajor> mat(1000,2000);              //声明一个 1000x2000 行优先的压缩的稀疏矩阵（元素类型：double）SparseVector<std::complex<float> > vec(1000);              //声明一个1000维的稀疏的列向量，元素类型为complex<float>SparseVector<double,RowMajor> vec(1000);                   //声明一个1000维的稀疏的行向量，元素类型为double

• 下面的代码，演示了如何访问稀疏矩阵和向量的维度：

// 参考链接：http://eigen.tuxfamily.org/dox/group__TutorialSparse.html#include <Eigen/Sparse>#include <iostream>using namespace std;using namespace Eigen;int main(int argc, char** argv){    SparseMatrix<double,RowMajor> mat(1000,2000);//创建一个行优先的，维度为1000x2000的稀疏矩阵    SparseVector<double,RowMajor> vec(1000);      //标准维度    cout << "mat.rows() = " << mat.rows() << endl;    cout << "mat.cols() = " << mat.cols() << endl;    cout << "mat.size() = " << mat.size() << endl;    cout << "vec.size() = " << vec.size() << endl;    //内/外维度    cout << "mat.innerSize() = " << mat.innerSize() << endl; //行优先，所以为列数    cout << "mat.outerSize() = " << mat.outerSize() << endl;    //非零元素个数    cout << "mat.nonZeros() = " << mat.nonZeros() << endl;    cout << "vec.nonZeros() = " << vec.nonZeros() << endl;    return 0;}

迭代所有的非零元素

• 使用coeffRef(i,j)函数可以随机访问稀疏对象的元素，但是用到的二叉搜索比较浪费时间。
• 大多数情况下，我们仅需迭代所有的非零元素。这时，可以先迭代外维度，然后迭代当前内维度向量的非零元素。非零元素的访问顺序和存储顺序相同？

// 参考链接：http://eigen.tuxfamily.org/dox/group__TutorialSparse.html#include <Eigen/Sparse>#include <iostream>using namespace std;using namespace Eigen;int main(int argc, char** argv){    SparseMatrix<double,RowMajor> mat(5,7);//创建一个行优先的，维度为1000x2000的稀疏矩阵    //随机访问（读/写）元素    cout << mat.coeffRef(0,0) << endl; //读取元素    mat.coeffRef(0,0) = 5;    cout << mat.coeffRef(0,0) << endl; //写入元素    //迭代访问稀疏矩阵    cout << "\n迭代访问稀疏矩阵的元素：" << endl;    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();   // 列标（此处等于k）            it.index(); // 内部索引，此处等于it.row()        }    cout << mat << endl;    //迭代访问稀疏向量    SparseVector<double, RowMajor> vec(5);    //vec.coeffRef(0,0) = 1;    for (SparseVector<double>::InnerIterator it(vec); it; ++it)    {        it.value(); // == vec[ it.index() ]        it.index(); //索引    }    cout << vec << endl;    return 0;}

填充稀疏矩阵

• 由于稀疏矩阵的特殊存储格式，添加新的非零元素时需要特别注意。例如，随机插入一个非零元素的复杂度是O(nnz)，nnz表示非零元素的个数。
• 最简单的保证性能的创建稀疏矩阵的方法是，先创建一个三元组列表，然后转换成SparseMatrix。如下所示：

// 参考链接：http://eigen.tuxfamily.org/dox/group__TutorialSparse.html#include <Eigen/Sparse>#include <iostream>using namespace std;using namespace Eigen;typedef Eigen::Triplet<double> T;typedef Eigen::SparseMatrix<double> SparseMatrixType;int main(int argc, char** argv){    int rows=10, cols = 10;    int estimation_of_entries = 10; //预计非零元素的个数    std::vector<T> tripletList;    tripletList.reserve(estimation_of_entries);    int j = 1; // 列标    for(int i=0; i<estimation_of_entries; i++){ //遍历行        int v_ij = i*i+1;        tripletList.push_back(T(i,j,v_ij));    }    SparseMatrixType mat(rows,cols);    mat.setFromTriplets(tripletList.begin(), tripletList.end()); //根据三元组列表生成稀疏矩阵    //     cout << mat << endl;    return 0;}

• 直接将非零元素插入到目标矩阵的方法更加高效，节省内存。如下所示：

// 参考链接：http://eigen.tuxfamily.org/dox/group__TutorialSparse.html#include <Eigen/Sparse>#include <iostream>using namespace std;using namespace Eigen;int main(int argc, char** argv){    int rows=10, cols = 10;    SparseMatrix<double> mat(rows,cols);         // 默认列优先    mat.reserve(VectorXi::Constant(cols,6)); //关键：为每一列保留6个非零元素空间    for(int i=0; i<3; i++){ //遍历行        for(int j=0;j<3; j++){            int v_ij = i+j+1;            mat.insert(i,j) = v_ij;                    // alternative: mat.coeffRef(i,j) += v_ij;        }    }    mat.makeCompressed(); //压缩剩余的空间    //     cout << mat << endl;    return 0;}

支持的操作和函数

• 由于稀疏矩阵的特殊存储格式，它没有密集矩阵那样的灵活性。Eigen的稀疏矩阵模块仅仅实现了密集矩阵API的一个子集。
• 下面使用sm表示稀疏矩阵，sv表示稀疏向量，dm表示密集矩阵，dv表示密集向量。

基本操作

• 稀疏矩阵支持大多数的逐元素的一元操作符和二元操作符。

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

• 注意：强制约束条件：它们的存储顺序必须匹配（相同）。比如，sm4 = sm1 + sm2 + sm3;，sm1，sm2，sm3必须都是行优先或列优先的；sm4没有要求。
• 计算$A^T+A$时，必须先将前一项生成一个兼容顺序的临时矩阵，然后再计算。如下所示：

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

• 二元操作符支持稀疏矩阵和密集矩阵的混合操作。

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

• 为了提升性能，稀疏矩阵加减法可以分成两步进行。好处是，完全发挥密集矩阵存储的高性能特性，仅在稀疏矩阵的非零元素上花费时间。仅如下所示。

dm2 = dm1;dm2 += sm1;

• 稀疏矩阵也支持转置，但没有原地计算的转置函数transposeInPlace()。

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

矩阵乘法

• Eigen支持不同的稀疏矩阵乘法运算。

稀疏矩阵和密集矩阵乘法：

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

对称稀疏矩阵和密集矩阵乘法：

dm2 = sm1.selfadjointView<>() * dm1;        // 当存储了A的所有元素dm2 = A.selfadjointView<Upper>() * dm1;     //仅当存储了A的上半部分dm2 = A.selfadjointView<Lower>() * dm1;     // 仅当存储了A的下半部分

稀疏矩阵间的乘法

• 稀疏矩阵之间的乘法有两种算法。默认的方法是一种保守的方法，保存可能出现的零。如下所示：

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

• 稀疏矩阵之间的乘法的第二种算法将会裁剪零或非常小的值，它通过prune()函数激活和控制。

sm3 = (sm1 * sm2).pruned();                  // 去除数值零sm3 = (sm1 * sm2).pruned(ref);               // 去除比ref小的元素sm3 = (sm1 * sm2).pruned(ref,epsilon);       // 去除比ref*epsilon小的元素

排列

• 稀疏矩阵也可以进行排列操作。

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

块操作

• 关于读操作，稀疏矩阵支持类似于密集矩阵相同的接口来访问块，列，行。
• 但由于性能的原因，稀疏子矩阵的写操作非常受限。当前仅列（列）优先稀疏矩阵的连续列（或行）可写。此外，在编译时必须知道这些信息。
• 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) = ...;

• 稀疏矩阵的2个方法SparseMatrixBase::innerVector() 和 SparseMatrixBase::innerVectors()，当稀疏矩阵是列优先存储方式时，这两个方法绑定到col/middleCols方法，当稀疏矩阵是行优先存储方式时，这两个方法绑定到row/middleRows方法。

三角和selfadjoint

• triangularView()函数可以用于解决矩阵的而三角部分，给定右边密集矩阵求解三角**。

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

• selfadjointView()函数支持不同的操作。

dm2 = sm1.selfadjointView<>() * dm1;        // 当存储了A的所有元素dm2 = A.selfadjointView<Upper>() * dm1;     //仅当存储了A的上半部分dm2 = A.selfadjointView<Lower>() * dm1;     // 仅当存储了A的下半部分

• 复制三角部分

sm2 = sm1.selfadjointView<Upper>();                               // 从上三角部分生成一个全的伴随矩阵sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>();      // 复制上三角部分到下三角部分

• 对称排列

PermutationMatrix<Dynamic,Dynamic> P = ...;sm2 = A.selfadjointView<Upper>().twistedBy(P);                                // 从A的上三角部分计算 P S P'，生成一个全矩阵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