scipy.linalg.lu_solve gives different results depending on the ordering of the input matrix (hence one of the results is wrong). For example, see the following Python commands: $ ipython -pylab ... In []: import scipy.linalg In []: a = rand(30, 5); A = dot(a.T, a); L = cholesky(A) ln []: b = rand(5); pivot = arange(5) In []: scipy.linalg.lu_solve((L, pivot), b) Out[]: array([ 0.05650387, -0.51047795, 1.15888144, -1.21440483, 0.18492365]) In []: scipy.linalg.lu_solve((L.copy(), pivot), b) Out[]: array([ 0.05650387, 0.01303743, 0.03212762, 0.03715339, 0.0364978 ]) Notice that the last two commands yield totally different results. Why? Inside the python code of lu_solve, depending on the order of L being 'C' or 'F', the function clapack.dgetrf or the function flapack.dgetrf is called. Usually, the two versions give the same interface, and the hope is that scipy will choose a more appropriate / faster one according to the order of L. However, the problem is that clapack.dgetrf and flapack.dgetrf don't take the same API, even though they both can deal with different matrix ordering: In []: from scipy.linalg.lapack import flapack, clapack In []: (flapack.dgetrf(L)[0] == flapack.dgetrf(L.copy())[0]).all() Out[]: True In []: (clapack.dgetrf(L)[0] == clapack.dgetrf(L.copy())[0]).all() Out[]: True In []: flapack.dgetrf(L)[0] Out[]: array([[ 9.76176212, 0. , 0. , 0. , 0. ], [ 0.70830302, 6.53864159, 0. , 0. , 0. ], [ 0.76641246, 0.47969069, 6.04051948, 0. , 0. ], [ 0.80153942, 0.47344331, 0.32977516, 5.89748659, -0. ], [ 0.7475908 , 0.3921305 , 0.36098544, 0.25955716, 5.41902844]]) In []: clapack.dgetrf(L)[0] Out[]: array([[ 9.76176212, 0. , 0. , 0. , 0. ], [ 6.91428558, 6.53864159, -0. , -0. , -0. ], [ 7.4815361 , 3.13652549, 6.04051948, -0. , -0. ], [ 7.82443714, 3.09567614, 1.99201329, 5.89748659, -0. ], [ 7.29780353, 2.56400081, 2.18053957, 1.5307349 , 5.41902844]]) This is reproducible with Fedora 13 on an Intel machine.
Fwiw, I can also reproduce the problem with Fedora 12.
Yeah this feels like something I need to have a discussion with upstream about. In the meantime do you have a suggestion on how to provide a workaround for Fedora users? Can I force lu_solve to use the dgetrf call that is working correctly? -jef
Oops, I said something wrong in this bug report (although the bug is real). dgetrf is used in lu_factor, and dgetrs is used in lu_solve (I said lu_solve calls dgetrf; that was wrong). Anyway, both lu_factor and lu_solve have this inconsistency behavior. It seems the Fortran lib (flapack) has the right behavior (matches the pydoc). I can think of two workarounds: a) (at package level) don't ship the clapack lib. thus scipy will be forced to only use the Fortran lib. OR b) (at user level) always use lu_factor and lu_solve in pairs. If one has to use only one of them, pass a copy of the numpy array with something like numpy.asarray(orig_array, order='F') a) is nicer since it doesn't have as severe a performance penalty as b), and it's transparent to the users. Thanks.
This message is a reminder that Fedora 13 is nearing its end of life. Approximately 30 (thirty) days from now Fedora will stop maintaining and issuing updates for Fedora 13. It is Fedora's policy to close all bug reports from releases that are no longer maintained. At that time this bug will be closed as WONTFIX if it remains open with a Fedora 'version' of '13'. Package Maintainer: If you wish for this bug to remain open because you plan to fix it in a currently maintained version, simply change the 'version' to a later Fedora version prior to Fedora 13's end of life. Bug Reporter: Thank you for reporting this issue and we are sorry that we may not be able to fix it before Fedora 13 is end of life. If you would still like to see this bug fixed and are able to reproduce it against a later version of Fedora please change the 'version' of this bug to the applicable version. If you are unable to change the version, please add a comment here and someone will do it for you. Although we aim to fix as many bugs as possible during every release's lifetime, sometimes those efforts are overtaken by events. Often a more recent Fedora release includes newer upstream software that fixes bugs or makes them obsolete. The process we are following is described here: http://fedoraproject.org/wiki/BugZappers/HouseKeeping
This also affects Fedora 14 and 15.
Ling, Have you brought this to the attention of upstream scipy? I'm hesitant on specifically not shipping the clapack to paper over the problem. And I don't know enough about the details for to have a constructive discussion about what the fix needs to be in the codebase. -jef
Fixed in upstream in release 0.10.0: http://projects.scipy.org/scipy/ticket/1458
This message is a notice that Fedora 15 is now at end of life. Fedora has stopped maintaining and issuing updates for Fedora 15. It is Fedora's policy to close all bug reports from releases that are no longer maintained. At this time, all open bugs with a Fedora 'version' of '15' have been closed as WONTFIX. (Please note: Our normal process is to give advanced warning of this occurring, but we forgot to do that. A thousand apologies.) Package Maintainer: If you wish for this bug to remain open because you plan to fix it in a currently maintained version, feel free to reopen this bug and simply change the 'version' to a later Fedora version. Bug Reporter: Thank you for reporting this issue and we are sorry that we were unable to fix it before Fedora 15 reached end of life. If you would still like to see this bug fixed and are able to reproduce it against a later version of Fedora, you are encouraged to click on "Clone This Bug" (top right of this page) and open it against that version of Fedora. Although we aim to fix as many bugs as possible during every release's lifetime, sometimes those efforts are overtaken by events. Often a more recent Fedora release includes newer upstream software that fixes bugs or makes them obsolete. The process we are following is described here: http://fedoraproject.org/wiki/BugZappers/HouseKeeping