A convolutional layer arrives in the literature as a triple sum. Zhao, Wang, Wang and Liu, A Faster Algorithm for Reducing the Computational Complexity of Convolutional Neural Networks (Algorithms, 2018), §2.1, defines one output cell:

               Q      M_w    N_w
  y[r,x,y]  =  ∑      ∑      ∑     W[r,q,u,v] · X[q, x+u, y+v]
              q=1    u=1    v=1
Symbol Range Meaning
r output filters selects the output feature map
q input channels selects the input feature map
u, v kernel rows, columns position inside the kernel window
x, y output positions position in the output feature map
Q, M_w, N_w — input channels, kernel height, kernel width

TL;DR — five additions to sw-MLPL, and what each one lets you write:

Keyword What it enables
windows(x, sizes) Every sliding neighborhood as an array — convolution, moving averages, stencils, any local operator, with no index arithmetic
trailing-axis broadcasting One small kernel meets a whole stack of patches without being copied out to match
reduce(:op, a, [2,3,4]) Several summations collapsed in one pass, in any order
reduce(:op, a, "channel") Axes chosen by meaning instead of position
grad through windows The convolution stops being something you compute and becomes something you train
Resource Link
sw-MLPL sw-ml-study/sw-mlpl
Playground mlpl.softwarewrighter.com — 0.22.0; everything here runs in the browser
The demos sw-ml-study/demo-ml-utils
The paper Zhao, Wang, Wang & Liu, Algorithms 11(10):159, 2018
Comments Discord

Two notes before the code, because both are places a transcription quietly stops matching its source. This is cross-correlation: X[q, x+u, y+v], not x-u, y-v, so the kernel is never flipped. Every mainstream framework does this and calls it convolution, and the native conv2d builtin agrees, but it is worth naming. And the paper counts from 1 while the array language counts from 0 — q = 1..Q is 0..Q-1 in code, and reduction axis numbers are 0-based too.

windows — neighborhoods become values

A 2D image with one neighborhood lifted out into a grid of patches shaped out_y by out_x by C by kh by kw

The hard term is X[q, x+u, y+v]. That index arithmetic says: for every output position, look at the small block of input around it. Written as loops, you manage the position counters and the offsets yourself, and the shape of the thing you are working on never appears in the source at all.

patches = windows(x, [3, 3])

A [channel, image_y, image_x] input becomes [out_y, out_x, channel, 3, 3]: a grid with one cell per output position, each cell holding an entire channel × 3 × 3 neighborhood. The offsets are gone from the source and have become structure in the array.

What it enables beyond convolution. Any operator defined on a neighborhood is now a two-step expression — take the windows, reduce them:

moving_avg = reduce(:add, windows(x, [5]), 1) / 5

Blur kernels, edge detectors, pooling, Game of Life neighbor counts, finite-difference stencils, moving-window statistics — all the same shape of statement. It windows the trailing axes and leaves earlier ones alone, so a channel or batch axis rides along untouched, and it takes an optional stride when you want the windows to skip rather than slide by one.

Broadcasting — one kernel, every patch

A single small kernel aligned against every patch in the stack and reused rather than copied

Now the multiply. A kernel is [channel, kernel_y, kernel_x]; the patches are [out_y, out_x, channel, kernel_y, kernel_x]. The same small kernel applies at every position, which is what “shared weights” means in a CNN.

Trailing-axis broadcasting lets those meet directly:

weighted = kernel * patches

The kernel’s three axes line up against the patches’ last three, and it is reused across the leading position axes rather than copied. The benefit is the copy that does not happen. Without it you must first replicate the kernel out to the full shape of the patch stack — for a modest [8,32,32] input that is over a million redundant values built and multiplied so the shapes would match. The equation says one kernel is applied everywhere; broadcasting is what lets the code say that too.

Multi-axis reduction — three sums, one statement

A rank-5 patch stack collapsing in one pass to a rank-2 out_y by out_x result

What remains is the summing, over channel, kernel row and kernel column at once.

One axis at a time, that is three nested calls, read inside-out and in reverse:

y = reduce(:add, reduce(:add, reduce(:add, weighted, 4), 3), 2)

The nesting also implies an order the mathematics does not have: ∑_q ∑_u ∑_v is three interchangeable sums, not three stacked passes. Reduction over a vector of axes says it once:

y = reduce(:add, weighted, [2, 3, 4])

One pass, rank drops by three, order irrelevant. Beyond convolution this is the ordinary case of summing a tensor down to what you actually want — totals across batch and spatial axes while keeping channels, and so on — without stacking calls to get there.

Named axes — meaning instead of position

Integers still make the reader do bookkeeping. Which axis was 2?

y = reduce(:add, weighted, "channel")

Several axes at once are named as one comma-separated string, "channel,kernel_y,kernel_x". The paper’s q, u, v are placeholders you decode from surrounding prose. channel, kernel_y, kernel_x are not. And unlike a comment, a label is checked: name an axis that does not exist and you get an error rather than a wrong answer. Labels also survive the operations that reshape and filter an array, so a name attached early still means something several steps later.

The layer

@formula "y[r,x,y] = ∑(q=1..Q) ∑(u=1..M_w) ∑(v=1..N_w) W[r,q,u,v] · X[q, x+u, y+v]"
def u:conv_layer(x, w, kh, kw) {
    "One convolutional layer: window, multiply, reduce.";
    p = label(windows(x, [kh, kw]), ["out_y", "out_x", "channel", "kernel_y", "kernel_x"]);
    reduce(:add, w * p, "channel,kernel_y,kernel_x")
}

Window, multiply, reduce. The annotation carries the source equation as data — readable at runtime, so the mathematics travels with the function instead of in a comment that drifts away from it.

The two lines name axes two different ways, which is worth flagging rather than glossing: label takes a bracketed list of names, while reduce takes them as one comma-separated string, and each rejects the other’s form. The labeling step is needed because only some names survive the window — a labeled [channel, image_y, image_x] input comes out as patches labeled image_y, image_x, channel with the two kernel axes unnamed, so the kernel axes have to be named before they can be reduced by name.

It is also the fast version

The reasonable worry is that this is a teaching toy you abandon for real work. Laying the windows out as a matrix and letting matmul contract them:

cols = reshape(windows(x, [kh, kw]), [oy * ox, c * kh * kw]);
y    = matmul(cols, reshape(kernel, [c * kh * kw]))

On an [8,32,32] input against [16,8,3,3] filters this runs in 1.198 ms against the hand-written native conv2d builtin’s 1.230 ms, and agrees with it exactly rather than within a tolerance. The readable spelling is not the slow one.

And then you can train it

Gradients from overlapping windows flowing back and accumulating into one shared input cell

Everything above computes a convolution. None of it learns one — and until now, nothing else in the language did either. The native conv2d is forward-only, so a kernel could be applied but never fitted; the array spelling had the same problem one level down, because windows had no gradient. A CNN whose filters cannot be learned is a filter bank.

Making windows differentiable is what closes that, and its backward pass is a nice piece of arithmetic in its own right. The forward direction gathers: each output position collects the neighborhood around it. Reversing that means each input cell must collect the gradient of every window it appeared in — and because the windows overlap, most cells appear in several. So the backward of a sliding-window gather is a scatter-add: contributions accumulate rather than overwrite. The overlap that makes a convolution a convolution is exactly what makes its gradient an accumulation.

Making that work meant a gradient for every step of the forward path, not just the window: the reduction had to reverse, and so did the broadcast — a kernel reused across a thousand positions collects gradient from all of them, so the backward of a broadcast is a sum over the axes it was stretched along.

The result is that the expression which reads like the paper’s equation is the expression you train. Both spellings are differentiable, so the matmul form is trainable too, and a reader who followed the forward arithmetic has already met every operation the backward pass reverses.

Why this was worth doing

conv2d already existed. None of these were needed to compute a convolution.

Four of them change how the layer reads. The offsets became an array. The three sums became one reduction. The axes acquired names. The kernel stopped being copied. Someone who has read the paper now recognizes the code, and someone who has read the code can reconstruct the paper — the notation and the executable form are the same object rather than two things you hope stay in sync.

The fifth changes what the layer is. A convolution you can only evaluate is a fixed filter, however elegantly it is written; a convolution you can differentiate is a layer that learns. The equation is no longer just something the code resembles — it is something the code does, and then improves at.