# [R] Avoiding loops to spare time and memory

David Brahm brahm at alum.mit.edu
Thu May 8 17:25:55 CEST 2003

```Dirk Enzmann <enzmann at kfn.uni-hannover.de> wrote:
> Is it possible to avoid the loop in the following function (or make the
> function otherwise more efficient)...

Prof Brian Ripley <ripley at stats.ox.ac.uk> replied:
> My guess is that most of the time is being spent in eigen(); if so you
> would have to use a different approach to gain much speed.

If M is an (p x q) matrix with q>>p, then at most the first p eigenvalues of
M'M are nonzero, and they can be found more efficiently using the left side of:
La.svd(M)\$d^2  =  eigen(crossprod(M))\$values[1:p]

I've re-written your function for the situation items >> cases (I leave the
opposite situation as an exercise).  I also threw in some sweep's, a zapsmall,
and a few other tricks.  A timing test shows it runs 57x faster for this
example (with items=500 and cases=10):
R> set.seed(1)
R> system.time(res1 <- RanEigen.orig(500,10,30))
[1] 229.80   0.07 230.38   0.00   0.00
R> set.seed(1)
R> system.time(res2 <- RanEigen(500,10,30))
[1] 4.00 0.02 4.07 0.00 0.00
R> all.equal(res1, res2)
[1] TRUE

and here's my function:

RanEigen <- function(items, cases, sample) {
if (items < cases) stop("I assume items >= cases")
EV <- matrix(0, sample, items)
for (i in 1:sample) {
X <- matrix(rnorm(cases*items), nrow=cases)
Xb <- sweep(X, 2, colMeans(X))
S <- crossprod(Xb)
Xc <- sweep(Xb, 2, 1/sqrt(diag(S)), "*")
EV[i, 1:cases] <- zapsmall(La.svd(Xc)\$d^2)  # The rest are zero
}
REigV <- cbind(1:items, colMeans(EV))
dimnames(REigV) <- list(rep('',items), c(' ','Eigenvalue'))
REigV
}

--
-- David Brahm (brahm at alum.mit.edu)

```