Convolutions & Computer Vision: Spatial Feature Hierarchies
Convolutions & Computer Vision: Spatial Feature Hierarchies
In Chapter 4: Neural Building Blocks, we explored fully connected linear layers.
What happens if you try to feed a standard high-resolution image into a fully connected layer? - A modest color photograph has resolution $224 \times 224 \times 3 \approx 150,528$ values. - Connecting this image to a single hidden layer of 1,024 neurons requires:
$$150,528 \times 1,024 \approx 154,140,672 \text{ parameters}$$
Over 154 million parameters for just one layer! The model would immediately overfit, consume gigabytes of VRAM, and completely ignore the spatial geometry of the image.
This chapter introduces the mathematical principles of 2D convolutions, the im2col matrix unrolling transformation, Winograd minimal filtering algebra, and how Deep Residual Networks (ResNet) revolutionized computer vision.
1. Core Intuition: A Sliding Magnifying Glass
Imagine inspecting a painting with a small $3\times3$ magnifying lens designed to detect sharp edges:
$$\mathbf{X}_{\text{patch}} = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 1 \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{K}_{\text{vertical}} = \begin{bmatrix} -1 & 0 & +1 \\ -1 & 0 & +1 \\ -1 & 0 & +1 \end{bmatrix} \implies (\mathbf{X} * \mathbf{K}) = \sum_{i,j} X_{ij} K_{ij} = +3$$
You slide this magnifying glass across every position in the image: 1. When the glass hovers over a flat region (all zeros), the output is zero. 2. When the glass crosses the boundary between dark and bright pixels, the product spikes to a high positive number.
Key Insight: Weight Sharing: Instead of learning a separate weight for every pixel coordinate, we learn one small $3\times3$ filter (9 weights) and slide it across every pixel. A feature useful in the top-left corner is equally useful in the bottom-right corner (Translation Invariance).
2. Mathematical Formalization of 2D Convolutions
Given an input feature map $\mathbf{X} \in \mathbb{R}^{C_{in} \times H \times W}$ and a collection of $C_{out}$ convolution filters $\mathbf{K} \in \mathbb{R}^{C_{out} \times C_{in} \times K_h \times K_w}$, the activation at output coordinate $(c_{out}, i, j)$ is:
$$\mathbf{Y}_{c_{out}, i, j} = \mathbf{b}_{c_{out}} + \sum_{c_{in}=1}^{C_{in}} \sum_{m=-k}^{k} \sum_{n=-k}^{k} \mathbf{X}_{c_{in}, \, i \cdot s + m, \, j \cdot s + n} \cdot \mathbf{K}_{c_{out}, c_{in}, m, n}$$
Where $s$ is the stride and padding $p$ determines boundary handling.
The Spatial Feature Hierarchy
As you stack convolutional layers: - Early Layers: Learn primitive local features (horizontal edges, diagonal gradients, color contrasts). - Middle Layers: Combine edges into textures and basic geometric motifs (corners, circles, honeycombs). - Deep Layers: Combine motifs into high-level semantic object parts (eyes, wheels, animal ears).
3. GPU Hardware Execution: The im2col Transformation
A naive C implementation of 2D convolution requires seven nested loops:
for (int b = 0; b < batch; b++)
for (int co = 0; co < c_out; co++)
for (int h = 0; h < out_h; h++)
for (int w = 0; w < out_w; w++)
for (int ci = 0; ci < c_in; ci++)
for (int kh = 0; kh < k_h; kh++)
for (int kw = 0; kw < k_w; kw++)
// compute...
Executing seven nested loops on a GPU results in abysmal memory bandwidth because threads read overlapping memory addresses incoherently.
The Solution: Matrix Multiplication via im2col
Chellapilla et al. (2006) introduced im2col (Image to Column):
(C_in · K_h · K_w) × (H_out · W_out)"] Filter["Filter Matrix:
(C_out) × (C_in · K_h · K_w)"] --> GEMM["High-Throughput GEMM:
Filters · Unrolled Columns"] Unroll --> GEMM GEMM --> Mat["Output Matrix:
(C_out) × (H_out · W_out)"] Mat -->|"Reshape"| Out["Output Tensor: (C_out, H_out, W_out)"]
By unrolling local $3\times3$ spatial patches into contiguous matrix columns, the entire convolution is transformed into a single, highly-optimized GEMM matrix multiplication.
4. Winograd Minimal Filtering: $2.25\times$ Arithmetic Reduction
For standard $3\times3$ convolutions, Lavin & Winograd (2016) demonstrated that computing outputs in the spatial domain performs redundant arithmetic.
Using polynomial residue number systems, Winograd $F(2 \times 2, 3 \times 3)$ computes a $2\times 2$ output block from a $4\times 4$ input tile and a $3\times 3$ filter using only 16 multiplications instead of 36:
$$\mathbf{Y} = \mathbf{A}^T \left[ (\mathbf{G} \mathbf{K} \mathbf{G}^T) \odot (\mathbf{B}^T \mathbf{X} \mathbf{B}) \right] \mathbf{A}$$
Where $\mathbf{B}^T, \mathbf{G}, \mathbf{A}^T$ are fixed transformation matrices consisting solely of small integer powers of 2 (computed via fast additions and bit shifts in GPU registers):
$$\mathbf{B}^T = \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 1 & 1 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 1 & 0 & -1 \end{bmatrix}, \quad \mathbf{G} = \begin{bmatrix} 1 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & \frac{1}{2} \\ \frac{1}{2} & -\frac{1}{2} & \frac{1}{2} \\ 0 & 0 & 1 \end{bmatrix}, \quad \mathbf{A}^T = \begin{bmatrix} 1 & 1 & 1 & 0 \\ 0 & 1 & -1 & -1 \end{bmatrix}$$
This yields an immediate $2.25\times$ reduction in multiplication FLOPs directly inside OpenCL compute workgroups.
5. Deep Residual Networks (ResNet): Solving Gradient Vanishing
In 2015, researchers observed a paradox: stacking more convolutional layers caused training accuracy to degrade rapidly, not because of overfitting, but because gradients could not flow back through 50 or 100 consecutive weight matrices.
He et al. (2015) introduced the Residual Identity Shortcut:
The Gradient Highway
Applying the chain rule to the residual block $\mathbf{y} = \mathcal{F}(\mathbf{x}) + \mathbf{x}$ reveals why this architecture changed deep learning forever:
$$\frac{\partial \mathcal{L}}{\partial \mathbf{x}} = \frac{\partial \mathcal{L}}{\partial \mathbf{y}} \cdot \left( \frac{\partial \mathcal{F}(\mathbf{x})}{\partial \mathbf{x}} + \mathbf{I} \right) = \frac{\partial \mathcal{L}}{\partial \mathbf{y}} \cdot \frac{\partial \mathcal{F}(\mathbf{x})}{\partial \mathbf{x}} + \mathbf{\frac{\partial \mathcal{L}}{\partial \mathbf{y}}}$$
Notice the term $+\frac{\partial \mathcal{L}}{\partial \mathbf{y}}$! Even if the weights in $\mathcal{F}(\mathbf{x})$ have near-zero gradients, the gradient flows unattenuated and completely intact directly back to the earliest layers through the identity shortcut.
6. Prototypical NetCL Implementation
NetCL provides hardware-accelerated 2D convolutions and complete ResNet models in netcl.nn:
import numpy as np
import netcl.autograd as ag
import netcl.nn as nn
from netcl.core.device import manager
from netcl.core.tensor import Tensor
q = manager.default("auto").queue
# 1. Standard 2D Convolutional Layer
# in_channels=3, out_channels=32, kernel_size=3, padding=1
conv = nn.Conv2d(q, in_channels=3, out_channels=32, kernel_size=3, padding=1)
# 2. Input image batch: (Batch=4, Channels=3, Height=32, Width=32)
images = Tensor.from_host(q, np.random.randn(4, 3, 32, 32).astype(np.float32))
with ag.Tape() as tape:
x_node = ag.tensor(images)
features = conv(x_node)
print(f"Output Feature Map Shape: {features.value.shape}")
assert features.value.shape == (4, 32, 32, 32)
# 3. Instantiate a complete ResNet-18 architecture in NetCL
from netcl.nn.resnet import ResNet18
resnet = ResNet18(q, num_classes=10)
print(f"ResNet-18 created with {len(list(resnet.parameters()))} parameter tensors.")
Related Documentation
- Concepts: Winograd Convolution: Mathematical derivation of Winograd transformation matrices and kernel implementation.
- Concepts: im2col: Memory layout and strided unrolling for GPU GEMM.
- Concepts: ResNet: Residual block variants and bottleneck architectures.
Next Steps in the Curriculum
Now that you understand 2D spatial feature extraction, how do we keep internal activations numerically stable across deep layers during high-throughput training?
Proceed to Chapter 6: Normalization & Regularization to explore Batch Normalization, Layer Normalization, and Dropout.