Re: matrix_type check

[Jaroslav Hajek]

On Fri, Apr 25, 2008 at 9:41 AM, David Bateman
<address@hidden> wrote:
> Jaroslav Hajek wrote:
>  > hello,
>  > please consider the attached changeset. It extends the current
>  > matrix_type check for positive definiteness
>  > (positive diagonal + symmetry) by also checking the necessary
>  > condition a(i,j)*a(j,i) < a(i,i) * a(j,j) for i != j,
>  > which can also be done in one pass through the matrix, and thus does
>  > not impose much additional cost,
>  > while greatly reducing the number of misclassified matrices (for
>  > example `a = rand(10); matrix_type (a+a')'
>  > is almost always wrong).
>  > A minor bug is also fixed: in current version, setting matrix_type via
>  > matrix_type (a, "type") always produces a warning.
>  > Few tests in matrix_type.cc relying on the old more optimistic
>  > behaviour are also fixed.
>  >
>  > cheers,
>  >
>  >
>  I haven't extensively checked the patch, but if it works then why not.

The property a(i,j)*a(j,i) < a(i,i)*a(j,j) (note that for complex
hermitian matrices both
sides are real) is a necessary (not sufficient) condition for positive
definiteness. This check can be done along with the check for
symmetry, so it costs only a little more than the old method. Cholesky
factorization is surely better but that is already an O(N^3) operation
(although the failure is likely to occur earlier in the process).


>  My positive definitebess test was alwaus a good enough test as the
>  Cholesky factorization would catch any misclassifications.. If this is
>  an accurate test for positive definiteness then does that mean that any
>  failure of the Cholesky factorization is due to a rank deficient matrix?

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)


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

>  About the only thing I disagree with in your changeset is the change
>
>  @@ -106,7 +101,7 @@
>         typ = MatrixType::Lower;
>        else if (hermitian)
>         typ = MatrixType::Hermitian;
>  -      else if (ncols == nrows)
>  +      else
>         typ = MatrixType::Full;
>      }
>    else
>
>  If the matrix is not square then the xGELSD should always be used and
>  the matrix should be flagged as singular to force this..
>

I think that the whole block (you see the closing brace of it) is
already inside an
`if (ncols == nrows)' clause, so that's why the check is redundant.
But I'll check that again.

regards,

>  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