# Create two (2x2) coviance matrices:
d <- 2 # The dimension
M <- 2 # The number of covariance matrices
Omega1 <- matrix(c(2, 0.5, 0.5, 2), nrow=d)
Omega2 <- matrix(c(1, -0.2, -0.2, 1), nrow=d)
# The decomposition with Omega1 as the first covariance matrix:
decomp1 <- diag_Omegas(Omega1, Omega2)
W <- matrix(decomp1[1:d^2], nrow=d, ncol=d) # Recover W
lambdas <- decomp1[(d^2 + 1):length(decomp1)] # Recover lambdas
tcrossprod(W) # = Omega1
W%*%tcrossprod(diag(lambdas), W) # = Omega2
# Reorder the covariance matrices in the decomposition so that now
# the first covariance matrix is Omega2:
decomp2 <- redecompose_Omegas(M=M, d=d, W=as.vector(W), lambdas=lambdas,
perm=2:1)
new_W <- matrix(decomp2[1:d^2], nrow=d, ncol=d) # Recover W
new_lambdas <- decomp2[(d^2 + 1):length(decomp2)] # Recover lambdas
tcrossprod(new_W) # = Omega2
new_W%*%tcrossprod(diag(new_lambdas), new_W) # = Omega1
Run the code above in your browser using DataLab