In your case, I think the right solution is to invert it by inv(cholfact(pd)). By calling cholfact, you are telling Julia that your matrix is positive definite and Julia will exploit that to give you a fast inverse which is also positive definite.
inv(factorize(Hermitian(pd))) is slow because it uses a factorization that only exploits symmetry (Bunch-Kaufman), but not positive definiteness (Cholesky). However, the Bunch-Kaufman factorization preserves symmetry and hence the result is positive definite. In contrast, when doing inv(pd) Julia has no idea that pd is positive definite or even symmetric and hence it defaults to use the LU factorization which won't preserve symmetry and therefore isposdef will return false. Hope it made sense. I'll probably have to write a section in the documentation about this soon. 2014-11-03 18:53 GMT-05:00 David van Leeuwen <[email protected]>: > Hello, > > I am struggling with the fact that covariance matrices computed from a > precision matrix aren't positive definite, according to `isposdef()` (they > should be according to the maths). > > It looks like the culprit is `inv(pd::Matrix)` which does not always > result in a positive definite matrix if `pd` is one. This is probably > because `inv()` is agnostic of the fact that the argument is positive > definite, and numerical details. > > Now I've tried to understand the support for special matrices, and I > believe that `inv(factorize(Hermitian(pd)))` is the proper way to do this. > Indeed the resulting matrix is positive definite. However, this > computation takes a lot longer than inv(), about 5--6 times as slow. I > would have expected that the extra symmetry would lead to a more efficient > matrix inversion. > > Is there something I'm doing wrong? > > Cheers, > > ---david >
