|
choroba has asked for the wisdom of the Perl Monks concerning the following question:
Spoiler alert: If you participate in The Weekly Challenge, don't read further if you haven't solved week 170 yet.
Kronecker product is a matrix operation that uses elements of one matrix to multiply the second matrix. Mohammad shows this example:
A = [ 1 2 ]
[ 3 4 ]
B = [ 5 6 ]
[ 7 8 ]
A x B = [ 1 x [ 5 6 ] 2 x [ 5 6 ] ]
[ [ 7 8 ] [ 7 8 ] ]
[ 3 x [ 5 6 ] 4 x [ 5 6 ] ]
[ [ 7 8 ] [ 7 8 ] ]
= [ 1x5 1x6 2x5 2x6 ]
[ 1x7 1x8 2x7 2x8 ]
[ 3x5 3x6 4x5 4x6 ]
[ 3x7 3x8 4x7 4x8 ]
= [ 5 6 10 12 ]
[ 7 8 14 16 ]
[ 15 18 20 24 ]
[ 21 24 28 32 ]
When I saw matrices, I immediately thought PDL. After finding a solution, I tried searching for existing solutions, and found the following at Rosetta Code:
#!/usr/bin/perl
use strict;
use warnings;
use PDL;
use PDL::NiceSlice;
sub kron{
my $A = shift;
my $B = shift;
my ($r0, $c0) = $A->dims;
my ($r1, $c1) = $B->dims;
my $kron = zeroes($r0 * $r1, $c0 * $c1);
for(my $i = 0; $i < $r0; ++$i){
for(my $j = 0; $j < $c0; ++$j){
$kron(
($i * $r1) : (($i + 1) * $r1 - 1),
($j * $c1) : (($j + 1) * $c1 - 1)
) .= $A($i,$j) * $B;
}
}
return $kron;
}
Wait, loops? The whole point of PDL is to hide loops. And indeed, my solution doesn't involve them. Also, a benchmark shows my solution is more than twice faster than the Rosetta Code one.
I'm far from an expert on PDL or matrices in general. But it seems we can easily multiply each element in a matrix by the same number, but it's not so easy to multiply them by different numbers. But we can multiply each element by a matrix, so we just need to "inflate" the matrix, so instead of
[ 1 2 ]
[ 3 4 ]
we'd have
[ 1 1 ] [ 2 2 ]
[ 1 1 ] [ 2 2 ]
[ 3 3 ] [ 4 4 ]
[ 3 3 ] [ 4 4 ]
That's what dummy does (two dummies, in fact, one in each dimension). Multiplying this with the second matrix gives us a result that has all the expected numbers, but a bit different dimensions.
[ 5 6 ] [ 10 12 ]
[ 7 8 ] [ 14 16 ]
[ 15 18 ] [ 20 24 ]
[ 21 24 ] [ 28 32 ]
Fortunately, PDL has all the needed methods to reshuffle the submatrices into the expected result.
sub kronecker_product {
my ($x, $y) = @_;
my ($x0, $x1) = $x->dims;
my ($y0, $y1) = $y->dims;
return (
$y * $x->dummy(0, $y0)->dummy(1, $y1)
)->xchg(1, 2)->reshape($x0 * $y0, $x1 * $y1)
}
It seems correct (it passes the test shown in the challenge and in the Wikipedia page), but my PDL is not so strong. Maybe it can be further improved?
map{substr$_->[0],$_->[1]||0,1}[\*||{},3],[[]],[ref qr-1,-,-1],[{}],[sub{}^*ARGV,3]
|
Re: Kronecker Product
by jwkrahn (Abbot) on Jun 21, 2022 at 08:12 UTC
|
$ perl -le'
my @A = ( [ 1, 2 ], [ 3, 4 ] );
print map "[ @$_ ]\n", @A;
my @B = ( [ 5, 6 ], [ 7, 8 ] );
print map "[ @$_ ]\n", @B;
my @A_B;
for my $x ( @A ) {
for my $y ( @B ) {
push @A_B, [ map( $x->[ 0 ] * $_, @$y ), map $x->[ 1 ] * $_, @
+$y ];
}
}
print map "[ @$_ ]\n", @A_B;
'
[ 1 2 ]
[ 3 4 ]
[ 5 6 ]
[ 7 8 ]
[ 5 6 10 12 ]
[ 7 8 14 16 ]
[ 15 18 20 24 ]
[ 21 24 28 32 ]
| [reply]
[d/l] |
|
|
It doesn't work for the larger example from Wikipedia:
my $X = pdl([1, -4, 7],
[-2, 3, 3]);
my $Y = pdl([8, -9, -6, 5],
[1, -3, -4, 7],
[2, 8, -8, -3],
[1, 2, -5, -1]);
my $XY =pdl([ 8, -9, -6, 5, -32, 36, 24, -20, 56, -63, -42, 35]
+,
[ 1, -3, -4, 7, -4, 12, 16, -28, 7, -21, -28, 49]
+,
[ 2, 8, -8, -3, -8, -32, 32, 12, 14, 56, -56, -21]
+,
[ 1, 2, -5, -1, -4, -8, 20, 4, 7, 14, -35, -7]
+,
[-16, 18, 12, -10, 24, -27, -18, 15, 24, -27, -18, 15]
+,
[ -2, 6, 8, -14, 3, -9, -12, 21, 3, -9, -12, 21]
+,
[ -4, -16, 16, 6, 6, 24, -24, -9, 6, 24, -24, -9]
+,
[ -2, -4, 10, 2, 3, 6, -15, -3, 3, 6, -15, -3]
+);
map{substr$_->[0],$_->[1]||0,1}[\*||{},3],[[]],[ref qr-1,-,-1],[{}],[sub{}^*ARGV,3]
| [reply]
[d/l]
[select] |
|
|
that's because he hardcoded $x->[0] and $x->[1]
here a generic version
use v5.12;
use warnings;
use Data::Dump qw/pp dd/;
my $X =
[
[1, -4, 7],
[-2, 3, 3]
];
my $Y =
[
[8, -9, -6, 5],
[1, -3, -4, 7],
[2, 8, -8, -3],
[1, 2, -5, -1]
];
pp $X;
pp $Y;
my $X_Y;
for my $x ( @$X ) {
for my $y ( @$Y ) {
push @$X_Y,
[
map {
my $xx = $_;
map { $xx * $_} @$y
} @$x
];
}
}
pp $X_Y;
[[1, -4, 7], [-2, 3, 3]]
[[8, -9, -6, 5], [1, -3, -4, 7], [2, 8, -8, -3], [1, 2, -5, -1]]
[
[8, -9, -6, 5, -32, 36, 24, -20, 56, -63, -42, 35],
[1, -3, -4, 7, -4, 12, 16, -28, 7, -21, -28, 49],
[2, 8, -8, -3, -8, -32, 32, 12, 14, 56, -56, -21],
[1, 2, -5, -1, -4, -8, 20, 4, 7, 14, -35, -7],
[-16, 18, 12, -10, 24, -27, -18, 15, 24, -27, -18, 15],
[-2, 6, 8, -14, 3, -9, -12, 21, 3, -9, -12, 21],
[-4, -16, 16, 6, 6, 24, -24, -9, 6, 24, -24, -9],
[-2, -4, 10, 2, 3, 6, -15, -3, 3, 6, -15, -3],
]
updates
cleaned up code
addded output | [reply]
[d/l]
[select] |
|
|
|
|
Re: Kronecker Product
by vr (Curate) on Jun 21, 2022 at 12:44 UTC
|
I have a nitpick and an improvement (well, in my view). The latter: the beauty of "array language" is that it allows one to be succinct. No need to bother with boring details of array shape (i.e. dimensions) explicitly: interpreter ain't stupid, it itself can figure out the info it needs. Size of leading dummy dimensions of $x, if "1" originally, will be stretched as required to fit shape of second array i.e. multiplicand i.e. $y. (Shape of which, in its turn, will be padded with dummy dimensions. So then both multiplicands have matching shapes of 4 dimensions.) Further, oh-please spare me the details that 2D shape of result will be pair of multiplications of widths and heights of original matrices. Can't tolerate too much maths for today already! I only know that 4D shape of multiplication result should be reduced: dims 0 and 2 clumped together (note implicit xchg (or mv) here) to produce 3D, then dims 1 and 2 also clumped to produce 2D. Will shape then be w1*w2 x h1*h2? Fine, but it's none of my business. So:
sub kronecker_product {
my ($x, $y) = @_;
return $x-> dummy( 0 )
-> dummy( 0 )
-> mult( $y, 0 )
-> clump( 0, 2 )
-> clump( 1, 2 )
}
Or:
use PDL::NiceSlice;
sub kronecker_product {
my ($x, $y) = @_;
( $x( *1, *1 ) * $y )
-> clump( 0, 2 )
-> clump( 1, 2 )
}
A nitpick: operation is non-commutative. Function args are $x, $y, in that order. There was no reason to write computation as $y * $x-> ...etc. Of course scalar multiplication of element by element of 4D arrays is commutative, and so result is correct either way. It's just my preference to keep written expression so it has no unjustified swaps of arguments as $x-> call() * $y. Sorry for nitpick.
| [reply]
[d/l]
[select] |
|
|
Thanks, that's exactly the reply I was waiting for.
I see no difference between your solution and mine in the benchmark, but I definitely understand your points and see the improvement in the code.
Are you aware of any reading on "array languages" that would be understandable to tyros? I understand your solution but I can't imagine how to get to it (to get mine, I used trial and error repeatedly).
map{substr$_->[0],$_->[1]||0,1}[\*||{},3],[[]],[ref qr-1,-,-1],[{}],[sub{}^*ARGV,3]
| [reply]
[d/l] |
|
|
But it's same trial and error on my part. I think if anything can count as systematic approach, as true textbooks, from very basics and up, -- that, for me, have been Learning J and J for C Programmers, books designed to teach "tyros" the "array language", their authors using a bit different perspectives.
| [reply] |
|
|
A nice intro to the concepts of array programming is https://en.wikipedia.org/wiki/Array_programming.
Don't feel bad about trial and error; my first go at this (and I have the dubious privilege of being current PDL maintainer) assumed (wrongly) that dupN would be suitable for both sides of the multiplication; furthermore, it took me a little while to figure out inflateN, even though, annoyingly, it is nearly identical to dupN except for swapping the two parts of the slice argument.
| [reply]
[d/l]
[select] |
|
|
I was inspired to copy the recently (as of 2.077) added dupN to make inflateN, which does the converse of dupN, suitable for Kronecker products.
sub PDL::inflateN {
my ($this, @times) = @_;
return $this->copy if !grep $_ != 1, @times;
my $sl = join ',', map "*$_,:", @times;
$this = $this->slice($sl);
$this = $this->clump($_, $_+1) for 0..$#times;
$this;
}
pdl> p $a = pdl '[[1 2][3 4]]'
[
[1 2]
[3 4]
]
pdl> p $b = pdl '[[5 6][7 8][9 10]]'
[
[ 5 6]
[ 7 8]
[ 9 10]
]
pdl> p $a->inflateN($b->dims)
[
[1 1 2 2]
[1 1 2 2]
[1 1 2 2]
[3 3 4 4]
[3 3 4 4]
[3 3 4 4]
]
pdl> p $b->dupN($a->dims)
[
[ 5 6 5 6]
[ 7 8 7 8]
[ 9 10 9 10]
[ 5 6 5 6]
[ 7 8 7 8]
[ 9 10 9 10]
]
sub kron { my ($x,$y) = @_; $x->inflateN($y->dims) * $y->dupN($x->dims
+) }
pdl> p kron($a,$b)
[
[ 5 6 10 12]
[ 7 8 14 16]
[ 9 10 18 20]
[15 18 20 24]
[21 24 28 32]
[27 30 36 40]
]
The above will broadcast nicely over any number of dimensions, not just 2. inflateN has been added and will be in the next PDL release (2.081). | [reply]
[d/l]
[select] |
|
|
I think dummy dimensions add only a few bytes to PDL-struct header, but your approach really physicalizes arrays from 2D to 4D (edit: well, it's still same 2D but squared dimensions for this example, sorry), which has grave consequences for performance. If no-one ever needs Kronecker product of 80x80 matrices, then count this comment as another nitpick :). + Note report of memory consumption assumes Win32, so disable/replace otherwise.
use strict;
use warnings;
use feature 'say';
use PDL;
use Time::HiRes 'time';
use constant DIM => 80;
my $x = random DIM, DIM;
my $y = random DIM, DIM;
report( 'choroba, vr', \&kronecker_product, $x, $y );
report( 'etj', \&kron, $x, $y );
sub report {
my ( $monk, $code, $x, $y ) = @_;
my $t = time;
my $k = $code-> ( $x, $y );
say "$monk:";
say 'time: ', time - $t;
say 'memory: ', mem();
}
sub mem {
qx{ typeperf "\\Process(perl)\\Working Set Peak" -sc 1 }
=~ /(\d+)\.\d+\"$/m;
( my $s = $1 ) =~ s/(\d{1,3}?)(?=(\d{3})+$)/$1,/g;
return $s
}
sub kronecker_product {
use PDL::NiceSlice;
my ( $x, $y ) = @_;
( $x( *1, *1 ) * $y )
-> clump( 0, 2 )
-> clump( 1, 2 )
}
sub PDL::dupN {
my ($this, @times) = @_;
return $this->copy if !grep $_ != 1, @times;
my $sl = join ',', map ":,*$_", @times;
$this = $this->slice($sl);
$this = $this->clump($_, $_+1) for 0..$#times;
$this;
}
sub PDL::inflateN {
my ($this, @times) = @_;
return $this->copy if !grep $_ != 1, @times;
my $sl = join ',', map "*$_,:", @times;
$this = $this->slice($sl);
$this = $this->clump($_, $_+1) for 0..$#times;
$this;
}
sub kron { my ($x,$y) = @_; $x->inflateN($y->dims) * $y->dupN($x->dims
+) }
__END__
choroba, vr:
time: 0.0958220958709717
memory: 355,741,696
etj:
time: 8.45677018165588
memory: 2,977,234,944
Edit 2. By the way, while we are at this example, the PDL::NiceSlice says
- arguments should not be quoted
- lone '*' inserts a dummy dimension of order 1
but I see that lone unquoted star(s) cause Perl to produce a warning:
Use of uninitialized value in concatenation (.) or string
so I had to insert a pair of very unsightly "1"'s, which of course devalues code nicety by half at least:). Is this a bug in PDL::NiceSlice, documentation flaw, or just me reading it incorrectly?
| [reply]
[d/l]
[select] |
|
|
|
|
