Convolution: local shared operators, shapes and gradients

Convolution makes two assumptions explicit: nearby measurements interact, and a detector should be reusable at different locations. Neither assumption says a classifier is invariant to every transformation. This lesson builds the operator, its gradient, and the design choices that remain useful in modern vision systems.

Why locality and sharing matter

LeCun and colleagues' 1998 Gradient-Based Learning Applied to Document Recognition is a primary account of convolutional document recognition, including LeNet-5 and deployed document systems. It is not evidence for the original lesson's imagined postal-service conversation or claim that nobody used neural vision for fourteen years. The durable historical lesson is architectural: learn local detectors jointly with the classifier, and reuse their weights across space.

For a 224×224 RGB image, flattening gives 150,528 inputs. A dense map to 1,024 units has 154,140,672 weights and 1,024 biases. FP32 weights alone occupy 616,562,688 bytes; including the1,024 biases gives616,566,784 bytes. Parameters, gradients and two Adam moments require roughly four parameter-sized arrays before activations and allocator overhead. This is a memory calculation, not a guarantee of out-of-memory on every laptop or a theorem that a billion labels are necessary.

A local unshared layer would already remove distant connections. Convolution adds sharing: one output-channel detector uses the same coefficients at every spatial position. A 7×7, 3-input/64-output convolution has 9,408 weights and 64 biases. Its output is a spatial map, not the dense layer's 1,024-vector, so the comparison illustrates different structural assumptions, not interchangeable layers at equal output shape. Shared parameters reduce statistical degrees of freedom, but generalization still depends on the data and training.

Prerequisites and outcomes. Use array slicing, scalar chain rules and shared-gradient accumulation. The CPU programs require NumPy, SciPy (correlate2d) and PyTorch; run the Python blocks in order so later blocks reuse their imports and definitions. By the end, compute the original 5×5/3×3 example, predict NCHW shapes including dilation/groups, implement forward and backward, and choose downsampling by its semantics rather than folklore.

Cross-correlation by hand

PyTorch2.9 Conv2d computes cross-correlation: the kernel is not spatially reversed. Mathematical convolution reverses it. A learned unconstrained filter family can absorb that reversal, but fixed filters, imported weights and backward implementations cannot ignore the convention.

For the original image and edge detector,

X = [[1,2,3,0,1],       K = [[1,0,-1],
     [0,1,2,3,1],            [1,0,-1],
     [3,1,0,2,2],            [1,0,-1]]
     [2,0,1,3,1],
     [1,3,2,1,0]]

At (0,0), row contributions are −2, −2, +3, giving −1. At (0,1), they are +2, −2, −1, again −1. The full output below corrects the old printed matrix; an output that merely has the right shape is not a numerical oracle.

import numpy as np
from scipy.signal import correlate2d
X = np.array([[1,2,3,0,1], [0,1,2,3,1], [3,1,0,2,2],
              [2,0,1,3,1], [1,3,2,1,0]], dtype=float)
K = np.array([[1,0,-1], [1,0,-1], [1,0,-1]], dtype=float)
Y = correlate2d(X, K, mode='valid')
expected = np.array([[-1,-1,1], [2,-6,-1], [3,-2,0]], dtype=float)
assert np.array_equal(Y, expected)
horizontal = np.array([[1,1,1], [0,0,0], [-1,-1,-1.]])
assert correlate2d(X, horizontal, mode='valid')[1,1] == 2
assert not np.array_equal(Y, correlate2d(X, K[::-1,::-1], mode='valid'))
print('original 5x5 output', Y.tolist())

The horizontal-kernel practice answer is 2: at (1,1), the upper row sums to 6 and lower row to 4. The asymmetric vertical filter supplies a negative control for accidentally reversing the kernel.

Shape arithmetic, padding and channel layout

For either spatial axis,

The effective span is $d(k-1)+1$. A three-point filter at dilation2 spans5 positions, not6. Count valid window starting locations, not just a vaguely halved image size. Padding is added on both sides; a remainder under stride means the final incomplete window is not evaluated.

Thus the original batches become (1,64,112,112), (1,128,28,28), and (8,16,28,28). With stride1 and odd effective kernel size, symmetric padding $p=d(k-1)/2$ preserves size. valid means no padding. PyTorch2.9's padding='same' supports stride1, whereas other frameworks can use SAME to mean ceil(input/stride); do not transfer string semantics across APIs. Even kernels are valid. Their alignment may require asymmetric padding; an explicit integer padding=1 does not secretly choose a different convention because the kernel is even. Track receptive-field centers when aligning segmentation masks.

PyTorch's logical shape is NCHW: batch, channels, height, width. TensorFlow commonly uses NHWC. Logical dimension order is distinct from memory layout: PyTorch can store a logical NCHW tensor in channels-last memory format. Performance depends on kernels and hardware, not a general claim that channels-first is faster. PIL images, float arrays, uint8 arrays and image loaders also differ in range; division by255 is not correct for every floating input. Inspect shape, dtype and range at the input boundary.

Multiple channels, groups and parameter accounting

For groups $g$, weights have shape $(C_{out},C_{in}/g,k_H,k_W)$ and both channel counts must be divisible by $g$. Each output sums over its connected input channels and spatial window, then adds one bias per output channel. With g1,32 inputs,64 outputs,k3, there are18,432 weights plus64 biases, not576 weights. A 1×1 convolution mixes channels at each location: 32→64 uses2,048 weights before biases; it is not a spatial no-op in feature space.

For depthwise convolution with multiplier1, $g=C_{in}=C_{out}=C$. A depthwise k×k layer followed by a pointwise C→C map uses $k^2C+C^2$ weights, instead of $k^2C^2$. This factorization constrains the operator and is not generally equivalent to an arbitrary dense convolution. With C32,k3, counts are1,312 versus9,216. More MACs per parameter can result from weight reuse: fewer parameters does not automatically mean lower latency.

A classifier typically chains convolution, nonlinearity and downsampling, then averages spatial features to (B,C) and maps them to class logits. Backpropagation through a spatial mean distributes each channel derivative by1/(HW), then activation derivatives and shared-kernel accumulation apply. A residual branch adds a shape-compatible shortcut: for y=x+F(x), the Jacobian is I+J_F and can still cancel. Shape-changing shortcuts need a projection. DL026 develops the full residual classifier and its post-activation qualification.

The following implementation retains the original multi-batch, multi-channel, stride and padding exercise, joining its formerly split function into a complete program. It deliberately supports square kernels, groups1 and dilation1; the general shape formula above is broader than this teaching implementation.

def conv2d_forward(X, W, b, stride=1, padding=0):
    X, W, b = np.asarray(X), np.asarray(W), np.asarray(b)
    if X.ndim != 4 or W.ndim != 4:
        raise ValueError('NCHW input and OIHW weights required')
    N, Ci, H, Wi = X.shape
    Co, Ck, kh, kw = W.shape
    if Ci != Ck or kh != kw or b.shape != (Co,):
        raise ValueError('incompatible channels, kernel or bias')
    if stride < 1 or padding < 0:
        raise ValueError('invalid stride or padding')
    Ho = (H + 2*padding-kh)//stride+1
    Wo = (Wi + 2*padding-kw)//stride+1
    if min(Ho, Wo) < 1:
        raise ValueError('kernel exceeds padded input')
    Xp = np.pad(X, ((0,0),(0,0),(padding,padding),(padding,padding)))
    Y = np.zeros((N,Co,Ho,Wo), dtype=np.result_type(X,W,b))
    for n in range(N):
        for oc in range(Co):
            for i in range(Ho):
                for j in range(Wo):
                    hs, ws = i*stride, j*stride
                    Y[n,oc,i,j] = (Xp[n,:,hs:hs+kh,ws:ws+kw]*W[oc]).sum()+b[oc]
    return Y
 
rng = np.random.default_rng(0)
image = rng.normal(size=(2,3,8,8))
weights = rng.normal(size=(4,3,3,3)); bias4 = np.zeros(4)
assert conv2d_forward(image, weights, bias4, 1, 1).shape == (2,4,8,8)
assert conv2d_forward(image, weights, bias4, 2, 1).shape == (2,4,4,4)
import torch
for stride in (1,2):
    got = conv2d_forward(image, weights, bias4, stride, 1)
    ref = torch.nn.functional.conv2d(torch.tensor(image), torch.tensor(weights),
                                    torch.tensor(bias4), stride=stride, padding=1)
    np.testing.assert_allclose(got, ref.numpy(), rtol=1e-12, atol=1e-12)
x2 = rng.normal(size=(1,2,4,4)); w2 = rng.normal(size=(2,2,3,3))
assert conv2d_forward(x2,w2,np.zeros(2)).shape == (1,2,2,2)
try:
    conv2d_forward(image.transpose(0,2,3,1), weights, bias4)
    raise AssertionError('wrong layout accepted')
except ValueError:
    pass

The inner sum collapses input-channel and spatial axes, not the output-channel axis. Cost is proportional to $NC_{out}H_{out}W_{out}C_{in}k_Hk_W$ multiply-accumulates. im2col unfolds windows into a matrix for GEMM but can materialize many repeated pixels. Frameworks select other algorithms too; the loop is an inspectable reference, not a throughput recommendation.

Receptive fields and downsampling are different questions

Begin with receptive-field size r1 and jump j1. For each layer, update $r'=r+(k_{effective}-1)j$ and $j'=sj$. At stride1/dilation1, L equal k×k layers give $1+L(k-1)$: two3×3 layers see5×5, three see7×7, five see11×11. The general recurrence also handles pooling and strided layers. Theoretical receptive field is possible geometric dependence, not a measured equal contribution from every input pixel.

At constant C and ignoring biases, two3×3 layers use18C² versus25C² for one5×5. Their intermediate nonlinearity changes the represented functions. This is a useful design tradeoff, not proof of a strict function-class superset at any width or guaranteed accuracy superiority. Without the nonlinearity their composition is linear, but its effective large kernel is constrained by the factorization and boundaries.

Strided convolution learns a weighted downsample. Max pooling selects a maximum; gradients go to selected maxima with a tie policy. Average pooling distributes gradients and averages away detail. They are not equivalent operators. Stride2 can save work relative to evaluating a stride1 convolution everywhere, but it is not universally four times faster. Both striding and pooling can alias high frequencies; consider low-pass filtering and evaluate shift sensitivity when alignment matters. Global average pooling removes the flatten-to-dense dependence on spatial size, but the upstream network still has minimum-size and padding assumptions.

For ResNet18's stem, k7,s2 gives112 resolution; the following k3,s2 max pool gives56. A later stride1 k3 convolution does not halve112 to56 by itself. The two named conv weight counts9,408 and36,864 sum46,272, but the intervening pooling must be included in the geometry.

Shared weights require summed gradients

For each output derivative G, accumulate $dW_{oc,ic,a,b} += G_{n,oc,i,j}X^p_{n,ic,is+a,js+b}$, $db_{oc} += G$, and $dX^p_{n,ic,is+a,js+b} += G W_{oc,ic,a,b}$. Overlaps require addition, not assignment. The input derivative is the adjoint of the forward linear map; a transposed-convolution API already implements its own index convention, so do not blindly flip a kernel again. Reduction over examples or pixels also changes gradient scale; no universal smaller-CNN-learning-rate rule follows from sharing.

The exact previously verified scalar-array program remains intact here:

import numpy as np
 
x = np.arange(1., 10.).reshape(3, 3)
k = np.array([[1., 0.], [0., -1.]])
 
def correlate(x, k):
    kh, kw = k.shape
    out = np.empty((x.shape[0]-kh+1, x.shape[1]-kw+1))
    for i in range(out.shape[0]):
        for j in range(out.shape[1]):
            out[i, j] = np.sum(x[i:i+kh, j:j+kw] * k)
    return out
 
y = correlate(x, k)
np.testing.assert_allclose(y, -4.)
dk = np.zeros_like(k); dx = np.zeros_like(x)
for i in range(y.shape[0]):
    for j in range(y.shape[1]):
        dk += y[i,j] * x[i:i+2, j:j+2]
        dx[i:i+2, j:j+2] += y[i,j] * k
 
def loss():
    return .5 * np.sum(correlate(x, k)**2)
 
for array, grad in ((x, dx), (k, dk)):
    for index in np.ndindex(array.shape):
        old = array[index]
        array[index] = old + 1e-5; plus = loss()
        array[index] = old - 1e-5; minus = loss()
        array[index] = old
        assert abs((plus-minus)/2e-5 - grad[index]) < 1e-7
print('output', y.tolist(), 'kernel gradient', dk.tolist())

Now the full solution to the original backward-pass stretch, including multi-channel input, bias, stride and padding. Its differentiation contract is finite real floating-point NumPy arrays with compatible dtypes; the tests use float64. zeros_like inherits the input dtype, so passing integer image/weight arrays does not create a valid floating-point gradient buffer. Convert such inputs before this exercise. Stride and padding are integers, and this reference does not implement complex-valued adjoints:

def conv2d_backward(X, W, G, stride=1, padding=0):
    b = np.zeros(W.shape[0], dtype=W.dtype)
    shape = conv2d_forward(X,W,b,stride,padding).shape
    if G.shape != shape:
        raise ValueError('upstream gradient shape mismatch')
    Xp = np.pad(X, ((0,0),(0,0),(padding,padding),(padding,padding)))
    dXp = np.zeros_like(Xp); dW = np.zeros_like(W)
    db = G.sum(axis=(0,2,3)); kh,kw = W.shape[-2:]
    for n in range(shape[0]):
        for oc in range(shape[1]):
            for i in range(shape[2]):
                for j in range(shape[3]):
                    hs,ws = i*stride,j*stride; g = G[n,oc,i,j]
                    dW[oc] += g*Xp[n,:,hs:hs+kh,ws:ws+kw]
                    dXp[n,:,hs:hs+kh,ws:ws+kw] += g*W[oc]
    dX = dXp if padding == 0 else dXp[:,:,padding:-padding,padding:-padding]
    return dX,dW,db
 
for stride,padding in ((1,0),(1,1),(2,1)):
    xx = rng.normal(size=(1,2,4,4)); ww = rng.normal(size=(2,2,3,3))
    bb = rng.normal(size=2)
    yy = conv2d_forward(xx,ww,bb,stride,padding); gg = rng.normal(size=yy.shape)
    dx,dw,db = conv2d_backward(xx,ww,gg,stride,padding)
    tx = torch.tensor(xx, requires_grad=True); tw = torch.tensor(ww, requires_grad=True)
    tb = torch.tensor(bb, requires_grad=True)
    ty = torch.nn.functional.conv2d(tx,tw,tb,stride=stride,padding=padding)
    (ty*torch.tensor(gg)).sum().backward()
    for a,b in ((dx,tx.grad.numpy()),(dw,tw.grad.numpy()),(db,tb.grad.numpy())):
        np.testing.assert_allclose(a,b,rtol=1e-12,atol=1e-12)
    for arr,grad in ((xx,dx),(ww,dw),(bb,db)):
        for idx in np.ndindex(arr.shape):
            old = arr[idx]; eps = 1e-5
            arr[idx] = old+eps; plus = (conv2d_forward(xx,ww,bb,stride,padding)*gg).sum()
            arr[idx] = old-eps; minus = (conv2d_forward(xx,ww,bb,stride,padding)*gg).sum()
            arr[idx] = old
            assert abs((plus-minus)/(2*eps)-grad[idx]) < 1e-8
print('multichannel forward/backward, stride and padding: pass')

Modern convolution: mechanisms rather than a winner list

The original survey's subjects remain useful, but its universal speedups, app-deployment stories and mixed benchmark configurations do not. Read each cited paper as a particular evaluated system.

• ConvNeXt (2022): modernizes ResNet using a modern training recipe, stage design, depthwise7×7 kernels, inverted bottlenecks, LayerNorm and GELU. Its87.8% headline is a scaled/pretrained result, not a22M-parameter model's generic score. Recipe and architecture ablations must be separated.
• ConvNeXt V2 (2023): co-designs a fully convolutional masked autoencoder and Global Response Normalization, which promotes competition across feature channels. The88.9% headline belongs to the650M-parameter Huge configuration, not every version or an automatic small-model gain.
• MobileNetV4 (2024): Universal Inverted Bottleneck unifies optional depthwise locations; Mobile MQA adds efficient attention in hybrid variants. The family includes convolutional and hybrid models. Its paper reports87.0% at3.8ms on a Pixel8 EdgeTPU for a specific model/protocol, not every phone camera or a CPU measurement here.
• EfficientNetV2 (2021): Fused-MBConv and progressive training address training efficiency; it is not merely EfficientNet's original compound scaling nor universally fastest below100M parameters.
• RepVGG and MobileOne: train linear convolution/normalization branches that can be folded into a simpler inference block. At evaluation, Conv+BN folds as $W'=\gamma W/\sqrt{v+\epsilon}$ and $b'=\beta+\gamma(b-\mu)/\sqrt{v+\epsilon}$. Pad a1×1 kernel to3×3 and sum compatible branches. This works only when no intervening nonlinearity prevents algebraic merging, with frozen BN statistics and matching stride/groups.
• RepLKNet: uses very large depthwise kernels and reparameterization to study large receptive fields; small3×3 kernels are not a universal optimum.
• FastViT: hybrid token mixing with structural reparameterization reduces memory-access costs. EfficientFormerV2 jointly searches architecture and size for mobile latency. Neither paper proves the specific commercial iOS deployment claimed by the old lesson.
• Swin partitions attention into shifted windows; MaxViT combines local-window and dilated-grid attention with convolutional blocks; Hiera simplifies a hierarchical transformer. Calling every one a pure CNN would erase their mechanism differences.

The durable decision variables are label/data scale, resolution, memory, actual deployment operators and measured latency. Convolution offers locality and compatible-grid translation equivariance. Boundary padding, stride and readout complicate invariance; there is no portable5–15-point or3–5× reward for choosing it.

Worked recall, debugging and practice

1. Stem shape and parameters: (1,3,224,224) through3→64,k7,s2,p3 becomes(1,64,112,112), with9,472 parameters including bias.
2. **Valid versus same:**32 with k5,p0,s1 becomes28; p2 gives32. For even kernels, specify asymmetric alignment rather than banning them.
3. Five3×3 layers: receptive field11 at stride1/dilation1. With intervening stride, use the jump recurrence; the shortcut formula no longer applies.
4. **Dense1000 versus conv1000 plus average:**150,529,000 dense parameters including bias versus148,000 conv parameters including bias for k7,3 input channels. This compares differently constrained maps to1000 features, not equal spatial computation.
5. Debug three failures: layout errors need NCHW/range assertions; shifted masks need padding/stride center tracing; conv bias immediately before train-mode BN is usually redundant because BN subtracts the channel mean. The issue is not a fixed accuracy penalty.
6. What changed in ConvNeXt? Both training and architecture. Identify AdamW/augmentation/schedule separately from depthwise7×7, normalization and bottleneck edits; no single universally dominant ingredient follows from a headline score.
7. Explain sharing without jargon: one small detector examines each nearby patch, using the same numbers every time. Training improves that detector from all its uses. This reduces separately learned coefficients but preserves a map of where detections occurred.
8. Backward stretch: use the complete accumulation solution above; deliberately replace += with = in a private experiment and an overlapping-window gradient check fails. A symmetric kernel alone cannot catch reversal.

Revisit shape arithmetic and the backward loop tomorrow; on day7 derive them before looking. The next lesson compares the actual LeNet, AlexNet and residual designs, rather than treating every CNN as the same recipe.

Sources and execution boundary

Primary reading: PyTorch2.9 Conv2d, LeCun et al.1998, VGG, ConvNeXt, ConvNeXtV2, MobileNetV4, EfficientNetV2, RepVGG, MobileOne, RepLKNet, FastViT, EfficientFormerV2.

CPU tests execute the displayed programs and compare framework derivatives and finite differences. They do not reproduce paper benchmarks, CIFAR training, pretrained weights, GPU or mobile timing. Course map · Next: CNN architectures.