Sci.linalg.lu permuation error

[bernhard.voigt at gmail.com]

Hi Luke,

you should send this to the scipy user list: scipy-user at scipy.org

Bernhard


On May 28, 10:44 am, Luke <hazelnu... at gmail.com> wrote:
> I'm trying to use Scipy's LU factorization.  Here is what I've got:
>
> from numpy import *
> import scipy as Sci
> import scipy.linalg
>
> A=array([[3., -2., 1., 0., 0.],[-1., 1., 0., 1., 0.],[4., 1., 0., 0.,
> 1.]])
> p,l,u=Sci.linalg.lu(A,permute_l = 0)
>
> now, according to the documentation:
> **************************************************
> Definition:     Sci.linalg.lu(a, permute_l=0, overwrite_a=0)
> Docstring:
>     Return LU decompostion of a matrix.
>
> Inputs:
>
>   a     -- An M x N matrix.
>   permute_l  -- Perform matrix multiplication p * l [disabled].
>
> Outputs:
>
>   p,l,u  -- LU decomposition matrices of a [permute_l=0]
>   pl,u   -- LU decomposition matrices of a [permute_l=1]
>
> Definitions:
>
>   a = p * l * u    [permute_l=0]
>   a = pl * u       [permute_l=1]
>
>   p   -  An M x M permutation matrix
>   l   -  An M x K lower triangular or trapezoidal matrix
>          with unit-diagonal
>   u   -  An K x N upper triangular or trapezoidal matrix
>          K = min(M,N)
> *********************************************
>
> So it would seem that from my above commands, I should have:
> A == dot(p,dot(l,u))
>
> but instead, this results in:
> In [21]: dot(p,dot(l,u))
> Out[21]:
> array([[-1.,  1.,  0.,  1.,  0.],
>        [ 4.,  1.,  0.,  0.,  1.],
>        [ 3., -2.,  1.,  0.,  0.]])
>
> Which isn't equal to A!!!
> In [23]: A
> Out[23]:
> array([[ 3., -2.,  1.,  0.,  0.],
>        [-1.,  1.,  0.,  1.,  0.],
>        [ 4.,  1.,  0.,  0.,  1.]])
>
> However, if I do use p.T instead of p:
> In [22]: dot(p.T,dot(l,u))
> Out[22]:
> array([[ 3., -2.,  1.,  0.,  0.],
>        [-1.,  1.,  0.,  1.,  0.],
>        [ 4.,  1.,  0.,  0.,  1.]])
>
> I get the correct answer.
>
> Either the documentation is wrong, or somehow Scipy is returning the
> wrong permutation matrix... anybody have any experience with this or
> tell me how to submit a bug report?
>
> Thanks.
> ~Luke