The singular value decomposition is X=UDV'
.
If you want to compute X'X
,
that would be (UDV')'(UDV')
, i.e.,
VDU'UDV'
, i.e., V D^2 V'
(U
is orthogonal and D
diagonal).
f <- function(x) {
s <- svd(x)
v <- s$v
d <- diag(s$d) # It is a vector: transform it to a diagonal matrix
v %*% d^2 %*% t(v)
}
x <- matrix( rnorm(200), nc=4 )
stopifnot( all( abs( f(x) - t(x) %*% x ) < 1e-12 ) )
To have the variance matrix, you still need to subtract the mean of all the columns, and divide by the number of observations.
stopifnot( all( abs(
var(x) -
f(scale(x, scale=FALSE)) / (nrow(x)-1)
) < 1e-12
) )