diff --git a/doc/modules/classes.rst b/doc/modules/classes.rst
index 6cd3128940d..495bcd223f6 100644
--- a/doc/modules/classes.rst
+++ b/doc/modules/classes.rst
@@ -639,6 +639,7 @@ From text
    linear_model.LassoLarsIC
    linear_model.LinearRegression
    linear_model.LogisticRegression
+   linear_model.LogisticRegressionCV
    linear_model.MultiTaskLasso
    linear_model.MultiTaskElasticNet
    linear_model.MultiTaskLassoCV
diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst
index 5a1b8580d87..55a8c5b8e13 100644
--- a/doc/modules/linear_model.rst
+++ b/doc/modules/linear_model.rst
@@ -649,10 +649,10 @@ rather than regression. Logistic regression is also known in the literature as
 logit regression, maximum-entropy classification (MaxEnt)
 or the log-linear classifier. In this model, the probabilities describing the possible outcomes of a single trial are modeled using a `logistic function <http://en.wikipedia.org/wiki/Logistic_function>`_.
 
-The implementation of logistic regression in scikit-learn can be accessed from 
-class :class:`LogisticRegression`. This 
+The implementation of logistic regression in scikit-learn can be accessed from
+class :class:`LogisticRegression`. This
 implementation can fit a multiclass (one-vs-rest) logistic regression with optional
-L2 or L1 regularization. 
+L2 or L1 regularization.
 
 As an optimization problem, binary class L2 penalized logistic regression minimizes
 the following cost function:
@@ -663,7 +663,14 @@ Similarly, L1 regularized logistic regression solves the following optimization
 
 .. math:: \underset{w, c}{min\,} \|w\|_1 + C \sum_{i=1}^n \log(\exp(- y_i (X_i^T w + c)) + 1) .
 
-L1 penalization yields sparse predicting weights.
+The solvers implemented in the class :class:`LogisticRegression`
+are "liblinear" (which is a wrapper around the C++ library,
+LIBLINEAR), "newton-cg" and "lbfgs".
+
+The lbfgs and newton-cg solvers only support L2 penalization and are found
+to converge faster for some high dimensional data. L1 penalization yields
+sparse predicting weights.
+
 For L1 penalization :func:`sklearn.svm.l1_min_c` allows to calculate
 the lower bound for C in order to get a non "null" (all feature weights to
 zero) model.
@@ -685,6 +692,12 @@ which is shipped with scikit-learn.
    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. In
+general the "newton-cg" and "lbfgs" solvers are found to be faster
+due to warm-starting. For the multiclass case, One-vs-All is used
+and an optimal C is obtained for each class.
+
 
 Stochastic Gradient Descent - SGD
 =================================
diff --git a/doc/tutorial/statistical_inference/supervised_learning.rst b/doc/tutorial/statistical_inference/supervised_learning.rst
index 7f54a1e92e9..3475542f718 100644
--- a/doc/tutorial/statistical_inference/supervised_learning.rst
+++ b/doc/tutorial/statistical_inference/supervised_learning.rst
@@ -370,8 +370,8 @@ function or **logistic** function:
     >>> logistic = linear_model.LogisticRegression(C=1e5)
     >>> logistic.fit(iris_X_train, iris_y_train)
     LogisticRegression(C=100000.0, class_weight=None, dual=False,
-              fit_intercept=True, intercept_scaling=1, penalty='l2',
-              random_state=None, tol=0.0001)
+              fit_intercept=True, intercept_scaling=1, max_iter=100,
+              penalty='l2', random_state=None, solver='liblinear', tol=0.0001)
 
 This is known as :class:`LogisticRegression`.
 
diff --git a/doc/whats_new.rst b/doc/whats_new.rst
index b23d961beed..b899171d588 100644
--- a/doc/whats_new.rst
+++ b/doc/whats_new.rst
@@ -21,6 +21,9 @@ New features
    - Add the :func:`metrics.label_ranking_average_precision_score` metrics. By
      `Arnaud Joly`_.
 
+   - Added :class:`linear_model.LogisticRegressionCV`. By
+     `Manoj Kumar`_, `Fabian Pedregosa`_, `Gael Varoquaux`_
+     and `Alexandre Gramfort`_.
 
 Enhancements
 ............
@@ -29,6 +32,9 @@ Enhancements
    - Add support for sample weights in scorer objects.  Metrics with sample
      weight support will automatically benefit from it.
 
+   - Added ``newton-cg`` and `lbfgs` solver support in
+     :class:`linear_model.LogisticRegression`. By `Manoj Kumar`_.
+
 
 Documentation improvements
 ..........................
diff --git a/sklearn/linear_model/__init__.py b/sklearn/linear_model/__init__.py
index 22665207d99..2a7c6b70280 100644
--- a/sklearn/linear_model/__init__.py
+++ b/sklearn/linear_model/__init__.py
@@ -22,7 +22,8 @@ from .sgd_fast import Hinge, Log, ModifiedHuber, SquaredLoss, Huber
 from .stochastic_gradient import SGDClassifier, SGDRegressor
 from .ridge import (Ridge, RidgeCV, RidgeClassifier, RidgeClassifierCV,
                     ridge_regression)
-from .logistic import LogisticRegression
+from .logistic import (LogisticRegression, LogisticRegressionCV,
+                       logistic_regression_path)
 from .omp import (orthogonal_mp, orthogonal_mp_gram, OrthogonalMatchingPursuit,
                   OrthogonalMatchingPursuitCV)
 from .passive_aggressive import PassiveAggressiveClassifier
@@ -48,6 +49,7 @@ __all__ = ['ARDRegression',
            'LinearRegression',
            'Log',
            'LogisticRegression',
+           'LogisticRegressionCV',
            'ModifiedHuber',
            'MultiTaskElasticNet',
            'MultiTaskElasticNetCV',
@@ -71,6 +73,7 @@ __all__ = ['ARDRegression',
            'lars_path',
            'lasso_path',
            'lasso_stability_path',
+           'logistic_regression_path',
            'orthogonal_mp',
            'orthogonal_mp_gram',
            'ridge_regression',
diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py
index 0e0b6aaa999..2162092affa 100644
--- a/sklearn/linear_model/logistic.py
+++ b/sklearn/linear_model/logistic.py
@@ -1,12 +1,565 @@
-# Authors: Fabian Pedregosa
-#          Alexandre Gramfort
-# License: 3-clause BSD
+"""
+Logistic Regression
+"""
+
+# Author: Gael Varoquaux <gael.varoquaux@normalesup.org>
+#         Fabian Pedregosa <f@bianp.net>
+#         Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
+#         Manoj Kumar <manojkumarsivaraj334@gmail.com>
+
+import numbers
+import warnings
 
 import numpy as np
+from scipy import optimize, sparse
 
-from .base import LinearClassifierMixin, SparseCoefMixin
+from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator
 from ..feature_selection.from_model import _LearntSelectorMixin
+from ..preprocessing import LabelEncoder
 from ..svm.base import BaseLibLinear
+from ..utils import check_array, check_consistent_length, compute_class_weight
+from ..utils.extmath import log_logistic, safe_sparse_dot
+from ..utils.optimize import newton_cg
+from ..utils.validation import as_float_array, DataConversionWarning
+from ..utils.fixes import expit
+from ..externals.joblib import Parallel, delayed
+from ..cross_validation import _check_cv
+from ..externals import six
+from ..metrics import SCORERS
+
+
+# .. some helper functions for logistic_regression_path ..
+def _intercept_dot(w, X, y):
+    """Computes y * np.dot(X, w).
+
+    It takes into consideration if the intercept should be fit or not.
+
+    Parameters
+    ----------
+    w : ndarray, shape (n_features,) or (n_features + 1,)
+        Coefficient vector.
+
+    X : {array-like, sparse matrix}, shape (n_samples, n_features)
+        Training data.
+
+    y : ndarray, shape (n_samples,)
+        Array of labels.
+    """
+    c = 0.
+    if w.size == X.shape[1] + 1:
+        c = w[-1]
+        w = w[:-1]
+
+    z = safe_sparse_dot(X, w) + c
+    return w, c, y * z
+
+
+def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None):
+    """Computes the logistic loss and gradient.
+
+    Parameters
+    ----------
+    w : ndarray, shape (n_features,) or (n_features + 1,)
+        Coefficient vector.
+
+    X : {array-like, sparse matrix}, shape (n_samples, n_features)
+        Training data.
+
+    y : ndarray, shape (n_samples,)
+        Array of labels.
+
+    alpha : float
+        Regularization parameter. alpha is equal to 1 / C.
+
+    sample_weight : ndarray, shape (n_samples,) optional
+        Array of weights that are assigned to individual samples.
+        If not provided, then each sample is given unit weight.
+
+    Returns
+    -------
+    out : float
+        Logistic loss.
+
+    grad : ndarray, shape (n_features,) or (n_features + 1,)
+        Logistic gradient.
+    """
+    _, n_features = X.shape
+    grad = np.empty_like(w)
+
+    w, c, yz = _intercept_dot(w, X, y)
+
+    if sample_weight is None:
+        sample_weight = np.ones(y.shape[0])
+
+    # Logistic loss is the negative of the log of the logistic function.
+    out = -np.sum(sample_weight * log_logistic(yz)) + .5 * alpha * np.dot(w, w)
+
+    z = expit(yz)
+    z0 = sample_weight * (z - 1) * y
+
+    grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w
+
+    # Case where we fit the intercept.
+    if grad.shape[0] > n_features:
+        grad[-1] = z0.sum()
+    return out, grad
+
+
+def _logistic_loss(w, X, y, alpha, sample_weight=None):
+    """Computes the logistic loss.
+
+    Parameters
+    ----------
+    w : ndarray, shape (n_features,) or (n_features + 1,)
+        Coefficient vector.
+
+    X : {array-like, sparse matrix}, shape (n_samples, n_features)
+        Training data.
+
+    y : ndarray, shape (n_samples,)
+        Array of labels.
+
+    alpha : float
+        Regularization parameter. alpha is equal to 1 / C.
+
+    sample_weight : ndarray, shape (n_samples,) optional
+        Array of weights that are assigned to individual samples.
+        If not provided, then each sample is given unit weight.
+
+    Returns
+    -------
+    out : float
+        Logistic loss.
+    """
+    w, c, yz = _intercept_dot(w, X, y)
+
+    if sample_weight is None:
+        sample_weight = np.ones(y.shape[0])
+
+    # Logistic loss is the negative of the log of the logistic function.
+    out = -np.sum(sample_weight * log_logistic(yz)) + .5 * alpha * np.dot(w, w)
+    return out
+
+
+def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None):
+    """Computes the logistic loss, gradient and the Hessian.
+
+    Parameters
+    ----------
+    w : ndarray, shape (n_features,) or (n_features + 1,)
+        Coefficient vector.
+
+    X : {array-like, sparse matrix}, shape (n_samples, n_features)
+        Training data.
+
+    y : ndarray, shape (n_samples,)
+        Array of labels.
+
+    alpha : float
+        Regularization parameter. alpha is equal to 1 / C.
+
+    sample_weight : ndarray, shape (n_samples,) optional
+        Array of weights that are assigned to individual samples.
+        If not provided, then each sample is given unit weight.
+
+    Returns
+    -------
+    out : float
+        Logistic loss.
+
+    grad : ndarray, shape (n_features,) or (n_features + 1,)
+        Logistic gradient.
+
+    Hs : callable
+        Function that takes the gradient as a parameter and returns the
+        matrix product of the Hessian and gradient.
+    """
+    n_samples, n_features = X.shape
+    grad = np.empty_like(w)
+    fit_intercept = grad.shape[0] > n_features
+
+    w, c, yz = _intercept_dot(w, X, y)
+
+    if sample_weight is None:
+        sample_weight = np.ones(y.shape[0])
+
+    # Logistic loss is the negative of the log of the logistic function.
+    out = -np.sum(sample_weight * log_logistic(yz)) + .5 * alpha * np.dot(w, w)
+
+    z = expit(yz)
+    z0 = sample_weight * (z - 1) * y
+
+    grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w
+
+    # Case where we fit the intercept.
+    if fit_intercept:
+        grad[-1] = z0.sum()
+
+    # The mat-vec product of the Hessian
+    d = sample_weight * z * (1 - z)
+    if sparse.issparse(X):
+        dX = safe_sparse_dot(sparse.dia_matrix((d, 0),
+                             shape=(n_samples, n_samples)), X)
+    else:
+        # Precompute as much as possible
+        dX = d[:, np.newaxis] * X
+
+    if fit_intercept:
+        # Calculate the double derivative with respect to intercept
+        # In the case of sparse matrices this returns a matrix object.
+        dd_intercept = np.squeeze(np.array(dX.sum(axis=0)))
+
+    def Hs(s):
+        ret = np.empty_like(s)
+        ret[:n_features] = X.T.dot(dX.dot(s[:n_features]))
+        ret[:n_features] += alpha * s[:n_features]
+
+        # For the fit intercept case.
+        if fit_intercept:
+            ret[:n_features] += s[-1] * dd_intercept
+            ret[-1] = dd_intercept.dot(s[:n_features])
+            ret[-1] += d.sum() * s[-1]
+        return ret
+
+    return out, grad, Hs
+
+
+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, copy=True,
+                             class_weight=None, dual=False, penalty='l2',
+                             intercept_scaling=1.):
+    """Compute a Logistic Regression model for a list of regularization
+    parameters.
+
+    This is an implementation that uses the result of the previous model
+    to speed up computations along the set of solutions, making it faster
+    than sequentially calling LogisticRegression for the different parameters.
+
+    Parameters
+    ----------
+    X : array-like or sparse matrix, shape (n_samples, n_features)
+        Input data.
+
+    y : array-like, shape (n_samples,)
+        Input data, target values.
+
+    Cs : int | array-like, shape (n_cs,)
+        List of values for the regularization parameter or integer specifying
+        the number of regularization parameters that should be used. In this
+        case, the parameters will be chosen in a logarithmic scale between
+        1e-4 and 1e4.
+
+    pos_class : int, None
+        The class with respect to which we perform a one-vs-all fit.
+        If None, then it is assumed that the given problem is binary.
+
+    fit_intercept : bool
+        Whether to fit an intercept for the model. In this case the shape of
+        the returned array is (n_cs, n_features + 1).
+
+    max_iter : int
+        Maximum number of iterations for the solver.
+
+    tol : float
+        Stopping criterion. For the newton-cg 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
+        Print convergence message if True.
+
+    solver : {'lbfgs', 'newton-cg', 'liblinear'}
+        Numerical solver to use.
+
+    coef : array-like, shape (n_features,), default None
+        Initialization value for coefficients of logistic regression.
+
+    copy : bool, default True
+        Whether or not to produce a copy of the data. Setting this to
+        True will be useful in cases, when logistic_regression_path
+        is called repeatedly with the same data, as y is modified
+        along the path.
+
+    class_weight : {dict, 'auto'}, optional
+        Over-/undersamples the samples of each class according to the given
+        weights. If not given, all classes are supposed to have weight one.
+        The 'auto' mode selects weights inversely proportional to class
+        frequencies in the training set.
+
+    dual : bool
+        Dual or primal formulation. Dual formulation is only implemented for
+        l2 penalty with liblinear solver. Prefer dual=False when
+        n_samples > n_features.
+
+    penalty : str, 'l1' or 'l2'
+        Used to specify the norm used in the penalization. The newton-cg and
+        lbfgs solvers support only l2 penalties.
+
+    intercept_scaling : float, default 1.
+        This parameter is useful only when the solver 'liblinear' is used
+        and self.fit_intercept is set to True. In this case, x becomes
+        [x, self.intercept_scaling],
+        i.e. a "synthetic" feature with constant value equals to
+        intercept_scaling is appended to the instance vector.
+        The intercept becomes intercept_scaling * synthetic feature weight
+        Note! the synthetic feature weight is subject to l1/l2 regularization
+        as all other features.
+        To lessen the effect of regularization on synthetic feature weight
+        (and therefore on the intercept) intercept_scaling has to be increased.
+
+    Returns
+    -------
+    coefs : ndarray, shape (n_cs, n_features) or (n_cs, n_features + 1)
+        List of coefficients for the Logistic Regression model. If
+        fit_intercept is set to True then the second dimension will be
+        n_features + 1, where the last item represents the intercept.
+
+    Cs : ndarray
+        Grid of Cs used for cross-validation.
+
+    Notes
+    -----
+    You might get slighly different results with the solver liblinear than
+    with the others since this uses LIBLINEAR which penalizes the intercept.
+    """
+    if isinstance(Cs, numbers.Integral):
+        Cs = np.logspace(-4, 4, Cs)
+
+    X = check_array(X, accept_sparse='csc', dtype=np.float64)
+    y = check_array(y, ensure_2d=False, copy=copy)
+    check_consistent_length(X, y)
+    n_classes = np.unique(y)
+
+    if pos_class is None:
+        if (n_classes.size > 2):
+            raise ValueError('To fit OvA, use the pos_class argument')
+        # np.unique(y) gives labels in sorted order.
+        pos_class = n_classes[1]
+
+    # If class_weights is a dict (provided by the user), the weights
+    # are assigned to the original labels. If it is "auto", then
+    # the class_weights are assigned after masking the labels with a OvA.
+    sample_weight = np.ones(X.shape[0])
+    le = LabelEncoder()
+
+    if isinstance(class_weight, dict):
+        if solver == "liblinear":
+            if n_classes.size == 2:
+                # Reconstruct the weights with keys 1 and -1
+                temp = {}
+                temp[1] = class_weight[pos_class]
+                temp[-1] = class_weight[n_classes[0]]
+                class_weight = temp.copy()
+            else:
+                raise ValueError("In LogisticRegressionCV the liblinear solver "
+                                 "cannot handle multiclass with class_weight "
+                                 "of type dict. Use the lbfgs, newton-cg "
+                                 "solvers or set class_weight='auto'")
+        else:
+            class_weight_ = compute_class_weight(class_weight, n_classes, y)
+            sample_weight = class_weight_[le.fit_transform(y)]
+
+    mask = (y == pos_class)
+    y[mask] = 1
+    y[~mask] = -1
+
+    # To take care of object dtypes
+    y = as_float_array(y, copy=False)
+    if class_weight == "auto":
+        class_weight_ = compute_class_weight(class_weight, [-1, 1], y)
+        sample_weight = class_weight_[le.fit_transform(y)]
+
+    if fit_intercept:
+        w0 = np.zeros(X.shape[1] + 1)
+    else:
+        w0 = np.zeros(X.shape[1])
+
+    if coef is not None:
+        # it must work both giving the bias term and not
+        if not coef.size in (X.shape[1], w0.size):
+            raise ValueError('Initialization coef is not of correct shape')
+        w0[:coef.size] = coef
+    coefs = list()
+
+    for C in Cs:
+        if solver == 'lbfgs':
+            func = _logistic_loss_and_grad
+            try:
+                out = optimize.fmin_l_bfgs_b(
+                    func, w0, fprime=None,
+                    args=(X, y, 1. / C, sample_weight),
+                    iprint=(verbose > 0) - 1, pgtol=tol, maxiter=max_iter)
+            except TypeError:
+                # old scipy doesn't have maxiter
+                out = optimize.fmin_l_bfgs_b(
+                    func, w0, fprime=None,
+                    args=(X, y, 1. / C, sample_weight),
+                    iprint=(verbose > 0) - 1, pgtol=tol)
+            w0 = out[0]
+            if out[2]["warnflag"] == 1:
+                warnings.warn("lbfgs failed to converge. Increase the number "
+                              "of iterations.")
+
+        elif solver == 'newton-cg':
+            grad = lambda x, *args: _logistic_loss_and_grad(x, *args)[1]
+            w0 = newton_cg(_logistic_loss_grad_hess, _logistic_loss, grad, w0,
+                           args=(X, y, 1. / C, sample_weight),
+                           maxiter=max_iter, tol=tol)
+        elif solver == 'liblinear':
+            lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=tol,
+                                    class_weight=class_weight, dual=dual,
+                                    penalty=penalty,
+                                    intercept_scaling=intercept_scaling)
+            lr.fit(X, y)
+            if fit_intercept:
+                w0 = np.concatenate([lr.coef_.ravel(), lr.intercept_])
+            else:
+                w0 = lr.coef_.ravel()
+        else:
+            raise ValueError("solver must be one of {'liblinear', 'lbfgs', "
+                             "'newton-cg'}, got '%s' instead" % solver)
+        coefs.append(w0)
+    return coefs, np.array(Cs)
+
+
+# helper function for LogisticCV
+def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10,
+                          scoring=None, fit_intercept=False,
+                          max_iter=100, tol=1e-4, class_weight=None,
+                          verbose=0, solver='lbfgs', penalty='l2',
+                          dual=False, copy=True, intercept_scaling=1.):
+    """Computes scores across logistic_regression_path
+
+    Parameters
+    ----------
+    X : {array-like, sparse matrix}, shape (n_samples, n_features)
+        Training data.
+
+    y : array-like, shape (n_samples,) or (n_samples, n_targets)
+        Target labels.
+
+    train : list of indices
+        The indices of the train set.
+
+    test : list of indices
+        The indices of the test set.
+
+    pos_class : int, None
+        The class with respect to which we perform a one-vs-all fit.
+        If None, then it is assumed that the given problem is binary.
+
+    Cs : list of floats | int
+        Each of the values in Cs describes the inverse of
+        regularization strength. If Cs is as an int, then a grid of Cs
+        values are chosen in a logarithmic scale between 1e-4 and 1e4.
+        If not provided, then a fixed set of values for Cs are used.
+
+    scoring : callable
+        For a list of scoring functions that can be used, look at
+        :mod:`sklearn.metrics`. The default scoring option used is
+        accuracy_score.
+
+    fit_intercept : bool
+        If False, then the bias term is set to zero. Else the last
+        term of each coef_ gives us the intercept.
+
+    max_iter : int
+        Maximum number of iterations for the solver.
+
+    tol : float
+        Tolerance for stopping criteria.
+
+    class_weight : {dict, 'auto'}, optional
+        Over-/undersamples the samples of each class according to the given
+        weights. If not given, all classes are supposed to have weight one.
+        The 'auto' mode selects weights inversely proportional to class
+        frequencies in the training set.
+
+    verbose : int
+        Amount of verbosity.
+
+    solver : {'lbfgs', 'newton-cg', 'liblinear'}
+        Decides which solver to use.
+
+    penalty : str, 'l1' or 'l2'
+        Used to specify the norm used in the penalization. The newton-cg and
+        lbfgs solvers support only l2 penalties.
+
+    dual : bool
+        Dual or primal formulation. Dual formulation is only implemented for
+        l2 penalty with liblinear solver. Prefer dual=False when
+        n_samples > n_features.
+
+    intercept_scaling : float, default 1.
+        This parameter is useful only when the solver 'liblinear' is used
+        and self.fit_intercept is set to True. In this case, x becomes
+        [x, self.intercept_scaling],
+        i.e. a "synthetic" feature with constant value equals to
+        intercept_scaling is appended to the instance vector.
+        The intercept becomes intercept_scaling * synthetic feature weight
+        Note! the synthetic feature weight is subject to l1/l2 regularization
+        as all other features.
+        To lessen the effect of regularization on synthetic feature weight
+        (and therefore on the intercept) intercept_scaling has to be increased.
+
+    Returns
+    -------
+    coefs : ndarray, shape (n_cs, n_features) or (n_cs, n_features + 1)
+        List of coefficients for the Logistic Regression model. If
+        fit_intercept is set to True then the second dimension will be
+        n_features + 1, where the last item represents the intercept.
+
+    Cs : ndarray
+        Grid of Cs used for cross-validation.
+
+    scores : ndarray, shape (n_cs,)
+        Scores obtained for each Cs.
+    """
+
+    log_reg = LogisticRegression(fit_intercept=fit_intercept)
+    log_reg._enc = LabelEncoder()
+    log_reg._enc.fit_transform([-1, 1])
+
+    X_train = X[train]
+    X_test = X[test]
+    y_train = y[train]
+    y_test = y[test]
+
+    if pos_class is not None:
+        mask = (y_test == pos_class)
+        y_test[mask] = 1
+        y_test[~mask] = -1
+
+    # To deal with object dtypes, we need to convert into an array of floats.
+    y_test = as_float_array(y_test, copy=False)
+
+    coefs, Cs = logistic_regression_path(X_train, y_train, Cs=Cs,
+                                         fit_intercept=fit_intercept,
+                                         solver=solver,
+                                         max_iter=max_iter,
+                                         class_weight=class_weight,
+                                         copy=copy, pos_class=pos_class,
+                                         tol=tol, verbose=verbose,
+                                         dual=dual, penalty=penalty,
+                                         intercept_scaling=intercept_scaling)
+
+    scores = list()
+
+    if isinstance(scoring, six.string_types):
+        scoring = SCORERS[scoring]
+    for w in coefs:
+        if fit_intercept:
+            log_reg.coef_ = w[np.newaxis, :-1]
+            log_reg.intercept_ = w[-1]
+        else:
+            log_reg.coef_ = w[np.newaxis, :]
+            log_reg.intercept_ = 0.
+        if scoring is None:
+            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)
 
 
 class LogisticRegression(BaseLibLinear, LinearClassifierMixin,
@@ -16,19 +569,25 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin,
     In the multiclass case, the training algorithm uses a one-vs.-all (OvA)
     scheme, rather than the "true" multinomial LR.
 
-    This class implements L1 and L2 regularized logistic regression using the
-    `liblinear` library. 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).
+    This class implements regularized logistic regression using the
+    `liblinear` library, newton-cg 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 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.
 
     Parameters
     ----------
-    penalty : string, 'l1' or 'l2'
-        Used to specify the norm used in the penalization.
+    penalty : str, 'l1' or 'l2'
+        Used to specify the norm used in the penalization. The newton-cg and
+        lbfgs solvers support only l2 penalties.
 
-    dual : boolean
-        Dual or primal formulation. Dual formulation is only
-        implemented for l2 penalty. Prefer dual=False when
+    dual : bool
+        Dual or primal formulation. Dual formulation is only implemented for
+        l2 penalty with liblinear solver. Prefer dual=False when
         n_samples > n_features.
 
     C : float, optional (default=1.0)
@@ -49,7 +608,7 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin,
         Note! the synthetic feature weight is subject to l1/l2 regularization
         as all other features.
         To lessen the effect of regularization on synthetic feature weight
-        (and therefore on the intercept) intercept_scaling has to be increased
+        (and therefore on the intercept) intercept_scaling has to be increased.
 
     class_weight : {dict, 'auto'}, optional
         Over-/undersamples the samples of each class according to the given
@@ -57,27 +616,34 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin,
         The 'auto' mode selects weights inversely proportional to class
         frequencies in the training set.
 
-    random_state: int seed, RandomState instance, or None (default)
+    max_iter : int
+        Useful only for the newton-cg and lbfgs solvers. Maximum number of
+        iterations taken for the solvers to converge.
+
+    random_state : int seed, RandomState instance, or None (default)
         The seed of the pseudo random number generator to use when
         shuffling the data.
 
-    tol: float, optional
+    solver : {'newton-cg', 'lbfgs', 'liblinear'}
+        Algorithm to use in the optimization problem.
+
+    tol : float, optional
         Tolerance for stopping criteria.
 
     Attributes
     ----------
-    `coef_` : array, shape = [n_classes, n_features]
+    `coef_` : array, shape (n_classes, n_features)
         Coefficient of the features in the decision function.
 
-    `intercept_` : array, shape = [n_classes]
+    `intercept_` : array, shape (n_classes,)
         Intercept (a.k.a. bias) added to the decision function.
         If `fit_intercept` is set to False, the intercept is set to zero.
 
     See also
     --------
-    SGDClassifier: incrementally trained logistic regression (when given
+    SGDClassifier : incrementally trained logistic regression (when given
         the parameter ``loss="log"``).
-    sklearn.svm.LinearSVC: learns SVM models using the same algorithm.
+    sklearn.svm.LinearSVC : learns SVM models using the same algorithm.
 
     Notes
     -----
@@ -99,12 +665,13 @@ class LogisticRegression(BaseLibLinear, 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):
+                 random_state=None, solver='liblinear', max_iter=100):
 
         super(LogisticRegression, self).__init__(
             penalty=penalty, dual=dual, loss='lr', tol=tol, C=C,
             fit_intercept=fit_intercept, intercept_scaling=intercept_scaling,
-            class_weight=class_weight, random_state=random_state)
+            class_weight=class_weight, random_state=random_state,
+            solver=solver, max_iter=max_iter)
 
     def predict_proba(self, X):
         """Probability estimates.
@@ -141,3 +708,289 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin,
             model, where classes are ordered as they are in ``self.classes_``.
         """
         return np.log(self.predict_proba(X))
+
+
+class LogisticRegressionCV(LogisticRegression, BaseEstimator,
+                           LinearClassifierMixin, _LearntSelectorMixin):
+    """Logistic Regression CV (aka logit, MaxEnt) classifier.
+
+    This class implements logistic regression using liblinear, newton-cg or
+    LBFGS optimizer. The newton-cg 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.
+
+    Parameters
+    ----------
+    Cs : list of floats | int
+        Each of the values in Cs describes the inverse of regularization
+        strength. If Cs is as an int, then a grid of Cs values are chosen
+        in a logarithmic scale between 1e-4 and 1e4.
+        Like in support vector machines, smaller values specify stronger
+        regularization.
+
+    fit_intercept : bool, default: True
+        Specifies if a constant (a.k.a. bias or intercept) should be
+        added the decision function.
+
+    class_weight : {dict, 'auto'}, optional
+        Over-/undersamples the samples of each class according to the given
+        weights. If not given, all classes are supposed to have weight one.
+        The 'auto' mode selects weights inversely proportional to class
+        frequencies in the training set.
+
+    cv : integer or cross-validation generator
+        The default cross-validation generator used is Stratified K-Folds.
+        If an integer is provided, then it is the number of folds used.
+        See the module :mod:`sklearn.cross_validation` module for the
+        list of possible cross-validation objects.
+
+    penalty : str, 'l1' or 'l2'
+        Used to specify the norm used in the penalization. The newton-cg and
+        lbfgs solvers support only l2 penalties.
+
+    dual : bool
+        Dual or primal formulation. Dual formulation is only implemented for
+        l2 penalty with liblinear solver. Prefer dual=False when
+        n_samples > n_features.
+
+    scoring : callabale
+        Scoring function to use as cross-validation criteria. For a list of
+        scoring functions that can be used, look at :mod:`sklearn.metrics`.
+        The default scoring option used is accuracy_score.
+
+    solver : {'newton-cg', 'lbfgs', 'liblinear'}
+        Algorithm to use in the optimization problem.
+
+    tol : float, optional
+        Tolerance for stopping criteria.
+
+    max_iter : int, optional
+        Maximum number of iterations of the optimization algorithm.
+
+    class_weight : {dict, 'auto'}, optional
+        Over-/undersamples the samples of each class according to the given
+        weights. If not given, all classes are supposed to have weight one.
+        The 'auto' mode selects weights inversely proportional to class
+        frequencies in the training set.
+
+    n_jobs : int, optional
+        Number of CPU cores used during the cross-validation loop. If given
+        a value of -1, all cores are used.
+
+    verbose : bool | int
+        Amount of verbosity.
+
+    refit : bool
+        If set to True, the scores are averaged across all folds, and the
+        coefs and the C that corresponds to the best score is taken, and a
+        final refit is done using these parameters.
+        Otherwise the coefs, intercepts and C that correspond to the
+        best scores across folds are averaged.
+
+    intercept_scaling : float, default 1.
+        This parameter is useful only when the solver 'liblinear' is used
+        and self.fit_intercept is set to True. In this case, x becomes
+        [x, self.intercept_scaling],
+        i.e. a "synthetic" feature with constant value equals to
+        intercept_scaling is appended to the instance vector.
+        The intercept becomes intercept_scaling * synthetic feature weight
+        Note! the synthetic feature weight is subject to l1/l2 regularization
+        as all other features.
+        To lessen the effect of regularization on synthetic feature weight
+        (and therefore on the intercept) intercept_scaling has to be increased.
+
+    Attributes
+    ----------
+    `coef_` : array, shape (1, n_features) or (n_classes, n_features)
+        Coefficient of the features in the decision function.
+
+        `coef_` is of shape (1, n_features) when the given problem
+        is binary.
+        `coef_` is readonly property derived from `raw_coef_` that
+        follows the internal memory layout of liblinear.
+
+    `intercept_` : array, shape (1,) or (n_classes,)
+        Intercept (a.k.a. bias) added to the decision function.
+        It is available only when parameter intercept is set to True
+        and is of shape(1,) when the problem is binary.
+
+    `Cs_` : array
+        Array of C i.e. inverse of regularization parameter values 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
+        during cross-validating across each fold and then across each Cs
+        after doing an OvA for the corresponding class.
+        Each dict value has shape (n_folds, len(Cs_), n_features) or
+        (n_folds, len(Cs_), n_features + 1) depending on whether the
+        intercept is fit or not.
+
+    `scores_` : dict
+        dict with classes as the keys, and the values as the
+        grid of scores obtained during cross-validating each fold, after doing
+        an OvA for the corresponding class.
+        Each dict value has shape (n_folds, len(Cs))
+
+    `C_` : array, shape (n_classes,) or (n_classes - 1,)
+        Array of C that maps to the best scores across every class. If refit is
+        set to False, then for each class, the best C is the average of the
+        C's that correspond to the best scores for each fold.
+
+    See also
+    --------
+    LogisticRegression
+
+    """
+
+    def __init__(self, Cs=10, fit_intercept=True, cv=None, dual=False,
+                 penalty='l2', scoring=None, solver='lbfgs', tol=1e-4,
+                 max_iter=100, class_weight=None, n_jobs=1, verbose=False,
+                 refit=True, intercept_scaling=1.):
+        self.Cs = Cs
+        self.fit_intercept = fit_intercept
+        self.cv = cv
+        self.dual = dual
+        self.penalty = penalty
+        self.scoring = scoring
+        self.tol = tol
+        self.max_iter = max_iter
+        self.class_weight = class_weight
+        self.n_jobs = n_jobs
+        self.verbose = verbose
+        self.solver = solver
+        self.refit = refit
+        self.intercept_scaling = 1.
+
+    def fit(self, X, y):
+        """Fit the model according to the given training data.
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix}, shape (n_samples, n_features)
+            Training vector, where n_samples in the number of samples and
+            n_features is the number of features.
+
+        y : array-like, shape (n_samples,)
+            Target vector relative to X.
+
+        Returns
+        -------
+        self : object
+            Returns self.
+        """
+        if self.solver != 'liblinear':
+            if self.penalty != 'l2':
+                raise ValueError("newton-cg and lbfgs solvers support only "
+                                 "l2 penalties.")
+            if self.dual:
+                raise ValueError("newton-cg and lbfgs solvers support only "
+                                 "the primal form.")
+
+        X = check_array(X, accept_sparse='csc', dtype=np.float64)
+        y = check_array(y, ensure_2d=False)
+
+        if y.ndim == 2 and y.shape[1] == 1:
+            warnings.warn(
+                "A column-vector y was passed when a 1d array was"
+                " expected. Please change the shape of y to "
+                "(n_samples, ), for example using ravel().",
+                DataConversionWarning
+                )
+            y = np.ravel(y)
+
+        check_consistent_length(X, y)
+
+        # init cross-validation generator
+        cv = _check_cv(self.cv, X, y, classifier=True)
+        folds = list(cv)
+
+        self._enc = LabelEncoder()
+        self._enc.fit(y)
+
+        labels = self.classes_
+        n_classes = len(labels)
+
+        if n_classes < 2:
+            raise ValueError("Number of classes have to be greater than one.")
+
+        if n_classes == 2:
+            # OvA in case of binary problems is as good as fitting
+            # the higher label
+            n_classes = 1
+            labels = labels[1:]
+
+        if self.class_weight and not(isinstance(self.class_weight, dict) or
+                                     self.class_weight == 'auto'):
+            raise ValueError("class_weight provided should be a "
+                             "dict or 'auto'")
+
+        fold_coefs_ = Parallel(n_jobs=self.n_jobs, verbose=self.verbose)(
+            delayed(_log_reg_scoring_path)(X, y, train, test,
+                                           pos_class=label,
+                                           Cs=self.Cs,
+                                           fit_intercept=self.fit_intercept,
+                                           penalty=self.penalty,
+                                           dual=self.dual,
+                                           solver=self.solver,
+                                           max_iter=self.max_iter,
+                                           tol=self.tol,
+                                           class_weight=self.class_weight,
+                                           verbose=max(0, self.verbose - 1),
+                                           scoring=self.scoring,
+                                           intercept_scaling=self.intercept_scaling)
+            for label in labels
+            for train, test in folds
+            )
+        coefs_paths, Cs, scores = zip(*fold_coefs_)
+
+        self.Cs_ = Cs[0]
+        coefs_paths = np.reshape(coefs_paths, (n_classes, len(folds),
+                                 len(self.Cs_), -1))
+        self.coefs_paths_ = dict(zip(labels, coefs_paths))
+        scores = np.reshape(scores, (n_classes, len(folds), -1))
+        self.scores_ = dict(zip(labels, scores))
+
+        self.C_ = list()
+        self.coef_ = list()
+        self.intercept_ = list()
+
+        for label in labels:
+            scores = self.scores_[label]
+            coefs_paths = self.coefs_paths_[label]
+
+            if self.refit:
+                best_index = scores.sum(axis=0).argmax()
+                C_ = self.Cs_[best_index]
+                self.C_.append(C_)
+                coef_init = np.mean(coefs_paths[:, best_index, :], axis=0)
+
+                w, _ = logistic_regression_path(
+                    X, y, pos_class=label, Cs=[C_], solver=self.solver,
+                    fit_intercept=self.fit_intercept, coef=coef_init,
+                    max_iter=self.max_iter, tol=self.tol,
+                    class_weight=self.class_weight,
+                    verbose=max(0, self.verbose - 1))
+                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.fit_intercept:
+                self.coef_.append(w[:-1])
+                self.intercept_.append(w[-1])
+            else:
+                self.coef_.append(w)
+                self.intercept_.append(0.)
+        self.C_ = np.asarray(self.C_)
+        self.coef_ = np.asarray(self.coef_)
+        self.intercept_ = np.asarray(self.intercept_)
+        return self
diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py
index f3d11ed0581..f547e6a9394 100644
--- a/sklearn/linear_model/tests/test_logistic.py
+++ b/sklearn/linear_model/tests/test_logistic.py
@@ -1,6 +1,8 @@
 import numpy as np
 import scipy.sparse as sp
+from scipy import linalg, optimize
 
+from sklearn.utils.testing import assert_almost_equal
 from sklearn.utils.testing import assert_array_equal
 from sklearn.utils.testing import assert_array_almost_equal
 from sklearn.utils.testing import assert_equal
@@ -8,15 +10,21 @@ from sklearn.utils.testing import assert_greater
 from sklearn.utils.testing import assert_raises
 from sklearn.utils.testing import assert_true
 from sklearn.utils.testing import raises
+from sklearn.utils.testing import ignore_warnings
 
-from sklearn.linear_model import logistic
-from sklearn import datasets
+from sklearn.linear_model.logistic import (
+    LogisticRegression,
+    logistic_regression_path, LogisticRegressionCV,
+    _logistic_loss_and_grad, _logistic_loss_grad_hess
+    )
+from sklearn.cross_validation import StratifiedKFold
+from sklearn.datasets import load_iris, make_classification
 
 X = [[-1, 0], [0, 1], [1, 1]]
 X_sp = sp.csr_matrix(X)
 Y1 = [0, 1, 1]
 Y2 = [2, 1, 0]
-iris = datasets.load_iris()
+iris = load_iris()
 
 
 def check_predictions(clf, X, y):
@@ -42,28 +50,26 @@ def test_predict_2_classes():
 
     Make sure it predicts the correct result on simple datasets.
     """
-    check_predictions(logistic.LogisticRegression(random_state=0), X, Y1)
-    check_predictions(logistic.LogisticRegression(random_state=0), X_sp, Y1)
+    check_predictions(LogisticRegression(random_state=0), X, Y1)
+    check_predictions(LogisticRegression(random_state=0), X_sp, Y1)
 
-    check_predictions(logistic.LogisticRegression(C=100, random_state=0),
-                      X, Y1)
-    check_predictions(logistic.LogisticRegression(C=100, random_state=0),
-                      X_sp, Y1)
+    check_predictions(LogisticRegression(C=100, random_state=0), X, Y1)
+    check_predictions(LogisticRegression(C=100, random_state=0), X_sp, Y1)
 
-    check_predictions(logistic.LogisticRegression(fit_intercept=False,
-                                                  random_state=0), X, Y1)
-    check_predictions(logistic.LogisticRegression(fit_intercept=False,
-                                                  random_state=0), X_sp, Y1)
+    check_predictions(LogisticRegression(fit_intercept=False,
+                                         random_state=0), X, Y1)
+    check_predictions(LogisticRegression(fit_intercept=False,
+                                         random_state=0), X_sp, Y1)
 
 
 def test_error():
     """Test for appropriate exception on errors"""
-    assert_raises(ValueError, logistic.LogisticRegression(C=-1).fit, X, Y1)
+    assert_raises(ValueError, LogisticRegression(C=-1).fit, X, Y1)
 
 
 def test_predict_3_classes():
-    check_predictions(logistic.LogisticRegression(C=10), X, Y2)
-    check_predictions(logistic.LogisticRegression(C=10), X_sp, Y2)
+    check_predictions(LogisticRegression(C=10), X, Y2)
+    check_predictions(LogisticRegression(C=10), X_sp, Y2)
 
 
 def test_predict_iris():
@@ -71,7 +77,7 @@ def test_predict_iris():
     n_samples, n_features = iris.data.shape
 
     target = iris.target_names[iris.target]
-    clf = logistic.LogisticRegression(C=len(iris.data)).fit(iris.data, target)
+    clf = LogisticRegression(C=len(iris.data)).fit(iris.data, target)
     assert_array_equal(np.unique(target), clf.classes_)
 
     pred = clf.predict(iris.data)
@@ -88,7 +94,7 @@ def test_sparsify():
     """Test sparsify and densify members."""
     n_samples, n_features = iris.data.shape
     target = iris.target_names[iris.target]
-    clf = logistic.LogisticRegression(random_state=0).fit(iris.data, target)
+    clf = LogisticRegression(random_state=0).fit(iris.data, target)
 
     pred_d_d = clf.decision_function(iris.data)
 
@@ -114,7 +120,7 @@ def test_inconsistent_input():
     y_ = np.ones(X_.shape[0])
     y_[0] = 0
 
-    clf = logistic.LogisticRegression(random_state=0)
+    clf = LogisticRegression(random_state=0)
 
     # Wrong dimensions for training data
     y_wrong = y_[:-1]
@@ -127,10 +133,7 @@ def test_inconsistent_input():
 
 def test_write_parameters():
     """Test that we can write to coef_ and intercept_"""
-    #rng = np.random.RandomState(0)
-    #X = rng.random_sample((5, 10))
-    #y = np.ones(X.shape[0])
-    clf = logistic.LogisticRegression(random_state=0)
+    clf = LogisticRegression(random_state=0)
     clf.fit(X, Y1)
     clf.coef_[:] = 0
     clf.intercept_[:] = 0
@@ -145,13 +148,277 @@ def test_nan():
     """
     Xnan = np.array(X, dtype=np.float64)
     Xnan[0, 1] = np.nan
-    logistic.LogisticRegression(random_state=0).fit(Xnan, Y1)
+    LogisticRegression(random_state=0).fit(Xnan, Y1)
+
+
+def test_consistency_path():
+    """Test that the path algorithm is consistent"""
+    rng = np.random.RandomState(0)
+    X = np.concatenate((rng.randn(100, 2) + [1, 1], rng.randn(100, 2)))
+    y = [1] * 100 + [-1] * 100
+    Cs = np.logspace(0, 4, 10)
+
+    f = ignore_warnings
+    # can't test with fit_intercept=True since LIBLINEAR
+    # penalizes the intercept
+    for method in ('lbfgs', 'newton-cg', 'liblinear'):
+        coefs, Cs = f(logistic_regression_path)(
+            X, y, Cs=Cs, fit_intercept=False, tol=1e-16, solver=method)
+        for i, C in enumerate(Cs):
+            lr = LogisticRegression(C=C, fit_intercept=False, tol=1e-16)
+            lr.fit(X, y)
+            lr_coef = lr.coef_.ravel()
+            assert_array_almost_equal(lr_coef, coefs[i], decimal=4)
+
+    # test for fit_intercept=True
+    for method in ('lbfgs', 'newton-cg', 'liblinear'):
+        Cs = [1e3]
+        coefs, Cs = f(logistic_regression_path)(
+            X, y, Cs=Cs, fit_intercept=True, tol=1e-4, solver=method)
+        lr = LogisticRegression(C=Cs[0], fit_intercept=True, tol=1e-4,
+                                intercept_scaling=10000)
+        lr.fit(X, y)
+        lr_coef = np.concatenate([lr.coef_.ravel(), lr.intercept_])
+        assert_array_almost_equal(lr_coef, coefs[0], decimal=4)
 
 
 def test_liblinear_random_state():
-    X, y = datasets.make_classification(n_samples=20)
-    lr1 = logistic.LogisticRegression(random_state=0)
+    X, y = make_classification(n_samples=20)
+    lr1 = LogisticRegression(random_state=0)
     lr1.fit(X, y)
-    lr2 = logistic.LogisticRegression(random_state=0)
+    lr2 = LogisticRegression(random_state=0)
     lr2.fit(X, y)
     assert_array_almost_equal(lr1.coef_, lr2.coef_)
+
+
+def test_logistic_loss_and_grad():
+    X_ref, y = make_classification(n_samples=20)
+    n_features = X_ref.shape[1]
+
+    X_sp = X_ref.copy()
+    X_sp[X_sp < .1] = 0
+    X_sp = sp.csr_matrix(X_sp)
+    for X in (X_ref, X_sp):
+        w = np.zeros(n_features)
+
+        # First check that our derivation of the grad is correct
+        loss, grad = _logistic_loss_and_grad(w, X, y, alpha=1.)
+        approx_grad = optimize.approx_fprime(
+            w, lambda w: _logistic_loss_and_grad(w, X, y, alpha=1.)[0], 1e-3
+            )
+        assert_array_almost_equal(grad, approx_grad, decimal=2)
+
+        # Second check that our intercept implementation is good
+        w = np.zeros(n_features + 1)
+        loss_interp, grad_interp = _logistic_loss_and_grad(
+            w, X, y, alpha=1.
+            )
+        assert_array_almost_equal(loss, loss_interp)
+
+        approx_grad = optimize.approx_fprime(
+            w, lambda w: _logistic_loss_and_grad(w, X, y, alpha=1.)[0], 1e-3
+            )
+        assert_array_almost_equal(grad_interp, approx_grad, decimal=2)
+
+
+def test_logistic_loss_grad_hess():
+    rng = np.random.RandomState(0)
+    n_samples, n_features = 50, 5
+    X_ref = rng.randn(n_samples, n_features)
+    y = np.sign(X_ref.dot(5 * rng.randn(n_features)))
+    X_ref -= X_ref.mean()
+    X_ref /= X_ref.std()
+    X_sp = X_ref.copy()
+    X_sp[X_sp < .1] = 0
+    X_sp = sp.csr_matrix(X_sp)
+    for X in (X_ref, X_sp):
+        w = .1 * np.ones(n_features)
+
+        # First check that _logistic_loss_grad_hess is consistent
+        # with _logistic_loss_and_grad
+        loss, grad = _logistic_loss_and_grad(w, X, y, alpha=1.)
+        loss_2, grad_2, hess = _logistic_loss_grad_hess(w, X, y, alpha=1.)
+        assert_array_almost_equal(grad, grad_2)
+
+        # Now check our hessian along the second direction of the grad
+        vector = np.zeros_like(grad)
+        vector[1] = 1
+        hess_col = hess(vector)
+
+        # Computation of the Hessian is particularly fragile to numerical
+        # errors when doing simple finite differences. Here we compute the
+        # grad along a path in the direction of the vector and then use a
+        # least-square regression to estimate the slope
+        e = 1e-3
+        d_x = np.linspace(-e, e, 30)
+        d_grad = np.array([
+            _logistic_loss_and_grad(w + t * vector, X, y, alpha=1.)[1]
+            for t in d_x
+            ])
+
+        d_grad -= d_grad.mean(axis=0)
+        approx_hess_col = linalg.lstsq(d_x[:, np.newaxis], d_grad)[0].ravel()
+
+        assert_array_almost_equal(approx_hess_col, hess_col, decimal=3)
+
+        # Second check that our intercept implementation is good
+        w = np.zeros(n_features + 1)
+        loss_interp, grad_interp = _logistic_loss_and_grad(
+            w, X, y, alpha=1.
+            )
+        loss_interp_2, grad_interp_2, hess = \
+            _logistic_loss_grad_hess(w, X, y, alpha=1.)
+        assert_array_almost_equal(loss_interp, loss_interp_2)
+        assert_array_almost_equal(grad_interp, grad_interp_2)
+
+
+def test_logistic_cv():
+    """test for LogisticRegressionCV object"""
+    n_samples, n_features = 50, 5
+    rng = np.random.RandomState(0)
+    X_ref = rng.randn(n_samples, n_features)
+    y = np.sign(X_ref.dot(5 * rng.randn(n_features)))
+    X_ref -= X_ref.mean()
+    X_ref /= X_ref.std()
+    lr_cv = LogisticRegressionCV(Cs=[1.], fit_intercept=False,
+                                 solver='liblinear')
+    lr_cv.fit(X_ref, y)
+    lr = LogisticRegression(C=1., fit_intercept=False)
+    lr.fit(X_ref, y)
+    assert_array_almost_equal(lr.coef_, lr_cv.coef_)
+
+    assert_array_equal(lr_cv.coef_.shape, (1, n_features))
+    assert_array_equal(lr_cv.classes_, [-1, 1])
+    assert_equal(len(lr_cv.classes_), 2)
+
+    coefs_paths = np.asarray(list(lr_cv.coefs_paths_.values()))
+    assert_array_equal(coefs_paths.shape, (1, 3, 1, n_features))
+    assert_array_equal(lr_cv.Cs_.shape, (1, ))
+    scores = np.asarray(list(lr_cv.scores_.values()))
+    assert_array_equal(scores.shape, (1, 3, 1))
+
+
+def test_logistic_cv_sparse():
+    X, y = make_classification(n_samples=50, n_features=5,
+                               random_state=0)
+    X[X < 1.0] = 0.0
+    csr = sp.csr_matrix(X)
+
+    clf = LogisticRegressionCV(fit_intercept=True)
+    clf.fit(X, y)
+    clfs = LogisticRegressionCV(fit_intercept=True)
+    clfs.fit(csr, y)
+    assert_array_almost_equal(clfs.coef_, clf.coef_)
+    assert_array_almost_equal(clfs.intercept_, clf.intercept_)
+    assert_equal(clfs.C_, clf.C_)
+
+
+def test_intercept_logistic_helper():
+    n_samples, n_features = 10, 5
+    X, y = make_classification(n_samples=n_samples, n_features=n_features,
+                               random_state=0)
+
+    # Fit intercept case.
+    alpha = 1.
+    w = np.ones(n_features + 1)
+    loss_interp, grad_interp, hess_interp = _logistic_loss_grad_hess(
+        w, X, y, alpha)
+
+    # Do not fit intercept. This can be considered equivalent to adding
+    # a feature vector of ones, i.e column of one vectors.
+    X_ = np.hstack((X, np.ones(10)[:, np.newaxis]))
+    loss, grad, hess = _logistic_loss_grad_hess(w, X_, y, alpha)
+
+    # In the fit_intercept=False case, the feature vector of ones is
+    # penalized. This should be taken care of.
+    assert_almost_equal(loss_interp + 0.5 * (w[-1] ** 2), loss)
+
+    # Check gradient.
+    assert_array_almost_equal(grad_interp[:n_features], grad[:n_features])
+    assert_almost_equal(grad_interp[-1] + alpha * w[-1], grad[-1])
+
+    rng = np.random.RandomState(0)
+    grad = rng.rand(n_features + 1)
+    hess_interp = hess_interp(grad)
+    hess = hess(grad)
+    assert_array_almost_equal(hess_interp[:n_features], hess[:n_features])
+    assert_almost_equal(hess_interp[-1] + alpha * grad[-1], hess[-1])
+
+
+def test_ova_iris():
+    """Test that our OvA implementation is correct using the iris dataset."""
+    train, target = iris.data, iris.target
+    n_samples, n_features = train.shape
+
+    # Use pre-defined fold as folds generated for different y
+    cv = StratifiedKFold(target, 3)
+    clf = LogisticRegressionCV(cv=cv)
+    clf.fit(train, target)
+
+    clf1 = LogisticRegressionCV(cv=cv)
+    target[target == 0] = 1
+    clf1.fit(train, target)
+
+    assert_array_equal(clf.scores_[2], clf1.scores_[2])
+    assert_array_equal(clf.intercept_[2:], clf1.intercept_)
+    assert_array_equal(clf.coef_[2][np.newaxis, :], clf1.coef_)
+
+    # Test the shape of various attributes.
+    assert_array_equal(clf.coef_.shape, (3, n_features))
+    assert_array_equal(clf.classes_, [0, 1, 2])
+    assert_equal(len(clf.classes_), 3)
+
+    coefs_paths = np.asarray(list(clf.coefs_paths_.values()))
+    assert_array_equal(coefs_paths.shape, (3, 3, 10, n_features + 1))
+    assert_array_equal(clf.Cs_.shape, (10, ))
+    scores = np.asarray(list(clf.scores_.values()))
+    assert_array_equal(scores.shape, (3, 3, 10))
+
+
+def test_logistic_regression_solvers():
+    X, y = make_classification(n_features=10, n_informative=5, random_state=0)
+    clf_n = LogisticRegression(solver='newton-cg', fit_intercept=False)
+    clf_n.fit(X, y)
+    clf_lbf = LogisticRegression(solver='lbfgs', fit_intercept=False)
+    clf_lbf.fit(X, y)
+    clf_lib = LogisticRegression(fit_intercept=False)
+    clf_lib.fit(X, y)
+    assert_array_almost_equal(clf_n.coef_, clf_lib.coef_, decimal=3)
+    assert_array_almost_equal(clf_lib.coef_, clf_lbf.coef_, decimal=3)
+    assert_array_almost_equal(clf_n.coef_, clf_lbf.coef_, decimal=3)
+
+
+def test_logistic_regression_solvers_multiclass():
+    X, y = make_classification(n_samples=20, n_features=20, n_informative=10,
+                               n_classes=3, random_state=0)
+    clf_n = LogisticRegression(solver='newton-cg', fit_intercept=False)
+    clf_n.fit(X, y)
+    clf_lbf = LogisticRegression(solver='lbfgs', fit_intercept=False)
+    clf_lbf.fit(X, y)
+    clf_lib = LogisticRegression(fit_intercept=False)
+    clf_lib.fit(X, y)
+    assert_array_almost_equal(clf_n.coef_, clf_lib.coef_, decimal=4)
+    assert_array_almost_equal(clf_lib.coef_, clf_lbf.coef_, decimal=4)
+    assert_array_almost_equal(clf_n.coef_, clf_lbf.coef_, decimal=4)
+
+
+def test_logistic_regressioncv_class_weights():
+    X, y = make_classification(n_samples=20, n_features=20, n_informative=10,
+                               n_classes=3, random_state=0)
+
+    # Test the liblinear fails when class_weight of type dict is
+    # provided, when it is multiclass
+    clf_lib = LogisticRegressionCV(class_weight={0: 0.1, 1: 0.2},
+                                   solver='liblinear')
+    assert_raises(ValueError, clf_lib.fit, X, y)
+
+    # Test for class_weight=auto
+    X, y = make_classification(n_samples=20, n_features=20, n_informative=10,
+                               random_state=0)
+    clf_lbf = LogisticRegressionCV(solver='lbfgs', fit_intercept=False,
+                                   class_weight='auto')
+    clf_lbf.fit(X, y)
+    clf_lib = LogisticRegressionCV(solver='liblinear', fit_intercept=False,
+                                   class_weight='auto')
+    clf_lib.fit(X, y)
+    assert_array_almost_equal(clf_lib.coef_, clf_lbf.coef_, decimal=4)
diff --git a/sklearn/svm/base.py b/sklearn/svm/base.py
index c7934c48381..8a92fe62c31 100644
--- a/sklearn/svm/base.py
+++ b/sklearn/svm/base.py
@@ -601,7 +601,8 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)):
     @abstractmethod
     def __init__(self, penalty='l2', loss='l2', dual=True, tol=1e-4, C=1.0,
                  multi_class='ovr', fit_intercept=True, intercept_scaling=1,
-                 class_weight=None, verbose=0, random_state=None):
+                 class_weight=None, verbose=0, random_state=None, max_iter=100,
+                 solver='liblinear'):
 
         self.penalty = penalty
         self.loss = loss
@@ -614,6 +615,8 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)):
         self.class_weight = class_weight
         self.verbose = verbose
         self.random_state = random_state
+        self.solver = solver
+        self.max_iter = max_iter
 
         # Check that the arguments given are valid:
         self._get_solver_type()
@@ -670,6 +673,19 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)):
         self : object
             Returns self.
         """
+
+        # Circular import, logistic_regression_path depends on LogisticRegression
+        # and hence BaseLibLinear.
+        from ..linear_model import logistic_regression_path
+
+        if self.solver != 'liblinear':
+            if self.penalty != 'l2':
+                raise ValueError("newton-cg and lbfgs solvers support only "
+                                 "l2 penalties.")
+            if self.dual:
+                raise ValueError("newton-cg and lbfgs solvers support only "
+                                 "the primal form.")
+
         self._enc = LabelEncoder()
         y_ind = self._enc.fit_transform(y)
         if len(self.classes_) < 2:
@@ -678,6 +694,7 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)):
 
         X = check_array(X, accept_sparse='csr', dtype=np.float64, order="C")
 
+        # Used in the liblinear solver.
         self.class_weight_ = compute_class_weight(self.class_weight,
                                                   self.classes_, y)
 
@@ -686,31 +703,67 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)):
                              "X has %s samples, but y has %s." %
                              (X.shape[0], y_ind.shape[0]))
 
-        liblinear.set_verbosity_wrap(self.verbose)
+        if self.solver not in ['liblinear', 'newton-cg', 'lbfgs']:
+            raise ValueError("Logistic Regression supports only liblinear,"
+                             " newton-cg and lbfgs solvers.")
 
-        rnd = check_random_state(self.random_state)
-        if self.verbose:
-            print('[LibLinear]', end='')
+        if self.solver == 'liblinear':
+            liblinear.set_verbosity_wrap(self.verbose)
 
-        # LibLinear wants targets as doubles, even for classification
-        y_ind = np.asarray(y_ind, dtype=np.float64).ravel()
-        raw_coef_ = liblinear.train_wrap(X, y_ind,
-                                         sp.isspmatrix(X),
-                                         self._get_solver_type(),
-                                         self.tol, self._get_bias(),
-                                         self.C, self.class_weight_,
-                                         rnd.randint(np.iinfo('i').max))
-        # Regarding rnd.randint(..) in the above signature:
-        # seed for srand in range [0..INT_MAX); due to limitations in Numpy
-        # on 32-bit platforms, we can't get to the UINT_MAX limit that
-        # srand supports
+            rnd = check_random_state(self.random_state)
+            if self.verbose:
+                print('[LibLinear]', end='')
+
+            # LibLinear wants targets as doubles, even for classification
+            y_ind = np.asarray(y_ind, dtype=np.float64).ravel()
+            raw_coef_ = liblinear.train_wrap(X, y_ind,
+                                             sp.isspmatrix(X),
+                                             self._get_solver_type(),
+                                             self.tol, self._get_bias(),
+                                             self.C, self.class_weight_,
+                                             rnd.randint(np.iinfo('i').max))
+            # Regarding rnd.randint(..) in the above signature:
+            # seed for srand in range [0..INT_MAX); due to limitations in Numpy
+            # on 32-bit platforms, we can't get to the UINT_MAX limit that
+            # srand supports
+
+            if self.fit_intercept:
+                self.coef_ = raw_coef_[:, :-1]
+                self.intercept_ = self.intercept_scaling * raw_coef_[:, -1]
+            else:
+                self.coef_ = raw_coef_
+                self.intercept_ = 0.
 
-        if self.fit_intercept:
-            self.coef_ = raw_coef_[:, :-1]
-            self.intercept_ = self.intercept_scaling * raw_coef_[:, -1]
         else:
-            self.coef_ = raw_coef_
-            self.intercept_ = 0.
+            if self.penalty != 'l2':
+                raise ValueError("newton-cg and lbfgs solvers support only "
+                                 "l2 penalties.")
+
+            n_tasks = len(self.classes_)
+            classes_ = self.classes_
+
+            if len(self.classes_) == 2:
+                n_tasks = 1
+                classes_ = classes_[1:]
+
+            self.coef_ = np.empty((n_tasks, X.shape[1]))
+            self.intercept_ = np.zeros(n_tasks)
+
+            for ind, class_ in enumerate(classes_):
+                coef_, _ = logistic_regression_path(
+                    X, y, pos_class=class_, Cs=[self.C],
+                    fit_intercept=self.fit_intercept,
+                    tol=self.tol, verbose=self.verbose,
+                    solver=self.solver, copy=True,
+                    max_iter=self.max_iter, class_weight=self.class_weight)
+
+                coef_ = coef_[0]
+                if self.fit_intercept:
+                    self.coef_[ind] = coef_[:-1]
+                    self.intercept_[ind] = coef_[-1]
+
+                else:
+                    self.coef_[ind] = coef_
 
         if self.multi_class == "crammer_singer" and len(self.classes_) == 2:
             self.coef_ = (self.coef_[1] - self.coef_[0]).reshape(1, -1)
diff --git a/sklearn/tests/test_common.py b/sklearn/tests/test_common.py
index ddc86c8afbc..920aab6bba2 100644
--- a/sklearn/tests/test_common.py
+++ b/sklearn/tests/test_common.py
@@ -319,7 +319,9 @@ def test_non_transformer_estimators_n_iter():
                 # These models are dependent on external solvers like
                 # libsvm and accessing the iter parameter is non-trivial.
                 if name in (['Ridge', 'SVR', 'NuSVR', 'NuSVC',
-                             'RidgeClassifier', 'SVC', 'RandomizedLasso']):
+                             'RidgeClassifier', 'SVC', 'RandomizedLasso',
+                             'LogisticRegressionCV', 'LogisticRegression',
+                             'LinearSVC']):
                     continue
 
                 # Tested in test_transformer_n_iter below
@@ -340,6 +342,7 @@ def test_transformer_n_iter():
         # Dependent on external solvers and hence accessing the iter
         # param is non-trivial.
         external_solver = ['Isomap', 'KernelPCA', 'LocallyLinearEmbedding',
-                           'RandomizedLasso']
+                           'RandomizedLasso', 'LogisticRegression',
+                           'LogisticRegressionCV', 'LinearSVC']
         if hasattr(estimator, "max_iter") and name not in external_solver:
             yield check_transformer_n_iter, name, estimator
diff --git a/sklearn/utils/estimator_checks.py b/sklearn/utils/estimator_checks.py
index 0ce6793d40f..874dae7c338 100644
--- a/sklearn/utils/estimator_checks.py
+++ b/sklearn/utils/estimator_checks.py
@@ -771,8 +771,9 @@ def check_cluster_overwrite_params(name, Clustering):
 
 
 def check_sparsify_multiclass_classifier(name, Classifier):
-    X = np.array([[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1]])
-    y = [1, 1, 1, 2, 2, 3]
+    X = np.array([[-2, -1], [-1, -1], [-1, -2], [1, 1], [1, 2], [2, 1],
+                  [-1, -2], [2, 2], [-2, -2]])
+    y = [1, 1, 1, 2, 2, 2, 3, 3, 3]
     est = Classifier()
 
     est.fit(X, y)
diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py
new file mode 100644
index 00000000000..c1f8403438d
--- /dev/null
+++ b/sklearn/utils/optimize.py
@@ -0,0 +1,176 @@
+"""
+Our own implementation of the Newton algorithm
+
+Unlike the scipy.optimize version, this version of the Newton conjugate
+gradient solver uses only one function call to retrieve the
+func value, the gradient value and a callable for the Hessian matvec
+product. If the function call is very expensive (e.g. for logistic
+regression with large design matrix), this approach gives very
+significant speedups.
+"""
+# This is a modified file from scipy.optimize
+# Original authors: Travis Oliphant, Eric Jones
+# Modifications by Gael Varoquaux
+# License: BSD
+
+import numpy as np
+import warnings
+from scipy.optimize.linesearch import line_search_wolfe2, line_search_wolfe1
+
+
+class _LineSearchError(RuntimeError):
+    pass
+
+
+def _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval,
+                         **kwargs):
+    """
+    Same as line_search_wolfe1, but fall back to line_search_wolfe2 if
+    suitable step length is not found, and raise an exception if a
+    suitable step length is not found.
+
+    Raises
+    ------
+    _LineSearchError
+        If no suitable step size is found
+
+    """
+    ret = line_search_wolfe1(f, fprime, xk, pk, gfk,
+                             old_fval, old_old_fval,
+                             **kwargs)
+
+    if ret[0] is None:
+        # line search failed: try different one.
+        ret = line_search_wolfe2(f, fprime, xk, pk, gfk,
+                                 old_fval, old_old_fval, **kwargs)
+
+    if ret[0] is None:
+        raise _LineSearchError()
+
+    return ret
+
+
+def newton_cg(func_grad_hess, func, grad, x0, args=(), eps=1e-4, tol=1e-4,
+              maxiter=100):
+    """
+    Minimization of scalar function of one or more variables using the
+    Newton-CG algorithm.
+
+    Parameters
+    ----------
+    func_grad_hess : callable
+        Should return the value of the function, the gradient, and a
+        callable returning the matvec product of the Hessian.
+
+    func : callable
+        Should return the value of the function.
+
+    grad : callable
+        Should return the function value and the gradient. This is used
+        by the linesearch functions.
+
+    x0 : float
+        Initial guess.
+
+    args: tuple, optional
+        Arguments passed to func_grad_hess, func and grad.
+
+    tol : float
+        Stopping criterion. The iteration will stop when
+        ``max{|g_i | i = 1, ..., n} <= tol``
+        where ``g_i`` is the i-th component of the gradient.
+
+    eps : float, optional
+        If fhess is approximated, use this value for the step size.
+
+    maxiter : int
+        Number of iterations.
+
+    Returns
+    -------
+    xk : float
+        Estimated minimum.
+    """
+    x0 = np.asarray(x0).flatten()
+    xk = x0
+    k = 1
+    old_fval = func(x0, *args)
+    old_old_fval = None
+
+    # Outer loop: our Newton iteration
+    while k <= maxiter:
+
+        # Compute a search direction pk by applying the CG method to
+        #  del2 f(xk) p = - fgrad f(xk) starting from 0.
+        fval, fgrad, fhess_p = func_grad_hess(xk, *args)
+
+        absgrad = np.abs(fgrad)
+        if np.max(absgrad) < tol:
+            break
+
+        maggrad = np.sum(absgrad)
+        eta = min([0.5, np.sqrt(maggrad)])
+        termcond = eta * maggrad
+        xsupi = np.zeros(len(x0), dtype=x0.dtype)
+        ri = fgrad
+        psupi = -ri
+        i = 0
+        dri0 = np.dot(ri, ri)
+
+        # Inner loop: solve the Newton update by conjugate gradient, to
+        # avoid inverting the Hessian
+        while np.sum(np.abs(ri)) > termcond:
+            Ap = fhess_p(psupi)
+            # check curvature
+            curv = np.dot(psupi, Ap)
+            if 0 <= curv <= 3*np.finfo(np.float64).eps:
+                break
+            elif curv < 0:
+                if (i > 0):
+                    break
+                else:
+                    # fall back to steepest descent direction
+                    xsupi = xsupi + dri0 / curv * psupi
+                    break
+            alphai = dri0 / curv
+            xsupi = xsupi + alphai * psupi
+            ri = ri + alphai * Ap
+            dri1 = np.dot(ri, ri)
+            betai = dri1 / dri0
+            psupi = -ri + betai * psupi
+            i = i + 1
+            dri0 = dri1          # update np.dot(ri,ri) for next time.
+
+        try:
+            alphak, fc, gc, old_fval, old_old_fval, gfkp1 = \
+                _line_search_wolfe12(func, grad, xk, xsupi, fgrad,
+                                     old_fval, old_old_fval, args=args)
+        except _LineSearchError:
+            warnings.warn('Line Search failed')
+            break
+
+        update = alphak * xsupi
+        xk = xk + update        # upcast if necessary
+        k += 1
+
+    if k > maxiter:
+        warnings.warn("newton-cg failed to converge. Increase the "
+                      "number of iterations.")
+    return xk
+
+
+###############################################################################
+# Tests
+
+if __name__ == "__main__":
+    A = np.random.normal(size=(10, 10))
+
+    def func(x):
+        Ax = A.dot(x)
+        return .5*(Ax).dot(Ax)
+
+    def func_grad_hess(x):
+        return func(x), A.T.dot(A.dot(x)), lambda x: A.T.dot(A.dot(x))
+
+    x0 = np.ones(10)
+    out = newton_cg(func_grad_hess, func, x0)