Gates → GPUs · Notes
Stanford CS149 · Lecture 10 companion

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.
Keep one question in mind: is this layer waiting for arithmetic units, or waiting for data? A matrix multiply can reuse each value many times if its implementation keeps that value nearby. An element-wise operation usually uses a value once, so running it as a separate pass may cost more in memory traffic than in arithmetic. The choices later in the lecture follow from that difference.
Part 1 · What a DNN computes

Layers are loops of multiply-adds

Just enough about neural networks to see the computation that has to run fast.

Section 1

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:

  1. Ask for less work. Choose a network topology that reaches the needed accuracy with fewer operations (Section 4).
  2. 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).
  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.

Section 2

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:

Four neurons, three inputs each. The weight matrix has 4 × 3 = 12 entries, plus one bias per neuron. Computing one input uses 12 weight-input products. If several input vectors are evaluated together, put them in columns: the layer becomes a matrix-matrix product, with opportunities to reuse each weight across the batch.

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.

Section 3

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.

Filter:

Click any output pixel. Its value is the sum of the nine input pixels under the window, each multiplied by the filter weight in the same position. Teal is positive, orange is negative. The input is a bright square on a dark background; blur softens its edges, the Sobel filters light up its left/right or top/bottom edges with opposite signs, and ReLU keeps only the positive responses.

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.

One possible block: conv → ReLU → pool. Under the diagram's same-size convolution and 2 × 2, stride-2 pool, W × H × K becomes (W/2) × (H/2) × K. The arrays between layers are activations. Many networks use other sequences or omit pooling.

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.

Section 4

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:

NetworkYearTop-1 accuracyParametersMADDs per image
VGG-16201471.5%138 M15 B
GoogLeNet201570%6.8 M1.5 B
ResNet-18201673%*11.7 M1.8 B
MobileNet-224201770.5%4.2 M0.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:

  1. Depthwise: C separate 3 × 3 × 1 filters, one per input channel. Spatial filtering, no channel mixing.
  2. Pointwise: K filters of size 1 × 1 × C. Channel mixing, no spatial extent.
Layer:
Output size H × W:
Input channels C:
Output channels K:
Filter size R × S:
The pointwise part dominates. The depthwise separable cost divided by the standard cost is exactly 1/K + 1/(R·S). With 3 × 3 filters and hundreds of channels that's about 1/9: 8–9 times less work. And the remaining work is almost all in the 1 × 1 convolutions.

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.

Operation count is not running time. A depthwise stage does little arithmetic per byte it reads, and some of its inputs may need to be fetched repeatedly. A GPU can therefore become limited by memory bandwidth or kernel overhead even after the multiply-add count falls. Count bytes and launches as well as arithmetic.
Part 2 · Making convolution fast

Approach 2: implement layers efficiently

Turn convolution into matrix multiply, then make matrix multiply use the memory hierarchy well.

Section 5

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?

Section 6

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"):

Click an output pixel or a row of the matrix.

im2col on a 4 × 4 image with zero padding. Each of the 16 output pixels gets a row of 9 values: its window, read left to right, top to bottom; positions outside the image are 0. Every interior input value appears in 9 different rows. Multiply this 16 × 9 matrix by the 9 filter weights and you get all 16 outputs.

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:

Conv layer → one GEMM. A[NPQ × RSC] holds one row per output pixel of every image; B[RSC × K] holds one column per filter; C = A·B is the output, one row per pixel and one column per output channel.

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.

Section 7

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.

Arithmetic:
Fast storage for the blocks:

A roofline for blocked GEMM. The machine is V100-like: 900 GB/s of memory bandwidth and either 15.7 TFLOP/s of fp32 or 125 TFLOP/s of fp16 tensor-core math. The model: three b × b blocks (A, B and C) must fit in the fast storage; a large 8192 × 8192 multiply is done block by block. The explorer counts reads and writes of C as well as the A/B traffic in the b/4 estimate above. Below the ridge point the multiply is bandwidth-bound; above it, compute-bound. (A teaching model: real kernels block at several levels at once and never get the full storage for one block.)

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:

  1. 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.
  2. 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.
  3. 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.

Section 8

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:

CaseOutputsOutput data (fp32)
N = 1 image, P = Q = 6464 × 64 × 128 × 1 = 524,2882 MB
N = 32 images, P = Q = 256256 × 256 × 128 × 32 = 268 M1 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.

Section 9

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.

Check it with your own numbers. Change any input or weight: Winograd's two outputs always match the direct ones. The slide prints the first product as (x0 − x1)·w0, which is a typo; switch to it and y0 comes out wrong unless x1 happens to equal x2.
Numerical results can change. Reordering floating-point arithmetic changes where rounding happens. Larger Winograd tiles can reduce multiplication count further but may magnify numerical error, especially at low precision. An implementation has to balance speed, temporary storage, and acceptable accuracy.

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.

Part 3 · Beyond one layer

Stop sending intermediate results through memory

Operator fusion, the fused-attention trick, and compilers that schedule for you.

Section 10

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.

Layer:
Math:

How much does fusion save in this model? The large layer has 32 images of 256 × 256, 128 input channels, 128 filters, and fp32 activations (about 1 GiB in and 1 GiB out of the convolution). The V100-like machine has 900 GB/s bandwidth. The tensor-core toggle raises the arithmetic ceiling while retaining fp32 activation traffic, isolating the resulting bottleneck shift. Each separate kernel is estimated by the larger of its arithmetic time and memory time; the model ignores launch costs and treats weight traffic as small. Orange marks an arithmetic limit and teal a memory limit.

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.

Exercisefrom the lecture
  1. 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 j and i by 2 and compute four tmp values). 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.

Section 11

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:

Chunk size:

Online softmax's running state. The running maximum m and exponential sum l are enough to update the normalization as chunks arrive. The explorer retains its eight example scores to display the final probabilities. A production fused-attention kernel instead consumes each chunk's weighted values immediately and keeps a running output accumulator.

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.

Keep small blocks and running state on chip. A block of Q, K, and V feeds two small matrix products and row reductions while the running output stays resident. The large N × N score and probability arrays are never stored. Q, K, V, and O still require storage proportional to N·d, and key/value blocks may be read again for other query blocks.

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).

Section 12

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

TechniqueExamplesWho does it
Better algorithms: better ML modelsdepth, filter width, number of filters per layer, stride; residual connections; depthwise separable convolutionsDesigned manually, and increasingly by automatic search over topologies
Software optimization: good schedules and kernelsloop blocking/tiling, fusion, implicit GEMM, WinogradLibrary authors and compilers, often with measured tuning
Approximation: compressing modelslower bit precision; automatic sparsification and pruning (not covered in this lecture)Tools plus retraining
Part 4 · Hardware

Hardware built for DNNs

Why GPUs are good at this, why they aren't ideal, and what specialized hardware changes.

Section 13

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.
Estimated control and delivery energy relative to arithmetic energy. These illustrative half-precision estimates, presented by Bill Dally in 2018, assign the overhead of instruction fetch, decode, control, and operand delivery as a percentage of useful arithmetic. An overhead of 2,000% means roughly twenty units of overhead energy for one unit of arithmetic energy in that model. A larger matrix instruction spreads those costs over more work; the percentages are not measurements of every GPU or workload.
Section 14

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 SMWhole GA100 GPU (108 SMs)Peak at 1.41 GHz
64 fp32 multiply-add ALUs6,912 ALUs6,912 × 2 × 1.41 G = 19.5 TFLOP/s
32 int32 ALUs3,456
4 tensor cores432 tensor cores312 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.

The lesson to carry forward: first decide how much work the network truly needs. Then keep reused values close to the arithmetic units, avoid storing intermediates that the next operation will consume immediately, and choose instructions and numeric formats suited to the calculation. Each step removes a different limit on speed.
Section 15

Check yourself

Questionsanswer before revealing
  1. 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.

  2. 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.

  3. 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.

  4. 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.

  5. 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.

  6. 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.

  7. 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.