scikit-learn · Micky774 · Jul 9, 2020 · Jul 9, 2020 · Jul 23, 2020 · Jul 24, 2020
diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst
@@ -974,36 +974,37 @@ small data-sets but for larger datasets its performance suffers. [9]_
 
 The following table summarizes the penalties supported by each solver:
 
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-|                              |                       **Solvers**                                        |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| **Penalties**                | **'liblinear'** | **'lbfgs'** | **'newton-cg'** | **'sag'** | **'saga'** |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| Multinomial + L2 penalty     |       no        |     yes     |       yes       |    yes    |    yes     |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| OVR + L2 penalty             |       yes       |     yes     |       yes       |    yes    |    yes     |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| Multinomial + L1 penalty     |       no        |     no      |       no        |    no     |    yes     |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| OVR + L1 penalty             |       yes       |     no      |       no        |    no     |    yes     |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| Elastic-Net                  |       no        |     no      |       no        |    no     |    yes     |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| No penalty ('none')          |       no        |     yes     |       yes       |    yes    |    yes     |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| **Behaviors**                |                                                                          |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| Penalize the intercept (bad) |       yes       |     no      |       no        |    no     |    no      |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| Faster for large datasets    |       no        |     no      |       no        |    yes    |    yes     |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
-| Robust to unscaled datasets  |       yes       |     yes     |       yes       |    no     |    no      |
-+------------------------------+-----------------+-------------+-----------------+-----------+------------+
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+|                              |                       **Solvers**                                                          |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| **Penalties**                | **'liblinear'** | **'lbfgs'** | **'newton-cg'** | **'sag'** | **'saga'** | **'trust-ncg'** |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| Multinomial + L2 penalty     |       no        |     yes     |       yes       |    yes    |    yes     |       yes       |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| OVR + L2 penalty             |       yes       |     yes     |       yes       |    yes    |    yes     |       yes       |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| Multinomial + L1 penalty     |       no        |     no      |       no        |    no     |    yes     |       no        |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| OVR + L1 penalty             |       yes       |     no      |       no        |    no     |    yes     |       no        |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| Elastic-Net                  |       no        |     no      |       no        |    no     |    yes     |       no        |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| No penalty ('none')          |       no        |     yes     |       yes       |    yes    |    yes     |       yes       |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| **Behaviors**                |                                                                          |                 |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| Penalize the intercept (bad) |       yes       |     no      |       no        |    no     |    no      |       no        |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| Faster for large datasets    |       no        |     no      |       no        |    yes    |    yes     |       yes       |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
+| Robust to unscaled datasets  |       yes       |     yes     |       yes       |    no     |    no      |       yes       |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+-----------------+
 
 The "lbfgs" solver is used by default for its robustness. For large datasets
 the "saga" solver is usually faster.
 For large dataset, you may also consider using :class:`SGDClassifier`
 with 'log' loss, which might be even faster but requires more tuning.
+The "trust-ncg" solver is useful for sparse datasets.
 
 .. topic:: Examples:
 

diff --git a/doc/whats_new/v1.2.rst b/doc/whats_new/v1.2.rst
@@ -142,6 +142,15 @@ Changelog
 - |Fix| Use dtype-aware tolerances for the validation of gram matrices (passed by users
   or precomputed). :pr:`22059` by :user:`Malte S. Kurz <MalteKurz>`.
 
+- |Enhancement| Added `trust-ncg` solver to `linear_model.LogisticRegression`.
+  This solver is fast for large datasets, similar to `saga`, while also being
+  robust to unscaled data, similar to `lbfgs`. The solver may require a higher
+  value of `max_iter` for proper convergence. It is compatible with both `ovr`
+  and `multinomial` options for `multi_class`, and can optimize with an `L2`
+  penalty.
+  :pr:`17877` by :user:`Rithvik Rao <rithvikrao>` and :user:`Ruby Werman <rubywerman>`,
+  and :pr:`22236` by :user:`Meekail Zain <micky774>`.
+
 - |Fix| Fixed an error in :class:`linear_model.LogisticRegression` with
   `solver="newton-cg"`, `fit_intercept=True`, and a single feature. :pr:`23608`
   by `Tom Dupre la Tour`_.

diff --git a/sklearn/linear_model/_logistic.py b/sklearn/linear_model/_logistic.py
@@ -43,7 +43,7 @@
 
 
 def _check_solver(solver, penalty, dual):
-    all_solvers = ["liblinear", "newton-cg", "lbfgs", "sag", "saga"]
+    all_solvers = ["liblinear", "newton-cg", "lbfgs", "sag", "saga", "trust-ncg"]
     if solver not in all_solvers:
         raise ValueError(
             "Logistic Regression supports only solvers in %s, got %s."
@@ -158,8 +158,8 @@ def _logistic_regression_path(
         Maximum number of iterations for the solver.
 
     tol : float, default=1e-4
-        Stopping criterion. For the newton-cg and lbfgs solvers, the iteration
-        will stop when ``max{|g_i | i = 1, ..., n} <= tol``
+        Stopping criterion. For the newton-cg, trust-ncg, and lbfgs solvers,
+        the iteration will stop when ``max{|g_i | i = 1, ..., n} <= tol``
         where ``g_i`` is the i-th component of the gradient.
 
     verbose : int, default=0
@@ -316,7 +316,7 @@ def _logistic_regression_path(
         w0 = np.zeros(n_features + int(fit_intercept), dtype=X.dtype)
         mask = y == pos_class
         y_bin = np.ones(y.shape, dtype=X.dtype)
-        if solver in ["lbfgs", "newton-cg"]:
+        if solver in ["lbfgs", "newton-cg", "trust-ncg"]:
             # HalfBinomialLoss, used for those solvers, represents y in [0, 1] instead
             # of in [-1, 1].
             mask_classes = np.array([0, 1])
@@ -333,7 +333,7 @@ def _logistic_regression_path(
             sample_weight *= class_weight_[le.fit_transform(y_bin)]
 
     else:
-        if solver in ["sag", "saga", "lbfgs", "newton-cg"]:
+        if solver in ["sag", "saga", "lbfgs", "newton-cg", "trust-ncg"]:
             # SAG, lbfgs and newton-cg multinomial solvers need LabelEncoder,
             # not LabelBinarizer, i.e. y as a 1d-array of integers.
             # LabelEncoder also saves memory compared to LabelBinarizer, especially
@@ -391,7 +391,7 @@ def _logistic_regression_path(
                 w0[:, : coef.shape[1]] = coef
 
     if multi_class == "multinomial":
-        if solver in ["lbfgs", "newton-cg"]:
+        if solver in ["lbfgs", "trust-ncg", "newton-cg"]:
             # scipy.optimize.minimize and newton-cg accept only ravelled parameters,
             # i.e. 1d-arrays. LinearModelLoss expects classes to be contiguous and
             # reconstructs the 2d-array via w0.reshape((n_classes, -1), order="F").
@@ -408,6 +408,17 @@ def _logistic_regression_path(
             func = loss.loss
             grad = loss.gradient
             hess = loss.gradient_hessian_product  # hess = [gradient, hessp]
+        # TODO: remove local LinearModelLoss after renaming `loss`
+        elif solver == "trust-ncg":
+            loss_ = LinearModelLoss(
+                base_loss=HalfMultinomialLoss(n_classes=classes.size),
+                fit_intercept=fit_intercept,
+            )
+            func = loss_.loss_gradient
+
+            def hessp(x, p, *args):
+                return loss_.gradient_hessian_product(x, *args)[1](p)
+
         warm_start_sag = {"coef": w0.T}
     else:
         target = y_bin
@@ -423,24 +434,45 @@ def _logistic_regression_path(
             func = loss.loss
             grad = loss.gradient
             hess = loss.gradient_hessian_product  # hess = [gradient, hessp]
+        # TODO: remove local LinearModelLoss after renaming `loss`
+        elif solver == "trust-ncg":
+            loss_ = LinearModelLoss(
+                base_loss=HalfBinomialLoss(), fit_intercept=fit_intercept
+            )
+            func = loss_.loss_gradient
+
+            def hessp(x, p, *args):
+                return loss_.gradient_hessian_product(x, *args)[1](p)
+
         warm_start_sag = {"coef": np.expand_dims(w0, axis=1)}
 
     coefs = list()
     n_iter = np.zeros(len(Cs), dtype=np.int32)
     for i, C in enumerate(Cs):
-        if solver == "lbfgs":
-            l2_reg_strength = 1.0 / C
-            iprint = [-1, 50, 1, 100, 101][
-                np.searchsorted(np.array([0, 1, 2, 3]), verbose)
-            ]
-            opt_res = optimize.minimize(
-                func,
-                w0,
-                method="L-BFGS-B",
-                jac=True,
-                args=(X, target, sample_weight, l2_reg_strength, n_threads),
-                options={"iprint": iprint, "gtol": tol, "maxiter": max_iter},
-            )
+        if solver in ["lbfgs", "trust-ncg"]:
+            if solver == "lbfgs":
+                l2_reg_strength = 1.0 / C
+                iprint = [-1, 50, 1, 100, 101][
+                    np.searchsorted(np.array([0, 1, 2, 3]), verbose)
+                ]
+                opt_res = optimize.minimize(
+                    func,
+                    w0,
+                    method="L-BFGS-B",
+                    jac=True,
+                    args=(X, target, sample_weight, l2_reg_strength, n_threads),
+                    options={"iprint": iprint, "gtol": tol, "maxiter": max_iter},
+                )
+            elif solver == "trust-ncg":
+                opt_res = optimize.minimize(
+                    func,
+                    w0,
+                    method=solver,
+                    jac=True,
+                    hessp=hessp,
+                    args=(X, target, sample_weight, 1.0 / C),
+                    options={"gtol": tol, "maxiter": max_iter},
+                )
             n_iter_i = _check_optimize_result(
                 solver,
                 opt_res,
@@ -517,7 +549,7 @@ def _logistic_regression_path(
 
         if multi_class == "multinomial":
             n_classes = max(2, classes.size)
-            if solver in ["lbfgs", "newton-cg"]:
+            if solver in ["lbfgs", "newton-cg", "trust-ncg"]:
                 multi_w0 = np.reshape(w0, (n_classes, -1), order="F")
             else:
                 multi_w0 = w0
@@ -842,7 +874,7 @@ class LogisticRegression(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
         Used when ``solver`` == 'sag', 'saga' or 'liblinear' to shuffle the
         data. See :term:`Glossary <random_state>` for details.
 
-    solver : {'newton-cg', 'lbfgs', 'liblinear', 'sag', 'saga'}, \
+    solver : {'newton-cg', 'lbfgs', 'liblinear', 'sag', 'saga', 'trust-ncg'}, \
             default='lbfgs'
 
         Algorithm to use in the optimization problem. Default is 'lbfgs'.
@@ -863,6 +895,7 @@ class LogisticRegression(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
            - 'liblinear'   -   ['l1', 'l2']
            - 'sag'         -   ['l2', 'none']
            - 'saga'        -   ['elasticnet', 'l1', 'l2', 'none']
+           - 'trust-ncg'   -   ['l2', 'none']
 
         .. note::
            'sag' and 'saga' fast convergence is only guaranteed on
@@ -881,6 +914,8 @@ class LogisticRegression(LinearClassifierMixin, SparseCoefMixin, BaseEstimator):
            SAGA solver.
         .. versionchanged:: 0.22
             The default solver changed from 'liblinear' to 'lbfgs' in 0.22.
+        .. versionchanged:: 1.1
+            Newton conjugate gradient trust-region (trust-ncg) solver.
 
     max_iter : int, default=100
         Maximum number of iterations taken for the solvers to converge.