Common Operations with RcppArmadillo

Overview

This post is a mind dump of working with the Armadillo library inside of R.

Software Requirements for RcppArmadillo

Unlike the majority of R packages that exist, the RcppArmadillo package requires a development toolchain to use. The toolchain varies from platform to platform. Please see one of the following guides for help:

Finally, please make sure to install the RcppArmadillo package via:

install.packages("RcppArmadillo")

Verify it works by running:

Rcpp::evalCpp("2 + 2", depends = "RcppArmadillo")
#> [1] 4

Common errors:

On macOS, if you receive an error such as -lgfortran and so on, then please make sure you installed the fortran binaries.

Using RcppArmadillo from a C++ script with sourceCpp()

The preferred method of compiling C++ code is to place it in a separate file and use sourceCpp("path/to/file.cpp") to compile and import C++ functions.

When writing C++ scripts that use armadillo, every C++ code file must have at the top:

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
  1. The #include<RcppArmadillo.h> statement provides the Rcpp.h and armadillo.h headers with the appropriate casting magic.
  2. The // [[Rcpp::depends(RcppArmadillo)]] statement adds onto the compiler’s search path the RcppArmadillo package directory where the armadillo.h header files are found. This approach takes advantage of Rcpp Attributes.

Optional After the Rcpp depends statement, you are able to use:

using namespace arma;

This avoids having code prefixed by a namespace. The downside to this approach is sometimes function names can overlap making for a fun evening of debugging.

Example

add_two.cpp

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

arma::vec add_two(arma::vec x){
  return x + 2;
}

add_two.R

# Load package
library("Rcpp")

# Compile C++ code and import function into R
sourceCpp("path/to/add_two.cpp")

# Example C++ function call from R
x = c(1,2,3)
add_two(1,2,3)

Examples of R syntax and conceptually corresponding Armadillo syntax

Based off of Armadillo’s syntax table with the left column containing R commands instead of Matlab/octave.

R Armadillo Notes
A[1, 1] A(0, 0) indexing in Armadillo starts at 0
A[n, p] A(n-1, p-1)
nrow(A) A.n_rows read only
ncol(A,2) A.n_cols
dim(Q)[3] Q.n_slices Q is a cube (3D array)
length(A) A.n_elem
A[, k] A.col(k) this is a conceptual example only; exact conversion from R to Armadillo syntax will require taking into account that indexing starts at 0
A[k, ] A.row(k)
A[, p:q] A.cols(p, q)
A[p:q,] A.rows(p, q)
Q[,,k] Q.slice(k) Q is a cube (3D array)
Q[,, t:u] Q.slices(t, u)
t(A) A.t() or trans(A) matrix transpose / Hermitian transpose (for complex matrices, the conjugate of each element is taken)
A = numeric(length(A)) A.zeros()
A = rep(1, length(A)) A.ones()
A = matrix(0, k, k) A = zeros(k,k)
matrix(1, k, k) A = ones(k,k)
A * B A % B element-wise multiplication
A / B A / B element-wise division
solve(A, B) solve(A,B) conceptually similar to inv(A)*B, but more efficient
A = A + 1 A++
A = A - 1 A–
matrix(c(1,2,3,4), nrow = 2) A << 1 << 2 << endr << 3 << 4 << endr; element initialisation, with special element endr indicating end of row
X = c(A) X = vectorise(A)
X = cbind(A, B) X = join_horiz(A,B)
X = rbind(A, B) X = join_vert(A,B)
print(A) Rcpp::Rcout << A << std::endl;

Row and Column Binds

When working with C++ objects, often there is a need to combine objects together. Traditionally, in R, there is an unlimited amount of times objects can be bound together. The same – unfortunately – cannot be said for C++. In particular, only two objects may be combined simultaneously. The objects can be joined together using either a row or a column matching scheme:

#include<RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// Performs the operation of rbind
// [[Rcpp::export]]
arma::mat bind_row_vectors(arma::vec A, arma::vec B){
  arma::mat out = join_vert(A, B); // same as join_cols
  return out;
}

// Performs the operation of cbind
// [[Rcpp::export]]
arma::mat bind_col_vectors(arma::vec A, arma::vec B){
  arma::mat out = join_horiz(A, B); // same as join_rows
  return out;
}
# Define two vectors
a = 1:3
b = 4:6

# Perform a row bind operation, e.g. rbind
bind_row_vectors(a, b)
##      [,1]
## [1,]    1
## [2,]    2
## [3,]    3
## [4,]    4
## [5,]    5
## [6,]    6
# Perform a column bind operation, e.g. cbind
bind_col_vectors(a, b)
##      [,1] [,2]
## [1,]    1    4
## [2,]    2    5
## [3,]    3    6

Subsetting with conditions

Accessing a portion of data or modifying values in place is possible by using one of the find*() variants:

  • which(X > 2): find(X > 2)
  • is.finite(X): find_finite(X)
    • Values excluding: -Inf, +Inf and NaN.
    • NB NA values in the integer case, e.g. ivec/imat, transform to INT_MIN. This method misses detecting this value. Use find(X == INT_MIND) to capture such values.
  • is.infinite(x) & is.nan(x) or !is.finite(x): find_nonfinite(X)
    • Values that contain: -Inf, +Inf and NaN.
  • seq_along(X)[!duplicated(X)]: find_unique(X)
    • Element index of the first instance of a value in a vector.
#include<RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// Obtains the subset of values matching condition
// [[Rcpp::export]]
arma::vec find_equal(arma::vec A, double b){
  arma::uvec idx = find(A == b); // Substitute == with >, <, >=, <=, !=
  arma::vec out = A.elem(idx);   // Retrieve elements from positional index
  return out;
}

// [[Rcpp::export]]
arma::vec find_and_replace(arma::vec A,
                           double find_val = 3.0,
                           double replace_val = 1.5) {
  arma::uvec idx = find(A >= find_val); // Substitute == with >, <, >=, <=, !=
  A.elem(idx).fill(replace_val);        // Retrieve elements from positional index
                                        // Fill with the appropriate value
  return A;
}
# Define a vector
A = c(1,4,3,3,4,2,3,2)
b = 3

# Locate elements
find_equal(A, b)
##      [,1]
## [1,]    3
## [2,]    3
## [3,]    3
find_and_replace(A, b)
##      [,1]
## [1,]  1.0
## [2,]  1.5
## [3,]  1.5
## [4,]  1.5
## [5,]  1.5
## [6,]  2.0
## [7,]  1.5
## [8,]  2.0