cf <- diag(3) cf[lower.tri(cf, diag = FALSE)] <- rnorm(3, sd=.2) c1 <- cf %*% t(cf) S <- diag(diag(cf)^2) A <- matrix(0,3,3) A[lower.tri(A, diag=FALSE)] <- cf[lower.tri(cf, diag = FALSE)] I <- diag(3) solve(I-A) %*% S %*% t(solve(I-A)) (A+I) %*% S %*% t(A+I) # more accurate