Week 5 seminar (complete)

Author

Joshua Loftus 

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\); 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.

  1. By linearity, \(\mathbb E[\hat R_n(f)] = \tfrac1n\sum_i \mathbb E\,\ell(y_i, f(x_i))\). Each \((x_i, y_i)\) has the distribution of \((X, Y)\) and \(f\) does not depend on the data, so each term is \(\mathbb E\,\ell(Y, f(X)) = R(f)\), and the average of \(n\) copies of \(R(f)\) is \(R(f)\). Used: that every row is distributed as the population. Not used: independence between rows; a fixed function’s training error is unbiased for its risk on any sample whose rows each have the population distribution. (Independence is what makes the standard error of a training or test score computable from the rows, and it is used in (d).)

  2. \(f_{\mathcal F} \in \mathcal F\) and \(\hat f\) minimizes \(\hat R_n\) over \(\mathcal F\), so \(\hat R_n(\hat f) \le \hat R_n(f_{\mathcal F})\) on every training set. Taking expectations, \(\mathbb E[\hat R_n(\hat f)] \le \mathbb E[\hat R_n(f_{\mathcal F})] = R(f_{\mathcal F})\), the last step by (a): \(f_{\mathcal F}\) is defined by the population and the class, not by the sample, so it is a fixed function.

  3. \(\hat f \in \mathcal F\) and \(f_{\mathcal F}\) minimizes \(R\) over \(\mathcal F\), so \(R(f_{\mathcal F}) \le R(\hat f)\) on every training set; take expectations. \(\blacksquare\) The right inequality used only the definitions of \(f_{\mathcal F}\) and \(\mathcal F\); it holds for any \(\hat f \in \mathcal F\), however it was chosen. The left one used (a), and so the sampling assumption, together with the fact that \(\hat f\) is the minimizer of \(\hat R_n\).

  4. Condition on the training data. Given them, \(\hat f\) is one fixed function, and the test rows are draws from the population that do not depend on it, so (a) gives \(\mathbb E[\hat R_{\text{test}}(\hat f) \mid \text{training data}] = R(\hat f)\). By the tower property, \(\mathbb E[\hat R_{\text{test}}(\hat f)] = \mathbb E\big[R(\hat f)\big]\). Independence of the test rows from the training data is what this step needs.

Common errors: applying (a) to \(\hat f\) itself (it is not a fixed function; that is the whole point); writing \(R(\hat f)\) as if it were a number (it is random, through \(\hat f\)); reading the theorem as “\(\hat R_n(\hat f) \le R(\hat f)\) on every sample” (it is a statement about averages over training sets).

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\). 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, the ratio that 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.

  1. \(Q_\lambda(\beta) = \mathbf y^T\mathbf y - 2\beta^T\mathbf X^T\mathbf y + \beta^T\mathbf X^T\mathbf X\beta + \lambda\beta^T\beta\), so \(\nabla Q_\lambda(\beta) = -2\mathbf X^T\mathbf y + 2\mathbf X^T\mathbf X\beta + 2\lambda\beta = 2\big[(\mathbf X^T\mathbf X + \lambda\mathbf I)\beta - \mathbf X^T\mathbf y\big]\). Setting it to zero gives \((\mathbf X^T\mathbf X + \lambda\mathbf I)\beta = \mathbf X^T\mathbf y\): the least-squares normal equations with \(\lambda\) added down the diagonal.

  2. For \(v \ne \mathbf 0\), \(v^T(\mathbf X^T\mathbf X + \lambda\mathbf I)v = \|\mathbf X v\|^2 + \lambda\|v\|^2 \ge \lambda\|v\|^2 > 0\). A positive definite matrix has no null vector, so it is invertible, and the stationary point is unique: \(\hat\beta_\lambda = (\mathbf X^T\mathbf X + \lambda\mathbf I)^{-1}\mathbf X^T\mathbf y\). Nothing here needed \(n \ge p\): when \(p > n\), \(\mathbf X^T\mathbf X\) has a null space of dimension at least \(p - n\) and least squares has infinitely many solutions, but \(\lambda\|v\|^2\) is positive on that null space too. \(\blacksquare\)

  3. The Hessian is \(2(\mathbf X^T\mathbf X + \lambda\mathbf I)\), with eigenvalues \(2(\mu_j + \lambda)\): smallest curvature \(2(\mu_p + \lambda) \ge 2\lambda > 0\), largest \(2(\mu_1 + \lambda)\). A positive definite Hessian everywhere makes \(Q_\lambda\) strictly convex, so its one stationary point is the global minimizer. The ratio of largest to smallest curvature falls from \(\mu_1/\mu_p\) (infinite when \(\mu_p = 0\), a direction of zero curvature, as in a regression with more columns than rows) to \((\mu_1 + \lambda)/(\mu_p + \lambda)\): the penalty rounds the valley, and the flat directions, which the data say nothing about, are the ones it shrinks hardest.

  4. In the orthonormal eigenvector basis \(u_1, \ldots, u_p\) of \(\mathbf X^T\mathbf X\), write \(\mathbf X^T\mathbf y = \sum_j c_j u_j\). Then \(\hat\beta_\lambda = \sum_j \frac{c_j}{\mu_j + \lambda} u_j\), so \(\|\hat\beta_\lambda\|^2 = \sum_j \frac{c_j^2}{(\mu_j + \lambda)^2} \le \frac{1}{\lambda^2}\sum_j c_j^2 = \|\mathbf X^T\mathbf y\|^2/\lambda^2\). As \(\lambda \to \infty\) the bound goes to zero: every coefficient is shrunk to nothing. At \(\lambda = 0\) with \(\mathbf X^T\mathbf X\) invertible, \(\hat\beta_0 = (\mathbf X^T\mathbf X)^{-1}\mathbf X^T\mathbf y\), least squares: \(\lambda\) interpolates between the least-squares fit and the zero fit, and each eigen-direction is shrunk by its own factor \(\mu_j/(\mu_j + \lambda)\), the flat ones most.

Common errors: a factor of \(2\) left on one term and not the other in (a) (the stationary point is then wrong by a factor); claiming invertibility from “\(\lambda > 0\)” without the quadratic form (the argument is that the form is positive on every nonzero \(v\)); in (d), bounding each coefficient by \(|c_j|/\lambda\) and stopping, which is fine, but then forgetting that \(\sum c_j^2 = \|\mathbf X^T\mathbf y\|^2\) needs the basis to be orthonormal.

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?

  1. Given the training data, each \(f_k\) is a fixed function. Take the expectation of one term of \(S_k\): what is the distribution of a selection row, given the training data?
  1. \(\min_k S_k \le S_j\) for every \(j\). Take conditional expectations on both sides, then choose \(j\).
  1. The minimum is \(r + 1\) only when every score is \(r + 1\). With what probability?

  1. Condition on the training data: \(f_k\) is then one fixed function. \(S_k\) is an average of \(m\) terms \(\ell(y_j, f_k(x_j))\), and each selection row \((x_j, y_j)\) is distributed as the population and independent of the training rows, so \(\mathbb E[\ell(y_j, f_k(x_j)) \mid \text{tr}] = \mathbb E\,\ell(Y, f_k(X)) = R(f_k)\) for each \(j\), and the average of \(m\) copies of \(R(f_k)\) is \(R(f_k)\). (If you have proved exercise 1(a), it is fine to say that it applies here, conditionally on the training data.)

  2. For every \(j\), \(S_{\hat k} = \min_k S_k \le S_j\), so \(\mathbb E[S_{\hat k} \mid \text{tr}] \le \mathbb E[S_j \mid \text{tr}] = R(f_j)\) by (a); as this holds for every \(j\), \(\mathbb E[S_{\hat k} \mid \text{tr}] \le \min_j R(f_j)\). And \(R(f_{\hat k}) \ge \min_j R(f_j)\) on every selection set, so \(\min_j R(f_j) \le \mathbb E[R(f_{\hat k}) \mid \text{tr}]\). It is exercise 1 with the finite class \(\{f_1, \ldots, f_K\}\) in place of \(\mathcal F\) and the selection score in place of \(\hat R_n\): choosing on a dataset is minimizing an empirical risk over a class, and the chosen score underestimates the chosen model’s risk for the same reason the training error underestimates the fitted model’s.

  3. The minimum of the two scores is \(r + 1\) only if both are, with probability \(\tfrac14\); otherwise it is \(r - 1\). So \(\mathbb E[S_{\hat k}] = \tfrac34(r - 1) + \tfrac14(r + 1) = r - \tfrac12\), while \(R(f_{\hat k}) = r\) whichever candidate wins: the leaderboard’s top score is half a unit too good on average. With \(K\) candidates, the minimum is \(r + 1\) with probability \(2^{-K}\), so \(\mathbb E[S_{\hat k}] = (r + 1) - 2(1 - 2^{-K}) = r - 1 + 2^{1 - K}\), which falls toward \(r - 1\) as \(K\) grows: in this setting, the more candidates are compared on the same rows, the more the winner’s score understates its risk, even though no candidate is any better. (With unequal risks or correlated scores the size of the gap is a different calculation; the direction is the same.)

  4. Condition on the training and selection data together: \(f_{\hat k}\) is then one fixed function, and the test rows are distributed as the population and independent of both, so each term of the test score has conditional expectation \(R(f_{\hat k})\) and the test score is unbiased for \(R(f_{\hat k})\). Nothing was chosen on the test rows; that is the rule, and it is why the vault is opened once.

Common errors: in (b), applying (a) to \(S_{\hat k}\) as if \(\hat k\) were fixed; in (c), forgetting that \(R(f_{\hat k}) = r\) exactly (both candidates have risk \(r\)), so the whole gap is in the score, not in the model.

Problem 2 (where the proof breaks). The optimism theorem does not cover every method. Its proof needs \(\hat f\) to be the minimizer of the training error over a class fixed before the data, and that fails for ridge regression at a fixed \(\lambda\), for a tree grown greedily, and for any knob chosen on the data (problem 1). The first two cases:

  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 the proof 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?

  1. \(R(\beta x) = (1 - \beta)^2\,\mathbb E X^2\). For the sample, minimize \(\sum_i (x_i - \beta x_i)^2 + \lambda\beta^2\) over the one number \(\beta\).
  1. The proof compares \(\hat R_n(\hat f)\) with \(\hat R_n(f_{\mathcal F})\). What does that comparison need \(\hat f\) to be?

  1. \(R(\beta x) = \mathbb E(X - \beta X)^2 = (1 - \beta)^2\,\mathbb E X^2\), minimized at \(\beta = 1\): \(f_{\mathcal F}(x) = x\) and \(R(f_{\mathcal F}) = 0\). For the sample, \(\tfrac{d}{d\beta}\big[\sum_i (1 - \beta)^2 x_i^2 + \lambda\beta^2\big] = 0\) gives \(\hat\beta_\lambda = \sum_i x_i^2 / (\sum_i x_i^2 + \lambda)\), strictly between \(0\) and \(1\). Its training error is \(\hat R_n(\hat\beta_\lambda) = (1 - \hat\beta_\lambda)^2 \tfrac1n\sum_i x_i^2 = \big(\tfrac{\lambda}{\sum_i x_i^2 + \lambda}\big)^2 \tfrac1n\sum_i x_i^2 > 0\) on every sample, so its expectation is positive while \(R(f_{\mathcal F}) = 0\). The hypothesis it violates: \(\hat f\) must minimize the training error \(\hat R_n\) over \(\mathcal F\), and \(\hat\beta_\lambda\) minimizes a different objective, the training error plus \(\lambda\beta^2\). The minimizer of \(\hat R_n\) over \(\mathcal F\) is \(\beta = 1\), with training error \(0\). (The right inequality still holds: \(\hat\beta_\lambda x \in \mathcal F\), so \(R(\hat\beta_\lambda) \ge 0\).)

  2. Yes: \(\mathcal F_t\) is a class fixed before the data, and \(\hat f\) minimizes \(\hat R_n\) over it, which is all the proof used. For \(t \ge 1\) the unconstrained minimizer \(\beta = 1\) is allowed, so \(\hat f(x) = x\) with training error \(0\); for \(t < 1\) the training error \((1 - \beta)^2 \tfrac1n\sum x_i^2\) decreases in \(\beta\) on \([-t, t]\), so \(\hat f(x) = tx\). Either way the theorem’s sandwich holds, with \(f_{\mathcal F_t}\) the same function (the population minimizer over \(\mathcal F_t\) is also \(\beta = \min(1, t)\)). The penalized and constrained forms trace out the same family of fits as \(\lambda\) and \(t\) vary, but the theorem is about what is minimized: the constrained fit minimizes the training error over \(\mathcal F_t\); the penalized fit minimizes the training error plus a penalty, which is not the training error.

  3. The proof needs \(\hat R_n(\hat f) \le \hat R_n(f)\) for every \(f \in \mathcal F\), which is what “\(\hat f\) minimizes \(\hat R_n\) over \(\mathcal F\)” says. The greedy tree chooses the first split without looking ahead, so another depth-\(2\) tree can have a smaller training error, and the comparison with \(f_{\mathcal F}\) is not available. The theorem is then silent: it does not say the greedy tree’s training error is above its risk, and in practice each split is chosen to reduce the training error, so the training error usually sits below the test error anyway, as it did on the lecture’s leaderboards. No guarantee is not the same as the opposite guarantee.

  4. The right one, \(R(f_{\mathcal F}) \le R(\hat f)\): it needs only \(\hat f \in \mathcal F\) and the definition of \(f_{\mathcal F}\). Whatever a method returns, its risk is at least the best risk its class allows; what the method and the sample decide is how far above.

Common errors: in (a), taking the expectation before computing the training error (the inequality is strict on every sample, so no expectation is needed to see it fail); in (c), concluding from “the theorem does not apply” that the opposite inequality holds.

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 = \|\mathbf y - \mathbf X\beta\|^2 + \lambda\|\beta\|^2\), the ridge objective.

  2. Write the least-squares normal equations for the regression of \(\tilde{\mathbf y}\) on \(\tilde{\mathbf X}\) and check that they are ridge regression’s, \((\mathbf X^T\mathbf X + \lambda\mathbf I)\beta = \mathbf X^T\mathbf y\). 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()?

  1. A stacked vector’s squared norm is the sum of the two blocks’ squared norms. (b) \(\tilde{\mathbf X}^T\tilde{\mathbf X}\) is a sum of two products, one per block.

  1. \(\tilde{\mathbf y} - \tilde{\mathbf X}\beta = \begin{pmatrix}\mathbf y - \mathbf X\beta \\ -\sqrt\lambda\,\beta\end{pmatrix}\), whose squared norm is \(\|\mathbf y - \mathbf X\beta\|^2 + \lambda\|\beta\|^2\).

  2. \(\tilde{\mathbf X}^T\tilde{\mathbf X} = \mathbf X^T\mathbf X + \lambda\mathbf I\) and \(\tilde{\mathbf X}^T\tilde{\mathbf y} = \mathbf X^T\mathbf y\), so the normal equations \(\tilde{\mathbf X}^T\tilde{\mathbf X}\beta = \tilde{\mathbf X}^T\tilde{\mathbf y}\) are exactly ridge regression’s. If \(\tilde{\mathbf X}v = \mathbf 0\) then the lower block gives \(\sqrt\lambda\,v = \mathbf 0\), so \(v = \mathbf 0\): full column rank, with no condition on \(n\). The ridge estimator is the least-squares estimator on a dataset with \(n + p\) rows for \(p\) columns, so always more rows than columns, which is another way to see why \(\mathbf X^T\mathbf X + \lambda\mathbf I\) is invertible.

  3. The \(j\)th extra row is an observation with input \(\mathbf e_j\) (one in position \(j\), zeros elsewhere) and outcome \(0\): it says “\(\beta_j\) should be \(0\)”, and \(\sqrt\lambda\) is how loudly it says it, since its squared residual is \(\lambda\beta_j^2\). With lm(): bind the \(p\) rows of \(\sqrt\lambda\,\mathbf I_p\) under X and \(p\) zeros under y, and fit without an intercept; the optional item at the end does it.

Common errors: forgetting the square root (the extra rows then contribute \(\lambda^2\|\beta\|^2\)); keeping an intercept in the augmented fit, which penalizes it through the zero outcomes.

Coding

1. The task and the split

The lecture’s task: predict a track’s popularity from ten audio features and its release year, squared loss (data note). Genre is left out today, so the numbers are those of the lecture’s rows without genre; an optional item at the end adds it.

feat <- c("danceability", "energy", "loudness", "speechiness", "acousticness",
          "instrumentalness", "liveness", "valence", "tempo", "duration_min", "release_year")
sp <- read_csv("https://ml4ds.com/files/spotify_tracks.csv", show_col_types = FALSE) |>
  distinct(track_id, .keep_all = TRUE) |>          # one row per track
  mutate(duration_min = duration_ms / 60000) |>
  select(track_popularity, all_of(feat))            # the outcome and the eleven inputs

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

set.seed(310)   # change this number: a different split, and different numbers below
id <- sample(nrow(sp))          # a random order of the 10,000 tracks
# sp_train <- sp[id[____], ]    # 5,000 tracks
# sp_tune <- ____               # 2,000
# sp_test <- ____               # 3,000
set.seed(310)   # change this number: a different split, and different numbers below
id <- sample(nrow(sp))
sp_train <- sp[id[1:5000], ]
sp_tune <- sp[id[5001:7000], ]
sp_test <- sp[id[7001:10000], ]

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

stopifnot(nrow(sp_train) == 5000, nrow(sp_tune) == 2000, nrow(sp_test) == 3000,
          length(unique(c(id[1:5000], id[5001:7000], id[7001:10000]))) == 10000)

Question 1. The file has 11,554 rows for 10,000 tracks, because 1,101 tracks are in more than one playlist, and distinct() kept one row per track before the split. 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?

A track in two playlists could land in the training set as one row and in the test set as the other, with the same audio features, the same year and the same popularity. A method that remembers its training rows, \(k\)-NN with small \(k\) or a forest, predicts that test row almost exactly, and the test score is too good. The rule: a row is fair game for scoring only if nothing about the model was chosen on it, and the model has seen this track. The unit of the split has to be the unit that repeats.

Step 2. The baselines on the tuning tracks, as in the lecture. score() builds one leaderboard row: a competitor’s MSE on the tuning tracks and its standard error. The mean is done; add the plane (least squares on every input) and \(k\) nearest neighbors on the standardized inputs, with \(k\) chosen on the tuning tracks from ks. scale() centers and scales each column; the tuning rows get the training means and standard deviations, so that both matrices are on the same scale.

y_train <- sp_train$track_popularity
y_tune <- sp_tune$track_popularity
score <- function(name, pred)
  tibble(competitor = name, tune_mse = mse(y_tune, pred), se = se_loss(y_tune, pred))
mu <- colMeans(sp_train[, feat])
s <- apply(sp_train[, feat], 2, sd)
X_train <- scale(sp_train[, feat], mu, s)
X_tune <- scale(sp_tune[, feat], mu, s)
ks <- c(1, 5, 10, 20, 50, 100, 200)
# fit_plane <- lm(____)
# knn_tune <- sapply(ks, function(k) ____)      # the tuning MSE at each k
# k_best <- ____
lb <- bind_rows(
  score("the mean", rep(mean(y_train), nrow(sp_tune)))
  # , score("least squares", ____)
  # , score("k-NN", ____)
)
lb
fit_plane <- lm(track_popularity ~ ., sp_train)
knn_tune <- sapply(ks, function(k) mse(y_tune, knn_predict(X_tune, X_train, y_train, k)))
k_best <- ks[which.min(knn_tune)]
rbind(k = ks, tune_mse = round(knn_tune))
         [,1] [,2] [,3] [,4] [,5] [,6] [,7]
k           1    5   10   20   50  100  200
tune_mse 1002  581  556  531  525  523  523
lb <- bind_rows(
  score("the mean", rep(mean(y_train), nrow(sp_tune))),
  score("least squares", predict(fit_plane, sp_tune)),
  score("k-NN", knn_predict(X_tune, X_train, y_train, k_best))
)
lb
# A tibble: 3 x 3
  competitor    tune_mse    se
  <chr>            <dbl> <dbl>
1 the mean          564.  12.7
2 least squares     526.  12.3
3 k-NN              523.  12.4

The plane and \(k\)-NN at its best \(k\) (100) are within a standard error or two of each other, both well below the mean.

2. Five structures, each tuned on the tuning 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 tuning 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; predict() on the tuning matrix then gives one column of predictions per \(\lambda\). Choose \(\lambda\), and count the coefficients that are exactly zero at it.

lasso_fit <- glmnet(X_train, y_train, standardize = FALSE)   # the columns already are
# lasso_tune <- apply(predict(lasso_fit, X_tune), 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", ____))
lasso_fit <- glmnet(X_train, y_train, standardize = FALSE)
lasso_tune <- apply(predict(lasso_fit, X_tune), 2, function(pred) mse(y_tune, pred))
lam_best <- lasso_fit$lambda[which.min(lasso_tune)]
b_lasso <- coef(lasso_fit, s = lam_best)
sum(b_lasso[-1] == 0)                       # zero coefficients (the intercept excluded)
[1] 0
rownames(b_lasso)[-1][b_lasso[-1] == 0]     # which ones
character(0)
lb <- bind_rows(lb, score("lasso", predict(lasso_fit, X_tune, s = lam_best)))

The chosen \(\lambda\) is 0.0095 on glmnet’s scale, near the small end of the path: with 5,000 tracks and 11 columns the data can afford every coefficient, 0 of the 11 are set to zero, and the lasso’s row on the leaderboard sits next to least squares’.

Step 4. Additive splines: lm() on a natural-spline basis of size \(J\) for each input, all in one least-squares fit, as in the lecture. spline_fit(J, data) writes the formula out. Choose \(J\) from Js.

spline_fit <- function(J, data)
  lm(track_popularity ~ ns(danceability, J) + ns(energy, J) + ns(loudness, J) + ns(speechiness, J) +
       ns(acousticness, J) + ns(instrumentalness, J) + ns(liveness, J) + ns(valence, J) +
       ns(tempo, J) + ns(duration_min, J) + ns(release_year, J), data)
Js <- c(2, 3, 4, 6)
# spline_tune <- sapply(Js, function(J) ____)
# J_best <- ____
# fit_spline <- spline_fit(J_best, sp_train)
# lb <- bind_rows(lb, score("additive splines", ____))
spline_fit <- function(J, data)
  lm(track_popularity ~ ns(danceability, J) + ns(energy, J) + ns(loudness, J) + ns(speechiness, J) +
       ns(acousticness, J) + ns(instrumentalness, J) + ns(liveness, J) + ns(valence, J) +
       ns(tempo, J) + ns(duration_min, J) + ns(release_year, J), data)
Js <- c(2, 3, 4, 6)
spline_tune <- sapply(Js, function(J) mse(y_tune, predict(spline_fit(J, sp_train), sp_tune)))
rbind(J = Js, tune_mse = round(spline_tune))
         [,1] [,2] [,3] [,4]
J           2    3    4    6
tune_mse  510  495  496  497
J_best <- Js[which.min(spline_tune)]
fit_spline <- spline_fit(J_best, sp_train)
lb <- bind_rows(lb, score("additive splines", predict(fit_spline, sp_tune)))

Step 5. A tree and a forest. tree_fit(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 (the slowest fit of the seminar, about half a minute).

tree_fit <- function(d, data)
  rpart(track_popularity ~ ., data, control = rpart.control(maxdepth = d, cp = 0, minbucket = 10))
depths <- c(2, 3, 4, 6, 8)
# tree_tune <- sapply(depths, function(d) ____)
# d_best <- ____
# fit_tree <- tree_fit(d_best, sp_train)
fit_rf <- randomForest(track_popularity ~ ., sp_train, ntree = 300)
# lb <- bind_rows(lb, score("tree", ____), score("forest", ____))
tree_fit <- function(d, data)
  rpart(track_popularity ~ ., data, control = rpart.control(maxdepth = d, cp = 0, minbucket = 10))
depths <- c(2, 3, 4, 6, 8)
tree_tune <- sapply(depths, function(d) mse(y_tune, predict(tree_fit(d, sp_train), sp_tune)))
rbind(depth = depths, tune_mse = round(tree_tune))
         [,1] [,2] [,3] [,4] [,5]
depth       2    3    4    6    8
tune_mse  518  513  520  533  558
d_best <- depths[which.min(tree_tune)]
fit_tree <- tree_fit(d_best, sp_train)
fit_rf <- randomForest(track_popularity ~ ., sp_train, ntree = 300)
lb <- bind_rows(lb, score("tree", predict(fit_tree, sp_tune)), score("forest", predict(fit_rf, sp_tune)))

Step 6. A network with one hidden layer of \(M\) units: nnet() on the standardized matrix, linout = TRUE for a numeric outcome, decay a small ridge penalty on the weights (fixed today), and maxit = 300 a fixed budget of iterations rather than a convergence test. The weights start at random values, so two fits of the same network differ slightly. Choose \(M\) from Ms.

net_fit <- function(M, X, y) nnet(X, y, size = M, decay = 0.1, linout = TRUE, maxit = 300, trace = FALSE)
Ms <- c(2, 4, 8, 16)
# net_tune <- sapply(Ms, function(M) ____)
# M_best <- ____
# fit_net <- net_fit(M_best, X_train, y_train)
# lb <- bind_rows(lb, score("network", ____))
net_fit <- function(M, X, y) nnet(X, y, size = M, decay = 0.1, linout = TRUE, maxit = 300, trace = FALSE)
Ms <- c(2, 4, 8, 16)
net_tune <- sapply(Ms, function(M) mse(y_tune, predict(net_fit(M, X_train, y_train), X_tune)))
rbind(M = Ms, tune_mse = round(net_tune))
         [,1] [,2] [,3] [,4]
M           2    4    8   16
tune_mse  496  508  512  531
M_best <- Ms[which.min(net_tune)]
fit_net <- net_fit(M_best, X_train, y_train)
lb <- bind_rows(lb, score("network", predict(fit_net, X_tune)))

The tuning leaderboard, with each competitor’s chosen knob:

lb |> arrange(tune_mse) |> mutate(across(c(tune_mse, se), round))
# A tibble: 8 x 3
  competitor       tune_mse    se
  <chr>               <dbl> <dbl>
1 additive splines      495    12
2 network               496    12
3 forest                499    12
4 tree                  513    12
5 k-NN                  523    12
6 least squares         526    12
7 lasso                 526    12
8 the mean              564    13
ggplot(lb, aes(tune_mse, fct_reorder(competitor, -tune_mse))) +
  geom_errorbar(aes(xmin = tune_mse - se, xmax = tune_mse + se), width = 0.3, colour = "grey60") +
  geom_point(size = 3) +
  labs(x = "MSE on the 2,000 tuning tracks", y = NULL, caption = "bars: one standard error")

c(k = k_best, lambda = signif(lam_best, 2), J = J_best, depth = d_best, M = M_best)
      k  lambda       J   depth       M 
1.0e+02 9.5e-03 3.0e+00 3.0e+00 2.0e+00

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

The three that can bend (the additive splines, the forest and the network) lead, about 30 below the plane, which is the structure’s worth on these data: a nonlinear year effect and loudness that a plane cannot follow. Among the three the gaps are a few units against bars of about 12, so this picture does not settle their order; the bars are each score’s own standard error, and the uncertainty in a difference of two scores is a different number (the paired comparison in part 3). The mean’s row and the plane’s are unbiased for their risks: both were fixed before the tuning tracks were looked at, and so is the forest’s, which has no knob. Every other row chose its knob on these same tracks and reports the minimum of several noisy scores, so each is biased downward (problem 1). The vault, next, scores every competitor on tracks nothing was chosen on.

3. The vault

Step 7. Refit every competitor on the training and tuning tracks together, each at the knob it chose, and predict the test tracks. The standardized matrices are rebuilt with the 7,000 rows’ means and standard deviations. \(k\)-NN has nothing to refit: its “model” is the rows it averages over, now 7,000 of them. Refit only the competitors you fitted in part 2; the solution skips the others.

sp_both <- bind_rows(sp_train, sp_tune)
y_both <- sp_both$track_popularity
y_test <- sp_test$track_popularity
mu_both <- colMeans(sp_both[, feat])
s_both <- apply(sp_both[, feat], 2, sd)
X_both <- scale(sp_both[, feat], mu_both, s_both)
X_test <- scale(sp_test[, feat], mu_both, s_both)
pred_test <- list(
  "the mean" = rep(mean(y_both), nrow(sp_test))
  # , "least squares" = ____
  # , "k-NN" = ____
  # , "lasso" = ____
  # , "additive splines" = ____
  # , "tree" = ____
  # , "forest" = ____
  # , "network" = ____
)
sp_both <- bind_rows(sp_train, sp_tune)
y_both <- sp_both$track_popularity
y_test <- sp_test$track_popularity
mu_both <- colMeans(sp_both[, feat])
s_both <- apply(sp_both[, feat], 2, sd)
X_both <- scale(sp_both[, feat], mu_both, s_both)
X_test <- scale(sp_test[, feat], mu_both, s_both)
fit_plane_both <- lm(track_popularity ~ ., sp_both)
lasso_both <- glmnet(X_both, y_both, lambda = lam_best, standardize = FALSE)
fit_spline_both <- spline_fit(J_best, sp_both)
pred_test <- list(
  "the mean" = rep(mean(y_both), nrow(sp_test)),
  "least squares" = predict(fit_plane_both, sp_test),
  "k-NN" = knn_predict(X_test, X_both, y_both, k_best),
  "lasso" = as.vector(predict(lasso_both, X_test)),
  "additive splines" = predict(fit_spline_both, sp_test)
)
# the competitors of step 5 and step 6, if they were fitted
if (exists("d_best")) pred_test$tree <- predict(tree_fit(d_best, sp_both), sp_test)
if (exists("fit_rf")) { fit_rf_both <- randomForest(track_popularity ~ ., sp_both, ntree = 300); pred_test$forest <- predict(fit_rf_both, sp_test) }
if (exists("M_best")) pred_test$network <- as.vector(predict(net_fit(M_best, X_both, y_both), X_test))

Checkpoint:

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

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

# vault <- imap_dfr(pred_test, ~ tibble(competitor = .y, test_mse = ____, se = ____)) |> arrange(test_mse)
# leaders <- ____                     # the top two of the tuning leaderboard, not of the vault
# d <- ____                           # their squared errors, track by track
# c(mean = mean(d), se = ____)
vault <- imap_dfr(pred_test, ~ tibble(competitor = .y, test_mse = mse(y_test, .x), se = se_loss(y_test, .x))) |>
  arrange(test_mse)
vault |> mutate(across(c(test_mse, se), round))
# A tibble: 8 x 3
  competitor       test_mse    se
  <chr>               <dbl> <dbl>
1 forest                499    10
2 network               499    10
3 additive splines      499    10
4 tree                  516    10
5 k-NN                  519    10
6 lasso                 534    10
7 least squares         534    10
8 the mean              561    10
leaders <- lb$competitor[order(lb$tune_mse)][1:2]
leaders
[1] "additive splines" "network"         
d <- (y_test - pred_test[[leaders[1]]])^2 - (y_test - pred_test[[leaders[2]]])^2
c(mean = mean(d), se = sd(d) / sqrt(length(d)))
     mean        se 
0.2503726 3.1953413

The tuning leaders were additive splines and network; in the vault their test scores differ by 0.3 with a paired standard error of 3.2. A paired standard error is the right one for a difference: both competitors are scored on the same tracks, and a hard track is hard for both, so most of each score’s own spread cancels in the difference. The pair was chosen before the vault was opened; choosing it from the vault’s own top two and then testing the difference would be choosing on the test set.

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

vault |> left_join(select(lb, competitor, tune_mse), by = "competitor") |>
  mutate(across(c(test_mse, se, tune_mse), round))
# A tibble: 8 x 4
  competitor       test_mse    se tune_mse
  <chr>               <dbl> <dbl>    <dbl>
1 forest                499    10      499
2 network               499    10      496
3 additive splines      499    10      495
4 tree                  516    10      513
5 k-NN                  519    10      523
6 lasso                 534    10      526
7 least squares         534    10      526
8 the mean              561    10      564

The three that bend stay the three leaders, and the vault does not separate them; the tree is alone in the middle; the plane, the lasso and \(k\)-NN stay a group behind it; the mean is last. The test scores sit at a different level from the tuning scores, the mean’s included (561 against 564), and the mean chose nothing: the two sets of tracks are different draws, and every refit used 7,000 rows instead of 5,000. So this comparison cannot isolate problem 1’s bias: the models and the evaluation samples both changed, and the bias is a statement about averages over draws of the tuning set; seeing it would take many splits. What the vault settles is the structure’s worth: the three that bend sit about 30 below the plane on both leaderboards (35 in the vault for its leader), scored on tracks nothing was chosen on, with standard errors that mean what they say.

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 = ____)
round(c(forest_train = mse(y_both, predict(fit_rf_both, sp_both)), forest_test = mse(y_test, predict(fit_rf_both, sp_test)),
        plane_train = mse(y_both, fitted(fit_plane_both)), plane_test = vault$test_mse[vault$competitor == "least squares"]))
forest_train  forest_test  plane_train   plane_test 
          95          499          527          534

The plane is covered: least squares is the empirical risk minimizer over a class fixed in advance, and its training error sits a little below its test error, the gap the theorem describes. The forest is not covered (each tree is grown greedily on a resample, and the ensemble is not a minimizer over any class), and its training error is about a fifth of its test error: it remembers the training tracks almost individually. Nothing in the theorem forbids a larger gap for a method it does not cover; it just says nothing about it. Either way the training error is the wrong number to compare models on, and the lecture’s hollow points were the same picture.

Optional, can be done at home: ridge two ways

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

lambda <- 100
X_c <- scale(X_train, scale = FALSE)        # centered columns, as in exercise 2
y_c <- y_train - mean(y_train)
# 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.

lambda <- 100
X_c <- scale(X_train, scale = FALSE)
y_c <- y_train - mean(y_train)
b_formula <- solve(t(X_c) %*% X_c + lambda * diag(ncol(X_c)), t(X_c) %*% y_c)
X_aug <- rbind(X_c, sqrt(lambda) * diag(ncol(X_c))); y_aug <- c(y_c, rep(0, ncol(X_c)))
b_lm <- coef(lm(y_aug ~ X_aug - 1))
max(abs(as.vector(b_formula) - b_lm))
[1] 7.993606e-14

The two agree to rounding. (glmnet(alpha = 0) fits ridge too, with \(\lambda\) on its own scale and the outcome rescaled internally; matching its numbers to these is fiddly and not the point.)

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: the path starts at the end where every coefficient is zero (the largest \(\lambda\); which side of the plot that is depends on your version of glmnet) and the coefficients leave zero one by one as \(\lambda\) falls. Where are the release year and loudness, the two inputs the lecture’s component plots bent most?

plot(lasso_fit, xvar = "lambda", label = TRUE)

first_in <- apply(as.matrix(coef(lasso_fit))[-1, ], 1, function(b) which(b != 0)[1])   # the step at which each enters
sort(first_in)
          energy     acousticness instrumentalness     duration_min 
               2                2                2                2 
    danceability         loudness         liveness            tempo 
               9               10               21               23 
     speechiness          valence     release_year 
              28               31               49

Energy and duration are among the first to enter, as in the lecture’s figure. Loudness enters at step 10 and the release year not until step 49 of 63, with a coefficient of -0.06 at the end of the path, although the year is the input the additive fit bends most. When an input enters depends on its linear association with popularity, and with the other inputs already in; the year’s effect is curved, so a line through it is nearly flat. Entering the path early is not the same as mattering.

Optional, can be done at home: add genre

The lecture’s leaderboard also used the playlist’s genre, a factor with six levels. Add it to the training and tuning tracks (playlist_genre in the file, matched by row after distinct()), refit the plane and the additive splines with + genre, and score them on the tuning tracks. Which methods take a factor as it is, and which need it as 0/1 columns (model.matrix()) first?

genre <- read_csv("https://ml4ds.com/files/spotify_tracks.csv", show_col_types = FALSE) |>
  distinct(track_id, .keep_all = TRUE) |> pull(playlist_genre) |> factor()
sp_train$genre <- genre[id[1:5000]]
sp_tune$genre <- genre[id[5001:7000]]
fit_plane_g <- lm(track_popularity ~ ., sp_train)
fit_spline_g <- lm(update(formula(fit_spline), . ~ . + genre), sp_train)
round(c(plane = mse(y_tune, predict(fit_plane_g, sp_tune)), splines = mse(y_tune, predict(fit_spline_g, sp_tune))))
  plane splines 
    515     486

lm(), rpart() and randomForest() take a factor in a formula as it is (lm builds the five 0/1 columns itself); glmnet(), nnet() and knn_predict() work on a numeric matrix, so the factor goes through model.matrix() first, and then the question is whether to standardize the 0/1 columns like the others. Genre takes about 10 off the plane’s score and about the same off the splines’.