[SciPy-dev] problem with linalg.cholesky?

[Andrew Jaffe]

Hi All,

OK, drilling down a bit further... The problem lies in the use of 
scipy.transpose() in _castCopyAndTranspose, since transpose() makes a 
contiguous, non-Fortran array into a non-contiguous fortran array (since 
it appears that at present it doesn't copy, just makes a new view). One 
solution is just to actually copy the array. The other function, 
_fastCopyAndTranspose(), seems to work fine in this situation, but I'm 
not sure how and why this differs from _castCopyAndTranspose (in 
particular, it doesn't seem to set the Fortran flag).

[However, the cholesky decomposition only makes sense for a symmetric 
(or Hermitian) array, so, if we're not changing the type, for real 
matrices really all we need to do is the cast and copy; for complex 
matrices we need to transpose (not hermitian conjugate since we really 
want to change the ordering), but this is equivalent to just taking the 
complex conjugate of every element, which must be faster than copying. 
If we are changing the type, I assume the speed difference isn't as 
important.]

Any ideas on what the best change(s) would be?

Andrew


> A bit more information on the problem:
> 
> Following a suggestion in the pickling thread going on here, ff I add in 
> lines to the cholesky_decomposition function in basic_lite.py to print 
> out the flags of the array, the following happens:
> 
> 
> input array flags {'WRITEABLE': True, 'UPDATEIFCOPY': False, 
> 'NOTSWAPPED': True, 'CONTIGUOUS': True, 'FORTRAN': False, 'ALIGNED': 
> True, 'OWNDATA': True}
> 
> cholesky: after _castCopyAndTranspose, flags= {'WRITEABLE': True, 
> 'UPDATEIFCOPY': False, 'NOTSWAPPED': True, 'CONTIGUOUS': False, 
> 'FORTRAN': True, 'ALIGNED': True, 'OWNDATA': True}
> 
> So, for some reason, _castCopyAndTranspose is somehow making the array 
> non-Contiguous, which is then a problem in the underlying lapack call.
> 
> What's going on here?
> 
> Also, an (unrelated?) point: I note that basic_lite.py exists in both 
> scipy/linalg and scipy/basic -- two different files with somewhat 
> different contents, imported in different ways! Which is correct? Should 
> they both exist? (They both give the same error for cholesky.) Is one a 
> leftover from an earlier installation?
> 
> Thanks,
> 
> Andrew
> 
> 
> 
> 
> Andrew Jaffe wrote:
> 
>>Hi All,
>>
>>One more data point: in addition to my problems on OS X, someone else 
>>has the same issue on Win XP... So it's not OS-specific, but it does 
>>occur only for some installations.
>>
>>Perhaps it depends on the underlying BLAS/LAPACK library?
>>
>>Andrew
>>
>>Andrew Jaffe wrote:
>>
>>
>>>hi all,
>>>
>>>(apologies that a similar question has appeared elsewhere...)
>>>
>>>In the newest incarnation of scipy_core, I am having trouble with the
>>>cholesky(a) routine. Here is some minimal code reproducing the bug  
>>>(on OS X)
>>>
>>>------------------------------------------------------------
>>
>>>from scipy import identity, __core_version__, Float64
>>
>>>import scipy.linalg as la
>>>print 'Scipy version: ', __core_version__
>>>i = identity(4, Float64)
>>>print 'identity matrix:'
>>>print i
>>>
>>>print 'about to get cholesky decomposition'
>>>c = la.cholesky(i)
>>>print c
>>>------------------------------------------------------------
>>>
>>>which gives
>>>
>>>------------------------------------------------------------
>>>Scipy version:  0.6.1
>>>identity matrix:
>>>[[ 1.  0.  0.  0.]
>>> [ 0.  1.  0.  0.]
>>> [ 0.  0.  1.  0.]
>>> [ 0.  0.  0.  1.]]
>>>about to get cholesky decomposition
>>>Traceback (most recent call last):
>>>  File "/Users/jaffe/Desktop/bad_cholesky.py", line 13, in ?
>>>    c = la.cholesky(i)
>>>  File "/Library/Frameworks/Python.framework/Versions/2.4/lib/ 
>>>python2.4/site-packages/scipy/linalg/basic_lite.py", line 117, in  
>>>cholesky_decomposition
>>>    results = lapack_routine('L', n, a, m, 0)
>>>lapack_lite.LapackError: Parameter a is not contiguous in  
>>>lapack_lite.dpotrf
>>>------------------------------------------------------------
>>>
>>>(The cholesky decomposition in this case should just be the matrix  
>>>itself; the same error occurs with a complex matrix.)
>>>
>>>
>>>Any ideas? Could this have anything to do with _CastCopyAndtranspose  
>>>in Basic_lite.py? (Since there are few other candidates for anything  
>>>that actually changes the matrix.)
>>>
>>>Thanks in advance,
>>>
>>>A
>>>
>>>
>>>______________________________________________________________________
>>>Andrew Jaffe                                    a.jaffe at imperial.ac.uk
>>>Astrophysics Group                                    +44 207 594-7526
>>>Blackett Laboratory, Room 1013                    FAX             7541
>>>Imperial College, Prince Consort Road
>>>London SW7 2AZ ENGLAND              http://astro.imperial.ac.uk/~jaffe