Interview: Eigen Matrix Library

·

·

What follows is an interview with Benoît Jacob and Gael Guennebaud, two developers of the open source Eigen project, which provides incredibly fast C++ matrix and vector numeric code.

OK, let’s start out with the basics. Could you introduce yourselves?

Benoit: I’m a Mathematics postdoc at the University of Toronto, coming from France. I’m working on C*-algebras, and teaching linear algebra. In addition to Eigen, I contribute to a few other free software projects, mostly within KDE. In the past I also contributed to Avogadro and a little bit to Open Babel.

Gael: I’m a French researcher in Computer Graphics currently working at INRIA Bordeaux. In particular, my research interests include real-time rendering and surface representations. My major contribution in the open-source world is by far what I did in Eigen, but I also contributed a bit to vcglib and MeshLab.

What is Eigen?

Benoit and Gael: Eigen is a free C++ template math library mainly focused on vectors, matrices, and linear algebra. It is a self-contained library covering a very broad range of use cases. For example, it covers both dense and sparse objects, and in the dense case, it covers both fixed-size and dynamic-size objects. Moreover it provides linear algebra algorithms, a geometry framework, etc. It has a very nice API for C++ programmers, and it embraces very high performance.

What drove you to create Eigen and Eigen 2?

Benoit: Eigen 1 was a very small project, 2500 LOC and just a few months of development. Its creation in 2006 was driven by the then-simple needs of some KDE and KOffice apps. However it quickly turned out that we had underestimated their needs, and Eigen 1 was insufficient. So in 2007, I started developing Eigen2. The aim was to create the linear algebra library that would satisfy all the needs of KDE and KOffice apps — a goal that, in retrospect, was very ambitious, and will only be reached with Eigen 2.1. After an initial experiment with TVMET‘s code, I decided to restart from scratch in August 2007 and quickly got a working implementation of expression templates. However, this early Eigen 2 was very small. Development speed really picked up when Gael joined in early 2008.

Gael: A bit more than a year ago, I became tired of going back and forth between my own fixed size vector/matrix classes and more generic linear algebra packages. So, I started looking at the other existing solutions without being excited by any of them. Since I have been using KDE for about 9 years, I was really curious to know what the KDE’s folks did in this area. At that time, it was exactly the start of Eigen2 which looked promising but the fact it was based on TVMET puzzled me. Eventually, Benoît had the great idea to restart the development of Eigen2 from scratch, and after one or two months he came up with a very lean design. Moreover his vision and feature plan for Eigen2 exactly matched my own, and being part of the KDE community was exciting too! At the beginning, I naively thought that after one or two months the job would be done! Instead, we started playing with exciting stuff like explicit vectorization, efficient matrix products, sparse matrix, etc.

Many people are familiar with other linear algebra and matrix libraries, including BLAS/LAPACK, Intel’s Math Kernel Library, and Apple’s vecLib framework. Can you explain how Eigen is different, besides being written in C++?

Benoit and Gael:

Giving a fair answer to that question would require a thorough comparison to all existing libraries that is obviously out the scope of this interview. Search for “C++ matrix library” to get an idea. For us, most important criteria include:

  • generality,
  • performance,
  • ease of use,
  • license policy.

Even more important is a sort of meta criterion which would glue together all these criteria. Indeed, what’s the point in having a very efficient library but which which cannot handle your specific needs, and/or that you cannot freely use/extend because of the license ? The central difference between Eigen and other libraries reside in this unmatched combination of advantages.

For instance, BLAS/LAPACK are certainly the most popular matrix/linear algebra packages. They can achieve very high performance, but only for relatively large and dense matrices, and through proprietary implementations such as Intel’s MKL or Apple’s vecLib. Goto BLAS too has a restrictive license. API-wise, BLAS/LAPACK are also very painful to use, though that limitation is somewhat overcome by the existence of many high-level wrappers (FLENS, matlab, etc.). In contrast, not only does Eigen provide a feature set to BLAS/LAPACK with better ease of use, but it also seamlessly supports small fixed size matrices, has experimental support for sparse matrices with optional back-ends, and includes other modules such as linear least squares and geometry.

As we said, expression templates bring unmatched expressiveness and performance for basic operations. It is therefore interesting to have a brief look at the few other C++ template libraries which are also based on expression templates. The most popular ones include TVMET, Blitz++, and Boost::uBLAS. However, in contrast to Eigen, none of them provide explicit vectorization nor include linear algebra algorithms. They also perform poorly for both small and large matrix sizes, and suffer from some shortcomings in the way expression templates have been implemented making them hard to use in practice. For example, contrary to Eigen, they do not decide automatically when to use lazy evaluation. Eigen handles most cases automatically while ultimately letting the user override the default behavior. For example, in a chained matrix product A*B*C, Eigen understands that lazy evaluation would be a bad idea so decides to evaluate A*B into a temporary, while Boost::uBLAS requires the user to do that explicitly. Another advantage of Eigen over these other C++ libraries, is that Eigen has still reasonable compilation times: we are very careful about that.

How does performance compare to other alternatives? Do you have some benchmarks?

Gael: Since good performance is one of our top priorities, we indeed regularly run several benchmarks. When speaking about performance it is important to distinguish between small fixed size and larger matrix which require very different kind of optimizations.

About the first category, Eigen includes meta-unrollers and an advanced explicit vectorization mechanism allowing us to generate machine code as good as hand written assembly code. Therefore, in all our benchmarks, we found Eigen outperforms all other libraries, or at least we made sure that was the case by fixing Eigen when needed 😉

For relatively large matrices, we have a comparison of several free and non-free libraries. These results demonstrate that Eigen outperforms all other Free libraries, including ATLAS. Compared to non-Free alternatives like Intel’s MKL and Goto, Eigen is usually faster for both Level-1 (vector-vector) and Level-2 (matrix-vector) equivalent BLAS operations. This is mainly due to expression templates, and a highly optimized matrix-vector product algorithm. However, we acknowledge that matrix-matrix operations are not yet as optimized as their respective BLAS Level-3 routines of MKL and Goto.

Regarding higher-level linear algebra algorithms like LU or Cholesky factorizations, we currently do not provide any block implementations for them. Therefore the equivalent LAPACK routines of vendor tuned implementations scale better for very large matrices. Compared to them, we also lack a support for multi-core computing. But nothing is engraved in stone, and we are confident that all these performance limitations for very large matrices will be overcome in the future, in the short term by adding an optional LAPACK backend to Eigen, and in the longer term, who knows 😉

What projects are using Eigen right now?

Benoit and Gael: Eigen is already used in a wide range of applications. Some of them include:

It is very interesting and motivating to see how many projects already switched, or are going to switched to Eigen, proving Eigen fills a real gap.

How would you write code with Eigen? Could you give me an example, say of multiplying the inverse of matrix A by matrix B?

Benoit: This operation amounts to solving the system AX=B which can be done for example using a LU decomposition. With Eigen you would do:

A.lu().solve(B, &X);

Now, depending on your context, you may prefer using another kind of decomposition. Here, Eigen also allows you to use a QR, LLt, LDLt, or SVD decomposition, each time with the same obvious syntax, for example:

A.svd().solve(B, &X);

If your matrices have small and fixed size, and especially if B is a vector, it can be faster to actually compute the inverse and the product:

A.inverse() * B;

Here is a full example program that you can compile. It works with a 1000×1000 matrix A and a 1000×500 matrix B, with double precision. I get a relative error of the order of 1e-13 which is pretty good. Eigen’s LU does full pivoting which guarantees that it works in all cases, and an user just reported speed to be on par with LAPACK. However, for this kind of large sizes, depending on the CPU cache size, we acknowledge that it would start to be beneficial to have a more cache-friendly implementation, typically working block-wise.

#include <Eigen/LU>// provides LU decomposition
#include <Eigen/Array>// provides random matrices

using namespace Eigen;
int main()
{
  MatrixXd A = MatrixXd::Random(1000, 1000);
  MatrixXd B = MatrixXd::Random(1000, 500);
  MatrixXd X;

  A.lu().solve(B, &X);

  std::cout << "relative error: " << (A*X - B).norm() / B.norm() << std::endl;
}

More examples are given in the API showcase.

Since Eigen is a template library, can you explain the license? Can you use it in proprietary projects?

Benoit: The license is dual LGPL3+/GPL2+. This means it’s LGPL3+, except that GPL2 projects can’t use LGPL3+ libraries so just for them we added the GPL2+ as an alternative. Yes, proprietary closed-source projects may use Eigen, without having to open any of their source code. The LGPL license allows it and we do have several closed-source users.

There used to be a problem with the older LGPL 2.1 and template libraries. Indeed, the phrasing of the LGPL 2.1 was not applicable to the case of C++ template libraries where the code is in header files and is included rather than linked to. Fortunately this issue was fixed in the LGPL3 which now explicitly adresses this issue in Section 3, in addition to having a much more generally applicable phrasing overall. See our Licensing FAQ.

What are some goals of Eigen moving forward? What kind of help do you need? What are some new features we might see in Eigen 2.1?

Benoit and Gael: Our goals for 2.1 are:

Our goals for 2.1 are:

  • Finish stabilizing the API (the API guarantee is only partial in 2.0)
  • Complete the Sparse module. One goal is to make it good enough for Step in KDE 4.3 and Krita.
  • Make sure all dense decompositions are world-class (LU is already good though we have improvements in mind, SVD is being completely rewritten, etc…)
  • Make fixed-size specializations of algorithms
  • Optimize the case of complex numbers
  • Vectorize more operations (e.g., exp() or sin() or log())
  • Investigate optionally leveraging new GCC 4.4 hotness: per-function optimization flags

We need a lot of help with all that, more details can be found in our To-do. If new contributors join the team, in the longer term, we could see some new modules covering statistics, fast fourier transform, non linear optimizations, etc.

We also welcome testing, especially on exotic (read: non-GNU or non-x86/x86-64) platforms. All bug reports are much appreciated.


Leave a Reply

Your email address will not be published. Required fields are marked *