Skip to content

Thought on crossprod

7 messages · Jonathan Rougier, Martin Maechler, Thomas Lumley +1 more

Original
Jonathan Rougier
# 
Hi everyone,

I do a lot of work with large variance matrices, and I like to use
"crossprod" for speed and to keep everything symmetric, i.e. I often
compute "crossprod(Q %*% t(A))" for "A %*% Sigma %*% t(A)", where
"Sigma" decomposes as "t(Q) %*% Q".  I notice in the code that
"crossprod", current definition
function (x, y = x) 
.Internal(crossprod(x, y))

does not recognise the case where y is the same as x as special.  I'd
like to extend crossprod to something like
function(x, y = NULL)
if (is.null(y))
  .Internal(crossprod.symm(x))
else
  .Internal(crossprod(x, y))

Before I start mesing up my R-1.4.1 installation, does anyone have any
thoughts on this?

Cheers, Jonathan.
Martin Maechler
# 
JonR> Hi everyone, I do a lot of work with large variance
    JonR> matrices, and I like to use "crossprod" for speed and
    JonR> to keep everything symmetric, i.e. I often compute
    JonR> "crossprod(Q %*% t(A))" for "A %*% Sigma %*% t(A)",
    JonR> where "Sigma" decomposes as "t(Q) %*% Q".  I notice in
    JonR> the code that "crossprod", current definition

    >> crossprod
    JonR> function (x, y = x) 
    JonR> .Internal(crossprod(x, y))

    JonR> does not recognise the case where y is the same as x as special.  I'd
    JonR> like to extend crossprod to something like

    >> crossprod
    JonR> function(x, y = NULL)
    JonR> if (is.null(y))
    JonR>    .Internal(crossprod.symm(x))
    JonR> else
    JonR>    .Internal(crossprod(x, y))

    JonR> Before I start mesing up my R-1.4.1 installation, does anyone have any
    JonR> thoughts on this?

Just for clarification:
The main point of your suggestion is to pass only half as much
data down to C, right?

Martin
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Thomas Lumley
#
On Fri, 15 Mar 2002, Martin Maechler wrote:

            

      
      
      
    

        
If this can be done with the BLAS functions then halving the necessary
memory bandwidth will also help -- that's one of the limiting factors in
ATLAS, for example.

	-thomas

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

  
  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

      
On Fri, 15 Mar 2002, Thomas Lumley wrote:

        

            

      
      
      
    

        
Looks like dsyrk can do this: despite the definition of *sy* BLAS routines
on my quick reference card it does not need a symmetric matrix, which is
why I overlooked it in earlier correspondence with Jonathan.  As from
1.5.0 dsyrk will be available in R.

  
    

      
      
    

  


  


      

    

    
2 days later

    

      
      

        


  

        

      Jonathan Rougier
    

    

      
      #
    

  


  

      
ripley@stats.ox.ac.uk 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.

Cheers, Jonathan.

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

      
On Mon, 18 Mar 2002, Jonathan Rougier wrote:

        
[...]

            

      
      
      
    

        
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).

Brian

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Brian Ripley
    

    

      
      #
    

  


  

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

        

            

      
      
      
    

        
Some data (on an old Sun):

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

without extra BLAS

            

      
      
      
    

        
[1] 7.21 0.02 7.31 0.00 0.00

with ATLAS

            

      
      
      
    

        
[1] 1.01 0.00 1.02 0.00 0.00

after adding the symmetric case

            

      
      
      
    

        
[1] 0.74 0.00 0.75 0.00 0.00

A <- matrix(rnorm(500*500), 5000)
without extra BLAS

            

      
      
      
    

        
[1] 0.83 0.00 0.85 0.00 0.00

with ATLAS

            

      
      
      
    

        
[1] 0.15 0.00 0.17 0.00 0.00

using symmetry

            

      
      
      
    

        
[1] 0.15 0.00 0.16 0.00 0.00

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