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).

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.

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.

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]
]
[download]

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
[download]

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?

your approach really physicalizes arrays from 2D to 4D

No, it does not. It uses clump and slice, both of which are affine transformations (as shown by the P2Child=>1 in https://github.com/PDLPorters/pdl/blob/master/Basic/Slices/slices.pd for the relevant operations). Physicalising only happens on output or use in calculations, and when loop fusion comes in, not even then. I'd love to hear of an approach that Kronecker-products m1,n1 and m2,n2 matrices without producing a m1*m2,n1*n2 in-memory object. While you're pondering that, please show some numbers demonstrating these "grave consequences" :-)

Is this a bug in PDL::NiceSlice, documentation flaw, or just me reading it incorrectly?

You're going to need to show some code that you're saying is behaving unexpectedly. This worked as described (with output):

pdl> p sequence(2)->(*1)
[
 [0]
 [1]
]
pdl> p sequence(2)->(*)
[
 [0]
 [1]
]
pdl> p sequence(2)->(*,*)
[
 [
  [0]
 ]
 [
  [1]
 ]
]
[download]

Naturally, you need to make sure a use PDL::NiceSlice is in scope. You will have noticed that generally I avoid using it in core PDL code, because it both slows down module-loading, and has limitations; see PDL::LinearAlgebra, where I had to remove a lot of it in favour of building arguments to slice that took account of adding optional parts for PDL::Complex.