I am examining the output of the prcomp() function in R for PCA in light of the singular value decomposition equation:

$X = U \cdot \Sigma \cdot V^{T}$, where:

$X$: is the standardized original data matrix

$\Sigma$: is the eigenvalues matrix given by the output of prcomp() as $sdev but squared root

$V$: is the right singular eigenvector matrix (wrongly called loadings) in the output of prcomp() as $rotation

$U$: is the left singular eigenvector is NOT given by the output of prcomp(), unlike svd() as $u, however prcomp() provides the raw component scores matrix ($U\cdot \Sigma$) as $x which is practically useful for plotting purposes.

Now suppose I want to compute the $U$ matrix from the output of the prcomp() function using the mtcars built-in dataset.

To achieve this mathematically, both sides of the equation need to be multiplied by $\Sigma^{-1} \cdot V$ this will result in identity matrices for both $V$ and $\Sigma$ on the right side leaving us with:

$U = X \cdot \Sigma^{-1} \cdot {V}$

Note: I used svd() function for the sake of validating the results

My trial in R

# the svd() method
A <- as.matrix(mtcars[,c(1:7,10,11)], row = 9, byrow = T)
S <- scale(A)
svd.pca <- svd(S)
svd.pca$u # left singular eigenvector matrix U

# the prcomp() method
prcomp.pca <- prcomp(A, center = TRUE,scale. = TRUE)

# to get inverse of a matrix in R we use the solve() function
diag(prcomp.pca$sdev) %*% (solve(diag(prcomp.pca$sdev))) # identity matrix for Sigma
prcomp.pca$rotation %*% t(prcomp.pca$rotation) # identity matrix for V

# now compute the U matrix from prcomp() output per the equation above for U
prcomp.pca$u <- S %*% (solve(diag(prcomp.pca$sdev))) %*% t(prcomp.pca$rotation)

round(svd.pca$u, 3) == round(prcomp.pca$u, 3) # FALSE, what am I missing?

So the computed $U$ matrix eventually didn't match the $U$ matrix of the svd() function, what am I missing here to obtain the $U$ matrix from the output of prcomp()?

  • 1
    $\begingroup$ Why not just inspect the source? The code is in stats::prcomp.default and is only a few dozen lines long. Note that without details of how the algorithm works, you cannot expect to recover $U$ exactly because it is not uniquely defined: in particular, any of its columns may be arbitrarily negated. $\endgroup$
    – whuber
    Oct 24 '20 at 15:37
  • 2
    $\begingroup$ @whuber, you mean stats:::prcomp.default 3 semicolons instead of 2. I am inspecting it right now. $\endgroup$
    – doctorate
    Oct 24 '20 at 16:11
  • 1
    $\begingroup$ @whuber, so from source code the difference was in the sigma part, now U matrix could be recovered. Thanks for the hint! $\endgroup$
    – doctorate
    Oct 24 '20 at 16:51

This can be achieved by computing the $\Sigma$ first from the prcomp() function to make it numerically equivalent to that of the svd()and it turned out after inspecting the source code of prcomp() this can be achieved by the following:

sigma <- prcomp.pca$sdev * sqrt(max(1, nrow(S) - 1))
round(sigma, 2) == round(svd.pca$d, 2) # TRUE 

diag(sigma) %*% solve(diag(sigma)) # identity matrix of Sigma
prcomp.pca$rotation %*% t(prcomp.pca$rotation) # identity matrix of V

prcomp.pca$u <- S %*% prcomp.pca$rotation %*% solve(diag(sigma))
round(svd.pca$u, 3) == round(prcomp.pca$u, 3) # TRUE

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.