Interview: Eigen Matrix Library

来源:互联网 发布:太极禅淘宝 编辑:程序博客网 时间:2024/06/03 19:56

http://www.macresearch.org/interview-eigen-matrix-library

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:

  • Computer graphics projects:
    • MeshLab
    • VcgLib
    • libmv
    • Krita
  • Robotics:
    • Yujin Robot
    • Willow Garage
  • KDE and KOffice are, of course historical users. More specifically:
    • the above-mentioned Krita
    • a little bit in KSpread
    • some screensavers, KGLlib, KGLEngine2d, SolidKreator. 
    • Kalzium's molecular viewer is using Eigen indirectly through Avogadro (see below).
    • In KDE 4.3, we have good hopes to add Step to this list.
  • Chemistry projects:
    • Avogadro
    • OpenBabel (in version 3)

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 1000x1000 matrix A and a 1000x500 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 matricesusing 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.

Comments

Comment viewing options

 
 
 
Select your preferred way to display the comments and click "Save settings" to activate your changes.

Benoit's Talk

Benoit also gave a recent talk about Eigen, which gives a few more code examples. His slides can be found here.

Sounds great!

Excellent. Over the years I have often searched for good and fast matrix / linear algebra C++ library, always hoping that something better had come along. But it hadn't. I ended up using ublas, which was okay, but not as fun to use or fast as it could have been. And (at least at the time) the Lapack bindings sort of existed, but were hidden away in some attachment of a Yahoo group and you had to integrate them yourself into the latest version if you really really wanted them.

Eigen sounds exactly like what I was hoping for. And the performance in those benchmarks seems pretty fantastic. Great job! I will definitely use this for my next C++ project.

Nice Interview and library!

Nice interview, Geoff, and it looks like a good library. I am particularly glad they didn't make the mistake of going full GPL. That would severely limit the project.

I spent several months writing a C++ program a few years back. I looked at all the options -- blitz++, pooma, ublas -- but none of them really were satisfactory, and I actually ended up rolling my own for the subset of linear algebra that I needed. If Eigen is to succeed, it needs to be maintained. That was one problem with very promising libraries like blitz++: there was no further development, so adopting them was a big risk.

I really hope Eigen succeeds.

Drew

PS A word of warning about C++ templates, and expression templates in particular: it all sounds great on paper, but in practice they can be tricky for a beginner. Often you will get 3 pages of error messages for a single error in your code. It can take some getting used to.

---------------------------
Drew McCormack
http://www.maccoremac.com
http://www.macanics.net
http://www.macresearch.org

Development & Error Messages

I think the Eigen development is very strong right now. Obviously the test with open source projects is if they can last 3-4 years or more. But the development mailing list is filled with newcomers and contributions, so it's off to a very good start.

As far as templates and error messages? Eigen really knows your concerns. Here's an example error message:

/Users/ghutchis/Devel/avogadro/libavogadro/src/extensions/orbitaldialog.cpp:282: instantiated from here
/usr/local/include/eigen2/Eigen/src/Core/Assign.h:404: error: ‘YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY’ is not a member of ‘Eigen::ei_static_assert’

So yes, I have about a page of error messages, but GCC spits up a helpful one at the end -- the line number of my mistake, plus Eigen's complaint that I need to cast() to go from a matrix of double to a matrix of float.

Eigen2 and Fink

FYI, eigen (version 2) is now available through Fink.


0 0