当前位置：首页 > news >正文

Eigen 中Sparse 模块的简单介绍和实战使用示例

news 2025/8/23 14:52:07

Eigen Sparse 模块详解

Eigen 的 稀疏矩阵（Sparse Matrix）模块 位于 #include <Eigen/Sparse>，主要用于存储和计算 大规模稀疏矩阵，避免使用 MatrixXd 那样的稠密存储浪费内存。

1. 存储结构

Eigen 默认采用 Compressed Column Storage (CCS)，也叫 列压缩格式，和 MATLAB/CSparse 类似。

稠密矩阵存储：MatrixXd 按行或列存储，每个元素都有。
稀疏矩阵存储：只存储 非零元及其行索引。

核心数据结构：

valuePtr[] ：非零元素值数组
innerIndexPtr[] ：非零元素对应的行索引
outerIndexPtr[] ：列起始位置索引

例如一个稀疏矩阵：

[ 1 0 0 2 ]
[ 0 3 0 0 ]
[ 4 0 5 6 ]

压缩存储：

values = [1,4,3,5,2,6]
rowIdx = [0,2,1,2,0,2]
colPtr = [0,2,3,5,6]   // 每一列的起始索引

2. 主要类

(1) `Eigen::SparseMatrix<T>`

用于存储稀疏矩阵，支持压缩和未压缩两种模式。

常见声明：

Eigen::SparseMatrix<double> A(m, n);  // m x n 稀疏矩阵

(2) `Eigen::SparseVector<T>`

一维稀疏向量，按 索引-值 方式存储。

(3) `Eigen::Triplet<T>`

用于 初始化稀疏矩阵 的三元组形式 (row, col, value)。
推荐：先收集 Triplet，再一次性构造 SparseMatrix，效率高。

(4) `Eigen::Map<SparseMatrix<T>>`

允许把外部内存映射为稀疏矩阵。

3. 构造和赋值

(1) 直接插入

Eigen::SparseMatrix<double> A(3,3);
A.insert(0,0) = 1.0;
A.insert(1,2) = 2.5;
A.insert(2,1) = -3.0;

插入时会触发动态分配，效率低，不推荐大量使用。

(2) 用 `Triplet` 批量初始化（推荐）

typedef Eigen::Triplet<double> T;
std::vector<T> tripletList;
tripletList.push_back(T(0,0,1.0));
tripletList.push_back(T(1,2,2.5));
tripletList.push_back(T(2,1,-3.0));Eigen::SparseMatrix<double> A(3,3);
A.setFromTriplets(tripletList.begin(), tripletList.end());

4. 基本操作

(1) 矩阵-向量乘法

Eigen::VectorXd x(3); x << 1,2,3;
Eigen::VectorXd y = A * x;

(2) 稀疏矩阵加法/乘法

Eigen::SparseMatrix<double> B(3,3);
Eigen::SparseMatrix<double> C = A + B;
Eigen::SparseMatrix<double> D = A * B;

(3) 转置

Eigen::SparseMatrix<double> At = A.transpose();

(4) 稀疏矩阵转稠密

Eigen::MatrixXd denseA = Eigen::MatrixXd(A);

(5) 稠密转稀疏

Eigen::MatrixXd M = Eigen::MatrixXd::Random(4,4);
Eigen::SparseMatrix<double> S = M.sparseView(); // 自动丢弃小元素

5. 稀疏求解器

Eigen 提供了 直接法 和 迭代法：

(1) 直接法（LU/Cholesky）

Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver;
solver.compute(A);          // A 必须对称正定
x = solver.solve(b);

(2) QR/LU 分解

Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
solver.compute(A);
x = solver.solve(b);

(3) 迭代法

共轭梯度 (ConjugateGradient) → 对称正定
BiCGSTAB → 一般稀疏矩阵
GMRES → 稀疏非对称系统

示例：

Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> cg;
cg.compute(A);
x = cg.solve(b);

6. 性能优化建议

批量构造：用 Triplet + setFromTriplets，避免频繁 insert。
压缩矩阵：A.makeCompressed()，提升运算效率。
避免拷贝：传递时尽量用引用。
数值稳定性：若矩阵条件数大，建议加预条件器（如 IncompleteCholesky）。
动态增长矩阵：如果要动态更新，先 A.reserve(nnz) 预分配空间。

7. 典型应用场景

大规模 SLAM / BA 后端优化（稀疏雅可比、Hessian）
PDE/有限元数值解
图优化问题（如 GTSAM、Ceres 后端里的稀疏矩阵运算）
图像分割 / 稀疏编码

8. 完整示例

#include <iostream>
#include <Eigen/Sparse>int main() {typedef Eigen::SparseMatrix<double> SpMat;typedef Eigen::Triplet<double> T;std::vector<T> tripletList;tripletList.push_back(T(0,0,4));tripletList.push_back(T(1,1,5));tripletList.push_back(T(2,2,6));tripletList.push_back(T(0,2,2));SpMat A(3,3);A.setFromTriplets(tripletList.begin(), tripletList.end());Eigen::VectorXd b(3);b << 1, 2, 3;Eigen::SimplicialLLT<SpMat> solver;solver.compute(A);Eigen::VectorXd x = solver.solve(b);std::cout << "Solution x = \n" << x << std::endl;
}

9.SLAM应用中的综合示例

一个完整的 Demo：手动构造 SLAM 后端的稀疏 H 矩阵（来自一组 2D 位置约束），并用 Eigen 稀疏求解器（SimplicialLLT）一次性解出增量。
示例问题：4 个 2D 位姿（只含位置，忽略旋转），带顺序里程计约束 + 1 个回环，外加对第一个位姿的先验以消除 gauge 自由度。

依赖：Eigen 3（#include <Eigen/Sparse>）。

// demo_slam_sparse.cpp
#include <Eigen/Sparse>
#include <iostream>
#include <vector>
#include <cmath>using SpMat = Eigen::SparseMatrix<double>;
using Triplet = Eigen::Triplet<double>;// 向 triplets 中添加 2x2 小块
void addBlock2x2(std::vector<Triplet>& T, int r, int c, const Eigen::Matrix2d& B) {T.emplace_back(r + 0, c + 0, B(0,0));T.emplace_back(r + 0, c + 1, B(0,1));T.emplace_back(r + 1, c + 0, B(1,0));T.emplace_back(r + 1, c + 1, B(1,1));
}// 将 2x1 向量累加到大 g 向量的对应位置
void addVec2(Eigen::VectorXd& g, int r, const Eigen::Vector2d& v) {g(r+0) += v(0);g(r+1) += v(1);
}int main() {// ========= 问题设置 =========// 4 个 2D 位姿（仅位置），真值为一个单位正方形的 4 个角点// x0=(0,0), x1=(1,0), x2=(1,1), x3=(0,1)const int N = 4;            // 位姿数量const int dim = 2;          // 每个位姿 2 维（x,y）const int nvar = N * dim;   // 总维度// 初始估计（故意加点偏差）std::vector<Eigen::Vector2d> x(N);x[0] = Eigen::Vector2d( 0.05, -0.02);x[1] = Eigen::Vector2d( 0.95,  0.10);x[2] = Eigen::Vector2d( 1.10,  1.05);x[3] = Eigen::Vector2d(-0.10,  0.95);// 观测（里程计 + 回环），模型: e = x_j - x_i - z_ijstruct Edge { int i,j; Eigen::Vector2d z; double w; }; // w=1/sigma^2std::vector<Edge> edges;// 顺序边：0-1: (1,0), 1-2: (0,1), 2-3: (-1,0), 3-0: (0,-1)（回环）double w = 100.0; // 信息权重（等价于 sigma=0.1，每维独立）edges.push_back({0,1, Eigen::Vector2d( 1.0,  0.0), w});edges.push_back({1,2, Eigen::Vector2d( 0.0,  1.0), w});edges.push_back({2,3, Eigen::Vector2d(-1.0,  0.0), w});edges.push_back({3,0, Eigen::Vector2d( 0.0, -1.0), w}); // 回环// 给 x0 加先验：x0_prior=(0,0)const Eigen::Vector2d x0_prior(0.0, 0.0);const double w_prior = 1e6; // 强先验，消除 gauge 自由度// ========= 构造正规方程 H dx = -g =========// 对于边 e_ij = xj - xi - z_ij，雅可比 J_i=-I, J_j=I// H 累加: H_ii += w*I, H_ij += -w*I, H_ji += -w*I, H_jj += w*I// g 累加: g_i += (-I)^T w e = -w*e, g_j += ( I)^T w e =  w*estd::vector<Triplet> triplets;Eigen::VectorXd g = Eigen::VectorXd::Zero(nvar);const Eigen::Matrix2d I2 = Eigen::Matrix2d::Identity();auto idx = [&](int k){ return k*dim; }; // 位姿 k 在大向量/矩阵中的起始索引for (const auto& e : edges) {int ri = idx(e.i);int rj = idx(e.j);// 计算残差 e = xj - xi - zEigen::Vector2d err = (x[e.j] - x[e.i]) - e.z;// H 的 4 个 2x2 blockaddBlock2x2(triplets, ri, ri,  e.w * I2);      // H_iiaddBlock2x2(triplets, ri, rj, -e.w * I2);      // H_ijaddBlock2x2(triplets, rj, ri, -e.w * I2);      // H_jiaddBlock2x2(triplets, rj, rj,  e.w * I2);      // H_jj// g 的两块addVec2(g, ri, -e.w * err); // g_iaddVec2(g, rj,  e.w * err); // g_j}// 先验: e_prior = x0 - x0_prior, J=I{int r0 = idx(0);addBlock2x2(triplets, r0, r0, w_prior * I2);                 // H_00 += wP * IEigen::Vector2d eprior = x[0] - x0_prior;addVec2(g, r0, w_prior * eprior);                            // g_0 += wP * (x0 - x0_prior)}// 组装稀疏矩阵 HSpMat H(nvar, nvar);H.setFromTriplets(triplets.begin(), triplets.end());H.makeCompressed();// RHS = -g  (Gauss-Newton: H dx = -g)Eigen::VectorXd rhs = -g;// ========= 求解 =========Eigen::SimplicialLLT<SpMat> solver; // H 应为对称正定solver.compute(H);if (solver.info() != Eigen::Success) {std::cerr << "Factorization failed.\n";return 1;}Eigen::VectorXd dx = solver.solve(rhs);if (solver.info() != Eigen::Success) {std::cerr << "Solve failed.\n";return 1;}// ========= 更新状态并打印 =========for (int k = 0; k < N; ++k) {x += dx(idx(k)+0);x += dx(idx(k)+1);}std::cout << "=== 解出的增量 dx ===\n" << dx.transpose() << "\n\n";std::cout << "=== 更新后的位姿（仅位置） ===\n";for (int k = 0; k < N; ++k) {std::cout << "x" << k << " =  << ", " << x << "]\n";}// 与真值比较std::vector<Eigen::Vector2d> gt = {{0,0},{1,0},{1,1},{0,1}};double rmse = 0.0;for (int k = 0; k < N; ++k) {rmse += (x[k] - gt[k]).squaredNorm();}rmse = std::sqrt(rmse / N);std::cout << "\nRMSE to ground truth = " << rmse << "\n";return 0;
}

代码细节要点

H 与 g 的装配：直接按 $\sum J^\top W J$ 、 $\sum J^\top W e$ 逐边累积；用 Triplet 高效构造稀疏矩阵。
先验：对 x0 加强先验（大权重）以固定参考系，避免奇异。
求解器：SimplicialLLT 适合 SPD；若矩阵可能不满足 SPD，可替换 SparseLDLT 或 SparseLU。
一次迭代即收敛：本例模型线性（雅可比常数），故一次 GN 迭代即可命中真值附近。