A Quick Introduction to the Eigen Library

Eigen is a C++ template library for linear algebra, widely used for its performance and convenience. This article overviews key Eigen features, focusing on those used in my blog post about automatic differentiation. It's designed for readers new to the library who want a quick start without diving deep into the documentation. For a more comprehensive tutorial on Eigen, please refer to the official documentation on dense matrices.

The Matrix Class

At the core of Eigen is the Matrix class, defined in the Eigen/Dense header, which represents a dense column-major matrix.

Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>

The template arguments are:

1. Scalar: The data type of the matrix elements (e.g., float, double, int)
2. RowsAtCompileTime: The number of rows (use -1 for dynamic size)
3. ColsAtCompileTime: The number of columns (use -1 for dynamic size)

For convenience, Eigen provides several type aliases:

• MatrixX<typename Scalar> : Equivalent to Matrix<typename Scalar, -1, -1> (dynamic size)
• MatrixXf, MatrixXd, MatrixXi : Dynamic-sized matrices of float, double, and int respectively

Basic Usage

Let's explore some fundamental operations with Eigen matrices:

// Fixed-size matrix
Eigen::Matrix<float, 2, 2> a{{1, 2}, {3, 4}};
std::cout << "a is \n" << a << "\n";

// Dynamic-sized matrix
Eigen::Matrix<float, -1, -1> b(4, 3);
b.resize(2, 2);
b.setZero();
b(0, 0) = 2;
b(1, 1) = 2;
std::cout << "b is \n" << b << "\n";

// Type casting
Eigen::MatrixX<double> c = b.cast<double>();

// Matrix multiplication
Eigen::MatrixX<float> d = a * b;
std::cout << "Eigen works: " << std::boolalpha 
            << (d == Eigen::MatrixX<float>{{2, 4}, {6, 8}}) 
            << "\n";

// Complex expression
Eigen::MatrixX<float> e = 
    1.23 * a * b + Eigen::MatrixX<float>::Constant(2, 2, 1.18);
std::cout << "The result is \n" << e << "\n";

// Lazy evaluation
auto f = 2 * (a + b) + 4.56 * (a - b);
Eigen::MatrixX<float> g = f;
std::cout << "f is \n" << g << "\n";

// Reduction methods
std::cout << "f.sum() is " << g.sum() << "\n";
std::cout << "f.prod() is " << g.prod() << "\n";
std::cout << "The mean of f is " << g.mean() << "\n";
std::cout << "The norm of f is " << g.norm() << "\n";

// Element-wise operations
Eigen::MatrixX<float> h = a.array() * b.array();
Eigen::MatrixX<float> i = h.array() + 7.89;

// Coefficient-wise functions
Eigen::MatrixX<float> j = i.cwiseAbs() + i.cwiseInverse();

This code demonstrates matrix creation, basic operations, lazy evaluation, reduction methods, and element-wise operations.

a is 
1 2
3 4
b is 
2 0
0 2
Eigen works: true
The result is 
 3.64   6.1
 8.56 11.02
f is 
 1.44 13.12
19.68 21.12
f.sum() is 55.36
f.prod() is 7852.63
The mean of f is 13.84
The norm of f is 31.7422

Efficient Function Arguments

When passing matrices to functions, it's important to consider efficiency. Here's an example of an inefficient approach:

void f(const Eigen::MatrixX<float>& m)
{
// some math
}

int main()
{
    Eigen::MatrixXf m = Eigen::MatrixXf::Random(4, 4);
    f(m.topLeftCorner(3,3));
}

In this case, topLeftCorner returns a type different of MatrixXf, causing an implicit copy.

A better alternative is to use the Eigen::Ref class, which can either represent a MatrixX or a block of it.

// A function that can modify the matrix passed to it
void f(Eigen::Ref<Eigen::MatrixX<float>> m)
{
    // some math
}

// A function that takes a const Ref 
void g(Eigen::Ref<const Eigen::MatrixX<float>> m)
{
    // some math
}

int main()
{
    Eigen::MatrixXf m = Eigen::MatrixXf::Random(4, 4);
    f(m);
    g(m.topLeftCorner(3,3));
    g(m + Eigen::MatrixXf::Identity(4, 4));
}

The first call to g avoids copying, while the second creates a temporary due to the mismatched memory layout.

For more detailed information and advanced usage, I highly recommend reading the Eigen documentation on this topic.