"""

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>

# Lars Buitinck

# Simon Wu <s8wu@uwaterloo.ca>

# Arthur Mensch <arthur.mensch@m4x.org

import numbers

import warnings

import numpy as np

from scipy import optimize, sparse

from scipy.special import expit, logsumexp

from joblib import Parallel, delayed, effective_n_jobs

from ._base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator

from ._sag import sag_solver

from ..preprocessing import LabelEncoder, LabelBinarizer

from ..svm._base import _fit_liblinear

from ..utils import check_array, check_consistent_length, compute_class_weight

from ..utils import check_random_state

from ..utils.extmath import (log_logistic, safe_sparse_dot, softmax,

squared_norm)

from ..utils.extmath import row_norms

from ..utils.optimize import _newton_cg, _check_optimize_result

from ..utils.validation import check_is_fitted, _check_sample_weight

from ..utils.validation import _deprecate_positional_args

from ..utils.multiclass import check_classification_targets

from ..utils.fixes import _joblib_parallel_args

from ..model_selection import check_cv

from ..metrics import get_scorer

_LOGISTIC_SOLVER_CONVERGENCE_MSG = (

"Please also refer to the documentation for alternative solver options:\n"

" https://scikit-learn.org/stable/modules/linear_model.html"

"#logistic-regression")

# .. 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 of shape (n_features,) or (n_features + 1,)

Coefficient vector.

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

y : ndarray of shape (n_samples,)

Array of labels.

Returns

-------

w : ndarray of shape (n_features,)

Coefficient vector without the intercept weight (w[-1]) if the

intercept should be fit. Unchanged otherwise.

c : float

The intercept.

yz : float

y * np.dot(X, w).

"""

c = 0.

if w.size == X.shape[1] + 1:

c = w[-1]

w = w[:-1]

z = safe_sparse_dot(X, w) + c

yz = y * z

return w, c, yz

def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None):

"""Computes the logistic loss and gradient.

Parameters

----------

w : ndarray of shape (n_features,) or (n_features + 1,)

Coefficient vector.

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

y : ndarray of shape (n_samples,)

Array of labels.

alpha : float

Regularization parameter. alpha is equal to 1 / C.

sample_weight : array-like of shape (n_samples,), default=None

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 of shape (n_features,) or (n_features + 1,)

Logistic gradient.

"""

n_samples, 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(n_samples)

# 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 of shape (n_features,) or (n_features + 1,)

Coefficient vector.

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

y : ndarray of shape (n_samples,)

Array of labels.

alpha : float

Regularization parameter. alpha is equal to 1 / C.

sample_weight : array-like of shape (n_samples,) default=None

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_grad_hess(w, X, y, alpha, sample_weight=None):

"""Computes the gradient and the Hessian, in the case of a logistic loss.

Parameters

----------

w : ndarray of shape (n_features,) or (n_features + 1,)

Coefficient vector.

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

y : ndarray of shape (n_samples,)

Array of labels.

alpha : float

Regularization parameter. alpha is equal to 1 / C.

sample_weight : array-like of shape (n_samples,) default=None

Array of weights that are assigned to individual samples.

If not provided, then each sample is given unit weight.

Returns

-------

grad : ndarray of 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])

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 grad, Hs

def _multinomial_loss(w, X, Y, alpha, sample_weight):

"""Computes multinomial loss and class probabilities.

Parameters

----------

w : ndarray of shape (n_classes * n_features,) or

(n_classes * (n_features + 1),)

Coefficient vector.

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

Y : ndarray of shape (n_samples, n_classes)

Transformed labels according to the output of LabelBinarizer.

alpha : float

Regularization parameter. alpha is equal to 1 / C.

sample_weight : array-like of shape (n_samples,)

Array of weights that are assigned to individual samples.

Returns

-------

loss : float

Multinomial loss.

p : ndarray of shape (n_samples, n_classes)

Estimated class probabilities.

w : ndarray of shape (n_classes, n_features)

Reshaped param vector excluding intercept terms.

Reference

---------

Bishop, C. M. (2006). Pattern recognition and machine learning.

Springer. (Chapter 4.3.4)

"""

n_classes = Y.shape[1]

n_features = X.shape[1]

fit_intercept = w.size == (n_classes * (n_features + 1))

w = w.reshape(n_classes, -1)

sample_weight = sample_weight[:, np.newaxis]

if fit_intercept:

intercept = w[:, -1]

w = w[:, :-1]

else:

intercept = 0

p = safe_sparse_dot(X, w.T)

p += intercept

p -= logsumexp(p, axis=1)[:, np.newaxis]

loss = -(sample_weight * Y * p).sum()

loss += 0.5 * alpha * squared_norm(w)

p = np.exp(p, p)

return loss, p, w

def _multinomial_loss_grad(w, X, Y, alpha, sample_weight):

"""Computes the multinomial loss, gradient and class probabilities.

Parameters

----------

w : ndarray of shape (n_classes * n_features,) or

(n_classes * (n_features + 1),)

Coefficient vector.

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

Y : ndarray of shape (n_samples, n_classes)

Transformed labels according to the output of LabelBinarizer.

alpha : float

Regularization parameter. alpha is equal to 1 / C.

sample_weight : array-like of shape (n_samples,)

Array of weights that are assigned to individual samples.

Returns

-------

loss : float

Multinomial loss.

grad : ndarray of shape (n_classes * n_features,) or \

(n_classes * (n_features + 1),)

Ravelled gradient of the multinomial loss.

p : ndarray of shape (n_samples, n_classes)

Estimated class probabilities

Reference

---------

Bishop, C. M. (2006). Pattern recognition and machine learning.

Springer. (Chapter 4.3.4)

"""

n_classes = Y.shape[1]

n_features = X.shape[1]

fit_intercept = (w.size == n_classes * (n_features + 1))

grad = np.zeros((n_classes, n_features + bool(fit_intercept)),

dtype=X.dtype)

loss, p, w = _multinomial_loss(w, X, Y, alpha, sample_weight)

sample_weight = sample_weight[:, np.newaxis]

diff = sample_weight * (p - Y)

grad[:, :n_features] = safe_sparse_dot(diff.T, X)

grad[:, :n_features] += alpha * w

if fit_intercept:

grad[:, -1] = diff.sum(axis=0)

return loss, grad.ravel(), p

def _multinomial_grad_hess(w, X, Y, alpha, sample_weight):

"""

Computes the gradient and the Hessian, in the case of a multinomial loss.

Parameters

----------

w : ndarray of shape (n_classes * n_features,) or

(n_classes * (n_features + 1),)

Coefficient vector.

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

Y : ndarray of shape (n_samples, n_classes)

Transformed labels according to the output of LabelBinarizer.

alpha : float

Regularization parameter. alpha is equal to 1 / C.

sample_weight : array-like of shape (n_samples,)

Array of weights that are assigned to individual samples.

Returns

-------

grad : ndarray of shape (n_classes * n_features,) or \

(n_classes * (n_features + 1),)

Ravelled gradient of the multinomial loss.

hessp : callable

Function that takes in a vector input of shape (n_classes * n_features)

or (n_classes * (n_features + 1)) and returns matrix-vector product

with hessian.

References

----------

Barak A. Pearlmutter (1993). Fast Exact Multiplication by the Hessian.

http://www.bcl.hamilton.ie/~barak/papers/nc-hessian.pdf

"""

n_features = X.shape[1]

n_classes = Y.shape[1]

fit_intercept = w.size == (n_classes * (n_features + 1))

# `loss` is unused. Refactoring to avoid computing it does not

# significantly speed up the computation and decreases readability

loss, grad, p = _multinomial_loss_grad(w, X, Y, alpha, sample_weight)

sample_weight = sample_weight[:, np.newaxis]

# Hessian-vector product derived by applying the R-operator on the gradient

# of the multinomial loss function.

def hessp(v):

v = v.reshape(n_classes, -1)

if fit_intercept:

inter_terms = v[:, -1]

v = v[:, :-1]

else:

inter_terms = 0

# r_yhat holds the result of applying the R-operator on the multinomial

# estimator.

r_yhat = safe_sparse_dot(X, v.T)

r_yhat += inter_terms

r_yhat += (-p * r_yhat).sum(axis=1)[:, np.newaxis]

r_yhat *= p

r_yhat *= sample_weight

hessProd = np.zeros((n_classes, n_features + bool(fit_intercept)))

hessProd[:, :n_features] = safe_sparse_dot(r_yhat.T, X)

hessProd[:, :n_features] += v * alpha

if fit_intercept:

hessProd[:, -1] = r_yhat.sum(axis=0)

return hessProd.ravel()

return grad, hessp

def _check_solver(solver, penalty, dual):

all_solvers = ['liblinear', 'newton-cg', 'lbfgs', 'sag', 'saga']

if solver not in all_solvers:

raise ValueError("Logistic Regression supports only solvers in %s, got"

" %s." % (all_solvers, solver))

all_penalties = ['l1', 'l2', 'elasticnet', 'none']

if penalty not in all_penalties:

raise ValueError("Logistic Regression supports only penalties in %s,"

" got %s." % (all_penalties, penalty))

if solver not in ['liblinear', 'saga'] and penalty not in ('l2', 'none'):

raise ValueError("Solver %s supports only 'l2' or 'none' penalties, "

"got %s penalty." % (solver, penalty))

if solver != 'liblinear' and dual:

raise ValueError("Solver %s supports only "

"dual=False, got dual=%s" % (solver, dual))

if penalty == 'elasticnet' and solver != 'saga':

raise ValueError("Only 'saga' solver supports elasticnet penalty,"

" got solver={}.".format(solver))

if solver == 'liblinear' and penalty == 'none':

raise ValueError(

"penalty='none' is not supported for the liblinear solver"

)

return solver

def _check_multi_class(multi_class, solver, n_classes):

if multi_class == 'auto':

if solver == 'liblinear':

multi_class = 'ovr'

elif n_classes > 2:

multi_class = 'multinomial'

else:

multi_class = 'ovr'

if multi_class not in ('multinomial', 'ovr'):

raise ValueError("multi_class should be 'multinomial', 'ovr' or "

"'auto'. Got %s." % multi_class)

if multi_class == 'multinomial' and solver == 'liblinear':

raise ValueError("Solver %s does not support "

"a multinomial backend." % solver)

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='auto',

random_state=None, check_input=True,

max_squared_sum=None, sample_weight=None,

l1_ratio=None):

"""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.

Note that there will be no speedup with liblinear solver, since it does

not handle warm-starting.

Read more in the :ref:`User Guide <logistic_regression>`.

Parameters

----------

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Input data.

y : array-like of shape (n_samples,) or (n_samples, n_targets)

Input data, target values.

pos_class : int, default=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 : int or array-like of shape (n_cs,), default=10

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.

fit_intercept : bool, default=True

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, default=100

Maximum number of iterations for the solver.

tol : float, default=1e-4

Stopping criterion. For the newton-cg and lbfgs solvers, the iteration

will stop when ``max{|g_i | i = 1, ..., n} <= tol``

where ``g_i`` is the i-th component of the gradient.

verbose : int, default=0

For the liblinear and lbfgs solvers set verbose to any positive

number for verbosity.

solver : {'lbfgs', 'newton-cg', 'liblinear', 'sag', 'saga'}, \

default='lbfgs'

Numerical solver to use.

coef : array-like of shape (n_features,), default=None

Initialization value for coefficients of logistic regression.

Useless for liblinear solver.

class_weight : dict or 'balanced', default=None

Weights associated with classes in the form ``{class_label: weight}``.

If not given, all classes are supposed to have weight one.

The "balanced" mode uses the values of y to automatically adjust

weights inversely proportional to class frequencies in the input data

as ``n_samples / (n_classes * np.bincount(y))``.

Note that these weights will be multiplied with sample_weight (passed

through the fit method) if sample_weight is specified.

dual : bool, default=False

Dual or primal formulation. Dual formulation is only implemented for

l2 penalty with liblinear solver. Prefer dual=False when

n_samples > n_features.

penalty : {'l1', 'l2', 'elasticnet'}, default='l2'

Used to specify the norm used in the penalization. The 'newton-cg',

'sag' and 'lbfgs' solvers support only l2 penalties. 'elasticnet' is

only supported by the 'saga' solver.

intercept_scaling : float, default=1.

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 equal 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.

multi_class : {'ovr', 'multinomial', 'auto'}, default='auto'

If the option chosen is 'ovr', then a binary problem is fit for each

label. For 'multinomial' the loss minimised is the multinomial loss fit

across the entire probability distribution, *even when the data is

binary*. 'multinomial' is unavailable when solver='liblinear'.

'auto' selects 'ovr' if the data is binary, or if solver='liblinear',

and otherwise selects 'multinomial'.

.. versionadded:: 0.18

Stochastic Average Gradient descent solver for 'multinomial' case.

.. versionchanged:: 0.22

Default changed from 'ovr' to 'auto' in 0.22.

random_state : int, RandomState instance, default=None

Used when ``solver`` == 'sag', 'saga' or 'liblinear' to shuffle the

data. See :term:`Glossary <random_state>` for details.

check_input : bool, default=True

If False, the input arrays X and y will not be checked.

max_squared_sum : float, default=None

Maximum squared sum of X over samples. Used only in SAG solver.

If None, it will be computed, going through all the samples.

The value should be precomputed to speed up cross validation.

sample_weight : array-like of shape(n_samples,), default=None

Array of weights that are assigned to individual samples.

If not provided, then each sample is given unit weight.

l1_ratio : float, default=None

The Elastic-Net mixing parameter, with ``0 <= l1_ratio <= 1``. Only

used if ``penalty='elasticnet'``. 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 of 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. For

``multiclass='multinomial'``, the shape is (n_classes, n_cs,

n_features) or (n_classes, n_cs, n_features + 1).

Cs : ndarray

Grid of Cs used for cross-validation.

n_iter : array of shape (n_cs,)

Actual number of iteration for each Cs.

Notes

-----

You might get slightly different results with the solver liblinear than

with the others since this uses LIBLINEAR which penalizes the intercept.

.. versionchanged:: 0.19

The "copy" parameter was removed.

"""

if isinstance(Cs, numbers.Integral):

Cs = np.logspace(-4, 4, Cs)

solver = _check_solver(solver, penalty, dual)

# Preprocessing.

if check_input:

X = check_array(X, accept_sparse='csr', dtype=np.float64,

accept_large_sparse=solver != 'liblinear')

y = check_array(y, ensure_2d=False, dtype=None)

check_consistent_length(X, y)

_, n_features = X.shape

classes = np.unique(y)

random_state = check_random_state(random_state)

multi_class = _check_multi_class(multi_class, solver, len(classes))

if pos_class is None and multi_class != 'multinomial':

if (classes.size > 2):

raise ValueError('To fit OvR, use the pos_class argument')

# np.unique(y) gives labels in sorted order.

pos_class = classes[1]

# If sample weights exist, convert them to array (support for lists)

# and check length

# Otherwise set them to 1 for all examples

sample_weight = _check_sample_weight(sample_weight, X,

dtype=X.dtype)

# If class_weights is a dict (provided by the user), the weights

# are assigned to the original labels. If it is "balanced", then

# the class_weights are assigned after masking the labels with a OvR.

le = LabelEncoder()

if isinstance(class_weight, dict) or multi_class == 'multinomial':

class_weight_ = compute_class_weight(class_weight,

classes=classes, y=y)

sample_weight *= class_weight_[le.fit_transform(y)]

# For doing a ovr, we need to mask the labels first. for the

# multinomial case this is not necessary.

if multi_class == 'ovr':

w0 = np.zeros(n_features + int(fit_intercept), dtype=X.dtype)

mask_classes = np.array([-1, 1])

mask = (y == pos_class)

y_bin = np.ones(y.shape, dtype=X.dtype)

y_bin[~mask] = -1.

# for compute_class_weight

if class_weight == "balanced":

class_weight_ = compute_class_weight(class_weight,

classes=mask_classes,

y=y_bin)

sample_weight *= class_weight_[le.fit_transform(y_bin)]

else:

if solver not in ['sag', 'saga']:

lbin = LabelBinarizer()

Y_multi = lbin.fit_transform(y)

if Y_multi.shape[1] == 1:

Y_multi = np.hstack([1 - Y_multi, Y_multi])

else:

# SAG multinomial solver needs LabelEncoder, not LabelBinarizer

le = LabelEncoder()

Y_multi = le.fit_transform(y).astype(X.dtype, copy=False)

w0 = np.zeros((classes.size, n_features + int(fit_intercept)),

order='F', dtype=X.dtype)

if coef is not None:

# it must work both giving the bias term and not

if multi_class == 'ovr':

if coef.size not in (n_features, w0.size):

raise ValueError(

'Initialization coef is of shape %d, expected shape '

'%d or %d' % (coef.size, n_features, w0.size))

w0[:coef.size] = coef

else:

# For binary problems coef.shape[0] should be 1, otherwise it

# should be classes.size.

n_classes = classes.size

if n_classes == 2:

n_classes = 1

if (coef.shape[0] != n_classes or

coef.shape[1] not in (n_features, n_features + 1)):

raise ValueError(

'Initialization coef is of shape (%d, %d), expected '

'shape (%d, %d) or (%d, %d)' % (

coef.shape[0], coef.shape[1], classes.size,

n_features, classes.size, n_features + 1))

if n_classes == 1:

w0[0, :coef.shape[1]] = -coef

w0[1, :coef.shape[1]] = coef

else:

w0[:, :coef.shape[1]] = coef

if multi_class == 'multinomial':

# scipy.optimize.minimize and newton-cg accepts only

# ravelled parameters.

if solver in ['lbfgs', 'newton-cg']:

w0 = w0.ravel()

target = Y_multi

if solver == 'lbfgs':

def func(x, *args): return _multinomial_loss_grad(x, *args)[0:2]

elif solver == 'newton-cg':

def func(x, *args): return _multinomial_loss(x, *args)[0]

def grad(x, *args): return _multinomial_loss_grad(x, *args)[1]

hess = _multinomial_grad_hess

warm_start_sag = {'coef': w0.T}

else:

target = y_bin

if solver == 'lbfgs':

func = _logistic_loss_and_grad

elif solver == 'newton-cg':

func = _logistic_loss

def grad(x, *args): return _logistic_loss_and_grad(x, *args)[1]

hess = _logistic_grad_hess

warm_start_sag = {'coef': np.expand_dims(w0, axis=1)}

coefs = list()

n_iter = np.zeros(len(Cs), dtype=np.int32)

for i, C in enumerate(Cs):

if solver == 'lbfgs':

iprint = [-1, 50, 1, 100, 101][

np.searchsorted(np.array([0, 1, 2, 3]), verbose)]

opt_res = optimize.minimize(

func, w0, method="L-BFGS-B", jac=True,

args=(X, target, 1. / C, sample_weight),

options={"iprint": iprint, "gtol": tol, "maxiter": max_iter}

)

n_iter_i = _check_optimize_result(

solver, opt_res, max_iter,

extra_warning_msg=_LOGISTIC_SOLVER_CONVERGENCE_MSG)

w0, loss = opt_res.x, opt_res.fun

elif solver == 'newton-cg':

args = (X, target, 1. / C, sample_weight)

w0, n_iter_i = _newton_cg(hess, func, grad, w0, args=args,

maxiter=max_iter, tol=tol)

elif solver == 'liblinear':

coef_, intercept_, n_iter_i, = _fit_liblinear(

X, target, C, fit_intercept, intercept_scaling, None,

penalty, dual, verbose, max_iter, tol, random_state,

sample_weight=sample_weight)

if fit_intercept:

w0 = np.concatenate([coef_.ravel(), intercept_])

else:

w0 = coef_.ravel()

elif solver in ['sag', 'saga']:

if multi_class == 'multinomial':

target = target.astype(X.dtype, copy=False)

loss = 'multinomial'

else:

loss = 'log'

# alpha is for L2-norm, beta is for L1-norm

if penalty == 'l1':

alpha = 0.

beta = 1. / C

elif penalty == 'l2':

alpha = 1. / C

beta = 0.

else: # Elastic-Net penalty

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,

verbose, random_state, False, max_squared_sum, warm_start_sag,

is_saga=(solver == 'saga'))

else:

raise ValueError("solver must be one of {'liblinear', 'lbfgs', "

"'newton-cg', 'sag'}, got '%s' instead" % solver)

if multi_class == 'multinomial':

n_classes = max(2, classes.size)

multi_w0 = np.reshape(w0, (n_classes, -1))

if n_classes == 2:

multi_w0 = multi_w0[1][np.newaxis, :]

coefs.append(multi_w0.copy())

else:

coefs.append(w0.copy())

n_iter[i] = n_iter_i

return np.array(coefs), np.array(Cs), n_iter

# 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, intercept_scaling=1.,

multi_class='auto', random_state=None,

max_squared_sum=None, sample_weight=None,

l1_ratio=None):

"""Computes scores across logistic_regression_path

Parameters

----------

X : {array-like, sparse matrix} of shape (n_samples, n_features)

Training data.

y : array-like of 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, default=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 : int or list of floats, default=10

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, default=None

A string (see model evaluation documentation) or

a scorer callable object / function with signature

``scorer(estimator, X, y)``. 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, default=False

If False, then the bias term is set to zero. Else the last

term of each coef_ gives us the intercept.

max_iter : int, default=100

Maximum number of iterations for the solver.

tol : float, default=1e-4

Tolerance for stopping criteria.

class_weight : dict or 'balanced', default=None

Weights associated with classes in the form ``{class_label: weight}``.

If not given, all classes are supposed to have weight one.

The "balanced" mode uses the values of y to automatically adjust

weights inversely proportional to class frequencies in the input data

as ``n_samples / (n_classes * np.bincount(y))``

Note that these weights will be multiplied with sample_weight (passed

through the fit method) if sample_weight is specified.

verbose : int, default=0

For the liblinear and lbfgs solvers set verbose to any positive

number for verbosity.

solver : {'lbfgs', 'newton-cg', 'liblinear', 'sag', 'saga'}, \

default='lbfgs'

Decides which solver to use.

penalty : {'l1', 'l2', 'elasticnet'}, default='l2'

Used to specify the norm used in the penalization. The 'newton-cg',

'sag' and 'lbfgs' solvers support only l2 penalties. 'elasticnet' is

only supported by the 'saga' solver.

dual : bool, default=False

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.

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.

multi_class : {'auto', 'ovr', 'multinomial'}, default='auto'

If the option chosen is 'ovr', then a binary problem is fit for each

label. For 'multinomial' the loss minimised is the multinomial loss fit

across the entire probability distribution, *even when the data is

binary*. 'multinomial' is unavailable when solver='liblinear'.

random_state : int, RandomState instance, default=None

Used when ``solver`` == 'sag', 'saga' or 'liblinear' to shuffle the

data. See :term:`Glossary <random_state>` for details.

max_squared_sum : float, default=None

Maximum squared sum of X over samples. Used only in SAG solver.

If None, it will be computed, going through all the samples.

The value should be precomputed to speed up cross validation.

sample_weight : array-like of shape(n_samples,), default=None

Array of weights that are assigned to individual samples.

If not provided, then each sample is given unit weight.

l1_ratio : float, default=None

The Elastic-Net mixing parameter, with ``0 <= l1_ratio <= 1``. Only

used if ``penalty='elasticnet'``. 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 of 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 of shape (n_cs,)

Scores obtained for each Cs.

n_iter : ndarray of shape(n_cs,)

Actual number of iteration for each Cs.

"""

X_train = X[train]

X_test = X[test]

y_train = y[train]

y_test = y[test]

if sample_weight is not None:

sample_weight = _check_sample_weight(sample_weight, X)

sample_weight = sample_weight[train]

coefs, Cs, n_iter = _logistic_regression_path(

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)

# The score method of Logistic Regression has a classes_ attribute.

if multi_class == 'ovr':

log_reg.classes_ = np.array([-1, 1])

elif multi_class == 'multinomial':

log_reg.classes_ = np.unique(y_train)

else:

raise ValueError("multi_class should be either multinomial or ovr, "

"got %d" % multi_class)

if pos_class is not None:

mask = (y_test == pos_class)

y_test = np.ones(y_test.shape, dtype=np.float64)

y_test[~mask] = -1.

scores = list()

scoring = get_scorer(scoring)

for w in coefs:

if multi_class == 'ovr':

w = w[np.newaxis, :]

if fit_intercept:

log_reg.coef_ = w[:, :-1]

log_reg.intercept_ = w[:, -1]

else:

log_reg.coef_ = w