library('expm')
set.seed(1)
controls <- matrix(nrow=50000,ncol=8)
colnames(controls) <- letters[1:8]
controls[,'a'] <- rnorm(50000,sd=3)
controls[,'b'] <- rnorm(50000,sd=2)
controls[,'c'] <- controls[,'a'] + rnorm(50000,sd=2)
controls[,'d'] <- controls[,'b'] + rnorm(50000)
controls[,'e'] <- rnorm(50000,sd=0.5) - controls[,'a']
controls[,'f'] <- rnorm(50000,sd=0.5) - controls[,'b']
controls[,'g'] <- rnorm(50000)
controls[,'h'] <- rnorm(50000)
controls <- apply(controls,2,scale,center=T,scale=T)

cases <- matrix(nrow=20000,ncol=8)
colnames(cases) <- letters[1:8]
cases[,'a'] <- rnorm(20000,sd=3)
cases[,'b'] <- rnorm(20000,sd=2)
cases[,'c'] <- cases[,'a'] + rnorm(20000,sd=2)
cases[,'d'] <- cases[,'b'] + rnorm(20000) + cases[,'a']
cases[,'e'] <- rnorm(20000,sd=0.5) - cases[,'a']
cases[,'f'] <- rnorm(20000,sd=0.5) - cases[,'b'] + cases[,'a']
cases[,'g'] <- rnorm(20000) + cases[,'a']
cases[,'h'] <- cases[,'a']  - rnorm(20000) + cases[,'b']
cases <- apply(cases,2,scale,center=T,scale=T)

R <- cov(controls)
S <- cov(cases)
round(R,2)
round(S,2)

sqrtinvR <- sqrtm(solve(R))
round(sqrtinvR %*% R %*% sqrtinvR,2) # Identity matrix
round(cov(controls %*% sqrtinvR),2) # Identity matrix, columns no longer correlated
round(sqrtinvR %*% S %*% sqrtinvR,2) # Why is cov(a,b) (eg. cell [1,2]) not close to 0?

round(cov(cases %*% sqrtinvR),2) # Done this way, the covariance between a and b (eg. cell [1,2]) is still not close to 0
round(cor(cases %*% sqrtinvR),2) # As a correlation it would be about 0.34