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
>
>
