import plotly.io as pio; pio.renderers.default = "notebook_connected"
Lecture 10: Optimization and Gradient Descent (CS 189, Fall 2026)
This notebook accompanies Lecture 10 and Chapter 7 of Bishop and Bishop, Deep Learning: Foundations and Concepts. It covers:
- simple one-dimensional optimization problems and convexity,
- the error surface of two small models (logistic regression and a non-linear sine regression model),
- the gradient of the error function,
- batch gradient descent and loss curves, and
- the Hessian, the quadratic approximation of the error function, and the eigenvectors of the Hessian at a stationary point.
The convergence analysis of gradient descent, momentum, and Adam are covered in the Lecture 11 notebook.
Sections used in lecture are marked in their headings. (Lecture Demo n) is run during the lecture from the matching "Demo" slide, and (Lecture Slido) produces the plots for a Slido question. Other sections produce plots that appear on the slides marked "Plots in demo notebook."
| Notebook section | Lecture 10 slide |
|---|---|
| Optimization Basics | Examples Optimization Problems (plots) |
| (Lecture Demo 1) The Error Surface | Demo: Understanding the Error Surface; Points on the Error Surface (plots) |
| (Lecture Slido) Match the Model with the Error Function | Slido: Match the model with the loss |
| (Lecture Demo 2) The Gradient of the Error Function | Demo: The Gradient of the Error Function |
| (Lecture Demo 3) Batch Gradient Descent | Demo: Batch Gradient Descent; Loss (Error) Curves |
| (Lecture Demo 4) Gradient Descent on a Parabola | Demo: Gradient Descent on a Parabola; the Slido question "Assuming a simple quadratic equation (e.g., x^2) is gradient descent guaranteed to converge?"; the result is the slide "What are the conditions for convergence?" |
| (Lecture Demo 5) Second-Order Structure of the Error Function | Hessian Exercise (plots) |
| (Lecture Demo 5) The Quadratic Approximation | Demo: Taylor Expansions and the Hessian; Eigenvectors of the Hessian (plots) |
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])
# Use the same scale for both weight axes on the contour plot. Otherwise one axis is
# stretched and the gradient no longer appears perpendicular to the contours.
# (The 3D scene takes no 2D axis, so the contour plot uses xaxis/yaxis.)
fig.update_yaxes(scaleanchor='x', scaleratio=1, constrain='domain', row=1, col=2)
fig.update_xaxes(constrain='domain', row=1, col=2)
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
Optimization Basics
Most learning problems can be written as the optimization problem $$ w^\star = \arg\min_{w \in \Theta} E(w), $$ where $E(w)$ is the error (objective) function and $\Theta$ is the constraint set. We begin with three one-dimensional examples. These are the plots on the slide "Examples Optimization Problems."
Problem 1. Unconstrained minimization of a quadratic: $$ \arg \min_{w \in \mathbb{R}} \; w^2 - 3w + 4. $$ Setting the derivative $2w - 3$ to zero gives the unique minimizer $w^\star = 3/2$, with minimum value $7/4$.
w = np.linspace(-2, 7, 100)
f = lambda w: w**2 - 3 * w + 4
fig = px.line(x=w, y=f(w), labels={'x': 'w', 'y': 'f(w)'})
fig.update_traces(line_width=5)
fig
Problem 2. The same objective restricted to the integers: $$ \arg \min_{w \in \mathbb{Z}} \; w^2 - 3w + 4. $$ The minimum value $2$ is attained at both $w = 1$ and $w = 2$, so the minimizer is not unique. Because the domain is discrete, the condition $f'(w) = 0$ can no longer be used to locate the solution.
w = np.arange(-2, 7)
f = lambda w: w**2 - 3 * w + 4
fig = px.scatter(x=w, y=f(w), labels={'x': 'w', 'y': 'f(w)'})
fig.update_traces(marker_size=7)
fig
Problem 3. Unconstrained minimization of a quartic: $$ \arg \min_{w \in \mathbb{R}} \; w^4 - 5 w^2 + w + 4. $$ The derivative $4w^3 - 10w + 1$ has three real roots, so the function has three stationary points: a local minimum at $w \approx 1.53$ ($f \approx -0.69$), a local maximum at $w \approx 0.10$, and the global minimum at $w \approx -1.63$ ($f \approx -3.86$). A method that only uses local information, such as gradient descent, may converge to either minimum depending on where it starts.
w = np.linspace(-3, 3, 100)
f = lambda w: w**4 - 5 * w**2 + w + 4
fig = px.line(x=w, y=f(w), labels={'x': 'w', 'y': 'f(w)'})
fig.update_traces(line_width=5)
fig.update_layout(margin=dict(l=20, r=20, t=20, b=20))
fig
# The stationary points are the real roots of f'(w) = 4w^3 - 10w + 1.
stationary = np.sort(np.roots([4, 0, -10, 1]).real)
for w_s in stationary:
print(f"w = {w_s:6.3f}, f(w) = {f(w_s):6.3f}, f''(w) = {12 * w_s**2 - 10:6.2f}")
w = -1.629, f(w) = -3.855, f''(w) = 21.84 w = 0.100, f(w) = 4.050, f''(w) = -9.88 w = 1.529, f(w) = -0.695, f''(w) = 18.04
Convexity
A set $\Theta$ is convex if for any $w_1, w_2 \in \Theta$ and any $0 \le t \le 1$, the point $t w_1 + (1-t) w_2$ also lies in $\Theta$. That is, the set contains the line segment between any two of its points.
A function $f$ defined on a convex set $\Theta$ is convex if for any $w_1, w_2 \in \Theta$ and any $0 \le t \le 1$, $$ f\big(t w_1 + (1-t) w_2\big) \le t f(w_1) + (1-t) f(w_2). $$ Geometrically, the secant line between any two points on the graph of $f$ lies on or above the graph. For a convex function every local minimum is a global minimum, and when $f$ is twice differentiable, convexity is equivalent to the Hessian being positive semi-definite everywhere.
The plots below show a secant line for the convex function of Problem 1 and for the non-convex function of Problem 3. In the second case the secant line lies below the graph near $w = 0$.
def plot_secant(w, f, x1, x2):
"""Plot a function together with the secant line between x1 and x2."""
fig = px.line(x=w, y=f(w), labels={'x': 'w', 'y': 'f(w)'})
fig.add_scatter(x=[x1, x2], y=[f(x1), f(x2)],
mode='markers+lines', name='Secant Line',
marker=dict(size=20, color='green'),
line=dict(color='green', dash="dash"))
fig.update_traces(line_width=5)
fig.update_layout(showlegend=False)
return fig
w = np.linspace(-2, 7, 100)
f = lambda w: w**2 - 3 * w + 4
fig_convex = plot_secant(w, f, -1, 5)
fig_convex.update_layout(title="Convex Function")
w = np.linspace(-3, 3, 100)
f = lambda w: w**4 - 5 * w**2 + w + 4
fig_nonconvex = plot_secant(w, f, -1, 1)
fig_nonconvex.update_layout(title="Non-convex Function")
fig = make_plot_grid([fig_convex, fig_nonconvex], 1, 2)
fig.update_layout(showlegend=False)
(Lecture Demo 1) The Error Surface
Each setting of the parameters $w$ corresponds to a point on the error surface $E(w)$ sitting over weight space (Bishop, Section 7.1). When there are only two parameters we can evaluate $E(w)$ on a grid and plot the surface directly. We do this for two models: logistic regression, whose error function is convex, and a non-linear regression model, whose error function is not.
Lecture: run from the slide "Demo: Understanding the Error Surface." The surface plots also appear on "Points on the Error Surface."
Logistic Regression on the Breast Cancer Data
We use two features of the breast cancer dataset in sklearn, "mean radius" and "mean texture", each standardized to zero mean and unit variance. The target is $t = 1$ for a benign tumor and $t = 0$ for a malignant tumor.
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)
display(cancer_df.head())
px.scatter(cancer_df, x='mean radius', y='mean texture', color='diagnosis', opacity=0.7,
color_discrete_map=DIAGNOSIS_COLORS)
The dataset: (569, 4)
| mean radius | mean texture | target | diagnosis | |
|---|---|---|---|---|
| 0 | 1.097064 | -2.073335 | 0 | malignant |
| 1 | 1.829821 | -0.353632 | 0 | malignant |
| 2 | 1.579888 | 0.456187 | 0 | malignant |
| 3 | -0.768909 | 0.253732 | 0 | malignant |
| 4 | 1.750297 | -1.151816 | 0 | malignant |
The two classes are approximately separated by a straight line, so a linear decision boundary is a reasonable model.
The Logistic Regression Model
The logistic regression model is $$ p(t=1 \mid x) = \sigma(w^\top x + b) = \frac{1}{1 + e^{-(w^\top x + b)}}. $$ So that the error surface can be plotted in two dimensions, we drop the bias $b$ and keep the two weights $w_1$ and $w_2$: $$ p(t=1 \mid x) = \sigma(w_1 x_1 + w_2 x_2). $$ Because the features are standardized, the data are centered at the origin, and a decision boundary that passes through the origin is a reasonable approximation.
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)
The following plots show the decision boundary and the predicted probability for two parameter settings.
guesses = np.array([[-1., -1.],
[-2., -4.]])
figs = [plot_lr_predictions(w, cancer_df) for w in guesses]
figs[1].update_traces(showlegend=False)
fig = make_plot_grid(figs, 1, 2)
fig.update_layout(height=600,
xaxis_range=figs[0].layout.xaxis.range, yaxis_range=figs[0].layout.yaxis.range,
xaxis2_range=figs[1].layout.xaxis.range, yaxis2_range=figs[1].layout.yaxis.range)
fig
The Cross-Entropy Error Surface
Let $y_n = \sigma(w^\top x_n)$ denote the predicted probability that $t_n = 1$. The average negative log-likelihood (cross-entropy) error function is
$$
E(w) = -\frac{1}{N} \sum_{n=1}^N \Big( t_n \ln y_n + (1 - t_n) \ln (1 - y_n) \Big).
$$
Writing $z_n = w^\top x_n$ and using $\ln \sigma(z) = z - \ln(1 + e^{z})$ and $\ln(1 - \sigma(z)) = -\ln(1 + e^{z})$, this simplifies to
$$
E(w) = \frac{1}{N} \sum_{n=1}^N \Big( \ln\big(1 + e^{z_n}\big) - t_n z_n \Big).
$$
The function $\ln(1 + e^{z})$ is evaluated with np.logaddexp(0, z), which avoids overflow for large $z$.
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)
We store the data, the grid of weights, and the quantities derived from them in a dictionary.
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
cancer['error'] = np.array([neg_log_likelihood(w) for w in cancer['ws']])
cancer['error'] = cancer['error'].reshape(cancer['w1'].shape)
fig = plot_loss(cancer['w1'], cancer['w2'], cancer['error'])
for i, g in enumerate(guesses):
fig.add_scatter3d(x=[g[0]], y=[g[1]], z=[neg_log_likelihood(g)],
mode='markers', marker=dict(size=10, color=colors[i + 2]),
name=f'w = {g}', legendgroup=str(i), row=1, col=1)
fig.add_scatter(x=[g[0]], y=[g[1]], mode='markers',
marker=dict(size=10, color=colors[i + 2]),
name=f'w = {g}', legendgroup=str(i), showlegend=False,
row=1, col=2)
fig
Grid Search
With the error evaluated on the grid, the simplest estimate of the minimizer is the grid point with the lowest error. This approach requires $K^D$ error evaluations for a grid with $K$ points per dimension in $D$ dimensions, so it is only practical for very small $D$.
best_ind = np.argmin(cancer['error'])
cancer['grid_best'] = cancer['ws'][best_ind]
fig = plot_loss(cancer['w1'], cancer['w2'], cancer['error'])
for i, g in enumerate(guesses):
fig.add_scatter3d(x=[g[0]], y=[g[1]], z=[neg_log_likelihood(g)],
mode='markers', marker=dict(size=10, color=colors[i + 2]),
name=f'w = {g}', legendgroup=str(i), row=1, col=1)
fig.add_scatter(x=[g[0]], y=[g[1]], mode='markers',
marker=dict(size=10, color=colors[i + 2]),
name=f'w = {g}', legendgroup=str(i), showlegend=False,
row=1, col=2)
label = f"Grid best [{cancer['grid_best'][0]:0.2f}, {cancer['grid_best'][1]:0.2f}]"
fig.add_scatter3d(x=[cancer['grid_best'][0]], y=[cancer['grid_best'][1]], z=[cancer['error'].min()],
mode='markers', marker=dict(size=10, color='red'),
name=label, legendgroup="Best", row=1, col=1)
fig.add_scatter(x=[cancer['grid_best'][0]], y=[cancer['grid_best'][1]],
mode='markers', marker=dict(size=10, color='red'),
name=label, legendgroup="Best", showlegend=False, row=1, col=2)
fig
plot_lr_predictions(cancer['grid_best'], cancer_df)
Squared Error for a Non-Linear Model
We now construct a less well-behaved error surface using a non-linear regression model with a squared error function. The data are generated from a noisy sine wave.
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']})
Because we generated the data, the true regression function is known: $$ h(x) = \sin(1.1 + 2.5 x), $$ and the targets are $y_n = h(x_n) + \epsilon_n$ with $\epsilon_n \sim \mathcal{N}(0, 0.5^2)$.
fig = px.scatter(sine_df, x='x', y='y')
fig.update_traces(marker_color='black')
data_trace = fig.data[0]
fig
The Sine Regression Model
We fit a model of the form $$ y(x, w) = \sin(w_0 + w_1 x), $$ with parameters $w = (w_0, w_1)$. The model is non-linear in both parameters: it is periodic in $w_0$ and oscillates faster as $w_1$ grows. As a result the error function is not convex.
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])
The following three parameter settings are evaluated on 100 test points and plotted against the data.
sine['guesses'] = np.array([[0, 2], [2, 3], [0, 3.5]])
sine['xhat'] = np.linspace(sine_df['x'].min(), sine_df['x'].max(), 100)
sine['pred_df'] = pd.DataFrame({'x': sine['xhat']})
for w in sine['guesses']:
sine['pred_df'][f'yhat(w={w})'] = sine_model(w, sine['xhat'])
sine['pred_df'].head()
| x | yhat(w=[0. 2.]) | yhat(w=[2. 3.]) | yhat(w=[0. 3.5]) | |
|---|---|---|---|---|
| 0 | 0.043371 | 0.086632 | 0.847619 | 0.151215 |
| 1 | 0.121225 | 0.240082 | 0.701797 | 0.411673 |
| 2 | 0.199080 | 0.387723 | 0.517864 | 0.641752 |
| 3 | 0.276935 | 0.525982 | 0.305809 | 0.824474 |
| 4 | 0.354790 | 0.651515 | 0.077147 | 0.946355 |
fig = go.Figure()
for i, w in enumerate(sine['pred_df'].columns[1:]):
fig.add_trace(go.Scatter(x=sine['pred_df']['x'], y=sine['pred_df'][w],
mode='lines', name=w,
line=dict(width=4, color=colors[i + 2])))
fig.add_trace(data_trace)
fig.update_layout(margin=dict(l=20, r=20, t=20, b=20),
xaxis_title='x', yaxis_title='y')
fig
None of these parameter settings fits the data well. To compare parameter settings quantitatively we need an error function.
The Mean Squared Error Surface
For a data set $\mathcal{D} = \{(x_n, y_n)\}_{n=1}^N$ we use the mean squared error $$ E(w) = \frac{1}{N} \sum_{n=1}^N \big(y_n - y(x_n, w)\big)^2 = \frac{1}{N} \sum_{n=1}^N \big(y_n - \sin(w_0 + w_1 x_n)\big)^2. $$
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)
We evaluate $E(w)$ on a grid over $(w_0, w_1)$ and plot the resulting surface. Unlike the logistic regression error, this surface has several basins separated by ridges.
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')
fig = plot_loss(sine['w0'], sine['w1'], sine['error'], names=names)
for i, w in enumerate(sine['guesses']):
fig.add_trace(go.Scatter3d(x=[w[0]], y=[w[1]], z=[sine_MSE(w)],
mode='markers', marker=dict(size=5, color=colors[i + 2]),
name=f'w=({w[0]}, {w[1]})', legendgroup=str(i)),
row=1, col=1)
fig.add_trace(go.Scatter(x=[w[0]], y=[w[1]],
mode='markers', marker=dict(size=20, color=colors[i + 2]),
name=f'w=({w[0]}, {w[1]})', legendgroup=str(i), showlegend=False),
row=1, col=2)
fig
Grid Search
As before, we take the grid point with the lowest error.
ind = np.argmin(sine['error'])
sine['grid_best'] = sine['ws'][ind, :]
sine['grid_best_error'] = sine['error'].flatten()[ind]
print(f"Best weights: {sine['grid_best']}, with error: {sine['grid_best_error']:.4f}")
Best weights: [1.29310345 2.44827586], with error: 0.2383
fig = plot_loss(sine['w0'], sine['w1'], sine['error'], names=names)
for i, w in enumerate(sine['guesses']):
fig.add_trace(go.Scatter3d(x=[w[0]], y=[w[1]], z=[sine_MSE(w)],
mode='markers', marker=dict(size=5, color=colors[i + 2]),
name=f'w=({w[0]}, {w[1]})', legendgroup=str(i)),
row=1, col=1)
fig.add_trace(go.Scatter(x=[w[0]], y=[w[1]],
mode='markers', marker=dict(size=20, color=colors[i + 2]),
name=f'w=({w[0]}, {w[1]})', legendgroup=str(i), showlegend=False),
row=1, col=2)
label = f"Grid best [{sine['grid_best'][0]:0.2f}, {sine['grid_best'][1]:0.2f}]"
fig.add_scatter3d(x=[sine['grid_best'][0]], y=[sine['grid_best'][1]], z=[sine['grid_best_error']],
mode='markers', marker=dict(size=10, color='red'),
name=label, legendgroup="Best", row=1, col=1)
fig.add_scatter(x=[sine['grid_best'][0]], y=[sine['grid_best'][1]],
mode='markers', marker=dict(size=20, color='red'),
name=label, legendgroup="Best", showlegend=False, row=1, col=2)
fig
The model at the best grid point is a much better fit to the data than the three initial guesses.
fig = go.Figure()
for i, w in enumerate(sine['pred_df'].columns[1:]):
fig.add_trace(go.Scatter(x=sine['pred_df']['x'], y=sine['pred_df'][w], mode='lines', name=w,
line=dict(width=4, color=colors[i + 2])))
fig.add_trace(data_trace)
fig.add_scatter(x=sine['xhat'], y=sine_model(sine['grid_best'], sine['xhat']),
mode='lines',
name=f"Grid best w=({sine['grid_best'][0]:.2f}, {sine['grid_best'][1]:.2f})",
line=dict(width=8, color='red'))
fig.update_layout(margin=dict(l=20, r=20, t=20, b=20),
xaxis_title='x', yaxis_title='y')
fig
(Lecture Slido) Match the Model with the Error Function
The following three cells produce the plots used in the Slido question that asks students to match each model with its error function. Each model is fit to one feature of the breast cancer data.
- Logistic regression with a single weight, $p(t=1 \mid x) = \sigma(w x)$, using the standardized mean radius as $x$. The error is a function of one parameter.
- Least squares regression through the origin, $y = w x$, predicting the standardized mean texture from the standardized mean radius. The error is a quadratic function of one parameter.
- Logistic regression with a bias, $p(t=1 \mid x) = \sigma(w_1 + w_2 x)$, with the feature shifted by $+3$. Because the shifted feature is not centered, the bias and the slope trade off against each other: many pairs $(w_1, w_2)$ with a similar decision boundary $x = -w_1 / w_2$ fit the data about equally well, so the error surface has a long, narrow valley.
A small amount of vertical jitter is added to the binary targets so that overlapping points are visible.
rng = np.random.default_rng(0)
x1 = cancer['x'][:, 0]
fig1 = px.scatter(x=x1, y=cancer['t'] + 0.02 * rng.normal(size=cancer['t'].shape),
color=cancer_df['diagnosis'], color_discrete_map=DIAGNOSIS_COLORS,
title="Cancer Data (Jittered)",
labels={'x': 'Mean Radius (standardized)', 'y': 't', 'color': 'diagnosis'})
ws = np.linspace(-20, 3, 100)
nll = np.array([neg_log_likelihood(np.array([w]), x1[:, None]) for w in ws])
fig2 = px.line(x=ws, y=nll, labels={'x': 'w', 'y': 'Negative Log-Likelihood'})
ind = np.argmin(nll)
best_w = ws[ind]
fig2.add_scatter(x=[best_w], y=[nll[ind]], mode='markers',
marker=dict(color='red', size=10), name="Best w")
xtest = np.linspace(x1.min(), x1.max(), 100)
fig1.add_scatter(x=xtest, y=sigmoid(best_w * xtest),
mode='lines', line=dict(color='black', width=4),
name=f'Logistic Model (w={best_w:.2f})')
fig = make_plot_grid([fig1, fig2], 1, 2)
fig.add_vline(x=best_w, line=dict(color='red', dash='dash'),
annotation_text=f"Best w={best_w:.2f}",
annotation_position="top right", row=1, col=2)
fig.update_layout(height=600)
x1, x2 = cancer['x'][:, 0], cancer['x'][:, 1]
fig1 = px.scatter(x=x1, y=x2, title="Cancer Data",
labels={'x': 'Mean Radius (standardized)', 'y': 'Mean Texture (standardized)'})
ws = np.linspace(-5, 5, 100)
sqloss = np.array([np.mean((x2 - w * x1) ** 2) for w in ws])
fig2 = px.line(x=ws, y=sqloss, labels={'x': 'w', 'y': 'Squared Loss'})
ind = np.argmin(sqloss)
best_w = ws[ind]
xtest = np.linspace(x1.min(), x1.max(), 100)
fig1.add_scatter(x=xtest, y=best_w * xtest,
mode='lines', line=dict(color='black', width=4),
name=f'Linear Model (w={best_w:.2f})')
fig2.add_scatter(x=[best_w], y=[sqloss[ind]], mode='markers',
marker=dict(color='red', size=10), name="Best w")
fig = make_plot_grid([fig1, fig2], 1, 2)
fig.add_vline(x=best_w, line=dict(color='red', dash='dash'),
annotation_text=f"Best w={best_w:.2f}",
annotation_position="top right", row=1, col=2)
fig.update_layout(height=600)
offset = 3
x1 = cancer['x'][:, 0] + offset
fig1 = px.scatter(x=x1, y=cancer['t'] + 0.02 * rng.normal(size=cancer['t'].shape),
color=cancer_df['diagnosis'], color_discrete_map=DIAGNOSIS_COLORS,
title="Cancer Data (Jittered, Shifted)",
labels={'x': 'Mean Radius (standardized + 3)', 'y': 't', 'color': 'diagnosis'})
w1, w2 = np.meshgrid(np.linspace(5, 15, 40), np.linspace(-8, 0, 40))
ws = np.stack([w1.flatten(), w2.flatten()]).T
# Features [1, x] so that w[0] is the bias and w[1] is the slope.
x_aug = np.stack([np.ones_like(x1), x1], axis=1)
nll = np.array([neg_log_likelihood(w, x_aug) for w in ws]).reshape(w1.shape)
fig2 = plot_loss(w1, w2, nll, names=('w1 (bias)', 'w2 (slope)'))
ind = np.argmin(nll)
best_w = ws[ind, :]
fig2.add_scatter3d(x=[best_w[0]], y=[best_w[1]], z=[nll.flatten()[ind]],
mode='markers', marker=dict(size=10, color='red'),
name=f'Best w=({best_w[0]:.2f}, {best_w[1]:.2f})', legendgroup='Best',
row=1, col=1)
fig2.add_scatter(x=[best_w[0]], y=[best_w[1]],
mode='markers', marker=dict(size=20, color='red'),
name=f'Best w=({best_w[0]:.2f}, {best_w[1]:.2f})', legendgroup='Best',
showlegend=False, row=1, col=2)
fig2.show()
xtest = np.linspace(x1.min(), x1.max(), 100)
fig1.add_scatter(x=xtest, y=sigmoid(best_w[1] * xtest + best_w[0]),
mode='lines', line=dict(color='black', width=4),
name=f'Logistic Model (w=({best_w[0]:.2f}, {best_w[1]:.2f}))')
fig1.update_layout(height=600)
(Lecture Demo 2) The Gradient of the Error Function
The gradient $\nabla E(w)$ is the vector of partial derivatives $\partial E / \partial w_j$. It points in the direction of the greatest rate of increase of $E$, and for a small step $\delta w$ the change in the error is $\delta E \simeq \delta w^\top \nabla E(w)$ (Bishop, Eq. 7.1). In the plots below the gradient is drawn as a vector field in weight space: cones under the surface on the left, and arrows on the contour plot on the right. The gradient lives in the $D$-dimensional weight space, not on the error surface itself.
Lecture: run from the slide "Demo: The Gradient of the Error Function," after the Slido question on the dimension of the gradient.
Gradient of the Cross-Entropy Error
For the logistic regression model without a bias, $y_n = \sigma(w^\top x_n)$, and the average negative log-likelihood is $$ E(w) = -\frac{1}{N} \sum_{n=1}^N \Big( t_n \ln \sigma(w^\top x_n) + (1 - t_n) \ln\big(1 - \sigma(w^\top x_n)\big) \Big). $$ The derivative of the logistic sigmoid satisfies $$ \frac{d\sigma}{dz} = \sigma(z)\big(1 - \sigma(z)\big). $$
Applying the chain rule, the $j^\text{th}$ component of the gradient is \begin{align*} \frac{\partial E}{\partial w_j} &= -\frac{1}{N} \sum_{n=1}^N \left( \frac{t_n}{\sigma(w^\top x_n)} \frac{\partial \sigma(w^\top x_n)}{\partial w_j} - \frac{1 - t_n}{1 - \sigma(w^\top x_n)} \frac{\partial \sigma(w^\top x_n)}{\partial w_j} \right) \\ &= -\frac{1}{N} \sum_{n=1}^N \Big( t_n \big(1 - \sigma(w^\top x_n)\big) - (1 - t_n)\, \sigma(w^\top x_n) \Big) x_{nj} \\ &= \frac{1}{N} \sum_{n=1}^N \big(\sigma(w^\top x_n) - t_n\big) x_{nj}. \end{align*} The gradient is the average of the prediction errors $y_n - t_n$ weighted by the inputs, the same form as for linear regression with a squared error (Bishop, Section 5.4.3).
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)
grad_NLL(np.array([-1., 2.]))
array([0.30254033, 0.44824536])
We evaluate the gradient at every grid point and plot the resulting vector field. The arrows are long where the surface is steep and short near the minimum.
(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)
fig = plot_gradient(cancer['w1'], cancer['w2'], cancer['error'],
cancer['dw1'], cancer['dw2'], scale=2)
fig
Gradient of the Squared Error for the Sine Model
For the sine model $y(x, w) = \sin(w_0 + w_1 x)$ with the mean squared error $$ E(w) = \frac{1}{N} \sum_{n=1}^N \big(y_n - \sin(w_0 + w_1 x_n)\big)^2, $$ the gradient is $$ \nabla E(w) = \begin{bmatrix} \partial E / \partial w_0 \\ \partial E / \partial w_1 \end{bmatrix}. $$ Writing $u_n = w_0 + w_1 x_n$ and applying the chain rule gives \begin{align*} \frac{\partial E}{\partial w_0} &= -\frac{2}{N} \sum_{n=1}^N \big(y_n - \sin u_n\big) \cos u_n, \\ \frac{\partial E}{\partial w_1} &= -\frac{2}{N} \sum_{n=1}^N \big(y_n - \sin u_n\big) \cos u_n \; x_n. \end{align*}
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])
grad_sine_MSE(np.array([0., 2.]))
array([-0.18866281, 0.66526642])
(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)
fig = plot_gradient(sine['w0'], sine['w1'], sine['error'],
sine['dw0'], sine['dw1'], scale=0.1, names=names)
fig
(Lecture Demo 3) Batch Gradient Descent
Gradient descent starts from an initial weight vector $w^{(0)}$ and repeatedly takes a small step in the direction of the negative gradient (Bishop, Eq. 7.16): $$ w^{(\tau)} = w^{(\tau-1)} - \eta \nabla E\big(w^{(\tau-1)}\big), $$ where $\eta > 0$ is the learning rate. Because $E$ is defined over the whole training set, each step requires a pass through all $N$ data points, and the method is therefore called batch gradient descent. The implementation below uses a fixed learning rate and stops when the update becomes smaller than a threshold $\epsilon$ or after a maximum number of iterations.
Lecture: run from the slide "Demo: Batch Gradient Descent." The loss curves below go with the slide "Loss (Error) Curves."
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)
Gradient Descent for Logistic Regression
Starting from $w^{(0)} = (0, 0)$ with $\eta = 1$, the iterates move downhill, perpendicular to the contours, and then along the floor of the valley towards the minimum. Progress along the valley is slow because the curvature in that direction is small, so after 100 iterations the final point is still short of the minimum $w^\star \approx (-3.96, -0.89)$ found below.
w0 = np.array([0., 0.])
path = gradient_descent(w0, grad_NLL, learning_rate=1, nepochs=100)
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)
Gradient Descent for the Sine Model
For the sine model we use a smaller learning rate, $\eta = 0.01$. The error surface is much more strongly curved in the $w_1$ direction than in the $w_0$ direction, because $\partial u_n / \partial w_1 = x_n$ and the inputs range up to $2.5\pi$.
w0 = np.array([2., 2.])
path = gradient_descent(w0, grad_sine_MSE, learning_rate=.01, nepochs=20)
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)
Loss Curves
In problems with more than two or three parameters the error surface cannot be visualized. Instead, we monitor the loss curve, which plots the error against the iteration number. A curve that decreases steadily indicates a suitable learning rate. A curve that oscillates or increases indicates that the learning rate is too large, and a curve that decreases very slowly suggests that it is too small.
The first example uses the sine model with $\eta = 0.05$.
w0 = np.array([2.5, 1.9])
path = gradient_descent(w0, grad_sine_MSE, learning_rate=.05, 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)
fig.show()
The second example returns to logistic regression, starting far from the minimum at $w^{(0)} = (-10, -5)$ with the much larger learning rate $\eta = 10$. The iteration is stable because the curvature of this error surface is small: near the minimum the largest eigenvalue of the Hessian is approximately $0.08$, so any $\eta < 2/\lambda_{\max} \approx 25$ is locally stable (this condition is derived in Lecture 11). The appropriate scale of $\eta$ is therefore problem dependent, and $0 < \eta < 1$ is only a common starting point.
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)
(Lecture Demo 4) Gradient Descent on a Parabola
Consider $f(x) = x^2 + 1$, with $f'(x) = 2x$. A gradient descent step is $$ x^{(\tau)} = x^{(\tau-1)} - 2\eta\, x^{(\tau-1)} = (1 - 2\eta)\, x^{(\tau-1)}, $$ so $x^{(\tau)} = (1 - 2\eta)^\tau x^{(0)}$. The iterates converge to the minimum at $x = 0$ if and only if $|1 - 2\eta| < 1$, that is $0 < \eta < 1$. With $\eta = 1.1$ the factor is $-1.2$, and the iterates alternate in sign and grow in magnitude. Gradient descent is therefore not guaranteed to converge, even on a simple quadratic.
The axes below are fixed, so the plot can be compared across learning rates. For $\eta$ slightly below 1 the iterates alternate in sign and shrink, at $\eta = 1$ they alternate between $\pm 5$ (the starting value $x^{(0)} = 5$), and for $\eta$ slightly above 1 they alternate and grow.
Lecture: run from the slide "Demo: Gradient Descent on a Parabola," around the Slido question "Assuming a simple quadratic equation (e.g., x^2) is gradient descent guaranteed to converge?" The plot with the default settings is the figure on the slide "What are the conditions for convergence?", which leads into the quadratic approximation below: near a minimum every smooth error function looks like this parabola, with the curvature given by the Hessian.
# Try values slightly below and above 1, e.g., 0.9, 0.95, 1.0, 1.05, 1.1.
learning_rate = 1.05
x0, nepochs = 5., 10
def grad_f(x):
return 2 * x
path = gradient_descent(x0, grad_f, learning_rate=learning_rate, nepochs=nepochs)
x = np.linspace(-13, 13, 200)
fig = px.line(x=x, y=x**2 + 1, labels={'x': 'x', 'y': 'f(x)'})
fig.update_traces(line_width=5)
fig.add_scatter(x=path, y=path**2 + 1, mode='markers+lines',
marker=dict(size=10, color='red'),
line=dict(color='black', width=2, dash="dash"),
name='Gradient Descent Path')
# Fixed axes so that the plot does not rescale when the learning rate changes.
# Points that leave this range (large learning rates) are not shown.
fig.update_layout(title=f"η = {learning_rate}, x(τ) = (1 − 2η)^τ x(0) with 1 − 2η = {1 - 2 * learning_rate:.2f}",
xaxis_range=[-13, 13], yaxis_range=[-5, 175], height=500)
fig
(Lecture Demo 5) Second-Order Structure of the Error Function
The Hessian $H = \nabla\nabla E(w)$ is the $D \times D$ matrix of second derivatives, $$ H_{jk} = \frac{\partial^2 E}{\partial w_j \, \partial w_k}. $$ It describes the local curvature of the error surface. Because mixed partial derivatives commute for smooth functions, $H$ is symmetric.
Lecture: the exercise below has the plots for the slide "Hessian Exercise." Demo 5 continues in "The Quadratic Approximation."
Exercise: A Simple Quadratic
For $E(w) = (w_0 - 1)^2 + (w_1 - 2)^2 + 1$ the gradient is $\nabla E(w) = \big(2(w_0 - 1),\, 2(w_1 - 2)\big)^\top$ and the Hessian is the constant matrix $H = 2I$. Both eigenvalues equal 2, so the stationary point $w^\star = (1, 2)$ is a minimum and the contours are circles.
import sympy as sp
w0_s, w1_s = sp.symbols('w0 w1')
E_ex = (w0_s - 1)**2 + (w1_s - 2)**2 + 1
display(sp.Matrix([sp.diff(E_ex, v) for v in (w0_s, w1_s)]))
display(sp.hessian(E_ex, (w0_s, w1_s)))
# Error surface and gradient field of the exercise function (plots for the "Hessian Exercise" slide).
ex_w0, ex_w1 = np.meshgrid(np.linspace(-1, 3, 13), np.linspace(0, 4, 13))
ex_E = (ex_w0 - 1)**2 + (ex_w1 - 2)**2 + 1
plot_gradient(ex_w0, ex_w1, ex_E, 2 * (ex_w0 - 1), 2 * (ex_w1 - 2),
scale=0.08, names=('w0', 'w1'))
The Hessian of the Logistic Regression Error
Starting from the gradient derived above, $$ \frac{\partial E}{\partial w_j} = \frac{1}{N} \sum_{n=1}^N \big(\sigma(w^\top x_n) - t_n\big) x_{nj}, $$ we differentiate with respect to $w_k$ and use $d\sigma/dz = \sigma(z)(1 - \sigma(z))$: \begin{align*} \frac{\partial^2 E}{\partial w_k \, \partial w_j} &= \frac{1}{N} \sum_{n=1}^N x_{nj} \frac{\partial}{\partial w_k} \sigma(w^\top x_n) \\ &= \frac{1}{N} \sum_{n=1}^N \sigma(w^\top x_n)\big(1 - \sigma(w^\top x_n)\big) x_{nj} x_{nk}. \end{align*} In matrix form, $H = \frac{1}{N} X^\top R X$ where $R$ is diagonal with $R_{nn} = y_n (1 - y_n)$. Since $R_{nn} > 0$, the Hessian is positive semi-definite for every $w$, which shows that the logistic regression error function is convex.
def hessian_NLL(w):
"""Hessian of the average negative log-likelihood of the logistic model."""
x = cancer['x']
p = logistic_model(w, x)
return (x * (p * (1 - p))[:, None]).T @ x / len(p)
The Hessian of the Sine Regression Error
With $u_n = w_0 + w_1 x_n$, the gradient of the mean squared error is $$ \frac{\partial E}{\partial w_0} = -\frac{2}{N} \sum_{n=1}^N (y_n - \sin u_n) \cos u_n, \qquad \frac{\partial E}{\partial w_1} = -\frac{2}{N} \sum_{n=1}^N (y_n - \sin u_n) \cos u_n \, x_n. $$ Differentiating again with the product rule, and defining $a_n = \cos^2 u_n + (y_n - \sin u_n) \sin u_n$, gives $$ \frac{\partial^2 E}{\partial w_0^2} = \frac{2}{N} \sum_{n=1}^N a_n, \qquad \frac{\partial^2 E}{\partial w_0 \, \partial w_1} = \frac{2}{N} \sum_{n=1}^N a_n x_n, \qquad \frac{\partial^2 E}{\partial w_1^2} = \frac{2}{N} \sum_{n=1}^N a_n x_n^2. $$ The term $(y_n - \sin u_n)\sin u_n$ can be negative, so $a_n$ can be negative and the Hessian need not be positive semi-definite. The check below finds points on the grid where the Hessian has a negative eigenvalue, so the error function is not convex.
def hessian_sine_MSE(w):
"""Hessian of the mean squared error of the sine model."""
x = sine['x']
y = sine['y']
u = w[0] + w[1] * x
a = np.cos(u) ** 2 + (y - np.sin(u)) * np.sin(u)
h00 = 2 * np.mean(a)
h01 = 2 * np.mean(a * x)
h11 = 2 * np.mean(a * x ** 2)
return np.array([[h00, h01], [h01, h11]])
# The smallest Hessian eigenvalue over the grid is negative, so E is not convex.
min_eig = min(np.linalg.eigvalsh(hessian_sine_MSE(w))[0] for w in sine['ws'])
print(f"Smallest Hessian eigenvalue on the grid: {min_eig:.2f}")
Smallest Hessian eigenvalue on the grid: -21.80
Checking the Derivation with Symbolic Differentiation
The sympy library can compute the same derivatives symbolically. We define the error as a symbolic sum over the data, compute its Hessian, and convert the result into a numerical function. On Colab, sympy may need to be upgraded first with !pip install --upgrade sympy.
import sympy as sp
w0_s, w1_s, x_s, y_s = sp.symbols('w0 w1 x y')
# Squared error for a single data point and its gradient and Hessian
E_n = (y_s - sp.sin(w0_s + w1_s * x_s))**2
display(sp.Matrix([sp.diff(E_n, v) for v in (w0_s, w1_s)]))
display(sp.hessian(E_n, (w0_s, w1_s)))
n_s = sp.symbols('n', integer=True, positive=True)
i_s = sp.Idx('i', n_s)
x_i, y_i = sp.IndexedBase('x'), sp.IndexedBase('y')
# Mean squared error over the data set as a symbolic sum
E_sym = sp.Sum((y_i[i_s] - sp.sin(w0_s + w1_s * x_i[i_s]))**2, (i_s, 0, n_s - 1)) / n_s
H_sym = sp.hessian(E_sym, (w0_s, w1_s))
H_fun = sp.lambdify((w0_s, w1_s, x_i, y_i, n_s), H_sym, modules='numpy')
def hessian_sine_MSE_sympy(w):
return np.array(H_fun(w[0], w[1], sine['x'], sine['y'], len(sine['x'])), dtype=float)
w_test = np.array([1.1, 2.5])
print(hessian_sine_MSE(w_test))
print(hessian_sine_MSE_sympy(w_test))
assert np.allclose(hessian_sine_MSE(w_test), hessian_sine_MSE_sympy(w_test))
[[ 0.9188191 3.45856751] [ 3.45856751 17.87279593]] [[ 0.9188191 3.45856751] [ 3.45856751 17.87279593]]
(Lecture Demo 5) The Quadratic Approximation
The second-order Taylor expansion of $E(w)$ around a point $\widehat{w}$ is (Bishop, Eq. 7.3) $$ E(w) \simeq E(\widehat{w}) + (w - \widehat{w})^\top b + \frac{1}{2} (w - \widehat{w})^\top H (w - \widehat{w}), $$ where $b = \nabla E(\widehat{w})$ and $H = \nabla\nabla E(\widehat{w})$. At a stationary point $w^\star$ the linear term vanishes. Let $H u_i = \lambda_i u_i$ with orthonormal eigenvectors $u_i$, and define $\alpha_i = u_i^\top (w - w^\star)$. Then (Bishop, Eq. 7.11) $$ E(w) \simeq E(w^\star) + \frac{1}{2} \sum_i \lambda_i \alpha_i^2. $$ The contours of constant error are ellipses whose axes are aligned with the eigenvectors $u_i$, with lengths proportional to $\lambda_i^{-1/2}$. If all eigenvalues are positive, $w^\star$ is a local minimum.
Lecture: run from the slide "Demo: Taylor Expansions and the Hessian." The eigenvector plot for the sine model is the figure on the slide "Eigenvectors of the Hessian."
def taylor_loss(w, w_hat, E_hat, g_hat, H_hat):
"""Second-order Taylor approximation of the error around w_hat."""
delta = w - w_hat
return E_hat + delta @ g_hat + 0.5 * delta @ H_hat @ delta
def plot_eigenvectors(w1, w2, taylor, w_star, lam, U, names=('w1', 'w2'), scaling=0.5):
"""Contours of a quadratic approximation with the Hessian eigenvectors drawn
from w_star with lengths proportional to 1/sqrt(lambda_i)."""
fig = go.Figure()
fig.add_contour(z=taylor.flatten(), x=w1.flatten(), y=w2.flatten(),
colorscale='viridis_r', opacity=0.5,
contours=dict(start=taylor.min(), end=taylor.max(),
size=(taylor.max() - taylor.min()) / 50),
colorbar=dict(x=1.05, y=0.3, len=0.75))
lengths = scaling / np.sqrt(np.abs(lam))
cx, cy = w_star
for i, (lam_i, u_i, s_i) in enumerate(zip(lam, U.T, lengths), start=1):
color = colors[i]
fig.add_scatter(x=[cx, cx + s_i * u_i[0]], y=[cy, cy + s_i * u_i[1]],
mode='lines+markers', line=dict(width=2, color=color),
name=f'u{i} (λ={lam_i:.3g})', legendgroup=f'u{i}')
fig.add_scatter(x=[cx, cx - s_i * u_i[0]], y=[cy, cy - s_i * u_i[1]],
mode='lines', line=dict(width=2, dash='dot', color=color),
legendgroup=f'u{i}', showlegend=False)
fig.update_layout(title="Quadratic Approximation of the Error",
xaxis_title=names[0], yaxis_title=names[1],
height=600, width=1200)
# Equal axis scales, so that orthogonal eigenvectors look orthogonal.
fig.update_xaxes(range=[w1.min(), w1.max()], constrain='domain')
fig.update_yaxes(range=[w2.min(), w2.max()], constrain='domain',
scaleanchor='x', scaleratio=1)
return fig
Logistic Regression
We expand around the minimum $w^\star$. The best grid point is only an approximation to the minimum, so we first refine it with gradient descent until the gradient is close to zero.
cancer['w_star'] = gradient_descent(cancer['grid_best'], grad_NLL,
learning_rate=10, nepochs=1000, epsilon=1e-10)[-1]
cancer['E_star'] = neg_log_likelihood(cancer['w_star'])
cancer['g_star'] = grad_NLL(cancer['w_star'])
cancer['H_star'] = hessian_NLL(cancer['w_star'])
cancer['lam'], cancer['U'] = np.linalg.eigh(cancer['H_star'])
print("w_star:", cancer['w_star'])
print("E(w_star):", cancer['E_star'])
print("gradient at w_star:", cancer['g_star'])
print("Hessian at w_star:\n", cancer['H_star'])
print("Eigenvalues:", cancer['lam'])
print("Eigenvectors (columns):\n", cancer['U'])
print("Condition number:", cancer['lam'][-1] / cancer['lam'][0])
w_star: [-3.96428592 -0.89302098] E(w_star): 0.2757570726465129 gradient at w_star: [7.69098239e-12 9.69307040e-13] Hessian at w_star: [[ 0.01430723 -0.00838548] [-0.00838548 0.07978486]] Eigenvalues: [0.0132504 0.0808417] Eigenvectors (columns): [[-0.99215133 -0.12504291] [-0.12504291 0.99215133]] Condition number: 6.101077926924686
Both eigenvalues are positive, so $w^\star$ is a minimum, as expected for a convex error function. The surface below overlays the quadratic approximation (orange) on the true error surface. The approximation is accurate near $w^\star$ and degrades farther away.
fig = plot_gradient(cancer['w1'], cancer['w2'], cancer['error'], cancer['dw1'], cancer['dw2'], scale=2)
cancer['taylor_loss'] = np.array([
taylor_loss(w, cancer['w_star'], cancer['E_star'], cancer['g_star'], cancer['H_star'])
for w in cancer['ws']]).reshape(cancer['w1'].shape)
tl = cancer['taylor_loss']
fig.add_surface(z=tl, x=cancer['w1'], y=cancer['w2'],
colorscale='plasma_r', opacity=0.5, showscale=False,
contours=dict(z=dict(show=True, highlightcolor="white",
start=tl.min(), end=tl.max(),
size=(tl.max() - tl.min()) / 50)),
row=1, col=1)
fig.update_layout(scene=dict(zaxis=dict(range=[0, 3])))
fig
The eigenvector with the smaller eigenvalue points along the long axis of the elliptical contours, the direction of lowest curvature.
plot_eigenvectors(cancer['w1'], cancer['w2'], cancer['taylor_loss'],
cancer['w_star'], cancer['lam'], cancer['U'])
Sine Regression
We repeat the analysis for the sine model, refining the best grid point to the nearby local minimum.
sine['w_star'] = gradient_descent(sine['grid_best'], grad_sine_MSE,
learning_rate=0.05, nepochs=5000, epsilon=1e-10)[-1]
sine['E_star'] = sine_MSE(sine['w_star'])
sine['g_star'] = grad_sine_MSE(sine['w_star'])
sine['H_star'] = hessian_sine_MSE(sine['w_star'])
sine['lam'], sine['U'] = np.linalg.eigh(sine['H_star'])
print("w_star:", sine['w_star'])
print("E(w_star):", sine['E_star'])
print("gradient at w_star:", sine['g_star'])
print("Hessian at w_star:\n", sine['H_star'])
print("Eigenvalues:", sine['lam'])
print("Eigenvectors (columns):\n", sine['U'])
print("Condition number:", sine['lam'][-1] / sine['lam'][0])
w_star: [1.11301782 2.48818099] E(w_star): 0.23391180154084013 gradient at w_star: [ 1.92756179e-09 -3.74951004e-10] Hessian at w_star: [[ 0.93080729 3.53734477] [ 3.53734477 18.42755312]] Eigenvalues: [ 0.24271707 19.11564335] Eigenvectors (columns): [[-0.98160118 0.19094271] [ 0.19094271 0.98160118]] Condition number: 78.75689883881263
Both eigenvalues are positive, so $w^\star$ is a local minimum. The error function is not convex, however: the Hessian is indefinite at other points of weight space (see the check above), and the quadratic approximation describes the surface only within the basin around $w^\star$. The two eigenvalues differ by a factor of about 80, so the contours near the minimum are long, narrow ellipses.
fig = plot_gradient(sine['w0'], sine['w1'], sine['error'], sine['dw0'], sine['dw1'],
scale=0.1, names=names)
sine['taylor_loss'] = np.array([
taylor_loss(w, sine['w_star'], sine['E_star'], sine['g_star'], sine['H_star'])
for w in sine['ws']]).reshape(sine['w1'].shape)
tl = sine['taylor_loss']
fig.add_surface(z=tl, x=sine['w0'], y=sine['w1'],
colorscale='plasma_r', opacity=0.5, showscale=False,
contours=dict(z=dict(show=True, highlightcolor="white",
start=tl.min(), end=tl.max(),
size=(tl.max() - tl.min()) / 50)),
row=1, col=1)
fig.update_layout(scene=dict(zaxis=dict(range=[0, 3])))
fig
plot_eigenvectors(sine['w0'], sine['w1'], sine['taylor_loss'],
sine['w_star'], sine['lam'], sine['U'], names=names)
Classifying Stationary Points
At a stationary point the quadratic approximation is $E(w^\star) + \frac{1}{2}\sum_i \lambda_i \alpha_i^2$, so the signs of the eigenvalues determine the type of the point: all positive gives a minimum, all negative a maximum, and mixed signs a saddle point. Moving along $u_i$ increases the error if $\lambda_i > 0$ and decreases it if $\lambda_i < 0$. The figure for the slides (slides/slides.ipynb, "Hessian directions at stationary points") draws the eigenvectors for each case. The helper below classifies a stationary point from its Hessian.
def classify_stationary_point(H, tol=1e-10):
lam = np.linalg.eigvalsh(H)
if np.all(lam > tol):
return 'minimum', lam
if np.all(lam < -tol):
return 'maximum', lam
if np.any(lam > tol) and np.any(lam < -tol):
return 'saddle point', lam
return 'degenerate (some eigenvalues are zero)', lam
for name, H_ex in [('minimum', np.diag([0.25, 1.0])), ('maximum', np.diag([-0.25, -1.0])),
('saddle', np.diag([0.25, -1.0])), ('sine model at w*', sine['H_star'])]:
print(f"{name:>18}: {classify_stationary_point(H_ex)}")
minimum: ('minimum', array([0.25, 1. ]))
maximum: ('maximum', array([-1. , -0.25]))
saddle: ('saddle point', array([-1. , 0.25]))
sine model at w*: ('minimum', array([ 0.24271707, 19.11564335]))
Next Lecture
The Lecture 11 notebook continues from this point. It uses the quadratic approximation developed above to analyze when gradient descent converges and how fast, and then introduces momentum, learning rate schedules, Adam, and stochastic gradient descent.