1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
## compute PCA of iris data set with svd() and sparsesvd()
library(sparsesvd)
library(Matrix)
data(iris)
M <- scale(as.matrix(iris[, 1:4]), scale=FALSE)
Ms <- Matrix(M) # not sparse, but a dMatrix
res1 <- svd(M)
res2 <- sparsesvd(Ms)
## check that eigenvalues are the same
print(res2$d, digits=3)
stopifnot(all.equal(res1$d, res2$d, tolerance=1e-12))
## these should be diagonal unit matrices
UtU <- abs(crossprod(res2$u, res1$u)) # diagonal entries may be 1 or -1
VtV <- abs(crossprod(res2$v, res1$v)) # (because sign of eigenvectors is arbitrary)
I1 <- diag(1, length(res1$d))
print(round(UtU, 12))
stopifnot(all.equal(UtU, I1, tolerance=1e-12))
print(round(VtV, 12))
stopifnot(all.equal(VtV, I1, tolerance=1e-12))
## check that SVD is reproducible
## (this is guaranteed by a deterministic RNG built into the SVDLIBC code)
for (i in 1:20) {
res <- sparsesvd(Ms)
if (!isTRUE(all.equal(res, res2))) stop("SVD not reproducible on iteration #", i)
}
|
