Multivariable, Vector and Matrix Calculus

Introduction

Data-science models rarely depend on one variable. A prediction may depend simultaneously on many features, parameters, embeddings, or intermediate representations. Multivariable calculus provides the language for describing how such quantities change. Vector and matrix calculus extend this language to the dimensions encountered in regression, neural networks, optimization, probabilistic modelling, and scientific computing.

The central objects are partial derivatives, gradients, directional derivatives, Jacobians, Hessians, and differentials. Together with the chain rule, they support sensitivity analysis, optimization, backpropagation, and automatic differentiation. Constrained optimization connects these ideas to practical problems such as resource allocation, portfolio construction, and regularized learning.

Learning Objectives

By the end of this chapter, you should be able to:

1. Compute and interpret partial derivatives and gradients.
2. Use directional derivatives to measure change along a specified direction.
3. Construct Jacobian and Hessian matrices and explain their dimensions.
4. Apply the multivariable chain rule to composite functions and computational graphs.
5. Solve basic unconstrained and constrained-extrema problems.
6. Derive common matrix-calculus results used in machine learning.
7. Explain reverse-mode automatic differentiation and its relation to backpropagation.
8. Identify numerical, modelling, and implementation limitations of derivative-based methods.

1. Functions of Several Variables

A scalar-valued function of n variables is a mapping

f:ℝⁿ → ℝ,

written as f(x₁, …, x_n). A vector-valued function maps an input vector to an output vector:

F:ℝⁿ → ℝ^m.

The input x = (x₁, …, x_n)^T is conventionally a column vector. A function is differentiable at a point if its local change can be approximated by a linear map. Specifically, for a small increment ∆x,

Here ∇f(x) is the gradient, ‖∆x‖ is Euclidean length, and the remainder term becomes negligible relative to ‖∆x‖ as the increment approaches zero. Differentiability is stronger than the mere existence of separate partial derivatives.

Partial derivatives

The partial derivative with respect to x_i holds all other variables constant:

where e_i is the i-th standard basis vector. Partial derivatives describe coordinate-wise sensitivity, but they should not be interpreted as the complete rate of change unless a direction is specified.

Example 1: A quadratic surface

Let

Then

At (x,y)=(1,2), the partial derivatives are 10 and 6. Increasing x alone by a small amount changes f approximately 10 times that amount; increasing y alone changes it approximately 6 times that amount.

2. Gradients and Directional Derivatives

For a differentiable scalar function, the gradient is the column vector of partial derivatives:

The gradient points in the direction of greatest local increase. Its magnitude is the maximum rate of increase per unit Euclidean distance.

For a nonzero direction vector v, define the unit direction u = v/‖v‖. The directional derivative is

If the direction is not normalized, ∇f(x)^T v gives the rate associated with the parameterization by v, not the rate per unit distance. This distinction is important in optimization and sensitivity analysis.

For Example 1, at (1,2), ∇f=(10,6)^T. In the direction v=(3,4)^T, the unit vector is (3/5,4/5), so

The steepest ascent direction is ∇f/‖∇f‖, while the steepest descent direction is its negative.

Geometric interpretation

A level set is a set on which f(x)=c. At a regular point where the gradient is nonzero, the gradient is perpendicular to the level set. This fact supports contour-based visualization and explains why gradient descent crosses level curves orthogonally in the infinitesimal limit.

3. Jacobians and Differentials

For F=(F₁,…,F_m)^T:ℝⁿ→ℝ^m, the Jacobian is the m×n matrix

Its rows are gradients of output components, and its columns describe how all outputs change when one input changes. The first-order approximation is

The dimensions must agree: an m×n Jacobian multiplying an n×1 increment produces an m×1 output increment.

Example 2: A nonlinear transformation

Let

Then

At (1,0), the Jacobian is [[0,1],[1,1]]. Consequently, a small change (∆x,∆y)^T produces approximately (∆y,∆x+∆y)^T.

The differential notation summarizes the same relationship:

For a scalar function, df = ∇f^Tdx. Differentials are especially useful for deriving matrix-calculus identities.

4. The Multivariable Chain Rule

Suppose y=F(x) and z=G(y), with compatible dimensions. If both maps are differentiable, then

The order matters. The inner Jacobian acts first, so it appears on the right.

For scalar functions, if z=g(u,v), u=u(x,y), and v=v(x,y), then

Example 3: A composed model

Let s=w^Tx+b and L=½(s−t)², where w∈ℝ^d, b and t are scalars. Since ∂L/∂s=s−t and ∂s/∂w=w? is not dimensionally meaningful as written, use differentials:

Therefore,

which gives

This derivation illustrates why differentials help prevent orientation errors.

5. Hessians and Second-Order Information

For a twice-differentiable scalar function, the Hessian is the matrix of second partial derivatives:

Under continuity of the second partial derivatives near the point, Clairaut’s theorem implies symmetry: ∂²f/∂x_i∂x_j = ∂²f/∂x_j∂x_i. A second-order Taylor approximation is

The quadratic form v^TH_fv describes curvature along direction v. A positive-definite Hessian has positive quadratic form for every nonzero v and indicates local convex curvature.

Classification of stationary points

A stationary point satisfies ∇f(x*)=0. If the Hessian is positive definite at the point, it is a strict local minimum; if negative definite, it is a strict local maximum. If it is indefinite, the point is a saddle. A singular or semidefinite Hessian gives an inconclusive second-order test and requires higher-order or other analysis.

For Example 1,

Its leading principal minors are 6 and 6×2−2×2=8, both positive, so H is positive definite. The unique stationary point solves 6x+2y=0 and 2x+2y=0, giving (0,0), a strict global minimum because the function is a positive-definite quadratic.

6. Extrema and Constrained Optimization

An unconstrained optimization problem has the form

A necessary condition for an interior differentiable local optimum is ∇f(x*)=0. It is not sufficient by itself; saddle points also have zero gradient.

A constrained problem restricts x to a feasible set. For an equality constraint g(x)=c, define the Lagrangian

At a regular constrained stationary point, where ∇g(x*)≠0, the Lagrange conditions are

They state that the objective gradient is normal to the constraint surface. Multiple equality constraints use ℒ=f+∑_k λ_k(g_k−c_k).

Worked example: constrained minimum

Minimize f(x,y)=x²+y² subject to x+y=1. The Lagrangian is

The stationary equations are 2x+λ=0, 2y+λ=0, and x+y=1. Thus x=y, and the constraint gives x=y=1/2. The minimum value is 1/2. Geometrically, this is the point on the line x+y=1 closest to the origin.

For inequality constraints, first-order optimality is expressed by the Karush–Kuhn–Tucker conditions. For a minimization problem with g_i(x)≤0, multipliers satisfy λ_i≥0, stationarity is ∇f+∑λ_i∇g_i=0, and complementary slackness is λ_i g_i(x)=0. Constraint qualification assumptions are needed for these conditions to be necessary.

7. Matrix Calculus for Data Science

Matrix calculus is not a single universal notation; numerator-layout and denominator-layout conventions differ. This chapter uses column gradients for scalar functions: ∇_x f is a column vector. Consistency is more important than a particular convention.

Useful identities, for constant compatible matrices and vectors, include:

The last identity follows from r=Ax−b and dr=A dx:

so the column gradient is 2A^Tr. Dimensions provide a useful check: if A is m×n, then A^T(Ax−b) is n×1, matching x.

Least squares and ridge regression

Given design matrix X∈ℝ^m×p, response y∈ℝ^m, and parameter β∈ℝ^p, the half-scaled least-squares objective is

Its gradient and Hessian are

If X^TX is invertible, the stationary solution is β̂=(X^TX)^-1X^Ty. In practice, solving the linear system is preferable to explicitly computing an inverse. With ridge penalty (λ/2)‖β‖², where λ>0,

The penalty improves conditioning and controls coefficient magnitude, but changes the estimand and introduces bias.

8. Optimization and Numerical Considerations

Gradient descent updates parameters by

where η_k>0 is the learning rate. For a quadratic function with Hessian eigenvalues between μ and L, fixed-step convergence is guaranteed under standard assumptions when 0<η<2/L; strong convexity, represented by μ>0, gives stronger convergence-rate results. Poor conditioning—large L/μ—causes slow progress and elongated optimization contours.

Newton’s method uses curvature:

provided the Hessian is nonsingular and the local model is appropriate. Implementations normally solve a linear system rather than form the inverse. Newton steps can be unstable for nonconvex objectives, indefinite Hessians, or poorly scaled data. Quasi-Newton methods approximate curvature, while stochastic gradient methods use minibatch estimates and trade exactness for scalability.

Finite differences can check derivatives. A central approximation is

Too large an h creates truncation error; too small an h creates cancellation and floating-point error. Finite differences are validation tools, not usually the preferred training mechanism.

9. Automatic Differentiation and Backpropagation

Automatic differentiation (AD) evaluates derivatives by decomposing a program into elementary operations and applying the chain rule. It is distinct from symbolic differentiation, which manipulates formulas, and numerical differentiation, which estimates derivatives from perturbed function values.

In forward mode, derivatives propagate with primal values from inputs toward outputs. It is efficient when the number of inputs is small relative to outputs. In reverse mode, a forward pass records intermediate values and a reverse pass propagates adjoints from outputs toward inputs. It is efficient for a scalar loss depending on many parameters, which is why it underlies neural-network backpropagation.

Consider a computational graph:

Figure: A computational graph and the direction of reverse-mode derivative propagation.

For L=½(s−t)², the reverse pass begins with ∂L/∂s=s−t. It then gives ∇_wL=(s−t)x and ∂L/∂b=s−t. Modern AD systems also support vector-Jacobian products and Jacobian-vector products without necessarily materializing a full Jacobian.

AD is exact up to floating-point arithmetic, but it does not remove numerical problems. Overflow, underflow, saturated activations, ill-conditioned operations, discontinuities, and memory requirements can still produce unreliable gradients. Gradient clipping, normalization, stable loss formulas, and careful initialization are engineering responses, not substitutes for mathematical understanding.

10. Applied Case: Training a Regularized Linear Predictor

Background and problem

A data team wants to predict a standardized continuous outcome from p correlated measurements. Correlation among predictors makes ordinary least squares unstable, and the team requires a reproducible, interpretable baseline.

Analysis

Let X contain the centered and scaled features, y the centered response, and β the coefficient vector. The team minimizes

with λ>0. Differentiation gives

Setting the gradient to zero produces the ridge system

Approach and solution

The team selects λ using cross-validation performed only on the training portion of the data. It solves the linear system using a numerically stable factorization, rather than explicitly calculating an inverse. A derivative check on a small synthetic dataset compares the analytic gradient with central finite differences.

Outcome

The ridge model gives more stable coefficients under resampling and improves validation error relative to the unregularized baseline. However, the selected λ and reported performance remain dependent on preprocessing, data splitting, and the possibility of distribution shift.

Learning points

1. Matrix dimensions expose errors before software execution.
2. Regularization changes both optimization curvature and statistical interpretation.
3. Validation must prevent leakage from test data.
4. A correct gradient does not guarantee a valid model or fair decision.

11. Limitations, Governance, and Best Practice

Derivative-based methods assume a meaningful local approximation. Nondifferentiable functions can sometimes be handled with subgradients, but ordinary Hessian reasoning may fail. Deep models are often nonconvex; a stationary point need not be globally optimal. Scaling features can substantially improve conditioning and optimization.

In high-impact applications, gradients and sensitivities should not be treated automatically as causal effects. A large derivative may reflect correlation, a chosen scale, or a local region unlike the deployment population. Models should document feature transformations, parameter choices, derivative checks, numerical precision, and failure modes. Sensitive attributes and their proxies require appropriate privacy, fairness, and governance review. Automatic differentiation can make a computation efficient, but it does not make the objective ethically appropriate or the data representative.

Best practice includes: state variable orientations explicitly; check dimensions; use stable linear algebra; standardize when justified; validate gradients on small cases; monitor exploding or vanishing gradients; distinguish local from global claims; and report assumptions alongside numerical results.

Chapter Summary

Multivariable calculus generalizes one-variable differentiation to functions of vectors and matrices. Partial derivatives describe coordinate-wise sensitivity, while the gradient gives the direction and magnitude of greatest local increase. Directional derivatives measure change along a normalized direction. The Jacobian is the linear approximation of a vector-valued map, and the Hessian records second-order curvature of a scalar function. The multivariable chain rule is the foundation of computational-graph differentiation and neural-network backpropagation.

Stationary-point tests use gradients and Hessians, while Lagrange multipliers and KKT conditions handle constraints under suitable regularity assumptions. Matrix calculus produces compact, dimensionally checkable derivatives for least squares, ridge regression, and many machine-learning objectives. Automatic differentiation applies the chain rule computationally and efficiently, especially in reverse mode for scalar losses. Nevertheless, numerical stability, conditioning, nonconvexity, data quality, fairness, and interpretation remain essential concerns.

Key Terms

Partial derivative: Derivative with respect to one variable while holding the others fixed.

Gradient: Column vector of first partial derivatives of a scalar function.

Directional derivative: Rate of change in a specified direction, usually normalized.

Level set: Set of points at which a function has a constant value.

Differentiability: Existence of a valid first-order linear approximation.

Jacobian: Matrix of first derivatives of a vector-valued function.

Differential: Linear first-order change associated with an infinitesimal increment.

Chain rule: Rule for differentiating a composition of functions.

Hessian: Matrix of second partial derivatives of a scalar function.

Positive definite: A matrix A satisfying v^TAv>0 for every nonzero v.

Stationary point: Point at which the gradient is zero.

Lagrangian: Objective augmented by weighted constraint functions.

KKT conditions: First-order conditions for constrained optimization, including inequalities.

Conditioning: Sensitivity of a numerical problem to perturbations.

Gradient descent: Iterative method that moves opposite the gradient.

Newton’s method: Optimization method using gradient and Hessian curvature.

Automatic differentiation: Chain-rule evaluation of derivatives through elementary operations.

Adjoint: Reverse-mode sensitivity propagated from an output toward inputs.

Backpropagation: Reverse-mode differentiation applied to layered computational models.

Regularization: Addition of a penalty or constraint to control complexity or ill-conditioning.

Review Questions

Short Answer Questions

1. Define the gradient of a scalar function and state its geometric meaning.
2. What are the dimensions of the Jacobian for a map from ℝⁿ to ℝ^m?
3. Why must a direction vector be normalized when computing a rate per unit distance?
4. State the second-order classification result for a positive-definite Hessian at a stationary point.
5. Distinguish automatic differentiation from finite-difference differentiation.

Long Answer Questions

1. Derive the gradient and Hessian of f(x)=½‖Ax−b‖², defining all dimensions and assumptions.
2. Explain the multivariable chain rule using Jacobian multiplication and discuss why matrix order matters.
3. Derive the ridge-regression estimator from its penalized objective and explain the effect of λ.
4. Compare forward-mode and reverse-mode automatic differentiation, including their computational trade-offs.
5. Explain Lagrange multipliers geometrically and solve a general quadratic-distance problem subject to one linear equality constraint.

Analytical / Application Questions

1. A model has loss L(w)=log(1+e^wTx). Derive ∇_wL and identify a numerically stable implementation concern.
2. For a twice-differentiable nonconvex objective with a singular Hessian at a stationary point, explain why the usual second-order test is inconclusive and propose further analyses.
3. Design a derivative-validation protocol for a neural-network loss using automatic differentiation and finite differences. Discuss step-size selection, floating-point error, and what the test cannot establish.

Multiple Choice Questions with Answers

Question 1

For a differentiable scalar function f, which expression gives the directional derivative in unit direction u?

A. ∇f(x)^Tu

B. f(x)u

C. ∇²f(x)u

D. ‖∇f(x)‖²u

Correct Answer: A. ∇f(x)^Tu

Question 2

If F maps ℝⁿ to ℝ^m, the Jacobian of F has which dimension?

A. n×m

B. m×n

C. m×m

D. n×n

Correct Answer: B. m×n

Question 3

For f(x,y)=x²+3xy, what is ∂f/∂y?

A. 2x+3y

B. x²+3x

C. 3x

D. 2y+3x

Correct Answer: C. 3x

Question 4

Under standard differentiability assumptions, the Jacobian of a composition G∘F is:

A. J_FJ_G

B. J_G+J_F

C. J_G(F(x))J_F(x)

D. J_G(F(x))−J_F(x)

Correct Answer: C. J_G(F(x))J_F(x)

Question 5

If the Hessian at a stationary point is positive definite, the point is a:

A. Strict local minimum

B. Strict local maximum

C. Saddle point

D. Point that must be globally maximal

Correct Answer: A. Strict local minimum

Question 6

For J(β)=½‖Xβ−y‖², the gradient with respect to β is:

A. Xβ−y

B. 2X(β−y)

C. X^T(Xβ−y)

D. XX^Tβ−y

Correct Answer: C. X^T(Xβ−y)

Question 7

In reverse-mode automatic differentiation, sensitivities generally propagate:

A. From inputs toward outputs only

B. From outputs toward inputs

C. Only between independent inputs

D. Without using the chain rule

Correct Answer: B. From outputs toward inputs

Question 8

For the constraint g(x)=c, the Lagrange first-order condition is:

A. ∇f=0 regardless of the constraint

B. ∇f+λ∇g=0 together with g=c

C. f+g=0

D. ∇f×∇g=0 only

Correct Answer: B. ∇f+λ∇g=0 together with g=c

Question 9

If A is m×n, what is the dimension of A^T(Ax−b) when x is n×1 and b is m×1?

A. m×1

B. n×1

C. m×n

D. n×m

Correct Answer: B. n×1

Question 10

Which statement about finite-difference derivatives is correct?

A. They are unaffected by floating-point precision.

B. Making the step size arbitrarily small always improves accuracy.

C. They estimate derivatives using function evaluations at perturbed inputs.

D. They require a Hessian to compute a first derivative.

Correct Answer: C. They estimate derivatives using function evaluations at perturbed inputs.