Repository: scikit-learn/scikit-learn
Pull request: [MRG] Add elastic net penalty to LogisticRegression
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/master/CONTRIBUTING.md#pull-request-checklist
-->

#### Reference Issues/PRs
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->

Closes #8288

#### What does this implement/fix? Explain your changes.

This PR allows the use of `penalty='elastic-net'` for the `LogisticRegression` class. Elastic net is already available in the saga solver, it's just not exposed yet.

#### Any other comments?

~~I set the `l1_ratio` parameter to 0.5 by default, because that's the default in `ElasticNet`. Should I set it to `None`? This would allow to raise an error if it's set to something else while `penalty != 'elastic-net'`.~~ 
EDIT: I set it to None and issue warning if penalty != elastic-net

#### TODO

- [X] expose `penalty='elastic-net'` in `LogisticRegression`
- [X] expose `penalty='elastic-net'` in `LogisticRegressionCV`
- [X] Update doc
- [X] write tests
- [X] Update examples

<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
linear classifier with elasticnet penalty
I think it's weird that we don't have a linear classifier with an elasticnet penalty (apart from SGDClassifier which is kinda hard to use). Maybe we should add that to the lbfgs solver in logistic regression?
glmnet does that afaik.

I'm not sure if I opened this issue before... I got a slight deja-vu.

Hints from issue discussion:
Are you suggesting doing the trick of splitting the weights into a positive and negative part and using the bounding box constraints? Otherwise l-bfgs is not an usual choice for sparse penalties.

On a related note, there are 4 orthogonal attributes of a classifier that are very difficult to find supported together:

1. Loss (e.g. log, squared hinge)
2. True multiclass (not ovr)
3. Elastic net regularization
4. Class/sample weights.

The only python solver I'm aware of that supports all of these is at the moment is SAGA in lightning.

It's actually frustrating in practice because i keep forgetting, and when I realize I need one of these features on a dataset, I need to rewrite a lot of code (esp. if using, e.g. LogisticRegressionCV before) and change the optimizer.
Can I work on this.
Could it be something like `sklearn.linear_model.RidgeClassifier` but with the elasticnet penalty?
I've not been following #8446. might it have elastic net?

On 3 Mar 2017 2:19 am, "Samuele Fiorini" <notifications@github.com> wrote:

> Could it be something like sklearn.linear_model.RidgeClassifier but with
> the elasticnet penalty?
>
> —
> You are receiving this because you are subscribed to this thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/scikit-learn/scikit-learn/issues/8288#issuecomment-283682213>,
> or mute the thread
> <https://github.com/notifications/unsubscribe-auth/AAEz6wXTPm4wyqn7xn8c2m09uswU7lHfks5rht4LgaJpZM4L22__>
> .
>

> apart from SGDClassifier which is kinda hard to use

We need to add a good stopping criterion for SGDClassifier. Right now it is a fixed `n_iter` instead of the usual `tol` + `max_iter` combo. Reusing the stopping criterion of the `MLPClassifer` would render the `SGDClassifier` much easier to use.
As mentioned by @jnothman, I have some hope that #8446 will enable us to tackle this.
given that #8446 is merged, this should be up for grabs, right?
My interpretation of @rhiever's results [here](https://arxiv.org/abs/1708.05070) is that the mixed penalty really helps - I don't know why else SGDClassifier would beat logistic regression so much. Though SGDClassifier is allowed to search over the loss as well, I guess?
@arthurmensch do you have any motivation to work on this? Shouldn't be too hard to expose, right?

For anyone else who wants to work on this: it's already implemented in SAGA, it's just not exposed in LogisticRegression....
If no one is currently working on that, I would like to try to do it.

Commit message: Fixed test

Patch to review:
diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst
index ab6b2994835f9..c8ff40a6ee1ae 100644
--- a/doc/modules/linear_model.rst
+++ b/doc/modules/linear_model.rst
@@ -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 the 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,8 +782,9 @@ 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="elastic-net"`.
 
 In a nutshell, the following table summarizes the solvers characteristics:
 
@@ -784,6 +795,7 @@ 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 penalty            no           no       no           no     yes
 ============================   ===========  =======  ===========  =====  ======
 Penalize the intercept (bad)   yes          no       no           no     no
 Faster for large datasets      no           no       no           yes    yes
@@ -830,13 +842,13 @@ with 'log' loss.
    :ref:`l1_feature_selection`.
 
 :class:`LogisticRegressionCV` implements Logistic Regression with
-builtin cross-validation to find out the optimal C parameter.
+builtin cross-validation to find out the optimal C and l1_ratio parameters.
 "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.
+multiclass case, if `multi_class` option is set to "ovr", an optimal C and
+l1_ratio is obtained for each class and if the `multi_class` option is set to
+"multinomial", an optimal C and l1_ratio is obtained by minimizing the
+cross-entropy loss.
 
 .. topic:: References:
 
diff --git a/doc/tutorial/statistical_inference/supervised_learning.rst b/doc/tutorial/statistical_inference/supervised_learning.rst
index 2c12f1038c285..4f050e045aa86 100644
--- a/doc/tutorial/statistical_inference/supervised_learning.rst
+++ b/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)
 
diff --git a/examples/linear_model/plot_logistic_l1_l2_sparsity.py b/examples/linear_model/plot_logistic_l1_l2_sparsity.py
index be63b144c260a..4eb09324c175f 100644
--- a/examples/linear_model/plot_logistic_l1_l2_sparsity.py
+++ b/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.
@@ -41,39 +42,56 @@
     # turn down tolerance for short training time
     clf_l1_LR = LogisticRegression(C=C, penalty='l1', tol=0.01)
     clf_l2_LR = LogisticRegression(C=C, penalty='l2', tol=0.01)
+    clf_EN_LR = LogisticRegression(C=C, l1_ratio=.3, penalty='elastic-net',
+                                   solver='saga', 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))
-
-    l1_plot = plt.subplot(3, 2, 2 * i + 1)
-    l2_plot = plt.subplot(3, 2, 2 * (i + 1))
+    print("{:<40} {:.2f}%".format("Sparsity with L1 penalty:", sparsity_l1_LR))
+    print("{:<40} {:.2f}%".format("Sparsity with L2 penalty:", sparsity_l2_LR))
+    print("{:<40} {:.2f}%".format("Sparsity with elastic-net penalty:",
+                                  sparsity_EN_LR))
+    print("{:<40} {:.2f}".format("Score with L1 penalty:",
+                                 clf_l1_LR.score(X, y)))
+    print("{:<40} {:.2f}".format("Score with L2 penalty:",
+                                 clf_l2_LR.score(X, y)))
+    print("{:<40} {:.2f}".format("Score with elastic-net penalty:",
+                                 clf_EN_LR.score(X, y)))
+
+    l1_plot = plt.subplot(3, 3, 3 * i + 1)
+    l2_plot = plt.subplot(3, 3, 3 * i + 2)
+    en_plot = plt.subplot(3, 3, 3 * i + 3)
     if i == 0:
         l1_plot.set_title("L1 penalty")
         l2_plot.set_title("L2 penalty")
+        en_plot.set_title("Elastic-net with l1_ratio = 0.3")
 
     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)
+    en_plot.imshow(np.abs(coef_EN_LR.reshape(8, 8)), interpolation='nearest',
+                   cmap='binary', vmax=1, vmin=0)
+    plt.text(-40, 3, "C = %.2f" % C)
 
     l1_plot.set_xticks(())
     l1_plot.set_yticks(())
     l2_plot.set_xticks(())
     l2_plot.set_yticks(())
+    en_plot.set_xticks(())
+    en_plot.set_yticks(())
 
 plt.show()
diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py
index 01a4f78ab0157..40f8a66bf30bf 100644
--- a/sklearn/linear_model/logistic.py
+++ b/sklearn/linear_model/logistic.py
@@ -436,7 +436,7 @@ def _check_solver(solver, penalty, dual):
         raise ValueError("Logistic Regression supports only solvers in %s, got"
                          " %s." % (all_solvers, solver))
 
-    all_penalties = ['l1', 'l2']
+    all_penalties = ['l1', 'l2', 'elastic-net']
     if penalty not in all_penalties:
         raise ValueError("Logistic Regression supports only penalties in %s,"
                          " got %s." % (all_penalties, penalty))
@@ -447,6 +447,10 @@ def _check_solver(solver, penalty, dual):
     if solver != 'liblinear' and dual:
         raise ValueError("Solver %s supports only "
                          "dual=False, got dual=%s" % (solver, dual))
+
+    if penalty == 'elastic-net' and solver != 'saga':
+        raise ValueError("Only 'saga' solver supports elastic-net penalty,"
+                         " got solver={}.".format(solver))
     return solver
 
 
@@ -473,14 +477,14 @@ def _check_multi_class(multi_class, solver, n_classes):
     return multi_class
 
 
-
 def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True,
                              max_iter=100, tol=1e-4, verbose=0,
                              solver='lbfgs', coef=None,
                              class_weight=None, dual=False, penalty='l2',
                              intercept_scaling=1., multi_class='warn',
                              random_state=None, check_input=True,
-                             max_squared_sum=None, sample_weight=None):
+                             max_squared_sum=None, sample_weight=None,
+                             l1_ratio=None):
     """Compute a Logistic Regression model for a list of regularization
     parameters.
 
@@ -549,9 +553,10 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True,
         l2 penalty with liblinear solver. Prefer dual=False when
         n_samples > n_features.
 
-    penalty : str, 'l1' or 'l2'
+    penalty : str, 'l1', 'l2', or 'elastic-net'
         Used to specify the norm used in the penalization. The 'newton-cg',
-        'sag' and 'lbfgs' solvers support only l2 penalties.
+        'sag' and 'lbfgs' solvers support only l2 penalties. 'elastic-net' is
+        only supported by the 'saga' solver.
 
     intercept_scaling : float, default 1.
         Useful only when the solver 'liblinear' is used
@@ -599,6 +604,13 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True,
         Array of weights that are assigned to individual samples.
         If not provided, then each sample is given unit weight.
 
+    l1_ratio : float or None, optional (default=None)
+        The elastic net mixing parameter, with ``0 <= l1_ratio <= 1``. Only
+        used if ``penalty='elastic-net'``. Setting ``l1_ratio=0`` is equivalent
+        to using ``penalty='l2'``, while setting ``l1_ratio=1`` is equivalent
+        to using ``penalty='l1'``. For ``0 < l1_ratio <1``, the penalty is a
+        combination of L1 and L2.
+
     Returns
     -------
     coefs : ndarray, shape (n_cs, n_features) or (n_cs, n_features + 1)
@@ -781,9 +793,13 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True,
             if penalty == 'l1':
                 alpha = 0.
                 beta = 1. / C
-            else:
+            elif penalty == 'l2':
                 alpha = 1. / C
                 beta = 0.
+            else:  # elastic-net penalty (alpha is 2-norm, beta is 1-norm)
+                alpha = (1. / C) * (1 - l1_ratio)
+                beta = (1. / C) * l1_ratio
+
             w0, n_iter_i, warm_start_sag = sag_solver(
                 X, target, sample_weight, loss, alpha,
                 beta, max_iter, tol,
@@ -815,7 +831,8 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10,
                           verbose=0, solver='lbfgs', penalty='l2',
                           dual=False, intercept_scaling=1.,
                           multi_class='warn', random_state=None,
-                          max_squared_sum=None, sample_weight=None):
+                          max_squared_sum=None, sample_weight=None,
+                          l1_ratio=None):
     """Computes scores across logistic_regression_path
 
     Parameters
@@ -877,9 +894,10 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10,
     solver : {'lbfgs', 'newton-cg', 'liblinear', 'sag', 'saga'}
         Decides which solver to use.
 
-    penalty : str, 'l1' or 'l2'
+    penalty : str, 'l1', 'l2', or 'elastic-net'
         Used to specify the norm used in the penalization. The 'newton-cg',
-        'sag' and 'lbfgs' solvers support only l2 penalties.
+        'sag' and 'lbfgs' solvers support only l2 penalties. 'elastic-net' is
+        only supported by the 'saga' solver.
 
     dual : bool
         Dual or primal formulation. Dual formulation is only implemented for
@@ -921,6 +939,13 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10,
         Array of weights that are assigned to individual samples.
         If not provided, then each sample is given unit weight.
 
+    l1_ratio : float or None, optional (default=None)
+        The elastic net mixing parameter, with ``0 <= l1_ratio <= 1``. Only
+        used if ``penalty='elastic-net'``. Setting ``l1_ratio=0`` is equivalent
+        to using ``penalty='l2'``, while setting ``l1_ratio=1`` is equivalent
+        to using ``penalty='l1'``. For ``0 < l1_ratio <1``, the penalty is a
+        combination of L1 and L2.
+
     Returns
     -------
     coefs : ndarray, shape (n_cs, n_features) or (n_cs, n_features + 1)
@@ -949,13 +974,13 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10,
         sample_weight = sample_weight[train]
 
     coefs, Cs, n_iter = logistic_regression_path(
-        X_train, y_train, Cs=Cs, fit_intercept=fit_intercept,
-        solver=solver, max_iter=max_iter, class_weight=class_weight,
-        pos_class=pos_class, multi_class=multi_class,
-        tol=tol, verbose=verbose, dual=dual, penalty=penalty,
-        intercept_scaling=intercept_scaling, random_state=random_state,
-        check_input=False, max_squared_sum=max_squared_sum,
-        sample_weight=sample_weight)
+        X_train, y_train, Cs=Cs, l1_ratio=l1_ratio,
+        fit_intercept=fit_intercept, solver=solver, max_iter=max_iter,
+        class_weight=class_weight, pos_class=pos_class,
+        multi_class=multi_class, tol=tol, verbose=verbose, dual=dual,
+        penalty=penalty, intercept_scaling=intercept_scaling,
+        random_state=random_state, check_input=False,
+        max_squared_sum=max_squared_sum, sample_weight=sample_weight)
 
     log_reg = LogisticRegression(solver=solver, multi_class=multi_class)
 
@@ -991,6 +1016,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10,
             scores.append(log_reg.score(X_test, y_test))
         else:
             scores.append(scoring(log_reg, X_test, y_test))
+
     return coefs, Cs, np.array(scores), n_iter
 
 
@@ -1002,25 +1028,27 @@ class LogisticRegression(BaseEstimator, LinearClassifierMixin,
     scheme if the 'multi_class' option is set to 'ovr', and uses the cross-
     entropy loss if the 'multi_class' option is set to 'multinomial'.
     (Currently the 'multinomial' option is supported only by the 'lbfgs',
-    'sag' and 'newton-cg' solvers.)
+    'sag', 'saga' and 'newton-cg' solvers.)
 
     This class implements regularized logistic regression using the
-    'liblinear' library, 'newton-cg', 'sag' and 'lbfgs' solvers. It can handle
-    both dense and sparse input. Use C-ordered arrays or CSR matrices
+    'liblinear' library, 'newton-cg', 'sag', 'saga' and 'lbfgs' solvers. It can
+    handle both dense and sparse input. Use C-ordered arrays or CSR matrices
     containing 64-bit floats for optimal performance; any other input format
     will be converted (and copied).
 
     The 'newton-cg', 'sag', and 'lbfgs' solvers support only L2 regularization
     with primal formulation. The 'liblinear' solver supports both L1 and L2
-    regularization, with a dual formulation only for the L2 penalty.
+    regularization, with a dual formulation only for the L2 penalty. The
+    elastic net regularization is only supported by the 'saga' solver.
 
     Read more in the :ref:`User Guide <logistic_regression>`.
 
     Parameters
     ----------
-    penalty : str, 'l1' or 'l2', default: 'l2'
+    penalty : str, 'l1', 'l2', or 'elastic-net', default: 'l2'
         Used to specify the norm used in the penalization. The 'newton-cg',
-        'sag' and 'lbfgs' solvers support only l2 penalties.
+        'sag' and 'lbfgs' solvers support only l2 penalties. 'elastic-net' is
+        only supported by the 'saga' solver.
 
         .. versionadded:: 0.19
            l1 penalty with SAGA solver (allowing 'multinomial' + L1)
@@ -1138,6 +1166,13 @@ class LogisticRegression(BaseEstimator, LinearClassifierMixin,
         context. ``-1`` means using all processors.
         See :term:`Glossary <n_jobs>` for more details.
 
+    l1_ratio : float or None, optional (default=None)
+        The elastic net mixing parameter, with ``0 <= l1_ratio <= 1``. Only
+        used if ``penalty='elastic-net'`. Setting ``l1_ratio=0`` is equivalent
+        to using ``penalty='l2'``, while setting ``l1_ratio=1`` is equivalent
+        to using ``penalty='l1'``. For ``0 < l1_ratio <1``, the penalty is a
+        combination of L1 and L2.
+
     Attributes
     ----------
 
@@ -1223,7 +1258,8 @@ class LogisticRegression(BaseEstimator, LinearClassifierMixin,
     def __init__(self, penalty='l2', dual=False, tol=1e-4, C=1.0,
                  fit_intercept=True, intercept_scaling=1, class_weight=None,
                  random_state=None, solver='warn', max_iter=100,
-                 multi_class='warn', verbose=0, warm_start=False, n_jobs=None):
+                 multi_class='warn', verbose=0, warm_start=False, n_jobs=None,
+                 l1_ratio=None):
 
         self.penalty = penalty
         self.dual = dual
@@ -1239,6 +1275,7 @@ def __init__(self, penalty='l2', dual=False, tol=1e-4, C=1.0,
         self.verbose = verbose
         self.warm_start = warm_start
         self.n_jobs = n_jobs
+        self.l1_ratio = l1_ratio
 
     def fit(self, X, y, sample_weight=None):
         """Fit the model according to the given training data.
@@ -1263,9 +1300,21 @@ def fit(self, X, y, sample_weight=None):
         -------
         self : object
         """
+        solver = _check_solver(self.solver, self.penalty, self.dual)
+
         if not isinstance(self.C, numbers.Number) or self.C < 0:
             raise ValueError("Penalty term must be positive; got (C=%r)"
                              % self.C)
+        if self.penalty == 'elastic-net':
+            if (not isinstance(self.l1_ratio, numbers.Number) or
+                    self.l1_ratio < 0 or self.l1_ratio > 1):
+                        raise ValueError("l1_ratio must be between 0 and 1;"
+                                         " got (l1_ratio=%r)" % self.l1_ratio)
+        elif self.l1_ratio is not None:
+            warnings.warn("l1_ratio parameter is only used when penalty is "
+                          "'elastic-net'. Got "
+                          "(penalty={})".format(self.penalty))
+            self.l1_ratio = None
         if not isinstance(self.max_iter, numbers.Number) or self.max_iter < 0:
             raise ValueError("Maximum number of iteration must be positive;"
                              " got (max_iter=%r)" % self.max_iter)
@@ -1273,8 +1322,6 @@ def fit(self, X, y, sample_weight=None):
             raise ValueError("Tolerance for stopping criteria must be "
                              "positive; got (tol=%r)" % self.tol)
 
-        solver = _check_solver(self.solver, self.penalty, self.dual)
-
         if solver in ['newton-cg']:
             _dtype = [np.float64, np.float32]
         else:
@@ -1348,13 +1395,12 @@ def fit(self, X, y, sample_weight=None):
         fold_coefs_ = Parallel(n_jobs=self.n_jobs, verbose=self.verbose,
                                prefer=prefer)(
             path_func(X, y, pos_class=class_, Cs=[self.C],
-                      fit_intercept=self.fit_intercept, tol=self.tol,
-                      verbose=self.verbose, solver=solver,
+                      l1_ratio=self.l1_ratio, fit_intercept=self.fit_intercept,
+                      tol=self.tol, verbose=self.verbose, solver=solver,
                       multi_class=multi_class, max_iter=self.max_iter,
                       class_weight=self.class_weight, check_input=False,
                       random_state=self.random_state, coef=warm_start_coef_,
-                      penalty=self.penalty,
-                      max_squared_sum=max_squared_sum,
+                      penalty=self.penalty, max_squared_sum=max_squared_sum,
                       sample_weight=sample_weight)
             for class_, warm_start_coef_ in zip(classes_, warm_start_coef))
 
@@ -1442,6 +1488,7 @@ class LogisticRegressionCV(LogisticRegression, BaseEstimator,
     of lbfgs optimizer. The newton-cg, sag and lbfgs solvers support only L2
     regularization with primal formulation. The liblinear solver supports both
     L1 and L2 regularization, with a dual formulation only for the L2 penalty.
+    Elastic net penalyy is only supported by the saga solver.
 
     For the grid of Cs values (that are set by default to be ten values in
     a logarithmic scale between 1e-4 and 1e4), the best hyperparameter is
@@ -1451,10 +1498,6 @@ class LogisticRegressionCV(LogisticRegression, BaseEstimator,
     present fit to be the coefficients got after convergence in the previous
     fit, so it is supposed to be faster for high-dimensional dense data.
 
-    For a multiclass problem, the hyperparameters for each class are computed
-    using the best scores got by doing a one-vs-rest in parallel across all
-    folds and classes. Hence this is not the true multinomial loss.
-
     Read more in the :ref:`User Guide <logistic_regression>`.
 
     Parameters
@@ -1485,9 +1528,10 @@ class LogisticRegressionCV(LogisticRegression, BaseEstimator,
         l2 penalty with liblinear solver. Prefer dual=False when
         n_samples > n_features.
 
-    penalty : str, 'l1' or 'l2'
+    penalty : str, 'l1', 'l2', or 'elastic-net', default: 'l2'
         Used to specify the norm used in the penalization. The 'newton-cg',
-        'sag' and 'lbfgs' solvers support only l2 penalties.
+        'sag' and 'lbfgs' solvers support only l2 penalties. 'elastic-net' is
+        only supported by the 'saga' solver.
 
     scoring : string, callable, or None
         A string (see model evaluation documentation) or
@@ -1589,6 +1633,13 @@ class LogisticRegressionCV(LogisticRegression, BaseEstimator,
         If None, the random number generator is the RandomState instance used
         by `np.random`.
 
+    l1_ratios : list of float or None, optional (default=None)
+        The list of elastic net mixing parameter, with ``0 <= l1_ratio <= 1``.
+        Only used if ``penalty='elastic-net'``. A value of 0 is equivalent to
+        using ``penalty='l2'``, while 1 is equivalent to using
+        ``penalty='l1'``. For ``0 < l1_ratio <1``, the penalty is a combination
+        of L1 and L2.
+
     Attributes
     ----------
     coef_ : array, shape (1, n_features) or (n_classes, n_features)
@@ -1607,6 +1658,9 @@ class LogisticRegressionCV(LogisticRegression, BaseEstimator,
         Array of C i.e. inverse of regularization parameter values used
         for cross-validation.
 
+    l1_ratios_ : array
+        Array of l1_ratios used for cross-validation.
+
     coefs_paths_ : array, shape ``(n_folds, len(Cs_), n_features)`` or \
                    ``(n_folds, len(Cs_), n_features + 1)``
         dict with classes as the keys, and the path of coefficients obtained
@@ -1632,6 +1686,12 @@ class LogisticRegressionCV(LogisticRegression, BaseEstimator,
         C's that correspond to the best scores for each fold.
         `C_` is of shape(n_classes,) when the problem is binary.
 
+    l1_ratio_ : array, shape (n_classes,) or (n_classes - 1,)
+        Array of l1_ratio that maps to the best scores across every class. If
+        refit is set to False, then for each class, the best l1_ratio is the
+        average of the l1_ratio's that correspond to the best scores for each
+        fold.  `l1_ratio_` is of shape(n_classes,) when the problem is binary.
+
     n_iter_ : array, shape (n_classes, n_folds, n_cs) or (1, n_folds, n_cs)
         Actual number of iterations for all classes, folds and Cs.
         In the binary or multinomial cases, the first dimension is equal to 1.
@@ -1659,7 +1719,7 @@ def __init__(self, Cs=10, fit_intercept=True, cv='warn', dual=False,
                  penalty='l2', scoring=None, solver='lbfgs', tol=1e-4,
                  max_iter=100, class_weight=None, n_jobs=None, verbose=0,
                  refit=True, intercept_scaling=1., multi_class='warn',
-                 random_state=None):
+                 random_state=None, l1_ratios=None):
         self.Cs = Cs
         self.fit_intercept = fit_intercept
         self.cv = cv
@@ -1676,6 +1736,7 @@ def __init__(self, Cs=10, fit_intercept=True, cv='warn', dual=False,
         self.intercept_scaling = intercept_scaling
         self.multi_class = multi_class
         self.random_state = random_state
+        self.l1_ratios = l1_ratios
 
     def fit(self, X, y, sample_weight=None):
         """Fit the model according to the given training data.
@@ -1705,6 +1766,21 @@ def fit(self, X, y, sample_weight=None):
         if not isinstance(self.tol, numbers.Number) or self.tol < 0:
             raise ValueError("Tolerance for stopping criteria must be "
                              "positive; got (tol=%r)" % self.tol)
+        if self.penalty == 'elastic-net':
+            if self.l1_ratios is None or len(self.l1_ratios) == 0 or any(
+                    (not isinstance(l1_ratio, numbers.Number) or l1_ratio < 0
+                     or l1_ratio > 1) for l1_ratio in self.l1_ratios):
+                raise ValueError("l1_ratios must be a list of numbers between "
+                                 "0 and 1; got (l1_ratios=%r)" %
+                                 self.l1_ratios)
+            l1_ratios_ = self.l1_ratios
+        else:
+            if self.l1_ratios is not None:
+                warnings.warn("l1_ratios parameter is only used when penalty "
+                              "is 'elastic-net'. Got (penalty={})".format(
+                                  self.penalty))
+
+            l1_ratios_ = [None]
 
         X, y = check_X_y(X, y, accept_sparse='csr', dtype=np.float64,
                          order="C",
@@ -1774,6 +1850,7 @@ def fit(self, X, y, sample_weight=None):
             prefer = 'threads'
         else:
             prefer = 'processes'
+
         fold_coefs_ = Parallel(n_jobs=self.n_jobs, verbose=self.verbose,
                                prefer=prefer)(
             path_func(X, y, train, test, pos_class=label, Cs=self.Cs,
@@ -1785,69 +1862,88 @@ def fit(self, X, y, sample_weight=None):
                       intercept_scaling=self.intercept_scaling,
                       random_state=self.random_state,
                       max_squared_sum=max_squared_sum,
-                      sample_weight=sample_weight
+                      sample_weight=sample_weight,
+                      l1_ratio=l1_ratio
                       )
             for label in iter_encoded_labels
-            for train, test in folds)
-
+            for train, test in folds
+            for l1_ratio in l1_ratios_)
+
+        # _log_reg_scoring_path will output different shapes depending on the
+        # multi_class param, so we need to reshape the outputs accordingly.
+        # After reshaping,
+        # - scores is of shape (n_classes X n_folds X n_features)
+        # - coefs_paths is of shape
+        #  (n_classes X n_folds X n_Cs . n_l1_ratios X n_features)
+        # - n_iter is of shape
+        #  (n_classes X n_folds X n_Cs . n_l1_ratios X n_features) or
+        #  (1 X n_folds X n_Cs . n_l1_ratios X n_features)
+        coefs_paths, Cs, scores, n_iter_ = zip(*fold_coefs_)
+        self.Cs_ = Cs[0]
         if multi_class == 'multinomial':
-            multi_coefs_paths, Cs, multi_scores, n_iter_ = zip(*fold_coefs_)
-            multi_coefs_paths = np.asarray(multi_coefs_paths)
-            multi_scores = np.asarray(multi_scores)
-
-            # This is just to maintain API similarity between the ovr and
-            # multinomial option.
-            # Coefs_paths in now n_folds X len(Cs) X n_classes X n_features
-            # we need it to be n_classes X len(Cs) X n_folds X n_features
-            # to be similar to "ovr".
-            coefs_paths = np.rollaxis(multi_coefs_paths, 2, 0)
-
-            # Multinomial has a true score across all labels. Hence the
-            # shape is n_folds X len(Cs). We need to repeat this score
-            # across all labels for API similarity.
-            scores = np.tile(multi_scores, (n_classes, 1, 1))
-            self.Cs_ = Cs[0]
-            self.n_iter_ = np.reshape(n_iter_, (1, len(folds),
-                                                len(self.Cs_)))
-
+            coefs_paths = np.reshape(
+                coefs_paths,
+                (len(folds),  len(l1_ratios_) * len(self.Cs_), n_classes, -1)
+            )
+            # equiv to coefs_paths = np.moveaxis(coefs_paths, (0, 1, 2, 3),
+            #                                                 (1, 2, 0, 3))
+            coefs_paths = np.swapaxes(coefs_paths, 0, 1)
+            coefs_paths = np.swapaxes(coefs_paths, 0, 2)
+            self.n_iter_ = np.reshape(
+                n_iter_,
+                (1, len(folds), len(self.Cs_) * len(l1_ratios_))
+            )
+            # repeat same scores accross all classes
+            scores = np.tile(scores, (n_classes, 1, 1))
         else:
-            coefs_paths, Cs, scores, n_iter_ = zip(*fold_coefs_)
-            self.Cs_ = Cs[0]
-            coefs_paths = np.reshape(coefs_paths, (n_classes, len(folds),
-                                                   len(self.Cs_), -1))
-            self.n_iter_ = np.reshape(n_iter_, (n_classes, len(folds),
-                                                len(self.Cs_)))
-
-        self.coefs_paths_ = dict(zip(classes, coefs_paths))
+            coefs_paths = np.reshape(
+                coefs_paths,
+                (n_classes, len(folds), len(self.Cs_) * len(l1_ratios_),
+                 -1)
+            )
+            self.n_iter_ = np.reshape(
+                n_iter_,
+                (n_classes, len(folds), len(self.Cs_) * len(l1_ratios_))
+            )
         scores = np.reshape(scores, (n_classes, len(folds), -1))
         self.scores_ = dict(zip(classes, scores))
+        self.coefs_paths_ = dict(zip(classes, coefs_paths))
 
         self.C_ = list()
+        self.l1_ratio_ = list()
         self.coef_ = np.empty((n_classes, X.shape[1]))
         self.intercept_ = np.zeros(n_classes)
-
-        # hack to iterate only once for multinomial case.
-        if multi_class == 'multinomial':
-            scores = multi_scores
-            coefs_paths = multi_coefs_paths
-
         for index, (cls, encoded_label) in enumerate(
                 zip(iter_classes, iter_encoded_labels)):
 
             if multi_class == 'ovr':
-                # The scores_ / coefs_paths_ dict have unencoded class
-                # labels as their keys
                 scores = self.scores_[cls]
                 coefs_paths = self.coefs_paths_[cls]
+            else:
+                # For multinomial, all scores are the same accross classes
+                scores = scores[0]
+                # coefs_paths will keep its original shape because
+                # logistic_regression_path expects it this way
 
             if self.refit:
+                # best_index is between 0 and (n_Cs . n_l1_ratios - 1)
+                # for example, with n_cs=2 and n_l1_ratios=3
+                # the layout of scores is
+                # [c1, c2, c1, c2, c1, c2]
+                #   l1_1 ,  l1_2 ,  l1_3
                 best_index = scores.sum(axis=0).argmax()
 
-                C_ = self.Cs_[best_index]
+                best_index_C = best_index % len(self.Cs_)
+                C_ = self.Cs_[best_index_C]
                 self.C_.append(C_)
+
+                best_index_l1 = best_index // len(self.Cs_)
+                l1_ratio_ = l1_ratios_[best_index_l1]
+                self.l1_ratio_.append(l1_ratio_)
+
                 if multi_class == 'multinomial':
-                    coef_init = np.mean(coefs_paths[:, best_index, :, :],
-                                        axis=0)
+                    coef_init = np.mean(coefs_paths[:, :, best_index, :],
+                                        axis=1)
                 else:
                     coef_init = np.mean(coefs_paths[:, best_index, :], axis=0)
 
@@ -1863,19 +1959,30 @@ def fit(self, X, y, sample_weight=None):
                     verbose=max(0, self.verbose - 1),
                     random_state=self.random_state,
                     check_input=False, max_squared_sum=max_squared_sum,
-                    sample_weight=sample_weight)
+                    sample_weight=sample_weight,
+                    l1_ratio=l1_ratio_)
                 w = w[0]
 
             else:
                 # Take the best scores across every fold and the average of all
                 # coefficients corresponding to the best scores.
                 best_indices = np.argmax(scores, axis=1)
-                w = np.mean([coefs_paths[i][best_indices[i]]
-                             for i in range(len(folds))], axis=0)
-                self.C_.append(np.mean(self.Cs_[best_indices]))
+                if self.multi_class == 'ovr':
+                    w = np.mean([coefs_paths[i, best_indices[i], :]
+                                 for i in range(len(folds))], axis=0)
+                else:
+                    w = np.mean([coefs_paths[:, i, best_indices[i], :]
+                                 for i in range(len(folds))], axis=0)
+
+                best_indices_C = best_indices % len(self.Cs_)
+                self.C_.append(np.mean(self.Cs_[best_indices_C]))
+
+                best_indices_l1 = best_indices // len(self.Cs_)
+                self.l1_ratio_.append(np.mean(l1_ratios_[best_indices_l1]))
 
             if multi_class == 'multinomial':
                 self.C_ = np.tile(self.C_, n_classes)
+                self.l1_ratio_ = np.tile(self.l1_ratio_, n_classes)
                 self.coef_ = w[:, :X.shape[1]]
                 if self.fit_intercept:
                     self.intercept_ = w[:, -1]
@@ -1885,6 +1992,9 @@ def fit(self, X, y, sample_weight=None):
                     self.intercept_[index] = w[-1]
 
         self.C_ = np.asarray(self.C_)
+        self.l1_ratio_ = np.asarray(self.l1_ratio_)
+        self.l1_ratios_ = np.asarray(l1_ratios_)
+
         return self
 
     def score(self, X, y, sample_weight=None):
diff --git a/sklearn/linear_model/sag.py b/sklearn/linear_model/sag.py
index 3e8861f26de83..42ae40f1a23e5 100644
--- a/sklearn/linear_model/sag.py
+++ b/sklearn/linear_model/sag.py
@@ -217,7 +217,7 @@ def sag_solver(X, y, sample_weight=None, loss='log', alpha=1., beta=0.,
     >>> clf.fit(X, y)
     ... #doctest: +NORMALIZE_WHITESPACE
     LogisticRegression(C=1.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='sag', tol=0.0001, verbose=0,
         warm_start=False)
diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py
index 4195405b86403..484e6b064ef8b 100644
--- a/sklearn/linear_model/tests/test_logistic.py
+++ b/sklearn/linear_model/tests/test_logistic.py
@@ -11,6 +11,8 @@
 from sklearn.metrics import log_loss
 from sklearn.metrics.scorer import get_scorer
 from sklearn.model_selection import StratifiedKFold
+from sklearn.model_selection import GridSearchCV
+from sklearn.model_selection import train_test_split
 from sklearn.preprocessing import LabelEncoder
 from sklearn.utils import compute_class_weight, _IS_32BIT
 from sklearn.utils.testing import assert_almost_equal
@@ -78,6 +80,7 @@ def test_predict_2_classes():
                                          random_state=0), X_sp, Y1)
 
 
+@pytest.mark.filterwarnings('ignore: Default solver will be changed')  # 0.22
 @pytest.mark.filterwarnings('ignore: Default multi_class will')  # 0.22
 def test_error():
     # Test for appropriate exception on errors
@@ -205,6 +208,7 @@ def test_multinomial_validation(solver):
     assert_raises(ValueError, lr.fit, [[0, 1], [1, 0]], [0, 1])
 
 
+@pytest.mark.filterwarnings('ignore: Default multi_class will')  # 0.22
 @pytest.mark.parametrize('LR', [LogisticRegression, LogisticRegressionCV])
 def test_check_solver_option(LR):
     X, y = iris.data, iris.target
@@ -224,7 +228,7 @@ def test_check_solver_option(LR):
     lr = LR(solver='liblinear', multi_class='multinomial')
     assert_raise_message(ValueError, msg, lr.fit, X, y)
 
-    # all solvers except 'liblinear'
+    # all solvers except 'liblinear' and 'saga'
     for solver in ['newton-cg', 'lbfgs', 'sag']:
         msg = ("Solver %s supports only l2 penalties, got l1 penalty." %
                solver)
@@ -236,6 +240,15 @@ def test_check_solver_option(LR):
         lr = LR(solver=solver, dual=True, multi_class='ovr')
         assert_raise_message(ValueError, msg, lr.fit, X, y)
 
+    # only saga supports elastic-net. We only test for liblinear because the
+    # error is raised before for the other solvers (solver %s supports only l2
+    # penalties)
+    for solver in ['liblinear']:
+        msg = ("Only 'saga' solver supports elastic-net penalty, got "
+               "solver={}.".format(solver))
+        lr = LR(solver=solver, penalty='elastic-net')
+        assert_raise_message(ValueError, msg, lr.fit, X, y)
+
 
 @pytest.mark.parametrize('model, params, warn_solver',
                          [(LogisticRegression, {}, True),
@@ -1379,6 +1392,220 @@ def test_warm_start_converge_LR():
     assert_allclose(lr_no_ws_loss, lr_ws_loss, rtol=1e-5)
 
 
+def test_elastic_net_coeffs():
+    # make sure elastic-net penalty gives different coefficients from l1 and l2
+    # with saga solver (l1_ratio different from 0 or 1)
+    X, y = make_classification(random_state=0)
+
+    C = 1 / .5
+    l1_ratio = .5
+    coeffs = list()
+    for penalty in ('elastic-net', 'l1', 'l2'):
+        lr = LogisticRegression(penalty=penalty, C=C, solver='saga',
+                                random_state=0, l1_ratio=l1_ratio)
+        lr.fit(X, y)
+        coeffs.append(lr.coef_)
+
+    elastic_net_coeffs, l1_coeffs, l2_coeffs = coeffs
+    # make sure coeffs differ by at least .1
+    assert not np.allclose(elastic_net_coeffs, l1_coeffs, rtol=0, atol=.1)
+    assert not np.allclose(elastic_net_coeffs, l2_coeffs, rtol=0, atol=.1)
+    assert not np.allclose(l2_coeffs, l1_coeffs, rtol=0, atol=.1)
+
+
+@pytest.mark.parametrize('C', [.001, .1, 1, 10, 100, 1000, 10000000])
+def test_elastic_net_l1_l2_equivalence(C):
+    # Make sure elastic-net is equivalent to l1 when l1_ratio=1 and l2 when
+    # l1_ratio=0.
+    X, y = make_classification(random_state=0)
+
+    lr = LogisticRegression(penalty='elastic-net', C=C, l1_ratio=0,
+                            solver='saga', random_state=0)
+    lr.fit(X, y)
+    elastic_net_l1_zero_coeffs = lr.coef_
+
+    lr = LogisticRegression(penalty='l2', C=C, solver='saga', random_state=0)
+    lr.fit(X, y)
+    l2_coeffs = lr.coef_
+
+    assert_array_almost_equal(elastic_net_l1_zero_coeffs, l2_coeffs)
+
+    lr = LogisticRegression(penalty='elastic-net', C=C, l1_ratio=1,
+                            solver='saga', random_state=0)
+    lr.fit(X, y)
+    elastic_net_l1_one_coeffs = lr.coef_
+
+    lr = LogisticRegression(penalty='l1', C=C, solver='saga', random_state=0)
+    lr.fit(X, y)
+    l1_coeffs = lr.coef_
+
+    assert_array_almost_equal(elastic_net_l1_one_coeffs, l1_coeffs)
+
+
+@pytest.mark.parametrize('C', [.001, 1, 100, 10000000])
+def test_elastic_net_vs_l1_l2(C):
+    # Make sure that elastic net with grid search on l1_ratio gives same or
+    # better results than just l1 or just l2.
+
+    X, y = make_classification(500, random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
+
+    cv = StratifiedKFold(5, random_state=0)
+    param_grid = {'l1_ratio': np.linspace(0, 1, 10)}
+
+    enet_clf = LogisticRegression(penalty='elastic-net', C=C, solver='saga',
+                                  random_state=0)
+    gs = GridSearchCV(enet_clf, param_grid, cv=cv, iid=False, refit=True)
+
+    l1_clf = LogisticRegression(penalty='l1', C=C, solver='saga',
+                                random_state=0)
+    l2_clf = LogisticRegression(penalty='l2', C=C, solver='saga',
+                                random_state=0)
+
+    for clf in (gs, l1_clf, l2_clf):
+        clf.fit(X_train, y_train)
+
+    assert gs.score(X_test, y_test) >= l1_clf.score(X_test, y_test)
+    assert gs.score(X_test, y_test) >= l2_clf.score(X_test, y_test)
+
+
+def test_elastic_net_CV_grid_search():
+    # make sure ElasticNetCV gives same best params as GridSearchCV
+
+    X, y = make_classification(random_state=0)
+    cv = StratifiedKFold(5, random_state=0)
+
+    l1_ratios = np.linspace(0, 1, 5)
+    Cs = np.logspace(-4, 4, 5)
+
+    lrcv = LogisticRegressionCV(penalty='elastic-net', Cs=Cs, solver='saga',
+                                cv=cv, l1_ratios=l1_ratios, random_state=0)
+    lrcv.fit(X, y)
+
+    param_grid = {'C': Cs, 'l1_ratio': l1_ratios}
+    lr = LogisticRegression(penalty='elastic-net', solver='saga',
+                            random_state=0)
+    gs = GridSearchCV(lr, param_grid, cv=cv)
+    gs.fit(X, y)
+
+    assert gs.best_params_['l1_ratio'] == lrcv.l1_ratio_[0]
+    assert gs.best_params_['C'] == lrcv.C_[0]
+
+
+def test_elastic_net_CV_multiclass_multinomial():
+    # make sure ElasticNetCV gives same best params as GridSearchCV with
+    # multinomial=True.
+
+    X, y = make_classification(n_samples=200, n_classes=3, n_informative=3,
+                               random_state=0)
+    cv = StratifiedKFold(5, random_state=0)
+
+    l1_ratios = np.linspace(0, 1, 5)
+    Cs = np.logspace(-4, 4, 5)
+
+    lrcv = LogisticRegressionCV(penalty='elastic-net', Cs=Cs, solver='saga',
+                                cv=cv, l1_ratios=l1_ratios, random_state=0,
+                                multi_class='multinomial')
+    lrcv.fit(X, y)
+
+    param_grid = {'C': Cs, 'l1_ratio': l1_ratios}
+    lr = LogisticRegression(penalty='elastic-net', solver='saga',
+                            random_state=0, multi_class='multinomial')
+    gs = GridSearchCV(lr, param_grid, cv=cv, iid=False)
+    gs.fit(X, y)
+
+    assert gs.best_params_['l1_ratio'] == lrcv.l1_ratio_[0]
+    assert gs.best_params_['C'] == lrcv.C_[0]
+
+
+def test_elastic_net_CV_multiclass_ovr():
+    # make sure ElasticNetCV gives similar predictions as GridSearchCV in a
+    # multiclass setting with ovr. We can't compare best_params like in the
+    # previous test because LogisticRegressionCV with multi_class='ovr' will
+    # have one C and one l1_param for each class, while LogisticRegression will
+    # share the parameters over the *n_classes* classifiers.
+
+    X, y = make_classification(n_samples=200, n_classes=3, n_informative=3,
+                               random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
+    cv = StratifiedKFold(5, random_state=0)
+
+    l1_ratios = np.linspace(0, 1, 5)
+    Cs = np.logspace(-4, 4, 5)
+
+    lrcv = LogisticRegressionCV(penalty='elastic-net', Cs=Cs, solver='saga',
+                                cv=cv, l1_ratios=l1_ratios, random_state=0,
+                                multi_class='ovr')
+    lrcv.fit(X_train, y_train)
+
+    param_grid = {'C': Cs, 'l1_ratio': l1_ratios}
+    lr = LogisticRegression(penalty='elastic-net', solver='saga',
+                            random_state=0, multi_class='ovr')
+    gs = GridSearchCV(lr, param_grid, cv=cv, iid=False)
+    gs.fit(X_train, y_train)
+
+    # Check that predictions are 80% the same
+    assert (lrcv.predict(X_train) == gs.predict(X_train)).mean() >= .8
+    assert (lrcv.predict(X_test) == gs.predict(X_test)).mean() >= .8
+
+
+@pytest.mark.parametrize('multi_class', ('ovr', 'multinomial'))
+def test_ElasticNetCV_no_refit(multi_class):
+    # Test ElasticNetCV when refit is False
+
+    n_classes = 3
+    n_features = 20
+    X, y = make_classification(n_samples=200, n_classes=n_classes,
+                               n_informative=n_classes, n_features=n_features,
+                               random_state=0)
+    cv = StratifiedKFold(5, random_state=0)
+
+    Cs = np.logspace(-4, 4, 3)
+    l1_ratios = np.linspace(0, 1, 2)
+
+    lrcv = LogisticRegressionCV(penalty='elastic-net', Cs=Cs, solver='saga',
+                                cv=cv, l1_ratios=l1_ratios, random_state=0,
+                                multi_class=multi_class, refit=False)
+    lrcv.fit(X, y)
+    assert lrcv.C_.shape == (n_classes,)
+    assert lrcv.l1_ratio_.shape == (n_classes,)
+    assert lrcv.coef_.shape == (n_classes, n_features)
+
+
+@pytest.mark.parametrize('l1_ratio', (-1, 2, None, 'something_wrong'))
+def test_l1_ratio_param(l1_ratio):
+
+    msg = "l1_ratio must be between 0 and 1; got (l1_ratio=%r)" % l1_ratio
+    assert_raise_message(ValueError, msg,
+                         LogisticRegression(penalty='elastic-net',
+                                            solver='saga',
+                                            l1_ratio=l1_ratio).fit, X, Y1)
+    if l1_ratio is not None:
+        msg = ("l1_ratio parameter is only used when penalty is 'elastic-net'."
+               " Got (penalty=l1)")
+        assert_warns_message(UserWarning, msg,
+                             LogisticRegression(penalty='l1', solver='saga',
+                                                l1_ratio=l1_ratio).fit, X, Y1)
+
+
+@pytest.mark.parametrize('l1_ratios', ([], [.5, 2], None, 'something_wrong'))
+def test_l1_ratios_param(l1_ratios):
+
+    msg = ("l1_ratios must be a list of numbers between 0 and 1; got "
+           "(l1_ratios=%r)" % l1_ratios)
+    assert_raise_message(ValueError, msg,
+                         LogisticRegressionCV(penalty='elastic-net',
+                                              solver='saga',
+                                              l1_ratios=l1_ratios, cv=2).fit,
+                         X, Y1)
+    if l1_ratios is not None:
+        msg = ("l1_ratios parameter is only used when penalty is "
+               "'elastic-net'. Got (penalty=l1)")
+        function = LogisticRegressionCV(penalty='l1', solver='saga',
+                                        l1_ratios=l1_ratios, cv=2).fit
+        assert_warns_message(UserWarning, msg, function, X, Y1)
+
+
 def test_logistic_regression_path_coefs_multinomial():
     # Make sure that the returned coefs by logistic_regression_path when
     # multi_class='multinomial' don't override each other (used to be a


Code review comments:
diff_hunk: @@ -1379,6 +1392,220 @@ def test_warm_start_converge_LR():
     assert_allclose(lr_no_ws_loss, lr_ws_loss, rtol=1e-5)
 
 
+def test_elastic_net_coeffs():
+    # make sure elastic-net penalty gives different coefficients from l1 and l2
+    # with saga solver (l1_ratio different from 0 or 1)
+    X, y = make_classification(random_state=0)
+
+    C = 1 / .5
+    l1_ratio = .5
+    coeffs = list()
+    for penalty in ('elastic-net', 'l1', 'l2'):
+        lr = LogisticRegression(penalty=penalty, C=C, solver='saga',
+                                random_state=0, l1_ratio=l1_ratio)
+        lr.fit(X, y)
+        coeffs.append(lr.coef_)
+
+    elastic_net_coeffs, l1_coeffs, l2_coeffs = coeffs
+    # make sure coeffs differ by at least .1
+    assert not np.allclose(elastic_net_coeffs, l1_coeffs, rtol=0, atol=.1)
+    assert not np.allclose(elastic_net_coeffs, l2_coeffs, rtol=0, atol=.1)
+    assert not np.allclose(l2_coeffs, l1_coeffs, rtol=0, atol=.1)
+
+
+@pytest.mark.parametrize('C', [.001, .1, 1, 10, 100, 1000, 10000000])
+def test_elastic_net_l1_l2_equivalence(C):
+    # Make sure elastic-net is equivalent to l1 when l1_ratio=1 and l2 when
+    # l1_ratio=0.
+    X, y = make_classification(random_state=0)
+
+    lr = LogisticRegression(penalty='elastic-net', C=C, l1_ratio=0,
+                            solver='saga', random_state=0)
+    lr.fit(X, y)
+    elastic_net_l1_zero_coeffs = lr.coef_
+
+    lr = LogisticRegression(penalty='l2', C=C, solver='saga', random_state=0)
+    lr.fit(X, y)
+    l2_coeffs = lr.coef_
+
+    assert_array_almost_equal(elastic_net_l1_zero_coeffs, l2_coeffs)
+
+    lr = LogisticRegression(penalty='elastic-net', C=C, l1_ratio=1,
+                            solver='saga', random_state=0)
+    lr.fit(X, y)
+    elastic_net_l1_one_coeffs = lr.coef_
+
+    lr = LogisticRegression(penalty='l1', C=C, solver='saga', random_state=0)
+    lr.fit(X, y)
+    l1_coeffs = lr.coef_
+
+    assert_array_almost_equal(elastic_net_l1_one_coeffs, l1_coeffs)
+
+
+@pytest.mark.parametrize('C', [.001, 1, 100, 10000000])
+def test_elastic_net_vs_l1_l2(C):
+    # Make sure that elastic net with grid search on l1_ratio gives same or
+    # better results than just l1 or just l2.
+
+    X, y = make_classification(500, random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
+
+    cv = StratifiedKFold(5, random_state=0)
+    param_grid = {'l1_ratio': np.linspace(0, 1, 10)}
+
+    enet_clf = LogisticRegression(penalty='elastic-net', C=C, solver='saga',
+                                  random_state=0)
+    gs = GridSearchCV(enet_clf, param_grid, cv=cv, iid=False, refit=True)
+
+    l1_clf = LogisticRegression(penalty='l1', C=C, solver='saga',
+                                random_state=0)
+    l2_clf = LogisticRegression(penalty='l2', C=C, solver='saga',
+                                random_state=0)
+
+    for clf in (gs, l1_clf, l2_clf):
+        clf.fit(X_train, y_train)
+
+    assert gs.score(X_test, y_test) >= l1_clf.score(X_test, y_test)
+    assert gs.score(X_test, y_test) >= l2_clf.score(X_test, y_test)
+
+
+def test_elastic_net_CV_grid_search():
+    # make sure ElasticNetCV gives same best params as GridSearchCV
path: sklearn/linear_model/tests/test_logistic.py
text: @user1:
LogisticRegressionCV with elastic net? There's also ElasticNetCV which is regression

diff_hunk: @@ -1607,6 +1658,9 @@ class LogisticRegressionCV(LogisticRegression, BaseEstimator,
         Array of C i.e. inverse of regularization parameter values used
         for cross-validation.
 
+    l1_ratios_ : array
+        Array of l1_ratios used for cross-validation.
+
     coefs_paths_ : array, shape ``(n_folds, len(Cs_), n_features)`` or \
path: sklearn/linear_model/logistic.py
text: @user1:
Shouldn't there be one for each `l1_ratio`? Surely l1_ratios should become another dimension. See also shape of `ElasticNetCV.mse_paths_` (which is quite different... and I'm not really sure why we're storing all the coefs; it seems expensive)

@user1:
Similarly, the shape mentioned in `scores_` below shoulld be updated.

@author:
I don't understand what you mean. `l1_ratios_` is just the 'processed' version of the `l1_ratios` parameter, so it's just an array of floats. Just like the `Cs_` attribute right above.

@user1:
But there should be a distinct vector of coefficients for each C and l1_ratio so the shapes documented here should be updated

@author:
Yes absolutely, sorry I totally misread your comment.
The new shape is ``(n_folds, len(Cs_) x len(l1_ratios), n_features)``, I assume this is more correct as adding a dimension could break users code?

diff_hunk: @@ -1379,6 +1392,220 @@ def test_warm_start_converge_LR():
     assert_allclose(lr_no_ws_loss, lr_ws_loss, rtol=1e-5)
 
 
+def test_elastic_net_coeffs():
+    # make sure elastic-net penalty gives different coefficients from l1 and l2
+    # with saga solver (l1_ratio different from 0 or 1)
+    X, y = make_classification(random_state=0)
+
+    C = 1 / .5
+    l1_ratio = .5
+    coeffs = list()
+    for penalty in ('elastic-net', 'l1', 'l2'):
+        lr = LogisticRegression(penalty=penalty, C=C, solver='saga',
+                                random_state=0, l1_ratio=l1_ratio)
+        lr.fit(X, y)
+        coeffs.append(lr.coef_)
+
+    elastic_net_coeffs, l1_coeffs, l2_coeffs = coeffs
+    # make sure coeffs differ by at least .1
+    assert not np.allclose(elastic_net_coeffs, l1_coeffs, rtol=0, atol=.1)
+    assert not np.allclose(elastic_net_coeffs, l2_coeffs, rtol=0, atol=.1)
+    assert not np.allclose(l2_coeffs, l1_coeffs, rtol=0, atol=.1)
+
+
+@pytest.mark.parametrize('C', [.001, .1, 1, 10, 100, 1000, 10000000])
+def test_elastic_net_l1_l2_equivalence(C):
+    # Make sure elastic-net is equivalent to l1 when l1_ratio=1 and l2 when
+    # l1_ratio=0.
+    X, y = make_classification(random_state=0)
+
+    lr = LogisticRegression(penalty='elastic-net', C=C, l1_ratio=0,
+                            solver='saga', random_state=0)
+    lr.fit(X, y)
+    elastic_net_l1_zero_coeffs = lr.coef_
+
+    lr = LogisticRegression(penalty='l2', C=C, solver='saga', random_state=0)
+    lr.fit(X, y)
+    l2_coeffs = lr.coef_
+
+    assert_array_almost_equal(elastic_net_l1_zero_coeffs, l2_coeffs)
+
+    lr = LogisticRegression(penalty='elastic-net', C=C, l1_ratio=1,
+                            solver='saga', random_state=0)
+    lr.fit(X, y)
+    elastic_net_l1_one_coeffs = lr.coef_
+
+    lr = LogisticRegression(penalty='l1', C=C, solver='saga', random_state=0)
+    lr.fit(X, y)
+    l1_coeffs = lr.coef_
+
+    assert_array_almost_equal(elastic_net_l1_one_coeffs, l1_coeffs)
+
+
+@pytest.mark.parametrize('C', [.001, 1, 100, 10000000])
+def test_elastic_net_vs_l1_l2(C):
+    # Make sure that elastic net with grid search on l1_ratio gives same or
+    # better results than just l1 or just l2.
+
+    X, y = make_classification(500, random_state=0)
+    X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)
+
+    cv = StratifiedKFold(5, random_state=0)
path: sklearn/linear_model/tests/test_logistic.py
text: @user1:
why are you doing this instead of cv=5? They are equivalent.