A More Complicated Example
I will admit, the above example is pretty contrived. Lets consider a more complicated example that was my acctual motivation for creating this package. Lets say we have a posterior distribution over covariance matricies and we want to visualize the posterior distribution of the first and second eigenvectors of these matrices. Importantly, we care about visualizing the “principle directions of variation” not just what the algorithm tells us are the first and second eigenvectors (this distnction should be clear in a second).
I start by simulating a distributio over covariance matricies. Here I am going to do it in a bit of a contrived way to show my point. I am going to simulate covariace matricies from a set of fixed eigenvectors and just two eigenvalues one that is bigger than the other. Yet I am going to randomly select which of the two eigenvectors gets the larger eigenvalue. In other words I have a bunch of covariance matricies which have the same rotation, but different scaling” along that rotation. Finally I add a bit of Wishart noise to ensure that the eigenvectors are not exactly the same just very very close to the true values - this just makes the plots prettier.
# Create a distribution of matricies that have identical eigenvectors but there
# is randomness in the eigenvalues and a tiny bit of added noise from a wishart
const <- 1/sqrt(2)
V <- cbind(c(const,const), c(const,-const))
Sigma <- array(0, dim=c(2, 2, 100))
for (i in 1:dim(Sigma)[3]){
Sigma[,,i] <- V %*% diag(sample(c(1, 1.3))) %*% t(V) + rWishart(1, 50, 0.002*diag(2))[,,1]
}Just as a reference here is a plot of the true eigenvectors shared by all matricies in the sample
We could take the eigen decomposition of each covariance matrix and then just plot the distribution over the first and second eigenvectors but there is are two problems:
- The ordering of the eigenvector is defined by the ranking of the eigenvalues.
- Transform is only defined up to the sign of the eignevectors.
The first problem implies that the “direction” reflected in the first eigenvalue of one sample might acctually correspond to the second eigenvalue in the second sample etc… Therefore we will need some method of matching “directions” between samples (get it?… we are going to do this matching with the Hungarian algoritym)
The second problem implies that we will have to be smart about how we do that matching, we acctually might need to match an eigenvector \(v\) in one sample with the eigenvector \(-u\) in a second (i.e., we may need to match to the negative of vector).
So lets do it the naive way and see what happens.
# Decompose each sample
Sigma.vectors <- array(0, dim=dim(Sigma))
for (i in 1:dim(Sigma)[3]){
decomp <- eigen(Sigma[,,i])
Sigma.vectors[,,i] <- decomp$vectors
}
# Plot distribution of the first two Eigenvectors
dat <- rbind(t(Sigma.vectors[,1,]),
t(Sigma.vectors[,2,]))
col <- c(rep("1", dim(Sigma)[3]), rep("2", dim(Sigma)[3]))
plot_arrows(dat, col)
Hard to interpret right? Weirdly there are eigenvectors going each of 4 directions for some reason with their lables all mixed up and some +/- weirdness happening.
So as I alluded to before, we will solve this using the Hungarian algoritym. We are going to pick a reference sample (in this case the first sample) and match the eigenvectors of each sample against the eigenvectors of the reference. In this case our “cost matrix” will be the pairwise euclidean distance between eigenvectors from the reference and eigenvectors of the sample to be matched. Notably, to deal with the sign issues (problem 2 above) we will use a trick, every eigenvector in the sample will be doubled (the second copy is the negation of the first). The following code does all this.
# Here is our distance function - calculates the pairwise distance
# between two sets of vectors.
cost_distance <- function(V1, V2){
n1 <- ncol(V1)
n2 <- ncol(V2)
D <- as.matrix(dist(t(cbind(V1, V2))))
return(D[1:n1, (n1+1):(n1+n2)])
}
# Make some empty containers to store results
Sigma.vectors <- array(0, dim=dim(Sigma))
# Initialize n=1 case - this is our reference against which
# we will match the rest of the samples
decomp <- eigen(Sigma[,,1])
Sigma.vectors[,,1] <- decomp$vectors
# Now match up each sample with the reference
for (i in 2:dim(Sigma)[3]){
decomp <- eigen(Sigma[,,i])
Sigma.vectors[,,i] <- decomp$vectors
# Here we add a second (negated) copy of the eigenvectors in the sample
Sigma.vectors.expanded <- cbind(Sigma.vectors[,,i], -Sigma.vectors[,,i])
# Here we compute the "cost matrix" (note this is rectangular)
# because we have the duplicated (negated) vectors from the sample
cost <- cost_distance(Sigma.vectors[,,1],Sigma.vectors.expanded)
# This is the key line where we solve the matching problem
matching <- RcppHungarian::HungarianSolver(cost)
# This is the key line were we reorder the eigenvectors / select if we wanted
# the negative or positive versions
Sigma.vectors[,,i] <- Sigma.vectors.expanded[,matching$pairs[,2]]
}
# Now we check our work visually
dat <- rbind(t(Sigma.vectors[,1,]),
t(Sigma.vectors[,2,]))
col <- c(rep("1", dim(Sigma)[3]), rep("2", dim(Sigma)[3]))
plot_arrows(dat, col)
Looks much better right? I woudl say this is about as good as this is going to get. Note, it may still not look identical to the original (the true eigenvectors) because we don’t know the true sign of the eigenvectors but the overall information conveyed by our result is identical to the original.