📚TheoryIntermediate

Matrix Calculus

Key Points

•
Matrix calculus extends ordinary calculus to functions whose inputs and outputs are vectors and matrices, letting you compute gradients, Jacobians, and Hessians systematically.
•
The gradient ∇f(x) is the direction of steepest increase of a scalar function and is essential for optimization algorithms like gradient descent.
•
The Jacobian J captures all first-order partial derivatives of a vector function and explains how small input changes propagate to outputs.
•
The Hessian H captures second-order curvature information and is critical for Newton-type optimization and understanding convexity.
•
Chain rules in matrix form let you combine derivatives of composed functions efficiently, which is exactly what backpropagation does.
•
Trace tricks, vectorization, and Kronecker products turn messy matrix derivatives into clean linear-algebra expressions.
•
Automatic differentiation computes exact derivatives via the chain rule; forward mode computes J·v (JVP) efficiently, while reverse mode computes Jᵀ·v (VJP) efficiently.
•
In ML, matrix calculus is the language of backprop: gradients of losses w.r.t. weights and activations follow directly from these rules.

Prerequisites

→Linear algebra (vectors, matrices, transpose, matrix multiplication) — Matrix calculus manipulates derivatives using linear-algebra operations and identities.
→Single-variable and multivariable calculus — Understanding partial derivatives, chain rule, and second derivatives is essential before generalizing to matrices.
→Optimization basics (gradient descent, Newton’s method) — Motivates why gradients and Hessians matter and how they are used.
→Numerical stability and floating-point arithmetic — Gradient computations require care with rounding, step sizes, and stable transformations (e.g., log-sum-exp).

Detailed Explanation

Tap terms for definitions

01Overview

Matrix calculus is the extension of differential calculus to functions whose inputs and outputs are vectors and matrices. Instead of single numbers, we differentiate objects with many components at once. The main tools are the gradient (for scalar-valued functions), the Jacobian (for vector-valued functions), and the Hessian (for second derivatives of scalar functions). These organize all partial derivatives into structured arrays that reflect how an output changes as each input component is perturbed. For practical work, matrix calculus comes with a set of algebraic identities—often expressed using traces, vectorization, and Kronecker products—that let you manipulate derivatives concisely and avoid indexing headaches. In modern machine learning, matrix calculus is indispensable. It powers backpropagation by systematically applying chain rules through compositions of linear and nonlinear functions. It also underlies optimization algorithms, from basic gradient descent to Newton’s method and quasi-Newton methods. With these tools, you can derive update rules for parameters, analyze curvature, and reason about stability and convergence. Combined with automatic differentiation (AD), matrix calculus turns derivative computation into a robust, programmable pipeline that is both exact and efficient.

02Intuition & Analogies

Imagine steering a drone in 3D space. Pushing the joystick a little forward changes your position; pushing it sideways changes it differently. The gradient is like a smart arrow that tells you which joystick direction makes your altitude (a scalar function) increase fastest. If your output is multi-dimensional—say, altitude, temperature, and battery drain—the Jacobian is like a panel of gauges, each showing how one output reacts to nudging each joystick axis. It’s a full sensitivity map. Curvature brings another layer: picture a hilly landscape. The Hessian is like a local terrain map that tells you how steepness itself changes as you step in different directions. In valleys it points you inward, on ridges it shows how the slope flips quickly from up to down. If the Hessian is positive definite, you’re at a bowl-shaped region—great for Newton’s method, which uses that curvature to jump toward the minimum. Matrix calculus identities are like shortcuts for avoiding messy coordinate-by-coordinate bookkeeping. Instead of differentiating each entry of a matrix product one at a time, you use a trace identity that turns the derivative into a tidy transpose. Vectorization and Kronecker products are like flattening and tiling tricks: they convert matrix problems into big linear systems where standard linear algebra does the heavy lifting. Automatic differentiation is your autopilot for derivatives. Forward mode tells you, given a push in input space (a direction v), how the outputs change (J·v). Reverse mode runs sensitivity backward: if you care about how a final scalar changes, it pushes that sensitivity back to all inputs (Jᵀ·v)—exactly what backprop does in neural nets.

03Formal Definition

Let f:

R^{n}

\to

R

be differentiable. The gradient of f at x

\in

R^{n}

is the column vector

\nabla

f(x)

\in

R^{n}

whose i-th entry is

\frac{\partial f}{\partial x _{i}}

. For a vector-valued function g:

R^{n}

\to

R^{m}

, the Jacobian

J_{g}

(x)

\in

R^{m \times n}

has entries [

J_{g}

(x)]_{ij} =

\frac{\partial g _{i}}{\partial x _{j}}

. If f:

R^{n}

\to

R

is twice differentiable, the Hessian

H_{f}

(x)

\in

R^{n \times n}

has entries [

H_{f}

(x)]_{ij} =

\frac{\partial ^{2} f}{\partial x _{i} \partial x _{j}}

. When f is sufficiently smooth (Clairaut’s theorem),

H_{f}

is symmetric. For matrix functions, we often represent differentials using the trace inner product: for scalar f(X), df =

⟨

\nabla_{X}

f, dX

⟩_{F}

tr

((

\nabla_{X}

f)^

⊤

dX). This avoids component-wise indices. The vectorization operator

vec

(

\cdot

) stacks a matrix into a long vector; the Kronecker product

\otimes

lets us translate between matrix equations and their vectorized forms, e.g.,

vec

(A X B) = (B^

⊤

\otimes

vec

(X). Chain rules generalize: if

y = g (x

) and f depends on y, then

\nabla_{x}

f = J_{g} (x

⊤

\nabla_{y}

f and, for second order,

\nabla_{x}^{2}

f = J_{g} (x

⊤

\nabla_{y}^{2}

f\,

J_{g}

(x) +

\sum_{k}

(

\nabla_{y}

f)_k

\nabla_{x}^{2}

g_{k}

(x). These formulas encode how sensitivities propagate through compositions, the backbone of backpropagation.

04When to Use

Use matrix calculus whenever your model parameters or variables are vectors/matrices and you need derivatives for optimization, sensitivity analysis, or uncertainty quantification.

Training ML models: Derive gradients of loss functions with respect to weights and biases in linear models, logistic regression, and neural networks. Reverse-mode (backprop) is ideal when outputs are few (often scalar loss) and inputs are many (parameters).
Algorithm design: Build Newton, Gauss–Newton, and quasi-Newton methods, which require Hessians or Hessian-vector products to exploit curvature for faster convergence.
Probabilistic modeling: Differentiate log-likelihoods, log-partition functions, and expectations. Use identities like d log det and trace rules for Gaussian models and graphical models.
Control and robotics: Compute Jacobians for kinematics and dynamics to linearize systems, perform state estimation, and propagate uncertainties.
Scientific computing: Sensitivity of PDE solvers or simulators with respect to parameters, where JVPs/VJPs enable efficient derivative propagation without forming dense Jacobians explicitly.
Implementation choice: Use forward-mode AD for small input dimension (cheap full Jacobian or J·v) and reverse-mode AD for small output dimension (cheap gradients and VJPs).

⚠️Common Mistakes

Shape/orientation confusion: Mixing row and column conventions leads to wrong transposes. Pick a consistent convention (column vectors by default) and check dimensions at every step.
Missing the symmetric factor: For f(x)=x^\top A x, the gradient is (A + A^\top)x, not Ax + x^\top A. If A is symmetric, it simplifies to 2Ax; don’t blindly write 2Ax when A is not symmetric.
Ignoring non-commutativity: Matrix multiplication does not commute. You cannot reorder terms inside derivatives; respect original order and apply trace tricks carefully.
Dropping chain-rule transposes: The gradient of a scalar function through y=g(x) is J_g(x)^\top \nabla_y f, not J_g(x) \nabla_y f. The transpose is essential to match shapes.
Confusing Jacobian with gradient: For f: \mathbb{R}^n \to \mathbb{R}, the Jacobian is a 1\times n row, while the gradient is an n\times 1 column. Many texts use different conventions—be explicit.
Numerical vs analytical mismatch: Finite-difference checks require the same vector/matrix shapes and careful step sizes; too-large steps give bias, too-small steps cause round-off error.
Overlooking stability: Expressions like log(\sum_i e^{z_i}) should use the log-sum-exp trick to avoid overflow; derivative code must mirror the numerically stable formula.
Forming dense Jacobians unnecessarily: For large problems, prefer JVP/VJP routines over materializing J; it saves memory and time.

Key Formulas

Gradient

\nabla f (x) = \frac{\partial f}{\partial x _{1}} ⋮ \frac{\partial f}{\partial x _{n}}

Explanation: The gradient stacks all first derivatives of a scalar function into a column vector. It points in the direction of steepest ascent.

Jacobian

J_{f} (x) = [\frac{\partial f _{i}}{\partial x _{j}}]_{i = 1.. m, j = 1.. n}

Explanation: The Jacobian collects partial derivatives for vector-valued functions. It maps input perturbations to output perturbations linearly at first order.

Hessian

H_{f} (x) = [\frac{\partial ^{2} f}{\partial x _{i} \partial x _{j}}]_{i, j = 1.. n}

Explanation: The Hessian captures curvature of a scalar function. It determines how the gradient changes with position.

Scalar Chain Rule

\frac{\partial f}{\partial x} = \frac{\partial f}{\partial y} \cdot \frac{\partial y}{\partial x}

Explanation: For scalar dependencies through an intermediate y, derivatives multiply according to the chain rule. Shapes must be aligned appropriately.

Matrix Chain Rule (Gradient)

\nabla_{x} f (g (x)) = J_{g} (x)^{⊤} \nabla_{y} f (y)

Explanation: When a scalar function depends on y=g(x), the gradient w.r.t. x is the Jacobian transpose times the gradient w.r.t. y. This is the core of backpropagation.

Quadratic Form Gradient

\nabla_{x} (x^{⊤} A x) = (A + A^{⊤}) x

Explanation: The gradient of a quadratic form is the symmetric part of A times x. If A is symmetric, this simplifies to 2Ax.

Trace Derivative 1

\frac{\partial o p er a t or nam e t r ( A ^{⊤} X )}{\partial X} = A

Explanation: Viewing df as tr( $G^{T}$ dX), this identity provides a clean way to read off the gradient with respect to X when $f = t r (A^{T}$ X).

Trace Derivative 2

\frac{\partial o p er a t or nam e t r ( A B )}{\partial A} = B^{⊤}

Explanation: Differentiating tr(AB)=tr(B A) with respect to A yields $B^{T}$ . It is often used to move derivatives through products.

Vec–Kronecker Identity

o p er a t or nam e v ec (A XB) = (B^{⊤} \otimes A) o p er a t or nam e v ec (X)

Explanation: Vectorizing a bilinear form turns it into a linear form using a Kronecker product. This is helpful to express matrix derivatives as large linear systems.

Log-Det Differential

d l o g d e tX = o p er a t or nam e t r (X^{- 1} d X)

Explanation: The differential of the log-determinant is the trace of $X^{- 1}$ times dX. It leads to gradients like ∂/∂X log det X = ( $X^{- 1}$ )^ $⊤$ .

Frobenius Inner Product

l an g l e A, B r an g l e_{F} = o p er a t or nam e t r (A^{⊤} B)

Explanation: This identity defines the inner product for matrices and underlies many trace-based derivative manipulations.

JVP & VJP

o p er a t or nam e J V P (v) = J (x) v, o p er a t or nam e V J P (u) = J (x)^{⊤} u

Explanation: Forward-mode AD computes J·v efficiently for a chosen direction v, while reverse-mode computes $J^{T}$ ·u efficiently for a chosen output weighting u.

Hessian–Vector Product

nab l a^{2} f (x); v = l e f t . f r a c d d e p s i l o n nab l a f (x + e p s i l o n v) r i g h t ∣_{e p s i l o n = 0}

Explanation: You can compute H·v via the directional derivative of the gradient, enabling curvature-aware methods without forming H explicitly.

Cost of Quadratic Form

T (n) = O (n^{2}) for x^{t o p} A x evaluation, T (n) = O (n^{2}) for (A + A^{t o p}) x

Explanation: Computing a quadratic form or its gradient uses matrix–vector products that scale quadratically with dimension n.

Complexity Analysis

Computational cost depends on both algebraic structure and the AD mode used. For a dense quadratic form f(x)=

x^{T}

A x with A∈

R^{n \times n}

, evaluating f requires one matrix–vector product Ax (O(

n^{2}

)) and a dot product

x^{T}

(Ax) (O(n)), yielding O(

n^{2}

) time and O(

n^{2}

) space if A is stored densely. Its gradient (A+

A^{T}

)x also costs O(

n^{2}

). The Hessian is constant (A+

A^{T}

), but explicitly storing it costs O(

n^{2}

). For a general vector function g:

R^{n}

→

R^{m}

, explicitly forming the Jacobian requires m·n partials. Finite-difference approximations need n function calls for forward differences (or 2n for central), each call costing

T_{f}

, giving O(n·

T_{f}

) to get one row/column pattern; a full Jacobian is O(n·

T_{g}

) or O(m·

T_{g}

) depending on how you sweep, and suffers from numerical error. Automatic differentiation is exact up to floating-point rounding and exploits program structure. Forward-mode propagates pairs (value, derivative) alongside computation and returns J·v in time O(

T_{p} ro g

), essentially a constant-factor overhead relative to the primal evaluation. To obtain the full Jacobian with forward mode, repeat with seed directions

e_{i}

, costing O(n·

T_{p} ro g

). Reverse-mode records a computation trace in the forward pass and then accumulates sensitivities backward; it computes

J^{T}

·u in O(

T_{p} ro g

). For a scalar output (

m = 1

), a single reverse sweep gives the full gradient in O(

T_{p} ro g

), while a full Jacobian would need m reverse sweeps, O(m·

T_{p} ro g

). Memory for reverse mode includes the tape of intermediates, often O(

T_{p} ro g

) in size; forward mode needs only local derivatives, typically O(n) or O(1) extra per variable depending on implementation. For large dense problems, avoid explicitly forming Jacobians/Hessians; prefer JVP/VJP and Hessian–vector products, reducing both time and memory overheads.

Code Examples

Analytic vs. numerical gradient and Hessian for a quadratic form

1 #include <iostream>
2 #include <vector>
3 #include <random>
4 #include <iomanip>
5 #include <cmath>
6 #include <cassert>
7 
8 // Compute f(x) = x^T A x + b^T x (A is n x n, x and b are n)
9 double f_quadratic(const std::vector<double>& A, const std::vector<double>& x, const std::vector<double>& b) {
10     size_t n = x.size();
11     // y = A x
12     std::vector<double> y(n, 0.0);
13     for (size_t i = 0; i < n; ++i)
14         for (size_t j = 0; j < n; ++j)
15             y[i] += A[i*n + j] * x[j];
16     // x^T y
17     double quad = 0.0;
18     for (size_t i = 0; i < n; ++i) quad += x[i] * y[i];
19     // b^T x
20     double lin = 0.0;
21     for (size_t i = 0; i < n; ++i) lin += b[i] * x[i];
22     return quad + lin;
23 }
24 
25 // Analytic gradient: ∇f = (A + A^T) x + b
26 std::vector<double> grad_analytic(const std::vector<double>& A, const std::vector<double>& x, const std::vector<double>& b) {
27     size_t n = x.size();
28     std::vector<double> g(n, 0.0);
29     for (size_t i = 0; i < n; ++i) {
30         for (size_t j = 0; j < n; ++j) {
31             g[i] += (A[i*n + j] + A[j*n + i]) * x[j];
32         }
33         g[i] += b[i];
34     }
35     return g;
36 }
37 
38 // Hessian: H = A + A^T (constant)
39 std::vector<double> hessian(const std::vector<double>& A) {
40     size_t n = static_cast<size_t>(std::sqrt((double)A.size()));
41     std::vector<double> H(n*n, 0.0);
42     for (size_t i = 0; i < n; ++i)
43         for (size_t j = 0; j < n; ++j)
44             H[i*n + j] = A[i*n + j] + A[j*n + i];
45     return H;
46 }
47 
48 // Central-difference numerical gradient
49 std::vector<double> grad_numeric(const std::vector<double>& A, const std::vector<double>& x, const std::vector<double>& b, double eps=1e-6) {
50     size_t n = x.size();
51     std::vector<double> g(n, 0.0);
52     std::vector<double> x1 = x, x2 = x;
53     for (size_t i = 0; i < n; ++i) {
54         x1[i] += eps; x2[i] -= eps;
55         double f1 = f_quadratic(A, x1, b);
56         double f2 = f_quadratic(A, x2, b);
57         g[i] = (f1 - f2) / (2.0 * eps);
58         x1[i] = x[i]; x2[i] = x[i];
59     }
60     return g;
61 }
62 
63 int main() {
64     std::mt19937 rng(123);
65     std::uniform_real_distribution<double> dist(-1.0, 1.0);
66 
67     size_t n = 5;
68     std::vector<double> A(n*n), x(n), b(n);
69     for (double &v : A) v = dist(rng);
70     for (double &v : x) v = dist(rng);
71     for (double &v : b) v = dist(rng);
72 
73     // Compute values
74     double fx = f_quadratic(A, x, b);
75     auto gA = grad_analytic(A, x, b);
76     auto gN = grad_numeric(A, x, b);
77     auto H  = hessian(A);
78 
79     // Report max absolute gradient error
80     double max_err = 0.0;
81     for (size_t i = 0; i < n; ++i) max_err = std::max(max_err, std::abs(gA[i] - gN[i]));
82 
83     std::cout << std::fixed << std::setprecision(6);
84     std::cout << "f(x) = " << fx << "\n";
85     std::cout << "Max |grad_analytic - grad_numeric| = " << max_err << "\n";
86 
87     // Sanity: check Hessian symmetry
88     double anti_sym = 0.0;
89     for (size_t i = 0; i < n; ++i)
90         for (size_t j = 0; j < n; ++j)
91             anti_sym = std::max(anti_sym, std::abs(H[i*n + j] - H[j*n + i]));
92     std::cout << "Max |H - H^T| = " << anti_sym << "\n";
93 
94     return 0;
95 }
96

We differentiate a quadratic objective f(x)=x^T A x + b^T x. The analytic gradient is (A+A^T)x + b and the Hessian is A+A^T. We verify the gradient by comparing it to a central-difference numerical estimate. The code also checks the Hessian’s symmetry. This example illustrates how matrix calculus converts multi-index derivatives into compact linear-algebra expressions and how to validate them numerically.

Time: O(n^2) to evaluate f and analytic ∇f, O(n^3) for numerical gradient via central differences (n perturbations × O(n^2) evaluation).Space: O(n^2) to store A and H; O(n) additional working memory.

Forward-mode automatic differentiation (JVP) with dual numbers

1 #include <iostream>
2 #include <vector>
3 #include <cmath>
4 #include <iomanip>
5 
6 // Minimal dual number for forward-mode AD (single directional derivative)
7 struct Dual {
8     double v; // value
9     double d; // derivative along a seeded direction
10     Dual(double v_=0.0, double d_=0.0) : v(v_), d(d_) {}
11 };
12 
13 // Overloaded operators for Dual
14 Dual operator+(const Dual& a, const Dual& b){ return Dual(a.v+b.v, a.d+b.d); }
15 Dual operator-(const Dual& a, const Dual& b){ return Dual(a.v-b.v, a.d-b.d); }
16 Dual operator*(const Dual& a, const Dual& b){ return Dual(a.v*b.v, a.v*b.d + a.d*b.v); }
17 Dual operator/(const Dual& a, const Dual& b){ return Dual(a.v/b.v, (a.d*b.v - a.v*b.d)/(b.v*b.v)); }
18 Dual operator-(const Dual& a){ return Dual(-a.v, -a.d); }
19 
20 // Elementary functions
21 Dual sin(const Dual& a){ return Dual(std::sin(a.v), std::cos(a.v)*a.d); }
22 Dual cos(const Dual& a){ return Dual(std::cos(a.v), -std::sin(a.v)*a.d); }
23 Dual exp(const Dual& a){ return Dual(std::exp(a.v), std::exp(a.v)*a.d); }
24 Dual log(const Dual& a){ return Dual(std::log(a.v), a.d / a.v); }
25 
26 // A vector-valued test function f: R^3 -> R^2
27 // y0 = x0 * x1 + sin(x2)
28 // y1 = exp(x0) + x1^2 - x2
29 std::vector<Dual> f(const std::vector<Dual>& x){
30     std::vector<Dual> y(2);
31     y[0] = x[0]*x[1] + sin(x[2]);
32     y[1] = exp(x[0]) + x[1]*x[1] - x[2];
33     return y;
34 }
35 
36 // Compute J*v at x using forward-mode AD by seeding the direction v
37 std::vector<double> jvp(const std::vector<double>& x0, const std::vector<double>& v){
38     size_t n = x0.size();
39     std::vector<Dual> x(n);
40     for(size_t i=0;i<n;++i) x[i] = Dual(x0[i], v[i]);
41     auto y = f(x);
42     std::vector<double> Jv(y.size());
43     for(size_t i=0;i<y.size();++i) Jv[i] = y[i].d;
44     return Jv;
45 }
46 
47 // Compute the full Jacobian by seeding basis vectors e_i (n forward sweeps)
48 std::vector<std::vector<double>> jacobian(const std::vector<double>& x0){
49     size_t n = x0.size();
50     auto y0 = f({Dual(x0[0],0), Dual(x0[1],0), Dual(x0[2],0)});
51     size_t m = y0.size();
52     std::vector<std::vector<double>> J(m, std::vector<double>(n, 0.0));
53     for(size_t j=0;j<n;++j){
54         std::vector<double> v(n, 0.0); v[j]=1.0; // seed e_j
55         auto Jv = jvp(x0, v);
56         for(size_t i=0;i<m;++i) J[i][j] = Jv[i];
57     }
58     return J;
59 }
60 
61 int main(){
62     std::vector<double> x0 = {0.3, -0.7, 1.2};
63     std::vector<double> v  = {1.0, 2.0, -0.5}; // direction for JVP
64 
65     auto Jv = jvp(x0, v);
66     auto J  = jacobian(x0);
67 
68     std::cout << std::fixed << std::setprecision(6);
69     std::cout << "J*v = [" << Jv[0] << ", " << Jv[1] << "]\n";
70     std::cout << "Jacobian at x0:\n";
71     for(const auto& row : J){
72         for(double a : row) std::cout << std::setw(12) << a;
73         std::cout << "\n";
74     }
75     return 0;
76 }
77

This forward-mode AD implementation uses dual numbers (value, derivative) to propagate derivatives through arithmetic and elementary functions automatically. Seeding a direction v in the input produces J·v (a Jacobian-vector product) at the output with one forward pass. Repeating with basis directions yields the full Jacobian. Forward mode is efficient when the input dimension is small or when you specifically need JVPs.

Time: O(T_prog) for one JVP (one forward sweep); O(n · T_prog) for the full Jacobian by repeating with n seed directions.Space: O(1) extra per active variable (stores value and one directional derivative).

Reverse-mode backprop for a linear layer: L = 1/2 ||W x − t||^2

1 #include <iostream>
2 #include <vector>
3 #include <random>
4 #include <cmath>
5 #include <iomanip>
6 
7 // Multiply W (m x n) by x (n)
8 std::vector<double> matvec(const std::vector<double>& W, size_t m, size_t n, const std::vector<double>& x){
9     std::vector<double> y(m, 0.0);
10     for(size_t i=0;i<m;++i)
11         for(size_t j=0;j<n;++j)
12             y[i] += W[i*n + j]*x[j];
13     return y;
14 }
15 
16 // Compute loss L = 1/2 ||y - t||^2
17 double loss(const std::vector<double>& y, const std::vector<double>& t){
18     double s=0.0; for(size_t i=0;i<y.size();++i){ double r=y[i]-t[i]; s+=0.5*r*r; } return s; }
19 
20 // Analytic gradients via matrix calculus rules:
21 // y = W x; r = y - t; L = 1/2 ||r||^2
22 // dL/dy = r
23 // dL/dx = W^T r
24 // dL/dW = r x^T
25 void analytic_grads(const std::vector<double>& W, size_t m, size_t n,
26                     const std::vector<double>& x, const std::vector<double>& t,
27                     std::vector<double>& gW, std::vector<double>& gx){
28     std::vector<double> y = matvec(W, m, n, x);
29     std::vector<double> r(m,0.0);
30     for(size_t i=0;i<m;++i) r[i] = y[i] - t[i]; // dL/dy
31     // dL/dx = W^T r
32     gx.assign(n,0.0);
33     for(size_t j=0;j<n;++j)
34         for(size_t i=0;i<m;++i)
35             gx[j] += W[i*n + j]*r[i];
36     // dL/dW = r x^T
37     gW.assign(m*n,0.0);
38     for(size_t i=0;i<m;++i)
39         for(size_t j=0;j<n;++j)
40             gW[i*n + j] = r[i]*x[j];
41 }
42 
43 // Finite-difference check for dL/dx
44 std::vector<double> numeric_grad_x(const std::vector<double>& W, size_t m, size_t n,
45                                    const std::vector<double>& x, const std::vector<double>& t, double eps=1e-6){
46     std::vector<double> gx(n,0.0), xp=x, xm=x;
47     for(size_t j=0;j<n;++j){
48         xp[j]+=eps; xm[j]-=eps;
49         double Lp = loss(matvec(W,m,n,xp), t);
50         double Lm = loss(matvec(W,m,n,xm), t);
51         gx[j] = (Lp-Lm)/(2.0*eps);
52         xp[j]=x[j]; xm[j]=x[j];
53     }
54     return gx;
55 }
56 
57 int main(){
58     size_t m=4, n=3;
59     std::mt19937 rng(42);
60     std::normal_distribution<double> N(0.0,1.0);
61     std::vector<double> W(m*n), x(n), t(m);
62     for(double& v:W) v=N(rng);
63     for(double& v:x) v=N(rng);
64     for(double& v:t) v=N(rng);
65 
66     std::vector<double> gW, gx;
67     analytic_grads(W,m,n,x,t,gW,gx);
68 
69     // Numerical check for gx
70     auto gx_num = numeric_grad_x(W,m,n,x,t);
71 
72     double max_err_x=0.0;
73     for(size_t j=0;j<n;++j) max_err_x = std::max(max_err_x, std::abs(gx[j]-gx_num[j]));
74 
75     std::cout << std::fixed << std::setprecision(6);
76     std::cout << "Max |dL/dx (analytic - numeric)| = " << max_err_x << "\n";
77 
78     // Print a few entries of dL/dW for illustration
79     std::cout << "Sample dL/dW entries:" << "\n";
80     for(size_t i=0;i<std::min<size_t>(m,2);++i){
81         for(size_t j=0;j<std::min<size_t>(n,3);++j){
82             std::cout << std::setw(12) << gW[i*n+j];
83         }
84         std::cout << "\n";
85     }
86 
87     return 0;
88 }
89

This example performs reverse-mode differentiation manually using matrix calculus identities for a linear layer with squared error. The forward pass computes y=W x and the loss L=1/2 ||y−t||^2. The backward pass applies chain rules to obtain dL/dy=r, dL/dx=W^T r, and dL/dW=r x^T. A finite-difference check validates dL/dx. This mirrors how backprop propagates sensitivities through linear layers in neural networks without materializing full Jacobians.

Time: O(mn) for matvec and gradient computations; O(n · (mn)) for numeric gradient check of x via central differences.Space: O(mn) to store W; O(m+n) for intermediates and gradients.

1	#include <iostream>
2	#include <vector>
3	#include <random>
4	#include <iomanip>
5	#include <cmath>
6	#include <cassert>
7
8	// Compute f(x) = x^T A x + b^T x (A is n x n, x and b are n)
9	double f_quadratic(const std::vector<double>& A, const std::vector<double>& x, const std::vector<double>& b) {
10	size_t n = x.size();
11	// y = A x
12	std::vector<double> y(n, 0.0);
13	for (size_t i = 0; i < n; ++i)
14	for (size_t j = 0; j < n; ++j)
15	y[i] += A[in + j] x[j];
16	// x^T y
17	double quad = 0.0;
18	for (size_t i = 0; i < n; ++i) quad += x[i] * y[i];
19	// b^T x
20	double lin = 0.0;
21	for (size_t i = 0; i < n; ++i) lin += b[i] * x[i];
22	return quad + lin;
23	}
24
25	// Analytic gradient: ∇f = (A + A^T) x + b
26	std::vector<double> grad_analytic(const std::vector<double>& A, const std::vector<double>& x, const std::vector<double>& b) {
27	size_t n = x.size();
28	std::vector<double> g(n, 0.0);
29	for (size_t i = 0; i < n; ++i) {
30	for (size_t j = 0; j < n; ++j) {
31	g[i] += (A[in + j] + A[jn + i]) * x[j];
32	}
33	g[i] += b[i];
34	}
35	return g;
36	}
37
38	// Hessian: H = A + A^T (constant)
39	std::vector<double> hessian(const std::vector<double>& A) {
40	size_t n = static_cast<size_t>(std::sqrt((double)A.size()));
41	std::vector<double> H(n*n, 0.0);
42	for (size_t i = 0; i < n; ++i)
43	for (size_t j = 0; j < n; ++j)
44	H[in + j] = A[in + j] + A[j*n + i];
45	return H;
46	}
47
48	// Central-difference numerical gradient
49	std::vector<double> grad_numeric(const std::vector<double>& A, const std::vector<double>& x, const std::vector<double>& b, double eps=1e-6) {
50	size_t n = x.size();
51	std::vector<double> g(n, 0.0);
52	std::vector<double> x1 = x, x2 = x;
53	for (size_t i = 0; i < n; ++i) {
54	x1[i] += eps; x2[i] -= eps;
55	double f1 = f_quadratic(A, x1, b);
56	double f2 = f_quadratic(A, x2, b);
57	g[i] = (f1 - f2) / (2.0 * eps);
58	x1[i] = x[i]; x2[i] = x[i];
59	}
60	return g;
61	}
62
63	int main() {
64	std::mt19937 rng(123);
65	std::uniform_real_distribution<double> dist(-1.0, 1.0);
66
67	size_t n = 5;
68	std::vector<double> A(n*n), x(n), b(n);
69	for (double &v : A) v = dist(rng);
70	for (double &v : x) v = dist(rng);
71	for (double &v : b) v = dist(rng);
72
73	// Compute values
74	double fx = f_quadratic(A, x, b);
75	auto gA = grad_analytic(A, x, b);
76	auto gN = grad_numeric(A, x, b);
77	auto H = hessian(A);
78
79	// Report max absolute gradient error
80	double max_err = 0.0;
81	for (size_t i = 0; i < n; ++i) max_err = std::max(max_err, std::abs(gA[i] - gN[i]));
82
83	std::cout << std::fixed << std::setprecision(6);
84	std::cout << "f(x) = " << fx << "\n";
85	std::cout << "Max \|grad_analytic - grad_numeric\| = " << max_err << "\n";
86
87	// Sanity: check Hessian symmetry
88	double anti_sym = 0.0;
89	for (size_t i = 0; i < n; ++i)
90	for (size_t j = 0; j < n; ++j)
91	anti_sym = std::max(anti_sym, std::abs(H[in + j] - H[jn + i]));
92	std::cout << "Max \|H - H^T\| = " << anti_sym << "\n";
93
94	return 0;
95	}
96