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 <andreasnoackjen...@gmail.com>
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 <ste...@karpinski.org>:
>
> So this works as desired if you just do inv(factorize(X))?
>>
>> On Wed, Nov 5, 2014 at 1:59 AM, Andreas Noack <
>> andreasnoackjen...@gmail.com> 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.h...@gmail.com>:
>>>
>>> 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 <
>>>> andreasnoackjen...@gmail.com>
>>>> > 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.vanleeu...@gmail.com>:
>>>> > >
>>>> > > 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
>>>>
>>>>
>>>
>>
>
