Home     |     .Net Programming    |     cSharp Home    |     Sql Server Home    |     Javascript / Client Side Development     |     Ajax Programming

Ruby on Rails Development     |     Perl Programming     |     C Programming Language     |     C++ Programming     |     IT Jobs

Python Programming Language     |     Laptop Suggestions?    |     TCL Scripting     |     Fortran Programming     |     Scheme Programming Language


 
 
Cervo Technologies
The Right Source to Outsource

MS Dynamics CRM 3.0

Python Programming Language

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

Hi Luke,

you should send this to the scipy user list: scipy-u@scipy.org

Bernhard

On May 28, 10:44 am, Luke <hazelnu@gmail.com> wrote:

Add to del.icio.us | Digg this | Stumble it | Powered by Megasolutions Inc