R devel - Oct 2002 - patch to mva:prcomp to use La.svd instead of svd (PR#2227)

Per the discussion about the problems with prcomp() when n << p,  which
boils down to a problem with svd() when n << p,
here is a patch to prcomp() which substitutes La.svd() instead of svd().

-Greg

(This is really a feature enhancement, but submitted to R-bugs to make sure
it doesn't get lost. )

*** R-1.6.0/src/library/mva/R/prcomp.R  Mon Aug 13 17:41:50 2001
--- R-1.6.0-GRW//src/library/mva/R/prcomp.R     Tue Oct 29 11:57:23 2002
***************
*** 1,18 ****
  prcomp <- function(x, retx = TRUE, center = TRUE, scale. = FALSE,
!                    tol = NULL) {
      x <- as.matrix(x)
      x <- scale(x, center = center, scale = scale.)
!     s <- svd(x, nu = 0)
      if (!is.null(tol)) {
          rank <- sum(s$d > (s$d[1]*tol))
          if (rank < ncol(x))
!             s$v <- s$v[, 1:rank, drop = FALSE]
      }
      s$d <- s$d / sqrt(max(1, nrow(x) - 1))
!     dimnames(s$v) <-
!         list(colnames(x), paste("PC", seq(len = ncol(s$v)), sep =
""))
!     r <- list(sdev = s$d, rotation = s$v)
!     if (retx) r$x <- x %*% s$v
      class(r) <- "prcomp"
      r
  }
--- 1,21 ----
  prcomp <- function (x, retx = TRUE, center = TRUE, scale. = FALSE,
!                     tol = NULL)
! {
      x <- as.matrix(x)
      x <- scale(x, center = center, scale = scale.)
!     s <- La.svd(x, nu = 0)
      if (!is.null(tol)) {
          rank <- sum(s$d > (s$d[1] * tol))
          if (rank < ncol(x))
!             s$vt <- s$vt[, 1:rank, drop = FALSE]
      }
      s$d <- s$d/sqrt(max(1, nrow(x) - 1))
!
!     dimnames(s$vt) <- list(paste("PC", seq(len = nrow(s$vt)), sep
= ""),
!                            colnames(x), )
!     r <- list(sdev = s$d, rotation = t(s$vt) )
!     if (retx)
!         r$x <- x %*% t(s$vt)
      class(r) <- "prcomp"
      r
  }


LEGAL NOTICE
Unless expressly stated otherwise, this message is confidential and may be
privileged. It is intended for the addressee(s) only. Access to this E-mail by
anyone else is unauthorized. If you are not an addressee, any disclosure or
copying of the contents of this E-mail or any action taken (or not taken) in
reliance on it is unauthorized and may be unlawful. If you are not an addressee,
please inform the sender immediately.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._