# Difference between revisions of "FAQ"

m (fix link) |
(→Is there a method to compute the (Moore-Penrose) pseudo inverse ?) |
||

Line 241: | Line 241: | ||

== Is there a method to compute the [http://en.wikipedia.org/wiki/Pseudo_inverse (Moore-Penrose) pseudo inverse] ? == | == Is there a method to compute the [http://en.wikipedia.org/wiki/Pseudo_inverse (Moore-Penrose) pseudo inverse] ? == | ||

− | + | Yes, using a complete orthogonal decomposition: | |

− | + | ||

+ | MatrixXd A = ...; | ||

+ | MatrixXd Ainv = A.completeOrthogonalDecomposition().pseudoInverse(); | ||

− | + | If you really want a SVD-based solution, then you can copy paste the following c++11 code snippet: | |

+ | |||

+ | template<typename MatType> | ||

+ | using PseudoInverseType = Eigen::Matrix<typename MatType::Scalar, MatType::ColsAtCompileTime, MatType::RowsAtCompileTime>; | ||

+ | |||

+ | template<typename MatType> | ||

+ | PseudoInverseType<MatType> pseudoInverse(const MatType &a, double epsilon = std::numeric_limits<double>::epsilon()) | ||

{ | { | ||

− | + | using WorkingMatType = Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, Eigen::Dynamic, 0, MatType::MaxRowsAtCompileTime, MatType::MaxColsAtCompileTime>; | |

− | + | Eigen::BDCSVD<WorkingMatType> svd(a, Eigen::ComputeThinU | Eigen::ComputeThinV); | |

− | + | svd.setThreshold(epsilon*std::max(a.cols(), a.rows())); | |

− | + | Eigen::Index rank = svd.rank(); | |

− | + | Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, MatType::RowsAtCompileTime, | |

− | + | 0, Eigen::BDCSVD<WorkingMatType>::MaxDiagSizeAtCompileTime, MatType::MaxRowsAtCompileTime> | |

− | + | tmp = svd.matrixU().leftCols(rank).adjoint(); | |

− | + | tmp = svd.singularValues().head(rank).asDiagonal().inverse() * tmp; | |

− | + | return svd.matrixV().leftCols(rank) * tmp; | |

− | + | } | |

− | See [http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 Bug 257] if there is any progress in natively providing this. | + | See [http://eigen.tuxfamily.org/bz/show_bug.cgi?id=257 Bug 257] if there is any progress in natively providing this SVD-based solution. |

There's no pseudo inverse for sparse matrices (yet?) in Eigen. | There's no pseudo inverse for sparse matrices (yet?) in Eigen. |

## Revision as of 15:25, 13 November 2019

## Contents

- 1 Licensing
- 2 Eigen and other libraries
- 3 Eigen and IDEs
- 4 Compilation
- 5 Runtime
- 6 Optimization
- 7 Vectorization
- 7.1 Which SIMD instruction sets are supported by Eigen?
- 7.2 How can I enable vectorization?
- 7.3 How can I disable vectorization?
- 7.4 I disabled vectorization, but I'm still getting annoyed about alignment issues!
- 7.5 How does vectorization depend on the compiler?
- 7.6 What can't be vectorized?
- 7.7 How can I check that vectorization is actually being used?

- 8 Design
- 9 Algorithms
- 10 Other languages

# Licensing

Since the 3.1.1 release, Eigen is licensed under the MPL2. We refer to the MPL2 FAQ for initial questions.

## Disabling non MPL2 features

Even though most of the Eigen code is under the MPL2, a few features are still under the LGPL license. This concerns: SimplicialCholesky, AMD ordering and constrained_cg. Non MPL2 features can be explicitly disabled by compiling with the EIGEN_MPL2_ONLY preprocessor token defined.

# Eigen and other libraries

## Should I use Eigen?

Probably, but check pit falls first.

## Why another matrix library? What is the need for Eigen?

First of all, see the Overview. No other library provides all of the features and benefits listed there.

The Eigen project started when some hackers from the large KDE meta-project realized the need for a single unified matrix library.

Some other libraries do satisfy very well certain specialized needs, but none is as versatile as Eigen, has such a nice API, etc.

The fact that so many projects are quickly adopting Eigen 2, shows that it fills a gap.

The state of existing matrix libraries before Eigen is that:

- some are Free Software
- some are fast
- some have a decent API
- some handle fixed-size matrices, some handle dynamic-size dense matrices, some handle sparse matrices
- some provide linear algebra algorithms (LU, QR, ...)
- some provide a geometry framework (quaternions, rotations...)

However Eigen is the first library to satisfy all these criteria.

## How does Eigen compare to BLAS/LAPACK?

Eigen covers many things that BLAS/LAPACK don't:

- Eigen handles fixed-size matrices and vectors, which are very widely used.
- Eigen has built-in support for sparse matrices and vectors.
- Eigen provides a lot of convenience features (see Geometry module, Array module, etc), which are also very widely used.

Using only one thread, Eigen compares very well performance-wise against the existing BLAS implementations. See the benchmark. It shows that:

- Eigen is faster than every Free BLAS, such as ATLAS or Boost::uBlas.
- Eigen is overall of comparable speed (faster or slower depending on what you do) to the best BLAS, namely Intel MKL and GOTO, both of which are non-Free.

However, currently Eigen parallelizes only general matrix-matrix products (bench), so it doesn't by itself take much advantage of parallel hardware.

Eigen has an incomparably better API than BLAS and LAPACK.

- See the API Showcase.
- For operations involving complex expressions, Eigen is inherently faster than any BLAS implementation because it can handle and optimize a whole operation globally -- while BLAS forces the programmer to split complex operations into small steps that match the BLAS fixed-function API, which incurs inefficiency due to introduction of temporaries. See for instance the benchmark result of a Y = a*X + b*Y operation which involves two calls to BLAS level1 routines while Eigen automatically generates a single vectorized loop.

Miscellaneous advantages (not specifically against BLAS/LAPACK):

- Eigen is only a compile-time dependency for your project. No need to redistribute, or ask your user to install, any library.
- Eigen is small, so it is feasible to include a copy of it in your own source tree, if you want to.
- Eigen is multi-platform, and is actually being used on a number of different operating systems, hardware platforms, and compilers.
- Eigen, compared to certain other C++ template libraries, is relatively easy on the compiler. Compilation times stay reasonable -- we are very careful about that.

# Eigen and IDEs

For details on how to correctly use Eigen with IDEs (such as Eclipse CDT), see IDEs.

# Compilation

## I need help with compiler errors!

- Did you forget to include a header? See Pit Falls.
- Did you check if you triggered a
*static assertion*? These are compile-time checks guarding from programming mistakes. Eigen has many of those. So even if you got a kilometer of compiler output, you might still find useful information from static assertion messages. Search for "static_assert". The static assertion messages themselves are UPPERCASE_SO_THEY_REALLY_STAND_OUT. - Did you forget the
`template`or`typename`keyword? This is a fairly obscure part of the C++ language, that programmers typically stumble upon when they are writing a template function which takes Eigen objects. See the page The template and typename keywords in C++ for some explanation. A typical error message from the GCC compiler is "expected primary-expression before ')' token". - Eigen is a C++ library and therefore there might be some incompatibility with some C99 or intrusive C header files, e.g.: complex.h, X.h, etc. In most cases, including Eigen first is enough to workaround the issue.

## Known MSVC issues

- MSVC 2010 sometime crashes when the "enable browse" compiler option (/FR) is activated.
- It has been reported that MSVC 2008 pro-sp1 generates wrong code for
`Matrix3d M; Matrix3d I = M * M.inverse();`

in 32 bits and Debug mode. and Eigen 3.2 This is specific to`Matrix3d`

. The workaround is is to call`.eval()`

:`M.inverse().eval()`

. (see bug 1442) - It has been reported that MSVC 2015 can be confused by the simple declaration of an Array of void* (e.g.,
`Eigen::Array< void*, 64, 1 > handles;`

), reporting errors in unrelated code (see bug 1440). Since Eigen's Array is not designed to hold pointers, the proper workaround is to use c++11 std::array.

# Runtime

## I need help with Assert crashes!

The asserts are there to protect you from later unexplained crashes due to bad memory accesses.

When you hit such an assert, rerun your program in a debugger and obtain a backtrace. Make sure that you have compiled your program with enough debugging info. This way, you will quickly be able to trace back to the root cause of the problem :)

The most dreaded assert is the "Unaligned array" assert. As you can see, that page is there to help you fix it. If however you are desperate about it, you can always get rid of it.

Other assertions are typically triggered when you have accessed coefficients with out-of-range indices in a matrix; or when you have mixed matrices of mismatched sizes.

# Optimization

## How do I get good performance?

There are many aspects to this question, but here are some points to get you started:

- Make sure you compile with optimization enabled. This can easily gain you a factor of ten or more.
- Enable vectorization, as described in How can I enable vectorization?.
- If your matrices are very small (size between 2 and 4), then using fixed-size matrices instead of dynamic-size matrices can get you a substantial speed-up.
- Sometimes, quite a lot of time can be spent by the creation of temporary objects to hold intermediate results. See the next question to track this down.
- Matrix multiplications are costly. See Writing Efficient Matrix Product Expressions for some advice.
- Define the NDEBUG macro when compiling; this disables some run-time checks which slightly speeds up your program. However, this should not be done prematurely, as it also hides otherwise hard-to-find programming mistakes (such as combining expressions of incompatible dimensions)!
- General programming advice also applies here. In particular, profile your code, find the bottleneck, and optimize that part.

## Where in my program are temporary objects created?

The Eigen library sometimes creates temporary matrices to hold intermediate results. This usually happens silently and may slow down your program, so it is useful to track down where temporary objects are created.

One possibility is to run your program under a debugger and set a break point which will be triggered when a temporary is created. For instance, you can set a break point in check_that_malloc_is_allowed() in Eigen/src/Core/util/Memory.h (this function is probably inlined if you compile with optimizations enabled).

Another possibility is to define the macro EIGEN_NO_MALLOC when compiling. This causes your program to abort whenever a temporary is created. More fine-grained checks are possible with the EIGEN_RUNTIME_NO_MALLOC macro. A minimal usage example follows:

#define EIGEN_RUNTIME_NO_MALLOC // Define this symbol to enable runtime tests for allocations #include <Eigen/Dense> int main(int argc, char** argv) { // It's OK to allocate here Eigen::MatrixXd A = Eigen::MatrixXd::Random(20, 20); Eigen::internal::set_is_malloc_allowed(false); // It's NOT OK to allocate here // An assertion will be triggered if an Eigen-related heap allocation takes place Eigen::internal::set_is_malloc_allowed(true); // It's OK to allocate again }

# Vectorization

## Which SIMD instruction sets are supported by Eigen?

Eigen supports SSE, AVX, AVX512, AltiVec/VSX (On Power7/8 systems in both little and big-endian mode), ARM NEON for 32 and 64-bit ARM SoCs, and now S390x SIMD (ZVector).

With SSE, at least SSE2 is required. SSE3, SSSE3 and SSE4 are optional, and will automatically be used if they are enabled.

Of course vectorization is not mandatory -- you can use Eigen on any CPU.

Note: For S390x SIMD, due to lack of hardware support for 32-bit vector float types, only 32-bit ints and 64-bit double support has been added.

## How can I enable vectorization?

You just need to tell your compiler to enable the corresponding instruction set, and Eigen will then detect it. If it is enabled by default, then you don't need to do anything. On GCC and clang you can simply pass -march=native to let the compiler enables all instruction set that are supported by your CPU.

On the x86 architecture, SSE is not enabled by default by most compilers. You need to enable SSE2 (or newer) manually. For example, with GCC, you would pass the -msse2 command-line option.

On the x86-64 architecture, SSE2 is generally enabled by default, but you can enable AVX and FMA for better performance

On PowerPC, you have to use the following flags: -maltivec -mabi=altivec, for AltiVec, or -mvsx for VSX-capable systems.

On 32-bit ARM NEON, the following: -mfpu=neon -mfloat-abi=softfp|hard, depending if you are on a softfp/hardfp system. Most current distributions are using a hard floating-point ABI, so go for the latter, or just leave the default and just pass -mfpu=neon.

On 64-bit ARM, SIMD is enabled by default, you don't have to do anything extra.

On S390X SIMD (ZVector), you have to use a recent gcc (version >5.2.1) compiler, and add the following flags: -march=z13 -mzvector.

## How can I disable vectorization?

You can disable Eigen's vectorization by defining the EIGEN_DONT_VECTORIZE preprocessor symbol.

If you also want to disable the "unaligned array" assertion or the 128bit alignment code, see the next entry below.

Also notice that your compiler may still be auto-vectorizing.

## I disabled vectorization, but I'm still getting annoyed about alignment issues!

For example, you're still getting the "unaligned array" assertion.

If you want to get rid of it, follow this link.

## How does vectorization depend on the compiler?

Eigen has its own vectorization system, it does not at all rely on the compiler to automatically vectorize. However it still needs some support from the compiler, in the form of *intrinsic functions* representing a single SIMD instruction each.

Eigen will automatically enable its vectorization if a supported SIMD instruction set and a supported compiler are detected. Otherwise, Eigen will automatically disable its vectorization and go on.

Eigen vectorization supports the following compilers:

- GCC 4.2 and newer,
- MSVC 2008 and newer,
- All other compilers (for example it works with clang and ICC).

Of course the reason why we "support all other compilers" is that so far we haven't seen other examples of compilers on which we should disable Eigen vectorization. If you know some, please let us know.

## What can't be vectorized?

SSE, AltiVec, NEON and S390x work with packets of 128 bits, or 16 bytes. This means 4 ints, or 4 floats, or 2 doubles, except for S390x which doesn't support 32-bit floats in hardware. Even though it is usually preferable to work on packets that are 128-bit aligned, starting from Eigen 3.3, Eigen will also vectorize unaligned vectors and matrices at a last resort. So only tiny vectors of matrices that are smaller than a packet won't be vectorized. This is for instance the case of Vector3f which is smaller than a packet. However, you may be able to use Vector4f class to perform Vector3f operations and make use of vectorization if you can carefully ensure that the last component is always zero. . Moreover, Matrix2f internally stores 4 floats, and vectorization will automatically take place if the expression can be seen as a 1D vector of 4 elements (e.g., adding two Matrix2f).

The vectorization of unaligned matrices or vectors can be disabled by defining EIGEN_UNALIGNED_VECTORIZE=0. In this case fixed-size vectors and matrices whose size in bytes is are not a multiple of 16 bytes won't be vectorized. Vectors and matrices of dynamic-sizes will still be vectorized by enforcing aligned stores at runtime.

The same remarks apply for AVX which operates on packets of 32 bytes for which 32 bytes alignment is required. Since AVX can also operate with 16 bytes packets, the optimal packet size and alignment is automatically figured-out for fixed-size vectors and matrices. For instance, 16bytes alignement and packets will be used for Vector2d, Vector4f, Vector6d, etc. For dynamic-sizes vectors and matrices, only 32 bytes packets will be used, meaning that VectorXf(4) won't be vectorized with AVX whereas it will be vectorized with SSE only.

Eigen also has some current limitations that can and will be overcome in the future. For instance, some advanced operations, such as visitors, aren't currently vectorized (this is on our to-do).

## How can I check that vectorization is actually being used?

First you can check that Eigen vectorization is enabled: the EIGEN_VECTORIZE preprocessor symbol is then defined.

Then, you may want to check the resulting assembly code. This is the best way to check that vectorization actually happened for a specific operation. Add some asm comments in your code around a line of code you're interested in, and tell your compiler to output assembly code. With GCC you could do:

Vector4f a, b; asm("#it begins here!") a += b; asm("#it ends here!")

See Studying assembly output for more information.

# Design

## Why Eigen's API is using signed integers for sizes, indices, etc.?

At a first glance it is indeed very tempting to use an unsigned type for values that are expected to be positive or zero. However, unsigned types relies on modulo 2^k arithmetic, which is definitely not what we want here. Moreover, in many places we can have negative indices and offsets. Mixing signed and unsigned integers is costly, error prone, and highly correlated to bugs. As a result, it is much simpler and safer to use signed integers everywhere, and reserve the usage of unsigned integers for bitwise manipulation and modulo 2^k arithmetic.

Hm... but the STL use unsigned integers everywhere. Right, but this has been acknowledge as a mistake by several members of the C++ standardization committee and C++ experts (Bjarne Stroustrup, Herb Sutter, Chandler Carruth, Eric Niebler, etc.). Sources:

- http://channel9.msdn.com/Events/GoingNative/2013/Interactive-Panel-Ask-Us-Anything (at 12mn, 42,mn and 63mn)
- http://ericniebler.com/2014/12/07/a-slice-of-python-in-c/#comment-92581

# Algorithms

## Is there a method to compute the (Moore-Penrose) pseudo inverse ?

Yes, using a complete orthogonal decomposition:

MatrixXd A = ...; MatrixXd Ainv = A.completeOrthogonalDecomposition().pseudoInverse();

If you really want a SVD-based solution, then you can copy paste the following c++11 code snippet:

template<typename MatType> using PseudoInverseType = Eigen::Matrix<typename MatType::Scalar, MatType::ColsAtCompileTime, MatType::RowsAtCompileTime>; template<typename MatType> PseudoInverseType<MatType> pseudoInverse(const MatType &a, double epsilon = std::numeric_limits<double>::epsilon()) { using WorkingMatType = Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, Eigen::Dynamic, 0, MatType::MaxRowsAtCompileTime, MatType::MaxColsAtCompileTime>; Eigen::BDCSVD<WorkingMatType> svd(a, Eigen::ComputeThinU | Eigen::ComputeThinV); svd.setThreshold(epsilon*std::max(a.cols(), a.rows())); Eigen::Index rank = svd.rank(); Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, MatType::RowsAtCompileTime, 0, Eigen::BDCSVD<WorkingMatType>::MaxDiagSizeAtCompileTime, MatType::MaxRowsAtCompileTime> tmp = svd.matrixU().leftCols(rank).adjoint(); tmp = svd.singularValues().head(rank).asDiagonal().inverse() * tmp; return svd.matrixV().leftCols(rank) * tmp; }

See Bug 257 if there is any progress in natively providing this SVD-based solution.

There's no pseudo inverse for sparse matrices (yet?) in Eigen.

## How to implement IIR filters?

BTK propose an experimental module for Eigen implementing Infinite Impulse Response filters for 1D signals there.

# Other languages

## Can I use Eigen with Python?

- sparray: a python binding of the Sparse module - alpha stage.
- pyeigen: bindings for MatrixXd.
- minieigen: bindings for several matrix and vector classes (Vector{2,3,6,X}, Vector{2,3,6,X}i, Matrix{3,6,X}, AlignedBox{2,3}, Quaternion); debian packages
- BTK proposes a fragment for SWIG to transform an Eigen matrix to a NumPy array.

## Can I use Eigen with R?

- The RcppEigen package provides bindings and more for R.

## Can I use Eigen with Java?

- jeigen, a Java wrapper for Eigen.