Eigen教程3 - 稀疏矩阵操作

博客新址: http://blog.xuezhisd.top
邮箱：xuezhisd@126.com

---

稀疏矩阵操作

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

• 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.14159265358979323846
typedef 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*sm1
sm1+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