Open mjhajharia opened 2 years ago
@sethaxen maybe I can start with correlation matrices and you can see covariance matrices
Is this for Cholesky factors? Or for the whole thing? We usually go with Cholesky-factor parameterizations, so to me, those are more important. But for symmetric positive definite, that's just an unconstraining transform from (0, inf) to (-inf, inf) for the diagonals.
The alternative is to take a symmetric matrix and push it through matrix exp. I have no clue as to what the Jacobian determinant looks like there.
I think correlation matrices should be broken out into a separate issue. Adam Haber has looked at those.
@sethaxen maybe I can start with correlation matrices and you can see covariance matrices
Works for me!
But for symmetric positive definite, that's just an unconstraining transform from (0, inf) to (-inf, inf) for the diagonals.
Agreed, we could just take the "best" Cholesky transform with the transform for the positive diagonals.
Another alternative is to use the "best" transform for orthonormal matrices for the eigenvectors and then a transform for the positive eigenvectors. I think I've seen other approaches in the literature before and will look into it.
The alternative is to take a symmetric matrix and push it through matrix exp. I have no clue as to what the Jacobian determinant looks like there.
From http://eprints.maths.manchester.ac.uk/2754/1/catalog.pdf, the eigendecomposition is still the preferred way to compute exp(A)
for symmetric A
. i.e. if A=U Λ U'
, then exp(A) = U exp(Λ) U'
. This is the composition of 3 bijective maps: eigendecomposition, diagonal exponential, then inverse of eigendecomposition, so we only need to work out the Jacobian correction from symmetric eigendecomposition. From Alan Edelman's lecture notes (Section 2.5 of http://web.mit.edu/18.325/www/handouts/handout3.pdf), the Jacobian correction for symmetric eigendecomposition is $\prod_{i < j} |\lambda_j - \lambda_i|$. So the correction for symmetric matrix exponential would be
$$\prod_{i=1}^n e^{\lambdai} \prod{j=i+1}^n \left| \frac{e^{\lambda_j} - e^{\lambda_i}}{\lambda_j - \lambda_i} \right |$$
I verified this using finite differences in Julia:
I'd expect ADing through this transform would be quite a bit more expensive than for the Cholesky transform, which only requires adding the unconstrained parameters behind the positive diagonals.
I think @adamhaber already has the correlation matrix and Cholesky for correlation matrix transforms.
@bob-carpenter I'm not certain if your comment is in reply to me, but I was only referring to symmetric positive definite matrices.
I verified this using finite differences in Julia:
Super cool! That correction is $\mathcal{O}(n^2)$. We have matrix_exp
implemented in Stan using the Al-Mohy and Higham algorithm.
f780af0 adds derivations for this section
The set of SPD matrices is not compact, so we need to add proper distributions. The Wisharts are conjugate for covariance in a multivariate normal model, so they're essentially like testing multivariate normal posteriors.