The LASSO:
least absolute shrinkage and selection operator
Pierre Ablin
12/03/2019
Parietal Tutorial
Overview
- Linear regression
- The Lasso model
- Non-smooth optimization: proximal operators
Supervised learning
y \in \mathbb{R}, \enspace \mathbf{x} \in \mathbb{R}^P, \enspace \varepsilon = \mathcal{N}(0, \sigma^2), \enspace f: \mathbb{R}^P \rightarrow \mathbb{R}
y = f(\mathbf{x}) + \varepsilon
GOAL: Learn \( f \) from some realizations of \( (\mathbf{x}, y) \)
Linear regression
y = f(\mathbf{x}) + \varepsilon
GOAL: Learn \( f \) from some realizations of \( (\mathbf{x}, y) \)
Assumption:
\(f\) is linear
f(\mathbf{x}) = \beta_1x_1 + \beta_2 x_2 + \cdots + \beta_P x_P = \mathbf{x}^{\top} \mathbf{\beta}
GOAL: Learn \( \beta \) from some realizations of \( (\mathbf{x}, y) \)
Maximum likelihood?
y = \mathbf{x}^{\top} \mathbf{\beta} + \varepsilon
For one sample, if we observe \( \mathbf{x}_{\text{obs}}\) and \( y_{\text{obs}} \) , we have:
p(y = y_{\text{obs}} | \mathbf{x} = \mathbf{x}_{\text{obs}}) = p( \mathbf{x}_{\text{obs}}^{\top} \mathbf{\beta} + \varepsilon = y_{\text{obs}}| \mathbf{x} = \mathbf{x}_{\text{obs}})
i.e. :
p(y = y_{\text{obs}} | \mathbf{x} = \mathbf{x}_{\text{obs}}) = p(\varepsilon = y_{\text{obs}} - \mathbf{x}_{\text{obs}}^{\top} \mathbf{\beta})
\\
= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp{(- \frac{( y_{\text{obs}} - \mathbf{x}_{\text{obs}}^{\top} \mathbf{\beta})^2}{2\sigma^2})}
Maximum likelihood?
For one sample, if we observe \( \mathbf{x}_{\text{obs}}\) and \( y_{\text{obs}} \) , we have:
p(y = y_{\text{obs}} | \mathbf{x} = \mathbf{x}_{\text{obs}}) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp{(- \frac{( y_{\text{obs}} - \mathbf{x}_{\text{obs}}^{\top} \mathbf{\beta})^2}{2\sigma^2})}
If we observe a whole dataset of \(N\) i.i.d. samples$$ (\mathbf{x}_1,y_1), \cdots, (\mathbf{x}_N, y_N) \enspace, $$
the likelihood of this observation is:
\mathcal{L}(\beta) = p(y_1 = y_1, \cdots, y_N = y_N | \mathbf{x}_1 = \mathbf{x}_1, \cdots, \mathbf{x}_N=\mathbf{x}_N)
\\
= p(y_1 = y_1| \mathbf{x}_1 = \mathbf{x}_1) \times \cdots \times p(y_N = y_N| \mathbf{x}_N = \mathbf{x}_N)
\\
maximum likelihood estimator:
\beta_{\text{MLE}} = \arg\max \mathcal{L}(\beta) = \arg\min \sum_{n=1}^N \frac12 (y_n - \mathbf{x}_n^{\top} \beta)^2
Matrix interlude
Define the \(\ell_2\) norm of a vector:
$$|| \mathbf{v}|| = \sqrt{\sum_{n=1}^N v_n^2}$$
And note in condensed matrix form:
$$ \mathbf{y} = [y_1,\cdots, y_N]^{\top} \in \mathbb{R}^N$$
X = [\mathbf{x}_1,\cdots,\mathbf{x}_N]^{\top}\in\mathbb{R}^{N \times P}
Then,
\sum_{n=1}^N (y_n - \mathbf{x}_n^{\top} \beta)^2 = ||\mathbf{y} - X \beta||^2
M.L.E.
\beta_{\text{MLE}} = \arg\min \frac12 ||\mathbf{y} - X \beta||^2
Pros:
- Consistant: as \( N \rightarrow \infty \), \(\beta_{\text{MLE}} \rightarrow \beta^* \)
- Comes from a clear statistical framework
- Can behave badly when the model is mispecified (e.g. \(f\) is not really linear)
- Lots of variance in \(\beta_{\text{MLE}}\) when \(N\) is not big enough
- Ill posed problem when \(N < P\) !
Cons:
What do we do when n < p?
Bet on sparsity / Ockham's razor: only a few coefficients in \( \beta \) play a role
The \( \ell_0 \) Norm:
\ell_0(\beta) = \text{Number of non-zero coefficients in }\beta
Example:
(1.1, 2.5, -6.3, 4.1) \rightarrow 4
\\
(1.1, 0, 0, 4.1) \rightarrow 2
\\
(0, 0, 0, 0) \rightarrow 0
Subset selection/matching pursuit
Idea: select only a few active coefficients with the \( \ell_0 \) norm.
\beta_{\text{subset}} = \arg\min \frac12 || \mathbf{y} - X\beta ||^2 \enspace \text{subject to} \enspace \ell_0(\beta) \leq t
Where \( t\) is an integer controlling the sparsity of the solution.
- Non-convexity
- Instability: adding a sample may completely change \( \beta \) and its support.
- NP-Hard (loooong to solve)
Problems:
The Lasso
Idea: relax the \(\ell_0\) norm to obtain a convex problem
\beta_{\text{lasso}} = \arg\min \frac12 || \mathbf{y} - X\beta ||^2 \enspace \text{subject to} \enspace \ell_1(\beta) \leq t
Where \( t\) is a level controlling the sparsity of the solution.
See notebook
\ell_1(\beta) = |\beta |_1 = |\beta_1| + \cdots + |\beta_P|
The Lasso: Lagrange formulation
\beta_{\text{lasso}} = \arg\min \frac12 || \mathbf{y} - X\beta ||^2 + \lambda |\beta|_1
Where \( \lambda \) is a level controlling the amount of sparsity.
- Relationship between \(\lambda\) and \(t\) depends on the data
- It promotes sparsity: there is a threshold \(\lambda_{\text{max}}\) such that \( \lambda >\) \(\lambda_{\text{max}}\) implies \(\beta_{\text{lasso}} = 0 \)
Derivation
\(\beta = 0\) is a solution to the lasso if:
0 \in \partial_{\beta=0}(\frac12 ||y - X \beta||^2 + \lambda |\beta|_1)
Where \(\partial \) is the subgradient.
\partial(\frac12 || y - X\beta||^2) = \{X^{\top}X\beta - X^{\top} y\}
\partial_{\beta=0}(\lambda|\beta|_1) = [-\lambda, \lambda]^P
So the condition writes:
\forall p, \enspace (X^{\top}y)_p \in [-\lambda , \lambda]
\lambda_{\text{max}} = |X^{\top}y|_{\infty}
Lasso is useful because...
- It performs at the same time model selection and estimation : it tells you which coefficients in \(\mathbf{x}\) are important.
Boston dataset:
Lasso is useful because...
- Leveraging sparsity enables fast solvers
If \(p = 1000 \) and \(n = 1000 \), computing \(X\beta\) takes \(\sim n \times p = 10^6 \) operations
Intuition:
But if we know that there are only \(10 \) non-zero coefficients in \(\beta\), it takes only \(\sim 10^4\) operations
The same reasoning applies to most quantities useful to estimate the model
Lasso estimation
how do we fit the model?
Gradient descent 101
Minimize \( f(\mathbf{x})\) where \(f\) is differentiable.
Iterate \(\mathbf{x}^{t+1} = \mathbf{x}^t - \eta \nabla f(\mathbf{x}^t) \)
Equivalent to minimizing a quadratic surrogate:
\mathbf{x}^{t+1} = \arg\min_{\mathbf{x}} f(\mathbf{x}^t) + \nabla f(\mathbf{x}^t)^{\top}(\mathbf{x} - \mathbf{x}^t) + \frac{1}{2\eta} ||\mathbf{x} - \mathbf{x}^t||_2^2
The simplest lasso:
\(P=1, N = 1, X = 1:\) can we minimize \(\frac12 (y - \beta)^2 + \lambda |\beta|\) ?
Yes: proximity operator
\text{prox}_{|\cdot |}(y, \lambda) = \arg\min_{\beta} \frac12(y-\beta)^2 + \lambda |\beta|
Soft-thresholding:
\text{prox}_{|\cdot |}(y, \lambda) =
- \(0\) if \(|y|\leq \lambda\)
- \(y- \lambda\) if \(y>\lambda\)
- \(y + \lambda\) if \(y < -\lambda\)
Soft thresholding
\text{st}(y, \lambda) =
- \(0\) if \(|y|\leq \lambda\)
- \(y- \lambda\) if \(y>\lambda\)
- \(y + \lambda\) if \(y < -\lambda\)
Separability:
\text{prox}_{|\cdot |_1}(\mathbf{z}, \lambda) = \arg\min_{\beta}\frac12 ||\mathbf{z} - \beta ||^2 + \lambda |\beta|_1 =(\text{st}(z_1, \lambda), \cdots, \text{st}(z_P, \lambda))
ISTA: Iterative soft threshdolding algorithm
\beta_{\text{lasso}} = \arg\min \frac12 || \mathbf{y} - X\beta ||^2 + \lambda |\beta|_1
Gradient of the smooth term:
\nabla_{\beta}(\frac12||\mathbf{y} - X \beta ||^2) = X^{\top}( X \beta - \mathbf{y})
\beta^{t+1} = \text{prox}_{|\cdot|_1}(\beta^{t} - \eta X^{\top}(X\beta^{t} - \mathbf{y}), \lambda\eta)
ISTA:
Comes from:
\mathbf{\beta}^{t+1} = \arg\min_{\beta} f(\mathbf{\beta}^t) + \nabla f(\beta^t)^{\top}(\beta - \beta^t) + \frac{1}{2\eta} ||\beta - \beta^t||_2^2 + \lambda |\beta|_1
ISTA: Iterative soft threshdolding algorithm
\beta^{t+1} = \text{prox}_{|\cdot|_1}(\beta^{t} - \eta X^{\top}(X\beta^{t} - \mathbf{y}), \lambda\eta)
- Some coefficients are \(0\) because of the prox
- One iteration takes \(O(\min(N, P) \times P)\)
Can be problematic for large \(P\) !
COORDInate descent
\beta^{t+1}_j = \text{prox}_{|\cdot|_1}(\beta^{t}_j - \eta X[j, :]^{\top}(X\beta^{t} - \mathbf{y}), \lambda\eta)
Idea : update only one coefficient of \( \beta \) at each iteration.
- If we know the residual \( \mathbf{r} =X\beta^{t} - \mathbf{y}\), each update is \(O(n) \) :)
Residual update:
After updating the coordinate \(j\):
\mathbf{r} \leftarrow \mathbf{r} + (\beta^{\text{new}}_j - \beta^{\text{old}}_j)X[j, :]
\( \beta^{t+1}_i = \beta^{t}_i \) for \( i \enspace \ne j \)
Thanks!
The LASSO
By Pierre Ablin
