Comment by imurray | Hacker News Reader

Comment by imurray | Hacker News Reader

Parent Full thread

I don't know APL, but from what I found quickly, APL's "conformability" rules[1] only seem to be a subset of the Yorick/NumPy broadcasting? Here's a Python example of broadcasting that APL doesn't seem to support(?): subtracting arrays of shape (3,1,7,5) and (1,4,1,5) to get an array of shape (3,4,7,5).

    from numpy.random import randn
    randn(3,1,7,5) - randn(1,4,1,5)

Inspired by other comments in this thread: The second line can also be run in current versions of Octave/Matlab as-is, BUT they didn't support it before NumPy, and neither did/does FORTRAN[2]. Mathworks resisted introducing this general form of broadcasting into Matlab for a long time. Octave then did it anyway (version 3.6.0 in 2012), and Matlab eventually followed (version R2016b). Before late 2016, Matlab code required `bsxfun(@minus, A, B)` instead of `A-B` to get broadcasting, and before bsxfun (introduced in R2007a) it was even more awkward (code often used `repmat`, or indexing tricks to explicitly expand arrays in memory to match shape).

[1] https://aplwiki.com/wiki/Conformability

[2] suggestion from 2022 to introduce broadcasting to Fortran: https://github.com/j3-fortran/fortran_proposals/issues/252

There was a little discussion about this just Monday in the array programming chat[0]. Not many people commented but it doesn't look like anyone was in favor. I think the concern is that while broadcasting 1s is an easy mechanism to understand, it changes the problem from "how do I match up array axes" to "how to I get 1s in the right places in the shape", which is no easier to do, and further from the programmer's intent. The modern APL way is to use the Rank operator, which can specify any broadcasting pattern. The example above would be f⍤1⍤2 1⍤2 on shapes 3 7 5 and 4 5, although in practice these sorts of alternating broadcasts don't show up much and even two copies of Rank is rare. One quote from the chat:

> a friend had written some code in—either numpy, or something else borrowing its model—and wanted to translate it to j. I could not figure out what she was trying to accomplish with the many seemingly-extraneous reshapes and transposes, until I found out what they did in numpy. The whole thing reduced to a single, trivial application of rank in j

[0] https://aplwiki.com/wiki/APL_Farm, main or "Array languages" channel

Thanks for the pointer. I can believe that a language that looks so different will find that different patterns and primitives are natural for it.

My experience from writing a lot of array-based code in NumPy/Matlab is that broadcasting absolutely has made it easier to write my code in those ecosystems. Axes of length 1 have often been in the right places already, or have been easy to insert. It's of course possible to create a big mess in any language; it seems likely that the NumPy code you saw could have been neater too.

In machine learning there can be many array dimensions floating around: batch-dims, sequence and/or channel-dims, weight matrices, and so on. It can be necessary to expand two or more dimensions, and/or line up dimensions quite carefully. Einops[1] has emerged from that community as a tool to succinctly express many operations that involve lots of array dimensions. You're likely to bump into more and more people who've used it, and again it seems there's some overlap with what Rank does. (And again, you'll see uses of Einops in the wild that are unnecessarily convoluted.)

[1] https://einops.rocks/ -- It works with all of the existing major array-based frameworks for Python (NumPy/PyTorch/Jax/etc), and the emerging array API standard for Python.

Of course the method with length-1 axes will look good if it's the only way to do broadcasting. I've hardly touched it (and never saw the NumPy code mentioned in the quote), but some people in the APL community have so I'm reporting on their thoughts. Although length-1 axes already being in the data sounds worrying: you have an array that doesn't depend on an axis, but there's an axis anyway to show where it would go if it did?

Named axis systems like einops are much better regarded. I've seen some implementations of similar functionality in APL, and various discussions of how it could be built into an array language (I have my own design where that's a core idea, but not using APL syntax).

> Although length-1 axes already being in the data sounds worrying: you have an array that doesn't depend on an axis, but there's an axis anyway to show where it would go if it did?

Example:

    A / A.sum(axis=3, keepdims=True)

Make all of the vectors along axis=3 sum up to one. There are other ways of doing it, but this way seems fairly clear to me. The shape of the denominator is the same as the numerator, except for a 1 in position 3. Unfortunately we have to specify `keepdims`, because the default of `False` removes the dimension being summed over, which doesn't work in general. `keepdims=True` is the behavior in Matlab/Octave, so the example becomes

    A ./ sum(A, 4)

with 4=3+1 because Matlab is 1-based like Fortran.

Reply on news.ycombinator.com