Re: matrix_type check

[Jaroslav Hajek]

On Fri, Apr 25, 2008 at 10:21 AM, David Bateman
<address@hidden> wrote:
> Jaroslav Hajek wrote:
>  > I'm not sure I understand. If this check succeeds, but Cholesky
>  > factorization fails,
>  > then the matrix is rank-deficient or has negative eigenvalues (yes,
>  > that can also happen).
>  > But it is much harder to create a matrix that is misclassified by this
>  > check - you can try.
>  > With the old check, just do
>  > A = rand(10); A = A + A';
>  > matrix_type (A)
>  > chol(A)
>  >
>  Currently the mldivide operator on a matrix flagged as positive definite
>  will do the following
>
>  * Attempt Cholesky factorization
>  * If that fails attempt LU factorization
>  * If that fails attempt a minimum norm solution
>
>  I was wondering if the second step should be removed. Of course, it was
>  a stupid question on my part... Your check can't see if there are
>  negative eigenvalues, and so the fallback to an LU factorization should
>  be kept.
>

Yes, I think so. Simply checking the diagonal for positivity
essentially means checking
the positive definiteness of projections on every 1-dimensional
subspace spanned by
a single column of identity matrix. The new check uses, in addition,
projections to 2-dimensional subspaces spanned by each pair of
columns. It is not certainly sufficient,
in fact I think that even (n-1)-dimensional subspaces would be insufficient.

>
>  >>  If so then the current factorization code should also be changed such
>  >>  that a failing Choleksy factorization falls back to a minimum norm
>  >>  solution rather than first trying an LU solution.
>  >>
>  >>
>  >
>  > Ah, now I get it. No, I don't think so. I think that the test can
>  > still pass even for a regular matrix with negative eigenvalues. It is
>  > an interesting question though - I'll try to research this a little,
>  > perhaps your guess is right.
>  >
>  I'm no longer sure unless you we can guarantee that a failing Cholesky
>  factorization is due to a rank deficient matrix rather than negative
>  eigenvalues. I suppose we should check for symmetric definite matrices
>  as well and use DSYTRF and ZHETRF to do the factorization.
>

I'm not sure how typical symmetric indefinite matrices are - I think I
have never used the Bunch-Kaufman factorization routines at all (but
that means only little given my short experience). Perhaps someone
other could comment on that. Certainly that would matrix division yet
smarter. Currently, "hermitian" in MatrixType is used for SPD/HPD,
thus one would need to add "hermitian_indef" or change "hermitian" to
"hermitian_posdef".


>
>
>  D.
>
>  --
>  David Bateman                                address@hidden
>  Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph)
>  Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob)
>  91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax)
>
>  The information contained in this communication has been classified as:
>
>  [x] General Business Information
>  [ ] Motorola Internal Use Only
>  [ ] Motorola Confidential Proprietary
>
>



-- 
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz