[R] crossprod vs %*% timing

[Robin Hankin]

  thanks everyone who replied.


Brian Ripley wrote:

>  t(a) %*% A %*% a is a quadratic form.  What varies `many many times'? If A
>  does not vary (often), you want to find B with B'B = A (e.g. via chol,

Yes indeed!   For me, "A" varies only very rarely, so it makes sense 
to use this fact.  However, keeping track of when it _does_ change is 
not straightforward, so I am now facing
a trade-off between clear understandable code, and fast code.

The time saving  appears to be negative for matrices under about 
50x50, but increases with the order
of the matrix.  With



  f1 <- function(a,X){ t(a) %*% X %*% a               }
  f2 <- function(a,X){ crossprod(t(crossprod(a,X)),a) }
  f3 <- function(a,X){  crossprod(a,X) %*% a           }
  f4 <- function(a,X.lower){jj <- crossprod(X.lower,a )
                             crossprod(jj,jj)}

timer <- function(n,r){
   a <- rnorm(n)
   X <- matrix(rnorm(n*n),n,n)
   X <- crossprod(X,X)
   X.lower <- t(chol(X))
   print(system.time( for(i in 1:r){ f1(a,X)}))
   print(system.time( for(i in 1:r){ f2(a,X)}))
   print(system.time( for(i in 1:r){ f3(a,X)}))
   print(system.time( for(i in 1:r){ f4(a,X)}))
}
print("n=50")
timer(50,100000)
print("n=500")
timer(500,1000)

I get


[1] "n=50"
[1] 8.62 0.07 8.49 0.00 0.00
[1] 3.40 0.01 3.59 0.00 0.00
[1] 1.90 0.02 1.79 0.00 0.00
[1] 2.15 0.01 2.09 0.00 0.00
[1] "n=500"
[1] 7.44 0.50 6.90 0.00 0.00
[1] 2.12 0.39 1.58 0.00 0.00
[1] 2.16 0.33 1.53 0.00 0.00
[1] 1.56 0.43 1.28 0.00 0.00


so it's not entirely clear that any single method dominates.



-- 
Robin Hankin
Uncertainty Analyst
Southampton Oceanography Centre
SO14 3ZH
tel +44(0)23-8059-7743
initialDOTsurname at soc.soton.ac.uk (edit in obvious way; spam precaution)