Soft-Thresholding Updates for Lasso GLMs

From a GLM Update to a One-Coordinate Problem

In a lasso GLM, a coefficient update must balance two forces: fitting the data and paying a penalty for nonzero coefficients. The key simplification is that, during a weighted least-squares step, all coefficients except one can be held fixed. The remaining one-dimensional minimization has a simple solution: a weighted partial-residual correlation is pulled toward zero by the penalty, and is set exactly to zero if it is not large enough.

Consider a GLM with linear predictor \(\eta_i=\beta_0+\sum_{j=1}^p x_{ij}\beta_j+o_i\), where \(o_i\) is a known offset. The intercept \(\beta_0\) is unpenalized, and the lasso penalty is \(\lambda\sum_{j=1}^p|\beta_j|\), with \(\lambda\geq 0\). Fisher scoring or IRLS supplies a working response \(z_i\) and positive working weight \(w_i\). With these held fixed, the penalized weighted least-squares surrogate is

Q(\beta_0,\beta)=\frac{1}{2n}\sum_{i=1}^n w_i \left(z_i-\beta_0-o_i-\sum_{k=1}^p x_{ik}\beta_k\right)^2 +\lambda\sum_{k=1}^p|\beta_k|.

The factor \(1/(2n)\) fixes the scaling used throughout this tutorial: the penalty \(\lambda\) is on the same scale as this average quadratic loss. The working response and weights are fixed only for the current weighted least-squares subproblem. A coordinate update minimizes this surrogate over one coefficient while holding the other coordinates fixed; it does not, in general, minimize the original GLM objective exactly.

Definition: For coordinate \(j\), define the partial residual excluding that coordinate by \(r_{ij}=z_i-\beta_0-o_i-\sum_{k\ne j}x_{ik}\beta_k\). Define its weighted partial-residual correlation and weighted squared norm by \(c_j=\frac{1}{n}\sum_{i=1}^n w_i x_{ij}r_{ij}\) and \(a_j=\frac{1}{n}\sum_{i=1}^n w_i x_{ij}^2\). When \(a_j>0\), the coordinate update is \(\beta_j^{\mathrm{new}}=S(c_j,\lambda)/a_j\), where \(S(c,\lambda)=\operatorname{sign}(c)(|c|-\lambda)_+\) and \((u)_+=\max(u,0)\).

Why Soft-Thresholding Is the Coordinate Minimizer

Holding the other terms fixed, the part of \(Q\) that depends on a candidate value \(b\) for \(\beta_j\) is

q_j(b)=\frac{1}{2n}\sum_{i=1}^n w_i(r_{ij}-x_{ij}b)^2+\lambda|b|.

Expanding the square and dropping terms that do not depend on \(b\) gives \(q_j(b)=\frac{a_j}{2}b^2-c_jb+\lambda|b|+\text{constant}\). This is a convex function when \(a_j>0\). Its minimum can be derived on the two half-lines \(b>0\) and \(b<0\), while also checking the possible minimum at \(b=0\).

Theorem: Suppose the working weights satisfy \(w_i>0\) and \(a_j>0\). The unique minimizer of \(q_j\) is \(S(c_j,\lambda)/a_j\). In particular, it equals zero exactly when \(|c_j|\leq\lambda\).

Proof. For \(b>0\), the derivative is \(a_jb-c_j+\lambda\). A stationary point in this region is \(b=(c_j-\lambda)/a_j\), which is positive exactly when \(c_j>\lambda\). For \(b<0\), the derivative is \(a_jb-c_j-\lambda\). Its stationary point is \(b=(c_j+\lambda)/a_j\), which is negative exactly when \(c_j<-\lambda\).

At zero, the subdifferential of \(|b|\) is the interval \([-1,1]\). Thus zero minimizes \(q_j\) exactly when \(0\in -c_j+\lambda[-1,1]\), or, equivalently, when \(|c_j|\leq\lambda\). If \(c_j>\lambda\), the positive stationary point is the minimizer; if \(c_j<-\lambda\), the negative stationary point is the minimizer; and otherwise the minimizer is zero. Since \(a_j>0\), the quadratic term is strictly convex, so this minimizer is unique. Combining the cases gives \(S(c_j,\lambda)/a_j\). \(\square\)

Worked Example: A Positive Update From Weighted Data

For one coordinate, suppose \(n=2\), the predictor values are \(x_{1j}=x_{2j}=1\), the working weights are \(w_1=1\) and \(w_2=3\), and the partial residuals are \(r_{1j}=1\) and \(r_{2j}=0.5\). Let \(\lambda=0.25\). Find the coordinate update.

First calculate the weighted squared norm:

a_j=\frac{1}{2}\{1(1)^2+3(1)^2\}=2.

The weighted partial-residual correlation is

c_j=\frac{1}{2}\{1(1)(1)+3(1)(0.5)\} =\frac{2.5}{2}=1.25.

Because \(c_j>\lambda\), soft-thresholding subtracts \(\lambda\) from the positive correlation. Therefore

\beta_j^{\mathrm{new}}=\frac{S(1.25,0.25)}{2} =\frac{1.25-0.25}{2}=0.5.

The update is positive, but smaller than the unpenalized weighted least-squares value \(c_j/a_j=1.25/2=0.625\). The penalty reduces the coefficient by \(0.125\).

Worked Example: When the Update Is Exactly Zero

Suppose \(n=2\), \(x_{1j}=1\), \(x_{2j}=2\), both weights are 1, and the partial residuals are \(r_{1j}=0.2\), \(r_{2j}=0.1\). Let \(\lambda=0.2\).

The weighted squared norm and correlation are

a_j=\frac{1}{2}(1^2+2^2)=2.5, \qquad c_j=\frac{1}{2}\{1(0.2)+2(0.1)\}=0.2.

Here \(|c_j|=\lambda\), so \(S(c_j,\lambda)=0\) and the update is \(\beta_j^{\mathrm{new}}=0/2.5=0\). Equality at the threshold still gives an exact zero. The coefficient does not become a small positive number: the penalty’s subgradient at zero permits the whole correlation to be balanced by the penalty.

Repeated Updates and Sparse Coefficients

Coordinate descent applies the one-coordinate minimization repeatedly, using the latest values of the other coefficients to form each partial residual. A change in one coefficient changes the partial residuals, and therefore the correlations used in subsequent updates. An update that sets a coefficient to zero removes that predictor’s current contribution to the linear predictor, but later updates of other coordinates can change the residual pattern. The algorithm must therefore revisit coordinates rather than treating a single sweep as final.

The next example uses a fixed weighted least-squares surrogate with no intercept and two predictors. It traces successive updates in the order \(\beta_1\), then \(\beta_2\). This isolates the coordinate mechanics: in a GLM, the same calculations occur within an IRLS weighted least-squares subproblem.

Worked Example: Two Repeated Coordinate Updates

Let \(n=2\), \(w_1=w_2=1\), \(z=(2,1)\), \(x_1=(1,0)\), and \(x_2=(1,1)\). There is no intercept or offset in this example. Take \(\lambda=0.1\) and start at \((\beta_1,\beta_2)=(0,0)\). Use the update formula first for coordinate 1 and then for coordinate 2.

Initially, the partial residual for coordinate 1 is \(r_1=z\). Its norm and correlation are \(a_1=(1/2)(1^2+0^2)=0.5\) and \(c_1=(1/2)(1\cdot2+0\cdot1)=1\). Thus \(\beta_1=S(1,0.1)/0.5=0.9/0.5=1.8\).

Now hold \(\beta_1=1.8\) fixed. The partial residual for coordinate 2 is \(r_2=z-x_1\beta_1=(0.2,1)\). Since \(a_2=(1/2)(1^2+1^2)=1\) and \(c_2=(1/2)(1\cdot0.2+1\cdot1)=0.6\), the update is \(\beta_2=S(0.6,0.1)/1=0.5\).

For a second pass, hold \(\beta_2=0.5\) fixed and update coordinate 1. Its partial residual is \(z-x_2\beta_2=(1.5,0.5)\), so \(c_1=(1/2)(1\cdot1.5+0\cdot0.5)=0.75\). With \(a_1=0.5\), the new value is \(\beta_1=S(0.75,0.1)/0.5=0.65/0.5=1.3\).

Next, with \(\beta_1=1.3\), coordinate 2 has partial residual \(z-x_1\beta_1=(0.7,1)\). Thus \(c_2=(1/2)(0.7+1)=0.85\), and \(\beta_2=S(0.85,0.1)/1=0.75\). The successive pairs are \((0,0)\), \((1.8,0.5)\), and \((1.3,0.75)\). Each update uses the latest partial residual, which is why the second pass differs from the first.

Result: Each exact coordinate update leaves the other coordinates fixed and minimizes the weighted least-squares surrogate over the selected coordinate. Consequently, the surrogate cannot increase at that update. If the selected coordinate has \(|c_j|\leq\lambda\), its minimizing value is exactly zero. These facts explain how repeated coordinate minimization can produce sparse estimates; they do not, by themselves, establish a particular convergence rate or guarantee for an entire GLM algorithm.

The nonincrease claim follows directly from minimization: if the old coordinate value is \(b_{\mathrm{old}}\) and the new value is \(b_{\mathrm{new}}\), then \(q_j(b_{\mathrm{new}})\leq q_j(b_{\mathrm{old}})\). Every term in the surrogate that does not depend on coordinate \(j\) is unchanged, so the full surrogate has the same nonincrease. This argument applies to an exact minimization of the current quadratic surrogate. When an outer IRLS step changes the working response or weights, the surrogate itself changes, so this comparison alone is not a proof that the original penalized likelihood decreases at every algorithmic step.

A GLM Coordinate Update in One IRLS Step

For a concrete link to GLMs, consider a Bernoulli response with a logistic link. At current linear predictors \(\eta_i\), the fitted means are \(\mu_i=1/(1+\exp(-\eta_i))\). The usual IRLS working quantities for a Bernoulli observation are \(w_i=\mu_i(1-\mu_i)\) and \(z_i=\eta_i+(y_i-\mu_i)/w_i\). The lasso coordinate formula applies to these working quantities, with any other coefficients and the intercept held at their current values.

Worked Example: One Penalized Logistic-Regression Surrogate Update

Take two observations with \(y=(0,1)\), one predictor \(x=(-1,1)\), no offset, and current intercept and slope both equal to zero. At these current values, \(\eta_i=0\) and \(\mu_i=0.5\) for each observation. Hence \(w_1=w_2=0.5(0.5)=0.25\). The working responses are

z_1=0+\frac{0-0.5}{0.25}=-2, \qquad z_2=0+\frac{1-0.5}{0.25}=2.

For the slope update, hold the intercept at zero. The partial residuals are therefore \((-2,2)\). With \(n=2\),

a_j=\frac{1}{2}\{0.25(-1)^2+0.25(1)^2\}=0.25, \qquad c_j=\frac{1}{2}\{0.25(-1)(-2)+0.25(1)(2)\}=0.5.

If \(\lambda=0.25\), the slope update is

\beta_j^{\mathrm{new}}=\frac{S(0.5,0.25)}{0.25} =\frac{0.25}{0.25}=1.

This is the slope minimizer for this weighted least-squares surrogate with the intercept held fixed. It is not, by itself, the exact minimizer of the original penalized logistic likelihood. After the inner coordinate updates, a GLM fitting procedure may refresh the working quantities and continue its outer iterations.

Common Mistakes and Practical Interpretation

Using the full residual instead of the partial residual. The correlation for coordinate \(j\) must exclude the current contribution \(x_{ij}\beta_j\). Otherwise the coefficient being updated is incorrectly treated as fixed inside its own residual.
Thresholding after dividing by the norm. With the stated scaling, the correct update is \(S(c_j,\lambda)/a_j\), not \(S(c_j/a_j,\lambda)/a_j\). The penalty threshold acts on the correlation before division by \(a_j\).
Assuming every small correlation gives a small nonzero estimate. If \(|c_j|\leq\lambda\), the minimizer is exactly zero. That thresholding is the source of lasso sparsity in the coordinate subproblem.
Calling a surrogate update an exact GLM likelihood update. The formula exactly minimizes one coordinate of the current penalized quadratic surrogate. The working weights and response depend on the current GLM fit and are generally recomputed in subsequent outer iterations.
Forgetting the loss scaling. Changing \(1/(2n)\) in the quadratic loss changes the scale of the correlation relative to \(\lambda\). The threshold and penalty must use a consistent objective convention.

In practical computation, the weighted correlation answers whether the current predictor has enough residual association to overcome the penalty. If not, its coefficient is set to zero for that update. If it does, the excess beyond the threshold determines the signed numerator, and the weighted predictor norm sets the scale of the resulting coefficient. Repeatedly recalculating these quantities is what lets the method adapt to the other coefficients while maintaining sparse solutions in each weighted least-squares subproblem.

Check Your Understanding

Use the stated objective scaling and coordinate-update formula to answer the following.

For \(a_j=2\), \(c_j=-0.8\), and \(\lambda=0.3\), what is the coordinate update?
State the necessary and sufficient condition for a coordinate minimizer to be zero.
Why does the partial residual used for coordinate \(j\) exclude \(x_{ij}\beta_j\)?
What does an exact coordinate update guarantee about the current weighted least-squares surrogate?
Why is one soft-thresholding update in an IRLS step not necessarily the exact minimizer of the original penalized GLM objective?