De-correlating multivariate normal data

A function for calculating the root \(\mathbf\Sigma^{\frac{1}{2}}\) of a covariance matrix, using a spectral decomposition:

makeRoot <- function(Sig){
  U <- eigen(Sig)
  RootSigma <- U$vectors %*% diag(sqrt(U$values)) %*% t(U$vectors)
  return(RootSigma)
}

An example covariance matrix \(\mathbf\Sigma\):

(mat <- matrix(c(1, 1.3, 1.3, 2), 2))

##      [,1] [,2]
## [1,]  1.0  1.3
## [2,]  1.3  2.0

cov2cor(mat)

##           [,1]      [,2]
## [1,] 1.0000000 0.9192388
## [2,] 0.9192388 1.0000000

And an example dataset \(\mathbf{X}\):

library("mvtnorm")
dat <- data.frame(rmvnorm(1000, sigma=mat))
head(dat)

##           X1          X2
## 1  1.2136601  1.49413055
## 2  0.4067333 -0.08562505
## 3 -0.3736126 -0.72329755
## 4 -0.2192473 -0.37922115
## 5 -0.8129728  0.13743406
## 6 -0.1860610 -0.67367826

cov(dat)

##          X1       X2
## X1 0.981142 1.256768
## X2 1.256768 1.905774

cor(dat)

##           X1        X2
## X1 1.0000000 0.9190802
## X2 0.9190802 1.0000000

De-correlate \(\mathbf{X}\) so as to get \(\mathbf{X}^* = \mathbf\Sigma^{-\frac{1}{2}} \mathbf{X}\):

DAT <- as.matrix(dat) %*% solve(makeRoot(cov(dat)))
cov(DAT)

##              [,1]         [,2]
## [1,] 1.000000e+00 2.269818e-16
## [2,] 2.269818e-16 1.000000e+00

cor(DAT)

##              [,1]         [,2]
## [1,] 1.000000e+00 2.269818e-16
## [2,] 2.269818e-16 1.000000e+00

Undo the de-correlation of \(\mathbf{X}\) (“re-correlate” \(\mathbf{X}^*\)):

dat2 <- DAT %*% makeRoot(cov(dat))
cov(dat2)

##          [,1]     [,2]
## [1,] 0.981142 1.256768
## [2,] 1.256768 1.905774

cor(dat2)

##           [,1]      [,2]
## [1,] 1.0000000 0.9190802
## [2,] 0.9190802 1.0000000