The problem

The problem, in short, is that we have an array of matrices and we want to apply a function "matrixwise" across the array. For example, let us consider the two-matrix case:

    my $A = [ [ 1, 2 ],
              [ 3, 4 ] ];   # matrix A

    my $B = [ [ 5, 6 ],
              [ 7, 8 ] ];   # matrix B

    my $mats = [ $A, $B ];  # array of matrices
[download]

Our goal is to compute the following matrix from $mats:

    # [ [ f(1,5), f(2,6) ],
    #   [ f(3,7), f(4,8) ] ]
[download]

for some given function f. (In the case of the knirirr's problem, as I understand it, f will compute the mean and standard deviation of its arguments.)

We can think about the array of matrices $mats as itself a 3-dimensional matrix. Its dimensions, most significant first, are (m,r,c), where m gives the position of the sub-matrix we want, r gives the row, and c gives the column:

    #                            c=0  c=1
    #                           ________
    #   Matrix A (m=0)         /   /   /
    #                    r=0  / 1 / 2 /.
    #                        /___/___/ .
    #                       /   /   /  .
    #                 r=1  / 3 / 4 /   .
    #                     /___/___/    .
    #                                  .
    #                           ________
    #   Matrix B (m=1)         /   /   /
    #                    r=0  / 5 / 6 /
    #                        /___/___/
    #                       /   /   /
    #                 r=1  / 7 / 8 /
    #                     /___/___/
    #
    #                   c=0  c=1
[download]

For example, the value at coordinate (0,0,0) is 1; and the value at (1,0,0) is 5.

A common solution

With this representation in mind, we can arrive at the following, typical approach to a solution:

    my $result;
    for my $r (0 .. @{$mats->[0]} - 1) {
        for my $c (0 .. @{$mats->[0][0]} - 1) {
            $result->[$r][$c] = f( map $_->[$r][$c], @$mats );
        }
    }
[download]

To see that it works, let us create a few functions to help us examine the results:

    use Data::Dumper;

    sub f {
        "f(" . join(", ", map dumpf($_), @_) . ")";
    }

    sub dumpf {
        local $Data::Dumper::Indent = 0;
        local $Data::Dumper::Terse = 1;
        local $_ = Dumper(@_);
        tr/\'//d;
        $_;
    }
[download]

Now we can see that the nested-for-loop solution does indeed compute the desired result:

    print dumpf($result), $/;
    # [ [f(1,5), f(2,6)],
    #   [f(3,7), f(4,8)] ]
[download]

One problem with the nested-for solution, however, is that it is hardcoded for the specifics of this particular, 3-dimensional case. Another problem is that the code has a high density of array operations, providing many opportunities to mix up index variables or make similar coding mistakes.

Nested-data helpers

In day-to-day coding, we frequently process arrays, and arrays of arrays, and AoAoAs, and similar structures. It might therefore be worthwhile to invest some time in helper functions that let us manipulate nested data more directly.

One such helper is the ever-handy mapat. It is similar to map except that it allows us to map at a particular depth of nesting in arbitrarily-nested arrays:

    sub mapat {
        my $depth = shift;
        my $f     = shift;
        $depth
            ? map [ map mapat($depth-1, $f, $_), @$_ ], @_
            : $f->(@_);
    }
[download]

To see how mapat works, let us target the various depths of nesting in our original 3-dimensional input matrix $mats:

    for (0..3) {
        print "mapat $_ f: ",
              dumpf(mapat($_, \&f, $mats)), $/;
    }
    # mapat 0 f: f([[[1,2],[3,4]],[[5,6],[7,8]]])
    # mapat 1 f: [f([[1,2],[3,4]]),f([[5,6],[7,8]])]
    # mapat 2 f: [[f([1,2]),f([3,4])],[f([5,6]),f([7,8])]]
    # mapat 3 f: [[[f(1),f(2)],[f(3),f(4)]],[[f(5),f(6)],[f(7),f(8)]]]
[download]

Note that mapat preserves the structure of its input down to the targeted depth. This is an important property that lets us process nested values "in structure," without having to tear down and build up arrays manually, like we did with the nested-for-loop solution.

While it is probably clear that mapat is a general and useful function for working with nested arrays, the function doesn't actually solve our original problem. Examining the output, we can see that mapat(2, ...) provides the correct output structure, but not the correct values.

The problem is that the elements passed to f are not grouped matrixwise. We are passing in rows extracted from our 3-dimensional array, and we ought to be passing in matrixwise, row-column stacks. That is, for each sub-matrix M and row R, we are passing to f all of the values whose coordinates are given by (m=M,r=R,c=*), where * matches any value. We ought to be passing in (m=*,r=R,c=C) for each valid row-column pair (R,C).

To provide the correct values to mapat, we can restructure $mats so that m becomes its least-significant dimension. In other words, we can convert $mats's dimensions from (m,r,c) to (r,c,m). Then mapat(2, ...) can be used to apply f to exactly the desired elements: (r=R,c=C,m=*).

We can use matrix transposition to perform the restructuring. Transposition swaps the order of the two most-significant dimensions of a matrix. It is commonly used in spreadsheets to convert rows into columns and vice versa. Here is how it works for the 2-dimensional case:

    # Dimensions
    #
    #   (r,c)        <=== transpose ===>    (c,r)
    #
    # Data
    #
    # [ [ 1, 2 ],                         [ [ 1, 3, 5 ], 
    #   [ 3, 4 ],    <=== transpose ===>    [ 2, 4, 6 ] ]
    #   [ 5, 6 ] ] 
    #
[download]

I usually implement transpose in terms of zip, which is a handy helper for functional programming. There are many other implementations lying around if you don't like this one:

    sub zip {
        my $len = min( map scalar @$_, @_ );
        map [ do { my $i=$_; map $_->[$i], @_ } ], 0..$len-1;
    }

    sub transpose {
        map [ zip(@$_) ], @_;
    }
[download]

With transpose in hand, we can restructure our original input matrix in two steps. First, we transpose at depth 0 to reorder the first two dimensions:

    #  0 1 2      0 1 2
    # (m,r,c) -> (r,m,c)  # via transpose(...)
    #
    # ([[[1,2],[3,4]],[[5,6],[7,8]]]
    #         -> [[[1,2],[5,6]],[[3,4],[7,8]]]
[download]

Then, we use mapat to transpose the result of the previous step at depth 1, in effect transposing the result's sub-matrices. This swaps the second and third dimensions:

    #  0 1 2      0 1 2
    # (r,m,c) -> (r,c,m)  # via mapat(1, \&transpose, ...)
    #
    # [[[1,2],[5,6]],[[3,4],[7,8]]]
    #         -> [[[1,5],[2,6]],[[3,7],[4,8]]]
[download]

Note that the output values are grouped matrixwise. This is exactly the grouping that we want.

The final step is to use mapat(2, ...) to apply f to the second "level" of our restructured matrix:

    mapat(2, \&f, [[[1,5],[2,6]],[[3,7],[4,8]]]);
    # [[f([1,5]),f([2,6])],[f([3,7]),f([4,8])]]
[download]

And that is our desired output.

Putting all of the above together results in our final function, matrixwise_map:

    sub matrixwise_map(&@) {
        my $f = shift;
        # force array context since we're always dealing w/ matricies
        ( mapat( 2, $f, mapat(1, \&transpose, transpose(@_)) ) )[0];
    }
[download]

It does exactly what we want, in one fell swoop:

    $result = matrixwise_map(\&f, $mats);

    print dumpf($result), $/;
    # [[f([1,5]),f([2,6])],[f([3,7]),f([4,8])]]
[download]

Back to our original problem

To return to the knirirr's original problem, here is how we might compute the matrixwise mean of the input matrices:

    use List::Util qw( sum );
    sub mean { sum(@_) / @_ }

    print dumpf( matrixwise_map {mean(@$_)} $mats ), $/;
    # [[3,4],[5,6]]
[download]

The matrixwise standard deviations can be computed similarly: Create a stdev function and matrixwise_map it over $mats. That's it!

What have we gained?

In solving a particular nested-array problem, we created two small functions mapat and transpose that we can use to solve many kinds of related problems. The mapat function lets us transform nested array data without disturbing its structure, and transpose lets us restructure dimensional data by rearranging dimensions. As an example of what these functions can do, we created matrixwise_map. It lets us apply functions matrixwise over an array of matrices.

Together or alone, all three of these functions are handy. If you find yourself writing a lot of code to process array data, consider whether you could benefit from functions like these or from problem-solving approaches similar to those we explored in this meditation.

Thanks!

Thanks for taking the time to read this meditation. If you can see any way to improve my writing or explanations, please let me know.

Cheers,
Tom