点击下方 卡片 ，关注 「3D视觉工坊」 公众号
选择 星标 ，干货第一时间送达

如需采购标定板，请联系

来自知乎|计算机视觉life整理， https://www.zhihu.com/question/28571059/answer/3057795309

作者 Sisyphus

欢迎来到C++数值优化的世界！Eigen能写的快和它的设计思路有关，涵盖了算法加速的几个方法。Eigen的设计思路，是把所有能优化的步骤放在编译时去优化。要想运行时变快，写Eigen代码的工程师需要显式地告诉Eigen矩阵的特性，告诉的越多，越能提供空间让编译器在编译时加速算法

比如一个矩阵乘法+

# include  Eigen:: Matrixxd Jacobian_i = Eigen:: MatrixXd:: Random(10, 10); Eigen:: Matrixxd Jacobian_j = Eigen:: Matrixxd:: Random(10, 10);
Eigen:: MatrixXd Hessian = Eigen:: Matrixxd:: Zero(10, 10);
Hessian += Jacobian_i. transpose() * Jacobian_j;

这个是一个稠密矩阵相乘。首先需要知道的是Eigen并不是一步一步地先做转置，再去乘，而是使用lazy evaluation的方法。Lazy evalution，是把计算本身尽可能放在最后做，减少内存访问。

在代码上来看，就是transpose（）函数并没有做转置，和operator*（）本身也没有做乘法，只是更新了一下flag，真正运行计算的是operator+=（）。从Eigen/src/Core/EigenBase.h可以看到，是在operator+=（）这边才真正去做内存读取和计算。

第二个可以用的优化技巧是改变内存分配的方式，使用Eigen时应该尽可能用静态内存+代替动态内存。在上述矩阵计算中，三个矩阵都是动态内存分配+，因为Eigen:：MatrixXd 是如下的缩写：

typedef MatrixXd Matrix

MatrixBase第二和第三个选项是行列的长度，有一项是Dynamic就会用动态内存分配。所以在已知矩阵大小时应尽可能声明大小，比如Matrix 。如果内存在整个程序中大小会变，但知道最大可能的大小，都可以告知Eigen，Eigen同样会选择用静态内存。

# include
Jacobian_i; Eigen:: Matrix« double, Eigen:: Dynamic, Eigen:: Dynamic, Eigen:: ColMajor, 10, 10>
Jacobian_j: Eigen:: Matrixcdouble, Eigen:: Dynamic, Eigen:: Dynamic, Eigen:: ColMajor, 10, 10>
Hessian =
Eigen:: Matrixcdouble, Eigen:: Dynamic, Eigen:: Dynamic, Eigen:: ColMajor, 10, 10>:: Zero();
Hessian += Jacobian_i. transpose()* Jacobian_j;

以上是普通稠密矩阵计算。如果矩阵本身有自身的性质，都可以通知Eigen，让Eigen用对应的加速方式。比如正定矩阵可以只用上三角进行计算，并且在求解时使用Eigen:：LLT这样又快又数值稳定的解法，详情见Linear algebra and decompositions。

总之Eigen的高效既来自于算法优化技巧，又来自于工程师本身对公式和内存计算的理解。

参考资料: Eigen源代码+:Files·master:libeigen/eigen

Eigen官方数程:The Matrix class

作者 朱小霖

感兴趣的朋友可以看一下这个pdf，是eigen的主要开发者对其功能和优化进行的总结：

https://download.tuxfamily.org/eigen/eigen_CGLibs_Giugno_Pisa_2013.pdf

这里也根据这个pdf简单说两句。

Eigen的优化主要可以看这个结构图： 从上到下，expression tree是通过c++ template的方法把计算提取成AST一样的东西，并惰性求值，从而减少了对中间变量的需求，这种优化被称为expression template。具体来说，就是对于这样一个计算

A = B + C - D;

按正常的C++，展开直接计算应该类似：

// tmp1 = B + C
Matrix tmp1(n);
for (int i=0; i  tmp1[i] = B[i] + C[i];
// tmp2 = tmp1 - D
Matrix tmp2(n);
for (int i=0; i  tmp2[i] = tmp1[i] - D[i];
// A = tmp2
for (int i=0; i  A[i] = tmp2[i];

而通过抽象成AST，把A和B, C, D通过树连接起来之后再对A的每个元素求值，就可以达到这样的效果：

for (int i=0; i  A[i] = B[i] + C[i] - D[i];

可以很明显对比出两者之间的区别。这里需要注意的是，他的AST用的是C template做的，可能因为Eigen这个库是只有头文件的。以及为了防止重复实现，使用了CRTP这种设计模式。不过这种优化仅限于element-wise的操作，矩阵乘法就不行了。所以eigen对于乘法做了2件事，第一是把常见的乘法操作提取出来，提取为： 单独进行优化。然后因为矩阵乘法需要得到中间变量，所以为了更好的利用expression template，eigen用tree optimizer把矩阵乘法运算往后移。tree optimizer主要就是这样的2条运算顺序调整：

// catch A * B + Y and builds Y' + A' * B'
TreeOpt,Y> > { … };
// catch X + A * B + Y and builds (X' + Y') + (A' * B')
TreeOpt >,Y> > { … };

这样就可以把运算

res -= m1 + m2 + m3*m4 + 2*m5 + m6*m7;

调整为

res -= m1 + m2 + 2*m5;
res -= m3*m4;
res -= m6*m7;

注意这其中的每一行可以理解为一个中间变量。

做完tree optimizer，就要开始利用硬件信息优化了，比如做向量化来利用SIMD，进行loop unrolling什么的。对于是否要进行这些优化，eigen对于运算记录了一个经验性的cost model（这里我说经验性的是因为我没看到他利用上了矩阵的形状，所以判断其可能认为所有加法cost都一样，如果有误请大家及时纠正），然后设置了阈值来判断优化是否有益。

还是希望有兴趣的同学去看一下pdf，应该比我这里说的详细很多。以上。

作者 积石息土

看了这里的答案，我也不知道是对还是错，但是我觉得很关键的一点没有人说，那就是向量化 （VECTORIZE）或者SIMD优化，而这正是eigen比你手写的代码快的关键因素。eigen有一个很关键的开关叫做 EIGEN_DONT_VECTORIZE ，默认情况下是关掉的，也就是说eigen在默认情况下是会选择使用simd优化的。我们来做一个测试，对比一下eigen和C语言基础库Morn里面的矩阵运算。Morn是我手写的，可以看做是一个普通程序员手写矩阵运算的正常水平。先看最常见的矩阵乘法+。测试代码如下：

# include "Eigen/Dense"
# include "morn_math.h"

# define T 1000

void test_mul(int K) 
{
  printf("matrix mul, size: %d*%d\n", к,к);
  Eigen:: Matrixxf m_a(K,K); 
  Eigen:: MatrixXf m_b(K,K);
  MMatrix * mat_a = mMatrixCreate(K,K); 
  MMatrix * mat_b = mMatrixCreate(K,к);
  
  for(int j=0;jfor(int i=0;i  {
    float v; 
    V=mRand(-1000,1000)/1000.8;
    m_a(j,i)=v; 
    mat_a->data[ jl[i]=V; 
    
    V=mRand(-1000,1000)/1000.0;
    m_b(j,i)=v; 
    mat_b->data[ jl[i]=v; 
    }
    Eigen:: MatrixXf m_mul; 
    mTimerBegin("eigen");
    for(int i=0;i      m_mul = m_a*m_b;
      mTimerEnd("eigen");
      
  MMatrix *mat_mul = mMatrixCreate();
  mTimerBegin("Morn");
  for(int i=0;i      mMatrixMul(mat_a,mat_b,mat_ mul);
  mTimerEnd("Morn");
  mMatrixRelease(mat_mul);
  
  mMlatrixRelease(mat_a);
  mMatrixRelease(mat_b);
}
int main()
{
  test_mul(18);
  test_mul(20);
  test_mul(50);
  test_mul(100);
  test_mul(200);
  test_mul(500);
  return 0;
}

以上分别测试了从10×10矩阵到500×500矩阵的乘法。测试时计时1000次。测试使用的编译器如下： 先看，当定义了EIGEN_DONT_VECTORIZE的结果。编译命令为： 测试结果为： 可以看到当EIGEN_DONT_VECTORIZE时，eigen并没有比Morn更快，甚至还比Morn稍慢一点。再测试一下默认的打开向量化时的结果，测试命令为：

测试结果为：

可以看到当使用eigen向量化运算的时候，eigen就明显快于Morn了，尤其对于大规模矩阵运算，有2到4倍的差距。

以上测试足以说明VECTORIZE对Eigen速度的提升效果。那么VECTORIZE到底是什么东西呢。在 eigen的"src/Core/util/ConfigureVectorization.h"文件里。

#if !(defined(EIGEN _DONT_VECTORIZE) Il defined(EIGEN_GPUCC))
  #if defined (EIGEN _SSE2_ON_NON_MSVC _BUT _NOT_OLD_GCC) || defined(EIGEN _SSE2_ON_MSVC_2008_OR_LATER)
    //Defines symbols for compile-time detection of which instructions are
    //used. 
    //EIGEN VECTORIZE YY is defined if and only if the instruction set y is used
    # define EIGEN VECTORIZE
    # define EIGEN VECTORIZE SSE
    # define EIGEN VECTORIZE SSE2

    //Detect sse3/ssse3/sse4:
    //gcc and icc defines _SSE3_, ..
    //there is no way to know about this on msvc. You can define EIGEN VECTORIZE SSE* if you
//want to force the use of those instructions with msvc.
  #ifdef _SSE3_
    #define EIGEN VECTORIZE_SSE3
  #endif
  #ifdef _SSSE3_
    #define EIGEN VECTORIZE_SSSE3
  #endif
  #ifdef _SSE4_1
    #define EIGEN_VECTORIZE_SSE4_1
  #endif
  #ifdef SSE4_2_
    #define EIGEN_VECTORIZE_SSE4_2
  #endif
  #ifdef _AVX_
    #ifndef EIGEN _USE_SYCL
      #define EIGEN VECTORIZE AVX
    #endif
    #define EIGEN_VECTORIZE_SSE3
    #define EIGEN_VECTORIZE_SSSE3
    #define EIGEN_VECTORIZE_SSE4_1
    #define EIGEN_VECTORIZE_SSE4_2




    

    #endif
    ...
  #endif
  ...
#endif

可以看到eigen的向量化就是使用sse、avx等等simd指令，而且eigen做的很好，已经覆盖了今天所有主流的simd指令集（以上摘抄的只是一小部分）。

但是对于一些不适用于SIMD优化的矩阵运算，eigen就没有这么快了，诸如矩阵转置，它打不打开EIGEN_DONT_VECTORIZE其实影响不大： 其它更多的一些测试：

矩阵求逆：

当eigen不使用向量化的时候：

当eigen使用向量化的时候： 计算行列式：

当eigen不使用向量化的时候：

当eigen使用向量化的时候： 求解线性方程组：

这里使用了eigen的QR分解和LU分解两种方法，Morn里的实现实际为LU分解。

当eigen不使用向量化的时候：

当eigen使用向量化的时候：

所以，结论是eigen实际上整合了各种平台下的各种优化手段（其实不限于SIMD优化，还有诸如GPU加速之类（但是一般人用的时候都默认不打开GPU加速）），做了很多繁杂的工作，因此才达到了优秀的性能表现。

作者 叶洪舟

以前也觉得 Eigen 很快，感觉和 mkl 可以媲美。但是昨天 debug 时一行一行看时间，发现这样一个矩阵乘法 要算将近一秒 ,而同样的计算用 MATLAB 只需要肉眼不可察觉的时间。

我 google 了一下，在 StackOverflow 上

(How to speed up Eigen library's matrix product?)

有人讨论说是 MATLAB 内部自动会调用多线程版的 mkl 里的矩阵乘法，而 Eigen 在通常状态下是单线程的，需要在编译时加上-fopen参数使用 openmp 开启多线程。我在自己的 Mac 上尝试了一下，在线程数为 4 的情况下，时间缩短为 0.57 秒，但相比 MATLAB 的肉眼不可察觉还是有差距。

我对计算机底层的硬件不是很熟，只是最近突然觉得处理起这种可以高度矢量化的问题时，(不是 非常 wisely 地) 使用Eigen 还是没有经过优化的 MATLAB 快。

P.S.: 最开始的单线程版本我使用的是 OS X 自带的 g++ 编译器，优化参数为 -Ofast.测试过加上 -msse2 没有明显区别。这个g++ 实际上是Apple 封装后的 clang++,很可惜的是不支持 openmp.所以为了开启多线程我安装了GNU g++ 编译器。其它编译参数不变。

更新

刚才又上 StackOverflow 上看了别人的几个帖子 (

Eigen vs Matlab: parallelized Matrix-Multiplication

)学习了一下，发现了几个新的 trick:

1.首先之前线程数我设置错误了。我的 CPU physically 应该只有两个核，四线程+是 hyper- thread 的结果，多出来的两个线程并不能加速。我把 OMP_NUM_THREADS 改为 2 以后，速度 变为 左右；2.如果计算矩阵乘法，使用 noalias() 是个不错的 trick,它可以避免生成 temporary 的矩阵存储中间结果。使用这个后缩短为 0.50s 左右；3.对 Eigen 3.3 或以上的版本可以加上 -mavx 和 -mfma 两个参数，进一步缩短为 0.40s 左右。这样总的下来提速了约 30%.

补充

经评论区提醒 Ofast 进行了一些 unsafe+ 的数学优化所以可能牺牲精度，在极端情况下甚至出错。参考如下 stackoverflow 链接的讨论：

stackoverflow.com/quest...

作者 丁一帆

eigen最好的地方不应该是head-only 可以直接扔到cuda里面骚操作嘛？

其他啥matlab做得到嘛-几乎天痛切换到cuda的诱惑比写simd啊avx啊高的不知道哪里去了

(然而windows下eigen 在cuda下面各种报错，我已经看腻了bug fix了)

作者 darren xu

见了好多答主都不开编译优化-_-|||。

VS下选Release，g++下打O3优化。

另外，多测几次取平均，开个10000的for循环，当然写c++时要有些技巧以免for循环也被编译优化 掉。

x86平台上该打的sse或avx优化要打上，arm平台上neon优化要打上。调用矢量化指令的前提是保 证数据内存对齐 (128bit或256bit对齐),否则也不会生效甚至挂掉。

英特尔cpu上matlab的一些矩阵操作大多都调用release版本的mkl库，Eigen也能调用。

至于openmp多线程优化，好像VS自带的1.0版本就是屎，一些操作可能数据同步就会占去大把时 间。在arm上，由于获取释放锁的开销比x86大得多，用了openmp可能会负优化。

作者 Carpathia

感觉Eigen相比于其他线性代数库最突出的特点是，“不是效率高的事情不做，不是最快的算法就不 提供”。

看下Eigen自己说的，

Eigen: What happens inside Eigen, on a simple example

总结来说Eigen做了这么几件事：

1.不用中间变量

2.小矩阵和大矩阵根据实际情况选用

3.特殊的赋值处理

4.矢量化SSE指令

市面上Eigen库最快，对于Armadillo+数学库的组合速度取决于数学库的速度，Armadillo本身数学库，就不用提了。

作者 生死看淡不服就干

Eigen库本身是一个模板类的数值运算库，在做矩阵运算+的时候一般能够取得比较好的结果.但是在 很多时候还有内存相关问题，内存分配的可以通过两种方式来解决：

· 减少动态内存的分配 · 自定义内存池

第一种方式是可以通过注意代码的实现方式，以及使用开发过程中的trick来解决，对原生库没有任何侵入性，对开发者需要相对深入了解Eigen库+本身.本文从这方面入手，详细介绍一些有效的优化手段

库

·在generalmatrixProcuct的具体实现中

使用了在栈上申请内存的宏定义函数*ei_declare_aligned_stack_constructed_variable.简化一下具体实现如下

#ifdef EIGEN ALLOCA