C++ 数值计算与线性代数：Eigen 与表达式模板

写一个矩阵库，最直观的做法是重载 operator+ 和 operator*，让代码像数学公式一样：

Matrix C = A + B * C_inv;   // 看起来很美好

问题在于，operator* 会先算出一个临时矩阵，operator+ 再基于它算出第二个临时矩阵。对一个 1000×1000 的 double 矩阵，每个临时对象就是 8MB，一次表达式产生数倍的内存分配与内存带宽消耗。数值计算是内存带宽受限的典型场景，多一次「写临时 + 读临时」就多两倍带宽，性能直接腰斩。

Eigen 的答案是把整个表达式延迟到赋值时一次性求值，用表达式模板（Expression Templates）在编译期把 A + B * C 组合成一棵表达式树，遍历时每个元素只计算一次，中间结果留在寄存器里。这是 C++ 模板元编程在数值计算领域最成功的应用之一。本文从这个问题出发，讲清 Eigen 的机制、选型与陷阱。

朴素运算符重载的代价

先量化问题。朴素实现下，D = A + B * C（均为 1000×1000）的执行序列是：

1. tmp1 = B * C：1000³ = 10⁹ 次乘加，写 8MB，读 16MB。
2. tmp1 再从内存读回，与 A 相加，写 8MB，读 16MB。
3. D = tmp1：又读 8MB、写 8MB。

总内存流量约 64MB，其中一半以上是临时对象造成的。如果用表达式模板融合成单次遍历：

for (i, j) D[i][j] = A[i][j] + dot(B.row(i), C.col(j));

中间结果全部在寄存器/栈上，内存流量降到「读 A、读 B、读 C、写 D」的必需量。矩阵乘法本身仍要 10⁹ 次乘加，但内存带宽压力大幅下降——对于大矩阵，这正是瓶颈所在。

这就是 Eigen 的核心价值：把「写起来像数学」与「跑起来无浪费」这两件通常矛盾的事统一起来。

Eigen 基础：Matrix、Vector 与尺寸选择

Eigen 的核心类型是 Matrix<Scalar, Rows, Cols>，Vector 是列数为 1 的特化。

#include <Eigen/Dense>
using namespace Eigen;

MatrixXd A(3, 3);          // 动态尺寸 double 矩阵
A << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;              // 逗号初始化

Vector3d v(1.0, 2.0, 3.0); // 固定尺寸 3 维向量
Matrix4f M = Matrix4f::Identity();   // 4x4 float 单位阵

尺寸参数的写法：

固定尺寸（编译期已知）与动态尺寸的性能差异显著：固定尺寸的矩阵存放在栈上，尺寸是编译期常量，循环可以完全展开、向量化，且没有堆分配。在图形学里，4×4 变换矩阵用 Matrix4f 比 MatrixXf 快得多，也避免了每帧数百次的小堆分配。

选择原则：

• 尺寸 ≤ 4 且固定 → 用固定尺寸（Matrix4f、Vector3d）。
• 尺寸较大或运行时决定 → 用动态尺寸（MatrixXd）。
• 尺寸不大但运行时才知（如 10~50） → 可以用 Matrix<double, N, N> 配合模板，或直接动态尺寸。

一个常见的错误是「所有矩阵都用 MatrixXd」，包括 3 维向量。这会让本该零开销的小向量运算变成堆分配。

存储顺序：列优先还是行优先

Eigen 默认列优先（column-major），与 Fortran、BLAS、LAPACK 一致。这对性能有实际影响：

Matrix<double, Dynamic, Dynamic, RowMajor> R(3, 3);   // 行优先

选择的关键是与访问模式匹配：如果算法主要按行遍历，行优先能获得更好的缓存局部性。但混用两种顺序会带来转换开销，一个项目内应尽量统一。与外部库（如 OpenCV 的 Mat、PyTorch 的张量）对接时，顺序不匹配会导致隐式的转置拷贝，往往是性能黑洞。

表达式模板：编译期的表达式树

Eigen 的表达式模板机制值得拆开看。当写 D = A + B 时，A + B 并不计算，而是返回一个类型 CwiseBinaryOp<scalar_sum_op<double>, MatrixXd, MatrixXd>——一个表达式对象，它持有 A、B 的引用，但没有任何数据。

auto expr = A + B;         // expr 的类型是表达式模板，不是 MatrixXd
// expr 里存的是对 A、B 的 const 引用，没有拷贝
MatrixXd D = expr;         // 赋值时才真正逐元素计算

多个运算嵌套时，表达式模板层层组合：

D = A + B * C;   // 类型形如 CwiseBinaryOp<sum, MatrixXd, Product<MatrixXd, MatrixXd>>

赋值给 MatrixXd 时，Eigen 遍历 D 的每个元素，对每个 (i,j) 计算整个表达式。对于 A + B*C，这不能融合成单循环（矩阵乘法需要内层求和），Eigen 会在内部展开为「先算矩阵乘法，再逐元素加」，但乘法的结果仍可能被优化掉中间存储——这正是 Eigen 在 noalias 与产品求值策略上的复杂之处。

对于逐元素运算（+、-、.*、sin、exp），表达式模板能做到完全的循环融合，这是它收益最大的场景：

// 单次遍历，无临时对象
D = (A.array() + 1.0).sin() + B.array().square();

这里的 .array() 把矩阵切换成「逐元素」语义（数组运算），sin、square 都是逐元素函数。矩阵语义与数组语义的区别是 Eigen 的一个关键设计：A * B 是矩阵乘法，A.array() * B.array() 是逐元素乘法。混用会得到意料之外的结果，.array() 与 .matrix() 是两者之间的转换。

表达式模板的本质是用类型编码计算结构，这与 https://plumephp.com/cpp-metaprogramming/ 里讨论的编译期计算同源：类型系统成了「计算图」的载体，编译器成了「求值器」。理解这一点，就能理解为什么 Eigen 的模板错误信息如此冗长——类型名里包含了整棵表达式树。

陷阱：aliasing、auto 与 eval()

表达式模板带来性能，也带来一串反直觉的陷阱。

陷阱一：aliasing（别名）

当赋值的目标与源出现在同一表达式里，且运算不是逐元素独立时，就地计算会读到已被覆盖的数据。

MatrixXd A = MatrixXd::Random(3, 3);
A = A.transpose();          // 危险！转置赋值给自己
A = A * A;                  // 危险！矩阵乘法不能就地

A = A.transpose() 在 Eigen 里会被自动检测并处理（Eigen 对已知的别名情况会先求值到临时变量），但 A = A * A 这类必须显式加 .eval() 或 noalias()：

A = (A * A).eval();         // 强制先算到临时，再赋值

反过来，当确定没有别名时可以显式 noalias() 消除临时：

C.noalias() = A * B;        // 告诉 Eigen：C 与 A、B 无重叠，可直接写入

noalias() 能省掉一次临时矩阵的分配与拷贝，但用错会得到静默的错误结果。只有在确认 C 与 A、B 不共享内存时才用。

陷阱二：auto 与表达式类型

auto 推导出的是表达式类型，不是矩阵类型：

auto expr = A + B;          // 类型是表达式对象，持有 A、B 的引用
MatrixXd D = expr;          // 此时才求值
// 若 A、B 在此之后被修改或销毁，expr 会悬垂或读到新值

auto bad = (A + B).eval();  // 正确：立刻求值成矩阵
auto worse = A.block(0, 0, 2, 2) + B.block(0, 0, 2, 2);  // 持有引用，注意生命周期

在函数返回值或长生命周期变量上用 auto 接表达式，是 Eigen 最隐蔽的 bug 来源。规则很简单：表达式只应作为临时值立即使用；要存储就用具体矩阵类型或 .eval()。

// 反例：返回表达式对象，引用局部变量 → 悬垂
auto make() { MatrixXd A(2,2); return A + A; }   // 返回表达式，A 已销毁！

// 正例
MatrixXd make() { MatrixXd A(2,2); return A + A; }  // 返回类型是 MatrixXd，会求值

陷阱三：跨表达式的隐式求值时机

Eigen 的求值策略（是否融合、是否分块）由内部启发式决定，不完全可预测。当性能不符合预期时，用 .eval() 或 .noalias() 显式控制：

// 大矩阵乘积，避免临时
C.noalias() = A * B;

// 逐元素链，通常无需干预，但可显式求值观察差异
D = (A.array() + B.array()).eval();

性能调优的第一步永远是测量。Eigen 本身有不错的向量化与分块，多数情况下不需要手工干预，但矩阵乘法这类核心运算值得对照 BLAS 基准验证。

分解与求解：选对算法比优化代码更重要

解线性方程组 Ax = b 时，直接求逆是最慢且最不数值稳定的做法。Eigen 提供一系列分解（decomposition），应按矩阵性质选择。

// 求解 Ax = b（推荐：用分解，不要用 A.inverse()）
VectorXd x = A.colPivHouseholderQr().solve(b);

// 或显式分解，复用于多个右端项
ColPivHouseholderQR<MatrixXd> qr(A);
VectorXd x1 = qr.solve(b1);
VectorXd x2 = qr.solve(b2);

选型表：

几条实用规则：

• 已知对称正定（如求解线性系统里的刚度矩阵、协方差）：用 LLT，速度是 PartialPivLU 的两倍。
• 最小二乘 min ||Ax - b||：用 ColPivHouseholderQR 或 HouseholderQR，不要用 A.transpose() * A 再解正规方程——那会把条件数平方，精度大幅下降。
• 求伪逆或判断秩：用 JacobiSVD（小矩阵）或 BDCSVD（大矩阵）。
• 不要 A.inverse()：求逆再乘是「更慢、更不稳、更不必要」的三重错误。

// 反例
VectorXd x = A.inverse() * b;          // 慢且不稳

// 正例
VectorXd x = A.partialPivLu().solve(b);

稀疏矩阵

当矩阵大部分元素为零（如有限元、图、推荐系统），必须用 SparseMatrix，否则内存与计算量都是灾难：

#include <Eigen/Sparse>
SparseMatrix<double> S(10000, 10000);
S.reserve(VectorXi::Constant(10000, 10));   // 预估每行非零元个数
S.insert(i, j) = value;                     // 或 setFromTriplets 批量构建
S.makeCompressed();

// 稀疏求解
SparseLU<SparseMatrix<double>> solver;
solver.compute(S);
VectorXd x = solver.solve(b);

稀疏矩阵的关键是存储格式（压缩稀疏列 CSC）、填充（fill-in）与排序（ordering）。reserve 预估非零元个数能避免反复重分配；setFromTriplets 适合从「三元组列表」批量构建。填充最小化（如 AMD、COLAMD 排序）对稀疏 Cholesky/LU 的性能影响巨大，Eigen 支持传入自定义排序。稀疏线性代数的完整方法学可以参看 HPC 稀疏线性代数 。

与外部内存互操作：Map 与 Ref

实际项目里，矩阵数据常常已经存在于别处：一段 C 数组、一个 std::vector、GPU 传来的缓冲区、另一个库的 Mat。Eigen 的 Map 能在不拷贝的前提下把这些内存当作 Eigen 矩阵使用。

double raw[9] = {1,2,3, 4,5,6, 7,8,9};
// 把 raw 解释为 3x3 列优先矩阵，零拷贝
Map<Matrix<double, 3, 3>> M(raw);

// 把 vector 当作动态矩阵
std::vector<double> buf(12);
Map<MatrixXd> N(buf.data(), 3, 4);   // 注意是列优先解释

Map 的关键参数是步长（stride）与存储顺序：

// 行优先解释
Map<Matrix<double, Dynamic, Dynamic, RowMajor>> R(raw, 3, 3);

// 带行步长的映射（数据不是紧密排列）
Map<MatrixXd, 0, Stride<Dynamic, 1>> S(ptr, rows, cols, Stride<Dynamic,1>(ld, 1));

陷阱在于 Map 默认按列优先解释，而 C 数组通常是行优先。把行优先的 C 数组直接 Map<MatrixXd> 会得到「转置后」的矩阵——数值上不报错，结果全错。对接外部数据时，务必确认对方的内存布局，必要时显式指定 RowMajor。

Ref<> 则是函数参数的最佳实践：它接受任何「行为像矩阵」的表达式，同时避免对临时对象取引用：

// 好：接受 MatrixXd、Map、block 等，且不会意外绑定临时
void normalize(Ref<VectorXd> v) {
    v /= v.norm();
}

// 差：只接受 VectorXd，传 block 会触发拷贝
void normalize_bad(VectorXd v);

用 Ref 而不是 const MatrixXd& 作为参数，能避免「传一个 block 或表达式时隐式构造临时矩阵」的开销。这是 Eigen 代码在 API 设计上的标准做法。

数值稳定性与精度

数值计算的 bug 往往不是崩溃，而是「结果看起来对但精度不够」。几个必须知道的稳定性问题。

条件数（condition number） 衡量问题对扰动的敏感度。条件数为 10¹⁵ 的矩阵，双精度（约 16 位有效数字）求解后可能一位有效数字都不剩。判断方法：

JacobiSVD<MatrixXd> svd(A);
double cond = svd.singularValues()(0) / svd.singularValues()(svd.singularValues().size()-1);

不要解正规方程。 求解最小二乘时，(AᵀA)x = Aᵀb 会把条件数平方——A 条件数 10⁸ 时，AᵀA 变成 10¹⁶，双精度直接失效。正确做法是用 QR 或 SVD：

// 反例：条件数被平方
VectorXd x = (A.transpose() * A).ldlt().solve(A.transpose() * b);

// 正例：QR 直接处理
VectorXd x = A.colPivHouseholderQr().solve(b);

避免「大数相减」。 两个相近的大数相减会造成灾难性抵消（catastrophic cancellation）。例如用 1 - cos(x) 计算小角度时，x 接近 0 会丢失全部有效数字，应改用 2 * sin(x/2)^2 或 expm1 这类数值稳定的形式。Eigen 提供了 stableNorm() 用于避免向量范数计算中的溢出与下溢。

累加顺序影响结果。 浮点加法不满足结合律，(a+b)+c 与 a+(b+c) 可能不同。大规模求和应使用 Kahan 求和或分块求和。Eigen 的 sum() 内部做了分块，比朴素循环更稳，但跨线程的归约顺序仍可能带来非确定性——并行数值程序的「结果不可复现」大多源于此。

单精度 vs 双精度 是显式的取舍：float 快一倍、省一半内存，但只有约 7 位有效数字。机器学习的推理常用 float 甚至 half，而传统科学计算几乎都用 double。选错精度的代价，通常远大于任何代码优化。

与 BLAS/LAPACK、SIMD 和多线程的协同

Eigen 的定位是「模板库」，它能把部分运算委托给高度优化的 BLAS/LAPACK，也能利用 SIMD 与多线程。

SIMD 向量化：Eigen 对逐元素运算与部分归约默认启用向量化（SSE/AVX/NEON）。启用方式取决于编译器选项，并可通过 EIGEN_MAX_ALIGN_* 控制对齐。核心矩阵乘法（GEMM）的向量化与分块策略是性能关键，与 https://plumephp.com/cpp-simd-vectorization-practice/ 里讨论的手工向量化思路一致——只是 Eigen 把它自动化了。

多线程：Eigen 默认单线程，开启并行需要定义宏并链接线程库：

#define EIGEN_DONT_PARALLELIZE       // 显式关闭
// 或编译时加 -fopenmp，Eigen 会在大型矩阵运算上并行

并行粒度由 Eigen 内部决定，通常对足够大的矩阵才启动。嵌套并行（外层 OpenMP 循环 + Eigen 内部并行）会导致线程超额订阅，反而变慢——多线程程序里用 Eigen 要小心这一点。

对接 BLAS/LAPACK：Eigen 可以链接 Intel MKL 或 OpenBLAS，把 A * B 这类大矩阵乘法交给厂商优化的 GEMM：

find_package(BLAS REQUIRED)
target_link_libraries(app PRIVATE Eigen3::Eigen ${BLAS_LIBRARIES})
# 或在编译时定义 EIGEN_USE_BLAS / EIGEN_USE_MKL_ALL

对大型稠密矩阵乘法，MKL 通常比 Eigen 自带的 GEMM 快（更好的分块与指令调度）；对中小矩阵与逐元素运算，Eigen 的自有实现往往更好（省去函数调用与边界检查开销）。选择依据是矩阵尺寸，不要一刀切。

如果计算要上 GPU，Eigen 的 Tensor 模块与 cuBLAS/cuSolver 可以协同，但要理解 CPU/GPU 之间的传输成本。GPU 加速的收益分析、profiling 方法与踩坑经验，可以对照 GPU profiling 与 Nsight 和 分布式训练 中关于大规模数值计算的讨论。整体性能优化方法论则与 https://plumephp.com/cpp-performance-optimization/ 相通：先测量、找瓶颈、再优化。

实践建议

1. 小尺寸用固定尺寸类型（Matrix4f、Vector3d），避免堆分配并启用完整展开。
2. 矩阵用 *，逐元素用 .array() *，不要混用两种语义。
3. 不要 A.inverse()，用分解的 .solve()；对称正定用 LLT，最小二乘用 QR。
4. 注意 aliasing：A = A.transpose()、A = A * A 需 .eval()；确定无重叠才用 .noalias()。
5. auto 不接表达式，要存储就用具体矩阵类型或 .eval()；表达式只作临时值。
6. 稀疏矩阵必须 reserve + makeCompressed，并关注填充与排序。
7. 存储顺序与访问模式匹配，与外部库对接时确认列/行优先是否一致。
8. 大矩阵乘法考虑接 BLAS，小矩阵留给 Eigen 自有实现。
9. 多线程程序里避免嵌套并行，Eigen 的内部并行与外层 OpenMP 会互相打架。
10. 一切以测量为准：先 profiling 找到真正的瓶颈，再决定是改算法还是改实现。
11. 接口参数用 Ref<> 而非 const MatrixXd&，避免调用方传入 block 或表达式时被隐式拷贝成临时矩阵。
12. 固定尺寸的小矩阵优先放栈上，既省去堆分配，也让编译器有机会完全展开循环并向量化。

Eigen 的设计哲学值得单独记一笔：它用模板元编程把「表达式的结构」编码进类型，从而在编译期完成循环融合与临时消除。这既是它高性能的根源，也是它编译错误难读、auto 行为反直觉的根源。用好 Eigen 的关键，是同时理解数学层面（选哪个分解、条件数多大）与类型层面（表达式何时求值、内存何时分配）这两件事——只懂前者会写出慢代码，只懂后者会写出错代码。