A simple array class with specialized linear algebra routines

[This article was first published on Rcpp Gallery, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Currie, Durban and Eilers write:

Data with an array structure are common in statistics, and the design or
regression matrix for analysis of such data can often be written as a
Kronecker product. Factorial designs, contingency tables and smoothing of
data on multidimensional grids are three such general classes of data and
models. In such a setting, we develop an arithmetic of arrays which allows
us to define the expectation of the data array as a sequence of nested
matrix operations on a coefficient array. We show how this arithmetic leads
to low storage, high speed computation in the scoring algorithm of the
generalized linear model.

For example, they show that if a design matrix X has the Kronecker
structure X = kronecker(Xd, ..., X2, X1) with X a partial model matrix
with n rows and c columns, linear functions X%*%theta of X and
a coefficient vector theta can be efficiently computed based only on the
partial model matrices where the entries in the vector X%*%theta (with
nd*...*n2*n1 entries) are the same as the entries in the array with
dimension c(n1, n2, ..., nd) returned by RH(Xd, ... , RH(X2, RH(X1,
Theta is an array with dimensions c(c1, c2, ..., cd) containing theta
and RH(X, A) – the “rotated H-transform” – is an operation generalizing
transposed pre-multiplication t(X %*% A) of a matrix A by a matrix X to
the case of higher dimensional array-valued A.

The code below implements a simple array class for numeric arrays and the
rotated H-transform in RcppArmadillo and compares the performance to both
the naive straight forward matrix multiplication based on the full model
matrix and an R-implementation of RH():

// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;

Offset and Array classes based on code by Romain Francois copied from
http://comments.gmane.org/gmane.comp.lang.r.rcpp/5932 on 2014-01-07.

class Offset{
     IntegerVector dim ;

     Offset( IntegerVector dim ) : dim(dim) {}

     int operator()( IntegerVector ind ){
         int ret = ind[0] ;
         int offset = 1 ;
         for(int d=1; d < dim.size(); d++) {
             offset = offset * dim[d-1] ; 
             ret = ret + ind[d] * offset ;
         return ret ;
     } ;
     IntegerVector getDims() const {
         return(dim) ;

} ;

class Array : public NumericVector {
     // NumericVector value;
     Offset dims ;

    Array( SEXP x) : NumericVector(x), 
                     dims( (IntegerVector)((RObject)x).attr("dim") ) {}
    Array( NumericVector x,  Offset d ): NumericVector(x), 
                                                dims(d) {}
    Array( Dimension d ): NumericVector( d ), dims( d ) {}

    IntegerVector getDims() const {
    NumericVector getValue()  const {
     inline double& operator()( IntegerVector ind) {
         int vecind = dims(ind);
         NumericVector value = this->getValue();  
         return value(vecind);
     } ;
     // change dims without changing order of elements (!= aperm)
     void resize(IntegerVector newdim) {
         int n = std::accumulate((this->getDims()).begin(), (this->getDims()).end(), 1, 
         int nnew = std::accumulate(newdim.begin(), newdim.end(), 1, 
         if(n != nnew)  stop("old and new old dimensions don't match.");
         this->dims = Offset(newdim);
     } ;
} ;

namespace Rcpp {
    // wrap(): converter from Array to an R array
    template <> SEXP wrap(const Array& A) {
        IntegerVector dims = A.getDims();
        //Dimension dims = A.getDims();
        Vector x = A;
        x.attr( "dim" ) = wrap(dims);
        return x; 

// [[Rcpp::export]]
Array rotate(Array A){
    Re-dimension an array from dim to c(dim[-1], dim[1]).
    Example: for a 2*3*4 array, indices 1:24 are shuffled into  
    1 3 5 ... 21 23 2 4 6 ... 20 22 24  
    i.e., a sequence of length prod(dims[-1])=12 from 1 to prod(dims)=24 
    ("baseseq") repeated twice ("space" = dims[1]) and shifted by 1 each time.
    IntegerVector dims = A.getDims() ;
    int ndims = dims.size() ;
    int space = dims[0] ;
    int length = std::accumulate(dims.begin(),dims.end(), 1, 
                                 std::multiplies()) / space ;
    IntegerVector baseseq = (seq_len(length) - 1) * space ;
    NumericVector old = A.getValue() ; 
    NumericVector ret(space*length) ; 
    for(int r=0; r < space; r++) {
        for(int j=0; j < length; j++){
            ret[ r * length + j ] = old[ baseseq[j] + r ] ;
        } ;
    } ;
    IntegerVector newdim(ndims) ;
    for(int d=0; d < ndims; d++){
        newdim(d) = dims[d+1] ;
    } ; 
    newdim[ndims-1] = dims[0] ;  
   Array rA = Array(ret, Offset(newdim)) ;
   return rA;

// [[Rcpp::export]]
Array RH(const arma::mat& X, Array A){
    Rotated H-transform of Array A by matrix X.
    H-transform generalizes premultiplication of A by X to array-valued A.
    For A with dimensions (c1, c2, ..., cd) and X with dim=(n, c1),
    H(X, A) is array(X*Aflat, dim=c(n, c2, ..., cd)), where Aflat is 
    array(A, c(c1, c2*c3*..*cd).
    IntegerVector olddims = A.getDims() ;
    int n = A.getDims()[0] ;
    int d = std::accumulate(A.getDims().begin(), A.getDims().end(), 1, 
                                 std::multiplies()) / n ;
    arma::mat Amod((A.getValue()).begin(), n, d, false) ;
    arma::vec tmp = vectorise(X * Amod);
    IntegerVector newdims = clone(A.getDims());
    newdims[0] = X.n_rows; 
    Array ret  = rotate(Array(as(wrap(tmp)), Offset(newdims)));
    return ret ; 

Set up test case:

Let’s look at a 3-dimensional example where X = X3 %x% X2 %x% X1 and
each X is a B-spline basis over seq(0, 1, len=n):


n1 <- 30; n2 <- 40; n3 <- 50
c1 <- 5; c2 <- 10; c3 <- 15
n <- n1*n2*n3
c <- c1*c2*c3

X1 <- bs(seq(0, 1, len=n1), df=c1)
X2 <- bs(seq(0, 1, len=n2), df=c2)
X3 <- bs(seq(0, 1, len=n3), df=c3)
X <- X3 %x% X2 %x% X1

theta_vec <- runif(c)
Theta <- array(theta_vec, dim=c(c1, c2, c3))

RH_r <- function(X, A){
    ## H-transform:
    A_flat <- array(A, dim=c(dim(A)[1], prod(dim(A)[-1])))
    ret <- array(X %*% A_flat, dim=c(nrow(X), dim(A)[-1]))
    ## Rotate:
    aperm(ret, c(2:length(dim(A)), 1))

Note that X is fairly large, with 6 × 104 rows and 750 columns.

Check correctness:

    array(X%*%theta_vec, dim=c(n1, n2, n3)),
    # RH(X3, RH(X2, RH(X1, Theta))):
           list(X3, X2, X1),

[1] TRUE

    array(X%*%theta_vec, dim=c(n1, n2, n3)),
    # RH(X3, RH(X2, RH(X1, Theta))):
           list(X3, X2, X1),

[1] TRUE

Check performance:

    array(X%*%theta_vec, dim=c(n1, n2, n3)),
           list(X3, X2, X1),
           list(X3, X2, X1),
    replications =  100)[,c(1,3:4)]

                                         test elapsed relative
1 array(X %*% theta_vec, dim = c(n1, n2, n3))  27.252  198.920
3   Reduce(RH, list(X3, X2, X1), Theta, TRUE)   0.137    1.000
2 Reduce(RH_r, list(X3, X2, X1), Theta, TRUE)   0.325    2.372

Note: An alternative version with proper formula notation can be found here.

To leave a comment for the author, please follow the link and comment on their blog: Rcpp Gallery.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)