NumPy のブロードキャスティングを PHP で模倣する

NumPy の配列には、ブロードキャスティングという仕組みがあります。これは、大きさの異なる多次元配列同士の算術演算に関する仕組みで、要素数が 1 の次元があった場合にそれを拡張し、二つの配列の要素数を揃えて演算結果を得るというものです。NumPy のマニュアルでは、以下のページに記載があります。また、そのページの末尾からリンクされている EricsBroadcastingDoc には具体例をともなう説明があります。
Broadcasting — NumPy v1.9 Manual

このブロードキャスティングを PHP で実装してみました。コードは以下のとおりです。関数名の bsxfun は Binary Singleton eXpansion FUNction の略です。ブロードキャスティング相当の機能は、MATLAB では bsxfun 関数で提供されており、今回はその関数名を使いました。

<?php
function bsxfun($f, $a, $b) {
    $da = depth($a);
    $db = depth($b);
    $a = newAxis($a, $db - $da);
    $b = newAxis($b, $da - $db);
    return bsxfunRec($f, $a, $b);
}

function depth($a) {
    return is_array($a) ? depth($a[0]) + 1 : 0;
}

function newAxis($a, $n) {
    for ($i = 0; $i < $n; ++$i) {
        $a = array($a);
    }
    return $a;
}

function bsxfunRec($f, $a, $b) {
    if (!is_array($a)) {
        return $f($a, $b);
    }
    $a = expandIfSingleton($a, count($b));
    $b = expandIfSingleton($b, count($a));
    return array_map(
        function ($a, $b) use ($f) { return bsxfunRec($f, $a, $b); },
        $a,
        $b);
}

function expandIfSingleton($a, $n) {
    if (count($a) == 1) {
        return array_fill(0, $n, $a[0]);
    }
    return $a;
}

簡単なテストプログラムでの実行例は次のとおりです。$a は 1 行 3 列の行列、$b は 3 行 1 列の行列を表しています。これらの配列に対して bsxfun 関数で要素ごとの積を求めると、3 行 3 列の行列が得られます。

$ cat bsxfun_test.php
<?php
require 'bsxfun.php';

$a = [1, 2, 3];
$b = [[1], [2], [3]];
$c = bsxfun(function ($a, $b) { return $a * $b; }, $a, $b);
echo json_encode($c) . "\n";

$ php bsxfun_test.php
[[1,2,3],[2,4,6],[3,6,9]]

実装した bsxfun 関数は、まず、引数の二つの配列の次元数をそれぞれ計算し (depth 関数)、次元数が異なる場合には、次元数の小さな側を配列で包んでいくことで、両者の次元数を揃えます (newAxis 関数)。この処理は、NumPy のマニュアルの以下の記述に対応します。この例では、一次元配列の Scale を 1x1x3 の三次元配列に変換してから後続の処理を行うようにします*1。

Arrays do not need to have the same number of dimensions. For example, if you have a 256x256x3 array of RGB values, and you want to scale each color in the image by a different value, you can multiply the image by a one-dimensional array with 3 values. Lining up the sizes of the trailing axes of these arrays according to the broadcast rules, shows that they are compatible:

Image  (3d array): 256 x 256 x 3
Scale  (1d array):             3
Result (3d array): 256 x 256 x 3

bsxfunRec 関数では、要素数が 1 の次元を拡張しながら (expandIfSingleton 関数)*2、再帰呼び出しによって各要素に $f を適用します。この部分は、マニュアルの以下の記述に対応します。

When either of the dimensions compared is one, the other is used. In other words, dimensions with size 1 are stretched or “copied” to match the other.
In the following example, both the A and B arrays have axes with length one that are expanded to a larger size during the broadcast operation:

A      (4d array):  8 x 1 x 6 x 1
B      (3d array):      7 x 1 x 5
Result (4d array):  8 x 7 x 6 x 5

なお、今回の実装では、入力の配列がブロードキャスティングの条件を満たさない場合を考慮していません。これは二つの場合があります。一つは、両者の要素数がどちらも 2 以上で、かつ、異なっている場合です。NumPy のブロードキャスティングでは、そのような場合にはエラーになると定められています。もう一つは、ジャグ配列が与えられた場合です。NumPy では、矩形配列の場合に限りブロードキャスティングが適用されるようです*3。

*1:今回の私の実装では、ここで実際に配列を作成して次元数を揃えていますが、本来そうする必要はありません。反復処理で実装することもできます。

*2:こちらも array_fill で要素数を揃えた配列を実際に作成していますが、本来その必要はありません。

*3:私は Python や NumPy には全く詳しくないのですが、NumPy の array では、ジャグ配列は dtype=object となり、矩形配列とは異なる内部構造を持つようです。