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