[R] multiplying a matrix by a vector

[Bert Gunter]

Berend, et. al.:

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

Nor do I. But, as Peter noted, my matrix multiplication solution
"should" be worst in terms of efficiency, and it clearly is. The other
two *should* be "similar", and they are. So your results are
consistent with what one should expect, whereas, as Peter noted,
Jeff's are not.

But as we all know, *actual* efficiency in use can be tricky to
determine (asymptopia is always questionable), as it may depend on
state, problem details, exact algorithms in use, etc. Which is why the
standard first principle in code optimization is: don't.

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 11:27 AM, Berend Hasselman <bhh at xs4all.nl> wrote:
>
>> 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
>
> ______________________________________________
> 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.