set.seed(1)
n <- 100
p <- 5
x <- matrix(rnorm(n * p), ncol = p)
y <- rnorm(n)
# standardize x variables
x_means <- colMeans(x)
x_sds <- apply(x, 2, function(x) sd(x) * sqrt((n-1) / n))
x_standardized <- (x - matrix(x_means, n, p, byrow = TRUE)) /
matrix(x_sds, n, p, byrow = TRUE)
# standardize y variables
y_mean <- mean(y)
y_sd <- sd(y) * sqrt((n-1) / n)
y_standardized <- (y - y_mean) / y_sdAppendix B — How we compute \(\lambda_{max}\)
Assume we are trying to build a linear model with the elastic net penalty. The elastic net model coefficients are defined as the solution to the following minimization problem:
\[ \min_{(\beta_0, \mathbf{\beta})} \quad \frac{1}{2n} \| \mathbf{y} - \beta_0 \mathbf{1} - \mathbf{X\beta} \|_2^2 + \lambda \left(\frac{1-\alpha}{2}\|\beta \|_2^2 + \alpha \|\beta \|_1 \right). \]
\(\lambda_{max}\) is defined to be the smallest value of the \(\lambda\) hyperparameter such that \(\mathbf{\beta} = \mathbf{0}\) is a solution to the problem above.
In practice, the computation of \(\lambda_{max}\) is quite involved and depends on many different function arguments. We shed some light on the computation here.
B.1 Standardized variables
First, assume that the response variable and input variables are all standardized to have mean 0 and variance 1. In that case, it can be shown that \(\hat\beta_0 = 0\), and so the minimization problem simplifies to
\[ \min_{\beta} \quad L(\beta) = \frac{1}{2n} \| \mathbf{y} - \mathbf{X}\beta \|_2^2 + \lambda \left(\frac{1-\alpha}{2}\|\beta \|_2^2 + \alpha \|\beta \|_1 \right). \]
The subgradient of the loss function \(L\) is
\[ \partial L = \frac{1}{n} \left(\mathbf{X}^\top \mathbf{X} \beta - \mathbf{X}^\top \mathbf{y} \right) + \lambda ((1-\alpha) \beta + \alpha \text{sgn}(\beta), \]
where \(\text{sgn}(\beta) \in \mathbb{R}^p\) is defined by
\[ (\text{sgn}(\beta))_j = \begin{cases} +1 &\text{if } \beta_j > 0, \\ -1 &\text{if } \beta_j < 0, \\ \in [-1, 1] &\text{if } \beta_j = 0. \end{cases} \]
A coefficient vector \(\hat\beta\) is an elastic net solution if and only if the zero vector \(\mathbf{0}\) belongs to the subgradient of \(L\) at \(\beta = \hat\beta\). Since \(\hat\beta = \mathbf{0}\) is a solution at \(\lambda = \lambda_{max}\), we have
\[ \begin{aligned} \mathbf{0} &\in \frac{1}{n} \left(\mathbf{X}^\top \mathbf{X} \mathbf{0} - \mathbf{X}^\top \mathbf{y} \right) + \lambda ((1-\alpha) \mathbf{0} + \alpha \text{sgn}( \mathbf{0}), \\ \mathbf{0} &\in -\frac{1}{n} \mathbf{X}^\top \mathbf{y} + \lambda \alpha \text{sgn}( \mathbf{0}), \\ - \frac{(\mathbf{X}^\top \mathbf{y})_j}{n} - \lambda \alpha &\leq 0 \leq - \frac{(\mathbf{X}^\top \mathbf{y})_j}{n} + \lambda \alpha \quad \text{for all } j = 1, \dots, p, \\ \lambda &\geq \frac{1}{\alpha} \max_j \left( \frac{|(\mathbf{X}^\top \mathbf{y})_j|}{n} \right). \end{aligned} \]
In the above, \((\mathbf{X}^\top \mathbf{y})_j \in \mathbb{R}^n\) refers to the \(j\)th element of \(\mathbf{X}^\top \mathbf{y}\). Thus, \(\lambda_{max} = \frac{1}{\alpha} \max_j \left( \frac{|(\mathbf{X}^\top \mathbf{y})_j|}{n} \right)\).
Let’s confirm this with a numerical example. First, generate some standardized data1:
We fit the data and pull out the \(\lambda_{max}\) value:
fit <- glmnet(x_standardized, y_standardized, standardize = FALSE)
fit$lambda[1]
#> [1] 0.1099471It matches what our math says it should be (by default, \(\alpha = 1\)):
max(abs(t(x_standardized) %*% y_standardized) / n)
#> [1] 0.1099471Let’s do the comparison again, but for \(\alpha = 0.2\):
alpha <- 0.2
fit <- glmnet(x_standardized, y_standardized, standardize = FALSE,
alpha = alpha)
fit$lambda[1]
#> [1] 0.5497355
max(abs(t(x_standardized) %*% y_standardized) / n) / alpha
#> [1] 0.5497355To avoid numerical issues arising from division by alpha, we set the minimum value of alpha to be \(10^{-3}\) for the purposes of computing \(\lambda_{max}\).
alpha <- 0
fit <- glmnet(x_standardized, y_standardized, standardize = FALSE,
alpha = alpha)
fit$lambda[1]
#> [1] 109.9471
max(abs(t(x_standardized) %*% y_standardized) / n) / (1e-3)
#> [1] 109.9471B.2 Unstandardized variables
Next, we assume that the response variable and input variables are not standardized to have mean 0 and variance 1, and glmnet is called with standardize = TRUE. In this case, glmnet standardized the response variable and input variables, then computes \(\lambda_{max}\) as in Section B.1. When reporting the lambda values in the returned object, the lambda sequence is returned to the original, unstandardized scale.
Again, we demonstrate this with an example. \(\lambda_{max}\) from the fit:
fit <- glmnet(x, y, standardize = TRUE)
fit$lambda[1]
#> [1] 0.1056963\(\lambda_{max}\) from the math:
# compute lambda_max on standardized scale, then return to original scale
max(abs(t(x_standardized) %*% y_standardized) / n) * y_sd
#> [1] 0.1056963B.3 Observation weights
Assume that we are fitting the elastic net for standardized variables with observation weights. In this case, the loss function is
\[ \min_{\beta} \quad L(\beta) = \frac{1}{2n} \sum_{i=1}^n w_i \left(y_i - \mathbf{x}_i^\top \beta \right)^2 + \lambda \left(\frac{1-\alpha}{2}\|\beta \|_2^2 + \alpha \|\beta \|_1 \right), \] where \(\mathbf{x}_j \in \mathbb{R}^p\) is the \(i\)th row of \(\mathbf{X}\) written as a column vector. If we define the diagonal matrix \(\mathbf{W} = \text{diag} (w_1, w_2, \dots, w_n)\), then we can rewrite the loss function as
\[ \min_{\beta} \quad L(\beta) = \frac{1}{2n} \| \mathbf{W}^{1/2}\mathbf{y} - \mathbf{W}^{1/2}\mathbf{X}\beta \|_2^2 + \lambda \left(\frac{1-\alpha}{2}\|\beta \|_2^2 + \alpha \|\beta \|_1 \right). \]
Performing the same computations as before, we find that
\[ \begin{aligned} \lambda_{max} &= \frac{1}{\alpha} \max_j \left( \frac{|((\mathbf{W}^{1/2}\mathbf{X})^\top (\mathbf{W}^{1/2}\mathbf{y}))_j|}{n} \right) \\ &= \frac{1}{\alpha} \max_j \left( \frac{|(\mathbf{X}^\top \mathbf{Wy})_j|}{n} \right). \end{aligned} \]
Let’s check with a numerical example. Recall that the presence of observation weights changes how the inputs to glmnet are standardized: it is the weighted mean and standard deviation that is scaled to 0 and 1 respectively. Let’s do the weighted standardization before comparing the \(\lambda_{max}\) values.
weights <- c(101, rep(1, times = n-1))
# get weighted feature means & standard deviations
x_w_means <- apply(x, 2, function(x) sum(weights * x) / sum(weights))
x_w_vars <- sapply(1:p, function(j) sum(weights * (x[, j] - x_w_means[j])^2) / sum(weights))
x_w_sds <- sqrt(x_w_vars)
# get transformed x matrix
x_w_std <- (x - matrix(x_w_means, n, p, byrow = TRUE)) /
matrix(x_w_sds, n, p, byrow = TRUE)
# transform y
y_w_mean <- sum(weights * y) / sum(weights)
y_w_sd <- sqrt(sum(weights * (y - y_w_mean)^2) / sum(weights))
y_w_std <- (y - y_w_mean) / y_w_sdHere is the \(\lambda_{max}\) from the fit:
fit <- glmnet(x_w_std, y_w_std, weights = weights, standardize = FALSE)
fit$lambda[1]
#> [1] 0.1395357And here is \(\lambda_{max}\) from the math:
max(abs(t(x_w_std) %*% (diag(weights) / sum(weights) * n) %*% y_w_std) / n)
#> [1] 0.1395357glmnetuses the 1/n rather than the 1/(n-1) formula for estimating the variance of the features and response variable.↩︎