Multivariable derivatives, matrix calculus, and automatic differentiation

A scalar training objective depends on many parameters, often through repeated computations and batches of records. Its derivative must account for every dependency, distinguish a partial observation from a whole-input approximation, and preserve array shapes. This lesson develops total differentials, Jacobians, Hessians and matrix gradients before implementing the chain rule as automatic differentiation.

Retrieval check: recall Module10 matrix multiplication, transposes and batch shapes, and Module16 local models, chain rules and nondifferentiable boundaries. The NumPy labs use the course’s pinned baseline; the scalar reverse engine needs only standard-library Python.

For a scalar f:Rⁿ→R on an open domain, the partial derivative with respect to coordinate i at a is the limit[f(a+h e_i)−f(a)]/h, where e_i is that coordinate’s unit vector. Other coordinates are held fixed. Partial existence only checks a coordinate line. It does not automatically certify continuity for simultaneous changes or a first-order approximation in every direction. A function’s multivariable domain still matters: a log of a coordinate combination or a zero denominator can exclude nearby paths.

For a vector v, the directional derivative is D_v f(a)=lim_{t→0}[f(a+tv)−f(a)]/t when this finite limit exists. A unit v describes rate per unit distance; a nonunit v also scales the speed along the path. The two-sided definition here allows both signs of t. One-sided directional derivatives at a boundary or for a nonsmooth function are differently qualified objects. The symbol D_v is not a guarantee that all these limits form one linear map in v.

A level set consists of inputs having the same scalar output. Where total differentiability holds, its gradient points in the direction of steepest first-order Euclidean increase and is perpendicular to differentiable level-set paths with nonzero tangent. These claims will follow from the differential and chain rule, rather than from a contour drawing alone. Different feature units change the Euclidean geometry unless a metric is specified; use the gradient under the chosen coordinate and inner-product convention.

Partials can miss discontinuity. Define f(x,y)=xy/(x²+y²) away from the origin and f(0,0)=0. Both coordinate-axis restrictions are zero, giving both partials zero at the origin. Along x=y=t≠0, the output is1/2, so continuity fails. Since differentiability implies continuity, those existing partials cannot establish total differentiability. Lab1 samples these paths to illustrate the already exact counterexample.

Even continuity and every directional derivative can be insufficient. Set g(x,y)=x³/(x²+y²) away from zero, g(0,0)=0. Its magnitude is≤|x|≤‖(x,y)‖, so it is continuous at zero. For nonzero v, the directional quotient is v₁³/(v₁²+v₂²). Axis rates are one and zero, but along v=(1,1) the rate is1/2 rather than the linear candidate one. The map of directions is not linear, so no single total derivative can represent it. Testing many straight lines still does not replace the uniform small-vector remainder requirement below.

A scalar f is differentiable at a when there exists a linear functional L with f(a+h)=f(a)+L(h)+r(h), where r(h)/‖h‖→0 as the whole vector h→0. The limit covers every small-vector path, not just coordinate directions. Every linear functional on finite-dimensional Euclidean space can be written gᵀh: expand h in coordinate basis, let g_i=L(e_i), and use linearity. The coefficients equal existing partials because restricting the remainder to one coordinate line recovers their limits.

We use the column gradient ∇f(a) of shape n×1, so the differential is df=∇f(a)ᵀdx, a scalar. The row derivative mapping inputs to a scalar is its transpose, shape1×n. A NumPy flat array of shape(n,) can represent those coefficients in code but is neither an explicit row nor column. Annotate the mathematical convention and each operation’s stored shape; otherwise a transpose error can hide when square examples happen to fit.

For a vector F:Rⁿ→Rᵐ, differentiability means F(a+h)=F(a)+Jh+r(h), with ‖r(h)‖/‖h‖→0. The Jacobian J has shape m×n: row i contains derivatives of output i with respect to all inputs. Its entries are J_{ij}=∂F_i/∂x_j. Thus dF=J dx. Componentwise scalar differentiability is equivalent in finite dimensions because a finite collection of vanishing remainders also vanishes in vector norm. The Jacobian is a linear map between input and output spaces, rather than a coordinate list with an arbitrary orientation.

Continuous partial derivatives in a neighbourhood supply a sufficient differentiability condition. For a scalar function, change coordinates one at a time and telescope the total increment. Apply the one-variable mean value theorem to each change. The resulting coefficients are coordinate partials at intermediate points approaching a. Subtract the partials at a; continuity makes every coefficient error≤ε for all sufficiently small h. The total remainder is at most εΣ|h_i|≤ε√n‖h‖, so its ratio vanishes. This proves the C¹ sufficient condition. It is not a necessary condition; a function can be differentiable at one point without having a continuous derivative throughout a neighbourhood.

Total differentiability implies every directional derivative equals ∇fᵀv and implies continuity: a bounded linear map sends h→0 to zero and the remainder does too. For unit v, Cauchy–Schwarz gives ∇fᵀv≤‖∇f‖, with equality in the gradient direction if the gradient is nonzero. If it is zero, every first-order unit-direction rate is zero; higher-order responses may remain. Along a differentiable level-set path c(t), the chain rule yields ∇f(c(t))ᵀc′(t)=0, proving the normal interpretation under the specified assumptions.

The differential is a mathematical first-order mapping, not a literal infinitesimal stored by Python. It can be applied to a specified finite perturbation as an approximation, with a remainder whose size must be controlled separately. Scaling coordinates or changing the metric changes which gradient vector represents the same differential. In a weighted inner product with positive-definite M, the representing gradient is M⁻¹ times the Euclidean gradient. The differential itself retains the same scalar action; geometric steepest directions depend on the chosen notion of length.

An explicit example separates differentiability at a point from C¹ regularity nearby. Define s(x)=x²sin(1/x) for nonzero x and s(0)=0. Its quotient at zero is x sin(1/x), bounded in magnitude by|x|, so s′(0)=0. For nonzero x, the rules give s′(x)=2x sin(1/x)−cos(1/x), which oscillates near zero. The function f(x,y)=s(x)+y² is differentiable at the origin because its magnitude is at most x²+y²=‖(x,y)‖², giving a zero linear map and a vanishing remainder ratio. Its derivative is not continuous nearby. Thus the neighbourhood C¹ condition is useful and sufficient, while demanding it as the definition would exclude valid pointwise differentiability.

For differentiable F:Rⁿ→Rᵐ and G:Rᵐ→Rᵖ with compatible local domains, the chain rule is J_{G∘F}(a)=J_G(F(a))J_F(a), shape(p×m)(m×n)=p×n. The inner map acts first on input increments, so its Jacobian is on the right. This order matches the matrix composition rule already established in Module10. A reversed product may be undefined or describe a different map even when square shapes permit it.

The proof uses local remainders. Write ΔF=J_Fh+r_F(h), with ‖r_F‖/‖h‖→0. In particular ‖ΔF‖/‖h‖ is bounded near zero. The outer increment is J_GΔF+r_G(ΔF). After substitution, J_G r_F is negligible relative to‖h‖; the second remainder is bounded by[‖r_G(ΔF)‖/‖ΔF‖][‖ΔF‖/‖h‖], tending to zero when ΔF≠0, and is zero directly when ΔF=0. Thus the product is the total derivative. This is the same safeguard against dividing by a vanished inner increment used for the scalar chain rule.

A Jacobian–vector product or JVP computes Jv for an input direction v, returning an output perturbation direction of length m. A vector–Jacobian product or VJP is traditionally written qᵀJ for an output row seed qᵀ. In our column convention its transpose is Jᵀq, an input adjoint vector of length n. These two equivalent notations must not be confused with multiplying J by an output seed. The seed specifies which output combination is being differentiated; a vector output does not have one unique gradient without that choice.

A rectangular example exposes shape mistakes that a square Jacobian can hide. Take F(x,y)=(xy,x²,y²) and L(u,v,w)=u+2v−3w. At(1,2), J_F has rows(2,1),(2,0),(0,4), so its shape is3×2. The output seed is q=(1,2,−3), length three. The valid pullback J_Fᵀq is(6,−11), length two, agreeing with derivatives of xy+2x²−3y². Trying J_Fq has incompatible inner dimensions. A direction v of length two instead produces three output rates through J_Fv. Writing shapes before numerical values makes the distinction visible without relying on an accidental dimension error later in a program.

Forward mode propagates an input tangent through each local JVP. Reverse mode propagates an output adjoint through each local transpose-Jacobian action. For a scalar output with many inputs, one reverse pass returns all input coefficients for seed one. Obtaining a full m×n Jacobian can require n input-basis forward passes or m output-basis reverse passes; neither mode automatically produces every entry in the cost of one scalar path. Special structure can alter practical costs, but the shape and seed contract remains.

For L(y) scalar and y=Ax, dy=A dx and dL=qᵀA dx, so the pullback is Aᵀq. For a composed nonlinear map, evaluate local Jacobians at the actual forward values before reversing them. Reusing a derivative from a different input or changing stored forward values before backward propagation can invalidate the result. Automatic differentiation systems store or recompute the required intermediate information; the chain rule explains why it is needed.

Shared inputs require sums of contributions. If a variable feeds two outputs or two parents, its differential influences each path, and the final scalar coefficient adds those influences. Reverse traversal must wait until all downstream contributions have reached a shared node before applying its local backward rule. A unique-node topological order ensures this schedule; treating a graph as a tree by recursively backpropagating without accumulation can lose paths or propagate a shared total multiple times.

For scalar f, the Hessian is the Jacobian of its column gradient, shape n×n. Entry H_{ij} is the derivative of gradient component i with respect to coordinate j. If second partials are continuous in a neighbourhood, mixed derivatives commute and H is symmetric. One justification considers the four-corner increment over a small coordinate rectangle and applies the mean value theorem in both orders. Dividing by the rectangle’s area and taking both widths to zero gives equal limits by mixed-partial continuity. State that hypothesis rather than assume every twice-written partial must agree.

Without continuity, mixed partials can differ. For f(x,y)=xy(x²−y²)/(x²+y²) away from zero and f(0,0)=0, f_x(0,y)=−y while f_y(x,0)=x. Hence at zero the derivative of f_x with respect to y is−1, while the derivative of f_y with respect to x is1. The smooth symmetry theorem’s neighbourhood hypothesis fails. This is not a contradiction; it identifies the boundary of a theorem that automatic symbolic rearrangement might hide.

For C² f, the second-order local model is f(a+h)=f(a)+∇f(a)ᵀh+(1/2)hᵀH(a)h+o(‖h‖²). To derive it, restrict to the scalar path g(t)=f(a+th). Chain rules give g′(t)=∇f(a+th)ᵀh and g″(t)=hᵀH(a+th)h. One-variable Taylor between zero and one provides an intermediate Hessian along that segment. Continuity of H makes its difference from H(a) vanish uniformly as h→0, so the quadratic remainder divided by‖h‖² tends to zero. This also shows why one fixed coordinate-axis expansion is weaker than the multivariable model.

Listing existing second coordinate partials is weaker than proving the gradient itself has a total derivative. The unequal-mixed-partial example above supplies such entries but does not supply the smooth second-order map required by this Taylor argument. Its failure should therefore be described as a failure of the hypotheses, rather than inserting a nonsymmetric entry table into a C² theorem. In numerical work, a nearly nonsymmetric estimate can also arise from finite differences or implementation errors even for a smooth objective. Identify the mathematical regularity first, then assess the estimator and array convention separately.

At an interior differentiable stationary point the gradient must be zero, by applying scalar Fermat along every coordinate. With C² regularity and positive-definite Hessian, the quadratic term dominates the smaller remainder and gives a strict local minimum; negative definite gives a maximum. An indefinite Hessian supplies positive and negative curvature directions and therefore a saddle. A positive-semidefinite Hessian with zero directions is inconclusive without more information. At zero, x⁴+y⁴ has a strict minimum, while x⁴−y⁴ is a saddle; both have zero Hessian. These are local classifications, with global optimisation developed next module.

For a symmetric Hessian, the directional second derivative along v is vᵀHv. Hessian–vector multiplication Hv also describes first-order change of the gradient along v. Computing such products can avoid storing all n² entries, which matters for large models. It still differentiates a gradient under sufficient smoothness; at a kink the classical Hessian may not exist. A finite difference of gradients is a separate estimator with its own step and rounding errors.

The positive-definite local-minimum test requires a uniform quadratic margin, not only positive values in a few sampled directions. The minimum of vᵀHv over the unit sphere is some λ>0 by compactness and positivity. Thus hᵀHh≥λ‖h‖² for all h. The little-o remainder can be bounded by λ‖h‖²/4 in a sufficiently small neighbourhood, leaving the stationary function increment at least λ‖h‖²/4 for every nonzero h there. This completes the domination argument behind the strict local test and explains why zero curvature directions prevent this certificate.

For vectors, collect the coefficient of dx in df=gᵀdx. For a matrix A, use the Frobenius inner product: df=Σ_{ij}(∇A f){ij}dA_{ij}=tr((∇A f)ᵀdA). The gradient has the same shape as A. Here tr sums a square matrix’s diagonal. The identity tr(BC)=tr(CB) for compatible rectangular factors follows by writing both asΣ{ij}B_{ij}C_{ji}; cyclic rearrangements can then collect a matrix differential while preserving valid dimensions.

For fixed b, f(x)=bᵀx gives df=bᵀdx and gradient b. For a general square A, f(x)=xᵀAx gives df=dxᵀAx+xᵀA dx=[(A+Aᵀ)x]ᵀdx. Thus the gradient is(A+Aᵀ)x, reducing to2Ax only when A is symmetric. The antisymmetric part contributes zero to the quadratic scalar and therefore no gradient. A formula2Ax for an arbitrary nonsymmetric A is not repaired merely because its output shape fits.

Loss reductions are mathematical definitions. For L=(1/B)Σr_i², gradient is(2/B)Xᵀr. If outputs have o columns and “mean” averages all B·o scalar entries, denominator is B·o rather than B. If it instead averages rows after summing output errors, denominator is B. These objectives differ by a factor o. Report the exact reduction; silently changing it alters derivatives and the effective scale of subsequent updates.

For matrix predictions XW with X:B×d,W:d×o,Y:B×o, write R=XW−Y. The half-squared Frobenius objective L=(1/2)‖R‖_F² has differential tr(RᵀdR), with dR=X dW+dX W−dY. Collecting coefficients gives ∇_W L=XᵀR of shape d×o, ∇_X L=R Wᵀ of shape B×d, and ∇_Y L=−R of shape B×o. Verify each by a differential or indexed sum before implementing it. The residual and gradient arrays need explicit output columns to prevent accidental broadcasting with a flat target.

Adding a shared bias b with mathematical shape o×1 gives prediction XW+1bᵀ. Its differential repeats db across B rows, so the bias gradient is Rᵀ1 for the half-squared objective, equivalently a sum of output adjoints over the batch. In NumPy b often has shape(o,), with backward sum(axis=0) returning that same stored shape. Broadcasting in the forward operation corresponds to summation in the backward operation because one parameter influences every replicated entry. Overwriting the bias with the last row’s contribution drops the others.

These derivations also apply to a two-layer network. With Z=XW₁+1b₁ᵀ,H=tanh(Z),P=HW₂+1b₂ᵀ, and output adjoint G=∂L/∂P, gradients are HᵀG for W₂ and row sums of G for b₂. Hidden adjoint is G W₂ᵀ; multiply elementwise by1−H² for Z’s adjoint D. Then W₁ gradient is XᵀD and b₁ gradient sums D over rows. Matrix products and elementwise products serve different roles and cannot be exchanged. The lab annotates every shape and verifies every parameter entry numerically at a smooth point.

For the lab’s batch of three records with two features, three hidden units and one output, W₁ is2×3, Z and H are3×3, W₂ is3×1, and P,G are3×1. The hidden pullback G W₂ᵀ is3×3; XᵀD is2×3, matching W₁. The stored biases have shapes(3,) and(1,), recovered by row sums. These dimensions deliberately make a missing output-weight transpose fail instead of silently multiplying square arrays. The loss is a mean over three scalar residual entries, so G=2R/3. Checking the stated loss before backward propagation is as necessary as checking the transposes: a correctly shaped gradient can still have the wrong reduction factor.

Automatic differentiation applies local derivative rules to the executed computation graph. It is distinct from finite differences, which perturb inputs and divide output differences, and from symbolic differentiation, which constructs algebraic expressions. Forward mode carries values and tangents; reverse mode first records forward values and then propagates adjoints in reverse dependency order. Both evaluate chain-rule derivatives of the selected mathematical primitive operations, with ordinary numerical rounding in those evaluations.

The literal floating-point input-output map includes discrete rounding behaviour; AD normally differentiates the intended real primitive formulas evaluated at stored values, rather than the staircase map of every rounding instruction. Thus “automatic” does not imply exact real arithmetic, eliminate conditioning, or establish differentiability across a branch. If a program selects a branch, the returned derivative describes that executed path under the primitive conventions. Thresholds and kinks need separate interpretation; a chosen ReLU backward value at zero remains a convention.

The lab engine stores parents and a local backward closure for each node. It first creates a topological list containing every reachable node once. Before a fresh backward pass it resets all adjoints, seeds the output, and traverses the reversed list. Each closure uses += for every parent edge, including when both operands are the same object. This ensures all downstream contributions are present before a shared node propagates its accumulated adjoint. Repeated passes with no reset would instead accumulate old and new requests, which is useful only when deliberately specified.

A seed other than one computes a scalar VJP scaled by that seed. At x=2, seed two returns74 and seed zero returns zero, while the ordinary derivative remains37. For multiple outputs, a vector seed specifies a weighted output combination. Gradients are therefore answers to a defined differentiated quantity, not context-free attributes permanently attached to a variable. Model state, graph reuse and output reduction must be included in that specification.

Gradient checks compare an analytic or AD result with central differences at selected smooth points and a range of meaningful step sizes. Use absolute and relative tolerances suitable for magnitude, and inspect shapes before comparing entries. An erroneous formula can agree at a symmetric or zero-residual point, so choose nontrivial values and perturb each parameter or independent directions. The network lab checks every small parameter entry, then explicitly exhibits a transpose fault, lost batch contribution and missing mean factor. Passing these checks supports that implementation on those examples; the differential derivation supplies the general formula.

Reverse mode needs stored or recomputed intermediate values and a sound graph schedule. Its practical memory and execution cost depend on the computation, not only the parameter count. Large systems can use checkpointing or Hessian–vector products, but those engineering choices retain the same chain and accumulation contracts. This foundation is sufficient to assess a manual backward pass and interpret what an automatic one is calculating before studying optimisation.

Lab1 · Vector derivatives and a partials counterexample

Download lab1_vector_derivatives.py

"""Finite differences check stated analytic vector derivatives at smooth points."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)

def F(z):
    x, y = z
    return np.array([x*x*y, np.sin(x)+y*y])

point = np.array([0.7, -0.4])
x, y = point
J = np.array([[2*x*y, x*x], [np.cos(x), 2*y]])
h = 1e-5
basis = np.eye(2)
numerical_J = np.column_stack([(F(point+h*e)-F(point-h*e))/(2*h) for e in basis])
print("Point / function value:", point, F(point))
print("Analytic Jacobian (outputs by inputs):\n", J)
print("Finite-difference Jacobian:\n", numerical_J)
print("Jacobian maximum error:", np.max(np.abs(J-numerical_J)))
assert J.shape == (2, 2) and np.allclose(J, numerical_J, atol=1e-9, rtol=0)

direction = np.array([0.6, 0.8])
directional = (F(point+h*direction)-F(point-h*direction))/(2*h)
print("JVP / directional difference:", J@direction, directional)
assert np.allclose(J@direction, directional, atol=1e-9, rtol=0)

def loss(z):
    return z[0]**2*z[1]+np.exp(z[1])
def grad(z):
    return np.array([2*z[0]*z[1], z[0]**2+np.exp(z[1])])
H = np.array([[2*y, 2*x], [2*x, np.exp(y)]])
numerical_grad = np.array([(loss(point+h*e)-loss(point-h*e))/(2*h) for e in basis])
numerical_H = np.column_stack([(grad(point+h*e)-grad(point-h*e))/(2*h) for e in basis])
print("Scalar gradient:", grad(point))
print("Hessian:\n", H)
print("Gradient / Hessian maximum errors:", np.max(np.abs(grad(point)-numerical_grad)), np.max(np.abs(H-numerical_H)))
assert np.allclose(grad(point), numerical_grad, atol=1e-9, rtol=0)
assert np.allclose(H, numerical_H, atol=1e-9, rtol=0)

# Axis partials alone do not imply continuity, hence do not imply differentiability.
def bad(z):
    a, b = z
    return a*b/(a*a+b*b) if a*a+b*b else 0.0
print("Origin axis samples / diagonal samples:")
for t in [1e-1, 1e-3, 1e-6]:
    print(t, bad([t, 0]), bad([0, t]), bad([t, t]))
    assert bad([t, 0]) == 0 and bad([0, t]) == 0 and bad([t, t]) == 0.5
print("Smooth derivative checks illustrate the analytic claims; axis samples cannot certify total differentiability.")

Point / function value: [ 0.7 -0.4] [-0.196     0.804218]
Analytic Jacobian (outputs by inputs):
 [[-0.56      0.49    ]
 [ 0.764842 -0.8     ]]
Finite-difference Jacobian:
 [[-0.56      0.49    ]
 [ 0.764842 -0.8     ]]
Jacobian maximum error: 1.702149532434305e-11
JVP / directional difference: [ 0.056    -0.181095] [ 0.056    -0.181095]
Scalar gradient: [-0.56     1.16032]
Hessian:
 [[-0.8      1.4    ]
 [ 1.4      0.67032]]
Gradient / Hessian maximum errors: 1.24074084340009e-11 1.5253354135325026e-11
Origin axis samples / diagonal samples:
0.1 0.0 0.0 0.5
0.001 0.0 0.0 0.5
1e-06 0.0 0.0 0.5
Smooth derivative checks illustrate the analytic claims; axis samples cannot certify total differentiability.

Derive J,∇f,H before running, then compare their shapes and numerical errors. Explain the JVP’s output dimension and the smoothness behind each check. The final axis/diagonal values illustrate a mathematical counterexample; they are not a sampling-based proof of continuity.

Lab2 · A minimal reverse-mode engine

Download lab2_reverse_engine.py

"""A tiny scalar reverse engine: unique-node order, reset, seed, and += adjoints."""
import math

class Value:
    def __init__(self, data, parents=(), label=""):
        self.data, self.grad = float(data), 0.0
        self.parents, self.label = parents, label
        self._backward = lambda: None

    def __add__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data+other.data, (self, other), "+")
        def backward():
            self.grad += out.grad
            other.grad += out.grad
        out._backward = backward
        return out
    __radd__ = __add__

    def __mul__(self, other):
        other = other if isinstance(other, Value) else Value(other)
        out = Value(self.data*other.data, (self, other), "*")
        def backward():
            self.grad += other.data*out.grad
            other.grad += self.data*out.grad
        out._backward = backward
        return out
    __rmul__ = __mul__

    def tanh(self):
        out = Value(math.tanh(self.data), (self,), "tanh")
        def backward():
            self.grad += (1-out.data*out.data)*out.grad
        out._backward = backward
        return out

    def backward(self, seed=1.0):
        visited, order = set(), []
        def visit(node):
            if node in visited:
                return
            visited.add(node)
            for parent in node.parents:
                visit(parent)
            order.append(node)
        visit(self)
        for node in order:
            node.grad = 0.0
        self.grad = float(seed)
        for node in reversed(order):
            node._backward()
        return order

x = Value(2, label="x")
u = x*x
u.label = "u=x*x"
square = u*u
square.label = "u*u"
z = square+u+x
z.label = "z"
order = z.backward()
print("Shared graph z=(x*x)^2+x*x+x at x=2")
for node in order:
    print(node.label, "value", node.data, "adjoint", node.grad)
assert z.data == 22 and x.grad == 37 and u.grad == 9
z.backward(seed=2)
print("Seed two input adjoint:", x.grad)
assert x.grad == 74
z.backward()
print("Fresh seeded pass resets rather than doubles:", x.grad)
assert x.grad == 37

a, b = Value(2), Value(-1)
mixed = a*b+a*a
mixed.backward()
print("Two-input value / gradient:", mixed.data, a.grad, b.grad)
assert mixed.data == 2 and (a.grad, b.grad) == (3, 2)

v = Value(0.4)
nonlinear = (v*v+0.5*v).tanh()
nonlinear.backward()
def reference(t):
    return math.tanh(t*t+0.5*t)
h = 1e-6
finite = (reference(0.4+h)-reference(0.4-h))/(2*h)
print("Nonlinear AD / finite difference:", v.grad, finite)
assert abs(v.grad-finite) < 1e-8
print("Overwriting a shared adjoint would discard a path; addition is part of the chain rule.")

Shared graph z=(x*x)^2+x*x+x at x=2
x value 2.0 adjoint 37.0
u=x*x value 4.0 adjoint 9.0
u*u value 16.0 adjoint 1.0
+ value 20.0 adjoint 1.0
z value 22.0 adjoint 1.0
Seed two input adjoint: 74.0
Fresh seeded pass resets rather than doubles: 37.0
Two-input value / gradient: 2.0 3.0 2.0
Nonlinear AD / finite difference: 1.1450754518266433 1.145075451791655
Overwriting a shared adjoint would discard a path; addition is part of the chain rule.

Trace the unique-node topological order and distinguish it from parent-edge contributions. Explain why x*x adds two contributions even though both parents reference the same object. Predict seed-two and reset results, then compare the nonlinear example with a finite difference at a smooth point.

Lab3 · Every two-layer parameter gradient

Download lab3_network_backward.py

"""Check every two-layer tanh-network parameter and expose shape/scaling faults."""
import numpy as np
np.set_printoptions(precision=6, suppress=True)
X = np.array([[0.2,-0.3], [0.5,0.7], [-0.4,0.1]])
Y = np.array([[0.1], [0.8], [-0.2]])
parameters = {
    "W1": np.array([[0.3,-0.2,0.1], [0.4,0.2,-0.5]]),
    "b1": np.array([0.1,-0.1,0.05]),
    "W2": np.array([[0.5], [-0.3], [0.2]]),
    "b2": np.array([0.02]),
}

def forward(p):
    H = np.tanh(X@p["W1"]+p["b1"])
    prediction = H@p["W2"]+p["b2"]
    residual = prediction-Y
    return np.mean(residual*residual), H, prediction, residual

loss, H, prediction, residual = forward(parameters)
# Mean of B*o entries, here B=3 and o=1; the exact denominator is explicit.
G = 2*residual/residual.size
D = (G@parameters["W2"].T)*(1-H*H)
gradients = {"W2": H.T@G, "b2": G.sum(axis=0), "W1": X.T@D, "b1": D.sum(axis=0)}
print("Loss / prediction:", loss, prediction.ravel())
print("Batch / input / hidden / output shapes:", X.shape, parameters["W1"].shape, H.shape, prediction.shape)
h = 1e-6
for name, array in parameters.items():
    numerical = np.zeros_like(array)
    for index in np.ndindex(array.shape):
        saved = array[index]
        array[index] = saved+h
        plus = forward(parameters)[0]
        array[index] = saved-h
        minus = forward(parameters)[0]
        array[index] = saved
        numerical[index] = (plus-minus)/(2*h)
    error = np.max(np.abs(numerical-gradients[name]))
    print(name, "shape", gradients[name].shape, "maximum derivative error", error)
    assert gradients[name].shape == array.shape and error < 1e-9

try:
    G@parameters["W2"]  # Missing transpose: (3,1) cannot multiply (3,1).
except ValueError:
    print("Missing transpose rejected by incompatible derivative shapes.")
wrong_bias = D[-1]  # Overwrite rather than accumulate contributions from all rows.
wrong_sum_scaling = 2*residual  # This differentiates a sum, not the stated mean.
print("Bias overwrite maximum error:", np.max(np.abs(wrong_bias-gradients["b1"])))
print("Missing mean factor scales output adjoint by:", residual.size)
assert not np.allclose(wrong_bias, gradients["b1"])
assert np.allclose(wrong_sum_scaling, residual.size*G)
print("Shape, shared contributions, and loss reduction are separate parts of the backward contract.")

Loss / prediction: 0.12469122832324815 [0.142506 0.231685 0.022003]
Batch / input / hidden / output shapes: (3, 2) (2, 3) (3, 3) (3, 1)
W1 shape (2, 3) maximum derivative error 2.7340560371236222e-11
b1 shape (3,) maximum derivative error 8.312184274217316e-12
W2 shape (3, 1) maximum derivative error 4.48385772955362e-12
b2 shape (1,) maximum derivative error 5.846850781310309e-12
Missing transpose rejected by incompatible derivative shapes.
Bias overwrite maximum error: 0.13066171364054682
Missing mean factor scales output adjoint by: 3
Shape, shared contributions, and loss reduction are separate parts of the backward contract.

Annotate forward and backward shapes and write the exact mean loss. Check each parameter entry, inspect the rejected transpose, and explain why overwriting a shared bias or omitting the mean factor violates the contract even when resulting arrays have plausible shapes.

Exercises1–12 are required. Optional Exercises13–14 add35 minutes beyond the twelve-hour core schedule.

Use gradient, Jacobian, Hessian and chain-rule selections in the authors’ Mathematics for Machine Learning companion and the primary MIT OpenCourseWare multivariable calculus course. The differential derivations and small reverse engine above provide the lesson’s own concrete shape and accumulation conventions.

Give a partial-existence counterexample and a C¹ differentiability proof. State column-gradient and output-by-input Jacobian conventions, derive a chain pullback and a second-order model, then recover quadratic and least-squares gradients through differentials.

Exit task: For X:B×d,W:d×o and half-squared residualR=XW−Y,deriveXᵀR,RWᵀ and the row-summed bias contribution. Extend through a tanh hidden layer, specifying elementwise factors and the exact mean reduction. Explain how a unique-node reverse schedule still accumulates both edges of x*x.

Ready to move on: you can derive and verify local vector/matrix sensitivities and interpret an automatic backward pass. Optimisation next studies what those derivatives do and do not imply about minimisers and update convergence. See the course overview.