C++ - Armadillo

Table of Contents

Introduction

Armadillo1 is an excellent C++ library, which provides efficient linear algebra related operations, based on a series of classes, e.g., for vectors, matrices, etc. Thanks to its friendly interfaces, similar to MATLAB, the library is very easy to use.

Armadillo is developped by Conrad Sanderson and Ryan Curtin, and released under Apache 2.0 license.

Installation

Matrix decomposition related operations are developped based on the existing third-party libraries, e.g., the open source OpenBLAS or LAPACK, or even the commercial MKL. In my opinion, openBLAS is the best choice, and can be installed in ArchLinux. 

pacman -S openblas

The latest codes resides in a git repository2. Based on the source codes can be compiled and installed as usual.

Besides, the library can be easily installed as system package. Taking ArchLinux for instance, armadillo can be installed by 

pacman -S armadillo

Usage

Note that the detailed usage of the libraries can be found in its online manual, and this post just provides a brief and incomplete summary as a quick lookup reference after my walking through the manual.

Typedef for basic elements

  • uword is a typedef for unsigned int.
  • sword is a typedef for (signed) int.
  • cx_float is a typedef for std::complex<float>.
  • cx_double is a typedef for std::complex<double>.

Fill form

The library defines a series of fill_form, which significantly facilitates the initialization of the data type classes.

  • fill::zeros sets all elements to 0.
  • fill::ones sets all elements to 1.
  • fill::eye initializes a matrix as an identity matrix.
  • fill::randu sets all elements to uniform distributed random values, i.e., \(\mathbb{U}(0, 1)\).
  • fill::randn sets all elements to standard Gaussian distributed random values, i.e., \(\mathbb{N}(0, 1)\).
  • fill::value(scalar) sets all elements to a specified scalar.
  • fill::none does not initialize the elements.

Classes

Matrix

Considering the type of component elements, matrix is defined as a template Mat<type>, where type can be signed or unsigned, real or complex, int or float or even double. A matrix is column-wise stored in the memory.

Typedef

Moreover, following types have been defined for convenience.

  • imat is a typedef for Mat<sword>.
  • umat is a typedef for Mat<uword>.
  • fmat is a typedef for Mat<float>.
  • mat or dmat are both typedefs for Mat<double>.
  • cx_fmat is a typedef for Mat<cx_float>.
  • cx_mat and cx_dmat are both typedefs for Mat<cx_double>.
Constructor

A rich number of constructors have been defined, e.g., taking mat for instance,

  • mat()
  • mat(n_rows, n_cols)
  • mat(n_rows, n_cols, fill_form)
  • mat::fixed<n_rows, n_cols> creates a matrix with fixed size.
  • mat::fixed<n_rows, n_cols>(fill_form)
  • mat(size(X))
  • mat(size(X), fill_form)
  • mat(mat)
  • mat(vec)
  • mat(rowvec)
  • mat(initializer_list)
  • mat(string)
  • mat(std::vector) creates a column vector, as a single-column matrix.
  • mat(sp_mat) creates a dense matrix based on a sparse matrix.
  • cx_mat(mat, mat) constructs a complex matrix based on a pair of real matrices, as real and imaginary parts respectively.
Submatrix view
Continuous
  • X.col(idx) / X.row(idx)
  • X.cols(idx1, idx2) / X.rows(idx1, idx2)
  • X.submat(row1, col1, row2, col2)
  • X(span(row1, row2), span(col1, col2))
  • X(row, col, size(n_rows, n_cols)) / X(row, col, size(Y))
  • X(span(row1, row2), col)
  • X(row, span(col1, col2))
  • X.head_cols(n_cols) / X.head_rows(n_rows)
  • X.tail_cols(n_cols) / X.tail_rows(n_rows)

span defines a continuous range, e.g.,

  • span(a, b) describes a contiguous range from \(a\) to \(b\), \(a < b\).
  • span(a) indicates a particular row/column/slice.
  • span() or span::all indicates the entire range.
Non-continuous
  • X.cols(idx_vec) / X.rows(idx_vec)
  • X.submat(row_vec, col_vec)
  • X(row_vec, col_vec)
Member functions
  • .eye() : Set the matrix to an identity matrix.
  • .diag(k = 0) : Return a column vector of a matrix.
  • .as_col() / .as_row() : Return a column/row vector by concatenating all the columns/rows of the matrix.
  • .each_col(idx_vec) / .each_row(idx_vec) / .each_col(lambda) / .each_row(lambda) : Apply a vector operation to each column/row of a matrix. The argument idx_vec is optional, and by default all columns/rows are used. If the lambda function lambda is specified, then the function must accept a reference to a Col<type> or Row<type> object with the same element type as the underlying matrix. Following operations have been supported.
    • +/-// : Addition/subtraction/division
    • % : Element-wise multiplication
    • = : Assignment (copy)
    • +=/-=//= : In-place addition/subtraction/division
    • %= : In-place element-wise multiplication
  • .st() : Return a transposed copy (only applicable for a complex matrix)
  • .t() : Return a Hermitian transposed copy3.
  • .i() : Return the inverse for a square matrix.
  • .is_square()
  • .is_vec() : Return true if the matrix is a column or row vector.
  • .is_colvec / .is_rowvec
  • .is_trimatu() / .is_trimatl()
  • .is_diagmat()

Vector

In the library, column vectors and row vectors are different, represented by Col<type> and Row<type> respectively.

Typedef
Column vector
  • icolvec and ivec are both typedefs for Col<sword>.
  • ucolvec and uvec are both typedefs for Col<uword>.
  • fcolvec and fvec are both typedefs for Col<float>.
  • dcolvec, dvec, colvec and vec are all typedefs for Col<double>.
  • cx_fcolvec and cx_fvec are both typedefs for Col<cx_float>.
  • cx_dcolvec, cx_dvec, cx_colvec and cx_vec are all typedefs for Col<cx_double>.
Row vector
  • irowvec is a typedef for Row<sword>.
  • urowvec is a typedef for Row<uword>.
  • frowvec is a typedef for Row<float>.
  • drowvec and rowvec are both typedefs for Row<double>.
Constructor

Taking vec/rowvec for instance, following constructs have been defined.

  • vec/rowvec()
  • vec/rowvec(n_elem)
  • vec/rowvec(n_elem, fill_form)
  • vec/rowvec::fixed<n_elem> creates a column/row vector with fixed size.
  • vec/rowvec::fixed<n_elem>(fill_form)
  • vec/rowvec(size(X))
  • vec/rowvec(size(X), fill_form)
  • vec/rowvec(vec/rowvec)
  • vec/rowvec(mat) creates a column/row vector based on a single-column/row matrix.
  • vec/rowvec(initializer_list)
  • vec/rowvec(string) creates a column/row vector based on a string with space as delimitor.
  • vec/rowvec(std::vector) creates a column/row vector based on a vector container in STL.
  • cx_vec/rowvec(vec/rowvec, vec/rowvec) constructs a complex column/row vector based on a pair of real column/row vectors, as real and imaginary parts respectively.
Subvector view
Continuous
  • V(span(idx1, idx2))
  • V.subvec(idx1, idx2) / V.subvec(idx, size(U))
  • V.head(n_elem) / V.tail(n_elem)
Non-continuous
  • X.elem(idx_vec) / X(idx_vec)

Cube

A cube is essentially a 3-dimension (3D) matrix, which is compriesd of a series of 2D matrices, a.k.a. slices. Cube class is also a template, i.e., Cube<type>.

Typedef
  • icube is a typedef for Cube<sword>.
  • ucube is a typedef for Cube<uword>.
  • fcube is a typedef for Cube<float>.
  • dcube and cube are both typedefs for Cube<double>.
  • cx_fcube is a typedef for Cube<cx_float>.
  • cx_dcube and cx_cube are both typedefs for Cube<cx_double>.
Constructor
  • cube()
  • cube(n_rows, n_cols, n_slices)
  • cube(n_rows, n_cols, n_slices, fill_form)
  • cube::fixed<n_rows, n_cols, n_slices> creates a cube with fixed size.
  • cube(size(X))
  • cube(size(X), fill_form)
  • cube(cube)
  • cx_cube(cube, cube) constructs a complex cube based on a pair of real cubes, as real and imaginary parts respectively.
Subcube view
Contiguous
  • Q.slice(idx) / Q.slices(idx1, idx2)
  • Q.col_as_mat(idx) / Q.row_as_mat(idx) returns a flattened matrix.
  • Q.row(idx) / Q.col(idx)
  • Q.rows(idx1, idx2) / Q.cols(idx1, idx2)
  • Q.subcube(row1, col1, slice1, row2, col2, slice2)
  • Q(span(row1, row2), span(col1, col2), span(slice1, slice2))
  • Q(row1, col1, slice1, size(n_rows, n_cols, n_slices)) / Q(row, col, slice, size(R))
  • Q.head_slices(n_slices) / Q.tail_slices(n_slices)
  • Q.tube(row, col)4 / Q.tube(row1, col1, row2, col2) / Q.tube(span(row1, row2), span(col1, col2)) / Q.tube(row, col, size(n_rows, n_cols))
Non-contiguous
  • Q.elem(idx_vec)
  • Q(idx_vec)
  • Q.slices(idx_vec)
Member functions
  • .each_slice(idx_vec) / .each_slice(lambda) / .each_slice(lambda, use_mp = false): Apply a matrix operation to each slice/matrix of a cube. The argument idx_vec is optional, and by default all slices are used. If the lambda function lambda is specified, then the function must accept a reference to a Mat<type> object with the same element type as the underlying cube. use_mp is a boolean indicator of OpenMP for multithreading. If it is set to true, then the lambda must be thread-safe, i.e., cannot write to variables outside of its scope.
  • .insert_slices(idx, X) : Insert X at idx th slice.
  • .insert_slices(idx, slice_num) : Insert slice_num slices at idx th slice, and the slice_num slices inserted are all set to 0.
  • .shed_slice(idx) : Remove idx th slice from the underlying cube.
  • .shed_slices(first, last) : Remove the slices between first th and last th slices.
  • .shed_slices(idx_vec) : Remove the slices with indices in idx_vec.

Field

According to the online manual, a field is a class to store arbitrary objects in matrix-like or cube-like layouts. Given a vector or a matrix or a cube, all the elements (e.g., scalars in a vector, vectors in a matrix, matrices in a cube) must have the identical type and equal size. But in a field, the elements just have the identical type, but their lengths or sizes can be different from one element to another.

Constructor

A field object can be created by any one of the constructors below, where object_type is the element type, which can be vector or matrix or cube.

  • field<object_type>()
  • field<object_type>(n_elem)
  • field<object_type>(n_rows, n_cols)
  • field<object_type>(n_rows, n_cols, n_slices)
  • field<object_type>(size(X))
  • field<object_type>(field<object_type>)
Subfield view
For a 2D field
  • F.row(idx)
  • F.col(idx)
  • F.rows(row1, row2)
  • F.cols(col1, col2)
  • F.subfield(row1, col1, row2, col2)
  • F(span(row1, row2), span(col1, col2))
  • F(row, col, size(G))
  • F(row, col, size(n_rows, n_cols))
For a 3D field
  • F.slice(idx)
  • F.slices(slice1, slice2)
  • F.subfield(row1, col1, slice1, row2, col2, slice2)
  • F(span(row1, row2), span(col1, col2), span(slice1, slice2))
  • F(row, col, slice, size(G))
  • F(row, col, slice, size(n_rows, n_cols, n_slices))

Length and size

Length and size related attributes of type uword, can be listed as below.

  • .n_elem : Total number of elements in a vector/matrix/cube/field
  • .n_cols : Number of columns in a vector/matrix/cube/field
  • .n_rows : Number of rows in a vector/matrix/cube/field
  • .n_slices : Number of slices in a cube/field
  • .n_nonzeros : Number of non-zero elements in a sparse matrix

The dimensions above are all read-only, but can be changed by

  • .set_size(n_elem) for a vector/field
  • .set_size(n_rows, n_cols) for a matrix/field
  • .set_size(n_rows, n_cols, n_slices) for a cube/field
  • .set_size(size(X)) for a vector/matrix/cube/field
  • .copy_size(X) set the size equal to that of object X.

Common member functions

  • .zeros() : Set all elements to 0.
  • .ones() : Set all elements to 1.
  • .randu() : Set all elements to uniform distributed random values, i.e., \(\mathbb{U}(0, 1)\).
  • .randn() : Set all elements to standard Gaussian distributed random values, i.e., \(\mathbb{N}(0, 1)\).
  • .fill(value) : Set all the elements to value.
  • .replace(old_value, new_value) : For all elements equal to old_value, set them to new_value.
  • .clean(threshold) : Replace each element with absolute value no larger than threshold by zero. For each complex element, its real and imaginary parts are individually replaced, irrespective of its absolute value and amplitude.
  • .reset() : Reset the size to zero, i.e., empty the object.
  • .set_size(n_elem) / .set_size(n_rows, n_cols) / .set_size(n_rows, n_cols, n_slices) / .set_size(size(X)) : Change the size of the object without preserving data.
  • .copy_size(X) : Set the size of the object equal to that of object X without preserving data.
  • .resize(n_elem) / .resize(n_rows, n_cols) / .resize(n_rows, n_cols, n_slices) / .resize(size(X)) : Grow/shrink the object while preserving the elements and their layout.
  • .reshape(n_rows, n_cols) / .reshape(n_rows, n_cols, n_slices) / .reshape(size(X)) : Recreate the object according to the given size using the elements in the previous version. The recreating is performed column by column.
  • .imbue(functor/lambda) : Fill/set with values yielded by a functor/lambda expression.
  • .transform(functor/lambda) : Transform each element using a functor/lambda expression.
  • .for_each(functor/lambda) : Process the object by passing each element to a functor/lambda expression.
  • .set_real(X) / .set_imag(X) : Set the real/imaginary part of an object. X must have the same size as the underlying object.
  • .insert_rows(idx, X) : Insert X at idx th row for a column-vector/matrix/cube.
  • .insert_rows(idx, row_num) : Insert row_num rows at idx th row for a column-vector/matrix/cube, and the row_num rows inserted are all set to 0.
  • .insert_cols(idx, X) : Insert X at idx th column for a row-vector/matrix/cube.
  • .insert_cols(idx, col_num) : Insert col_num columns at idx th column for a column-vector/matrix/cube, and the col_num columns inserted are all set to 0.
  • .shed_row(idx) : Remove idx th row from the underlying column-vector/matrix/cube.
  • .shed_rows(first, last) : Remove the rows between first th and last th rows from the underlying column-vector/matrix/cube.
  • .shed_col(idx) : Remove idx th column from the underlying row-vector/matrix/cube.
  • .shed_cols(first, last) : Remove the columns between first th and last th columns from the underlying row-vector/matrix/cube.
  • .swap_rows(idx1, idx2) : Swap 2 rows of the underlying column-vector/matrix.
  • .swap_cols(idx1, idx2) : Swap 2 columns of the underlying row-vector/matrix.
  • .min() / .max() : Return the extreme5 value of a matrix or a cube.
  • .index_min() / .index_max() : Return the linear index of the extreme5 value of a matrix or a cube.
  • .in_range : Check the validity of a given location, e.g., return false for an empty object, out-of-bounds.
  • .is_empty() : Check the object whether it is empty.
  • .is_sorted(sort_direction = "ascend", dim = 0) : Check a vector or a matrix.
  • .print() / .print(header) / .print(stream) / .print(stream, header) : Print the content of the object.
  • .save / .load

Notes: Functions .zeros, .ones, .randu, .randn, .eye above can also take dimensions as arguments to set the length and size of the object. In this case, it combines the operations of argument-free provoking and .set_size.

Overloaded operators

A rich number of operators have been overloaded for column/row vector, matrix and cube classes, e.g.,

  • * : Matrix multiplication of two objects, but not applicable to cube, unless multiplying by a scalar.
  • % : Element-wise multiplication, a.k.a. Schur product.
  • / : Element-wise division of an object by another object or a scalar.
  • == / != : Element-wise equality/non-equality evaluation
  • > / < / >= / <= : Element-wise comparsion
  • && / || : Element-wise logical AND/OR evaluation

Element access

  • (n) : Access the \(n\) th element of a vector/matrix/cube or \(n\) th object of a field.
  • .at(n) / [n] : Similar to (n), but without a bound check6.
  • (i, j) : Access the element in \(i\) th row and \(j\) th column of a matrix, or the object in \(i\) th row and \(j\) th column of a 2D field.
  • .at(i, j) : Similar to (i, j), but without a bound check.
  • (i, j, k) : Access the element in \(i\) th row, \(j\) th column, and \(k\) th slice; or the object in \(i\) th row, \(j\) th column, and \(k\) th slice for a 3D field.
  • .at(i, j, k) : Similar to (i, j, k), but without bound check.

Note that, 2D index using rectangle bracket (e.g., \([i, j]\), \([i, j, k]\)) does not work correctly. We should use \((i, j)\), \((i, j, k)\) instead.

Iterators

In armadillo, matrices/vectors/cubes are somewhat similar to the containers in C++ STL, and consequently can be accessed via iterators.

Member functions
  • .begin() : Iterator referring to the first element.
  • .end() : Iterator referring to the past-the-end element.
  • .begin_col(idx) / .begin_row(idx) : Iterator referring to idx th column/row of a row/column vector or a matrix.
  • .end_col(idx) / .end_row(idx) : Iterator referring to the past-the-end of idx th column/row of a row/column vector or a matrix.
  • .begin_slice(idx) : Iterator referring to idx th slice of a cube.
  • .end_slice(idx) : Iterator referring to the past-the-end of idx th slice of a cube.
Iterator types
For read/write access
  • mat::iterator / vec::iterator / rowvec::iterator / cube::iterator
  • mat::col_iterator / vec::col_iterator / rowvec::col_iterator
  • mat::row_iterator / vec::row_iterator / rowvec::row_iterator
  • cube::slice_iterator
For read-only access
  • mat::const_iterator / vec::const_iterator / rowvec::const_iterator / cube::const_iterator
  • mat::const_col_iterator / vec::const_col_iterator / rowvec::const_col_iterator
  • mat::const_row_iterator / vec::const_row_iterator / rowvec::const_row_iterator
  • cube::const_slice_iterator

Container-compatible functions

  • .front() / .back() : Access the first/last element in a vector.
  • .clear() : Empty the object.
  • .empty() : Return true/false for an empty/non-empty object.
  • .size() : Total number of elements in the object, i.e., .n_elem.

Generation

  • ones(size(X)) / zeros(size(X))
  • randu / randi / randn / randg

Matrix

  • eye(n_rows, n_cols) / eye(size(X))
  • ones(n_rows, n_cols) / zeros(n_rows, n_cols)
  • toeplitz / circ_toeplitz
  • diagmat(X, k=0)
  • trimatu(X, k) / trimatl(X, k)

Vector

  • ones(n_elem) / zeros(n_elem)
  • linspace(start, end, N=100)
  • logspace(a, b, N=50): Generate a logarithmically spaced vector from \(10^a\) to \(10^b\) (included).
  • regspace(start, end) / regspace(start, step, end)
  • randperm(N) / randperm(N, M) : Generate a vector with random permutation of integers from 0 to N-1.
  • diagvec(X, k=0)

Cube

  • ones(n_rows, n_cols, n_slices) / zeros(n_rows, n_cols, n_slices)

Functions

  • abs(X) / arg(X) : Per-element absolute values or phase angles (in radians)
  • conj(X) : Per-element complex conjugate
  • accu(X) : Sum all elements.
  • sum(X, dim=0)
  • all(V) / all(X, dim) : Check whether all the elements (in dimension dim) are non-zero.
  • any(V) / any(X, dim) : Check whether there is any non-zero element (dimension dim) exist.
  • conv_to<type>::from(X) : Type convertion/cast
  • diff(V, k=0) / diff(X, k, dim=0)
  • real(X) / imag(X)
  • min / max
  • norm(X, p=2)
  • pow
  • shuffle(X, dim=0)
  • sort / sort_index
  • unique(X)
  • vectorise(X, dim=0) : Concatenate all the columns (dim=0) or rows (dim=1).
  • fft / ifft
  • cor / cov : Correlation/covariance
  • Statistics: mean, median, stddev, var, range, hist / histc, quantile.

Matrix

  • as_scalar(expression) : Convert the expression (resulting a 1x1 matrix) to a scalar.
  • cond(X) / rcond(X) : Calculate the (reciprocal) condition number of a matrix.
  • det(X) / log_det(X) / log_det_sympd(X)
  • trace(X)
  • expmat(X) / expmat_sym(X) : Matrix exponential of a square (symmetric/Hermitian) matrix.
  • trans(X) / strans(X)
  • fliplr(X) / flipud(X)
  • inv / inv_sympd / pinv
  • intersect(A, B)
  • join_rows(A, B, ...) / join_horiz(A, B, ...)
  • join_cols(A, B, ...) / join_vert(A, B, ...)
  • kron(A, B)
  • fft2 / ifft2
  • normalise(X, p=2, dim=0) : Return the normalized version with unit p norm along the dimension dim.
  • powmat(X, n) / sqrtmat(X) / sqrtmat_sympd(X)
  • rank(X)
  • trimatu_ind(size(X), k) / trimatl_ind(size(X), k) : Return a column vector containing the indices of elements that form the upper/lower triangle part of matrix.
  • orth(X) : Find the orthonormal basis of the range space of the matrix.
  • conv2(a, b) : 2D convolution
  • Decomposition: chol, eig_sym / eig_gen, lu, qr, svd, svd_econ.

Vector

  • dot(a, b) : \(a \cdot b\)
  • cdot(a, b) : \(a^{*} \cdot b\)
  • norm_dot(a, b) : \(\dfrac{a \cdot b}{\|a\|\|b\|}\)
  • normalise(V, p=2) : Return the normalized version with unit p norm.
  • roots(P) : Calculate the complex roots of a polynormial function (represented by P).
  • conv(a, b) : 1D convolution

Cube

  • join_slices(A, B)

Miscellaneous

Constant

The constants are stored in Datum<type> class, where type can be float or double.

  • datum is a typedef for Datum<double>.
  • fdatum is a typedef for Datum<float>.
Expression Constant
datum::pi \(\pi\)
datum::tau \(2\pi\)
datum::inf \(\infty\)
datum::nan NaN7
datum::eps \(\epsilon\)8
datum::e \(e\)
datum::sqrt2 \(\sqrt{2}\)
datum::log_min=/=datum::log_max Log of minimum/maximum value (type and machine dependent)
datum::gratio Golden ratio
datum::F Faraday constant (in coulombs)
datum::h=/=datum::h_bar Planck constant (\(h\)) / reduced Planck constant (\(\frac{h}{2\pi}\))
datum::c_0 Speed of light in vacuum (m/s)

Timer

Class wall_clock is simple timer to measure the elapsed seconds. It has 2 member functions.

  • .tic() starts the timer.
  • .toc() returns the number of seconds since the last call to .tic().

Footnotes:

1

https://arma.sourceforge.net

2

https://gitlab.com/conradsnicta/armadillo-code.git

3

Reduce to transposition for real matrices.

4

Tube is a subcube across all the slices.

5

For a complex matrix/cube, absolute values are compared.

6

For this reason, it is not recommended.

7

Not a number (NaN) does not equal anything, even itself.

8

Machine epsilon, i.e., difference between 1 and the next representable value.