[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: matrix_type check
From: |
Jaroslav Hajek |
Subject: |
Re: matrix_type check |
Date: |
Fri, 25 Apr 2008 09:58:12 +0200 |
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