the same trick applies to all other solvers and I think we should implement it for those as well. I'd be happy to have this in just for the default solver, though.

hm right now subtracting the mean even in the sparse case... need to benchmark/think about the sparse case a bit more

This is actually a regression because this used to work with liblinear, which was the default solver in 0.21.

In the binary case, liblinear is slightly worse than the preconditioned solution by default, but much better than current master.

If we want to scale the data, I think it should be by default (either with a parameter or not), otherwise one will still get the warning anyway so it does not solve anything.

If we start scaling the data for the user then I am wondering about the following. One will not be able to robust scale the data since that we will reapply a standard scaling on the top. Does the optimization problem to be solved will be affected and does it matter?

@glemaitre just to be clear, this PR doesn't scale the data. It solves the original optimization problem, just better and faster (if I scaled the data, the loss on the unscaled data would be high, it's not). That was @GaelVaroquaux's suggestion and I agree that it's more in line with our general principles (and it addresses your concern).

Can we pass the X_pre, X_mean, and X_scale into _logistic_loss_and_grad and friends and use this to compute the loss and gradients? (This is to try to avoid scaling sparse data)

Right now I'm doing that for X_scale (it's the only way to make this work), but not X_mean. But should be possible.

Also, it does feel weird to have it for one solver and not the others. But I guess the docs can explain that.

@adrinjalali as I said above, I intend to remove the parameter everywhere again and it's just for easy comparison.

removing from the milestone (this one got complicated)

@tomMoral the synthetic data I mentioned to trigger conditioning issues for benchopt is here: #15583 (comment)

Hm the svd on subsampled data idea sounds good... damn, I haven't looked at this in a while...

Regarding #15583 (comment), it was pointed out to me that we might be favoring the version without preconditioning because the first and last feature are actually much more important given how we set up the problem.
If we adjust the coefficents to reflect the change of units we did for the features, the results are different:

from time import time
import numpy as np
from sklearn.datasets import make_low_rank_matrix
from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import scale

n_samples, n_features = 1000, 10000
max_iter = 1000
print(f"n_samples={n_samples}, n_features={n_features}")

rng = np.random.RandomState(0)
w_true = rng.randn(n_features)

X = make_low_rank_matrix(n_samples, n_features, random_state=rng)
# X = rng.randn(n_samples, n_features)
X[:, 0] *= 1e3
X[:, -1] *= 1e3
w_true[0] /= 1e3
w_true[-1] /= 1e3

z = X @ w_true + 1
z += 1e-1 * rng.randn(n_samples)

# Balanced binary classification problem
y = (z > np.median(z)).astype(np.int32)

for C in [1e3, 1, 1e-3]:
    print(f"\nC={C}")
    print("Without data preconditioning:")
    clf = LogisticRegression(precondition=False, C=C,
                             max_iter=max_iter)
    tic = time()
    clf.fit(X, y)
    duration = time() - tic
    print(f"n_iter={clf.n_iter_[0]},"
          f" obj={clf.objective_value_:.6f},"
          f" duration={duration:.3f}s")

    print("Without data preconditioning with scaling:")
    clf_scaled = LogisticRegression(precondition=False, C=C,
                                    max_iter=max_iter)

    print("With data preconditioning:")
    clf_pre = LogisticRegression(precondition=True, C=C,
                                 max_iter=max_iter)
    tic = time()
    clf_pre.fit(X, y)
    duration_pre = time() - tic
    print(f"n_iter={clf_pre.n_iter_[0]},"
          f" obj={clf_pre.objective_value_:.6f},"
          f" duration={duration_pre:.3f}s")
    print(f"speedup (pre): {duration/duration_pre:.1f}x")

n_samples=1000, n_features=10000

C=1000.0
Without data preconditioning:
n_iter=809, obj=367.675372, duration=5.530s
Without data preconditioning with scaling:
With data preconditioning:
n_iter=144, obj=367.674329, duration=1.083s
speedup (pre): 5.1x

C=1
Without data preconditioning:
n_iter=24, obj=675.094608, duration=0.196s
Without data preconditioning with scaling:
With data preconditioning:
n_iter=31, obj=675.094611, duration=0.295s
speedup (pre): 0.7x

C=0.001
Without data preconditioning:
n_iter=8, obj=692.023023, duration=0.083s
Without data preconditioning with scaling:
With data preconditioning:
n_iter=23, obj=692.024219, duration=0.249s
speedup (pre): 0.3x

Note: as discussed in #18795, it seems that when features are proprocessed with PolynomialCountSketch(degree=2, n_components=1000) or more we get another real life case where our linear model solvers struggle with data-conditioning issues (although I did not investigate in details).

I also wonder if multiclass problems with many classes including some rare classes could not cause another family of data-related conditioning issues (for LogisticRegression with the multinomial loss).

This seems a bit stalled.

The findings of our study confirm that preconditioning least-squares problems is hard and
that at present there is no single approach that works well for all problems.

Gould, N.I., & Scott, J.A. (2017). The State-of-the-Art of Preconditioners for Sparse Linear Least-Squares Problems. ACM Transactions on Mathematical Software (TOMS), 43, 1 - 35. https://doi.org/10.1145/3014057 PDF

Edit: More on that matter in

• Preconditioned Conjugate Gradient Methods in Truncated Newton Frameworks for Large-scale Linear Classification (2018)
• A Study on Truncated Newton Methods for Linear Classification (2022)

It is also a bit ironic to remove the "normalize" parameter from so many linear models and to try to reinvent it for logistic regression 😄

Note that the glum authors mention that intern standardization is beneficial for the optimizer, see https://glum.readthedocs.io/en/latest/background.html#Standardization and https://glum.readthedocs.io/en/latest/motivation.html#software-optimizations.

See #15556 (comment).