Backpropagation

July 2026 – Vladislav Kruglikov

Recently, I implemented a small automatic differentiation engine in C++. In this post, I want to share some of the ideas behind it.

Gradient

For a scalar-valued function of several variables, the gradient at a point is a vector that points in the direction of steepest increase. Its magnitude is the maximum rate of increase at that point. Taking the gradient at every point gives a gradient vector field. At any particular point, one vector from that field describes the local direction and rate.

For a two-dimensional matrix $X$ and a scalar loss $L$ , every element $X_{ij}$ can affect the loss, so $\partial L / \partial X$ contains one partial derivative for each matrix element. Conceptually, we can flatten the independent scalar entries of $X$ into a vector. By definition, their derivatives form an ordinary gradient vector. The matrix-shaped gradient is a view of that vector, described by a shape and strides over its storage, and does not need to be contiguous.

Backpropagation therefore starts from a scalar objective, usually the loss. A vector-valued output has a Jacobian, not one gradient vector, and the vector itself does not define a unique meaning of «better». We first need a preference that combines its components into a scalar.

That scalar objective can encode the preference as a weighted combination. During backpropagation, its upstream gradient supplies those weights, reducing the output Jacobian to the single gradient direction relevant to the loss.

Jacobian

The gradient is the Jacobian of a scalar-valued function, up to whether derivatives are written as a row or a column.

Let $J_{ab}$ be the Jacobian of a vector $\mathbf{a}$ with respect to a vector $\mathbf{b}$ . Its entry in row $i$ and column $j$ is

(J_{ab})_{ij} = \frac{\partial a_i}{\partial b_j}

Therefore, the full Jacobian is

J_{ab} = \frac{\partial \mathbf{a}}{\partial \mathbf{b}} = \begin{bmatrix} \frac{\partial a_1}{\partial b_1} & \frac{\partial a_1}{\partial b_2} & \cdots & \frac{\partial a_1}{\partial b_n} \\ \frac{\partial a_2}{\partial b_1} & \frac{\partial a_2}{\partial b_2} & \cdots & \frac{\partial a_2}{\partial b_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial a_m}{\partial b_1} & \frac{\partial a_m}{\partial b_2} & \cdots & \frac{\partial a_m}{\partial b_n} \end{bmatrix}

If $\mathbf{c}$ depends on $\mathbf{a}$ and $\mathbf{a}$ depends on $\mathbf{b}$ , the Jacobians compose by matrix multiplication

J_{cb} = J_{ca}J_{ab}

To prove it, expand entry $(i,j)$ of the matrix product using $(J_{ab})_{ij}=\frac{\partial a_i}{\partial b_j}$

\begin{aligned} (J_{ca}J_{ab})_{ij} &= \sum_k (J_{ca})_{ik}(J_{ab})_{kj} \\ &= \sum_k \frac{\partial c_i}{\partial a_k} \frac{\partial a_k}{\partial b_j} \\ &= \frac{\partial c_i}{\partial b_j} \\ &= (J_{cb})_{ij} \end{aligned}

Every entry of $J_{ca}J_{ab}$ is therefore equal to the corresponding entry of $J_{cb}$ .

Now consider a sequence of dependencies

\mathbf{a}\to\mathbf{b}\to\mathbf{c}\to\mathbf{d}\to l

The Jacobian from $\mathbf{a}$ to $l$ is the product of the Jacobians along the path. Applying the same rule to the rightmost pair after each equals sign gives

\begin{aligned} J_{ld}J_{dc}J_{cb}J_{ba} &= J_{ld}J_{dc}J_{ca} \\ &= J_{ld}J_{da} \\ &= J_{la} \end{aligned}

If $l$ is scalar and $\mathbf{x}\in\mathbb{R}^n$ , the final Jacobian $J_{lx}$ has shape $1\times n$ . Its entries are $\frac{\partial l}{\partial x_i}$ , so this row vector contains the gradient information for $\mathbf{x}$ . Under the column-gradient convention,

\nabla_{\mathbf{x}}l = J_{lx}^\mathsf{T}

If $\mathbf{x}$ is a parameter vector, gradient descent updates it as

\mathbf{x}\leftarrow\mathbf{x}-\eta\nabla_{\mathbf{x}}l

Computing the gradient this way requires materializing every intermediate Jacobian, which can be a huge matrix.

These Jacobians are often sparse because one input usually affects only part of the output, leaving most entries equal to zero.

Consider a dependency chain $a \to b \to c \to d \to l$ , where $l$ is a scalar. Suppose we already know $J_{ld}$ . This is the upstream gradient. Because $l$ is scalar, $J_{ld}$ contains exactly one derivative for each element of $d$ . Mathematically it has shape $1\times |d|$ when $d$ is flattened, but in code those entries are stored or viewed with the same tensor shape as $d$ . How do we compute $J_{lc}$ ?

J_{lc} = J_{ld}J_{dc}

There are two approaches.

Materialize $J_{dc}$ and multiply $J_{ld}$ by it. This is memory-heavy because $J_{dc}$ can be huge and mostly sparse.
Compute the vector-Jacobian product directly, without forming $J_{dc}$ .

Automatic differentiation

Matrix multiplication can be decomposed into scalar multiplications and additions, often executed as fused multiply-add operations. Therefore, if we know how to backpropagate through scalar multiplication and addition, we can also backpropagate through matrix multiplication.

Doing this one scalar at a time would be extremely inefficient. Computing each output element would require traversing the graph and making many function calls, creating far more overhead than useful work. A much better approach is to represent the entire matrix multiplication as a single operation and execute it with a CUDA kernel built using CUTLASS. We can do exactly that while preserving automatic differentiation.

Deriving gradients by brute force

A general way to derive the backward pass for any operation is to start with every possible path from an input to the loss. Most of those paths usually have zero influence, so the full expression quickly collapses into something simple.

Consider a matrix multiplication

Y = XW

Let the upstream gradient be $G = \frac{\partial L}{\partial Y}$ , with the same shape as $Y$ . For one element $X_{ab}$ , the chain rule initially sums over every element of $Y$

\frac{\partial L}{\partial X_{ab}} = \sum_i \sum_j \frac{\partial Y_{ij}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{ij}}

Now substitute $Y_{ij} = \sum_k X_{ik}W_{kj}$ into the full expression

\frac{\partial L}{\partial X_{ab}} = \sum_i \sum_j \frac{\partial}{\partial X_{ab}} \left(\sum_k X_{ik}W_{kj}\right) \frac{\partial L}{\partial Y_{ij}}

First, move the derivative inside the sum over $k$

\frac{\partial L}{\partial X_{ab}} = \sum_i \sum_j \left( \sum_k \frac{\partial \left(X_{ik}W_{kj}\right)}{\partial X_{ab}} \right) \frac{\partial L}{\partial Y_{ij}}

Because $W$ does not depend on $X$ , we can move it outside the derivative

\frac{\partial L}{\partial X_{ab}} = \sum_i \sum_j \sum_k W_{kj} \frac{\partial X_{ik}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{ij}}

Next, split the sum over $i$ into the case $i=a$ and all cases where $i\ne a$

\frac{\partial L}{\partial X_{ab}} = \sum_j \sum_k W_{kj} \frac{\partial X_{ak}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{aj}} + \sum_{i\ne a}\sum_j\sum_k W_{kj} \frac{\partial X_{ik}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{ij}}

Every derivative in the second term is zero because $X_{ik}$ and $X_{ab}$ are different elements when $i\ne a$

\frac{\partial L}{\partial X_{ab}} = \sum_j \sum_k W_{kj} \frac{\partial X_{ak}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{aj}}

Now split the sum over $k$ into the case $k=b$ and all cases where $k\ne b$

\frac{\partial L}{\partial X_{ab}} = \sum_j W_{bj} \frac{\partial X_{ab}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{aj}} + \sum_j\sum_{k\ne b} W_{kj} \frac{\partial X_{ak}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{aj}}

Every derivative in the second term is zero because $X_{ak}$ and $X_{ab}$ are different elements when $k\ne b$ . Only the left term remains

\frac{\partial L}{\partial X_{ab}} = \sum_j W_{bj} \frac{\partial X_{ab}}{\partial X_{ab}} \frac{\partial L}{\partial Y_{aj}}

Finally, $\frac{\partial X_{ab}}{\partial X_{ab}}=1$ , so

\frac{\partial L}{\partial X_{ab}} = \sum_j W_{bj}\frac{\partial L}{\partial Y_{aj}}

This is already very close to matrix multiplication. To produce element $(a,b)$ , we need the dot product of row $a$ and column $b$ . Move the upstream gradient to the left, then rewrite $W_{bj}$ as $(W^\mathsf{T})_{jb}$

\frac{\partial L}{\partial X_{ab}} = \sum_j \frac{\partial L}{\partial Y_{aj}} (W^\mathsf{T})_{jb}

Now the indices have exactly the form of matrix multiplication. We take row $a$ of $\frac{\partial L}{\partial Y}$ and column $b$ of $W^\mathsf{T}$ . Therefore

\frac{\partial L}{\partial X} = GW^\mathsf{T}

We can therefore compute the gradient with a fast CUDA GEMM instead of traversing a graph of scalar multiplications and additions.

Gradient with Respect to the Weights

The gradient with respect to the weights follows exactly the same process. Start by summing every possible effect of $W_{ab}$ on the output

\frac{\partial L}{\partial W_{ab}} = \sum_i \sum_j \frac{\partial Y_{ij}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ij}}

Substitute $Y_{ij} = \sum_k X_{ik}W_{kj}$ into the full expression

\frac{\partial L}{\partial W_{ab}} = \sum_i \sum_j \frac{\partial}{\partial W_{ab}} \left(\sum_k X_{ik}W_{kj}\right) \frac{\partial L}{\partial Y_{ij}}

Move the derivative inside the sum over $k$

\frac{\partial L}{\partial W_{ab}} = \sum_i \sum_j \left( \sum_k \frac{\partial\left(X_{ik}W_{kj}\right)}{\partial W_{ab}} \right) \frac{\partial L}{\partial Y_{ij}}

Because $X$ does not depend on $W$ , we can move it outside the derivative

\frac{\partial L}{\partial W_{ab}} = \sum_i\sum_j\sum_k X_{ik} \frac{\partial W_{kj}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ij}}

First, split the sum over $j$ into the case $j=b$ and all cases where $j\ne b$

\frac{\partial L}{\partial W_{ab}} = \sum_i\sum_k X_{ik} \frac{\partial W_{kb}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ib}} + \sum_i\sum_{j\ne b}\sum_k X_{ik} \frac{\partial W_{kj}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ij}}

Every derivative in the second term is zero because $W_{kj}$ and $W_{ab}$ are different elements when $j\ne b$ . Only the left term remains

\frac{\partial L}{\partial W_{ab}} = \sum_i\sum_k X_{ik} \frac{\partial W_{kb}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ib}}

Now split the sum over $k$ into the case $k=a$ and all cases where $k\ne a$

\frac{\partial L}{\partial W_{ab}} = \sum_i X_{ia} \frac{\partial W_{ab}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ib}} + \sum_i\sum_{k\ne a} X_{ik} \frac{\partial W_{kb}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ib}}

Every derivative in the second term is zero because $W_{kb}$ and $W_{ab}$ are different elements when $k\ne a$ . Only the left term remains

\frac{\partial L}{\partial W_{ab}} = \sum_i X_{ia} \frac{\partial W_{ab}}{\partial W_{ab}} \frac{\partial L}{\partial Y_{ib}}

Finally, $\frac{\partial W_{ab}}{\partial W_{ab}}=1$ , so

\frac{\partial L}{\partial W_{ab}} = \sum_i X_{ia}\frac{\partial L}{\partial Y_{ib}}

Again, this is close to matrix multiplication. Rewrite $X_{ia}$ as $(X^\mathsf{T})_{ai}$

\frac{\partial L}{\partial W_{ab}} = \sum_i (X^\mathsf{T})_{ai} \frac{\partial L}{\partial Y_{ib}}

This is row $a$ of $X^\mathsf{T}$ multiplied by column $b$ of the upstream gradient. Therefore

\frac{\partial L}{\partial W} = X^\mathsf{T}G

Gradient of ReLU

ReLU is elementwise, so its backward rule is the same for a vector, matrix, or N-dimensional tensor. Conceptually, flatten the logical elements into a vector and apply the scalar rule independently to each element

y_i = \max(0, x_i)

\frac{\partial L}{\partial x_i} = \begin{cases} \frac{\partial L}{\partial y_i}, & x_i > 0 \\ 0, & x_i \le 0 \end{cases}

At $x_i=0$ , the implementation chooses the conventional subgradient $0$ . The result can then be viewed with the original shape. When the storage is contiguous, the implementation can process it as a flat contiguous buffer. Other tensor views may require respecting their strides.

Semantics of the backward pass

Each operation has a forward pass and a backward pass. The forward pass computes output tensors from input tensors and remembers the values needed to compute the local derivatives later.

Backpropagation starts from a scalar loss, whose derivative with respect to itself is $1$ .

The backward pass topologically sorts the computation graph and traverses it in reverse, from the result through the operations that produced it and toward the input tensors. Each operation consumes its accumulated upstream gradient, computes contributions to the gradients of its inputs, and accumulates those contributions.

Reverse topological order ensures that an operation is processed only after all of its downstream consumers have contributed. Its upstream gradient is therefore complete before the backward pass moves further toward the inputs.