Yes and no. The problem is that small floating point noise destroys symmetry/Hermitianity and therefore factorize concludes that the matrix is not positive definite (I know about the other definition). If you construct your positive definite matrix by A'A, then Julia makes it exactly symmetric/Hermitian, but if you do e.g. A'*Diagonal(rand(size(A,1)))*A then your matrix is still positive definite in infinite precision, but it is almost never exactly symmetric/Hermitian in finite precision.
Hence it is a good idea to use inv(cholfact(X)) whenever you know that your matrix should be considered positive definite. This is also much faster as it bypasses all the checks in factorize. 2014-11-04 20:22 GMT-05:00 Stefan Karpinski <[email protected]>: > So this works as desired if you just do inv(factorize(X))? > > On Wed, Nov 5, 2014 at 1:59 AM, Andreas Noack < > [email protected]> 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 <[email protected]>: >> >> 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 >>> >>> >> >
