Matrix Calculus for Deep Learning, Without the Fog

Every few months I need the same deep learning derivatives again. Not all of them. The same small cluster: affine maps, element-wise activations, Jacobians, chain rules, and the softmax-cross-entropy shortcut. I keep ending up back at Parr and Howard’s excellent article¹. This post is my own rewrite of a large chunk of that material, organized the way I wish my notes looked when I am half-awake and trying to remember why a gradient has the shape it has.

The thread running through all of it is simple:

• Geometry: a Jacobian is the best local linear map near the current point.
• Algebra: backprop is repeated composition of those local maps.
• Probability: the famous output-layer gradients are just “predicted distribution minus what reality said.”

The One Convention That Matters

With that convention, there are only a few cases to remember:

Most confusion in matrix calculus is not real math confusion. It is convention confusion. Once you settle the orientation, the rest is bookkeeping.

A Jacobian Is a Local Matrix

The Jacobian is not a bag of partial derivatives. It is the matrix that best describes how a nonlinear function behaves near the current point:

That is the entire geometric story. The Jacobian is the local linear approximation. If you zoom in far enough, a smooth nonlinear function looks linear, and the Jacobian is the matrix of that local linear map.

The columns of $J_f(x)$ are especially useful to keep in mind. The $j$-th column tells you what happens to the $j$-th basis direction $e_j$. If you know where the basis goes, you know where any vector goes, because any vector is a weighted sum of basis vectors.

That widget is literally the mental model I use for Jacobians. A matrix is a story about basis vectors. A Jacobian is the same story, but only locally.

Finite differences are a good way to keep yourself honest². If you ever derive a Jacobian on paper and want a quick check, approximate each column numerically:

import numpy as np

def jacobian_fd(f, x, eps=1e-6):
    x = np.asarray(x, dtype=float)
    y = np.asarray(f(x), dtype=float)
    J = np.zeros((y.size, x.size))
    for j in range(x.size):
        step = np.zeros_like(x)
        step[j] = eps
        J[:, j] = (f(x + step) - f(x - step)) / (2 * eps)
    return J

def f(x):
    x1, x2 = x
    return np.array([
        x1 * x2,
        np.sin(x1) + x2**2,
    ])

x = np.array([1.2, -0.7])
J_analytic = np.array([
    [x[1], x[0]],
    [np.cos(x[0]), 2 * x[1]],
])

np.testing.assert_allclose(
    jacobian_fd(f, x),
    J_analytic,
    rtol=1e-5,
    atol=1e-5,
)

Four Local Maps Worth Memorizing

Backprop is easier when you stop memorizing giant end-to-end formulas and memorize the local pieces instead.

1. Sum Reduction

If

then the gradient with respect to $x$ is just the all-ones vector:

Each coordinate contributes equally to the sum. Nothing subtle is happening.

2. Dot Product

If

then

This is one of the most reused facts in neural networks. A neuron pre-activation is just a dot product plus a bias term.

3. Affine Map

If

then the Jacobian of $z$ with respect to $x$ is simply

That is unsurprising: affine maps are already linear, up to a constant shift. The bias disappears in the derivative because constants do not change under infinitesimal perturbations.

4. Element-Wise Nonlinearity

If

where $\phi$ acts coordinate-wise, then the Jacobian is diagonal:

Why diagonal? Because each output coordinate depends only on the matching input coordinate. There is no cross-talk.

Put together, those four pieces cover a surprising amount of deep learning practice.

The Chain Rules That Power Backprop

Here is the version of the chain rule that shows up most often in neural nets.

Suppose

where $L$ is a scalar loss. Then

This is the backprop pattern in one line. A gradient flowing backward through a vector-valued function gets multiplied by the transpose of the local Jacobian.

Why the transpose? Because the Jacobian maps forward perturbations $\Delta x$ into forward perturbations $\Delta y$. Gradients live in the dual direction: they tell you how a scalar changes with respect to those perturbations. Backward flow is the transpose action.

For two vector-valued functions,

the forward Jacobian chain rule is

So the whole picture is:

• forward pass: multiply Jacobians in forward order
• backward pass: multiply transposed Jacobians in reverse order

If the function in the middle is element-wise, its Jacobian is diagonal, so the backward step becomes an element-wise multiplication:

where $\odot$ is the Hadamard product.

That is why backprop code is full of expressions like grad *= relu_prime(z) or dz = da * sigmoid_prime(z).

The easiest way to make that pattern concrete is to run one tiny network forward and then walk the scalar loss backward through each local map:

One Neuron, End to End

Take the smallest nontrivial neural building block:

with sigmoid

Assume a binary label $y \in \{0,1\}$ and binary cross-entropy loss

The derivative with respect to the activation is

The sigmoid derivative is

Multiply them:

That simplification is worth remembering. It is the binary cousin of the softmax-cross-entropy shortcut.

Once you have $\frac{\partial L}{\partial z}$, the rest is immediate:

This is both algebraically useful and probabilistically intuitive.

• If $y = 1$ but the neuron predicts $a = 0.1$, then $a-y = -0.9$. Gradient descent subtracts a negative quantity, pushing $z$ upward and making the positive class more likely.
• If $y = 0$ but the neuron predicts $a = 0.9$, then $a-y = 0.9$. Gradient descent pushes $z$ downward.

The sign is not arbitrary. It is the correction needed to move the model’s probability toward the observed outcome.

import numpy as np

x = np.array([1.5, -2.0, 0.7])
w = np.array([0.4, -0.3, 0.2])
b = -0.1
y = 1.0

z = w @ x + b
a = 1.0 / (1.0 + np.exp(-z))
loss = -(y * np.log(a) + (1.0 - y) * np.log(1.0 - a))

dz = a - y
dw = dz * x
db = dz

def loss_from_w(w_):
    z_ = w_ @ x + b
    a_ = 1.0 / (1.0 + np.exp(-z_))
    return -(y * np.log(a_) + (1.0 - y) * np.log(1.0 - a_))

eps = 1e-6
dw_fd = np.zeros_like(w)
for i in range(w.size):
    step = np.zeros_like(w)
    step[i] = eps
    dw_fd[i] = (loss_from_w(w + step) - loss_from_w(w - step)) / (2 * eps)

np.testing.assert_allclose(dw, dw_fd, rtol=1e-5, atol=1e-5)

Softmax + Cross-Entropy: The Famous Collapse

For multiclass classification, the output layer is usually

followed by cross-entropy with one-hot target vector $y$:

If the true class is $t$, this becomes

That sentence already tells you the probabilistic meaning. Cross-entropy is the surprisal of the probability assigned to the event that actually happened. Low probability on the true class means high loss.

The full softmax Jacobian is

Component-wise, that means

The diagonal terms are positive because raising a class’s own logit raises its probability. The off-diagonal terms are negative because probability mass is conserved: if one class goes up, others must come down.

The miracle is what happens after you compose that Jacobian with cross-entropy. The gradient with respect to logits collapses to

No giant expression. No explicit Jacobian matrix in your implementation. Just predicted probabilities minus target probabilities.

This formula is easier to remember if you read it probabilistically:

• $p$ is what the model currently believes.
• $y$ is what reality said.
• $p-y$ is the discrepancy.

Gradient descent acts on that discrepancy.

If the model gave too much mass to the wrong classes, those coordinates are positive and get pushed down. If it gave too little mass to the true class, that coordinate is negative and gets pushed up.

One clean derivation starts from

where the division is element-wise, and then applies the scalar-after-vector chain rule:

Carrying out the algebra gives the compact result $p-y$. In practice, most libraries fuse softmax and cross-entropy into one numerically stable operation because that compact form is both cheaper and more stable.

import numpy as np

def softmax(z):
    shifted = z - np.max(z)
    exp_z = np.exp(shifted)
    return exp_z / exp_z.sum()

z = np.array([2.0, -0.5, 0.2])
y = np.array([1.0, 0.0, 0.0])

def cross_entropy_from_logits(z_):
    p_ = softmax(z_)
    return -np.sum(y * np.log(p_))

p = softmax(z)
grad = p - y

eps = 1e-6
grad_fd = np.zeros_like(z)
for i in range(z.size):
    step = np.zeros_like(z)
    step[i] = eps
    grad_fd[i] = (
        cross_entropy_from_logits(z + step)
        - cross_entropy_from_logits(z - step)
    ) / (2 * eps)

np.testing.assert_allclose(grad, grad_fd, rtol=1e-5, atol=1e-5)

A Cheat Sheet I Actually Remember

If I had to compress the whole page into one sentence, it would be this:

Deep learning derivatives look hard only until you see them as local linear maps plus probability corrections.

That is what a Jacobian is. That is what the chain rule composes. That is what the output-layer gradients are doing.

Further Reading

If you want the denser reference version, go straight to Parr and Howard’s original write-up³. If you want a practical habit that pays off immediately, keep finite-difference checks around in small NumPy scratch files. They are fast, they are humbling, and they save a lot of algebraic overconfidence.