[Rd] Thought on crossprod

[Prof Brian D Ripley]

On Mon, 18 Mar 2002 ripley@stats.ox.ac.uk wrote:

> On Mon, 18 Mar 2002, Jonathan Rougier wrote:
>
> [...]
>
> > Sorry -- I've been away for the weekend.  My intention in the original
> > post was not to save on memory but simply to do half as many
> > inner-products (it seems inelegant not to exploit the symmetry!).  I had
> > not thought about losing the benefit of an optimised BLAS (I am not
> > using one myself at present).  I take it from this correspondence that,
> > for now, my best approach might be to write my own "crossprod.symm" in
> > C, but that this might actually slow me down in the future.
>
> You really should investigate optimized BLASs before worrying about this.
>
> I think the best approach is for us to modify crossprod to recognize the
> symmetric case in 1.5.0 and call the appropriate BLAS routine.
> Or, at least do so for real (not complex) inputs on IEEE-compliant
> machines (which we believe all current R platforms are).

Some data (on an old Sun):

A <- matrix(rnorm(500*500), 500)

without extra BLAS
> system.time(B <- crossprod(A, A))
[1] 7.21 0.02 7.31 0.00 0.00

with ATLAS
> system.time(B <- crossprod(A,A))
[1] 1.01 0.00 1.02 0.00 0.00

after adding the symmetric case
> system.time(B <- crossprod(A))
[1] 0.74 0.00 0.75 0.00 0.00

A <- matrix(rnorm(500*500), 5000)
without extra BLAS
> system.time(B <- crossprod(A,A))
[1] 0.83 0.00 0.85 0.00 0.00

with ATLAS
>  system.time(B <- crossprod(A,A))
[1] 0.15 0.00 0.17 0.00 0.00

using symmetry
>  system.time(B <- crossprod(A))
[1] 0.15 0.00 0.16 0.00 0.00

So the potential gains are minor compared to using an optimized BLAS.

-- 
Brian D. Ripley,                  ripley@stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._