Regularized GLMs
Jonathan Taylor
Winter 2022
Penalized regression
- The objective function for a GLM (as a minimization) can be written as
\[
\Lambda(X\beta) - \beta'(X'Y) = - \log L(\beta|Y) = \ell(\beta|Y).
\]
- With some abuse of notation
\[
\Lambda(\alpha) = \sum_{i=1}^n \Lambda(\alpha_i)
\]
and \(\Lambda\) is the log MGF of \(Y\)’s as a function of the natural parameter.
Penalized regression
\[
\Lambda(X\beta) - \beta'(X'Y) + {\cal P}(\beta)
\]
Three examples we’ve seen:
Gaussian \(\Lambda(\alpha) = \alpha^2/2\) (or scale objective \(\sigma^2\));
Logistic \(\Lambda(\alpha) = \log(1 + e^{\alpha})\).
Poisson log-linear \(\Lambda(\alpha) = e^{\alpha}\).
Ridge regression
- Ridge regression adds a quadratic pnalty
\[
{\cal P}(\beta) = \frac{\lambda}{2} \|\beta\|_2^2 = \frac{\lambda}{2} \sum_{j=1}^p \beta_j^2.
\]
\[
\hat{\beta}_{\lambda} = \text{argmin}_{\beta} \Lambda(X\beta) -\beta'(X'Y) + \frac{\lambda}{2} \|\beta\|^2_2
\]
- Features are usually centered and scaled so coefficients are unitless.
“Bayesian” interpretation
- If we were to put a \(N(0, \lambda^{-1} I)\) prior on \(\beta\) its log posterior would be (up to normalizing constant depending only on \((X,Y)\))
\[
-(X\beta)'Y + \Lambda(X\beta) + \frac{\lambda}{2} \|\beta\|^2_2
\]
- Ridge estimator yields MAP (Maximum A Posteriori) estimator (i.e. the posterior mode).
Easy to add weights
\[
{\cal P}(\beta) = \frac{1}{2} \sum_{j=1}^p \lambda_j \beta_j^2
\]
unpenalized variables have \(\lambda_j=0\).
General quadratics also fine
\[
{\cal P}(\beta) = \frac{1}{2}\beta'Q\beta
\]
- Corresponds to \(N(0,Q^{-1})\) prior.
Overall objective
\[
\beta \mapsto \Lambda(X\beta) - \beta'(X'Y) + {\cal P}(\beta) = \ell(\beta) + \frac{1}{2} \beta'Q\beta
\]
Solving for the ridge estimator
\[
\begin{aligned}
\nabla {\cal P}(\beta) &= Q \beta \\
\nabla^2 {\cal P}(\beta) &= Q
\end{aligned}
\]
Newton-Raphson / Fisher scoring
\[
\hat{\beta}^{(t)} = \hat{\beta}^{(t-1)} - (X'W^{(t-1)}X + Q)^{-1}
\left(X'(E_{\hat{\beta}^{(t-1)}}(Y) - Y) + Q \hat{\beta}^{(t-1)}\right)
\]
\[
\hat{\beta}^{(t)} = \hat{\beta}^{(t-1)} - (X'W^{(t-1)}X + Q)^{-1}
\left(X'W^{(t-1)}g'(\hat{\mu}^{(t-1)})(\hat{\mu}^{(t-1)} - Y) + Q \hat{\beta}^{(t-1)}\right)
\]
- Not the only way to solve this problem…
Why do we regularize?
- Consider Gaussian case with \(Q=\lambda I\)
\[
\hat{\beta} = (X'X+\lambda I)^{-1}X'Y \sim N\left((X'X+\lambda I)^{-1}X'X\beta, \sigma^2 (X'X+\lambda I)^{-1}X'X(X'X+\lambda I)^{-1}\right)
\]
- Writing \(X=UDV'\) we see
\[
\begin{aligned}
(X'X+\lambda I)^{-1}X'X &= \sum_{i=1}^p \frac{d^2_i}{d^2_i+\lambda} V_iV_i' \\
(X'X+\lambda I)^{-1}X'X(X'X+\lambda I)^{-1} &= \sum_{i=1}^p \frac{d^2_i}{(d^2_i+\lambda)^2} V_iV_i' \\
\end{aligned}
\]
Bias
- Writing \(\beta=\sum_j \alpha_j V_j\):
\[
\|E[\hat{\beta}]-\beta\|^2_2 = \lambda^2 \sum_{i=1}^p \frac{\alpha_i^2}{(d_i^2 + \lambda)^2} \\
\]
Variance
\[
\text{Tr}(\text{Var}(\hat{\beta})) = \sum_{i=1}^p \frac{d^2_i}{(d^2_i+\lambda)^2}
\leq \sum_{i=1}^p \frac{1}{d_i^2}
\]
MSE
\[
MSE(\hat{\beta}) = E[\|\hat{\beta}-\beta\|^2_2] = \|E[\hat{\beta}]-\beta\|^2_2 + \text{Tr}(\text{Var}(\hat{\beta}))
\]
- Assuming \(\beta \neq 0\)
\[\frac{d}{d\lambda}\left(\|E[\hat{\beta}]-\beta\|^2_2 + \text{Tr}(\text{Var}(\hat{\beta}))\right)\biggl|_{\lambda=0} < 0
\]
Ridge solution path for SPAM data
## Loading required package: Matrix
## Loaded glmnet 4.0-2
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Cross-validation to choose \(\lambda\)
Re-solve at \(\lambda_{\min}\)
Ridge solution path for South African heart data
Cross-validation to choose \(\lambda\)
Re-solve at \(\lambda_{1SE}\)
LASSO
- The LASSO adds a penalty to the objective
\[
{\cal P}(\beta) = \lambda \|\beta\|_1 = \sum_{j=1}^p \lambda_j |\beta_j|.
\]
- The indices for which \(\lambda_j=0\) are unpenalized.
The LASSO picture
Why LASSO?
There are several sparsity inducing penalties out there: SCAD, MCP, etc.
But we know less about the “solutions” compared to LASSO because LASSO is a convex problem.
Solvable even in very high dimensions. (“same” for SCAD, MCP, etc.)
KKT conditions give an exact description of the solution however we decide to solve it (not same for SCAD, MCP, etc), at least under uniqueness.
Will the LASSO solution be unique?
In broad generality, yes. Uniqueness of the LASSO.
Solving the LASSO
Many algorithms out there. Many of you will end up taking EE364 (Convex Optimization) where you’ll see things like this in more detail.
We’ll consider two simple algorithms: coordinate descent, proximal gradient descent.
Everyone has their own way of talking about LASSO. These notes are just my way.
Coordinate descent
If the penalty \({\cal P}\) were smooth, then we could use something like Newton-Raphson as we did for fitting the GLM.
The penalty is differentiable everywhere except points where one of the \(\beta_j\)’s is 0: \(\implies\) not smooth.
However, the penalty is separable meaning
\[
{\cal P}(\beta) = \sum_{i=1}^p p_i(\beta_i).
\]
Block separability allows for a partition of \(\{1,\dots,p\}= \cup_{g \in G} g\) – group LASSO an example of this.
When the penalty is separable, the problem can be solved by successively solving low-dimensional problems.
Basic iteration
- Each iteration \(t\) chooses a coordinate \(j\) to update. Let the current position be \(\hat{\beta}^t\) define
\[
\hat{\beta}^{(t+1)}_j =
\text{argmin}_{\beta_j \in \mathbb{R}} \left[\beta_j \mapsto \Lambda (X_{-j}\hat{\beta}^{(t)}_{-j} + X_j \beta_j) - \beta_j(X_j'Y) + p_j(\beta_j) \right]
\]
Many space and time-saving tricks can be exploited here.
When solution is sparse, this can find the solution quickly!
A descent method: objective is monotone decreasing across iterations.
For convex loss and separable penalty this converges to a minimizer: Tseng
Refs: glmnet,
Coordinate descent and Gibbs sampler
Suppose we have a prior with independent coordinates with negative log density \({\cal P}(\beta)\).
The Gibbs sampler considers same “objective” that a coordinate descent algorithm does.
Instead of minimizing this function, it draws from a density proportional to its negative exponential.
Each iteration \(t\) chooses a coordinate \(j\) to update. Let the current position be \(\hat{\beta}^t\) draw \(\hat{\beta}^{(t+1)}_j\) from a density proportional to:
\[
\beta_j \mapsto \exp\left(-\left[\Lambda (X_{-j}\hat{\beta}^{(t)}_{-j} + X_j \beta_j) - \beta_j(X_j'Y) + p_j(\beta_j)\right] \right)
\]
Gibbs: sample
Coordinate Descent: optimize
Logistic regression with offset
\[
\beta \mapsto -(X\beta + Z)'Y + \sum_{i=1}^n \log(1 + \exp(X_i'\beta + Z_i))
\]
\[
X' \left(Y - \frac{\exp(X\widehat{\beta} + Z)}{1 + \exp(X\widehat{\beta} + Z)} \right) = 0
\]
- Generally speaking offset \(Z\) shifts introduces affine predictor \(\beta \mapsto X\beta + Z\).
Example with offset
## [,1]
## [1,] 4.493916e-09
## [2,] -5.790705e-08
- At iteration \(t\), for each \(j\), if
\[
\begin{aligned}
\lambda &> \left|\frac{\partial}{\partial \beta_j} \left[\Lambda(X_{-j} \widehat{\beta}_{-j} + X_j \beta_j) - \beta_j \cdot (X_j'Y)\right] \biggl|_{\beta_j=0} \right| \\
&= \left| X_j'\left( \pi(\hat{\beta}_{-j}, 0_j) - Y \right) \right|
\end{aligned}
\]
then \[
0 = \text{argmin}_{\beta_j \in \mathbb{R}} \left[\beta_j \mapsto \Lambda (X_{-j}\hat{\beta}^{(t)}_{-j} + X_j \beta_j) - \beta_j(X_j'Y) + p_j(\beta_j) \right]
\]
- Easy to check whether a coefficient should stay or be assigned to 0…
LASSO solution path for SPAM data
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
Cross-validation to choose \(\lambda\)
Re-solve at \(\lambda_{\min}\)
## [1] 53
Predictions at new X
## [1] 0.5748062 0.9800819 0.9999740 0.7488135 0.7487251
LASSO solution path for South African heart data
Cross-validation to choose \(\lambda\)
Re-solve at \(\lambda_{1SE}\)
## [1] 5
Proximal gradient descent
\[
\text{argmin}_{\beta} \frac{L}{2}\|Z-\beta\|^2_2 + {\cal P}(\beta)
\]
is soft-thresholding, most importantly it is fast
\[
\hat{\beta}_i = \max(|Z_i| - \lambda_i / L, 0) \cdot \text{sign}(Z_i)
\]
Basic update
- Given current state \(\hat{\beta}^{(t)}\) define
\[
Q^L(\beta;\hat{\beta}^{(t)}) = \ell(\hat{\beta}^{(t)}) + \nabla \ell(\hat{\beta}^{(t)})'(\beta - \hat{\beta}^{(t)}) + \frac{L}{2} \|\beta - \hat{\beta}^{(t)}\|^2_2 + {\cal P}(\beta)
\]
- The next iterate is defined as
\[
\hat{\beta}^{(t+1)} = \text{argmin}_{\beta} Q^L(\beta; \hat{\beta}^{(t)}) =
\text{argmin}_{\beta} \frac{L}{2} \|\beta - \hat{\beta}^{(t)} + L^{-1} \nabla \ell(\hat{\beta}^{(t)})\|^2_2 + {\cal P}(\beta).
\]
Quadratic approximation similar to Newton-Raphson.
If \(L\) is chosen so that \(L \cdot I \geq \nabla^2 \ell(\beta)\) for all \(\beta\) then this is an example of Majorization-Minimization (MM) method.
Stepsize \(L\) only needs to be chosen to be a descent – backtracking typically used.
There are several accelerated proximal gradient descent algorithms out there, some local to us: NESTA, TFOCS.
Proximal algorithm is key ingredient to all of these methods.
Majorization-Minimization
\[
\text{minimize}_{\beta} f(\beta)
\]
- We assume we have a function \(M(\beta;\beta')\) such that
\[
\begin{aligned}
Q(\beta;\beta) &= f(\beta) \\
Q(\beta;\beta') & \geq f(\beta) \qquad \forall \beta, \beta'
\end{aligned}
\]
\[
\hat{\beta}^{(t+1)} = \text{argmin}_{\beta} Q(\beta;\beta^{(t)})
\]
- Then, \(f(\hat{\beta}^{(t+1)}) \leq f(\hat{\beta}^{(t)})\). Yields a descent method!
Gradient descent and Langevin sampler
\[
\hat{\beta}^{(t+1)} = \hat{\beta}^{(t)} + \alpha \nabla \left[ \ell(\hat{\beta}^{(t)}) + {\cal P}(\hat{\beta}^{(t)}) \right] + \sqrt{2 \alpha} \cdot \xi^{(t)}
\]
with \(\xi^{(t)} \sim N(0,I)\).
Proximal map
\[
Z \mapsto \text{argmin}_{\beta} \frac{1}{2}\|\beta-Z\|^2_2 + {\cal P}(\beta)
\]
is called the proximal operator of \({\cal P}\).
Proximal gradient descent methods rely on \({\cal P}\) being fast.
Each step above evaluates the proximal map of \(L^{-1} \cdot {\cal P}\).
Reference: Parikh and Boyd
Solving the LASSO
\[
\text{minimize}_{\beta} \frac{1}{2} \|Y - X\beta\|^2_2
\]
subject to \({\cal P}(\beta) \leq C\).
- Proximal gradient descent replaces soft-thresholding with map
\[
Z \mapsto \text{argmin}_{\beta:\|\beta\|_1 \leq C} \frac{L}{2} \|\beta-Z\|^2_2.
\]
- Solving LASSO in this form will repeatedly project
\[
\hat{\beta}^{(t)} - \frac{1}{L} \nabla \ell(\hat{\beta}^{(t)}).
\]
onto the convex set \(\left\{\beta: \|\beta\|_1 \leq C \right\}.\)
KKT conditions
- How will we know we have solved
\[
\text{minimize}_{\beta} \ell(\beta) + {\cal P}(\beta) \ ?
\]
\[
\nabla \ell (\hat{\beta}) + \hat{u} = 0 \iff X'Y = X' \nabla \Lambda(X\hat{\beta}) + \hat{u}
\]
where
\[
\hat{u} \in \partial {\cal P}(\hat{\beta}).
\]
The set \(\partial {\cal P}(\beta)\) is called the subdifferential of \({\cal P}\) at \(\beta\).
If \({\cal P}\) is differentiable at \(\beta\), then the subdifferential is the derivative.
Can also deduce stationarity conditions from proximal gradient algorithm.
Stationarity from proximal gradient
- If \(\hat{\beta}\) is a statitionary point of proximal gradient with step size \(L^{-1}\), then
\[
\hat{\beta} = S_{\lambda L^{-1}}\left(\hat{\beta} - L^{-1} \nabla \ell(\hat{\beta})\right).
\]
- If \(\hat{\beta}_i=0\) then
\[
|L^{-1}\nabla \ell(\hat{\beta})_i| \leq \lambda L^{-1} \ \implies \ |\nabla \ell(\hat{\beta})_i| \leq \lambda
\]
- If \(\hat{\beta}_i \neq 0\) then
\[
\hat{\beta}_i - L^{-1} \nabla \ell(\hat{\beta})_i = \hat{\beta}_i + \text{sign}(\hat{\beta}_i) \cdot \lambda L^{-1}
\ \implies \ -\nabla\ell(\hat{\beta})_i = \lambda \cdot \text{sign}(\hat{\beta}_i)
\]
\[
{\cal P}(\beta) = \sup_{u: |u_i| \leq \lambda_i} u'\beta.
\]
- For penalties of the form
\[
{\cal P}(\beta) = \sup_{u \in K} u'\beta,
\]
the subdifferential is
\[
u \in \partial {\cal P}(\beta) \iff \beta \in N_u(K)
\]
with \(N_u(K)\) the set of normal directions at \(u \in K\).
\[
\left\{u: u'\beta = {\cal P}(\beta) \right\}.
\]
Normal cones of a cube
\[
\hat{u} \in \partial{\cal P}(\hat{\beta}) \iff \|u_{-E}/\lambda_{-E}\|_{\infty} \leq 1, u_E = \lambda_E s_E.
\]
- In other words: active signs are fixed, inactive score is not too big.
## [1] 1.322557
## [1] 1.573924 1.736504 0.000000
## [,1]
## [1,] -1.32255694
## [2,] -1.32255695
## [3,] 0.02511821
## [1] 1.322557
## [1] 1.57392444 1.73650397 0.02511821
Let’s perturb \(Y\) a little and see what happens.
## [1] 1.570243 1.739898 0.000000
## [,1]
## [1,] -1.32255592
## [2,] -1.32255695
## [3,] 0.03358971
## [1] 1.57024314 1.73989803 0.03358971
## [1] 1.57392444 1.73650397 0.02511821
Take-away
- The pair \((\hat{\beta}_E, u_{-E})\) parameterize the sufficient statistic \(X'Y\).
KKT conditions: deeper dive
Restricted problem
\[
\begin{aligned}
\nabla \ell(\hat{\beta})_E &= - \lambda_E s_E \\
\|\nabla \ell(\hat{\beta})_{-E}/\lambda_{-E}\|_{\infty} &\leq 1 \\
\end{aligned}
\]
- Having fixed \(E\) we can define restricted objective
\[
\beta_E \overset{\ell^E}{\mapsto} \ell(\beta_E, 0)
\]
\[
\nabla \ell^E(\hat{\beta}_E) = -\lambda_E s_E
\]
with constraint \(\text{sign}(\hat{\beta}_E) = s_E\).
Gaussian case
- Let’s take one Newton-Raphson step from \(\hat{\beta}_E\):
\[
\bar{\beta}_E = \hat{\beta}_E + (X_E'X_E)^{-1} \lambda_E s_E
\]
or
\[
\hat{\beta}_E = \bar{\beta}_E - (X_E'X_E)^{-1} \lambda_E s_E
\]
In the Gaussian case this reads \[
\text{sign}((X_E'X_E)^{-1}(X_E'Y - \lambda_E s_E)) = s_E.
\]
Allowing some \(\lambda_i\)’s to be zero changes things a little – only have conditions for \(\lambda_i > 0\).
Inactive block
- Let’s look at our other set of KKT equations.
\[
\begin{aligned}
\nabla \ell(\hat{\beta}_E, 0)_{-E} &\approx \nabla \ell(\bar{\beta}_E, 0)_{-E} + \nabla^2 \ell(\bar{\beta}_E,0)[-E,E](\hat{\beta}_E - \bar{\beta}_E) \\
& \approx \nabla \ell(\bar{\beta}_E, 0)_{-E} - \nabla^2 \ell(\bar{\beta}_E,0)_{-E} \nabla^2 \ell^E(\bar{\beta}_E)^{-1} \lambda_E s_E
\end{aligned}
\]
\[
\|\nabla \ell(\hat{\beta}_E, 0)_{-E}/\lambda_{-E}\|_{\infty} \leq 1.
\]
\[
\nabla^2 \ell(\bar{\beta}_E,0)_{-E} \nabla^2 \ell^E(\bar{\beta}_E)^{-1} \lambda_E s_E
\]
will behave as a fixed vector (with signs fixed).
Gaussian case
\[
\|(X_{-E}'(Y - P_EY) + X_{-E}'X_E(X_E'X_E)^{-1} \lambda_Es_E)/\lambda_{-E} \|_{\infty} \leq 1.
\]
Model selection consistency
Now suppose that model is correct with parameter \(\beta^*\) such that \(\beta^*_{-A}=0\).
When will the LASSO recover \(A\) with correct signs \(\text{sign}(\beta^*_A)\)? (For simplicity let’s now set all \(\lambda_i\)’s equal).
References (among others): Zhao and Yu, Candes and Plan
\[
\text{sign}(\bar{\beta}_A - \lambda \nabla \ell^2(\bar{\beta}_A) s^*_A) = s^*_A
\]
with \(s^*_A=\text{sign}(\beta^*_A)\).
- If \(\beta_A^*\) are far away enough from 0 (beta-min condition) then subtracting the second term will not change signs of \(\bar{\beta}_A\).
- Let’s look at the inactive block. We saw that if the active set is going to be \(A\) and signs \(s_A^*\) then
\[
\|(\nabla \ell(\bar{\beta}_A, 0)_{-A} - I_A \lambda s_A^*) / \lambda\|_{\infty} = \|(\nabla \ell(\bar{\beta}_A, 0)_{-A}/\lambda - I_A s_A^*) / \|_{\infty} \leq 1.
\]
- For this to happen, essentially we need the irrepresentable condition (Zhao and Yu)
\[
\max_{\|u_A\|_{\infty}\leq 1} \|I_A u_A\|_{\infty} \leq 1
\]
How large should \(\lambda\) be?
\[
\lambda \approx C \cdot E[\|\nabla \ell(\beta_A^*,0)\|_{\infty}], \qquad C > 1.
\]
- Such choice often guarantees LASSO is consistent. (Large literature on this, c.f. Negabhan et al.)