[SciPy-user] bug in cho_factor?

[Warren Weckesser]

It appears that cho_factor is a wrapper for the lapack routine potrf.
This function doesn't zero out the irrelevant part of the matrix; it
appears to leave it untouched.  In your case, the upper triangular
part of your answer is correct; the subdiagonal part is irrelevant
(and is apparently left over from the input matrix).

Apparently this is not a problem, since the output of cho_factor is
intended to be used in cho_solve (as the docstring states).  Presumably
cho_solve ignores the irrelevant parts of the input matrix.

Warren

On Sat, Jul 19, 2008 at 11:29 AM, Emanuele Olivetti <emanuele at relativita.com>
wrote:

> In scipy v0.7.0.dev4543 (svn) the documentation of
> scipy.linalg.cho_factor says that cho_factor returns a tuple
> made of:
> ----
>    c : array, shape (M, M)
>        Upper- or lower-triangular Cholesky factor of A
>    lower : array, shape (M, M)
>        Flag indicating whether the factor is lower or upper triangular
> ----
> But the following simple example shows that 'c' is not triangular
> and the result of cho_factor() is pretty different from that of
> cholesky()!
> ----
> import numpy as N
> import scipy.linalg as SL
> A = N.array([[2,1],[1,2]])
> c_wrong,lower = SL.cho_factor(A)
> print c_wrong
> c_correct = SL.cholesky(A)
> print c_correct
> print c_wrong==c_correct
> ----
> The output is:
> ----
> [[ 1.41421356  0.70710678]
>  [ 1.          1.22474487]]
> [[ 1.41421356  0.70710678]
>  [ 0.          1.22474487]]
> [[ True  True]
>  [False  True]]
> ----
> Why c_wrong[1,0] is not zero? And why is it 1.0 ?
>
> I went mad tracking this issue in my code 8-|
>
> Regards,
>
> Emanuele
>
> P.S.: there is a type in cho_factor documentation: 'lower' is
> not an 'array' but an 'int'.
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20080719/0cb6493f/attachment.html>