In utils.jl, there is a mistake in the formula below when it reconstructs the covariance matrix after eigen decomposition:
|
return Symmetric(eigvect' * Diagonal(eigval) * eigvect) |
The correct order should be:
eigvect * Diagonal(eigval) * eigvect'
i.e., the first matrix is not transposed while the last is. To quickly check,
using LinearAlgebra;
A = Symmetric(rand(3,3));
eigval, eigvect = eigen(A);
norm(eigvect * Diagonal(eigval) * eigvect' - A) < 1e-12 # true
norm(eigvect' * Diagonal(eigval) * eigvect - A) < 1e-12 # false