These are the outputs for my run (adding a print "perl version $]";print "PDL version $PDL::VERSION"; at the beginning to print the versions):

Mathematically speaking, a [3x3] x [3x2] is the right balance of dimensions for a successful matrix multiplication (resulting in a [3x2] matrix). The use of PDL::Matrix used to be the way to make PDL match mathematical order for the dimensions, but maybe they've changed that between 2.019 and 2.057. (I don't think I can test, because the Strawberry Perl PDL setup is rather dependent on the libraries that Strawberry bundles, and I don't know if upgrading PDL will work... maybe I'll try in a sandbox to avoid ruining my "main" installation).

I tried getting rid of the use PDL::Matrix, but then that gets rid of the mpdl command, causing errors during perl's compile stage, and I'm not sure what hoops I would need to jump through to make my script compatible with both 2.019 and 2.057.

Okay, 5.30.2 & PDL 2.021 and 5.32.1 & PDL 2.025 (the last pre-bundled Perl+PDL from Strawberry) both worked like my 5.30.0 & PDL 2.019.

But I was able to install PDL 5.080 with Perl 5.32.1, to get slightly newer PDL but slightly older Perl than you. But it showed the same dimensional difference.

So I did some debug to make it work without PDL::Matrix

#!perl -l

use 5.012; # //, strict, say
use warnings;

use PDL 2.057;
#use PDL::Matrix;

print "perl version $]";
print "PDL version $PDL::VERSION";

print my $red_in  = (pdl [ # reference
    [419 , 17 ], #  p01
    [772 , 13 ], #  p02
    [348 , 102], #   p03
    [720 , 95 ], #  p04
    [58  , 238], #   p05
    [218 , 236], #   p06
    [55  , 333], #   p07
    [58  , 391], #   p08
    [60  , 498], #   p09
    [58  , 554], #   p10
    [1061, 553], #   p11
    [60  , 892], #   p12
])->append(ones(1))->transpose;
print my $red_out = (pdl [ # target_41
    [487, 466], # p01
    [531, 630], # p02
    [439, 447], # p03
    [486, 618], # p04
    [340, 328], # p05
    [358, 401], # p06
    [297, 343], # p07
    [274, 349], # p08
    [227, 366], # p09
    [205, 374], # p10
    [313, 809], # p11
    [77 , 420], # p12
])->append(ones(1))->transpose;

print my $blue_in = (pdl [ # reference
    [438,331],  # q01
    [432,392],  # q02
    [436,493],  # q03
    [436,552],  # q04
    [439,602],  # q05
])->append(ones(1))->transpose;
print my $blue_out_for_checking = (pdl [ # target_41
    [339,508],  # q01
    [314,518],  # q02
    [265,533],  # q03
    [245,537],  # q04
    [228,544],  # q05
])->append(ones(1))->transpose;

$\ = '';
print for
    "sx  = ", my $sx  = sum( $red_in->transpose->slice(0) ),          
+          "\n",
    "sy  = ", my $sy  = sum( $red_in->transpose->slice(1) ),          
+          "\n",
    "s1  = ", my $s1  = sum( $red_in->transpose->slice(2) ),          
+          "\n",
    "sxx = ", my $sxx = sum( $red_in->transpose->slice(0)**2 ),       
+          "\n",
    "syy = ", my $syy = sum( $red_in->transpose->slice(1)**2 ),       
+          "\n",
    "sxy = ", my $sxy = sum( $red_in->transpose->slice(0)*$red_in->tra
+nspose->slice(1) ),  "\n",
    "sv  = ", my $sv  = sum( $red_out->transpose->slice(0) ),         
+          "\n",
    "sw  = ", my $sw  = sum( $red_out->transpose->slice(1) ),         
+          "\n",
    "sxv = ", my $sxv = sum( $red_in->transpose->slice(0)*$red_out->tr
+anspose->slice(0) ), "\n",
    "sxw = ", my $sxw = sum( $red_in->transpose->slice(0)*$red_out->tr
+anspose->slice(1) ), "\n",
    "syv = ", my $syv = sum( $red_in->transpose->slice(1)*$red_out->tr
+anspose->slice(0) ), "\n",
    "syw = ", my $syw = sum( $red_in->transpose->slice(1)*$red_out->tr
+anspose->slice(1) ), "\n",
    ;

$\ = "\n";

print "G => ", my $G = pdl
[
    [ $sxx, $sxy, $sx ],
    [ $sxy, $syy, $sy ],
    [ $sx,  $sy,  $s1 ],
];

print "Q => ", my $Q = pdl
[
    [ $sxv, $sxw ],
    [ $syv, $syw ],
    [ $sv,  $sw  ],
];


print "P => ", my $P = inv($G) x $Q;

# print "P->transpose", $P->transpose;
# print "P->append->transpose", $P->append(zeros(1))->transpose;

print "abcdef001 => ", my $abcdef001 = $P->append(zeros(1))->transpose
+->set(2,2,1);

print "blue out calc => ", my $blue_out_calc = $abcdef001 x $blue_in;

print "error => ", $blue_out_calc - $blue_out_for_checking; # no more 
+than one pixel off

__END__

perl version 5.032001
PDL version 2.080

[
 [ 419  772  348  720   58  218   55   58   60   58 1061   60]
 [  17   13  102   95  238  236  333  391  498  554  553  892]
 [   1    1    1    1    1    1    1    1    1    1    1    1]
]


[
 [487 531 439 486 340 358 297 274 227 205 313  77]
 [466 630 447 618 328 401 343 349 366 374 809 420]
 [  1   1   1   1   1   1   1   1   1   1   1   1]
]


[
 [438 432 436 436 439]
 [331 392 493 552 602]
 [  1   1   1   1   1]
]


[
 [339 314 265 245 228]
 [508 518 533 537 544]
 [  1   1   1   1   1]
]

sx  = 3887
sy  = 3922
s1  = 12
sxx = 2604611
syy = 2052390
sxy = 929565
sv  = 4034
sw  = 5551
sxv = 1608891
sxw = 2354880
syv = 945962
syw = 1755275
G =>
[
 [2604611  929565    3887]
 [ 929565 2052390    3922]
 [   3887    3922      12]
]

Q =>
[
 [1608891 2354880]
 [ 945962 1755275]
 [   4034    5551]
]

P =>
[
 [0.11504241 0.44420653]
 [-0.4325134 0.11994661]
 [ 440.26231  279.49489]
]

abcdef001 =>
[
 [0.11504241 -0.4325134  440.26231]
 [0.44420653 0.11994661  279.49489]
 [         0          0          1]
]

blue out calc =>
[
 [ 347.48895  320.41538  277.19169   251.6734  230.39286]
 [ 513.75967  518.41118  532.30261  539.37946  546.70941]
 [         1          1          1          1          1]
]

error =>
[
 [   8.488949   6.4153774   12.191694   6.6734037   2.3928612]
 [  5.7596723  0.41117608 -0.69739064   2.3794591    2.709409]
 [          0           0           0           0           0]
]
[download]

This shows the 12-ish pixel error, but will hopefully work with your 2.057 (and definitely worked with my 2.080).

Minimal pair:

perl -MPDL -MPDL::Matrix -e '$m = mpdl [[2,0],[0,2]]; $v = vpdl [1,2];
+ print $m x $v'

[
 [2]
 [4]
]
[download]

vs.

perl -MPDL -MPDL::Matrix -e '$m = mpdl [[2,0],[0,2]]; $v = vpdl [1,2];
+ print inv($m) x $v'
Dim mismatch in matmult of [2x2] x [2x1]: 2 != 1 at /usr/lib64/perl5/v
+endor_perl/PDL/Primitive.pm line 274.
    PDL::matmult(PDL=SCALAR(0x5600a82130e0), PDL::Matrix=SCALAR(0x5600
+a8202ba0), PDL=SCALAR(0x5600a8d718c8)) called at /usr/lib64/perl5/ven
+dor_perl/PDL/Primitive.pm line 31
    PDL::__ANON__(PDL=SCALAR(0x5600a82130e0), PDL::Matrix=SCALAR(0x560
+0a8202ba0), "") called at -e line 1
[download]

This looks like a bug in PDL::Matrix, specifically with inv.

As now noted on pryrt's excellent bug report (thank you), this will still work as expected if you switch your inv($m) to $m->inv. This will then engage PDL::Matrix's overload of x, which is what you need here.

I would strongly encourage you not to use PDL::Matrix, because somewhat similar to PDL::Complex, it is full of gotchas. Even though after some head-scratching it became clear what exactly was happening here, it's not particularly clear how to fix it - a keyhole change to the Perl function PDL::MatrixOps::inv to make it ensure the output was the right class would resolve this exact issue. But since that uses lu_backsub for its work, and that is an eye-wateringly complicated Perl function as well, it will not be easy. That function's output variable comes from some PDL::Slices operations.

If someone were minded to at least use git bisect to identify the commit that broke it, that would be very helpful. Even better would be a pull-request that fixes this (obviously with a test, in t/matrix.t), that will use our amazing CI to check it doesn't break anything else in main PDL or known downstreams. A top tip would be that make coretest runs t/matrix.t, with a very minimal build otherwise that allows a very quick dev cycle.

As a hopefully practical (as well as better-performance, because it uses LAPACK) alternative, have you taken a look at PDL::LinearAlgebra?

kikichiyo,

Thanks for that SSCCE.

I included it in a bug report

The way to avoid needing mpdl is to transpose your data, and use pdl, and I would recommend doing so. But as observed elsewhere, there are almost certainly out-of-the-box solutions to this particular problem.