Andreas, Thanks for the detailed explanations, that's great to know.
It would be useful indeed to include some more details on the matrix hierarchy in the documentation. I had been looking for a type indicating positive definiteness, but the closest I could find was Hermitian (Cholesky is not mentioned in 1.18.2, and there is no constructor exported for the user). `chol(a)` does the factorization, but does not bless the result to be of type Cholesky, probably for a good reason. Indeed inv(factorize(pd)) is slightly faster than inv(pd) (for very large matrices the differences becomes more noticeable), but more importantly for my application: it retains positive definiteness in the result! Thanks agin, ---david On Wednesday, November 5, 2014 1:59:53 AM UTC+1, Andreas Noack wrote: > > factorize(Matrix) is a full service check, but if you specify structure as > David did with Hermitian, factorize dispatches to an appropriate method. In > this case it is Bunch-Kaufman. E.g. > > julia> A = randn(3,3);A = A'A; > > julia> typeof(factorize(A)) > Cholesky{Float64} (constructor with 1 method) > > julia> typeof(factorize(Hermitian(A))) > BunchKaufman{Float64} (constructor with 1 method) > > 2014-11-04 19:31 GMT-05:00 Tim Holy <tim....@gmail.com <javascript:>>: > >> 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 <andreasno...@gmail.com >> <javascript:>> >> > 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 <david.va...@gmail.com >> <javascript:>>: >> > > >> > > 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 >> >> >