On Wednesday, November 05, 2014 01:19:55 AM Stefan Karpinski wrote: > I know that factorize checks a bunch of properties of the matrix to be > factorized – is positive definiteness not something that it checks? Should > it?
Positive definiteness is not a quick check. For example, the matrix `ones(2,2)` is symmetric and has all positive entries but is not positive-definite. You have to finish computing the Cholesky factorization before you can be sure it's positive definite, at which point you should of course just keep that factorization. --Tim > > On Tue, Nov 4, 2014 at 5:59 PM, Andreas Noack <[email protected]> > wrote: > > 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
