r11 <- matrix(c(1,.4,.4,1),2,2) r12 <- matrix(c(.4,.2,.5,.3),byrow=T,ncol=2) r22 <- matrix(c(1,.2,.2,1),2,2) (r11.l <- eigen(r11)$val) (e11 <- eigen(r11)$vec) r11.m12 <- (1/sqrt(r11.l[1])) * (e11[,1] %*% t(e11[,1])) + (1/sqrt(r11.l[2])) * (e11[,2] %*% t(e11[,2])) solve(r11.m12 %*% r11.m12)