Week 5

Throughout, \(\bfX\in\Rbb^{n\times d}\) has rank \(r\geq1\) and compact SVD \(\bfX=\sum_{i=1}^r\sigma_i\bfu_i\bfv_i^\top\), and \(\bfy\in\Rbb^n\). Let \(\lambda>0\), and let \(\widehat{\bftheta}_\lambda\) denote the ridge regression estimator \[ \widehat{\bftheta}_\lambda =\argmin_{\bftheta\in\Rbb^d}\Big\{\|\bfy-\bfX\bftheta\|_2^2+\lambda\|\bftheta\|_2^2\Big\}. \]

Problem 1

The ridge solution \(\widehat{\bftheta}_\lambda=(\bfX^\top\bfX+\lambda\bfI_d)^{-1}\bfX^\top\bfy\) requires solving a \(d\times d\) linear system. When there are many more features than observations (\(d\gg n\)), there is a better way! In this problem you will find a second formula that only requires an \(n\times n\) system.

Show that \(\bfX\bfX^\top+\lambda\bfI_n\) is invertible.
Show that \((\bfX^\top\bfX+\lambda\bfI_d)\bfX^\top=\bfX^\top(\bfX\bfX^\top+\lambda\bfI_n)\).
Conclude that \[ \widehat{\bftheta}_\lambda=\bfX^\top(\bfX\bfX^\top+\lambda\bfI_n)^{-1}\bfy. \]
A study measures the expression of \(d=10{,}000\) genes in each of \(n=100\) patients. Solving a \(k\times k\) linear system takes roughly \(k^3\) arithmetic operations. Roughly how many times faster is the formula of part 3 than the one from lecture?

Solution sketch

Since \(\bfX\bfX^\top+\lambda\bfI_n\) is square, it suffices to show its null space is trivial. Suppose \((\bfX\bfX^\top+\lambda\bfI_n)\bfz=\mathbf{0}\). Multiplying on the left by \(\bfz^\top\) gives \[ 0=\bfz^\top\bfX\bfX^\top\bfz+\lambda\bfz^\top\bfz=\|\bfX^\top\bfz\|_2^2+\lambda\|\bfz\|_2^2\geq\lambda\|\bfz\|_2^2. \] Since \(\lambda>0\), this forces \(\|\bfz\|_2=0\), so \(\bfz=\mathbf{0}\).
Both sides equal \(\bfX^\top\bfX\bfX^\top+\lambda\bfX^\top\).
Multiply the identity of part 2 on the left by \((\bfX^\top\bfX+\lambda\bfI_d)^{-1}\) and on the right by \((\bfX\bfX^\top+\lambda\bfI_n)^{-1}\), both of which exist by part 1 and the fact that \(\bfX^\top\bfX+\lambda\bfI_d\) is invertible for every \(\lambda>0\): \[ \bfX^\top(\bfX\bfX^\top+\lambda\bfI_n)^{-1}=(\bfX^\top\bfX+\lambda\bfI_d)^{-1}\bfX^\top. \] Apply both sides to \(\bfy\) and use the ridge solution.
The lecture formula solves a \(10{,}000\times10{,}000\) system, about \((10^4)^3=10^{12}\) operations; part 3 solves a \(100\times100\) system, about \((10^2)^3=10^6\) operations. That is \((d/n)^3=100^3=10^6\) times fewer. (Forming \(\bfX\bfX^\top\) costs about \(n^2d=10^8\) operations, still far less than \(d^3\).)

Problem 2

It is tempting to choose \(\lambda\) by making the training error \(\|\bfy-\bfX\widehat{\bftheta}_\lambda\|_2^2\) as small as possible. This problem is a guided derivation of training error always prefers \(\lambda=0\), which shows that this never works. Let \(\bfP=\sum_{i=1}^r\bfu_i\bfu_i^\top=\bfU_r\bfU_r^\top\) be the orthogonal projection onto \(\Ccal(\bfX)\).

Use ridge regression filters the singular values to show that the fitted values are \[ \bfX\widehat{\bftheta}_\lambda=\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\,\bfu_i. \]
Show that \[ \bfy-\bfX\widehat{\bftheta}_\lambda =\sum_{i=1}^r\frac{\lambda}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\,\bfu_i+(\bfy-\bfP\bfy). \]
Conclude that \[ \|\bfy-\bfX\widehat{\bftheta}_\lambda\|_2^2 =\sum_{i=1}^r\Big(\frac{\lambda}{\sigma_i^2+\lambda}\Big)^2(\bfu_i^\top\bfy)^2+\|\bfy-\bfP\bfy\|_2^2, \] and that the training error never decreases as \(\lambda\) grows. Hint: What is \(\bfu_i^\top(\bfy-\bfP\bfy)\)?
Find the limits of the training error as \(\lambda\to0^+\) and as \(\lambda\to\infty\), and interpret each.

Solution sketch

Apply \(\bfX\) to the formula for \(\widehat{\bftheta}_\lambda\) and use \(\bfX\bfv_i=\sigma_i\bfu_i\): \[ \bfX\widehat{\bftheta}_\lambda =\sum_{i=1}^r\frac{\sigma_i}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\,\sigma_i\bfu_i =\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\,\bfu_i. \]
Write \(\bfy=\bfP\bfy+(\bfy-\bfP\bfy)=\sum_{i=1}^r(\bfu_i^\top\bfy)\bfu_i+(\bfy-\bfP\bfy)\) and subtract part 1, using \(1-\frac{\sigma_i^2}{\sigma_i^2+\lambda}=\frac{\lambda}{\sigma_i^2+\lambda}\).
The vectors \(\bfu_1,\ldots,\bfu_r\) are orthonormal, and \(\bfy-\bfP\bfy\perp\Ccal(\bfX)\), so \(\bfu_i^\top(\bfy-\bfP\bfy)=0\) for every \(i\). So the \(r+1\) terms in part 2 are mutually orthogonal, and the Pythagorean theorem gives the formula. Each \(\lambda/(\sigma_i^2+\lambda)=1-\sigma_i^2/(\sigma_i^2+\lambda)\) is positive and increases with \(\lambda\), so each term never decreases, and neither does the sum.
As \(\lambda\to0^+\), every factor \(\lambda/(\sigma_i^2+\lambda)\to0\), so the training error tends to \(\|\bfy-\bfP\bfy\|_2^2\), the least squares residual. As \(\lambda\to\infty\), every factor tends to \(1\), so the training error tends to \(\sum_{i=1}^r(\bfu_i^\top\bfy)^2+\|\bfy-\bfP\bfy\|_2^2=\|\bfP\bfy\|_2^2+\|\bfy-\bfP\bfy\|_2^2=\|\bfy\|_2^2\), the error of the model \(\bftheta=\mathbf{0}\). Minimizing training error therefore always picks least squares, since least squares is defined as the minimizer of the training error; judging \(\lambda\) honestly requires data not used in the fit.

Problem 3

Ridge regression was defined with a penalty, but it can also be read as least squares with a budget on the size of the coefficients. Let \(\bftheta^\dagger=\sum_{i=1}^r\frac{\bfu_i^\top\bfy}{\sigma_i}\bfv_i\) be the minimum-norm least squares solution.

Use ridge regression filters the singular values to show that \[ \|\widehat{\bftheta}_\lambda\|_2^2=\sum_{i=1}^r\frac{\sigma_i^2}{(\sigma_i^2+\lambda)^2}(\bfu_i^\top\bfy)^2. \] Conclude that \(\|\widehat{\bftheta}_\lambda\|_2\) never increases as \(\lambda\) grows, and that \(\|\widehat{\bftheta}_\lambda\|_2\leq\|\bftheta^\dagger\|_2\).
Let \(t=\|\widehat{\bftheta}_\lambda\|_2\). Show that \(\widehat{\bftheta}_\lambda\) solves the constrained problem \[ \min_{\bftheta\in\Rbb^d}\|\bfy-\bfX\bftheta\|_2^2 \quad\text{subject to}\quad \|\bftheta\|_2\leq t. \] Hint: If \(\|\bftheta\|_2\leq t\), compare the ridge objective at \(\bftheta\) and at \(\widehat{\bftheta}_\lambda\).
Combine parts 1 and 2 with Problem 2 to describe, in a sentence or two, what happens to the budget \(t\) and to the training error as \(\lambda\) increases from \(0\) to \(\infty\).

Solution sketch

The vectors \(\bfv_1,\ldots,\bfv_r\) are orthonormal, so the squared norm of \(\sum_i c_i\bfv_i\) is \(\sum_ic_i^2\); apply this with \(c_i=\frac{\sigma_i}{\sigma_i^2+\lambda}(\bfu_i^\top\bfy)\). Each coefficient \(\sigma_i^2/(\sigma_i^2+\lambda)^2\) decreases as \(\lambda\) grows, so the sum never increases. Since \((\sigma_i^2+\lambda)^2>\sigma_i^4\), each term is at most \((\bfu_i^\top\bfy)^2/\sigma_i^2\), and summing these gives \(\|\bftheta^\dagger\|_2^2\) by the same orthonormality argument.
The point \(\widehat{\bftheta}_\lambda\) is feasible, since \(\|\widehat{\bftheta}_\lambda\|_2=t\). Let \(\bftheta\) be any point with \(\|\bftheta\|_2\leq t\). Since \(\widehat{\bftheta}_\lambda\) minimizes the ridge objective, \[ \|\bfy-\bfX\bftheta\|_2^2+\lambda\|\bftheta\|_2^2 \geq\|\bfy-\bfX\widehat{\bftheta}_\lambda\|_2^2+\lambda t^2. \] Rearranging, \(\|\bfy-\bfX\bftheta\|_2^2\geq\|\bfy-\bfX\widehat{\bftheta}_\lambda\|_2^2+\lambda(t^2-\|\bftheta\|_2^2)\geq\|\bfy-\bfX\widehat{\bftheta}_\lambda\|_2^2\). So no feasible point fits better than \(\widehat{\bftheta}_\lambda\).
As \(\lambda\) increases, the budget \(t=\|\widehat{\bftheta}_\lambda\|_2\) shrinks, from at most \(\|\bftheta^\dagger\|_2\) (approached as \(\lambda\to0^+\)) down toward \(0\) (as \(\lambda\to\infty\)). A smaller budget allows a worse fit, which is why the training error increases with \(\lambda\), as found in Problem 2. Large \(\lambda\) means a tight budget on the coefficients.

Problem 4

For least squares with \(\operatorname{rank}(\bfX)=d\), the fit uses \(d\) parameters. Ridge regression also has \(d\) coefficients, but implicity adds the information that we should be ‘shrinking’ them. In theory, therefore, it should count as fitting “fewer” parameters. Let’s make that precise! \[ \bfH_\lambda=\bfX(\bfX^\top\bfX+\lambda\bfI_d)^{-1}\bfX^\top\in\Rbb^{n\times n}, \] so that the ridge fitted values are \(\widehat{\bfy}=\bfX\widehat{\bftheta}_\lambda=\bfH_\lambda\bfy\).

Show that \[ \bfH_\lambda=\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda}\,\bfu_i\bfu_i^\top. \]
Let’s define the effective degrees of freedom for ridge regression as \(\operatorname{df}(\lambda)=\operatorname{tr}(\bfH_\lambda)\). Show that \[ \operatorname{df}(\lambda)=\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda}. \] Moreover, show that \(\operatorname{df}\) is strictly decreasing in \(\lambda\), and that \(\operatorname{df}(\lambda)\to r\) as \(\lambda\to0^+\) and \(\operatorname{df}(\lambda)\to0\) as \(\lambda\to\infty\). What matrix does \(\bfH_\lambda\) approach as \(\lambda\to0^+\)?
Suppose \(d=3\) and \(\bfX\) has singular values \(10\), \(1\), and \(0.1\). Compute \(\operatorname{df}(1)\). Roughly how many parameters is ridge regression with \(\lambda=1\) “really” fitting, and which directions \(\bfv_i\) does it keep?

Solution sketch

Write \(\bfX^\top\bfX=\sum_{i=1}^d\sigma_i^2\bfv_i\bfv_i^\top\) (with \(\sigma_i=0\) for \(i>r\)), so \((\bfX^\top\bfX+\lambda\bfI_d)^{-1}=\sum_{i=1}^d\frac{1}{\sigma_i^2+\lambda}\bfv_i\bfv_i^\top\). Then, using \(\bfX\bfv_i=\sigma_i\bfu_i\), \[ \bfH_\lambda =\sum_{i=1}^d\frac{1}{\sigma_i^2+\lambda}(\bfX\bfv_i)(\bfX\bfv_i)^\top =\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda}\,\bfu_i\bfu_i^\top. \]
The trace is linear and \(\operatorname{tr}(\bfu_i\bfu_i^\top)=\bfu_i^\top\bfu_i=1\), which gives the formula. Each term \(\sigma_i^2/(\sigma_i^2+\lambda)\) is strictly decreasing in \(\lambda\) because \(\sigma_i>0\), tends to \(1\) as \(\lambda\to0^+\), and tends to \(0\) as \(\lambda\to\infty\); there are \(r\) terms. As \(\lambda\to0^+\), \(\bfH_\lambda\to\sum_{i=1}^r\bfu_i\bfu_i^\top=\bfU_r\bfU_r^\top\), the orthogonal projection onto \(\Ccal(\bfX)\), whose trace is \(r\). So \(\operatorname{df}\) interpolates between the number of parameters least squares fits (namely \(r\), which is \(d\) when \(\operatorname{rank}(\bfX)=d\)) and \(0\).
\(\operatorname{df}(1)=\frac{100}{101}+\frac{1}{2}+\frac{0.01}{1.01}=\frac{100}{101}+\frac12+\frac{1}{101}=1.5\). Ridge keeps \(\bfv_1\) almost entirely (filter factor \(\approx0.99\)), keeps half of \(\bfv_2\), and essentially discards \(\bfv_3\) (filter factor \(\approx0.01\)), so it fits about one and a half parameters rather than three.