Skip to content
Prev 31777 / 398506 Next

Avoiding loops to spare time and memory

Dirk Enzmann <enzmann at kfn.uni-hannover.de> wrote:
Prof Brian Ripley <ripley at stats.ox.ac.uk> replied:
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
}