Lasso Regularization (L1 Regularization) › Matrix Notation

We can rewrite the LASSO optimization problem in matrix notation, which provides a more compact and mathematically elegant representation. This notation is particularly useful for understanding the relationship between LASSO and other optimization problems, and it's the form commonly used in machine learning libraries like scikit-learn.

Let's start with our familiar summation form and transform it step by step:

$$\min_{\beta} \left\{ \sum_{i=1}^n \left( y_i - \beta_0 - \sum_{j=1}^p \beta_j x_{ij} \right)^2 + \lambda \sum_{j=1}^p |\beta_j| \right\}$$

Step 1: Define the matrices and vectors

First, let's define our data in matrix form:

• $y$: An $n \times 1$ vector containing all target values:

$$y = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix}$$

• $X$: An $n \times (p+1)$ matrix containing all features (including a column of ones for the intercept):

$$X = \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{bmatrix}$$

• $\beta$: A $(p+1) \times 1$ vector containing all coefficients (including intercept):

$$\beta = \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_p \end{bmatrix}$$

Step 2: Express predictions in matrix form

The predicted values for all observations can be written as:

$$\hat{y} = X\beta$$

This single matrix multiplication $X\beta$ computes all predictions simultaneously. Row $i$ of $X$ multiplied by $\beta$ gives:

$$1 \cdot \beta_0 + x_{i1} \cdot \beta_1 + x_{i2} \cdot \beta_2 + \cdots + x_{ip} \cdot \beta_p = \beta_0 + \sum_{j=1}^p \beta_j x_{ij}$$

This is exactly the prediction for observation $i$:

$$\hat{y}_i = \beta_0 + \sum_{j=1}^p \beta_j x_{ij}$$

Step 3: Express residuals in matrix form

The residuals (errors) for all observations are:

$$r = y - X\beta = \begin{bmatrix} y_1 - \hat{y}_1 \\ y_2 - \hat{y}_2 \\ \vdots \\ y_n - \hat{y}_n \end{bmatrix}$$

Step 4: Express the sum of squared errors in matrix form

The sum of squared errors becomes:

$$\text{SSE} = \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n r_i^2 = r^T r = (y - X\beta)^T (y - X\beta)$$

Let's verify this step by step:

$$r^T r = \begin{bmatrix} r_1 & r_2 & \cdots & r_n \end{bmatrix} \begin{bmatrix} r_1 \\ r_2 \\ \vdots \\ r_n \end{bmatrix} = r_1^2 + r_2^2 + \cdots + r_n^2 = \sum_{i=1}^n r_i^2$$

Since $r_i = y_i - \hat{y}_i$, we have $\sum_{i=1}^n r_i^2 = \sum_{i=1}^n (y_i - \hat{y}_i)^2$, thus:

$$\sum_{i=1}^n r_i^2 = \sum_{i=1}^n (y_i - \hat{y}_i)^2$$

Step 5: Express the L1 penalty in matrix form

The L1 penalty term can be written using the L1 norm:

$$\lambda \sum_{j=1}^p |\beta_j| = \lambda \|\beta\|_1$$

where $\|\beta\|_1 = \sum_{j=1}^p |\beta_j|$ is the L1 norm of the coefficient vector. This means that the L1 penalty is the sum of the absolute values of the coefficients, or otherwise known as the "Manhattan distance" between the coefficients and the origin. It is the same as summing up all the features:

$$\|\beta\|_1 = \sum_{j=1}^p |\beta_j|$$

Final matrix form:

What we have now is a compact way of writing the LASSO optimization problem using matrix notation.

$$\min_{\beta} \left\{ (y - X\beta)^T (y - X\beta) + \lambda \|\beta\|_1 \right\}$$

Scikit-learn uses a slightly different but equivalent formulation that includes a normalization factor:

$$\min_{\beta} \left\{ \frac{1}{2n} \|y - X\beta\|_2^2 + \alpha \|\beta\|_1 \right\}$$

where:

• $\alpha$: regularization parameter in scikit-learn, related to traditional $\lambda$ by $\alpha = \frac{\lambda}{2n}$
• $n$: number of observations in the dataset
• $\|y - X\beta\|_2^2$: squared L2 (Euclidean) norm of residuals, equivalent to $\sum_{i=1}^n (y_i - \hat{y}_i)^2$ (the SSE)
• $\|\beta\|_1$: L1 (Manhattan) norm of coefficients, equal to $\sum_{j=1}^p |\beta_j|$ (sum of absolute values, intercept not penalized)

The expression $\|y - X\beta\|_2^2$ uses the L2 norm (also called the Euclidean norm):

$$\|y - X\beta\|_2 = \sqrt{\sum_{i=1}^n (y_i - \hat{y}_i)^2} = \sqrt{\sum_{i=1}^n (y_i - (X\beta)_i)^2}$$

where:

• $\|y - X\beta\|_2$: L2 norm of the residual vector (Euclidean distance)
• $y_i$: actual target value for observation $i$
• $\hat{y}_i$: predicted value for observation $i$, equal to $(X\beta)_i$
• $n$: number of observations

When we square the L2 norm, we get:

$$\|y - X\beta\|_2^2 = \sum_{i=1}^n (y_i - \hat{y}_i)^2 = (y - X\beta)^T (y - X\beta)$$

This is exactly our sum of squared errors (SSE)! The L2 norm squared gives us the same quantity as our matrix form $(y - X\beta)^T (y - X\beta)$. The L2 or Euclidean norm is simply the square root of the sum of the squares of the elements of the vector.

We use the L2 norm notation because it is mathematically equivalent to the matrix form $(y - X\beta)^T (y - X\beta)$, but more concise and widely recognized. The L2 norm is standard in optimization and statistics, making formulas easier to read and interpret. This notation also highlights the connection to Euclidean distance, which helps build intuition about how the loss function measures the overall error.

We now got the left side of the equation, let's look at the right side. The expression $\|\beta\|_1$ is the L1 norm of the coefficient vector:

$$\|\beta\|_1 = \sum_{j=1}^p |\beta_j| = |\beta_1| + |\beta_2| + \cdots + |\beta_p|$$

where:

• $\|\beta\|_1$: L1 norm (Manhattan distance) of the coefficient vector
• $|\beta_j|$: absolute value of the $j$-th coefficient
• $p$: number of features (excluding intercept)

This is exactly our L1 penalty term! The L1 norm is simply the sum of the absolute values of the elements of the vector.

Note: The intercept $\beta_0$ is typically not included in the L1 penalty, so the sum runs from $j=1$ to $j=p$, not from $j=0$ to $j=p$.

We know that the L1 norm offers several important advantages: it encourages sparsity by driving some coefficients exactly to zero, which enables automatic feature selection. Further, it has a clear geometric interpretation: it measures the "Manhattan distance" (the sum of absolute values) from the origin in coefficient space, and is a convex function, which ensures that the optimization problem remains tractable and that efficient algorithms can find a global minimum.

The Normalization Factor: $\frac{1}{2n}$Link Copied

Let's take a closer look at the normalization factor $\frac{1}{2n}$ that appears in front of the L2 norm in the scikit-learn formulation. In scikit-learn, the regularization strength is controlled by the parameter alpha ($\alpha$), which is related to the traditional $\lambda$ parameter by $\alpha = \frac{\lambda}{2n}$.

This means that the scikit-learn objective function is normalized by both the number of samples and a factor of 2, making the loss function consistent with the mean squared error and simplifying gradient calculations.

This formulation is sometimes referred to as the "normalized" or "scaled" LASSO objective, and is the standard approach used in most modern machine learning libraries.

The factor $\frac{1}{2n}$ serves several important purposes:

1. Sample size normalization: Dividing by $n$ makes the loss function independent of sample size. Without this, larger datasets would have larger loss values simply because they have more observations.

2. Gradient scaling: The factor $\frac{1}{2}$ simplifies the gradient calculation. Let's take the derivative of $\frac{1}{2n}\|y - X\beta\|_2^2$ with respect to $\beta$. To see how we arrive at this gradient, let's rewrite the squared L2 norm and take the derivative step by step:

   $$\frac{\partial}{\partial \beta} \left[ \frac{1}{2n}\|y - X\beta\|_2^2 \right] = \frac{\partial}{\partial \beta} \left[ \frac{1}{2n} (y - X\beta)^T (y - X\beta) \right]$$

   Let's see how we can rewrite the term $(y - X\beta)^T (y - X\beta)$ in a more manageable form. The term $(y - X\beta)^T (y - X\beta)$ is a quadratic form that represents the sum of squared errors in matrix notation. To make it easier to take derivatives and understand the structure, we can rewrite it using the distributive property of matrix multiplication:

   $$(y - X\beta)^T (y - X\beta) = y^T y - y^T X\beta - (X\beta)^T y + (X\beta)^T (X\beta)$$

   Notice that $(X\beta)^T y$ is a scalar and is equal to its transpose $y^T X\beta$, so the two middle terms combine to $-2\beta^T X^T y$. The last term, $(X\beta)^T (X\beta)$, can be rewritten as $\beta^T X^T X \beta$. Putting it all together, we get:

   $$(y - X\beta)^T (y - X\beta) = y^T y - 2\beta^T X^T y + \beta^T X^T X \beta$$

   Taking the derivative with respect to $\beta$:

   • The derivative of $y^T y$ with respect to $\beta$ is 0 (since it does not depend on $\beta$)
   • The derivative of $-2\beta^T X^T y$ with respect to $\beta$ is $-2 X^T y$
   • The derivative of $\beta^T X^T X \beta$ with respect to $\beta$ is $2 X^T X \beta$

   So, the derivative is:

   $$\frac{\partial}{\partial \beta} (y - X\beta)^T (y - X\beta) = -2 X^T y + 2 X^T X \beta = 2 X^T (X\beta - y) = -2 X^T (y - X\beta)$$

   where:

   • $X^T$: transpose of the feature matrix
   • $y$: target vector
   • $\beta$: coefficient vector

   Including the $\frac{1}{2n}$ factor:

   $$\frac{\partial}{\partial \beta} \left[ \frac{1}{2n}(y - X\beta)^T (y - X\beta) \right] = \frac{1}{2n} \cdot 2 X^T (X\beta - y) = \frac{1}{n} X^T (X\beta - y)$$

   Or, equivalently (by factoring out the negative sign):

   $$= -\frac{1}{n} X^T (y - X\beta)$$

   where:

   • $\frac{1}{n}$: normalization factor (average over all observations)
   • $X^T (X\beta - y)$: gradient of the squared error term with respect to $\beta$

   The gradient $\frac{1}{n} X^T (X\beta - y)$ points in the direction of steepest ascent of the loss function. In gradient descent for minimization, we move in the opposite direction by subtracting this gradient from the current $\beta$. Without the $\frac{1}{2}$ factor, we would have an extra factor of 2 in the gradient.

3. Consistency with statistical literature: This normalization makes the objective function consistent with the mean squared error (MSE), since $\frac{1}{n}\|y - X\beta\|_2^2$ is exactly the MSE.

The Regularization Parameter: $\alpha$Link Copied

In scikit-learn, the regularization parameter is called alpha ($\alpha$) instead of lambda ($\lambda$). The relationship between the two formulations is:

$$\alpha = \frac{\lambda}{2n}$$

where:

• $\alpha$: scikit-learn regularization parameter
• $\lambda$: traditional LASSO regularization parameter
• $n$: number of observations in the dataset

This means:

• Small $\alpha$: Less regularization, coefficients closer to OLS solution
• Large $\alpha$: More regularization, more coefficients driven to zero

Let's put it all together with a concrete example. Suppose we have:

• $n = 100$ observations
• $p = 5$ features (plus intercept)
• $\alpha = 0.1$

The scikit-learn objective function becomes:

$$\min_{\beta} \left\{ \frac{1}{200} \|y - X\beta\|_2^2 + 0.1 \|\beta\|_1 \right\}$$

where:

• $\frac{1}{200}$: normalization factor ($\frac{1}{2n}$ where $n = 100$)
• $0.1$: regularization parameter $\alpha$
• $\|y - X\beta\|_2^2$: squared L2 norm of residuals
• $\|\beta\|_1$: L1 norm of coefficients

Expanding this:

$$\min_{\beta} \left\{ \frac{1}{200} \sum_{i=1}^{100} (y_i - \hat{y}_i)^2 + 0.1 \sum_{j=1}^5 |\beta_j| \right\}$$

where:

• $i = 1, 2, \ldots, 100$: index over observations
• $j = 1, 2, \ldots, 5$: index over features
• $y_i$: actual target value for observation $i$
• $\hat{y}_i$: predicted value for observation $i$
• $\beta_j$: coefficient for feature $j$

This is equivalent to:

$$\min_{\beta} \left\{ \frac{1}{200} \sum_{i=1}^{100} \left( y_i - \beta_0 - \sum_{j=1}^5 \beta_j x_{ij} \right)^2 + 0.1 (|\beta_1| + |\beta_2| + |\beta_3| + |\beta_4| + |\beta_5|) \right\}$$

where:

• $\beta_0$: intercept term (not penalized)
• $x_{ij}$: value of feature $j$ for observation $i$
• $|\beta_j|$: absolute value of coefficient $j$

Matrix notation has several important practical benefits. First, matrix operations are highly optimized in modern numerical libraries. Second, expressing the objective in terms of norms directly connects LASSO to the broader field of convex optimization. Third, this standardized form is what most machine learning libraries, such as scikit-learn, use in their implementations. Finally, this framework is flexible: it extends naturally to other regularization techniques like Ridge regression and Elastic Net, which will be discussed in later sections.

By understanding the matrix formulation, you will be better equipped to interpret results from scikit-learn, select appropriate values for the alpha parameter, recognize the connections between different regularization methods, and even implement custom optimization algorithms if your application requires it.

After all of this math, let's take a look at a visualization that demonstrates how L1 regularization works by showing how coefficients change as the regularization parameter increases. This is called a "coefficient path" plot.

Recall that scikit-learn's Lasso solves the objective function $\tfrac{1}{2n}\,\lVert y - X\beta \rVert_2^2 + \alpha\,\lVert \beta \rVert_1$, where $\alpha = \tfrac{\lambda}{2n}$. This objective function defines what the algorithm is optimizing, and the choice of $\alpha$ directly controls the strength of regularization. In the plot below, we vary $\alpha$ to show its effect on the coefficients. The underlying implementation typically uses an algorithm such as the Iterative Shrinkage-Thresholding Algorithm (ISTA) or coordinate descent to efficiently solve this objective, but the objective function itself is what determines the solution.