[Scipy-tickets] [SciPy] #1458: lu_factor gives different results, depending on whether the input is C or Fortran contiguous (also affects lu_solve)

SciPy Trac scipy-tickets@scipy....
Thu Jun 9 09:49:23 CDT 2011


#1458: lu_factor gives different results, depending on whether the input is C or
Fortran contiguous (also affects lu_solve)
--------------------------+-------------------------------------------------
 Reporter:  mike.wimmer   |       Owner:  somebody
     Type:  defect        |      Status:  new     
 Priority:  normal        |   Milestone:          
Component:  scipy.linalg  |     Version:  0.9.0   
 Keywords:                |  
--------------------------+-------------------------------------------------

Old description:

> scipy.inalg.lu_factor returns different results when called with a C or
> Fortran contiguous array. The behavior can easily be seen for the
> following small example:
>
> import numpy as np[[BR]]
> import scipy.linalg as la[[BR]]
>

> A=np.array([[2,1],[0,1]])[[BR]]
> la.lu_factor(A)[[BR]]
> la.lu_factor(np.asfortranarray(A))
>
> which for me yields the output
>
> (array([[ 2. ,  0.5],[[BR]]
>        [ 0. ,  1. ]]), array([0, 1], dtype=int32))
> (array([[ 2.,  1.],[[BR]]
>        [ 0.,  1.]]), array([0, 1], dtype=int32))
>
> The matrices lu are obviously different, and only the second output (the
> one from la.lu_factor(np.asfortranarray(A))) are of the form as described
> by the documentation of lu_factor, with unit triangular L in the strictly
> lower and U in the upper triangle of lu.
>
> I suspect that lu_factor acutally factors the transposed system when a C
> contiguous array is used.
>
> In fact, the error cancels out when lu_solve is called:
>
> lu = la.lu_factor(A)[[BR]]
> la.lu_solve(lu, [1,2])[[BR]]
> lu = la.lu_factor(np.asfortranarray(A))[[BR]]
> la.lu_solve(lu, [1,2])[[BR]]
>
> give identical results:
>
> array([-0.5,  2. ])[[BR]]
> array([-0.5,  2. ])
>
> I suspect again this is because lu_solve actually solves the transposed
> system, if it gets as an input a C contiguous LU factorization. It is
> however easy to mess up lu_solve, too:
>
> lu = la.lu_factor(A)[[BR]]
> la.lu_solve((np.asfortranarray(lu[0]), lu[1]), [1,2])
>
> gives the wrong result
>
> array([ 0.,  2.])
>
> It seems to me that lu_factor and lu_solve use the way the array is
> stored as a flag in whether to internally transpoe the system or not.
>
> The function lu() is not affected by this problem.
>
> In my opinion this behavior is not appropriate: Results shouldn't depend
> on what kind of internal memory layout one has.
>
> I see two ways to fix this: Either always return the same result - this
> will require internally converting C to Fortran contiguous. Or change the
> documentation of lu_factor and lu_solve to indicate that the return value
> of lu_factor may have different meanings depending on the memory layout
> of the input, and that the return value of lu_factor should not be
> touched before passing it to lu_solve.

New description:

 scipy.inalg.lu_factor returns different results when called with a C or
 Fortran contiguous array. The behavior can easily be seen for the
 following small example:
 {{{
 import numpy as np
 import scipy.linalg as la


 A=np.array([[2,1],[0,1]])
 la.lu_factor(A)
 la.lu_factor(np.asfortranarray(A))
 }}}
 which for me yields the output
 {{{
 (array([[ 2. ,  0.5],
        [ 0. ,  1. ]]), array([0, 1], dtype=int32))
 (array([[ 2.,  1.],
        [ 0.,  1.]]), array([0, 1], dtype=int32))
 }}}
 The matrices lu are obviously different, and only the second output (the
 one from la.lu_factor(np.asfortranarray(A))) are of the form as described
 by the documentation of lu_factor, with unit triangular L in the strictly
 lower and U in the upper triangle of lu.

 I suspect that lu_factor acutally factors the transposed system when a C
 contiguous array is used.

 In fact, the error cancels out when lu_solve is called:
 {{{
 lu = la.lu_factor(A)
 la.lu_solve(lu, [1,2])
 lu = la.lu_factor(np.asfortranarray(A))
 la.lu_solve(lu, [1,2])
 }}}
 give identical results:
 {{{
 array([-0.5,  2. ])
 array([-0.5,  2. ])
 }}}
 I suspect again this is because lu_solve actually solves the transposed
 system, if it gets as an input a C contiguous LU factorization. It is
 however easy to mess up lu_solve, too:
 {{{
 lu = la.lu_factor(A)
 la.lu_solve((np.asfortranarray(lu[0]), lu[1]), [1,2])
 }}}
 gives the wrong result
 {{{
 array([ 0.,  2.])
 }}}
 It seems to me that lu_factor and lu_solve use the way the array is stored
 as a flag in whether to internally transpoe the system or not.

 The function lu() is not affected by this problem.

 In my opinion this behavior is not appropriate: Results shouldn't depend
 on what kind of internal memory layout one has.

 I see two ways to fix this: Either always return the same result - this
 will require internally converting C to Fortran contiguous. Or change the
 documentation of lu_factor and lu_solve to indicate that the return value
 of lu_factor may have different meanings depending on the memory layout of
 the input, and that the return value of lu_factor should not be touched
 before passing it to lu_solve.

--

Comment(by warren.weckesser):

 Tweak markup in the description.

-- 
Ticket URL: <http://projects.scipy.org/scipy/ticket/1458#comment:1>
SciPy <http://www.scipy.org>
SciPy is open-source software for mathematics, science, and engineering.


More information about the Scipy-tickets mailing list