Sci.linalg.lu permuation error
Luke
hazelnusse at gmail.com
Mon May 28 04:44:59 EDT 2007
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
More information about the Python-list
mailing list