473,544 Members | 2,340 Online

# Sci.linalg.lu permuation error

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.linal g.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

May 28 '07 #1
1 1620
Hi Luke,

you should send this to the scipy user list: sc********@scip y.org

Bernhard
On May 28, 10:44 am, Luke <hazelnu...@gma il.comwrote:
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.linal g.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

May 29 '07 #2

This thread has been closed and replies have been disabled. Please start a new discussion.

### Similar topics

 2 4363 by: AIM | last post by: Error in msvc in building inheritance.obj to build hello.pyd Hello, I am trying to build the boost 1.31.0 sample extension hello.cpp. I can not compile the file inheritance.cpp because the two files containing some templates: adjacency_list.hpp and mem_fn.hpp can not compile. Does anyone have any solutions? 2 9726 by: I. Myself | last post by: And it has to run on Windows, so it can't use xplt. I would prefer that it use the simplest multi-dimensional model, z = k + a*x1 + b*x2 + c*x3 + d*x4 Anyone have such a thing? Thanks, Mitchell Timin 3 1725 by: Sophiecarr | last post by: Hi, I am trying to use DLookup in relation to input from an option group. The option group has 3 options (Course Title), in a 'Core Modules' table, I have set the primary field to be the 'Course Title' and have two further fields containing core modules. I am aiming to update a text box with these core modules once the user clicks on their... 5 1640 by: punk86 | last post by: I have this login codings. which it have a parsing error at line 9 which starts from the \$query. can anyone help? i have problems with parsing errors. :(