There are several chapters of Trefethen, or Demmel, on these exact topics. Just have to translate their pseudocode or Matlab addenda into Julia.
On Tuesday, November 4, 2014 5:41:55 PM UTC-8, Stefan Karpinski wrote: > > Ah, good to know. There's so much depth here. This could be a chapter or > two of a book on effective numerical analysis in Julia :-) > > On Wed, Nov 5, 2014 at 2:33 AM, Andreas Noack <[email protected] > <javascript:>> wrote: > >> 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] >> <javascript:>>: >> >> 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] >>> <javascript:>> 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] <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 < >>>>> [email protected] <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 < >>>>> [email protected] <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 >>>>> >>>>> >>>> >>> >> >
