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 unitdiagonal
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 1 1520
Hi Luke,
you should send this to the scipy user list: sc********@scipy.org
Bernhard
On May 28, 10:44 am, Luke <hazelnu...@gmail.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.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 unitdiagonal
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
This discussion thread is closed Replies have been disabled for this discussion. Similar topics
2 posts
views
Thread by AIM 
last post: by

2 posts
views
Thread by I. Myself 
last post: by
            