博客新址: http://blog.xuezhisd.top
邮箱：xuezhisd@126.com
稀疏矩阵操作
操作和求解稀疏问题需要的模块：
- SparseCore
- SparseMatrix 和 SparseVector 类，基本线性代数（包括三角求解器）
- SparseCholesky
- 稀疏LLT和LDLTCholesky分解，解决稀疏正定问题。
- SparseLU
- SparseQR
- 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();
矩阵乘法
稀疏矩阵和密集矩阵乘法：
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
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)
