scikit-learn · lesteve · Jun 24, 2022 · Jun 15, 2022 · Jun 15, 2022 · Jun 16, 2022
diff --git a/doc/whats_new/v1.2.rst b/doc/whats_new/v1.2.rst
@@ -142,6 +142,11 @@ Changelog
   `solver="newton-cg"`, `fit_intercept=True`, and a single feature. :pr:`23608`
   by `Tom Dupre la Tour`_.
 
+- |API| The default value for the `solver` parameter in
+  :class:`linear_model.QuantileRegressor` will change from `"interior-point"`
+  to `"highs"` in version 1.4.
+  :pr:`23637` by :user:`Guillaume Lemaitre <glemaitre>`.
+
 :mod:`sklearn.metrics`
 ......................
 

diff --git a/examples/linear_model/plot_quantile_regression.py b/examples/linear_model/plot_quantile_regression.py
@@ -111,13 +111,20 @@
 #
 # We will use the quantiles at 5% and 95% to find the outliers in the training
 # sample beyond the central 90% interval.
+from sklearn.utils.fixes import sp_version, parse_version
+
+# This is line is to avoid incompatibility if older SciPy version.
+# You should use `solver="highs"` with recent version of SciPy.
+solver = "highs" if sp_version >= parse_version("1.6.0") else "interior-point"
+
+# %%
 from sklearn.linear_model import QuantileRegressor
 
 quantiles = [0.05, 0.5, 0.95]
 predictions = {}
 out_bounds_predictions = np.zeros_like(y_true_mean, dtype=np.bool_)
 for quantile in quantiles:
-    qr = QuantileRegressor(quantile=quantile, alpha=0)
+    qr = QuantileRegressor(quantile=quantile, alpha=0, solver=solver)
     y_pred = qr.fit(X, y_normal).predict(X)
     predictions[quantile] = y_pred
 
@@ -179,7 +186,7 @@
 predictions = {}
 out_bounds_predictions = np.zeros_like(y_true_mean, dtype=np.bool_)
 for quantile in quantiles:
-    qr = QuantileRegressor(quantile=quantile, alpha=0)
+    qr = QuantileRegressor(quantile=quantile, alpha=0, solver=solver)
     y_pred = qr.fit(X, y_pareto).predict(X)
     predictions[quantile] = y_pred
 
@@ -250,7 +257,7 @@
 from sklearn.metrics import mean_squared_error
 
 linear_regression = LinearRegression()
-quantile_regression = QuantileRegressor(quantile=0.5, alpha=0)
+quantile_regression = QuantileRegressor(quantile=0.5, alpha=0, solver=solver)
 
 y_pred_lr = linear_regression.fit(X, y_pareto).predict(X)
 y_pred_qr = quantile_regression.fit(X, y_pareto).predict(X)

diff --git a/sklearn/linear_model/_quantile.py b/sklearn/linear_model/_quantile.py
@@ -44,10 +44,16 @@ class QuantileRegressor(LinearModel, RegressorMixin, BaseEstimator):
     solver : {'highs-ds', 'highs-ipm', 'highs', 'interior-point', \
             'revised simplex'}, default='interior-point'
         Method used by :func:`scipy.optimize.linprog` to solve the linear
-        programming formulation. Note that the highs methods are recommended
-        for usage with `scipy>=1.6.0` because they are the fastest ones.
-        Solvers "highs-ds", "highs-ipm" and "highs" support
-        sparse input data and, in fact, always convert to sparse csc.
+        programming formulation.
+
+        From `scipy>=1.6.0`, it is recommended to use the highs methods because
+        they are the fastest ones. Solvers "highs-ds", "highs-ipm" and "highs"
+        support sparse input data and, in fact, always convert to sparse csc.
+
+        From `scipy>=1.11.0`, "interior-point" is not available anymore.
+
+        .. versionchanged:: 1.4
+           The default of `solver` will change to `"highs"` in version 1.4.
 
     solver_options : dict, default=None
         Additional parameters passed to :func:`scipy.optimize.linprog` as
@@ -91,7 +97,10 @@ class QuantileRegressor(LinearModel, RegressorMixin, BaseEstimator):
     >>> rng = np.random.RandomState(0)
     >>> y = rng.randn(n_samples)
     >>> X = rng.randn(n_samples, n_features)
-    >>> reg = QuantileRegressor(quantile=0.8).fit(X, y)
+    >>> # the two following lines are optional in practice
+    >>> from sklearn.utils.fixes import sp_version, parse_version
+    >>> solver = "highs" if sp_version >= parse_version("1.6.0") else "interior-point"
+    >>> reg = QuantileRegressor(quantile=0.8, solver=solver).fit(X, y)
     >>> np.mean(y <= reg.predict(X))
     0.8
     """
@@ -102,7 +111,7 @@ def __init__(
         quantile=0.5,
         alpha=1.0,
         fit_intercept=True,
-        solver="interior-point",
+        solver="warn",
         solver_options=None,
     ):
         self.quantile = quantile
@@ -166,14 +175,14 @@ def fit(self, X, y, sample_weight=None):
                 f"The argument fit_intercept must be bool, got {self.fit_intercept}"
             )
 
-        if self.solver not in (
-            "highs-ds",
-            "highs-ipm",
-            "highs",
-            "interior-point",
-            "revised simplex",
-        ):
-            raise ValueError(f"Invalid value for argument solver, got {self.solver}")
+        if self.solver == "warn":
+            warnings.warn(
+                "The default solver will change from 'interior-point' to 'highs' in "
+                "version 1.4. Set `solver='highs'` or to the desired solver to silence "
+                "this warning.",
+                FutureWarning,
+            )
+            solver = "interior-point"
         elif self.solver in (
             "highs-ds",
             "highs-ipm",
@@ -183,8 +192,23 @@ def fit(self, X, y, sample_weight=None):
                 f"Solver {self.solver} is only available "
                 f"with scipy>=1.6.0, got {sp_version}"
             )
+        else:
+            solver = self.solver
+
+        if solver not in (
+            "highs-ds",
+            "highs-ipm",
+            "highs",
+            "interior-point",
+            "revised simplex",
+        ):
+            raise ValueError(f"Invalid value for argument solver, got {solver}")
+        elif solver == "interior-point" and sp_version >= parse_version("1.11.0"):
+            raise ValueError(
+                f"Solver {solver} is not anymore available in SciPy >= 1.11.0."
+            )
 
-        if sparse.issparse(X) and self.solver not in ["highs", "highs-ds", "highs-ipm"]:
+        if sparse.issparse(X) and solver not in ["highs", "highs-ds", "highs-ipm"]:
             raise ValueError(
                 f"Solver {self.solver} does not support sparse X. "
                 "Use solver 'highs' for example."
@@ -200,7 +224,7 @@ def fit(self, X, y, sample_weight=None):
             )
 
         # make default solver more stable
-        if self.solver_options is None and self.solver == "interior-point":
+        if self.solver_options is None and solver == "interior-point":
             solver_options = {"lstsq": True}
         else:
             solver_options = self.solver_options
@@ -243,7 +267,7 @@ def fit(self, X, y, sample_weight=None):
             c[0] = 0
             c[n_params] = 0
 
-        if self.solver in ["highs", "highs-ds", "highs-ipm"]:
+        if solver in ["highs", "highs-ds", "highs-ipm"]:
             # Note that highs methods always use a sparse CSC memory layout internally,
             # even for optimization problems parametrized using dense numpy arrays.
             # Therefore, we work with CSC matrices as early as possible to limit
@@ -268,7 +292,7 @@ def fit(self, X, y, sample_weight=None):
             c=c,
             A_eq=A_eq,
             b_eq=b_eq,
-            method=self.solver,
+            method=solver,
             options=solver_options,
         )
         solution = result.x

diff --git a/sklearn/linear_model/tests/test_quantile.py b/sklearn/linear_model/tests/test_quantile.py
@@ -23,6 +23,11 @@ def X_y_data():
     return X, y
 
 
+@pytest.fixture
+def default_solver():
+    return "highs" if sp_version >= parse_version("1.6.0") else "interior-point"
+
+
 @pytest.mark.parametrize(
     "params, err_msg",
     [
@@ -40,6 +45,10 @@ def X_y_data():
         ),
     ],
 )
+@pytest.mark.filterwarnings(
+    # FIXME (1.4): remove once we changed the default solver to "highs"
+    "ignore:The default solver will change from 'interior-point'"
+)
 def test_init_parameters_validation(X_y_data, params, err_msg):
     """Test that invalid init parameters raise errors."""
     X, y = X_y_data
@@ -85,11 +94,13 @@ def test_too_new_solver_methods_raise_error(X_y_data, solver):
         [0.5, 100, 2, 0],
     ],
 )
-def test_quantile_toy_example(quantile, alpha, intercept, coef):
+def test_quantile_toy_example(quantile, alpha, intercept, coef, default_solver):
     # test how different parameters affect a small intuitive example
     X = [[0], [1], [1]]
     y = [1, 2, 11]
-    model = QuantileRegressor(quantile=quantile, alpha=alpha).fit(X, y)
+    model = QuantileRegressor(
+        quantile=quantile, alpha=alpha, solver=default_solver
+    ).fit(X, y)
     assert_allclose(model.intercept_, intercept, atol=1e-2)
     if coef is not None:
         assert_allclose(model.coef_[0], coef, atol=1e-2)
@@ -99,13 +110,15 @@ def test_quantile_toy_example(quantile, alpha, intercept, coef):
 
 
 @pytest.mark.parametrize("fit_intercept", [True, False])
-def test_quantile_equals_huber_for_low_epsilon(fit_intercept):
+def test_quantile_equals_huber_for_low_epsilon(fit_intercept, default_solver):
     X, y = make_regression(n_samples=100, n_features=20, random_state=0, noise=1.0)
     alpha = 1e-4
     huber = HuberRegressor(
         epsilon=1 + 1e-4, alpha=alpha, fit_intercept=fit_intercept
     ).fit(X, y)
-    quant = QuantileRegressor(alpha=alpha, fit_intercept=fit_intercept).fit(X, y)
+    quant = QuantileRegressor(
+        alpha=alpha, fit_intercept=fit_intercept, solver=default_solver
+    ).fit(X, y)
     assert_allclose(huber.coef_, quant.coef_, atol=1e-1)
     if fit_intercept:
         assert huber.intercept_ == approx(quant.intercept_, abs=1e-1)
@@ -114,26 +127,26 @@ def test_quantile_equals_huber_for_low_epsilon(fit_intercept):
 
 
 @pytest.mark.parametrize("q", [0.5, 0.9, 0.05])
-def test_quantile_estimates_calibration(q):
+def test_quantile_estimates_calibration(q, default_solver):
     # Test that model estimates percentage of points below the prediction
     X, y = make_regression(n_samples=1000, n_features=20, random_state=0, noise=1.0)
     quant = QuantileRegressor(
         quantile=q,
         alpha=0,
-        solver_options={"lstsq": False},
+        solver=default_solver,
     ).fit(X, y)
     assert np.mean(y < quant.predict(X)) == approx(q, abs=1e-2)
 
 
-def test_quantile_sample_weight():
+def test_quantile_sample_weight(default_solver):
     # test that with unequal sample weights we still estimate weighted fraction
     n = 1000
     X, y = make_regression(n_samples=n, n_features=5, random_state=0, noise=10.0)
     weight = np.ones(n)
     # when we increase weight of upper observations,
     # estimate of quantile should go up
     weight[y > y.mean()] = 100
-    quant = QuantileRegressor(quantile=0.5, alpha=1e-8, solver_options={"lstsq": False})
+    quant = QuantileRegressor(quantile=0.5, alpha=1e-8, solver=default_solver)
     quant.fit(X, y, sample_weight=weight)
     fraction_below = np.mean(y < quant.predict(X))
     assert fraction_below > 0.5
@@ -146,7 +159,7 @@ def test_quantile_sample_weight():
     reason="The `highs` solver is available from the 1.6.0 scipy version",
 )
 @pytest.mark.parametrize("quantile", [0.2, 0.5, 0.8])
-def test_asymmetric_error(quantile):
+def test_asymmetric_error(quantile, default_solver):
     """Test quantile regression for asymmetric distributed targets."""
     n_samples = 1000
     rng = np.random.RandomState(42)
@@ -171,7 +184,7 @@ def test_asymmetric_error(quantile):
     model = QuantileRegressor(
         quantile=quantile,
         alpha=0,
-        solver="highs",
+        solver=default_solver,
     ).fit(X, y)
     # This test can be made to pass with any solver but in the interest
     # of sparing continuous integration resources, the test is performed
@@ -206,7 +219,7 @@ def func(coef):
 
 
 @pytest.mark.parametrize("quantile", [0.2, 0.5, 0.8])
-def test_equivariance(quantile):
+def test_equivariance(quantile, default_solver):
     """Test equivariace of quantile regression.
 
     See Koenker (2005) Quantile Regression, Chapter 2.2.3.
@@ -223,7 +236,7 @@ def test_equivariance(quantile):
     )
     # make y asymmetric
     y += rng.exponential(scale=100, size=y.shape)
-    params = dict(alpha=0, solver_options={"lstsq": True, "tol": 1e-10})
+    params = dict(alpha=0, solver=default_solver)
     model1 = QuantileRegressor(quantile=quantile, **params).fit(X, y)
 
     # coef(q; a*y, X) = a * coef(q; y, X)
@@ -252,6 +265,7 @@ def test_equivariance(quantile):
     assert_allclose(model2.coef_, np.linalg.solve(A, model1.coef_), rtol=1e-5)
 
 
+@pytest.mark.filterwarnings("ignore:`method='interior-point'` is deprecated")
 def test_linprog_failure():
     """Test that linprog fails."""
     X = np.linspace(0, 10, num=10).reshape(-1, 1)
@@ -275,12 +289,14 @@ def test_linprog_failure():
 )
 @pytest.mark.parametrize("solver", ["highs", "highs-ds", "highs-ipm"])
 @pytest.mark.parametrize("fit_intercept", [True, False])
-def test_sparse_input(sparse_format, solver, fit_intercept):
+def test_sparse_input(sparse_format, solver, fit_intercept, default_solver):
     """Test that sparse and dense X give same results."""
     X, y = make_regression(n_samples=100, n_features=20, random_state=1, noise=1.0)
     X_sparse = sparse_format(X)
     alpha = 1e-4
-    quant_dense = QuantileRegressor(alpha=alpha, fit_intercept=fit_intercept).fit(X, y)
+    quant_dense = QuantileRegressor(
+        alpha=alpha, fit_intercept=fit_intercept, solver=default_solver
+    ).fit(X, y)
     quant_sparse = QuantileRegressor(
         alpha=alpha, fit_intercept=fit_intercept, solver=solver
     ).fit(X_sparse, y)
@@ -289,3 +305,26 @@ def test_sparse_input(sparse_format, solver, fit_intercept):
         assert quant_sparse.intercept_ == approx(quant_dense.intercept_)
         # check that we still predict fraction
         assert 0.45 <= np.mean(y < quant_sparse.predict(X_sparse)) <= 0.55
+
+
+# TODO (1.4): remove this test in 1.4
+def test_warning_new_default(X_y_data):
+    """Check that we warn about the new default solver."""
+    X, y = X_y_data
+    model = QuantileRegressor()
+    with pytest.warns(FutureWarning, match="The default solver will change"):
+        model.fit(X, y)
+
+
+def test_error_interior_point_future(X_y_data, monkeypatch):
+    """Check that we will raise a proper error when requesting
+    `solver='interior-point'` in SciPy >= 1.11.
+    """
+    X, y = X_y_data
+    import sklearn.linear_model._quantile
+
+    with monkeypatch.context() as m:
+        m.setattr(sklearn.linear_model._quantile, "sp_version", parse_version("1.11.0"))
+        err_msg = "Solver interior-point is not anymore available in SciPy >= 1.11.0."
+        with pytest.raises(ValueError, match=err_msg):
+            QuantileRegressor(solver="interior-point").fit(X, y)
diff --git a/sklearn/tests/test_docstring_parameters.py b/sklearn/tests/test_docstring_parameters.py
@@ -20,6 +20,7 @@
 from sklearn.utils.estimator_checks import _enforce_estimator_tags_y
 from sklearn.utils.estimator_checks import _enforce_estimator_tags_x
 from sklearn.utils.estimator_checks import _construct_instance
+from sklearn.utils.fixes import sp_version, parse_version
 from sklearn.utils.deprecation import _is_deprecated
 from sklearn.datasets import make_classification
 from sklearn.linear_model import LogisticRegression
@@ -266,6 +267,11 @@ def test_fit_docstring_attributes(name, Estimator):
     if Estimator.__name__ in ("KMeans", "MiniBatchKMeans"):
         est.set_params(n_init="auto")
 
+    # TODO(1.4): TO BE REMOVED for 1.4 (avoid FutureWarning)
+    if Estimator.__name__ == "QuantileRegressor":
+        solver = "highs" if sp_version >= parse_version("1.6.0") else "interior-point"
+        est.set_params(solver=solver)
+
     # In case we want to deprecate some attributes in the future
     skipped_attributes = {}