R Under development (unstable) (2026-01-25 r89330 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2026 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require("irlba")
Loading required package: irlba
Loading required package: Matrix
> 
> # prcomp convenience function
> x  <- matrix(rnorm(200), nrow=20)
> p1 <- prcomp_irlba(x, n=3)
> p2 <- prcomp(x, tol=0.7)
> if (!isTRUE(all.equal(p1$sdev[1:2], p2$sdev[1:2])))
+ {
+   stop("Failed basic prcomp test")
+ }
> 
> s <- summary(p1)
> 
> # scaling bug identified in issue #21
> normalize_signs <- function(X, Y) {
+   for (i in 1:ncol(X)) {
+     if (sign(X[1, i]) != sign(Y[1, i])) {
+       Y[, i] <- -Y[, i]
+     }
+   }
+   return(Y)
+ }
> 
> all.equal_pca <- function(X, Y) {
+   Y <- normalize_signs(X, Y)
+   return(all.equal(X, Y, check.attributes=F, tolerance=1e-4))
+ }
> 
> set.seed(1)
> X <- matrix(rnorm(2000), ncol=40)
> M <- 5 # number of PCA components
> centers <- colMeans(X)
> sds <- apply(X, 2, sd)
> rms <- apply(X, 2, function(x) sqrt(sum(x^2) / (length(x) - 1)))
> Xc <- sweep(X, 2, centers, `-`)
> Xs <- sweep(X, 2, sds, `/`)
> Xcs <- sweep(Xc, 2, sds, `/`)
> Xrms <- sweep(X, 2, rms, `/`)
> 
> # unscaled
> scaled <- FALSE
> centered <- FALSE
> pca <- prcomp(X, center=centered, scale.=scaled)
> sv <- svd(X)
> svir <- irlba(X, nv=M, nu=M)
> pcair <- prcomp_irlba(X, n=M, center=centered, scale.=scaled)
> Xpca <- predict(pca)[, 1:M]
> Xsvl <- sv$u[, 1:M] %*% diag(sv$d[1:M])
> Xsvr <- X %*% sv$v[, 1:M]
> Xsvirl <- svir$u %*% diag(svir$d)
> Xsvirr <- X %*% svir$v
> Xpcair <- predict(pcair)
> Xpcair2 <- X %*% pcair$rotation
> 
> if (! isTRUE(all.equal_pca(Xsvl, Xsvr)) &&
+      isTRUE(all.equal_pca(Xpca, Xsvl)) &&
+      isTRUE(all.equal_pca(Xsvirl, Xsvirr)) &&
+      isTRUE(all.equal_pca(Xpca, Xsvirl)) &&
+      isTRUE(all.equal_pca(Xpcair, Xpcair2)) &&
+      isTRUE(all.equal_pca(Xpca, Xpcair)) &&
+      isTRUE(all.equal_pca(Xpcair, Xsvirl)))
+ {
+   stop("failed unscaled, uncentered prcomp")
+ }
> 
> # scaled, uncentered
> scaled <- TRUE
> centered <- FALSE
> pca <- prcomp(X, center=centered, scale.=scaled)
> sv <- svd(Xrms)
> svir <- irlba(X, nv=M, nu=M, scale=rms)
> pcair <- prcomp_irlba(X, n=M, center=centered, scale.=scaled)
> 
> Xpca <- predict(pca)[, 1:M]
> Xsvl <- sv$u[, 1:M] %*% diag(sv$d[1:M])
> Xsvr <- Xrms %*% sv$v[, 1:M]
> Xsvirl <- svir$u %*% diag(svir$d)
> Xsvirr <- Xrms %*% svir$v
> Xpcair <- predict(pcair)
> Xpcair2 <- Xrms %*% pcair$rotation
> 
> if (! isTRUE(all.equal_pca(Xsvl, Xsvr)) &&
+       isTRUE(all.equal_pca(Xpca, Xsvl)) &&
+       isTRUE(all.equal_pca(Xsvirl, Xsvirr)) &&
+       isTRUE(all.equal_pca(Xpca, Xsvirl)) &&
+       isTRUE(all.equal_pca(Xpcair, Xpcair2)) &&
+       isTRUE(all.equal_pca(Xpca, Xpcair)) &&
+       isTRUE(all.equal_pca(Xpcair, Xsvirl)))
+ {
+   stop("failed scaled, uncentered prcomp")
+ }
> 
> 
> # issue #25 prcomp_irlba regression (error in scale. handling)
> set.seed(1)
> x <- matrix(rnorm(100), 10)
> p <- prcomp_irlba(x, 3, scale.=TRUE, fastpath=FALSE)
> p <- prcomp_irlba(x, 3, scale.=TRUE, fastpath=TRUE)
> 
> 
> # issue #32 (and also issue #25 redux) more checks for proper
> # variance proportion computation
> library(irlba)
> set.seed(1)
> x  <- matrix(rnorm(200), nrow=20)
> n  <- 3
> s1 <- summary(prcomp_irlba(x, n=n, center=TRUE, scale.=FALSE))
> s2 <- summary(prcomp(x, tol=0.7, center=TRUE, scale.=FALSE))
> if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) &&
+              all.equal(s1$importance, s2$importance[, 1:n])))
+ {
+   stop("center=TRUE scale.=FALSE prcomp variance computation")
+ }
> 
> s1 <- summary(prcomp_irlba(x, n=3, center=TRUE, scale.=TRUE))
> s2 <- summary(prcomp(x, tol=0.8, center=TRUE, scale.=TRUE))
> if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) &&
+              all.equal(s1$importance, s2$importance[, 1:n])))
+ {
+   stop("center=TRUE scale.=TRUE prcomp variance computation")
+ }
> 
> s1 <- summary(prcomp_irlba(x, n=3, center=FALSE, scale.=TRUE))
> s2 <- summary(prcomp(x, tol=0.8, center=FALSE, scale.=TRUE))
> if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) &&
+              all.equal(s1$importance, s2$importance[, 1:n])))
+ {
+   stop("center=FALSE scale.=TRUE prcomp variance computation")
+ }
> 
> s1 <- summary(prcomp_irlba(x, n=3, center=FALSE, scale.=FALSE))
> s2 <- summary(prcomp(x, tol=0.7, center=FALSE, scale.=FALSE))
> if (! isTRUE(all.equal(s1$sdev, s2$sdev[1:n]) &&
+              all.equal(s1$importance, s2$importance[, 1:n])))
+ {
+   stop("center=FALSE, scale.=FALSE prcomp variance computation")
+ }
> 
> proc.time()
   user  system elapsed 
   1.29    0.21    1.48