Strange behaviour of \ w.r.t. chol

[Marco Caliari]

Dear David,

matrix B is not singular, as you can see from its eigenvalues (and from 
Matlab). Moreover, the Cholesky factorization can be successfully 
computed, as you can see running chol(B). So, if the test, even if 
weak, detects B to be 'Positive Definite' and if chol(B) can be computed, 
why det(B)=0 and B\y is not solved using Cholesky?

Best regards,

Marco

On Tue, 19 Feb 2008, David Bateman wrote:

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