[SciPy-User] Permutation convention for LU decomposition

Sebastian Walter sebastian.walter@gmail....
Sun Oct 31 07:18:32 CDT 2010


On Sun, Oct 31, 2010 at 2:42 AM, Jason Grout
<jason-sage@creativetrax.com> wrote:
> I notice that in Lapack, Matlab, and Mathematica, the LU decomposition
> routine for a matrix A returns a P, L, and U matrices so that:
>
> PA=LU
>
> In scipy, however, the LU decomposition routine gives three matrices so
> that:
>
> A=PLU
>
> (i.e., the P matrix is the inverse of the P matrix returned by the other
> software)

P should be orthogonal and thus P^{-1} = P^T
so PA = LU or A = P L U shouldn't be a big issue.

>
> I'm curious why this design decision was made.  Was there precedent with
> other software to make it A=PLU instead of PA=LU?  Is it more natural in
> applications?  I realize it's just a convention.

Algorithms like dgetrf (in LAPACK)
do not return a matrix P but an IPIV array that stores _successive_
row interchanges.
E.g.:

In [1]: import numpy

In [2]: import scipy.linalg

In [3]: A = numpy.random.rand(3,3)

In [4]: print scipy.linalg.lapack.clapack.dgetrf(A.copy())
(array([[ 0.73664967,  0.2875308 ,  0.67223657],
       [ 0.82181513,  0.62191285, -0.00461938],
       [ 0.13967213,  0.42672366, -0.02721027]]), array([1, 1, 2],
dtype=int32), 0)

In [5]: print scipy.linalg.lu(A.copy())
(array([[ 0.,  1.,  0.],
       [ 1.,  0.,  0.],
       [ 0.,  0.,  1.]]), array([[ 1.        ,  0.        ,  0.        ],
       [ 0.82181513,  1.        ,  0.        ],
       [ 0.13967213,  0.42672366,  1.        ]]), array([[ 0.73664967,
 0.2875308 ,  0.67223657],
       [ 0.        ,  0.62191285, -0.00461938],
       [ 0.        ,  0.        , -0.02721027]]))


The array([1, 1, 2], dtype=int32) are the pivots, typically denoted
IPIV in LAPACK. I.e., you first interchange row 0 <-> 1, then  1 <-> 1
and finally 2<->2. In total you have interchanged
row 0 and 1. Row 2 remains unchanged.

So, arguably the natural representation of the permutation is IPIV and
not the permutation matrix P.
I'm not familiar with the  scipy.linalg.lu implementation. If it is
based on dgetref it probably computes P from IPIV.


Sebastian


>
> Thanks,
>
> Jason
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>


More information about the SciPy-User mailing list