Week 4 seminar

Author

Joshua Loftus

Note

This is the student version of the notebook. The code you write goes where a line says # ... or has a ____. A code chunk with a blank starts with #| eval: false, so the file renders as it stands; delete that line once you have filled the blank. Working through the chunks one at a time in RStudio (Ctrl/Cmd+Enter) is unaffected. Answers are revealed in the seminar.

library(tidyverse)       # should've already install.packages("tidyverse")
library(gapminder)       # might need to install.packages("gapminder")
library(palmerpenguins)  # might need to install.packages("palmerpenguins")
theme_set(theme_minimal())

Before the seminar

The lecture proved one result, the descent lemma, which you prove again with the slides closed (exercise 1). It stated the ERM form of the bias–variance decomposition and left its two-line proof to you (exercise 2), set the Newton-Raphson exercise on the logistic log-likelihood (exercise 3), and stated no free lunch (exercise 4, a proof on a finite domain). Do all four on paper before the seminar; the seminar’s problems are variations on them. You may use, without proof: the fundamental theorem of calculus along a line, \(f(y) - f(x) = \int_0^1 \nabla f(x + t(y - x))^T (y - x)\, dt\) for a differentiable \(f\); the Cauchy-Schwarz inequality \(|a^T b| \le \|a\|\,\|b\|\); and the tower property, \(\mathbb E[W] = \mathbb E\big[\mathbb E[W \mid X]\big]\). Throughout, assume that every expectation that appears exists.

WarningHints

Each exercise has hints, collapsed below it. Do not open a hint unless you have been stuck for at least a few minutes, and open them in order. Exercise 1’s proof is on the lecture slides. Solutions to exercises 2 to 4 are shown in the seminar if you ask.

Exercise 1 (descent lemma). Let \(f: \mathbb R^p \to \mathbb R\) be differentiable with an \(L\)-Lipschitz gradient: \(\|\nabla f(x) - \nabla f(y)\| \le L\,\|x - y\|\) for all \(x, y\), for some \(L > 0\). Fix \(x\), a step size \(\gamma > 0\), and write \(y = x - \gamma\nabla f(x)\) and \(v = y - x\).

  1. Write \(f(y) - f(x)\) as an integral along the segment from \(x\) to \(y\), and split it into \(\nabla f(x)^T v\) plus a remainder.

  2. Bound the remainder by \(\tfrac{L}{2}\|v\|^2\).

  3. Substitute \(v = -\gamma\nabla f(x)\) and show that \(f(y) \le f(x) - \gamma\big(1 - \tfrac12 L\gamma\big)\|\nabla f(x)\|^2\). For which \(\gamma\) is the right-hand side below \(f(x)\), and what does \(\gamma = 1/L\) give?

  4. The lemma says each step decreases \(f\). Give a function with an \(L\)-Lipschitz gradient for which the iterates \(x_{t+1} = x_t - \nabla f(x_t)/L\) never converge, and say in one sentence what the lemma does and does not guarantee.

  1. Apply the fundamental theorem of calculus along the line to \(f\), with \(x_t = x + tv\). To split, add and subtract \(\nabla f(x)^T v\) inside the integral.
  1. Cauchy-Schwarz on \([\nabla f(x_t) - \nabla f(x)]^T v\), then the Lipschitz bound, then \(\|x_t - x\| = t\|v\|\). What is \(\int_0^1 t\, dt\)?
  1. Think of a function with no minimum at all. Its gradient is as Lipschitz as you like.

Exercise 2 (bias–variance decomposition, ERM form). Risk \(R(f) = \mathbb E\,\ell(Y, f(X))\). Let \(f^\star\) minimize \(R\) over all predictors, let \(\mathcal F\) be a class of predictors with a best member \(f_{\mathcal F} = \arg\min_{f \in \mathcal F} R(f)\), and let \(\hat f \in \mathcal F\) be the predictor an algorithm returns from the sample (the empirical risk minimizer, or anything else in \(\mathcal F\)).

  1. Show that \(R(\hat f) - R(f_{\mathcal F}) \ge 0\) and \(R(f_{\mathcal F}) - R(f^\star) \ge 0\).

  2. Show that \(R(\hat f) = R(f^\star) + \big[R(\hat f) - R(f_{\mathcal F})\big] + \big[R(f_{\mathcal F}) - R(f^\star)\big]\), and say which term is random and which is fixed once \(\mathcal F\) is chosen.

  3. For squared loss, use week 2’s optimal model theorem to write the approximation term as a squared distance.

  4. The smallest class there is, \(\mathcal F = \{\text{constant functions}\}\), squared loss, \(\hat f = \bar y\) from an i.i.d. sample of size \(n\). Find \(f_{\mathcal F}\), show that the approximation error is \(\mathrm{Var}(f^\star(X))\), and show that the expected estimation error is \(\mathrm{Var}(Y)/n\).

  1. What does “\(f_{\mathcal F}\) minimizes \(R\) over \(\mathcal F\)” say about \(R(\hat f)\), given that \(\hat f \in \mathcal F\)?
  1. The optimal model theorem: \(R(f) - R(f^\star) = \mathbb E[(f^\star(X) - f(X))^2]\) for every \(f\). Apply it to \(f = f_{\mathcal F}\).
  1. For a constant \(c\), \(R(c) = \mathbb E(Y - c)^2 = \mathrm{Var}(Y) + (c - \mathbb E Y)^2\). Then the optimal model theorem once more, and \(\mathrm{Var}(Y) = \mathbb E[\mathrm{Var}(Y \mid X)] + \mathrm{Var}(\mathbb E[Y \mid X])\).

Exercise 3 (score equations, Hessian, Newton-Raphson). Logistic regression with design matrix \(\mathbf X\) (\(n\) rows \(\mathbf x_i^T\)), labels \(y_i \in \{0, 1\}\), and \(p_i(\beta) = 1/(1 + e^{-\mathbf x_i^T\beta})\).

  1. Show that the log-likelihood can be written \(\ell(\beta) = \sum_i \big[y_i\,\mathbf x_i^T\beta - \log(1 + e^{\mathbf x_i^T\beta})\big]\), and that its gradient is \(\nabla\ell(\beta) = \sum_i \mathbf x_i\,(y_i - p_i(\beta)) = \mathbf X^T(\mathbf y - \mathbf p(\beta))\).

  2. Differentiate once more to get the Hessian \(\nabla^2\ell(\beta) = -\sum_i p_i(1 - p_i)\,\mathbf x_i\mathbf x_i^T = -\mathbf X^T\mathbf W\mathbf X\) with \(\mathbf W = \mathrm{diag}(p_i(1 - p_i))\).

  3. Show that \(v^T\nabla^2\ell(\beta)\,v \le 0\) for every \(v\), so \(\ell\) is concave, and say when the inequality is strict for every \(v \ne 0\).

  4. Newton-Raphson finds a root of \(\nabla\ell\) by solving the linear approximation \(\nabla\ell(\beta_t) + \nabla^2\ell(\beta_t)(\beta - \beta_t) = \mathbf 0\) for \(\beta\). Write the update \(\beta_{t+1}\) in terms of \(\mathbf X\), \(\mathbf W_t\), \(\mathbf y\) and \(\mathbf p_t\).

  1. \(\log p_i = -\log(1 + e^{-\mathbf x_i^T\beta})\) and \(\log(1 - p_i) = -\mathbf x_i^T\beta - \log(1 + e^{-\mathbf x_i^T\beta})\). For the gradient, \(\tfrac{d}{dz}\log(1 + e^z) = e^z/(1 + e^z)\), which is \(p_i\) at \(z = \mathbf x_i^T\beta\).
  1. \(\tfrac{\partial}{\partial\beta} p_i(\beta) = p_i(1 - p_i)\,\mathbf x_i\). (c) \(v^T\mathbf x_i\mathbf x_i^T v = (\mathbf x_i^T v)^2\).

Exercise 4 (no free lunch on a finite domain). A finite set \(\mathcal X\) has \(N\) points. A labeling is a function \(f: \mathcal X \to \{0, 1\}\). A method sees the labels on a fixed training set \(S \subset \mathcal X\) of \(n < N\) points and returns a classifier \(\hat h\), which may be any function of those \(n\) labels. Its error on the points it has not seen is \[ \mathrm{err}(\hat h, f) = \frac{1}{N - n} \sum_{x \notin S} \mathbb 1\big[\hat h(x) \ne f(x)\big]. \]

  1. Draw \(f\) uniformly at random from all \(2^N\) labelings. Show that \(\mathbb E\big[\mathrm{err}(\hat h, f)\big] = \tfrac12\) for every method.

  2. Two methods, \(A\) and \(B\). Suppose \(A\)’s error on the unseen points is smaller than \(B\)’s for some labeling. Show that \(B\)’s is smaller than \(A\)’s for another labeling.

  3. What would have to be true of the world for some method to do better than \(\tfrac12\) on average? Where does that assumption sit in the ERM form of exercise 2?

Uniform over all \(2^N\) labelings means that the \(N\) labels \(f(x)\) are independent fair coin flips.

Fix \(x \notin S\). \(\hat h(x)\) depends only on the labels on \(S\), and \(f(x)\) is independent of them. Condition on the labels on \(S\): what is \(P\big(\hat h(x) \ne f(x)\big)\)?

  1. Average \(\mathrm{err}(\hat h_A, f) - \mathrm{err}(\hat h_B, f)\) over all \(2^N\) labelings.

Problems

Problem 1 (three observations). Three observations, \((x_i, y_i) = (-1, 0), (0, 1), (1, 3)\), a line \(\beta_0 + \beta_1 x\), and the objective \(\hat R_n(\beta) = \tfrac13\sum_i (y_i - \beta_0 - \beta_1 x_i)^2\).

  1. Compute the gradient of \(\hat R_n\) at \(\beta = (0, 0)\) and take one gradient-descent step with \(\gamma = 0.1\).

  2. At what \(\beta\) does the iteration stop moving, and which equations does that point satisfy?

  3. The curvatures of \(\hat R_n\) are the eigenvalues of \(2\mathbf X^T\mathbf X/n\). Compute them. What is the largest step size for which gradient descent converges from every start, and is \(\gamma = 0.5\) safe?

  4. Which coordinate converges more slowly with \(\gamma = 0.5\), and by what factor does its error shrink per step?

Problem 2 (one direction, then two). One direction with curvature \(a > 0\): \(f(\beta) = \tfrac12 a\beta^2\).

  1. Show that gradient descent with step \(\gamma\) gives \(\beta_t = (1 - \gamma a)^t\beta_0\), and find all \(\gamma\) for which it converges.

  2. Which \(\gamma\) reaches the minimum in one step?

  3. Two directions with curvatures \(10\) and \(1\) and one shared step size. Which \(\gamma\) makes the slower of the two directions shrink fastest, and what is the shrinkage factor per step then?

Problem 3 (a random row). \(\hat R_n(\beta) = \tfrac1n\sum_{i=1}^n \ell_i(\beta)\). Pick an index \(I\) uniformly at random from \(\{1, \dots, n\}\), with the data and \(\beta\) held fixed.

  1. Show that \(\mathbb E_I[\nabla\ell_I(\beta)] = \nabla\hat R_n(\beta)\). Which assumption about where the data came from did you use?

  2. A batch of \(b\) indices drawn independently with replacement, gradient averaged: give its mean and its variance in terms of \(\sigma^2(\beta) = \mathrm{Var}_I[\nabla\ell_I(\beta)]\) (one coordinate at a time is fine).

  3. If the batch is drawn without replacement, is the variance larger or smaller, and what is it when \(b = n\)?

  4. The descent lemma says a full gradient step of size \(1/L\) decreases \(\hat R_n\). Does a step of size \(1/L\) along \(\nabla\ell_I(\beta)\) decrease \(\hat R_n\)? Give a reason or a counterexample.

Problem 4 (the valley with shrinking steps). Gradient descent on the lecture’s valley, \(f(\beta) = \tfrac12(\beta_1^2 + 10\beta_2^2)\), from \(\beta_0 = (-4, 1.2)\) with the step sizes \(\gamma_t = 0.3 \cdot 0.6^t\), \(t = 1, 2, \dots\) (panel C of the lecture).

  1. Show that each coordinate obeys \(\beta_{j, t} = \beta_{j, 0}\prod_{s=1}^t (1 - a_j\gamma_s)\) with \(a_1 = 1\), \(a_2 = 10\).

  2. For the first coordinate, compute the first four factors and say what \(\beta_{1, t}\) converges to (two decimals are enough). Is it the minimizer?

  3. Compute \(\sum_{t \ge 1}\gamma_t\). Which of the two conditions on step sizes from the lecture does this schedule violate, and what does the sum say about how far the iterate can travel?

Coding

1. Gradient descent, written once

Life expectancy on \(\log_{10}\) GDP per capita in 2007, as in week 2. The input is centered; step 4 shows why.

g07 <- gapminder |> filter(year == 2007) |> mutate(lgdp = log10(gdpPercap))
y <- g07$lifeExp
X <- cbind(1, g07$lgdp - mean(g07$lgdp))
n <- nrow(X)

Step 1. Write the objective and its gradient as functions of X, y and beta: \(\hat R_n(\beta) = \tfrac1n\|\mathbf y - \mathbf X\beta\|^2\) and \(\nabla\hat R_n(\beta) = -\tfrac2n\mathbf X^T(\mathbf y - \mathbf X\beta)\).

ls_loss <- function(X, y, beta) {
  # ...
}
ls_grad <- function(X, y, beta) {
  # ...
}

Step 2. The loop, as a function you will reuse: start from beta0, take steps steps of size gamma along the negative gradient, and record the loss and the gradient norm after each step. It takes the loss and gradient functions as arguments, so that the same loop fits a different model when those two change.

gradient_descent <- function(X, y, loss, grad, gamma, steps, beta0 = rep(0, ncol(X))) {
  # ...
  list(beta = as.vector(beta), losses = losses, grad_norms = grad_norms)
}

Checkpoint (run it; if it stops, fix step 1 or 2 before going on):

chk <- gradient_descent(X, y, ls_loss, ls_grad, gamma = 0.3, steps = 50)
stopifnot(length(chk$losses) == 50, all(diff(chk$losses) <= 1e-8),
          max(abs(chk$beta - coef(lm(y ~ X - 1)))) < 0.01)

Step 3. Run it with gamma = 0.3 for 50 steps and compare beta with coef(lm(y ~ X - 1)). Then run gamma = 0.05 and gamma = 1.05, and plot the three loss curves against the step (log scale on the vertical axis) and the three gradient norms.

fits <- list(`0.05` = ____, `0.3` = ____, `1.05` = ____)
# a data frame with one row per (gamma, step), then two plots
# ...

Question 1. The curvatures here, the eigenvalues of \(2\mathbf X^T\mathbf X/n\), are \(2\) and about \(0.69\). Using problem 2, explain what each of the three runs did, and what the gradient norm shows that the loss does not.

Step 4. Now uncenter the input, X_raw <- cbind(1, g07$lgdp). The curvatures are now about \(30.6\) and \(0.045\). Choose a safe step (say \(\gamma = 1/30\)), run 300 steps, and compare with lm again.

X_raw <- cbind(1, g07$lgdp)
# ...

Question 2. The fitted line is the same as before (the slope is unchanged; only the intercept’s meaning moved). Why does gradient descent now need thousands of steps for the same problem, and what does that say about scaling inputs before descent?

Step 5. Logistic regression by the same loop. Chinstrap against Adelie from flipper length (standardized), as in the lecture’s penguin fit. Write the negative log-likelihood per bird, \(\hat R_n(\beta) = -\ell(\beta)/n\), and its gradient \(-\tfrac1n\mathbf X^T(\mathbf y - \mathbf p(\beta))\), and run gradient_descent with gamma = 4 for 40 steps. Compare with glm, and check the score equations at your answer: \(\mathbf X^T(\mathbf y - \hat{\mathbf p})\) should be zero.

pg <- palmerpenguins::penguins |>   # modeldata also has a penguins; be explicit
  filter(species != "Gentoo") |>
  drop_na(flipper_length_mm) |>
  mutate(chinstrap = as.numeric(species == "Chinstrap"),
         z = as.numeric(scale(flipper_length_mm)))
X_pg <- cbind(1, pg$z)
y_pg <- pg$chinstrap
logit_loss <- function(X, y, beta) {
  # ...
}
logit_grad <- function(X, y, beta) {
  # ...
}
# fit_logit <- gradient_descent(____)
# ...

Checkpoint:

stopifnot(max(abs(fit_logit$beta - coef(glm(chinstrap ~ z, family = binomial, data = pg)))) < 1e-3)

Question 3. Exactly two things changed between the least-squares run and this one. Name them. Why is a step of \(4\) safe here when \(1.05\) was not for least squares?

Step 6. The gradient without the calculus. Write numeric_grad(loss, X, y, beta, h = 1e-6), which returns the vector whose \(j\)th entry is the central difference \[ \frac{\hat R_n(\beta + h\,e_j) - \hat R_n(\beta - h\,e_j)}{2h}, \] with \(e_j\) the \(j\)th coordinate vector. Check it against ls_grad and logit_grad at some \(\beta\) of your choice, then fit the penguin model again with gradient_descent, passing a gradient function that calls numeric_grad on logit_loss.

numeric_grad <- function(loss, X, y, beta, h = 1e-6) {
  # ...
}
# fit_num <- gradient_descent(X_pg, y_pg, logit_loss, grad = ____, gamma = 4, steps = 40)

Checkpoint:

stopifnot(max(abs(numeric_grad(logit_loss, X_pg, y_pg, c(0.5, -0.3)) - logit_grad(X_pg, y_pg, c(0.5, -0.3)))) < 1e-6)

Question. To fit a different model with what you now have, say a Poisson regression or a neural network, what would you have to write, and what would it cost you per step?

2. Many solutions, and stopping early

Step 1. More coefficients than observations: \(n = 100\), \(p = 400\). Implement mini-batch stochastic gradient descent for least squares from beta <- rep(0, p): in each epoch, shuffle the rows (sample(n)) and, for each batch of 10 rows, take a step of size gamma = 0.005 along the batch’s gradient. Record after each epoch the training MSE and the distance \(\|\beta_t - \beta\|\) to the true coefficients.

set.seed(1)
n_big <- 100
p_big <- 400
X_big <- matrix(rnorm(n_big * p_big), n_big, p_big)
beta_true <- rpois(p_big, 1.5)
y_big <- X_big %*% beta_true + rnorm(n_big, sd = 2)
epochs <- 40
batch_size <- 10
gamma <- 0.005
beta_sgd <- rep(0, p_big)
sgd_path <- tibble(epoch = 1:epochs, train_mse = NA_real_, dist_to_truth = NA_real_)
for (e in 1:epochs) {
  # ...
  sgd_path$train_mse[e] <- ls_loss(X_big, y_big, beta_sgd)
  sgd_path$dist_to_truth[e] <- sqrt(sum((beta_sgd - beta_true)^2))
}
# ...

Question 4. The training MSE goes to zero. Why is that possible, and why should it worry you? The distance to the true coefficients levels off at about \(33\), from a starting distance of \(38\): what is it that the algorithm cannot learn from these data?

Step 2. The lecture’s wide regression: \(n = 80\), \(p = 60\), noisy, and a held-out set of 2,000 rows from the same model. Run full-batch gradient descent from zero for 300 steps with step 0.5 / L, where L is twice the largest eigenvalue of \(\mathbf X^T\mathbf X/n\), and record the training and the held-out MSE after every step. Plot both against the step (log scale on the horizontal axis), and find the step at which the held-out MSE is smallest.

set.seed(12)
n_w <- 80
p_w <- 60
X_w <- matrix(rnorm(n_w * p_w), n_w, p_w)
beta_w <- rnorm(p_w, sd = 0.5)
y_w <- X_w %*% beta_w + rnorm(n_w, sd = 3)
X_held <- matrix(rnorm(2000 * p_w), 2000, p_w)
y_held <- X_held %*% beta_w + rnorm(2000, sd = 3)
# L <- ____
beta_gd <- rep(0, p_w)
es_path <- tibble(step = 1:300, training = NA_real_, held_out = NA_real_)
# ...

Question 5. The training error falls at every step, as the descent lemma says it must. Where would you have stopped if you had only seen it? You chose the number of steps by looking at the held-out MSE: can you now report that held-out MSE as an honest estimate of the chosen model’s error? What would you need instead?

Optional, can be done at home: which of the many solutions

Compare the stochastic gradient descent estimate from part 2, step 1, with week 2’s minimum-norm least-squares solution, MASS::ginv(X_big) %*% y_big, and compute the norms of the two.

# ...