Implementing ConvNets from Scratch in C++

01 // THE LANDSCAPE IN 2013

It is September 2013. AlexNet is still the paper everyone is talking about. It won ImageNet by a large margin, and convolutional neural networks suddenly looked practical rather than academic. I wanted to understand the mechanics: what the layers computed, how gradients flowed, and where the runtime went.

Most vision code I saw was still SIFT, HOG, Haar cascades, or some hand-written OpenCV pipeline. Caffe was early. Theano existed, but it felt far away from normal C++ code. Torch7 was Lua. TensorFlow and PyTorch were not options yet.

So I wrote a small one myself. Raw C++. No autograd. No layer framework. Just arrays, loops, a little matrix code, and the backprop equations from LeCun's papers.

"The useful work was not just getting MNIST accuracy. It was learning how to trust gradients by checking them numerically."

02 // ARCHITECTURE

The architecture I chose is essentially LeNet-5, Yann LeCun's 1998 design for handwritten digit recognition. It's small enough to train on a single CPU core in reasonable time, but deep enough to demonstrate the core ideas: local receptive fields, weight sharing, spatial downsampling, and learned hierarchical features. The full pipeline is Input → Conv → ReLU → Pool → Conv → ReLU → Pool → FC → FC → Softmax.

The first convolutional layer applies 6 filters of size 5×5 to the 28×28 grayscale input, producing 6 feature maps of size 24×24. Max-pooling with a 2×2 window and stride 2 halves each spatial dimension to 12×12. The second conv layer applies 16 filters of size 5×5 across all 6 input channels, yielding 16 feature maps of size 8×8, which pool down to 4×4. The flattened 256-dimensional vector feeds through two fully connected layers (120 and 84 units) before a 10-way softmax output for digit classification.

In C++, each layer is a struct holding its weights, biases, and cached activations. I kept the implementation explicit: one forward and backward function per layer type. That made it easier to compare analytic gradients against finite differences and isolate shape bugs.

03 // FORWARD PASS

The convolution operation is conceptually a sliding dot product. For a single input channel and a single filter, the output feature map is computed by sliding the kernel across the input and computing the inner product at each position:

For multi-channel inputs (e.g., the second conv layer which takes 6 input feature maps), the summation extends over the channel dimension as well. Each of the 16 output filters has a 5×5×6 weight tensor — 150 parameters per filter — and the dot products are accumulated across all input channels before adding the bias.

The naïve implementation is a six-deep nested loop: output filter, output row, output column, input channel, kernel row, kernel column. On a 2013-era Intel i5, this takes about 12ms per image for the first conv layer alone. The standard optimization is the im2col trick: rearrange the input patches into columns of a matrix, then express the entire convolution as a single matrix multiplication.

Here F is the number of filters, C is the number of input channels, K is the kernel size, and HW is the number of output spatial positions. The im2col matrix has redundant data — each input pixel appears in up to K² columns — but the tradeoff is worth it because we can now call an optimized BLAS sgemm routine instead of hand-written loops. On the CPU, this alone gave a 4× speedup.

ReLU is trivially max(0, x)applied elementwise. Max-pooling takes the maximum value in each 2×2 window, also recording the index of the maximum for use during backpropagation (the "switches").

04 // BACKPROPAGATION

The loss function is cross-entropy between the softmax output and the one-hot target label. For a single training example with true class c, the loss is:

The gradient of the loss with respect to the pre-softmax logits is ∂L/∂z_k = ŷ_k − y_k, where y is the one-hot target. That compact result becomes the seed for the rest of backpropagation. Each downstream layer is a local Jacobian-vector product.

For a fully connected layer with y = W·x + b, the gradients are straightforward:

where δ = ∂L/∂y is the upstream gradient. The convolutional layer is trickier. The gradient with respect to the weights is itself a convolution — you correlate the input with the upstream gradient. The gradient with respect to the inputrequires a "full convolution" — convolving the upstream gradient with the 180°-rotated filter, padded to reconstruct the input spatial dimensions:

In practice, both operations are implemented via im2col and matrix multiplication, just with different rearrangements. The max-pool backward pass is sparse: gradients flow only through the positions that were selected as maxima during the forward pass, zeroing out all other positions. ReLU's backward pass multiplies the upstream gradient by 1 where the input was positive and 0 where it was negative.

Getting this right took careful testing. The most common bug was transposing a matrix that shouldn't have been transposed, or vice versa. I wrote a numerical gradient checker: perturb each weight by ε = 1e-5 and compare the finite-difference gradient to the analytic one. When the relative error stayed below 1e-7 across parameters, I trusted the implementation enough to train.

05 // SGD & MOMENTUM

With correct gradients in hand, the simplest optimizer is vanilla stochastic gradient descent. For each mini-batch of B training examples, we compute the average gradient and update:

Learning rate η is the single most important hyperparameter. Too large and the loss diverges; too small and training stalls. In 2013, the standard approach was to start with η = 0.01 and manually decay it by 10× when the validation loss plateaued — no cosine annealing, no warmup schedules, no learning rate finders. Just stare at the loss curve and adjust.

Momentum improved convergence and made the loss curve less noisy. Instead of stepping directly in the gradient direction, I kept a velocity vector that accumulated recent gradients:

Weight initialization was an unsolved problem in 2013. Xavier initialization had been published in 2010, but Kaiming initialization (which accounts for ReLU's zero-mean asymmetry) wouldn't arrive until He et al. 2015. I used the common heuristic of the era: sample each weight from N(0, 0.01)and initialize biases to zero. This worked for five layers but would catastrophically fail for deeper networks — a problem I didn't fully appreciate until later.

06 // CUDA ACCELERATION

The CPU implementation trained on MNIST at about 45 seconds per epoch with batch size 32. That was fine for learning, but slow enough to be annoying. CUDA 5.5 was available, and the GTX 680 had enough cores to make the experiment worth trying.

The key insight is that convolution via im2col reduces to a General Matrix Multiply (GEMM), and NVIDIA provides a heavily optimized GEMM implementation in cuBLAS. Calling cublasSgemm instead of my hand-written CPU loops gave an immediate 30× speedup on the convolution layers. The remaining bottleneck was memory transfer — copying input batches from host RAM to device GDDR5.

Writing custom CUDA kernels for ReLU and max-pooling taught me the fundamentals of GPU programming: thread blocks, warps, shared memory tiling, and memory coalescing. The most important lesson was that GPUs are not just "fast CPUs" — they require fundamentally different data layout. Strided access patterns that cost nothing on a CPU (because the cache line handles it) can destroy GPU throughput by causing uncoalesced global memory reads.

With cuBLAS for GEMM and small custom kernels for the rest, MNIST training dropped to under 1 second per epoch. That made iteration faster enough to run real experiments: learning-rate sweeps, initialization checks, and ablations on pooling and hidden-layer size.

07 // RESULTS & MNIST

After four weeks of implementation and debugging, the network was training and converging. On the MNIST test set (10,000 images), the final results told a clear story:

The 1.3% gap between the ConvNet and SVM+HOG might seem small, but on MNIST it represents roughly 130 fewer misclassifications out of 10,000 test images. More importantly, the ConvNet required zero feature engineering. The HOG pipeline needed careful tuning of cell size, block normalization, and histogram bins. The ConvNet learned everything from pixels.

The most fascinating result was visualizing the learned first-layer filters. After training, the 6 filters in Conv1 had spontaneously organized into edge and corner detectors — horizontal edges, vertical edges, diagonal edges, and simple curvature detectors. These are qualitatively identical to the Gabor-like filters that neuroscientists had observed in the primary visual cortex (V1) of cats and primates decades earlier. The network had independently rediscovered the same feature representation that biological visual systems evolved over millions of years.

"When your randomly-initialized 5×5 filters converge to edge detectors without any supervision beyond digit labels, you know the representation learning hypothesis is real. The features are not prescribed — they are discovered."

08 // LOOKING BACK

Writing this from scratch was useful because it removed the mystery. After implementing forward, backward, update, gradient checking, and a few CUDA kernels, I understood where the numbers came from and how to debug a failing training run.

The frameworks that followed solved all of these problems. Caffe gave us declarative layer definitions in protobuf. Torch7 brought dynamic computation in Lua. Theano compiled symbolic graphs to optimized CUDA. Then TensorFlow (2015) industrialized the computation graph, and PyTorch (2016) made dynamic graphs the default. Today, you can train a ResNet-50 on ImageNet in under an hour with eight lines of PyTorch. The entire forward-backward-update loop that took me weeks to implement is a single call to loss.backward().

But I believe the raw math still matters. When your model doesn't converge, the debugging process is the same as it was in 2013: check your gradients numerically, visualize your activations, inspect your weight distributions. The tools are better, but the mental model is identical. Understanding that a convolution is a structured matrix multiplication, that backpropagation through a conv layer is a full convolution with a rotated kernel, that momentum is an exponential moving average of gradients — these are the invariants that survive every framework migration.

Frameworks and hardware will keep changing. The useful part of this exercise is simpler: once you have implemented the chain rule by hand, the abstractions are easier to trust and easier to debug.