# ST443 Machine Learning and data mining
Week 7

Philipp Sterzinger

p.sterzinger@lse.ac.uk

27 November 2025

# Homework 4

## Question 1

> Regularised least squares. Consider the linear regression model:

> \[ y = X \beta + \epsilon \] where \(y = [y_1, \ldots, y_n]^\top \in \Re^n\), \(X = [x_1,\ldots, x_n] \in \Re^{n \times p}\), \(\beta = [\beta_1,\ldots, \beta_p]^\top \in \Re^p\) and \(\epsilon = [\epsilon_1,\ldots, \epsilon_n]^\top \in \Re^n\) is a vector of vector of i.i.d. random variables with zero mean and variance \(\sigma^2\).

1. Ordinary least squares. The ordinary least squares (OLS) parameter estimate \(\hat{\beta}^{\text{OLS}}\) is defined as a minimiser of the following optimisation problem: \[ \min_{\beta \in \Re^p} L(\beta), \quad L(\beta) = \| y - X\beta \|_2^2 \,. \] Show that the OLS estimator is given by \[ \hat{\beta}^{\text{OLS}} = (X^\top X)^{-1}X^\top y \,. \] whenever \(X\) has full column rank (i.e. predictors are not linearly dependent). Does an OLS estimator always exist and when it exists under what condition is it unique?

## Question 1.a

### Show

- If \(\text{rank}(X) = p\), \[\hat{\beta}^{\text{OLS}} = (X^\top X)^{-1}X^\top y\]

### Know

- If \(\text{rank}(X) = p\), then \(X^\top X\) is invertible

- \[ \nabla_x (x^\top a) = a \]

- \[ \nabla_x (x^\top A x) = Ax + A^\top x \]

## Question 1.a

\[ \begin{aligned} L(\beta) &= \| y - X\beta \|_2^2 \\ &= (y - X\beta)^\top (y - X\beta) \\ &= y^\top y - y^\top X\beta - \beta^\top X^\top y - \beta^\top (X^\top X) \beta \\ \end{aligned} \] Differentiating yields: \[ \begin{aligned} \nabla_\beta L(\beta) = -2 X^\top y - (X^\top X) \beta - (X^\top X)^\top \beta &= -2X^\top y - 2 (X^\top X) \beta \overset{!}{=} 0 \\ \iff (X^\top X) \beta &= X^\top y \\ \iff \beta &= (X^\top X)^{-1} X^\top y \end{aligned} \] Checking for optimality: \[ \nabla_\beta \nabla^\top_\beta L(\beta) = - 2 (X^\top X) \preceq 0, \quad \forall \beta \,, \] So \(L(\beta)\) is convex and \(\hat{\beta}^{\text{OLS}}\) is a minimum (and unique).

- Q: Why is \((X^\top X)\) invertible?

## Question 1.a

### Invertibility of \(X^\top X\)

> If \(X^\top X\) has full rank, then it is invertible

## Question 1.a

We know that \(X\) has full column rank, i.e. \[ X v = 0 \implies v = 0 \]

What this says is:

\[ Xv = \begin{bmatrix} {\mid} & {\mid} & \cdots & {\mid} \\ x_1 & x_2 & \cdots & x_p \\ {\mid} & {\mid} & \cdots & {\mid} \\ \end{bmatrix} \begin{bmatrix} v_1 \\ \vdots \\ v_p \end{bmatrix} = v_1\cdot \begin{bmatrix} {\mid} \\ x_1 \\ {\mid} \end{bmatrix} + v_2 \cdot\begin{bmatrix} {\mid} \\ x_2 \\ {\mid} \end{bmatrix} + \cdots + v_p\cdot\begin{bmatrix} {\mid} \\ x_p \\ {\mid} \end{bmatrix} \] is just a linear combination of the columns of \(X\).

- If \(X\) has full column rank, they are linearly independent, so no linear combination of \(p-1\) columns can achieve the other one
- Equivalently, \(Xv = 0\) must imply that \(v = 0\)

## Question 1.a

### Invertibility of \(X^\top X\)

How can we use this find the rank of \(X^\top X\)?

Assume there is a \(v: (X^\top X) v = 0\). Then,

\[ \begin{aligned} v^\top (X^\top X) v = v^\top 0 = 0 \,, \end{aligned} \]

But

\[ v^\top (X^\top X)v = (Xv)^\top (Xv) =\tilde{v}^\top \tilde{v} = 0 \iff \tilde{v} = 0 \,. \]

By construction

\[ \tilde{v} = Xv = 0 \implies v = 0 \,. \]

So \(X^\top X\) must have full rank, therefore it is invertible.

## Question 1.a

### Invertibility of \(X^\top X\)

What happens if \(X\) has linearly dependent columns?

Code

    library(plotly)

    vertices <- matrix(c(
      0,0,0,
      1,0,0,
      1,1,0,
      0,1,0,
      0,0,1,
      1,0,1,
      1,1,1,
      0,1,1
    ), ncol = 3, byrow = TRUE)

    edges <- list(
      c(1,2), c(2,3), c(3,4), c(4,1),
      c(5,6), c(6,7), c(7,8), c(8,5),
      c(1,5), c(2,6), c(3,7), c(4,8)
    )

    edge_lines <- function(V) {
      xs <- ys <- zs <- c()
      for(e in edges) {
        i <- e[1]; j <- e[2]
        xs <- c(xs, V[i,1], V[j,1], NA)
        ys <- c(ys, V[i,2], V[j,2], NA)
        zs <- c(zs, V[i,3], V[j,3], NA)
      }
      list(x = xs, y = ys, z = zs)
    }

    A_full <- matrix(c(
      1.2, 0.3, 0.2,
      0.1, 1.1, -0.4,
      0.2, 0.5, 1.3
    ), 3, 3, byrow = TRUE)

    A_rank_def <- matrix(c(
      1,0,0,
      0,1,0,
      0,0,0
    ), 3, 3, byrow = TRUE)

    V_full     <- vertices %*% t(A_full)
    V_rank_def <- vertices %*% t(A_rank_def)

    orig <- edge_lines(vertices)
    full <- edge_lines(V_full)
    flat <- edge_lines(V_rank_def)

    v1 <- c(0.2, 0.4, 0.1)
    v2 <- c(0.2, 0.4, 0.9)

    v1_full     <- as.vector(A_full %*% v1)
    v2_full     <- as.vector(A_full %*% v2)
    v1_rank_def <- as.vector(A_rank_def %*% v1)
    v2_rank_def <- as.vector(A_rank_def %*% v2)

    vector_arrow <- function(fig, vec, scene = NULL,
                             line_color = NULL, cone_color = NULL,
                             opacity = 1) {
      fig <- fig %>%
        add_trace(
          x = c(0, vec[1]), y = c(0, vec[2]), z = c(0, vec[3]),
          type = "scatter3d", mode = "lines",
          showlegend = FALSE,
          scene = scene,
          opacity = opacity,
          line = if (is.null(line_color)) list() else list(color = line_color)
        )
      fig %>%
        add_trace(
          type = "cone",
          x = vec[1], y = vec[2], z = vec[3],
          u = vec[1], v = vec[2], w = vec[3],
          showscale = FALSE,
          sizemode = "absolute",
          sizeref = 0.2,
          anchor = "tip",
          scene = scene,
          opacity = opacity,
          colorscale = if (is.null(cone_color)) NULL else list(
            c(0, cone_color),
            c(1, cone_color)
          ),
          inherit = FALSE
        )
    }

    connector <- function(fig, from, to, scene = NULL) {
      fig %>%
        add_trace(
          x = c(from[1], to[1]),
          y = c(from[2], to[2]),
          z = c(from[3], to[3]),
          type = "scatter3d",
          mode = "lines",
          showlegend = FALSE,
          scene = scene,
          line = list(color = "grey", dash = "dot")
        )
    }

    p1 <- plot_ly(
      x = orig$x, y = orig$y, z = orig$z,
      type = "scatter3d", mode = "lines",
      showlegend = FALSE
    )
    p1 <- vector_arrow(p1, v1)
    p1 <- vector_arrow(p1, v2)

    p2 <- plot_ly(
      x = full$x, y = full$y, z = full$z,
      type = "scatter3d", mode = "lines",
      showlegend = FALSE,
      scene = "scene2"
    )
    p2 <- vector_arrow(p2, v1, "scene2",
                       line_color = "grey",
                       cone_color = "grey",
                       opacity = 0.5)
    p2 <- vector_arrow(p2, v2, "scene2",
                       line_color = "grey",
                       cone_color = "grey",
                       opacity = 0.5)
    p2 <- vector_arrow(p2, v1_full, "scene2")
    p2 <- vector_arrow(p2, v2_full, "scene2")
    p2 <- connector(p2, v1, v1_full, "scene2")
    p2 <- connector(p2, v2, v2_full, "scene2")

    p3 <- plot_ly(
      x = flat$x, y = flat$y, z = flat$z,
      type = "scatter3d", mode = "lines",
      showlegend = FALSE,
      scene = "scene3"
    )
    p3 <- vector_arrow(p3, v1, "scene3",
                       line_color = "grey",
                       cone_color = "grey",
                       opacity = 0.5)
    p3 <- vector_arrow(p3, v2, "scene3",
                       line_color = "grey",
                       cone_color = "grey",
                       opacity = 0.5)
    p3 <- vector_arrow(p3, v1_rank_def, "scene3")
    p3 <- vector_arrow(p3, v2_rank_def, "scene3")
    p3 <- connector(p3, v1, v1_rank_def, "scene3")
    p3 <- connector(p3, v2, v2_rank_def, "scene3")

    fig <- subplot(p1, p2, p3, nrows = 1) %>%
      layout(
        scene  = list(aspectmode = "cube"),
        scene2 = list(aspectmode = "cube"),
        scene3 = list(aspectmode = "cube"),
        annotations = list(
          list(
            x = 0.16, y = 1.05, text = "Original Domain",
            xref = "paper", yref = "paper",
            showarrow = FALSE
          ),
          list(
            x = 0.50, y = 1.05, text = "Invertible map",
            xref = "paper", yref = "paper",
            showarrow = FALSE
          ),
          list(
            x = 0.84, y = 1.05, text = "Rank-deficient map",
            xref = "paper", yref = "paper",
            showarrow = FALSE
          )
        )
      )

    fig

## Question 1.a

### Uniqueness of \(\hat{\beta}^{\text{OLS}}\)

Suppose that there exists another solution \(\hat{b}\) to the least-squares problem.

Since \(L(\beta)\) is convex, \(\hat{b}\) must be an interior point solution satisfying the FOC, i.e.

\[ (X^\top X)\hat{b} = X^\top y \iff \hat b = (X^\top X)^{-1}X^\top y = \hat{\beta}^{\text{OLS}} \,, \]

so the solution is unique.

## Question 1.a

### What happens if \(X\) is rank deficient

Assume a minimiser \(\hat \beta\) exists.

Since \(X\) is rank deficient, there exists a \(v \neq 0: Xv = 0\).

Then

\[ L(\hat \beta + v) = \| y - X(\hat \beta + v) \|^2_2 = \| y - X \hat \beta - Xv \|_2^2 = \| y - X \hat \beta \|_2^2 = L(\hat{\beta}) \,, \]

so \(\hat \beta + v\) is also a minimiser.

## Question 1.a

### What happens if \(X\) is rank deficient

To understand existence, recall that we seek to minimise

\[ \| y - X \beta \|_2^2 \,, \]

i.e. we seek the \(\beta \in \Re^p\), such that \(X \beta\), which lives in the column span of \(X\), has smallest distance to \(y\).

## Question 1.a

### What happens if \(X\) is rank deficient

Recall that orthogonal projections minimise distance: Projection matrix onto the column span of \(X\) is a matrix \(P\) such that:

1. \(Py \in \text{span}(X)\)
2. \(Py \perp y - Py\)

\[\begin{equation} \begin{aligned} \|y - \eta\|_2^2 &= \|(Py + (I-P)y) - \eta\|_2^2 \\ &= \|(Py - \eta) + (I-P)y\|_2^2 \\ &= \|Py - \eta\|_2^2 + \|(I-P)y\|_2^2 + 2\langle Py - \eta,\,(I-P)y\rangle \\ &= \|Py - \eta\|_2^2 + \|(I-P)y\|_2^2 \\ &\ge \|(I-P)y\|_2^2\,, \end{aligned} \end{equation}\] with equality if \(\eta = Py\).

## Question 1.a

### What happens if \(X\) is rank deficient

Code

    library(jsonlite)

    a <- 0.03
    b <- -0.02
    y <- c(0.7, 0.8, 1.2)

    n <- c(-a, -b, 1)
    n_hat <- n / sqrt(sum(n^2))
    y_proj <- as.numeric(y - sum(y * n_hat) * n_hat)

    u <- c(1, 0, a)
    e1 <- u / sqrt(sum(u^2))

    N <- 61
    L <- 1.5
    t_vals <- seq(-L, L, length.out = N)

    P <- t(sapply(t_vals, function(t) y_proj + t * e1))
    dist <- sqrt(rowSums((P - matrix(y, nrow = N, ncol = 3, byrow = TRUE))^2))
    mid_idx <- ceiling(N / 2)

    B <- max(c(abs(P[,1]), abs(P[,2]), abs(y[1]), abs(y[2]),
               abs(y_proj[1]), abs(y_proj[2]))) + 0.5
    xs <- seq(-B, B, length.out = 30)
    ys <- seq(-B, B, length.out = 30)
    zmat <- outer(xs, ys, function(x, yv) a * x + b * yv)
    z_flat <- as.numeric(t(zmat))

    payload <- list(
      xs    = xs,
      ys    = ys,
      z     = z_flat,
      nx    = length(xs),
      ny    = length(ys),
      y     = y,
      yProj = y_proj,
      Px    = P[,1],
      Py    = P[,2],
      Pz    = P[,3],
      tVals = t_vals,
      dist  = as.numeric(dist),
      midIdx = mid_idx - 1
    )

    json <- toJSON(payload, auto_unbox = TRUE, null = "null", digits = I(16))
    cat(sprintf('<script id="projection-data" type="application/json">%s</script>', json))

# Question 1.b.

2. Ridge regression. The ridge regression parameter estimatê \(\hat{\beta}^{\text{ridge}}\) is defined as a minimiser of the following optimisation problem: \[ \min_{\beta \in \Re^p} L(\beta) + \lambda \|y\|_2^2 = \| y - X\beta \|_2^2 + \lambda \| \beta\|_2^2 \,, \] where \(\lambda \geq 0\) is a regularisation hyperparameter. Assume that \(\lambda > 0\), (if \(\lambda = 0\), then we have the OLS problem).

1. Parameter vector. Show that \[ \hat{\beta}^{\text{ridge}} = (X^\top X + \lambda I_p)^{-1}X^\top y \,, \] where \(I_p\) is a \(p \times p\) identity matrix.

## Question 1.b.1

### Show

- \(\hat{\beta}^{\text{ridge}} = (X^\top X + \lambda I_p)^{-1}X^\top y\)

### Know

- \[ \nabla_x (x^\top a) = a \]

- \[ \nabla_x (x^\top A x) = Ax + A^\top x \]

## Question 1.b.1

\[ \begin{aligned} L(\beta) + \lambda \|\beta \|_2^2 = y^\top y - 2 y^\top X \beta + \beta^\top (X^\top X) \beta + \lambda \beta^\top \beta \,. \end{aligned} \] Differentiating yields: \[ \begin{aligned} \nabla_\beta \left\{ L(\beta) + \lambda \|\beta \|_2^2\right\} &= - 2 X^\top y + 2 (X^\top X) \beta + 2\lambda \beta \overset{!}{=} 0\\ \iff (X^\top X) \beta + \lambda \beta &= X^\top y \\ \iff (X^\top X + \lambda I_p ) \beta&= X^\top y \\ \iff \beta &= (X^\top X + \lambda I_p )^{-1} X^\top y \,. \end{aligned} \] Checking optimality: \[ \begin{aligned} \nabla_\beta \nabla^\top_\beta \left\{ L(\beta) + \lambda \|\beta \|_2^2\right\} &= -2(X^\top X + \lambda I_p) \prec 0, \quad \forall \beta \in \Re^p \,. \end{aligned} \] So the ridge penalised least-squares problem is strictly convex even if \(X\) does not have full column rank. Thus, the ridge estimator exists and is unique regardless of the rank of \(X\).

## Question 1.b.1

### Invertibility of \((X^\top X + \lambda I_p)\)

Assume there is a \(v \in \Re^p: (X^\top X + \lambda I_p)v = 0\). Then: \[ \begin{aligned} v^\top (X^\top X + \lambda I_p) v &= v^\top (X^\top X)v + \lambda v^\top v\\ &= \tilde v^\top \tilde v + \lambda v^\top v, \quad (\tilde v = Xv) \\ &\geq v^\top v \lambda \,. \end{aligned} \]

Hence \((X^\top X + \lambda I_p)v = 0\) must imply that \(v = 0\) for \(\lambda >0\), regardless of the rank of \(X\).

## Question 1.b.1

### Invertibility of \((X^\top X + \lambda I_p)\)

Intuition: Ridge augments the data to (potentially) increase the rank of \(X\):

\[ \begin{aligned} L(\beta) + \lambda \|\beta\|_2^2 &= \| y - X\beta \|_2^2 + \lambda \|\beta\|_2^2 \\ &= \left\| \begin{bmatrix} y \\ 0 \end{bmatrix} - \begin{bmatrix} X \\ \sqrt{\lambda} I_p \end{bmatrix} \beta \right\|_2^2 \,, \end{aligned} \]

## Question 1.b.

2. Covariance of the parameter vector. Assume that we observe \(X\) but not \(y\). Hence, we may regard \(\hat{\beta}^{\text{Ridge}}\) as a random vector conforming to \[ \hat{\beta}^{\text{ridge}} = (X^\top X + \lambda I_p)^{-1}X^\top y \,. \] Consider the eigen-decomposition \[ X^{\top}X = Q D Q^{\top} \] where \(Q^{\top}Q = QQ^{\top} = I_p\) and \(D\) is the diagonal matrix with diagonal elements \(d_1,\ldots,d_p\).

Show that \[ \operatorname{Cov}\bigl[\hat{\beta}^{\text{Ridge}}\bigr] = \sigma^2 Q D^* Q^{\top} \,, \] where \(D^*\) is a diagonal matrix with diagonal elements \(d_1/(d_1+\lambda)^2,\ldots,d_p/(d_p+\lambda)^2\).

Remark: covariance of a random column vector \(Z\) is defined by \[ \operatorname{Cov}[Z] = \mathrm{E}\bigl[(Z-\mathrm{E}[Z])(Z-\mathrm{E}[Z])^{\top}\bigr]\,. \]

## Question 1.b.2

### Show

- \(\text{Cov}(\hat{\beta}^{\text{ridge}}) = \sigma^2 Q D^* Q^{\top}\)

### Know

- \(\text{Cov}[Z] = \mathrm{E}[(Z-\mathrm{E}[Z])(Z-\mathrm{E}[Z])^{\top}]\)

- \(X^\top X = QDQ^\top\)

- \(QQ^\top = Q^\top Q = I_p\)

## Question 1.b.2

### Simplifying variance-covariance matrix

\[ \begin{aligned} \text{Cov}[Z] &= \mathrm{E}[(Z-\mathrm{E}[Z])(Z-\mathrm{E}[Z])^{\top}] \\ &= \mathrm{E} [ZZ^\top] - \mathrm{E}[Z \mathrm{E}[Z]^\top] - \mathrm{E}[\mathrm{E}[Z]Z^\top] + \mathrm{E}[\mathrm{E}[Z]\mathrm{E}[Z]^\top] \\ &= \mathrm{E} [ZZ^\top] - \mathrm{E}[Z]\mathrm{E}[Z]^\top \,. \end{aligned} \]

## Question 1.b.2

### Calculating quantities

Using that

\[ y = X\beta + \epsilon, \quad \mathrm{E}[y] = X\beta, \quad \mathrm{E}[yy^\top] = \sigma^\top X\beta \beta^\top X^\top \] \[ \begin{aligned} \mathrm{E}[\hat{\beta}^{\text{ridge}}(\hat{\beta}^{\text{ridge}})^\top] &= \mathrm{E}[(X^\top X + \lambda I_p)^{-1}X^\top y((X^\top X + \lambda I_p)^{-1}X^\top y)^\top] \\ &= \mathrm{E}[(X^\top X + \lambda I_p)^{-1}X^\top yy^\top X (X^\top X + \lambda I_p)^{-1}] \\ &= (X^\top X + \lambda I_p)^{-1}X^\top\mathrm{E}[yy^\top]X (X^\top X + \lambda I_p)^{-1} \\ &= \sigma^2 (X^\top X + \lambda I_p)^{-1}X^\top X (X^\top X + \lambda I_p)^{-1} \\ & +(X^\top X + \lambda I_p)^{-1}X^\top X \beta \beta^\top X^\top X (X^\top X + \lambda I_p)^{-1} \,. \end{aligned} \]

\[ \mathrm{E}[\hat{\beta}^{\text{ridge}}] = (X^\top X + \lambda I_p)^{-1}X^\top \mathrm{E}(y) = (X^\top X + \lambda I_p)^{-1}X^\top X\beta \,. \]

\[ \mathrm{E}[\hat{\beta}^{\text{ridge}}] \mathrm{E}[\hat{\beta}^{\text{ridge}}]^\top = (X^\top X + \lambda I_p)^{-1}X^\top X\beta \beta^\top X^\top X (X^\top X + \lambda I_p)^{-1} \,. \]

## Question 1.b.2

### Putting it together

\[ \begin{aligned} \text{Cov}(\hat{\beta}^{\text{ridge}}(\hat{\beta}^{\text{ridge}})^\top) &= \mathrm{E}[\hat{\beta}^{\text{ridge}}(\hat{\beta}^{\text{ridge}})^\top] - \mathrm{E}[\hat{\beta}^{\text{ridge}}] \mathrm{E}[\hat{\beta}^{\text{ridge}}]^\top \\ &= \sigma^2 (X^\top X + \lambda I_p)^{-1}X^\top X (X^\top X + \lambda I_p)^{-1} \end{aligned} \]

### Shortcut

- \(\text{var}(A y) = A \text{var}(y) A^\top\) for random vector \(y\), nonrandom matrix \(A\)

\[ \begin{aligned} \text{var}(\hat{\beta}^{\text{ridge}}) &= \text{var}\left( (X^\top X + \lambda I_p)^{-1} X^\top y \right) \\ &= (X^\top X + \lambda I_p)^{-1} X^\top\text{var}(y)\left[(X^\top X + \lambda I_p)^{-1} X^\top \right]^\top \\ &= \sigma^2 (X^\top X + \lambda I_p)^{-1} X^\top X (X^\top X + \lambda I_p)^{-1} \end{aligned} \]

## Question 1.b.2

### SVD

Using that \(X = QDQ^\top\), with \(QQ^\top = Q^\top Q = I_p\), we have

\[ X^\top X = (QDQ^\top)^\top QDQ^\top = QDQ^\top QDQ^\top = QDDQ^\top \,. \]

## Question 1.b.2

### SVD

To figure out the inverse of \(X^\top X + \lambda I_p\), we note that \(Q^{-1} = Q^\top\) \[ X^\top X + \lambda I_p = QDQ^\top + \lambda QQ^\top = Q(D + \lambda I_p) Q^\top \,. \]

\[ D + \lambda I_p = \begin{bmatrix} D_1 + \lambda & 0 & \cdots & 0 \\ 0 & D_2 + \lambda & \ddots & \vdots \\ \vdots & 0 & \ddots & \vdots \\ 0 & \cdots & \cdots & D_p + \lambda \end{bmatrix} \,, \] so that

\[ (D + \lambda I_p)^{-1} = \begin{bmatrix}\frac{1}{D_1+\lambda} & 0 & \cdots & 0 \\ 0 & \frac{1}{D_2+\lambda} & \ddots & \vdots \\ \vdots & 0 & \ddots & \vdots \\ 0 & \cdots & \cdots & \frac{1}{D_p+\lambda} \end{bmatrix} \,. \]

## Question 1.b.2

### SVD

\[ \begin{aligned} (X^\top X + \lambda I_p)^{-1} &= (QDDQ^\top + \lambda QQ^\top)^{-1} \\ &= [Q(D + \lambda I_p) Q^\top]^{-1} \\ &= (Q^T)^{-1} (D + \lambda I_p)^{-1} Q^{-1} \\ &= Q (D + \lambda I_p)^{-1} Q^\top \,, \end{aligned} \]

where we used that for two invertible matrices \(A,B\),

\[ (AB)^{-1} = B^{-1}A^{-1} \,. \]

## Question 1.b.2

### Simplifying

\[ \begin{aligned} \sigma^2 (X^\top X + \lambda I_p)^{-1} X^\top X (X^\top X + \lambda I_p)^{-1} &= \sigma^2 Q (D + \lambda I_p)^{-1} Q^\top Q D Q^\top Q (D + \lambda I_p)^{-1}Q^\top \\ &= \sigma^2 Q(D + \lambda I_p)^{-1} D (D + \lambda I_p)^{-1}Q^\top \\ &= \sigma^2 QD^*Q^\top \\ \end{aligned} \] with \[ D^* = \begin{bmatrix}\frac{D_1}{(D_1+\lambda)^2} & 0 & \cdots & 0 \\ 0 & \frac{D_2}{(D_2+\lambda)^2} & \ddots & \vdots \\ \vdots & 0 & \ddots & \vdots \\ 0 & \cdots & \cdots & \frac{D_p}{(D_p+\lambda)^2} \end{bmatrix} \,. \]

## Question 1.b.3

> Using the assertions in (b.1) and (b.2), discuss how \(\hat \beta^{\text{ridge}}\) behaves as \(\lambda\) grows large.

### Shrinkage

\[ \begin{aligned} \hat{\beta}^{\text{ridge}} = (X^\top X + \lambda I_p)^{-1}X^\top y = Q(D + \lambda I_p)^{-1} Q^\top X^\top y \,. \end{aligned} \]

If \(\lambda \to \infty\), \(D_i / (D_i + \lambda)^2 \to 0\), so \(\hat \beta^{\text{ridge}} \to 0\)

### Bias

\[ \begin{aligned} \mathrm{E}[\hat{\beta}^{\text{ridge}}] = Q (D + \lambda I_p)^{-1} Q^\top Q D Q^\top \beta = Q (D + \lambda I_p)^{-1} D Q^\top \beta \end{aligned} \] If \(\lambda = 0\), Bias is zero, if \(\lambda \to \infty\), bias is \(-\beta\).

### Variance

As \(\lambda \to \infty\), variance goes to zero (Classic bias variance tradeoff)

## Learnings

- Be comfortable with basic matrix calculus
- Understand limitations of OLS when \(X\) is collinear
- Comparison of Ridge to OLS in linear model
  - Understand existence, shrinkage, bias-variance (why we like ridge for prediction)

## Check your understanding

You observe a linear model \(y = X\beta + \varepsilon\) with \(X \in \Re^{n \times p}\), and \(X\) is known to be rank-deficient (\(\operatorname{rank}(X) < p\)).

1. Does an OLS minimiser of \(\|y - X\beta\|_2^2\) always exist?

Answer

Yes, a minimiser of \(\|y - X\beta\|_2^2\) always exists; geometrically, we are projecting onto the closed subspace \(\operatorname{col}(X)\)

2. Is it necessarily unique?

Answer

It is not unique when \(\operatorname{rank}(X) < p\). There are infinitely many minimisers; if \(\hat\beta\) is one minimiser and \(v \neq 0: Xv = 0\), then \(\hat\beta + v\) is also a minimiser.

3. How does ridge regression change the answer to (2)?

Answer

Ridge regression minimises \(\|y - X\beta\|_2^2 + \lambda\|\beta\|_2^2\) with \(\lambda > 0\). The Hessian is \(2(X^\top X + \lambda I_p)\), which is positive definite for any \(\lambda > 0\), so the solution is unique even when \(X\) is rank-deficient.

## Check your understanding

Let \(X^\top X = Q D Q^\top\) with \(D = \operatorname{diag}(d_1,\dots,d_p)\) and \(Q^\top Q = I_p\). Consider the ridge estimator \[ \hat\beta^{\text{ridge}} = (X^\top X + \lambda I_p)^{-1} X^\top y. \]

Suppose \(n \ll p\), and the columns of \(X\) are highly collinear. You care only about minimising test error, not about unbiasedness of \(\hat\beta\).

Give two clear reasons why ridge regression can outperform OLS in this regime, and one situation where making \(\lambda\) too large will start to harm predictive performance.

Answer

- OLS has very high variance and is unstable when \(n \ll p\) and predictors are collinear; small perturbations in the data can cause large changes in \(\hat\beta\). Ridge stabilises the estimator by shrinking coefficients, reducing variance of predictions.
- Ridge effectively downweights directions in predictor space that are poorly informed by the data (small eigenvalues of \(X^\top X\)), which reduces overfitting to noise and typically improves generalisation.
- If \(\lambda\) is made too large, the fitted model is overly shrunk toward zero, underfitting the true signal. In this case, the bias in predictions becomes dominant and test error increases, even though the variance is small.

## Question 1.c

> Penalised least squares analysis. To gain some insights about variable selection procedures, consider the case when \(X\) has orthonormal columns, i.e. \(X^\top X = I_p\)

1. Separability. Consider the following penalised least squares problem: find \(\hat\beta\) that is a minimiser of the following optimisation problem: \[ \min_{\beta \in \Re^p} \frac{1}{2} \| y - X\beta \|_2^2+ \sum_{j=1}^p c_j(|\beta_j|,\lambda) \,, \] where each \(c_j : \Re_+ \times \Re_+^m \to \Re_+\) is an increasing function. Define \[ z = X^\top y = (z_1, \dots, z_p)^\top \,. \] Show that the penalised least squares problem is equivalent to finding \(\hat\beta\) that is a minimiser of the following optimisation problem: \[ \min_{\beta \in \Re^p} \frac{1}{2} \sum_{j=1}^p (z_j - \beta_j)^2 + \sum_{j=1}^p c_j(|\beta_j|,\lambda) \,. \]

## Question 1.c.1

### Show

- penalised least squares is separable componentwise

### Know

- \(z = X^\top y\)

- \(X^\top X = I_p\)

## Question 1.c.1

\[ \begin{aligned} \| y - X\beta \|_2^2 &= (y - X\beta)^\top (y - X\beta) \\ &= y^\top y - 2y^\top X\beta + \beta^\top (X^\top X) \beta \\ &= y^\top y - 2z^\top \beta + \beta^\top \beta \\ % &= y^\top I_py - 2z^\top \beta + \beta^\top \beta \\ % &= y^\top y - 2z^\top \beta + \beta^\top \beta \\ &= z^\top z - 2z^\top \beta + \beta^\top \beta + (y^\top y - z^\top z) \\ &= \|z - \beta \|_2^2 + C, \quad C \perp \beta \end{aligned} \]

## Question 1.c.1

Moreover, show that the last optimisation problem is equivalent to componentwise solving of the following penalised least squares problem: find \(\hat\theta\) that is a minimiser of the optimisation problem: \[ \min_{\theta \in \Re} \frac{1}{2} (z - \theta)^2 + c(|\theta|,\lambda) \,, \] where for each \(j \in \{1,\dots,p\}\), \(z \equiv z_j\) and \(c \equiv c_j\).

\[ \begin{aligned} \min_{\beta \in \Re^p} \, \frac{1}{2} \|z - \beta\|^2 + \sum_{j = 1}^p c(|\beta_i|,\lambda) &= \min_{\beta \in \Re^p} \,\sum_{j = 1}^p \frac{1}{2} (z_j - \beta_j)^2 + \sum_{j = 1}^p c(|\beta_j|,\lambda)\\ &= \min_{\beta \in \Re^p}\sum_{j = 1}^p \left(\frac{1}{2} (z_j - \beta_j)^2 + c(|\beta_j|,\lambda) \right) \\ &= \sum_{j = 1}^p \min_{\beta_j \in \Re} \, \left(\frac{1}{2} (z_j - \beta_j)^2 + c(|\beta_j|,\lambda) \right) \,. \end{aligned} \]

## Question 1.c.1

### Intuition separability

We can think of orthonormal designs as independence of features, in that no feature holds any information about any other feature.

In that case, it makes sense that a change in the \(j\)th feature (i.e. column of \(X\)), should only affect our choice of the \(j\)th component of \(\beta\).

Mathematically, the cross term

\[ \beta^\top (X^\top X) \beta = \sum_{i = 1}^p \sum_{j = 1}^p \beta_i \beta_j x^\top_i x_j \] captures all the dependencies of the components of \(\beta\).

Statistically, we can think about the case where all the rows of \(X\) are i.i.d. with variance-covariance matrix \(I_p\).

In that case, we cannot learn anything about the signal by considering dependence of features, so we should ignore it.

\[ \hat{\beta}^{\text{OLS}} = (X^\top X)^{-1}X^\top (X\beta + \epsilon) = \beta + X^\top \epsilon, \quad \text{var}(X^\top \epsilon) = \sigma^2 I_p \,. \]

## Question 1.c.2

2. For the \(\ell_2\) penalty, the optimisation problem is to find \(\hat{\theta}^{\text{ridge}}\) that minimises \[ \frac{1}{2} (z - \theta)^2 + \frac{\lambda}{2} \theta^2 \,, \] over \(\theta \in \Re\).

Code

    lambda <- 0.1
    z <- 3

    fun <- function(theta){
      0.5 * (z - theta)^2 + 0.5 * lambda * theta^2
    }

    thetas <- seq(from = -10, to = 10, length.out = 1000)

    dfun <- function(theta){-(z - theta)  + lambda * theta}

    par(mfrow = c(1,2))

    plot(thetas, fun(thetas), type = "l", linewidth = 2, main = "Ridge loss")

    plot(thetas, dfun(thetas), type = "l", linewidth = 2, main = "Ridge loss derivative")
    abline(h = 0, col = "red", linetype = "dash")

\[ \frac{\partial}{\partial \theta} \left(\frac{1}{2}(z - \theta)^2 + \frac{\lambda}{2}\theta^2 \right) = \theta - z + \lambda \theta \overset{!}{=}0 \implies \hat{\theta}^{\text{ridge}} = \frac{z}{1 + \lambda} \,. \]

## Question 1.c.2

### Ridge shrinkage

Code

    zs <- seq(from = -5, to = 5, length.out = 1000)

    lambda <- 1

    ridge <- function(z){
      z / (1 + lambda)
    }

    plot(zs, ridge(zs), type = "l", col = "blue", lwd = 2, main = "Ridge shrinkage", xlab = "z", ylab = expression(hat(theta)))
    lines(zs,zs, col = "blue", lwd = 2, lty = "dashed")

## Question 1.c.3

3. For \(\ell_1\) penalty, the optimisation problem is to find \(\hat{\theta}^{\text{soft}}\) that minimises \[ L(\theta) = \frac{1}{2}(z - \theta)^2 + \lambda |\theta| \,. \]

## Question 1.c.3

Problem:

- At \(\theta = 0\), the derivative is not defined
- Need to deal with the case \(\theta = 0\) separately

Code

    lambda <- 3.
    z <- 5

    fun <- function(theta){
      0.5 * (z - theta)^2 + lambda * abs(theta)
    }

    thetas <- seq(from = -10, to = 10, length.out = 1000)

    dfun <- function(theta){-(z - theta) + lambda * sign(theta)}

    par(mfrow = c(1,2))

    plot(thetas, fun(thetas), type = "l", linewidth = 2, main = "Lasso loss")
    abline(v = 0, col = "red", linetype = "dash")
    plot(thetas, dfun(thetas), type = "l", linewidth = 2, main = "Lasso loss derivative")
    abline(v = 0, col = "red", linetype = "dash")

## Question 1.c.3

### When is \(\hat{\theta}^{\text{soft}} = 0\) optimal?

If for all \(\theta \in \Re \setminus \{0\}\):

\[ \begin{aligned} L(\theta) &\geq L(0) \\ \frac{1}{2} (z - \theta)^2 + \lambda |\theta| &\geq \frac{1}{2}z^2 \\ \frac{1}{2}z^2 - \theta z + \theta^2+ \lambda |\theta| &\geq \frac{1}{2}z^2 \\ - \text{sign}(\theta) z |\theta| + |\theta|^2+ \lambda |\theta| &\geq 0 \\ - \text{sign}(\theta) z + |\theta| + \lambda &\geq 0 \\ \end{aligned} \]

## Question 1.c.3

### When is \(\hat{\theta}^{\text{soft}} = 0\) optimal?

\[ |\theta| + \lambda \geq \text{sign}(\theta) z \]

1. \(\lim\limits_{|\theta| \to \infty} (- \text{sign}(\theta) z + |\theta| + \lambda) = \infty > 0\)

2. \(\lim\limits_{|\theta| \to 0^+} (- \text{sign}(\theta) z + |\theta| + \lambda) = \lim\limits_{|\theta| \to 0^+} - z + \theta + \lambda \geq 0\) so that \(z \leq \lambda\)

3. \(\lim\limits_{|\theta| \to 0^-} (- \text{sign}(\theta) z + |\theta| + \lambda) = \lim\limits_{|\theta| \to 0^+} z - \theta + \lambda \geq 0\) so that \(z \geq -\lambda\)

Hence, if \(\hat{\theta}^{\text{soft}} = 0\) is optimal, then \(z \in [-\lambda, \lambda]\).

Therefore, whenever \(|z| > \lambda\), \(\hat{\theta}^{\text{soft}} = 0\) is not optimal

In this case, \(L(\theta)\) is differentiable and we can find solutions as usual

## Question 1.c.3

### \(\hat{\theta}^{\text{soft}} \neq 0\)

- If \(|z| > \lambda\), \(\hat{\theta}^{\text{soft}} = 0\) is not optimal

- \(L(\theta)\) is convex, differentiable, so solution must be characterised by FOCs

\[ \frac{\partial}{\partial \theta} L(\theta) = \begin{cases} -(z - \theta) + \lambda & \text{if } \theta > 0, \\ -(z - \theta) - \lambda & \text{if } \theta < 0 \,. \end{cases} \]

1. Assume that \(\hat{\theta}^{\text{soft}} > 0\) is optimal:

\[ \frac{\partial}{\partial \theta} \left(\frac{1}{2} (z - \theta)^2 + \lambda |\theta| \right) = \frac{\partial}{\partial \theta} \left(\frac{1}{2} (z - \theta)^2 + \lambda \theta\right) = -(z - \theta) + \lambda \overset{!}{=} 0 \] so that \[ \hat{\theta}^{\text{soft}} = z - \lambda > 0 \iff z > \lambda \,. \]

## Question 1.c.3

### \(\hat{\theta}^{\text{soft}} \neq 0\)

2. Assume that \(\hat{\theta}^{\text{soft}} < 0\) is optimal:

\[ \frac{\partial}{\partial \theta} \left(\frac{1}{2} (z - \theta)^2 + \lambda |\theta|\right) = \frac{\partial}{\partial \theta} \left( \frac{1}{2} (z - \theta)^2 - \lambda \theta\right) = -(z - \theta) - \lambda \overset{!}{=} 0 \] so that \[ \hat{\theta}^{\text{soft}} = z + \lambda < 0 \iff z < -\lambda \,. \]

## Question 1.c.3

- We found three candidate minimisers \(\hat{\theta}^{\text{soft}}\) and necessary conditions on \(z, \lambda\) for optimality
- We need to verify that these are also sufficient
- We check the three cases
  1. \(z > \lambda\)
  2. \(z < - \lambda\)
  3. \(|z| \leq \lambda\)

## Question 1.c.3

### \(z > \lambda\)

- For \(\theta > 0\), \[ \frac{\partial}{\partial \theta}L(\theta) = -(z - \theta) + \lambda = (\lambda - z) + \theta \,, \] is strictly increasing, negative for \(\theta < z - \lambda\), positive for \(\theta > z - \lambda\).

- For \(\theta < 0\), \[ \frac{\partial}{\partial \theta}L(\theta) = -(z - \theta) - \lambda = \theta - (z - \lambda) < 0\,, \] is strictly decreasing, and the minimum is at the boundary \(\theta = 0\).

- For \(\theta = 0\), \[ L(z-\lambda) - L(0) = \frac{1}{2}\lambda^2 + \lambda (z - \lambda) - \frac{1}{2} z^2 = -\frac{1}{2}(z - \lambda)^2 < 0 \,. \]

- If \(z > \lambda\), then \(\hat \theta ^{\text{soft}} = z - \lambda\) is the unique minimiser of \(L(\theta)\)

## Question 1.c.3

### \(|z| \leq \lambda\)

- Have already shown that for any \(\theta \in \Re \setminus \{0\}\), \(L(\theta) > L(0)\)
- Could also check derivatives on both half lines \((-\infty, 0)\), \((0, \infty)\) and conclude solution must be on boundary (i.e. \(\theta = 0\))

### \(z < - \lambda\)

- Similar to case \(z > \lambda\)

In conclusion \[ \begin{aligned} \hat{\theta}^{\text{soft}} &= \begin{cases} z + \lambda & \text{if } z < - \lambda \\ 0 & \text{if } |z| \leq \lambda \\ z - \lambda & \text{if } z > \lambda \end{cases} \quad = \mathbb{1} \{|z| > \lambda \} \text{sign}(z)(|z|-\lambda) \,. \end{aligned} \]

## Question 1.c.3

In conclusion \[ \begin{aligned} \hat{\theta}^{\text{soft}} &= \begin{cases} z + \lambda & \text{if } z < - \lambda \\ 0 & \text{if } |z| \leq \lambda \\ z - \lambda & \text{if } z > \lambda \end{cases} \\ &= \mathbb{1} \{|z| > \lambda \} \text{sign}(z)(|z|-\lambda) \,. \end{aligned} \]

Code

    zs <- seq(from = -5, to = 5, length.out = 1000)

    lambda <- 1

    st <- function(z){
      as.numeric(abs(z) > lambda) * sign(z) * (abs(z) - lambda)
    }

    plot(zs, st(zs), type = "l", col = "blue", lwd = 2, main = "Soft thresholding", xlab = "z", ylab = expression(hat(theta)))
    lines(zs, zs, col = "blue", lwd = 1, lty = "dashed")
    abline(v = -lambda, col = "grey", lwd = 2, lty = "dashed")
    abline(v = lambda, col = "grey", lwd = 2, lty = "dashed")

## Question 1.c.4

4. For the hard-threshold penalty, the optimisation problem is to find \(\hat{\theta}^{\text{hard}}\) that minimises

\[ L(\theta) = \frac{1}{2} (z - \theta)^2 + \frac{\lambda^2}{2} \mathbb{1}\{ |\theta| > 0 \} \]

Code

    lambda <- 3.
    z <- 3

    fun <- function(theta){
      0.5 * (z - theta)^2 + lambda^2 / 2 * as.numeric(abs(theta) > 0)
    }

    thetas <- seq(from = -10, to = 10, length.out = 1000)

    dfun <- function(theta){-(z - theta) + lambda * sign(theta)}

    par(mfrow = c(1,2))

    plot(thetas, fun(thetas), type = "l", linewidth = 2, main = "Hard threshold loss")
    abline(v = 0, col = "red", linetype = "dash")
    points(0, 0.5*z^2, col = "black", pch = 20)
    plot(thetas, dfun(thetas), type = "l", linewidth = 2, main = "Hard threshold loss derivative")
    abline(v = 0, col = "red", linetype = "dash")

## Question 1.c.4

### When is \(\hat{\theta}^{\text{hard}} = 0\) optimal?

If for all \(\theta \in \Re \setminus \{0\}\):

\[ \begin{aligned} L(\theta) &\geq L(0) \\ \frac{1}{2} (z - \theta)^2 + \frac{\lambda^2}{2} \mathbb{1}\{|\theta| > 0 \} & \geq \frac{1}{2}z^2 \\ \frac{1}{2}z^2 - \theta z + \theta^2+ \frac{\lambda^2}{2} & \geq \frac{1}{2}z^2 \\ - z \theta + \frac{1}{2}\theta^2 + \frac{1}{2}\lambda^2 &\geq 0 \\ \theta^2 - 2z\theta + \lambda^2 &\geq 0 \\ \theta^2 - 2z\theta + z^2 + \lambda^2 - z^2 &\geq 0 \\ (\theta - z)^2 +\lambda^2 - z^2 &\geq 0 \\ \end{aligned} \]

- \((\theta - z)^2 \geq 0\) and equal iff \(z = \theta\)
- In that case, need \(\lambda^2 \geq z^2\) or, equivalently \(|z| \leq \lambda\)

## Question 1.c.4

### \(\hat{\theta}^{\text{hard}} \neq 0\)

\[ L(\theta) = (z - \theta)^2 \,, \]

which is minimised at \(\hat{\theta}^{\text{hard}} = z\).

### Checking sufficiency

Similar to c.3

- In conclusion

\[ \hat{\theta}^{\text{hard}} = \begin{cases} 0 & \text{if } |z| \leq \lambda\,, \\ z & \text{if } |z| > \lambda \,. \end{cases} \]

## Question 1.c.4

### Hard threshold

Code

    zs <- seq(from = -5, to = 5, length.out = 1000)

    lambda <- 1

    ht <- function(z){
      ifelse(abs(z)<= lambda, 0, z)
    }

    hts <- ht(zs)

    plot(zs[hts < 0], hts[hts<0], type = "l", col = "blue", lwd = 2, main = "Hard thresholding", xlim = c(-5,5), ylim = c(-5,5), xlab = "z", ylab = expression(hat(theta)))
    lines(zs[hts == 0], hts[hts == 0], type = "l", col = "blue", lwd = 2)
    lines(zs[hts > 0], hts[hts>0], type = "l", col = "blue", lwd = 2)
    lines(zs, zs, col = "blue", lwd = 1, lty = "dashed")
    abline(v = -lambda, col = "grey", lwd = 2, lty = "dashed")
    abline(v = lambda, col = "grey", lwd = 2, lty = "dashed")

## Question 1.c.4

### \(k\)-best subset selection

- Given a set \(S \subset \{1, 2, \ldots, n\}: |S| = k\), find \(\hat{\beta}_S\) is the \(k\)-vector that minimises \[ \|y - X_{S} \hat{\beta}_{S} \|_2^2 = \| Z_S - X_S^\top X_S \beta \|_2^2 + C = \| Z_S - \beta_S \|_2^2 \,, \] where \(X_S\) is \(X\) with only columns from \(S\) and \(Z_S = X_S^\top y\). Hence \(\hat \beta_S = Z_S\)

- Overall loss for \(S\): \[ \| y - X_S\hat{\beta}_S \|_2^2 = \| y - X_SX_S^\top y \| = \| y\|_2^2 - \| X_S^\top y\|_2^2 = C - \|Z_S \|_2^2 \]

- Note that \(Z_S\) is the \(k\)-subvector of \(Z\) with all the entries selected by \(S\)

\[ \| y - X_S\hat{\beta}_S \|_2^2 = C - \sum_{i \in S} z_i^2 \]

and \(k\)-best subset picks the \(k\) entries for which \(|z_i|\) is largest

## Question 1.c.4

### \(k\)-best subset selection vs. Hard threshold

- \(k\)-best subset: Pick \(k\) entries for which \(|z_i|\) is largest

- Hard threshold: \[ \hat{\theta}_i^{\text{hard}} = \begin{cases} 0 & \text{if } |z_i| \leq \lambda\,, \\ z_i & \text{if } |z_i| > \lambda \,. \end{cases} \]

- Order \(z_{(n)} \geq z_{(n-1)} \geq \ldots z_{(n - k)} \geq z_{(n - k - 1)} \geq \ldots\)

- If we choose \(z_{(n-k)} > \lambda \geq z_{(n-k-1)}\), we get \(k\)-nonzero entries and the solution is equivalent

## Question 1.c.5

### Elastic net

For the elastic net penalty, the optimisation problem is to find \(\hat{\theta}^{\text{en}}\) that minimises

\[ L(\theta) = \frac{1}{2}(z - \theta)^2 + \lambda_1 \theta^2 + \lambda_2| \theta | \,. \]

Code

    lambda <- 2.
    lambda2 <- .5
    z <- 3

    fun <- function(theta){
      0.5 * (z - theta)^2 + lambda * abs(theta) + lambda2 * theta^2
    }

    thetas <- seq(from = -5, to = 5, length.out = 1000)

    dfun <- function(theta){-(z - theta) + lambda * sign(theta) + lambda2 * theta}

    par(mfrow = c(1,2))

    plot(thetas, fun(thetas), type = "l", linewidth = 2, main = "Elastic net loss")
    abline(v = 0, col = "red", linetype = "dash")
    plot(thetas, dfun(thetas), type = "l", linewidth = 2, main = "Elastic net loss derivative")
    abline(v = 0, col = "red", linetype = "dash")

## Question 1.c.5

### When is \(\hat{\theta}^{\text{en}} = 0\) optimal?

If for all \(\theta \in \Re\): \[ \begin{aligned} L(\theta) &\geq L(0) \\ \frac{1}{2} (z - \theta)^2 + \lambda_1 \theta^2 + \lambda_2 |\theta| &\geq \frac{1}{2} z^2 \\ \frac{1}{2} z^2 - \theta z + \frac{1}{2} \theta^2 + \lambda_1 \theta^2 + \lambda_2 |\theta| &\geq \frac{1}{2} z^2 \\ \left(\frac{1}{2} + \lambda_1 \right)(|\theta|)^2 + |\theta|(\lambda_2 - \text{sign}(\theta)z) &\geq 0 \\ \left(\frac{1}{2} + \lambda_1 \right) |\theta| + \lambda_2 - \text{sign}(\theta)z &\geq 0 \\ \left(\frac{1}{2} + \lambda_1 \right) |\theta| + \lambda_2 &\geq \text{sign}(\theta) z \,. \end{aligned} \]

## Question 1.c.5

### When is \(\hat{\theta}^{\text{en}} = 0\) optimal?

1. \(\forall \theta > 0\):

\[ \left(\frac{1}{2} + \lambda_1 \right) \theta + \lambda_2 \geq z \implies z \leq \lambda_2 \]

2. \(\forall \theta < 0\):

\[ -\left(\frac{1}{2} + \lambda_1 \right) \theta + \lambda_2 \geq -z \implies z \geq -\lambda_2 \]

In conclusion, \(\hat{\theta}^{\text{en}} = 0\) is optimal whenever \(|z| \leq \lambda_2\)

## Question 1.c.5

### \(\hat{\theta}^{\text{en}} \neq 0\)

1. \(\hat{\theta}^{\text{en}} > 0\)

\[ \frac{\partial}{\partial \theta} L(\theta) = \frac{\partial}{\partial \theta} \left(\frac{1}{2}(z - \theta)^2 +\lambda_1 \theta^2 + \lambda_2 \theta \right) = - (z - \theta) + 2 \lambda_1 \theta + \lambda_2 \overset{!}{=} 0 \]

\[ \hat{\theta}^{\text{en}} = \frac{z - \lambda_2}{1 + 2 \lambda_2} > 0 \iff z > \lambda_2 \,. \]

2. \(\hat{\theta}^{\text{en}} < 0\)

\[ \frac{\partial}{\partial \theta} L(\theta) = \frac{\partial}{\partial \theta} \left(\frac{1}{2}(z - \theta)^2 +\lambda_1 \theta^2 - \lambda_2 \theta \right) = - (z - \theta) + 2 \lambda_1 \theta - \lambda_2 \overset{!}{=} 0 \]

\[ \hat{\theta}^{\text{en}} = \frac{z + \lambda_2}{1 + 2 \lambda_2} < 0 \iff z <- \lambda_2 \,. \]

## Question 1.c.5

### Checking sufficiency

Similar to c.3

### Elastic net shrinkage

\[ \hat{\theta}^{\text{en}} = \frac{\mathbb{1}\{|z| > \lambda_2 \} \text{sign}(z)(|z|-\lambda_2)}{1 + 2 \lambda_1} \,. \]

Code

    zs <- seq(from = -5, to = 5, length.out = 1000)

    lambda1 <- .25
    lambda2 <- .5

    en <- function(z){
      as.numeric(abs(z) > lambda1) * sign(z) * (abs(z) - lambda1) / (1 + 2 * lambda2)
    }

    plot(zs, en(zs), type = "l", col = "blue", lwd = 2, main = "Elastic net shrinkage", xlab = "z", ylab = expression(hat(theta)))
    lines(zs, zs, col = "blue", lwd = 1, lty = "dashed")
    abline(v = -lambda1, col = "grey", lwd = 2, lty = "dashed")
    abline(v = lambda1, col = "grey", lwd = 2, lty = "dashed")

## Learnings

- Understand implications of orthonormal designs
- Understand shrinkage properties of different penalties
- Deal with nondifferentiable, convex functions

## Check your understanding

> Why can ridge beat thresholding rules for prediction?

Assume the orthonormal design setting with scalar problems \[ \hat\theta(z) = \arg\min_{\theta \in \mathbb{R}} \left\{ \frac{1}{2}(z-\theta)^2 + c(|\theta|,\lambda)\right\}, \] and compare ridge to soft/hard thresholding when you only care about test MSE, not sparsity.

Answer

- Stability: Ridge uses a smooth map \(z \mapsto \hat\theta^{\text{ridge}}(z) = z/(1+\lambda)\), so small changes in \(z\) give small changes in \(\hat\theta\). Thresholding rules introduce a discontinuity at \(|z| = \lambda\), so tiny perturbations around the threshold can flip a coefficient from zero to nonzero (or vice versa), which increases prediction variance.
- Better use of many weak signals. With many small but real effects, ridge shrinks all of them instead of zeroing them out. Thresholding rules aggressively delete coefficients below a cut-off, which can throw away cumulative signal and hurt predictive performance even if the resulting model is sparser.

## Check your understanding

> Why is soft thresholding often preferred to hard thresholding for prediction?

Both rules threshold \(z\), but in different ways. Explain why the soft rule is usually safer for test error.

Answer

- Soft is continuous, hard is not. Soft thresholding \[ \hat\theta^{\text{soft}}(z) = \mathbb{1}\{|z|>\lambda\}\,\text{sign}(z)\,(|z| - \lambda) \] is continuous in \(z\). Hard thresholding \[ \hat\theta^{\text{hard}}(z) = \begin{cases} 0 & \text{if } |z|\le\lambda,\\ z & \text{if } |z|>\lambda \end{cases} \] has a jump at \(|z|=\lambda\). This jump makes fits highly sensitive to noise near the threshold and inflates test error.
- Convex vs non-convex. Soft thresholding comes from a convex \(\ell_1\) penalty, giving a unique global minimiser and well-behaved optimisation. Hard thresholding corresponds to a non-convex penalty; in general designs this leads to multiple local minima and more unstable model selection, which again tends to worsen predictive risk.

## Check your understanding

> Bias–variance tradeoff of ridge in the orthonormal case

Assume \(z = \beta + \epsilon\) with \(\epsilon \sim (0,\sigma^2)\) and ridge estimator \[ \hat\theta^{\text{ridge}}(z) = \frac{z}{1+\lambda}. \] Describe how bias and variance depend on \(\lambda\), and why a nonzero \(\lambda\) can reduce test error.

Answer

- Bias. \[ \mathbb{E}[\hat\theta^{\text{ridge}}] = \frac{\beta}{1+\lambda}, \quad \text{bias}(\hat\theta^{\text{ridge}}) = -\frac{\lambda}{1+\lambda}\,\beta. \] Bias grows in magnitude as \(\lambda\) increases (stronger shrinkage toward \(0\)).
- Variance. \[ \operatorname{Var}(\hat\theta^{\text{ridge}}) = \frac{1}{(1+\lambda)^2}\,\operatorname{Var}(z) = \frac{\sigma^2}{(1+\lambda)^2}, \] which decreases as \(\lambda\) increases; large \(\lambda\) makes the estimator very low variance.
- Tradeoff and test MSE. The risk (per coordinate) is \[ \mathbb{E}\big[(\hat\theta^{\text{ridge}} - \beta)^2\big] = \underbrace{\frac{\lambda^2}{(1+\lambda)^2}\beta^2}_{\text{bias}^2} + \underbrace{\frac{\sigma^2}{(1+\lambda)^2}}_{\text{variance}}. \] As \(\lambda\) increases, variance drops faster than bias initially, so test MSE can decrease compared to OLS (\(\lambda=0\)). If \(\lambda\) is pushed too large, bias dominates and test error increases again: that is the classical bias–variance tradeoff.