Eigen 完全教程：从入门到精通

版本: Eigen 5.0.1
难度: 初级到高级
目标读者: C++ 开发者、科研人员、工程师
编译标准: Eigen 5.0.x 最低要求 C++14，本教程示例默认按 C++17 组织

快速导航

• 初学者：按 安装篇 → 基础篇 → 矩阵基础篇 → 高级操作篇 → 线性代数篇 的顺序学习
• 有经验者：可直接跳转到 线性代数篇、几何变换篇、稀疏矩阵篇 等专题
• 迁移用户：重点阅读附录：迁移变更汇总，了解 Eigen 5.0/5.0.1 的兼容性变化与推荐写法

关于本教程

本教程旨在为 C++ 开发者、科研人员和工程师提供一份全面、系统、实用的 Eigen 学习指南。从基础的矩阵操作到高级的性能优化，从线性代数算法到几何变换应用，本教程涵盖了 Eigen 库的核心功能和最佳实践。

教程特色

• 循序渐进：从基础概念到高级应用，适合不同水平的读者
• 实战导向：每个章节都包含丰富的代码示例和实战案例
• 版本更新：针对 Eigen 5.0 的最新特性进行更新
• 问题导向：提供常见问题解答和调试技巧

参考文档

Eigen 官方文档

一、安装篇：环境配置

1.1 快速安装

Eigen 是纯头文件库，安装完成后无需额外编译库文件。

# Linux
sudo apt-get install libeigen3-dev

# macOS
brew install eigen

# Windows (vcpkg)
vcpkg install eigen3

# 或下载源码
git clone https://gitlab.com/libeigen/eigen.git
cd eigen && git checkout 5.0.1

安装完成后的最小验证

安装完成后先运行一个最小示例，确认头文件路径和编译选项正确。

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Matrix2d A;
    A << 1, 2,
         3, 4;

    std::cout << "A =\n" << A << "\n";
    std::cout << "det(A) = " << A.determinant() << "\n";
    return 0;
}

# Linux/macOS
g++ -std=c++17 -I/path/to/eigen -O2 -o eigen_check eigen_check.cpp

# Windows (MinGW)
g++ -std=c++17 -I"C:\path\to\eigen" -O2 -o eigen_check.exe eigen_check.cpp

预期现象：

• 程序成功编译
• 终端输出矩阵 A
• 输出 det(A) = -2

常见失败原因：

• -I 路径没有指向 Eigen 头文件根目录
• C++ 标准过低（Eigen 5.0.x 最低要求 C++14）
• Windows 下路径或引号写法不正确

1.2 CMake 集成

下面给出两个完整可运行的最小 CMake 示例。

方式1：系统已安装 Eigen

cmake_minimum_required(VERSION 3.16)
project(MyProject LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

add_executable(myapp main.cpp)

find_package(Eigen3 5.0 REQUIRED NO_MODULE)
target_link_libraries(myapp PRIVATE Eigen3::Eigen)

注意：find_package(Eigen3 ...) 的可用性取决于系统是否正确安装了 Eigen 5.x 的 CMake 配置文件。并非所有发行包都提供该配置；若此方式失败，可使用下面的 FetchContent 方式。

方式2：使用 FetchContent 自动下载

cmake_minimum_required(VERSION 3.16)
project(MyProject LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

include(FetchContent)
FetchContent_Declare(
    eigen
    GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
    GIT_TAG 5.0.1
)
FetchContent_MakeAvailable(eigen)

add_executable(myapp main.cpp)
target_link_libraries(myapp PRIVATE Eigen3::Eigen)

对应的最小 main.cpp

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Vector3d v(1, 2, 3);
    std::cout << "v = " << v.transpose() << "\n";
    return 0;
}

构建命令

cmake -S . -B build
cmake --build build
./build/myapp

命令行直接编译示例：

g++ -std=c++17 -I/path/to/eigen -O3 -march=native -o myapp myapp.cpp

1.3 头文件模块说明

选择指南：

• 最小依赖：仅使用<Eigen/Core>
• 通用场景：使用<Eigen/Dense>（包含大部分功能）
• 几何应用：使用<Eigen/Geometry>
• 大规模计算：使用<Eigen/Sparse>
• 单独功能：按需包含特定模块（如<Eigen/QR>）

对应官方文档：Getting Started

二、基础篇：Eigen简介与快速入门

覆盖 Eigen 的整体概览、最小示例、命名空间使用规范。

完成安装验证后，继续阅读：

• 矩阵基础篇：矩阵/向量声明、初始化、运算与存储
• 高级操作篇：块操作、Array、Map、Ref 等进阶主题

2.1 什么是Eigen？

Eigen 是一个高性能的C++模板库，专注于线性代数、矩阵和向量运算。核心特点：

Matrix vs Array：Eigen提供两种核心类型

• Matrix：用于线性代数运算（矩阵乘法、分解、求解线性方程等）
• Array：用于逐元素运算（sin、exp、逐元素乘法等）

A * B 在 Matrix 语义下表示矩阵乘法，在 Array 语义下表示逐元素乘法。

2.2 快速入门示例

#include <iostream>
#include <Eigen/Dense>

int main() {
    // 使用固定大小的矩阵（尺寸在编译期确定，通常更容易被优化）
    Eigen::Matrix3d A;
    A << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;
    
    // 使用动态大小矩阵（尺寸在运行时确定）
    Eigen::MatrixXd B(3, 3);
    B.setRandom();  // 随机初始化

    // 已知尺寸时，也可继续使用固定大小矩阵
    Eigen::Matrix3d B_fixed = Eigen::Matrix3d::Random();
    
    // 向量操作
    Eigen::Vector3d v(1, 2, 3);
    
    // 矩阵运算
    Eigen::MatrixXd C = A * B;        // 固定大小 × 动态大小 -> 结果也更自然地写成动态矩阵
    Eigen::Matrix3d C_fixed = A * B_fixed; // 固定大小 × 固定大小
    Eigen::Vector3d y = A * v;        // 矩阵-向量乘法
    double dot = v.dot(v);            // 点积
    
    std::cout << "A =\n" << A << "\n\n";
    std::cout << "A * B 的尺寸: " << C.rows() << " x " << C.cols() << "\n";
    std::cout << "A * v =\n" << y << "\n";
    std::cout << "v · v = " << dot << "\n";
    
    return 0;
}

编译命令

未完成环境配置的，先完成 安装篇 的最小验证。

# Linux/macOS
g++ -std=c++17 -I/path/to/eigen -O2 -o eigen_intro eigen_intro.cpp

# Windows (MinGW)
g++ -std=c++17 -I"C:\eigen" -O2 -o eigen_intro.exe eigen_intro.cpp

2.3 命名空间使用指南

// 方式1：完整命名空间（推荐，避免命名冲突）
Eigen::Matrix3d A;
Eigen::VectorXd v;

// 方式2：函数内局部using（推荐）
void my_function() {
    using namespace Eigen;
    Matrix3d A;  // 仅在此函数内有效
}

// 方式3：使用using声明
using Eigen::Matrix3d;
Matrix3d A;

最佳实践：头文件中必须使用完整命名空间，源文件中可使用局部 using namespace。示例代码中用 using namespace Eigen; 的情况均为示例场景，不照搬到公共头文件中。

后续章节会逐步引入更多类型、表达式和分解器，保持类型和命名空间书写清晰有助于减少混淆。

常见问题 (FAQ)

Q: Eigen与MATLAB相比有什么优势？

A:

• 性能: C++编译代码执行效率远高于MATLAB解释执行
• 部署: 可生成独立可执行文件，无需运行时环境
• 集成: 易于集成到现有C++项目中
• 成本: 完全开源免费（MPL2许可证）

Q: 固定大小矩阵为什么比动态大小矩阵更容易被优化？

A: 主要有两个原因：

1. 栈分配 vs 堆分配：固定大小矩阵（如 Matrix4d）的数据是对象内部的静态数组，直接在栈上分配，没有堆内存分配开销；动态大小矩阵（如 MatrixXd）在运行时分配堆内存。
2. 向量化与对齐：当固定大小矩阵的尺寸是 16 字节的整数倍时（如 Vector2d、Matrix4d），Eigen 申请 SIMD 要求的对齐（16/32/64 字节），利用对齐的 SSE/AVX 指令高效读写。这类为“固定大小可向量化“（fixed-size vectorizable）类型。

固定大小可向量化类型在 new 分配、STL 容器、按值传递等场景下需要注意内存对齐问题。详见 调试与排错篇 9.2 节。

Q: Eigen支持GPU计算吗？

A: Eigen 很早就支持在 CUDA 内核中使用固定大小的矩阵、向量和数组类型；这并不是 5.0 才新增的能力。对于动态大小矩阵，通常更适合结合 cuBLAS 等库，或采用 CPU/GPU 混合计算方案。

练习题

1. 基础练习: 创建一个3×3单位矩阵，计算其行列式和逆矩阵。
2. 思考练习: 解释为什么固定大小矩阵(Matrix3d)通常比动态大小矩阵(MatrixXd)更容易被优化？（提示：从内存分配和向量化两个角度思考）
3. 衔接练习: 思考下面三个概念分别会在后续哪一章深入展开：
   • Matrix 与 Array 的区别
   • 可逆矩阵与线性方程组求解
   • Map 与 Ref 的用途

小结

• Eigen 是一个面向线性代数的高性能 C++ 模板库
• 支持固定大小和动态大小对象
• Matrix 与 Array 的语义不同，不能混为一谈
• 教程示例默认使用 C++17
• 头文件中应避免 using namespace Eigen

下一步：矩阵基础篇

对应官方文档：Getting Started

三、矩阵基础篇：声明与基本运算

前置：已完成安装，了解 Matrix 与 Array 的基本区别。

3.1 矩阵和向量的声明

Eigen 中，矩阵表示线性变换、系数表或二维数据；列向量表示坐标、状态、参数或观测值；行向量常用于一行数据、转置结果或统计场景。

后续章节默认以列向量作为主要约定；如果没有特别说明，VectorXd 一般指动态大小的列向量。

Eigen提供了丰富的类型别名，简化矩阵和向量的声明。

类型命名规则

Matrix[尺寸][数据类型]

尺寸：X = 动态大小, N = 固定大小 N (如 2, 3, 4)
数据类型：d = double, f = float, i = int, cd = complex<double>

常用类型速查表

声明示例

// 固定大小矩阵（尺寸在编译期确定，通常更容易被优化）
Eigen::Matrix3d A;           // 3×3双精度矩阵
Eigen::Vector3d v1;          // 3维列向量

// 动态大小矩阵（尺寸在运行时确定）
Eigen::MatrixXd D(10, 10);   // 10×10双精度矩阵
Eigen::VectorXd v2(100);     // 100维动态列向量

// 特殊矩阵
Eigen::Matrix3d I = Eigen::Matrix3d::Identity();  // 单位矩阵
Eigen::MatrixXd Z = Eigen::MatrixXd::Zero(5, 5);  // 零矩阵
Eigen::MatrixXd R = Eigen::MatrixXd::Random(3, 3);// 随机矩阵

3.2 矩阵初始化方法

逗号初始化适合小型固定大小对象；循环赋值适合按规则填充；预定义函数适合零矩阵、单位矩阵、常量矩阵、随机矩阵等场景。

// 方法1：逗号初始化
Eigen::Matrix3d A;
A << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;

// 方法2：逐个元素赋值
Eigen::MatrixXd B(3, 3);
for (int i = 0; i < 3; ++i)
    for (int j = 0; j < 3; ++j)
        B(i, j) = i * 3 + j + 1;

// 方法3：预定义函数
Eigen::Matrix3d C = Eigen::Matrix3d::Zero();      // 全零
Eigen::Matrix3d D = Eigen::Matrix3d::Identity();  // 单位矩阵
Eigen::Matrix3d E = Eigen::Matrix3d::Constant(5); // 全为5
Eigen::Matrix3d F = Eigen::Matrix3d::Random();    // 随机值

3.3 基础矩阵运算

A * B 是矩阵乘法，A.array() * B.array() 是逐元素乘法。涉及 sqrt、exp、逐元素乘除法时，使用 Array 视图。

Eigen::Matrix3d A, B;
A << 2, -1, 0,
     -1, 2, -1,
     0, -1, 2;                    // 可逆矩阵
B << 9, 8, 7,
     6, 5, 4,
     3, 2, 1;

// 算术运算
Eigen::Matrix3d C = A + B;        // 矩阵加法
Eigen::Matrix3d AB = A * B;       // 矩阵乘法
Eigen::Matrix3d scaled = A * 2.0; // 数乘

// 转置
Eigen::Matrix3d At = A.transpose();

// 矩阵属性
double det = A.determinant();       // 行列式
double tr = A.trace();              // 迹
double norm = A.norm();             // Frobenius范数
Eigen::Matrix3d Ainv = A.inverse(); // 逆矩阵（仅对可逆方阵有意义）

// 逐元素运算（需转换为Array）
Eigen::Matrix3d coeff_product = A.array() * B.array(); // 逐元素乘法

// 逐元素平方根要求元素非负；这里单独构造一个满足前提的示例矩阵
Eigen::Matrix3d P;
P << 1, 4, 9,
     16, 25, 36,
     49, 64, 81;
Eigen::Matrix3d sqrtP = P.array().sqrt();

3.4 向量运算

dot() 是点积，返回标量；cross() 限于三维向量；normalized() 返回新向量（不修改原对象），normalize() 原地修改。

Eigen::Vector3d v1(1, 2, 3);
Eigen::Vector3d v2(4, 5, 6);

// 基本运算
double dot = v1.dot(v2);          // 点积
Eigen::Vector3d cross = v1.cross(v2);  // 叉积（仅3D向量）
double norm = v1.norm();          // 欧几里得范数

// normalized() 返回新的归一化向量，不修改 v1 本身
Eigen::Vector3d normalized = v1.normalized();

// 向量元素运算（此时 v1 仍是原始值 {1, 2, 3}）
v1.sum();                         // 所有元素之和
v1.prod();                        // 所有元素之积
v1.mean();                        // 平均值
v1.minCoeff();                    // 最小值
v1.maxCoeff();                    // 最大值

// normalize() 直接修改原向量
v1.normalize();  // 此后 v1 变为单位向量

3.5 数据类型与存储

block、Map、Ref、切片等功能与对象尺寸是否已知、内存是否连续、存储顺序密切相关。

• 小矩阵优先用固定大小类型
• 尺寸运行时确定时用动态大小类型
• Eigen 默认列优先（ColMajor）

固定大小 vs 动态大小

存储顺序

// Eigen默认列优先存储（与MATLAB、Fortran一致）
Eigen::Matrix<double, 3, 4, Eigen::ColMajor> A;  // 列优先（默认）
Eigen::Matrix<double, 3, 4, Eigen::RowMajor> B;  // 行优先（与C一致）

// 访问元素
A(i, j);  // 第i行第j列（从0开始）

对应官方文档：The Matrix class | Matrix and vector arithmetic

四、高级操作篇：块操作与内存映射

本篇介绍块操作、Array、Map 和 Ref，关注表达式语义、内存布局、视图与拷贝。

重点：

1. 区分取数据（视图）和复制数据（拷贝）
2. 行向量与列向量的类型差异
3. Map 与 Ref 的适用场景
4. 零拷贝写法与产生临时对象的写法

4.1 块操作与子矩阵

4.1.1 先理解：视图(view) 与拷贝(copy)

很多“取子矩阵”的操作默认返回的是表达式视图，而不是立即复制一份新矩阵。
这里要把“表达式”和“视图”区分开来看：

• block()、row()、col()、head()、segment() 这类 API，通常返回的是引用原对象数据的视图表达式
• 赋值给 MatrixXd、VectorXd 等具体类型时发生拷贝
• 修改视图会影响原矩阵；修改拷贝不会影响原矩阵
• 但并不是所有 Eigen 表达式都等于“原地引用原数据”；有些表达式只是延迟求值的组合表达式

初学时可以先把这类块操作理解为“默认是视图，不是副本”，但不要把所有“表达式对象”都等同理解成“直接别名”。

Eigen::MatrixXd A(6, 6);
A.setRandom();

// 方式1：复制出一个新的矩阵
Eigen::MatrixXd block_copy = A.block(1, 1, 3, 3);

// 方式2：直接操作原矩阵中的那一块
A.block(0, 0, 2, 2) = Eigen::Matrix2d::Identity();

规则：

• MatrixXd x = ...; 往往表示“我要一个独立副本”
• A.block(...) = ...; 往往表示“我直接修改原矩阵中的局部区域”

4.1.2 常见块操作

Eigen::MatrixXd A(6, 6);
A.setRandom();

// 复制出一个3×3子块
Eigen::MatrixXd block_copy = A.block(1, 1, 3, 3);  // 从(1,1)开始取3×3块

// 直接修改原矩阵中的块
A.block(0, 0, 2, 2) = Eigen::Matrix2d::Identity();

// 行和列操作：注意行向量和列向量的类型不同
Eigen::RowVectorXd row0 = A.row(0);
Eigen::VectorXd col1 = A.col(1);

// 特定块操作
Eigen::MatrixXd top_left = A.topLeftCorner(3, 3);
Eigen::MatrixXd bottom_right = A.bottomRightCorner(3, 3);
Eigen::MatrixXd top_rows = A.topRows(3);
Eigen::MatrixXd left_cols = A.leftCols(3);

// 向量特定操作
Eigen::VectorXd v(10);
Eigen::VectorXd head_copy = v.head(3);         // 前3个元素（复制）
Eigen::VectorXd tail_copy = v.tail(3);         // 后3个元素（复制）
Eigen::VectorXd segment_copy = v.segment(2, 4);// 从索引2开始的4个元素（复制）

// slicing / indexing API 在 Eigen 3.4 已引入；Eigen 5.0 中需要注意 last/all 的命名空间变化
Eigen::VectorXd even = v(Eigen::seq(0, Eigen::placeholders::last, 2));  // 每2个元素取一个

注意事项：

• A.row(i) 返回的是行表达式，语义上更接近行向量
• A.col(j) 返回的是列表达式，语义上更接近列向量
• 若只是“读一块数据”，可以直接使用表达式；若要长期保存或脱离原矩阵使用，再显式复制

4.2 Array 类：逐元素运算

(Matrix 与 Array 的核心差别见 基础篇 和 矩阵基础篇。)

Eigen::ArrayXXd A(3, 3);
A << 1, 2, 3, 4, 5, 6, 7, 8, 9;

// 逐元素运算
Eigen::ArrayXXd B = A * 2;      // 逐元素乘法
Eigen::ArrayXXd C = A.sqrt();   // 平方根
Eigen::ArrayXXd D = A.exp();    // 指数
Eigen::ArrayXXd E = A.log();    // 对数
Eigen::ArrayXXd F = A.pow(2);   // 平方
Eigen::ArrayXXd G = A.abs();    // 绝对值

// 比较运算
Eigen::Array<bool, 3, 3> mask = A > 5;

// Matrix与Array互转
Eigen::MatrixXd M = A.matrix();     // Array -> Matrix
Eigen::ArrayXXd A2 = M.array();     // Matrix -> Array

// 链式操作
Eigen::MatrixXd result = (M.array() * 2).exp().matrix();

4.3 Map 类：零拷贝内存映射

Map 将现有内存缓冲区直接映射为 Eigen 对象，无需数据拷贝。适合以下场景：

• 已有外部内存（C 数组、std::vector、设备/共享内存缓冲区）
• 将其“看成”一个 Eigen 矩阵或向量
• 希望避免额外复制

可以把 Map 理解为：
“不是自己的矩阵数据，只是临时包装成 Eigen 类型来用。”

// 从C数组映射
double raw_data[9] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
Eigen::Map<Eigen::Matrix3d> mat(raw_data);  // 列优先映射

// 从std::vector映射
std::vector<double> vec = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
Eigen::Map<Eigen::MatrixXd> mat_from_vec(vec.data(), 3, 4);

// 修改原始数据
mat(0, 0) = 100;  // 直接修改raw_data[0]

// 行优先映射
Eigen::Map<Eigen::Matrix<double, 3, 3, Eigen::RowMajor>> mat_row(raw_data);

// 步长映射
double stride_data[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
Eigen::Map<Eigen::VectorXd, 0, Eigen::InnerStride<2>> strided_vec(stride_data, 5);

注意事项：

• Map 对象不拥有数据，必须确保原始内存在 Map 生命周期内有效
• Map<Eigen::Matrix3d> 默认按 列主序 解释内存；也就是说，上面的 raw_data 会按列填入矩阵，而不是按很多 C/C++ 初学者直觉中的“逐行填入”
• 需要按行解释时显式使用 RowMajor
• 如果数据未正确对齐，可能影响 SIMD 优化，某些场景下还可能触发对齐相关问题（详见 第九章：调试与排错篇 的 9.2 错误3：内存对齐错误）
• 只读映射使用 Map<const MatrixXd>
• Map 更强调“把外部内存包装成 Eigen 对象”，而不是“通用函数参数适配”

4.4 Ref 类：函数参数传递

给函数传递矩阵/向量参数时有四种可选写法。先给出决策框架，再解释 Ref 的原理。

4.4.1 四种参数写法的对比

4.4.2 MatrixXd& 为什么不够用？

直接用 MatrixXd& 作为参数有一个容易被忽略的限制：它只能绑定到类型完全匹配的左值。以下调用全部编译失败：

void foo(const Eigen::MatrixXd& mat);

Eigen::MatrixXd A(4, 4);
Eigen::Matrix4d B;               // 固定大小矩阵——NOT MatrixXd
Eigen::Map<Eigen::MatrixXd> M;   // Map 类型——NOT MatrixXd

foo(A);                          // OK
foo(A.block(1, 1, 2, 2));       // 编译错误：block 返回的是 Block<...> 表达式，不是 MatrixXd
foo(B);                          // 编译错误：Matrix4d ≠ MatrixXd
foo(M);                          // 编译错误：Map<MatrixXd> ≠ MatrixXd

这就是为什么 Eigen 提供了 Ref——它通过内部的类型擦除机制，可以绑定到任何布局兼容的矩阵表达式上。

4.4.3 Ref 的工作原理与限制

Ref<MatrixXd> 在构造时做两件事：

1. 检查传入表达式的内存布局是否兼容（步长、连续性、对齐），不兼容则编译报错
2. 包装——如果兼容，Ref 内部持有指向原数据的指针，零拷贝

// Ref 可以接受多种来源
void process(Eigen::Ref<Eigen::MatrixXd> mat) {
    mat(0, 0) = 999;  // 直接修改原始数据（零拷贝视图）
}

Eigen::MatrixXd A(4, 4);
process(A);                              // OK：完整矩阵
process(A.block(1, 1, 3, 3));           // OK：块表达式（布局兼容时）
process(A.topLeftCorner(3, 3));         // OK：角块表达式

Eigen::Matrix4d B;
process(B);                              // OK：固定大小矩阵（布局兼容时）

Ref 的局限：它不是万能适配器。以下情况会出问题：

// 布局不兼容
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> A_row(4, 4);
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor> A_col(4, 4);
A_row.setRandom();
A_col.setRandom();

// 以下可能在编译期或运行时报错（RowMajor vs ColMajor 的内步长可能不兼容）
// process(A_row.block(0, 0, 3, 3));

4.4.4 关键区分：可写 Ref vs 只读 Ref<const>

一个重要细节：Ref<const MatrixXd> 比 Ref<MatrixXd> 更灵活。

原因是 const 版本不要求可写语义，对传入表达式的内存布局约束更宽松。大多数只需读数据的函数用 Ref<const MatrixXd>。

// 只读函数：推荐 Ref<const MatrixXd>
double compute_norm(Eigen::Ref<const Eigen::MatrixXd> mat) {
    return mat.norm();  // 不修改 mat，只读取
}

// 可写函数：用 Ref<MatrixXd>
void zero_first_row(Eigen::Ref<Eigen::MatrixXd> mat) {
    mat.row(0).setZero();  // 需要写入
}

4.4.5 决策树：到底用哪个？

需要修改传入矩阵？
├── 是 → 调用方可能传 block/Map/固定大小矩阵？
│   ├── 是 → 用 Ref<MatrixXd>
│   └── 否 → 用 MatrixXd&（最简单，零开销）
└── 否（只读）→ 调用方可能传 block/Map/固定大小矩阵？
    ├── 是 → 用 Ref<const MatrixXd>（最推荐）
    ├── 否 → 用 const MatrixXd&（最简单，零开销）
    └── 表达式类型极度不确定 → 用 template<typename Derived>

实际中的推荐默认值：

• 只读参数：优先 Ref<const MatrixXd>。开销极小，兼容性远超 const MatrixXd&
• 可写参数：优先 Ref<MatrixXd>，除非调用方永远只传 MatrixXd 本身
• 泛型代码：用模板参数 template<typename Derived> + const MatrixBase<Derived>&

4.4.6 完整示例

#include <Eigen/Dense>
#include <iostream>

// 只读参数：接受 block、Map、固定大小等各种来源
double frobenius_squared(Eigen::Ref<const Eigen::MatrixXd> mat) {
    return mat.squaredNorm();
}

// 可写参数：接受 block、Map、固定大小等各种来源
void add_identity(Eigen::Ref<Eigen::MatrixXd> mat) {
    mat.diagonal().array() += 1.0;
}

// 泛型版本：接受任意表达式（包括 A + B 这类中间表达式）
template<typename Derived>
double trace_abs(const Eigen::MatrixBase<Derived>& mat) {
    return mat.template cast<double>().cwiseAbs().sum();
}

int main() {
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(3, 3);

    // 直接传 MatrixXd —— 所有写法都 OK
    frobenius_squared(A);
    add_identity(A);

    // 传 block —— MatrixXd& 会失败，Ref 能接住
    frobenius_squared(A.block(0, 0, 2, 2));
    add_identity(A.block(0, 0, 2, 2));  // 直接修改 A 的子块

    // 传固定大小矩阵
    Eigen::Matrix3d B = Eigen::Matrix3d::Random();
    frobenius_squared(B);
    add_identity(B);

    // 传表达式 —— 只有模板能接住
    // frobenius_squared(A * 2.0);  // 编译错误：A * 2.0 是临时表达式，无法绑定到 Ref
    double t = trace_abs(A * 2.0);  // 模板 OK

    return 0;
}

注意：A * 2.0 这类表达式是纯代数表达式，不存储连续内存，因此无法绑定到 Ref（即使是 Ref<const>）。如果需要接收这种中间表达式，只能用模板参数。

Map vs Ref 选择指南：

可以记一个简单规则：

• 已有一块内存，要把它“看成矩阵“：优先想 Map
• 在写函数参数，希望兼容常见 Eigen 对象：优先想 Ref（只读用 Ref<const>，可写用 Ref）
• 想接受非常多种表达式类型（包括中间表达式）：优先考虑模板

4.5 Broadcasting（广播）

Broadcasting 是 Eigen 中一种将向量自动“扩展“以匹配矩阵维度的机制，用于逐元素运算。它类似于 NumPy 的 broadcasting 概念，但 Eigen 的机制更显式。

基本概念

对矩阵每行加不同值或每列乘不同系数时，broadcasting 提供简洁写法：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::MatrixXd mat(3, 4);
    mat << 1, 2, 3, 4,
           5, 6, 7, 8,
           9, 10, 11, 12;
    
    Eigen::VectorXd v(3);
    v << 1, 2, 3;
    
    // 方式1：使用 colwise/rowwise + 显式广播
    Eigen::MatrixXd result = mat.colwise() + v;  // v 被"广播"到每一列
    // 等价于：result.col(j) = mat.col(j) + v
    std::cout << "mat.colwise() + v:\n" << result << "\n\n";
    
    // 方式2：逐行运算
    Eigen::RowVectorXd r(4);
    r << 10, 20, 30, 40;
    result = mat.rowwise() + r;  // r 被"广播"到每一行
    std::cout << "mat.rowwise() + r:\n" << result << "\n\n";
    
    // 方式3：结合 Array 使用（更多元素级操作）
    result = (mat.array().colwise() * v.array()).matrix();
    std::cout << "mat.colwise() * v:\n" << result << "\n\n";
    
    return 0;
}

输出：

mat.colwise() + v:
 2  3  4  5
 7  8  9 10
12 13 14 15

mat.rowwise() + r:
11 22 33 44
15 26 37 48
19 30 41 52

mat.colwise() * v:
 1  2  3  4
10 12 14 16
27 30 33 36

常用 Broadcasting 操作：

注意事项：

• colwise() 搭配列向量使用，rowwise() 搭配行向量使用
• Broadcasting 操作返回的是表达式而非副本，可高效链式调用
• 维度必须兼容：行数匹配 colwise，列数匹配 rowwise

4.6 Reshape（重塑）

reshaped() 以不同的行/列布局重新解释同一个矩阵，不拷贝数据。

#include <Eigen/Dense>
#include <iostream>

int main() {
    // 创建一个 2x6 矩阵
    Eigen::MatrixXd A(2, 6);
    A << 1, 2, 3, 4, 5, 6,
         7, 8, 9, 10, 11, 12;
    
    std::cout << "原始矩阵 A (2x6):\n" << A << "\n\n";
    
    // 重塑为 4x3（仍然是同一块内存）
    Eigen::MatrixXd B = A.reshaped(4, 3);
    std::cout << "重塑后 B (4x3):\n" << B << "\n\n";
    
    // 重塑为列向量（固定形状）
    Eigen::VectorXd v = A.reshaped();
    std::cout << "重塑为向量:\n" << v.transpose() << "\n\n";
    
    // 原地修改重塑后的视图会影响原矩阵
    A.reshaped(3, 4).row(0).setConstant(99);
    std::cout << "修改重塑视图后的 A:\n" << A << "\n\n";
    
    // 固定大小重塑
    Eigen::Matrix<double, 3, 4> C = A.reshaped<Eigen::FixedRows>(3, 4);
    std::cout << "固定大小重塑:\n" << C << "\n";
    
    return 0;
}

关键点：

• reshaped() 返回的依然是视图，不拷贝数据
• 如要独立副本，将结果赋值给具体类型（如上面的 Eigen::MatrixXd B）
• 元素总数必须不变：rows * cols == original.rows() * original.cols()
• 数据按列优先顺序重新排列

reshaped() vs resize()：

4.7 STL 互操作

Eigen 提供了与 C++ 标准库的兼容接口，便于在 STL 算法和容器中使用 Eigen 对象。

迭代器

#include <Eigen/Dense>
#include <algorithm>
#include <iostream>

int main() {
    Eigen::VectorXd v(5);
    v << 5, 3, 1, 4, 2;
    
    // 使用 STL 算法操作 Eigen 向量
    std::sort(v.begin(), v.end());
    std::cout << "排序后: " << v.transpose() << "\n";  // 1 2 3 4 5
    
    // 查找元素
    auto it = std::find(v.begin(), v.end(), 3);
    if (it != v.end())
        std::cout << "找到 3，索引: " << (it - v.begin()) << "\n";
    
    // 逐行排序矩阵
    Eigen::MatrixXd M(3, 4);
    M.setRandom();
    for (int i = 0; i < M.rows(); ++i)
        std::sort(M.row(i).begin(), M.row(i).end());
    std::cout << "逐行排序后的 M:\n" << M << "\n";
    
    return 0;
}

Range-for 遍历

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::MatrixXd M(3, 4);
    M.setRandom();
    
    // 遍历所有元素（列优先顺序）
    std::cout << "所有元素（列优先）: ";
    for (auto val : M.reshaped()) {
        std::cout << val << " ";
    }
    std::cout << "\n";
    
    return 0;
}

STL 容器中的 Eigen 类型

#include <Eigen/Dense>
#include <vector>
#include <map>

// 动态大小类型可直接放入 STL 容器（无对齐问题）
std::vector<Eigen::MatrixXd> matrices;
std::map<int, Eigen::VectorXd> vectors;

// 固定大小可向量化类型需要 aligned_allocator
// 详见第 9 章：调试与排错篇
std::vector<Eigen::Vector4d, Eigen::aligned_allocator<Eigen::Vector4d>> fixed_vectors;

对应官方文档：Broadcasting | Reshape and Slicing | STL iterators and algorithms | Block operations | The Array class and coefficient-wise operations | Map class

五、线性代数篇：矩阵分解与求解

5.1 线性方程组求解

Eigen提供多种求解器应对不同类型的线性系统。

密集矩阵求解器选择指南

线性方程组求解示例

#include <Eigen/Dense>
#include <iostream>

int main() {
    // 系统 Ax = b
    Eigen::Matrix3d A;
    A << 2, -1, 0,
         -1, 2, -1,
         0, -1, 2;
    
    Eigen::Vector3d b(1, 2, 3);
    
    // ========== 方法1：Cholesky分解（已知A为对称正定时首选）==========
    Eigen::LLT<Eigen::Matrix3d> llt(A);
    Eigen::Vector3d x1 = llt.solve(b);
    
    // ========== 方法2：LU分解（一般可逆方阵的常见选择）==========
    Eigen::PartialPivLU<Eigen::Matrix3d> lu(A);
    Eigen::Vector3d x2 = lu.solve(b);
    
    // ========== 方法3：列主元QR（更稳妥的通用方案）==========
    Eigen::ColPivHouseholderQR<Eigen::Matrix3d> qr(A);
    Eigen::Vector3d x3 = qr.solve(b);
    
    // ========== 不推荐作为默认写法：显式求逆 ==========
    // 对于线性系统，通常优先使用 solve()，而不是 A.inverse() * b
    Eigen::Vector3d x_inverse = A.inverse() * b;
    
    // ========== 验证解的正确性 ==========
    double residual = (A * x3 - b).norm();
    std::cout << "QR 残差范数: " << residual << "\n";
    std::cout << "逆矩阵写法的残差范数: " << (A * x_inverse - b).norm() << "\n";
    
    // ========== 多次求解（分解复用）==========
    // 当需要求解多个右端项时，先分解再求解更高效
    Eigen::MatrixXd B(3, 5);
    B.setRandom();
    Eigen::MatrixXd X = lu.solve(B);  // 同时求解5个右端项
    
    return 0;
}

最小二乘问题

#include <Eigen/Dense>
#include <iostream>

int main() {
    // 超定系统：寻找x使得 ||Ax - b||最小
    Eigen::MatrixXd A(10, 3);
    A.setRandom();
    
    Eigen::VectorXd b(10);
    b.setRandom();
    
    // ========== 方法1：正规方程（最快，但数值稳定性较差）==========
    Eigen::VectorXd x1 = (A.transpose() * A).ldlt().solve(A.transpose() * b);
    
    // ========== 方法2：列主元QR（稳妥、常用）==========
    Eigen::VectorXd x2 = A.colPivHouseholderQr().solve(b);
    
    // ========== 方法3：CompleteOrthogonalDecomposition（推荐默认选择）==========
    // 对秩亏、欠定、最小范数解等情况更稳妥
    Eigen::VectorXd x3 = A.completeOrthogonalDecomposition().solve(b);
    
    // ========== 方法4：SVD（最稳健，但通常更慢）==========
    // 适合病态问题，或需要奇异值/奇异向量
    // 在 Eigen 5.x 中，thin/full U/V 的运行时选项已弃用，推荐使用编译时模板参数
    Eigen::VectorXd x4 = A.bdcSvd<Eigen::ComputeThinU | Eigen::ComputeThinV>().solve(b);
    
    // 计算残差
    std::cout << "ColPivHouseholderQR 残差: "
              << (A * x2 - b).norm() << "\n";
    std::cout << "CompleteOrthogonalDecomposition 残差: "
              << (A * x3 - b).norm() << "\n";
    std::cout << "SVD 残差: "
              << (A * x4 - b).norm() << "\n";
    
    return 0;
}

5.2 特征值分解

#include <Eigen/Dense>
#include <iostream>
#include <complex>

int main() {
    Eigen::Matrix3d A;
    A << 4, 2, 1,
         2, 5, 3,
         1, 3, 6;
    
    // ========== 自伴矩阵（对称/厄米特）的特征值分解 ==========
    // 使用SelfAdjointEigenSolver（更快更稳定）
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(A);
    
    if (eigensolver.info() != Eigen::Success) {
        std::cerr << "特征值分解失败!\n";
        return 1;
    }
    
    std::cout << "特征值:\n" << eigensolver.eigenvalues() << "\n\n";
    std::cout << "特征向量（列向量）:\n" << eigensolver.eigenvectors() << "\n";
    
    // 验证：A * v = lambda * v
    Eigen::Vector3d v = eigensolver.eigenvectors().col(0);
    double lambda = eigensolver.eigenvalues()(0);
    std::cout << "\n验证: ||Av - λv|| = " << (A * v - lambda * v).norm() << "\n";
    
    // ========== 一般矩阵的特征值分解 ==========
    Eigen::Matrix3d B;
    B << 1, 2, 3,
         4, 5, 6,
         7, 8, 10;
    
    Eigen::EigenSolver<Eigen::Matrix3d> es(B);
    if (es.info() == Eigen::Success) {
        // 特征值可能是复数
        std::cout << "\n特征值（可能复数）:\n" << es.eigenvalues() << "\n";
        
        // 提取实部（如果矩阵有复数特征值）
        Eigen::VectorXd real_eigenvalues = es.eigenvalues().real();
    }
    
    // ========== 广义特征值问题：Ax = λBx ==========
    Eigen::Matrix3d M;
    M.setRandom();
    M = M * M.transpose();  // 确保正定
    
    Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::Matrix3d> ges(A, M);
    std::cout << "\n广义特征值:\n" << ges.eigenvalues() << "\n";
    
    return 0;
}

5.3 奇异值分解(SVD)

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::MatrixXd A(5, 4);
    A.setRandom();
    
    // ========== BDCSVD（分治SVD，大型矩阵更快）==========
    // 在 Eigen 5.x 中，thin/full U/V 的运行时选项已弃用，推荐使用编译时模板参数
    Eigen::BDCSVD<Eigen::MatrixXd, Eigen::ComputeThinU | Eigen::ComputeThinV> svd1(A);

    // ========== 检查分解是否成功 ==========
    if (svd1.info() != Eigen::Success) {
        std::cerr << "SVD分解失败！可能是数值问题。\n";
        return -1;
    }
    
    // ========== JacobiSVD（Jacobi方法，更精确）==========
    // 在 Eigen 5.x 中，thin/full U/V 的运行时选项已弃用，推荐使用编译时模板参数
    Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::ComputeThinU | Eigen::ComputeThinV> svd2(A);
    
    // 同样可以检查分解状态
    if (svd2.info() == Eigen::Success) {
        std::cout << "JacobiSVD分解成功\n";
    }
    
    // 获取结果
    Eigen::VectorXd singular_values = svd1.singularValues();
    Eigen::MatrixXd U = svd1.matrixU();
    Eigen::MatrixXd V = svd1.matrixV();
    
    std::cout << "奇异值:\n" << singular_values << "\n\n";
    
    // ========== 矩阵伪逆（Moore-Penrose逆）==========
    // 方式1：推荐使用 SVD 的 solve() 方法获得伪逆解
    Eigen::MatrixXd I = Eigen::MatrixXd::Identity(A.rows(), A.rows());
    Eigen::MatrixXd A_pinv = svd1.solve(I);  // 利用 SVD 构造伪逆
    
    // 方式2：手动构造（仅作学习参考）
    // double tolerance = 1e-10 * singular_values(0);
    // Eigen::VectorXd singular_values_inv = singular_values;
    // for (int i = 0; i < singular_values.size(); ++i) {
    //     if (singular_values(i) > tolerance)
    //         singular_values_inv(i) = 1.0 / singular_values(i);
    //     else
    //         singular_values_inv(i) = 0;
    // }
    // Eigen::MatrixXd A_pinv = V * singular_values_inv.asDiagonal() * U.transpose();
    
    // ========== 矩阵秩 ==========
    int rank = svd1.rank();
    std::cout << "矩阵秩: " << rank << "\n";
    
    // ========== 矩阵条件数 ==========
    double cond = singular_values(0) / singular_values(singular_values.size() - 1);
    std::cout << "条件数: " << cond << "\n";
    
    // ========== 低秩近似 ==========
    // 保留前k个奇异值
    int k = 2;
    Eigen::MatrixXd A_approx = U.leftCols(k) * 
                               singular_values.head(k).asDiagonal() * 
                               V.leftCols(k).transpose();
    
    double approx_error = (A - A_approx).norm();
    std::cout << "低秩近似误差: " << approx_error << "\n";
    
    return 0;
}

info() 方法说明（SVD 分解仅返回 Success 或 NoConvergence，为完整参考，此处列出所有分解可能返回的状态）：

最佳实践：在生产代码中始终检查 info() 返回值，特别是在处理用户输入或不确定的矩阵时。

5.4 矩阵分解详解与选择

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Matrix3d A;
    A << 4, 2, 1,
         2, 5, 3,
         1, 3, 6;
    
    // ========== Cholesky分解：A = LL^T（A必须正定）==========
    Eigen::LLT<Eigen::Matrix3d> llt(A);
    Eigen::Matrix3d L = llt.matrixL();
    std::cout << "Cholesky L:\n" << L << "\n";
    std::cout << "验证 L*L^T:\n" << L * L.transpose() << "\n\n";

    // ========== 原地分解（Inplace Decomposition）：适合已分配好内存的场景 ==========
    Eigen::Matrix3d M = A;
    Eigen::LLT<Eigen::Ref<Eigen::Matrix3d>> llt_inplace(M);  // 在 M 的内存上原地分解
    // 注意：M 会被修改，原矩阵内容不再可用
    
    // ========== LDL^T分解（适合对称正定，也常用于某些半定/边界场景）==========
    Eigen::LDLT<Eigen::Matrix3d> ldlt(A);
    
    // ========== LU分解：PA = LU ==========
    Eigen::PartialPivLU<Eigen::Matrix3d> plu(A);
    Eigen::Matrix3d L_lu = Eigen::Matrix3d::Identity();
    L_lu.triangularView<Eigen::StrictlyLower>() = plu.matrixLU();
    Eigen::Matrix3d U_lu = plu.matrixLU().triangularView<Eigen::Upper>();
    std::cout << "LU分解:\nL =\n" << L_lu << "\nU =\n" << U_lu << "\n\n";
    
    // ========== QR分解：A = QR ==========
    Eigen::HouseholderQR<Eigen::Matrix3d> qr(A);
    Eigen::Matrix3d Q = qr.householderQ();
    Eigen::Matrix3d R = qr.matrixQR().triangularView<Eigen::Upper>();
    std::cout << "QR分解:\nQ =\n" << Q << "\nR =\n" << R << "\n";
    
    return 0;
}

对应官方文档：Linear algebra and decompositions | Catalogue of dense decompositions

六、几何变换篇：旋转与变换

本篇介绍Eigen的几何变换功能，包括旋转矩阵、四元数、欧拉角和仿射变换。

阅读提示：

• AngleAxis、Quaternion、Translation 这类类型更偏向“抽象变换表示”
• Transform、Affine3d、Isometry3d 更偏向“变换矩阵表示”
• 对单个旋转进行表达时，四元数和轴角通常更自然；对多个点批量变换时，旋转矩阵和变换矩阵往往更直接

6.1 旋转表示方式对比

三维空间中的旋转有多种数学表示方式，各有优劣：

6.2 四元数详解

什么是四元数？

四元数是一种扩展复数的数学概念，由爱尔兰数学家Hamilton于1843年发明。一个四元数表示为：

q = w + xi + yj + zk

其中 w 是实部（标量部分），(x, y, z) 是虚部（向量部分），i, j, k 满足：

• i² = j² = k² = ijk = -1

为什么四元数适合表示旋转？

1. 紧凑性：仅需4个数（vs 旋转矩阵9个数）
2. 无奇异性：不存在万向节锁问题
3. 插值友好：球面线性插值(SLERP)可平滑过渡
4. 数值稳定：归一化即可恢复正交性

四元数与旋转的关系：

绕单位向量 n 旋转角度 θ 的四元数为：

q = cos(θ/2) + sin(θ/2)(nₓi + nᵧj + n_zk)

#include <Eigen/Geometry>
#include <iostream>
#include <cmath>

int main() {
    // ========== 四元数创建方式 ==========
    
    // 方式1：直接赋值
    // 注意：Eigen中四元数构造函数参数顺序为 (w, x, y, z)
    // 但 coeffs() 返回的顺序是 (x, y, z, w)
    Eigen::Quaterniond q1(1, 0, 0, 0);  // 单位四元数 (w=1, x=0, y=0, z=0)，表示无旋转
    std::cout << "单位四元数 coeffs(): " << q1.coeffs().transpose() << "\n";
    // 输出: 单位四元数 coeffs(): 0 0 0 1  (顺序为 x, y, z, w)
    std::cout << "四元数元素: w=" << q1.w() << ", x=" << q1.x() 
              << ", y=" << q1.y() << ", z=" << q1.z() << "\n";
    // 输出: 四元数元素: w=1, x=0, y=0, z=0
    
    // 方式2：从轴角构造（绕Z轴旋转45度）
    // 注意：轴向量必须是单位向量；这里 (0,0,1) 已经是单位向量
    Eigen::AngleAxisd rotation_vector(EIGEN_PI / 4, Eigen::Vector3d(0, 0, 1));
    Eigen::Quaterniond q2(rotation_vector);
    std::cout << "绕Z轴45度: " << q2.coeffs().transpose() << "\n";
    // 输出: 绕Z轴45度: 0 0 0.382683 0.92388
    
    // 方式3：从旋转矩阵构造
    Eigen::Matrix3d R = rotation_vector.toRotationMatrix();
    Eigen::Quaterniond q3(R);
    
    // 方式4：从两个向量构造（从v1旋转到v2的最短旋转）
    Eigen::Vector3d v1(1, 0, 0), v2(0, 1, 0);
    Eigen::Quaterniond q4 = Eigen::Quaterniond::FromTwoVectors(v1, v2);
    std::cout << "X轴到Y轴的旋转: " << q4.coeffs().transpose() << "\n";
    // 输出: X轴到Y轴的旋转: 0 0 0.707107 0.707107
    
    // ========== 四元数基本操作 ==========
    
    // 归一化（确保表示有效旋转）
    q2.normalize();
    
    // 获取旋转矩阵
    Eigen::Matrix3d R_from_q = q2.toRotationMatrix();
    std::cout << "\n旋转矩阵:\n" << R_from_q << "\n";
    // 输出:
    // 旋转矩阵:
    //  0.707107 -0.707107         0
    //  0.707107  0.707107         0
    //         0         0         1
    
    // 旋转向量
    Eigen::Vector3d v(1, 0, 0);
    Eigen::Vector3d v_rotated = q2 * v;  // 等价于 q2.toRotationMatrix() * v
    std::cout << "旋转后的向量: " << v_rotated.transpose() << "\n";
    // 输出: 旋转后的向量: 0.707107 0.707107 0
    
    // ========== 四元数复合（旋转串联）==========
    
    // 先绕X轴90度，再绕Z轴90度
    Eigen::Quaterniond qx(Eigen::AngleAxisd(EIGEN_PI / 2, Eigen::Vector3d::UnitX()));
    Eigen::Quaterniond qz(Eigen::AngleAxisd(EIGEN_PI / 2, Eigen::Vector3d::UnitZ()));
    
    // 四元数乘法：q_combined = qz * qx 表示先qx后qz
    Eigen::Quaterniond q_combined = qz * qx;
    
    // 验证：旋转X轴
    Eigen::Vector3d x_axis(1, 0, 0);
    Eigen::Vector3d result = q_combined * x_axis;
    std::cout << "复合旋转后的X轴: " << result.transpose() << "\n";
    // 输出: 复合旋转后的X轴: 0 1 0
    
    // ========== 四元数插值（SLERP）==========
    
    // 球面线性插值：在两个旋转之间平滑过渡
    Eigen::Quaterniond q_start = Eigen::Quaterniond::Identity();
    Eigen::Quaterniond q_end(Eigen::AngleAxisd(EIGEN_PI / 2, Eigen::Vector3d::UnitZ()));
    
    // t=0时为q_start，t=1时为q_end，t=0.5为中间旋转
    for (double t : {0.0, 0.25, 0.5, 0.75, 1.0}) {
        Eigen::Quaterniond q_interp = q_start.slerp(t, q_end);
        Eigen::AngleAxisd aa(q_interp);
        std::cout << "t=" << t << ": 旋转角度 " << aa.angle() * 180 / EIGEN_PI << "度\n";
    }
    // 输出:
    // t=0: 旋转角度 0度
    // t=0.25: 旋转角度 22.5度
    // t=0.5: 旋转角度 45度
    // t=0.75: 旋转角度 67.5度
    // t=1: 旋转角度 90度
    
    return 0;
}

四元数使用注意事项：

1. 必须归一化：数值计算误差可能导致四元数不再单位化，定期调用normalize()
2. 乘法顺序：q1 * q2 表示先应用 q2，再应用 q1（与矩阵乘法一致）
3. 双倍覆盖：q 和 -q 表示相同的旋转，比较时需注意

6.3 欧拉角

什么是欧拉角？

欧拉角使用三个角度描述三维旋转，由数学家欧拉于1775年提出。常见的约定包括：

万向节锁问题：

当中间轴旋转±90°时，第一轴和第三轴重合，丢失一个自由度。

#include <Eigen/Geometry>
#include <iostream>
#include <cmath>

int main() {
    // ========== 欧拉角创建 ==========
    
    // 定义欧拉角（弧度）
    double yaw = EIGEN_PI / 6;    // 30度，绕Z轴
    double pitch = EIGEN_PI / 4;  // 45度，绕Y轴
    double roll = EIGEN_PI / 3;   // 60度，绕X轴
    
    // 从欧拉角创建旋转矩阵（ZYX顺序，内旋）
    // 内旋：依次绕物体自身的Z、Y、X轴旋转
    Eigen::Matrix3d R;
    R = Eigen::AngleAxisd(yaw, Eigen::Vector3d::UnitZ()) *
        Eigen::AngleAxisd(pitch, Eigen::Vector3d::UnitY()) *
        Eigen::AngleAxisd(roll, Eigen::Vector3d::UnitX());
    
    std::cout << "旋转矩阵（ZYX顺序）:\n" << R << "\n\n";
    
    // ========== 从旋转矩阵提取欧拉角 ==========
    
    // eulerAngles(2, 1, 0) 表示按 Z-Y-X 顺序提取
    // 注意：返回值会遵循 Eigen 当前实现的规范化约定；
    // 不同版本之间，等价旋转可能对应不同但等价的角度表示
    Eigen::Vector3d euler = R.eulerAngles(2, 1, 0);
    std::cout << "提取的欧拉角（弧度）: " << euler.transpose() << "\n";
    std::cout << "提取的欧拉角（度）: " 
              << (euler * 180 / EIGEN_PI).transpose() << "\n\n";
    // 输出:
    // 提取的欧拉角（弧度）: 0.523599 0.785398 1.0472
    // 提取的欧拉角（度）: 30 45 60
    
    // ========== 万向节锁演示 ==========
    
    std::cout << "===== 万向节锁演示 =====\n";
    
    // 当pitch = ±90°时，发生万向节锁
    double pitch_locked = EIGEN_PI / 2;
    Eigen::Matrix3d R_locked;
    R_locked = Eigen::AngleAxisd(yaw, Eigen::Vector3d::UnitZ()) *
               Eigen::AngleAxisd(pitch_locked, Eigen::Vector3d::UnitY()) *
               Eigen::AngleAxisd(roll, Eigen::Vector3d::UnitX());
    
    // 尝试提取欧拉角
    Eigen::Vector3d euler_locked = R_locked.eulerAngles(2, 1, 0);
    std::cout << "pitch=90°时提取的欧拉角: " << (euler_locked * 180 / EIGEN_PI).transpose() << "\n";
    // 这里的结果不应被当作唯一“标准答案”；
    // 在奇异位置附近，不同但等价的欧拉角表示都可能出现
    
    // ========== 避免万向节锁的建议 ==========
    
    std::cout << "\n避免万向节锁的方法:\n";
    std::cout << "1. 使用四元数代替欧拉角进行旋转计算\n";
    std::cout << "2. 限制pitch角度范围（如-89°到89°）\n";
    std::cout << "3. 在需要用户输入时才使用欧拉角，内部计算用四元数\n";
    
    return 0;
}

Eigen 5.0兼容性说明：Eigen 5.0 中欧拉角的返回值形式更加规范（canonical），因此与旧版本相比，返回值的角度区间或具体数值表示可能发生变化；但 eulerAngles(a0, a1, a2) 的轴顺序语义并没有改变。如需跨版本兼容，不应把某个固定角度区间当作协议，应把结果视为“当前版本下的一种等价表示”。

6.4 仿射变换

什么是变换矩阵？

变换矩阵是4×4的齐次变换矩阵，综合表示旋转和平移：

T = | R₃ₓ₃  t₃ₓ₁ |
    | 0₁ₓ₃   1    |

其中 R 是3×3旋转矩阵，t 是3×1平移向量

为什么使用齐次坐标？

• 统一表示旋转和平移
• 变换复合简化为矩阵乘法
• 便于处理刚体运动链

在实际使用中，还要区分三类对象，它们在同一个变换矩阵下的变换方式不同：

• 点（point）：T * point，既受旋转影响，也受平移影响
• 方向向量（direction vector）：T.rotation() * vector，只受旋转影响，平移对其无意义
• 法向量（normal vector）：在一般仿射变换（含非均匀缩放）下，法向量不能简单用 T 或 T.rotation() 变换。正确变换为 (T^{-1})^T，即逆转置。对于 Affine3d/Isometry3d，实践中常用 T.rotation().inverse().transpose() * normal

为什么法向量需要特殊处理？ 考虑一个平面经过非均匀缩放后，其法向量方向不会简单与顶点同步缩放。对于纯刚体变换（旋转+平移，Isometry3d），法向量与方向向量的变换相同，因为旋转矩阵是正交矩阵（R^{-1} = R^T），逆转置等于自身。

// 区分三种向量的变换方式
Eigen::Affine3d T = Eigen::Translation3d(1, 2, 3)
                  * Eigen::AngleAxisd(EIGEN_PI / 4, Eigen::Vector3d::UnitZ())
                  * Eigen::Scaling(2.0, 1.0, 1.0);  // X 方向缩放 2 倍

Eigen::Vector3d point(1, 0, 0);
Eigen::Vector3d dir(0, 1, 0);
Eigen::Vector3d normal(1, 0, 0);  // 假设为某平面的法向量

// 点：直接乘变换矩阵
Eigen::Vector3d transformed_point = T * point;

// 方向向量：只用旋转部分（去掉平移）
Eigen::Vector3d transformed_dir = T.rotation() * dir;

// 法向量：用逆转置（对于含非均匀缩放的仿射变换）
Eigen::Vector3d transformed_normal =
    T.rotation().inverse().transpose() * normal;
// 注意：若为纯刚体变换（Isometry3d），rotation().inverse().transpose() == rotation()，
// 此时法向量与方向向量的变换相同

#include <Eigen/Geometry>
#include <iostream>
#include <cmath>

int main() {
    // ========== 创建变换矩阵 ==========
    
    // 方式1：分别设置旋转和平移
    Eigen::Affine3d T1 = Eigen::Affine3d::Identity();
    T1.rotate(Eigen::AngleAxisd(EIGEN_PI / 4, Eigen::Vector3d::UnitZ()));
    T1.pretranslate(Eigen::Vector3d(1, 2, 3));
    
    std::cout << "变换矩阵T1:\n" << T1.matrix() << "\n\n";
    // 输出:
    // 变换矩阵T1:
    //  0.707107 -0.707107         0         1
    //  0.707107  0.707107         0         2
    //         0         0         1         3
    //         0         0         0         1
    
    // 方式2：链式构造（推荐）
    // 注意：乘法顺序从右到左执行
    Eigen::Affine3d T2 = 
        Eigen::Translation3d(1, 2, 3) *                    // 后平移
        Eigen::AngleAxisd(EIGEN_PI / 4, Eigen::Vector3d::UnitZ());  // 先旋转
    
    // ========== 应用变换 ==========
    
    Eigen::Vector3d point(1, 0, 0);
    Eigen::Vector3d transformed = T2 * point;
    
    std::cout << "原始点: " << point.transpose() << "\n";
    std::cout << "变换后: " << transformed.transpose() << "\n\n";
    // 输出:
    // 原始点: 1 0 0
    // 变换后: 1.70711 2.70711 3
    
    // ========== 变换复合 ==========
    
    Eigen::Affine3d T_a = Eigen::Translation3d(1, 0, 0) * Eigen::Affine3d::Identity();
    Eigen::Affine3d T_b = Eigen::Translation3d(0, 1, 0) * Eigen::Affine3d::Identity();
    
    // T_combined = T_b * T_a 表示先T_a后T_b
    Eigen::Affine3d T_combined = T_b * T_a;
    
    Eigen::Vector3d origin(0, 0, 0);
    std::cout << "原点经T_a后: " << (T_a * origin).transpose() << "\n";
    std::cout << "再经T_b后: " << (T_combined * origin).transpose() << "\n\n";
    // 输出:
    // 原点经T_a后: 1 0 0
    // 再经T_b后: 1 1 0
    
    // ========== 逆变换 ==========
    
    Eigen::Affine3d T_inv = T2.inverse();
    Eigen::Vector3d back = T_inv * transformed;
    
    std::cout << "逆变换还原: " << back.transpose() << "\n";
    std::cout << "还原误差: " << (back - point).norm() << "\n\n";
    // 输出:
    // 逆变换还原: 1 0 0
    // 还原误差: 0
    
    // ========== 提取旋转和平移 ==========
    
    Eigen::Matrix3d R = T2.rotation();
    Eigen::Vector3d t = T2.translation();
    
    std::cout << "旋转部分:\n" << R << "\n";
    std::cout << "平移部分: " << t.transpose() << "\n";
    
    return 0;
}

6.5 缩放变换

#include <Eigen/Geometry>
#include <iostream>

int main() {
    // ========== 各向同性缩放 ==========
    Eigen::UniformScaling<double> uniform_scale(2.0);  // 各方向放大2倍
    
    // ========== 各向异性缩放 ==========
    Eigen::DiagonalMatrix<double, 3> aniso_scale(2, 3, 4);  // X放大2倍，Y放大3倍，Z放大4倍
    
    // ========== 完整变换链 ==========
    // 变换顺序：缩放 → 旋转 → 平移（从右到左）
    Eigen::Quaterniond q(Eigen::AngleAxisd(EIGEN_PI / 4, Eigen::Vector3d::UnitZ()));
    
    Eigen::Affine3d full_transform = 
        Eigen::Translation3d(1, 0, 0) *    // 最后平移
        q *                                 // 然后旋转
        Eigen::Scaling(2.0, 3.0, 4.0);      // 先缩放
    
    Eigen::Vector3d p(1, 1, 1);
    Eigen::Vector3d result = full_transform * p;
    
    std::cout << "完整变换结果: " << result.transpose() << "\n";
    // 输出: 完整变换结果: -0.828427 4.24264 4
    
    return 0;
}

6.6 实战：三维刚体变换

#include <Eigen/Geometry>
#include <iostream>
#include <vector>
#include <cmath>

// 3自由度机械臂的正向运动学
class RobotArm3DOF {
private:
    std::vector<double> link_lengths_;  // 各连杆长度
    
public:
    RobotArm3DOF() : link_lengths_({1.0, 0.8, 0.5}) {}
    
    // 计算末端执行器位置
    Eigen::Vector3d forwardKinematics(const std::vector<double>& joint_angles) {
        if (joint_angles.size() != 3) {
            throw std::invalid_argument("需要3个关节角度");
        }
        
        // 从基坐标系开始
        Eigen::Affine3d T = Eigen::Affine3d::Identity();
        
        for (size_t i = 0; i < 3; ++i) {
            // 旋转关节（绕Z轴）
            T.rotate(Eigen::AngleAxisd(joint_angles[i], Eigen::Vector3d::UnitZ()));
            // 沿X轴平移到下一个关节
            T.translate(Eigen::Vector3d(link_lengths_[i], 0, 0));
        }
        
        return T.translation();
    }
    
    // 计算完整的末端位姿
    Eigen::Affine3d forwardKinematicsFull(const std::vector<double>& joint_angles) {
        Eigen::Affine3d T = Eigen::Affine3d::Identity();
        
        for (size_t i = 0; i < 3; ++i) {
            T.rotate(Eigen::AngleAxisd(joint_angles[i], Eigen::Vector3d::UnitZ()));
            T.translate(Eigen::Vector3d(link_lengths_[i], 0, 0));
        }
        
        return T;
    }
    
    // 打印工作空间边界（简化版）
    void printWorkspaceBoundary() {
        std::cout << "工作空间半径范围: [" 
                  << link_lengths_[2] << ", " 
                  << (link_lengths_[0] + link_lengths_[1] + link_lengths_[2]) 
                  << "]\n";
    }
};

int main() {
    RobotArm3DOF arm;
    arm.printWorkspaceBoundary();
    
    // 设置关节角度（弧度）
    std::vector<double> angles = {EIGEN_PI / 6, EIGEN_PI / 4, EIGEN_PI / 3};
    
    std::cout << "\n关节角度（度）:\n";
    for (size_t i = 0; i < angles.size(); ++i) {
        std::cout << "  关节" << i+1 << ": " << angles[i] * 180 / EIGEN_PI << "°\n";
    }
    
    // 计算末端位置
    Eigen::Vector3d end_pos = arm.forwardKinematics(angles);
    std::cout << "\n末端执行器位置: " << end_pos.transpose() << "\n";
    
    // 计算末端姿态
    Eigen::Affine3d end_pose = arm.forwardKinematicsFull(angles);
    Eigen::Vector3d end_orientation = end_pose.rotation() * Eigen::Vector3d::UnitX();
    std::cout << "末端执行器朝向: " << end_orientation.transpose() << "\n";
    
    // 示例输出:
    // 工作空间半径范围: [0.5, 2.3]
    // 
    // 关节角度（度）:
    //   关节1: 30°
    //   关节2: 45°
    //   关节3: 60°
    // 
    // 末端执行器位置: 1.0458 1.43934 0
    // 末端执行器朝向: -0.258819 0.965926 0
    
    return 0;
}

6.7 几何变换常见问题

Q: 四元数和旋转矩阵如何选择？

A:

• 使用四元数：姿态估计、动画插值、存储/传输
• 使用旋转矩阵：矩阵运算密集、需要直接访问旋转后的坐标轴

Q: 变换矩阵乘法顺序如何理解？

A: 变换从右向左执行。T = T1 * T2 表示先应用 T2，再应用 T1。

Q: 如何处理坐标系转换？

A: 使用变换矩阵的逆。如果 T_A_B 表示从B到A的变换，则 T_B_A = T_A_B.inverse()。

对应官方文档：Geometry | Space transformations

七、稀疏矩阵篇：大规模计算

7.1 稀疏矩阵基础

什么是稀疏矩阵？

稀疏矩阵是指大部分元素为零的矩阵。在实际应用中，许多大规模矩阵（如有限元刚度矩阵、社交网络邻接矩阵）的非零元素占比通常小于1%。

存储方式对比：

选择指南：

下面这些阈值只适合作为粗略经验，不是通用定律。是否该用稀疏矩阵，除了非零元素占比，还取决于：

• 主要操作是矩阵-向量乘、矩阵分解，还是频繁随机访问
• 稀疏结构是否规则（如带状、块状、图结构）
• 是否有适合的稀疏求解器与预处理器
• 实际平台上的基准测试结果

非零元素占比明显较高  → 往往更适合先考虑稠密矩阵
非零元素占比明显较低  → 往往更适合先考虑稀疏矩阵
中间情况              → 以实际运算模式和基准测试为准

7.2 稀疏矩阵存储格式

Eigen支持两种主要的稀疏矩阵存储格式：

1. 列压缩存储（CSC，默认）

以下为简化示意图（Eigen 实际在非压缩模式下使用 4 个紧凑数组：Values、InnerIndices、OuterStarts、InnerNNZs）。

原始矩阵:          CSC存储:
| 1 0 2 |         values:  [1, 4, 2, 5, 3]
| 0 0 3 |   →     indices: [0, 1, 0, 1, 1]  (行索引)
| 4 5 0 |         starts:  [0, 2, 4, 5]     (每列起始位置)

2. 行压缩存储（CSR）

原始矩阵:          CSR存储:
| 1 0 2 |         values:  [1, 2, 3, 4, 5]
| 0 0 3 |   →     indices: [0, 2, 2, 0, 1]  (列索引)
| 4 5 0 |         starts:  [0, 2, 3, 5]     (每行起始位置)

#include <Eigen/Sparse>
#include <iostream>

int main() {
    // ========== 创建稀疏矩阵 ==========
    
    // 默认使用列压缩存储（CSC）
    Eigen::SparseMatrix<double> sp_col_major(1000, 1000);
    
    // 行压缩存储（CSR）
    Eigen::SparseMatrix<double, Eigen::RowMajor> sp_row_major(1000, 1000);
    
    // ========== 高效填充方式：三元组列表 ==========
    
    // 预估非零元素数量（优化内存分配）
    // 这里的 VectorXi::Constant(1000, 3) 表示：按列预估，每一列大约有 3 个非零元素
    sp_col_major.reserve(Eigen::VectorXi::Constant(1000, 3));
    
    std::vector<Eigen::Triplet<double>> triplets;
    triplets.reserve(3000);  // 预分配内存
    
    // 填充三对角矩阵
    for (int i = 0; i < 1000; ++i) {
        triplets.emplace_back(i, i, 2.0);       // 主对角线
        if (i > 0)
            triplets.emplace_back(i, i-1, -1.0); // 下对角线
        if (i < 999)
            triplets.emplace_back(i, i+1, -1.0); // 上对角线
    }
    
    // 批量设置（高效）
    sp_col_major.setFromTriplets(triplets.begin(), triplets.end());
    
    // 压缩存储（移除多余空间）
    sp_col_major.makeCompressed();
    
    std::cout << "矩阵维度: " << sp_col_major.rows() << " x " << sp_col_major.cols() << "\n";
    std::cout << "非零元素数: " << sp_col_major.nonZeros() << "\n";
    std::cout << "稀疏密度: " 
              << 100.0 * sp_col_major.nonZeros() / (sp_col_major.rows() * sp_col_major.cols()) 
              << "%\n";
    
    // 输出:
    // 矩阵维度: 1000 x 1000
    // 非零元素数: 2998
    // 稀疏密度: 0.2998%
    
    // ========== 内存占用对比 ==========
    
    // 稠密矩阵内存: 1000 * 1000 * 8 bytes = 8 MB
    // 稀疏矩阵内存: 2998 * (8 + 4 + 间接开销) ≈ 48 KB
    std::cout << "\n内存节省: " 
              << (1.0 - (double)sp_col_major.nonZeros() / (1000 * 1000)) * 100 
              << "%\n";
    // 输出: 内存节省: 99.7002%
    
    return 0;
}

7.2.1 稀疏矩阵运算

#include <Eigen/Sparse>
#include <iostream>

int main() {
    // 创建两个稀疏矩阵
    Eigen::SparseMatrix<double> A(4, 4);
    std::vector<Eigen::Triplet<double>> triplets;
    
    triplets.emplace_back(0, 0, 1);
    triplets.emplace_back(1, 1, 2);
    triplets.emplace_back(2, 2, 3);
    triplets.emplace_back(3, 3, 4);
    triplets.emplace_back(0, 1, 0.5);
    triplets.emplace_back(1, 0, 0.5);
    
    A.setFromTriplets(triplets.begin(), triplets.end());
    
    // ========== 基本运算 ==========
    
    // 稀疏矩阵 + 稀疏矩阵
    Eigen::SparseMatrix<double> B = A + A;
    
    // 稀疏矩阵 * 标量
    Eigen::SparseMatrix<double> C = 2.0 * A;
    
    // 稀疏矩阵 * 稠密向量
    Eigen::VectorXd v(4);
    v << 1, 2, 3, 4;
    Eigen::VectorXd result = A * v;
    
    std::cout << "A * v = " << result.transpose() << "\n";
    // 输出: A * v = 2 4.5 9 16
    
    // 稀疏矩阵 * 稀疏矩阵
    Eigen::SparseMatrix<double> D = A * A;
    
    // ========== 访问元素 ==========
    
    // 方式1：迭代非零元素
    std::cout << "\n非零元素:\n";
    for (int k = 0; k < A.outerSize(); ++k) {
        for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it) {
            std::cout << "  A(" << it.row() << "," << it.col() << ") = " << it.value() << "\n";
        }
    }
    
    // 方式2：检查特定元素（通常比稠密矩阵随机访问更慢）
    // 实际代价与矩阵是否压缩、存储顺序以及访问模式有关
    double val = A.coeff(0, 0);  // 返回A(0,0)，不存在则返回0
    
    // 方式3：修改元素（需要非压缩状态）
    A.coeffRef(2, 3) = 1.5;  // 插入新元素
    
    return 0;
}

7.2.2 稀疏直接求解器

直接求解器对中小规模问题（通常 < 10000 阶）效率高且结果可靠。以下示例展示主要直接求解器的用法：

#include <Eigen/Sparse>
#include <iostream>

int main() {
    // 创建大规模稀疏矩阵（二维拉普拉斯算子）
    int n = 50;  // 50x50网格
    Eigen::SparseMatrix<double> A(n * n, n * n);
    std::vector<Eigen::Triplet<double>> triplets;
    triplets.reserve(5 * n * n);
    
    for (int i = 0; i < n * n; ++i) {
        triplets.emplace_back(i, i, 4.0);  // 主对角线
        if (i % n > 0) triplets.emplace_back(i, i - 1, -1.0);     // 左邻居
        if (i % n < n - 1) triplets.emplace_back(i, i + 1, -1.0); // 右邻居
        if (i >= n) triplets.emplace_back(i, i - n, -1.0);        // 上邻居
        if (i < n * (n - 1)) triplets.emplace_back(i, i + n, -1.0); // 下邻居
    }
    A.setFromTriplets(triplets.begin(), triplets.end());
    A.makeCompressed();
    
    Eigen::VectorXd b = Eigen::VectorXd::Ones(n * n);
    Eigen::VectorXd x;
    
    std::cout << "矩阵大小: " << n * n << " x " << n * n << "\n";
    std::cout << "非零元素: " << A.nonZeros() << "\n\n";
    
    // SimplicialLLT：对称正定矩阵，Cholesky分解
    Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> llt;
    llt.compute(A);
    if (llt.info() == Eigen::Success) {
        x = llt.solve(b);
        std::cout << "SimplicialLLT 残差: " << (A * x - b).norm() / b.norm() << "\n";
    }
    
    // SimplicialLDLT：更稳定，可处理半正定
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt;
    ldlt.compute(A);
    if (ldlt.info() == Eigen::Success) {
        x = ldlt.solve(b);
        std::cout << "SimplicialLDLT 残差: " << (A * x - b).norm() / b.norm() << "\n";
    }
    
    // SparseLU：通用求解器，非对称方阵
    Eigen::SparseLU<Eigen::SparseMatrix<double>> sparse_lu;
    sparse_lu.compute(A);
    if (sparse_lu.info() == Eigen::Success) {
        x = sparse_lu.solve(b);
        std::cout << "SparseLU 残差: " << (A * x - b).norm() / b.norm() << "\n";
    }
    
    // SparseQR：适用于矩形矩阵和最小二乘问题
    Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> sparse_qr;
    sparse_qr.compute(A);
    if (sparse_qr.info() == Eigen::Success) {
        x = sparse_qr.solve(b);
        std::cout << "SparseQR 残差: " << (A * x - b).norm() / b.norm() << "\n";
    }
    
    return 0;
}

迭代求解器和预处理技术见 7.3 节，适合大规模问题。

7.2.3 稀疏求解器选择指南

决策流程：

矩阵是否对称正定？
  ├─ 是 → 规模 < 10000？
  │        ├─ 是 → SimplicialLLT
  │        └─ 否 → ConjugateGradient + 预处理
  └─ 否 → 需要最小二乘解？
           ├─ 是 → SparseQR
           └─ 否 → 规模 < 10000？
                    ├─ 是 → SparseLU
                    └─ 否 → BiCGSTAB + 预处理

7.3 迭代求解器

对于大规模稀疏线性方程组，迭代求解器通常比直接求解器更高效。

#include <Eigen/Sparse>
#include <Eigen/IterativeLinearSolvers>
#include <iostream>

int main() {
    // 创建大规模稀疏矩阵（示例：二维拉普拉斯算子）
    int n = 100;
    Eigen::SparseMatrix<double> A(n * n, n * n);
    std::vector<Eigen::Triplet<double>> triplets;
    
    for (int i = 0; i < n * n; ++i) {
        triplets.push_back(Eigen::Triplet<double>(i, i, 4.0));
        if (i % n > 0) triplets.push_back(Eigen::Triplet<double>(i, i - 1, -1.0));
        if (i % n < n - 1) triplets.push_back(Eigen::Triplet<double>(i, i + 1, -1.0));
        if (i >= n) triplets.push_back(Eigen::Triplet<double>(i, i - n, -1.0));
        if (i < n * (n - 1)) triplets.push_back(Eigen::Triplet<double>(i, i + n, -1.0));
    }
    A.setFromTriplets(triplets.begin(), triplets.end());
    
    Eigen::VectorXd b = Eigen::VectorXd::Ones(n * n);
    Eigen::VectorXd x;
    
    // ========== 共轭梯度法（CG）==========
    // 适用于对称正定矩阵，收敛最快
    Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> cg;
    cg.compute(A);
    cg.setTolerance(1e-10);
    cg.setMaxIterations(1000);
    x = cg.solve(b);
    
    std::cout << "共轭梯度法:\n";
    std::cout << "  迭代次数: " << cg.iterations() << "\n";
    std::cout << "  估计误差: " << cg.error() << "\n";
    std::cout << "  状态: " << (cg.info() == Eigen::Success ? "成功" : "失败") << "\n\n";
    
    // ========== 最小二乘共轭梯度法（LeastSquaresConjugateGradient）==========
    // 适用于稀疏最小二乘问题
    Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<double>> lscg;
    lscg.compute(A);
    x = lscg.solve(b);
    
    std::cout << "LeastSquaresConjugateGradient:\n";
    std::cout << "  迭代次数: " << lscg.iterations() << "\n\n";
    
    // ========== 双共轭梯度稳定法（BiCGSTAB）==========
    // 适用于非对称矩阵
    Eigen::BiCGSTAB<Eigen::SparseMatrix<double>> bicg;
    bicg.compute(A);
    bicg.setTolerance(1e-10);
    x = bicg.solve(b);
    
    std::cout << "BiCGSTAB:\n";
    std::cout << "  迭代次数: " << bicg.iterations() << "\n";
    std::cout << "  估计误差: " << bicg.error() << "\n\n";
    
    // ========== 带预处理的迭代求解器 ==========
    // 预处理可以显著加速收敛
    using Preconditioner = Eigen::IncompleteCholesky<double>;
    Eigen::ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower | Eigen::Upper, Preconditioner> cg_precond;
    cg_precond.compute(A);
    x = cg_precond.solve(b);
    
    std::cout << "预处理CG:\n";
    std::cout << "  迭代次数: " << cg_precond.iterations() << " (预处理后)\n";
    
    return 0;
}

迭代求解器选择指南

预处理技术

// 常用预处理类型
using DiagonalPrecond = Eigen::DiagonalPreconditioner<double>;           // 对角预处理
using IncompleteLUT = Eigen::IncompleteLUT<double>;                      // 不完全LU分解
using IncompleteCholesky = Eigen::IncompleteCholesky<double>;            // 不完全Cholesky分解

// 使用预处理
Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, IncompleteLUT> solver;
solver.preconditioner().setDroptol(0.001);  // 设置丢弃容差
solver.compute(A);
x = solver.solve(b);

常见问题 (FAQ)

Q: 求解器返回info() != Success怎么办？

A: 检查以下几点：

1. 矩阵是否奇异（行列式为零）
2. 对于Cholesky，矩阵是否正定
3. 数值是否溢出/下溢
4. 尝试更稳定的求解器（如SVD）

Q: 如何选择合适的求解器？

A: 参考上方 7.2.3 节的决策流程。

Q: 特征值计算结果与MATLAB不同？

A: 特征向量的符号和顺序可能不同，这是正常的。验证：

// 检查 A*v = lambda*v
assert((A * v - lambda * v).norm() < 1e-10);

练习题

1. 线性系统: 实现一个函数，使用LU分解求解多个右端项的线性系统。
2. PCA实现: 给定数据矩阵X，实现主成分分析（PCA），返回主成分和方差解释率。
3. 矩阵平方根: 利用特征值分解计算正定矩阵的平方根。

对应官方文档：Sparse matrix manipulations

八、优化篇：性能调优技巧

8.1 编译优化选项

# 推荐的通用高性能编译配置
g++ -std=c++17 -O3 -march=native -DEIGEN_NO_DEBUG -DNDEBUG \
    -I/path/to/eigen -o optimized_app app.cpp

# 如需启用多线程并行，添加 -fopenmp（详见 8.4 节）
# g++ -std=c++17 -O3 -march=native -fopenmp -DEIGEN_NO_DEBUG -DNDEBUG \
#     -I/path/to/eigen -o optimized_app app.cpp

# 注意：-ffast-math 可能导致数值精度问题
# 仅在对精度要求不高的场景使用，例如图形渲染
# g++ ... -ffast-math ...  # 不推荐用于科学计算

关于-ffast-math的警告：

-ffast-math选项会启用激进的浮点优化，可能导致：

• 违反IEEE 754浮点标准
• NaN和Inf检查失效
• 数值精度降低
• 某些数学函数结果不准确

建议：

• 科学计算和数值分析：不要使用-ffast-math
• 图形渲染、游戏开发：可以谨慎使用
• 如需高性能，优先考虑-O3 -march=native -ffp-contract=fast
• -fopenmp：仅在矩阵规模大且多核场景下启用（详见 8.4 节）

向量化指令集

8.2 内存管理最佳实践

固定大小 vs 动态大小

#include <Eigen/Dense>
#include <chrono>

void benchmark() {
    const int N = 1000000;
    
    // 固定大小（栈分配）- 在这类小矩阵场景下通常更快
    auto start = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < N; ++i) {
        Eigen::Matrix4d A = Eigen::Matrix4d::Random();
        Eigen::Matrix4d B = Eigen::Matrix4d::Random();
        Eigen::Matrix4d C = A * B;
    }
    auto fixed_time = std::chrono::high_resolution_clock::now() - start;
    
    // 动态大小（堆分配）- 在这类小矩阵场景下通常较慢
    start = std::chrono::high_resolution_clock::now();
    for (int i = 0; i < N; ++i) {
        Eigen::MatrixXd A = Eigen::MatrixXd::Random(4, 4);
        Eigen::MatrixXd B = Eigen::MatrixXd::Random(4, 4);
        Eigen::MatrixXd C = A * B;
    }
    auto dynamic_time = std::chrono::high_resolution_clock::now() - start;
    
    // 对这类小尺寸矩阵，固定大小常常会明显更快；
    // 具体收益依赖尺寸、编译器、目标平台和表达式形态
}

避免动态内存分配

// 不推荐：运行时确定大小导致堆分配
template<typename T>
T compute(int n) {
    Eigen::VectorXd v(n);  // 堆分配
    // ...
}

// 推荐：使用模板参数在编译时确定大小
template<int N>
Eigen::Matrix<double, N, 1> compute() {
    Eigen::Matrix<double, N, 1> v;  // 栈分配
    // ...
}

// 或者使用“动态大小 + 编译期最大上限”
// 这种类型仍然是动态尺寸，只是给出一个编译期上限，便于优化小尺寸场景
template<int MaxN>
using VectorMax = Eigen::Matrix<double, Eigen::Dynamic, 1, 0, MaxN, 1>;

void compute(int n) {
    assert(n <= 100 && "n 超过编译期最大上限");
    VectorMax<100> v(n);  // 运行时长度为 n，最大不超过 100
}

预分配内存

// 动态矩阵预分配
Eigen::MatrixXd A;
A.resize(1000, 1000);  // 一次性分配

// 多次操作时避免重复分配
Eigen::MatrixXd B(1000, 1000), C(1000, 1000);
B.setRandom();
C.setRandom();

for (int i = 0; i < 1000; ++i) {
    // A.resize(...)  // 不要这样做！
    A.noalias() = B * C;  // 当 A 与 B/C 不重叠时，可用 noalias() 避免不必要的临时对象
}

8.3 表达式模板与临时对象

Eigen的表达式模板是其高性能的核心机制，通过延迟求值避免创建临时对象。

什么是表达式模板？

Eigen使用表达式模板（Expression Templates）技术来优化运算。传统方式会为每个中间结果创建临时对象，而表达式模板通过延迟求值，将整个表达式编译成一个高效的循环。

传统方式（无表达式模板）：

// 假设没有表达式模板
MatrixXd D = A + B + C;
// 实际执行：
// 1. temp1 = A + B;  // 创建临时矩阵
// 2. D = temp1 + C;  // 再创建临时矩阵
// 结果：多次内存分配、多次遍历数据

Eigen的表达式模板：

MatrixXd D = A + B + C;
// 实际执行：
// 编译器生成类似以下的代码：
for (int i = 0; i < D.size(); ++i) {
    D(i) = A(i) + B(i) + C(i);
}
// 结果：零临时对象、单次遍历、编译时优化

表达式模板的工作原理

// Eigen内部实现（简化）
template<typename Lhs, typename Rhs>
struct MatrixSum {
    const Lhs& lhs;
    const Rhs& rhs;
    
    double operator()(int i, int j) const {
        return lhs(i, j) + rhs(i, j);
    }
};

// A + B 返回的是 MatrixSum<MatrixXd, MatrixXd>，而不是新的矩阵
auto sum = A + B;  // sum是一个表达式对象，不计算结果

MatrixXd C = sum;  // 此时才真正计算（赋值时触发求值）

何时会创建临时对象？

情况1：显式转换类型

MatrixXd A, B, C;
VectorXd v;

// 会创建临时对象
MatrixXd temp = A * B;        // 显式创建临时对象
v = (A * B).col(0);           // A*B的结果需要先计算出来

// 避免不必要的临时对象
C.noalias() = A * B;          // 当C与A/B不重叠时可使用noalias()

情况2：表达式过于复杂

// 复杂表达式不一定更慢，但当存在子表达式复用、别名风险或可读性问题时，
// 显式拆分通常更容易分析和调优
MatrixXd result = A * B + C * D * E + F.transpose() * G;

// 如果需要复用中间结果，可以显式求值
MatrixXd temp1 = (C * D).eval();
MatrixXd temp2 = temp1 * E;
MatrixXd result2 = A * B + temp2 + F.transpose() * G;

情况3：使用auto导致的陷阱

MatrixXd A, B;

// 错误：auto保存的是表达式对象，不是结果
auto sum = A + B;  // sum是表达式对象，不拥有数据

// 如果A或B在sum求值前被修改，结果可能与"保存当时结果"的直觉不一致；
// 如果它们已经销毁，则还可能出现悬空引用问题
A(0, 0) = 999;     // 修改A
MatrixXd C = sum;  // 这里得到的是"此刻再求值"的结果，而不是定义sum时就保存好的结果

// 更危险的场景：操作数是临时对象
auto bad = (A + B) + MatrixXd::Ones(3, 3);
// (A + B) 和 MatrixXd::Ones(3, 3) 作为临时对象在表达式结束后被销毁，
// 此时 bad 内部引用的是已销毁的对象（悬空引用），后续求值导致未定义行为

// 正确做法
MatrixXd sum = A + B;  // 显式指定类型，立即求值

noalias()优化

MatrixXd A, B, C;

// 普通赋值（可能创建临时对象）
C = A * B;  // 编译器可能创建临时对象，因为C可能与A或B重叠

// 使用noalias()（承诺C不与A或B重叠）
C.noalias() = A * B;  // 零临时对象，性能最优

// 注意：如果C确实与A或B重叠，使用noalias()会导致错误结果
C = A * C;            // 正确：编译器会处理重叠
C.noalias() = A * C;  // 错误：结果不正确！

性能优化技巧

技巧1：链式操作

// 推荐：链式操作，单次遍历
MatrixXd D = (A + B).cwiseProduct(C);

// 不推荐：分步操作，多次遍历
MatrixXd temp = A + B;
MatrixXd D = temp.cwiseProduct(C);

技巧2：使用in-place操作

MatrixXd A, B;

// 推荐：原地操作
A += B;
A *= 2.0;

// 不推荐：创建新对象
A = A + B;
A = A * 2.0;

技巧3：预分配结果矩阵

MatrixXd A(1000, 1000), B(1000, 1000);
A.setRandom();
B.setRandom();

// 推荐：预分配
MatrixXd C(1000, 1000);
C.noalias() = A * B;

// 不推荐：动态分配
MatrixXd C = A * B;  // 可能触发重新分配

常见陷阱与解决方案

8.4 多线程并行计算

关于 OpenMP

OpenMP（Open Multi-Processing）是一种面向共享内存多处理器系统的并行编程模型。它通过编译器指令（#pragma omp）在 C/C++/Fortran 代码中标记可并行区域，由编译器和运行时库自动将工作分配到多个线程。

Eigen 内部利用 OpenMP 对部分运算进行自动并行化——无需修改算法代码，只需在编译时开启 OpenMP 支持即可。

何时启用 OpenMP

满足以下条件时，启用 OpenMP 能带来显著的性能提升：

• 矩阵规模较大（矩阵乘法约 128×128 以上，LU 分解需更大规模）
• 机器有多个物理核心（通常 ≥ 4 核）
• 程序的计算瓶颈确实在矩阵运算上（而非 I/O 或其他逻辑）

以下情况不应启用：

• 矩阵规模很小（< 128×128）：线程调度开销超过并行收益
• 单核或双核机器
• 程序已在应用层使用 OpenMP 手动并行：与 Eigen 内部线程竞争导致性能下降
• 程序需要在嵌套并行场景运行（Eigen 不支持嵌套 OpenMP）

线程模型

Eigen 使用 OpenMP 的默认线程池。线程数由以下优先级决定（高到低）：

1. Eigen::setNbThreads(n) — Eigen API 显式设置
2. omp_set_num_threads(n) — OpenMP API 设置
3. OMP_NUM_THREADS 环境变量
4. 系统默认（通常等于逻辑核心数）

查询当前线程数：Eigen::nbThreads()

重要警告：OpenMP 报告的“核心数“通常是逻辑核心数（含超线程），而 Eigen 的矩阵乘法内核已几乎占满单个物理核心的算力。在超线程开启的机器上，线程数应设为物理核心数而非逻辑核心数，否则性能会因缓存污染和调度开销显著下降（可能慢数倍）。

构建环境要求

macOS 注意：Apple Clang 默认不带 OpenMP 运行时，需通过 Homebrew 安装 libomp 并添加 -Xpreprocessor -fopenmp -lomp。

CMake 配置：

find_package(OpenMP)
if(OpenMP_CXX_FOUND)
    target_link_libraries(myapp PRIVATE OpenMP::OpenMP_CXX)
endif()

并行化支持范围

线程安全注意事项

• Matrix::Random() / setRandom() 不是线程安全的（基于 std::rand）。多线程中生成随机矩阵需使用 C++11 <random> 引擎自行填充
• 自定义标量类型在 OpenMP 环境下不应抛出异常，否则导致未定义行为
• 应用层若已使用 OpenMP，应在 Eigen 侧禁用并行：Eigen::setNbThreads(1) 或定义 EIGEN_DONT_PARALLELIZE

与 -fopenmp 编译选项的关系

-fopenmp 是编译器和链接器选项，Eigen 仅在检测到 _OPENMP 宏已定义（即编译器开启了 OpenMP）时才会启用内部并行代码路径。不加 -fopenmp 编译的程序中，Eigen 的所有运算均为单线程执行。

禁用 Eigen 内部多线程

// 运行时禁用
Eigen::setNbThreads(1);

// 或编译时全局禁用
#define EIGEN_DONT_PARALLELIZE

并行化示例

#include <Eigen/Dense>
#include <iostream>

int main() {
    // 查看 Eigen 使用的线程数
    std::cout << "Eigen 线程数: " << Eigen::nbThreads() << "\n";

    // 显式设置线程数（设为物理核心数）
    Eigen::setNbThreads(4);

    Eigen::MatrixXd A = Eigen::MatrixXd::Random(2000, 2000);
    Eigen::MatrixXd B = Eigen::MatrixXd::Random(2000, 2000);

    Eigen::MatrixXd C = A * B;  // 自动并行化

    return 0;
}

# 编译（GCC/Linux）
g++ -std=c++17 -O3 -march=native -fopenmp \
    -I/path/to/eigen -o parallel_demo parallel_demo.cpp

# 编译（macOS + Homebrew libomp）
g++ -std=c++17 -O3 -march=native -Xpreprocessor -fopenmp -lomp \
    -I/path/to/eigen -o parallel_demo parallel_demo.cpp

# 运行时控制线程数
OMP_NUM_THREADS=4 ./parallel_demo

8.5 向量化优化

#include <Eigen/Dense>
#include <iostream>

int main() {
    // 检查向量化支持
    #ifdef __AVX2__
    std::cout << "编译时启用了AVX2\n";
    #endif
    
    #ifdef __SSE4_2__
    std::cout << "编译时启用了SSE4.2\n";
    #endif
    
    // 确保数据对齐以获得最佳向量化性能
    // Eigen的固定大小矩阵自动对齐
    Eigen::Matrix4d A;  // 自动16/32字节对齐
    
    // 动态矩阵也自动对齐
    Eigen::MatrixXd B(100, 100);
    
    // 使用Map时确保对齐
    alignas(32) double data[100];  // C++11对齐
    Eigen::Map<Eigen::VectorXd> v(data, 100);
    
    return 0;
}

8.6 性能分析技巧

#include <Eigen/Dense>
#include <chrono>

class Timer {
    std::chrono::high_resolution_clock::time_point start;
    const char* name;
public:
    Timer(const char* n) : name(n) {
        start = std::chrono::high_resolution_clock::now();
    }
    ~Timer() {
        auto end = std::chrono::high_resolution_clock::now();
        auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
        std::cout << name << ": " << duration.count() << " μs\n";
    }
};

void benchmark_operations() {
    const int N = 1000;
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(N, N);
    Eigen::MatrixXd B = Eigen::MatrixXd::Random(N, N);
    Eigen::VectorXd v = Eigen::VectorXd::Random(N);
    
    {
        Timer t("矩阵乘法");
        Eigen::MatrixXd C = A * B;
    }
    
    {
        Timer t("矩阵-向量乘法");
        Eigen::VectorXd w = A * v;
    }
    
    {
        Timer t("LU分解");
        Eigen::PartialPivLU<Eigen::MatrixXd> lu(A);
    }
    
    {
        Timer t("Cholesky分解");
        Eigen::LLT<Eigen::MatrixXd> llt(A.transpose() * A);
    }
}

常见问题 (FAQ)

Q: 如何检查Eigen是否使用了向量化？

A: 定义EIGEN_VECTORIZE宏查看：

#include <Eigen/Core>
#ifdef EIGEN_VECTORIZE
std::cout << "向量化已启用\n";
#endif

Q: 程序崩溃提示“unaligned memory access“

A: 这通常是结构体中包含固定大小可向量化 Eigen 类型导致的：

// 错误
struct MyData {
    Eigen::Vector4d v;  // 需要 16/32 字节对齐
    int x;
};

// 正确
struct MyData {
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW  // 重载 new 以正确对齐
    Eigen::Vector4d v;
    int x;
};

详细讨论：对齐问题涉及四种根源（类成员、STL 容器、按值传递、GCC 栈对齐）和多种解决方案。完整的分类讲解、C++17 与旧标准的差异、验证方法，见 调试与排错篇 9.2 节。

Q: 如何禁用多线程？

A:

Eigen::setNbThreads(1);
// 或编译时定义
// -DEIGEN_DONT_PARALLELIZE

8.7 CUDA支持

Eigen 早已支持在 CUDA 内核中使用固定大小的矩阵和向量；这并不是 Eigen 5.0 才新增的能力。

使用示例：

#include <Eigen/Core>

__global__ void matrix_multiply_kernel(float* result, const float* a, const float* b, int n) {
    // 在CUDA内核中使用Eigen固定大小矩阵
    Eigen::Map<const Eigen::Matrix3f> mat_a(a);
    Eigen::Map<const Eigen::Matrix3f> mat_b(b);
    Eigen::Map<Eigen::Matrix3f> mat_result(result);
    
    mat_result = mat_a * mat_b;
}

// 注意：
// 1. 适合固定大小矩阵（如 Matrix3f、Vector4d 等）
// 2. 动态大小对象通常不适合直接在 kernel 中这样使用
// 3. 需要使用 Eigen::Map 映射设备内存
// 4. 在 .cu 文件中通常需要关闭主机端 SIMD 向量化，并注意索引类型与编译器兼容性

编译配置：

find_package(CUDA REQUIRED)
include_directories(${EIGEN_INCLUDE_DIR})

cuda_add_executable(my_cuda_app kernel.cu)

限制：

• 仅支持固定大小类型
• 不支持动态内存分配
• 部分高级功能不可用

练习题

1. 性能测试: 比较固定大小和动态大小矩阵在不同尺寸下的性能差异。
2. 内存分析: 使用noalias()优化矩阵链式乘法，测量性能提升。
3. 并行优化: 测试不同线程数对大型矩阵乘法的影响，找到最优线程数。

对应官方文档：Using Eigen in multi-threaded applications

九、调试与排错篇

9.1 常见编译错误

错误1：头文件未找到

fatal error: Eigen/Core: No such file or directory

解决方案：

# 检查Eigen安装路径
g++ -I/usr/include/eigen3 ...  # 系统安装路径
# 或
g++ -I/path/to/eigen ...       # 手动安装路径

错误2：C++标准不满足

error: Eigen requires at least c++14 support

解决方案：

g++ -std=c++14 ...  # 或 -std=c++17

错误3：模块头文件缺失或包含不完整

error: ‘JacobiSVD’ is not a member of ‘Eigen’

或：

error: incomplete type ‘Eigen::SelfAdjointEigenSolver<...>’ used in nested name specifier

常见原因：

• 只包含了 <Eigen/Core>，但实际使用了分解、特征值或几何相关功能
• 头文件包含层级不完整

解决方案：

# 若使用大多数稠密矩阵功能，直接包含：
#include <Eigen/Dense>

# 若希望按模块最小化包含，可按需使用：
#include <Eigen/SVD>
#include <Eigen/Eigenvalues>
#include <Eigen/Geometry>

说明： 这类错误比线程链接问题更常见，也更符合 Eigen 初学者的实际踩坑路径。

9.2 运行时错误

错误1：矩阵尺寸不匹配

Eigen::MatrixXd A(3, 3), B(4, 4);
Eigen::MatrixXd C = A * B;  // 运行时错误（调试模式）

调试方法：

// 使用断言检查
assert(A.cols() == B.rows() && "矩阵维度不匹配");

// 或使用Eigen的调试模式（默认开启）
// 编译时添加 -DEIGEN_NO_DEBUG 禁用（发布模式）

错误2：奇异矩阵

Eigen::Matrix2d A;
A << 1, 2,
     2, 4;  // 奇异矩阵（行列式为0）
Eigen::Matrix2d A_inv = A.inverse();  // 结果可能为Inf/NaN

调试方法：

// 检查行列式
double det = A.determinant();
if (std::abs(det) < 1e-10) {
    std::cerr << "警告：矩阵接近奇异，行列式 = " << det << "\n";
    // 使用伪逆或其他方法
}

// 或使用条件数（这里直接使用编译时模板参数写法）
Eigen::JacobiSVD<Eigen::Matrix2d, Eigen::ComputeFullU | Eigen::ComputeFullV> svd(A);
const auto singular_values = svd.singularValues();
double sigma_max = singular_values(0);
double sigma_min = singular_values(singular_values.size() - 1);

if (sigma_min <= 1e-15) {
    std::cerr << "警告：矩阵奇异或接近奇异，最小奇异值 = " << sigma_min << "\n";
} else {
    double cond = sigma_max / sigma_min;
    if (cond > 1e10) {
        std::cerr << "警告：矩阵病态，条件数 = " << cond << "\n";
    }
}

错误3：内存对齐错误

前置概念：什么是“固定大小可向量化“（fixed-size vectorizable）类型？

Eigen 中 “固定大小可向量化”（fixed-size vectorizable） 类型的判断标准：

该类型在编译期大小固定，且大小是 16 字节的整数倍。

因为 SIMD 指令（SSE/AVX/AVX512）以 128 bit（16 字节）、256 bit（32 字节）或 512 bit（64 字节）的数据包为单位工作，这些数据包在对齐到数据包大小时读写效率最高。Eigen 会为满足条件的类型申请对应的对齐（16/32/64 字节），从而充分利用 SIMD 向量化。

典型的“固定大小可向量化“类型包括：

注意：动态大小类型（如 VectorXd、MatrixXd）在堆上分配自己的系数数组，会自行处理绝对对齐，因此不受下文讨论的问题影响。

问题根源：

Eigen 使用 SIMD 指令（SSE/AVX）来优化性能。固定大小可向量化类型必须在正确对齐的内存地址上创建，否则 SIMD 指令会触发段错误（segfault）。

Eigen 通常会自动处理这些对齐问题——它通过 alignas 关键字声明对齐要求，并重载 operator new 来返回对齐的指针。但在少数边界情况下，这些对齐设置会被绕过，从而导致崩溃。

实际错误信息：

程序因对齐问题崩溃时，通常出现以下断言失败信息：

my_program: path/to/eigen/Eigen/src/Core/DenseStorage.h:44:
Eigen::internal::matrix_array<T, Size, MatrixOptions, Align>::internal::matrix_array()
[with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:
Assertion (reinterpret_cast<size_t>(array) & (sizemask)) == 0 && "this assertion
is explained here: ... **** READ THIS WEB PAGE !!! ****"' failed.

定位问题根源：

用调试器获取 backtrace 定位触发断言的代码位置：

$ gdb ./my_program          # 启动 GDB
> run                       # 运行程序
...                         # 复现崩溃
> bt                        # 获取 backtrace

四大根源与对应解决方案

对齐问题的根源可分为四类，下面逐一说明。

C++17 用户注意：如果你使用 C++17 且编译器足够新（GCC >= 7, Clang >= 5, MSVC >= 19.12），编译器会自动处理大部分对齐问题。在这种情况下，EIGEN_MAKE_ALIGNED_OPERATOR_NEW 宏为空操作，STL 容器的对齐也由编译器保证。但如果你需要兼容旧编译器或旧标准，请继续阅读下文。

根源 1：结构体/类中包含固定大小可向量化 Eigen 成员

类中包含固定大小可向量化 Eigen 类型成员，通过 new 动态分配时，返回的指针可能不满足对齐要求。

// 问题代码
class Foo {
    Eigen::Vector4d v;  // 需要 32 字节对齐（AVX 下）
    int x;
};

Foo *foo = new Foo;  // foo 指针可能未对齐，导致 foo->v 也未对齐

原因：Eigen::Vector4d 本身有对齐的 operator new，但 Foo 没有。当 new Foo 时，用的是 Foo 的 operator new（默认无对齐），返回的 foo 指针可能不是 32 字节对齐的。成员 v 的对齐是相对于 Foo 起始地址的——如果 foo 不对齐，foo->v 也不会对齐。

解决方案 A：使用 EIGEN_MAKE_ALIGNED_OPERATOR_NEW（推荐）

class Foo {
    Eigen::Vector4d v;
    int x;
public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW  // 重载 operator new，保证返回对齐指针
};

Foo *foo = new Foo;  // 安全

解决方案 B：条件性对齐（模板参数决定是否需对齐）

template<int N> class Foo {
    typedef Eigen::Matrix<float, N, 1> Vector;
    enum { NeedsToAlign = (sizeof(Vector) % 16) == 0 };
    Vector v;
public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign)
};

Foo<4> *foo4 = new Foo<4>;  // 128bit 对齐
Foo<3> *foo3 = new Foo<3>;  // 默认系统对齐（Vector3f 为 12 字节，不需要特殊对齐）

解决方案 C：禁用对齐（性能有损失，但不影响正确性）

class Foo {
    Eigen::Matrix<double, 4, 1, Eigen::DontAlign> v;
};

关于成员排列顺序：Eigen 成员不需要放在类的最前面——对齐是自动处理的。但从节省内存的角度，建议将对齐要求高的成员放在前面，避免编译器插入过多的填充字节。例如 AVX 下：

class Foo {
    Eigen::Vector4d v;  // 32 字节对齐
    int x;              // 后面可能有填充
};
// 优于
class Foo {
    int x;              // 后面编译器会插入 24 字节填充
    Eigen::Vector4d v;
};

根源 2：STL 容器或手动内存分配

将固定大小可向量化 Eigen 类型（或包含它们的类）放入 STL 容器时，默认的 std::allocator 不保证对齐。

// 问题代码
std::vector<Eigen::Vector2d> my_vector;
std::map<int, Eigen::Matrix4f> my_map;

解决方案 A：使用 Eigen::aligned_allocator

// vector
std::vector<Eigen::Vector2d, Eigen::aligned_allocator<Eigen::Vector2d>> my_vector;

// map（注意：第三个参数 std::less 是默认值，但必须显式写出才能指定第四个参数）
std::map<int, Eigen::Vector4d, std::less<int>,
         Eigen::aligned_allocator<std::pair<const int, Eigen::Vector4d>>> my_map;

解决方案 B：特化 std::vector（C++98/03 环境）

#include <Eigen/StdVector>

// 针对特定类型特化 vector
EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(Eigen::Matrix2d)

// 之后可以直接使用
std::vector<Eigen::Matrix2d> vec;  // 自动使用对齐分配

注意：特化必须在所有用到 std::vector<该类型> 的代码之前声明，否则编译器会用默认分配器编译那些实例。

手动内存分配：使用 std::make_shared 或 std::allocate_shared 时，同样需要传入对齐分配器。

根源 3：按值传递 Eigen 对象

在 C++17 之前，将固定大小可向量化 Eigen 对象按值传递给函数不仅低效，还可能导致崩溃：

// 问题代码（pre-C++17）
void my_function(Eigen::Vector2d v);  // 按值传递，可能违反对齐

struct Foo { Eigen::Vector2d v; };
void my_function(Foo v);              // 同样的问题

解决方案：改用 const 引用传递

void my_function(const Eigen::Vector2d& v);
void my_function(const Foo& v);

注意：函数返回 Eigen 对象按值是没问题的，只有参数传递时存在此问题。C++17 及以上标准中，编译器已能正确处理此场景。

根源 4：编译器对栈对齐的错误假设（GCC on Windows）

这是仅见于 Windows 上 GCC（MinGW、TDM-GCC）的历史问题（GCC 4.5 后已修复）。

现象：在一个普通函数中声明局部 Eigen 变量，也会触发对齐断言：

void foo() {
    Eigen::Quaternionf q;  // 可能崩溃
}

原因：GCC 默认假设栈已 16 字节对齐，但在 Windows 上栈仅保证 4 字节对齐。当函数从其他线程或被其他编译器编译的二进制调用时，栈对齐可能被破坏。

解决方案（任选其一）：

// 方案 A：单函数标记
__attribute__((force_align_arg_pointer)) void foo() {
    Eigen::Quaternionf q;
}

# 方案 B：全局编译选项（告知 GCC 栈只保证 4 字节对齐）
g++ -mincoming-stack-boundary=2 ...

# 方案 C：全局编译选项（等同于给所有函数加 force_align_arg_pointer）
g++ -mstackrealign ...

如果不需要最优向量化，如何彻底禁用对齐？

三种方式（按影响范围递增）：

1. 单个类型禁用：使用 Eigen::DontAlign 模板参数

   Eigen::Matrix<double, 4, 1, Eigen::DontAlign> v;
   

2. 降低静态对齐阈值：定义 EIGEN_MAX_STATIC_ALIGN_BYTES 为 0（禁用所有 16 字节及以上静态对齐）或 16（仅禁用 32/64 字节对齐）。注意这会破坏 ABI 兼容性。

3. 完全禁用向量化：同时定义 EIGEN_DONT_VECTORIZE 和 EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT。这保留了对齐机制（维持 ABI 兼容），但关闭了向量化。

如何验证代码的对齐安全性？

代码有单元测试时，可通过链接仅返回 8 字节对齐缓冲区的自定义 malloc 库暴露对齐缺陷。编译时定义 EIGEN_MALLOC_ALREADY_ALIGNED=0。

总结：最佳实践

9.3 错误处理机制

Eigen提供多层错误处理机制，从编译时检查到运行时验证。

编译时检查（静态断言）

// 维度检查在编译时进行
Eigen::Matrix<double, 3, 4> A;
Eigen::Matrix<double, 4, 3> B;
auto C = A * B;  // 编译通过：3x4 * 4x3 = 3x3
// auto D = A + B;  // 编译错误：维度不匹配

运行时检查（调试模式）

// 调试模式下，Eigen会检查运行时错误
Eigen::MatrixXd X(3, 3), Y(4, 4);
// auto Z = X + Y;  // 调试模式触发断言

// 禁用调试检查（发布模式）
// 编译选项：-DEIGEN_NO_DEBUG -DNDEBUG

分解运算状态检查

// 所有分解类都提供info()方法
Eigen::LLT<Eigen::MatrixXd> llt(A);
if (llt.info() == Eigen::NumericalIssue) {
    // 矩阵非正定
    std::cerr << "Cholesky分解失败：矩阵非正定\n";
}

Eigen::BDCSVD<Eigen::MatrixXd, Eigen::ComputeThinU | Eigen::ComputeThinV> svd(A);
if (svd.info() != Eigen::Success) {
    std::cerr << "SVD分解失败\n";
}

数值检查方法

Eigen::MatrixXd A = Eigen::MatrixXd::Random(100, 100);

// 检查NaN
if (A.hasNaN()) {
    std::cerr << "矩阵包含NaN\n";
}

// 检查Inf
if (!A.allFinite()) {
    std::cerr << "矩阵包含Inf\n";
}

// 检查条件数（判断是否病态）
Eigen::JacobiSVD<Eigen::MatrixXd> svd(A);
const auto singular_values = svd.singularValues();
double sigma_max = singular_values(0);
double sigma_min = singular_values(singular_values.size() - 1);

if (sigma_min <= 1e-15) {
    std::cerr << "矩阵奇异或接近奇异，最小奇异值 = " << sigma_min << "\n";
} else {
    double cond = sigma_max / sigma_min;
    if (cond > 1e10) {
        std::cerr << "矩阵病态，条件数 = " << cond << "\n";
    }
}

错误处理最佳实践

9.4 调试技巧

打印矩阵信息

#include <Eigen/Dense>
#include <iostream>

template<typename Derived>
void print_matrix_info(const Eigen::MatrixBase<Derived>& mat, const std::string& name) {
    std::cout << "=== " << name << " ===\n";
    std::cout << "尺寸: " << mat.rows() << " x " << mat.cols() << "\n";
    std::cout << "内容:\n" << mat << "\n";
    std::cout << "范数: " << mat.norm() << "\n";
    std::cout << "最小值: " << mat.minCoeff() << "\n";
    std::cout << "最大值: " << mat.maxCoeff() << "\n";
    std::cout << "==================\n\n";
}

// 使用
Eigen::MatrixXd A = Eigen::MatrixXd::Random(5, 5);
print_matrix_info(A, "矩阵A");

检查数值稳定性

bool check_numerical_stability(const Eigen::MatrixXd& A) {
    // 检查NaN
    if (A.hasNaN()) {
        std::cerr << "错误：矩阵包含NaN\n";
        return false;
    }
    
    // 检查Inf
    if (!A.allFinite()) {
        std::cerr << "错误：矩阵包含Inf\n";
        return false;
    }
    
    // 检查数值范围
    double max_abs = A.cwiseAbs().maxCoeff();
    if (max_abs > 1e100) {
        std::cerr << "警告：矩阵元素过大 (" << max_abs << ")\n";
    }
    if (max_abs < 1e-100 && max_abs > 0) {
        std::cerr << "警告：矩阵元素过小 (" << max_abs << ")\n";
    }
    
    return true;
}

9.5 性能瓶颈定位

#include <Eigen/Dense>
#include <chrono>
#include <map>
#include <string>

// 简单的性能分析器
class Profiler {
    static inline std::map<std::string, std::pair<int, double>> stats;
public:
    static void record(const std::string& name, double ms) {
        stats[name].first++;
        stats[name].second += ms;
    }
    
    static void report() {
        std::cout << "\n=== 性能报告 ===\n";
        for (const auto& [name, data] : stats) {
            auto [count, total] = data;
            std::cout << name << ": " << count << " 次, 平均 " 
                      << total / count << " ms, 总计 " << total << " ms\n";
        }
    }
};

#define PROFILE_SCOPE(name) \
    auto _start = std::chrono::high_resolution_clock::now(); \
    struct _profiler_##__LINE__ { \
        const char* _name; \
        ~_profiler_##__LINE__() { \
            auto _end = std::chrono::high_resolution_clock::now(); \
            auto _dur = std::chrono::duration_cast<std::chrono::microseconds>(_end - _start).count() / 1000.0; \
            Profiler::record(_name, _dur); \
        } \
    } _prof_##__LINE__{name};

// 使用示例
void my_function() {
    PROFILE_SCOPE("my_function");
    
    PROFILE_SCOPE("矩阵创建");
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(1000, 1000);
    
    PROFILE_SCOPE("矩阵运算");
    Eigen::MatrixXd B = A * A.transpose();
}

9.6 错误处理最佳实践

分解运算的错误检查

Eigen的分解类提供info()方法检查计算是否成功：

#include <Eigen/Dense>
#include <iostream>

int main() {
    Eigen::Matrix3d A;
    A << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;
    
    // 特征值分解
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(A.transpose() * A);
    if (eigensolver.info() != Eigen::Success) {
        std::cerr << "特征值分解失败!\n";
        return -1;
    }
    std::cout << "特征值:\n" << eigensolver.eigenvalues() << "\n";
    
    // LU分解求解线性方程
    Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 1);
    Eigen::FullPivLU<Eigen::MatrixXd> lu(A);
    if (!lu.isInvertible()) {
        std::cerr << "矩阵不可逆!\n";
        return -1;
    }
    Eigen::VectorXd x = lu.solve(B);
    if (lu.info() != Eigen::Success) {
        std::cerr << "求解失败!\n";
        return -1;
    }
    
    // SVD分解：
    // - BDCSVD：通常更适合较大矩阵
    // - JacobiSVD：通常更适合较小矩阵或更强调精度的场景
    // Eigen 5.x 推荐使用编译时模板参数指定 thin/full U/V 选项
    Eigen::BDCSVD<Eigen::MatrixXd, Eigen::ComputeThinU | Eigen::ComputeThinV> svd(A);
    if (svd.info() != Eigen::Success) {
        std::cerr << "SVD分解失败!\n";
        return -1;
    }
    
    return 0;
}

安全的矩阵求逆

#include <Eigen/Dense>
#include <stdexcept>

Eigen::MatrixXd safe_inverse(const Eigen::MatrixXd& A, double threshold = 1e-10) {
    if (A.rows() != A.cols()) {
        throw std::invalid_argument("矩阵必须是方阵");
    }
    
    // 这里计算 thin U/V，是因为后续需要调用 solve() 来构造稳定的逆或伪逆
    // Eigen 5.x 推荐使用编译时模板参数指定 thin/full U/V 选项
    Eigen::BDCSVD<Eigen::MatrixXd, Eigen::ComputeThinU | Eigen::ComputeThinV> svd(A);
    double cond = svd.singularValues()(0) / 
                  svd.singularValues()(svd.singularValues().size() - 1);
    
    if (cond > 1.0 / threshold) {
        std::cerr << "警告：矩阵病态，条件数 = " << cond << "\n";
        return svd.solve(Eigen::MatrixXd::Identity(A.rows(), A.cols()));
    }
    
    return A.inverse();
}

数值稳定性检查

bool check_matrix_valid(const Eigen::MatrixXd& A, const std::string& name = "矩阵") {
    if (A.hasNaN()) {
        std::cerr << "错误：" << name << "包含NaN值\n";
        return false;
    }
    if (!A.allFinite()) {
        std::cerr << "错误：" << name << "包含Inf值\n";
        return false;
    }
    return true;
}

// 使用示例
Eigen::MatrixXd A = Eigen::MatrixXd::Random(3, 3);
Eigen::MatrixXd B = A.inverse();
if (!check_matrix_valid(B, "逆矩阵")) {
    // 处理错误
}

常见问题 (FAQ)

Q: 调试模式下程序很慢，如何优化？

A: 发布模式编译：

g++ -O3 -DNDEBUG -DEIGEN_NO_DEBUG ...

Q: 如何启用Eigen的详细调试信息？

A:

#define EIGEN_RUNTIME_NO_MALLOC  // 检测堆分配
#define EIGEN_INITIALIZE_MATRICES_BY_NAN  // 初始化矩阵为NaN

Q: 程序在矩阵运算时崩溃，如何定位？

A:

1. 确保使用调试模式编译（不加-DNDEBUG）
2. 使用Eigen::Matrix的.data()检查原始内存
3. 使用Valgrind检测内存错误：

valgrind --tool=memcheck --leak-check=full ./myapp

练习题

1. 调试练习: 编写一个包含常见错误的程序，练习使用调试技巧定位问题。
2. 性能分析: 使用提供的Profiler类分析一个复杂算法的性能瓶颈。

对应官方文档：Unaligned array assertion | Structures Having Eigen Members | Using STL Containers with Eigen

十、实战案例篇

通过多个独立案例展示 Eigen 在统计学习、几何计算、状态估计和图像处理中的实际用法。

阅读说明

1. 每节示例均为独立程序，包含各自的 main()，不拼接为单个源文件编译。
2. 示例默认按 C++17 编写（部分代码使用结构化绑定等语法）；需统一到 C++14 时，手动改写对应示例。
3. 偏向“教学演示”而不是“生产级实现”，重点是展示 Eigen 的用法与建模思路。
4. 阅读时结合前文相关章节：矩阵基础、线性代数、几何变换与调试篇。
5. 不同案例中的数据组织方式：有的示例使用“每行一个样本”，有的使用“每列一个点”，还有的直接把矩阵当作二维网格或图像。

10.0 案例路线图

按需选择：

• 侧重数据分析：先读 10.1 和 10.2
• 侧重几何/机器人：先读 10.3
• 侧重控制与估计：先读 10.4
• 侧重矩阵操作与图像卷积：先读 10.5

10.1 主成分分析（PCA）

前置：

• 矩阵与向量基本运算
• 协方差矩阵
• 特征值分解 / SelfAdjointEigenSolver

数据约定：

• 本节示例采用“每行一个样本”的数据组织方式

**编译**：

• 推荐使用 C++17
• 本节代码为独立示例程序

#include <Eigen/Dense>
#include <iostream>
#include <vector>

/**
 * PCA实现
 * @param X 数据矩阵（每行一个样本）
 * @param n_components 保留的主成分数量
 * @return 包含投影矩阵和方差解释率的pair
 */
std::pair<Eigen::MatrixXd, Eigen::VectorXd> pca(
    const Eigen::MatrixXd& X, 
    int n_components) {
    
    // 1. 数据中心化
    Eigen::VectorXd mean = X.colwise().mean();
    Eigen::MatrixXd X_centered = X.rowwise() - mean.transpose();
    
    // 2. 计算协方差矩阵
    Eigen::MatrixXd cov = (X_centered.transpose() * X_centered) / (X.rows() - 1);
    
    // 3. 特征值分解
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> solver(cov);
    if (solver.info() != Eigen::Success) {
        throw std::runtime_error("PCA特征值分解失败");
    }
    
    // 4. 按特征值降序排序
    // SelfAdjointEigenSolver返回升序排列的特征值
    // 需要降序时，同时反转特征值和对应的特征向量列
    Eigen::VectorXd eigenvalues = solver.eigenvalues().reverse();
    Eigen::MatrixXd eigenvectors = solver.eigenvectors().rowwise().reverse();
    
    // 5. 选择前n_components个主成分
    Eigen::MatrixXd components = eigenvectors.leftCols(n_components);
    
    // 6. 计算方差解释率
    double total_var = eigenvalues.sum();
    Eigen::VectorXd explained_variance_ratio = eigenvalues.head(n_components) / total_var;
    
    return {components, explained_variance_ratio};
}

// 使用示例
int main() {
    // 生成示例数据
    Eigen::MatrixXd X(100, 3);
    X.setRandom();
    // 添加一些相关性
    X.col(1) = X.col(0) * 0.8 + X.col(1) * 0.2;
    
    auto [components, variance_ratio] = pca(X, 2);
    
    std::cout << "主成分:\n" << components << "\n\n";
    std::cout << "方差解释率: " << variance_ratio.transpose() << "\n";
    std::cout << "累计解释率: " << variance_ratio.sum() * 100 << "%\n";
    
    // 数据投影
    Eigen::VectorXd mean = X.colwise().mean();
    Eigen::MatrixXd X_projected = (X.rowwise() - mean.transpose()) * components;
    std::cout << "\n投影后数据维度: " << X_projected.rows() << " x " << X_projected.cols() << "\n";
    
    return 0;
}

10.2 线性回归与多项式拟合

前置：

• 矩阵乘法与转置
• 最小二乘问题
• QR 分解
• CompleteOrthogonalDecomposition

数据约定：

• 线性回归部分采用“每行一个样本”的数据组织方式
• 多项式拟合部分使用一维向量构造范德蒙德矩阵

**编译**：

• 推荐使用 C++17
• 本节代码为独立示例程序

#include <Eigen/Dense>
#include <iostream>
#include <vector>

/**
 * 多元线性回归
 * 求解: y = X * beta + epsilon
 * 最小化: ||y - X*beta||^2
 */
class LinearRegression {
public:
    Eigen::VectorXd coefficients;
    double intercept = 0;
    
    void fit(const Eigen::MatrixXd& X, const Eigen::VectorXd& y) {
        // 添加偏置列（全1）
        Eigen::MatrixXd X_with_bias(X.rows(), X.cols() + 1);
        X_with_bias.col(0).setOnes();
        X_with_bias.rightCols(X.cols()) = X;
        
        // 使用列主元QR求解：作为常见、稳妥的默认选择
        // 如果问题可能秩亏，或需要更稳妥地处理欠定/最小范数解，
        // 可以考虑 completeOrthogonalDecomposition()
        coefficients = X_with_bias.colPivHouseholderQr().solve(y);
        
        intercept = coefficients(0);
        coefficients = coefficients.tail(X.cols());
    }
    
    Eigen::VectorXd predict(const Eigen::MatrixXd& X) const {
        return X * coefficients + Eigen::VectorXd::Constant(X.rows(), intercept);
    }
    
    double score(const Eigen::MatrixXd& X, const Eigen::VectorXd& y) const {
        Eigen::VectorXd y_pred = predict(X);
        double ss_res = (y - y_pred).squaredNorm();
        double ss_tot = (y - y.mean() * Eigen::VectorXd::Ones(y.size())).squaredNorm();
        return 1.0 - ss_res / ss_tot;  // R^2 score
    }
};

/**
 * 多项式拟合
 * 拟合: y = a0 + a1*x + a2*x^2 + ... + an*x^n
 */
Eigen::VectorXd polynomial_fit(const Eigen::VectorXd& x, 
                               const Eigen::VectorXd& y, 
                               int degree) {
    // 构建范德蒙德矩阵
    Eigen::MatrixXd V(x.size(), degree + 1);
    for (int i = 0; i <= degree; ++i) {
        V.col(i) = x.array().pow(i);
    }
    
    // 求解最小二乘问题
    // 这里继续使用列主元QR；若更关注秩亏情形的稳健性，
    // 可改为 V.completeOrthogonalDecomposition().solve(y)
    return V.colPivHouseholderQr().solve(y);
}

// 使用示例
int main() {
    // 线性回归示例
    std::cout << "=== 线性回归示例 ===\n";
    
    // 生成数据: y = 2*x1 + 3*x2 + 1 + noise
    Eigen::MatrixXd X(100, 2);
    X.setRandom();
    Eigen::VectorXd y = 2 * X.col(0) + 3 * X.col(1) + 
                        Eigen::VectorXd::Ones(100) + 
                        0.1 * Eigen::VectorXd::Random(100);
    
    LinearRegression model;
    model.fit(X, y);
    
    std::cout << "系数: " << model.coefficients.transpose() << "\n";
    std::cout << "截距: " << model.intercept << "\n";
    std::cout << "R^2 score: " << model.score(X, y) << "\n";
    
    // 多项式拟合示例
    std::cout << "\n=== 多项式拟合示例 ===\n";
    
    Eigen::VectorXd x_vals(10);
    x_vals << 0, 1, 2, 3, 4, 5, 6, 7, 8, 9;
    Eigen::VectorXd y_vals = x_vals.array().square() + 2 * x_vals.array() + 1;
    
    Eigen::VectorXd coeffs = polynomial_fit(x_vals, y_vals, 2);
    std::cout << "多项式系数 (从低次到高次): " << coeffs.transpose() << "\n";
    
    return 0;
}

10.3 三维刚体变换与ICP算法

前置：

• 刚体变换、旋转矩阵与四元数
• 奇异值分解（SVD）
• 点云配准的基本概念

数据约定：

• 本节点云采用“每列一个点”的组织方式，即矩阵形状为 3 x N

实现说明：

• 本节中的 ICP 为教学简化版
• 最近邻搜索使用暴力法，仅用于展示 Eigen 在线性代数部分的用法
• 生产环境通常需要 KD-tree 等更高效的数据结构
• 本节重点在于说明：如何用 Eigen 表达刚体变换、去中心化、协方差矩阵与 SVD 求解

**编译**：

• 推荐使用 C++17
• 本节代码为独立示例程序

#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <iostream>
#include <vector>

using namespace Eigen;

/**
 * 计算两组点之间的刚性变换（SVD方法）
 * 求解 R, t 使得: target = R * source + t
 */
std::pair<Matrix3d, Vector3d> rigid_transform_3D(
    const MatrixXd& source, 
    const MatrixXd& target) {
    
    // 计算质心
    Vector3d centroid_src = source.rowwise().mean();
    Vector3d centroid_tgt = target.rowwise().mean();
    
    // 去质心
    MatrixXd src_demean = source.colwise() - centroid_src;
    MatrixXd tgt_demean = target.colwise() - centroid_tgt;
    
    // 计算协方差矩阵
    Matrix3d H = src_demean * tgt_demean.transpose();
    
    // SVD分解（这里直接使用编译时模板参数写法）
    JacobiSVD<Matrix3d, ComputeFullU | ComputeFullV> svd(H);
    Matrix3d U = svd.matrixU();
    Matrix3d V = svd.matrixV();
    
    // 计算旋转矩阵
    Matrix3d R = V * U.transpose();
    
    // 处理反射情况（确保R是旋转矩阵而非反射矩阵）
    // SVD 求刚体变换时，UV^T 可能产生反射矩阵（det = -1），
    // 此时需将最小奇异值对应的 V 列取反，修正为旋转矩阵
    if (R.determinant() < 0) {
        V.col(2) *= -1;
        R = V * U.transpose();
    }
    
    // 计算平移
    Vector3d t = centroid_tgt - R * centroid_src;
    
    return {R, t};
}

/**
 * 简单ICP（迭代最近点）算法实现
 */
class SimpleICP {
public:
    struct Result {
        Matrix3d R;
        Vector3d t;
        double rmse;
        int iterations;
    };
    
    Result align(const MatrixXd& source, const MatrixXd& target, 
                 int max_iterations = 50, double tolerance = 1e-6) {
        MatrixXd src = source;
        Matrix3d R_total = Matrix3d::Identity();
        Vector3d t_total = Vector3d::Zero();
        
        double prev_error = std::numeric_limits<double>::max();
        
        for (int iter = 0; iter < max_iterations; ++iter) {
            // 步骤1: 找到最近点对（简化版本：假设已知对应关系）
            // 实际应用中需要KD树等数据结构
            MatrixXd matched_target = findNearestNeighbors(src, target);
            
            // 步骤2: 计算刚性变换
            auto [R, t] = rigid_transform_3D(src, matched_target);
            
            // 步骤3: 应用变换
            src = (R * src).colwise() + t;
            
            // 累计变换
            R_total = R * R_total;
            t_total = R * t_total + t;
            
            // 计算误差
            double error = (src - matched_target).colwise().norm().mean();
            
            if (std::abs(prev_error - error) < tolerance) {
                return {R_total, t_total, error, iter + 1};
            }
            prev_error = error;
        }
        
        return {R_total, t_total, prev_error, max_iterations};
    }
    
private:
    MatrixXd findNearestNeighbors(const MatrixXd& src, const MatrixXd& tgt) {
        // 暴力搜索最近邻点（仅用于演示）
        // 生产环境推荐使用nanoflann或FLANN库加速
        MatrixXd result(3, src.cols());
        for (int i = 0; i < src.cols(); ++i) {
            double min_dist = std::numeric_limits<double>::max();
            int min_idx = 0;
            for (int j = 0; j < tgt.cols(); ++j) {
                double dist = (src.col(i) - tgt.col(j)).norm();
                if (dist < min_dist) {
                    min_dist = dist;
                    min_idx = j;
                }
            }
            result.col(i) = tgt.col(min_idx);
        }
        return result;
    }
};

// 使用示例
int main() {
    // 创建源点云
    MatrixXd source(3, 100);
    source.setRandom();
    
    // 应用已知变换生成目标点云
    AngleAxisd rotation(EIGEN_PI / 6, Vector3d::UnitZ());
    Vector3d translation(1, 2, 3);
    
    MatrixXd target = (rotation.toRotationMatrix() * source).colwise() + translation;
    
    // 添加噪声
    target += 0.01 * MatrixXd::Random(3, 100);
    
    // 使用ICP恢复变换
    SimpleICP icp;
    auto result = icp.align(source, target);
    
    std::cout << "=== ICP配准结果 ===\n";
    std::cout << "迭代次数: " << result.iterations << "\n";
    std::cout << "最终RMSE: " << result.rmse << "\n\n";
    
    std::cout << "估计的旋转矩阵:\n" << result.R << "\n\n";
    std::cout << "估计的平移向量: " << result.t.transpose() << "\n\n";
    
    // 验证
    std::cout << "真实旋转:\n" << rotation.toRotationMatrix() << "\n";
    std::cout << "真实平移: " << translation.transpose() << "\n";
    
    return 0;
}

10.4 卡尔曼滤波器

前置：

• 矩阵乘法、转置与逆
• 线性状态空间模型
• 协方差与高斯噪声的基本概念

实现说明：

• 本节使用的是线性卡尔曼滤波教学示例
• 更新步骤中使用了 LDLT 分解求解代替显式求逆（与教程第 5 章的 solve 优先原则一致），对初学者也更安全
• 本节重点不是构建工业级滤波器，而是展示如何用 Eigen 组织状态向量、系统矩阵、协方差矩阵与更新公式

**编译**：

• 推荐使用 C++17
• 本节代码为独立示例程序

#include <Eigen/Dense>
#include <iostream>

using namespace Eigen;

/**
 * 线性卡尔曼滤波器实现
 * 
 * 状态方程: x_k = F * x_{k-1} + B * u_k + w_k
 * 观测方程: z_k = H * x_k + v_k
 * 
 * w_k ~ N(0, Q)  过程噪声
 * v_k ~ N(0, R)  观测噪声
 */
class KalmanFilter {
public:
    KalmanFilter(int state_dim, int measurement_dim) :
        n(state_dim), m(measurement_dim) {
        
        x = VectorXd::Zero(n);      // 状态估计
        P = MatrixXd::Identity(n, n); // 协方差矩阵
        F = MatrixXd::Identity(n, n); // 状态转移矩阵
        B = MatrixXd::Zero(n, 1);   // 控制矩阵（n×1，默认为零——不使用控制输入）
        H = MatrixXd::Zero(m, n);   // 观测矩阵
        Q = MatrixXd::Identity(n, n) * 0.01; // 过程噪声
        R = MatrixXd::Identity(m, m) * 0.1;  // 观测噪声
    }
    
    // 预测步骤
    void predict(const VectorXd& u = VectorXd()) {
        // 状态预测
        x = F * x;
        if (u.size() > 0) {
            x += B * u;
        }
        
        // 协方差预测
        P = F * P * F.transpose() + Q;
    }
    
    // 更新步骤
    void update(const VectorXd& z) {
        // 计算卡尔曼增益
        MatrixXd S = H * P * H.transpose() + R;
        // K = P * H^T * S^{-1}，等价于求解 K^T = S^{-1} * (H * P)
        // 使用 LDLT 分解代替显式求逆（S 为对称正定，数值更稳定）
        MatrixXd K = S.ldlt().solve(H * P).transpose();
        
        // 状态更新
        VectorXd y = z - H * x;  // 观测残差
        x = x + K * y;
        
        // 协方差更新（Joseph形式，数值更稳定）
        MatrixXd I_KH = MatrixXd::Identity(n, n) - K * H;
        P = I_KH * P * I_KH.transpose() + K * R * K.transpose();
    }
    
    VectorXd getState() const { return x; }
    MatrixXd getCovariance() const { return P; }
    
public:
    MatrixXd F, B, H, Q, R;
    
private:
    int n, m;
    VectorXd x;
    MatrixXd P;
};

// 使用示例：一维运动跟踪
int main() {
    // 状态: [位置, 速度]
    // 观测: [位置]
    KalmanFilter kf(2, 1);
    
    // 设置系统矩阵
    double dt = 0.1;  // 时间步长
    kf.F << 1, dt,
            0, 1;
    kf.H << 1, 0;
    kf.Q << 0.01, 0,
            0, 0.001;
    kf.R << 0.1;
    
    // 模拟真实轨迹和观测
    double true_pos = 0;
    double true_vel = 2;
    
    std::cout << "=== 卡尔曼滤波跟踪 ===\n";
    std::cout << "时间\t观测位置\t估计位置\t估计速度\n";
    
    for (int i = 0; i < 20; ++i) {
        // 真实状态更新
        true_pos += true_vel * dt;
        
        // 生成带噪声的观测
        double measurement = true_pos + 0.1 * ((double)rand() / RAND_MAX - 0.5);
        
        // 卡尔曼滤波
        kf.predict();
        kf.update(VectorXd::Constant(1, measurement));
        
        VectorXd state = kf.getState();
        std::cout << i * dt << "\t" << measurement << "\t\t" 
                  << state(0) << "\t\t" << state(1) << "\n";
    }
    
    return 0;
}

10.5 图像处理 - 高斯滤波与边缘检测

前置：

• 矩阵块操作
• Array / Matrix 的逐元素运算
• 二维卷积与 Sobel 算子的基本概念

数据约定：

• 本节直接使用 MatrixXd 表示灰度图像
• 图像矩阵的每个元素代表一个像素值

实现说明：

• 本节实现为简化教学版本
• 未包含 padding、边界处理、SIMD 优化和实际图像 I/O
• 重点在于展示如何使用 Eigen 进行卷积和梯度计算
• 本节尤其适合和第四章一起阅读，以体会块操作与逐元素运算在二维数据上的配合方式

**编译**：

• 推荐使用 C++17
• 本节代码为独立示例程序

#include <Eigen/Dense>
#include <iostream>

using namespace Eigen;

/**
 * 创建高斯核
 */
MatrixXd gaussian_kernel(int size, double sigma) {
    MatrixXd kernel(size, size);
    int center = size / 2;
    double sum = 0;
    
    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j) {
            int x = i - center;
            int y = j - center;
            kernel(i, j) = std::exp(-(x*x + y*y) / (2 * sigma * sigma));
            sum += kernel(i, j);
        }
    }
    
    return kernel / sum;  // 归一化
}

/**
 * 2D卷积（简化实现，无填充）
 */
MatrixXd convolve2d(const MatrixXd& input, const MatrixXd& kernel) {
    int output_rows = input.rows() - kernel.rows() + 1;
    int output_cols = input.cols() - kernel.cols() + 1;
    
    MatrixXd output(output_rows, output_cols);
    
    for (int i = 0; i < output_rows; ++i) {
        for (int j = 0; j < output_cols; ++j) {
            output(i, j) = (input.block(i, j, kernel.rows(), kernel.cols()).array() 
                           * kernel.array()).sum();
        }
    }
    
    return output;
}

/**
 * Sobel边缘检测
 */
std::pair<MatrixXd, MatrixXd> sobel_edge_detection(const MatrixXd& image) {
    // Sobel算子
    MatrixXd sobel_x(3, 3);
    sobel_x << -1, 0, 1,
               -2, 0, 2,
               -1, 0, 1;
    
    MatrixXd sobel_y(3, 3);
    sobel_y << -1, -2, -1,
                0,  0,  0,
                1,  2,  1;
    
    MatrixXd grad_x = convolve2d(image, sobel_x);
    MatrixXd grad_y = convolve2d(image, sobel_y);
    
    return {grad_x, grad_y};
}

// 使用示例
int main() {
    // 创建测试图像（简单的边缘）
    MatrixXd image(10, 10);
    image.setZero();
    image.rightCols(5).setOnes();  // 右半边为1
    
    std::cout << "=== 原始图像 ===\n" << image << "\n\n";
    
    // 高斯滤波
    MatrixXd gaussian = gaussian_kernel(3, 1.0);
    std::cout << "=== 高斯核 ===\n" << gaussian << "\n\n";
    
    MatrixXd blurred = convolve2d(image, gaussian);
    std::cout << "=== 高斯滤波后 ===\n" << blurred << "\n\n";
    
    // Sobel边缘检测
    auto [grad_x, grad_y] = sobel_edge_detection(image);
    MatrixXd magnitude = (grad_x.array().square() + grad_y.array().square()).sqrt();
    
    std::cout << "=== 边缘强度 ===\n" << magnitude << "\n";
    
    return 0;
}

附录：迁移变更汇总（Eigen 3.4 → 5.0/5.0.1）

本附录汇总了从 Eigen 3.4 迁移到 Eigen 5.0/5.0.1 时需要注意的所有变化，按主题组织，标注每个变更的类型（破坏性 / 行为语义变化 / 推荐用法变化）。

版本号语义变更

Eigen 5.0 开始采用语义化版本控制（Semantic Versioning）。

• 旧习惯：WORLD.MAJOR.MINOR
• 新习惯：MAJOR.MINOR.PATCH
• 历史兼容宏 EIGEN_WORLD_VERSION 仍保持为 3

原则：对外写 Eigen 5.0.1，不把 5.0.1 理解成 3.5.1。

输出当前版本：

#include <Eigen/Core>
#include <iostream>

int main() {
    std::cout << "Eigen version: "
              << EIGEN_MAJOR_VERSION << "."
              << EIGEN_MINOR_VERSION << "."
              << EIGEN_PATCH_VERSION << "\n";
}

类型：信息性变更

C++ 标准要求

建议：编译时至少使用 -std=c++14，教程与工程实践中更推荐 -std=c++17。

set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)

类型：破坏性变更

Eigen::all / Eigen::last 命名空间迁移

出于与其他项目的命名冲突考虑，Eigen 5.0 将 all 和 last 从 Eigen:: 移至 Eigen::placeholders::。

旧写法（3.4）：

using namespace Eigen;
Eigen::VectorXd v2 = v(seq(0, last, 2));

Eigen 5.0 推荐写法：

using namespace Eigen;
using namespace Eigen::placeholders;

Eigen::VectorXd v2 = v(seq(0, last, 2));

完整的 placeholders 命名空间包含：all、last、lastp1、end、lastN。

此外，Eigen 还提供了 Eigen::indexing 子命名空间作为快捷方式：

using namespace Eigen::indexing;  // 导入所有 slicing 相关符号

类型：破坏性变更

SVD 计算选项：运行时 → 编译时

Eigen 5.0 中，SVD 的 thin/full U/V 运行时选项被弃用，推荐改用编译时模板参数。

旧写法：

Eigen::BDCSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);

推荐写法：

Eigen::BDCSVD<Eigen::MatrixXd, Eigen::ComputeThinU | Eigen::ComputeThinV> svd(A);
// JacobiSVD 同理：
Eigen::JacobiSVD<Eigen::MatrixXd, Eigen::ComputeThinU | Eigen::ComputeThinV> svd2(A);

说明：旧代码短期内可能还能工作，但应尽早迁移；新代码直接使用编译时模板参数。

类型：推荐用法变化（旧写法已弃用）

欧拉角返回值更 canonical

Eigen 5.0 调整了欧拉角的返回形式，使其更接近规范表示。

• 等价旋转的数值表示可能与旧版本不同
• 角度范围可能变化
• eulerAngles(a0, a1, a2) 的轴顺序语义没有改变

建议：

• 不要依赖某个固定角度区间作为协议
• 若用于内部旋转计算，优先使用四元数或旋转矩阵

Eigen::Matrix3d R = ...;
Eigen::Vector3d ea = R.eulerAngles(2, 1, 0);  // ZYX 顺序
// 结果应视为"当前版本下的一种等价表示"

类型：行为语义变化

随机数生成行为变化

MatrixXd::Random() 在不同 Eigen 版本中可能得到不同序列。

建议：

• 不要把 Random() 的具体输出当作回归测试基准
• 若需要可复现实验，使用自己的随机数引擎与固定种子，再写入 Eigen 容器

类型：行为语义变化

Eigen::aligned_allocator 不再继承 std::allocator

原因与标准库分配器模型变化有关。依赖旧继承关系的代码（如自定义容器封装、STL 容器适配）可能需要调整。

对于普通用户，对齐用法本身不变：

// 这仍然是正确的写法
std::vector<Eigen::Vector4d, Eigen::aligned_allocator<Eigen::Vector4d>> v;

类型：破坏性变更

Eigen BLAS 返回类型变更

Eigen BLAS 的返回类型从 int 改为 void，以提高与其他 BLAS 实现的兼容性。

影响范围：BLAS/LAPACK 接口兼容层、与外部线性代数库的集成代码。

类型：破坏性变更

LGPL 代码移除

所有 LGPL 许可代码已移除（如 ConstrainedConjugateGradient）。这些本属于“unsupported“模块，使用不广泛。

类型：破坏性变更（移除功能）

直接包含内部头文件将被阻止

例如 ../src/... 这类内部路径不再允许直接包含。

类型：破坏性变更

以下内容不是 Eigen 5.0 新增

为避免误解，下面这些能力在 Eigen 3.4 中已经存在：

• slicing / indexing API
• C++ 模板别名（如 Vector<float, 3> 写法）
• bfloat16
• 在 CUDA kernel 中使用固定大小 Eigen 类型
• SVDBase::info()
• 若干稀疏矩阵与子块/视图相关接口

迁移检查清单

从 Eigen 3.4 迁移到 5.0/5.0.1 时，逐项检查：

• 编译选项已升级到 C++14 或更高
• 不再使用 Eigen::all / Eigen::last（改用 Eigen::placeholders::）
• 不再使用 SVD 的运行时 thin/full 选项（改用编译时模板参数）
• 不依赖 Random() 的固定输出序列
• 不依赖某种固定形式的欧拉角返回值
• 不直接包含 Eigen 内部头文件
• 不依赖 aligned_allocator 的旧继承关系

总结

升级到 Eigen 5.0 时，最值得关注的不是“多了哪些新函数“，而是：

1. 代码是否满足 C++14+
2. 是否依赖旧版 all/last 命名空间
3. 是否仍在使用旧式 SVD 运行时选项
4. 是否依赖旧的欧拉角或随机数行为
5. 是否有 STL / allocator / BLAS 兼容层代码需要适配

参考：Eigen 5.0 CHANGELOG | 开始学习

附录：API 快速索引

C. API 快速索引

矩阵创建与初始化

算术运算

块操作

矩阵分解与求解

几何变换

Array 类（逐元素运算）

D. 常见错误速查表

E. 快速参考表

常用类型别名

常用函数

F. 编译选项速查

# 说明：Eigen 5.0.x 的最低要求是 C++14，本教程中的示例默认按 C++17 组织

# 基础编译
g++ -std=c++17 -I/path/to/eigen -O2 program.cpp

# 通用高性能编译
g++ -std=c++17 -I/path/to/eigen -O3 -march=native \
    -DEIGEN_NO_DEBUG -DNDEBUG program.cpp

# 多线程并行（详见优化篇 8.4 节）
# g++ -std=c++17 -I/path/to/eigen -O3 -march=native -fopenmp \
#     -DEIGEN_NO_DEBUG -DNDEBUG program.cpp

# 调试编译
g++ -std=c++17 -I/path/to/eigen -O0 -g program.cpp

G. 推荐学习资源

1. 官方文档: https://eigen.tuxfamily.org/
2. 官方教程: https://eigen.tuxfamily.org/dox/GettingStarted.html
3. API参考: https://eigen.tuxfamily.org/dox/modules.html
4. GitLab仓库: https://gitlab.com/libeigen/eigen

总结

本教程系统介绍了 Eigen 库的各个方面：

• 基础篇：了解 Eigen 的核心特性和设计理念
• 安装篇：掌握多平台安装和 CMake 集成方法
• 矩阵基础篇：熟练进行矩阵与向量操作
• 进阶篇：掌握块操作、广播、Map/Ref 等高级用法
• 线性代数篇：学习矩阵分解与求解的一切
• 几何变换篇：掌握旋转、四元数、仿射变换
• 稀疏矩阵篇：处理大规模稀疏计算
• 优化篇：学会性能调优技巧
• 调试篇：具备问题排查能力
• 实战篇：能够应用于实际项目

1. 熟练运用Eigen进行科学计算
2. 根据场景选择合适的算法和数据结构
3. 优化代码性能，发挥硬件潜力
4. 独立解决开发中遇到的问题

本教程基于Eigen 5.0.1版本编写，以官方文档为准。