Adding a converged lasso by coordinate descent next to the as-submitted one

• Status: Accepted
• Date recorded: 2026-10-06
• Decision in one line: The 2026 analyses add a lasso solved to convergence by cyclic coordinate descent, verified against scikit-learn, while the as-submitted single-step lasso stays wherever the report's numbers are reproduced.

Context

The notebook's LassoRegression was meant to be iterative soft-thresholding (ISTA). Two quirks change what it computes. Every iteration restarts from the zero vector Ao, so the ten repetitions all return the same single thresholded gradient step. The survivors are also scaled by (1 / 1 + thr), which Python reads as 1 + thr, so they are scaled up instead of down. The port reproduces both on purpose, because the report's numbers depend on them.

The 2026 bias-variance analysis needs a lasso that actually solves its optimisation problem. With the single step, the estimates are heavily shrunk even at ρ = 0, so a bias-variance curve would mostly describe the bug.

Decision

Add lassoCoordinateDescent alongside the original. It is labelled "converged lasso (2026)" wherever it appears and is used only in the 2026 analyses on /bias-variance and /recovery. It does not replace the as-submitted lasso on /retrieve, /tune or /pcr.

For each pixel v it minimises

with the update $a_j \leftarrow S\big(c_j - \sum_{k\ne j} G_{jk} a_k,\ N\rho\big) / G_{jj}$, where $G = D^\top D$ and $c = D^\top x_v$. This is the problem the 2021 code was approximating, because its threshold N ρ × step is a proximal step on the same objective. It is also the problem scikit-learn's Lasso(alpha=ρ, fit_intercept=False) solves, which makes it checkable.

Options considered

1. Fix the original ISTA loop, carrying the iterate forward and correcting the scale. With the notebook's step size of $1 / (1.1\,\lVert D^\top D\rVert_F)$, convergence is slow and needs many hundreds of iterations per pixel.
2. Accelerated proximal gradient (FISTA). Faster than ISTA, but it still needs a step size and a stopping rule tuned for this problem.
3. Cyclic coordinate descent with a shared Gram matrix. Each pixel has only six coefficients and every pixel shares DᵀD, so each sweep is cheap and the updates are exact.
4. Call a library. There is no maintained lasso for the browser, and compiling scikit-learn to WebAssembly would be heavy for one function.
5. A closed form. The lasso has one only for orthogonal designs, and these time courses are correlated, with |r| up to 0.77 between TC4, TC5 and TC6.

Why

Coordinate descent is the standard lasso solver (it is what scikit-learn uses), it has no step size to tune, and with six coefficients per pixel it converges in a handful of sweeps. Solving exactly the objective the 2021 code intended keeps the comparison honest. The difference between the two lassos is then the solver alone, with the same ρ meaning the same penalty.

What happened

• On the 2021 data at ρ = 0.05, 0.40 and 0.625 the coefficients match scikit-learn (tolerance 1e-14) to better than 1e-9, including the exact number of non-zero coefficients. At ρ = 0 it equals least squares, and at ρ = 1 every coefficient is zero. Coordinate descent converged in every fit of the 50-realisation bias-variance run, and /bias-variance reports this for every run.
• At the 2021 penalty ρ = 0.625 the two solvers give clearly different maps (some coefficients differ by more than 0.05).
• With 50 realisations, the converged lasso's error on the spatial maps, relative to the noise-free target, is lowest at ρ = 0.40 (0.437, 95% interval 0.436 to 0.439). The as-submitted lasso is best at ρ = 0.25 (0.726), and its bias² stays above 0.59 at every ρ.
• Across 100 realisations at ρ = 0.625, the converged lasso recovers the spatial maps better than the as-submitted one. The aligned $\Sigma c_S$ rises from 4.82 to 5.15, a paired difference of 0.32 (95% interval 0.31 to 0.33) that holds in all 100 realisations.

Not everything improved. On the time courses at ρ = 0.625 the converged lasso scores slightly lower than the as-submitted one ($\Sigma c_T$ 5.419 against 5.435). Both solvers also break down in the same way when temporal noise grows. With the penalty fixed at 0.625 and temporal noise variance 2, the mean recovery of the maps falls to about 0.1, because standardising X shrinks the signal below the fixed threshold.

What I'd change

• Tune ρ for the converged solver directly, by cross-validation over pixels or by the bias-variance target, instead of reusing the value chosen for the single-step solver.
• Scale the penalty with an estimate of the noise level, so the same rule works when the noise changes.
• Keep a converged reference implementation from the start, and test any hand-written solver against it before tuning anything with it.