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 at
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 at 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 at 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