Bug 623000 - Inconsistent C/Fortran API causing wrong behavior of scipy.linalg.lu_solve
Inconsistent C/Fortran API causing wrong behavior of scipy.linalg.lu_solve
Status: CLOSED WONTFIX
Product: Fedora
Classification: Fedora
Component: scipy (Show other bugs)
15
i686 Linux
low Severity medium
: ---
: ---
Assigned To: Jef Spaleta
Fedora Extras Quality Assurance
:
Depends On:
Blocks:
  Show dependency treegraph
 
Reported: 2010-08-10 22:59 EDT by Ling Li
Modified: 2012-08-07 15:33 EDT (History)
3 users (show)

See Also:
Fixed In Version:
Doc Type: Bug Fix
Doc Text:
Story Points: ---
Clone Of:
Environment:
Last Closed: 2012-08-07 15:33:10 EDT
Type: ---
Regression: ---
Mount Type: ---
Documentation: ---
CRM:
Verified Versions:
Category: ---
oVirt Team: ---
RHEL 7.3 requirements from Atomic Host:
Cloudforms Team: ---


Attachments (Terms of Use)

  None (edit)
Description Ling Li 2010-08-10 22:59:04 EDT
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.
Comment 1 Ling Li 2010-08-10 23:00:58 EDT
Fwiw, I can also reproduce the problem with Fedora 12.
Comment 2 Jef Spaleta 2010-08-12 12:42:08 EDT
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
Comment 3 Ling Li 2010-08-12 22:31:12 EDT
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.
Comment 4 Bug Zapper 2011-06-01 07:41:41 EDT
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
Comment 5 Ling Li 2011-06-01 12:33:58 EDT
This also affects Fedora 14 and 15.
Comment 6 Jef Spaleta 2011-06-17 08:37:50 EDT
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
Comment 7 pav 2012-06-13 10:08:21 EDT
Fixed in upstream in release 0.10.0: 
http://projects.scipy.org/scipy/ticket/1458
Comment 8 Fedora End Of Life 2012-08-07 15:33:13 EDT
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

Note You need to log in before you can comment on or make changes to this bug.