[SciPy-User] Permutation convention for LU decomposition
Charles R Harris
charlesr.harris@gmail....
Sun Oct 31 13:13:11 CDT 2010
On Sun, Oct 31, 2010 at 6:18 AM, Sebastian Walter <
sebastian.walter@gmail.com> wrote:
> 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.
>
>
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
More information about the SciPy-User
mailing list