Compute Classical Principal Components via SVD or Eigen

Compute classical principal components (PC) via SVD (svd or eigenvalue decomposition (eigen) with non-trivial rank determination.

classPC(x, scale = FALSE, center = TRUE, signflip = TRUE,
        via.svd = n > p, scores = FALSE)

Arguments

x

a numeric matrix.

scale

logical indicating if the matrix should be scaled; it is mean centered in any case (via scale(*, scale=scale)c

center

logical or numeric vector for “centering” the matrix.

signflip

logical indicating if the sign(.) of the loadings should be determined should flipped such that the absolutely largest value is always positive.

via.svd

logical indicating if the computation is via SVD or Eigen decomposition; the latter makes sense typically only for n <= p.

scores

logical indicating

Author

Valentin Todorov; efficiency tweaks by Martin Maechler

Value

a list with components

rank

the (numerical) matrix rank of x; an integer number, say \(k\), from 0:min(dim(x)). In the \(n > p\) case, it is rankMM(x).

eigenvalues

the \(k\) eigenvalues, in the \(n > p\) case, proportional to the variances.

loadings

the loadings, a \(p \times k\) matrix.

scores

if the scores argument was true, the \(n \times k\) matrix of scores, where \(k\) is the rank above.

center

a numeric \(p\)-vector of means, unless the center argument was false.

scale

if the scale argument was not false, the scale used, a \(p\)-vector.

See also

In spirit very similar to R's standard prcomp and princomp, one of the main differences being how the rank is determined via a non-trivial tolerance.

Examples

set.seed(17)
x <- matrix(rnorm(120), 10, 12) # n < p {the unusual case}
pcx  <- classPC(x)
(k <- pcx$rank) # = 9  [after centering!]
#> [1] 9
pc2  <- classPC(x, scores=TRUE)
pcS  <- classPC(x, via.svd=TRUE)
all.equal(pcx, pcS, tol = 1e-8)
#> [1] TRUE
## TRUE: eigen() & svd() based PC are close here
pc0 <- classPC(x, center=FALSE, scale=TRUE)
pc0$rank # = 10  here *no* centering (as E[.] = 0)
#> [1] 10

## Loadings are orthnormal:
zapsmall( crossprod( pcx$loadings ) )
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#>  [1,]    1    0    0    0    0    0    0    0    0
#>  [2,]    0    1    0    0    0    0    0    0    0
#>  [3,]    0    0    1    0    0    0    0    0    0
#>  [4,]    0    0    0    1    0    0    0    0    0
#>  [5,]    0    0    0    0    1    0    0    0    0
#>  [6,]    0    0    0    0    0    1    0    0    0
#>  [7,]    0    0    0    0    0    0    1    0    0
#>  [8,]    0    0    0    0    0    0    0    1    0
#>  [9,]    0    0    0    0    0    0    0    0    1

## PC Scores are roughly orthogonal:
S.S <- crossprod(pc2$scores)
print.table(signif(zapsmall(S.S), 3), zero.print=".")
#>       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]  
#>  [1,] 35.700      .      .      .      .      .      .      .      .
#>  [2,]      . 21.600      .      .      .      .      .      .      .
#>  [3,]      .      . 17.300      .      .      .      .      .      .
#>  [4,]      .      .      . 14.300      .      .      .      .      .
#>  [5,]      .      .      .      .  6.120      .      .      .      .
#>  [6,]      .      .      .      .      .  5.450      .      .      .
#>  [7,]      .      .      .      .      .      .  2.300      .      .
#>  [8,]      .      .      .      .      .      .      .  1.390      .
#>  [9,]      .      .      .      .      .      .      .      .  0.493
stopifnot(all.equal(pcx$eigenvalues, diag(S.S)/k))

## the usual n > p case :
pc.x <- classPC(t(x))
pc.x$rank # = 10, full rank in the n > p case
#> [1] 10

cpc1 <- classPC(cbind(1:3)) # 1-D matrix
stopifnot(cpc1$rank == 1,
          all.equal(cpc1$eigenvalues, 1),
          all.equal(cpc1$loadings, 1))