help-octave
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

chol(a)


From: Salvatore Filippone
Subject: chol(a)
Date: Tue, 15 Nov 2011 11:38:41 +0100

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


reply via email to

[Prev in Thread] Current Thread [Next in Thread]