Week 5 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(splines)         # ns(); ships with R
library(glmnet)          # might need to install.packages("glmnet")
library(rpart)           # might need to install.packages("rpart")
library(randomForest)    # might need to install.packages("randomForest")
library(nnet)            # one-hidden-layer networks; ships with R
theme_set(theme_minimal())
mse <- function(y, pred) mean((y - pred)^2)
se_loss <- function(y, pred) sd((y - pred)^2) / sqrt(length(y))   # standard error of an MSE
# k-NN regression with several standardized inputs, as in the lecture: squared distances by
# matrix algebra, then the mean outcome of the k nearest training rows.
knn_predict <- function(Xnew, Xtrain, ytrain, k) {
  d2 <- outer(rowSums(Xnew^2), rowSums(Xtrain^2), "+") - 2 * Xnew %*% t(Xtrain)
  apply(d2, 1, function(d) mean(ytrain[order(d)[1:k]]))
}

Before the seminar

The lecture stated one result, the optimism theorem, and left its proof to you (exercise 1). It used ridge regression without deriving its estimator (exercise 2). Do both on paper before the seminar; the seminar’s problems are variations on them. You may use, without proof: linearity of expectation; the tower property, \(\mathbb E[W] = \mathbb E\big[\mathbb E[W \mid V]\big]\); that for a symmetric matrix \(\mathbf A\) the gradient of \(\beta^T\mathbf A\beta\) is \(2\mathbf A\beta\) and its Hessian is \(2\mathbf A\) (week 4); and that a symmetric matrix has real eigenvalues and an orthonormal basis of eigenvectors. Throughout, assume that every expectation that appears exists and every minimizer named 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. Solutions are shown in the seminar if you ask.

Exercise 1 (optimism theorem). Training data \((x_i, y_i)\), \(i = 1, \ldots, n\), independent draws from the population of \((X, Y)\); a loss \(\ell\); the risk \(R(f) = \mathbb E\,\ell(Y, f(X))\) and the training error \(\hat R_n(f) = \tfrac1n\sum_{i=1}^n \ell(y_i, f(x_i))\). A class \(\mathcal F\) is fixed before the data; \(f_{\mathcal F}\) minimizes \(R\) over \(\mathcal F\); \(\hat f\) minimizes \(\hat R_n\) over \(\mathcal F\) (the empirical risk minimizer). Expectations are over the training data.

  1. Let \(f\) be any fixed function, chosen without looking at the data. Show that \(\mathbb E[\hat R_n(f)] = R(f)\). Which part of “independent draws from the population” did you use, and which part did you not use?

  2. Show that \(\hat R_n(\hat f) \le \hat R_n(f_{\mathcal F})\) for every training set, and deduce \(\mathbb E[\hat R_n(\hat f)] \le R(f_{\mathcal F})\).

  3. Show that \(R(f_{\mathcal F}) \le R(\hat f)\) for every training set, and deduce \(R(f_{\mathcal F}) \le \mathbb E[R(\hat f)]\). Of the two inequalities in the theorem, which used the sampling assumption and which used only the definitions?

  4. A test set \((x'_j, y'_j)\), \(j = 1, \ldots, m\), of independent draws from the population, independent of the training data, scores \(\hat f\) by \(\hat R_{\text{test}}(\hat f) = \tfrac1m\sum_j \ell(y'_j, \hat f(x'_j))\). Show that \(\mathbb E[\hat R_{\text{test}}(\hat f)] = \mathbb E[R(\hat f)]\): the test score lands, on average, on the right end of the theorem.

  1. \(\hat R_n(f)\) is an average of \(n\) terms; take the expectation of each term. What is the distribution of \((x_i, y_i)\)?
  1. \(\hat f\) minimizes \(\hat R_n\) over \(\mathcal F\), and \(f_{\mathcal F}\) is a member of \(\mathcal F\). Then apply (a) to \(f_{\mathcal F}\): is it a fixed function?
  1. Condition on the training data. Given them, \(\hat f\) is a fixed function and the test rows are still draws from the population; apply (a), then the tower property.

Exercise 2 (ridge regression’s normal equations). \(\mathbf X\) is \(n \times p\) with centered columns and \(\mathbf y\) is centered, so there is no intercept to fit (the lecture’s convention: the intercept is not penalized). For \(\lambda > 0\) the ridge estimator minimizes \[ Q_\lambda(\beta) = \|\mathbf y - \mathbf X\beta\|^2 + \lambda\|\beta\|^2. \]

  1. Expand \(Q_\lambda\), compute its gradient, and show that a stationary point satisfies \((\mathbf X^T\mathbf X + \lambda\mathbf I)\beta = \mathbf X^T\mathbf y\).

  2. Show that \(\mathbf X^T\mathbf X + \lambda\mathbf I\) is positive definite, hence invertible, for every \(\lambda > 0\), whatever \(n\) and \(p\) are. Conclude that \(\hat\beta_\lambda = (\mathbf X^T\mathbf X + \lambda\mathbf I)^{-1}\mathbf X^T\mathbf y\) is the unique stationary point, even when \(p > n\).

  3. Compute the Hessian of \(Q_\lambda\). In week 4’s language, what are its smallest and largest curvatures, in terms of the eigenvalues \(\mu_1 \ge \dots \ge \mu_p \ge 0\) of \(\mathbf X^T\mathbf X\)? Why is the stationary point the global minimizer, and what does the penalty do to the ratio of largest to smallest curvature that week 4 showed sets the speed of gradient descent?

  4. Show that \(\|\hat\beta_\lambda\| \le \|\mathbf X^T\mathbf y\|/\lambda\), so \(\hat\beta_\lambda \to \mathbf 0\) as \(\lambda \to \infty\). What is \(\hat\beta_\lambda\) at \(\lambda = 0\) when \(\mathbf X^T\mathbf X\) is invertible?

  1. \(\|\mathbf y - \mathbf X\beta\|^2 = \mathbf y^T\mathbf y - 2\beta^T\mathbf X^T\mathbf y + \beta^T\mathbf X^T\mathbf X\beta\), and \(\|\beta\|^2 = \beta^T\mathbf I\beta\). Differentiate term by term.
  1. For any \(v\), \(v^T\mathbf X^T\mathbf X v = \|\mathbf X v\|^2\). What is \(v^T(\mathbf X^T\mathbf X + \lambda\mathbf I)v\) when \(v \ne \mathbf 0\)?

(c), (d) \(\mathbf X^T\mathbf X\) and \(\mathbf X^T\mathbf X + \lambda\mathbf I\) have the same eigenvectors; the eigenvalues of the second are \(\mu_j + \lambda\). Write \(\mathbf X^T\mathbf y\) in the eigenvector basis.

Problems

Problem 1 (the selection winner’s score). \(K\) candidate predictors \(f_1, \ldots, f_K\) have been fitted on the training rows, say by \(K\) methods or one method at \(K\) values of its knob. A selection set of \(m\) rows, independent draws from the population and independent of the training rows, scores them: \(S_k = \tfrac1m\sum_j \ell(y_j, f_k(x_j))\). The winner is \(\hat k = \arg\min_k S_k\), and \(S_{\hat k}\) is its score, the number at the top of the leaderboard.

  1. Show that \(\mathbb E[S_k \mid \text{training data}] = R(f_k)\) for each \(k\).

  2. Show that \(\mathbb E[S_{\hat k} \mid \text{training data}] \le \min_k R(f_k) \le \mathbb E[R(f_{\hat k}) \mid \text{training data}]\). Which argument from exercise 1 is this, with what in place of \(\mathcal F\) and of \(\hat R_n\)?

  3. Two candidates with the same risk \(r\), whose scores are independent and each equal to \(r - 1\) or \(r + 1\) with probability \(\tfrac12\). Compute \(\mathbb E[S_{\hat k}]\) and \(\mathbb E[R(f_{\hat k})]\). Now \(K\) such candidates: what is \(\mathbb E[S_{\hat k}]\), and what happens as \(K\) grows?

  4. The winner is then scored once on the test set. Why does that score not have the same problem?

Problem 2 (where the proof breaks). The lecture said the optimism theorem covers least squares and ridge or lasso with a fixed radius, and not ridge at a fixed \(\lambda\), a greedily grown tree, or a knob chosen on the data.

  1. A population in which \(Y = X\) exactly (no noise), \(\mathbb E X^2 > 0\), squared loss, and the class \(\mathcal F = \{x \mapsto \beta x : \beta \in \mathbb R\}\) of lines through the origin. Find \(f_{\mathcal F}\) and \(R(f_{\mathcal F})\). For a sample \(x_1, \ldots, x_n\) (not all zero) and a fixed \(\lambda > 0\), find the ridge estimate \(\hat\beta_\lambda\) and its training error \(\hat R_n(\hat\beta_\lambda)\), and show that \(\mathbb E[\hat R_n(\hat\beta_\lambda)] > R(f_{\mathcal F})\): the left inequality fails. Which hypothesis of the theorem does \(\hat\beta_\lambda\) violate?

  2. Replace the penalty by the constraint \(|\beta| \le t\) for a fixed \(t > 0\): \(\mathcal F_t = \{x \mapsto \beta x : |\beta| \le t\}\), and \(\hat f\) the minimizer of \(\hat R_n\) over \(\mathcal F_t\). Does the theorem apply? What is \(\hat f\) here, for \(t \ge 1\) and for \(t < 1\)?

  3. A tree of depth \(2\) grown greedily (the lecture’s splitting) is a member of the class of all depth-\(2\) trees on the same inputs. Which line of your proof of exercise 1 is not available for it? Does that mean its training error must sit above its risk on average, or only that the theorem does not say?

  4. Which of the theorem’s two inequalities holds for every \(\hat f \in \mathcal F\), whatever algorithm produced it?

Problem 3 (ridge as least squares on more rows). Stack \(p\) extra rows under the data: \(\tilde{\mathbf X} = \begin{pmatrix}\mathbf X \\ \sqrt\lambda\,\mathbf I_p\end{pmatrix}\), an \((n + p) \times p\) matrix, and \(\tilde{\mathbf y} = \begin{pmatrix}\mathbf y \\ \mathbf 0_p\end{pmatrix}\).

  1. Show that \(\|\tilde{\mathbf y} - \tilde{\mathbf X}\beta\|^2 = Q_\lambda(\beta)\), the ridge objective of exercise 2.

  2. Write week 2’s normal equations for the regression of \(\tilde{\mathbf y}\) on \(\tilde{\mathbf X}\) and check that they are exercise 2’s. Show that \(\tilde{\mathbf X}\) has full column rank for every \(\lambda > 0\), even when \(p > n\).

  3. In words: what does each of the \(p\) extra rows say, and what does \(\sqrt\lambda\) do? How would you fit ridge regression with nothing but lm()?

Coding

1. The task and the split

The lecture’s task: predict a track’s popularity from ten audio features, its release year and its genre, squared loss (data note). The three-way split with the lecture’s seed, so that your numbers match its leaderboards.

sp_file <- "https://ml4ds.com/files/spotify_tracks.csv"
# the course file; a local copy is used if one is beside this notebook's folder
for (p in c("../../../files/spotify_tracks.csv", "../../../../files/spotify_tracks.csv")) if (file.exists(p)) sp_file <- p
sp <- read_csv(sp_file, show_col_types = FALSE) |>
  distinct(track_id, .keep_all = TRUE) |>          # one row per track
  mutate(duration_min = duration_ms / 60000, genre = factor(playlist_genre))
feat <- c("danceability", "energy", "loudness", "speechiness", "acousticness",
          "instrumentalness", "liveness", "valence", "tempo", "duration_min", "release_year")
f_plane <- reformulate(c(feat, "genre"), "track_popularity")   # every input, as a formula

Step 1. Split the 10,000 tracks at random into 5,000 training, 2,000 selection and 3,000 test tracks. The test tracks are the vault: nothing before part 3 touches them.

set.seed(5)
sp_id <- sample(nrow(sp))        # a random order of the 10,000 tracks
# sp_train <- sp[sp_id[____], ]  # 5,000 tracks
# sp_sel <- ____                 # 2,000
# sp_test <- ____                # 3,000

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

stopifnot(nrow(sp_train) == 5000, nrow(sp_sel) == 2000, nrow(sp_test) == 3000,
          length(intersect(sp_train$track_id, c(sp_sel$track_id, sp_test$track_id))) == 0,
          length(intersect(sp_sel$track_id, sp_test$track_id)) == 0)

Question 1. The file has 11,554 rows for 10,000 tracks, because 1,101 tracks are in more than one playlist, and the split above is by track. If you had skipped distinct() and split the rows, what could happen to a track, and which of the lecture’s rules would the test score then break?

Step 2. The baselines on the selection tracks, as in the lecture. score() builds one leaderboard row: a competitor’s selection MSE and its standard error. The mean is done; add the plane (least squares on every input) and \(k\) nearest neighbors on the eleven numeric inputs, standardized, with \(k\) chosen on the selection tracks from ks. design() builds a standardized matrix of inputs, with genre as five 0/1 columns when asked, each column centered and scaled by the reference set’s mean and standard deviation.

y_tr <- sp_train$track_popularity
y_se <- sp_sel$track_popularity
score <- function(name, pred_sel)
  tibble(competitor = name, sel_mse = mse(y_se, pred_sel), se = se_loss(y_se, pred_sel))
design <- function(d, ref, with_genre = TRUE) {
  f <- reformulate(if (with_genre) c(feat, "genre") else feat, "track_popularity")
  X <- model.matrix(f, d)[, -1]
  R <- model.matrix(f, ref)[, -1]
  scale(X, colMeans(R), apply(R, 2, sd))
}
ks <- c(1, 5, 10, 20, 50, 100, 200)
# fit_plane <- lm(____)
X_tr <- design(sp_train, sp_train, with_genre = FALSE)
X_se <- design(sp_sel, sp_train, with_genre = FALSE)
# knn_sel <- sapply(ks, function(k) ____)      # the selection MSE at each k
# k_best <- ____
lb <- bind_rows(
  score("the training mean", rep(mean(y_tr), nrow(sp_sel)))
  # , score("least squares, every input", ____)
  # , score(paste("k-NN, numeric inputs, k =", k_best), ____)
)
lb

2. Five structures, each tuned on the selection tracks

The pattern is the same for every family: fit on the training tracks at each value of the knob, score each fit on the selection tracks, keep the best. Each step adds one row to lb.

Step 3. The lasso. glmnet() fits the whole path of \(\lambda\) values at once, on the standardized matrix with the genre columns; predict() on the selection matrix then gives one column of predictions per \(\lambda\). Choose \(\lambda\), and count the coefficients that are exactly zero at it.

Xg_tr <- design(sp_train, sp_train)
Xg_se <- design(sp_sel, sp_train)
lasso_fit <- glmnet(Xg_tr, y_tr, standardize = FALSE)   # the columns already are
# lasso_sel <- apply(predict(lasso_fit, Xg_se), 2, function(pred) ____)
# lam_best <- lasso_fit$lambda[____]
# coefficients at lam_best: coef(lasso_fit, s = lam_best); how many are zero?
# ...
# lb <- bind_rows(lb, score("lasso, every input", ____))

Step 4. Additive splines: lm() on a natural-spline basis of size \(J\) for each numeric input, plus genre, as in the lecture. gam_formula(J) writes the formula. Choose \(J\) from Js.

gam_formula <- function(J) reformulate(c(sprintf("ns(%s, df = %d)", feat, J), "genre"), "track_popularity")
Js <- c(2, 3, 4, 6)
# gam_sel <- sapply(Js, function(J) ____)
# J_best <- ____
# fit_gam <- ____
# lb <- bind_rows(lb, score(paste("additive splines, J =", J_best), ____))

Step 5. A tree and a forest. tree_at(d, data) grows a tree of depth \(d\); choose the depth from depths. The forest has no knob today: randomForest() with its defaults, 300 trees, seeded so that everyone grows the same forest (it takes a minute).

tree_at <- function(d, data) rpart(f_plane, data, control = rpart.control(maxdepth = d, cp = 0, minbucket = 10))
depths <- c(2, 3, 4, 6, 8)
# tree_sel <- sapply(depths, function(d) ____)
# d_best <- ____
# fit_tree <- tree_at(d_best, sp_train)
set.seed(1)
fit_rf <- randomForest(f_plane, sp_train, ntree = 300)
# lb <- bind_rows(lb, score(paste("tree, depth", d_best), ____), score("forest, defaults", ____))

Step 6. A network with one hidden layer of \(M\) units, nnet() on the standardized matrix (linout = TRUE for a numeric outcome; decay is a small ridge penalty on the weights, fixed today). net_at() seeds each fit, since the optimization starts from random weights. Choose \(M\) from Ms.

net_at <- function(M, X, y) {
  set.seed(10)
  nnet(X, y, size = M, decay = 0.1, linout = TRUE, maxit = 300, trace = FALSE)
}
Ms <- c(2, 4, 8, 16)
# net_sel <- sapply(Ms, function(M) ____)
# M_best <- ____
# fit_net <- ____
# lb <- bind_rows(lb, score(paste("network, M =", M_best), ____))

The selection leaderboard:

lb |> arrange(sel_mse) |> mutate(across(c(sel_mse, se), round))
ggplot(lb, aes(sel_mse, fct_reorder(competitor, -sel_mse))) +
  geom_errorbar(aes(xmin = sel_mse - se, xmax = sel_mse + se), width = 0.3, colour = "grey60") +
  geom_point(size = 3) +
  labs(x = "MSE on the 2,000 selection tracks", y = NULL, caption = "bars: one standard error")

Question 2. Which competitors lead, and by how much against the plane with the same inputs? The bars are each score’s own standard error: can they tell the leaders apart? Which rows of this leaderboard are unbiased estimates of their competitor’s risk, and which are not (problem 1)?

3. The vault

Step 7. Refit every competitor on the training and selection tracks together, each at the knob it chose, and predict the test tracks. The standardized matrices are rebuilt with the 7,000 rows as the reference. \(k\)-NN has nothing to refit: its “model” is the rows it averages over, now 7,000 of them.

sp_both <- bind_rows(sp_train, sp_sel)
y_both <- sp_both$track_popularity
y_te <- sp_test$track_popularity
X_both <- design(sp_both, sp_both, with_genre = FALSE); X_te <- design(sp_test, sp_both, with_genre = FALSE)
Xg_both <- design(sp_both, sp_both); Xg_te <- design(sp_test, sp_both)
pred_te <- list(
  "the training mean" = rep(mean(y_both), nrow(sp_test))
  # , "least squares, every input" = ____
  # , "k-NN, numeric inputs" = ____
  # , "lasso" = ____
  # , "additive splines" = ____
  # , "tree" = ____
  # , "forest" = ____
  # , "network" = ____
)

Checkpoint:

stopifnot(all(sapply(pred_te, length) == nrow(sp_test)))

Step 8. Open the vault once: every competitor’s test MSE with its standard error, in one table. Then the paired comparison of the two leaders: the per-track difference of their squared errors, its mean and its standard error.

# vault <- imap_dfr(pred_te, ~ tibble(competitor = .y, test_mse = ____, se = ____)) |> arrange(test_mse)
# the two leaders, paired on the 3,000 test tracks:
# d <- ____
# c(mean = mean(d), se = ____)

Question 3. Compare the vault with the selection leaderboard: the order, and each competitor’s test score against its selection score. Problem 1 says a tuned competitor’s selection score is too good on average. Can you see that here? What else changed between the two leaderboards?

Question 4. The forest’s training error against its test error, and the plane’s: compute both and compare the gaps. Which competitor does the optimism theorem cover, and does the uncovered one behave any differently?

# c(forest_train = ____, forest_test = ____, plane_train = ____, plane_test = ____)

Optional, can be done at home: ridge three ways

Problem 3 says ridge regression is least squares on an augmented dataset. Check it on the Spotify training matrix: fit ridge with \(\lambda = 100\) by the formula of exercise 2, by lm() on the augmented rows, and compare the two coefficient vectors.

lambda <- 100
X_c <- scale(Xg_tr, scale = FALSE)          # centered columns, as in exercise 2
y_c <- y_tr - mean(y_tr)
# b_formula <- ____
# X_aug <- rbind(X_c, ____); y_aug <- c(y_c, ____)
# b_lm <- coef(lm(y_aug ~ X_aug - 1))
# max(abs(b_formula - b_lm))

solve(A, b) solves \(\mathbf A\beta = \mathbf b\). diag(p) is \(\mathbf I_p\). The augmented fit has no intercept: - 1 in the formula.

Optional, can be done at home: the lasso path on the tracks

The lecture drew the lasso’s solution path on these tracks. Draw it from your lasso_fit with plot(lasso_fit, xvar = "lambda", label = TRUE), and find the step of the path at which each input enters as \(\lambda\) falls (large \(\lambda\) is on the right in this plot). Where are the release year and loudness, the two inputs the lecture’s component plots bent most?