Eigen C++线性代数运算慢于NumPy的原因及优化咨询

问题描述 #

我用GCC -O3 编译了一段Eigen的C++矩阵求逆代码，结果比等价的NumPy np.linalg.inv(a) 慢3-4倍，原本以为Eigen性能更好，肯定哪里操作错了，求优化方法。

C++测试代码 #

#include <iostream>
#include <iomanip>
#include <chrono>
#include <Eigen/LU>
#include <Eigen/Dense>

using Eigen::MatrixXd;

int main()
{
  const int max_rep = 5;
  const int N=1000;
  MatrixXd inv_m(N,N);

  std::cout << "**** Eigen Library -- matrix inversion speed test ****" << std::endl;
  std::cout << "Matrix size = " << N << std::endl;
  
  for(int i=1;i<=max_rep;i++) {
    MatrixXd m = MatrixXd::Random(N,N);
    auto start = std::chrono::system_clock::now();
    inv_m = m.inverse();
    auto stop = std::chrono::system_clock::now();
    std::cout << "Test no." << std::setw(3) << i;
    std::cout << ", CPU time = " 
          << std::chrono::duration<double>{stop - start}.count()
              << " seconds." << std::endl;
  }  
  return 0;
}

Python测试代码 #

import numpy as np
import time

for n in range(500, 5001, 500):
    a=np.random.rand(n,n)
    t=time.time();
    b=np.linalg.inv(a)
    t=time.time()-t;
    print(f'Size {n=:5d}, {t=:.2f}')

NumPy配置信息 #

In [26]: np.show_config()
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
accelerate_info:
  NOT AVAILABLE
atlas_3_10_blas_threads_info:
  NOT AVAILABLE
atlas_3_10_blas_info:
  NOT AVAILABLE
atlas_blas_threads_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
blas_info:
    libraries = ['blas', 'blas']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    include_dirs = ['/usr/local/include', '/usr/include']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    libraries = ['blas', 'blas']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    include_dirs = ['/usr/local/include', '/usr/include']
    language = c
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
  NOT AVAILABLE
openblas_clapack_info:
  NOT AVAILABLE
flame_info:
  NOT AVAILABLE
atlas_3_10_threads_info:
  NOT AVAILABLE
atlas_3_10_info:
  NOT AVAILABLE
atlas_threads_info:
  NOT AVAILABLE
atlas_info:
  NOT AVAILABLE
lapack_info:
    libraries = ['lapack', 'lapack']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    language = f77
lapack_opt_info:
    libraries = ['lapack', 'lapack', 'blas', 'blas']
    library_dirs = ['/usr/lib/x86_64-linux-gnu']
    language = c
    define_macros = [('NO_ATLAS_INFO', 1), ('HAVE_CBLAS', None)]
    include_dirs = ['/usr/local/include', '/usr/include']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_KNM,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL

优化方案 #

1. 让Eigen链接系统优化的BLAS/LAPACK库 #

NumPy默认调用系统的BLAS/LAPACK实现（即便不是MKL/OpenBLAS，系统版本也经过架构优化），而Eigen默认使用纯C++实现，大矩阵运算性能差距明显。

• 在代码开头添加宏定义，强制Eigen使用外部BLAS/LAPACK：

  #define EIGEN_USE_BLAS
  #define EIGEN_USE_LAPACKE
  #include <Eigen/LU>
  #include <Eigen/Dense>
  

• 编译时链接系统BLAS和LAPACK库：

  g++ -O3 -march=native -I/usr/include/eigen3 your_code.cpp -o eigen_inv -lblas -llapack
  

2. 添加针对性编译优化选项 #

除了-O3，补充以下选项最大化CPU利用：

• -march=native：让编译器针对你的CPU生成最优指令（匹配你的AVX2、FMA3等扩展）
• -ffast-math：启用快速数学优化（若对数值精度要求极高需谨慎）
  完整编译命令：

g++ -O3 -march=native -ffast-math -I/usr/include/eigen3 your_code.cpp -o eigen_inv -lblas -llapack

3. 消除不必要的矩阵拷贝 #

原代码每次循环都重新分配MatrixXd m，增加内存开销，可复用内存：

int main()
{
  const int max_rep = 5;
  const int N=1000;
  MatrixXd inv_m(N,N);
  MatrixXd m(N,N); // 提前分配内存

  std::cout << "**** Eigen Library -- matrix inversion speed test ****" << std::endl;
  std::cout << "Matrix size = " << N << std::endl;
  
  for(int i=1;i<=max_rep;i++) {
    m = MatrixXd::Random(N,N); // 赋值而非重新分配
    auto start = std::chrono::system_clock::now();
    inv_m = m.inverse();
    auto stop = std::chrono::system_clock::now();
    std::cout << "Test no." << std::setw(3) << i;
    std::cout << ", CPU time = " 
          << std::chrono::duration<double>{stop - start}.count()
              << " seconds." << std::endl;
  }  
  return 0;
}

也可以直接使用LU分解类避免临时对象：

Eigen::PartialPivLU<MatrixXd> lu(m);
inv_m = lu.inverse();

4. 启用多线程优化 #

Eigen支持OpenMP多线程，编译时开启即可利用多核CPU：

• 代码开头添加宏：

  #define EIGEN_USE_OPENMP
  

• 编译命令添加-fopenmp：

  g++ -O3 -march=native -ffast-math -fopenmp -I/usr/include/eigen3 your_code.cpp -o eigen_inv -lblas -llapack -lgomp
  

5. 考虑单精度运算（若允许） #

如果业务场景不需要双精度，改用MatrixXf代替MatrixXd，单精度运算速度远快于双精度，能进一步拉开性能差距。

内容的提问来源于stack exchange，提问作者Mico Mrkaic