[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)
More information about the R-help
mailing list