Gradient Descent and Numerical Optimization › Implementation in NumPy
Machine Learning from ScratchAll chapters Gradient Descent and Numerical Optimization › Implementation in NumPySection view · Sections in this chapter7 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
We now implement least-squares fitting with four synthetic features and an explicit intercept. The design matrix contains a leading column of ones, so the parameter vector has five entries: the intercept followed by four slopes. Including that intercept is essential to this comparison: centering the features can then be absorbed into the intercept without changing the family of predictions. We compare gradient descent with NumPy's least-squares implementation, using the same design matrix and objective.
PythonIn [19]Expand
import numpy as np
rng = np.random.default_rng(7)
def make_regression(n=200, noise=0.5, seed=0):
"""Four synthetic features and an explicit, unscaled intercept column."""
rng = np.random.default_rng(seed)
features = rng.normal(size=(n, 4))
slopes = np.array([2.0, -1.0, 0.5, 0.0])
y = features @ slopes + noise * rng.normal(size=n)
X = np.column_stack([np.ones(n), features])
beta = np.concatenate([[0.0], slopes])
return X, y, beta
X_raw, y, beta_true = make_regression()
def sse_objective(theta, X, y):
"""Half the mean squared error for this linear design."""
r = y - X @ theta
return 0.5 * np.mean(r**2)
def sse_gradient(theta, X, y):
"""Analytic gradient of half the mean squared error."""
r = y - X @ theta
return -X.T @ r / len(y)
def full_batch_gd(X, y, eta=0.1, steps=200, tol=1e-10):
"""Full gradients; count gradient checks and reject nonfinite iterates."""
theta = np.zeros(X.shape[1])
hist = [sse_objective(theta, X, y)]
gradient_checks = 0
for _ in range(steps):
g = sse_gradient(theta, X, y)
gradient_checks += 1
if not np.all(np.isfinite(g)):
raise FloatingPointError(
"Nonfinite gradient; inspect the learning rate."
)
if np.linalg.norm(g) < tol:
break
candidate = theta - eta * g
with np.errstate(over="ignore", invalid="ignore"):
value = sse_objective(candidate, X, y)
if not np.all(np.isfinite(candidate)) or not np.isfinite(value):
raise FloatingPointError(
"Nonfinite iterate; reduce the learning rate."
)
theta = candidate
hist.append(value)
return theta, np.array(hist), gradient_checksConsoleOut [20]Expand
design matrix: 200 rows, 5 columns feature stds (intercept excluded): [0.973 0.946 1.027 1.042] true coefficients: [ 0. 2. -1. 0.5 0. ] initial objective: 2.5285
The objective is half the mean squared error, , where is the row count, contains the targets, and contains the predictions. Differentiating gives ; the factor of two from the square cancels the one-half in the objective. Duplicating every row leaves this average and its gradient unchanged, although changing the data distribution can still change the curvature and a suitable learning rate. The corresponding sum of squared errors is .
Before interpreting a fit, compare the analytic gradient with a finite difference. This is the numerical check introduced in the formula section; it can expose sign errors or missing factors without relying on the eventual training loss.
PythonIn [21]Expand
def finite_difference_gradient(f, theta, h=1e-6):
"""Central difference estimate, one coordinate at a time."""
theta = np.asarray(theta, dtype=float)
grad = np.zeros_like(theta)
for j in range(len(theta)):
e = np.zeros_like(theta)
e[j] = h
grad[j] = (f(theta + e) - f(theta - e)) / (2 * h)
return grad
theta_test = np.random.default_rng(11).normal(size=X_raw.shape[1])
g_analytic = sse_gradient(theta_test, X_raw, y)
g_numeric = finite_difference_gradient(
lambda t: sse_objective(t, X_raw, y), theta_test
)
max_rel_err = np.max(
np.abs(g_analytic - g_numeric) / np.maximum(1e-8, np.abs(g_numeric))
)
assert max_rel_err < 1e-5ConsoleOut [22]Expand
analytic gradient: [ 0.192819 -0.607492 1.941206 -1.096276 -0.270025] numeric gradient: [ 0.192819 -0.607492 1.941206 -1.096276 -0.270025] Gradient check passed with max relative error 1.86e-09
With the gradient check passed, we fit the model and compare with np.linalg.lstsq, an established least-squares solver. This comparison checks the fitted coefficients and objective for this full-rank design.
PythonIn [23]Expand
# Closed-form reference
theta_ols = np.linalg.lstsq(X_raw, y, rcond=None)[0]
theta_gd, hist_gd, used = full_batch_gd(X_raw, y, eta=0.1, steps=2000)
coef_diff = np.max(np.abs(theta_ols - theta_gd))ConsoleOut [24]Expand
parameter updates: 252; full-gradient checks: 253 stopping reason: gradient tolerance returned gradient norm: 9.357e-11 OLS: [-0.06483 1.96236 -0.99078 0.51876 -0.00477] GD: [-0.06483 1.96236 -0.99078 0.51876 -0.00477] max abs difference: 7.83e-11 objective OLS 0.09391853 vs GD 0.09391853
The four feature columns are approximately standard normal, and the leading intercept column remains one. We multiply the first feature by 200 to change its units while leaving the available predictions unchanged: its slope can be divided by 200. A fixed learning rate of 0.1 becomes unstable on this design. We show only six such updates to keep the demonstration finite, then use for a longer stable run on the same badly scaled problem.
PythonIn [25]Expand
X_bad = X_raw.copy()
X_bad[:, 1] = X_bad[:, 1] * 200.0
# Condition number of X^T X is the square of the design matrix condition number.
cond_raw = np.linalg.cond(X_raw)
cond_bad = np.linalg.cond(X_bad)
hessian_bad = X_bad.T @ X_bad / len(X_bad)
lambda_max_bad = np.linalg.eigvalsh(hessian_bad)[-1]
unsafe_eta = 0.1
theta_unsafe, hist_unsafe, checks_unsafe = full_batch_gd(
X_bad, y, eta=unsafe_eta, steps=6
)
eta_bad = 1.0 / lambda_max_bad
theta_bad, hist_bad, used_bad = full_batch_gd(X_bad, y, eta=eta_bad, steps=2000)
assert unsafe_eta > 2.0 / lambda_max_bad
assert hist_unsafe[-1] > hist_unsafe[0]
assert np.all(np.isfinite(hist_bad))ConsoleOut [26]Expand
cond(X_raw) = 1.2, cond(X^T X) = 1.4 cond(X_bad) = 207.1, cond(X^T X) = 42870.3 Well-scaled objective: 0.093919; gradient checks: 253 Unsafe eta=0.1: objective after 6 updates = 1.612e+43 Bad-design stability ceiling: 5.281e-05 unsafe stopping reason: deliberate six-update budget unsafe returned gradient norm: 1.105e+24 Stable eta=2.641e-05: objective = 0.616614; gradient checks: 2000 safe badly scaled stopping reason: update budget safe badly scaled returned gradient norm: 9.840e-01
The six large-step updates increase the objective sharply because 0.1 exceeds this design's stability ceiling. The second badly scaled run remains finite at , but makes slow progress in the shallow directions. Standardizing the four feature columns lets us use a much larger stable step. We leave the intercept column unchanged because its standard deviation is zero.
PythonIn [27]Expand
mu = X_bad[:, 1:].mean(axis=0)
sd = X_bad[:, 1:].std(axis=0)
assert np.all(sd > 0)
X_scaled = np.column_stack([np.ones(len(X_bad)), (X_bad[:, 1:] - mu) / sd])
theta_scaled, hist_scaled, used_scaled = full_batch_gd(
X_scaled, y, eta=0.1, steps=2000
)
cond_scaled = np.linalg.cond(X_scaled)
theta_scaled_in_bad_units = np.concatenate(
[
[theta_scaled[0] - mu @ (theta_scaled[1:] / sd)],
theta_scaled[1:] / sd,
]
)
scaled_reference = np.linalg.lstsq(X_scaled, y, rcond=None)[0]
assert np.allclose(
X_bad @ theta_scaled_in_bad_units, X_scaled @ theta_scaled, atol=1e-10
)
assert np.allclose(X_raw @ theta_ols, X_scaled @ scaled_reference, atol=1e-10)
assert np.isclose(
hist_scaled[-1], sse_objective(theta_ols, X_raw, y), atol=1e-10
)ConsoleOut [28]Expand
Standardized: cond(X)=1.1, J=0.093919, gradient checks=234 standardized stopping reason: gradient tolerance standardized returned gradient norm: 9.087e-11 Badly scaled, stable step: cond(X)=207.1, J=0.616614, gradient checks=2000 Same OLS optimum: 0.093919 Prediction-map max error: 8.88e-16
With an explicit intercept, standardization preserves the set of possible linear predictions. To map standardized slopes back, divide each by its feature standard deviation; subtract the mean-weighted slopes from the standardized intercept. The assertions above verify that this mapping preserves predictions and that the two reference fits attain the same least-squares optimum. The finite-budget badly scaled gradient-descent run can still be farther from that optimum. All 200 rows here are a synthetic training set; in an evaluated model, estimate scaling statistics on its training split and reuse them unchanged for validation and test data.
Learning rate from the Hessian spectrum
For this full-rank quadratic, the Hessian spectrum gives the stable interval . With our half-mean-squared-error objective the Hessian is , where is the design and is its largest eigenvalue. The next cell computes the eigenvalues and prints , where is the smallest eigenvalue. This rate minimizes the worst-case Euclidean parameter-error contraction for a fixed step on this quadratic.
The update budget trades computation against closeness to the optimizer. Stopping early can also change a model's generalization behavior; its relation to regularization depends on the model and update rule. Tolerances need to match the units and scale of the parameterization and objective. A small absolute gradient threshold is an algorithmic condition, not a universal accuracy measure.
PythonIn [29]Expand
def hessian_spectrum(X):
"""Eigenvalues of the least-squares Hessian X^T X / n."""
H = X.T @ X / len(X)
return np.sort(np.linalg.eigvalsh(H))
ev_scaled = hessian_spectrum(X_scaled)
ev_bad = hessian_spectrum(X_bad)
eta_opt_scaled = 2.0 / (ev_scaled[0] + ev_scaled[-1])
eta_opt_bad = 2.0 / (ev_bad[0] + ev_bad[-1])ConsoleOut [30]Expand
scaled eigenvalues: [0.9371 0.9771 1. 1.0073 1.0785] kappa 1.15, eta_max 1.8545, eta_opt 0.9923 unscaled eigenvalues: [8.800000e-01 8.900000e-01 1.090000e+00 1.190000e+00 3.787083e+04] kappa 42870.3, eta_max 0.000053, eta_opt 0.000053
The badly scaled problem has both a small stable-rate ceiling and a large eigenvalue ratio. At the optimal fixed rate, the worst-case parameter-error contraction is , where . Its factor is therefore much closer to one than the standardized problem's factor. This compares worst-case relative error across starting directions; a particular run also depends on which eigendirections are present in its initial error.
The eigenvalues make the conditioning gap visible. The standardized design has a spectrum clustered near one, while the badly scaled design spans several orders of magnitude. That spread worsens the best worst-case fixed-rate contraction on this quadratic.
VisualizationOut [31]Expand
Stochastic and mini-batch methods
The next implementation reshuffles the training rows each epoch and updates using the mean gradient of the current batch. We use a nonempty design, nonnegative integer epoch budgets, positive integer batch sizes, a positive initial rate, and nonnegative decay. We count per-example gradient contributions used for updates, and separately note that objective evaluation and diagnostic checks also cost work. A final short mini-batch is divided by its actual number of rows.
PythonIn [32]Expand
def sgd(X, y, eta0=0.05, epochs=30, decay=None, seed=0, batch_size=None):
"""SGD with optional mini-batch and learning-rate decay.
Returns parameters, one objective per epoch, and example contributions.
"""
rng = np.random.default_rng(seed)
n, p = X.shape
theta = np.zeros(p)
b = 1 if batch_size is None else batch_size
hist, evals = [], 0
for epoch in range(epochs):
order = rng.permutation(n)
for start in range(0, n, b):
idx = order[start : start + b]
eta = eta0 if decay is None else eta0 / (1 + decay * evals / n)
g = -X[idx].T @ (y[idx] - X[idx] @ theta) / len(idx)
theta = theta - eta * g
evals += len(idx)
hist.append(sse_objective(theta, X, y))
return theta, np.array(hist), evals
theta_sgd, hist_sgd, evals_sgd = sgd(
X_scaled, y, eta0=0.1, epochs=40, decay=1.0, batch_size=1, seed=3
)
theta_mb, hist_mb, evals_mb = sgd(
X_scaled, y, eta0=0.1, epochs=40, decay=1.0, batch_size=32, seed=3
)
full_batch_evals = len(X_scaled) * (len(hist_scaled) - 1)ConsoleOut [33]Expand
SGD : J = 0.093933, update contributions = 8000 mini-batch: J = 0.104950, update contributions = 8000 full batch: J = 0.093919, update contributions = 46600 SGD and mini-batch stopping reason: completed the 40-epoch budget SGD returned full-gradient norm: 5.436e-03 mini-batch returned full-gradient norm: 1.485e-01
These are results for the stated data, learning rates, sampling seed, and finite budgets. They do not prove that one variant always reaches a target objective with fewer operations. The decay rule changes the step according to the number of examples already used for updates; it is a choice for this illustration, not a general convergence certificate.
The horizontal axis below counts per-example gradient contributions in parameter updates. For full-batch descent it increases by per update; the stochastic curves are recorded after each epoch of contributions. This omits the extra cost of evaluating full objectives, terminal gradient checks, and implementation overhead, so it is not a wall-clock comparison.
PythonIn [34]Expand
eval_sgd = np.arange(1, len(hist_sgd) + 1) * len(X_scaled)
eval_mb = np.arange(1, len(hist_mb) + 1) * len(X_scaled)
eval_full = np.arange(len(hist_scaled)) * len(X_scaled)VisualizationOut [35]Expand
Optional: Newton on least squares
On the full-rank least-squares design, the Hessian is constant and positive definite. We can therefore demonstrate the full Newton step directly. The implementation below is specific to this quadratic and does not add a line search or a logistic-loss example. General nonquadratic objectives need the safeguards discussed above.
PythonIn [36]Expand
def newton_ls(X, y, steps=5):
"""Newton on full-rank least squares, with one history entry per update."""
H = X.T @ X / len(X)
theta = np.zeros(X.shape[1])
hist = [sse_objective(theta, X, y)]
for _ in range(steps):
g = sse_gradient(theta, X, y)
if np.linalg.norm(g) < 1e-10:
break
theta = theta - np.linalg.solve(H, g)
hist.append(sse_objective(theta, X, y))
return theta, np.array(hist)
theta_newton, hist_newton = newton_ls(X_bad, y, steps=5)ConsoleOut [37]Expand
Newton on ill-conditioned LS: J = 0.0939185272, steps used 1 Newton objective history: [2.52846336 0.09391853] Newton solution: [-0.065 0.01 -0.991 0.519 -0.005] Gradient norm at the solution: 3.03e-14 Newton stopping reason: gradient tolerance coefficients match OLS: True
The printed gradient norm and reference comparison quantify the numerical error. One update reaches the least-squares solution to the reported precision on this design. The exact-arithmetic result is independent of the condition number, while the accuracy of a floating-point linear solve can still depend on scaling and conditioning.
Key Parameters
Below are the main parameters that affect how gradient-based optimization behaves, with guidance on how to select them.
etascales the negative gradient. For the positive-definite quadratics here, convergence from every start requires , and minimizes the worst-case fixed-rate parameter contraction. General objectives need an appropriate step-selection rule.stepscaps the update count. More iterations can reduce optimization error in a stable run; they do not repair an unstable rate. Early stopping can regularize some fitting procedures, but the optimization examples here do not evaluate that effect on held-out data.tolis the absolute gradient-norm threshold used by this illustrative implementation. Interpret it in the current objective and parameter units, and report the returned residual and budget.batch_sizeis the number of examples used per update. It changes gradient noise, work per update, and updates per epoch; 32 is the setting demonstrated here.decaycontrols the displayed schedule , where counts per-example gradient contributions already used. The code uses this example-count clock, which differs from counting parameter updates.
Before increasing the budget or loosening a tolerance, check the gradient calculation and the response to the learning rate. A decay schedule and a batch size should be assessed together under a stated budget and sampling rule.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!