help-octave
[Top][All Lists]
Advanced

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

Re: chol problem in 2.9.5


From: David Bateman
Subject: Re: chol problem in 2.9.5
Date: Thu, 25 May 2006 20:59:05 +0200
User-agent: Mozilla Thunderbird 1.0.6-7.6.20060mdk (X11/20050322)

Dmitri A. Sergatskov wrote:
> octave:1> a=randn(2)
> a =
> 
>   -0.94562   1.05567
>   -0.69245   1.99836
> 
> octave:2> x=a'*a
> x =
> 
>    1.3737  -2.3820
>   -2.3820   5.1079
> 
> octave:3> chol(x)
> ans =
> 
>    1.17204  -2.03238
>    0.00000   0.98860
> Upper Triangular
> 
> octave:4> chol(x)
> error: chol: matrix not positive definite
> octave:4>
> ---
> 
> So after calling chol a second time on the same matrix it gives an error.
> That seems breaks chol for good for remaining of the octave session:

This is the old octave-forge chol.cc function being used here. You can
tell since it prints "Upper Triangular". So strictly speaking this is a
bug in octave-forge.

The fact is in the octave CVS I submitted a patch to correctly determine
is a matrix is positive definite and use the LAPACK cholesky
factorization functions. In addition to special cases for upper/lower
triangular functions using the lapack back substitution code rather than
the home rolled version in the octave-forge code, so the version in
octave's core should now be significantly better than the octave-forge
forge.

This completely removes the need for the octave-forge code, and so I
disabled it in octave-forge (though I didn't delete it is a good example
of how to write a user type). It also means I don't see much point in
tracking down this bug...

In 2.9.5+ I get

octave:1> a = [-0.94562   1.05567;-0.69245   1.99836];
octave:2>  x=a'*a
x =

  1.3737  -2.3820
  -2.3820  5.1079

octave:3> chol(x)
ans =

  1.17204  -2.03237
  0.00000  0.98861

octave:4> chol(x)
ans =

  1.17204  -2.03237
  0.00000  0.98861

octave:5> matrix_type(x)
ans = Positive Definite
octave:6> matrix_type(chol(x))
ans = Upper

Regards
David


reply via email to

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