Efficiently evaluating deep neural networks
A trained neural network takes an input, passes it through a sequence of layers, and produces an output. Many of those layers repeat the same multiply-and-add operation millions of times. That sounds ideal for parallel hardware, yet a fast chip can still spend much of its time moving data instead of calculating. This lecture follows the computation from a single neuron to a whole network and asks, at each step, how to do less work or make better use of the hardware.
- You'll be able to
- Explain what a fully connected layer, a convolution layer, ReLU and pooling compute, and count their cost. Explain why depthwise separable convolutions (MobileNet) are cheaper. Map a convolution onto a matrix multiplication, explicitly (im2col) or implicitly, and explain the cost of each. Explain how blocking raises arithmetic intensity. Explain how operator fusion saves memory traffic, including the chunked-softmax trick behind fused attention. Explain why tensor cores and DNN accelerators exist.
- Background
- Lecture 3 (arithmetic intensity, bandwidth-bound execution) and Lecture 7 (GPUs, CUDA, shared memory). Unit 5 of these notes covers the roofline model, tiling and number formats in more depth.
- Source
- Follows the structure and examples of Stanford CS149 Parallel Computing, Fall 2023, Lecture 10, "Efficiently Evaluating DNNs" (Kayvon Fatahalian and Kunle Olukotun). Lectures 8–9 aren't covered here yet. The explanations, diagrams and explorers are original study material, not a copy of the slides.
The efficiency challenge
A network's topology is the choice, order, and shape of its layers. VGG, ResNet, Inception, and MobileNet make different choices, so they demand different amounts and kinds of work. The same network may also run on a server GPU, a phone CPU, or a dedicated neural processor. An implementation that is quick on one device may be slow on another. We need to understand both the work the model asks for and the machine that will perform it.
The lecture groups the answers into three strategies, which structure this page:
- Ask for less work. Choose a network topology that reaches the needed accuracy with fewer operations (Section 4).
- Use data well. Arrange each layer so values are reused near the arithmetic units, and combine neighboring operations so intermediate results need fewer trips to memory (Parts 2 and 3).
- Make each operation cheaper. Use an appropriate numeric precision and hardware that performs many multiply-adds with little instruction overhead (Parts 3 and 4).
Here evaluation, or inference, means applying a network whose weights have already been learned to a new input. We will follow the forward calculation only. Training also performs this forward pass, then computes gradients and updates the weights; it brings additional work and memory needs.
Neurons, layers and matrix-vector products
Start with the lecture's small expression, max(max(0, a·b + c·d) + e·f + g·h, i·j). To evaluate it, multiply pairs of numbers, add some results, and take two maxima. A neural network is a much larger graph of operations of the same general kind. The graph matters because it tells us which results can be computed independently and which must wait for earlier results.
One familiar building block is a neuron. It multiplies each input xi by a weight wi, adds the products and a bias b, then applies a function f:
out = f( Σi xi·wi + b ) // f is often ReLU: f(v) = max(0, v)
The weights and bias are the parameters found during training; they stay fixed during this evaluation. The function f makes the layer nonlinear, allowing a stack of layers to represent more than one large weighted sum. A common choice is ReLU (rectified linear unit): negative inputs become zero, and positive inputs pass through.
A fully connected layer has several neurons, each reading every input value. Put one neuron's weights in each row of a matrix. Multiplying that matrix by the input vector computes all the weighted sums at once; then add the biases and apply f to each result:
The number of weights is inputs × outputs. A 224 × 224 RGB image contains 150,528 input values; connecting all of them to just 1,000 neurons would require 150,528,000 weights, before biases. Storing and using that many weights in an early image layer is costly. The next layer type keeps the useful local structure of an image without connecting every pixel to every output.
Convolution: local connections with shared weights
The first layer need not look at the whole image to recognize a small feature. Two choices make it much cheaper:
- Local connections. An edge or corner can be detected from nearby pixels. Each output can therefore read a small input window, such as 3 × 3, instead of the entire image.
- Shared weights. The same kind of edge might occur anywhere. Reuse one set of weights at every location. Moving the input feature then moves its response in the output; this is translation equivariance. The shared weights are a filter, or kernel, and applying it across the image is called a convolution layer.
Here is the lecture's 3 × 3 example on a single-channel image. The input array includes a one-pixel border whose values have been chosen beforehand, so every 3 × 3 window is in bounds. At each output position, the code multiplies nine pixels by the same nine weights and adds the products. With every weight equal to 1/9, the result is an average, or box blur:
int WIDTH = 1024, HEIGHT = 1024; float input[(WIDTH+2) * (HEIGHT+2)]; // initialized image plus a 1-pixel border float output[WIDTH * HEIGHT]; float weights[] = {1/9., 1/9., 1/9., 1/9., 1/9., 1/9., 1/9., 1/9., 1/9.}; for (int j=0; j<HEIGHT; j++) for (int i=0; i<WIDTH; i++) { float tmp = 0.f; for (int jj=0; jj<3; jj++) for (int ii=0; ii<3; ii++) tmp += input[(j+jj)*(WIDTH+2) + (i+ii)] * weights[jj*3 + ii]; output[j*WIDTH + i] = tmp; }
The filter changes the question asked of each window. Sobel filters, for example, respond to changes in brightness across a horizontal or vertical direction. They highlight edges instead of blurring them. Try the choices below. In machine-learning libraries, this sliding dot product is usually called “convolution” even though the filter is not reversed as it would be in the mathematical definition.
Convolution layers in a real network
An RGB input has three channels; deeper layers may have dozens or hundreds. To produce one output channel, one filter looks across all input channels: on a 64-channel input, a 3 × 3 filter has 3 × 3 × 64 weights. Apply K different filters and the output has K channels. AlexNet's first layer, for example, uses 96 filters of size 11 × 11 × 3.
One familiar network pattern applies ReLU after convolution, then pooling. A 2 × 2 max pool keeps the largest value in each nonoverlapping 2 × 2 region of a channel. With stride 2 and even dimensions, it halves the width and height. The diagram assumes the convolution uses padding that preserves the spatial dimensions; other padding and stride choices give different sizes.
We can now count the convolution's work. Let its output have H × W positions and K channels; let the input have C channels and each filter have R × S spatial weights. Each output value sums C·R·S products:
multiply-adds = H · W · K · C · R · S // every output pixel, every filter, every weight parameters = K · C · R · S (+ K biases) // weights reused at every position
The output dimensions already reflect the chosen stride and padding. Notice the difference between the two counts: a larger image increases the number of multiply-adds, but it does not create more filter weights. That repeated use of the same weights is a source of reuse for a fast implementation.
Approach 1: design topologies that need less work
The first optimization happens before implementation: choose a network that asks for less work. The lecture compares several ImageNet classifiers from 2014–2017. Top-1 accuracy is the fraction of images for which the model's highest-scoring label is correct; MADDs counts multiply-adds for one image. These published figures are useful for scale, though evaluation details differ among models:
| Network | Year | Top-1 accuracy | Parameters | MADDs per image |
|---|---|---|---|---|
| VGG-16 | 2014 | 71.5% | 138 M | 15 B |
| GoogLeNet | 2015 | 70% | 6.8 M | 1.5 B |
| ResNet-18 | 2016 | 73%* | 11.7 M | 1.8 B |
| MobileNet-224 | 2017 | 70.5% | 4.2 M | 0.6 B |
* 10-crop result; ResNet's single-crop results are similar to the other networks'.
VGG-16 and MobileNet-224 land in roughly the same accuracy range in this table, while their listed multiply-add counts differ by about 25×: 15 B / 0.6 B. That is a change in the work requested, before any faster chip or kernel is considered. The accuracy figures are not a controlled head-to-head comparison; for example, the ResNet-18 entry uses ten crops per image. ResNet's residual connections provide a route around layers during training, making very deep models easier to optimize. They are one example of how network structure changes both accuracy and computational cost.
MobileNet's trick: depthwise separable convolution
MobileNet's key change is easier to see with the count from Section 3. At every output position, a standard R × S convolution uses R·S·C·K multiply-adds: each of K filters spans all C input channels. Depthwise separable convolution breaks the spatial and channel work into two stages:
- Depthwise: C separate 3 × 3 × 1 filters, one per input channel. Spatial filtering, no channel mixing.
- Pointwise: K filters of size 1 × 1 × C. Channel mixing, no spatial extent.
Per output position, the two stages cost R·S·C + C·K multiply-adds. Divide by the standard cost R·S·C·K and you get the ratio in the caption: 1/K + 1/(R·S). For a 3 × 3 filter and a large K, the cost approaches one ninth of the original. MobileNet's paper reports that about 95% of its multiply-adds are in the remaining 1 × 1 convolutions. A 1 × 1 convolution can be viewed as a matrix multiplication over pixels and channels, so the next part can use decades of work on fast matrix multiplication.
The direct implementation: seven nested loops
Before optimizing convolution, write down the work it must perform. The following direct version processes a batch, meaning several images evaluated together. It uses stride 1 and assumes the input already has a suitable border, so every window stays in bounds and the output has IN_H × IN_W positions:
float input[BATCH][IN_H+FILTER_Y-1][IN_W+FILTER_X-1][IN_DEPTH]; // padded input float output[BATCH][IN_H][IN_W][NUM_FILTERS]; // output activations float weights[NUM_FILTERS][FILTER_Y][FILTER_X][IN_DEPTH]; float biases[NUM_FILTERS]; for (int img=0; img<BATCH; img++) for (int j=0; j<IN_H; j++) for (int i=0; i<IN_W; i++) for (int f=0; f<NUM_FILTERS; f++) { float tmp = biases[f]; for (int kk=0; kk<IN_DEPTH; kk++) // sum over input channels for (int jj=0; jj<FILTER_Y; jj++) // spatial convolution (Y) for (int ii=0; ii<FILTER_X; ii++) // spatial convolution (X) tmp += weights[f][jj][ii][kk] * input[img][j+jj][i+ii][kk]; output[img][j][i][f] = tmp; }
The outer four loops select an image, output row, output column, and filter. The inner three visit all input channels and spatial positions of that filter. For each output, tmp starts at biases[f] and accumulates one product per filter weight. The code describes the answer clearly, but it leaves performance decisions to the compiler and cache. Neighboring output windows overlap, and the next filter reads the same input window again. Keeping those values near the arithmetic units could avoid many slow reloads. How can we arrange this reuse systematically, across many layer shapes?
Convolution as matrix multiplication ("explicit GEMM")
GEMM stands for general matrix-matrix multiplication. Its fast implementations already know how to reuse data from registers and caches, so we would like to express convolution as a GEMM. Start with one input channel and one 3 × 3 filter. Copy the nine input values needed by each output pixel into one row of a new matrix. The filter's nine weights form a column. Their matrix product gives every output pixel. This rearrangement is called im2col ("image to columns"):
Now give every output pixel one row, including pixels from other images in the batch. Put each filter in a separate column of the weight matrix. An input with C channels contributes R·S·C values to each row, since a filter spans all channels. Using NVIDIA's symbols, let the batch contain N images, the convolution output have P × Q positions, and there be K filters of size R × S × C:
The difficulty is the intermediate matrix A. An interior input value appears in about R·S overlapping windows, so A can be about R·S times as large as the original activations. It must also be written out and read back by GEMM. In the lecture's example, 32 images of 256 × 256 pixels with 128 fp32 channels occupy 1 GiB; their 3 × 3 im2col matrix occupies about 9 GiB. That is a large temporary buffer and extra memory traffic, though the total traffic increase depends on the cache behavior, weights, and output size.
A 1 × 1 convolution has no overlapping spatial window to copy. Conceptually, the activation tensor is already an (N·P·Q) × C matrix, and the filters form a C × K matrix. With a compatible memory layout, GEMM can use that view without a separate im2col buffer. This is one reason MobileNet's many 1 × 1 convolutions are attractive to optimized libraries.
Making GEMM fast: blocking
Writing a matrix product as a loop is easy; making it reuse data is the hard part. Start with a straightforward loop that computes one element of C at a time:
float A[M][K], B[K][N], C[M][N]; // compute C += A * B #pragma omp parallel for for (int j=0; j<M; j++) for (int i=0; i<N; i++) for (int k=0; k<K; k++) C[j][i] += A[j][k] * B[k][i];
Different values of j can be handled in parallel, and the innermost loop forms a dot product. But the loop order gives a large matrix little reliable reuse of B: as k changes, B[k][i] jumps between rows, and another output element may revisit it only much later. A can sometimes remain cached across neighboring outputs, and a compiler can hold C[j][i] in a register. If we pessimistically charge two fresh fp32 operand loads to every multiply-add, that is 2 floating-point operations for 8 bytes, or 0.25 FLOPs per byte. Actual traffic depends on the cache and compiler; the point is that the code does not deliberately organize reuse.
Blocking, or tiling, makes that reuse deliberate. Instead of finishing one output at a time, multiply a small block of A by a small block of B and accumulate into a small block of C. Choose the blocks so they stay in a fast storage level while many operations use them. This sketch assumes each matrix dimension is a multiple of its block size; production code also handles smaller edge blocks:
for (int jblock=0; jblock<M; jblock+=BLOCKSIZE_J) for (int iblock=0; iblock<N; iblock+=BLOCKSIZE_I) for (int kblock=0; kblock<K; kblock+=BLOCKSIZE_K) // multiply one block of A by one block of B, accumulate into one block of C for (int j=0; j<BLOCKSIZE_J; j++) for (int i=0; i<BLOCKSIZE_I; i++) for (int k=0; k<BLOCKSIZE_K; k++) C[jblock+j][iblock+i] += A[jblock+j][kblock+k] * B[kblock+k][iblock+i];
For square b × b blocks, a step reads b² values from A and b² from B, then performs b³ multiply-adds. Each loaded input can contribute to b outputs before the block is replaced. Counting two FLOPs per multiply-add and four bytes per fp32 input gives 2b³ / (8b²) = b/4 FLOPs per byte of A/B traffic. This simple estimate leaves out loading and storing C and the traffic at other memory levels, but shows why larger reusable blocks raise arithmetic intensity.
Should BLOCKSIZE simply be as large as possible? The roofline's ridge point is the intensity at which the compute and bandwidth ceilings meet; reaching it is a useful target. Beyond that point, larger blocks do not raise this model's peak throughput. They can also crowd out other data, leave too few blocks to occupy all cores, or waste work at matrix edges. The best size depends on the machine and matrix shape. Notice that faster tensor-core arithmetic moves the ridge point to the right: faster math needs more data reuse to stay busy.
Blocking at every level, and vectorizing the block
A real implementation usually blocks at several levels: larger tiles for cache or GPU shared memory, smaller ones for a nearer cache, and tiny accumulators in registers. Each level tries to use a value several times before fetching it from the next, slower level.
Inside the innermost block the code must also use SIMD. The lecture shows three ways, and the right one depends on the block's shape:
- Vectorize the i loop (across a row of C). Broadcast one element of A to all lanes (
splat), multiply by a vector of B's row, accumulate into a vector of C. Good spatial locality in B, but the working set grows by the SIMD width. Needs a wide i dimension. - Vectorize the k loop (a dot product). If the i dimension is small, pre-transpose the block of B so that a row of A and a "row" of Bᵀ are both contiguous, then use SIMD dot products.
- Compute a SIMD_WIDTH × SIMD_WIDTH tile of C in registers, with pre-transposed A and C blocks, so the innermost operations are independent of each other (no serial dependency on one accumulator).
Now look back at MobileNet's layer list: the matrix shapes change from layer to layer (112 × 112 images with 32 channels at the start, 7 × 7 images with 1,024 channels at the end). Different layers of one network may need different schedules, which is exactly why writing DNN libraries is hard work.
Implicit GEMM, and giving the GPU enough work
Explicit GEMM gets a regular matrix operation but first writes the large im2col matrix to memory. Implicit GEMM keeps the matrix as a virtual view: when the GEMM kernel needs an entry of A, it calculates which image, location, and input channel that entry represents. It reads the original activation at that address instead of reading a stored copy of A.
- A naive version can calculate the addresses for every individual multiply-add. It avoids the large intermediate but still misses opportunities for reuse.
- A tiled version forms just the needed sub-block in on-chip shared memory, then uses it for many multiply-adds. The full im2col matrix is never stored off chip. Overlapping tiles may still reread some original activations, so the exact DRAM traffic depends on tiling and cache behavior.
The idea is close to Lecture 7's shared-memory convolution tile: bring in a small working set, reuse it, then move to the next tile. NVIDIA's CUTLASS provides building blocks for these kernels, including tiled matrix operations and iterators that map matrix coordinates back to convolution tensors. A library can tune the tile size and instruction shape for the particular layer and GPU.
A GPU needs a lot of parallel work
Good reuse alone is not enough on a GPU: many blocks must be ready at once to occupy its 80 streaming multiprocessors. The lecture compares two V100 convolution workloads, both with 128 output filters. The number of outputs gives a rough sense of how much independent work the kernel can divide into blocks:
| Case | Outputs | Output data (fp32) |
|---|---|---|
| N = 1 image, P = Q = 64 | 64 × 64 × 128 × 1 = 524,288 | 2 MB |
| N = 32 images, P = Q = 256 | 256 × 256 × 128 × 32 = 268 M | 1 GB |
In the lecture's measurements, effective convolution throughput rises from roughly 20 to 85 TFLOP/s as the images or batch grow, then levels off. Those figures describe convolution implementations using the GPU's tensor-core math; 85 TFLOP/s is above the V100's ordinary fp32 ALU peak. More input channels also create more work and potential reuse per output. A 1 × 1 filter has less spatial reuse than a larger filter, so its measured rate may be lower even though it performs fewer total operations. Batching often raises throughput by providing more work and reusing weights across images, but a service may wait longer to gather a batch, increasing request latency.
Algorithmic improvements: FFT and Winograd
Blocking preserves the same mathematical work while arranging data movement better. A different algorithm can sometimes reach the same convolution result with fewer multiplications. Two examples are useful because their trade-offs differ:
- FFT convolution. Transform the input and filter into the frequency domain, multiply matching frequency values, then transform the result back. The transforms cost roughly O(n log n) for a spatial region of n values, so their setup pays off mainly for sufficiently large filters or regions. Contributions from input channels can be summed before the inverse transform.
- Winograd convolution. Algebraically reorganize a small-filter convolution so a group of neighboring outputs shares intermediate products. It saves multiplications but introduces extra additions and transformations.
Here is the lecture's one-dimensional example. A three-weight filter slides across four input values to produce two outputs: y0 = x0w0 + x1w1 + x2w2, and y1 = x1w0 + x2w1 + x3w2. Computing those sums directly takes six multiplications and four additions. Winograd F(2, 3) forms four different products and recombines them:
m1 = (x0 − x2) · w0 m2 = (x1 + x2) · (w0 + w1 + w2)/2 // the filter factors depend only on w: m3 = (x2 − x1) · (w0 − w1 + w2)/2 // precompute them once per filter m4 = (x1 − x3) · w2 y0 = m1 + m2 + m3 y1 = m2 − m3 − m4
The input transformations and output combinations add work, but the four filter-dependent expressions can be prepared once and reused at many positions. This form uses four runtime products instead of six; the corresponding two-dimensional 2 × 2 output tile uses 16 instead of 36. That is 2.25× fewer multiplications, not 2.25× faster execution. Modern hardware may perform a multiply and an add together, so the extra additions, data movement, and rounding behavior all matter when choosing an algorithm.
Libraries pick the algorithm for you
This is why a library offers more than one convolution kernel. NVIDIA's cuDNN, Intel's oneDNN, and other vendor libraries provide implementations based on direct loops, GEMM, FFT, or Winograd. The fastest choice depends on dimensions, precision, device, workspace available for temporary data, and accuracy requirements. A library may use heuristics or benchmark candidates for a particular workload. Once a layer is fast, however, the memory traffic between layers becomes the next problem.
Memory traffic between operations: fusion
Imagine a convolution followed by scale and bias on each output value, then a 2 × 2 max pool. If each is a separate GPU kernel, convolution writes a full activation tensor to memory, scale-and-bias reads and rewrites it, and pooling reads it again. In the large example below, that tensor is about 1 GiB. The arithmetic in scale-and-bias is tiny compared with these transfers.
Fusion computes a consumer's result while its producer's values are still nearby. The first step is simple: apply scale and bias before storing each convolution result. The following line illustrates that first step only; it still writes the full-size output:
...
output[img][j][i][f] = tmp*scale + bias; // scale+bias fused: applied before the store
To fuse pooling too, the kernel must keep four scaled convolution results from one 2 × 2 region, take their maximum, and store only that pooled value. The full convolution output need not become a separate memory-resident array. Input and output activations still move, and implementations may have other traffic for weights or temporary state, but the large intermediate passes can disappear.
With ordinary fp32 arithmetic, a 3 × 3 convolution may dominate the sequence, so separate element-wise passes look relatively cheap. Make the convolution faster with tensor cores, or use a 1 × 1 filter that does less arithmetic per output, and those same memory passes take a much larger share of total time. Improving one stage can expose the cost of the next; fusion removes the exposed intermediate transfers.
- How would you fuse the 2 × 2 max pool into the conv loop nest?
Show one answer
Reorder the loops so each iteration produces a 2 × 2 group of output pixels for a filter (step
jandiby 2 and compute fourtmpvalues). Take the max of the four (after scale and bias) and store one value into a half-size output array. The full-size conv output is never written to memory. In a blocked GPU kernel, make each tile's width and height even, so every pooling window lies inside one tile.
A fusion trick for attention
Fusion can reach beyond two neighboring layers. In a transformer, attention lets each token in a sequence use information from other tokens. Each token has a query vector, a key vector, and a value vector. For N tokens and vectors of width d, collect them into three N × d matrices Q, K, and V. The simplified calculation for one attention head is:
S = Q Kᵀ // N × N: how much each token attends to each other token P = softmax(S) // softmax along each row of S O = P V // N × d
Each entry of S compares one query with one key. A row of P converts those scores into weights that sum to one, and the corresponding output row is a weighted sum of V. Real attention may scale the scores and mask disallowed positions; we omit both here to focus on data movement. The problem is that S and P each have N × N entries. A conventional sequence of three operations writes S, rereads it for softmax, writes P, then rereads P for the final multiply. It seems difficult to fuse these steps because a softmax weight depends on the maximum and sum of every score in its row.
First look at one score row x. Subtract its largest entry before exponentiating; this prevents overflow and does not change the normalized answer:
m(x) = maxi xi // subtract the max so exp never overflows
f(x) = [ e^(x1 − m(x)), …, e^(xN − m(x)) ]
l(x) = Σi f(x)i
softmax(x) = f(x) / l(x)
The key observation is that the maximum m and exponential sum l can be updated chunk by chunk. Suppose we have processed chunk 1 and kept its two numbers. Chunk 2 has its own maximum and sum. If the combined maximum changes, the old sum must be expressed relative to the new maximum before we add the new sum:
m(x) = max( m(x(1)), m(x(2)) ) l(x) = e^(m(x(1)) − m(x))·l(x(1)) + e^(m(x(2)) − m(x))·l(x(2))
For example, if the old maximum is 3 and the new one is 5, every old exponential is now measured relative to a value two larger, so its contribution is multiplied by e−2. That lets us find the final denominator with only two running numbers. It does not by itself produce every probability: to emit the full softmax row afterward, we would need to retain or recompute the scores. Fused attention avoids producing that row at all. Step through the running maximum and sum:
Fused attention
Now keep a third running quantity: a weighted value sum for each query row. For a score chunk with new combined maximum m, rescale the old weighted sum by e(old m − m), add this chunk's e(score − m)·V contributions, and divide the final weighted sum by the final l. This is the same rescaling used for the denominator, applied to vector-valued data.
A fused kernel handles a block Qi of query rows against successive blocks Kj and Vj. It forms a small score block QiKᵀj, updates the row maxima and sums, and folds that block's weighted V values into the running output. Once every key block has been visited, it writes Oi. S and P never need to be complete arrays in off-chip memory.
At N = 4,096, just one fp32 score matrix has 4,096² ≈ 16.8 million values, or 64 MiB. Multiple heads and sequences multiply that demand. Fusing attention removes those N² score intermediates from memory, at the cost of extra rescaling work and carefully sized on-chip tiles. For sufficiently long sequences, saving that memory traffic can outweigh the extra arithmetic. This is the central idea of FlashAttention (Dao et al., 2022).
Fusion in frameworks, and lower precision
Writing a custom kernel for every useful combination of layers would be an enormous task. Early libraries therefore provided a few popular fused patterns. More flexible systems take a graph of operations and generate a kernel for a supported fused subgraph. The cuDNN backend API is one route; compilers such as XLA, TVM, and Triton-based systems offer others. They can avoid large intermediate memory transfers, though compiling, launching, and selecting a kernel still have costs, and not every sequence can be fused profitably.
Halide's compute_at gives a concrete way to think about a schedule: compute a producer's values inside the consumer's loop, near the point where they are needed, instead of first storing the producer's entire output. The arithmetic meaning stays the same, but the order and location of computation change. That is the common thread between a hand-written fused kernel and a compiler-generated one.
Low precision
Data movement also falls when values use fewer bits. Many implementations use 16-bit floating point or 8-bit integers for some weights and activations; binary networks are a more extreme research direction. If values remain in the narrower format throughout a transfer, halving their width can halve those data bytes and storage. Real savings depend on conversions, packing, metadata, and whether the hardware has fast arithmetic for that format. Lower precision may also change the network's answers, so accuracy must be checked. Unit 5 explains fp16, bf16, int8, and quantization in more detail.
Summary of the optimization techniques
| Technique | Examples | Who does it |
|---|---|---|
| Better algorithms: better ML models | depth, filter width, number of filters per layer, stride; residual connections; depthwise separable convolutions | Designed manually, and increasingly by automatic search over topologies |
| Software optimization: good schedules and kernels | loop blocking/tiling, fusion, implicit GEMM, Winograd | Library authors and compilers, often with measured tuning |
| Approximation: compressing models | lower bit precision; automatic sparsification and pruning (not covered in this lecture) | Tools plus retraining |
Why a GPU is a good platform, and why it isn't the ideal one
Why use a GPU? Large convolution and matrix-multiply layers offer many independent outputs, regular control flow, and opportunities to reuse inputs. GPUs supply many arithmetic units, substantial memory bandwidth, and libraries such as cuDNN that know how to schedule these layers.
Why specialize further? A GPU is programmable across many workloads. Its warp instructions already share control across lanes, but fetching and decoding those instructions, moving operands through registers, and scheduling warps still cost energy. When most useful work is a repeating matrix operation, a specialized unit can amortize that control over a much larger piece of math.
Watch the amount of arithmetic represented by one instruction grow:
- a scalar fused multiply-add
a·x + b: one multiplication and one addition, counted as 2 floating-point operations; - a four-component dot product accumulated into a result: four multiplications and four additions, counted as 8 operations;
- a 4 × 4 matrix multiply-accumulate
A·B + C: 64 multiply-adds, counted as 128 operations.
Tensor cores, and dedicated DNN accelerators
NVIDIA's tensor cores apply this idea to small matrix tiles. A warp-level matrix instruction asks them to multiply blocks and accumulate the result, such as D = A·B + C. The A100 example in the lecture uses fp16 inputs with fp32 accumulation for the tensor-core figure. The table compares its theoretical peak rates with ordinary fp32 units:
| Per SM | Whole GA100 GPU (108 SMs) | Peak at 1.41 GHz |
|---|---|---|
| 64 fp32 multiply-add ALUs | 6,912 ALUs | 6,912 × 2 × 1.41 G = 19.5 TFLOP/s |
| 32 int32 ALUs | 3,456 | |
| 4 tensor cores | 432 tensor cores | 312 TFLOP/s (fp16 in, fp32 accumulate) |
Divide the tensor-core peak by 432 cores and 1.41 billion clocks per second: the idealized count is about 512 FLOPs, or 256 multiply-adds per tensor core per clock. That is the arithmetic in an 8 × 4 matrix multiplied by a 4 × 8 matrix and accumulated into an 8 × 8 result. The table's 312 versus 19.5 TFLOP/s ratio is about 16, but it compares mixed-precision tensor math with ordinary fp32 math, not two ways to compute every problem at the same precision. A layer reaches neither peak automatically: its matrix tiles, supported formats, data supply, and numerical tolerance all matter.
Going all the way: accelerators
A DNN accelerator pushes specialization further. Examples in the lecture include Google's TPU, Apple's Neural Engine, phone NPUs, Graphcore's IPU, and Cerebras's wafer-scale engine. Their designs differ, but many devote substantial area to matrix arithmetic and on-chip storage, with control designed around larger units of work. The trade-off is less flexibility than a general CPU or GPU. GPUs themselves have also moved in this direction by adding tensor cores.
Check yourself
- A convolution layer produces a 56 × 56 × 128 output from 64 input channels with 3 × 3 filters. How many multiply-adds and parameters (ignoring biases)?
Show answer
Multiply-adds: 56 × 56 × 128 × 64 × 9 = 231,211,008 (about 231 M). Parameters: 128 × 64 × 9 = 73,728. The parameter count doesn't depend on the image size.
- By what factor does a depthwise separable convolution reduce the work of a 3 × 3 convolution with 256 output channels?
Show answer
The cost ratio is 1/K + 1/(R·S) = 1/256 + 1/9 ≈ 0.115, so about 8.7× less work.
- Why does explicit im2col increase DRAM traffic, and how does implicit GEMM avoid it?
Show answer
Each interior input value appears in about R × S im2col rows, making the temporary matrix much larger than the input; writing and rereading that matrix adds off-chip traffic. Implicit GEMM forms only the sub-block currently needed in on-chip storage, so there is no full im2col buffer. Other tiles may still reread some original input values.
- In the blocked-GEMM model, three b × b fp32 blocks fit in 96 KB. What's the largest b, and the approximate arithmetic intensity?
Show answer
3 × b² × 4 bytes ≤ 98,304 gives b² ≤ 8,192, so the largest integer is b = 90 (90² = 8,100). Counting only reads of the A and B blocks, intensity is approximately b/4 ≈ 22 FLOPs per byte, above this teaching model's fp32 ridge point of about 17. Real kernels need room for more state and have additional traffic.
- Use Winograd F(2, 3) with x = (1, 2, 3, 4) and w = (1, 1, 1). What are m1…m4 and the outputs?
Show answer
m1 = (1 − 3)·1 = −2; m2 = (2 + 3)·3/2 = 7.5; m3 = (3 − 2)·1/2 = 0.5; m4 = (2 − 4)·1 = −2. y0 = −2 + 7.5 + 0.5 = 6 = 1 + 2 + 3; y1 = 7.5 − 0.5 + 2 = 9 = 2 + 3 + 4.
- A row is processed in two chunks. Chunk 1 has m = 3 and l = 2.0. Chunk 2 has m = 5 and l = 1.5. What are the combined m and l?
Show answer
m = 5. l = e3 − 5 × 2.0 + e0 × 1.5 = 0.1353 × 2.0 + 1.5 ≈ 1.771.
- Why does fusing an element-wise operation into a convolution matter more when the convolution runs on tensor cores?
Show answer
The element-wise pass costs the same memory time either way (read and write the whole tensor). Tensor cores shrink the conv's compute time, so the unfused memory passes become a larger fraction of the total, and can dominate it.
Based on Stanford CS149 (Fall 2023), Lecture 10, "Efficiently Evaluating DNNs", by Kayvon Fatahalian and Kunle Olukotun. The explanations, diagrams and interactive models on this page are original study material; the network statistics, V100/A100 figures and energy estimates are as presented in the lecture, which cites Howard et al. 2017 (MobileNet), Canziani et al. 2017, NVIDIA, and Bill Dally. Fused attention follows Dao et al., "FlashAttention" (2022). For the course's own slides and videos, see the CS149 course site.