import plotly.io as pio; pio.renderers.default = "notebook_connected"
Lecture 11: Convergence, Momentum, Adam, and Stochastic Gradient Descent (CS 189, Fall 2026)
This notebook accompanies Lecture 11 and Chapter 7 of Bishop and Bishop, Deep Learning: Foundations and Concepts. It continues from the Lecture 10 notebook and covers:
- the convergence of gradient descent on a quadratic error surface,
- momentum,
- adaptive learning rates (AdaGrad, RMSProp, and Adam),
- learning rate schedules and warmup, and
- stochastic and mini-batch gradient descent.
Sections used in lecture are marked in their headings. (Lecture Demo n) is run during the lecture from the matching "Demo" slide.
| Notebook section | Lecture 11 slide |
|---|---|
| Recap: The Quadratic Approximation | Recap slides; the eigenvector plots ("Eigenvectors of the Hessian") are Lecture Demo 5 in the Lecture 10 notebook |
| (Lecture Demo 1) Convergence of Gradient Descent on a Quadratic | Demo: Oscillating convergence; Issues with Gradient Descent (plots) |
| (Lecture Demo 2) Gradient Descent with Momentum | Demo: Momentum |
| (Lecture Demo 3) Adaptive Learning Rates | Demo: Adam |
| (Lecture Demo 4) Learning Rate Schedules | Demo: Learning Rate Schedules |
| (Lecture Demo 5) Stochastic Gradient Descent | Demo: Stochastic Gradient Descent |
| (Lecture Demo 6) Mini-Batch Stochastic Gradient Descent | Demo: Mini-batch Gradient Descent; How do I Pick the Batch Size? |
import numpy as np
import pandas as pd
import plotly.express as px
import plotly.graph_objects as go
from plotly import figure_factory as ff
from plotly.subplots import make_subplots
colors = px.colors.qualitative.Plotly
px.defaults.width = 800
Plotting Code
The functions below produce the visualizations used throughout the notebook. They are not part of the course material, but they must be run before the rest of the notebook.
def make_plot_grid(figs, rows, cols):
"""Arrange several Plotly figures in a grid of subplots.
The traces and the 2D axis titles of each figure are copied into the grid.
Layout elements such as shapes (e.g., add_vline) and annotations are not
copied and must be added to the returned figure directly.
"""
def get_trace_type(fig):
for trace in fig.data:
if trace.type == 'surface':
return 'surface'
if trace.type in ('scatter3d', 'mesh3d', 'cone'):
return 'scene'
return 'xy'
types = [get_trace_type(fig) for fig in figs]
specs = [[{'type': ty} for ty in types[i:i + cols]]
for i in range(0, len(figs), cols)]
fig_grid = make_subplots(rows=rows, cols=cols, specs=specs,
subplot_titles=[fig.layout.title.text or '' for fig in figs])
for i, (fig, ty) in enumerate(zip(figs, types)):
row, col = i // cols + 1, i % cols + 1
fig_grid.add_traces(fig.data, rows=row, cols=col)
if ty == 'xy':
fig_grid.update_xaxes(title_text=fig.layout.xaxis.title.text, row=row, col=col)
fig_grid.update_yaxes(title_text=fig.layout.yaxis.title.text, row=row, col=col)
return fig_grid
# Fixed colors for the true class and fixed symbols for the predicted class.
DIAGNOSIS_COLORS = {'benign': 'blue', 'malignant': 'red'}
PREDICTION_SYMBOLS = {'Predicted 1 (benign)': 'cross', 'Predicted 0 (malignant)': 'circle-open'}
def plot_lr_predictions(w, cancer_df):
"""Plot the data, the decision boundary, and the predicted probability surface
of the logistic regression model p(t=1|x) = sigmoid(w^T x).
In this dataset t=1 denotes a benign tumor and t=0 a malignant tumor. Color shows the
true class (benign blue, malignant red). The symbol shows the predicted class:
+ for a prediction of 1 (benign) and o for a prediction of 0 (malignant).
"""
df = cancer_df.copy()
p = logistic_model(w, df[['mean radius', 'mean texture']].values)
df['prediction'] = np.where(p > 0.5, 'Predicted 1 (benign)', 'Predicted 0 (malignant)')
fig = px.scatter(df, x='mean radius', y='mean texture',
symbol='prediction', color='diagnosis',
symbol_map=PREDICTION_SYMBOLS, color_discrete_map=DIAGNOSIS_COLORS,
category_orders={'diagnosis': list(DIAGNOSIS_COLORS),
'prediction': list(PREDICTION_SYMBOLS)})
for (i, t) in enumerate(fig.data):
t.legendgroup = str(i)
# Decision boundary: w_1 x_1 + w_2 x_2 = 0
xs = np.linspace(df['mean radius'].min(), df['mean radius'].max(), 100)
decision_boundary = -(w[0] * xs) / w[1]
fig.add_scatter(x=xs, y=decision_boundary, mode='lines',
name='Decision Boundary', legendgroup='Decision Boundary',
line=dict(color='black', width=2, dash='dash'))
# Probability surface p(t=1|x)
ys = np.linspace(df['mean texture'].min(), df['mean texture'].max(), 100)
X, Y = np.meshgrid(xs, ys)
Z = logistic_model(w, np.c_[X.ravel(), Y.ravel()]).reshape(X.shape)
# Light red where the model predicts malignant (p near 0), white at the decision
# boundary (p = 0.5), and light blue where it predicts benign (p near 1).
fig.add_contour(x=xs, y=ys, z=Z, zmin=0, zmax=1,
colorscale=[[0, '#fddbc7'], [0.5, '#f7f7f7'], [1, '#d1e5f0']],
contours=dict(start=0.1, end=0.9, size=0.1), line=dict(width=0.5, color='gray'),
name='p(t=1|x)',
colorbar=dict(title='p(t=1|x)', x=1.05, y=0.3, len=0.75))
fig.update_layout(title=f'w=({w[0]:.2f}, {w[1]:.2f})',
xaxis_range=[xs.min(), xs.max()], yaxis_range=[ys.min(), ys.max()],
xaxis_title='Mean Radius (standardized)',
yaxis_title='Mean Texture (standardized)',
width=800, height=600)
return fig
def plot_loss(w1, w2, error, ncontours=50, names=('w1', 'w2')):
"""Plot the error surface (left) and its contours (right) on a grid of weights.
Args:
w1, w2: Meshgrid arrays of the two weights.
error: Array of error values with the same shape as w1 and w2.
ncontours: Number of contour levels.
names: Axis labels for the two weights.
"""
surf_fig = go.Figure()
surf_fig.add_surface(z=error, x=w1, y=w2,
colorscale='Viridis_r', opacity=0.7, showscale=False,
contours=dict(z=dict(show=True, highlightcolor="white",
start=error.min(), end=error.max(),
size=(error.max() - error.min()) / ncontours)))
surf_fig.update_layout(title="Error Surface")
contour_fig = go.Figure()
contour_fig.add_contour(x=w1.flatten(), y=w2.flatten(), z=error.flatten(),
colorscale='Viridis_r', opacity=0.7,
contours=dict(start=error.min(), end=error.max(),
size=(error.max() - error.min()) / ncontours),
colorbar=dict(x=1.05, y=0.35, len=0.75))
contour_fig.update_layout(title="Error Contours")
fig = make_plot_grid([surf_fig, contour_fig], 1, 2).update_layout(height=800)
fig.update_layout(scene=dict(xaxis_title=names[0], yaxis_title=names[1],
zaxis_title='Error', aspectmode='cube'))
fig.update_layout(xaxis_range=[w1.min(), w1.max()], yaxis_range=[w2.min(), w2.max()],
xaxis_title=names[0], yaxis_title=names[1])
return fig
def plot_gradient(w1, w2, error, dw1, dw2, scale=1.0, names=('w1', 'w2')):
"""Plot the error surface and contours together with the gradient field.
The cones (left) and arrows (right) show the gradient vector at each grid
point. The gradient points in the direction of steepest ascent, so gradient
descent moves in the opposite direction.
"""
fig = plot_loss(w1, w2, error, names=names)
fig.add_trace(
go.Cone(
x=w1.flatten(), y=w2.flatten(), z=np.zeros_like(error).flatten(), # ground plane
u=dw1.flatten(), v=dw2.flatten(), w=np.zeros_like(error).flatten(), # no vertical component
sizeref=2, anchor="tail", showscale=False
), 1, 1)
quiver = ff.create_quiver(
x=w1.flatten(), y=w2.flatten(), u=dw1.flatten(), v=dw2.flatten(),
line_width=2, line_color="white",
scale=scale, arrow_scale=.2, showlegend=False)
fig.add_traces(quiver.data, rows=1, cols=2)
return fig
def add_solution_path(fig, errors, ws):
"""Overlay an optimization path on a figure created by plot_loss or plot_gradient."""
s = np.linspace(0, 1, len(ws))
fig.add_scatter3d(x=ws[:, 0], y=ws[:, 1], z=errors, marker_color=s, marker_size=5,
mode='lines+markers', line=dict(color='black', width=2), opacity=0.5,
name='Optimization Path', legendgroup='Optimization Path',
row=1, col=1)
fig.add_scatter(x=ws[:, 0], y=ws[:, 1], marker_color=s,
mode='lines+markers', line=dict(color='black', width=2), opacity=0.5,
name='Optimization Path', legendgroup='Optimization Path',
showlegend=False,
row=1, col=2)
fig.add_scatter3d(x=[ws[-1, 0]], y=[ws[-1, 1]], z=[errors[-1]],
mode='markers', marker=dict(color='red', size=10),
name='Final Solution', legendgroup='Final Solution',
row=1, col=1)
fig.add_scatter(x=[ws[-1, 0]], y=[ws[-1, 1]],
mode='markers', marker=dict(color='red', size=20),
name='Final Solution', legendgroup='Final Solution',
showlegend=False,
row=1, col=2)
return fig
def make_loss_curve(path, error_func):
"""Plot the error as a function of the iteration number for a path of weights."""
errors = [error_func(w) for w in path]
fig = px.line(x=np.arange(len(errors)), y=errors,
labels={'x': 'Iteration (gradient steps)', 'y': 'Error'})
fig.update_traces(line_width=4)
fig.update_layout(margin=dict(l=20, r=20, t=20, b=20))
return fig
Setup: Models, Error Functions, and Gradients
This section recreates the two problems from the Lecture 10 notebook: logistic regression on two standardized features of the breast cancer data, and non-linear regression with the model $y(x, w) = \sin(w_0 + w_1 x)$. The derivations of the error functions and gradients are given in that notebook. Each gradient function accepts optional data arguments, so the same function can be used for batch, stochastic, and mini-batch gradient descent.
Logistic Regression
from sklearn import datasets
from sklearn.preprocessing import StandardScaler
cancer_dict = datasets.load_breast_cancer(as_frame=True)
cancer_df = pd.DataFrame(cancer_dict.data, columns=cancer_dict.feature_names)
cancer_df = cancer_df[['mean radius', 'mean texture']].copy()
# Target t=1 denotes a benign tumor and t=0 a malignant tumor.
cancer_df['target'] = cancer_dict.target.values
cancer_df['diagnosis'] = cancer_dict.target_names[cancer_df['target']]
# Standardize each feature to zero mean and unit variance.
scaler = StandardScaler()
cancer_df[['mean radius', 'mean texture']] = scaler.fit_transform(
cancer_df[['mean radius', 'mean texture']])
print("The dataset:", cancer_df.shape)
The dataset: (569, 4)
cancer = dict()
cancer['x'] = cancer_df[['mean radius', 'mean texture']].values
cancer['t'] = cancer_df['target'].values.astype(float)
cancer['npts'] = 30
cancer['w1'], cancer['w2'] = np.meshgrid(
np.linspace(-10, 1, cancer['npts']),
np.linspace(-5, 1.3, cancer['npts']))
cancer['ws'] = np.stack([cancer['w1'].flatten(), cancer['w2'].flatten()]).T
def sigmoid(z):
"""Logistic sigmoid function."""
return 1 / (1 + np.exp(-z))
def logistic_model(w, x):
"""Logistic regression model p(t=1|x) = sigmoid(w^T x) (no bias term)."""
return sigmoid(x @ w)
def neg_log_likelihood(w, x=None, t=None):
"""Average negative log-likelihood (cross-entropy) of the logistic model.
Uses the identity -[t ln sigmoid(z) + (1-t) ln(1-sigmoid(z))] = ln(1+e^z) - t z,
evaluated with the numerically stable np.logaddexp(0, z) = ln(1+e^z).
Defaults to the full training set.
"""
x = cancer['x'] if x is None else x
t = cancer['t'] if t is None else t
z = x @ w
return np.mean(np.logaddexp(0, z) - t * z)
def grad_NLL(w, x=None, t=None):
"""Gradient of the average negative log-likelihood of the logistic model.
Accepts the full data set (default), a mini-batch, or a single data point.
"""
x = cancer['x'] if x is None else np.atleast_2d(x)
t = cancer['t'] if t is None else np.atleast_1d(t)
y = logistic_model(w, x)
return np.mean((y - t).reshape(-1, 1) * x, 0)
cancer['error'] = np.array([neg_log_likelihood(w) for w in cancer['ws']])
cancer['error'] = cancer['error'].reshape(cancer['w1'].shape)
best_ind = np.argmin(cancer['error'])
cancer['grid_best'] = cancer['ws'][best_ind]
(cancer['dw1'], cancer['dw2']) = np.array([grad_NLL(w) for w in cancer['ws']]).T
cancer['dw1'] = cancer['dw1'].reshape(cancer['w1'].shape)
cancer['dw2'] = cancer['dw2'].reshape(cancer['w1'].shape)
Sine Regression
sine = dict()
np.random.seed(42)
sine['n'] = 200
sine['x'] = np.sort(np.random.rand(sine['n']) * 2.5 * np.pi) # sorted for plotting
sine['y'] = np.sin(1.1 + 2.5 * sine['x']) + 0.5 * np.random.randn(sine['n'])
sine_df = pd.DataFrame({'x': sine['x'], 'y': sine['y']})
def sine_model(w, x):
"""Non-linear regression model y(x, w) = sin(w_0 + w_1 x)."""
return np.sin(w[0] + x * w[1])
def sine_MSE(w, x=None, y=None):
"""Mean squared error of the sine model. Defaults to the full training set."""
x = sine['x'] if x is None else x
y = sine['y'] if y is None else y
return np.mean((y - sine_model(w, x)) ** 2)
def grad_sine_MSE(w, x=None, y=None):
"""Gradient of the mean squared error of the sine model.
Accepts the full data set (default), a mini-batch, or a single data point.
"""
x = sine['x'] if x is None else np.atleast_1d(x)
y = sine['y'] if y is None else np.atleast_1d(y)
u = w[0] + w[1] * x
r = y - np.sin(u)
grad_w0 = -2 * np.mean(r * np.cos(u))
grad_w1 = -2 * np.mean(r * np.cos(u) * x)
return np.array([grad_w0, grad_w1])
sine['npts'] = 30
sine['w0'], sine['w1'] = np.meshgrid(
np.linspace(-1.5, 3, sine['npts']), np.linspace(1, 4, sine['npts']))
# Stack the grid into an array of weight vectors (w_0, w_1).
sine['ws'] = np.stack([sine['w0'].flatten(), sine['w1'].flatten()]).T
sine['error'] = np.array([sine_MSE(w) for w in sine['ws']]).reshape(sine['w0'].shape)
names = ('w0', 'w1')
ind = np.argmin(sine['error'])
sine['grid_best'] = sine['ws'][ind, :]
(sine['dw0'], sine['dw1']) = np.array([grad_sine_MSE(w) for w in sine['ws']]).T
sine['dw0'] = sine['dw0'].reshape(sine['w1'].shape)
sine['dw1'] = sine['dw1'].reshape(sine['w1'].shape)
Batch Gradient Descent
Batch gradient descent with a fixed learning rate $\eta$ (Bishop, Eq. 7.16): $$ w^{(\tau)} = w^{(\tau-1)} - \eta \nabla E\big(w^{(\tau-1)}\big). $$
def gradient_descent(w_0, gradient, learning_rate=1, nepochs=10, epsilon=1e-6):
"""Batch gradient descent with a fixed learning rate.
Args:
w_0: Initial weight vector.
gradient: Function that returns the gradient of the error at w.
learning_rate: Learning rate (step size) eta.
nepochs: Maximum number of iterations.
epsilon: Stop when the norm of the update falls below this threshold.
Returns:
path: Array containing the weight vector at each iteration.
"""
w_old = np.asarray(w_0, dtype=float)
path = [w_old]
for _ in range(nepochs):
w = w_old - learning_rate * gradient(w_old)
path.append(w)
if np.linalg.norm(w - w_old) < epsilon:
break
w_old = w
return np.array(path)
As a reminder, gradient descent on the logistic regression error with $\eta = 10$, starting from $w^{(0)} = (-10, -5)$:
w0 = np.array([-10., -5.])
path = gradient_descent(w0, grad_NLL, learning_rate=10, nepochs=100)
fig_loss = make_loss_curve(path, neg_log_likelihood)
fig_loss.show()
errors = [neg_log_likelihood(w) for w in path]
fig = plot_gradient(cancer['w1'], cancer['w2'], cancer['error'],
cancer['dw1'], cancer['dw2'], scale=2)
add_solution_path(fig, errors, path)
Recap: The Quadratic Approximation
The Hessian, the quadratic approximation of the error function, and the eigenvectors of the Hessian at a stationary point are developed in the Lecture 10 notebook. That notebook also shows that gradient descent on the parabola $f(x) = x^2 + 1$ converges only if $0 < \eta < 1 = 2/f''(x)$; this notebook generalizes that condition to error functions of many variables. We use the quadratic approximation: near a stationary point $w^\star$ with Hessian eigenpairs $(\lambda_i, u_i)$, $$ E(w) \simeq E(w^\star) + \frac{1}{2} \sum_i \lambda_i \alpha_i^2, \qquad \alpha_i = u_i^\top (w - w^\star). $$ The eigenvector plots on the slides "Eigenvectors of the Hessian" are produced by Lecture Demo 5 ("The Quadratic Approximation") in the Lecture 10 notebook.
(Lecture Demo 1) Convergence of Gradient Descent on a Quadratic
Consider gradient descent on the quadratic approximation. In the coordinates $\alpha_i$ the gradient is $\partial E / \partial \alpha_i = \lambda_i \alpha_i$, so each step gives (Bishop, Eqs. 7.26 to 7.29) $$ \alpha_i^{(\tau)} = (1 - \eta \lambda_i)\, \alpha_i^{(\tau-1)} = (1 - \eta \lambda_i)^\tau \alpha_i^{(0)}. $$ The components evolve independently. The iterates converge to $w^\star$ if and only if $|1 - \eta \lambda_i| < 1$ for every $i$, which requires all $\lambda_i > 0$ and $$ 0 < \eta < \frac{2}{\lambda_{\max}}. $$ When $1 - \eta\lambda_i < 0$ the component $\alpha_i$ changes sign at every step, which produces oscillation across the valley.
For every stable learning rate $\eta < 2/\lambda_{\max}$, the factor of the component along $u_{\min}$ satisfies $|1 - \eta\lambda_{\min}| > 1 - 2\lambda_{\min}/\lambda_{\max}$ (Bishop, Eq. 7.30). When the condition number $\kappa = \lambda_{\max}/\lambda_{\min}$ is large, this bound is close to 1, so that component converges slowly whatever $\eta$ we choose. The $u_{\min}$ component is not always the slowest one: for $\eta$ close to $2/\lambda_{\max}$ the factor $|1 - \eta\lambda_{\max}|$ is also close to 1. The step size that minimizes the largest factor is $\eta^\star = 2/(\lambda_{\max} + \lambda_{\min})$, which gives the rate $(\kappa - 1)/(\kappa + 1)$ in every direction.
Lecture: the runs below give the plots on "Issues with Gradient Descent (GD)" ($\eta = 0.1$ and $\eta = 0.9$) and on "Demo: Oscillating convergence".
We study the poorly conditioned quadratic $$ E(w) = \frac{1}{2} w^\top H w, \qquad H = \begin{bmatrix} 0.2 & 0 \\ 0 & 2 \end{bmatrix}, $$ with minimum at $w^\star = 0$. The eigenvectors are the coordinate axes, $\lambda_{\min} = 0.2$, $\lambda_{\max} = 2$, and $\kappa = 10$. Gradient descent converges only for $\eta < 2/\lambda_{\max} = 1$.
H_quad = np.array([[0.2, 0.], [0., 2.]])
def quad(w):
"""Poorly conditioned quadratic E(w) = 0.5 w^T H w (w may be a batch of rows)."""
return 0.5 * np.sum((w @ H_quad) * w, axis=-1)
def grad_quad(w):
"""Gradient of the quadratic, H w."""
return w @ H_quad
w1, w2 = np.meshgrid(np.linspace(-5, 5, 30), np.linspace(-5, 5, 30))
ws = np.hstack([w1.reshape(-1, 1), w2.reshape(-1, 1)])
error = quad(ws).reshape(w1.shape)
contour = go.Contour(x=w1.flatten(), y=w2.flatten(), z=error.flatten(), colorscale='Viridis_r',
contours=dict(start=0, end=20, size=.5))
go.Figure(data=contour)
def plot_quad_path(path, title=None):
"""Plot an optimization path on the contours of the quadratic."""
fig = go.Figure()
fig.add_trace(contour)
fig.add_scatter(x=path[:, 0], y=path[:, 1],
mode='lines+markers', line=dict(color='black', width=2),
marker=dict(size=10, color='black', symbol="arrow-bar-up", angleref="previous"),
name='Optimization Path', showlegend=False)
fig.update_layout(margin=dict(l=5, r=5, t=40 if title else 5, b=5), title=title)
return fig
def report_factors(eta):
"""Print the per-step contraction factor 1 - eta * lambda_i for each eigen-direction."""
for lam_i in np.diag(H_quad):
print(f"lambda = {lam_i:4.1f}: 1 - eta*lambda = {1 - eta * lam_i:6.3f}")
Small learning rate. Starting at $(-4, 0)$, only the low-curvature component is non-zero. With $\eta = 0.1$ its factor is $1 - 0.1 \times 0.2 = 0.98$, so progress along the valley is slow.
eta = 0.1
report_factors(eta)
path = gradient_descent(np.array([-4., 0.]), grad_quad, learning_rate=eta, nepochs=50)
plot_quad_path(path, f"eta = {eta}, start (-4, 0)")
lambda = 0.2: 1 - eta*lambda = 0.980 lambda = 2.0: 1 - eta*lambda = 0.800
Learning rate at the stability limit. Starting at $(-4, -2)$ with $\eta = 1 = 2/\lambda_{\max}$, the factor along $w_2$ is exactly $-1$. The $w_2$ component alternates between $-2$ and $+2$ indefinitely while the $w_1$ component converges with factor $0.8$.
eta = 1.0
report_factors(eta)
path = gradient_descent(np.array([-4., -2.]), grad_quad, learning_rate=eta, nepochs=50)
plot_quad_path(path, f"eta = {eta}, start (-4, -2)")
lambda = 0.2: 1 - eta*lambda = 0.800 lambda = 2.0: 1 - eta*lambda = -1.000
Slightly smaller learning rate. With $\eta = 0.9$ the factor along $w_2$ is $-0.8$, which gives a damped oscillation across the valley, and the factor along $w_1$ is $0.82$.
eta = 0.9
report_factors(eta)
path = gradient_descent(np.array([-4., -2.]), grad_quad, learning_rate=eta, nepochs=50)
plot_quad_path(path, f"eta = {eta}, start (-4, -2)")
lambda = 0.2: 1 - eta*lambda = 0.820 lambda = 2.0: 1 - eta*lambda = -0.800
Optimal fixed learning rate. With $\eta^\star = 2/(\lambda_{\max} + \lambda_{\min}) = 2/2.2 \approx 0.91$, both factors have magnitude $(\kappa - 1)/(\kappa + 1) = 9/11 \approx 0.82$. No fixed learning rate does better on this problem, and the rate approaches 1 as $\kappa$ grows. This limitation motivates momentum and adaptive learning rates.
eta = 2 / (H_quad.max() + H_quad[H_quad > 0].min())
report_factors(eta)
path = gradient_descent(np.array([-4., -2.]), grad_quad, learning_rate=eta, nepochs=50)
plot_quad_path(path, f"eta = {eta:.3f}, start (-4, -2)")
lambda = 0.2: 1 - eta*lambda = 0.818 lambda = 2.0: 1 - eta*lambda = -0.818
(Lecture Demo 2) Gradient Descent with Momentum
Momentum adds a fraction $\mu$ of the previous update to the current one (Bishop, Eq. 7.31): $$ \Delta w^{(\tau-1)} = -\eta \nabla E\big(w^{(\tau-1)}\big) + \mu \Delta w^{(\tau-2)}, \qquad w^{(\tau)} = w^{(\tau-1)} + \Delta w^{(\tau-1)}, $$ with $0 \le \mu < 1$. Substituting for $\Delta w^{(\tau-2)}$ repeatedly, with $\Delta w = 0$ before the first step, gives $$ \Delta w^{(\tau-1)} = -\eta \sum_{k=0}^{\tau-1} \mu^k \, \nabla E\big(w^{(\tau-1-k)}\big), $$ so each step is an exponentially weighted sum of past gradients, and a gradient from $k$ steps ago has weight $\mu^k$. If the gradient is approximately constant (low curvature), the sum is approximately $\nabla E/(1 - \mu)$ and the effective learning rate increases to $\eta/(1 - \mu)$. If the gradient changes sign at every step (high curvature), the sum is approximately $\nabla E/(1 + \mu)$ and the effective learning rate decreases to $\eta/(1 + \mu)$.
Deep learning libraries such as PyTorch maintain a velocity $v^{(\tau)} = \mu v^{(\tau-1)} + \nabla E(w^{(\tau-1)})$ and update $w^{(\tau)} = w^{(\tau-1)} - \eta v^{(\tau)}$. With $\Delta w^{(\tau-1)} = -\eta v^{(\tau)}$ this is the same algorithm when $\eta$ is fixed. The implementation below uses Bishop's form and stores $\Delta w$ in the variable dw.
Lecture: run from the slide "Demo: Momentum".
def gd_momentum(w_0, gradient, learning_rate=1, nepochs=10, epsilon=1e-6, momentum=0.9):
"""Gradient descent with momentum, Bishop Eq. 7.31.
Args:
w_0: Initial weight vector.
gradient: Function that returns the gradient of the error at w.
learning_rate: Learning rate eta.
nepochs: Maximum number of iterations.
epsilon: Stop when the norm of the update falls below this threshold.
momentum: Momentum parameter mu, with 0 <= mu < 1.
Returns:
path: Array containing the weight vector at each iteration.
"""
w_old = np.asarray(w_0, dtype=float)
path = [w_old]
dw = np.zeros_like(w_old) # previous update Delta w
for _ in range(nepochs):
dw = momentum * dw - learning_rate * gradient(w_old)
w = w_old + dw
path.append(w)
if np.linalg.norm(w - w_old) < epsilon:
break
w_old = w
return np.array(path)
A quick check that the velocity form used by PyTorch gives the same iterates as Bishop's form for a fixed learning rate:
def gd_velocity(w_0, gradient, learning_rate=1, nepochs=10, momentum=0.9):
"""Momentum in the velocity form: v <- mu v + grad, w <- w - eta v."""
w = np.asarray(w_0, dtype=float)
path, v = [w], np.zeros_like(w)
for _ in range(nepochs):
v = momentum * v + gradient(w)
w = w - learning_rate * v
path.append(w)
return np.array(path)
p1 = gd_momentum(np.array([-4., -2.]), grad_quad, learning_rate=0.5, momentum=0.8, nepochs=20, epsilon=0)
p2 = gd_velocity(np.array([-4., -2.]), grad_quad, learning_rate=0.5, momentum=0.8, nepochs=20)
print("max difference between the two forms:", np.abs(p1 - p2).max())
max difference between the two forms: 0.0
We apply momentum to the toy quadratic function from the previous section, starting at $(-4, -2)$ with $\eta = 1$ and $\mu = 0.4$. Without momentum this learning rate lies exactly on the stability limit.
eta, mu = 1.0, 0.4
path = gd_momentum(np.array([-4., -2.]), grad_quad, learning_rate=eta, momentum=mu, nepochs=50)
plot_quad_path(path, f"Momentum: eta = {eta}, mu = {mu}")
Momentum on the logistic regression error, starting from $w^{(0)} = (0, 0)$. Along the long, flat valley of this error surface the gradient points in a consistent direction, so momentum accelerates the iterates; they overshoot the minimum and then return. Set momentum=0 to compare with plain gradient descent at the same learning rate.
w0 = np.array([0., 0.])
path = gd_momentum(w0, grad_NLL, learning_rate=1., momentum=0.9, nepochs=100)
fig_loss = make_loss_curve(path, neg_log_likelihood)
fig_loss.show()
errors = [neg_log_likelihood(w) for w in path]
fig = plot_gradient(cancer['w1'], cancer['w2'], cancer['error'],
cancer['dw1'], cancer['dw2'], scale=2)
add_solution_path(fig, errors, path)
Momentum can also be applied to the non-convex sine regression problem. Set momentum=0 to compare with plain gradient descent at the same learning rate.
w0 = np.array([1.5, 2.])
path = gd_momentum(w0, grad_sine_MSE, learning_rate=.1, momentum=0.5, nepochs=50)
fig_loss = make_loss_curve(path, sine_MSE)
fig_loss.show()
errors = [sine_MSE(w) for w in path]
fig = plot_gradient(sine['w0'], sine['w1'], sine['error'],
sine['dw0'], sine['dw1'], scale=.1, names=names)
add_solution_path(fig, errors, path)
(Lecture Demo 3) Adaptive Learning Rates
The optimal learning rate depends on the local curvature, which can differ between directions in weight space. Adaptive methods maintain a separate effective learning rate for each parameter $w_i$ (Bishop, Section 7.3.3). Let $g_i = \partial E / \partial w_i$ denote the current gradient component.
AdaGrad accumulates the sum of squared gradients (Eqs. 7.39 and 7.40): $$ r_i^{(\tau)} = r_i^{(\tau-1)} + g_i^2, \qquad w_i^{(\tau)} = w_i^{(\tau-1)} - \frac{\eta}{\sqrt{r_i^{(\tau)}} + \delta}\, g_i. $$ AdaGrad therefore takes smaller steps along parameters whose gradients have been large in the past. Because $r_i$ only grows, the effective learning rate decreases monotonically and can become too small late in training.
RMSProp replaces the sum with an exponentially weighted moving average (Eqs. 7.41 and 7.42): $$ r_i^{(\tau)} = \beta r_i^{(\tau-1)} + (1 - \beta) g_i^2, \qquad w_i^{(\tau)} = w_i^{(\tau-1)} - \frac{\eta}{\sqrt{r_i^{(\tau)}} + \delta}\, g_i. $$
Adam combines RMSProp with a moving average of the gradient itself, which acts as momentum (Eqs. 7.43 to 7.47): \begin{align*} s_i^{(\tau)} &= \beta_1 s_i^{(\tau-1)} + (1 - \beta_1) g_i, & r_i^{(\tau)} &= \beta_2 r_i^{(\tau-1)} + (1 - \beta_2) g_i^2, \\ \widehat{s}_i^{(\tau)} &= \frac{s_i^{(\tau)}}{1 - \beta_1^\tau}, & \widehat{r}_i^{(\tau)} &= \frac{r_i^{(\tau)}}{1 - \beta_2^\tau}, \\ w_i^{(\tau)} &= w_i^{(\tau-1)} - \eta \frac{\widehat{s}_i^{(\tau)}}{\sqrt{\widehat{r}_i^{(\tau)}} + \delta}. \end{align*} The factors $1/(1 - \beta^\tau)$ correct the bias towards zero caused by initializing $s_i^{(0)} = r_i^{(0)} = 0$. In all three methods $\delta$ is a small constant for numerical stability. Bishop gives $\beta_2 = 0.99$ as a typical value. The original Adam paper (Kingma and Ba, 2014) and the PyTorch default use $\beta_1 = 0.9$ and $\beta_2 = 0.999$, which we use below.
Note that Adam's first moment $s$ is an exponential moving average: the weights $(1 - \beta_1)\beta_1^k$ on past gradients sum to one, so for a constant gradient $s$ approaches $\nabla E$ itself. Adam also divides by $\sqrt{\widehat{r}_i}$, so its step in each coordinate has magnitude of about $\eta$ whatever the scale of the gradient; the first step is exactly $-\eta\, \mathrm{sign}(g_i)$ up to $\delta$.
Lecture: run from the slide "Demo: Adam".
def gd_adam(w_0, gradient, learning_rate=lambda t: 0.1, nepochs=100,
delta=1e-8, beta1=0.9, beta2=0.999):
"""Batch gradient descent with the Adam update.
Args:
w_0: Initial weight vector.
gradient: Function that returns the gradient of the error at w.
learning_rate: Function of the number of completed steps (0, 1, 2, ...) that returns eta.
nepochs: Number of iterations.
delta: Small constant added to the denominator for numerical stability.
beta1: Decay rate of the moving average of the gradient (first moment).
beta2: Decay rate of the moving average of the squared gradient (second moment).
Returns:
path: Array containing the weight vector at each iteration.
"""
w = np.asarray(w_0, dtype=float)
path = [w]
s = np.zeros_like(w) # first moment (momentum)
r = np.zeros_like(w) # second moment (RMSProp)
for t in range(1, nepochs + 1):
g = gradient(w)
s = beta1 * s + (1 - beta1) * g
r = beta2 * r + (1 - beta2) * g**2
s_hat = s / (1 - beta1**t) # bias correction
r_hat = r / (1 - beta2**t)
w = w - learning_rate(t - 1) * s_hat / (np.sqrt(r_hat) + delta) # eta^(tau-1), as in the other methods
path.append(w)
return np.array(path)
w0 = np.array([-1, 2.5])
path = gd_adam(w0, grad_sine_MSE, learning_rate=lambda t: 0.1, nepochs=100)
fig_loss = make_loss_curve(path, sine_MSE)
fig_loss.show()
errors = [sine_MSE(w) for w in path]
fig = plot_gradient(sine['w0'], sine['w1'], sine['error'],
sine['dw0'], sine['dw1'], scale=.1, names=names)
add_solution_path(fig, errors, path)
Separating the Update Rule from the Optimization Loop
The algorithms above differ only in how the gradient is turned into an update. We therefore write each update rule as an optimizer object with a step(w, g) method, which returns the new weight vector given the current weights and gradient, and write a single optimization loop. Each optimizer takes the learning rate as a function of the number of completed steps $t = 0, 1, 2, \dots$, so the step that produces $w^{(\tau)}$ uses $\eta^{(\tau-1)}$ (Bishop, Eq. 7.35), as torch.optim.lr_scheduler does. MomentumOptimizer applies the current learning rate to the new gradient only, as in Bishop's form; PyTorch instead multiplies the whole velocity by the current learning rate, so the two differ while the learning rate is changing. The same optimizer objects are used with stochastic gradient descent below. This mirrors the design of torch.optim.
class VanillaOptimizer:
"""Gradient descent update w <- w - eta_t g."""
def __init__(self, learning_rate=lambda t: 0.1):
self.learning_rate = learning_rate
self.t = 0
def step(self, w, g):
w = w - self.learning_rate(self.t) * g
self.t += 1
return w
class MomentumOptimizer:
"""Momentum: v <- mu v - eta_t g, w <- w + v (v stores the previous update Delta w)."""
def __init__(self, learning_rate=lambda t: 0.1, momentum=0.9):
self.learning_rate = learning_rate
self.momentum = momentum
self.t = 0
self.v = None
def step(self, w, g):
if self.v is None:
self.v = np.zeros_like(w)
self.v = self.momentum * self.v - self.learning_rate(self.t) * g
self.t += 1
return w + self.v
class AdaGradOptimizer:
"""AdaGrad: r <- r + g^2, w <- w - eta_t g / (sqrt(r) + delta)."""
def __init__(self, learning_rate=lambda t: 0.1, delta=1e-8):
self.learning_rate = learning_rate
self.delta = delta
self.t = 0
self.r = None
def step(self, w, g):
if self.r is None:
self.r = np.zeros_like(w)
self.r = self.r + g**2
w = w - self.learning_rate(self.t) * g / (np.sqrt(self.r) + self.delta)
self.t += 1
return w
class RMSPropOptimizer:
"""RMSProp: r <- beta r + (1 - beta) g^2, w <- w - eta_t g / (sqrt(r) + delta)."""
def __init__(self, learning_rate=lambda t: 0.01, beta=0.9, delta=1e-8):
self.learning_rate = learning_rate
self.beta = beta
self.delta = delta
self.t = 0
self.r = None
def step(self, w, g):
if self.r is None:
self.r = np.zeros_like(w)
self.r = self.beta * self.r + (1 - self.beta) * g**2
w = w - self.learning_rate(self.t) * g / (np.sqrt(self.r) + self.delta)
self.t += 1
return w
class AdamOptimizer:
"""Adam with bias-corrected first and second moment estimates."""
def __init__(self, learning_rate=lambda t: 0.1, beta1=0.9, beta2=0.999, delta=1e-8):
self.learning_rate = learning_rate
self.beta1 = beta1
self.beta2 = beta2
self.delta = delta
self.t = 0
self.s = None
self.r = None
def step(self, w, g):
if self.s is None:
self.s = np.zeros_like(w) # first moment (momentum)
self.r = np.zeros_like(w) # second moment (RMSProp)
eta = self.learning_rate(self.t) # schedule at the number of completed steps
self.t += 1
self.s = self.beta1 * self.s + (1 - self.beta1) * g
self.r = self.beta2 * self.r + (1 - self.beta2) * g**2
s_hat = self.s / (1 - self.beta1**self.t) # bias correction
r_hat = self.r / (1 - self.beta2**self.t)
return w - eta * s_hat / (np.sqrt(r_hat) + self.delta)
def gd(w_0, gradient, optimizer, nepochs=100, epsilon=1e-8):
"""Batch gradient descent with a pluggable update rule.
Args:
w_0: Initial weight vector.
gradient: Function that returns the gradient of the error at w.
optimizer: Object with a step(w, g) method that returns the updated weights.
nepochs: Maximum number of iterations.
epsilon: Stop when the norm of the update falls below this threshold.
Returns:
path: Array containing the weight vector at each iteration.
"""
w = np.asarray(w_0, dtype=float)
path = [w]
for _ in range(nepochs):
w = optimizer.step(w, gradient(w))
path.append(w)
if np.linalg.norm(w - path[-2]) < epsilon:
break
return np.array(path)
The loss curves below compare the five update rules on the sine regression problem from the same initial point. The learning rates are chosen separately for each method, because the scale of the update differs between methods.
w0 = np.array([-1, 2.5])
optimizers = {
'Gradient descent': VanillaOptimizer(learning_rate=lambda t: 0.1),
'Momentum': MomentumOptimizer(learning_rate=lambda t: 0.01, momentum=0.9),
'AdaGrad': AdaGradOptimizer(learning_rate=lambda t: 0.1),
'RMSProp': RMSPropOptimizer(learning_rate=lambda t: 0.1),
'Adam': AdamOptimizer(learning_rate=lambda t: 0.1),
}
fig = go.Figure()
paths = {}
for name, opt in optimizers.items():
paths[name] = gd(w0, grad_sine_MSE, opt, nepochs=100)
fig.add_scatter(y=[sine_MSE(w) for w in paths[name]], mode='lines', name=name,
line=dict(width=3))
fig.update_layout(xaxis_title='Iteration', yaxis_title='Error',
margin=dict(l=20, r=20, t=20, b=20))
fig
path = paths['Adam']
errors = [sine_MSE(w) for w in path]
fig = plot_gradient(sine['w0'], sine['w1'], sine['error'],
sine['dw0'], sine['dw1'], scale=.1, names=names)
add_solution_path(fig, errors, path)
(Lecture Demo 4) Learning Rate Schedules
In practice a larger learning rate is used at the start of training and is then reduced, so that $\eta$ becomes a function of the iteration index $\tau$ (Bishop, Eq. 7.35): $$ w^{(\tau)} = w^{(\tau-1)} - \eta^{(\tau-1)} \nabla E\big(w^{(\tau-1)}\big). $$ Each schedule below is implemented as a function that returns a callable $\tau \mapsto \eta^{(\tau)}$. Here $\tau$ counts iterations (gradient steps), and $T$ is the total number of iterations.
Lecture: the comparison at the end of this section is the figure on "Demo: Learning Rate Schedules".
Linear Schedule
$$ \eta^{(\tau)} = \left(1 - \frac{\tau}{T}\right) \eta_{\text{start}} + \frac{\tau}{T} \eta_{\text{end}}, $$ after which the learning rate is held at $\eta_{\text{end}}$ (Bishop, Eq. 7.36).
def linear_learning_rate(initial_lr, T, end_lr=0):
"""Linear decay from initial_lr to end_lr over T iterations, then constant."""
return lambda t: initial_lr * (1 - min(1, t / T)) + end_lr * min(1, t / T)
def plot_schedule(lr, T=100):
fig = px.line(x=np.arange(T), y=[lr(t) for t in range(T)], markers=True,
labels={'x': 'Iteration', 'y': 'Learning Rate'})
fig.update_traces(line_width=4)
fig.update_layout(margin=dict(l=20, r=20, t=20, b=20))
return fig
plot_schedule(linear_learning_rate(initial_lr=0.1, T=100, end_lr=0.01))
Power-Law Schedule
$$ \eta^{(\tau)} = \eta_{\text{start}} \left(1 + \frac{\tau}{s}\right)^{-c} $$ (Bishop, Eq. 7.37), where $s$ sets the time scale of the decay and $c$ its exponent.
def power_law_learning_rate(initial_lr, s=1, c=1):
"""Power-law decay: initial_lr * (1 + t/s)^(-c)."""
return lambda t: initial_lr * (1 + t / s) ** (-c)
plot_schedule(power_law_learning_rate(initial_lr=0.1, s=10, c=1))
Exponential and Step Decay
The exponential schedule (Bishop, Eq. 7.38) is $$ \eta^{(\tau)} = \eta_{\text{start}}\, c^{\tau/s}, \qquad 0 < c < 1. $$ The slides use the common variant step decay, which holds the learning rate constant for $s$ iterations at a time: $\eta^{(\tau)} = \eta_{\text{start}}\, c^{\lfloor \tau/s \rfloor}$.
def exponential_learning_rate(initial_lr, c=0.9, s=10, staircase=False):
"""Exponential decay initial_lr * c^(t/s), or step decay c^floor(t/s) if staircase."""
if staircase:
return lambda t: initial_lr * c ** np.floor(t / s)
return lambda t: initial_lr * c ** (t / s)
plot_schedule(exponential_learning_rate(initial_lr=0.1, c=0.9, s=10, staircase=True))
Cosine Schedule with Linear Warmup
The cosine schedule decays from $\eta_{\max}$ to $\eta_{\min}$ over $T$ iterations: $$ \eta^{(\tau)} = \eta_{\min} + \frac{1}{2}\left(\eta_{\max} - \eta_{\min}\right)\left(1 + \cos\left(\frac{\pi \tau}{T}\right)\right). $$ Large models are commonly trained with a short warmup phase before the decay, during which the learning rate increases linearly from 0 to $\eta_{\max}$ over the first $T_w$ iterations, $\eta^{(\tau)} = \eta_{\max}\, \tau / T_w$ for $\tau < T_w$, and then follows the cosine schedule over the remaining $T - T_w$ iterations. Warmup avoids large, poorly directed steps at the start of training, when the optimizer state (for example the moving averages in Adam) has not yet adapted.
def cosine_learning_rate(T, max_lr, min_lr=0):
"""Cosine decay from max_lr to min_lr over T iterations."""
return lambda t: min_lr + 0.5 * (max_lr - min_lr) * (1 + np.cos(np.pi * min(t, T) / T))
def warmup_cosine_learning_rate(T, T_warmup, max_lr, min_lr=0):
"""Linear warmup over T_warmup iterations followed by cosine decay to min_lr at T."""
cosine = cosine_learning_rate(T - T_warmup, max_lr, min_lr)
return lambda t: max_lr * t / T_warmup if t < T_warmup else cosine(t - T_warmup)
plot_schedule(warmup_cosine_learning_rate(T=100, T_warmup=10, max_lr=0.1, min_lr=0.01))
Comparing Schedules
T = 100
schedules = {
'Linear': linear_learning_rate(initial_lr=0.1, T=T, end_lr=0.01),
'Power law': power_law_learning_rate(initial_lr=0.1, s=10, c=1),
'Step decay': exponential_learning_rate(initial_lr=0.1, c=0.9, s=5, staircase=True),
'Cosine': cosine_learning_rate(T=T, max_lr=0.1, min_lr=0.01),
'Warmup + cosine': warmup_cosine_learning_rate(T=T, T_warmup=10, max_lr=0.1, min_lr=0.01),
}
t = np.arange(T)
fig = go.Figure()
for name, lr in schedules.items():
fig.add_scatter(x=t, y=[lr(tau) for tau in t], mode='lines+markers', name=name,
line=dict(width=4))
fig.update_layout(margin=dict(l=20, r=20, t=20, b=20),
xaxis_title='Iteration', yaxis_title='Learning Rate', height=600)
(Lecture Demo 5) Stochastic Gradient Descent
The error function is an average over the training data, $$ E(w) = \frac{1}{N} \sum_{n=1}^N E_n(w), $$ so each batch gradient step evaluates $N$ per-example gradients. For a model with $W$ parameters each of these costs $O(W)$, so a batch step costs $O(NW)$. Stochastic gradient descent (Bishop, Section 7.2.3) replaces $\nabla E$ with the gradient of the error on a single data point $n$: $$ w^{(\tau)} = w^{(\tau-1)} - \eta \nabla E_n\big(w^{(\tau-1)}\big). $$ If $n$ is drawn uniformly from $\{1, \dots, N\}$, the stochastic gradient is an unbiased estimate of the batch gradient: $$ \mathbb{E}_n\big[\nabla E_n(w)\big] = \frac{1}{N} \sum_{n=1}^N \nabla E_n(w) = \nabla E(w). $$ Bishop defines $E(w) = \sum_n E_n(w)$ as a sum rather than an average. The two conventions differ by the constant factor $N$, which is absorbed into the learning rate.
The slide "Stochastic Gradient Descent" samples a random index at each step. In practice the data are shuffled at the start of each epoch and then visited in order, so each point is used exactly once per epoch ("Stochastic Gradient Descent (Shuffling)"). The implementation below runs for a fixed number of epochs and uses a seeded random number generator so that the results are reproducible.
def sgd(w_0, x, t, gradient, optimizer, nepochs=2, seed=0):
"""Stochastic gradient descent with one data point per step.
Args:
w_0: Initial weight vector.
x, t: Training inputs and targets.
gradient: Function gradient(w, x_n, t_n) returning the gradient on the given data.
optimizer: Object with a step(w, g) method that returns the updated weights.
nepochs: Number of passes over the data.
seed: Seed for the random shuffling of the data.
Returns:
path: Array containing the weight vector after each step.
"""
rng = np.random.default_rng(seed)
w = np.asarray(w_0, dtype=float)
path = [w]
for _ in range(nepochs):
for n in rng.permutation(len(t)): # reshuffle every epoch
w = optimizer.step(w, gradient(w, x[n], t[n]))
path.append(w)
return np.array(path)
SGD with Adam on the sine regression problem. Each epoch consists of $N = 200$ steps, each using a single data point. The path is noisy because each step follows the gradient of a single term.
w0 = np.array([-1, 2.5])
path = sgd(w0, sine['x'], sine['y'], grad_sine_MSE,
AdamOptimizer(learning_rate=lambda t: 0.1), nepochs=2)
errors = [sine_MSE(w) for w in path]
fig = plot_gradient(sine['w0'], sine['w1'], sine['error'],
sine['dw0'], sine['dw1'], scale=.1, names=names)
fig = add_solution_path(fig, errors, path)
fig.show()
fig_loss = make_loss_curve(path, sine_MSE)
fig_loss.show()
SGD with Adam and a cosine learning rate schedule on the logistic regression problem, for one epoch of $N = 569$ steps. The decreasing learning rate reduces the noise in the iterates as training proceeds. This is the figure on "Demo: Stochastic Gradient Descent".
w0 = np.array([-1., 2.])
nepochs = 1
coslr = cosine_learning_rate(T=nepochs * len(cancer['t']), max_lr=0.1, min_lr=0.01)
path = sgd(w0, cancer['x'], cancer['t'], grad_NLL,
AdamOptimizer(learning_rate=coslr), nepochs=nepochs)
errors = [neg_log_likelihood(w) for w in path]
fig = plot_gradient(cancer['w1'], cancer['w2'], cancer['error'],
cancer['dw1'], cancer['dw2'], scale=1)
fig = add_solution_path(fig, errors, path)
fig.show()
fig_loss = make_loss_curve(path, neg_log_likelihood)
fig_loss.show()
(Lecture Demo 6) Mini-Batch Stochastic Gradient Descent
A mini-batch $\mathcal{B}$ of $B$ data points gives the gradient estimate $$ \nabla E_{\mathcal{B}}(w) = \frac{1}{B} \sum_{n \in \mathcal{B}} \nabla E_n(w), $$ which is also unbiased and costs $O(BW)$. For $B = 1$ this is SGD, and for $B = N$ it is batch gradient descent. If $\sigma$ is the standard deviation of a single-point gradient component, the standard deviation of the mini-batch estimate is approximately $\sigma/\sqrt{B}$ (Bishop, Section 7.2.4 and Exercise 7.8). Increasing the batch size by a factor of 100 therefore reduces the noise only by a factor of 10.
Mini-batch SGD shuffles the data at the start of each epoch and then takes successive blocks of $B$ points (Bishop, Algorithm 7.2). Each epoch therefore has $\lceil N/B \rceil$ steps (for the Slido question, $N/2$ steps when $B = 2$ and $N$ is even). With $N = 200$ and $B = 32$ each epoch has $\lceil 200/32 \rceil = 7$ steps, the last of which uses the remaining 8 points.
def mb_sgd(w_0, x, t, gradient, optimizer, batch_size=32, nepochs=2, seed=0):
"""Mini-batch stochastic gradient descent.
Args:
w_0: Initial weight vector.
x, t: Training inputs and targets.
gradient: Function gradient(w, x_batch, t_batch) returning the average gradient on the batch.
optimizer: Object with a step(w, g) method that returns the updated weights.
batch_size: Number of data points B in each mini-batch.
nepochs: Number of passes over the data.
seed: Seed for the random shuffling of the data.
Returns:
path: Array containing the weight vector after each step.
"""
rng = np.random.default_rng(seed)
w = np.asarray(w_0, dtype=float)
path = [w]
for _ in range(nepochs):
order = rng.permutation(len(t)) # reshuffle every epoch
for start in range(0, len(t), batch_size):
batch = order[start:start + batch_size]
w = optimizer.step(w, gradient(w, x[batch], t[batch]))
path.append(w)
return np.array(path)
w0 = np.array([-1, 2.5])
path = mb_sgd(w0, sine['x'], sine['y'], grad_sine_MSE,
AdamOptimizer(learning_rate=lambda t: 0.1),
batch_size=32, nepochs=10)
errors = [sine_MSE(w) for w in path]
fig = plot_gradient(sine['w0'], sine['w1'], sine['error'],
sine['dw0'], sine['dw1'], scale=.1, names=names)
fig = add_solution_path(fig, errors, path)
fig.show()
fig_loss = make_loss_curve(path, sine_MSE)
fig_loss.show()
Mini-batch SGD with Adam on the logistic regression problem, with $B = 32$ and 5 epochs. This is the figure on "Demo: Mini-batch Gradient Descent".
w0 = np.array([0., -4.])
path = mb_sgd(w0, cancer['x'], cancer['t'], grad_NLL,
AdamOptimizer(learning_rate=lambda t: 0.2),
batch_size=32, nepochs=5)
errors = [neg_log_likelihood(w) for w in path]
fig = plot_gradient(cancer['w1'], cancer['w2'], cancer['error'],
cancer['dw1'], cancer['dw2'], scale=2)
fig = add_solution_path(fig, errors, path)
fig.show()
fig_loss = make_loss_curve(path, neg_log_likelihood)
fig_loss.show()