[mlakolar]

Summary

Two changes to the scaled lasso of Sun and Zhang (2012):

1. scaledLasso! solves the square-root lasso in one run. The scaled lasso minimizes ‖r‖²/(2nσ) + σ/2 + λ ∑ ωⱼ|βⱼ|; the minimum over σ is at σ̂ = ‖r‖/√n, which leaves the square-root lasso ‖r‖/√n + λ ∑ ωⱼ|βⱼ|. scaledLasso! now solves that by one run of coordinate descent on CDSqrtLassoLoss (penalty √n λ). Before, it alternated lasso fits with updates of σ̂ until σ̂ changed by less than optTol (1e-2 by default), which stopped short of the solution and refitted a lasso at every iteration. Of IterLassoOptions it now uses only optionsCD (the other fields still drive feasibleLasso!).
2. CDSqrtQuadraticLoss(A, b, c; check=true), the square-root lasso loss √(x'Ax + 2b'x + c) from sufficient statistics (A = X'X/n, b = -X'y/n, c = y'y/n gives ‖y - Xx‖/√n). It caches Ax and the quadratic, so a coordinate update costs O(p) instead of O(n) and a coordinate that does not move O(1); the step is the exact minimizer, with the closed form of CDSqrtLassoLoss. scaledLasso!(x, f::CDSqrtQuadraticLoss, λ, ω, options) solves the scaled lasso from it. With A = S, b = -S[:, j], c = S[j, j] and ω[j] = Inf (an infinite weight keeps a coordinate at zero) it is the regression of variable j on the others in the precision matrix estimator of Sun and Zhang (2013); check = false skips the O(p²) checks of A when many losses share it.

Momentum and Anderson acceleration move the affine cache Ax; a hook, _cacheMoved!, computes the quadratic again after a move.

sol = scaledLasso(X, y, λ, ω)                       # one run of coordinate descent
f = CDSqrtQuadraticLoss(X'X / n, -X'y / n, dot(y, y) / n)
sol = scaledLasso!(SparseIterate(p), f, λ, ω)       # from sufficient statistics
# variable j on the others, from a covariance matrix S checked once
ω = sqrt.(diag(S)); ω[j] = Inf
sol = scaledLasso!(SparseIterate(p), CDSqrtQuadraticLoss(S, -S[:, j], S[j, j]; check=false), λ0, ω)

Benchmarks

benchmark/scaledlasso_bench.jl on USC CARC Discovery (one CPU, AMD EPYC 7513, Julia 1.13.2; run 20261023T1216Z-cd-scaledlasso; tables in benchmark/README.md):

• Regression. The new scaledLasso was 1.2–5.4× faster than the default of 0.6.0 in 14 of 15 cases (3% slower in the last, n = 5000, p = 500, ρ = 0.9), whose coefficients were off by 2e-5 to 3e-2 (the new ones by at most 6e-6); against 0.6.0 at optTol = 1e-6, 1.2–19× faster.
• Sufficient statistics were 3.6–20× faster than the data when n > p (forming A costs O(np²) once), and slower when p ≥ n (up to 25× at p = 5000, n = 200).
• The precision matrix estimator of Sun and Zhang, all p columns from the shared S: 2–90 ms for p = 200–1000, 3.4–62× faster than the new scaledLasso on copies of the data, 7–73× faster than 0.6.0's default.

Against CovSel's exact homotopy (slassoinv in CovSel #16, which follows the lasso path of each column to the scaled-lasso fixed point exactly), on its 30 speed problems (AR(1), band, Erdős–Rényi; p = 100–2000; n = p/2 and 2p; same λ₀): column by column with CDSqrtQuadraticLoss(...; check=false) took 0.9–2.1× the time of slassoinv with the default options, 0.6–1.6× with order = :cyclic (faster than the homotopy for Erdős–Rényi graphs and larger p), and 1.0–2.8× at optTol = 1e-10; the estimates agree to 1e-6 at the default tolerance and to 1e-9 at optTol = 1e-10. Before check = false, the O(p²) checks of S in every loss made it 26–76× slower at p = 2000. So CovSel could call CoordinateDescent for this estimator at about the same speed, trading the exact solution for one within the tolerance.

Behavior changes

• scaledLasso! returns the solution of the jointly convex problem (to the tolerance of coordinate descent); results differ from 0.6.0 by up to the old optTol. The fields maxIter, optTol, initProcedure, sinit, σinit of IterLassoOptions no longer affect it, and it no longer warns about the outer iteration.
• The internal _findInitSigma! is removed.
• Version 0.6.1.

Relation to #23

Draft #23 (CoordinateDescent 0.7 on ProximalOPT) deletes the engine files this PR touches (cd_differentiable_function.jl, anderson.jl, momentum.jl) and rewrites lasso.jl. Porting it means a ProximalOPT model for the square root of a quadratic (with the cached quadratic recomputed after momentum moves, as _cacheMoved! does here), and scaledLasso! on ProximalOPT's square-root least squares model.

Tests

CDSqrtQuadraticLoss against CDSqrtLassoLoss (same solution, with momentum on and off, the cyclic order with Anderson acceleration, and the path from zero), its value and gradient, infinite weights, a coordinate the loss does not depend on, check = false, input checks; scaledLasso at the optimality conditions of the lasso at λσ̂; from sufficient statistics against the data; a precision-matrix column against the regression on the other columns; argument conversion. The old outer-iteration test now checks that maxIter/optTol do not matter. 1181 tests pass on Julia 1.13 and 1.10; examples/04_sqrt_scaled_feasible.jl updated and runs.

🤖 Generated with Claude Code

[codecov]

Codecov Report

❌ Patch coverage is 97.26027% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.81%. Comparing base (1277dab) to head (b31ae4a).

Files with missing lines	Patch %	Lines
src/cd_differentiable_function.jl	96.42%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #24      +/-   ##
==========================================
- Coverage   97.88%   97.81%   -0.08%     
==========================================
  Files           9        9              
  Lines        1322     1370      +48     
==========================================
+ Hits         1294     1340      +46     
- Misses         28       30       +2     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.