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