Peter
Peter Campbell Smith

Primorials and
Kronecker products

Weekly challenge 170 — 20 June 2022

Week 170 - 20 Jun 2022

Task 2

Task — Kronecker product

You are given 2 matrices, @A and @B. Write a script to calculate the Kronecker Product of the matrices as described on this Wikipedia page.

Examples


Example 1:
A = [ 1 2 ]
    [ 3 4 ]

B = [ 5 6 ]
    [ 7 8 ]

A ⊗ B = [ 1 × [ 5 6 ]   2 × [ 5 6 ] ]
        [     [ 7 8 ]       [ 7 8 ] ]
        [ 3 × [ 5 6 ]   4 × [ 5 6 ] ]
        [     [ 7 8 ]       [ 7 8 ] ]

      = [ 1×5 1×6 2×5 2×6 ]
        [ 1×7 1×8 2×7 2×8 ]
        [ 3×5 3×6 4×5 4×6 ]
        [ 3×7 3×8 4×7 4×8 ]

      = [  5  6 10 12 ]
        [  7  8 14 16 ]
        [ 15 18 20 24 ]
        [ 21 24 28 32 ]

Analysis

Leopold Kronecker
Leopold Kronecker

You come across Leopold Kronecker's musings all over the place, and I see that at least a dozen mathematical 'things' are named after him. Fortunately, this one isn't too hard to follow, thanks to someone giving the key to solving it in Wikipedia, albeit in very small print.

As with some other recent challenges, it takes more lines code to display the data neatly than to perform the required calculation.

I used two helper functions. First is quorem which takes two arguments, a dividend and a divisor, and returns the quotient and the remainder as per primary school division. The second is show which simply displays a matrix with the columns neatly lined up. I use that to display the two input matrices ($A and $B) and the resulting product matrix ($K).

Wikipedia kindly gives us the formula for each element K(i, j) as a function of two elements of A and B. The calculation of which elements of A and B are needed involves some modular arithmetic, which is where quorem comes in.

The show function makes two passes over the matrix to be displayed. In the first pass it determines which element is the widest, for which Perl's length function is useful as it intrinsically copes with minus signs. The cell width is then made one space wider than the widest number for a pleasing result.

I chose to store the 3 matrices as 2-dimensional arrays - actually arrays of references, so that element X(i, j) is $X->[$i]->[$j]. This nicely parallels the mathematical representation.

Try it 

Try running the script with any input:



example: [1, 2], [3, 4]



example: [5, 6], [7, 8]

Script


#!/usr/bin/perl

# Peter Campbell Smith - 2022-06-20
# PWC 170 task 2

use v5.28;
use strict;
use warnings;
use utf8;
binmode(STDOUT, ':utf8');

my ($test, $A, $B, $K, $m, $n, $p, $q, $i, $j, $i_quo_p, $i_rem_p, $j_quo_q, $j_rem_q);

for $test (1, 2) {
    
    if ($test == 1) {   # Mohammad's example

        $A = [[ 1, 2 ],
              [ 3, 4 ]];

        $B = [[ 5, 6 ],
              [ 7, 8 ]];
              
    } else {   # Wikipedia's 2nd example

        $A = [[ 1, -4, 7],
              [-2,  3, 3]];
              
        $B = [[8, -9, -6,  5],
              [1, -3, -4,  7],
              [2,  8, -8, -3],
              [1,  2, -5, -1]];
    }
    
    # display input
    show('A', $A);
    show('B', $B);
          
    $m = scalar(@$A);           # no of rows in A (0 to $m - 1)
    $n = scalar(@{$A->[0]});    # no of columns in A
    $p = scalar(@$B);           # no of rows in B 
    $q = scalar(@{$B->[0]});    # no of columns in B

    # loop over the rows and cols of product matrix $K ($m * $p rows, $n * $q cols)
    for $i (0 .. $m * $p - 1) {
        for $j (0 .. $n * $q - 1) {
            
            # get the quotient and remainder on dividing i by p and j by q
            ($i_quo_p, $i_rem_p) = quorem($i, $p);
            ($j_quo_q, $j_rem_q) = quorem($j, $q);
            
            # Wikipedia formula for element i, j of K
            $K->[$i]->[$j] = $A->[$i_quo_p]->[$j_quo_q] * $B->[$i_rem_p]->[$j_rem_q];
        }
    }
    
    # display the result
    show('A × B', $K);
    print qq[-----\n\n];
}
        
sub quorem {

    my ($dividend, $divisor, $quotient, $remainder);

    # returns the integer quotient and remainder from dividend / divisor (both integers > 0)
    $dividend = shift;
    $divisor = shift;
    
    $quotient = int($dividend / $divisor);
    $remainder = $dividend - $quotient * $divisor;
    
    return ($quotient, $remainder);
}

sub show {
    
    my ($caption, $M, $margin, $i, $j, $this_width, $cell_width);
    
    # pretty print matrix
    ($caption, $M) = @_;
    
    $caption .= ' = ';
    $margin = length($caption);
    
    # first pass to get width of cells
    $cell_width = 0;
    for $i (0 .. scalar(@$M) - 1) {
        for $j (0 .. scalar(@{$M->[0]}) - 1) {
            $this_width = length($M->[$i]->[$j]);
            $cell_width = $this_width if $this_width > $cell_width;
        }
    }
    $cell_width += 1;
    
    # second pass to print out - loop over rows
    for $i (0 .. scalar(@$M) - 1) {
        print qq{$caption\[};       # first line
        $caption = ' ' x $margin;   # subsequent lines
        
        # and over columns, padding each value to $cell_width
        for $j (0 .. scalar(@{$M->[0]}) - 1) {
            printf("%${cell_width}s", $M->[$i]->[$j]);
        }
        print qq[ \]\n];
    }
    print qq[\n];
}

Output


A = [ 1 2 ]
    [ 3 4 ]

B = [ 5 6 ]
    [ 7 8 ]

A ⊗ B = [  5  6 10 12 ]
        [  7  8 14 16 ]
        [ 15 18 20 24 ]
        [ 21 24 28 32 ]

-----

A = [  1 -4  7 ]
    [ -2  3  3 ]

B = [  8 -9 -6  5 ]
    [  1 -3 -4  7 ]
    [  2  8 -8 -3 ]
    [  1  2 -5 -1 ]

A ⊗ B = [   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 ]

-----