Week 6: Regularization for Model Selection

DSAN 5300: Statistical Learning
Spring 2025, Georgetown University

Author
Affiliation

Jeff Jacobs

Published

Tuesday, February 18, 2025

Open slides in new tab →

Schedule

Today’s Planned Schedule:

Start End Topic
Lecture 6:30pm 7:00pm Extended Recap / Clarification →
7:00pm 7:20pm “Manual” Model Selection: Subsets →
7:20pm 7:40pm Key Regularization Building Block: Lp Norm →
7:40pm 8:00pm Regularized Regression Intro →
Break! 8:00pm 8:10pm
8:10pm 8:50pm Basically Lasso is the Coolest Thing Ever →
8:50pm 9:00pm Scary-Looking But Actually-Fun W07 Preview →

Roadmap

  • [Week 4] Oh no! When we go beyond linear models, we have to worry about overfitting!
    • New goal! Maximize generalizability rather than accuracy
    • Evaluate models on unseen test data rather than training data
  • [Week 5] Cross-Validation (CV) as a tool for Model Assessment: For more complex, non-linear models, is there some way we can try to… “foresee” how well a trained model will generalize?
    • Answer: Yes! Cross-validation!
  • [Week 6] Regularization as a tool for Model Selection: Now that we have a method (CV) for imperfectly measuring “generalizability”, is there some way we can try to… allow models to optimize CV but penalize them for unnecessary complexity?
    • Answer: Yes! Regularization methods like LASSO and Elastic Net!

[Reminder (W04)] New Goal: Generalizability

Goal 2.0: Statistical Learning

Find…

  • A function y^=f(x)
  • That best predicts Y for given values of X
  • For data that has not yet been observed! 😳❓

Clarification: Target Diagrams

Low Variance High Variance
Low Bias
High Bias
Figure 1: Adapted from Fortmann-Roe (2012), “Understanding the Bias-Variance Tradeoff”

Why Was This Helpful for 5100?

  • Law of Large Numbers:
    • Avg(many sample means s) true mean μ
  • θ^ unbiased estimator for θ:
    • Avg(Estimates θ^) true θ

The Low Bias, High Variance case

Relevance for CV Error

  • In Goal 2.0 world, we choose models on the basis of estimated test error (before, with Goal 1.0, we only used e.g. MSE, RSS, R2, which was fine for linear regression)
  • Data D = single realization of DGP (for 5300, only relevance is why we don’t look at test set)
  • [DTrain|DTest] = random permutation of D
  • Bullseye on target = true test error
    (We could compute this, but then we’d have to end the study, collect more data… better alternative on next slide!)
  • Darts thrown around bullseye = estimated test errors (CV fold errors!)
    • They don’t hit bullseye because we’re inferring DGP from from sample
  • True test error = f(D)=f([DTrain|DTest])
  • Validation error = f(DTrain)=f([DSubTr|DVal]),
    • Validation error is an estimate, using a smaller sample DTrain drawn from the same distribution (DGP) as true test error!

True Test Error vs. CV Error

Note the icons! Test set = Lake monster: pulling out of water to evaluate kills it 😵

True Test Error εTest=ErrTest

Data D “arises” out of (unobservable) DGP

Randomly chop D into [DTrainDTest]

ErrTestTest error,no cap=f(DTrainfitMθevalDTestThis kills monster 😢)

Issue: can only be evaluated once, ever 😱

Validation Set Error εVal=ε^Test=Err^Test

DGPD; Randomly chop into [DTrainDTest]

Leave DTest alone until end of study

Randomly chop DTrain into [DSubTrDVal]

Err^TestTest error,capping a bit=f(DSubTrfitMθevalDValMonster still alive!)

K-Fold Cross-Validation Error ε(K)=ε^Test=E^rrTest

DGPD; Randomly chop into [DTrainDTest]; Leave DTest for end of study

Randomly chop DTrain into [DTrFold(1)|DTrFold(2)||DTrFold(K)]

For i{1,,K}:

     εValFold(i)=f(DTrFold(i)fitMθevalDTrFold(i))

E^rrTestTest error,less cap!=1Ki=1KεValFold(i)  (…monster still alive, even after all that!)

General Issue with CV: It’s… Halfway There

CV plots will often look like (complexity on x-axis and CV error on y-axis):

Code
library(tidyverse) |> suppressPackageStartupMessages()
library(latex2exp) |> suppressPackageStartupMessages()
cpl_label <- TeX("$M_0$")
sim1k_delta_df <- tibble(
    complexity=1:7,
    cv_err=c(8, 2, 1, 1, 1, 1, 2),
    label=c("","",TeX("$M_3$"),"","",TeX("$M_6$"),"")
)
sim1k_delta_df |> ggplot(aes(x=complexity, y=cv_err, label=label)) +
  geom_line(linewidth=1) +
  geom_point(size=(2/3)*g_pointsize) +
  geom_text(vjust=-0.7, size=10, parse=TRUE) +
  scale_x_continuous(
    breaks=seq(from=1,to=7,by=1)
  ) +
  theme_dsan(base_size=22) +
  labs(
    title="Generic CV Error Plot",
    x = "Complexity",
    y = "CV Error"
  )

  • We “know” M3 preferable to M6 (same error yet, less overfitting) “1SE rule”
  • But… heuristic optimal! What are we gaining/losing as we move M6M3?
  • Enter REGULARIZATION!

CV Now Goes Into Your Toolbox

(We will take it back out later, I promise!)

Model Selection

  • [Week 5 ✅] We have a metric (CV error) for evaluating different models w.r.t. Goal 2.0…
  • [Week 6 so far] It gets us halfway to what we want, by showing us a “basin” of models with low CV error
  • [Now] Statistically-principled approach to overfitting that
    • Shows us why the 1SE rule “works”
    • Quantifies tradeoff: if I only have enough data to estimate βj for J0<J feats, which J0 should I choose?

Optimal but Infeasible: Best Subset Selection

Algorithm: Best Subset Selection

Let M0 be null model: Predicts y^(xi)=y for any xi

For k=1,2,,J:

     Fit all (Jk) possible models with k predictors, Mk is model with lowest RSS

Choose from M0,,MJ using CV or heuristics: AIC, BIC, adjusted R2

Feasible but Suboptimal: Stepwise Selection

Algorithm: Forward Stepwise Selection

Let M0 be null model: Predicts y^(xi)=y for any xi

For k=0,1,,J1:

     Fit Jk models, each adds single feature to Mk; call “best” model Mk+1

Choose from M0,,MJ using CV error (or heuristics: AIC, BIC, adjusted R2)

Algorithm: Backward Stepwise Selection

Let MJ be full model: Contains all features

For k=J,J1,,1:

     Fit k models, each removes single feature from Mk; call “best” model Mk1

Choose from M0,,MJ using CV error (or heuristics: AIC, BIC, adjusted R2)

Stepwise Selection Algorithms are Greedy

  • Like a mouse who chases closest cheese path with most cheese
  • Can get “trapped” in sub-optimal model, if (e.g.) feature is in M4 but not in M1,M2,M3!

k Mk (Best Subset) Mk (Forward Stepwise)
1 rating rating
2 rating, income rating, income
3 rating, income, student rating, income, student
4 cards, income, student, limit rating, income, student, limit

Regularization for Automatic Model Selection

  • Ridge Regression
  • Lasso
  • Elastic Net

Key Building Block: Lp Norms

  • Technically you have seen these before: distance metrics!

sim(u,v)vupLp norm of uv=(i=1n|viui|p)1/p

  • Euclidean distance is L2 norm: if u=(0,0) and v=(3,4),

vu2=(|30|2+|40|2)1/2=251/2=25=5

  • We’ll use even simpler form: distance of coefficients β from zero (so, u=0)

βp=(j=1J|βj|p)1/p

L1 and L2 Norms

  • L1 norm has a nice closed form expression as a sum (efficient, vectorizable!):

β1=(j=1J|βj|1)1/1=j=1J|βj|

  • L2 norm almost a similarly nice expression, besides this zigzag line thing ( )

β2=(j=1J|βj|2)1/2=j=1Jβj2🤨

  • Can always convert bound on true Euclidean distance like β210 into bound on squared Euclidean distance like β22100. Hence we’ll use squared Lp norm:

β22=((j=1J|βj|p)1/2)2=j=1Jβj2💆

Different Norms Different Distances from 0

  • Unit Disk in L2: All points v=(vx,vy)R2 such that

v21v221

  • Unit Disk in L1: All points v=(vx,vy)R2 such that

v11

Regularized Regression (Finally!)

General Form:

βreg=argminβ[1Ni=1N(y^i(β)yi)2MSE from before+λPenalty forβ22Dist from 0]

Three Main Types of Regularized Regression

Ridge Regression:

βridge=argminβ[1Ni=1N(y^i(β)yi)2+λβ22]

LASSO:

βlasso=argminβ[1Ni=1N(y^i(β)yi)2+λβ1]

Elastic Net:

βEN=argminβ[1Ni=1N(y^i(β)yi)2+λ2β22+λ1β1]

(Does anyone recognize λ from Lagrange multipliers?)

Top Secret Equivalent Forms

Ridge Regression:

βridge=arg{minβ[1Ni=1N(y^i(β)yi)2] subject to β22s}

LASSO:

βlasso=arg{minβ[1Ni=1N(y^i(β)yi)2] subject to β1s}

Elastic Net:

βEN=arg{minβ[1Ni=1N(y^i(β)yi)2] subject to β22s2β1s1}

The Key Plot

Code
library(tidyverse) |> suppressPackageStartupMessages()
library(latex2exp) |> suppressPackageStartupMessages()
library(ggforce) |> suppressPackageStartupMessages()
library(patchwork) |> suppressPackageStartupMessages()
# Bounding the space
xbound <- c(-1, 1)
ybound <- c(0, 1.65)
stepsize <- 0.05
dx <- 0.605
dy <- 1.6
# The actual function we're plotting contours for
b_inter <- 1.5
my_f <- function(x,y) 8^(b_inter*(x-dx)*(y-dy) - (x-dx)^2 - (y-dy)^2)
x_vals <- seq(from=xbound[1], to=xbound[2], by=stepsize)
y_vals <- seq(from=ybound[1], to=ybound[2], by=stepsize)
data_df <- expand_grid(x=x_vals, y=y_vals)
data_df <- data_df |> mutate(
  z = my_f(x, y)
)
# Optimal beta df
beta_opt_df <- tibble(
  x=121/200, y=8/5, label=c(TeX("$\\beta^*_{OLS}$"))
)
# Ridge optimal beta
ridge_opt_df <- tibble(
  x=0.111, y=0.998, label=c(TeX("$\\beta^*_{ridge}$"))
)
# Lasso diamond
lasso_df <- tibble(x=c(1,0,-1,0,1), y=c(0,1,0,-1,0), z=c(1,1,1,1,1))
lasso_opt_df <- tibble(x=0, y=1, label=c(TeX("$\\beta^*_{lasso}$")))

# And plot
base_plot <- ggplot() +
  geom_contour_filled(
    data=data_df, aes(x=x, y=y, z=z),
    alpha=0.8, binwidth = 0.04, color='black', linewidth=0.65
  ) +
  # y-axis
  geom_segment(aes(x=0, xend=0, y=-Inf, yend=Inf), color='white', linewidth=0.5, linetype="solid") +
  # Unconstrained optimal beta
  geom_point(data=beta_opt_df, aes(x=x, y=y), size=2) +
  geom_label(
    data=beta_opt_df, aes(x=x, y=y, label=label),
    hjust=-0.45, vjust=0.65, parse=TRUE, alpha=0.9
  ) +
  scale_fill_viridis_d(option="C") +
  #coord_equal() +
  labs(
    #title = "Model Selection: Ridge vs. Lasso Constraints",
    x = TeX("$\\beta_1$"),
    y = TeX("$\\beta_2$")
  )
ridge_plot <- base_plot +
  geom_circle(
    aes(x0=0, y0=0, r=1, alpha=I(0.1), linetype="circ", color='circ'), fill=NA, linewidth=0.5
  )
  # geom_point(
  #   data=data.frame(x=0, y=0), aes(x=x, y=y),
  #   shape=21, size=135.8, color='white', stroke=1.2, linestyle="dashed"
  # )
lasso_plot <- ridge_plot +
  geom_polygon(
    data=lasso_df, aes(x=x, y=y, linetype="diamond", color="diamond"),
    fill='white',
    alpha=0.5,
    linewidth=1
  ) +
  # Ridge beta
  geom_point(data=ridge_opt_df, aes(x=x, y=y), size=2) +
  geom_label(
    data=ridge_opt_df, aes(x=x, y=y, label=label),
    hjust=2, vjust=-0.15, parse=TRUE, alpha=0.9
  ) +
  # Lasso beta
  geom_point(data=lasso_opt_df, aes(x=x, y=y), size=2) +
  geom_label(
    data=lasso_opt_df, aes(x=x, y=y, label=label),
    hjust=-0.75, vjust=-0.15, parse=TRUE, alpha=0.9
  ) +
  ylim(ybound[1], ybound[2]) +
  # xlim(xbound[1], xbound[2]) +
  scale_linetype_manual("Line", values=c("diamond"="solid", "circ"="dashed"), labels=c("a","b")) +
  scale_color_manual("Color", values=c("diamond"="white", "circ"="white"), labels=c("c","d")) +
  # scale_fill_manual("Test") +
  # x-axis
  geom_segment(aes(x=-Inf, xend=Inf, y=0, yend=0), color='white') +
  theme_dsan(base_size=16) +
  coord_fixed() +
  theme(
    legend.position = "none",
    axis.line = element_blank(),
    axis.ticks = element_blank()
  )
lasso_plot

Bayesian Interpretation

  • Belief A: Most/all of the features you included have important effect on Y
  • A Gaussian prior on βj, μ=0
  • If XN(μ,σ2), pdf of X is

fX(x)=12πσ2exp[12(xμσ)2]

Code
library(tidyverse) |> suppressPackageStartupMessages()
library(latex2exp) |> suppressPackageStartupMessages()
prior_labs <- labs(
  x = TeX("$\\beta_j$"),
  y = TeX("$f(\\beta_j)$")
)
ggplot() +
  stat_function(fun=dnorm, linewidth=1) +
  xlim(-3, 3) +
  theme_dsan(base_size=28) +
  prior_labs

  • Gaussian prior Ridge Regression!
    (High complexity penalty λ low σ2)
  • Belief B: Only a few of the features you included have important effect on Y
  • B Laplacian prior on βj, μ=0
  • If XL(μ,b), pdf of X is

f(x)=12bexp[|xμb|]

Code
library(tidyverse) |> suppressPackageStartupMessages()
library(latex2exp) |> suppressPackageStartupMessages()
library(extraDistr) |> suppressPackageStartupMessages()
ggplot() +
  stat_function(fun=dlaplace, linewidth=1) +
  xlim(-3, 3) +
  theme_dsan(base_size=28) +
  prior_labs

  • Laplacian prior Lasso!
    (High complexity penalty λ low b)

Ok So… How Do We Find These Magical λ Values?

  • Cross-Validation!
  • (…that’s, uh, that’s it! That’s the whole slide!)
  • (But let’s look at what we find when we use CV…)

Varying λ Tradeoff Curve!

Ridge Regression Case:

Increasing λ in minβ[RSS+λβ22]
Decreasing s in minβ[RSS] s.t. β22s

Only Change: L1 instead of L2 Norm 🤯

Lasso Case:

Increasing λ in minβ[RSS+λβ1]
Decreasing s in minβ[RSS] s.t. β1s

Week 7 Preview: Linear Functions

  • You’ve already seen polynomial regression:

Y=β0+β1X+β2X2+β3X3++βdXd

Y=β0+β1cos(πX)+β2sin(πX)++β2d1cos(πdX)+β2dsin(πdX)

New Regression Just Dropped

Piecewise regression:

Choose K cutpoints c1,,cK

Let Ck(X)=1[ck1X<ck], (c0, cK+1)

Y=β0+β1C1(X)+β2C2(X)++βKCK(X)

Decomposing Fancy Regressions into Core “Pieces”

  • Q: What do all these types of regression have in common?

  • A: They can all be written in the form

    Y=β0+β1b1(X)+β2b2(X)++βdbd(X)

  • Where b() is called a basis function

    • Linear (d=1): b1(X)=X
    • Polynomial: bj(X)=Xj
    • Piecewise: bj(X)=1[cj1X<cj]