[R] A slightly unorthodox matrix product.

[Berry, Charles]

> On Aug 4, 2018, at 11:43 AM, Jeff Newmiller <jdnewmil using dcn.davis.ca.us> wrote:
> 
> Sometimes a good old for loop performs best, even if it doesn't look sexy:
> 
> 

Fair enough, but a vectorized solution beats them all (see below).

Also,

[SNIP]


> # Charles
> ans1b <- function( a, b )
> {

The lapply you put here was from Eric's solution:

>  xxx <- lapply( seq.int( nrow( A ) )
>               , function( i ) {
>                    A[ i, ] %o% B[ i, ]
>                 }


This is what I had in mind:

ans1b.corrected <- function( a, b ) {
  yyy <- sapply( seq.int( nrow( a ) )
                 , function( i ) a[ i, ] %o% b[ i, ]
                 , simplify = "array"
  )
  zzz <- aperm( yyy, c( 3, 1, 2 ) )
  zzz
}

On my system it is slower than a for loop but a lot faster than your benchmark showed with the superfluous code from Eric's solution.

For speed, a vectorized solution is faster than a for loop by a factor of 3 on my laptop:

ans0 <- function(A,B){
  nca <- ncol(A)
  ncb <- ncol(B)
  j.index <- rep(1:nca, times = ncb)
  k.index <- rep(1:nca, each = ncb)
  res <- array(A[, j.index] * B[, k.index], c(nrow(A), nca, ncb))
  res
  }


> microbenchmark(
+   res0 <- ans0(A, B),
+   res1 <- ans1(A, B),
+   res1b <- ans1b.corrected(A, B),
+   res2 <- ans2(A, B),
+   res3 <- ans3(A, B)
+ )
Unit: microseconds
                           expr     min       lq      mean   median       uq     max neval   cld
             res0 <- ans0(A, B)  13.281  18.4960  21.52723  19.9905  23.4750  61.556   100 a    
             res1 <- ans1(A, B) 353.121 369.8635 409.77788 381.5840 444.3290 701.256   100     e
 res1b <- ans1b.corrected(A, B)  82.816  89.4185 101.52321  95.4275 107.1700 217.357   100   c  
             res2 <- ans2(A, B)  49.674  54.4825  61.78278  58.7540  65.5265 172.625   100  b   
             res3 <- ans3(A, B) 317.772 342.4220 392.25065 360.4675 436.2125 602.346   100    d 
> 


FWIW, if there are dimnames on A and B, sapply( row names(A), ..., simplify="array") preserves them without further ado.

Chuck