✔️ Linting Passed

All linting checks passed. Your pull request is in excellent shape! ☀️

_{Generated for commit: 537b41a. Link to the linter CI: here}

A little benchmark on the leukemia dataset

1. Full Lasso Path with 100 steps

n_samples = 72, n_features = 7129

Starting path computation...
Lasso without screening - time: 0.18 s
Lasso with screening    - time: 0.10 s
Lasso without screening - time: 1.66 s
Lasso with screening    - time: 0.53 s
Lasso without screening - time: 9.75 s
Lasso with screening    - time: 1.62 s

Code

# https://github.com/mathurinm/celer/blob/main/examples/plot_leukemia_path.py
import time
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

from sklearn.linear_model import lasso_path
from sklearn.datasets import fetch_openml

print(__doc__)

print("Loading data...")
dataset = fetch_openml("leukemia")
X = np.asfortranarray(dataset.data.astype(float))
y = 2 * ((dataset.target != "AML") - 0.5)
n_samples = len(y)

y -= np.mean(y)
y /= np.std(y)

print("Starting path computation...")
alpha_max = np.max(np.abs(X.T.dot(y))) / n_samples

n_alphas = 100
alphas = alpha_max * np.geomspace(1, 0.01, n_alphas)

tols = [1e-2, 1e-3, 1e-4]
results = np.zeros([2, len(tols)])
for tol_ix, tol in enumerate(tols):
    t0 = time.time()
    _, coefs, gaps = lasso_path(
        X, y, alphas=alphas, tol=tol, max_iter=10_000, do_screening=False)
    results[0, tol_ix] = time.time() - t0
    print('Lasso without screening - time: %.2f s' % results[0, tol_ix])

    t0 = time.time()
    _, coefs, dual_gaps = lasso_path(
        X, y, alphas=alphas, tol=tol, max_iter=10_000, do_screening=True)
    results[1, tol_ix] = time.time() - t0
    print('Lasso with screening    - time: %.2f s' % results[1, tol_ix])

df = pd.DataFrame(results.T, columns=["no screening", "with screening"])
df.index = [str(tol) for tol in tols]
df.plot.bar(rot=0)
plt.xlabel("stopping tolerance")
plt.ylabel("path computation time (s)")
plt.tight_layout()
plt.show()

2. Single alpha

Code

# For single values of alpha
import time
import pandas as pd
import numpy as np
import seaborn as sns

from sklearn.linear_model import lasso_path
from sklearn.datasets import fetch_openml


print("Loading data...")
dataset = fetch_openml("leukemia")
X = np.asfortranarray(dataset.data.astype(float))
y = 2 * ((dataset.target != "AML") - 0.5)
n_samples = len(y)

y -= np.mean(y)
y /= np.std(y)

print("Starting path computation...")
alpha_max = np.max(np.abs(X.T.dot(y))) / n_samples

n_alphas = 1

tols = [1e-2, 1e-3, 1e-4]
results = []  # scale_alpha_max, tol, do_screening, time
penalty_scales = [1e-1, 1e-2, 1e-5]

for alpha_idx, scale_alpha_max in enumerate(penalty_scales):
    alpha = [alpha_max * scale_alpha_max]
    print(f"alpha = alpha_max * {scale_alpha_max}")
    for tol_ix, tol in enumerate(tols):
        t0 = time.time()
        _, coefs, gaps = lasso_path(
            X, y, alphas=alpha, tol=tol, max_iter=10_000, do_screening=False)
        t = time.time() - t0
        results += [[scale_alpha_max, tol, False, t]] 
        print('  Lasso without screening - time: %.2f s' % t)

        t0 = time.time()
        _, coefs, dual_gaps = lasso_path(
            X, y, alphas=alpha, tol=tol, max_iter=10_000, do_screening=True)
        t = time.time() - t0
        results += [[scale_alpha_max, tol, True, t]] 
        print('  Lasso with screening    - time: %.2f s' % t)

df = pd.DataFrame(
    results,
    columns=["scale_alpha_max", "tol", "do_screening", "time [s]"],
)
df

sns.catplot(
    df, x="tol", y="time [s]", hue="do_screening", col="scale_alpha_max", kind="bar", errorbar=None, sharey=False)

cc @mathurinm @agramfort @josephsalmon

Based on the build failure of some CI runners, it seems that this code requires bumping up the minimum Cython version. I think it's fine because it's a not runtime dependency.

Any idea why this PR would lead to a change of predicted values in this case?

2117 Examples
2118 --------
2119 >>> from sklearn.linear_model import LassoCV
2120 >>> from sklearn.datasets import make_regression
2121 >>> X, y = make_regression(noise=4, random_state=0)
2122 >>> reg = LassoCV(cv=5, random_state=0).fit(X, y)
2123 >>> reg.score(X, y)
2124 0.9993
2125 >>> reg.predict(X[:1,])
Expected:
    array([-78.4951])
Got:
    array([-79.4755331])

Maybe we have nearly tied choices for the optimal value of alpha on this data?

Any idea why this PR would lead to a change of predicted values in this case?

The solution is within the uncertainty of tol. One of the main changes is that this PR adds an additional convergence check before the main loop: if the provided (warm started) coefficients already satisfy the stopping criterion, then stop. In main, it always calculates one complete cycle over coordinates/features.

Also: could you update the docstrings and the user guide to reference the gap safe screening rule paper?

First of all, @ogrisel, thanks for your review, much appreciated.

From a cursory glance at the diff and the paper, this PR implements some basic active set based on the gap safe rule written in the Celer paper, however it does not seem to implement the precise active set construction described in Algorithm 4. On top of this, the Celer project implements other things such as the "dual extrapolation" from https://arxiv.org/abs/1907.05830 which probably amplifies the expected speed-ups. So I suspect that the performance of this PR is still far from the potential gains one could obtain from the full-fledged Celer implementation. Have you tried to run some benchmark against it to assess the remaining difference?

And I have nowhere claimed to have added working sets or dual extrapolation, see PR description. Again, this is (under 200 lines of code) was (now over 500 loc) a pretty neat and small PR, introducing screening rules for a single Cython function that already gives a massive (but not maximum) speedup.

The neat thing is that it is minimal:

• Add check of stopping criterion before entering the main loop.
• Whenever the stopping criterion is calulcated, we already have the dual gap. Therefore we then screen variables. (This step could later be fine tuned).

I don't mean to block this PR for that because it seems to already provide a net improvement over main, and I am fine with proceeding iteratively, but we should not imply in the doc or changelog that scikit-learn implements all the things described in the Celer paper (and in even less so if you consider the dual extrapolation paper on top).

It nowhere does!

Also: could you update the docstrings and the user guide to reference the gap safe screening rule paper?

This took me a whole day, that's longer than the implementation of the screening rule itself!

And I have nowhere claimed to have added working sets or dual extrapolation, see PR description.

Sorry, I did not mean to imply that you made that claim in the PR, but I was just pointing that we should not leave the impression that we implement the full Celer algorithm when referencing the Celer paper in the docstring of the estimators (or the user guide).

Task list:

• remove attribute self._do_screening
• merge with main after Cython update MNT update Cython 3.0.10 to 3.1.2 #31905
  Solved/Circumvented with 46d3a09 and e71e89b
• if ENH/FIX stopping criterion for coordinate descent gap <= tol #31906 gets merged first (preferred way), update to tol <=

@mathurinm Thanks for your review. Very valuable!

If CD pipelines are green this is ready for merge.

Note that I removed all the private do_screening things. Only in enet_path one can pass do_screening as private (and undocumented) variable (in params).

Were are back at a neat PR: Without the addition to the user guide the changed lines of code are: +278 -78