[R] A slightly unorthodox matrix product.

[Jeff Newmiller]

Sometimes a good old for loop performs best, even if it doesn't look sexy:

########
A <- matrix( 1:12, nrow = 3 )
B <- matrix( 1:15, nrow = 3 )

library(abind)
# Eric
ans1 <- function( a, b ) {
   xxx <- lapply( seq.int( nrow( A ) )
                , function( i ) {
                     A[ i, ] %o% B[ i, ]
                  }
                )
   yyy <- do.call( abind, c( xxx, list( along = 3 ) ) )
   zzz <- aperm( yyy, c( 3, 1, 2 ) )
   zzz
}
# Charles
ans1b <- function( a, b ) {
   xxx <- lapply( seq.int( nrow( A ) )
                , function( i ) {
                     A[ i, ] %o% B[ i, ]
                  }
                )
   yyy <- sapply( seq.int( nrow( a ) )
                , function( i ) a[ i, ] %o% b[ i, ]
                , simplify = "array"
                )
   zzz <- aperm( yyy, c( 3, 1, 2 ) )
   zzz
}
# Jeff #1
ans2 <- function( a, b ) {
   zzz <- array( rep( NA, nrow( a ) * ncol( a ) * ncol( b ) )
               , dim = c( nrow( a ), ncol( a ), ncol( b ) )
               )
   jseq <- seq.int( ncol( a ) )
   kseq <- seq.int( ncol( b ) )
   for ( i in seq.int( nrow( a ) ) ) {
     zzz[ i, jseq, kseq ] <- outer( a[ i, ], b[ i, ] )
   }
   zzz
}
# Jeff #2
ans3 <- function( a, b ) {
   idxs <- expand.grid( i = seq.int( nrow( a ) )
                      , j = seq.int( ncol( a ) )
                      , k = seq.int( ncol( b ) )
                      )
   ij <- as.matrix( idxs[ , c( "i", "j" ) ] )
   ik <- as.matrix( idxs[ , c( "i", "k" ) ] )
   array( a[ ij ] * b[ ik ]
        , dim = c( nrow( a ), ncol( a ), ncol( b ) )
        )
}

library(microbenchmark)

microbenchmark( res1 <- ans1( A, B )
               , res1b <- ans1b( A, B )
               , res2 <- ans2( A, B )
               , res3 <- ans3( A, B )
               )
#> Unit: microseconds
#>                  expr     min       lq      mean   median       uq
#>    res1 <- ans1(A, B) 660.489 688.3460 4199.5385 742.5505 805.1860
#>  res1b <- ans1b(A, B) 224.436 236.2250  427.4806 246.3240 269.6425
#>    res2 <- ans2(A, B)  91.538  96.9075  287.9596 102.0335 110.8825
#>    res3 <- ans3(A, B) 508.642 528.9700  860.6295 563.5470 619.5285
#>        max neval
#>  344769.27   100
#>   17062.63   100
#>   18212.11   100
#>   23041.89   100
all( res1 == res2 )
#> [1] TRUE
all( res1 == res1b )
#> [1] TRUE
all( res1 == res3 )
#> [1] TRUE
res3
#> , , 1
#>
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    4    7   10
#> [2,]    4   10   16   22
#> [3,]    9   18   27   36
#>
#> , , 2
#>
#>      [,1] [,2] [,3] [,4]
#> [1,]    4   16   28   40
#> [2,]   10   25   40   55
#> [3,]   18   36   54   72
#>
#> , , 3
#>
#>      [,1] [,2] [,3] [,4]
#> [1,]    7   28   49   70
#> [2,]   16   40   64   88
#> [3,]   27   54   81  108
#>
#> , , 4
#>
#>      [,1] [,2] [,3] [,4]
#> [1,]   10   40   70  100
#> [2,]   22   55   88  121
#> [3,]   36   72  108  144
#>
#> , , 5
#>
#>      [,1] [,2] [,3] [,4]
#> [1,]   13   52   91  130
#> [2,]   28   70  112  154
#> [3,]   45   90  135  180

#' Created on 2018-08-04 by the [reprex package](http://reprex.tidyverse.org) (v0.2.0).
########

On Sat, 4 Aug 2018, Berry, Charles wrote:

>
>
>> On Aug 4, 2018, at 10:01 AM, Eric Berger <ericjberger using gmail.com> wrote:
>>
>> Hi Rolf,
>> A few edits because (i) nrow(a) should be nrow(A) and (ii) you have
>> calculated C[j,k,i] = A[i,j]*B[i,k], (iii) minor style change on lapply.
>>
>> library(abind)
>> xxx <- lapply(1:nrow(A),function(i){A[i,]%o%B[i,]})
>> yyy <- do.call(abind,c(xxx,list(along=3)))
>
> Or use the simplify="array" gambit in sapply:
>
> yyy <- sapply(1:nrow(A), function(i) A[i,] %o% B[i,], simplify="array")
>
>> zzz <- aperm(yyy,c(3,1,2))
>>
>
> HTH,
>
> Chuck
>
> ______________________________________________
> R-help using r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

---------------------------------------------------------------------------
Jeff Newmiller                        The     .....       .....  Go Live...
DCN:<jdnewmil using dcn.davis.ca.us>        Basics: ##.#.       ##.#.  Live Go...
                                       Live:   OO#.. Dead: OO#..  Playing
Research Engineer (Solar/Batteries            O.O#.       #.O#.  with
/Software/Embedded Controllers)               .OO#.       .OO#.  rocks...1k