Math for Data Science · MDS-4
The math behind gradients, Jacobians, and backpropagation — built for machine learning.
Multivariable and matrix calculus is the engine room of modern machine learning: every gradient step, every backpropagation pass, every optimizer rests on it. This course builds the geometry of functions of several variables, then the partial derivatives, gradients, Jacobians, and Hessians that power them, and finishes with the chain rule and matrix-calculus identities used to derive gradients for real models. You will verify every formula numerically in NumPy/autograd so the math is never a black box.
Course Outline
Every lesson is self-contained: a full explanation, worked step-by-step examples (in Python and R/Julia), common mistakes, a try-it-yourself, and an interactive quiz. Jump to a unit:
A scalar field takes a point in n-dimensional space and returns a single number — the foundation of every loss function in ML.
A scalar field is a function f: R^n -> R that assigns one real number to each input point. In machine learning the loss L(theta) is exactly such a field: it maps a parameter vector to a single scalar error we want to minimize. We evaluate it by substituting every coordinate at once, e.g. f(x,y) = x^2 + 3y^2 at (1,2) gives 1 + 12 = 13. Scalar fields have no inherent direction in their output, but their rate of change has direction — captured later by the gradient. Understanding them as height surfaces over an input plane is the mental model that makes contour plots, gradients, and gradient descent intuitive. Every optimization problem is a search for the lowest point of a scalar field.
Worked Example 1
Problem. Evaluate f(x,y) = x^2 + 3*y^2 at the point (2, -1).
Answer. f(2,-1) = 7
Worked Example 2
Problem. For f(x,y,z) = x*y + z^2, compute f(1, 4, -2).
Answer. f(1,4,-2) = 8
Worked Example 3
Problem. A squared-error loss L(w1,w2) = (w1 - 3)^2 + (w2 + 1)^2 is a scalar field on R^2. Find its minimum value and where it occurs.
Answer. Minimum value 0 at (w1,w2) = (3,-1)
Problem. Evaluate f(x,y) = 2*x^2 - x*y + y^2 at (1, 2), then state whether the output is a scalar or a vector.
Solution. 2*x^2 = 2*1 = 2; -x*y = -1*2 = -2; y^2 = 4. Sum: 2 - 2 + 4 = 4. The output 4 is a single scalar, so f is a scalar field f: R^2 -> R. Answer: 4 (scalar).
A level set collects all inputs giving the same output height; drawing several of them yields a contour map of the surface.
The graph of f: R^2 -> R is a surface z = f(x,y). A level set is {(x,y) : f(x,y) = c} — every input that produces the same height c. In 2D these are contour lines, exactly like a topographic map. A contour plot draws several level sets at evenly spaced c values. Two facts make this the workhorse of optimization intuition: contour spacing encodes steepness (tightly packed lines mean a steep region, a large gradient), and the gradient is always perpendicular to the contour through a point. For a convex bowl like x^2 + 3y^2 the contours are nested ellipses shrinking toward the minimum; for a saddle x^2 - y^2 they form hyperbolas. Reading a loss landscape's contours tells you where gradient descent will move and how fast.
Worked Example 1
Problem. Describe the level set f(x,y) = x^2 + y^2 = 4.
Answer. A circle of radius 2 centered at the origin
Worked Example 2
Problem. For f(x,y) = x^2 + 3*y^2, identify the shape of the level sets and which axis they stretch along.
Answer. Nested ellipses, elongated along the x-axis
Worked Example 3
Problem. For the saddle g(x,y) = x^2 - y^2, what curve is the level set g = 0?
Answer. Two crossing lines y = x and y = -x
Problem. Sketch (in words) the level sets of f(x,y) = y - x^2 for c = -1, 0, 1, and identify the family of curves.
Solution. Set y - x^2 = c, i.e. y = x^2 + c. Each level set is an upward parabola shifted vertically: c=-1 gives y=x^2-1, c=0 gives y=x^2, c=1 gives y=x^2+1. Answer: a family of vertically stacked parabolas y = x^2 + c.
A multivariable limit exists only if the function approaches the same value along every path to the point — far stricter than the 1D two-sided rule.
For f: R^n -> R, the limit as a point approaches a equals L only if f tends to L along every possible approach path. In one dimension we check left and right; in higher dimensions there are infinitely many directions and curved paths, so finding two paths with different limits proves the limit does NOT exist. f is continuous at a if the limit there equals f(a). This matters for ML because activation functions and losses must be continuous (and usually differentiable) for gradients to be meaningful; a discontinuity breaks gradient descent. A classic failure is f(x,y) = xy/(x^2+y^2) at the origin: along y=0 it is 0, but along y=x it is 1/2, so no limit exists and the function cannot be made continuous there.
Worked Example 1
Problem. Show that lim of f(x,y) = (x^2 - y^2)/(x^2 + y^2) at (0,0) does not exist.
Answer. Limit does not exist (path-dependent: 1 vs -1)
Worked Example 2
Problem. Evaluate lim of f(x,y) = (x^2*y)/(x^2 + y^2) at (0,0).
Answer. Limit = 0 (exists)
Worked Example 3
Problem. Is f(x,y) = (xy)/(x^2 + y^2) (with f(0,0)=0) continuous at the origin?
Answer. Not continuous at the origin
Problem. Determine whether lim of f(x,y) = (x*y^2)/(x^2 + y^4) exists at (0,0).
Solution. Try straight lines y = m*x: f = (x*m^2*x^2)/(x^2 + m^4*x^4) = (m^2*x^3)/(x^2(1+m^4 x^2)) -> 0. But try the parabola x = y^2: f = (y^2*y^2)/(y^4 + y^4) = y^4/(2y^4) = 1/2. Lines give 0 but the parabola gives 1/2, so the limit does not exist.
A vector-valued function maps one parameter to a point in R^n, tracing a curve whose velocity is its component-wise derivative.
A vector-valued function r: R -> R^n sends a single parameter t to a point r(t) = (x(t), y(t), ...), tracing a curve through space. Its derivative r'(t) is taken component by component and gives the velocity (tangent) vector pointing along the direction of travel; its magnitude is the speed. These curves are how we describe paths — for example the trajectory of parameters during training, or a line search direction in optimization. In ML, composing a scalar loss with a parametric path, L(r(t)), and differentiating via the chain rule is exactly how directional derivatives and line searches are computed. Mastering component-wise differentiation now sets up the Jacobian (the multi-output generalization) in Unit 3.
Worked Example 1
Problem. For r(t) = (cos t, sin t), find the velocity r'(t) and its speed.
Answer. r'(t) = (-sin t, cos t), constant speed 1
Worked Example 2
Problem. For r(t) = (t^2, 3*t, t^3), compute r'(1).
Answer. r'(1) = (2, 3, 3)
Worked Example 3
Problem. A parameter path is r(t) = (1 - t)*(0,0) + t*(4,2) for t in [0,1]. Find r'(t) and interpret it.
Answer. r'(t) = (4, 2); a straight line in the direction (4,2)
Problem. For r(t) = (e^t, t^2, ln(t)) with t > 0, find r'(t) and evaluate at t = 1.
Solution. Differentiate componentwise: d/dt e^t = e^t; d/dt t^2 = 2t; d/dt ln t = 1/t. So r'(t) = (e^t, 2t, 1/t). At t=1: (e^1, 2, 1) = (e, 2, 1) approx (2.718, 2, 1).
Open sets and neighborhoods formalize 'wiggle room' around a point — the setting in which derivatives and limits are well defined.
A neighborhood of a point a is a small open ball {x : ||x - a|| < r} of all points within distance r. A set is open if every one of its points has such a ball entirely inside the set — no boundary included. This matters because partial derivatives, gradients, and limits require approaching a point from all directions, which is only possible at interior points of an open domain. The domain of f is the set of inputs where f is defined; for f = ln(x*y) the domain excludes points where x*y <= 0. In optimization, working on an open domain guarantees we can take a small gradient step in any direction without leaving the feasible region, which is why unconstrained problems assume open domains.
Worked Example 1
Problem. Find the domain of f(x,y) = sqrt(9 - x^2 - y^2).
Answer. Domain = {(x,y) : x^2 + y^2 <= 9}, the disk of radius 3
Worked Example 2
Problem. Find the domain of f(x,y) = ln(x - y).
Answer. Domain = {(x,y) : x > y}, an open half-plane
Worked Example 3
Problem. Is the set {(x,y) : x^2 + y^2 < 1} open? Is {(x,y) : x^2 + y^2 <= 1} open?
Answer. The open disk is open; the closed disk is not open
Problem. Find the domain of f(x,y) = 1 / (x^2 + y^2 - 4) and state whether it is open.
Solution. The denominator must be nonzero: x^2 + y^2 - 4 != 0, i.e. x^2 + y^2 != 4. The domain is all of R^2 except the circle of radius 2. Removing a circle (a closed curve) from the open plane leaves an open set, since every remaining point has a small ball avoiding the circle. Answer: domain = R^2 minus the circle x^2+y^2=4; it is open.
A quadratic form x^T A x encodes a bowl, dome, or saddle; the sign pattern of A's eigenvalues decides which.
A quadratic form is Q(x) = x^T A x with A a symmetric matrix; it is the degree-two heart of every local approximation of a loss (the Hessian term in a Taylor expansion). Its geometry is dictated entirely by the eigenvalues of A. If all eigenvalues are positive, A is positive definite and Q is a bowl with a unique minimum at the origin — the ideal, convex case for optimization. All negative gives a dome (maximum). Mixed signs give a saddle: up in some directions, down in others, where gradient descent stalls. Diagonal A makes this transparent: Q = lambda1 x1^2 + lambda2 x2^2 + .... Recognizing the bowl/saddle structure of the quadratic approximation tells you whether a critical point of a loss is a true minimum.
Worked Example 1
Problem. Write Q(x,y) = 2*x^2 + 3*y^2 in the form x^T A x and identify the matrix A.
Answer. A = [[2,0],[0,3]]; eigenvalues 2,3 > 0 so a bowl
Worked Example 2
Problem. For Q(x,y) = x^2 + 4*x*y + y^2, find A and classify the surface.
Answer. A = [[1,2],[2,1]], eigenvalues 3 and -1: a saddle
Worked Example 3
Problem. Classify Q(x,y) = -x^2 - 2*y^2 and evaluate it at (1,1).
Answer. Negative definite (dome/maximum); Q(1,1) = -3
Problem. Express Q(x,y) = 3*x^2 - 2*x*y + 3*y^2 as x^T A x, find the eigenvalues of A, and classify the surface.
Solution. a11=3, a22=3, cross term -2xy means a12 = -1, so A = [[3,-1],[-1,3]]. Eigenvalues: det(A - lambda I) = (3-lambda)^2 - 1 = 0 -> 3-lambda = +/-1 -> lambda = 2, 4. Both positive, so A is positive definite and Q is a bowl with a minimum at the origin.
Pick a two-variable loss-like function such as f(x,y) = x^2 + 3*y^2 and a saddle such as g(x,y) = x^2 - y^2. Use NumPy meshgrid to evaluate each on a grid and Matplotlib to draw filled contour plots side by side. Annotate which function is a bowl (a convex minimum) and which is a saddle, and explain how the spacing of contour lines signals steepness.
Deliverable · A notebook or script producing both contour plots plus a short markdown note classifying each surface and reading off its steepest region.
1. A function f: R^3 -> R takes how many inputs and produces how many outputs?
Answer C. The notation R^3 -> R means three real inputs mapped to a single real scalar output.
2. On a contour plot, what does it mean when level curves are packed very close together?
Answer B. Closely spaced contours mean the height changes a lot over a small distance, indicating a steep region.
3. Which expression is a quadratic form?
Answer C. A quadratic form is x^T A x, a homogeneous degree-two polynomial encoded by a matrix A.
4. For a multivariable limit to exist at a point, the function must approach the same value...
Answer C. Unlike one dimension with two sides, higher dimensions have infinitely many approach paths; all must agree.
5. A vector-valued function r(t) = (cos t, sin t) describes what object?
Answer B. Mapping a single parameter t to a point in R^2 traces a curve; here it is the unit circle.
I can describe a function of several variables and sketch its level sets and contours.
I can reason about limits and continuity for functions of two or more variables.
I can recognize the surface shape implied by a quadratic form.
A partial derivative measures how f changes when you nudge one input while freezing all the others — the atomic operation of every gradient.
The partial derivative of f with respect to x, written df/dx or f_x, is the ordinary derivative of f treating every other variable as a constant. Geometrically it is the slope of the surface along a slice parallel to the x-axis. For f(x,y) = x^2*y, holding y fixed gives f_x = 2*x*y; holding x fixed gives f_y = x^2. These partials are the building blocks of the gradient and therefore of backpropagation: every weight's gradient in a network is ultimately a partial derivative of the loss with respect to that weight, with all other weights momentarily frozen. Mastering the freeze-the-rest mechanic is what lets you differentiate losses with thousands of parameters without confusion. Second partials (f_xx, f_xy) extend this to curvature, previewing the Hessian.
Worked Example 1
Problem. For f(x,y) = x^2*y + 3*y, compute f_x and f_y.
Answer. f_x = 2*x*y, f_y = x^2 + 3
Worked Example 2
Problem. For f(x,y) = sin(x*y), find f_x and f_y at (pi, 1).
Answer. f_x(pi,1) = -1, f_y(pi,1) = -pi
Worked Example 3
Problem. For f(x,y) = e^(x^2 + y), compute the second partials f_xx and f_xy.
Answer. f_xx = (2 + 4*x^2)*e^(x^2+y), f_xy = 2*x*e^(x^2+y)
Problem. For f(x,y) = x/y + y^2, compute f_x and f_y, then evaluate both at (2, 1).
Solution. f_x: y constant, d/dx(x/y) = 1/y, d/dx(y^2) = 0, so f_x = 1/y. f_y: x constant, d/dy(x/y) = -x/y^2, d/dy(y^2) = 2y, so f_y = -x/y^2 + 2y. At (2,1): f_x = 1/1 = 1; f_y = -2/1 + 2 = 0. Answer: f_x = 1, f_y = 0.
Stacking all partial derivatives into one vector gives the gradient — the single object gradient descent follows.
The gradient of a scalar field f: R^n -> R is the column vector of its partial derivatives, grad f = (f_x1, f_x2, ..., f_xn). It is a function of position: at each point it returns a vector. Two facts make it central to ML. First, the gradient points in the direction of steepest increase of f, and its negative is the steepest-decrease direction that gradient descent steps along. Second, its magnitude ||grad f|| is the maximum rate of change at that point. In backpropagation the entire purpose is to compute grad L with respect to every parameter; once we have it, the update theta <- theta - eta*grad L reduces the loss. Always order the components to match the variable order so dot products with directions and updates stay shape-consistent.
Worked Example 1
Problem. Compute the gradient of f(x,y) = x^2 + 3*y^2 and evaluate it at (1, 2).
Answer. grad f(1,2) = (2, 12)
Worked Example 2
Problem. For f(x,y,z) = x*y + y*z + z*x, find grad f at (1, 1, 1).
Answer. grad f(1,1,1) = (2, 2, 2)
Worked Example 3
Problem. At which point is the gradient of f(x,y) = x^2 - 4*x + y^2 + 2*y equal to the zero vector?
Answer. grad f = 0 at (2, -1)
Problem. Find grad f for f(x,y) = ln(x^2 + y^2) and evaluate at (1, 2).
Solution. Let u = x^2 + y^2. f = ln(u). f_x = (1/u)*2x = 2x/(x^2+y^2); f_y = 2y/(x^2+y^2). At (1,2): u = 1+4 = 5, so f_x = 2/5, f_y = 4/5. grad f(1,2) = (2/5, 4/5) = (0.4, 0.8).
The directional derivative measures the rate of change of f along an arbitrary direction, computed as the gradient dotted with a unit vector.
The directional derivative D_u f at a point is the rate at which f changes as you move in the direction of a unit vector u. The key formula is D_u f = grad f . u, the dot product of the gradient with the (normalized) direction. If u is not unit length you must normalize it first, or the rate is scaled wrongly. Because grad f . u = ||grad f|| cos(theta), the rate is largest when u aligns with the gradient (theta = 0) — recovering the steepest-ascent fact — and zero when u is perpendicular to the gradient (moving along a contour). In ML this is how a line search evaluates the slope of the loss along a chosen update direction, and how we confirm the negative gradient is the locally fastest way down.
Worked Example 1
Problem. For f(x,y) = x^2 + y^2, find the directional derivative at (1,1) in the direction v = (3,4).
Answer. D_u f = 14/5 = 2.8
Worked Example 2
Problem. For f(x,y) = x*y, compute the directional derivative at (2,3) toward (1,0).
Answer. D_u f = 3
Worked Example 3
Problem. For f(x,y) = x^2 - y^2 at (1,1), find the direction of zero change and confirm with the dot product.
Answer. Direction (1,1)/sqrt(2) gives D_u f = 0
Problem. For f(x,y) = 3*x + 4*y^2, find the directional derivative at (1, 1) in the direction (1, 2).
Solution. grad f = (3, 8y); at (1,1) = (3, 8). Normalize (1,2): ||(1,2)|| = sqrt(5), u = (1/sqrt5, 2/sqrt5). D_u f = (3,8).(1/sqrt5, 2/sqrt5) = (3 + 16)/sqrt5 = 19/sqrt5 ≈ 8.497.
The tangent plane is the best flat approximation to a surface at a point, built directly from the gradient.
At a point (a,b) the tangent plane to z = f(x,y) is z = f(a,b) + f_x(a,b)(x-a) + f_y(a,b)(y-b). It is the multivariable analog of the tangent line and the first-order Taylor approximation: it matches the surface's height and both slopes at the contact point. Compactly, f(a + dx) ≈ f(a) + grad f(a) . dx. This linearization is the conceptual core of gradient-based optimization: each gradient step assumes the loss is locally well-approximated by its tangent plane, so moving along -grad f decreases the linear model and (for a small enough step) the true loss. Forward-mode automatic differentiation literally propagates this linear approximation. The error of the approximation is governed by the second derivatives (the Hessian).
Worked Example 1
Problem. Find the tangent plane to f(x,y) = x^2 + y^2 at (1, 2).
Answer. z = 5 + 2*(x-1) + 4*(y-2)
Worked Example 2
Problem. Use linear approximation to estimate f(1.1, 1.9) for f(x,y) = x^2 + y^2 near (1, 2).
Answer. f(1.1,1.9) ≈ 4.8
Worked Example 3
Problem. Find the tangent plane to f(x,y) = x*e^y at (2, 0).
Answer. z = x + 2*y
Problem. Find the tangent plane to f(x,y) = ln(x) + y^2 at (1, 2), then approximate f(1.05, 2).
Solution. f(1,2) = ln(1) + 4 = 4. f_x = 1/x -> 1; f_y = 2y -> 4. Plane: z = 4 + 1*(x-1) + 4*(y-2). Approximate at (1.05, 2): dx=(0.05, 0), z ≈ 4 + 1*0.05 + 4*0 = 4.05. (True: ln(1.05)+4 ≈ 4.0488.)
Among all directions you could step, the gradient is the one that increases f fastest — the reason its negative drives descent.
The directional derivative D_u f = grad f . u = ||grad f|| cos(theta) is maximized when cos(theta) = 1, i.e. when u points exactly along grad f. So the gradient direction is steepest ascent, with rate of increase equal to ||grad f||; the opposite direction, -grad f, is steepest descent. This is the single most important fact for optimization: gradient descent updates parameters by theta <- theta - eta*grad L precisely because -grad L is the locally fastest way to reduce the loss. Perpendicular to the gradient the rate is zero — that is why level sets (contours) are orthogonal to the gradient. The magnitude ||grad f|| also tells you how aggressive a step you can take before the linear approximation breaks down.
Worked Example 1
Problem. For f(x,y) = x^2 + y^2 at (3, 4), give the direction of steepest ascent and the maximum rate of increase.
Answer. Direction (0.6, 0.8); max rate 10
Worked Example 2
Problem. For the loss L(w1,w2) = (w1-1)^2 + (w2+2)^2 at (0,0), give the steepest-descent direction.
Answer. Descent direction (1, -2)/sqrt5 (i.e. -grad L = (2,-4))
Worked Example 3
Problem. One gradient-descent step on L(w1,w2)=(w1-1)^2+(w2+2)^2 from (0,0) with learning rate 0.1.
Answer. New point (0.2, -0.4), loss fell 5 -> 3.2
Problem. For f(x,y) = 4*x + 3*y at the point (0,0), give the steepest-ascent unit direction and the maximum rate of increase.
Solution. grad f = (4, 3) everywhere (linear function). ||grad f|| = sqrt(16+9) = 5. Steepest-ascent unit direction = (4,3)/5 = (0.8, 0.6); maximum rate of increase = 5.
Gradient descent repeatedly steps downhill along the negative gradient to find a minimum of a loss surface.
Gradient descent is the workhorse optimizer of machine learning. Starting from an initial parameter vector theta_0, it iterates theta_{k+1} = theta_k - eta*grad L(theta_k), where eta is the learning rate. Each step moves opposite the gradient — the steepest-descent direction — by an amount proportional to the slope's magnitude, so steps shrink automatically as you approach a flat minimum where grad L ≈ 0. The learning rate is decisive: too small and convergence crawls; too large and the iterates oscillate or diverge. On a convex (bowl) loss the method converges to the global minimum; on non-convex losses it finds a local minimum or stalls near saddle points. Backpropagation exists precisely to supply grad L cheaply at every iteration of this loop.
Worked Example 1
Problem. Minimize f(x) = x^2 (one variable) from x0 = 4 with eta = 0.1. Do two steps.
Answer. x1 = 3.2, x2 = 2.56 (f decreasing toward 0)
Worked Example 2
Problem. For f(x,y) = x^2 + 10*y^2, take one step from (1, 1) with eta = 0.1 and note the zig-zag risk.
Answer. New point (0.8, -1.0); steep y-direction overshoots (ill-conditioned)
Worked Example 3
Problem. Show that with eta = 1.0 on f(x) = x^2, gradient descent from x0 = 1 fails to converge.
Answer. Diverges/oscillates between 1 and -1: eta too large
Problem. Minimize f(x,y) = (x-2)^2 + (y-3)^2 from (0,0) with eta = 0.5. Take one step and give the new loss.
Solution. grad f = (2(x-2), 2(y-3)) = (-4, -6) at (0,0). Step: (0,0) - 0.5*(-4,-6) = (2, 3). New loss = (2-2)^2 + (3-3)^2 = 0 — it reached the exact minimum in one step because eta = 0.5 is exactly 1/(2*curvature) for this isotropic bowl.
Automatic differentiation computes exact gradients from code, letting you verify hand derivations and scale to large models.
Hand-deriving gradients is essential for understanding but error-prone at scale, so practitioners rely on automatic differentiation (autograd, JAX, PyTorch). Autodiff records the operations of f as a computation graph and applies the chain rule mechanically, returning the exact gradient (not a finite-difference estimate). Two checks anchor the workflow. First, a finite-difference approximation, (f(x+h) - f(x-h))/(2h) per coordinate, gives a numerical gradient you can compare to the analytic one to catch bugs to within ~1e-6. Second, autograd's grad/jacobian functions give the analytic answer directly. In practice you always assert your hand gradient matches autograd before trusting it. This lesson connects everything: partials -> gradient -> the tooling that powers backpropagation in real frameworks.
Worked Example 1
Problem. Hand-compute grad f for f(x,y) = x^2*y + sin(x*y), then state what autograd should return at (1, 2).
Answer. grad f(1,2) ≈ (3.168, 0.584)
Worked Example 2
Problem. Estimate df/dx of f(x,y) = x^2*y at (3, 2) by central finite differences with h = 1e-4.
Answer. Numerical f_x ≈ 12.0, matching analytic 12
Worked Example 3
Problem. Explain why we use central differences (f(x+h)-f(x-h))/(2h) rather than forward differences for the gradient check.
Answer. Central differences have O(h^2) error, giving a tighter, more reliable check
Problem. For f(x,y) = x*y^2, compute the analytic gradient at (2,3) and verify the x-partial by a central finite difference with h = 1e-4 (you may approximate).
Solution. Analytic: f_x = y^2 = 9; f_y = 2*x*y = 12. So grad f(2,3) = (9, 12). Central difference for f_x: (f(2+h,3) - f(2-h,3))/(2h) = ((2.0001)*9 - (1.9999)*9)/(2e-4) = (18.0009 - 17.9991)/2e-4 = 0.0018/0.0002 = 9.0. Matches the analytic f_x = 9.
For f(x,y) = x^2*y + sin(x*y), compute both partial derivatives by hand, then implement f in Python and use autograd (or jax.grad) to compute the gradient automatically. Evaluate both at three test points and assert they agree to within 1e-6. Then run five steps of gradient descent from a starting point and print how f decreases.
Deliverable · A script that prints hand-derived and autograd gradients side by side with a passing numeric assertion, plus a short log of five descent steps showing the loss dropping.
1. When computing the partial derivative of f with respect to x, the other variables are treated as...
Answer B. A partial derivative differentiates with respect to one variable while holding all others constant.
2. The gradient vector of f points in the direction of...
Answer C. The gradient points toward the greatest rate of increase; its negative is the steepest descent direction.
3. In gradient descent, the update rule for parameters theta is approximately:
Answer B. We step opposite the gradient, scaled by the learning rate, to reduce the loss.
4. A directional derivative in the direction of a unit vector u equals:
Answer B. The directional derivative is grad(f) dot u, projecting the gradient onto the chosen direction.
5. What does the tangent plane to a surface at a point represent?
Answer B. The tangent plane is the multivariable analog of the tangent line: the local linear approximation.
I can compute partial derivatives and assemble them into a gradient vector.
I can compute a directional derivative and explain why the gradient is steepest ascent.
I can take a gradient-descent step and verify the gradient numerically with autograd.
A vector-to-vector function bundles several scalar fields into one map, the structure of every neural-network layer.
A function f: R^n -> R^m takes an n-vector input and returns an m-vector output; it is just m scalar component functions f_1, ..., f_m stacked together, each depending on all n inputs. A linear layer y = W x + b is exactly such a map (W is m-by-n). Thinking in these terms is essential because a deep network is a composition of vector-to-vector functions, and its training requires differentiating that composition. The right derivative object for an R^n -> R^m map is no longer a single gradient vector but a matrix — the Jacobian — collecting the partials of every output with respect to every input. This lesson sets up that structure: identify inputs (n), outputs (m), and the component functions before differentiating.
Worked Example 1
Problem. Identify n, m, and the component functions of f(x,y) = (x+y, x*y, x^2).
Answer. n=2, m=3; components f_1=x+y, f_2=x*y, f_3=x^2
Worked Example 2
Problem. Evaluate the linear map f(x) = W x with W = [[1,2],[0,3],[4,0]] at x = (1, -1).
Answer. f(1,-1) = (-1, -3, 4)
Worked Example 3
Problem. A network layer is a(x) = sigma(W x + b) with elementwise sigma. State its input/output dimensions if W is 4x3.
Answer. a: R^3 -> R^4 (n=3, m=4)
Problem. For f(x,y,z) = (x*y, y+z, x*z, z^2), state n and m and list the component functions.
Solution. Inputs x, y, z so n = 3. Four outputs so m = 4, giving f: R^3 -> R^4. Components: f_1 = x*y, f_2 = y+z, f_3 = x*z, f_4 = z^2.
The Jacobian stacks the gradients of every output component into one m-by-n matrix — the multi-output derivative.
For f: R^n -> R^m, the Jacobian J is the m-by-n matrix whose (i,j) entry is df_i/dx_j: row i is the gradient (transposed) of output component f_i, and column j collects how every output responds to input x_j. It is the best linear approximation of f near a point, generalizing the single-variable derivative and the gradient. For a linear map f(x) = W x the Jacobian is simply W. The Jacobian is the central object of backpropagation: the chain rule multiplies the Jacobians of successive layers, and reverse-mode AD forms vector-Jacobian products with them. Always get the shape right first (m rows, n columns); a shape mismatch is the most common error in matrix-calculus derivations.
Worked Example 1
Problem. Compute the Jacobian of f(x,y) = (x^2 + y, x*y, sin(x)).
Answer. J = [[2x, 1], [y, x], [cos x, 0]]
Worked Example 2
Problem. Evaluate the Jacobian of f(x,y) = (x^2 + y, x*y, sin x) at (0, 1).
Answer. J(0,1) = [[0,1],[1,0],[1,0]]
Worked Example 3
Problem. For the linear map f(x) = W x with W = [[2,0,1],[ -1,3,0]], give the Jacobian.
Answer. J = W = [[2,0,1],[-1,3,0]]
Problem. Find the Jacobian of f(x,y) = (e^x*y, x + y^2) and evaluate it at (0, 1).
Solution. m=2, n=2 -> 2x2. Row 1 (e^x*y): [d/dx, d/dy] = [e^x*y, e^x]. Row 2 (x + y^2): [1, 2y]. So J = [[e^x*y, e^x],[1, 2y]]. At (0,1): e^0=1, so J = [[1, 1],[1, 2]].
Near a point, a vector function is well-approximated by its value plus the Jacobian times the displacement.
The first-order (linear) approximation of f: R^n -> R^m near a point a is f(a + dx) ≈ f(a) + J(a) dx, where J(a) is the Jacobian evaluated at a. This is the vector generalization of the tangent-line/tangent-plane idea: J(a) dx is the linear map that best predicts how the output moves when the input is nudged by dx. It underlies forward-mode automatic differentiation, which propagates a small input perturbation through each layer by multiplying by that layer's Jacobian. In optimization and sensitivity analysis it answers 'if I perturb the inputs slightly, how do the outputs change?' The approximation error grows with the step size and is controlled by second derivatives, so it is trustworthy only for small dx.
Worked Example 1
Problem. For f(x,y) = (x^2, x*y) linearize near (1, 2) and estimate f(1.1, 2.1).
Answer. f(1.1,2.1) ≈ (1.2, 2.3)
Worked Example 2
Problem. Use the Jacobian to find how the output of f(x,y) = (x+y, x-y) changes for dx = (0.05, -0.02).
Answer. delta_f ≈ (0.03, 0.07)
Worked Example 3
Problem. For f(x) = (sin x, cos x) near x = 0, give the linear approximation and estimate f(0.1).
Answer. f(0.1) ≈ (0.1, 1.0)
Problem. For f(x,y) = (x*y, y^2) linearize near (2, 1) and estimate f(2.1, 0.9).
Solution. f(2,1) = (2, 1). J = [[y, x],[0, 2y]] = [[1, 2],[0, 2]] at (2,1). dx = (0.1, -0.1). J dx = (1*0.1 + 2*(-0.1), 0*0.1 + 2*(-0.1)) = (-0.1, -0.2). f ≈ (2,1) + (-0.1,-0.2) = (1.9, 0.8). (True: (1.89, 0.81).)
The Hessian gathers all second partial derivatives of a scalar function, encoding its local curvature.
For a scalar field f: R^n -> R, the Hessian H is the n-by-n matrix with entry H[i,j] = d^2 f / (dx_i dx_j). It is the Jacobian of the gradient, and it describes the curvature of f. For functions with continuous second partials, Clairaut's theorem guarantees mixed partials are equal (f_xy = f_yx), so the Hessian is symmetric. The Hessian is the quadratic term in the second-order Taylor expansion f(a+dx) ≈ f(a) + grad f(a).dx + (1/2) dx^T H dx, which is exactly the local bowl/saddle quadratic form from Unit 1. In optimization it powers Newton's method and the second-derivative test, and its eigenvalues reveal whether a critical point is a minimum, maximum, or saddle.
Worked Example 1
Problem. Compute the Hessian of f(x,y) = x^2 + x*y + 2*y^2.
Answer. H = [[2, 1], [1, 4]]
Worked Example 2
Problem. Find the Hessian of f(x,y) = x^3 + x*y^2 at (1, 2).
Answer. H(1,2) = [[6, 4], [4, 2]]
Worked Example 3
Problem. Verify symmetry of the Hessian for f(x,y) = e^(x*y).
Answer. f_xy = f_yx = (1 + x*y)*e^(x*y); Hessian symmetric
Problem. Compute the Hessian of f(x,y) = x^2*y + y^3 at (1, 1).
Solution. f_x = 2xy, f_y = x^2 + 3y^2. f_xx = 2y; f_xy = 2x; f_yy = 6y. At (1,1): f_xx = 2, f_xy = 2, f_yy = 6. H(1,1) = [[2, 2],[2, 6]] (symmetric).
The eigenvalues of the Hessian classify a critical point: all positive means a minimum, all negative a maximum, mixed a saddle.
At a critical point (where grad f = 0) the Hessian's eigenvalues determine the local shape. If all eigenvalues are positive the Hessian is positive definite and the point is a local minimum; all negative gives a local maximum; mixed signs give a saddle; a zero eigenvalue makes the test inconclusive. This is the second-derivative test, the multivariable generalization of f'' > 0. A function is convex if its Hessian is positive semidefinite everywhere, which is the gold standard in optimization: for a convex loss every local minimum is the global minimum, so gradient descent cannot get trapped. For 2x2 Hessians a quick proxy is the determinant test: D = f_xx*f_yy - f_xy^2 > 0 with f_xx > 0 means a minimum; D < 0 means a saddle.
Worked Example 1
Problem. Classify the critical point of f(x,y) = x^2 + x*y + 2*y^2 at the origin using the Hessian.
Answer. Local minimum at the origin (positive definite Hessian)
Worked Example 2
Problem. Classify the critical point of g(x,y) = x^2 - y^2 at the origin.
Answer. Saddle point (eigenvalues +2 and -2)
Worked Example 3
Problem. Find and classify the critical point of f(x,y) = x^2 + 4*y^2 - 2*x - 8*y.
Answer. Local minimum at (1, 1)
Problem. Classify the critical point of f(x,y) = x^2 - 4*x + y^2 + 2*x*y at its critical point.
Solution. grad f = (2x - 4 + 2y, 2y + 2x) = 0. From the second equation x = -y; substitute into the first: 2(-y) - 4 + 2y = -4 ≠ 0, so there is NO solution — the gradient is never zero, meaning no critical point. Checking the Hessian H = [[2,2],[2,2]] confirms it is singular (eigenvalues 0 and 4, positive semidefinite but degenerate), consistent with a function having no isolated critical point (a parabolic trough heading to -infinity along one direction).
Saddle points are critical points that trap naive optimizers; high-dimensional loss landscapes are full of them.
A saddle point has a zero gradient but an indefinite Hessian — positive curvature along some directions and negative along others — so it is a minimum in some slices and a maximum in others. Gradient descent slows dramatically near saddles because the gradient is tiny, yet they are not minima. In high-dimensional deep-learning losses, critical points are overwhelmingly saddles rather than bad local minima, which is why momentum, stochastic noise, and second-order information help escape them: a direction of negative curvature (an eigenvector of a negative Hessian eigenvalue) is a guaranteed descent direction even when the gradient nearly vanishes. Understanding curvature, summarized by the Hessian's eigenvalues, explains both why optimization stalls and how to design optimizers that don't.
Worked Example 1
Problem. Show that f(x,y) = x^2 - y^2 has a saddle at the origin and give a descent direction.
Answer. Saddle at origin; descend along the y-axis (negative-curvature eigenvector)
Worked Example 2
Problem. Explain why gradient descent stalls near the saddle of f(x,y) = x^2 - y^2 if started near (epsilon, 0).
Answer. Near-zero gradient on the y=0 line makes steps vanish; the optimizer stalls without noise
Worked Example 3
Problem. Classify the critical point of the monkey-saddle-like f(x,y) = x^2 + y^2 - 3*x^2*y... evaluate curvature at the origin only.
Answer. Origin is a local minimum (H = 2I, positive definite)
Problem. For f(x,y) = x*y, find the critical point, build the Hessian, and classify it.
Solution. grad f = (y, x) = 0 only at the origin. f_xx = 0, f_yy = 0, f_xy = 1, so H = [[0,1],[1,0]]. Eigenvalues solve lambda^2 - 1 = 0 -> lambda = +1, -1. Mixed signs (indefinite) -> the origin is a saddle point. A descent direction is the eigenvector of -1, namely (1,-1)/sqrt2, along which f = -t^2/2 decreases.
For f(x,y) = (x^2 + y, x*y, sin(x)) compute the 3x2 Jacobian, and for g(x,y) = x^2 + x*y + 2*y^2 compute the 2x2 Hessian, both by hand. Verify each with autograd's jacobian and hessian. Then compute the Hessian's eigenvalues with NumPy and state whether g has a minimum, maximum, or saddle at the origin.
Deliverable · A script that prints the hand and autodiff Jacobian/Hessian, asserts they match, and reports the Hessian eigenvalues with a one-line critical-point classification.
1. For a function f: R^4 -> R^2, the Jacobian matrix has dimensions:
Answer B. The Jacobian is m x n: outputs (2) as rows and inputs (4) as columns, so 2x4.
2. The Hessian matrix collects which derivatives of a scalar function?
Answer B. The Hessian is the matrix of all second-order partial derivatives, describing curvature.
3. If the Hessian at a critical point is positive definite (all eigenvalues > 0), the point is a:
Answer C. Positive definite curvature means the surface bows upward in every direction: a local minimum.
4. A saddle point occurs when the Hessian's eigenvalues are...
Answer C. Mixed-sign eigenvalues mean the function curves up in some directions and down in others: a saddle.
5. Why is convexity desirable in optimization?
Answer B. For convex functions every local minimum is global, so gradient descent cannot get stuck in a worse minimum.
I can assemble the Jacobian of a vector-valued function and use it to linearize.
I can build the Hessian and use its eigenvalues to classify a critical point.
I can explain how curvature and convexity affect optimization.
The multivariable chain rule multiplies Jacobians along a composition — the single identity that backpropagation is built on.
If z = g(y) and y = f(x), then the composition h(x) = g(f(x)) has Jacobian J_h(x) = J_g(f(x)) * J_f(x): you multiply the outer function's Jacobian (evaluated at the intermediate value) by the inner function's Jacobian. For a scalar output composed through one intermediate variable, this collapses to the familiar dz/dx = (dz/dy)(dy/dx), and when y depends on x through several paths you SUM over the paths: dz/dt = sum_i (dz/dy_i)(dy_i/dt). This single rule is exactly what a neural network needs: a network is a chain of layer functions, and the gradient of the loss is the product of all the layer Jacobians. Getting the multiplication order and shapes right is the whole game of backprop.
Worked Example 1
Problem. Let z = y^2 with y = 3*x + 1. Find dz/dx at x = 1.
Answer. dz/dx = 6y = 24 at x=1
Worked Example 2
Problem. Let z = x^2 + y^2 with x = cos t, y = sin t. Find dz/dt.
Answer. dz/dt = 0
Worked Example 3
Problem. Compose linear maps z = B y, y = A x with A = [[1,2],[0,1]], B = [[2,0],[1,3]]. Find the Jacobian dz/dx.
Answer. dz/dx = B*A = [[2,4],[1,5]]
Problem. Let z = e^u with u = x*y. Using the chain rule, find dz/dx and dz/dy, then evaluate at (1, 2).
Solution. dz/du = e^u. u_x = y, u_y = x. So dz/dx = e^u*y, dz/dy = e^u*x. At (1,2): u = 2, e^2 ≈ 7.389. dz/dx = 7.389*2 ≈ 14.78; dz/dy = 7.389*1 ≈ 7.389.
A computation graph turns a formula into a DAG of operations, making the chain rule a mechanical walk over edges.
A computation graph is a directed acyclic graph whose nodes are elementary operations (add, multiply, exp, ...) and whose edges carry intermediate values from inputs to the final output. Building it for an expression exposes every dependency path, and the chain rule becomes a bookkeeping rule over edges: the derivative of the output with respect to a node is the sum, over all outgoing edges, of (local edge derivative) times (downstream derivative). Each edge stores a local gradient — the partial derivative of the child node with respect to the parent. Frameworks like PyTorch and JAX build this graph automatically during the forward pass and then traverse it to differentiate. Drawing the graph by hand for a small expression is the clearest way to see why backprop visits nodes in reverse topological order.
Worked Example 1
Problem. Draw the computation graph (in words) for f = (x + y) * z and list the intermediate nodes.
Answer. Nodes: a = x+y, f = a*z; edges x->a, y->a, a->f, z->f
Worked Example 2
Problem. For f = (x + y)*z, give the local gradient on each edge.
Answer. x->a:1, y->a:1, a->f:z, z->f:(x+y)
Worked Example 3
Problem. Use the graph of f = (x+y)*z to compute df/dx and df/dz at (x,y,z) = (1,2,3).
Answer. df/dx = 3, df/dz = 3 at (1,2,3)
Problem. Build the computation graph for f = x * sin(y) and give df/dx and df/dy at (2, 0).
Solution. Nodes: s = sin(y) (edge y->s with local grad cos(y)); f = x*s (edges x->f local grad s, s->f local grad x). df/dx = s = sin(y); at (2,0) = sin0 = 0. df/dy = (df/ds)*(ds/dy) = x*cos(y); at (2,0) = 2*cos0 = 2. Answer: df/dx = 0, df/dy = 2.
Both modes apply the chain rule, but in opposite directions — and the choice decides whether differentiating is cheap or expensive.
Automatic differentiation comes in two modes. Forward mode propagates derivatives from inputs to outputs, computing one column of the Jacobian (the sensitivity to one input) per pass; it is cheap when there are few inputs. Reverse mode propagates derivatives from outputs back to inputs, computing one row of the Jacobian (the sensitivity of one output) per pass; it is cheap when there are few outputs. Machine learning has a scalar loss (one output) and millions of parameters (many inputs), so reverse mode — backpropagation — wins overwhelmingly: it gets the entire gradient in a single backward pass, whereas forward mode would need one pass per parameter. The cost rule of thumb: forward mode scales with the number of inputs, reverse mode with the number of outputs.
Worked Example 1
Problem. A function has 1,000,000 inputs and 1 output. How many passes does each mode need to get the full gradient?
Answer. Forward: 1,000,000 passes; Reverse: 1 pass
Worked Example 2
Problem. A function f: R^2 -> R^5 (2 inputs, 5 outputs). Which mode computes the full Jacobian with fewer passes?
Answer. Forward mode (2 passes vs 5)
Worked Example 3
Problem. For z = sin(x*y) with x, y inputs, sketch one forward-mode pass tracking dx (set dx=1, dy=0) to get dz/dx.
Answer. dz/dx = y*cos(x*y) (single forward pass with dx=1, dy=0)
Problem. A neural network loss has 10 million parameters and outputs a single scalar loss. State which AD mode to use and how many backward passes are needed for the full gradient.
Solution. Use reverse mode (backpropagation). Because there is exactly one output (the scalar loss), reverse mode obtains the gradient with respect to all 10 million parameters in a single backward pass. Forward mode would require 10 million passes (one per parameter), so it is infeasible here.
Backpropagation is reverse-mode AD on a network: a forward pass caches values, then a backward pass propagates gradients output-to-input.
Backpropagation computes the gradient of a scalar loss with respect to every parameter in two sweeps. The forward pass evaluates the network layer by layer, caching each intermediate value. The backward pass starts with dL/dL = 1 at the output and walks the computation graph in reverse topological order; at each node it multiplies the incoming upstream gradient by the node's local gradient (the chain rule) and accumulates contributions where a value fanned out. The quantity flowing backward into a node is often called its 'delta' or adjoint. Because every parameter's gradient is obtained in this one backward pass, the cost is roughly the same as the forward pass regardless of how many parameters there are — the property that makes training deep networks feasible. Reuse of cached forward values is what keeps it efficient.
Worked Example 1
Problem. For L = (a - t)^2 with a = w*x, do the backward pass to get dL/dw at w=2, x=3, t=5.
Answer. dL/dw = 6
Worked Example 2
Problem. For L = (sigma(w*x) - t)^2 with sigma the logistic function, give dL/dw symbolically.
Answer. dL/dw = 2*(a - t)*a*(1 - a)*x, where a = sigma(w*x)
Worked Example 3
Problem. Compute dL/dw numerically for Example 2 at w=0, x=1, t=1.
Answer. dL/dw = -0.25
Problem. For L = (a - t)^2 with a = relu(w*x), w=2, x=3, t=5, do the backward pass for dL/dw (relu'(z)=1 for z>0).
Solution. Forward: z = w*x = 6 > 0 so a = relu(6) = 6; L = (6-5)^2 = 1. Backward: dL/da = 2(a-t) = 2(1) = 2. da/dz = relu'(6) = 1. dz/dw = x = 3. dL/dw = 2*1*3 = 6.
A single dense layer's parameter gradients follow directly from the matrix chain rule, with neat transpose patterns.
Consider one layer z = W x + b, a = sigma(z), feeding a scalar loss L. Given the upstream gradient dL/da, backprop produces the layer's gradients with the matrix chain rule. First dL/dz = dL/da (elementwise-times) sigma'(z), the 'delta' for this layer. Then the parameter gradients are dL/dW = delta * x^T (an outer product giving a matrix the same shape as W) and dL/db = delta, and the gradient passed to the previous layer is dL/dx = W^T * delta. The transposes are not arbitrary: they are exactly what makes the shapes line up (W is m-by-n, so dL/dW must be m-by-n, forcing the outer product delta x^T). These three formulas are the entire backward pass of a dense layer and recur in every framework.
Worked Example 1
Problem. For z = W x + b with no activation and L scalar, given delta = dL/dz, give dL/dW, dL/db, dL/dx.
Answer. dL/dW = delta x^T, dL/db = delta, dL/dx = W^T delta
Worked Example 2
Problem. With x = (1, 2), W = [[1,0],[0,1]], delta = dL/dz = (3, 4), compute dL/dW and dL/dx.
Answer. dL/dW = [[3,6],[4,8]], dL/dx = (3, 4)
Worked Example 3
Problem. Add a sigmoid: a = sigma(z), with dL/da = (0.5, -0.5) and z = (0, 0). Find delta = dL/dz.
Answer. delta = dL/dz = (0.125, -0.125)
Problem. For a layer z = W x with W = [[2,1]], x = (1,3) (so z is scalar), and dL/dz = 5, compute dL/dW and dL/dx.
Solution. z = 2*1 + 1*3 = 5 (1x1). dL/dW = delta x^T = 5 * (1,3) = [[5, 15]] (same 1x2 shape as W). dL/dx = W^T delta = (2,1)^T * 5 = (10, 5). Answer: dL/dW = [[5, 15]], dL/dx = (10, 5).
Reverse mode never forms a Jacobian explicitly; it computes the vector-Jacobian product v^T J, which is all backprop ever needs.
The core operation of reverse-mode AD is the vector-Jacobian product (VJP): given an upstream gradient row vector v = dL/dy, a layer y = f(x) returns v^T J_f = dL/dx without ever materializing the full Jacobian J. This is crucial because J can be enormous (outputs times inputs), but the VJP only needs an efficient rule for 'multiply a vector by J^T'. For a linear layer y = W x the VJP is just W^T v; for an elementwise activation it is v times the diagonal of derivatives. Backprop chains VJPs from the loss backward through every layer. Because each layer supplies a cheap VJP rule rather than a dense Jacobian, reverse mode stays memory- and compute-efficient even for networks with billions of parameters — the reason ML standardizes on it.
Worked Example 1
Problem. For y = W x with W = [[1,2,3],[4,5,6]] (2x3) and upstream v = (1, 1), compute the VJP v^T J = dL/dx.
Answer. dL/dx = (5, 7, 9)
Worked Example 2
Problem. For an elementwise y_i = x_i^2 with upstream v = (2, 3) at x = (1, 4), give the VJP.
Answer. dL/dx = (4, 24)
Worked Example 3
Problem. Explain why computing v^T J is cheaper than forming J for a layer with 10^4 inputs and 10^4 outputs.
Answer. VJP avoids the 10^8-entry matrix, using an O(n) rule instead
Problem. For y = W x with W = [[2,0],[1,3]] and upstream gradient v = (1, 2), compute the vector-Jacobian product dL/dx.
Solution. For a linear layer J = W, so the VJP is W^T v. W^T = [[2,1],[0,3]]. W^T v = (2*1 + 1*2, 0*1 + 3*2) = (4, 6). Answer: dL/dx = (4, 6) — obtained without ever forming a separate Jacobian.
A finite-difference gradient check is the standard sanity test confirming a hand-derived backward pass is correct.
Even correct-looking backprop code hides sign and transpose bugs, so we validate it with a numerical gradient check. For each parameter we perturb it by a tiny h and form the central difference (L(theta + h) - L(theta - h)) / (2h), comparing this numerical gradient to the analytic one from backprop. We score agreement with the relative error |g_analytic - g_numeric| / (|g_analytic| + |g_numeric| + eps); values below about 1e-5 to 1e-7 indicate a correct implementation, while anything near 1e-2 signals a bug. The central difference is used (not forward) because its error is O(h^2), giving a tight check. This is mandatory discipline before trusting a manual backward pass and is the bridge from the math of this course to reliable ML code.
Worked Example 1
Problem. For L(w) = w^2 with analytic gradient 2w, verify at w = 3 using central differences with h = 1e-4.
Answer. Numeric 6.0 matches analytic 6 (check passes)
Worked Example 2
Problem. Compute the relative error between g_analytic = 6.0 and g_numeric = 6.00003.
Answer. Relative error ≈ 2.5e-6 (passes)
Worked Example 3
Problem. Backprop returns -6 but the numerical gradient is +6 at w=3. Diagnose the bug.
Answer. Sign error in backprop (relative error 1.0); fix the missing minus sign
Problem. For L(w) = 3*w^2 with analytic gradient 6w, do a central-difference check at w = 2 with h = 1e-3 and report the relative error.
Solution. Analytic: 6*2 = 12. L(2.001) = 3*(2.001)^2 = 3*4.004001 = 12.012003; L(1.999) = 3*(1.999)^2 = 3*3.996001 = 11.988003. Numeric = (12.012003 - 11.988003)/(2e-3) = 0.024/0.002 = 12.0. Relative error = |12 - 12|/(12+12) = 0. The check passes.
Build a one-hidden-layer network y = W2 @ sigma(W1 @ x) with a squared-error loss in NumPy. Derive the gradients of the loss with respect to W1 and W2 using the multivariable chain rule, implement them as a manual backward pass, and confirm each against a finite-difference numerical gradient check to within 1e-5. Explain why reverse mode is the right choice here.
Deliverable · A NumPy script with a manual forward and backward pass plus a gradient-check function that prints the max relative error for each weight matrix.
1. The multivariable chain rule combines the derivatives of composed functions by...
Answer B. For composed maps, the Jacobian of the composition is the product of the individual Jacobians.
2. Backpropagation is essentially which form of automatic differentiation?
Answer C. Backprop is reverse-mode AD: gradients flow from the scalar loss backward to the parameters.
3. Reverse-mode AD is most efficient when a function has...
Answer A. A scalar loss with millions of parameters is the canonical case; reverse mode gets all gradients in one pass.
4. A computation graph is best described as a:
Answer B. It is a DAG: nodes are operations, edges carry intermediate results, with no cycles.
5. Why do we run a numerical gradient check during backprop development?
Answer B. Finite-difference checks validate that the hand-derived backward pass matches the true gradient.
I can apply the multivariable chain rule across a computation graph.
I can explain why reverse-mode (backprop) is efficient for scalar losses.
I can derive and numerically gradient-check the gradients of a small network.
A double integral sums a function's values over a 2D region, giving the signed volume under its surface.
The double integral of f(x,y) over a region R, written ∫∫_R f dA, is the limit of Riemann sums of f over shrinking area elements dA = dx dy — the signed volume between the surface z = f(x,y) and the xy-plane. Over a rectangle [a,b]x[c,d] the bounds are constant; over a general region one variable's limits depend on the other (e.g. y runs from g1(x) to g2(x)). In data science the double integral is how we compute total probability and expectations from a joint density, and how we normalize 2D distributions. The key skill is describing the region correctly as bounds before integrating; a sketch of R prevents the most common errors. The integral of f = 1 over R simply returns the area of R.
Worked Example 1
Problem. Evaluate ∫∫_R x*y dA over the rectangle R = [0,1]x[0,2].
Answer. 1
Worked Example 2
Problem. Evaluate ∫∫_R (x + y) dA over R = [0,2]x[1,3].
Answer. 12
Worked Example 3
Problem. Evaluate ∫∫_R x dA over the triangle R bounded by y=0, x=1, y=x.
Answer. 1/3
Problem. Evaluate ∫∫_R y dA over R = [0,1]x[0,1].
Solution. Inner over y: ∫_0^1 y dy = [y^2/2]_0^1 = 1/2 (the x-integration multiplies a constant). Outer over x: ∫_0^1 (1/2) dx = 1/2. Answer: 1/2.
Fubini's theorem lets us evaluate a double integral as two single integrals, and the order can be swapped when it helps.
An iterated integral evaluates a double integral one variable at a time: integrate the inner variable first (treating the outer as constant), then the outer. Fubini's theorem says that for a well-behaved f the result is independent of order: ∫∫ f dx dy = ∫∫ f dy dx. Choosing the easier order can drastically simplify a problem — sometimes one order is elementary while the other has no closed form. For a general region, swapping the order requires re-describing the region's bounds (redraw it and read the limits the other way), which is the trickiest step. This skill matters in probability when marginalizing a joint density: you integrate out one variable, and picking the convenient order avoids intractable integrals.
Worked Example 1
Problem. Evaluate ∫_0^1 ∫_0^2 (x + y) dy dx, then redo as dx dy to confirm Fubini.
Answer. 3 (order-independent)
Worked Example 2
Problem. Swap the order of ∫_0^1 ∫_x^1 f dy dx into dx dy form (describe the new limits).
Answer. ∫_0^1 ∫_0^y f dx dy
Worked Example 3
Problem. Use a convenient order to evaluate ∫_0^1 ∫_y^1 e^(x^2) dx dy (the dx integral has no elementary form).
Answer. (e - 1)/2 ≈ 0.859
Problem. Evaluate ∫_0^2 ∫_0^1 x^2*y dx dy.
Solution. Inner over x: ∫_0^1 x^2*y dx = y*[x^3/3]_0^1 = y/3. Outer over y: ∫_0^2 (y/3) dy = (1/3)*[y^2/2]_0^2 = (1/3)*2 = 2/3. Answer: 2/3.
A triple integral extends integration to 3D regions, computing volumes, masses, and 3D expectations.
The triple integral ∫∫∫_E f dV integrates over a solid region E with volume element dV = dx dy dz. Evaluated as three nested single integrals, the innermost limits may depend on the two outer variables and the middle on the outermost. Setting f = 1 returns the volume of E. Triple integrals compute mass (density integrated over a solid), centers of mass, and — most relevant to data science — probabilities and expectations under three-dimensional joint densities. As with double integrals, the hard part is describing the region as a nested set of bounds; a clear ordering (e.g. z from a surface to a surface, then y, then x) keeps the setup tractable. The mechanics are just repeated single-variable integration with care about which variables are held constant.
Worked Example 1
Problem. Evaluate ∫∫∫_E 1 dV over the box E = [0,1]x[0,2]x[0,3] (its volume).
Answer. 6
Worked Example 2
Problem. Evaluate ∫_0^1 ∫_0^1 ∫_0^1 x*y*z dz dy dx.
Answer. 1/8
Worked Example 3
Problem. Find the volume of the tetrahedron E under x + y + z <= 1 in the first octant.
Answer. Volume = 1/6
Problem. Evaluate ∫_0^2 ∫_0^1 ∫_0^3 (x + z) dz dy dx.
Solution. Inner over z: ∫_0^3 (x+z) dz = 3x + [z^2/2]_0^3 = 3x + 9/2. Over y (0 to 1): multiply by (1-0)=1, still 3x + 9/2. Over x: ∫_0^2 (3x + 9/2) dx = [3x^2/2 + (9/2)x]_0^2 = 6 + 9 = 15. Answer: 15.
Substituting new coordinates rescales the area/volume element by the absolute Jacobian determinant.
When you change variables from (x,y) to (u,v) via a map x=x(u,v), y=y(u,v), the area element transforms as dx dy = |det J| du dv, where J is the Jacobian matrix of partials [[x_u, x_v],[y_u, y_v]]. The absolute value of its determinant measures how the transformation stretches or shrinks area locally — the multivariable generalization of the single-variable u-substitution factor du/dx. This is the bridge back to Unit 3's Jacobian: there it linearized a map; here its determinant rescales volume. It is essential in probability for the change-of-variables formula for densities (a transformed random variable's density picks up a |det J^{-1}| factor) and underlies normalizing flows in modern generative models. Always take the absolute value — negative determinants flip orientation but area is positive.
Worked Example 1
Problem. For the polar map x = r*cos t, y = r*sin t, compute the Jacobian determinant.
Answer. det J = r, giving dx dy = r dr dtheta
Worked Example 2
Problem. For the linear map u = 2x, v = 3y, find the factor relating du dv to dx dy.
Answer. dx dy = (1/6) du dv
Worked Example 3
Problem. Use polar coordinates to evaluate ∫∫_D (x^2 + y^2) dA over the unit disk D.
Answer. pi/2
Problem. For the map u = x + y, v = x - y, find |det J| (with x,y as functions of u,v) and the relation between dx dy and du dv.
Solution. Solve: x = (u+v)/2, y = (u-v)/2. J = [[x_u, x_v],[y_u, y_v]] = [[1/2, 1/2],[1/2, -1/2]]. det J = (1/2)(-1/2) - (1/2)(1/2) = -1/4 - 1/4 = -1/2, so |det J| = 1/2. Therefore dx dy = (1/2) du dv.
Curvilinear coordinates match the symmetry of a region, replacing messy bounds with simple ones plus a Jacobian factor.
When a region has circular or spherical symmetry, switching coordinates simplifies both the bounds and the integrand. Polar (2D): x=r cos t, y=r sin t, dA = r dr dt. Cylindrical (3D): add z unchanged, dV = r dr dt dz. Spherical (3D): x=rho sin phi cos t, y=rho sin phi sin t, z=rho cos phi, with dV = rho^2 sin phi drho dphi dt. The extra factors (r, r, rho^2 sin phi) are the Jacobian determinants from the previous lesson. These coordinates turn a disk, cylinder, or ball — awkward in Cartesian limits — into simple rectangular ranges. In statistics this is exactly why a 2D Gaussian is integrated in polar form: the r^2 in the exponent and the r Jacobian combine into an elementary integral, proving the normalization constant.
Worked Example 1
Problem. Find the area of the disk of radius 2 using polar coordinates.
Answer. 4pi
Worked Example 2
Problem. Find the volume of a cylinder of radius 1 and height 3 using cylindrical coordinates.
Answer. 3*pi
Worked Example 3
Problem. Find the volume of a ball of radius 1 using spherical coordinates.
Answer. 4*pi/3
Problem. Use polar coordinates to evaluate ∫∫_D dA where D is the disk x^2 + y^2 <= 9 (confirm it equals the area).
Solution. In polar, D is 0<=r<=3, 0<=t<=2pi, with dA = r dr dt. ∫_0^{2pi}∫_0^3 r dr dt = ∫_0^{2pi} [r^2/2]_0^3 dt = ∫_0^{2pi} (9/2) dt = (9/2)(2pi) = 9pi. This matches pi*r^2 = pi*9 = 9pi, the area of a radius-3 disk.
Expectations and probabilities are integrals of quantities weighted by a density — the link from calculus to statistics and ML.
For a continuous random vector with joint density p(x,y), probability is the integral of p over a region and the expectation of any g is E[g] = ∫∫ g(x,y) p(x,y) dA. A valid density is nonnegative and integrates to 1 over its support — that normalization is itself a double integral. Marginalizing integrates out a variable, p_X(x) = ∫ p(x,y) dy, and conditional densities, moments (mean, variance), and the log-likelihoods minimized in ML all reduce to such integrals. This lesson ties the unit's machinery to data science: every expectation in a loss, every normalizing constant in a model, every Bayesian marginal is a multiple integral, often evaluated numerically (grid sums or Monte Carlo) when no closed form exists. Change of variables (the Jacobian) is what lets us transform densities cleanly.
Worked Example 1
Problem. Verify p(x,y) = 4*x*y on [0,1]x[0,1] is a valid density (integrates to 1).
Answer. Integral = 1, so p is a valid density
Worked Example 2
Problem. For that density p = 4xy, compute E[x] = ∫∫ x p dA.
Answer. E[x] = 2/3
Worked Example 3
Problem. Find the marginal p_X(x) of p(x,y) = 4xy on the unit square.
Answer. p_X(x) = 2x for x in [0,1]
Problem. For the density p(x,y) = 4xy on [0,1]x[0,1], compute E[x*y] = ∫∫ x*y*p dA.
Solution. E[xy] = ∫_0^1∫_0^1 xy*4xy dx dy = ∫_0^1∫_0^1 4x^2 y^2 dx dy. Inner over x: 4y^2*[x^3/3]_0^1 = (4/3)y^2. Outer over y: ∫_0^1 (4/3)y^2 dy = (4/3)(1/3) = 4/9. Answer: E[xy] = 4/9.
Take a 2D standard normal density and verify by numerical double integration (e.g. NumPy grid sum or scipy.integrate.dblquad) that it integrates to 1 over a large region. Then compute the expectation of x^2 + y^2 under this density numerically and compare it to the analytic value of 2. Briefly note how the Jacobian determinant would appear if you switched to polar coordinates.
Deliverable · A script reporting the numerical integral of the density (close to 1.0), the numerical expectation of x^2+y^2, and a short note on the polar Jacobian factor r.
1. An iterated double integral is evaluated by integrating...
Answer B. We evaluate the inner integral first (treating the outer variable as constant), then the outer integral.
2. When changing variables in a multiple integral, you must multiply by the absolute value of the:
Answer C. The Jacobian determinant rescales area/volume elements under the coordinate transformation.
3. Converting a 2D integral to polar coordinates introduces what extra factor?
Answer A. The area element dx dy becomes r dr dtheta; the Jacobian of the polar map is r.
4. The expectation of a function g under a probability density p is computed as:
Answer B. Expectation is the integral of g(x) weighted by the density p(x) over all outcomes.
5. A valid joint probability density must integrate to what value over its entire support?
Answer C. Total probability is 1, so any valid density integrates to one over its full domain.
I can set up and evaluate double and triple integrals over given regions.
I can apply a change of variables using the Jacobian determinant.
I can interpret a multiple integral as an expectation over a probability density.
A vector field attaches a vector to every point in space, modeling flows, forces, and the gradient of a loss.
A vector field F: R^n -> R^n assigns a vector to each point, e.g. F(x,y) = (P(x,y), Q(x,y)) in 2D. You visualize it by drawing arrows whose direction and length vary with position — a velocity field of a fluid, a force field, or, crucially for ML, the negative-gradient field -grad L whose arrows show which way gradient descent flows at each point. A streamline is a curve everywhere tangent to the field, the path a particle follows. Distinguishing the two main behaviors — spreading out (captured by divergence) versus swirling (captured by curl) — is the goal of this unit. Recognizing that the gradient of a scalar field is itself a vector field connects optimization to the divergence/curl machinery developed next.
Worked Example 1
Problem. Describe the field F(x,y) = (x, y) and evaluate it at (1, 2) and (-1, 0).
Answer. A radial 'source' field; F(1,2)=(1,2), F(-1,0)=(-1,0)
Worked Example 2
Problem. Describe the rotational field F(x,y) = (-y, x) and evaluate at (1, 0) and (0, 1).
Answer. Counterclockwise rotation; F(1,0)=(0,1), F(0,1)=(-1,0)
Worked Example 3
Problem. Show that F(x,y) = grad f for f(x,y) = x^2 + y^2 is a vector field and identify it.
Answer. F = (2x, 2y), the radial gradient field of x^2+y^2
Problem. For F(x,y) = (y, -x), evaluate the field at (1,1) and (0,2) and describe the overall flow.
Solution. F(1,1) = (1, -1) (points down-right); F(0,2) = (2, 0) (points right). Tracing these, the arrows circulate clockwise around the origin — F = (y, -x) is a clockwise rotation field (the opposite sense to (-y, x)).
Divergence is a scalar measuring how much a vector field spreads out from a point — net outflow per unit volume.
The divergence of F = (P, Q) in 2D (or (P,Q,R) in 3D) is div F = dP/dx + dQ/dy (+ dR/dz), the dot product of the del operator with F, nabla . F. It is a scalar field: positive where the field acts as a source (arrows spreading out), negative at a sink (arrows converging), and zero for an incompressible field. Intuitively it measures net outflow per unit volume around a point. In physics it appears in the continuity and Maxwell equations; in data science the divergence of a probability flow appears in continuity equations for diffusion models and in the change-of-density formula for continuous normalizing flows, where the log-density evolves by the integral of -div of the velocity field. Computing it is just summing the diagonal partials.
Worked Example 1
Problem. Compute div F for the radial field F(x,y) = (x, y).
Answer. div F = 2
Worked Example 2
Problem. Compute div F for the rotation field F(x,y) = (-y, x).
Answer. div F = 0 (divergence-free)
Worked Example 3
Problem. Compute div F for F(x,y,z) = (x^2, x*y, z) at (1, 2, 3).
Answer. div F = 3x + 1; at (1,2,3) it is 4
Problem. Compute the divergence of F(x,y) = (x^2*y, y^2) and evaluate at (2, 1).
Solution. P = x^2*y -> dP/dx = 2xy. Q = y^2 -> dQ/dy = 2y. div F = 2xy + 2y. At (2,1): 2*2*1 + 2*1 = 4 + 2 = 6. Answer: div F = 2xy + 2y, equal to 6 at (2,1).
Curl measures the local rotation or swirl of a vector field; it is a vector in 3D and a scalar in 2D.
The curl of a vector field measures how much it circulates locally. In 3D, curl F = nabla x F is a vector whose direction is the rotation axis and whose magnitude is twice the local angular speed. In 2D, F = (P, Q) has a scalar curl dQ/dx - dP/dy, the z-component of the 3D curl. A field with zero curl everywhere is called irrotational; a key theorem says a field is conservative (a gradient field) exactly when its curl is zero on a simply connected domain — directly relevant to optimization, since every gradient field -grad L is curl-free. Computing 2D curl is a single subtraction of cross-partials; the 3D version is the determinant-style expansion of nabla x F. Curl and divergence together fully characterize a smooth field's local behavior.
Worked Example 1
Problem. Compute the 2D curl of the rotation field F(x,y) = (-y, x).
Answer. curl F = 2
Worked Example 2
Problem. Compute the 2D curl of the radial field F(x,y) = (x, y).
Answer. curl F = 0 (irrotational)
Worked Example 3
Problem. Compute the 2D curl of F(x,y) = (x*y, x^2) and evaluate at (1, 2).
Answer. curl F = x; at (1,2) it is 1
Problem. Compute the 2D curl of F(x,y) = (3*y, -2*x) and state whether the field rotates.
Solution. P = 3y -> dP/dy = 3. Q = -2x -> dQ/dx = -2. curl = dQ/dx - dP/dy = -2 - 3 = -5. The curl is a nonzero negative constant, so the field rotates clockwise everywhere with constant strength 5.
A line integral sums a field's effect along a curve; for conservative fields it depends only on the endpoints.
A line integral ∫_C F . dr accumulates the component of a vector field along a curve C — for a force field it is the work done moving along the path. Parametrizing C by r(t), t in [a,b], it becomes ∫_a^b F(r(t)) . r'(t) dt, a single integral. A field F is conservative if F = grad(phi) for some scalar potential phi; then the Fundamental Theorem of Line Integrals gives ∫_C F . dr = phi(end) - phi(start), independent of the path. Equivalently, conservative fields are exactly the curl-free fields (on simply connected domains), and their loop integrals vanish. This matters in ML because -grad L is conservative with potential L: the change in loss along any parameter path equals the difference in loss values, not the route taken — the very property gradient descent exploits.
Worked Example 1
Problem. Evaluate ∫_C F . dr for F = (y, x) along r(t) = (t, t^2), t in [0,1].
Answer. ∫_C F . dr = 1
Worked Example 2
Problem. Show F = (y, x) is conservative and find its potential, then redo Example 1 via endpoints.
Answer. phi = x*y; path integral = phi(1,1) - phi(0,0) = 1
Worked Example 3
Problem. Evaluate the loop integral of the conservative field F = (y, x) around any closed curve.
Answer. 0 (closed-loop integral of a conservative field)
Problem. Evaluate ∫_C F . dr for F = (2x, 2y) along the straight path r(t) = (t, t), t in [0,1], and verify with a potential.
Solution. r'(t) = (1,1); F(r(t)) = (2t, 2t); F . r' = 2t + 2t = 4t. ∫_0^1 4t dt = [2t^2]_0^1 = 2. Check: F = grad(x^2 + y^2), potential phi = x^2+y^2. phi(1,1) - phi(0,0) = 2 - 0 = 2. Matches.
Grad, div, and curl are three uses of the single del operator nabla, mapping between scalar and vector fields.
The del operator nabla = (d/dx, d/dy, d/dz) is a symbolic vector of partial derivatives, and the three core operators are different products with it. Gradient: nabla f turns a scalar field into a vector field pointing uphill. Divergence: nabla . F dots del with a vector field to yield a scalar (spreading). Curl: nabla x F crosses del with a vector field to yield a vector (rotation). Two identities are foundational: curl(grad f) = 0 (gradient fields are irrotational, so every loss-gradient field is curl-free) and div(curl F) = 0. The Laplacian nabla^2 f = div(grad f) = f_xx + f_yy + ... is the divergence of the gradient and appears in diffusion, smoothing, and the regularization terms of many ML models. Knowing which operator maps scalar<->vector keeps derivations type-correct.
Worked Example 1
Problem. For f(x,y) = x^2*y, compute grad f, then div(grad f) (the Laplacian).
Answer. grad f = (2xy, x^2); Laplacian = 2y
Worked Example 2
Problem. Verify curl(grad f) = 0 for f(x,y) = x^2*y in 2D (scalar curl).
Answer. curl(grad f) = 0 (gradient fields are irrotational)
Worked Example 3
Problem. Compute the Laplacian of f(x,y) = sin(x) + y^3.
Answer. Laplacian = -sin(x) + 6*y
Problem. For f(x,y) = e^x*cos(y), compute the Laplacian nabla^2 f.
Solution. f_x = e^x cos y -> f_xx = e^x cos y. f_y = -e^x sin y -> f_yy = -e^x cos y. nabla^2 f = f_xx + f_yy = e^x cos y - e^x cos y = 0. (This f is harmonic — its Laplacian vanishes everywhere.)
Gradient fields, divergence, and curl underpin optimization, diffusion models, physics-informed networks, and more.
Vector calculus is not just physics background; it is active in modern ML. The negative-gradient field -grad L is the conservative vector field that gradient descent follows, with the loss as its potential. Divergence drives continuity equations: continuous normalizing flows and diffusion/score-based generative models track how a probability density changes via the divergence of a learned velocity or score field (the log-density change integrates -div). The Laplacian appears in graph learning (the graph Laplacian), image smoothing, and regularizers. Physics-informed neural networks embed div, curl, and Laplacian operators directly into their loss to enforce PDE constraints like conservation and incompressibility. Recognizing these operators lets you read and design such models. The throughline: scalar potentials and their gradient fields connect every optimization problem to this vector-calculus vocabulary.
Worked Example 1
Problem. Explain why the gradient-descent flow field -grad L is conservative and what its potential is.
Answer. It is grad(-L); potential -L, hence conservative/curl-free
Worked Example 2
Problem. A 2D incompressible flow has div v = 0. Check whether v = (-y, x) qualifies.
Answer. Yes: div v = 0, so it is incompressible
Worked Example 3
Problem. In a continuous normalizing flow, the log-density changes by -div(velocity). For velocity v = (a*x, a*y), find the instantaneous rate of log-density change.
Answer. d/dt log p = -2a
Problem. For a normalizing-flow velocity v(x,y) = (x^2, -y), find div v and the instantaneous log-density rate -div v at the point (1, 5).
Solution. div v = d/dx(x^2) + d/dy(-y) = 2x - 1. At (1,5): 2*1 - 1 = 1. The log-density rate is -div v = -1, so log p is decreasing at rate 1 there (the field is locally expanding, thinning the density).
Define a 2D vector field such as F(x,y) = (-y, x) and a radial field G(x,y) = (x, y) in NumPy. Compute the divergence and curl of each by hand and verify with finite-difference approximations on a grid. Then use Matplotlib quiver to plot both fields and explain which one rotates and which one spreads out.
Deliverable · A script printing hand-computed and numerical divergence/curl for both fields, plus quiver plots with a short note identifying rotation versus radial flow.
1. A vector field assigns to each point in space a...
Answer B. A vector field maps each point to a vector, such as a velocity, force, or flow direction.
2. The divergence of a vector field measures...
Answer B. Divergence is a scalar capturing net outflow per unit volume: positive at sources, negative at sinks.
3. The curl of a vector field measures...
Answer C. Curl quantifies how much the field swirls around a point; it is zero for irrotational fields.
4. A conservative vector field is one that can be written as the...
Answer B. Conservative fields equal grad(phi) for some potential phi, making line integrals path-independent.
5. For a conservative field, the line integral between two points depends on...
Answer B. Path independence means only the start and end potentials matter, not the route between them.
I can describe a vector field and compute its divergence and curl.
I can evaluate a line integral and recognize a conservative field.
I can relate gradient, divergence, and curl to data-science and physics settings.
Before any matrix-calculus derivation, fix whether derivatives come out as rows or columns — the layout convention that keeps shapes consistent.
Matrix calculus has two equally valid conventions for arranging the partials of a vector with respect to a vector. In numerator layout the derivative dy/dx of an m-vector y by an n-vector x is m-by-n (outputs index rows) — this matches the Jacobian. In denominator layout it is n-by-m (the transpose), so the gradient of a scalar comes out as a column. Mixing them is the single biggest source of transpose bugs in ML derivations. The practical fix is to pick one convention, keep it for the whole derivation, and sanity-check shapes at every step: a gradient must have the same shape as the parameter it differentiates. This lesson is the contract that makes the identities in the rest of the unit reliable.
Worked Example 1
Problem. For y in R^3 and x in R^2, give the shape of dy/dx in numerator and denominator layout.
Answer. Numerator: 3x2; Denominator: 2x3
Worked Example 2
Problem. For a scalar loss L and parameter vector w in R^n, what shape must dL/dw be to update w?
Answer. dL/dw is an n-vector (column), matching w
Worked Example 3
Problem. In numerator layout the derivative of a scalar L by a column vector x is a row 1xn. How do you recover the usual gradient column?
Answer. Transpose the 1xn row to get the n x 1 gradient column
Problem. You differentiate a scalar loss L by a weight matrix W of shape 4x3. In any consistent convention, what shape must dL/dW be, and why?
Solution. dL/dW must have the same shape as W, namely 4x3. The reason is the update rule W <- W - eta*(dL/dW): adding/subtracting requires identical shapes, and each entry (dL/dW)[i,j] is the partial of the scalar L with respect to W[i,j]. So the gradient of a scalar by a matrix is a matrix of the same shape.
Two identities — d(a^T x)/dx = a and d(x^T A x)/dx = (A + A^T)x — are the atoms of nearly every ML gradient.
Most loss gradients reduce to two building blocks. The linear form a^T x has gradient a: differentiating a sum a_1 x_1 + ... + a_n x_n with respect to x_j gives a_j, so the whole gradient is the constant vector a. The quadratic form x^T A x has gradient (A + A^T) x, which simplifies to 2 A x when A is symmetric — exactly the Hessian-times-position structure from Unit 1. A quick way to see the quadratic result: x^T A x = sum_{i,j} A_ij x_i x_j, and d/dx_k picks up sum_j A_kj x_j (from i=k) plus sum_i A_ik x_i (from j=k), i.e. row k of A plus column k of A, which is the k-th entry of (A + A^T)x. These two identities, combined with the chain rule, derive least squares, ridge, and many more.
Worked Example 1
Problem. Differentiate f(x) = a^T x with a = (2, -1, 3) with respect to x.
Answer. grad f = a = (2, -1, 3)
Worked Example 2
Problem. For symmetric A = [[2,1],[1,3]], find the gradient of x^T A x at x = (1, 1).
Answer. grad = 2 A x = (6, 8)
Worked Example 3
Problem. For non-symmetric A = [[1,2],[0,3]], find the gradient of x^T A x at x = (1, 1).
Answer. grad = (A + A^T)x = (4, 8)
Problem. For A = [[4,2],[2,1]] (symmetric) and a = (1, -3), find the gradient of f(x) = x^T A x - a^T x at x = (1, 2).
Solution. grad = 2 A x - a (linear and quadratic identities combined). A x = (4*1+2*2, 2*1+1*2) = (8, 4); 2 A x = (16, 8). Subtract a = (1,-3): grad = (16-1, 8-(-3)) = (15, 11). Answer: grad f(1,2) = (15, 11).
Squared error, ridge, and logistic losses have compact gradients built from the linear/quadratic identities and the chain rule.
Training reduces to differentiating a loss with respect to parameters, and a handful of patterns recur. Mean-squared error per example, (y_hat - y)^2 with y_hat = w^T x, has gradient 2(y_hat - y)x by the chain rule (outer 2*residual times inner x). Ridge regularization adds lambda*||w||^2, contributing 2*lambda*w. Logistic regression's binary cross-entropy with sigmoid produces the famously clean gradient (sigma(w^T x) - y)x — the residual times the input, identical in form to linear regression. Recognizing these closed forms lets you implement and check gradients quickly, and reveals why so many models share the 'prediction-minus-target times feature' structure. Each is a direct application of the linear-form derivative composed with an outer scalar derivative via the chain rule.
Worked Example 1
Problem. For the per-example squared loss L = (w^T x - y)^2, derive dL/dw.
Answer. dL/dw = 2*(w^T x - y)*x
Worked Example 2
Problem. Add ridge: L = (w^T x - y)^2 + lambda*||w||^2. Give dL/dw.
Answer. dL/dw = 2*(w^T x - y)*x + 2*lambda*w
Worked Example 3
Problem. For logistic regression L = -[y*log(p) + (1-y)*log(1-p)] with p = sigma(w^T x), show dL/dw = (p - y)x.
Answer. dL/dw = (p - y)*x
Problem. For L = (w^T x - y)^2 with w = (1, 2), x = (3, 1), y = 4, compute the numeric gradient dL/dw.
Solution. Prediction w^T x = 1*3 + 2*1 = 5. Residual r = 5 - 4 = 1. dL/dw = 2*r*x = 2*1*(3,1) = (6, 2). Answer: dL/dw = (6, 2).
Composed matrix operations differentiate by multiplying Jacobians, with transposes placed so shapes align.
The matrix chain rule is the multivariable chain rule with explicit attention to layout. For a composition where a scalar loss L depends on z = W x through z, the gradient with respect to W uses dL/dW = (dL/dz) x^T and the gradient passed back is dL/dx = W^T (dL/dz) — the same outer-product and transpose patterns seen in Unit 4's layer derivation, now stated as identities. The governing principle: each link contributes its local Jacobian, and the transposes are forced by requiring the result to match the differentiated object's shape. For elementwise maps the 'Jacobian' is diagonal, so its action is an elementwise (Hadamard) product, not a matrix multiply. Mastering where transposes go turns any layered computation into a mechanical gradient derivation.
Worked Example 1
Problem. For z = W x and scalar L, given g = dL/dz, derive dL/dx and dL/dW.
Answer. dL/dx = W^T g, dL/dW = g x^T
Worked Example 2
Problem. For a = sigma(z) elementwise with z = W x, given dL/da, find dL/dz.
Answer. dL/dz = dL/da ⊙ sigma'(z) (elementwise)
Worked Example 3
Problem. Numerically: W = [[1,2],[3,4]], x = (1,1), g = dL/dz = (1, 0). Compute dL/dW and dL/dx.
Answer. dL/dW = [[1,1],[0,0]], dL/dx = (1, 2)
Problem. For z = W x with W = [[2,0],[1,3]], x = (1, 2), and upstream g = dL/dz = (1, 1), compute dL/dW and dL/dx.
Solution. dL/dW = g x^T = (1,1)^T (1,2) = [[1*1, 1*2],[1*1, 1*2]] = [[1,2],[1,2]] (shape of W, 2x2). dL/dx = W^T g = [[2,1],[0,3]] (1,1) = (2+1, 0+3) = (3, 3). Answer: dL/dW = [[1,2],[1,2]], dL/dx = (3, 3).
Softmax followed by cross-entropy yields the elegant gradient p - y, the cornerstone of classifier training.
Softmax turns a score vector z into probabilities p_i = exp(z_i)/sum_j exp(z_j). Its Jacobian has entries dp_i/dz_j = p_i(delta_ij - p_j) — diagonal terms p_i(1-p_i), off-diagonal -p_i p_j. Cross-entropy loss for a one-hot label y is L = -sum_i y_i log(p_i). The remarkable result, after composing the cross-entropy derivative with the softmax Jacobian, is that everything telescopes to dL/dz = p - y: the predicted distribution minus the true one. This is why frameworks fuse softmax and cross-entropy into one op — it is numerically stabler and the gradient is trivially p - y. Knowing this derivation explains the universal 'prediction minus target' gradient form and lets you implement a classifier's backward pass in a single line.
Worked Example 1
Problem. Compute softmax of z = (1, 2, 3).
Answer. p ≈ (0.090, 0.245, 0.665)
Worked Example 2
Problem. With p ≈ (0.090, 0.245, 0.665) and one-hot label y = (0, 0, 1), compute dL/dz.
Answer. dL/dz ≈ (0.090, 0.245, -0.335)
Worked Example 3
Problem. Verify the softmax Jacobian diagonal entry for class i: dp_i/dz_i = p_i(1 - p_i) using p_1 ≈ 0.090.
Answer. dp_1/dz_1 ≈ 0.082
Problem. For logits z = (2, 0, 0) and true class 0 (y = (1,0,0)), compute the softmax p and the gradient dL/dz = p - y.
Solution. exp(z) = (e^2, 1, 1) ≈ (7.389, 1, 1), sum ≈ 9.389. p ≈ (0.787, 0.1065, 0.1065). Gradient dL/dz = p - y = (0.787-1, 0.1065-0, 0.1065-0) ≈ (-0.213, 0.107, 0.107). The negative entry on the true class pulls its logit up, as expected.
A compact table of the recurring identities lets you assemble any ML gradient by composition.
Rather than re-deriving from scratch, practitioners keep a cheat sheet of vetted identities and combine them with the chain rule. The core entries (numerator-consistent): d(a^T x)/dx = a^T; d(x^T a)/dx = a^T; d(x^T A x)/dx = x^T(A + A^T), or 2x^T A for symmetric A; d(A x)/dx = A; d(||x||^2)/dx = 2x^T; d(softmax+CE)/dz = (p - y)^T. To use the sheet, decompose a loss into these primitives, differentiate each, and chain them with correct transposes, checking shapes at every step. This lesson is meta-skill: the value is not memorizing but knowing the small set of patterns that compose into every gradient you will meet, and always verifying the assembled result numerically against autograd before trusting it.
Worked Example 1
Problem. Use the cheat sheet to differentiate f(x) = a^T x + x^T A x (symmetric A) with respect to x.
Answer. grad f = a + 2 A x
Worked Example 2
Problem. Differentiate the regularized objective f(w) = ||w||^2 with respect to w, and state the identity used.
Answer. grad f = 2 w
Worked Example 3
Problem. Assemble the gradient of f(x) = (A x)^T (A x) using cheat-sheet identities.
Answer. grad f = 2 A^T A x
Problem. Use cheat-sheet identities to find the gradient of f(w) = b^T w + ||w||^2 at w = (1, 2) with b = (3, -1).
Solution. d(b^T w)/dw = b = (3, -1). d(||w||^2)/dw = 2 w = (2, 4). Sum: grad f = b + 2w = (3+2, -1+4) = (5, 3). Answer: grad f(1,2) = (5, 3).
The full least-squares gradient 2 X^T (X w - y) follows in two lines from the chain rule and quadratic-form identity.
Linear regression minimizes L(w) = ||X w - y||^2 over weights w, with X the n-by-d design matrix. Write the residual r = X w - y, so L = r^T r = ||r||^2. By the squared-norm identity dL/dr = 2 r, and since r is linear in w with Jacobian X (dr/dw = X), the matrix chain rule gives dL/dw = X^T (2 r) = 2 X^T (X w - y). Setting this to zero yields the normal equations X^T X w = X^T y. This single derivation showcases the whole unit: a quadratic objective, the chain rule, transposes placed for shape consistency, and a result you verify against autograd. Adding ridge lambda||w||^2 simply appends 2 lambda w, giving the regularized gradient 2 X^T(Xw - y) + 2 lambda w.
Worked Example 1
Problem. Derive dL/dw for L = ||X w - y||^2 step by step.
Answer. dL/dw = 2 X^T (X w - y)
Worked Example 2
Problem. Set the gradient to zero to obtain the normal equations.
Answer. X^T X w = X^T y
Worked Example 3
Problem. Compute the numeric gradient for X = [[1,0],[0,1],[1,1]], y = (1, 2, 2), w = (0, 0).
Answer. dL/dw = (-6, -8)
Problem. For X = [[1,2],[3,4]], y = (1, 0), w = (1, 1), compute the least-squares gradient 2 X^T (X w - y).
Solution. X w = (1*1+2*1, 3*1+4*1) = (3, 7). Residual r = X w - y = (3-1, 7-0) = (2, 7). X^T r = (1*2 + 3*7, 2*2 + 4*7) = (2+21, 4+28) = (23, 32). Gradient = 2*(23, 32) = (46, 64). Answer: dL/dw = (46, 64).
Using matrix-calculus identities, derive the gradient of the least-squares loss L(w) = ||X w - y||^2 with respect to w, arriving at 2 X^T (X w - y). Implement L and its analytic gradient in NumPy, generate random X, y, and w, and confirm the analytic gradient matches autograd to within 1e-6. Then derive and verify the gradient of softmax cross-entropy for a single example.
Deliverable · A script printing the analytic least-squares gradient versus autograd with a passing assertion, plus the softmax cross-entropy gradient verified the same way.
1. Using the standard identity, the derivative of the linear form a^T x with respect to x is:
Answer B. Differentiating a^T x with respect to x yields the constant vector a.
2. For a symmetric matrix A, the gradient of the quadratic form x^T A x with respect to x is:
Answer C. In general (A + A^T) x; for symmetric A this is 2 A x.
3. The gradient of the least-squares loss ||X w - y||^2 with respect to w is:
Answer B. Applying the chain rule to the squared residual gives 2 X^T (X w - y).
4. Why must you fix a layout convention (numerator vs denominator) before doing matrix calculus?
Answer B. Layout decides the orientation of derivatives; mixing conventions produces transposed, mismatched shapes.
5. Softmax is used in classification to...
Answer B. Softmax exponentiates and normalizes scores into nonnegative values summing to one: a distribution.
I can apply matrix-calculus identities for linear and quadratic forms with a consistent layout.
I can derive the gradient of a common loss such as least squares or cross-entropy.
I can verify a matrix-calculus derivation numerically with autograd.
Where this leads
Part of the open-source Crunch Academy tech-education path. Pair it with the data-science and ML tracks when you're ready.
All Math for Data Science courses Crunch Academy home