[R] fast repeated matrix operations

[Paul Rathouz]

Hi -- I am wondering if there is a suite of R functions or an R package for the following types of purpose:

Suppose I had a nxpxp array (say, A) and I wanted to invert each pxp "layer" in the array (i.e., solve(A[i,,]) for i in 1:n).  Or, suppose I had a nxpxq array (A) and a nxqxr array (B) and I wanted to multiply the pxq matrix from each layer of the first array with the corresponding qxr layer of the second array, i.e., A[i,,]%*%B[i,,] for i in 1:n.

Is there a computationally efficient way to do this? I can of course use a for() loop.  It can also be done, I think, with the tensorA package.  But, a quick test on inverse ( inv.tensor() ) showed that a for() loop around solve() is faster than inv.tensor() from tensorA. But, i am concerned that a for() loop is not very efficient.

A long while back, I designed some R functions that called fortran code to do the looping.  Those functions are rudimentary, but I could develop them further.  A quick test suggested they were very much faster than a for() loop around solve(). 

Finally, I am not just interested in solve() and matrix multiplication, but other matrix operations as well.

Before going further, I am in search of existing solutions. -- pr

Paul Rathouz, PhD
Professor and Chair
Department of Biostatistics & Medical Informatics
University of Wisconsin School of Medicine and Public Health
K6/446 CSC, Box 4675
600 Highland Avenue
Madison, WI 53792-4675
608.263.1706