[R] multiplying a matrix by a vector

Berend Hasselman bhh at xs4all.nl
Fri Nov 4 19:27:55 CET 2016


> On 4 Nov 2016, at 16:41, Jeff Newmiller <jdnewmil at dcn.davis.ca.us> wrote:
> 
> Sara wins on memory use.
> 
> Rui wins on speed.
> 
> Bert wins on clarity.
> 
> library(microbenchmark)
> 
> N <- 1000
> x <- matrix( runif( N*N ), ncol=N )
> y <- seq.int( N )
> 
> microbenchmark( { t( y * t(x) ) }
>              , { x %*% diag( y ) }
>              , { sweep( x, 2, y, `*` ) }
>              )
> Unit: milliseconds
>                        expr       min        lq    median        uq      max neval
>         {     t(y * t(x)) }  6.659562  7.475414  7.871341  8.182623 47.01105 100
>       {     x %*% diag(y) }  9.859292 11.014021 11.281334 11.733825 48.79463 100
> {     sweep(x, 2, y, `*`) } 16.535938 17.682175 18.283572 18.712342 55.47159 100


I get different results with R3.2.2 on Mac OS X (using reference BLAS).


library(rbenchmark)

N <- 1000
x <- matrix( runif( N*N ), ncol=N )
y <- seq.int( N )

benchmark(  t( y * t(x) ), x %*% diag( y ), sweep( x, 2, y, `*` ) ,
           columns=c("test","elapsed","relative"), replications=10
         )

#                  test elapsed relative
# 3 sweep(x, 2, y, `*`)   0.132    1.000
# 1         t(y * t(x))   0.189    1.432
# 2       x %*% diag(y)   7.928   60.061

library(microbenchmark)
microbenchmark( { t( y * t(x) ) }
             , { x %*% diag( y ) }
             , { sweep( x, 2, y, `*` ) },
             times=10, unit="s"
             )

# Unit: seconds
#                         expr        min         lq       mean     median uq        max neval cld
#          {     t(y * t(x)) } 0.01099410 0.01132096 0.01597058 0.01332414  0.01430672 0.04447432    10  a
#        {     x %*% diag(y) } 0.76588887 0.78581717 0.79878255 0.79811626  0.80607877 0.82903945    10   b
#  {     sweep(x, 2, y, `*`) } 0.01177478 0.01246777 0.01409457 0.01314718   0.01600818 0.01802171    10  a


sweep appears to very good.

I don't quite understand why I get a very different ranking.

Berend



More information about the R-help mailing list