[MRG] Add elastic net penalty to LogisticRegression (scikit-learn#11646)

* First draft on elasticnet penaly for LogisticRegression * Some basic tests * Doc update * First draft for LogisticRegressionCV. It seems to be working for binary classification and for multiclass when multi_class='ovr'. I'm having a hard time figuring out the intricacies of multi_class='multinomial'. * Changed default to None for l1_ratio. added warning message is user sets l1_ratio while penalty is not elastic-net * Some more doc * Updated example to plot elastic net sparsity * Fixed flake8 * Fixed test by not modifying attribute in fit * Fixed doc issues * WIP * Partially fixed logistic_reg_CV for multinomial. Also added some comments that are hopefully clear. Still need to fix refit=False * Fixed doc issue * WIP * Fixed test for refit=False in LogisticRegressionCV * Fixed Python 2 numpy version issue * minor doc updates * Weird doc error... * Added test to ensure that elastic net is at least as good as L1 or L2 once l1_ratio has been optimized with grid search Also addressed minor reviews * Fixed test * addressed comments * Added back ignore warning on tests * Added a functional test * Scale data in test... Now failing * elastic-net --> elasticnet * Updated doc for some attributes and checked their shape in tests * Added l1_ratio dimension to coefs_paths and scores attr * improve example + fix test * FIX incorrect lagged_update in SAGA * Add non-regression test for SAGA's bug * FIX flake8 and warning * Re fixed warning * Updated some tests * Addressed comments * more comments and added dimension to LogisticRegressionCV.n_iter_ attribute * Updated whatsnew for 0.21 * better doc shape looks * Fixed whatnew entry after merges * Added dot * Addressed comments + standardized optional default param docstrings * Addessed comments * use swapaxes instead of unsupported moveaxis (hopefully fixes tests)

Author: NicolasHug

doc/modules/linear_model.rst

@@ -338,7 +338,7 @@ the algorithm to fit the coefficients.
 
 .. _elastic_net:
 
-Elastic Net
+Elastic-Net
 ===========
 :class:`ElasticNet` is a linear regression model trained with L1 and L2 prior
 as regularizer. This combination allows for learning a sparse model where
@@ -390,7 +390,7 @@ the duality gap computation used for convergence control.
 
 .. _multi_task_elastic_net:
 
-Multi-task Elastic Net
+Multi-task Elastic-Net
 ======================
 
 The :class:`MultiTaskElasticNet` is an elastic-net model that estimates sparse
@@ -730,7 +730,7 @@ or the log-linear classifier. In this model, the probabilities describing the po
 
 The implementation of logistic regression in scikit-learn can be accessed from
 class :class:`LogisticRegression`. This implementation can fit binary, One-vs-
-Rest, or multinomial logistic regression with optional L2 or L1
+Rest, or multinomial logistic regression with optional L2, L1 or Elastic-Net
 regularization.
 
 As an optimization problem, binary class L2 penalized logistic regression
@@ -739,12 +739,22 @@ minimizes the following cost function:
 .. math:: \min_{w, c} \frac{1}{2}w^T w + C \sum_{i=1}^n \log(\exp(- y_i (X_i^T w + c)) + 1) .
 
 Similarly, L1 regularized logistic regression solves the following
-optimization problem
+optimization problem:
 
 .. math:: \min_{w, c} \|w\|_1 + C \sum_{i=1}^n \log(\exp(- y_i (X_i^T w + c)) + 1).
 
+Elastic-Net regularization is a combination of L1 and L2, and minimizes the
+following cost function:
+
+.. math:: \min_{w, c} \frac{1 - \rho}{2}w^T w + \rho \|w\|_1 + C \sum_{i=1}^n \log(\exp(- y_i (X_i^T w + c)) + 1),
+
+where :math:`\rho` controls the strengh of L1 regularization vs L2
+regularization (it corresponds to the `l1_ratio` parameter).
+
 Note that, in this notation, it's assumed that the observation :math:`y_i` takes values in the set
-:math:`{-1, 1}` at trial :math:`i`.
+:math:`{-1, 1}` at trial :math:`i`. We can also see that Elastic-Net is
+equivalent to L1 when :math:`\rho = 1` and equivalent to L2 when
+:math:`\rho=0`.
 
 The solvers implemented in the class :class:`LogisticRegression`
 are "liblinear", "newton-cg", "lbfgs", "sag" and "saga":
@@ -772,10 +782,12 @@ than other solvers for large datasets, when both the number of samples and the
 number of features are large.
 
 The "saga" solver [7]_ is a variant of "sag" that also supports the
-non-smooth `penalty="l1"` option. This is therefore the solver of choice
-for sparse multinomial logistic regression.
+non-smooth `penalty="l1"`. This is therefore the solver of choice for sparse
+multinomial logistic regression. It is also the only solver that supports
+`penalty="elasticnet"`.
 
-In a nutshell, the following table summarizes the penalties supported by each solver:
+In a nutshell, the following table summarizes the penalties supported by
+each solver:
 
 +------------------------------+-----------------+-------------+-----------------+-----------+------------+
 | | **Solvers** |
@@ -790,6 +802,8 @@ In a nutshell, the following table summarizes the penalties supported by each so
 +------------------------------+-----------------+-------------+-----------------+-----------+------------+
 | OVR + L1 penalty | yes | no | no | no | yes |
 +------------------------------+-----------------+-------------+-----------------+-----------+------------+
+| Elastic-Net | no | no | no | no | yes |
++------------------------------+-----------------+-------------+-----------------+-----------+------------+
 | **Behaviors** | |
 +------------------------------+-----------------+-------------+-----------------+-----------+------------+
 | Penalize the intercept (bad) | yes | no | no | no | no |
@@ -799,8 +813,8 @@ In a nutshell, the following table summarizes the penalties supported by each so
 | Robust to unscaled datasets | yes | yes | yes | no | no |
 +------------------------------+-----------------+-------------+-----------------+-----------+------------+
 
-The "saga" solver is often the best choice but requires scaling. The "liblinear" solver is
-used by default for historical reasons.
+The "saga" solver is often the best choice but requires scaling. The
+"liblinear" solver is used by default for historical reasons.
 
 For large dataset, you may also consider using :class:`SGDClassifier`
 with 'log' loss.
@@ -838,14 +852,11 @@ with 'log' loss.
  thus be used to perform feature selection, as detailed in
  :ref:`l1_feature_selection`.
 
-:class:`LogisticRegressionCV` implements Logistic Regression with
-builtin cross-validation to find out the optimal C parameter.
-"newton-cg", "sag", "saga" and "lbfgs" solvers are found to be faster
-for high-dimensional dense data, due to warm-starting. For the
-multiclass case, if `multi_class` option is set to "ovr", an optimal C
-is obtained for each class and if the `multi_class` option is set to
-"multinomial", an optimal C is obtained by minimizing the cross-entropy
-loss.
+:class:`LogisticRegressionCV` implements Logistic Regression with built-in
+cross-validation support, to find the optimal `C` and `l1_ratio` parameters
+according to the ``scoring`` attribute. The "newton-cg", "sag", "saga" and
+"lbfgs" solvers are found to be faster for high-dimensional dense data, due
+to warm-starting (see :term:`Glossary <warm_start>`).
 
 .. topic:: References:
 

doc/tutorial/statistical_inference/supervised_learning.rst

@@ -183,6 +183,7 @@ Linear models: :math:`y = X\beta + \epsilon`
  [ 0.30349955 -237.63931533 510.53060544 327.73698041 -814.13170937
  492.81458798 102.84845219 184.60648906 743.51961675 76.09517222]
 
+
  >>> # The mean square error
  >>> np.mean((regr.predict(diabetes_X_test) - diabetes_y_test)**2)
  ... # doctest: +ELLIPSIS
@@ -378,7 +379,7 @@ function or **logistic** function:
  ... multi_class='multinomial')
  >>> log.fit(iris_X_train, iris_y_train) # doctest: +NORMALIZE_WHITESPACE
  LogisticRegression(C=100000.0, class_weight=None, dual=False,
- fit_intercept=True, intercept_scaling=1, max_iter=100,
+ fit_intercept=True, intercept_scaling=1, l1_ratio=None, max_iter=100,
  multi_class='multinomial', n_jobs=None, penalty='l2', random_state=None,
  solver='lbfgs', tol=0.0001, verbose=0, warm_start=False)
 

doc/whats_new/v0.21.rst

@@ -22,6 +22,9 @@ random sampling procedures.
 
 - Decision trees and derived ensembles when both `max_depth` and
  `max_leaf_nodes` are set. |Fix|
+- :class:`linear_model.LogisticRegression` and
+ :class:`linear_model.LogisticRegressionCV` with 'saga' solver. |Fix|
+
 
 Details are listed in the changelog below.
 
@@ -146,6 +149,15 @@ Support for Python 3.4 and below has been officially dropped.
  affects all ensemble methods using decision trees.
  :pr:`12344` by :user:`Adrin Jalali <adrinjalali>`.
 
+:mod:`sklearn.linear_model`
+...........................
+
+- |Feature| :class:`linear_model.LogisticRegression` and
+ :class:`linear_model.LogisticRegressionCV` now support Elastic-Net penalty,
+ with the 'saga' solver. :issue:`11646` by :user:`Nicolas Hug <NicolasHug>`.
+
+- |Fix| Fixed a bug in the 'saga' solver where the weights would not be
+ correctly updated in some cases. :issue:`11646` by `Tom Dupre la Tour`_.
 
 Multiple modules
 ................

examples/linear_model/plot_logistic_l1_l2_sparsity.py

@@ -4,10 +4,11 @@
 ==============================================
 
 Comparison of the sparsity (percentage of zero coefficients) of solutions when
-L1 and L2 penalty are used for different values of C. We can see that large
-values of C give more freedom to the model. Conversely, smaller values of C
-constrain the model more. In the L1 penalty case, this leads to sparser
-solutions.
+L1, L2 and Elastic-Net penalty are used for different values of C. We can see
+that large values of C give more freedom to the model. Conversely, smaller
+values of C constrain the model more. In the L1 penalty case, this leads to
+sparser solutions. As expected, the Elastic-Net penalty sparsity is between
+that of L1 and L2.
 
 We classify 8x8 images of digits into two classes: 0-4 against 5-9.
 The visualization shows coefficients of the models for varying C.
@@ -35,45 +36,55 @@
 # classify small against large digits
 y = (y > 4).astype(np.int)
 
+l1_ratio = 0.5 # L1 weight in the Elastic-Net regularization
+
+fig, axes = plt.subplots(3, 3)
 
 # Set regularization parameter
-for i, C in enumerate((1, 0.1, 0.01)):
+for i, (C, axes_row) in enumerate(zip((1, 0.1, 0.01), axes)):
  # turn down tolerance for short training time
  clf_l1_LR = LogisticRegression(C=C, penalty='l1', tol=0.01, solver='saga')
  clf_l2_LR = LogisticRegression(C=C, penalty='l2', tol=0.01, solver='saga')
+ clf_en_LR = LogisticRegression(C=C, penalty='elasticnet', solver='saga',
+ l1_ratio=l1_ratio, tol=0.01)
  clf_l1_LR.fit(X, y)
  clf_l2_LR.fit(X, y)
+ clf_en_LR.fit(X, y)
 
  coef_l1_LR = clf_l1_LR.coef_.ravel()
  coef_l2_LR = clf_l2_LR.coef_.ravel()
+ coef_en_LR = clf_en_LR.coef_.ravel()
 
  # coef_l1_LR contains zeros due to the
  # L1 sparsity inducing norm
 
  sparsity_l1_LR = np.mean(coef_l1_LR == 0) * 100
  sparsity_l2_LR = np.mean(coef_l2_LR == 0) * 100
+ sparsity_en_LR = np.mean(coef_en_LR == 0) * 100
 
  print("C=%.2f" % C)
- print("Sparsity with L1 penalty: %.2f%%" % sparsity_l1_LR)
- print("score with L1 penalty: %.4f" % clf_l1_LR.score(X, y))
- print("Sparsity with L2 penalty: %.2f%%" % sparsity_l2_LR)
- print("score with L2 penalty: %.4f" % clf_l2_LR.score(X, y))
+ print("{:<40} {:.2f}%".format("Sparsity with L1 penalty:", sparsity_l1_LR))
+ print("{:<40} {:.2f}%".format("Sparsity with Elastic-Net penalty:",
+ sparsity_en_LR))
+ print("{:<40} {:.2f}%".format("Sparsity with L2 penalty:", sparsity_l2_LR))
+ print("{:<40} {:.2f}".format("Score with L1 penalty:",
+ clf_l1_LR.score(X, y)))
+ print("{:<40} {:.2f}".format("Score with Elastic-Net penalty:",
+ clf_en_LR.score(X, y)))
+ print("{:<40} {:.2f}".format("Score with L2 penalty:",
+ clf_l2_LR.score(X, y)))
 
- l1_plot = plt.subplot(3, 2, 2 * i + 1)
- l2_plot = plt.subplot(3, 2, 2 * (i + 1))
  if i == 0:
- l1_plot.set_title("L1 penalty")
- l2_plot.set_title("L2 penalty")
-
- l1_plot.imshow(np.abs(coef_l1_LR.reshape(8, 8)), interpolation='nearest',
- cmap='binary', vmax=1, vmin=0)
- l2_plot.imshow(np.abs(coef_l2_LR.reshape(8, 8)), interpolation='nearest',
- cmap='binary', vmax=1, vmin=0)
- plt.text(-8, 3, "C = %.2f" % C)
-
- l1_plot.set_xticks(())
- l1_plot.set_yticks(())
- l2_plot.set_xticks(())
- l2_plot.set_yticks(())
+ axes_row[0].set_title("L1 penalty")
+ axes_row[1].set_title("Elastic-Net\nl1_ratio = %s" % l1_ratio)
+ axes_row[2].set_title("L2 penalty")
+
+ for ax, coefs in zip(axes_row, [coef_l1_LR, coef_en_LR, coef_l2_LR]):
+ ax.imshow(np.abs(coefs.reshape(8, 8)), interpolation='nearest',
+ cmap='binary', vmax=1, vmin=0)
+ ax.set_xticks(())
+ ax.set_yticks(())
+
+ axes_row[0].set_ylabel('C = %s' % C)
 
 plt.show()