#206 failure in MatrixOps: decomp & backsub

critical
closed
nobody
other (94)
5
2008-10-28
2008-07-21
W_Luis
No

The following programs fails to solve a simple 2x2 system of equations; I believe something is wrong with either lu_decomp and/or lu_backsub from PDL::MatrixOps. (or am I making a something wrong?). Maybe it is related with amd64.

Program:

#!/usr/bin/perl -w
use strict;
use PDL;
use PDL::NiceSlice;
my $m=pdl([[2,1],[1,2]]); #make a simple matrix
my ($lu,$perm,$par)=lu_decomp($m); #decompose
my $x=pdl([[1],[0]]); #make a simple column vector
#back substitute directly and transpossing
my $y=lu_backsub($lu,$perm,$par,$x);
my $y1=lu_backsub($lu,$perm,$par,$x->transpose)->transpose;
print "Solve \$m x \$y= \$x\n";
print "Original matrix times result: ", $m x $y;
print "Original matrix times alternative result: ", $m x $y1;
print "One of the results should have yielded the original \$x: $x";

Output:

Solve $m x $y= $x
Original matrix times result:
[
[ 1]
[0.5]
]
Original matrix times alternative result:
[
[ 0.38888889]
[-0.13888889]
]
One of the results should have yielded the original $x:
[
[1]
[0]
]

The output from perldl -V is attached

Discussion

1 2 > >> (Page 1 of 2)
  • W_Luis
    W_Luis
    2008-07-21

    output from perldl -V

     
    Attachments
  • Chris Marshall
    Chris Marshall
    2008-07-27

    Logged In: YES
    user_id=44920
    Originator: NO

    The signature of the lu_backsub() routine is:

    (lu(m,m); perm(m); b(m))

    Which means that the calling args for the two tests
    should be

    (lu(2,2); perm(2); x(1,2))

    and

    (lu(2,2); perm(2); x(2,1))

    neither of which matches the required signature.
    I tried using $x = pdl([1,0]) instead with call
    of

    (lu(2,2); perm(2); x(2))

    which does match and got for the result

    perldl> $y = lu_backsub($lu,$perm,$d,$x)

    perldl> p $y

    [
    [ 0.66666667 -0.33333333]
    [ 0.66666667 -0.33333333]
    ]

    perldl> p $y x $m

    [
    [1 0]
    [1 0]
    ]

    which does give the correct solution
    although it looks like there is a problem
    with the threading.

    I tried padding the $lu and $perm
    arrays with extra singleton dimension
    which should trivially thread and got
    the incorrect answer in case #2 above.

    When I used

    $lu = $lu->dummy(-1,2)-sever
    $perm = $perm->dummy(-1,2)->sever
    $x = $x->dummy(-1,2)->sever

    then I got the repeated "correct"
    answer as above (i.e. it has two
    identical solutions but the solutions
    are 2x2 rather than 2x1).

    SUMMARY: the lu_decomp() appears to work correctly.
    The lu_backsub() gives correct (but redundant) solns
    if the arguments have the correct dimensionality.
    There may be a problem with trivial threading or
    the threading in general. Needs more looking at...

     
  • Chris Marshall
    Chris Marshall
    2008-07-27

    lu_backsub() pp implementation and test

     
    Attachments
  • Chris Marshall
    Chris Marshall
    2008-07-27

    Logged In: YES
    user_id=44920
    Originator: NO

    I implemented a PDL::PP version of the NR in C back substitution code (see attached Inline::Pdlpp code) which appeared to work for the 2x2 problem of this bug report. When I tried it on a problem with a nontrivial permutation, it does not work. No more time to work on this but maybe this can help someone finish this off. The PP version of the code, lu_backsub2(), seems to handle the loop/thread dimensions but something is breaking when a permutation is obtained in the lu_decomp().
    File Added: lu_back.pl

     
  • Chris Marshall
    Chris Marshall
    2008-07-29

    Logged In: YES
    user_id=44920
    Originator: NO

    I took another look at the NR C code I translated into a PP lu_backup2() routine. The LU decomp still looks solid. The problem remaining appears to be the back substitution. I think it may be tied into row/col-major issues and dimension ordering---something that can be tricky to keep straight at times.

     
  • Chris Marshall
    Chris Marshall
    2008-09-18

    Another alternative here would be to implement the GSL versions of these functions. That might be a better alternative in the long run.

     
  • W_Luis
    W_Luis
    2008-09-18

    I actually solved my problem using the Slatec routines (geco and gesl), but the bug in MatrixOps still worries me.

     
  • Chris Marshall
    Chris Marshall
    2008-10-06

    This needs to be resolved for PDL-2.4.4 if possible.

     
  • Chris Marshall
    Chris Marshall
    2008-10-06

    • priority: 5 --> 7
     
  • Chris Marshall
    Chris Marshall
    2008-10-13

    Bug fixed in CVS.
    Thanks for reporting the problem!

     
1 2 > >> (Page 1 of 2)