Armadillo

• Developed by Dr. Conrad Sanderson and Dr Ryan Curtin

• Template based linear algebra library with high level syntax (like R or Matlab)

• Heavy lifting is (mostly) handled by LAPACK (i.e. benefits from OpenBLAS)

• Supports vectors, matrices, and cubes in dense or sparse format

• Some builtin expression optimization via template meta-programming

• Header only or shared library versions available

Basic types

Armadillo has 4 basic (dense) templated types:

arma::Col<type>,
arma::Row<type>, 
arma::Mat<type>, 
arma::Cube<type>

These types can be specialized using one of the following data types:

float, double,
std::complex<float>, std::complex<double>, 
short, int, long, 
unsigned short, unsigned int, unsigned long

typedef Shortcuts

For convenience the following typedefs are defined:

• Vectors:

arma::vec     = arma::colvec     =  arma::Col<double>
arma::dvec    = arma::dcolvec    =  arma::Col<double>
arma::fvec    = arma::fcolvec    =  arma::Col<float>
arma::cx_vec  = arma::cx_colvec  =  arma::Col<cx_double>
arma::cx_dvec = arma::cx_dcolvec =  arma::Col<cx_double>
arma::cx_fvec = arma::cx_fcolvec =  arma::Col<cx_float>
arma::uvec    = arma::ucolvec    =  arma::Col<uword>
arma::ivec    = arma::icolvec    =  arma::Col<sword>

• Matrices

arma::mat     = arma::Mat<double>
arma::dmat    = arma::Mat<double>
arma::fmat    = arma::Mat<float>
arma::cx_mat  = arma::Mat<cx_double>
arma::cx_dmat = arma::Mat<cx_double>
arma::cx_fmat = arma::Mat<cx_float>
arma::umat    = arma::Mat<uword>
arma::imat    = arma::Mat<sword>

RcppArmadillo

• Written and maintained by Dirk Eddelbuettel, Romain Francois, Doug Bates and Binxiang Ni

• Provides the header only version of Armadillo along with additional wrappers

  • Wrappers provide easy conversion between Rcpp types and Armadillo types
  • Enables use of Rcpp attributes and related tools

• Include the following in your C++ code:

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

Example Program

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

// [[Rcpp::export]]
arma::mat test_randu(int n, int m) {
  arma::mat A = arma::randu<arma::mat>(n,m);
  return A;
}

test_randu(4,5)

          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.1667928 0.6575326 0.5676687 0.6074325 0.2540894
[2,] 0.3053677 0.6481933 0.9097105 0.7308671 0.9109742
[3,] 0.3001430 0.7471385 0.6986718 0.9681653 0.3233815
[4,] 0.8806327 0.5012938 0.3641878 0.6065834 0.2796219

test_randu(3,1)

          [,1]
[1,] 0.9545359
[2,] 0.6865738
[3,] 0.7338450

arma class attributes

Element access

For an arma::vec v,

For an arma::mat m,

Element access - Cubes

For an arma::cube c,

Data Organization

fastLm example

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

// [[Rcpp::export]]
Rcpp::List fastLm(const arma::mat& X, const arma::colvec& y) {
    int n = X.n_rows, k = X.n_cols;
        
    arma::colvec coef = arma::solve(X, y);    // fit model y ~ X
    arma::colvec res  = y - X*coef;           // residuals

    // std.errors of coefficients
    double s2 = std::inner_product(res.begin(), res.end(), res.begin(), 0.0)/(n - k);
                                                        
    arma::colvec std_err = arma::sqrt(s2 * arma::diagvec(arma::pinv(arma::trans(X)*X)));

    return Rcpp::List::create(
      Rcpp::Named("coefficients") = coef,
      Rcpp::Named("stderr")       = std_err,
      Rcpp::Named("df.residual")  = n - k
    );
}

Multivariate Normal Distribution - Review

For an \(n\)-dimension multivate normal distribution with covariance \(\boldsymbol{\Sigma}\) can be written as

\[ \underset{n \times 1}{\boldsymbol{Y}} \sim N(\underset{n \times 1}{\boldsymbol{\mu}}, \; \underset{n \times n}{\boldsymbol{\Sigma}}) \]

where \(\big\{\boldsymbol{\Sigma}\big\}_{ij} = \rho_{ij} \sigma_i \sigma_j\)

\[ \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix} \sim N\left( \begin{pmatrix} \mu_1\\ \mu_2\\ \vdots\\ \mu_n \end{pmatrix}, \, \begin{pmatrix} \rho_{11}\sigma_1\sigma_1 & \rho_{12}\sigma_1\sigma_2 & \cdots & \rho_{1n}\sigma_1\sigma_n \\ \rho_{21}\sigma_2\sigma_1 & \rho_{22}\sigma_2\sigma_2 & \cdots & \rho_{2n}\sigma_2\sigma_n\\ \vdots & \vdots & \ddots & \vdots \\ \rho_{n1}\sigma_n\sigma_1 & \rho_{n2}\sigma_n\sigma_2 & \cdots & \rho_{nn}\sigma_n\sigma_n \\ \end{pmatrix} \right) \]

Multivariate Normal Distribution - Sampling

To generate draws from an \(n\)-dimensional multivate normal with mean \(\underset{n \times 1}{\boldsymbol{\mu}}\) and covariance matrix \(\underset{n \times n}{\boldsymbol{\Sigma}}\),

• Find a matrix \(\underset{n \times n}{\boldsymbol{L}}\) such that \(\boldsymbol{\Sigma} = \boldsymbol{L}\,\boldsymbol{L}^t\)

  • most often we use \(\boldsymbol{L} = \text{Chol}(\boldsymbol{\Sigma})\) where \(\boldsymbol{L}\) is a lower triangular matrix.

• Draw \(n\) iid unit normals, \(N(0,1)\), as \(\underset{n \times 1}{\boldsymbol{z}}\)

• Obtain multivariate normal draws using \[ \underset{n \times 1}{\boldsymbol{y}} = \underset{n \times 1}{\boldsymbol{\mu}} + \underset{n \times n}{\boldsymbol{L}} \, \underset{n \times 1}{\boldsymbol{z}} \]

Multivariate Normal Distribution - Density

For the \(n\) dimensional multivate normal given on the last slide, its density is given by

\[ f\big(\boldsymbol{Y} | \boldsymbol{\mu}, \boldsymbol{\Sigma}\big) = (2\pi)^{-n/2} \; \det(\boldsymbol{\Sigma})^{-1/2} \; \exp{\left(-\frac{1}{2} \underset{1 \times n}{(\boldsymbol{Y}-\boldsymbol{\mu})'} \underset{n \times n}{\boldsymbol{\Sigma}^{-1}} \underset{n \times 1}{(\boldsymbol{Y}-\boldsymbol{\mu})}\right)} \]

and its log density is given by

\[ \log f\big(\boldsymbol{Y} | \boldsymbol{\mu}, \boldsymbol{\Sigma}\big) = -\frac{n}{2} \log 2\pi - \frac{1}{2} \log \det(\boldsymbol{\Sigma}) - \frac{1}{2} \underset{1 \times n}{(\boldsymbol{Y}-\boldsymbol{\mu})'} \underset{n \times n}{\boldsymbol{\Sigma}^{-1}} \underset{n \times 1}{(\boldsymbol{Y}-\boldsymbol{\mu})} \]