[SciPy-user] bug in cho_factor?
Emanuele Olivetti
emanuele@relativita....
Sat Jul 19 12:50:47 CDT 2008
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@relativita.com <mailto:emanuele@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@scipy.org <mailto:SciPy-user@scipy.org>
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>
More information about the SciPy-user
mailing list