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 分解方法选择指南

矩阵类型推荐分解方法特点
一般方阵PartialPivLUA.lu()快速,适度稳定
一般方阵FullPivLUA.fullPivLu()慢,最稳定
对称正定LLTA.llt()最快,要求 SPD
对称正定LDLTA.ldlt()快,更稳定
非方阵ColPivHouseholderQRA.colPivHouseholderQr()列主元 QR
非方阵JacobiSVDA.jacobiSvd()最稳定,最慢
稀疏矩阵SparseLUA.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)
SVDJacobiSVD<MatrixXd> svd(A)

参考资料

评论

留言经审核后显示,欢迎留下想法。