Re: chol problem in 2.9.5

[David Bateman]

Quentin Spencer wrote:
> David Bateman wrote:
> 
>>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...
>>  
> 
> 
> I hadn't realized there were redundant versions of chol, but I can 
> verify that both RPMs have a version:
> /home/qspencer> rpm -ql octave | grep chol.oct
> /usr/libexec/octave/2.9.5/oct/i686-redhat-linux-gnu/chol.oct
> /usr/libexec/octave/2.9.5/oct/i686-redhat-linux-gnu/spchol.oct
> /home/qspencer> rpm -ql octave-forge | grep chol.oct
> /usr/libexec/octave/site/oct/api-v18/i686-redhat-linux-gnu/octave-forge/chol.oct
> 
> So, I could fix this by just excluding chol.oct from the octave-forge 
> RPM. Is there any downside to doing this with octave 2.9.5?
> 


Yes there is a downside, in that the octave-forge version of chol flags
the matrix as upper triangular in a special user type and uses only
back/forward subsititution for "r = chol(A); x = r \ (r' \ b)", rather
than a full LU solve (twice). So even though the LAPACK
back-substitution code isn't used in this version it should still be
much faster for the solve than the octave 2.9.5 version. However, the
octave 2.9.5+ version should beat it, as it automatically recognizes r
as being triangular and uses LAPACK's back/forward substitution code in
the solve.

D.