第 6 章 稀疏矩阵
稀疏矩阵操作
压缩稀疏行（格式）
求解稀疏线性系统
无矩阵求解器
稀疏矩阵： 涵盖稀疏矩阵操作、压缩稀疏行格式、求解稀疏线性系统以及矩阵自由求解器四个部分，每个部分都附有详细的中文注释说明，所有示例均使用 C++ 语言，并利用了 Eigen 库（当然也可以不用 Eigen，只是这里统一采用 Eigen 进行说明）。
6.1 稀疏矩阵操作
在本节中，我们介绍如何利用 Eigen 构造、初始化和遍历稀疏矩阵。Eigen 中的稀疏矩阵主要采用 Triplet 列表方式进行初始化，下面的例子中构造了一个 5×5 的稀疏矩阵，并通过迭代器输出其中的非零元素。
#include <iostream>
#include <vector>
#include <Eigen/Sparse> // 引入 Eigen 稀疏矩阵模块
usingnamespace Eigen;
usingnamespace std;
intmain(){
// 定义一个 5x5 的稀疏矩阵，元素类型为 double
SparseMatrix<double>spMat(5,5);
// 使用 Triplet 列表存储非零元素
typedef Triplet<double> T;
std::vector<T> tripletList;
// 向 tripletList 中添加非零元素，格式为 Triplet<行, 列, 数值>
tripletList.push_back(T(0,0,10));
tripletList.push_back(T(1,2,20));
tripletList.push_back(T(2,4,30));
tripletList.push_back(T(3,1,40));
tripletList.push_back(T(4,3,50));
// 利用 tripletList 构造稀疏矩阵
spMat.setFromTriplets(tripletList.begin(), tripletList.end());
// 遍历稀疏矩阵中所有的非零元素
// 注意：如果 spMat 是列主序（默认），outerSize() 返回的是列数；如果是行主序，则返回行数
for(int k =0; k < spMat.outerSize();++k){
// 内部迭代器遍历第 k 列（或第 k 行，取决于存储方式）上的所有非零元素
for(SparseMatrix<double>::InnerIterator it(spMat, k); it;++it){
// 输出元素的行号、列号以及对应的数值
cout <<"A("<< it.row()<<","<< it.col()<<") = "<< it.value()<< endl;
}
}
return0;
}6.2 压缩稀疏行（CSR）格式
压缩稀疏行格式（CSR）通常使用行主序存储。Eigen 默认是列主序，但我们可以通过指定存储方式为 RowMajor 来实现行主序存储。下面的示例展示如何构造一个行主序的稀疏矩阵，并访问其底层的 CSR 数据结构（行指针、列索引、非零值数组）。
#include <iostream>
#include <vector>
#include <Eigen/Sparse> // 引入稀疏矩阵模块
usingnamespace Eigen;
usingnamespace std;
intmain(){
// 定义一个 5x5 的行主序稀疏矩阵，RowMajor 表示按行存储（符合 CSR 格式）
SparseMatrix<double, RowMajor>spMat(5,5);
// 同样使用 Triplet 列表来初始化非零元素
typedef Triplet<double> T;
std::vector<T> tripletList;
tripletList.push_back(T(0,0,10));
tripletList.push_back(T(1,2,20));
tripletList.push_back(T(2,4,30));
tripletList.push_back(T(3,1,40));
tripletList.push_back(T(4,3,50));
spMat.setFromTriplets(tripletList.begin(), tripletList.end());
// Eigen 中行主序稀疏矩阵底层存储即为 CSR 格式，下面获取其内部数据：
// outerIndexPtr() 返回行指针数组（每一行非零元素在内存中的起始位置索引）
// innerIndexPtr() 返回列索引数组（存储每个非零元素的列索引）
// valuePtr() 返回非零元素数组（存储每个非零元素的具体数值）
cout <<"outerIndexPtr (行指针数组): ";
for(int i =0; i <= spMat.rows();++i){
cout << spMat.outerIndexPtr()[i]<<" ";
}
cout << endl;
cout <<"innerIndexPtr (列索引数组): ";
for(int i =0; i < spMat.nonZeros();++i){
cout << spMat.innerIndexPtr()[i]<<" ";
}
cout << endl;
cout <<"valuePtr (非零值数组): ";
for(int i =0; i < spMat.nonZeros();++i){
cout << spMat.valuePtr()[i]<<" ";
}
cout << endl;
return0;
}6.3 求解稀疏线性系统
在本节中，我们介绍如何使用 Eigen 的稀疏矩阵求解器来求解线性系统 Ax=b
（此处假设 A 为对称正定矩阵）。我们采用 SimplicialLLT 分解方法进行求解。
#include <iostream>
#include <vector>
#include <Eigen/Sparse> // 引入稀疏矩阵模块
#include <Eigen/SparseCholesky> // 引入稀疏 Cholesky 分解模块
usingnamespace Eigen;
usingnamespace std;
intmain(){
// 定义一个 4x4 的对称正定稀疏矩阵 A
SparseMatrix<double>A(4,4);
typedef Triplet<double> T;
std::vector<T> tripletList;
// 构造一个简单的对称正定矩阵 A：
// 例如： A = [ 4 -1 0 0
// -1 4 -1 0
// 0 -1 4 -1
// 0 0 -1 4 ]
tripletList.push_back(T(0,0,4));
tripletList.push_back(T(0,1,-1));
tripletList.push_back(T(1,0,-1));
tripletList.push_back(T(1,1,4));
tripletList.push_back(T(1,2,-1));
tripletList.push_back(T(2,1,-1));
tripletList.push_back(T(2,2,4));
tripletList.push_back(T(2,3,-1));
tripletList.push_back(T(3,2,-1));
tripletList.push_back(T(3,3,4));
// 利用 tripletList 初始化矩阵 A
A.setFromTriplets(tripletList.begin(), tripletList.end());
// 定义右侧向量 b
VectorXd b(4);
b <<15,10,10,15;
// 使用 SimplicialLLT 进行稀疏 Cholesky 分解，求解 Ax = b
SimplicialLLT<SparseMatrix<double>> solver;
solver.compute(A);
if(solver.info()!= Success){
cout <<"矩阵分解失败"<< endl;
return-1;
}
// 求解线性系统
VectorXd x = solver.solve(b);
if(solver.info()!= Success){
cout <<"求解失败"<< endl;
return-1;
}
// 输出求得的解向量 x
cout <<"解向量 x = \n"<< x << endl;
return0;
}6.4 矩阵自由求解器
矩阵自由求解器（Matrix Free Solvers）通常不需要显式存储整个矩阵，而是通过定义矩阵-向量乘法操作来实现迭代求解。下面给出一个共轭梯度（Conjugate Gradient, CG）方法的示例，其中矩阵 A 仅通过一个函数 A_mul 实现其与向量的乘法。
#include <iostream>
#include <Eigen/Dense> // 引入 Eigen Dense 模块用于向量操作
#include <cmath> // 用于数学函数 sqrt
usingnamespace Eigen;
usingnamespace std;
// 定义一个矩阵向量乘法函数，计算 y = A * x
// 本例中 A 定义为上述 4x4 矩阵
voidA_mul(const VectorXd &x, VectorXd &y){
int n = x.size();
y.resize(n);
y.setZero();// 将 y 初始化为零向量
for(int i =0; i < n;++i){
// 计算对角线元素的贡献
y(i)+=4*x(i);
// 如果存在左边元素，则添加对应的乘积（下三角元素）
if(i -1>=0){
y(i)+=-1*x(i -1);
}
// 如果存在右边元素，则添加对应的乘积（上三角元素）
if(i +1< n){
y(i)+=-1*x(i +1);
}
}
}
// 矩阵自由的共轭梯度求解器
// 参数：b 为右侧向量，tol 为收敛容限，maxIter 为最大迭代次数
VectorXd matrixFreeCG(const VectorXd &b,double tol =1e-6,int maxIter =1000){
int n = b.size();
// 初始猜测解设为零向量
VectorXd x =VectorXd::Zero(n);
// 初始残差 r = b - A*x，因 x = 0 所以 r = b
VectorXd r = b;
// 初始搜索方向设为 r
VectorXd p = r;
// 计算初始残差平方
double rsold = r.dot(r);
for(int i =0; i < maxIter;++i){
VectorXd Ap(n);
A_mul(p, Ap);// 计算 A 乘以搜索方向 p
// 计算步长 alpha = (r^T * r) / (p^T * A * p)
double alpha = rsold / p.dot(Ap);
// 更新解向量： x = x + alpha * p
x = x + alpha * p;
// 更新残差： r = r - alpha * A * p
r = r - alpha * Ap;
// 计算新的残差平方
double rsnew = r.dot(r);
// 检查是否收敛，如果残差的 2 范数小于 tol，则认为已收敛
if(sqrt(rsnew)< tol){
cout <<"共轭梯度法在 "<< i +1<<" 次迭代后收敛."<< endl;
break;
}
// 更新搜索方向： p = r + (rsnew/rsold) * p
p = r +(rsnew / rsold)* p;
// 更新 rsold
rsold = rsnew;
}
return x;
}
intmain(){
int n =4;
// 定义右侧向量 b
VectorXd b(n);
b <<15,10,10,15;
// 使用矩阵自由共轭梯度求解器求解 Ax = b
VectorXd x =matrixFreeCG(b);
// 输出求得的解向量 x
cout <<"解向量 x = \n"<< x << endl;
return0;
}说明
- 矩阵自由求解器思想
：与传统方法需要先构造并存储整个稀疏矩阵不同，此方法仅需定义一个“矩阵-向量乘法”操作函数
A_mul，适用于大规模问题或隐式定义的算子。 - 共轭梯度法
：适用于对称正定矩阵，在每次迭代中仅需一次矩阵向量乘法，非常适合矩阵自由的求解场景。
以上便是第6章关于稀疏矩阵的完整教程。每个代码示例都配有详细的中文注释，帮助读者理解稀疏矩阵的构造、存储格式以及如何利用不同的求解器求解稀疏线性系统。
扫一扫
4283
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
