Eigenvectors are only defined up to multiplication by a scalar, and, consequently, the sign of eigenvectors can flip arbitrarily between different calls to eigenvectors_sym even if the input matrices are “close” to one another. That’s of course not a problem in theory, but I was wondering how this might affect HMC in practice because it relies on the unconstrained space being smooth.
Suppose we have vector[n] f = U * z, where matrix [n, n] U = eigenvectors_sym(A) and z is another vector. Then if the sign of a column of U flips, the sign of the corresponding element of z should also flip to preserve f. But z and z after flipping the sign of an element may be far from one another in the unconstrained space.
I’ve found normalizing the eigenvector matrix by ensuring elements in the first row are all positive to accelerate sampling in a toy example (sampling a GP using the eigen rather than Cholesky decomposition). Is this sort of normalization common practice? Are there better approaches to resolve the sign ambiguity?