Third Party Libraries

VCPKG & CMake

Introduction

  • VCPKG is a package manager for C++ that simplifies the process of installing and managing libraries.
  • Provides a consistent, cross-platform build system for C++ libraries.
  • Prevent version conflicts and simplify the build process.

Install and use packages

For macOS, install pkg-config using Homebrew first:

brew install pkg-config

Then follow the instructions on vcpkg documentation to install and use packages.

Eigen

Introduction

Eigen is a C++ template library for linear algebra: matrices, vectors, numerical solvers, and related algorithms.

  • Versatile: supports dense and sparse matrices, and supports both fixed and dynamic sizes.
  • Fast: optimized for performance and efficiency.
  • Easy to use: provides a simple and intuitive API for linear algebra operations.

Quick Start

You can install Eigen either through vcpkg or download the source code then use add_subdirectory in your CMakeLists.txt.

#include <iostream>
#include <Eigen/Dense>
 
using Eigen::MatrixXd;
 
int main()
{
    MatrixXd m(2,2);
    m(0,0) = 3;
    m(1,0) = 2.5;
    m(0,1) = -1;
    m(1,1) = m(1,0) + m(0,1);
    std::cout << m << std::endl;
}

The Matrix class

The three mandatory template parameters(inside the angle brackets) of Matrix are:

  • Scalar: the type of the scalar elements in the matrix. (e.g. float, double, int)
  • int RowsAtCompileTime: the number of rows of the matrix.
  • int ColsAtCompileTime: the number of columns of the matrix.

For example:

// A 3x3 matrix of floats
Matrix<float, 3, 3> matrix;
// A matrix with a dynamic number of rows and columns
Matrix<int, Dynamic, Dynamic> matrix;

Type aliases

Eigen provides some type aliases for common matrices:1

  • MatrixNt for Matrix<type, N, N>
    • Example: MatrixXi for Matrix<int, Dynamic, Dynamic>
  • MatrixXNt for Matrix<type, Dynamic, N>
    • Example: MatrixX3i for Matrix<int, Dynamic, 3>
  • MatrixNXt for Matrix<type, N, Dynamic>
    • Example: Matrix4Xd for Matrix<double, 4, Dynamic>

Vector and RowVector

Vector by default is a column vector.

  • VectorNt for Matrix<type, N, 1>
    • Example: Vector2f for Matrix<float, 2, 1>
  • RowVectorNt for Matrix<type, 1, N>
    • Example: RowVector3d for Matrix<double, 1, 3>

Constructors

  • Matrix3f a;
    • Creates a 3x3 matrix, all the elements are uninitialized.
  • MatrixXf b
    • Creates a matrix with dynamic size, its size is currently 0-by-0, memory is not allocated yet.
  • MatrixXf c(10, 10)
    • Creates a 10x10 matrix, all the elements are uninitialized.
  • Vector2f d;
    • Creates a 2x1 vector, all the elements are uninitialized.

Here’s the way to assign values during declaration:

Matrix<int, 5, 1> b {1, 2, 3, 4, 5};
Matrix<int, 1, 5> c = {1, 2, 3, 4, 5};

MatrixXi a {      // construct a 2x2 matrix
      {1, 2},     // first row
      {3, 4}      // second row
};

Matrix<double, 2, 3> b {
      {2, 3, 4},
      {5, 6, 7},
};

Accessing elements

The primary coefficient accessors and mutators in Eigen are the overloaded parenthesis operators.

int main()
{
    Eigen::MatrixXd m(2,2);
    m(0,0) = 3;
    m(1,0) = 2.5;
    m(0,1) = -1;
    m(1,1) = m(1,0) + m(0,1);
    std::cout << "Here is the matrix m:\n" << m << std::endl;
    Eigen::VectorXd v(2);
    v(0) = 4;
    v(1) = v(0) - 1;
    std::cout << "Here is the vector v:\n" << v << std::endl;
}

Comma-initialization

Matrix and vector coefficients can be conveniently set using the so-called comma-initializer syntax.

  • To use this syntax, the matrix must have a fixed size or already be resized to the desired size before assigning values.
Matrix<int, 2, 3> a;
a << 1, 2, 3,
     4, 5, 6;
std::cout << a << std::endl;

Resizing and copying

.resize() can be used to change the size of a matrix or vector.

  • After resizing, the matrix/vector will be uninitialized.
MatrixXd a;
a.resize(4, 4); // resizes a to 4x4

Matrix3d b;
b.resize(4, 3); // resizes b to 4x3

b = a; // copies the elements of a into b
  • The copy in Eigen is a deep copy.

Slice

  • Eigen::seq(start, end, step): generates a sequence of indices starting from start to end with a step size of step.
    • step is optional, the default step size is 1.
  • Eigen::seqN(start, length, increment): generates a sequence of indices starting from start to start + length * increment with a step size of increment.
  • Eigen::indexing::all: generates a sequence of all indices.
Eigen::MatrixXf m(3, 3);
m << 1, 2, 3,
    4, 5, 6,
    7, 8, 9;
std::cout << m(Eigen::indexing::all, Eigen::seq(1, 2)) << std::endl;

Arithmetic operations

Scalar with matrix/vector:

  • +, -, *, / are applied to each element of the matrix/vector.

Matrix/vector with matrix/vector:

  • +, - are applied element-wise.
  • *, / are matrix multiplication and division.
  • dot(): dot product.
  • cross(): cross product.

Array in Eigen

The Array class provides general-purpose arrays, as opposed to the Matrix class which is intended for linear algebra.

  • Take the same template parameters as the Matrix class.
  • The arithmetic operations between Array and Matrix/Vector are element-wise.
  • .matrix() can convert an Array to a Matrix.
  • .array() can convert a Matrix/Vector to an Array.

Broadcasting

Sometimes we want to perform an operation between a matrix and a vector,2 it’s commonly used in numpy or matlab.

Suppose we want to do the following operation, each column of the matrix is added by the corresponding vector.

\[ \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ \end{bmatrix} + \begin{bmatrix} 1 \\ 2 \\ \end{bmatrix} \]

MatrixXf a(2, 3);
a << 1, 2, 3,
     4, 5, 6;
VectorXf b(2);
b << 1, 2;
std::cout << a.colwise() + b << std::endl;

Optional: Tensor in Eigen

Tensor is not officially supported in Eigen, but you can use the Eigen::Tensor class to define a tensor.

The template parameters are:

  • Scalar: the type of the scalar elements in the tensor.
  • Rank: the number of dimensions of the tensor.
Tensor<float, 3> a(10, 10, 10); // a 10x10x10 tensor
Tensor<float, 4> b(10, 3, 5, 7); // a 10x3x5x7 tensor

To perform product between tensors, you must use define the contraction dimension.

Eigen::Tensor<float, 3> tensorA(2, 3, 4); // Shape: (2, 3, 4)
Eigen::Tensor<float, 2> tensorB(4, 5);   // Shape: (4, 5)

// Contract along the 3rd axis of tensorA and 1st axis of tensorB
Eigen::array<Eigen::IndexPair<int>, 1> contraction_dims = { Eigen::IndexPair<int>(2, 0) };

// Perform contraction
Eigen::Tensor<float, 3> result = tensorA.contract(tensorB, contraction_dims);

Footnotes

  1. N can be fixed number or X for dynamic size, t can be one of i(int), f(float), d(double), cf(complex float), cd(complex double)↩︎

  2. Indeed, it can be done between a matrix with more dimensions and a matrix with fewer dimensions.↩︎