> I've been trying to find a way to improve (if
> possible) the efficiency of the code, especially when the number of
> components to calculate is higher than 100.
> pca.nipals <- function(X, ncomp, iter = 50, toler =
sqrt(.Machine$double.eps))
> # X...data matrix, ncomp...number of components,
> # iter...maximal number of iterations per component,
> # toler...precision tolerance for calculation of components
> {
>
> Xh <- scale(X, center = TRUE, scale = FALSE)
> nr <- 0
> T <- NULL; P <- NULL # matrix of scores and loadings
> for (h in 1:ncomp)
> {
> th <- Xh[, 1]
> ende <- FALSE
> while (!ende)
> {
> nr <- nr + 1
> ph <- t(crossprod(th, Xh) * as.vector(1 /
> crossprod(th)))
> ph <- ph * as.vector(1/sqrt(crossprod(ph)))
> thnew <- t(tcrossprod(t(ph), Xh) *
> as.vector(1/(crossprod(ph))))
> prec <- crossprod(th-thnew)
> th <- thnew
> if (prec <= (tol^2)) ende <- TRUE
> if (it <= nr) ende <- TRUE # didn't converge
> }
>
> Xh <- Xh - tcrossprod(th)
> T <- cbind(T, th); P <- cbind(P, ph)
> nr <- 0
> }
> list(T = T, P = P)
> }
First a disclaimer: I've just had a quick eyeball of your code; I haven't
run it, so this is speculative.
You've named the tolerance variable as 'toler', but in the line where you
come to check it, it is called tol:
if (prec <= (tol^2)) ende <- TRUE
I suspect that this means you have 'tol' as a global variable somewhere,
which may or may not be the number you want. Even if the correct variable
is being found, I suspect that you want prec <= tol, rather than the
square of it (you are probably wasting time on iterations calculating
excessive decimal places).
Also, it may be slightly faster to use a for loop instead of the while
loop, as so:
for(j in 1:iter)
{
#your calculations
if(prec <=tol) break
}
Regards,
Richie.
Mathematical Sciences Unit
HSL
------------------------------------------------------------------------
ATTENTION:
This message contains privileged and confidential inform...{{dropped:20}}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
