Matrix Calculus: Gradients, Jacobians & Backprop

New to derivatives? A derivative just measures "if I nudge this input a tiny bit, how much does the output change, and in which direction?" A partial derivative is the same question asked about one input at a time, while every other input is frozen in place — "if I only nudge x (not y), how does the output move?" That's the entire prerequisite for this lesson; nothing more advanced is assumed.

Beginner: the gradient of a function with respect to a vector is just a vector holding every one of those partial derivatives — one number per input, all collected together. The gradient always points in the direction of steepest increase, which is exactly why gradient descent moves in the opposite direction.

Intermediate: the Jacobian generalizes this to functions that output a whole vector, not just a single number: it's a matrix where row i is the gradient of output i with respect to every input. The Hessian (section 1.16) is the Jacobian of the gradient itself — the matrix of all second derivatives.

Advanced: the chain rule, written in matrix form, says the Jacobian of a composition of functions is the product of their individual Jacobians. This one sentence is, precisely and without exaggeration, what backpropagation is: applying the chain rule layer by layer through a neural network, where "layer by layer" means "Jacobian by Jacobian."

"Gradient checking" — comparing an analytic gradient formula against a numeric finite-difference approximation — is a standard debugging technique whenever you implement backpropagation by hand.

#include <cstdio>
#include <vector>

using Matrix = std::vector<std::vector<double>>;
using Vector = std::vector<double>;

Matrix A = {{2., 0.}, {0., 5.}};

Vector matvec(const Matrix& A, const Vector& x) {
    Vector out(A.size(), 0.0);
    for (size_t i = 0; i < A.size(); ++i)
        for (size_t j = 0; j < x.size(); ++j) out[i] += A[i][j] * x[j];
    return out;
}

double dot(const Vector& a, const Vector& b) {
    double s = 0.0;
    for (size_t i = 0; i < a.size(); ++i) s += a[i] * b[i];
    return s;
}

double f(const Vector& x) {
    return dot(x, matvec(A, x));
}

Vector analyticGrad(const Vector& x) {
    Vector Ax = matvec(A, x);
    for (double& v : Ax) v *= 2.0;   // the identity above, since A is symmetric
    return Ax;
}

Vector numericGrad(const Vector& x, double eps = 1e-6) {
    Vector g(x.size(), 0.0);
    for (size_t i = 0; i < x.size(); ++i) {
        Vector xPlus = x, xMinus = x;
        xPlus[i] += eps;
        xMinus[i] -= eps;
        g[i] = (f(xPlus) - f(xMinus)) / (2 * eps);
    }
    return g;
}

int main() {
    Vector x = {1.5, -0.7};
    Vector ag = analyticGrad(x);
    Vector ng = numericGrad(x);
    std::printf("analytic: [%.6f, %.6f]\n", ag[0], ag[1]);
    std::printf("numeric:  [%.6f, %.6f]\n", ng[0], ng[1]);  // should match closely
    return 0;
}