RcppArmadillo Package in R (Example)
As a data scientist, operations that fall under linear algebra are an everyday part of the job. With package RcppArmadillo, the functionality of the rich C++ Armadillo library for linear algebra becomes available in the R programming language. We show you some examples of the functionality of RcppArmadillo in the following.
This post contains the following structure:
Let’s get right to it.
RcppArmadillo
R package Rcpp allows you to integrate C++ code into R, see our blog post. For C++, there exist various libraries. One particularly useful library for data scientists is the Armadillo library (Sanderson 2010) for linear algebra and scientific computing. Luckily, R package RcppArmadillo (Eddelbuettel, Sanderson 2014) allows us to use the rich functionality of Armadillo in R.
A further advantage of RcppArmadillo is that the syntax is pretty straightforward and similar to R. On the Armadillo webpage, you can find a documentation of the different Armadillo classes and functions. In the following, we make some examples to show you the functionality of the package.
Example 1: Generalized Least Squares
We load the Rcpp and RcppArmadillo package.
if (!require('Rcpp', quietly = TRUE)) { install.packages('Rcpp') } library('Rcpp') # Load package 'Rcpp' if (!require('RcppArmadillo', quietly = TRUE)) { install.packages('RcppArmadillo') } library('RcppArmadillo') # Load package 'RcppArmadillo'
Let’s program the standard generalized least squares (GLS) solution of the fixed effects, given as \(\hat{\boldsymbol{\beta}} = (\boldsymbol{X}^{\top} \boldsymbol{V}^{-1} \boldsymbol{X})^{-1} \boldsymbol{X}^{\top} \boldsymbol{V}^{-1} \boldsymbol{y} \), where \(\hat{\boldsymbol{\beta}}\) is a \(p\)-vector of fixed effects, \(\boldsymbol{X}\) is a \(n\times p\)-dimensional model matrix of full rank, \(\boldsymbol{V}\) is a \(n\times n\)-covariance matrix, and \(\boldsymbol{y}\) is a \(n\)-vector of observations.
In R, the solution can be written as follows.
GLS_r <- function ( X, V, y) { V_inv <- solve(V) solve(t(X) %*% V_inv %*% X) %*% t(X) %*% V_inv %*% y }
We can write a C++ function for it in R as follows.
cppFunction(depends = "RcppArmadillo", ' arma::vec GLS_cpp( arma::mat X, arma::mat V, arma::vec y ) { arma::mat V_inv = inv(V); arma::vec res = inv(X.t() * V_inv * X) * X.t() * V_inv * y; return res; } ' )
In the first row, you see that we explicitly include the RcppArmadillo dependency in the code. The Armadillo vectors and matrices are defined with the arma::mat or arma::vec expression. With these arma expressions, we can conduct the matrix calculations visible in function GLS_cpp(). When comparing the R and C++ code, we see that the code looks very similar. That’s a big advantage of RcppArmadillo, programming the matrix calculations is pretty straightforward and similar to R.
To illustrate the functions, we use the example application presented in the book Faraway (2004) Linear Models with R, pages 97-98. You can take a look at the book to learn more about the data. In brief, we use the Longley regression data and regress the number of employed persons in 1947-1962 (dependent variable) on the gross national product (GNP) and the number of persons 14 years and older (independent variables). Our model input is defined as follows.
X <- model.matrix(lm(Employed ~ GNP + Population, data = longley)) V <- diag(16) V <- 0.31041^abs(row(V)-col(V)) y <- as.matrix(longley$Empl, ncol = 1)
We apply functions GLS_r() and GLS_cpp() to the data.
GLS_r (X, V, y) # [,1] # (Intercept) 94.8988949 # GNP 0.0673895 # Population -0.4742741 GLS_cpp(X, V, y) # [,1] # [1,] 94.8988949 # [2,] 0.0673895 # [3,] -0.4742741
The functions return the GLS solution of \(\hat{\boldsymbol{\beta}}\). Note that usually, vectors like \(\hat{\boldsymbol{\beta}}\) are defined as column-vectors (like a matrix with only one column).
As always, when integrating C++ with R, we are interested in whether this will save us any computation time.
library(microbenchmark) microbenchmark::microbenchmark("R" = GLS_r (X, V, y), "Cpp" = GLS_cpp(X, V, y), times = 10^{5}) # Unit: microseconds # expr min lq mean median uq max neval # R 53.2 56.500 73.56015 58.201 63.702 116693.400 1e+05 # Cpp 7.9 8.601 12.16419 10.201 11.001 8617.502 1e+05
Take a look at the difference in the mean/median time of the 10^5 repetitions. It’s huge!
Example 2: Build a Block-Diagonal Matrix From List/Array
We make some more practical examples with RcppArmadillo.
As an exercise, assume we have a list or an array of covariance matrices and we want to transform these into a block-diagonal matrix. We generate the following example 2-by-2 matrices for a list of dimension \(D=3\).
V_list <- list() D = 3 for (d in 1:D) { V_list[[d]] <- matrix(c(3, rep(0.6 * sqrt(3 * 4), 2), 4), ncol = 2, byrow = TRUE) } V_array <- array(unlist(V_list), dim = c(2, 2, D))
The last line of code transforms this list into an array of same dimension.
Now, in R we can use function bdiag() from the Matrix package which returns a sparse block-diagonal matrix of V_list.
Matrix::bdiag(V_list) # 6 x 6 sparse Matrix of class "dgCMatrix" # # [1,] 3.000000 2.078461 . . . . # [2,] 2.078461 4.000000 . . . . # [3,] . . 3.000000 2.078461 . . # [4,] . . 2.078461 4.000000 . . # [5,] . . . . 3.000000 2.078461 # [6,] . . . . 2.078461 4.000000
As an exercise, we want to define a C++ function which returns the same output.
cppFunction(depends = "RcppArmadillo", ' arma::mat cpp_bdiag(arma::cube V) { // From cube V, we get the number of slices, rows and columns int n_slices = V.n_slices; int n_rows = V.n_rows; int n_cols = V.n_cols; Rcout << "n_slices: " << n_slices << "" << std::endl; Rcout << "n_rows: " << n_rows << "" << std::endl; Rcout << "n_cols: " << n_cols << "" << std::endl; // Initialize the block-diagonal matrix V_bdiag arma::mat A = arma::mat(n_rows, n_cols); arma::mat V_bdiag = arma::repmat(A, n_slices, n_slices); Rcout << "A:\\n" << A << "" << std::endl; Rcout << "V_bdiag:\\n" << V_bdiag << "" << std::endl; // Fill in the values of cube V into V_bdiag for (int i = 0; i < n_slices; ++i) { int a = 2 * i; int b = a + 1; V_bdiag.submat(a, a, b, b) = V.slice(i); } return V_bdiag; } ' )
In this exercise, we learn several things. Function cpp_bdiag() takes a arma::cube object as input. This corresponds to a 3-dimensional array in R. Command V.n_slices gives us the number of slices of cube V. Remember that in C++, we have to initialize all object before any loop. With arma::mat(n_rows, n_cols), we initialize an empty matrix with specified number of rows and columns.
With V_bdiag.submat(a, a, b, b), we index a sub-matrix of V_bdiag. The Rcout commands can be used to print intermediate results or complete objects in the console when conducting the function. It is very useful to make sure you know what is happening between the lines. The output of cpp_bdiag() is:
# n_slices: 3 # n_rows: 2 # n_cols: 2 # A: # 0 0 # 0 0 # # V_bdiag: # 0 0 0 0 0 0 # 0 0 0 0 0 0 # 0 0 0 0 0 0 # 0 0 0 0 0 0 # 0 0 0 0 0 0 # 0 0 0 0 0 0 # # [,1] [,2] [,3] [,4] [,5] [,6] # [1,] 3.000000 2.078461 0.000000 0.000000 0.000000 0.000000 # [2,] 2.078461 4.000000 0.000000 0.000000 0.000000 0.000000 # [3,] 0.000000 0.000000 3.000000 2.078461 0.000000 0.000000 # [4,] 0.000000 0.000000 2.078461 4.000000 0.000000 0.000000 # [5,] 0.000000 0.000000 0.000000 0.000000 3.000000 2.078461 # [6,] 0.000000 0.000000 0.000000 0.000000 2.078461 4.000000
Function cpp_bdiag() returns a dense matrix. If you work with large matrices, using arma::sp_mat instead of arma::mat can be useful, it is for sparse matrices. Try it!
Example 3: Some More Function Examples for Matrices and Vectors
We want to make some additional examples to see how the functions and classes in RcppArmadillo work.
cppFunction(depends = "RcppArmadillo", ' int cpp_example3() { // Form a 2-by-2 matrix of 0s and fill some values arma::mat Mat1 = arma::mat(2, 2); Mat1(0, 1) = 2; Mat1(1, 1) = 4; Rcout << "Mat1:\\n" << Mat1 << "" << std::endl; Rcout << "The diagonal of Mat1:\\n" << Mat1.diag() << "" << std::endl; Rcout << "First column of Mat1:\\n" << Mat1.col(0) << "" << std::endl; // Form a 2-vector of 0s arma::vec x = arma::vec(2); Rcout << "x:\\n" << x << "" << std::endl; // Form a 2-vector of 1s arma::vec y = arma::ones<arma::vec>(2); Rcout << "y:\\n" << y << "" << std::endl; Rcout << "Return the unique elements of y: " << arma::unique(y) << "" << std::endl; // Form a 2-vector of uniform random values arma::vec z = arma::vec(2, arma::fill::randu); Rcout << "z:\\n" << z << "" << std::endl; return 0; } ' )
Function cpp_example3() returns the results.
set.seed(3) cpp_example3() # Mat1: # 0 2.0000 # 0 4.0000 # # The diagonal of Mat1: # 0 # 4.0000 # # First column of Mat1: # 0 # 0 # # x: # 0 # 0 # # y: # 1.0000 # 1.0000 # # Return the unique elements of y: 1.0000 # # z: # 0.1680 # 0.8075 # # [1] 0
We see different versions of initializing matrices or vectors, e.g. with random numbers. Again, the Rcout command comes in handy to show us what the return of the single functions is.
Video & Further Resources
A list of the RcppArmadillo functions and classes is given here. For general information on Rcpp, we also recommend you to take a look at the official webpage.
Furthermore, take a look at the following tutorial by Dirk Eddelbuettel held at the virtual useR! 2020 conference, which was published on the R Consortium Youtube channel.
We have some more blog posts which could be of interest to you.
- sugar Functions in Rcpp Package in R (Example)
- Construct 3D Array in Rcpp Package in R (2 Examples)
- Speed Up Loop Using Rcpp Package in R
- R programming Overview
We presented you some examples of R package RcppArmadillo which integrates C++ library Armadillo into R. If you have questions or want to comment on the port, let us know below.
This page was created in collaboration with Anna-Lena Wölwer. Have a look at Anna-Lena’s author page to get more information about her academic background and the other articles she has written for Statistics Globe.
4 Comments. Leave new
Dear Joachim ,
I would like to express my special thanks for all examples and tutorials you provided.
I am writing to ask your help for solving my issue occuring when I am running my code given below.
When I want to run my function called “my.cpp” in RStudio using command Rcpp::sourceCpp(‘E:/cpp files/my.cpp’) it appears a new box shown “R Session Aborted R encountered a fatal error. The Session was terminated. Start New Sesssion.”
My code is given below.
Could you please help me to sort out this issue ?
#####################################################################
#include
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat my(int p,int q, arma::mat Y, arma::vec mu, arma::mat sigma, arma::mat lambda) {
arma::vec x0(q);
arma::vec m(q);
arma::mat Omega = sigma + lambda*lambda.t();
arma::mat Omegainv = inv(Omega);
arma::mat I = I.eye( p, p );
int j=0;
arma::mat Delta = I – lambda.t()*Omegainv*lambda;
m = lambda.t()*Omegainv*(Y.col(j) – mu);
x0 = arma::randn(q);
x0 = m + arma::chol(Delta)*x0;
return(x0);
}
#####################################################################
n<-2;
p<-2;
q<-2;
Y<-matrix( c(rnorm(n), rnorm(n) ), nrow=2, ncol = n);
M<-c(0,0);
S<-matrix(c(1,.5,.5,1), nrow=p, ncol = p);
L<-matrix(c(0,0,2,2), nrow=p, ncol = q);
my(n, p, q, Y, M, S, L)
Hello Mahdi,
I was on vacation. That’s why I couldn’t respond sooner. Do you still need help?
Best,
Cansu
Hi Cansu,
Thanks for your contact.
Yes, I am waiting for your response.
I appreciate if you could help me for sorting out my problem.
Best, Mahdi.
Hello Mahdi,
I have just checked your code. Unfortunately, I am not familiar with this content. I advise you to post it on our Facebook discussion group. Someone there could help.
Best,
Cansu