Convolutional Networks, and Why They Suit Images

import numpy as np

image = np.array([
    [0, 0, 0, 9, 9, 9],
    [0, 0, 0, 9, 9, 9],
    [0, 0, 0, 9, 9, 9],
    [0, 0, 0, 9, 9, 9],
    [0, 0, 0, 9, 9, 9],
    [0, 0, 0, 9, 9, 9],
], dtype=float)

# A vertical-edge detector: dark on the left, bright on the right.
kernel = np.array([
    [-1, 0, 1],
    [-1, 0, 1],
    [-1, 0, 1],
], dtype=float)


def convolve(image, kernel):
    kh, kw = kernel.shape
    out_h, out_w = image.shape[0] - kh + 1, image.shape[1] - kw + 1
    out = np.zeros((out_h, out_w))
    for row in range(out_h):
        for column in range(out_w):
            patch = image[row:row + kh, column:column + kw]
            out[row, column] = (patch * kernel).sum()
    return out


result = convolve(image, kernel)
print('image')
print(image.astype(int))
print()
print('kernel')
print(kernel.astype(int))
print()
print('output')
print(result.astype(int))
print()
print('The edge is at column 2 of the output and nowhere else. The')
print('kernel has nine weights and it found that edge at every one')
print('of the four rows, because the SAME nine weights were')
print('applied at all four positions. That is weight sharing, and')
print('it is the entire reason this works on images.')

import numpy as np


def convolve(image, kernel):
    kh, kw = kernel.shape
    out = np.zeros((image.shape[0] - kh + 1, image.shape[1] - kw + 1))
    for row in range(out.shape[0]):
        for column in range(out.shape[1]):
            out[row, column] = (image[row:row + kh,
                                      column:column + kw] * kernel).sum()
    return out


# A bright square on a dark background.
image = np.zeros((7, 7))
image[2:5, 2:5] = 9.0

kernels = {
    'vertical edge': np.array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]], float),
    'horizontal edge': np.array([[-1, -1, -1], [0, 0, 0], [1, 1, 1]], float),
    'blur': np.full((3, 3), 1 / 9),
    'sharpen': np.array([[0, -1, 0], [-1, 5, -1], [0, -1, 0]], float),
    'identity': np.array([[0, 0, 0], [0, 1, 0], [0, 0, 0]], float),
}

print('input')
print(image.astype(int))
for name, kernel in kernels.items():
    print()
    print(f'{name}:')
    print(convolve(image, kernel).round(1))
print()
print('Five kernels, one image, five completely different outputs.')
print('In a real network nobody writes these down - they are')
print('weights, and the first layer of a trained image model')
print('reliably learns something very close to the edge detectors')
print('above.')

import pandas as pd


def output_size(input_size, kernel, stride=1, padding=0):
    return (input_size + 2 * padding - kernel) // stride + 1


rows = []
for kernel, stride, padding, note in [
        (3, 1, 0, 'shrinks by 2 every layer'),
        (3, 1, 1, '"same" padding - size preserved'),
        (5, 1, 2, 'same, with a wider view'),
        (3, 2, 1, 'halves the size - the usual downsample'),
        (7, 2, 3, 'halves it, ResNet stem'),
        (1, 1, 0, 'pointwise - mixes channels, not pixels'),
        (2, 2, 0, 'max pooling, the classic')]:
    rows.append({'kernel': kernel, 'stride': stride, 'padding': padding,
                 '224 ->': output_size(224, kernel, stride, padding),
                 '32 ->': output_size(32, kernel, stride, padding),
                 'note': note})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('The second row is the one to default to: kernel 3, stride 1,')
print('padding 1 keeps the spatial size exactly, so you can stack')
print('as many as you like without doing arithmetic. Downsampling')
print('is then a deliberate separate step, which is how every')
print('modern architecture is organised.')

import pandas as pd

rows = []
for size, channels_in, channels_out in [(32, 3, 16), (64, 16, 32),
                                        (224, 3, 64), (112, 64, 128)]:
    pixels_in = size * size * channels_in
    pixels_out = size * size * channels_out

    dense = pixels_in * pixels_out + pixels_out
    conv = (3 * 3 * channels_in * channels_out) + channels_out

    rows.append({'input': f'{channels_in}x{size}x{size}',
                 'output channels': channels_out,
                 'fully connected': f'{dense:,}',
                 '3x3 convolution': f'{conv:,}',
                 'ratio': f'{dense / conv:,.0f}x'})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('A fully-connected layer mapping one 224x224 colour image to')
print('64 channels at the same resolution needs 483 BILLION')
print('weights - about two terabytes in float32, for one layer.')
print('The convolution needs 1,792, which is 7 kilobytes. That is')
print('not an optimisation, it is the difference between a model')
print('that exists and one that does not.')
print()
print('And the parameter count is independent of the image size:')
print('the same 1,792 weights work on 32x32 and on 4000x3000.')

import torch
import torch.nn as nn

torch.manual_seed(0)

layer = nn.Conv2d(in_channels=3, out_channels=16, kernel_size=3, padding=1)
batch = torch.randn(8, 3, 32, 32)
output = layer(batch)

print('input shape ', tuple(batch.shape), ' (batch, channels, height, width)')
print('output shape', tuple(output.shape))
print()
weight, bias = layer.weight, layer.bias
print('weight shape', tuple(weight.shape),
      '= (out channels, in channels, kernel, kernel)')
print('bias shape  ', tuple(bias.shape))
print('parameters  ', weight.numel() + bias.numel(),
      '= 16 x 3 x 3 x 3 + 16')
print()
print('Each of the 16 output channels has its own 3x3x3 kernel: it')
print('looks at a 3x3 neighbourhood across ALL THREE input')
print('channels at once. "16 channels" means sixteen different')
print('learned feature detectors, each producing its own 32x32 map')
print('of where it fired.')

import pandas as pd

# Receptive field of a stack: r = r + (k - 1) * jump, jump = jump * stride
rows = []
field, jump, size = 1, 1, 224
for index, (kernel, stride, name) in enumerate([
        (7, 2, 'conv 7x7 /2'), (3, 2, 'maxpool 3x3 /2'),
        (3, 1, 'conv 3x3'), (3, 1, 'conv 3x3'),
        (3, 2, 'conv 3x3 /2'), (3, 1, 'conv 3x3'),
        (3, 2, 'conv 3x3 /2'), (3, 1, 'conv 3x3'),
        (3, 2, 'conv 3x3 /2'), (3, 1, 'conv 3x3')], 1):
    field = field + (kernel - 1) * jump
    jump = jump * stride
    size = (size + 2 * (kernel // 2) - kernel) // stride + 1
    rows.append({'layer': f'{index}. {name}',
                 'feature map': f'{size}x{size}',
                 'receptive field': f'{field}x{field}',
                 'covers': f'{min(100.0, 100 * field / 224):.0f}% of a 224px image'})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('Layer 1 sees a 7x7 patch - 3% of the image, enough for an')
print('edge and nothing more. Layer 10 sees 195x195, 87% of it,')
print('which is enough for "a dog". Depth is how the receptive')
print('field grows; stride is what makes it grow fast. Notice the')
print('jumps: crossing a stride-2 layer doubles how much each')
print('later 3x3 covers, so the field grows by 32 pixels a layer')
print('near the end against 8 near the start.')

import numpy as np

feature_map = np.array([
    [1, 3, 2, 1],
    [4, 9, 1, 0],
    [2, 1, 5, 7],
    [1, 0, 2, 3],
], dtype=float)


def max_pool(x, size=2):
    h, w = x.shape[0] // size, x.shape[1] // size
    return np.array([[x[r * size:(r + 1) * size,
                        c * size:(c + 1) * size].max()
                      for c in range(w)] for r in range(h)])


def average_pool(x, size=2):
    h, w = x.shape[0] // size, x.shape[1] // size
    return np.array([[x[r * size:(r + 1) * size,
                        c * size:(c + 1) * size].mean()
                      for c in range(w)] for r in range(h)])


print('feature map')
print(feature_map.astype(int))
print()
print('max pool 2x2')
print(max_pool(feature_map).astype(int))
print()
print('average pool 2x2')
print(average_pool(feature_map).round(2))
print()

# The invariance claim, tested rather than asserted.
shifted = np.roll(feature_map, 1, axis=1)
print('after shifting the map one pixel right:')
print(shifted.astype(int))
print('max pool  ', max_pool(shifted).astype(int).ravel(),
      ' (before:', max_pool(feature_map).astype(int).ravel(), ')')
print()
print('Note what that shows: pooling is only invariant to shifts')
print('INSIDE a window. A one-pixel shift that moves a value')
print('across a window boundary changes the output. Pooling buys')
print('a little tolerance, not translation invariance, and the')
print('difference is why heavily-pooled networks still need')
print('augmentation.')

import torch
import torch.nn as nn

torch.manual_seed(0)
body = nn.Sequential(nn.Conv2d(3, 32, 3, stride=2, padding=1), nn.ReLU(),
                     nn.Conv2d(32, 64, 3, stride=2, padding=1), nn.ReLU())
head = nn.Sequential(nn.AdaptiveAvgPool2d(1), nn.Flatten(), nn.Linear(64, 10))

for size in [32, 64, 224]:
    image = torch.randn(1, 3, size, size)
    features = body(image)
    logits = head(features)
    shape = f'{size}x{size}'
    print(f'input {shape:<9} -> features {str(tuple(features.shape)[1:]):<15}'
          f' -> logits {tuple(logits.shape)}')
print()
print('Three different input sizes, one set of weights, the same')
print('ten logits out. A Flatten followed by a Linear would have')
print('failed on two of the three, because its input size is')
print('baked into the weight matrix. This is why global average')
print('pooling replaced it.')

import numpy as np
import torch
import torch.nn as nn
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

torch.manual_seed(0)
torch.set_num_threads(4)

digits = load_digits()
X = digits.images[:, None, :, :] / 16.0           # (n, 1, 8, 8)
y = digits.target
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3,
                                          random_state=0, stratify=y)
X_tr = torch.tensor(X_tr, dtype=torch.float32)
X_te = torch.tensor(X_te, dtype=torch.float32)
y_tr = torch.tensor(y_tr, dtype=torch.long)
y_te = torch.tensor(y_te, dtype=torch.long)
print(f'{len(X_tr)} train, {len(X_te)} test, images {tuple(X_tr.shape[1:])}')


def dense_model():
    return nn.Sequential(nn.Flatten(), nn.Linear(64, 64), nn.ReLU(),
                         nn.Linear(64, 10))


def conv_model():
    return nn.Sequential(
        nn.Conv2d(1, 16, 3, padding=1), nn.ReLU(),
        nn.Conv2d(16, 32, 3, padding=1), nn.ReLU(),
        nn.MaxPool2d(2),                            # 8x8 -> 4x4
        nn.AdaptiveAvgPool2d(1), nn.Flatten(),
        nn.Linear(32, 10))


def train(build, epochs=40, rate=0.01, seed=0):
    torch.manual_seed(seed)
    model = build()
    optimiser = torch.optim.AdamW(model.parameters(), lr=rate)
    loss_fn = nn.CrossEntropyLoss()
    for _ in range(epochs):
        model.train()
        order = torch.randperm(len(X_tr))
        for start in range(0, len(X_tr), 64):
            index = order[start:start + 64]
            optimiser.zero_grad()
            loss = loss_fn(model(X_tr[index]), y_tr[index])
            loss.backward()
            optimiser.step()
    model.eval()
    with torch.no_grad():
        train_acc = (model(X_tr).argmax(1) == y_tr).float().mean().item()
        test_acc = (model(X_te).argmax(1) == y_te).float().mean().item()
    parameters = sum(p.numel() for p in model.parameters())
    return parameters, train_acc, test_acc


print()
print('model           parameters  train acc  test acc')
for name, build in [('fully connected', dense_model), ('convolutional', conv_model)]:
    parameters, train_acc, test_acc = train(build)
    print(f'{name:<15} {parameters:10,}  {train_acc:9.4f}  {test_acc:8.4f}')
print()
print('A dead heat: 0.9722 each. On 8x8 digits that are already')
print('centred and cropped, the dense model has nothing to gain -')
print('every digit is in the same place, so weight sharing buys no')
print('information. The advantage is invisible until the')
print('assumption behind it is tested.')
print()

# The claim everyone makes about convolutions, tested: are they
# more robust to a shifted input?
def shift_right(images, pixels=1):
    out = torch.zeros_like(images)
    out[:, :, :, pixels:] = images[:, :, :, :-pixels]
    return out


print('and on test images shifted one pixel to the right:')
for name, build in [('fully connected', dense_model), ('convolutional', conv_model)]:
    torch.manual_seed(0)
    model = build()
    optimiser = torch.optim.AdamW(model.parameters(), lr=0.01)
    loss_fn = nn.CrossEntropyLoss()
    for _ in range(40):
        model.train()
        order = torch.randperm(len(X_tr))
        for start in range(0, len(X_tr), 64):
            index = order[start:start + 64]
            optimiser.zero_grad()
            loss_fn(model(X_tr[index]), y_tr[index]).backward()
            optimiser.step()
    model.eval()
    with torch.no_grad():
        clean = (model(X_te).argmax(1) == y_te).float().mean().item()
        shifted = (model(shift_right(X_te)).argmax(1) == y_te).float().mean().item()
    print(f'{name:<15} clean {clean:.4f}  shifted {shifted:.4f}  '
          f'lost {clean - shifted:.4f}')
print()
print('Both models lose accuracy on a shift that no human would')
print('even notice, and the convolutional one loses less. "Shift')
print('invariant" is an overstatement of what a convolution gives')
print('you - it shares weights, so a shifted feature is still')
print('detected, but pooling boundaries and the classifier head')
print('still see a different input.')

import pandas as pd

# VGG's claim: two stacked 3x3 layers have the receptive field of
# one 5x5 and fewer parameters.
channels = 64
one_5x5 = 5 * 5 * channels * channels
two_3x3 = 2 * (3 * 3 * channels * channels)
print(f'one 5x5 layer, 64->64 channels : {one_5x5:,} weights')
print(f'two 3x3 layers, 64->64 channels: {two_3x3:,} weights')
print(f'both see a 5x5 region; the pair is {1 - two_3x3 / one_5x5:.0%} smaller')
print('and has a non-linearity in the middle, which the single')
print('5x5 layer does not.')
print()

# MobileNet's claim: split a convolution into a per-channel spatial
# part and a 1x1 channel-mixing part.
rows = []
for cin, cout in [(32, 64), (64, 128), (256, 512)]:
    standard = 3 * 3 * cin * cout
    separable = 3 * 3 * cin + cin * cout
    rows.append({'channels': f'{cin} -> {cout}',
                 'standard 3x3': f'{standard:,}',
                 'depthwise separable': f'{separable:,}',
                 'reduction': f'{standard / separable:.1f}x'})
print(pd.DataFrame(rows).to_string(index=False))
print()
print('Roughly an 8-9x reduction, and it grows with the channel')
print('count. This single substitution is what put convolutional')
print('networks on phones.')

import pandas as pd

rows = []
for name, size, channels in [('conv1', 112, 64), ('conv2', 56, 128),
                             ('conv3', 28, 256), ('conv4', 14, 512),
                             ('conv5', 7, 512)]:
    activations = size * size * channels
    weights = 3 * 3 * channels * channels
    rows.append({'layer': name, 'map': f'{channels}x{size}x{size}',
                 'activations per image': f'{activations:,}',
                 'weights in the layer': f'{weights:,}'})
table = pd.DataFrame(rows)
print(table.to_string(index=False))
print()
print('conv1 holds 802,816 activations per image and 36,864')
print('weights. conv5 holds 25,088 activations and 2.4 million')
print('weights. They are inverses of each other, so "reduce the')
print('model size" and "stop running out of memory" are two')
print('different jobs on two different layers.')

import numpy as np
import torch
import torch.nn as nn

torch.manual_seed(0)

image = np.random.default_rng(0).normal(0, 1, (1, 1, 6, 6)).astype(np.float32)
layer = nn.Conv2d(1, 1, 3, padding=0, bias=False)
kernel = layer.weight.detach().numpy()[0, 0]

framework = layer(torch.tensor(image)).detach().numpy()[0, 0]

by_hand = np.zeros((4, 4), dtype=np.float32)
for row in range(4):
    for column in range(4):
        patch = image[0, 0, row:row + 3, column:column + 3]
        by_hand[row, column] = (patch * kernel).sum()

print('framework')
print(framework.round(4))
print()
print('by hand')
print(by_hand.round(4))
print()
print('max difference', float(np.abs(framework - by_hand).max()))
print('agree:', bool(np.allclose(framework, by_hand, atol=1e-5)))
print()
print('They agree, which tells you something that trips up anyone')
print('coming from a signal-processing background: what deep')
print('learning calls "convolution" is really cross-correlation.')
print('A true convolution flips the kernel first. Since the kernel')
print('is learned, the flip makes no difference to what the')
print('network can represent - but it matters the moment you')
print('compare against scipy.signal.convolve2d, which does flip.')