Hi Luke, you should send this to the scipy user list: [EMAIL PROTECTED]
Bernhard On May 28, 10:44 am, Luke <[EMAIL PROTECTED]> 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 -- http://mail.python.org/mailman/listinfo/python-list