Eigen的速度为什么这么快？

2015-03-06财经

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

比如一个矩阵乘法

#include <Eigen/Dense> 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 < double , Dynamic , Dynamic , ColMajor >

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

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

静态内存分配不但让我们节省了new/delete的开销，还给Eigen内部继续优化提供了可能。Eigen内置Single-Instruction-Multiple-Data( SIMD )指令集，对稠密矩阵有很好的优化，如果能触发CPU SIMD的指令，能收获成倍的计算效率。

最后是内存的读写。Lazy evaluation的另一个特点是不生成中间变量，减少内存搬运次数，而Eigen为了防止矩阵覆盖自己，对矩阵-矩阵乘法会生成一个中间变量。如果我们知道等式左右两边没有相同的项，则可以通知Eigen去取消中间变量。

Hessian . noalias () += Jacobian_i . transpose () * Jacobian_j ;

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

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

参考资料：

Eigen源代码：Files · master · libeigen / eigen

Eigen官方教程：The Matrix class