# [R] crossprod vs %*% timing

Robin Hankin rksh at soc.soton.ac.uk
Thu Oct 7 17:12:01 CEST 2004

```  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

 "n=50"
 8.62 0.07 8.49 0.00 0.00
 3.40 0.01 3.59 0.00 0.00
 1.90 0.02 1.79 0.00 0.00
 2.15 0.01 2.09 0.00 0.00
 "n=500"
 7.44 0.50 6.90 0.00 0.00
 2.12 0.39 1.58 0.00 0.00
 2.16 0.33 1.53 0.00 0.00
 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)

```