Re: matrix_type check

[David Bateman]

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

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

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