class: center, middle, inverse, title-slide # Armadillo ### Colin Rundel ### 2019-03-19 --- exclude: true ```r knitr::opts_chunk$set(cache=TRUE) ``` --- ## Armadillo <br/> <img src="imgs/arma.png" width="66%" style="display: block; margin: auto;" /> <br/> * 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: .center-column[ ```c++ arma::Col<type>, arma::Row<type>, arma::Mat<type>, arma::Cube<type> ``` ] These types can be specialized using one of the following data types: .center-column[ ```c++ float, double, std::complex<float>, std::complex<double>, short, int, long, unsigned short, unsigned int, unsigned long ``` ] --- ## typedef Shortcuts For convenience the following typedefs have been defined: * Vectors: .center-column[ ```c++ 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 .center-column[ ```c++ 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 <br/> * Requirements - include the following in your C++ code .center-column[ ```c++ // [[Rcpp::depends(RcppArmadillo)]] #include <RcppArmadillo.h> ``` ] --- ## Example Program ```cpp // [[Rcpp::depends(RcppArmadillo)]] #include <RcppArmadillo.h> // [[Rcpp::export]] void test_randu() { arma::mat A = arma::randu<arma::mat>(4,5); Rcpp::Rcout << A << "\n"; } ``` ```r test_randu() ``` ``` ## 0.2784 0.1081 0.1102 0.9330 0.0741 ## 0.0976 0.4692 0.4514 0.3798 0.4071 ## 0.2892 0.7233 0.9054 0.2519 0.6407 ## 0.2051 0.4325 0.4188 0.2002 0.7786 ``` --- ## Class attributes | Attribute | Description | |----------------|---------------------------------------------------------------------------| | `.n_rows` | number of rows; present in Mat, Col, Row, Cube, field and SpMat | | `.n_cols` | number of columns; present in Mat, Col, Row, Cube, field and SpMat | | `.n_elem` | total number of elements; present in Mat, Col, Row, Cube, field and SpMat | | `.n_slices` | number of slices; present in Cube and field | --- ```cpp // [[Rcpp::depends(RcppArmadillo)]] #include <RcppArmadillo.h> // [[Rcpp::export]] void test_attr(arma::mat m) { Rcpp::Rcout << "m.n_rows = " << m.n_rows << "\n"; Rcpp::Rcout << "m.n_cols = " << m.n_cols << "\n"; Rcpp::Rcout << "m.n_elem = " << m.n_elem << "\n"; } ``` .columns[ .col[ ```r test_attr(matrix(0, 3, 3)) ``` ``` ## m.n_rows = 3 ## m.n_cols = 3 ## m.n_elem = 9 ``` ```r test_attr(matrix(1, 4, 5)) ``` ``` ## m.n_rows = 4 ## m.n_cols = 5 ## m.n_elem = 20 ``` ] .col[ ```r test_attr(1:10) ``` ``` ## Error in test_attr(1:10): Not a matrix. ``` ```r test_attr(as.matrix(1:10)) ``` ``` ## m.n_rows = 10 ## m.n_cols = 1 ## m.n_elem = 10 ``` ] ] --- ## Element access For an `arma::vec v`, | Call | Description | |-------------|-----------------------------------------------------| | `v(i)` | Access the `i`-th element with bounds checking | | `v.at(i)` | Access the `i`-th element without bounds checking | | `v[i]` | Access the `i`-th element without bounds checking | For an `arma::mat m`, | Call | type | Description | |-------------|---------|-----------------------------------------------------| | `m(i)` | `mat` | Access the `i`-th element, treating object as flat and in column major order | | `m(i,j)` | `mat` | Access the element in `i`-th row and `j`-th column with bounds checking | | `m.at(i,j)` | `mat` | Access the element in `i`-th row and `j`-th column without bounds checking | --- ## Element access - Cubes For an `arma::cube c`, | Call | type | Description | |-------------|---------|-----------------------------------------------------| | `c(i)` | `cube` | Access the i-th element, treating object as flat and in column major order | | `c(i,j,k)` | `cube` | Access the element in `i`-th row, `j`-th column, and `k`-th slice with bounds checking | | `c.at(i,j,k)`| `cube` | Access the element in `i`-th row, `j`-th column, and `k`-th slice without bounds checking | --- ## Data Organization <img src="imgs/colmajor.jpg" width="66%" style="display: block; margin: auto;" /> .footnote[[Image source](https://scc.ustc.edu.cn/zlsc/tc4600/intel/2017.0.098/mkl/common/mkl_userguide/GUID-ABCC618B-43C4-4DCD-ADA2-6F061B5116CD.html) ] --- ```cpp // [[Rcpp::depends(RcppArmadillo)]] #include <RcppArmadillo.h> // [[Rcpp::export]] void test_order(arma::mat m) { for(int i=0; i!=m.n_elem; ++i) { Rcpp::Rcout << m(i) << " "; } Rcpp::Rcout << "\n"; } ``` ```r m = matrix(1:9, 3, 3) c(m) ``` ``` ## [1] 1 2 3 4 5 6 7 8 9 ``` ```r test_order(m) ``` ``` ## 1 2 3 4 5 6 7 8 9 ``` ```r c(t(m)) ``` ``` ## [1] 1 4 7 2 5 8 3 6 9 ``` ```r test_order(t(m)) ``` ``` ## 1 4 7 2 5 8 3 6 9 ``` --- ## `fastLm` example ```cpp // [[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 ); } ``` --- .smaller[ ```r library(dplyr) n=1e5 d = tibble( x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n), x4 = rnorm(n), x5 = rnorm(n), ) %>% mutate( y = 3 + x1 - x2 + 2*x3 -2*x4 + 3*x5 - rnorm(n) ) res = bench::press( size = c(100, 1000, 10000, 100000), { d = d[seq_len(size),] X = model.matrix(y ~ ., d) y = as.matrix(d$y) bench::mark( lm(y~., data=d), lm.fit(X,y), .lm.fit(X,y), fastLm(X,y), check = FALSE ) } ) ``` ``` ## Running with: ## size ``` ``` ## 1 100 ``` ``` ## 2 1000 ``` ``` ## 3 10000 ``` ``` ## 4 100000 ``` ] --- <img src="Lec11_files/figure-html/unnamed-chunk-12-1.png" style="display: block; margin: auto;" /> --- class: middle # MVN Example