In machine learning, parameter estimation, also known as model fitting, is formulated as an
optimization problem:
\[
\begin{align*}
\boldsymbol{\theta}^* &\in \operatorname*{arg\,min}_{\boldsymbol{\theta} \in \Theta} \mathcal{L}(\boldsymbol{\theta}), \\\\
\Theta &\subseteq \mathbb{R}^d,
\end{align*}
\]
where \(\mathcal{L}\) is the loss function (or objective function), \(\Theta\) is the
parameter space, and \(d\) is the number of parameters. Since a global minimizer is in general NP-hard to
certify, our practical target is a local minimum.
Definition: Local Minimum
A point \(\boldsymbol{\theta}^* \in \Theta\) is a local minimum of \(\mathcal{L}\) if
\[
\exists\, \delta \gt 0 \text{ such that } \mathcal{L}(\boldsymbol{\theta}^*) \leq \mathcal{L}(\boldsymbol{\theta})
\quad \text{for all } \boldsymbol{\theta} \in \Theta \text{ with } \|\boldsymbol{\theta} - \boldsymbol{\theta}^*\| \lt \delta.
\]
If the inequality is strict for \(\boldsymbol{\theta} \neq \boldsymbol{\theta}^*\), we call
\(\boldsymbol{\theta}^*\) a strict local minimum.
For a twice continuously differentiable objective, local minima can be characterized analytically. Write
\(\mathbf{g}(\boldsymbol{\theta}) = \nabla \mathcal{L}(\boldsymbol{\theta})\) for the
gradient
(the column-vector case of the matrix definition) and
\(\mathbf{H}(\boldsymbol{\theta}) = \nabla^2 \mathcal{L}(\boldsymbol{\theta})\) for the
Hessian, a symmetric matrix because continuous mixed partials commute, a fact we
do not prove on this page.
Theorem: Second-Order Optimality Conditions
Let \(\mathcal{L}\) be \(C^2\) on the interior of \(\Theta\) and let \(\boldsymbol{\theta}^*\)
lie in the interior of \(\Theta\).
(Necessary) If \(\boldsymbol{\theta}^*\) is a local minimum, then
\(\mathbf{g}(\boldsymbol{\theta}^*) = \mathbf{0}\) and \(\mathbf{H}(\boldsymbol{\theta}^*) \succeq 0\)
(positive semidefinite).
(Sufficient) If \(\mathbf{g}(\boldsymbol{\theta}^*) = \mathbf{0}\) and
\(\mathbf{H}(\boldsymbol{\theta}^*) \succ 0\) (positive definite),
then \(\boldsymbol{\theta}^*\) is a strict local minimum.
Proof sketch:
Both directions rest on the second-order Taylor expansion, obtained from the \(k = 1\) case of
the
multivariate Taylor theorem
by using continuity of the second partials to replace each integral remainder
\(\int_0^1 (1-t)\,\partial_{ij}\mathcal{L}(\boldsymbol{\theta}^* + t\mathbf{h})\,dt\) by
\(\tfrac{1}{2}\partial_{ij}\mathcal{L}(\boldsymbol{\theta}^*) + o(1)\),
\[
\begin{align*}
\mathcal{L}(\boldsymbol{\theta}^* + \mathbf{h})
&= \mathcal{L}(\boldsymbol{\theta}^*)
+ \mathbf{g}(\boldsymbol{\theta}^*)^\top \mathbf{h} \\\\
&\quad + \tfrac{1}{2}\, \mathbf{h}^\top \mathbf{H}(\boldsymbol{\theta}^*) \mathbf{h}
+ o(\|\mathbf{h}\|^2)
\quad \text{as } \|\mathbf{h}\| \to 0.
\end{align*}
\]
(1) Necessity. Suppose \(\boldsymbol{\theta}^*\) is a local minimum. If
\(\mathbf{g}(\boldsymbol{\theta}^*) \neq \mathbf{0}\), taking
\(\mathbf{h} = -\varepsilon\, \mathbf{g}(\boldsymbol{\theta}^*)\) gives
\(\mathcal{L}(\boldsymbol{\theta}^* + \mathbf{h}) - \mathcal{L}(\boldsymbol{\theta}^*)
= -\varepsilon \|\mathbf{g}(\boldsymbol{\theta}^*)\|^2 + O(\varepsilon^2)\), which is negative for small
\(\varepsilon \gt 0\). This contradicts the local-minimum property, so
\(\mathbf{g}(\boldsymbol{\theta}^*) = \mathbf{0}\). With this, for any unit vector \(\mathbf{v}\) and
small \(\varepsilon \gt 0\),
\(\mathcal{L}(\boldsymbol{\theta}^* + \varepsilon \mathbf{v}) - \mathcal{L}(\boldsymbol{\theta}^*)
= \tfrac{\varepsilon^2}{2}\, \mathbf{v}^\top \mathbf{H}(\boldsymbol{\theta}^*) \mathbf{v} + o(\varepsilon^2) \geq 0\)
forces \(\mathbf{v}^\top \mathbf{H}(\boldsymbol{\theta}^*) \mathbf{v} \geq 0\) for all \(\mathbf{v}\),
that is, \(\mathbf{H}(\boldsymbol{\theta}^*) \succeq 0\).
(2) Sufficiency. Suppose \(\mathbf{g}(\boldsymbol{\theta}^*) = \mathbf{0}\) and
\(\mathbf{H}(\boldsymbol{\theta}^*) \succ 0\). By positive-definiteness,
\(\lambda_{\min}(\mathbf{H}(\boldsymbol{\theta}^*)) \gt 0\). Then for all sufficiently small
\(\mathbf{h} \neq \mathbf{0}\),
\(\tfrac{1}{2}\mathbf{h}^\top \mathbf{H}(\boldsymbol{\theta}^*) \mathbf{h} + o(\|\mathbf{h}\|^2)
\geq \tfrac{1}{2} \lambda_{\min}\, \|\mathbf{h}\|^2 + o(\|\mathbf{h}\|^2) \gt 0\), so
\(\boldsymbol{\theta}^*\) is a strict local minimum.
Convex vs. non-convex landscapes
The sufficient condition (\(\mathbf{H} \succ 0\)) certifies only a strict local minimum. A
non-convex objective typically has many stationary points (points where
\(\mathbf{g} = \mathbf{0}\)): local minima (where necessarily \(\mathbf{H} \succeq 0\)), local maxima
(where necessarily \(\mathbf{H} \preceq 0\)), and saddle points, which are neither and
include every stationary point where \(\mathbf{H}\) is indefinite. First-order methods can converge to
stationary points of any type and cannot distinguish them without further information. This is one
reason the convexity structure studied in the next section is so valuable. Under convexity, stationarity
already certifies global optimality.
Convexity
We often design training objectives to be convex, because in a convex optimization problem every
local minimum is automatically a global minimum. Convexity therefore removes the principal obstacle of non-convex
optimization. The underlying geometry is that of a
convex set: a subset
\(\mathcal{S} \subseteq \mathbb{R}^n\) that contains the line segment joining any two of its points.
Definition: Convex Function
Let \(\mathcal{S} \subseteq \mathbb{R}^n\) be a convex set. A function \(f: \mathcal{S} \to \mathbb{R}\) is a
convex function if for every \(\mathbf{x}, \mathbf{y} \in \mathcal{S}\) and every
\(\lambda \in [0,1]\),
\[
f\!\left(\lambda \mathbf{x} + (1-\lambda)\mathbf{y}\right)
\leq \lambda f(\mathbf{x}) + (1-\lambda) f(\mathbf{y}).
\]
It is strictly convex if the inequality is strict whenever \(\mathbf{x} \neq \mathbf{y}\)
and \(\lambda \in (0,1)\).
In one dimension, for twice-differentiable functions, convexity is captured entirely by the sign of the second
derivative. We begin with a lemma that links convexity to the first-order Taylor inequality. This lemma is the
technical bridge used in both directions of the next theorem.
Lemma: First-Order Characterization of Convexity
Let \(f: I \to \mathbb{R}\) be differentiable on an open interval \(I \subseteq \mathbb{R}\). Then \(f\) is
convex on \(I\) if and only if
\[
f(y) \geq f(x) + f'(x)(y - x) \quad \text{for all } x, y \in I. \tag{*}
\]
Proof:
(\(\Rightarrow\)) Assume \(f\) is convex. Fix \(x, y \in I\) and \(\lambda \in (0,1)\), and
let \(z_\lambda = \lambda y + (1-\lambda)x = x + \lambda(y-x)\). By convexity,
\[
\begin{align*}
f(z_\lambda) &\leq \lambda f(y) + (1-\lambda) f(x), \\\\
\text{that is,} \quad \frac{f(x + \lambda(y-x)) - f(x)}{\lambda} &\leq f(y) - f(x).
\end{align*}
\]
Letting \(\lambda \downarrow 0\) and using differentiability gives \(f'(x)(y-x) \leq f(y) - f(x)\), which is
(*).
(\(\Leftarrow\)) Assume (*) holds. Fix \(a, b \in I\), \(\lambda \in [0,1]\), and set
\(z = \lambda a + (1-\lambda) b\). Applying (*) at the base point \(z\) to both \(a\) and \(b\):
\[
\begin{align*}
f(a) &\geq f(z) + f'(z)(a - z), \\\\
f(b) &\geq f(z) + f'(z)(b - z).
\end{align*}
\]
Multiplying the first inequality by \(\lambda\), the second by \((1-\lambda)\), and adding:
\[
\begin{align*}
\lambda f(a) + (1-\lambda) f(b)
&\geq f(z) + f'(z)\bigl[\lambda (a-z) + (1-\lambda)(b-z)\bigr].
\end{align*}
\]
The bracketed term equals \(\lambda a + (1-\lambda) b - z = 0\) by the definition of \(z\), so
\(\lambda f(a) + (1-\lambda) f(b) \geq f(z)\), which is convexity.
Theorem: Second-Order Condition for Convexity (Univariate)
Let \(f: I \to \mathbb{R}\) be twice differentiable on an open interval \(I\). Then \(f\) is convex on
\(I\) if and only if \(f''(x) \geq 0\) for all \(x \in I\).
Proof:
(\(\Leftarrow\)) Suppose \(f'' \geq 0\) on \(I\). By
Taylor's theorem, Lagrange form,
for any \(x, y \in I\) with
\(x \neq y\) (the case \(x = y\) is trivial) there exists \(c\) between \(x\) and \(y\) such that
\[
\begin{align*}
f(y) &= f(x) + f'(x)(y - x) + \tfrac{1}{2} f''(c)(y - x)^2 \\\\
&\geq f(x) + f'(x)(y - x),
\end{align*}
\]
since the quadratic remainder is non-negative. This is (*) of the preceding lemma, hence \(f\) is convex.
(\(\Rightarrow\)) Suppose \(f\) is convex. By the lemma, (*) holds. Fix \(a \lt b\) in \(I\)
and apply (*) in both directions:
\[
\begin{align*}
f(b) &\geq f(a) + f'(a)(b - a) \\\\
&\Longrightarrow f'(a) \leq \frac{f(b) - f(a)}{b - a},
\end{align*}
\]
\[
\begin{align*}
f(a) &\geq f(b) + f'(b)(a - b) \\\\
&\Longrightarrow f'(b) \geq \frac{f(b) - f(a)}{b - a}.
\end{align*}
\]
Hence \(f'(a) \leq f'(b)\) whenever \(a \lt b\), so \(f'\) is monotonically non-decreasing. For a
differentiable \(f'\), this is equivalent to \(f'' \geq 0\).
Convex functions appear throughout machine learning. For example, the
cross-entropy loss for
classification is convex in the
predicted probability vector, and the ReLU activation \(\phi(x) = \max(0, x)\) used in
neural networks is convex on
\(\mathbb{R}\). (Note that cross-entropy composed with a deep network is convex in the probabilities but
generally non-convex in the network parameters. This failure of convexity is the origin of the non-convex
optimization landscape that dominates deep learning.)
For the multivariate case, we reduce to the univariate theorem via restriction to lines.
Theorem: Second-Order Condition for Convexity (Multivariate)
Let \(f: \mathcal{S} \to \mathbb{R}\) be \(C^2\) on an open convex set \(\mathcal{S} \subseteq \mathbb{R}^n\).
Then:
\(f\) is convex on \(\mathcal{S}\) \(\iff\) \(\nabla^2 f(\mathbf{x}) \succeq 0\) for all \(\mathbf{x} \in \mathcal{S}\).
If \(\nabla^2 f(\mathbf{x}) \succ 0\) for all \(\mathbf{x} \in \mathcal{S}\), then \(f\) is strictly
convex (the converse fails in general, since \(f(x) = x^4\) is strictly convex while \(f''(0) = 0\)).
Proof:
For \(\mathbf{x}, \mathbf{y} \in \mathcal{S}\), choose an open interval \(J \supseteq [0,1]\) with
\(\mathbf{x} + t(\mathbf{y}-\mathbf{x}) \in \mathcal{S}\) for all \(t \in J\) (possible because the
segment is compact and \(\mathcal{S}\) is open), and define the one-variable restriction
\(\varphi: J \to \mathbb{R}\) by
\[
\varphi(t) = f\!\left(\mathbf{x} + t(\mathbf{y} - \mathbf{x})\right).
\]
Since \(f \in C^2\) and the map \(t \mapsto \mathbf{x} + t(\mathbf{y}-\mathbf{x})\) is smooth, \(\varphi\) is
\(C^2\) on \(J\) by the chain rule.
A direct computation gives
\[
\varphi''(t) = (\mathbf{y} - \mathbf{x})^\top\, \nabla^2 f\!\left(\mathbf{x} + t(\mathbf{y}-\mathbf{x})\right)\, (\mathbf{y} - \mathbf{x}).
\]
Equivalence of convexity on restrictions. \(f\) is convex on \(\mathcal{S}\) if and only if
every such \(\varphi\) is convex on \(J\). The "only if" direction is immediate from the definition,
and the "if" direction follows by taking \(t = \lambda\) in
\(\varphi(\lambda) \leq \lambda \varphi(1) + (1-\lambda)\varphi(0)\).
(1) By the univariate theorem, each \(\varphi\) is convex \(\iff\) \(\varphi''(t) \geq 0\) on
\(J\). Varying \(\mathbf{x}, \mathbf{y}\), and using that \(\mathbf{y} - \mathbf{x}\) ranges over a
ball around \(\mathbf{0}\) because \(\mathcal{S}\) is open while the inequality is homogeneous in
\(\mathbf{v}\), this translates to
\(\mathbf{v}^\top \nabla^2 f(\mathbf{p}) \mathbf{v} \geq 0\) for all \(\mathbf{p} \in \mathcal{S}\) and all
\(\mathbf{v} \in \mathbb{R}^n\), that is, \(\nabla^2 f \succeq 0\) on \(\mathcal{S}\).
(2) If \(\nabla^2 f \succ 0\) on \(\mathcal{S}\), then \(\varphi''(t) \gt 0\) for
\(\mathbf{x} \neq \mathbf{y}\), so the remainder in the Taylor argument for the univariate theorem
is strictly positive for \(y \neq x\), the inequalities of the lemma become strict, and each
\(\varphi\) is strictly convex, which gives strict convexity of \(f\).
Why Hessian structure matters for ML
The Hessian governs more than convexity. Its spectrum controls, for a positive-definite Hessian,
the condition number \(\kappa = \lambda_{\max} / \lambda_{\min}\), which governs
the worst-case rate of gradient-based methods on strongly convex problems. We develop that theme
formally in the convergence analysis later in this
section. In deep learning, the Hessian's effective rank, spectral gap, and curvature distribution
at initialization have become active research objects. They are studied for their influence on
generalization and trainability, and they inform the choice of preconditioner.
Gradient Descent
By the first-order Taylor expansion,
\[
\mathcal{L}(\boldsymbol{\theta} + \mathbf{d})
\approx \mathcal{L}(\boldsymbol{\theta}) + \nabla \mathcal{L}(\boldsymbol{\theta})^\top \mathbf{d}
\quad (\|\mathbf{d}\| \text{ small}).
\]
To make \(\mathcal{L}\) decrease fastest, we want \(\mathbf{d}\) to minimize
\(\nabla \mathcal{L}(\boldsymbol{\theta})^\top \mathbf{d}\) subject to \(\|\mathbf{d}\| = 1\). Assuming
\(\nabla \mathcal{L}(\boldsymbol{\theta}) \neq \mathbf{0}\), the Cauchy-Schwarz inequality
\(|\mathbf{a}^\top\mathbf{b}| \leq \|\mathbf{a}\|\,\|\mathbf{b}\|\), with equality exactly when the two vectors
are linearly dependent, which we take as known, shows that
\(\nabla \mathcal{L}(\boldsymbol{\theta})^\top \mathbf{d} \geq -\|\nabla \mathcal{L}(\boldsymbol{\theta})\|\)
for every unit vector \(\mathbf{d}\), with equality only at
\[
\mathbf{d}^\star = -\frac{\nabla \mathcal{L}(\boldsymbol{\theta})}{\|\nabla \mathcal{L}(\boldsymbol{\theta})\|},
\]
the direction of steepest descent. (At a stationary point
\(\nabla \mathcal{L}(\boldsymbol{\theta}) = \mathbf{0}\), the first-order model predicts no change in any
direction, and the method halts.)
Gradient descent iterates along the steepest-descent direction with a chosen step size \(\eta \gt 0\):
\[
\boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} - \eta\, \nabla \mathcal{L}(\boldsymbol{\theta}^{(k)}).
\]
Under suitable regularity (for example, a Lipschitz-continuous gradient with constant \(L\), called
\(L\)-smoothness, together with \(\eta \leq 1/L\) and \(\mathcal{L}\) bounded below), the gradient
norms \(\|\nabla \mathcal{L}(\boldsymbol{\theta}^{(k)})\|\) tend to zero.
Algorithm 1: GRADIENT_DESCENTInput: objective \(\mathcal{L}\), tolerance \(\epsilon\), learning rate \(\eta\); Output: \(\boldsymbol{\theta}\) with \(\|\nabla \mathcal{L}(\boldsymbol{\theta})\| \lt \epsilon\); begin
\(k \leftarrow 0\); choose an initial point \(\boldsymbol{\theta}^{(0)}\); repeat
\(\mathbf{d}^{(k)} \leftarrow -\nabla \mathcal{L}(\boldsymbol{\theta}^{(k)})\);
\(\boldsymbol{\theta}^{(k+1)} \leftarrow \boldsymbol{\theta}^{(k)} + \eta\, \mathbf{d}^{(k)}\);
\(k \leftarrow k + 1\); until \(\|\nabla \mathcal{L}(\boldsymbol{\theta}^{(k)})\| \lt \epsilon\);
Output \(\boldsymbol{\theta}^{(k)}\); end
Gradient descent is called a first-order method because it uses only the gradient. Choosing
\(\eta\) optimally, rather than by trial-and-error, leads to the line search machinery. The
Armijo and
Wolfe conditions are developed
together with Newton's method. Analyzing how fast gradient
descent converges leads to the contraction-mapping analysis at the
end of this section.
Stochastic Gradient Descent
In stochastic optimization, the objective is an expectation:
\[
\mathcal{L}(\boldsymbol{\theta}) = \mathbb{E}_{q(z)}\bigl[\ell(\boldsymbol{\theta}, z)\bigr],
\]
where \(z\) is a random input (a training example, or a latent noise term) whose distribution \(q(z)\) does not
depend on \(\boldsymbol{\theta}\). At each iteration, stochastic gradient descent (SGD) samples
\(z^{(k)} \sim q\) and updates
\[
\boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} - \eta^{(k)}\, \nabla_{\boldsymbol{\theta}}\,
\ell(\boldsymbol{\theta}^{(k)}, z^{(k)}).
\]
By the Leibniz rule (interchange of expectation and gradient, valid under mild regularity since \(q(z)\) does
not depend on \(\boldsymbol{\theta}\)),
\(\mathbb{E}_{z^{(k)}}[\nabla_{\boldsymbol{\theta}}\, \ell(\boldsymbol{\theta}, z^{(k)})] = \nabla \mathcal{L}(\boldsymbol{\theta})\),
so SGD follows the full gradient in expectation. Classical convergence theory, which goes back to Robbins and
Monro and is not reproduced on this page, works with a diminishing step-size schedule
\(\sum_k \eta^{(k)} = \infty\) and \(\sum_k (\eta^{(k)})^2 \lt \infty\), under which, given suitable smoothness
and boundedness assumptions on \(\ell\) and with \(\mathcal{L}\) bounded below, the gradients
\(\nabla \mathcal{L}(\boldsymbol{\theta}^{(k)})\) tend to zero almost surely.
In practice, computing the full-batch gradient on every step is prohibitive when the dataset size
\(N\) is large. The mini-batch approach is the standard compromise. At each step, we
draw a random subset \(B \subset \{1, \ldots, N\}\) of size \(|B|\), typically \(32\), \(64\), or
\(128\), and use the batch-averaged gradient
\[
\hat{\mathbf{g}}^{(k)} = \frac{1}{|B|} \sum_{i \in B} \nabla_{\boldsymbol{\theta}}\, \ell(\boldsymbol{\theta}^{(k)}, z_i).
\]
This is an unbiased estimator of \(\nabla \mathcal{L}\) with lower variance than a single-sample
gradient. Mini-batching also maps cleanly onto GPU parallelism. Both properties make it the default in
modern deep learning.
Algorithm 2: MINI_BATCH_SGDInput: dataset \(X\), objective \(\ell\), tolerance \(\epsilon\), batch size \(|B|\), learning rate \(\eta\), max_epoch; Output: approximate stationary point \(\boldsymbol{\theta}^*\); begin
\(k \leftarrow 0\); \(e \leftarrow 0\); choose an initial point \(\boldsymbol{\theta}^{(0)}\); repeat
Shuffle \(X\) randomly;
Partition \(X\) into mini-batches \(\{B_1, B_2, \ldots, B_m\}\) of size \(|B|\) (the last possibly smaller); for each mini-batch \(B_j\):
\(\hat{\mathbf{g}}^{(k)} \leftarrow \frac{1}{|B_j|} \sum_{i \in B_j} \nabla_{\boldsymbol{\theta}}\, \ell(\boldsymbol{\theta}^{(k)}, z_i)\);
\(\boldsymbol{\theta}^{(k+1)} \leftarrow \boldsymbol{\theta}^{(k)} - \eta\, \hat{\mathbf{g}}^{(k)}\);
\(k \leftarrow k + 1\); end for
\(e \leftarrow e + 1\); until \(\|\hat{\mathbf{g}}^{(k-1)}\| \lt \epsilon\) or \(e \geq\) max_epoch;
Output \(\boldsymbol{\theta}^{(k)}\); end
One full pass over the shuffled dataset is called an epoch. Shuffling between epochs
is important. Without it, the same mini-batches recur in the same order every epoch, so the gradient
noise becomes a fixed periodic pattern rather than fresh randomness, and it is this randomness that is
often credited with helping SGD escape shallow basins.
A sample implementation of mini-batch SGD for linear regression follows.
import numpy as np
# Fixed parameters for mini-batch SGD
MAX_EPOCH = 10000
BATCH_SIZE = 64
TOLERANCE = 1e-3
LEARNING_RATE = 0.01
# Sample data dimensions
N_SAMPLES = 10000
N_FEATURES = 3
# Mini-batch Stochastic Gradient Descent
def mini_batch_sgd(X, y):
n_samples = X.shape[0]
theta = np.random.randn(X.shape[1]) * 0.01 # initial theta
for i in range(MAX_EPOCH):
# Shuffle data
indices = np.random.permutation(n_samples)
x_shuffled = X[indices]
y_shuffled = y[indices]
for j in range(0, n_samples, BATCH_SIZE):
end = j + BATCH_SIZE
x_batch = x_shuffled[j:end]
y_batch = y_shuffled[j:end]
# Gradient of MSE on the mini-batch
g = grad_f(theta, x_batch, y_batch)
# Update parameters
theta -= LEARNING_RATE * g
# Convergence check on the last mini-batch gradient
if np.linalg.norm(g) < TOLERANCE:
print(f"Converged in {i + 1} epochs.")
break
return theta
# Gradient of the MSE objective for linear regression
def grad_f(theta, x_batch, y_batch):
return (1 / x_batch.shape[0]) * x_batch.T @ ((x_batch @ theta) - y_batch)
# Main
if __name__ == "__main__":
# Generate sample data
X = np.random.rand(N_SAMPLES, N_FEATURES)
true_theta = np.array([5, -1, -9])
y = X @ true_theta + np.random.randn(X.shape[0]) * 0.1 # additive noise
# Run mini-batch SGD
estimated_theta = mini_batch_sgd(X, y)
# Relative error
relative_error = np.linalg.norm(estimated_theta - true_theta) / np.linalg.norm(true_theta)
print("True Parameters: ", true_theta)
print("Estimated Parameters:", estimated_theta)
print(f"Relative Error: {relative_error:.6f}")
Sub-gradient Descent
Many ML objectives are not differentiable everywhere, notably those using \(\ell_1\) regularization,
ReLU, or hinge loss. Gradient descent is not directly applicable at the kinks. When the objective is convex, as
for a convex loss of a linear model plus an \(\ell_1\) penalty, the remedy is to generalize the gradient to a
set-valued notion.
Definition: Subgradient
Let \(f\) be a convex function on a convex set \(\operatorname{dom}(f) \subseteq \mathbb{R}^n\) and
let \(\mathbf{x} \in \operatorname{dom}(f)\). A vector
\(\mathbf{g} \in \mathbb{R}^n\) is a subgradient of \(f\) at \(\mathbf{x}\) if
\[
f(\mathbf{z}) \geq f(\mathbf{x}) + \mathbf{g}^\top (\mathbf{z} - \mathbf{x})
\quad \text{for all } \mathbf{z} \in \operatorname{dom}(f). \tag{\dagger}
\]
Definition: Subdifferential
The subdifferential of \(f\) at \(\mathbf{x}\), denoted \(\partial f(\mathbf{x})\),
is the set of all subgradients of \(f\) at \(\mathbf{x}\):
\[
\partial f(\mathbf{x}) = \{\, \mathbf{g} \in \mathbb{R}^n : (\dagger) \text{ holds for all } \mathbf{z} \,\}.
\]
We say \(f\) is subdifferentiable at \(\mathbf{x}\) if
\(\partial f(\mathbf{x}) \neq \emptyset\). For a convex function, \(\partial f(\mathbf{x})\) is
closed and convex, being the intersection over \(\mathbf{z}\) of the closed half-spaces in
\(\mathbf{g}\) cut out by (\(\dagger\)). At every interior point of \(\operatorname{dom}(f)\) it
is also non-empty, and there \(\partial f(\mathbf{x}) = \{\nabla f(\mathbf{x})\}\) precisely when
\(f\) is differentiable at \(\mathbf{x}\). These two facts rest on a separating-hyperplane
argument that we do not reproduce
on this page.
Example: absolute value.
For \(f(x) = |x|\) on \(\mathbb{R}\),
\[
\partial f(x) = \begin{cases}
\{-1\} & \text{if } x \lt 0,\\\\
[-1, 1] & \text{if } x = 0,\\\\
\{+1\} & \text{if } x \gt 0.
\end{cases}
\]
At \(x \neq 0\), \(f\) is differentiable and the subdifferential collapses to the ordinary derivative. At
\(x = 0\), every slope \(g \in [-1, 1]\) satisfies \(|z| \geq 0 + g \cdot z\) for all \(z\), giving a full
interval of subgradients.
Sub-gradient descent selects any
\(\mathbf{g}^{(k)} \in \partial \mathcal{L}(\boldsymbol{\theta}^{(k)})\) and updates
\[
\boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} - \eta^{(k)}\, \mathbf{g}^{(k)}.
\]
Because \(-\mathbf{g}^{(k)}\) need not be a descent direction, the objective values need not decrease
monotonically, and because the chosen subgradient need not vanish at the optimum (the subdifferential there
contains \(\mathbf{0}\) but is in general larger), a diminishing step-size schedule is essential for
convergence. For a convex objective with a minimizer and bounded subgradients, the best objective value found in
the first \(k\) steps approaches the optimum at rate \(O(1/\sqrt{k})\) under a suitable step-size choice, slower
than gradient descent's \(O(1/k)\) under \(L\)-smoothness. The gap is the price of handling non-smoothness.
Why this matters for modern ML
The subdifferential underpins the practical use of non-smooth regularizers. Lasso's \(\ell_1\)
penalty, the SVM hinge loss, and the ReLU activation are all analyzed through their
subdifferentials. The proximal gradient method and ISTA/FISTA
algorithms explicitly exploit this structure to achieve sparse solutions. Automatic
differentiation frameworks such as PyTorch and JAX resolve the ambiguity at kinks by a fixed
convention, which for ReLU is a canonical element of the subdifferential (\(0\) for
\(\mathrm{ReLU}'(0)\) in both PyTorch and JAX). The convention is mathematically arbitrary, since
any element of the subdifferential is admissible.