Coordinate Descent for Lasso

Why Update One Coefficient at a Time?

The lasso estimates regression coefficients by balancing fit to the observed response against a penalty on coefficient size. With several predictors, the coefficients are coupled through the residual sum of squares, so minimizing over all coefficients at once can be inconvenient. Coordinate descent instead holds all but one coefficient fixed, minimizes the objective for the remaining coefficient, and repeats this process through the predictors.

For a fixed penalty, each update is a one-dimensional minimization. Its solution is a soft-thresholding rule: an association with the response must exceed the penalty before it produces a nonzero coefficient. The update also depends on the predictor's squared length, so the familiar soft-thresholding formula without a denominator applies directly when the predictors are standardized to a particular scale.

Definition: For a response vector $y$, predictor columns $x_1,\ldots,x_p$, and coefficient vector $\beta$, use the lasso objective $$ Q(\beta)=\frac{1}{2n}\left\|y-\sum_{k=1}^p x_k\beta_k\right\|^2+\lambda\sum_{k=1}^p|\beta_k|, $$ where $n$ is the number of observations and $\lambda\geq 0$. The intercept, if present, is not penalized. For coefficient $j$, its partial residual is $r_j=y-\sum_{k\ne j}x_k\beta_k$: it includes the response and all other predictors' current contributions, but excludes predictor $j$.

In this tutorial, a centered predictor is called standardized when $x_j^Tx_j/n=1$; the response is centered as well, so the intercept can be treated separately. This convention fixes the scale used in the update. The objective's factor $1/(2n)$ also fixes the meaning of $\lambda$: changing the objective's scaling changes the corresponding threshold.

Deriving the Coordinate Update

To update $\beta_j$, hold the other coefficients fixed. The part of the objective that depends on a candidate value $b$ for $\beta_j$ is $$ \frac{1}{2n}\|r_j-x_jb\|^2+\lambda|b|. $$ Expanding the squared norm gives a quadratic term, a linear term measuring the association between $x_j$ and its partial residual, and the absolute-value penalty.

\frac{1}{2n}\|r_j-x_jb\|^2+\lambda|b| =\text{constant}+\frac{z_j}{2}b^2-c_jb+\lambda|b|, \qquad z_j=\frac{x_j^Tx_j}{n},\quad c_j=\frac{x_j^Tr_j}{n}.

Assume $z_j>0$, so the predictor is not a zero column. For $b>0$, the derivative of the terms depending on $b$ is $z_jb-c_j+\lambda$, and its zero is $b=(c_j-\lambda)/z_j$, which is valid in this region only when $c_j>\lambda$. For $b<0$, the derivative is $z_jb-c_j-\lambda$; its zero is $b=(c_j+\lambda)/z_j$, valid when $c_j<-\lambda$. If $-\lambda\leq c_j\leq\lambda$, the objective is minimized at the kink $b=0$: the left and right slopes there are respectively $-c_j-\lambda\leq0$ and $-c_j+\lambda\geq0$. These cases yield the update below.

Coordinate-descent update: Define $S(c,\lambda)=\operatorname{sign}(c)(|c|-\lambda)_+$, where $(a)_+=\max(a,0)$. Then $$ \beta_j\leftarrow \frac{S(c_j,\lambda)}{z_j} =\frac{S\!\left(x_j^Tr_j/n,\lambda\right)}{x_j^Tx_j/n}. $$ When $x_j^Tx_j/n=1$, this reduces to $\beta_j\leftarrow S(c_j,\lambda)$. The coefficient is zero whenever $|c_j|\leq\lambda$.

This derivation explains why the rule uses a partial-residual association rather than an unadjusted association with $y$: the contributions of the other predictors must already be accounted for. It also explains why a nonstandardized predictor requires division by $z_j$. Each exact coordinate update minimizes the objective over that coordinate with the others held fixed, so it cannot increase the objective.

Worked Example: Deriving a One-Coordinate Update

A centered, standardized predictor has $x_j^Tx_j/n=1$. At the current values of all other coefficients, suppose $x_j^Tr_j/n=0.74$, and let $\lambda=0.29$. Find the updated coefficient and determine what a negative association of the same magnitude would produce.

Method. Use the partial-residual coordinate update, with $z_j=1$. The positive association exceeds the penalty, so $$ S(0.74,0.29)=\operatorname{sign}(0.74)(|0.74|-0.29)_+=1(0.45)=0.45. $$ Thus the updated coefficient is $\beta_j=0.45/1=0.45$.

If the association is $-0.74$, then $$ S(-0.74,0.29)=(-1)(0.74-0.29)=-0.45, $$ so the updated coefficient is $-0.45$. In contrast, if the association were $0.20$, its magnitude would be below $0.29$, and the update would be zero. The calculation avoids the mistake of shrinking every coefficient by a fixed amount while ignoring its sign or setting a coefficient to zero without comparing the association to the threshold.

Check. The positive and negative cases have equal magnitudes and opposite signs, and the smaller association is thresholded to zero, as required by the update.

Successive Updates and Partial Residuals

A complete coordinate-descent pass updates each coefficient in a specified order, using the most recent values of the other coefficients. Thus, after changing one coefficient, the partial residual for the next coefficient generally changes too. The method is iterative: later passes revisit earlier coefficients after their associated partial residuals have changed.

Worked Example: Two Successive Passes With Correlated Predictors

Consider four centered observations and two standardized predictors. Let $$ x_1=(1,1,-1,-1)^T,\quad x_3=(1,-1,1,-1)^T,\quad x_2=0.6x_1+0.8x_3=(1.4,-0.2,0.2,-1.4)^T. $$ The columns $x_1$ and $x_3$ are orthogonal, each has squared length $4$, and $x_2$ also has squared length $4$. Therefore $x_1^Tx_2/4=0.6$, and both fitted predictors $x_1,x_2$ are standardized under the convention $x_j^Tx_j/n=1$. Let $x_4=(1,-1,-1,1)^T$ and $y=1.2x_1+0.8x_3+0.5x_4=(2.5,-0.1,-0.9,-1.5)^T$. Use $\lambda=0.5$, start with $(\beta_1,\beta_2)=(0,0)$, and update predictor 1 before predictor 2.

Method. Since $n=4$, the raw associations are $$ \frac{x_1^Ty}{4}=1.2,\qquad \frac{x_2^Ty}{4}=0.6(1.2)+0.8(0.8)=1.36. $$ The initial partial residual for either predictor is $y$. The first update is therefore $$ \beta_1=S(1.2,0.5)=0.70. $$ Now form predictor 2's partial residual, $r_2=y-x_1\beta_1$. Its association is $$ \frac{x_2^Tr_2}{4} =\frac{x_2^Ty}{4}-\frac{x_2^Tx_1}{4}\beta_1 =1.36-(0.6)(0.70)=0.94. $$ Hence the first update of predictor 2 is $\beta_2=S(0.94,0.5)=0.44$. The coefficients after the first pass are $(0.70,0.44)$.

For the second pass, update predictor 1 using $r_1=y-x_2\beta_2$: $$ \frac{x_1^Tr_1}{4} =1.20-(0.6)(0.44)=0.936,\qquad \beta_1=S(0.936,0.5)=0.436. $$ Then update predictor 2 using the new value of $\beta_1$: $$ \frac{x_2^Tr_2}{4} =1.36-(0.6)(0.436)=1.0984,\qquad \beta_2=S(1.0984,0.5)=0.5984. $$ After two passes, the coefficients are $(0.436,0.5984)$. These are intermediate iterates, not a claim that the algorithm has converged.

Common mistake avoided. The value $1.36$ is the raw association $x_2^Ty/4$, before removing predictor 1's contribution. It is correct to subtract $0.6\beta_1$ from that value to obtain the partial-residual association. If the value had instead been stated as $x_2^Tr_2/4$, it would already include that adjustment, and subtracting $0.6\beta_1$ again would be wrong.

Check. Each update uses the other predictor's current coefficient, and the second pass uses the newly updated $\beta_1=0.436$. The resulting partial-residual associations, $0.936$ and $1.0984$, each exceed $0.5$, so both updates are positive.

Scaling, Residual Bookkeeping, and Stopping

Standardization makes the denominator $z_j$ equal to one, but the general update should be retained whenever the predictor columns do not have equal squared lengths. For efficient computation, an implementation may keep the full residual $e=y-\sum_kx_k\beta_k$. The partial residual for coordinate $j$ is then $r_j=e+x_j\beta_j$, because the full residual subtracts the current contribution of $x_j$. After updating $\beta_j$ to a new value, the full residual can be updated by subtracting $x_j$ times the coefficient change.

Worked Example: A Nonstandardized Predictor and an Inactive Coordinate

For one coordinate, suppose $n=5$, $x_j^Tx_j=20$, and the current partial residual satisfies $x_j^Tr_j=6.5$. For a second standardized coordinate, suppose its partial-residual association is $0.42$. Use $\lambda=0.30$ to find both updates.

Method. For the first coordinate, $$ z_j=\frac{20}{5}=4,\qquad c_j=\frac{6.5}{5}=1.30. $$ Its soft-thresholded association is $S(1.30,0.30)=1.00$. The coordinate update is therefore $$ \beta_j=\frac{1.00}{4}=0.25. $$ For the second coordinate, standardization gives $z=1$. Since $0.42>0.30$, its update is $S(0.42,0.30)=0.12$.

Common mistake avoided. Omitting the first coordinate's denominator would give $1.00$, four times the correct update. Another error would be applying the denominator to the second coordinate despite its stated standardization. The denominator is $x_j^Tx_j/n$ for the coordinate being updated; it is not automatically one for every column.

Check. The first coordinate's partial association divided by its squared-length factor is $1.30/4=0.325$ before thresholding, and the penalty reduces the numerator by $0.30$, giving $1.00/4=0.25$. The standardized coordinate has an association above the threshold and remains nonzero.

In practice, coordinate descent cycles through the coefficients until successive coefficient vectors or objective values change by less than a chosen tolerance. A small change in one coordinate alone is not enough to establish that the full cycle has stabilized. The penalty convention must also remain fixed: for the objective $Q$ used here, the threshold is $\lambda$; writing the residual term without the factor $1/(2n)$ would require a correspondingly rescaled penalty to represent the same optimization problem.

Common Exam Traps

Using the raw response association after other coefficients have been fitted. The update uses $x_j^Tr_j/n$, not $x_j^Ty/n$, unless all other coefficients are currently zero. The partial residual removes the other predictors' current contributions.
Subtracting an adjustment twice. If an association is already calculated using $r_j$, do not subtract the other predictors' contributions from it again.
Forgetting the predictor-length denominator. The update is $S(c_j,\lambda)/z_j$. The simpler formula applies only when $z_j=1$.
Using a different penalty scaling without changing the threshold. The threshold follows from the exact objective convention, including its factor $1/(2n)$.
Assuming an update is the final fit. Each coordinate update is conditional on the current values of the other coefficients. Further passes may change it.

Key takeaway: For each lasso coordinate, first construct the partial residual that excludes that predictor, then calculate its association and apply soft thresholding. Divide by $x_j^Tx_j/n$; omit this division only for a predictor standardized to squared length $n$.

Check Your Understanding

For each problem, use the objective scaling and standardization stated here. The final answer is provided for checking.

A standardized predictor has partial-residual association $0.82$, and $\lambda=0.30$. Find its coordinate update. (Answer: $0.52$.)
A standardized predictor has partial-residual association $-0.90$, and $\lambda=0.25$. Find its coordinate update, including its sign. (Answer: $-0.65$.)
Two standardized predictors have correlation $0.50$. Starting at zero, their raw associations with $y$ are $1.10$ and $1.00$, respectively. With $\lambda=0.40$, update predictor 1 and then predictor 2. What are the coefficients after this pass? (Answer: $(0.70,0.25)$.)
A standardized predictor's current full residual is $e$, its coefficient is $0.40$, and $x_j^Te/n=-0.12$. Find its partial-residual association and its update when $\lambda=0.30$. (Answer: partial-residual association $0.28$; updated coefficient $0$.)
A predictor has $x_j^Tx_j/n=2.5$, partial-residual association $1.40$, and $\lambda=0.40$. Compute its update and state why applying the standardized formula directly would be incorrect. (Answer: $0.40$; the standardized formula would give $1.00$ because it omits division by $2.5$.)

What You'll Learn