Strange behaviour of \ w.r.t. chol

[David Bateman]

Marco Caliari wrote:
> Hi.
>
> Please consider the following script
>
> n = 11;
> x = linspace(0,10,n)';
> m = 7;
> A = (x * ones (1, m+1)) .^ (ones (n, 1) * (m : -1 : 0));
> B = A'*A;
> matrix_type(B)
> issymmetric(B)
> eig(B)
> det(B)
> y = A'*sin(x);
> B\y
> R = chol(B);
> R\(R'\y)
>
> I get det(B)=0, even if matrix_type(B)='Positive Definite' and no 
> eigenvalue is 0. Moreover, I would expect that \ uses the Cholesky 
> factorization, but the result of B\y is different from R\(R'\y). Matlab 
> gives det(B) = 5.32e+32, a warning of bad scaled matrix, but the same 
> results. Octave has no problem with m = 6.
> Finally, a minor thing: why B=matrix_type(B,matrix_type(B)) gives a 
> warning (invalid matrix type)?
> I tried Octave 3.0.0 with atlas_sse2 libs (atlas_sse2) and Octave 2.9.12 
> with atlas_base libs.
>
> Best regards,
>
> Marco
>
>   
Marco,

I've explained this issue on the lists previous, so perhaps the
documentation of matrix_type should be updated. The test for positive
definiteness is a weak test. It just tests for probably positive
definiteness, looking to see if the matrix is hermitian with a positive
real diagonal. A full test for positive definite matrices would take as
long as the factorization. How the result of matrix_type is then used is
that the matrix being flagged as probably positive definite a cholesky
factorization is attempted, if that fails the factorization falls back
to an LU factorization and if that fails the factorization falls back to
xGELSD in the case of a full matrix and a Dulmange-Mendhelsohn QR
blocked factorization in the case of a sparse matrix.

If you run matrix_type(B) after the above, you'll find its now marked as
being singular. How the factorization is really done (ignoring the
output of matrix_type) is documented in sections 18.1 and 20.2 of the
manual.

D.







This is explained in the documentation of matrix_type

-- 
David Bateman                                David.Bateman at motorola.com
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