Sorry I missed your intent on the sapply.
Slick work on the vectorizing, but for the future reference it was
slightly buggy:
#######
A <- matrix( 1:12, nrow = 3 )
B <- matrix( 1:15, nrow = 3 )
# for loop
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
}
# fast but buggy
ans0 <- function( A, B ) {
nca <- ncol( A )
ncb <- ncol( B )
j.index <- rep( seq.int( nca ), times = ncb)
k.index <- rep( seq.int( nca ), each = ncb)
res <- array( A[ , j.index ] * B[ , k.index ]
, c( nrow( A ), nca, ncb )
)
res
}
# bugfixed
ans0b <- function( A, B ) {
nca <- ncol( A )
ncb <- ncol( B )
j.index <- rep( seq.int( nca ), times = ncb )
k.index <- rep( seq.int( ncb ), each = nca )
res <- array( A[ , j.index ] * B[ , k.index ]
, c( nrow( A ), nca, ncb )
)
res
}
library(microbenchmark)
microbenchmark( res2 <- ans2( A, B )
, res0b <- ans0b( A, B )
, res0 <- ans0( A, B )
)
#> Unit: microseconds
#> expr min lq mean median uq max
#> res2 <- ans2(A, B) 84.987 87.8185 270.2153 96.4315 99.4175 17531.77
#> res0b <- ans0b(A, B) 17.940 19.2055 126.8974 20.8800 22.2865 10616.36
#> res0 <- ans0(A, B) 18.041 19.1670 126.1183 20.5530 21.9545 10532.44
#> neval
#> 100
#> 100
#> 100
all( res2 == res0 )
#> [1] FALSE
all( res2 == res0b )
#> [1] TRUE
#' 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 11:43 AM, Jeff Newmiller <jdnew...@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
---------------------------------------------------------------------------
Jeff Newmiller The ..... ..... Go Live...
DCN:<jdnew...@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
______________________________________________
R-help@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.
