C++ Eigen库详解:高性能线性代数计算

次浏览

概述

Eigen 是一个高质量的 C++ 模板库,用于线性代数运算、矩阵和向量操作、数值计算及相关算法。它具有以下特点:

特性 说明
高性能 表达式模板、惰性求值、SIMD 优化
易用 直观的 API,类似 MATLAB 的语法
Header-only 无需编译,只需包含头文件
跨平台 支持 Windows、Linux、macOS
开源 MPL2 许可证,商用友好

一、安装与配置

1.1 安装方式

方式一:包管理器(推荐)

1
2
3
4
5
6
7
8
# Ubuntu/Debian
sudo apt-get install libeigen3-dev

# macOS
brew install eigen

# Fedora
sudo dnf install eigen3-devel

方式二:源码安装

1
2
3
4
5
6
7
8
# 下载
git clone https://gitlab.com/libeigen/eigen.git
cd eigen

# 安装到系统目录
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX=/usr/local
sudo make install

方式三:直接包含(Header-only)

1
2
3
4
# 只需将 Eigen 目录放入项目中
wget https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz
tar -xzf eigen-3.4.0.tar.gz
# 将 eigen-3.4.0/Eigen 复制到项目中

1.2 CMake 配置

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
cmake_minimum_required(VERSION 3.10)
project(eigen_demo)

# 方式一:使用系统安装的 Eigen
find_package(Eigen3 3.4 REQUIRED)

# 方式二:使用本地 Eigen 目录
# include_directories(${CMAKE_SOURCE_DIR}/eigen)

add_executable(demo main.cpp)
target_link_libraries(demo Eigen3::Eigen)

1.3 快速开始

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
#include <iostream>
#include <Eigen/Dense>

int main() {
    // 定义一个 3x3 矩阵
    Eigen::Matrix3d m;
    m << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;
    
    std::cout << "Matrix m:\n" << m << std::endl;
    std::cout << "m * m:\n" << m * m << std::endl;
    
    return 0;
}

二、矩阵与向量基础

2.1 矩阵类型定义

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
#include <Eigen/Dense>

// 固定大小矩阵(编译时确定大小)
Eigen::Matrix2d m2d;           // 2x2 double 矩阵
Eigen::Matrix3f m3f;           // 3x3 float 矩阵
Eigen::Matrix4d m4d;           // 4x4 double 矩阵
Eigen::Matrix<double, 5, 5> m5; // 5x5 double 矩阵

// 动态大小矩阵(运行时确定大小)
Eigen::MatrixXd mxd;           // 动态大小 double 矩阵
Eigen::MatrixXf mxf;           // 动态大小 float 矩阵

// 向量(列向量)
Eigen::Vector3d v3d;           // 3维 double 列向量
Eigen::VectorXd vxd;           // 动态大小 double 列向量

// 行向量
Eigen::RowVector3d rv3d;       // 3维 double 行向量
Eigen::RowVectorXd rvxd;       // 动态大小 double 行向量

2.2 矩阵初始化

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <iostream>
#include <Eigen/Dense>

int main() {
    // 1. 默认初始化(未定义值)
    Eigen::Matrix3d m1;
    
    // 2. 零矩阵
    Eigen::Matrix3d m2 = Eigen::Matrix3d::Zero();
    
    // 3. 单位矩阵
    Eigen::Matrix3d m3 = Eigen::Matrix3d::Identity();
    
    // 4. 常数矩阵
    Eigen::Matrix3d m4 = Eigen::Matrix3d::Constant(2.5);  // 所有元素为 2.5
    
    // 5. 随机矩阵
    Eigen::Matrix3d m5 = Eigen::Matrix3d::Random();
    
    // 6. 逗号初始化
    Eigen::Matrix3d m6;
    m6 << 1, 2, 3,
          4, 5, 6,
          7, 8, 9;
    
    // 7. 从数组初始化
    double data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
    Eigen::Map<Eigen::Matrix3d> m7(data);
    
    // 8. 动态大小矩阵
    Eigen::MatrixXd m8(3, 4);  // 3行4列
    m8 = Eigen::MatrixXd::Zero(3, 4);
    
    std::cout << "Identity matrix:\n" << m3 << std::endl;
    
    return 0;
}

2.3 元素访问

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <iostream>
#include <Eigen/Dense>

int main() {
    Eigen::Matrix3d m;
    m << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;
    
    // 使用 () 运算符(从0开始)
    std::cout << "m(0,0) = " << m(0, 0) << std::endl;  // 1
    std::cout << "m(1,2) = " << m(1, 2) << std::endl;  // 6
    
    // 修改元素
    m(0, 0) = 10;
    m(2, 2) = 30;
    
    // 向量访问
    Eigen::Vector3d v(1, 2, 3);
    std::cout << "v[0] = " << v[0] << std::endl;    // 1 (使用 [])
    std::cout << "v(1) = " << v(1) << std::endl;    // 2 (使用 ())
    
    // 行列操作
    m.row(0) = Eigen::Vector3d(100, 100, 100);  // 设置第一行
    m.col(1) = Eigen::Vector3d(0, 0, 0);        // 设置第二列
    
    std::cout << "Modified m:\n" << m << std::endl;
    
    return 0;
}

三、矩阵运算

3.1 基本运算

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include <iostream>
#include <Eigen/Dense>

int main() {
    Eigen::Matrix2d A, B;
    A << 1, 2,
         3, 4;
    B << 5, 6,
         7, 8;
    
    // 加减法
    std::cout << "A + B:\n" << A + B << std::endl;
    std::cout << "A - B:\n" << A - B << std::endl;
    
    // 标量运算
    std::cout << "2 * A:\n" << 2 * A << std::endl;
    std::cout << "A / 2:\n" << A / 2.0 << std::endl;
    
    // 矩阵乘法
    std::cout << "A * B:\n" << A * B << std::endl;
    
    // 逐元素运算(使用 array)
    std::cout << "A.cwiseProduct(B):\n" << A.cwiseProduct(B) << std::endl;  // 逐元素乘法
    std::cout << "A.array() * B.array():\n" << A.array() * B.array() << std::endl;
    
    // 转置
    std::cout << "A^T:\n" << A.transpose() << std::endl;
    
    // 共轭(复数矩阵)
    // A.conjugate()
    
    // 逆矩阵
    std::cout << "A^(-1):\n" << A.inverse() << std::endl;
    
    // 行列式
    std::cout << "det(A) = " << A.determinant() << std::endl;
    
    // 迹(对角线元素和)
    std::cout << "trace(A) = " << A.trace() << std::endl;
    
    // 范数
    std::cout << "Frobenius norm: " << A.norm() << std::endl;
    std::cout << "L1 norm: " << A.lpNorm<1>() << std::endl;
    std::cout << "L-inf norm: " << A.lpNorm<Eigen::Infinity>() << std::endl;
    
    return 0;
}

3.2 向量运算

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <iostream>
#include <Eigen/Dense>

int main() {
    Eigen::Vector3d a(1, 2, 3);
    Eigen::Vector3d b(4, 5, 6);
    
    // 点积
    double dot = a.dot(b);
    std::cout << "a · b = " << dot << std::endl;  // 32
    
    // 叉积(仅适用于3维向量)
    Eigen::Vector3d cross = a.cross(b);
    std::cout << "a × b = " << cross.transpose() << std::endl;  // [-3, 6, -3]
    
    // 范数
    std::cout << "L2 norm: " << a.norm() << std::endl;
    std::cout << "Squared norm: " << a.squaredNorm() << std::endl;
    std::cout << "L1 norm: " << a.lpNorm<1>() << std::endl;
    
    // 归一化
    Eigen::Vector3d a_normalized = a.normalized();  // 返回归一化向量
    a.normalize();  // 原地归一化
    
    // 向量外积(得到矩阵)
    Eigen::Matrix3d outer = a * b.transpose();
    std::cout << "Outer product:\n" << outer << std::endl;
    
    return 0;
}

3.3 块操作

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
#include <iostream>
#include <Eigen/Dense>

int main() {
    Eigen::Matrix4d m;
    m << 1,  2,  3,  4,
         5,  6,  7,  8,
         9,  10, 11, 12,
         13, 14, 15, 16;
    
    // 动态大小块
    std::cout << "Block (2x2 starting at (1,1)):\n" 
              << m.block(1, 1, 2, 2) << std::endl;
    
    // 固定大小块(更高效)
    std::cout << "Block<2,2>:\n" << m.block<2, 2>(1, 1) << std::endl;
    
    // 行块
    std::cout << "Rows 1-2:\n" << m.middleRows(1, 2) << std::endl;
    
    // 列块
    std::cout << "Cols 1-2:\n" << m.middleCols(1, 2) << std::endl;
    
    // 角块
    std::cout << "Top-left 2x2:\n" << m.topLeftCorner(2, 2) << std::endl;
    std::cout << "Top-right 2x2:\n" << m.topRightCorner(2, 2) << std::endl;
    std::cout << "Bottom-left 2x2:\n" << m.bottomLeftCorner(2, 2) << std::endl;
    std::cout << "Bottom-right 2x2:\n" << m.bottomRightCorner(2, 2) << std::endl;
    
    // 对角线
    std::cout << "Diagonal: " << m.diagonal().transpose() << std::endl;
    
    // 修改块
    m.block<2, 2>(0, 0) = Eigen::Matrix2d::Zero();
    std::cout << "After modification:\n" << m << std::endl;
    
    return 0;
}

四、解线性方程组

4.1 直接求解

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include <iostream>
#include <Eigen/Dense>

int main() {
    // 求解 Ax = b
    Eigen::Matrix3d A;
    A << 1, 2, 3,
         4, 5, 6,
         7, 8, 10;  // 注意:系数矩阵不能奇异
    
    Eigen::Vector3d b(3, 3, 4);
    
    // 方法1:直接求逆(不推荐,效率低)
    Eigen::Vector3d x1 = A.inverse() * b;
    std::cout << "Using inverse: " << x1.transpose() << std::endl;
    
    // 方法2:LU 分解(推荐)
    Eigen::Vector3d x2 = A.lu().solve(b);
    std::cout << "Using LU: " << x2.transpose() << std::endl;
    
    // 方法3:PartialPivLU(部分主元 LU)
    Eigen::PartialPivLU<Eigen::Matrix3d> lu(A);
    Eigen::Vector3d x3 = lu.solve(b);
    std::cout << "Using PartialPivLU: " << x3.transpose() << std::endl;
    
    // 方法4:FullPivLU(完全主元 LU,更稳定)
    Eigen::FullPivLU<Eigen::Matrix3d> full_lu(A);
    Eigen::Vector3d x4 = full_lu.solve(b);
    std::cout << "Using FullPivLU: " << x4.transpose() << std::endl;
    
    return 0;
}

4.2 根据矩阵特性选择分解方法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#include <iostream>
#include <Eigen/Dense>

int main() {
    Eigen::MatrixXd A(4, 4);
    A << 4, 1, 0, 0,
         1, 4, 1, 0,
         0, 1, 4, 1,
         0, 0, 1, 4;
    Eigen::Vector4d b(1, 2, 3, 4);
    
    // 对称正定矩阵:使用 LLT(Cholesky 分解)
    // 要求:A 必须是对称正定的
    Eigen::LLT<Eigen::MatrixXd> llt(A);
    if (llt.info() == Eigen::Success) {
        Eigen::VectorXd x = llt.solve(b);
        std::cout << "LLT solution: " << x.transpose() << std::endl;
    }
    
    // 另一种:LDLT 分解(更稳定)
    Eigen::LDLT<Eigen::MatrixXd> ldlt(A);
    if (ldlt.info() == Eigen::Success) {
        Eigen::VectorXd x = ldlt.solve(b);
        std::cout << "LDLT solution: " << x.transpose() << std::endl;
    }
    
    // 非方阵:QR 分解
    Eigen::MatrixXd B(4, 3);
    B << 1, 2, 3,
         4, 5, 6,
         7, 8, 9,
         10, 11, 12;
    Eigen::Vector4d b2(1, 2, 3, 4);
    
    Eigen::VectorXd x_qr = B.colPivHouseholderQr().solve(b2);
    std::cout << "QR solution (least squares): " x_qr.transpose() << std::endl;
    
    return 0;
}

4.3 分解方法选择指南

矩阵类型 推荐分解 方法 特点
一般方阵 PartialPivLU A.lu() 快速,适度稳定
一般方阵 FullPivLU A.fullPivLu() 慢,最稳定
对称正定 LLT A.llt() 最快,要求 SPD
对称正定 LDLT A.ldlt() 快,更稳定
非方阵 ColPivHouseholderQR A.colPivHouseholderQr() 列主元 QR
非方阵 JacobiSVD A.jacobiSvd() 最稳定,最慢
稀疏矩阵 SparseLU A.sparseLU() 稀疏 LU

五、特征值与奇异值分解

5.1 特征值分解

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#include <iostream>
#include <Eigen/Dense>

int main() {
    // 对称矩阵的特征值分解
    Eigen::Matrix3d A;
    A << 2, 1, 0,
         1, 3, 1,
         0, 1, 2;
    
    // SelfAdjointEigenSolver 用于对称/Hermitian 矩阵
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigensolver(A);
    
    if (eigensolver.info() != Eigen::Success) {
        std::cerr << "Eigenvalue decomposition failed!" << std::endl;
        return -1;
    }
    
    // 特征值(从小到大排列)
    std::cout << "Eigenvalues: " << eigensolver.eigenvalues().transpose() << std::endl;
    
    // 特征向量(列向量)
    std::cout << "Eigenvectors:\n" << eigensolver.eigenvectors() << std::endl;
    
    // 验证: A * v = λ * v
    Eigen::Vector3d v0 = eigensolver.eigenvectors().col(0);
    double lambda0 = eigensolver.eigenvalues()(0);
    std::cout << "Verification A*v = " << (A * v0).transpose() << std::endl;
    std::cout << "Verification λ*v = " << (lambda0 * v0).transpose() << std::endl;
    
    // 一般矩阵的特征值分解
    Eigen::Matrix3d B;
    B << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;
    
    Eigen::EigenSolver<Eigen::Matrix3d> solver(B);
    std::cout << "\nGeneral eigenvalues: " << solver.eigenvalues().transpose() << std::endl;
    std::cout << "General eigenvectors:\n" << solver.eigenvectors() << std::endl;
    
    return 0;
}

5.2 奇异值分解(SVD)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#include <iostream>
#include <Eigen/Dense>

int main() {
    Eigen::Matrix3d A;
    A << 1, 2, 3,
         4, 5, 6,
         7, 8, 9;
    
    // JacobiSVD(最稳定)
    Eigen::JacobiSVD<Eigen::Matrix3d> svd(A, 
        Eigen::ComputeFullU | Eigen::ComputeFullV);
    
    std::cout << "Singular values: " << svd.singularValues().transpose() << std::endl;
    std::cout << "Matrix U:\n" << svd.matrixU() << std::endl;
    std::cout << "Matrix V:\n" << svd.matrixV() << std::endl;
    
    // 验证: A = U * Σ * V^T
    Eigen::Matrix3d reconstructed = svd.matrixU() * 
        svd.singularValues().asDiagonal() * svd.matrixV().transpose();
    std::cout << "Reconstructed A:\n" << reconstructed << std::endl;
    
    // BDCSVD(大规模矩阵更快)
    Eigen::BDCSVD<Eigen::MatrixXd> bdcsvd(A);
    std::cout << "BDCSVD singular values: " << bdcsvd.singularValues().transpose() << std::endl;
    
    // 计算伪逆(Moore-Penrose 伪逆)
    double threshold = 1e-6 * svd.singularValues()(0);
    Eigen::Vector3d sig_inv = (svd.singularValues().array() > threshold)
        .select(svd.singularValues().array().inverse(), 0);
    Eigen::Matrix3d pseudo_inverse = svd.matrixV() * 
        sig_inv.asDiagonal() * svd.matrixU().transpose();
    std::cout << "Pseudo-inverse:\n" << pseudo_inverse << std::endl;
    
    return 0;
}

5.3 SVD 应用:最小二乘问题

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include <iostream>
#include <Eigen/Dense>

int main() {
    // 过定方程组 Ax = b(方程数 > 未知数)
    Eigen::MatrixXd A(5, 3);
    A << 1, 0, 1,
         0, 1, 1,
         1, 1, 0,
         1, 1, 1,
         0, 0, 1;
    
    Eigen::Vector5d b(1, 2, 3, 4, 5);
    
    // 方法1:正规方程 A^T * A * x = A^T * b
    Eigen::VectorXd x1 = (A.transpose() * A).ldlt().solve(A.transpose() * b);
    std::cout << "Normal equations: " << x1.transpose() << std::endl;
    
    // 方法2:QR 分解
    Eigen::VectorXd x2 = A.colPivHouseholderQr().solve(b);
    std::cout << "QR decomposition: " << x2.transpose() << std::endl;
    
    // 方法3:SVD(最稳定)
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, Eigen::ComputeThinU | Eigen::ComputeThinV);
    Eigen::VectorXd x3 = svd.solve(b);
    std::cout << "SVD solution: " << x3.transpose() << std::endl;
    
    // 残差
    std::cout << "Residual norm: " << (A * x3 - b).norm() << std::endl;
    
    return 0;
}

六、稀疏矩阵

6.1 创建稀疏矩阵

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
#include <iostream>
#include <Eigen/Sparse>

int main() {
    // 定义稀疏矩阵类型
    typedef Eigen::SparseMatrix<double> SpMat;
    typedef Eigen::Triplet<double> T;
    
    // 创建 4x4 稀疏矩阵
    SpMat A(4, 4);
    
    // 使用三元组列表填充
    std::vector<T> triplets;
    triplets.push_back(T(0, 0, 1.0));
    triplets.push_back(T(1, 1, 2.0));
    triplets.push_back(T(2, 2, 3.0));
    triplets.push_back(T(3, 3, 4.0));
    triplets.push_back(T(0, 1, 0.5));
    triplets.push_back(T(1, 0, 0.5));
    
    A.setFromTriplets(triplets.begin(), triplets.end());
    
    std::cout << "Sparse matrix A:\n" << Eigen::MatrixXd(A) << std::endl;
    std::cout << "Non-zeros: " << A.nonZeros() << std::endl;
    
    // 遍历非零元素
    for (int k = 0; k < A.outerSize(); ++k) {
        for (SpMat::InnerIterator it(A, k); it; ++it) {
            std::cout << "A(" << it.row() << "," << it.col() << ") = " 
                      << it.value() << std::endl;
        }
    }
    
    return 0;
}

6.2 稀疏矩阵求解

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
#include <iostream>
#include <Eigen/Sparse>

int main() {
    typedef Eigen::SparseMatrix<double> SpMat;
    typedef Eigen::Triplet<double> T;
    
    // 创建三对角稀疏矩阵
    int n = 100;
    SpMat A(n, n);
    std::vector<T> triplets;
    
    for (int i = 0; i < n; i++) {
        triplets.push_back(T(i, i, 2.0));
        if (i > 0) triplets.push_back(T(i, i-1, -1.0));
        if (i < n-1) triplets.push_back(T(i, i+1, -1.0));
    }
    A.setFromTriplets(triplets.begin(), triplets.end());
    
    // 右侧向量
    Eigen::VectorXd b = Eigen::VectorXd::Constant(n, 1.0);
    
    // 使用 SparseLU 求解
    Eigen::SparseLU<SpMat> solver;
    solver.analyzePattern(A);  // 分析稀疏结构
    solver.factorize(A);       // 数值分解
    
    if (solver.info() != Eigen::Success) {
        std::cerr << "Decomposition failed!" << std::endl;
        return -1;
    }
    
    Eigen::VectorXd x = solver.solve(b);
    
    if (solver.info() != Eigen::Success) {
        std::cerr << "Solve failed!" << std::endl;
        return -1;
    }
    
    std::cout << "Solution (first 5 elements): " 
              << x.head(5).transpose() << std::endl;
    
    // 验证
    std::cout << "Residual norm: " << (A * x - b).norm() << std::endl;
    
    return 0;
}

七、几何模块

7.1 旋转与变换

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
#include <iostream>
#include <Eigen/Geometry>

int main() {
    // 旋转矩阵
    Eigen::Matrix3d rotation_matrix;
    rotation_matrix = Eigen::AngleAxisd(M_PI / 4, Eigen::Vector3d::UnitZ());
    std::cout << "Rotation matrix (45° around Z):\n" << rotation_matrix << std::endl;
    
    // 旋转向量(轴角表示)
    Eigen::AngleAxisd rotation_vector(M_PI / 3, Eigen::Vector3d(1, 1, 1).normalized());
    std::cout << "Angle: " << rotation_vector.angle() << std::endl;
    std::cout << "Axis: " << rotation_vector.axis().transpose() << std::endl;
    
    // 四元数
    Eigen::Quaterniond q(rotation_matrix);
    std::cout << "Quaternion: " q.coeffs().transpose() << std::endl;
    
    // 四元数转旋转矩阵
    Eigen::Matrix3d m = q.toRotationMatrix();
    
    // 欧拉角(ZYX顺序)
    Eigen::Vector3d euler = rotation_matrix.eulerAngles(2, 1, 0); // ZYX
    std::cout << "Euler angles (ZYX): " << euler.transpose() << std::endl;
    
    // 旋转向量
    Eigen::Vector3d v(1, 0, 0);
    Eigen::Vector3d v_rotated = rotation_matrix * v;
    std::cout << "Original: " << v.transpose() << std::endl;
    std::cout << "Rotated: " << v_rotated.transpose() << std::endl;
    
    return 0;
}

7.2 欧式变换矩阵

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <iostream>
#include <Eigen/Geometry>

int main() {
    // 3D 欧式变换矩阵(4x4 齐次坐标)
    Eigen::Isometry3d T = Eigen::Isometry3d::Identity();
    
    // 设置旋转(绕 Z 轴旋转 45°)
    T.rotate(Eigen::AngleAxisd(M_PI / 4, Eigen::Vector3d::UnitZ()));
    
    // 设置平移
    T.pretranslate(Eigen::Vector3d(1, 2, 3));
    
    std::cout << "Transform matrix:\n" << T.matrix() << std::endl;
    
    // 提取旋转和平移
    Eigen::Matrix3d R = T.rotation();
    Eigen::Vector3d t = T.translation();
    
    std::cout << "Rotation:\n" << R << std::endl;
    std::cout << "Translation: " << t.transpose() << std::endl;
    
    // 变换点
    Eigen::Vector3d p(1, 0, 0);
    Eigen::Vector3d p_transformed = T * p;
    std::cout << "Original point: " << p.transpose() << std::endl;
    std::cout << "Transformed point: " << p_transformed.transpose() << std::endl;
    
    // 变换的逆
    Eigen::Isometry3d T_inv = T.inverse();
    std::cout << "Inverse transform:\n" << T_inv.matrix() << std::endl;
    
    return 0;
}

八、性能优化技巧

8.1 使用固定大小矩阵

1
2
3
4
5
// 推荐:编译时确定大小
Eigen::Matrix3d m1;  // 快,在栈上分配

// 仅在大小未知时使用
Eigen::MatrixXd m2(3, 3);  // 慢,在堆上分配

8.2 避免临时对象

1
2
3
4
5
6
7
8
Eigen::MatrixXd A, B, C, D;

// 慢:产生临时对象
Eigen::MatrixXd result = A * B + C * D;

// 快:使用 noalias()(如果结果不与操作数重叠)
Eigen::MatrixXd result;
result.noalias() = A * B + C * D;

8.3 表达式模板

1
2
3
4
5
6
7
// Eigen 使用惰性求值,以下代码不会产生临时矩阵
Eigen::VectorXd v1, v2, v3, v4;
Eigen::VectorXd result = v1 + v2 + v3 + v4;  // 一次遍历完成

// 强制求值
auto expr = v1 + v2;  // 此时未计算
Eigen::VectorXd sum = expr;  // 此时计算

8.4 启用 SIMD

1
2
3
4
5
// CMake 中启用优化
// target_compile_options(demo PRIVATE -O3 -march=native)

// 或者使用 OpenMP
// target_link_libraries(demo Eigen3::Eigen OpenMP::OpenMP_CXX)

九、实战示例:最小二乘拟合

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
#include <iostream>
#include <vector>
#include <Eigen/Dense>

// 多项式拟合:y = a0 + a1*x + a2*x^2 + ...
Eigen::VectorXd polynomial_fit(const std::vector<double>& x_data,
                                const std::vector<double>& y_data,
                                int degree) {
    int n = x_data.size();
    Eigen::MatrixXd A(n, degree + 1);
    Eigen::VectorXd b(n);
    
    // 构建范德蒙德矩阵
    for (int i = 0; i < n; i++) {
        double x = 1.0;
        for (int j = 0; j <= degree; j++) {
            A(i, j) = x;
            x *= x_data[i];
        }
        b(i) = y_data[i];
    }
    
    // 使用 SVD 求解最小二乘
    Eigen::JacobiSVD<Eigen::MatrixXd> svd(A, 
        Eigen::ComputeThinU | Eigen::ComputeThinV);
    return svd.solve(b);
}

int main() {
    // 生成测试数据:y = 1 + 2x + 0.5x^2 + 噪声
    std::vector<double> x_data, y_data;
    for (int i = 0; i < 20; i++) {
        double x = i * 0.5;
        double y = 1 + 2*x + 0.5*x*x + ((rand() % 100) - 50) * 0.01;
        x_data.push_back(x);
        y_data.push_back(y);
    }
    
    // 二次多项式拟合
    Eigen::VectorXd coeffs = polynomial_fit(x_data, y_data, 2);
    
    std::cout << "Fitted polynomial coefficients:\n";
    std::cout << "y = " << coeffs(0) << " + " 
              << coeffs(1) << "x + " 
              << coeffs(2) << "x^2" << std::endl;
    
    return 0;
}

十、总结

10.1 常用头文件

1
2
3
4
5
6
7
#include <Eigen/Core>        // 核心模块
#include <Eigen/Dense>       // 稠密矩阵
#include <Eigen/Geometry>    // 几何模块
#include <Eigen/Sparse>      // 稀疏矩阵
#include <Eigen/Eigenvalues> // 特征值
#include <Eigen/SVD>         // 奇异值分解
#include <Eigen/QR>          // QR 分解

10.2 快速参考

操作 代码
创建矩阵 Matrix3d m; MatrixXd m(rows, cols)
初始化 m << 1,2,3; m.setZero(); m.setIdentity()
元素访问 m(i,j); m.row(i); m.col(j)
矩阵运算 m+n; m-n; m*n; m.transpose(); m.inverse()
解方程 A.lu().solve(b); A.ldlt().solve(b)
特征值 SelfAdjointEigenSolver<Matrix3d> es(A)
SVD JacobiSVD<MatrixXd> svd(A)

参考资料

使用 Hugo 构建
主题 StackJimmy 设计