Eigen

Very useful library for Transformations. First interacted through F1TENTH, getting better foundation through Visual SLAM project.

Eigen is a special library built with pure header files (this is the amazing part!). This means you can only locate its header files, not binary files like .so or .a. When you use it, you only need to import Eigen’s header file. You don’t need to link the library file (because it doesn’t have any library files). Now let’s write a piece of code below to actually practice the use of Eigen

Because the Eigen library only has header files, we don’t need to link the program to the library with the target_link_libraries statement.

Resources

Fundamentally understanding how Eigen works

Installation

sudo apt install libeigen3-dev

CMake

In CMake:

find_package(Eigen3 REQUIRED)
include_directories(${EIGEN3_INCLUDE_DIR})

Notice that we don’t need add_library nor target_link_libraries. This is because Eigen only has header files (no .cpp files), so we don’t need to link the program. You can just run them in the include!!! ahh I see.

Eigen is ROW-major by default

This distinction is important to make once you get into the territory of storing large values in a particular dimension. Source

Types of Matrices in Eigen

  1. Fixed-size matrices: Dimensions determined at compile time. They allow for optimizations such as stack allocation and loop unrolling.
Matrix3d m;
 
// Under the hood
Matrix<double, 3,3> m;
  1. Dynamic-size matrices: Dimensions determined at runtime. These matrices offer flexibility, but with potential additional overhead due to heap allocation.
MatrixXd m(3,3);
 
// Under the hood
Matrix<double, Dynamic, Dynamic>;
m.resize(1,9);
  1. Hybrid matrices: One dimension is determined at compile time, the other at runtime. They offer a balance between the flexibility of dynamic-size matrices and the performance optimizations of fixed-size matrices.
Matrix<double, 3, Dynamic> m(3,3);

Are Eigen matrices heap or stack allocated?

Found the answer in this stack overflow thread. Seems like dynamic-size matrices are always heap allocated. Fixed-size matrices are oftentimes stack allocated, but sometimes also heap allocated.

Matrices value types

Matrix3

  • Matrix3d is a 3x3 of type double (the d stands for double)
  • Matrix3f is a 3x3 of type float (the f stands for float)
  • Matrix3i is a 3x3 of type integer (the i stands for integers)

You can apply the same logic for dynamic Matrices MatrixX

  • MatrixXd is of type double
  • MatrixXf is of type float (the f stands for float)
  • MatrixXi is of type integer

Vectors (they’re just 1d matrices)

Vector3d v;
 
// Under the hood
Matrix<double, 3, 1> v;
 
// Dynamic vector shape
VectorXd v(3);
 
// For vectors, you can actually directly set the values in the constructor
Vector3d v(1,0,0);

There also exists row vectors ( instead of )

  • RowVector3d
  • RowVectorXd

Matrix Initialization

What is the preferred way to initialize Eigen matrices? Depends on what you need. Don’t need to initialize matrix to 0 if you will overwrite it right after.

Rule of Thumb

As a general rule of thumb, if it is a dynamic matrix, you need the two extra arguments in front, which specify the dimensions of the matrix.

  1. Using Comma Initializer List: Eigen provides a comma-initializer list feature to fill a matrix or vector with coefficients:
Eigen::MatrixXd mat(3,3);
Eigen::Matrix3d mat;
 
mat << 1, 2, 3,
       4, 5, 6,
       7, 8, 9;
  1. Using Zero Method: If you want to initialize a matrix to be filled with zeros:
Eigen::Matrix3d mat = Eigen::Matrix3d::Zero();
Eigen::MatrixXd mat = Eigen::MatrixXd::Zero(3,3);
Eigen::Matrix<int, 2, 3> mat = Eigen::Matrix<int, 2, 3>::Zero();
 
// You can also do it as a method
mat.setZero();
  1. Using Constant Method: To fill a matrix with a constant value:
Eigen::Matrix3d mat = Eigen::Matrix3d::Constant(1.23);
Eigen::MatrixXd mat = Eigen::MatrixXd::Constant(3,3,1.23);
 
// You can also do it as a method
mat.setConstant(1.23);
  1. Using Identity Method: To create an identity matrix:
Eigen::Matrix3d mat = Eigen::MatrixXd::Identity();
Eigen::MatrixXd mat = Eigen::MatrixXd::Identity(3,3);
 
// You can also do it as a method
mat.setIdentity();
  1. Using Random Method: To fill a matrix with random values:
Eigen::Matrix3d mat = Eigen::MatrixXd::Random();
Eigen::MatrixXd mat = Eigen::MatrixXd::Random(3,3);
 
// You can also do it as a method
mat.setRandom();

Const matrix

https://stackoverflow.com/questions/25999407/initialize-a-constant-eigen-matrix-in-a-header-file

There are at least two possibilities. The first one is using the comma initializer features of Eigen:

Eigen::Matrix3d A((Eigen::Matrix3d() << 1, 2, 3, 4, 5, 6, 7, 8, 9).finished());
  • This is row-wise

The second is using the Matrix3d(const double*) constructor which copies data from a raw pointer. In this case, the values must be provided in the same order than the storage order of the destination, so column-wise in most cases:

const double B_data[] = {1, 4, 7, 2, 5, 8, 3, 6, 9};
Eigen::Matrix3d B(B_data);

This is stored Column-Major!

This is because the default storage order is column first.

More

Basic matrix operations

Matrix3d m1, m2;
m1 + m2; // element-wise addition
m1 - m2; // element-wise subtraction
m1 * m2; // DIFFERENT: matrix multiplication
m1 / m2; // element-wise division
 
m1.cwiseProduct(m2); // element-wise multiplication

Make sure both matrices are of the SAME type!

In Eigen, for performance reasons, you must explicitly convert the matrix type. And if you forget to do this, Eigen will (not very friendly) prompt you with a very long “YOU MIXED DIFFERENT NUMERIC TYPES …” compilation error.

matrix_33.transpose();
matrix_33.sum();
matrix_33.trace();
10 * matrix_33;
matrix_33.inverse();
matrix_33.determinant();

Casting values

matrix_23.cast<double>()
Solving

Solving things (covered in Visual SLAM book), I don’t quite understand this yet. Okay I understand, it’s just solving for

  • There is a solution when is invertible,
MatrixXd matrix_NN = MatrixXd::Random(50, 50);
matrix_NN = matrix_NN * matrix_NN.transpose(); // guarantee semi-positive definite
 
VectorXd v_Nd = VectorXd::Random(50);

There are 3 ways to solve for x

  1. use the inverse (SLOW)
MatrixXd x = matrix_NN.inverse() * v_Nd;
cout << x.transpose() << endl;
  1. QR Decomposition
x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
cout << x.transpose() << endl;
  1. Cholesky decomposition (works for positive definite matrices)
x = matrix_NN.ldlt().solve(v_Nd);
cout << x.transpose() << endl;

CheatSheet

Taken from https://zhangboyi.gitlab.io/post/2018-06-28-eigen-cheat-sheet/

// These are the typedefs provided by Eigen;
 
// fixed size
typedef Matrix<float, 4, 4> Matrix4f;
typedef Matrix<int, 1, 2> RowVector2i;
 
// dynamic size
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<int, Dynamic, 1> VectorXi; // column vector
 
 
// The matrix is stored column-major
 
// constructors
Matrix3f a;
MatrixXf b;
MatrixXf a(10,15);
VectorXf b(30);
 
// constructors to initialize value of small fixed-size vectors
Vector4d c(5.0, 6.0, 7.0, 8.0);
 
// coefficient accesors
MatrixXd m(2,2);
m(0,0) = 3;
m(1,0) = 2.5;
 
VectorXd v(2);
v(0) = 4;
v(1) = v(0) - 1;
 
// comma-initialization
Matrix3f m;
m << 1, 2, 3,
     4, 5, 6,
     7, 8, 9;
std::cout << m;
/* prints:
1 2 3
4 5 6
7 8 9
*/
 
// resizing
MatrixXd m(2,5);
m.resize(4,3);
std::cout << "The matrix m is of size " << m.rows() << "x" << m.cols() << std::endl;
std::cout << "It has " << m.size() << " coefficients" << std::endl;
/* prints:
The matrix m is of size 4x3
It has 12 coefficients
*/
 
// matrix addition, subtraction, and multiplication
Matrix2d a;
a << 1, 2,
     3, 4;
MatrixXd b(2,2);
b << 2, 3,
     1, 4;
a + b;
a - b;
a += b;
a * b;
a = a * a; // no aliasing issue
// scalar multiplication
a * 2.5   // a + 2.5 will fail
    
// transpose
a.transpose();
a.transposeInPlace();  // a = a.transpose() would give aliasing issue
 
// vector dot product and cross product
Vector3d v(1,2,3);
Vector3d w(0,1,2);
v.dot(w);
v.cross(w);
 
// arithmetic reduction operations
a.sum();		// 10
a.prod();		// 24
a.mean();		// 2.5
a.minCoeff();	// 1
a.maxCoeff();	// 4
a.trace();		// 5
std::ptrdiff_t i, j;
a.minCoeff(&i, &j); // 0, 0
 
// coefficient-wize operation
// array is like matrix, but is for coefficient-wize operations
MatrixXf m(2,2);
MatrixXf n(2,2);
m << 1,2,
     3,4;
n << 5,6,
     7,8;
m.array() * n.array();
/* result:
5  12
21 32
*/
m.cwiseProduct(n);
/* result:
5  12
21 32
*/
m.array() + 4;
/* result:
5 6
7 8
*/
 
// Interfacing with raw buffers: the Map class
Map<MatrixXf> mf(pf,rows,columns);  // pf is a float* pointing to the array of memory
Map<const Vector4i> mi(pi);
 
int array[8];
for(int i = 0; i < 8; ++i) array[i] = i;
cout << "Column-major:\n" << Map<Matrix<int,2,4> >(array) << endl;
cout << "Row-major:\n" << Map<Matrix<int,2,4,RowMajor> >(array) << endl;
cout << "Row-major using stride:\n" << Map<Matrix<int,2,4>, Unaligned, Stride<1,4> >(array) << endl;
/* prints:
Column-major:
0 2 4 6
1 3 5 7
Row-major:
0 1 2 3
4 5 6 7
Row-major using stride:
0 1 2 3
4 5 6 7
*/
 
// Reductions, visitors and broadcasting
// https://eigen.tuxfamily.org/dox/group__TutorialReductionsVisitorsBroadcasting.html

Summary

Common Issues

http://library.isr.ist.utl.pt/docs/roswiki/eigen(2f)Troubleshooting.html

When using fixed-size vectorizable Eigen types in STL containers, you must use an aligned allocator. Eigen provides one:

// Instead of
std::vector<Vector4f> v; // WRONG
// Use
std::vector<Vector4f, Eigen::aligned_allocator<Vector4f>> m;
  • umm, even this is not perfectly right?

http://library.isr.ist.utl.pt/docs/roswiki/eigen(2f)Troubleshooting.html