[SciPy-User] Permutation convention for LU decomposition

[Charles R Harris]

On Sun, Oct 31, 2010 at 6:18 AM, Sebastian Walter <
sebastian.walter at gmail.com> wrote:

> On Sun, Oct 31, 2010 at 2:42 AM, Jason Grout
> <jason-sage at 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.
>
>
The matrix form of P could be a problem with really big matrices. Same with
the orthogonal matrices in SVD and QR, which are stored internally as
Householder reflections. I don't see an easy way around this without
complicating matters at the user level.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20101031/f7eb3190/attachment.html>