Pooling Operations

Introduction

§10.4 closed with a promise: “max-pooling, average-pooling, and even modern alternatives like strided convolutions are variations on a single theme — a stride-S window reducer.” This section pays that off. We will show that the whole family of pooling operators is specified by two choices: the window geometry (which inherits the $O = \lfloor (I + 2P - K)/S \rfloor + 1$ formula from §10.3), and a reducer function applied inside each window. Swap the reducer and you traverse the whole design space: max, average, L^p, stochastic, fractional, mixed, and eventually the all-convolutional replacement.

The Core Insight: Pooling is a deliberately parameter-free, non-linear downsampler. It trades a tiny amount of expressiveness for three concrete wins: smaller feature maps, approximate translation invariance, and a bigger effective receptive field for every downstream layer. The famous “complex cells” of Hubel & Wiesel that §10.1 introduced map almost directly onto this operator.

You already met pooling informally — the animation in §10.2's CNN2DConvolutionVisualizer halves a 7×7 feature map with a 2×2 max-pool, and every CNN in Chapter 11 (LeNet, AlexNet, VGG) will lean on it. Here we open the hood: what each variant does mathematically, how gradients route through it, and why modern architectures have been steadily replacing it with strided convolutions.

Learning Objectives

After working through this section, you will be able to:

1. Implement max-pool and average-pool from scratch in NumPy and reproduce PyTorch's output byte-for-byte.
2. Predict output shapes using $O = \lfloor (I + 2P - K)/S \rfloor + 1$ for any kernel, stride, and padding.
3. Derive the backward pass for both max and average pool and verify it against torch.autograd.
4. Explain Global Average Pooling and why it replaced fully-connected heads in GoogLeNet, ResNet, and every modern classifier.
5. Use nn.AdaptiveAvgPool2d to make a network accept variable input sizes.
6. Argue — with citations — why ResNet and ConvNeXt eliminated max-pool in favour of stride-2 convolutions.

The Four Jobs Pooling Does

Before any maths, let us be honest about why pooling exists. It is not a decorative layer: it does four concrete jobs, and if you remove it you must replace it with something that does the same jobs. Identifying these jobs lets us reason about what the modern stride-2 conv actually replaces.

Max Pooling

The most common pooling variant, and the one you will see in AlexNet, VGG, and every textbook diagram. For a window $W_{ij}$ of size $K \times K$ centred at output position $(i,j)$:

$y[i,j] = \max_{(m,n) \in W_{ij}} x[m,n]$

Intuition: “report the loudest activation in this neighbourhood.” Under a ReLU activation (§5.3), pixel values represent how strongly a learned feature fires — so max-pool asks whether the feature is present, not exactly where. This is Boureau, Ponce & LeCun (2010) [3]'s formal argument for why max-pool works well on sparse codes: when most cells are zero, the max of a window effectively detects presence.

Hand-Computation on a 4×4 Feature Map

Let us do a full worked example on the smallest non-trivial input. The 4×4 grid below is exactly the one you will see animated in the Pooling Playground further down, so you can match each hand-computed value against the live widget.

Result: [[6, 4], [8, 6]] — a 2×2 output from a 4×4 input, exactly as the master formula predicts: $(4 - 2)/2 + 1 = 2$.

Max Pool in Pure Python

     c0  c1  c2  c3
r0    1   3   2   4
r1    5   6   1   2
r2    7   2   3   1
r3    4   8   5   6

     c0   c1
r0    0    0
r1    0    0  (before loop)

[[6., 4.],
 [8., 6.]]

1import numpy as np
2
3# 4x4 feature map — same values shown in the Pooling Visualizer above
4x = np.array([
5    [1, 3, 2, 4],
6    [5, 6, 1, 2],
7    [7, 2, 3, 1],
8    [4, 8, 5, 6],
9], dtype=float)
10
11def max_pool2d(x, K, S):
12    """2-D max-pool with kernel size K and stride S, no padding."""
13    H, W = x.shape
14    out_H = (H - K) // S + 1
15    out_W = (W - K) // S + 1
16    y = np.zeros((out_H, out_W))
17    for i in range(out_H):
18        for j in range(out_W):
19            window = x[i*S : i*S + K,
20                       j*S : j*S + K]
21            y[i, j] = window.max()
22    return y
23
24print(max_pool2d(x, K=2, S=2))
25# [[6. 4.]
26#  [8. 6.]]

Max Pool in PyTorch

[[[[6., 4.],
    [8., 6.]]]]

1import torch
2import torch.nn as nn
3import torch.nn.functional as F
4
5# Same 4x4 feature map, now as an (N, C, H, W) tensor
6x = torch.tensor([[[
7    [1., 3., 2., 4.],
8    [5., 6., 1., 2.],
9    [7., 2., 3., 1.],
10    [4., 8., 5., 6.],
11]]])  # shape: (1, 1, 4, 4)
12
13# Option A — the Module form (preferred inside nn.Sequential)
14pool = nn.MaxPool2d(kernel_size=2, stride=2)
15y1 = pool(x)
16
17# Option B — the functional form (preferred when you don't need state)
18y2 = F.max_pool2d(x, kernel_size=2, stride=2)
19
20print(y1.squeeze())
21# tensor([[6., 4.],
22#         [8., 6.]])
23
24assert torch.equal(y1, y2)                  # same numbers
25assert sum(p.numel() for p in pool.parameters()) == 0   # ZERO learnable params

Average Pooling

The other classical reducer. For the same window $W_{ij}$:

$y[i,j] = \frac{1}{|W_{ij}|} \sum_{(m,n) \in W_{ij}} x[m,n]$

LeCun et al.'s original LeNet-5 [1] actually used a learnable scaled average pool (one coefficient and one bias per feature map), but the modern convention settled on plain unweighted averaging. Average-pool preserves the total energy of the window and, unlike max, routes gradient to every input cell — a property that will matter a great deal in the Backward Pass section below.

Average Pool in Pure Python

[[3.75, 2.25],
 [5.25, 3.75]]

1import numpy as np
2
3x = np.array([
4    [1, 3, 2, 4],
5    [5, 6, 1, 2],
6    [7, 2, 3, 1],
7    [4, 8, 5, 6],
8], dtype=float)
9
10def avg_pool2d(x, K, S):
11    """Identical to max_pool2d but swaps .max() for .mean()."""
12    H, W = x.shape
13    out_H = (H - K) // S + 1
14    out_W = (W - K) // S + 1
15    y = np.zeros((out_H, out_W))
16    for i in range(out_H):
17        for j in range(out_W):
18            window = x[i*S : i*S + K, j*S : j*S + K]
19            y[i, j] = window.mean()
20    return y
21
22print(avg_pool2d(x, K=2, S=2))
23# [[3.75 2.25]
24#  [5.25 3.75]]

Average Pool in PyTorch

1import torch
2import torch.nn as nn
3import torch.nn.functional as F
4
5x = torch.tensor([[[
6    [1., 3., 2., 4.],
7    [5., 6., 1., 2.],
8    [7., 2., 3., 1.],
9    [4., 8., 5., 6.],
10]]])
11
12# Module form
13avg = nn.AvgPool2d(kernel_size=2, stride=2)
14y_avg = avg(x)
15
16# Functional form
17y_func = F.avg_pool2d(x, kernel_size=2, stride=2)
18
19print(y_avg.squeeze())
20# tensor([[3.7500, 2.2500],
21#         [5.2500, 3.7500]])
22
23# Important subtlety: padding + count_include_pad
24x_padded_demo = torch.tensor([[[[1., 2.], [3., 4.]]]])
25print(F.avg_pool2d(x_padded_demo, 2, 1, padding=1,
26                   count_include_pad=True))   # divides by 4 (incl. zeros)
27print(F.avg_pool2d(x_padded_demo, 2, 1, padding=1,
28                   count_include_pad=False))  # divides only by real cells

The Master Pooling Formula

Pooling is a convolution that reduces instead of weight-summing, so it inherits §10.3's master formula. Set dilation to 1 and you have:

$O = \left\lfloor \frac{I + 2P - K}{S} \right\rfloor + 1$

with the same symbol dictionary: $I$ = input size, $K$ = kernel size, $P$ = padding per side, $S$ = stride. The formula applies identically to max, avg, L^p, stochastic and mixed pooling — only the reducer changes, never the geometry.

1import torch
2import torch.nn.functional as F
3import math
4
5def pool_out(I, K, S, P):
6    """Master pooling formula (same as conv with dilation=1, §10.3)."""
7    return (I + 2*P - K) // S + 1
8
9# Verify against PyTorch for 20 random valid configurations
10torch.manual_seed(0)
11for _ in range(20):
12    I = torch.randint(4, 65, (1,)).item()
13    K = torch.randint(2, 8, (1,)).item()
14    S = torch.randint(1, K+1, (1,)).item()
15    P = torch.randint(0, K//2 + 1, (1,)).item()
16    if 2*P >= K:  # PyTorch rejects pad >= kernel/2
17        continue
18
19    x = torch.randn(1, 1, I, I)
20    y = F.max_pool2d(x, kernel_size=K, stride=S, padding=P)
21
22    predicted = pool_out(I, K, S, P)
23    actual = y.shape[-1]
24    assert predicted == actual, (I, K, S, P, predicted, actual)
25    print(f"I={I:2d}  K={K}  S={S}  P={P}  →  {predicted}×{predicted}  ✓")

Interactive: Pooling Playground

Switch between max and average modes, step through the windows, and watch how the output builds up cell by cell. The 4×4 grid matches the hand-computation above; the 6×6 grid lets you see a longer sliding pattern.

Max vs Average — When to Use Which

The choice is not arbitrary; there is a real theoretical and empirical literature. Boureau, Ponce & LeCun (2010) [3] analysed pooling on sparse codes under simple Bernoulli activation models, and Scherer, Müller & Behnke (2010) [4] ran a direct bake-off on object-recognition benchmarks. Their findings cohere:

Worked Example — Continuing the 7×7 from §10.2

§10.2 animated a 7×7 feature map sliding under a 3×3 kernel inside CNN2DConvolutionVisualizer. We reuse exactly that matrix here so you can track specific cells through the conv → pool pipeline end-to-end.

     c0 c1 c2 c3 c4 c5 c6
r0    3  3  2  1  0  2  1
r1    0  0  1  3  1  0  2
r2    3  1  2  2  3  1  0
r3    2  0  0  2  2  3  1
r4    2  0  0  0  1  2  2
r5    1  3  2  1  0  1  3
r6    0  2  1  3  2  0  1

     c0  c1  c2
r0    3   3   2
r1    3   2   3
r2    3   2   2

      c0    c1    c2
r0   1.50  1.75  0.75
r1   1.50  1.50  2.25
r2   1.50  0.75  1.00

1import numpy as np
2
3# 7x7 feature map — SAME matrix used in the §10.2 CNN2DConvolutionVisualizer
4x = np.array([
5    [3, 3, 2, 1, 0, 2, 1],
6    [0, 0, 1, 3, 1, 0, 2],
7    [3, 1, 2, 2, 3, 1, 0],
8    [2, 0, 0, 2, 2, 3, 1],
9    [2, 0, 0, 0, 1, 2, 2],
10    [1, 3, 2, 1, 0, 1, 3],
11    [0, 2, 1, 3, 2, 0, 1],
12], dtype=float)
13
14# 2x2 max pool, stride 2 — the 'VGG/AlexNet halving block'
15def max_pool2d(x, K, S):
16    H, W = x.shape
17    out_H = (H - K) // S + 1
18    out_W = (W - K) // S + 1
19    y = np.zeros((out_H, out_W))
20    for i in range(out_H):
21        for j in range(out_W):
22            y[i, j] = x[i*S:i*S+K, j*S:j*S+K].max()
23    return y
24
25print("max-pool 2x2 stride 2 →")
26print(max_pool2d(x, 2, 2))
27# [[3. 3. 2.]
28#  [3. 2. 3.]
29#  [3. 2. 2.]]
30
31# Same function, same input, with .mean() instead of .max():
32def avg_pool2d(x, K, S):
33    H, W = x.shape
34    out_H = (H - K) // S + 1
35    out_W = (W - K) // S + 1
36    y = np.zeros((out_H, out_W))
37    for i in range(out_H):
38        for j in range(out_W):
39            y[i, j] = x[i*S:i*S+K, j*S:j*S+K].mean()
40    return y
41
42print("avg-pool 2x2 stride 2 →")
43print(avg_pool2d(x, 2, 2))
44# [[1.50 1.75 0.75]
45#  [1.50 1.50 2.25]
46#  [1.50 0.75 1.00]]

Global Average Pooling

If you set the pooling window equal to the full spatial extent of the feature map, the output is one number per channel — a C-dimensional vector. This is Global Average Pooling (GAP), introduced by Lin, Chen & Yan in Network-in-Network (2014) [7] and immediately adopted by GoogLeNet (2015) [8] and then every ResNet-style architecture [9].

$y_c = \frac{1}{H \cdot W} \sum_{i=1}^{H} \sum_{j=1}^{W} x_{c,i,j}$

Why GAP replaced the giant FC head

The AlexNet (2012) head looked like Flatten(6144) → FC(4096) → FC(4096) → FC(1000) — roughly 58 million parameters just in the head. GoogLeNet swapped this for GAP(1024) → FC(1000), i.e. ~1 million parameters. Same accuracy on ImageNet, with vastly less overfitting risk. Four structural wins:

1. Zero parameters in the pool itself.
2. Works for any input size. GAP always produces a C-dim vector regardless of H and W, so the classifier downstream is size-agnostic.
3. Acts as a structural regulariser. Lin et al. (2014) [7] note that GAP forces each feature map to correspond to a class, because the classifier can only average over it.
4. Localisation comes for free. Class-Activation Maps (Zhou et al. 2016, [Ref 9 of §10.6 foreshadow]) use the fact that GAP is a weighted sum to highlight which spatial cells drove the prediction.

tensor([[3.75, 2.00, 5.00]])

1import torch
2import torch.nn as nn
3import torch.nn.functional as F
4
5# Imagine the final feature map of a classifier: 3 channels, 4x4 spatial
6x = torch.tensor([[
7    # channel 0 — same 4x4 matrix we've been using
8    [[1., 3., 2., 4.],
9     [5., 6., 1., 2.],
10     [7., 2., 3., 1.],
11     [4., 8., 5., 6.]],
12    # channel 1 — uniform response
13    [[2., 2., 2., 2.],
14     [2., 2., 2., 2.],
15     [2., 2., 2., 2.],
16     [2., 2., 2., 2.]],
17    # channel 2 — uniform high response
18    [[5., 5., 5., 5.],
19     [5., 5., 5., 5.],
20     [5., 5., 5., 5.],
21     [5., 5., 5., 5.]],
22]])  # shape: (1, 3, 4, 4)
23
24# Three equivalent ways to spell GAP — pick the one that reads best
25gap_module = nn.AdaptiveAvgPool2d(output_size=1)       # (1) the canonical module
26y1 = gap_module(x).flatten(1)                          # (1, 3)
27
28y2 = F.adaptive_avg_pool2d(x, 1).flatten(1)            # (2) functional
29y3 = x.mean(dim=(2, 3))                                # (3) pure tensor ops, same result
30
31print(y1)   # tensor([[3.7500, 2.0000, 5.0000]])
32assert torch.allclose(y1, y2)
33assert torch.allclose(y1, y3)

Adaptive Pooling

Adaptive pooling is the generalisation that enables GAP: you specify the output size, and PyTorch picks the kernel and stride schedule for you. This removes the hard-coded 7×7 feature-map assumption that plagued early CNN code and is the reason torchvision models accept any reasonable input resolution.

For an input of height $H_{\text{in}}$ and a target height $H_{\text{out}}$, PyTorch's rule for output cell $i$ is:

$\text{start}_i = \left\lfloor \frac{i \cdot H_{\text{in}}}{H_{\text{out}}} \right\rfloor, \quad \text{end}_i = \left\lceil \frac{(i+1) \cdot H_{\text{in}}}{H_{\text{out}}} \right\rceil$

Two consequences: bins can have different widths when $H_{\text{in}}$ is not a multiple of $H_{\text{out}}$, and adjacent bins can overlap by one cell. Adaptive pool is therefore not a plain “pick the right kernel” reduction.

1import torch
2import torch.nn.functional as F
3
4# Feed the SAME network feature sizes that don't divide evenly into 7
5# AdaptiveAvgPool2d has to invent a variable-size binning schedule on the fly.
6for H in [7, 8, 10, 14, 17, 224]:
7    x = torch.arange(H*H, dtype=torch.float32).view(1, 1, H, H)
8    y = F.adaptive_avg_pool2d(x, 7)
9    print(f"input {H:3d}x{H:3d}  →  output {y.shape[-2]}x{y.shape[-1]}")
10# input   7x  7  →  output 7x7
11# input   8x  8  →  output 7x7
12# input  10x 10  →  output 7x7
13# input  14x 14  →  output 7x7
14# input  17x 17  →  output 7x7
15# input 224x224  →  output 7x7
16
17# The bin boundaries for output cell i are:
18#   start = floor(i  *  H_in / H_out)
19#   end   = ceil((i+1) * H_in / H_out)
20# Bins can have different widths when H_in is not a multiple of H_out.
21def adaptive_bins(H_in, H_out):
22    import math
23    return [(math.floor(i*H_in/H_out), math.ceil((i+1)*H_in/H_out))
24            for i in range(H_out)]
25
26print(adaptive_bins(10, 7))
27# [(0, 2), (1, 3), (2, 5), (4, 6), (5, 8), (7, 9), (8, 10)]

The Backward Pass — How Gradients Route

Pooling has no parameters, but it still has a non-trivial backward pass: an upstream gradient $\partial L / \partial y$ must be routed back to the input cells. How it is routed depends entirely on the reducer.

Max pool — the argmax router

The derivative of $y = \max(x_1, \ldots, x_{K^2})$ with respect to $x_k$ is 1 if $x_k$ is the max, and 0 otherwise (technically a sub-gradient since max is non-smooth at ties, but ties are measure-zero in practice). Therefore:

$\frac{\partial L}{\partial x[m,n]} = \frac{\partial L}{\partial y[i,j]} \quad \text{if } (m,n) = \arg\max_{W_{ij}} x, \text{ else } 0$

The consequence is striking: in a window of size $K \times K$, exactly one cell receives gradient and $K^2 - 1$ receive none. For our 4×4 example with 2×2 pools, 12 of 16 input cells get nothing.

Average pool — the uniform splitter

Average pool is a linear operator, so its derivative is constant: each input cell of each window gets $1/K^2$ of the upstream gradient. No zero-gradient cells, no dead zones.

$\frac{\partial L}{\partial x[m,n]} = \frac{1}{K^2} \sum_{W_{ij} \ni (m,n)} \frac{\partial L}{\partial y[i,j]}$

Backward Pass in Pure Python

[[0. 0. 0. 2.]
 [0. 1. 0. 0.]
 [0. 0. 0. 0.]
 [0. 3. 0. 4.]]

[[0.25 0.25 0.50 0.50]
 [0.25 0.25 0.50 0.50]
 [0.75 0.75 1.00 1.00]
 [0.75 0.75 1.00 1.00]]

1import numpy as np
2
3# FORWARD — same 4x4 input, 2x2 stride-2 max pool
4x = np.array([
5    [1, 3, 2, 4],
6    [5, 6, 1, 2],
7    [7, 2, 3, 1],
8    [4, 8, 5, 6],
9], dtype=float)
10
11def max_pool_forward(x, K, S):
12    H, W = x.shape
13    out_H = (H - K) // S + 1
14    out_W = (W - K) // S + 1
15    y = np.zeros((out_H, out_W))
16    argmax = np.zeros((out_H, out_W, 2), dtype=int)   # cache WHERE the max was
17    for i in range(out_H):
18        for j in range(out_W):
19            window = x[i*S:i*S+K, j*S:j*S+K]
20            y[i, j] = window.max()
21            # np.unravel_index turns a flat argmax into (row, col) within the window
22            di, dj = np.unravel_index(window.argmax(), window.shape)
23            argmax[i, j] = [i*S + di, j*S + dj]        # store in INPUT coordinates
24    return y, argmax
25
26# BACKWARD — route upstream gradient to the argmax cell of each window
27def max_pool_backward(dL_dy, argmax, x_shape):
28    dL_dx = np.zeros(x_shape)
29    out_H, out_W = dL_dy.shape
30    for i in range(out_H):
31        for j in range(out_W):
32            r, c = argmax[i, j]
33            dL_dx[r, c] += dL_dy[i, j]   # += handles the case two windows share the same max
34    return dL_dx
35
36# AVERAGE POOL — gradient splits uniformly over all K^2 cells
37def avg_pool_backward(dL_dy, x_shape, K, S):
38    dL_dx = np.zeros(x_shape)
39    out_H, out_W = dL_dy.shape
40    for i in range(out_H):
41        for j in range(out_W):
42            dL_dx[i*S:i*S+K, j*S:j*S+K] += dL_dy[i, j] / (K * K)
43    return dL_dx
44
45# Let the upstream gradient be [[1,2],[3,4]] to make tracing easy
46y, argmax = max_pool_forward(x, 2, 2)
47dL_dy = np.array([[1, 2], [3, 4]], dtype=float)
48
49print("max-pool backward dL/dx:")
50print(max_pool_backward(dL_dy, argmax, x.shape))
51# [[0. 0. 0. 2.]
52#  [0. 1. 0. 0.]
53#  [0. 0. 0. 0.]
54#  [0. 3. 0. 4.]]
55
56print("avg-pool backward dL/dx:")
57print(avg_pool_backward(dL_dy, x.shape, 2, 2))
58# [[0.25 0.25 0.50 0.50]
59#  [0.25 0.25 0.50 0.50]
60#  [0.75 0.75 1.00 1.00]
61#  [0.75 0.75 1.00 1.00]]

Verification Against torch.autograd

1import torch
2import torch.nn.functional as F
3import numpy as np
4
5x_t = torch.tensor([
6    [1., 3., 2., 4.],
7    [5., 6., 1., 2.],
8    [7., 2., 3., 1.],
9    [4., 8., 5., 6.],
10], requires_grad=True).unsqueeze(0).unsqueeze(0)   # (1,1,4,4)
11
12# MAX-POOL — let autograd compute the gradient
13y = F.max_pool2d(x_t, kernel_size=2, stride=2)
14dL_dy = torch.tensor([[[[1., 2.], [3., 4.]]]])
15y.backward(dL_dy)
16print("autograd (max):")
17print(x_t.grad.squeeze())
18# tensor([[0., 0., 0., 2.],
19#         [0., 1., 0., 0.],
20#         [0., 0., 0., 0.],
21#         [0., 3., 0., 4.]])
22
23# AVG-POOL — reset grads first
24x_t.grad = None
25y = F.avg_pool2d(x_t, kernel_size=2, stride=2)
26y.backward(dL_dy)
27print("autograd (avg):")
28print(x_t.grad.squeeze())
29# tensor([[0.2500, 0.2500, 0.5000, 0.5000],
30#         [0.2500, 0.2500, 0.5000, 0.5000],
31#         [0.7500, 0.7500, 1.0000, 1.0000],
32#         [0.7500, 0.7500, 1.0000, 1.0000]])
33
34# Both match the NumPy reference implementation byte-for-byte.

Variants — Stochastic, L^p, Mixed, Fractional

The 2013–2015 literature explored many alternative reducers. None dethroned max and avg as defaults, but several are useful tools and they make the “single family, different reducer” framing concrete.

1import numpy as np
2
3# Stochastic pooling (Zeiler & Fergus, 2013) — sample a cell with probability
4# proportional to its value, not argmax
5def stochastic_pool2d(x, K, S, rng):
6    H, W = x.shape
7    out_H = (H - K) // S + 1
8    out_W = (W - K) // S + 1
9    y = np.zeros((out_H, out_W))
10    for i in range(out_H):
11        for j in range(out_W):
12            w = x[i*S:i*S+K, j*S:j*S+K].flatten()
13            probs = w / w.sum()           # requires non-negative inputs (e.g., post-ReLU)
14            y[i, j] = rng.choice(w, p=probs)
15    return y
16
17# L^p pooling (Sermanet, Chintala, LeCun 2013) — interpolates between avg (p=1) and max (p→∞)
18def lp_pool2d(x, K, S, p):
19    H, W = x.shape
20    out_H = (H - K) // S + 1
21    out_W = (W - K) // S + 1
22    y = np.zeros((out_H, out_W))
23    for i in range(out_H):
24        for j in range(out_W):
25            w = x[i*S:i*S+K, j*S:j*S+K]
26            y[i, j] = (np.mean(w**p)) ** (1/p)
27    return y
28
29# Mixed pooling (Yu, Wang, Chen, Wei 2014) — random (or learnable) convex mix
30def mixed_pool2d(x, K, S, lam):
31    H, W = x.shape
32    out_H = (H - K) // S + 1
33    out_W = (W - K) // S + 1
34    y = np.zeros((out_H, out_W))
35    for i in range(out_H):
36        for j in range(out_W):
37            w = x[i*S:i*S+K, j*S:j*S+K]
38            y[i, j] = lam * w.max() + (1.0 - lam) * w.mean()
39    return y

Pool vs Strided Convolution

The central architectural question of the late 2010s: if pooling just halves the spatial resolution, why not use a stride-2 convolution that halves it and learns what to keep? Springenberg, Dosovitskiy, Brox & Riedmiller (2015) [10] made this case explicitly in “Striving for Simplicity: The All Convolutional Net” and showed that replacing every max-pool with a stride-2 conv reaches equal or better accuracy on CIFAR-10 and ImageNet.

1import torch
2import torch.nn as nn
3
4# Same goal — halve a (C=64, 32x32) feature map in both dimensions.
5pool = nn.MaxPool2d(kernel_size=2, stride=2)           # 0 parameters
6sconv = nn.Conv2d(64, 64, kernel_size=3,
7                  stride=2, padding=1)                 # 64·64·3·3 + 64 = 37,000 parameters
8
9x = torch.randn(1, 64, 32, 32)
10print("pool  output:", pool(x).shape)    # torch.Size([1, 64, 16, 16])
11print("sconv output:", sconv(x).shape)   # torch.Size([1, 64, 16, 16])
12
13pool_params = sum(p.numel() for p in pool.parameters())
14sconv_params = sum(p.numel() for p in sconv.parameters())
15print(f"pool params:  {pool_params:>6d}")   #       0
16print(f"sconv params: {sconv_params:>6d}")  #  36,928
17
18# Same spatial downsampling, but sconv LEARNS what to keep.
19# Springenberg et al. (2015) showed that replacing every max-pool with
20# a stride-2 conv reaches equal or better accuracy on CIFAR-10 / ImageNet.

Why Pooling Declined in Modern CNNs

Post-2016, every flagship architecture has moved max-pool out of the spine:

1. ResNet (He et al. 2016) [9] uses stride-2 3×3 conv inside every transition block and keeps only a single MaxPool2d(3, 2, 1) after the stem. Residual connections mean the signal can skip past the learned downsampler, so the network never has to choose between preserving information and reducing resolution.
2. DenseNet (2017) uses average pool in transition blocks, not max.
3. Vision Transformer (Dosovitskiy et al. 2021) [14] contains no pooling at all. Tokenisation via a 16×16 stride-16 patch embedding simultaneously reduces resolution and learns a linear projection.
4. ConvNeXt (Liu et al. 2022) [15] — the explicitly ResNet-modernising architecture — uses stride-2 depthwise convs for downsampling inside stages and a “patchify stem” of stride-4 conv up front. Zero max-pool.

What has survived everywhere is global average pool as the head — cheap, learnable-free, size-agnostic, class-activation-map-friendly.

Takeaway: Pooling in the body of the network is an optional inductive-bias choice that trades parameters for regularisation. Pooling at the head (GAP) is nearly universal and there is no sign of it going away.

Design Patterns in Real Architectures

Summary

Commit these to memory

1. Output-shape formula: $O = \lfloor (I + 2P - K)/S \rfloor + 1$ (same as conv with dilation = 1).
2. Max-pool gradient rule: one cell per window receives the full upstream gradient; the rest receive zero.
3. Avg-pool gradient rule: every cell receives $1/K^2$ of the upstream gradient.
4. GAP = AdaptiveAvgPool2d(1) = .mean(dim=(2,3)). Three spellings for the modern CNN classifier head.
5. Stride-2 conv can replace max-pool at the cost of $K^2 C^2$ parameters per layer. Modern architectures prefer this when data is plentiful.

Exercises

Conceptual

1. Given a (1, 3, 28, 28) input, what is the output shape after nn.MaxPool2d(kernel_size=3, stride=2, padding=1)? Show your use of the master formula.
2. Explain in one sentence why max-pool's backward pass creates zero-gradient cells but avg-pool's does not.
3. A classifier is trained on 224×224 images but deployed on 384×384 images. Which single layer change lets the model work at the new resolution without retraining? (Hint: not the conv layers.)
4. You observe that every cell of a particular feature map eventually becomes an argmax of at least one window during training. Is max-pool still problematic for gradient flow here? Why or why not?
5. Why did Springenberg et al. (2015) [10] get a parameter increase when replacing pool with stride-2 conv, yet their total model was still competitive? (Hint: compare to the parameters they could remove elsewhere.)

Coding

1. Extend max_pool2d to support padding. Pad with $-\infty$ (not 0) so the max is unaffected by the padded cells, and verify against F.max_pool2d(x, K, S, padding=P).
2. Implement max_unpool2d: given the argmax indices cached during the forward pass, reconstruct a sparse tensor where each argmax cell holds the corresponding output value. This is the core operation of the Zeiler & Fergus (2014) deconvnet visualisation. Check your implementation against PyTorch's nn.MaxUnpool2d.
3. Replace every max-pool in a small VGG-style network with a stride-2 3×3 conv (Springenberg et al. [10]). Train both on CIFAR-10 for 10 epochs and compare final accuracy and training curves.
4. Implement the bin schedule in adaptive_bins(H_in, H_out) and feed your output into a manual NumPy avg-pool. Verify the result matches F.adaptive_avg_pool2d byte-for-byte on 20 random input sizes.

Challenge

Reproduce the GoogLeNet head switch. Take a small pre-trained CNN with a dense classifier head (e.g. Flatten → FC(4096) → FC(num_classes)), replace the head with AdaptiveAvgPool2d(1) → Flatten → FC(C, num_classes), and fine-tune for a few epochs. Report (a) the parameter reduction, (b) the change in validation accuracy, and (c) a visualisation of the Class Activation Map produced by the GAP head using the technique of Zhou et al. (2016). This single experiment is the clearest demonstration of why GAP replaced dense heads across nearly every vision architecture from 2015 onwards.

References

1. LeCun, Y., Bottou, L., Bengio, Y., & Haffner, P. (1998). Gradient-Based Learning Applied to Document Recognition. Proceedings of the IEEE 86(11), 2278–2324. DOI:10.1109/5.726791. (First large-scale use of subsampling layers in LeNet-5.)
2. Goodfellow, I., Bengio, Y., & Courville, A. (2016). Deep Learning, §9.3 “Pooling”. MIT Press. https://www.deeplearningbook.org/ (Canonical textbook treatment of pooling as approximate invariance.)
3. Boureau, Y.-L., Ponce, J., & LeCun, Y. (2010). A Theoretical Analysis of Feature Pooling in Visual Recognition. Proc. 27th ICML, 111–118. (Bernoulli-activation analysis of max vs avg — cited for the sparse-features argument.)
4. Scherer, D., Müller, A., & Behnke, S. (2010). Evaluation of Pooling Operations in Convolutional Architectures for Object Recognition. ICANN, LNCS 6354, 92–101. Springer. (Empirical NORB / CIFAR comparison.)
5. Zeiler, M. D., & Fergus, R. (2013). Stochastic Pooling for Regularization of Deep Convolutional Neural Networks. ICLR. arXiv:1301.3557.
6. Sermanet, P., Chintala, S., & LeCun, Y. (2013). Convolutional Neural Networks Applied to House Numbers Digit Classification. ICPR. arXiv:1204.3968. (Introduces L^p-pooling in the SVHN classifier.)
7. Lin, M., Chen, Q., & Yan, S. (2014). Network In Network. ICLR. arXiv:1312.4400. (Introduces Global Average Pooling.)
8. Szegedy, C., Liu, W., Jia, Y., Sermanet, P., Reed, S., Anguelov, D., Erhan, D., Vanhoucke, V., & Rabinovich, A. (2015). Going Deeper with Convolutions (GoogLeNet). CVPR. arXiv:1409.4842. (GAP replacing giant FC head.)
9. He, K., Zhang, X., Ren, S., & Sun, J. (2016). Deep Residual Learning for Image Recognition. CVPR. arXiv:1512.03385. (Stride-2 conv + GAP head in ResNet.)
10. Springenberg, J. T., Dosovitskiy, A., Brox, T., & Riedmiller, M. (2015). Striving for Simplicity: The All Convolutional Net. ICLR workshop. arXiv:1412.6806. (Empirical case that stride-2 conv ≥ max-pool.)
11. Graham, B. (2014). Fractional Max-Pooling. arXiv:1412.6071.
12. Yu, D., Wang, H., Chen, P., & Wei, Z. (2014). Mixed Pooling for Convolutional Neural Networks. RSKT, LNCS 8818. Springer. (Random and learnable convex combination of max and average.)
13. Ronneberger, O., Fischer, P., & Brox, T. (2015). U-Net: Convolutional Networks for Biomedical Image Segmentation. MICCAI. arXiv:1505.04597. (Max-pool in the encoder path, inverted by transposed convolution in the decoder — motivating §10.6.)
14. Dosovitskiy, A. et al. (2021). An Image is Worth 16×16 Words: Transformers for Image Recognition at Scale (ViT). ICLR. arXiv:2010.11929. (First flagship vision architecture with no pooling at all.)
15. Liu, Z., Mao, H., Wu, C.-Y., Feichtenhofer, C., Darrell, T., & Xie, S. (2022). A ConvNet for the 2020s (ConvNeXt). CVPR. arXiv:2201.03545. (Modern ResNet-style CNN with stride-2 depthwise downsampling, no max-pool.)
16. PyTorch documentation, torch.nn.MaxPool2d, AvgPool2d, AdaptiveAvgPool2d, AdaptiveMaxPool2d, FractionalMaxPool2d, functional.max_pool2d, functional.avg_pool2d, functional.adaptive_avg_pool2d. https://pytorch.org/docs/stable/nn.html (Authoritative reference forcount_include_pad, ceil_mode, and adaptive-pool bin schedule.)

In the next section we tackle the inverse of every downsampling operation we have seen so far — transposed convolution. It is the upsampler that lets decoders (U-Net), generators (DCGAN), and super-resolution networks grow feature maps back up to image resolution.