[R] multiplying a matrix by a vector

[Bert Gunter]

OK, just for fun, I added a couple of *obvious* (or should have been)
approaches  and got results that exactly agreed with Peter D's
comments: all O(n^2) results were essentially the same, and my matrix
multiplication O(n^3) solution was **far** worse. So I'm not sure
quite where Jeff's result came from.

microbenchmark( { t( y * t(x) ) }
                ,{ x %*% diag( y ) }
                ,{ sweep( x, 2, y, `*` ) }
                ,{vapply(seq_along(y),function(i)x[,i]*y[i],.1*seq_along(y))}
                ,{scale(x,center = FALSE,scale = 1/y)}
)

Unit: milliseconds

 expr        min         lq      mean
                                                         {     t(y *
t(x)) }   8.163510   9.880709  17.29899
                                                       {     x %*%
diag(y) } 637.639318 665.391992 686.12427
                                                 {     sweep(x, 2, y,
`*`) }   9.383453  10.497793  20.81000
 {     vapply(seq_along(y), function(i) x[, i] * y[i], 0.1 *
seq_along(y)) }   9.980690  11.748857  18.05533
                               {     scale(x, center = FALSE, scale =
1/y) }  10.427894  12.342690  19.57287

    median        uq       max neval
  11.72617  13.36316 100.00340   100
 679.38066 699.03779 827.17972   100
  13.11326  15.18546  95.41382   100
  13.47730  14.50548  98.80840   100
  14.08266  15.70222  99.52050   100


Note that scale() is basically a wrapper for sweep() with the api the
OP requested, so it would seem to be the "natural" choice. Also, my
matrix solution was an awful suggestion and should be tossed intp the
trash.

Cheers,
Bert
Bert Gunter

"The trouble with having an open mind is that people keep coming along
and sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )


On Fri, Nov 4, 2016 at 1:52 PM, Stefan Evert <stefanML at collocations.de> wrote:
>
>> On 4 Nov 2016, at 17:35, peter dalgaard <pdalgd at gmail.com> wrote:
>>
>> Notice though, that Bert loses (or _should_ lose) for larger values of N, since that method involves O(N^3) operations whereas the other two are O(N^2). I am a bit surprised that sweep() is so inefficient even at N=1000.
>
> I also got different results on Mac OS X (with R 3.3.1 and vecLib BLAS):
>
> Unit: milliseconds
>                         expr      min       lq     mean   median       uq       max neval
>          {     t(y * t(x)) } 18.46369 19.61954 25.77548 21.24242 22.72943  88.03469   100
>        {     x %*% diag(y) } 28.12786 29.87109 37.83814 31.59671 32.77839 101.68553   100
>  {     sweep(x, 2, y, `*`) } 11.19871 12.76442 21.48670 14.51618 15.70665 100.51444   100
>
> Note that sweep() is considerably _faster_ than the other two methods (which I found quite surprising).
>
> Of course, if you care about speed (and memory efficiency), you could also
>
>         library(wordspace)
>         scaleMargins(x, cols=y)
>
> Best,
> Stefan
>
> ______________________________________________
> 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.