473,544 Members | 2,340 Online
Bytes | Software Development & Data Engineering Community
+ Post

Home Posts Topics Members FAQ

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. :( <?php session_start (); $username = $_POST;
0
7781
jinu1996
by: jinu1996 | last post by:
In today's digital age, having a compelling online presence is paramount for businesses aiming to thrive in a competitive landscape. At the heart of this digital strategy lies an intricately woven tapestry of website design and digital marketing. It's not merely about having a website; it's about crafting an immersive digital experience that...
0
7717
tracyyun
by: tracyyun | last post by:
Dear forum friends, With the development of smart home technology, a variety of wireless communication protocols have appeared on the market, such as Zigbee, Z-Wave, Wi-Fi, Bluetooth, etc. Each protocol has its own unique characteristics and advantages, but as a user who is planning to build a smart home system, I am a bit confused by the...
0
5928
agi2029
by: agi2029 | last post by:
Let's talk about the concept of autonomous AI software engineers and no-code agents. These AIs are designed to manage the entire lifecycle of a software development project—planning, coding, testing, and deployment—without human intervention. Imagine an AI that can take a project description, break it down, write the code, debug it, and then...
1
5306
isladogs
by: isladogs | last post by:
The next Access Europe User Group meeting will be on Wednesday 1 May 2024 starting at 18:00 UK time (6PM UTC+1) and finishing by 19:30 (7.30PM). In this session, we are pleased to welcome a new presenter, Adolph Dupré who will be discussing some powerful techniques for using class modules. He will explain when you may want to use classes...
0
3427
by: TSSRALBI | last post by:
Hello I'm a network technician in training and I need your help. I am currently learning how to create and manage the different types of VPNs and I have a question about LAN-to-LAN VPNs. The last exercise I practiced was to create a LAN-to-LAN VPN between two Pfsense firewalls, by using IPSEC protocols. I succeeded, with both firewalls in...
0
3421
by: adsilva | last post by:
A Windows Forms form does not have the event Unload, like VB6. What one acts like?
1
1848
by: 6302768590 | last post by:
Hai team i want code for transfer the data from one system to another through IP address by using C# our system has to for every 5mins then we have to update the data what the data is updated we have to send another system
1
993
muto222
by: muto222 | last post by:
How can i add a mobile payment intergratation into php mysql website.
0
677
bsmnconsultancy
by: bsmnconsultancy | last post by:
In today's digital era, a well-designed website is crucial for businesses looking to succeed. Whether you're a small business owner or a large corporation in Toronto, having a strong online presence can significantly impact your brand's success. BSMN Consultancy, a leader in Website Development in Toronto offers valuable insights into creating...

By using Bytes.com and it's services, you agree to our Privacy Policy and Terms of Use.

To disable or enable advertisements and analytics tracking please visit the manage ads & tracking page.