[SciPy-user] bug in cho_factor?

[Emanuele Olivetti]

Thanks a lot for the explanation. I believe cho_factor's docstring
should be updated in order to mention these facts. It is definitely
unexpected that the result of the two decompositions are different
and can cause problems like I had (a couple of hours spent). A
clear "Warning" should fit. Consider that U,lower=cho_factor(A)
outputs an U that does not satisfy A==N.dot(U.T,U) !!

Do you know why cho_factor does not zeros out the matrix?
Is it for performance reasons? Otherwise there are no reasons
for using cho_factor() instead of cholesky().

Best Regards,

Emanuele


Warren Weckesser wrote:
> Emanuele,
>
> After a closer look at the code, it appears that the only difference
> between cholesky() and cho_factor() is that cho_factor() does not
> set the irrelevant terms to zero.  So, if you are only going to
> use the result to call cho_solve(), you can use cho_factor(), but
> if you need the correct square triangular matrix with zeros in the
> subdiagonal, use cholesky().
>
> Warren
>
>
> On Sat, Jul 19, 2008 at 11:29 AM, Emanuele Olivetti
> <emanuele at relativita.com <mailto: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 <mailto:SciPy-user at scipy.org>
>     http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>