Re: chol(a)

[c.]

On 15 Nov 2011, at 11:38, Salvatore Filippone wrote:

> Hello there,
> While checking exercises for some students, I came across a rather
> surprising (for me) behaviour.
> Consider the following SPD matrix:
> octave:1> a=[12,2,3,4;
>> 2,14,5,3;
>> 3,5,16,6;
>> 4,3,6,16]
> a =
> 
>   12    2    3    4
>    2   14    5    3
>    3    5   16    6
>    4    3    6   16
> 
> If I perform the Cholesky factorization on the lower or upper part, I
> get the expected results;
> octave:2> chol(a)
> ans =
> 
>   3.46410   0.57735   0.86603   1.15470
>   0.00000   3.69685   1.21725   0.63117
>   0.00000   0.00000   3.71057   1.14045
>   0.00000   0.00000   0.00000   3.60107
> 
> octave:4> chol(a,'lower')
> ans =
> 
>   3.46410   0.00000   0.00000   0.00000
>   0.57735   3.69685   0.00000   0.00000
>   0.86603   1.21725   3.71057   0.00000
>   1.15470   0.63117   1.14045   3.60107
> 
> Now, it is my understanding that internally Octave uses LAPACK; if I
> ask for the factorization of an upper triangle in LAPACK, the code
> will only look at the upper part of the input.
> This is consistent with what I get if I pass triu(a) as in
>> chol(triu(a))
> ans =
> 
>   3.46410   0.57735   0.86603   1.15470
>   0.00000   3.69685   1.21725   0.63117
>   0.00000   0.00000   3.71057   1.14045
>   0.00000   0.00000   0.00000   3.60107
> 
> In LAPACK, I can do a perfectly symmetric choice by passing only the
> lower triangle, but Octave in this case seems to only factor the
> diagonal:
> octave:5> chol(tril(a),'lower')
> ans =
> 
>   3.46410   0.00000   0.00000   0.00000
>   0.00000   3.74166   0.00000   0.00000
>   0.00000   0.00000   4.00000   0.00000
>   0.00000   0.00000   0.00000   4.00000
> 
> 
> Having been an LAPACK user for a long time this is a bit surprising.
> Note that "that other package" (i.e. Matlab) does what I expect it to
> do, i.e. with choll(tril(a),'lower') produces the transpose of the
> output from chol(triu(a))
> 
> Bug or feature? If the latter, then maybe the documentation ought to
> be more specific.
> 
> Thanks a lot
> Salvatore

Dear Salvatore,

based on the documentation, I'm not sure that chol is meant to work as you 
expect, 
on the other hand matlab seems to do what you expect.

Incompatibilities with matlab are usually regarded as bugs in Octave, so I'd 
recommend
reporting this as a bug on http://savannah.gnu.org/bugs/?group=octave

Thanks,
Carlo