From 6e4f02c4273b2d4f11fec5e3de460ffe45b1c82e Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Fri, 26 Jul 2013 11:24:16 +0200 Subject: [PATCH 01/51] Implementation of logistic_regression_path. This is mostly the code from Gael's logistic_cv branch. My work is only some minor cleanup some added tests and benchmarking. --- sklearn/linear_model/__init__.py | 5 +- sklearn/linear_model/logistic.py | 450 +++++++++++++++++++- sklearn/linear_model/tests/test_logistic.py | 118 +++++ sklearn/utils/optimize.py | 103 +++++ 4 files changed, 673 insertions(+), 3 deletions(-) create mode 100644 sklearn/utils/optimize.py 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..0b7e5a994c4 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -2,11 +2,327 @@ # Alexandre Gramfort # License: 3-clause BSD -import numpy as np +""" +Logistic Regression +""" -from .base import LinearClassifierMixin, SparseCoefMixin +# Author: Gael Varoquaux +# Fabian Pedregosa +# Alexandre Gramfort + +import numbers +>>>>>>> Implementation of logistic_regression_path. +import numpy as np +from scipy import optimize, sparse + +from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator from ..feature_selection.from_model import _LearntSelectorMixin from ..svm.base import BaseLibLinear +from ..utils import as_float_array +from ..cross_validation import check_cv +from ..utils.optimize import newton_cg + + +# .. some helper functions for logistic_regression_path .. +def _phi(t, copy=True): + # helper function: return 1. / (1 + np.exp(-t)) + if copy: + t = np.copy(t) + t *= -1. + t = np.exp(t, t) + t += 1 + t = np.reciprocal(t, t) + return t + + +def _logistic_loss_and_grad(w, X, y, alpha): + # the logistic loss and its gradient + z = X.dot(w) + yz = y * z + out = np.empty(yz.shape, yz.dtype) + idx = yz > 0 + out[idx] = np.log(1 + np.exp(-yz[idx])) + out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) + out = out.sum() + .5 * alpha * w.dot(w) + + z = _phi(yz, copy=False) + z0 = (z - 1) * y + grad = X.T.dot(z0) + alpha * w + return out, grad + + +def _logistic_loss(w, X, y, alpha): + # the logistic loss and + z = X.dot(w) + yz = y * z + out = np.empty(yz.shape, yz.dtype) + idx = yz > 0 + out[idx] = np.log(1 + np.exp(-yz[idx])) + out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) + out = out.sum() + .5 * alpha * w.dot(w) + #print 'Loss %r' % out + return out + + +def _logistic_loss_grad_hess(w, X, y, alpha): + # the logistic loss, its gradient, and the matvec application of the + # Hessian + z = X.dot(w) + yz = y * z + out = np.empty(yz.shape, yz.dtype) + idx = yz > 0 + out[idx] = np.log(1 + np.exp(-yz[idx])) + out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) + out = out.sum() + .5 * alpha * w.dot(w) + + z = _phi(yz, copy=False) + z0 = (z - 1) * y + grad = X.T.dot(z0) + alpha * w + + # The mat-vec product of the Hessian + d = z * (1 - z) + if sparse.issparse(X): + def Hs(s): + ret = d * X.dot(s) + return X.T.dot(ret) + alpha * s + else: + # Precompute as much as possible + d = np.sqrt(d, d) + # XXX: how to do this with sparse matrices? + dX = d[:, np.newaxis] * X + def Hs(s): + ret = dX.T.dot(dX.dot(s)) + ret += alpha * s + return ret + #print 'Loss/grad/hess %r, %r' % (out, grad.dot(grad)) + return out, grad, Hs + + +def _logistic_loss_and_grad_intercept(w_c, X, y, alpha): + w = w_c[:-1] + c = w_c[-1] + + z = X.dot(w) + z += c + yz = y * z + out = np.empty(yz.shape, yz.dtype) + idx = yz > 0 + out[idx] = np.log(1 + np.exp(-yz[idx])) + out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) + out = out.sum() + .5 * alpha * w.dot(w) + + z = _phi(yz, copy=False) + z0 = (z - 1) * y + grad = np.empty_like(w_c) + grad[:-1] = X.T.dot(z0) + alpha * w + grad[-1] = z0.sum() + return out, grad + + +def _logistic_loss_intercept(w_c, X, y, alpha): + w = w_c[:-1] + c = w_c[-1] + + z = X.dot(w) + z += c + yz = y * z + out = np.empty(yz.shape, yz.dtype) + idx = yz > 0 + out[idx] = np.log(1 + np.exp(-yz[idx])) + out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) + out = out.sum() + .5 * alpha * w.dot(w) + + #print 'Loss %r' % out + return out + + +def _logistic_loss_grad_hess_intercept(w_c, X, y, alpha): + w = w_c[:-1] + c = w_c[-1] + + z = X.dot(w) + z += c + yz = y * z + out = np.empty(yz.shape, yz.dtype) + idx = yz > 0 + out[idx] = np.log(1 + np.exp(-yz[idx])) + out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) + out = out.sum() + .5 * alpha * w.dot(w) + + z = _phi(yz, copy=False) + z0 = (z - 1) * y + grad = np.empty_like(w_c) + grad[:-1] = X.T.dot(z0) + alpha * w + z0_sum = z0.sum() + grad[-1] = z0_sum + # The mat-vec product of the Hessian + d = z * (1 - z) + if sparse.issparse(X): + def Hs(s): + ret = np.empty_like(s) + ret[:-1] = d * X.dot(s[:-1]) + ret[:-1] += X.T.dot(ret[:-1]) + alpha * s + ret[-1] = z0_sum * s[-1] + return + else: + # Precompute as much as possible + d = np.sqrt(d, d) + dX = d[:, np.newaxis] * X + #print 'Loss/grad/hess %r, %r' % (out, grad.dot(grad)) + def Hs(s): + ret = np.empty_like(s) + ret[:-1] = dX.T.dot(dX.dot(s[:-1])) + ret[:-1] += alpha * s[:-1] + # XXX: I am not sure that this last line of the Hessian is right + # Without the intercept the Hessian is right, though + ret[-1] = z0_sum * s[-1] + return ret + + return out, grad, Hs + +def logistic_regression_path(X, y, Cs=10, fit_intercept=True, + max_iter=100, gtol=1e-4, verbose=0, + method='trust-ncg', callback=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. + + 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 : array-like or integer of 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. + + fit_intercept : boolean + 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 : integer + Maximum number of iterations for the solver. + + gtol : float + Stopping criterion. The iteration will stop when + ``max{|g_i | i = 1, ..., n} <= gtol`` + where ``g_i`` is the i-th component of the gradient. Only used + by the methods 'lbfgs' and 'trust-ncg' + + verbose: int + Print convergence message if True. + + method : {'lbfgs', 'newton-cg', 'liblinear'} + Numerical solver to use. + + callback : callable + Function to be called before and after the fit of each regularization + parameter. Must have the signature callback(w, X, y, alpha). + + + Returns + ------- + coefs: array 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 seconds dimension will be + n_features + 1, where the last item represents the intercept. + + + Notes + ----- + You might get slighly different results with the solver trust-ncg than + with the others since this uses LIBLINEAR penalizes the intercept. + """ + if isinstance(Cs, numbers.Integral): + Cs = np.logspace(-4, 4, Cs) + Cs = np.sort(Cs) + y = np.sign(y - np.asarray(y).mean()) + X = as_float_array(X, copy=False) + if not (np.unique(y).size == 2): + raise NotImplementedError('logistic_regression_path is currently only ' + 'implemented for the binary class case') + if fit_intercept: + w0 = np.zeros(X.shape[1] + 1) + func = _logistic_loss_and_grad_intercept + else: + w0 = np.zeros(X.shape[1]) + func = _logistic_loss_and_grad + coefs = list() + + for C in Cs: + if callback is not None: + callback(w0, X, y, 1. / C) + if method == 'lbfgs': + out = optimize.fmin_l_bfgs_b( + func, w0, fprime=None, + args=(X, y, 1./C), + iprint=verbose > 0, pgtol=gtol, maxiter=max_iter) + w0 = out[0] + elif method == 'newton-cg': + if fit_intercept: + func_grad_hess = _logistic_loss_grad_hess_intercept + func = _logistic_loss_intercept + else: + func_grad_hess = _logistic_loss_grad_hess + func = _logistic_loss + + w0 = newton_cg(func_grad_hess, func, w0, args=(X, y, 1./C), + maxiter=max_iter) + elif method == 'liblinear': + lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=gtol) + lr.fit(X, y) + if fit_intercept: + w0 = np.concatenate([lr.coef_.ravel(), lr.intercept_]) + else: + w0 = lr.coef_.ravel() + else: + raise ValueError("method must be one of {'trust-ncg', 'lbfgs', " + "'newton-cg'}") + if callback is not None: + callback(w0, X, y, 1. / C) + coefs.append(w0) + return coefs, Cs + + +# helper function for LogisticCV +def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, + fit_intercept=False, + max_iter=100, gtol=1e-4, + tol=1e-4, verbose=0): + log_reg = LogisticRegression(fit_intercept=fit_intercept) + log_reg._enc = LabelEncoder() + log_reg._enc.fit_transform([-1, 1]) + + coefs, Cs = logistic_regression_path(X[train], y[train], Cs=Cs, + fit_intercept=fit_intercept, + solver=solver, + max_iter=max_iter, + gtol=gtol, tol=tol, verbose=verbose) + scores = list() + X_test = X[test] + y_test = y[test] + 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, :] + 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, @@ -141,3 +457,133 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, model, where classes are ordered as they are in ``self.classes_``. """ return np.log(self.predict_proba(X)) + + +class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, + _LearntSelectorMixin): + """Logistic Regression (aka logit, MaxEnt) classifier. + + This class implements L2 regularized logistic regression using and + LBFGS optimizer. + + Parameters + ---------- + Cs : list of floats, integer + + fit_intercept : bool, default: True + Specifies if a constant (a.k.a. bias or intercept) should be + added the decision function. + + tol: float, optional + Tolerance for stopping criteria. + + Attributes + ---------- + `coef_` : array, shape = [n_classes-1, n_features] + Coefficient of the features in the decision function. + + `coef_` is readonly property derived from `raw_coef_` that \ + follows the internal memory layout of liblinear. + + `intercept_` : array, shape = [n_classes-1] + Intercept (a.k.a. bias) added to the decision function. + It is available only when parameter intercept is set to True. + + See also + -------- + LogisticRegression + + """ + + + def __init__(self, Cs=10, fit_intercept=True, cv=None, scoring=None, + solver='newton', tol=1e-4, gtol=1e-4, max_iter=100, + n_jobs=1, verbose=False): + self.Cs = Cs + self.fit_intercept = fit_intercept + self.cv = cv + self.scoring = scoring + self.tol = tol + self.gtol = gtol + self.max_iter = max_iter + self.n_jobs = n_jobs + self.verbose = verbose + self.solver = solver + + 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. + """ + self._enc = LabelEncoder() + y = self._enc.fit_transform(y) + if len(self.classes_) != 2: + raise ValueError("LogisticRegressionCV works only on 2 " + "class problems. Please use " + "OneVsOneClassifier or OneVsRestClassifier") + + if X.shape[0] != y.shape[0]: + raise ValueError("X and y have incompatible shapes.\n" + "X has %s samples, but y has %s." % + (X.shape[0], y.shape[0])) + + # Transform to [-1, 1] classes, as y is [0, 1] + y = 2 * y + y -= 1 + + # init cross-validation generator + cv = check_cv(self.cv, X, y, classifier=True) + folds = list(cv) + + fold_coefs_ = Parallel(n_jobs=self.n_jobs, verbose=self.verbose)( + delayed(_log_reg_scoring_path)(X, y, train, test, + Cs=self.Cs, + fit_intercept=self.fit_intercept, + solver=self.solver, + max_iter=self.max_iter, + gtol=self.gtol, tol=self.tol, + verbose=max(0, self.verbose - 1), + scoring=self.scoring, + ) + for train, test in folds + ) + coefs_paths, Cs, scores = zip(*fold_coefs_) + self.Cs_ = Cs[0] + self.coefs_paths_ = coefs_paths + self.scores_ = np.array(scores) + best_index = self.scores_.sum(axis=0).argmax() + self.C_ = self.Cs_[best_index] + coef_init = np.mean([c[best_index] for c in coefs_paths], axis=0) + w = logistic_regression_path(X, y, C=[self.C_], + fit_intercept=self.fit_intercept, + w0=coef_init, + solver=self.solver, + max_iter=self.max_iter, + gtol=self.gtol, tol=self.tol, + verbose=max(0, self.verbose-1), + ) + w = w[0] + if self.fit_intercept: + self.coef_ = w[np.newaxis, :-1] + self.intercept_ = w[-1] + else: + self.coef_ = w[np.newaxis, :] + self.intercept_ = 0 + return self + + @property + def classes_(self): + return self._enc.classes_ + diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index f3d11ed0581..1c78ea60ed8 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -1,5 +1,6 @@ import numpy as np import scipy.sparse as sp +from scipy import linalg, optimize from sklearn.utils.testing import assert_array_equal from sklearn.utils.testing import assert_array_almost_equal @@ -148,6 +149,33 @@ def test_nan(): logistic.LogisticRegression(random_state=0).fit(Xnan, Y1) +def test_consistency_path(): + """Test that the path algorithm is consistent""" + Cs = np.logspace(0, 4, 10) + # can't test with fit_intercept=True since LIBLINEAR + # penalizes the intercept + for method in ('lbfgs', 'newton-cg', 'liblinear'): + coefs, Cs = logistic.logistic_regression_path( + X, Y1, Cs=Cs, fit_intercept=False, gtol=1e-16, method=method) + for i, C in enumerate(Cs): + lr = logistic.LogisticRegression( + C=C,fit_intercept=False, tol=1e-16) + lr.fit(X, Y1) + lr_coef = lr.coef_.ravel() + assert_array_almost_equal(lr_coef, coefs[i], decimal=1) + + # test for fit_intercept=True + for method in ('lbfgs', 'newton-cg', 'liblinear'): + Cs = [1e3] + coefs, Cs = logistic.logistic_regression_path( + X, Y1, Cs=Cs, fit_intercept=True, gtol=1e-16, method=method) + lr = logistic.LogisticRegression( + C=Cs[0], fit_intercept=True, tol=1e-16) + lr.fit(X, Y1) + lr_coef = np.concatenate([lr.coef_.ravel(), lr.intercept_]) + assert_array_almost_equal(lr_coef, coefs[0], decimal=1) + + def test_liblinear_random_state(): X, y = datasets.make_classification(n_samples=20) lr1 = logistic.LogisticRegression(random_state=0) @@ -155,3 +183,93 @@ def test_liblinear_random_state(): lr2 = logistic.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 = datasets.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._logistic_loss_and_grad(w, X, y, alpha=1.) + approx_grad = optimize.approx_fprime(w, + lambda w: logistic._logistic_loss_and_grad(w, X, y, + alpha=1.)[0], + 1e-3 + ) + np.testing.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._logistic_loss_and_grad_intercept(w, + X, y, alpha=1.) + np.testing.assert_allclose(loss, loss_interp) + + approx_grad = optimize.approx_fprime(w, + lambda w: + logistic._logistic_loss_and_grad_intercept(w, X, y, + alpha=1.)[0], + 1e-3 + ) + np.testing.assert_array_almost_equal(grad_interp, approx_grad, decimal=2) + + +def test__logistic_loss_grad_hess(): + n_samples, n_features = 100, 5 + X_ref = np.random.randn(n_samples, n_features) + y = np.sign(X_ref.dot(5 * np.random.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._logistic_loss_and_grad(w, X, y, alpha=1.) + loss_2, grad_2, hess = logistic._logistic_loss_grad_hess(w, X, y, + alpha=1.) + np.testing.assert_array_almost_equal(grad, grad_2) + # XXX: we should check a few simple properties of our problem, such + # as the fact that if X=0, the problem is alpha * ||w||**2, so we + # know the hessian + + # 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._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() + + np.testing.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._logistic_loss_and_grad_intercept(w, + X, y, alpha=1.) + loss_interp_2, grad_interp_2, hess = \ + logistic._logistic_loss_grad_hess_intercept(w, + X, y, alpha=1.) + np.testing.assert_allclose(loss_interp, loss_interp_2) + np.testing.assert_allclose(grad_interp, grad_interp_2) \ No newline at end of file diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py new file mode 100644 index 00000000000..fe3041b7735 --- /dev/null +++ b/sklearn/utils/optimize.py @@ -0,0 +1,103 @@ +""" +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 +from scipy.optimize.linesearch import line_search_BFGS + +def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, + maxiter=100, disp=False): + """ + Minimization of scalar function of one or more variables using the + Newton-CG algorithm. + + func: callable + Should return the value of the function, the gradient, and a + callable returning the matvec product of the Hessian + """ + avextol = xtol + + x0 = np.asarray(x0).flatten() + xtol = len(x0)*avextol + update = [2*xtol] + xk = x0 + k = 0 + old_fval = None + + # Outer loop: our Newton iteration + while (np.sum(np.abs(update)) > xtol) and (k < maxiter): + # Compute a search direction pk by applying the CG method to + # del2 f(xk) p = - grad f(xk) starting from 0. + fval, grad, fhess_p = func_grad_hess(xk, *args) + maggrad = np.sum(np.abs(grad)) + eta = min([0.5, np.sqrt(maggrad)]) + termcond = eta * maggrad + xsupi = np.zeros(len(x0), dtype=x0.dtype) + ri = grad + 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: + 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. + + if old_fval is None: + old_fval = fval + alphak, fc, gc, old_fval = line_search_BFGS(func, xk, xsupi, grad, + old_fval, args=args) + + update = alphak * xsupi + xk = xk + update # upcast if necessary + k += 1 + + return xk + + +############################################################################### +# Tests + +if __name__ == "__main__": + A = np.random.normal(size=(10, 10)) + + def func(x): + print 'Call to f: x %r' % x + Ax = A.dot(x) + return .5*(Ax).dot(Ax) + + def func_grad_hess(x): + print 'Call to f_g_h: x %r' % 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) \ No newline at end of file From bb74497c67e00354ba5dcf63be5d7711b7b591bb Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Mon, 5 Aug 2013 11:28:56 +0200 Subject: [PATCH 02/51] Take into account @agramfort's comments. Only cosmetic changes. --- sklearn/linear_model/logistic.py | 55 +++++++++------------ sklearn/linear_model/tests/test_logistic.py | 15 +++--- 2 files changed, 30 insertions(+), 40 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 0b7e5a994c4..6ab6da360ee 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -39,7 +39,7 @@ def _logistic_loss_and_grad(w, X, y, alpha): # the logistic loss and its gradient z = X.dot(w) yz = y * z - out = np.empty(yz.shape, yz.dtype) + out = np.empty_like(yz) idx = yz > 0 out[idx] = np.log(1 + np.exp(-yz[idx])) out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) @@ -55,7 +55,7 @@ def _logistic_loss(w, X, y, alpha): # the logistic loss and z = X.dot(w) yz = y * z - out = np.empty(yz.shape, yz.dtype) + out = np.empty_like(yz) idx = yz > 0 out[idx] = np.log(1 + np.exp(-yz[idx])) out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) @@ -69,7 +69,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): # Hessian z = X.dot(w) yz = y * z - out = np.empty(yz.shape, yz.dtype) + out = np.empty_like(yz) idx = yz > 0 out[idx] = np.log(1 + np.exp(-yz[idx])) out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) @@ -81,19 +81,17 @@ def _logistic_loss_grad_hess(w, X, y, alpha): # The mat-vec product of the Hessian d = z * (1 - z) + d = np.sqrt(d, out=d) if sparse.issparse(X): - def Hs(s): - ret = d * X.dot(s) - return X.T.dot(ret) + alpha * s + dX = sparse.dia_matrix((d, 0), shape=(d.size, d.size)).dot(X) else: # Precompute as much as possible - d = np.sqrt(d, d) - # XXX: how to do this with sparse matrices? dX = d[:, np.newaxis] * X - def Hs(s): - ret = dX.T.dot(dX.dot(s)) - ret += alpha * s - return ret + + def Hs(s): + ret = dX.T.dot(dX.dot(s)) + ret += alpha * s + return ret #print 'Loss/grad/hess %r, %r' % (out, grad.dot(grad)) return out, grad, Hs @@ -105,7 +103,7 @@ def _logistic_loss_and_grad_intercept(w_c, X, y, alpha): z = X.dot(w) z += c yz = y * z - out = np.empty(yz.shape, yz.dtype) + out = np.empty_like(yz) idx = yz > 0 out[idx] = np.log(1 + np.exp(-yz[idx])) out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) @@ -126,7 +124,7 @@ def _logistic_loss_intercept(w_c, X, y, alpha): z = X.dot(w) z += c yz = y * z - out = np.empty(yz.shape, yz.dtype) + out = np.empty_like(yz) idx = yz > 0 out[idx] = np.log(1 + np.exp(-yz[idx])) out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) @@ -143,7 +141,7 @@ def _logistic_loss_grad_hess_intercept(w_c, X, y, alpha): z = X.dot(w) z += c yz = y * z - out = np.empty(yz.shape, yz.dtype) + out = np.empty_like(yz) idx = yz > 0 out[idx] = np.log(1 + np.exp(-yz[idx])) out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) @@ -157,26 +155,20 @@ def _logistic_loss_grad_hess_intercept(w_c, X, y, alpha): grad[-1] = z0_sum # The mat-vec product of the Hessian d = z * (1 - z) + d = np.sqrt(d, out=d) if sparse.issparse(X): - def Hs(s): - ret = np.empty_like(s) - ret[:-1] = d * X.dot(s[:-1]) - ret[:-1] += X.T.dot(ret[:-1]) + alpha * s - ret[-1] = z0_sum * s[-1] - return + dX = sparse.dia_matrix((d, 0), shape=(d.size, d.size)).dot(X) else: # Precompute as much as possible - d = np.sqrt(d, d) dX = d[:, np.newaxis] * X - #print 'Loss/grad/hess %r, %r' % (out, grad.dot(grad)) - def Hs(s): - ret = np.empty_like(s) - ret[:-1] = dX.T.dot(dX.dot(s[:-1])) - ret[:-1] += alpha * s[:-1] - # XXX: I am not sure that this last line of the Hessian is right - # Without the intercept the Hessian is right, though - ret[-1] = z0_sum * s[-1] - return ret + def Hs(s): + ret = np.empty_like(s) + ret[:-1] = dX.T.dot(dX.dot(s[:-1])) + ret[:-1] += alpha * s[:-1] + # XXX: I am not sure that this last line of the Hessian is right + # Without the intercept the Hessian is right, though + ret[-1] = z0_sum * s[-1] + return ret return out, grad, Hs @@ -495,7 +487,6 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, """ - def __init__(self, Cs=10, fit_intercept=True, cv=None, scoring=None, solver='newton', tol=1e-4, gtol=1e-4, max_iter=100, n_jobs=1, verbose=False): diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 1c78ea60ed8..e0fbd295031 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -202,13 +202,13 @@ def test__logistic_loss_and_grad(): alpha=1.)[0], 1e-3 ) - np.testing.assert_array_almost_equal(grad, approx_grad, decimal=2) + 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._logistic_loss_and_grad_intercept(w, X, y, alpha=1.) - np.testing.assert_allclose(loss, loss_interp) + assert_array_almost_equal(loss, loss_interp) approx_grad = optimize.approx_fprime(w, lambda w: @@ -216,7 +216,7 @@ def test__logistic_loss_and_grad(): alpha=1.)[0], 1e-3 ) - np.testing.assert_array_almost_equal(grad_interp, approx_grad, decimal=2) + assert_array_almost_equal(grad_interp, approx_grad, decimal=2) def test__logistic_loss_grad_hess(): @@ -236,7 +236,7 @@ def test__logistic_loss_grad_hess(): loss, grad = logistic._logistic_loss_and_grad(w, X, y, alpha=1.) loss_2, grad_2, hess = logistic._logistic_loss_grad_hess(w, X, y, alpha=1.) - np.testing.assert_array_almost_equal(grad, grad_2) + assert_array_almost_equal(grad, grad_2) # XXX: we should check a few simple properties of our problem, such # as the fact that if X=0, the problem is alpha * ||w||**2, so we # know the hessian @@ -261,8 +261,7 @@ def test__logistic_loss_grad_hess(): d_grad -= d_grad.mean(axis=0) approx_hess_col = linalg.lstsq(d_x[:, np.newaxis], d_grad)[0].ravel() - np.testing.assert_array_almost_equal(approx_hess_col, hess_col, - decimal=3) + 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) @@ -271,5 +270,5 @@ def test__logistic_loss_grad_hess(): loss_interp_2, grad_interp_2, hess = \ logistic._logistic_loss_grad_hess_intercept(w, X, y, alpha=1.) - np.testing.assert_allclose(loss_interp, loss_interp_2) - np.testing.assert_allclose(grad_interp, grad_interp_2) \ No newline at end of file + assert_array_almost_equal(loss_interp, loss_interp_2) + assert_array_almost_equal(grad_interp, grad_interp_2) \ No newline at end of file From 91ed4bb376275ff403e5f4d9e84643d422e666fa Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Mon, 5 Aug 2013 11:50:26 +0200 Subject: [PATCH 03/51] Some fixes for LogisticCV object --- sklearn/linear_model/logistic.py | 23 ++++++++++++--------- sklearn/linear_model/tests/test_logistic.py | 13 +++++++++++- 2 files changed, 25 insertions(+), 11 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 6ab6da360ee..3bac93d9167 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -17,8 +17,10 @@ from scipy import optimize, sparse 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 as_float_array +from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv from ..utils.optimize import newton_cg @@ -174,7 +176,7 @@ def _logistic_loss_grad_hess_intercept(w_c, X, y, alpha): def logistic_regression_path(X, y, Cs=10, fit_intercept=True, max_iter=100, gtol=1e-4, verbose=0, - method='trust-ncg', callback=None): + method='liblinear', callback=None): """ Compute a Logistic Regression model for a list of regularization parameters. @@ -277,8 +279,8 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, else: w0 = lr.coef_.ravel() else: - raise ValueError("method must be one of {'trust-ncg', 'lbfgs', " - "'newton-cg'}") + raise ValueError("method must be one of {'liblinear', 'lbfgs', " + "'newton-cg'}, got '%s' instead" % method) if callback is not None: callback(w0, X, y, 1. / C) coefs.append(w0) @@ -289,16 +291,16 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, fit_intercept=False, max_iter=100, gtol=1e-4, - tol=1e-4, verbose=0): + tol=1e-4, verbose=0, method='liblinear'): log_reg = LogisticRegression(fit_intercept=fit_intercept) log_reg._enc = LabelEncoder() log_reg._enc.fit_transform([-1, 1]) coefs, Cs = logistic_regression_path(X[train], y[train], Cs=Cs, fit_intercept=fit_intercept, - solver=solver, + method=method, max_iter=max_iter, - gtol=gtol, tol=tol, verbose=verbose) + gtol=gtol, verbose=verbose) scores = list() X_test = X[test] y_test = y[test] @@ -488,7 +490,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, """ def __init__(self, Cs=10, fit_intercept=True, cv=None, scoring=None, - solver='newton', tol=1e-4, gtol=1e-4, max_iter=100, + solver='newton-cg', tol=1e-4, gtol=1e-4, max_iter=100, n_jobs=1, verbose=False): self.Cs = Cs self.fit_intercept = fit_intercept @@ -519,6 +521,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Returns self. """ self._enc = LabelEncoder() + X = as_float_array(X, copy=False) y = self._enc.fit_transform(y) if len(self.classes_) != 2: raise ValueError("LogisticRegressionCV works only on 2 " @@ -542,7 +545,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, delayed(_log_reg_scoring_path)(X, y, train, test, Cs=self.Cs, fit_intercept=self.fit_intercept, - solver=self.solver, + method=self.solver, max_iter=self.max_iter, gtol=self.gtol, tol=self.tol, verbose=max(0, self.verbose - 1), @@ -560,9 +563,9 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, w = logistic_regression_path(X, y, C=[self.C_], fit_intercept=self.fit_intercept, w0=coef_init, - solver=self.solver, + method=self.solver, max_iter=self.max_iter, - gtol=self.gtol, tol=self.tol, + gtol=self.gtol, verbose=max(0, self.verbose-1), ) w = w[0] diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index e0fbd295031..db683dab86c 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -271,4 +271,15 @@ def test__logistic_loss_grad_hess(): logistic._logistic_loss_grad_hess_intercept(w, X, y, alpha=1.) assert_array_almost_equal(loss_interp, loss_interp_2) - assert_array_almost_equal(grad_interp, grad_interp_2) \ No newline at end of file + assert_array_almost_equal(grad_interp, grad_interp_2) + +def test_logistic_cv(): + # test for LogisticRegressionCV object + n_samples, n_features = 100, 5 + X_ref = np.random.randn(n_samples, n_features) + y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) + X_ref -= X_ref.mean() + X_ref /= X_ref.std() + lr_cv = logistic.LogisticRegressionCV() + lr_cv.fit(X_ref, y) + # TODO: do something From 4d2e962fbe5b1062f164a2fa9ae8ff901426a413 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 09:34:02 +0200 Subject: [PATCH 04/51] Docstring of LogisticRegressionCV --- sklearn/linear_model/logistic.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 3bac93d9167..535a1ccd432 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -471,6 +471,9 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, tol: float, optional Tolerance for stopping criteria. + scoring: callabale + Scoring function to use as cross-validation criteria. + Attributes ---------- `coef_` : array, shape = [n_classes-1, n_features] From 9139307b42b970a7b04df728e79789a90df355dc Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 09:49:38 +0200 Subject: [PATCH 05/51] Docstring --- sklearn/linear_model/logistic.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 535a1ccd432..05f85fad3c7 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -463,6 +463,10 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Parameters ---------- Cs : list of floats, integer + Each of the values in Cs describes the inverse of regularization + strength and must be a positive float. + 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 From 4cb2d9348f73a8bfa55823ffaa60327196686692 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 10:12:40 +0200 Subject: [PATCH 06/51] Refactor and some bug fixing --- sklearn/linear_model/logistic.py | 49 ++++++++++++--------- sklearn/linear_model/tests/test_logistic.py | 4 +- 2 files changed, 31 insertions(+), 22 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 05f85fad3c7..b4ace77b285 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -23,6 +23,7 @@ from ..utils import as_float_array from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv from ..utils.optimize import newton_cg +from ..externals import six # .. some helper functions for logistic_regression_path .. @@ -176,7 +177,8 @@ def _logistic_loss_grad_hess_intercept(w_c, X, y, alpha): def logistic_regression_path(X, y, Cs=10, fit_intercept=True, max_iter=100, gtol=1e-4, verbose=0, - method='liblinear', callback=None): + solver='liblinear', callback=None, + coef=None): """ Compute a Logistic Regression model for a list of regularization parameters. @@ -215,13 +217,15 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, verbose: int Print convergence message if True. - method : {'lbfgs', 'newton-cg', 'liblinear'} + solver : {'lbfgs', 'newton-cg', 'liblinear'} Numerical solver to use. callback : callable Function to be called before and after the fit of each regularization parameter. Must have the signature callback(w, X, y, alpha). + coef: array-lime, shape (n_features,) + Initialization value for coefficients of logistic regression. Returns ------- @@ -250,18 +254,24 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, else: w0 = np.zeros(X.shape[1]) func = _logistic_loss_and_grad + + 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 callback is not None: callback(w0, X, y, 1. / C) - if method == 'lbfgs': + if solver == 'lbfgs': out = optimize.fmin_l_bfgs_b( func, w0, fprime=None, - args=(X, y, 1./C), + args=(X, y, 1. / C), iprint=verbose > 0, pgtol=gtol, maxiter=max_iter) w0 = out[0] - elif method == 'newton-cg': + elif solver == 'newton-cg': if fit_intercept: func_grad_hess = _logistic_loss_grad_hess_intercept func = _logistic_loss_intercept @@ -271,7 +281,7 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, w0 = newton_cg(func_grad_hess, func, w0, args=(X, y, 1./C), maxiter=max_iter) - elif method == 'liblinear': + elif solver == 'liblinear': lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=gtol) lr.fit(X, y) if fit_intercept: @@ -279,8 +289,8 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, else: w0 = lr.coef_.ravel() else: - raise ValueError("method must be one of {'liblinear', 'lbfgs', " - "'newton-cg'}, got '%s' instead" % method) + raise ValueError("solver must be one of {'liblinear', 'lbfgs', " + "'newton-cg'}, got '%s' instead" % solver) if callback is not None: callback(w0, X, y, 1. / C) coefs.append(w0) @@ -298,7 +308,7 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, coefs, Cs = logistic_regression_path(X[train], y[train], Cs=Cs, fit_intercept=fit_intercept, - method=method, + solver=method, max_iter=max_iter, gtol=gtol, verbose=verbose) scores = list() @@ -478,6 +488,9 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, scoring: callabale Scoring function to use as cross-validation criteria. + solver: {'newton-cg', 'lbfgs', 'liblinear'} + Algorithm to use in the optimization problem. + Attributes ---------- `coef_` : array, shape = [n_classes-1, n_features] @@ -541,7 +554,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, (X.shape[0], y.shape[0])) # Transform to [-1, 1] classes, as y is [0, 1] - y = 2 * y + y *= 2 y -= 1 # init cross-validation generator @@ -567,20 +580,16 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, best_index = self.scores_.sum(axis=0).argmax() self.C_ = self.Cs_[best_index] coef_init = np.mean([c[best_index] for c in coefs_paths], axis=0) - w = logistic_regression_path(X, y, C=[self.C_], - fit_intercept=self.fit_intercept, - w0=coef_init, - method=self.solver, - max_iter=self.max_iter, - gtol=self.gtol, - verbose=max(0, self.verbose-1), - ) + w = logistic_regression_path( + X, y, Cs=[self.C_], fit_intercept=self.fit_intercept, + coef=coef_init, solver=self.solver, max_iter=self.max_iter, + gtol=self.gtol, verbose=max(0, self.verbose-1)) w = w[0] if self.fit_intercept: - self.coef_ = w[np.newaxis, :-1] + self.coef_ = w[:-1] self.intercept_ = w[-1] else: - self.coef_ = w[np.newaxis, :] + self.coef_ = w self.intercept_ = 0 return self diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index db683dab86c..c84a1642544 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -156,7 +156,7 @@ def test_consistency_path(): # penalizes the intercept for method in ('lbfgs', 'newton-cg', 'liblinear'): coefs, Cs = logistic.logistic_regression_path( - X, Y1, Cs=Cs, fit_intercept=False, gtol=1e-16, method=method) + X, Y1, Cs=Cs, fit_intercept=False, gtol=1e-16, solver=method) for i, C in enumerate(Cs): lr = logistic.LogisticRegression( C=C,fit_intercept=False, tol=1e-16) @@ -168,7 +168,7 @@ def test_consistency_path(): for method in ('lbfgs', 'newton-cg', 'liblinear'): Cs = [1e3] coefs, Cs = logistic.logistic_regression_path( - X, Y1, Cs=Cs, fit_intercept=True, gtol=1e-16, method=method) + X, Y1, Cs=Cs, fit_intercept=True, gtol=1e-16, solver=method) lr = logistic.LogisticRegression( C=Cs[0], fit_intercept=True, tol=1e-16) lr.fit(X, Y1) From 582d7827dfd64366331cebf6956e70c705e5711b Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 10:14:51 +0200 Subject: [PATCH 07/51] Docstring --- sklearn/linear_model/logistic.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index b4ace77b285..76c8fe35ccb 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -472,16 +472,19 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Parameters ---------- - Cs : list of floats, integer + Cs: list of floats, integer Each of the values in Cs describes the inverse of regularization strength and must be a positive float. Like in support vector machines, smaller values specify stronger regularization. - fit_intercept : bool, default: True + fit_intercept: bool, default: True Specifies if a constant (a.k.a. bias or intercept) should be added the decision function. + max_iter: integer, optional + Maximum number of iterations of the optimization algorithm. + tol: float, optional Tolerance for stopping criteria. @@ -583,7 +586,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, w = logistic_regression_path( X, y, Cs=[self.C_], fit_intercept=self.fit_intercept, coef=coef_init, solver=self.solver, max_iter=self.max_iter, - gtol=self.gtol, verbose=max(0, self.verbose-1)) + gtol=self.gtol, verbose=max(0, self.verbose - 1)) w = w[0] if self.fit_intercept: self.coef_ = w[:-1] From 8cdbc6f3565c29c99ba5196b4094236db4e4f740 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 11:08:50 +0200 Subject: [PATCH 08/51] FIX missing import --- sklearn/linear_model/logistic.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 76c8fe35ccb..360e0635225 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -24,6 +24,7 @@ from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv from ..utils.optimize import newton_cg from ..externals import six +from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. From 0171d80ffb31a7257e120ffd49598b0aca92bd25 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 11:38:02 +0200 Subject: [PATCH 09/51] FIX: bug in LogisticRegressionCV --- sklearn/linear_model/logistic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 360e0635225..b0e2bced6e5 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -588,7 +588,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, X, y, Cs=[self.C_], fit_intercept=self.fit_intercept, coef=coef_init, solver=self.solver, max_iter=self.max_iter, gtol=self.gtol, verbose=max(0, self.verbose - 1)) - w = w[0] + w = w[0][0] if self.fit_intercept: self.coef_ = w[:-1] self.intercept_ = w[-1] From 668e0d15554dbd03f8f9c7b40b43009ba83113cb Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 11:39:31 +0200 Subject: [PATCH 10/51] Make tests deterministic --- sklearn/linear_model/tests/test_logistic.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index c84a1642544..0063739d8cf 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -276,6 +276,7 @@ def test__logistic_loss_grad_hess(): def test_logistic_cv(): # test for LogisticRegressionCV object n_samples, n_features = 100, 5 + np.random.seed(0) X_ref = np.random.randn(n_samples, n_features) y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) X_ref -= X_ref.mean() From af6c7ed5e5b89d7f021dc5850ed97b17ed592d19 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 11:53:54 +0200 Subject: [PATCH 11/51] FIX: coef.shape --- sklearn/linear_model/logistic.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index b0e2bced6e5..a4cc9b38a64 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -588,10 +588,10 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, X, y, Cs=[self.C_], fit_intercept=self.fit_intercept, coef=coef_init, solver=self.solver, max_iter=self.max_iter, gtol=self.gtol, verbose=max(0, self.verbose - 1)) - w = w[0][0] + w = w[0][0][:, np.newaxis].T if self.fit_intercept: - self.coef_ = w[:-1] - self.intercept_ = w[-1] + self.coef_ = w[:, :-1] + self.intercept_ = w[:, -1] else: self.coef_ = w self.intercept_ = 0 From 9fc0618c50d6e3f0b927f3f831dc88f1120bed3e Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 12:07:13 +0200 Subject: [PATCH 12/51] Just to be sure --- sklearn/linear_model/logistic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index a4cc9b38a64..7a12d794b4c 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -594,7 +594,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, self.intercept_ = w[:, -1] else: self.coef_ = w - self.intercept_ = 0 + self.intercept_ = 0. return self @property From 48a862886eee3d20628b96c2e9cf77a22a2d3428 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 13:47:28 +0200 Subject: [PATCH 13/51] Add test --- sklearn/linear_model/tests/test_logistic.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 0063739d8cf..603071f80e7 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -281,6 +281,8 @@ def test_logistic_cv(): y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) X_ref -= X_ref.mean() X_ref /= X_ref.std() - lr_cv = logistic.LogisticRegressionCV() + lr_cv = logistic.LogisticRegressionCV(Cs=[.1]) lr_cv.fit(X_ref, y) - # TODO: do something + lr = logistic.LogisticRegression(C=.1) + lr.fit(X_ref, y) + assert_array_almost_equal(lr.coef_, lr_cv.coef_, decimal=2) From 8cfa8ad6ccb47d4e0c46b924e23a7b519866098c Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Tue, 17 Sep 2013 17:59:20 +0200 Subject: [PATCH 14/51] BUG when fit_intercept=False --- sklearn/linear_model/logistic.py | 1 + sklearn/linear_model/tests/test_logistic.py | 4 ++-- sklearn/utils/optimize.py | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 7a12d794b4c..d3a1ac0ef37 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -323,6 +323,7 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, 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: diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 603071f80e7..44c289a50fc 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -281,8 +281,8 @@ def test_logistic_cv(): y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) X_ref -= X_ref.mean() X_ref /= X_ref.std() - lr_cv = logistic.LogisticRegressionCV(Cs=[.1]) + lr_cv = logistic.LogisticRegressionCV(Cs=[1.], fit_intercept=False) lr_cv.fit(X_ref, y) - lr = logistic.LogisticRegression(C=.1) + lr = logistic.LogisticRegression(C=1, fit_intercept=False) lr.fit(X_ref, y) assert_array_almost_equal(lr.coef_, lr_cv.coef_, decimal=2) diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index fe3041b7735..ac65e4d70f5 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -74,9 +74,9 @@ def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, if old_fval is None: old_fval = fval + alphak, fc, gc, old_fval = line_search_BFGS(func, xk, xsupi, grad, old_fval, args=args) - update = alphak * xsupi xk = xk + update # upcast if necessary k += 1 From 1443465c56daebc84cfb335688e4d19a6698de45 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Wed, 18 Sep 2013 10:41:32 +0200 Subject: [PATCH 15/51] Fallback for failing line search --- sklearn/utils/optimize.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index ac65e4d70f5..a280a7ef63f 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -14,7 +14,8 @@ significant speedups. # License: BSD import numpy as np -from scipy.optimize.linesearch import line_search_BFGS +import warnings +from scipy.optimize.linesearch import line_search_BFGS, line_search_wolfe2 def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, maxiter=100, disp=False): @@ -77,6 +78,15 @@ def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, alphak, fc, gc, old_fval = line_search_BFGS(func, xk, xsupi, grad, old_fval, args=args) + if alphak is None: + # line search failed + out = line_search_wolfe2( + lambda x: func(x, *args), + lambda x: func_grad_hess(x, *args)[1], xk, xsupi) + alphak, fc, gc = out[0], out[1], out[2] + warnings.warn( + 'Failed to find a suitable descent direction, the algorithm' + + 'will now terminate') update = alphak * xsupi xk = xk + update # upcast if necessary k += 1 From 228f5b043fd19c4ffcff9c05c77b46381a44325a Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Wed, 18 Sep 2013 10:48:18 +0200 Subject: [PATCH 16/51] Remove warning (not needed any more) --- sklearn/utils/optimize.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index a280a7ef63f..b69fcd58b37 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -84,9 +84,6 @@ def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, lambda x: func(x, *args), lambda x: func_grad_hess(x, *args)[1], xk, xsupi) alphak, fc, gc = out[0], out[1], out[2] - warnings.warn( - 'Failed to find a suitable descent direction, the algorithm' + - 'will now terminate') update = alphak * xsupi xk = xk + update # upcast if necessary k += 1 From 174205e261b5bb7ee26e47a4e01ad5438c0a9882 Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Wed, 18 Sep 2013 11:31:52 +0200 Subject: [PATCH 17/51] Compatibility for old scipy --- sklearn/linear_model/logistic.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index d3a1ac0ef37..38b24700357 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -267,10 +267,17 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, if callback is not None: callback(w0, X, y, 1. / C) if solver == 'lbfgs': - out = optimize.fmin_l_bfgs_b( - func, w0, fprime=None, - args=(X, y, 1. / C), - iprint=verbose > 0, pgtol=gtol, maxiter=max_iter) + try: + out = optimize.fmin_l_bfgs_b( + func, w0, fprime=None, + args=(X, y, 1. / C), + iprint=verbose > 0, pgtol=gtol, 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), + iprint=verbose > 0, pgtol=gtol) w0 = out[0] elif solver == 'newton-cg': if fit_intercept: From 02c4a289fa4d70ce548f3ed74e91c9cfc95eb1fc Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Wed, 18 Sep 2013 12:05:28 +0200 Subject: [PATCH 18/51] iterate some more on line search --- sklearn/linear_model/logistic.py | 5 ++- sklearn/utils/optimize.py | 67 ++++++++++++++++++++++++-------- 2 files changed, 53 insertions(+), 19 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 38b24700357..3320c38fea7 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -67,7 +67,6 @@ def _logistic_loss(w, X, y, alpha): #print 'Loss %r' % out return out - def _logistic_loss_grad_hess(w, X, y, alpha): # the logistic loss, its gradient, and the matvec application of the # Hessian @@ -283,11 +282,13 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, if fit_intercept: func_grad_hess = _logistic_loss_grad_hess_intercept func = _logistic_loss_intercept + grad = lambda x, *args: _logistic_loss_and_grad_intercept(x, *args)[1] else: func_grad_hess = _logistic_loss_grad_hess func = _logistic_loss + grad = lambda x, *args: _logistic_loss_and_grad(x, *args)[1] - w0 = newton_cg(func_grad_hess, func, w0, args=(X, y, 1./C), + w0 = newton_cg(func_grad_hess, func, grad, w0, args=(X, y, 1./C), maxiter=max_iter) elif solver == 'liblinear': lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=gtol) diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index b69fcd58b37..a6d639c5ae4 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -15,9 +15,39 @@ significant speedups. import numpy as np import warnings -from scipy.optimize.linesearch import line_search_BFGS, line_search_wolfe2 +from scipy.optimize.linesearch import line_search_wolfe2, line_search_wolfe1 -def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, +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=(), xtol=1e-5, eps=1e-4, maxiter=100, disp=False): """ Minimization of scalar function of one or more variables using the @@ -30,22 +60,23 @@ def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, avextol = xtol x0 = np.asarray(x0).flatten() - xtol = len(x0)*avextol - update = [2*xtol] + xtol = len(x0) * avextol + update = [2 * xtol] xk = x0 k = 0 - old_fval = None + old_fval = func(x0, *args) + old_old_fval = None # Outer loop: our Newton iteration while (np.sum(np.abs(update)) > xtol) and (k < maxiter): # Compute a search direction pk by applying the CG method to - # del2 f(xk) p = - grad f(xk) starting from 0. - fval, grad, fhess_p = func_grad_hess(xk, *args) - maggrad = np.sum(np.abs(grad)) + # del2 f(xk) p = - fgrad f(xk) starting from 0. + fval, fgrad, fhess_p = func_grad_hess(xk, *args) + maggrad = np.sum(np.abs(fgrad)) eta = min([0.5, np.sqrt(maggrad)]) termcond = eta * maggrad xsupi = np.zeros(len(x0), dtype=x0.dtype) - ri = grad + ri = fgrad psupi = -ri i = 0 dri0 = np.dot(ri, ri) @@ -62,6 +93,7 @@ def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, if (i > 0): break else: + # fall back to steepest descent direction xsupi = xsupi + dri0 / curv * psupi break alphai = dri0 / curv @@ -76,14 +108,15 @@ def newton_cg(func_grad_hess, func, x0, args=(), xtol=1e-5, eps=1e-4, if old_fval is None: old_fval = fval - alphak, fc, gc, old_fval = line_search_BFGS(func, xk, xsupi, grad, - old_fval, args=args) - if alphak is None: - # line search failed - out = line_search_wolfe2( - lambda x: func(x, *args), - lambda x: func_grad_hess(x, *args)[1], xk, xsupi) - alphak, fc, gc = out[0], out[1], out[2] + + 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 From 999c23285900c7452f783bfa7c81a8b44cf443ef Mon Sep 17 00:00:00 2001 From: Fabian Pedregosa Date: Wed, 18 Sep 2013 12:07:58 +0200 Subject: [PATCH 19/51] cosmetic --- sklearn/linear_model/logistic.py | 2 +- sklearn/utils/optimize.py | 4 ---- 2 files changed, 1 insertion(+), 5 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 3320c38fea7..453763035c9 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -11,7 +11,7 @@ Logistic Regression # Alexandre Gramfort import numbers ->>>>>>> Implementation of logistic_regression_path. + import numpy as np from scipy import optimize, sparse diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index a6d639c5ae4..acabcad3639 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -105,10 +105,6 @@ def newton_cg(func_grad_hess, func, grad, x0, args=(), xtol=1e-5, eps=1e-4, i = i + 1 dri0 = dri1 # update np.dot(ri,ri) for next time. - if old_fval is None: - old_fval = fval - - try: alphak, fc, gc, old_fval, old_old_fval, gfkp1 = \ _line_search_wolfe12(func, grad, xk, xsupi, fgrad, From 12adea826bbfe29a97d29b66b832f89bc38b33c7 Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Sun, 16 Feb 2014 11:26:04 +0530 Subject: [PATCH 20/51] COSMIT --- sklearn/linear_model/tests/test_logistic.py | 97 ++++++++++----------- 1 file changed, 45 insertions(+), 52 deletions(-) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 44c289a50fc..540b7323822 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -10,14 +10,17 @@ from sklearn.utils.testing import assert_raises from sklearn.utils.testing import assert_true from sklearn.utils.testing import raises -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_and_grad_intercept, + _logistic_loss_grad_hess, _logistic_loss_grad_hess_intercept) +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): @@ -43,28 +46,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(): @@ -72,7 +73,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) @@ -89,7 +90,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) @@ -115,7 +116,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] @@ -131,7 +132,7 @@ def test_write_parameters(): #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 @@ -146,7 +147,7 @@ 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(): @@ -155,11 +156,10 @@ def test_consistency_path(): # can't test with fit_intercept=True since LIBLINEAR # penalizes the intercept for method in ('lbfgs', 'newton-cg', 'liblinear'): - coefs, Cs = logistic.logistic_regression_path( + coefs, Cs = logistic_regression_path( X, Y1, Cs=Cs, fit_intercept=False, gtol=1e-16, solver=method) for i, C in enumerate(Cs): - lr = logistic.LogisticRegression( - C=C,fit_intercept=False, tol=1e-16) + lr = LogisticRegression(C=C, fit_intercept=False, tol=1e-16) lr.fit(X, Y1) lr_coef = lr.coef_.ravel() assert_array_almost_equal(lr_coef, coefs[i], decimal=1) @@ -167,26 +167,25 @@ def test_consistency_path(): # test for fit_intercept=True for method in ('lbfgs', 'newton-cg', 'liblinear'): Cs = [1e3] - coefs, Cs = logistic.logistic_regression_path( + coefs, Cs = logistic_regression_path( X, Y1, Cs=Cs, fit_intercept=True, gtol=1e-16, solver=method) - lr = logistic.LogisticRegression( - C=Cs[0], fit_intercept=True, tol=1e-16) + lr = LogisticRegression(C=Cs[0], fit_intercept=True, tol=1e-16) lr.fit(X, Y1) lr_coef = np.concatenate([lr.coef_.ravel(), lr.intercept_]) assert_array_almost_equal(lr_coef, coefs[0], decimal=1) 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 = datasets.make_classification(n_samples=20) + X_ref, y = make_classification(n_samples=20) n_features = X_ref.shape[1] X_sp = X_ref.copy() @@ -196,26 +195,21 @@ def test__logistic_loss_and_grad(): w = np.zeros(n_features) # First check that our derivation of the grad is correct - loss, grad = logistic._logistic_loss_and_grad(w, X, y, alpha=1.) + loss, grad = _logistic_loss_and_grad(w, X, y, alpha=1.) approx_grad = optimize.approx_fprime(w, - lambda w: logistic._logistic_loss_and_grad(w, X, y, - alpha=1.)[0], - 1e-3 - ) + 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._logistic_loss_and_grad_intercept(w, + loss_interp, grad_interp = _logistic_loss_and_grad_intercept(w, X, y, alpha=1.) assert_array_almost_equal(loss, loss_interp) approx_grad = optimize.approx_fprime(w, - lambda w: - logistic._logistic_loss_and_grad_intercept(w, X, y, - alpha=1.)[0], - 1e-3 - ) + lambda w: _logistic_loss_and_grad_intercept( + w, X, y, alpha=1.)[0], 1e-3) assert_array_almost_equal(grad_interp, approx_grad, decimal=2) @@ -233,9 +227,8 @@ def test__logistic_loss_grad_hess(): # First check that _logistic_loss_grad_hess is consistent # with _logistic_loss_and_grad - loss, grad = logistic._logistic_loss_and_grad(w, X, y, alpha=1.) - loss_2, grad_2, hess = logistic._logistic_loss_grad_hess(w, X, y, - alpha=1.) + 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) # XXX: we should check a few simple properties of our problem, such # as the fact that if X=0, the problem is alpha * ||w||**2, so we @@ -253,8 +246,8 @@ def test__logistic_loss_grad_hess(): e = 1e-3 d_x = np.linspace(-e, e, 30) d_grad = np.array([ - logistic._logistic_loss_and_grad( - w + t*vector, X, y, alpha=1.)[1] + _logistic_loss_and_grad( + w + t * vector, X, y, alpha=1.)[1] for t in d_x ]) @@ -265,14 +258,14 @@ def test__logistic_loss_grad_hess(): # Second check that our intercept implementation is good w = np.zeros(n_features + 1) - loss_interp, grad_interp = logistic._logistic_loss_and_grad_intercept(w, + loss_interp, grad_interp = _logistic_loss_and_grad_intercept(w, X, y, alpha=1.) loss_interp_2, grad_interp_2, hess = \ - logistic._logistic_loss_grad_hess_intercept(w, - X, y, alpha=1.) + _logistic_loss_grad_hess_intercept(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 = 100, 5 @@ -281,8 +274,8 @@ def test_logistic_cv(): y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) X_ref -= X_ref.mean() X_ref /= X_ref.std() - lr_cv = logistic.LogisticRegressionCV(Cs=[1.], fit_intercept=False) + lr_cv = LogisticRegressionCV(Cs=[1.], fit_intercept=False) lr_cv.fit(X_ref, y) - lr = logistic.LogisticRegression(C=1, fit_intercept=False) + lr = LogisticRegression(C=1, fit_intercept=False) lr.fit(X_ref, y) assert_array_almost_equal(lr.coef_, lr_cv.coef_, decimal=2) From 53e259333d6c49d1aa55252b614919d9d6e43610 Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Mon, 26 May 2014 09:45:09 +0530 Subject: [PATCH 21/51] Replaced Label Encoder with Label Binarizer --- sklearn/linear_model/logistic.py | 21 ++++++--------------- 1 file changed, 6 insertions(+), 15 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 453763035c9..3e0508e7cf3 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -17,9 +17,9 @@ from scipy import optimize, sparse from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator from ..feature_selection.from_model import _LearntSelectorMixin -from ..preprocessing import LabelEncoder +from ..preprocessing import LabelEncoder, LabelBinarizer from ..svm.base import BaseLibLinear -from ..utils import as_float_array +from ..utils import as_float_array, check_arrays from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv from ..utils.optimize import newton_cg @@ -553,23 +553,14 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, self : object Returns self. """ - self._enc = LabelEncoder() - X = as_float_array(X, copy=False) - y = self._enc.fit_transform(y) + X, y = check_arrays(X, y, copy=False) + self._lb = LabelBinarizer(neg_label=-1, pos_label=1) + y = np.squeeze(self._lb.fit_transform(y)) if len(self.classes_) != 2: raise ValueError("LogisticRegressionCV works only on 2 " "class problems. Please use " "OneVsOneClassifier or OneVsRestClassifier") - if X.shape[0] != y.shape[0]: - raise ValueError("X and y have incompatible shapes.\n" - "X has %s samples, but y has %s." % - (X.shape[0], y.shape[0])) - - # Transform to [-1, 1] classes, as y is [0, 1] - y *= 2 - y -= 1 - # init cross-validation generator cv = check_cv(self.cv, X, y, classifier=True) folds = list(cv) @@ -608,5 +599,5 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, @property def classes_(self): - return self._enc.classes_ + return self._lb.classes_ From b088c2add71e27954f97e950ec4ce1dc283ca008 Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Sat, 31 May 2014 23:16:52 +0530 Subject: [PATCH 22/51] Replaced helper function _phi by special.expit --- sklearn/linear_model/logistic.py | 25 +++++++------------------ 1 file changed, 7 insertions(+), 18 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 3e0508e7cf3..8dee42423ad 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -13,7 +13,7 @@ Logistic Regression import numbers import numpy as np -from scipy import optimize, sparse +from scipy import optimize, sparse, special from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator from ..feature_selection.from_model import _LearntSelectorMixin @@ -28,17 +28,6 @@ from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. -def _phi(t, copy=True): - # helper function: return 1. / (1 + np.exp(-t)) - if copy: - t = np.copy(t) - t *= -1. - t = np.exp(t, t) - t += 1 - t = np.reciprocal(t, t) - return t - - def _logistic_loss_and_grad(w, X, y, alpha): # the logistic loss and its gradient z = X.dot(w) @@ -49,7 +38,7 @@ def _logistic_loss_and_grad(w, X, y, alpha): out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) out = out.sum() + .5 * alpha * w.dot(w) - z = _phi(yz, copy=False) + z = special.expit(yz) z0 = (z - 1) * y grad = X.T.dot(z0) + alpha * w return out, grad @@ -78,7 +67,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) out = out.sum() + .5 * alpha * w.dot(w) - z = _phi(yz, copy=False) + z = special.expit(yz) z0 = (z - 1) * y grad = X.T.dot(z0) + alpha * w @@ -112,7 +101,7 @@ def _logistic_loss_and_grad_intercept(w_c, X, y, alpha): out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) out = out.sum() + .5 * alpha * w.dot(w) - z = _phi(yz, copy=False) + z = special.expit(yz) z0 = (z - 1) * y grad = np.empty_like(w_c) grad[:-1] = X.T.dot(z0) + alpha * w @@ -150,7 +139,7 @@ def _logistic_loss_grad_hess_intercept(w_c, X, y, alpha): out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) out = out.sum() + .5 * alpha * w.dot(w) - z = _phi(yz, copy=False) + z = special.expit(yz) z0 = (z - 1) * y grad = np.empty_like(w_c) grad[:-1] = X.T.dot(z0) + alpha * w @@ -584,11 +573,11 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, best_index = self.scores_.sum(axis=0).argmax() self.C_ = self.Cs_[best_index] coef_init = np.mean([c[best_index] for c in coefs_paths], axis=0) - w = logistic_regression_path( + w, _ = logistic_regression_path( X, y, Cs=[self.C_], fit_intercept=self.fit_intercept, coef=coef_init, solver=self.solver, max_iter=self.max_iter, gtol=self.gtol, verbose=max(0, self.verbose - 1)) - w = w[0][0][:, np.newaxis].T + w = w[0][:, np.newaxis].T if self.fit_intercept: self.coef_ = w[:, :-1] self.intercept_ = w[:, -1] From 137d40bad609cdfdb8f8d7daf731997936b54134 Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Sun, 1 Jun 2014 13:24:10 +0530 Subject: [PATCH 23/51] Do away with intercept helper functions --- sklearn/linear_model/logistic.py | 179 ++++++++------------ sklearn/linear_model/tests/test_logistic.py | 11 +- 2 files changed, 71 insertions(+), 119 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 8dee42423ad..2699e22793d 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -20,6 +20,7 @@ from ..feature_selection.from_model import _LearntSelectorMixin from ..preprocessing import LabelEncoder, LabelBinarizer from ..svm.base import BaseLibLinear from ..utils import as_float_array, check_arrays +from ..utils.extmath import log_logistic, safe_sparse_dot from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv from ..utils.optimize import newton_cg @@ -30,140 +31,101 @@ from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. def _logistic_loss_and_grad(w, X, y, alpha): # the logistic loss and its gradient - z = X.dot(w) + fit_intercept = False + c = 0 + _, n_features = X.shape + grad = np.empty_like(w) + + # the fit_intercept case + if w.size == n_features + 1: + fit_intercept = True + c = w[-1] + w = w[:-1] + + z = safe_sparse_dot(X, w) + z += c yz = y * z - out = np.empty_like(yz) - idx = yz > 0 - out[idx] = np.log(1 + np.exp(-yz[idx])) - out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) - out = out.sum() + .5 * alpha * w.dot(w) + + # Logistic loss is the negative of the log of the logistic function. + out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) z = special.expit(yz) z0 = (z - 1) * y - grad = X.T.dot(z0) + alpha * w + + grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w + if fit_intercept: + grad[-1] = z0.sum() return out, grad def _logistic_loss(w, X, y, alpha): + + # For the fit_intercept case. + c = 0 + if w.size == X.shape[1] + 1: + c = w[-1] + w = w[:-1] + # the logistic loss and - z = X.dot(w) + z = safe_sparse_dot(X, w) + z += c yz = y * z - out = np.empty_like(yz) - idx = yz > 0 - out[idx] = np.log(1 + np.exp(-yz[idx])) - out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) - out = out.sum() + .5 * alpha * w.dot(w) - #print 'Loss %r' % out + + # Logistic loss is the negative of the log of the logistic function. + out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) return out + def _logistic_loss_grad_hess(w, X, y, alpha): # the logistic loss, its gradient, and the matvec application of the # Hessian - z = X.dot(w) + + n_samples, n_features = X.shape + fit_intercept = False + c = 0 + grad = np.empty_like(w) + + if w.size == n_features + 1: + fit_intercept = True + c = w[-1] + w = w[:-1] + + z = safe_sparse_dot(X, w) + z += c yz = y * z - out = np.empty_like(yz) - idx = yz > 0 - out[idx] = np.log(1 + np.exp(-yz[idx])) - out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) - out = out.sum() + .5 * alpha * w.dot(w) + + # Logistic loss is the negative of the log of the logistic function. + out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) z = special.expit(yz) z0 = (z - 1) * y - grad = X.T.dot(z0) + alpha * w + grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w + if fit_intercept: + z0_sum = np.sum(z0) + grad[-1] = np.sum(z0) # The mat-vec product of the Hessian d = z * (1 - z) d = np.sqrt(d, out=d) if sparse.issparse(X): - dX = sparse.dia_matrix((d, 0), shape=(d.size, d.size)).dot(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 def Hs(s): - ret = dX.T.dot(dX.dot(s)) - ret += alpha * s + ret = np.empty_like(s) + ret[:n_features] = dX.T.dot(dX.dot(s[:n_features])) + ret[:n_features] += alpha * s[:n_features] + if fit_intercept: + # XXX: Is this right? + ret[-1] = z0_sum * s[-1] return ret #print 'Loss/grad/hess %r, %r' % (out, grad.dot(grad)) return out, grad, Hs -def _logistic_loss_and_grad_intercept(w_c, X, y, alpha): - w = w_c[:-1] - c = w_c[-1] - - z = X.dot(w) - z += c - yz = y * z - out = np.empty_like(yz) - idx = yz > 0 - out[idx] = np.log(1 + np.exp(-yz[idx])) - out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) - out = out.sum() + .5 * alpha * w.dot(w) - - z = special.expit(yz) - z0 = (z - 1) * y - grad = np.empty_like(w_c) - grad[:-1] = X.T.dot(z0) + alpha * w - grad[-1] = z0.sum() - return out, grad - - -def _logistic_loss_intercept(w_c, X, y, alpha): - w = w_c[:-1] - c = w_c[-1] - - z = X.dot(w) - z += c - yz = y * z - out = np.empty_like(yz) - idx = yz > 0 - out[idx] = np.log(1 + np.exp(-yz[idx])) - out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) - out = out.sum() + .5 * alpha * w.dot(w) - - #print 'Loss %r' % out - return out - - -def _logistic_loss_grad_hess_intercept(w_c, X, y, alpha): - w = w_c[:-1] - c = w_c[-1] - - z = X.dot(w) - z += c - yz = y * z - out = np.empty_like(yz) - idx = yz > 0 - out[idx] = np.log(1 + np.exp(-yz[idx])) - out[~idx] = (-yz[~idx] + np.log(1 + np.exp(yz[~idx]))) - out = out.sum() + .5 * alpha * w.dot(w) - - z = special.expit(yz) - z0 = (z - 1) * y - grad = np.empty_like(w_c) - grad[:-1] = X.T.dot(z0) + alpha * w - z0_sum = z0.sum() - grad[-1] = z0_sum - # The mat-vec product of the Hessian - d = z * (1 - z) - d = np.sqrt(d, out=d) - if sparse.issparse(X): - dX = sparse.dia_matrix((d, 0), shape=(d.size, d.size)).dot(X) - else: - # Precompute as much as possible - dX = d[:, np.newaxis] * X - def Hs(s): - ret = np.empty_like(s) - ret[:-1] = dX.T.dot(dX.dot(s[:-1])) - ret[:-1] += alpha * s[:-1] - # XXX: I am not sure that this last line of the Hessian is right - # Without the intercept the Hessian is right, though - ret[-1] = z0_sum * s[-1] - return ret - - return out, grad, Hs - def logistic_regression_path(X, y, Cs=10, fit_intercept=True, max_iter=100, gtol=1e-4, verbose=0, solver='liblinear', callback=None, @@ -239,10 +201,8 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, 'implemented for the binary class case') if fit_intercept: w0 = np.zeros(X.shape[1] + 1) - func = _logistic_loss_and_grad_intercept else: w0 = np.zeros(X.shape[1]) - func = _logistic_loss_and_grad if coef is not None: # it must work both giving the bias term and not @@ -255,6 +215,7 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, if callback is not None: callback(w0, X, y, 1. / C) if solver == 'lbfgs': + func = _logistic_loss_and_grad try: out = optimize.fmin_l_bfgs_b( func, w0, fprime=None, @@ -268,17 +229,9 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, iprint=verbose > 0, pgtol=gtol) w0 = out[0] elif solver == 'newton-cg': - if fit_intercept: - func_grad_hess = _logistic_loss_grad_hess_intercept - func = _logistic_loss_intercept - grad = lambda x, *args: _logistic_loss_and_grad_intercept(x, *args)[1] - else: - func_grad_hess = _logistic_loss_grad_hess - func = _logistic_loss - grad = lambda x, *args: _logistic_loss_and_grad(x, *args)[1] - - w0 = newton_cg(func_grad_hess, func, grad, w0, args=(X, y, 1./C), - maxiter=max_iter) + 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), maxiter=max_iter) elif solver == 'liblinear': lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=gtol) lr.fit(X, y) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 540b7323822..c6633071117 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -12,8 +12,7 @@ from sklearn.utils.testing import raises from sklearn.linear_model.logistic import (LogisticRegression, logistic_regression_path, LogisticRegressionCV, - _logistic_loss_and_grad, _logistic_loss_and_grad_intercept, - _logistic_loss_grad_hess, _logistic_loss_grad_hess_intercept) + _logistic_loss_and_grad, _logistic_loss_grad_hess) from sklearn.datasets import load_iris, make_classification X = [[-1, 0], [0, 1], [1, 1]] @@ -203,12 +202,12 @@ def test__logistic_loss_and_grad(): # Second check that our intercept implementation is good w = np.zeros(n_features + 1) - loss_interp, grad_interp = _logistic_loss_and_grad_intercept(w, + 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_intercept( + lambda w: _logistic_loss_and_grad( w, X, y, alpha=1.)[0], 1e-3) assert_array_almost_equal(grad_interp, approx_grad, decimal=2) @@ -258,10 +257,10 @@ def test__logistic_loss_grad_hess(): # Second check that our intercept implementation is good w = np.zeros(n_features + 1) - loss_interp, grad_interp = _logistic_loss_and_grad_intercept(w, + loss_interp, grad_interp = _logistic_loss_and_grad(w, X, y, alpha=1.) loss_interp_2, grad_interp_2, hess = \ - _logistic_loss_grad_hess_intercept(w, X, y, alpha=1.) + _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) From 4fe04e20c277e9b8499a965e3221ad1a704ffd2f Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Tue, 3 Jun 2014 17:26:36 +0530 Subject: [PATCH 24/51] Refactor fit_intercept case --- sklearn/linear_model/logistic.py | 55 ++++++++++++-------------------- 1 file changed, 20 insertions(+), 35 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 2699e22793d..133b12e9693 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -29,22 +29,25 @@ from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. -def _logistic_loss_and_grad(w, X, y, alpha): - # the logistic loss and its gradient - fit_intercept = False - c = 0 - _, n_features = X.shape - grad = np.empty_like(w) +def _intercept_dot(w, X, y): - # the fit_intercept case - if w.size == n_features + 1: - fit_intercept = True + c = None + if w.size == X.shape[1] + 1: c = w[-1] w = w[:-1] z = safe_sparse_dot(X, w) - z += c - yz = y * z + if c is not None: + z += c + return w, c, y*z + + +def _logistic_loss_and_grad(w, X, y, alpha): + # the logistic loss and its gradient + _, n_features = X.shape + grad = np.empty_like(w) + + w, c, yz = _intercept_dot(w, X, y) # Logistic loss is the negative of the log of the logistic function. out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) @@ -53,23 +56,14 @@ def _logistic_loss_and_grad(w, X, y, alpha): z0 = (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w - if fit_intercept: + if c is not None: grad[-1] = z0.sum() return out, grad -def _logistic_loss(w, X, y, alpha): +def _logistic_loss(w, X, y, alpha, fit_intercept=False): - # For the fit_intercept case. - c = 0 - if w.size == X.shape[1] + 1: - c = w[-1] - w = w[:-1] - - # the logistic loss and - z = safe_sparse_dot(X, w) - z += c - yz = y * z + w, c, yz = _intercept_dot(w, X, y) # Logistic loss is the negative of the log of the logistic function. out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) @@ -81,18 +75,9 @@ def _logistic_loss_grad_hess(w, X, y, alpha): # Hessian n_samples, n_features = X.shape - fit_intercept = False - c = 0 grad = np.empty_like(w) - if w.size == n_features + 1: - fit_intercept = True - c = w[-1] - w = w[:-1] - - z = safe_sparse_dot(X, w) - z += c - yz = y * z + w, c, yz = _intercept_dot(w, X, y) # Logistic loss is the negative of the log of the logistic function. out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) @@ -100,7 +85,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): z = special.expit(yz) z0 = (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w - if fit_intercept: + if c is not None: z0_sum = np.sum(z0) grad[-1] = np.sum(z0) @@ -118,7 +103,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): ret = np.empty_like(s) ret[:n_features] = dX.T.dot(dX.dot(s[:n_features])) ret[:n_features] += alpha * s[:n_features] - if fit_intercept: + if c is not None: # XXX: Is this right? ret[-1] = z0_sum * s[-1] return ret From 2e3b92f25d6affc47fcd356be1b0129ba123f3db Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Wed, 4 Jun 2014 23:35:49 +0530 Subject: [PATCH 25/51] FIX: Fixed hessian value for intercept --- sklearn/linear_model/logistic.py | 27 +++++++++++++++++++++------ 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 133b12e9693..308a88e61a3 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -30,7 +30,9 @@ from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. def _intercept_dot(w, X, y): + """ + """ c = None if w.size == X.shape[1] + 1: c = w[-1] @@ -43,7 +45,9 @@ def _intercept_dot(w, X, y): def _logistic_loss_and_grad(w, X, y, alpha): - # the logistic loss and its gradient + """ + + """ _, n_features = X.shape grad = np.empty_like(w) @@ -62,7 +66,9 @@ def _logistic_loss_and_grad(w, X, y, alpha): def _logistic_loss(w, X, y, alpha, fit_intercept=False): + """ + """ w, c, yz = _intercept_dot(w, X, y) # Logistic loss is the negative of the log of the logistic function. @@ -71,6 +77,10 @@ def _logistic_loss(w, X, y, alpha, fit_intercept=False): def _logistic_loss_grad_hess(w, X, y, alpha): + """ + + + """ # the logistic loss, its gradient, and the matvec application of the # Hessian @@ -91,7 +101,6 @@ def _logistic_loss_grad_hess(w, X, y, alpha): # The mat-vec product of the Hessian d = z * (1 - z) - d = np.sqrt(d, out=d) if sparse.issparse(X): dX = safe_sparse_dot(sparse.dia_matrix((d, 0), shape=(n_samples, n_samples)), X) @@ -99,15 +108,21 @@ def _logistic_loss_grad_hess(w, X, y, alpha): # Precompute as much as possible dX = d[:, np.newaxis] * X + if c is not None: + # Calculate the double derivative with respect to intecept. + dd_intercept = dX.sum(axis=0) + def Hs(s): ret = np.empty_like(s) - ret[:n_features] = dX.T.dot(dX.dot(s[:n_features])) + ret[:n_features] = X.T.dot(dX.dot(s[:n_features])) ret[:n_features] += alpha * s[:n_features] + if c is not None: - # XXX: Is this right? - ret[-1] = z0_sum * s[-1] + ret[:n_features] += s[-1] * dd_intercept + ret[-1] = dd_intercept.dot(s[:n_features]) + ret[-1] += z0_sum * s[-1] return ret - #print 'Loss/grad/hess %r, %r' % (out, grad.dot(grad)) + return out, grad, Hs From 200b925786d3aaf23cf8ed4a96cb5452e417aded Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Wed, 4 Jun 2014 23:59:51 +0530 Subject: [PATCH 26/51] DOC: Add docs for helper functions --- sklearn/linear_model/logistic.py | 57 ++++++++++++++++++++++++++++++-- 1 file changed, 54 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 308a88e61a3..9ab1266d0d2 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -31,7 +31,19 @@ from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. def _intercept_dot(w, X, y): """ + Computes y * np.dot(w, X), taking 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 = None if w.size == X.shape[1] + 1: @@ -46,7 +58,21 @@ def _intercept_dot(w, X, y): def _logistic_loss_and_grad(w, X, y, alpha): """ + 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 + Inverse of the cross_validation parameter. """ _, n_features = X.shape grad = np.empty_like(w) @@ -67,7 +93,21 @@ def _logistic_loss_and_grad(w, X, y, alpha): def _logistic_loss(w, X, y, alpha, fit_intercept=False): """ + 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 + Inverse of the cross_validation parameter. """ w, c, yz = _intercept_dot(w, X, y) @@ -78,12 +118,22 @@ def _logistic_loss(w, X, y, alpha, fit_intercept=False): def _logistic_loss_grad_hess(w, X, y, alpha): """ + 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 + Inverse of the cross_validation parameter. """ - # the logistic loss, its gradient, and the matvec application of the - # Hessian - n_samples, n_features = X.shape grad = np.empty_like(w) @@ -194,6 +244,7 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, if isinstance(Cs, numbers.Integral): Cs = np.logspace(-4, 4, Cs) Cs = np.sort(Cs) + y = np.sign(y - np.asarray(y).mean()) X = as_float_array(X, copy=False) if not (np.unique(y).size == 2): From ddb3395bff5c069cfcdeae960e18683a058bc1b9 Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Thu, 5 Jun 2014 11:59:31 +0530 Subject: [PATCH 27/51] ENH: LogisticRegressionCV can now handle sparse matrices --- sklearn/linear_model/logistic.py | 11 +++++++---- sklearn/linear_model/tests/test_logistic.py | 21 +++++++++++++++++++-- sklearn/utils/optimize.py | 2 -- 3 files changed, 26 insertions(+), 8 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 9ab1266d0d2..044838594c0 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -19,7 +19,7 @@ from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator from ..feature_selection.from_model import _LearntSelectorMixin from ..preprocessing import LabelEncoder, LabelBinarizer from ..svm.base import BaseLibLinear -from ..utils import as_float_array, check_arrays +from ..utils import atleast2d_or_csc, as_float_array, check_arrays from ..utils.extmath import log_logistic, safe_sparse_dot from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv @@ -159,8 +159,9 @@ def _logistic_loss_grad_hess(w, X, y, alpha): dX = d[:, np.newaxis] * X if c is not None: - # Calculate the double derivative with respect to intecept. - dd_intercept = dX.sum(axis=0) + # Calculate the double derivative with respect to intecept + # 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) @@ -246,7 +247,8 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, Cs = np.sort(Cs) y = np.sign(y - np.asarray(y).mean()) - X = as_float_array(X, copy=False) + X = atleast2d_or_csc(X, dtype=np.float64) + X, y = check_arrays(X, y, copy=False) if not (np.unique(y).size == 2): raise NotImplementedError('logistic_regression_path is currently only ' 'implemented for the binary class case') @@ -546,6 +548,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, self : object Returns self. """ + X = atleast2d_or_csc(X, dtype=np.float64) X, y = check_arrays(X, y, copy=False) self._lb = LabelBinarizer(neg_label=-1, pos_label=1) y = np.squeeze(self._lb.fit_transform(y)) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index c6633071117..e55046ade2e 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -183,7 +183,7 @@ def test_liblinear_random_state(): assert_array_almost_equal(lr1.coef_, lr2.coef_) -def test__logistic_loss_and_grad(): +def test_logistic_loss_and_grad(): X_ref, y = make_classification(n_samples=20) n_features = X_ref.shape[1] @@ -212,7 +212,7 @@ def test__logistic_loss_and_grad(): assert_array_almost_equal(grad_interp, approx_grad, decimal=2) -def test__logistic_loss_grad_hess(): +def test_logistic_loss_grad_hess(): n_samples, n_features = 100, 5 X_ref = np.random.randn(n_samples, n_features) y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) @@ -278,3 +278,20 @@ def test_logistic_cv(): lr = LogisticRegression(C=1, fit_intercept=False) lr.fit(X_ref, y) assert_array_almost_equal(lr.coef_, lr_cv.coef_, decimal=2) + + +def test_logistic_cv_sparse(): + X, y = make_classification(n_samples=100, n_features=5) + X[X < 1.0] = 0.0 + csr = sp.csr_matrix(X) + csc = sp.csc_matrix(X) + + for fit_intercept in [True, False]: + clf = LogisticRegressionCV(fit_intercept=fit_intercept) + clf.fit(X, y) + for data in [csr, csc]: + clfs = LogisticRegressionCV(fit_intercept=fit_intercept) + clfs.fit(data, y) + assert_array_almost_equal(clfs.coef_, clf.coef_) + assert_array_almost_equal(clfs.intercept_, clf.intercept_) + assert_equal(clfs.C_, clf.C_) diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index acabcad3639..0c43673c0ff 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -127,12 +127,10 @@ if __name__ == "__main__": A = np.random.normal(size=(10, 10)) def func(x): - print 'Call to f: x %r' % x Ax = A.dot(x) return .5*(Ax).dot(Ax) def func_grad_hess(x): - print 'Call to f_g_h: x %r' % x return func(x), A.T.dot(A.dot(x)), lambda x: A.T.dot(A.dot(x)) x0 = np.ones(10) From b8c6b4c5c25662a61210e455b0a92cc3cec0361d Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Thu, 5 Jun 2014 14:30:32 +0530 Subject: [PATCH 28/51] TST: Add tests to explicitly check hessian, loss and gradient for fit_intercept --- sklearn/linear_model/logistic.py | 2 +- sklearn/linear_model/tests/test_logistic.py | 36 ++++++++++++++++++++- 2 files changed, 36 insertions(+), 2 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 044838594c0..1fa2272b9d2 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -171,7 +171,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): if c is not None: ret[:n_features] += s[-1] * dd_intercept ret[-1] = dd_intercept.dot(s[:n_features]) - ret[-1] += z0_sum * s[-1] + ret[-1] += d.sum() * s[-1] return ret return out, grad, Hs diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index e55046ade2e..ebf6e425374 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -2,6 +2,7 @@ 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 @@ -281,7 +282,8 @@ def test_logistic_cv(): def test_logistic_cv_sparse(): - X, y = make_classification(n_samples=100, n_features=5) + X, y = make_classification(n_samples=100, n_features=5, + random_state=np.random.RandomState(0)) X[X < 1.0] = 0.0 csr = sp.csr_matrix(X) csc = sp.csc_matrix(X) @@ -295,3 +297,35 @@ def test_logistic_cv_sparse(): 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 = 100, 5 + X, y = make_classification(n_samples=n_samples, n_features=n_features, + random_state=np.random.RandomState(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(100)[:, 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_equal(loss_interp + 0.5 * (w[-1]**2), loss) + + # Check gradient. + assert_array_equal(grad_interp[:n_features], grad[:n_features]) + assert_equal(grad_interp[-1] + alpha * w[-1], grad[-1]) + + np.random.seed(0) + grad = np.random.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]) From a4d9f77805beec418891e54beb43e5fc6c783452 Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Fri, 6 Jun 2014 02:08:23 +0530 Subject: [PATCH 29/51] More docs and tests for LogisticRegressionCV --- sklearn/linear_model/logistic.py | 63 ++++++++++++++++++++- sklearn/linear_model/tests/test_logistic.py | 24 ++++++++ 2 files changed, 84 insertions(+), 3 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 1fa2272b9d2..b509c1a8e34 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -306,6 +306,38 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, fit_intercept=False, max_iter=100, gtol=1e-4, tol=1e-4, verbose=0, method='liblinear'): + """ + 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 values + + train : list of indices + The indices of the train set + + test : list of indices + The indices of the test set + + fit_intercept : bool + If False, then the bias term is set to zero. + + max_iter : int + Maximum no. of iterations for the solver. + + gtol : float + Stopping criteria + + verbose : int + Amount of verbosity + + method : {'lbfgs', 'newton-cg', 'liblinear'} + Decides which solver to use. + """ log_reg = LogisticRegression(fit_intercept=fit_intercept) log_reg._enc = LabelEncoder() log_reg._enc.fit_transform([-1, 1]) @@ -487,10 +519,16 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Specifies if a constant (a.k.a. bias or intercept) should be added the decision function. + 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. + max_iter: integer, optional Maximum number of iterations of the optimization algorithm. - tol: float, optional + gtol: float, optional Tolerance for stopping criteria. scoring: callabale @@ -499,18 +537,37 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, solver: {'newton-cg', 'lbfgs', 'liblinear'} Algorithm to use in the optimization problem. + verbose : bool or integer + Amount of verbosity. + + n_jobs : integer, optional + Number of CPU cores used during the cross-validation loop. If given + a value of -1, all cores are used. + Attributes ---------- - `coef_` : array, shape = [n_classes-1, n_features] + `coef_` : array, shape = (n_classes-1, n_features) Coefficient of the features in the decision function. `coef_` is readonly property derived from `raw_coef_` that \ follows the internal memory layout of liblinear. - `intercept_` : array, shape = [n_classes-1] + `intercept_` : array, shape = (n_classes-1) Intercept (a.k.a. bias) added to the decision function. It is available only when parameter intercept is set to True. + `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 + 1) or + (n_folds, len(Cs_), n_features + 1) + path of coefficients obtained during cross-validating across each + fold and then across each Cs. + + `scores_` : array, shape = [n_folds, len(Cs_)] + grid of scores obtained during cross-validating each fold. + See also -------- LogisticRegression diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index ebf6e425374..ae380888bcc 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -329,3 +329,27 @@ def test_intercept_logistic_helper(): 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_error_multitask(): + # Right now LogisticRegressionCV cannot handle multiple classes. + X, y = make_classification(n_samples=100, n_features=50, + n_informative=20, n_classes=3) + assert_raises(ValueError, LogisticRegressionCV().fit, X, y) + + +def test_shape_attributes_logregcv(): + n_samples, n_features = 100, 100 + X, y = make_classification(n_samples=n_samples, n_features=n_features, + n_informative=20) + clf = LogisticRegressionCV(cv=3) + clf.fit(X, y) + + assert_array_equal(clf.coef_.shape, (1, n_features)) + assert_array_equal(clf.classes_, [0, 1]) + assert_equal(len(clf.classes_), 2) + + coefs_paths = np.asarray(clf.coefs_paths_) + assert_array_equal(coefs_paths.shape, (3, 10, n_features + 1)) + assert_array_equal(clf.Cs_.shape, (10, )) + assert_array_equal(clf.scores_.shape, (3, 10)) From 7c935a32358839588b9509d6b71efda466ffae88 Mon Sep 17 00:00:00 2001 From: Manoj-Kumar-S Date: Sat, 7 Jun 2014 01:58:00 +0530 Subject: [PATCH 30/51] FIX: Doctests --- sklearn/linear_model/logistic.py | 136 ++++++++++---------- sklearn/linear_model/tests/test_logistic.py | 7 +- 2 files changed, 71 insertions(+), 72 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index b509c1a8e34..22c04b76aa9 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -1,7 +1,3 @@ -# Authors: Fabian Pedregosa -# Alexandre Gramfort -# License: 3-clause BSD - """ Logistic Regression """ @@ -9,11 +5,12 @@ Logistic Regression # Author: Gael Varoquaux # Fabian Pedregosa # Alexandre Gramfort +# Manoj Kumar import numbers import numpy as np -from scipy import optimize, sparse, special +from scipy import optimize, sparse from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator from ..feature_selection.from_model import _LearntSelectorMixin @@ -21,6 +18,7 @@ from ..preprocessing import LabelEncoder, LabelBinarizer from ..svm.base import BaseLibLinear from ..utils import atleast2d_or_csc, as_float_array, check_arrays from ..utils.extmath import log_logistic, safe_sparse_dot +from ..utils.fixes import expit from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv from ..utils.optimize import newton_cg @@ -30,19 +28,19 @@ from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. def _intercept_dot(w, X, y): - """ - Computes y * np.dot(w, X), taking into consideration - if the intercept should be fit or not. + """Computes y * np.dot(w, X). + + It takes into consideration if the intercept should be fit or not. Parameters ---------- - w : ndarray, shape = (n_features,) or (n_features + 1,) + 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) + y : ndarray, shape (n_samples,) Array of labels """ c = None @@ -57,18 +55,17 @@ def _intercept_dot(w, X, y): def _logistic_loss_and_grad(w, X, y, alpha): - """ - Computes the logistic loss and gradient. + """Computes the logistic loss and gradient. Parameters ---------- - w : ndarray, shape = (n_features,) or (n_features + 1,) + 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) + y : ndarray, shape (n_samples,) Array of labels. alpha : float @@ -82,7 +79,7 @@ def _logistic_loss_and_grad(w, X, y, alpha): # Logistic loss is the negative of the log of the logistic function. out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) - z = special.expit(yz) + z = expit(yz) z0 = (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w @@ -92,12 +89,11 @@ def _logistic_loss_and_grad(w, X, y, alpha): def _logistic_loss(w, X, y, alpha, fit_intercept=False): - """ - Computes the logistic loss and gradient. + """Computes the logistic loss. Parameters ---------- - w : ndarray, shape = (n_features,) or (n_features + 1,) + w : ndarray, shape (n_features,) or (n_features + 1,) Coefficient vector X : {array-like, sparse matrix}, shape (n_samples, n_features) @@ -117,12 +113,11 @@ def _logistic_loss(w, X, y, alpha, fit_intercept=False): def _logistic_loss_grad_hess(w, X, y, alpha): - """ - Computes the logistic loss, gradient and the Hessian. + """Computes the logistic loss, gradient and the Hessian. Parameters ---------- - w : ndarray, shape = (n_features,) or (n_features + 1,) + w : ndarray, shape (n_features,) or (n_features + 1,) Coefficient vector X : {array-like, sparse matrix}, shape (n_samples, n_features) @@ -142,7 +137,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): # Logistic loss is the negative of the log of the logistic function. out = -np.sum(log_logistic(yz)) + .5 * alpha * np.dot(w, w) - z = special.expit(yz) + z = expit(yz) z0 = (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w if c is not None: @@ -159,7 +154,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): dX = d[:, np.newaxis] * X if c is not None: - # Calculate the double derivative with respect to intecept + # 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))) @@ -178,11 +173,10 @@ def _logistic_loss_grad_hess(w, X, y, alpha): def logistic_regression_path(X, y, Cs=10, fit_intercept=True, - max_iter=100, gtol=1e-4, verbose=0, + max_iter=100, tol=1e-4, verbose=0, solver='liblinear', callback=None, coef=None): - """ - Compute a Logistic Regression model for a list of regularization + """Compute a Logistic Regression model for a list of regularization parameters. This is an implementation that uses the result of the previous model @@ -210,11 +204,10 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, max_iter : integer Maximum number of iterations for the solver. - gtol : float + tol : float Stopping criterion. The iteration will stop when - ``max{|g_i | i = 1, ..., n} <= gtol`` - where ``g_i`` is the i-th component of the gradient. Only used - by the methods 'lbfgs' and 'trust-ncg' + ``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. @@ -231,12 +224,11 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, Returns ------- - coefs: array of shape (n_cs, n_features) or (n_cs, n_features + 1) + 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 seconds dimension will be + fit_intercept is set to True then the second dimension will be n_features + 1, where the last item represents the intercept. - Notes ----- You might get slighly different results with the solver trust-ncg than @@ -244,7 +236,6 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, """ if isinstance(Cs, numbers.Integral): Cs = np.logspace(-4, 4, Cs) - Cs = np.sort(Cs) y = np.sign(y - np.asarray(y).mean()) X = atleast2d_or_csc(X, dtype=np.float64) @@ -273,20 +264,20 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, out = optimize.fmin_l_bfgs_b( func, w0, fprime=None, args=(X, y, 1. / C), - iprint=verbose > 0, pgtol=gtol, maxiter=max_iter) + iprint=verbose > 0, 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), - iprint=verbose > 0, pgtol=gtol) + iprint=verbose > 0, pgtol=tol) w0 = out[0] 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), maxiter=max_iter) elif solver == 'liblinear': - lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=gtol) + lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=tol) lr.fit(X, y) if fit_intercept: w0 = np.concatenate([lr.coef_.ravel(), lr.intercept_]) @@ -304,10 +295,9 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, # helper function for LogisticCV def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, fit_intercept=False, - max_iter=100, gtol=1e-4, - tol=1e-4, verbose=0, method='liblinear'): - """ - Computes scores across logistic_regression_path + max_iter=100, tol=1e-4, + verbose=0, method='liblinear'): + """Computes scores across logistic_regression_path Parameters ---------- @@ -315,7 +305,7 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, Training data. y : array-like, shape (n_samples,) or (n_samples, n_targets) - Target values + Target labels train : list of indices The indices of the train set @@ -323,14 +313,25 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, test : list of indices The indices of the test set + Cs: list of floats, integer + Each of the values in Cs describes the inverse of + regularization strength and must be a positive float. + 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. + 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 no. of iterations for the solver. - gtol : float - Stopping criteria + tol : float + Tolerance for stopping criteria. verbose : int Amount of verbosity @@ -346,7 +347,7 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, fit_intercept=fit_intercept, solver=method, max_iter=max_iter, - gtol=gtol, verbose=verbose) + tol=tol, verbose=verbose) scores = list() X_test = X[test] y_test = y[test] @@ -423,10 +424,10 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, 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. @@ -525,34 +526,36 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, See the module :mod:`sklearn.cross_validation` module for the list of possible cross-validation objects. - max_iter: integer, optional - Maximum number of iterations of the optimization algorithm. - - gtol: float, optional - Tolerance for stopping criteria. - scoring: callabale - Scoring function to use as cross-validation criteria. + 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. - verbose : bool or integer - Amount of verbosity. + tol: float, optional + Tolerance for stopping criteria. + + max_iter: integer, optional + Maximum number of iterations of the optimization algorithm. n_jobs : integer, optional Number of CPU cores used during the cross-validation loop. If given a value of -1, all cores are used. + verbose : bool or integer + Amount of verbosity. + Attributes ---------- - `coef_` : array, shape = (n_classes-1, n_features) + `coef_` : array, shape (n_classes-1, n_features) Coefficient of the features in the decision function. `coef_` is readonly property derived from `raw_coef_` that \ follows the internal memory layout of liblinear. - `intercept_` : array, shape = (n_classes-1) + `intercept_` : array, shape (n_classes-1) Intercept (a.k.a. bias) added to the decision function. It is available only when parameter intercept is set to True. @@ -560,12 +563,12 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Array of C i.e inverse of regularization parameter values used for cross-validation. - `coefs_paths_` : array, shape = (n_folds, len(Cs_), n_features + 1) or + `coefs_paths_` : array, shape (n_folds, len(Cs_), n_features + 1) or (n_folds, len(Cs_), n_features + 1) path of coefficients obtained during cross-validating across each fold and then across each Cs. - `scores_` : array, shape = [n_folds, len(Cs_)] + `scores_` : array, shape (n_folds, len(Cs_)) grid of scores obtained during cross-validating each fold. See also @@ -575,14 +578,13 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, """ def __init__(self, Cs=10, fit_intercept=True, cv=None, scoring=None, - solver='newton-cg', tol=1e-4, gtol=1e-4, max_iter=100, + solver='newton-cg', tol=1e-4, max_iter=100, n_jobs=1, verbose=False): self.Cs = Cs self.fit_intercept = fit_intercept self.cv = cv self.scoring = scoring self.tol = tol - self.gtol = gtol self.max_iter = max_iter self.n_jobs = n_jobs self.verbose = verbose @@ -593,11 +595,11 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Parameters ---------- - X : {array-like, sparse matrix}, shape = [n_samples, n_features] + 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] + y : array-like, shape (n_samples,) Target vector relative to X Returns @@ -624,7 +626,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, fit_intercept=self.fit_intercept, method=self.solver, max_iter=self.max_iter, - gtol=self.gtol, tol=self.tol, + tol=self.tol, verbose=max(0, self.verbose - 1), scoring=self.scoring, ) @@ -640,7 +642,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, w, _ = logistic_regression_path( X, y, Cs=[self.C_], fit_intercept=self.fit_intercept, coef=coef_init, solver=self.solver, max_iter=self.max_iter, - gtol=self.gtol, verbose=max(0, self.verbose - 1)) + tol=self.tol, verbose=max(0, self.verbose - 1)) w = w[0][:, np.newaxis].T if self.fit_intercept: self.coef_ = w[:, :-1] diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index ae380888bcc..38b498be2db 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -157,7 +157,7 @@ def test_consistency_path(): # penalizes the intercept for method in ('lbfgs', 'newton-cg', 'liblinear'): coefs, Cs = logistic_regression_path( - X, Y1, Cs=Cs, fit_intercept=False, gtol=1e-16, solver=method) + X, Y1, 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, Y1) @@ -168,7 +168,7 @@ def test_consistency_path(): for method in ('lbfgs', 'newton-cg', 'liblinear'): Cs = [1e3] coefs, Cs = logistic_regression_path( - X, Y1, Cs=Cs, fit_intercept=True, gtol=1e-16, solver=method) + X, Y1, Cs=Cs, fit_intercept=True, tol=1e-16, solver=method) lr = LogisticRegression(C=Cs[0], fit_intercept=True, tol=1e-16) lr.fit(X, Y1) lr_coef = np.concatenate([lr.coef_.ravel(), lr.intercept_]) @@ -230,9 +230,6 @@ def test_logistic_loss_grad_hess(): 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) - # XXX: we should check a few simple properties of our problem, such - # as the fact that if X=0, the problem is alpha * ||w||**2, so we - # know the hessian # Now check our hessian along the second direction of the grad vector = np.zeros_like(grad) From 51a1831aaa05c34dab40284d748ea35f23f82c9b Mon Sep 17 00:00:00 2001 From: MechCoder Date: Sat, 7 Jun 2014 11:19:23 +0000 Subject: [PATCH 31/51] TST: Improved tests --- sklearn/linear_model/tests/test_logistic.py | 25 ++++++++++----------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 38b498be2db..f88834e82fe 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -214,7 +214,7 @@ def test_logistic_loss_and_grad(): def test_logistic_loss_grad_hess(): - n_samples, n_features = 100, 5 + n_samples, n_features = 50, 5 X_ref = np.random.randn(n_samples, n_features) y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) X_ref -= X_ref.mean() @@ -265,7 +265,7 @@ def test_logistic_loss_grad_hess(): def test_logistic_cv(): # test for LogisticRegressionCV object - n_samples, n_features = 100, 5 + n_samples, n_features = 50, 5 np.random.seed(0) X_ref = np.random.randn(n_samples, n_features) y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) @@ -279,7 +279,7 @@ def test_logistic_cv(): def test_logistic_cv_sparse(): - X, y = make_classification(n_samples=100, n_features=5, + X, y = make_classification(n_samples=50, n_features=5, random_state=np.random.RandomState(0)) X[X < 1.0] = 0.0 csr = sp.csr_matrix(X) @@ -297,7 +297,7 @@ def test_logistic_cv_sparse(): def test_intercept_logistic_helper(): - n_samples, n_features = 100, 5 + n_samples, n_features = 10, 5 X, y = make_classification(n_samples=n_samples, n_features=n_features, random_state=np.random.RandomState(0)) @@ -309,16 +309,16 @@ def test_intercept_logistic_helper(): # 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(100)[:, np.newaxis])) + 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_equal(loss_interp + 0.5 * (w[-1]**2), loss) + assert_almost_equal(loss_interp + 0.5 * (w[-1]**2), loss) # Check gradient. - assert_array_equal(grad_interp[:n_features], grad[:n_features]) - assert_equal(grad_interp[-1] + alpha * w[-1], grad[-1]) + assert_array_almost_equal(grad_interp[:n_features], grad[:n_features]) + assert_almost_equal(grad_interp[-1] + alpha * w[-1], grad[-1]) np.random.seed(0) grad = np.random.rand(n_features + 1) @@ -330,15 +330,14 @@ def test_intercept_logistic_helper(): def test_error_multitask(): # Right now LogisticRegressionCV cannot handle multiple classes. - X, y = make_classification(n_samples=100, n_features=50, - n_informative=20, n_classes=3) + X, y = make_classification(n_samples=10, n_features=20, n_informative=10, + n_classes=3) assert_raises(ValueError, LogisticRegressionCV().fit, X, y) def test_shape_attributes_logregcv(): - n_samples, n_features = 100, 100 - X, y = make_classification(n_samples=n_samples, n_features=n_features, - n_informative=20) + n_samples, n_features = 10, 10 + X, y = make_classification(n_samples=n_samples, n_features=n_features) clf = LogisticRegressionCV(cv=3) clf.fit(X, y) From f9871cd824727ebad95bdc01e29f41329f3355a6 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Mon, 23 Jun 2014 20:23:15 +0530 Subject: [PATCH 32/51] ENH: Added one-vs-all fit in case of multi-class data --- sklearn/linear_model/logistic.py | 161 ++++++++++++++------ sklearn/linear_model/tests/test_logistic.py | 26 +++- 2 files changed, 132 insertions(+), 55 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 22c04b76aa9..ef77e9d0bbd 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -172,10 +172,10 @@ def _logistic_loss_grad_hess(w, X, y, alpha): return out, grad, Hs -def logistic_regression_path(X, y, Cs=10, fit_intercept=True, +def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, max_iter=100, tol=1e-4, verbose=0, solver='liblinear', callback=None, - coef=None): + coef=None, copy=False): """Compute a Logistic Regression model for a list of regularization parameters. @@ -197,6 +197,10 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, 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 : boolean Whether to fit an intercept for the model. In this case the shape of the returned array is (n_cs, n_features + 1). @@ -222,6 +226,12 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, coef: array-lime, shape (n_features,) Initialization value for coefficients of logistic regression. + copy: boolean + 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. + Returns ------- coefs: ndarray, shape (n_cs, n_features) or (n_cs, n_features + 1) @@ -237,12 +247,20 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, if isinstance(Cs, numbers.Integral): Cs = np.logspace(-4, 4, Cs) - y = np.sign(y - np.asarray(y).mean()) X = atleast2d_or_csc(X, dtype=np.float64) - X, y = check_arrays(X, y, copy=False) - if not (np.unique(y).size == 2): - raise NotImplementedError('logistic_regression_path is currently only ' - 'implemented for the binary class case') + X, y = check_arrays(X, y, copy=copy) + + if pos_class is None: + n_classes = np.unique(y) + if not (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] + + mask = (y == pos_class) + y[mask] = 1 + y[~mask] = -1 + if fit_intercept: w0 = np.zeros(X.shape[1] + 1) else: @@ -293,8 +311,8 @@ def logistic_regression_path(X, y, Cs=10, fit_intercept=True, # helper function for LogisticCV -def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, - fit_intercept=False, +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, verbose=0, method='liblinear'): """Computes scores across logistic_regression_path @@ -313,6 +331,10 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, 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, integer Each of the values in Cs describes the inverse of regularization strength and must be a positive float. @@ -339,18 +361,33 @@ def _log_reg_scoring_path(X, y, train, test, Cs=10, scoring=None, method : {'lbfgs', 'newton-cg', 'liblinear'} Decides which solver to use. """ + log_reg = LogisticRegression(fit_intercept=fit_intercept) log_reg._enc = LabelEncoder() log_reg._enc.fit_transform([-1, 1]) - coefs, Cs = logistic_regression_path(X[train], y[train], Cs=Cs, + X_train = X[train] + X_test = X[test] + y_train = y[train] + y_test = y[test] + + if pos_class is not None: + # In order to avoid a copy in y, mask test and train separately + mask = (y_train == pos_class) + y_train[mask] = 1 + y_train[~mask] = -1 + mask = (y_test == pos_class) + y_test[mask] = 1 + y_test[~mask] = -1 + + coefs, Cs = logistic_regression_path(X_train, y_train, Cs=Cs, fit_intercept=fit_intercept, solver=method, max_iter=max_iter, tol=tol, verbose=verbose) + scores = list() - X_test = X[test] - y_test = y[test] + if isinstance(scoring, six.string_types): scoring = SCORERS[scoring] for w in coefs: @@ -549,27 +586,38 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Attributes ---------- - `coef_` : array, shape (n_classes-1, n_features) + `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 (n_classes-1) + `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. + 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 + 1) or + `coefs_paths_` : array, shape (n_folds, len(Cs_), n_features) or (n_folds, len(Cs_), n_features + 1) - path of coefficients obtained during cross-validating across each - fold and then across each Cs. + 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_` : array, shape (n_folds, len(Cs_)) - grid of scores obtained during cross-validating each fold. + `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)) See also -------- @@ -600,7 +648,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, n_features is the number of features. y : array-like, shape (n_samples,) - Target vector relative to X + Target vector relative to X. Returns ------- @@ -609,19 +657,23 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, """ X = atleast2d_or_csc(X, dtype=np.float64) X, y = check_arrays(X, y, copy=False) - self._lb = LabelBinarizer(neg_label=-1, pos_label=1) - y = np.squeeze(self._lb.fit_transform(y)) - if len(self.classes_) != 2: - raise ValueError("LogisticRegressionCV works only on 2 " - "class problems. Please use " - "OneVsOneClassifier or OneVsRestClassifier") # init cross-validation generator cv = check_cv(self.cv, X, y, classifier=True) folds = list(cv) + self.classes_ = labels = np.unique(y) + n_classes = len(labels) + + if n_classes == 2: + # OvA in case of binary problems is as good as fitting + # the higher label + n_classes = 1 + labels = labels[1:] + 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, method=self.solver, @@ -630,29 +682,42 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, verbose=max(0, self.verbose - 1), scoring=self.scoring, ) + for label in labels for train, test in folds ) coefs_paths, Cs, scores = zip(*fold_coefs_) + self.Cs_ = Cs[0] - self.coefs_paths_ = coefs_paths - self.scores_ = np.array(scores) - best_index = self.scores_.sum(axis=0).argmax() - self.C_ = self.Cs_[best_index] - coef_init = np.mean([c[best_index] for c in coefs_paths], axis=0) - w, _ = logistic_regression_path( - X, y, Cs=[self.C_], fit_intercept=self.fit_intercept, - coef=coef_init, solver=self.solver, max_iter=self.max_iter, - tol=self.tol, verbose=max(0, self.verbose - 1)) - w = w[0][:, np.newaxis].T - if self.fit_intercept: - self.coef_ = w[:, :-1] - self.intercept_ = w[:, -1] - else: - self.coef_ = w - self.intercept_ = 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] + 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, copy=True, + verbose=max(0, self.verbose - 1)) + w = w[0] + if self.fit_intercept: + self.coef_.append(w[:-1]) + self.intercept_.append(w[-1]) + else: + self.coef_.append(w) + self.intercept_.append(0.) + self.coef_ = np.asarray(self.coef_) + self.intercept_ = np.asarray(self.intercept_) return self - - @property - def classes_(self): - return self._lb.classes_ - diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index f88834e82fe..c2bc93ce724 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -328,24 +328,36 @@ def test_intercept_logistic_helper(): assert_almost_equal(hess_interp[-1] + alpha * grad[-1], hess[-1]) -def test_error_multitask(): +def test_multitask(): # Right now LogisticRegressionCV cannot handle multiple classes. X, y = make_classification(n_samples=10, n_features=20, n_informative=10, n_classes=3) - assert_raises(ValueError, LogisticRegressionCV().fit, X, y) + clf = LogisticRegressionCV(cv=3) + clf.fit(X, y) + + assert_array_equal(clf.coef_.shape, (3, 20)) + 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, 20 + 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_shape_attributes_logregcv(): - n_samples, n_features = 10, 10 - X, y = make_classification(n_samples=n_samples, n_features=n_features) + n_samples, n_features = 10, 20 clf = LogisticRegressionCV(cv=3) + X, y = make_classification(n_samples=n_samples, n_features=n_features) clf.fit(X, y) assert_array_equal(clf.coef_.shape, (1, n_features)) assert_array_equal(clf.classes_, [0, 1]) assert_equal(len(clf.classes_), 2) - coefs_paths = np.asarray(clf.coefs_paths_) - assert_array_equal(coefs_paths.shape, (3, 10, n_features + 1)) + coefs_paths = np.asarray(list(clf.coefs_paths_.values())) + assert_array_equal(coefs_paths.shape, (1, 3, 10, n_features + 1)) assert_array_equal(clf.Cs_.shape, (10, )) - assert_array_equal(clf.scores_.shape, (3, 10)) + scores = np.asarray(list(clf.scores_.values())) + assert_array_equal(scores.shape, (1, 3, 10)) From d4b544e0246bc73d6478df9ccee4a496a456734a Mon Sep 17 00:00:00 2001 From: MechCoder Date: Tue, 24 Jun 2014 12:39:37 +0530 Subject: [PATCH 33/51] ENH: Added refit parameter --- sklearn/linear_model/logistic.py | 46 ++++++++++++++++++++++++-------- 1 file changed, 35 insertions(+), 11 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index ef77e9d0bbd..2bf0b5c9b9b 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -584,6 +584,13 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, verbose : bool or integer 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 and C that corresponds to the best score across + each fold is taken and retuned after averaging. + Attributes ---------- `coef_` : array, shape (1, n_features) or @@ -619,6 +626,9 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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. + See also -------- LogisticRegression @@ -627,7 +637,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, def __init__(self, Cs=10, fit_intercept=True, cv=None, scoring=None, solver='newton-cg', tol=1e-4, max_iter=100, - n_jobs=1, verbose=False): + n_jobs=1, verbose=False, refit=True): self.Cs = Cs self.fit_intercept = fit_intercept self.cv = cv @@ -637,6 +647,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, self.n_jobs = n_jobs self.verbose = verbose self.solver = solver + self.refit = refit def fit(self, X, y): """Fit the model according to the given training data. @@ -701,17 +712,30 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, for label in labels: scores = self.scores_[label] coefs_paths = self.coefs_paths_[label] - 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, copy=True, - verbose=max(0, self.verbose - 1)) - w = w[0] + 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, copy=True, + verbose=max(0, self.verbose - 1)) + w = w[0] + + else: + # Take the best scores across every fold and the average of all + # coefficients coressponding 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]) From 191119b74acd2f2520d2390af857c9ed02db9725 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Thu, 26 Jun 2014 14:04:31 +0530 Subject: [PATCH 34/51] TST: Tests to verify OvA behavior --- sklearn/linear_model/tests/test_logistic.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index c2bc93ce724..d634fdc3831 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -14,6 +14,7 @@ from sklearn.utils.testing import raises 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]] @@ -361,3 +362,21 @@ def test_shape_attributes_logregcv(): assert_array_equal(clf.Cs_.shape, (10, )) scores = np.asarray(list(clf.scores_.values())) assert_array_equal(scores.shape, (1, 3, 10)) + + +def test_ova_iris(): + # Test that our OvA implementation is correct using the iris dataset. + train, target = iris.data, iris.target + + # 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_) From 60f454f8ece5026c0cb72dd3ab258eb6545ac94d Mon Sep 17 00:00:00 2001 From: MechCoder Date: Thu, 3 Jul 2014 15:28:47 +0200 Subject: [PATCH 35/51] FIX: PEP8 and other cosmits --- sklearn/linear_model/logistic.py | 80 +++++++++------------ sklearn/linear_model/tests/test_logistic.py | 61 ++++++++-------- 2 files changed, 65 insertions(+), 76 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 2bf0b5c9b9b..7eaea7aefde 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -14,9 +14,9 @@ from scipy import optimize, sparse from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator from ..feature_selection.from_model import _LearntSelectorMixin -from ..preprocessing import LabelEncoder, LabelBinarizer +from ..preprocessing import LabelEncoder from ..svm.base import BaseLibLinear -from ..utils import atleast2d_or_csc, as_float_array, check_arrays +from ..utils import atleast2d_or_csc, check_arrays from ..utils.extmath import log_logistic, safe_sparse_dot from ..utils.fixes import expit from ..externals.joblib import Parallel, delayed @@ -69,7 +69,7 @@ def _logistic_loss_and_grad(w, X, y, alpha): Array of labels. alpha : float - Inverse of the cross_validation parameter. + Regularization parameter. alpha is equal to 1 / C. """ _, n_features = X.shape grad = np.empty_like(w) @@ -88,7 +88,7 @@ def _logistic_loss_and_grad(w, X, y, alpha): return out, grad -def _logistic_loss(w, X, y, alpha, fit_intercept=False): +def _logistic_loss(w, X, y, alpha): """Computes the logistic loss. Parameters @@ -99,11 +99,11 @@ def _logistic_loss(w, X, y, alpha, fit_intercept=False): X : {array-like, sparse matrix}, shape (n_samples, n_features) Training data - y : ndarray, shape (n_samples) + y : ndarray, shape (n_samples,) Array of labels. alpha : float - Inverse of the cross_validation parameter. + Regularization parameter. alpha is equal to 1 / C. """ w, c, yz = _intercept_dot(w, X, y) @@ -123,11 +123,11 @@ def _logistic_loss_grad_hess(w, X, y, alpha): X : {array-like, sparse matrix}, shape (n_samples, n_features) Training data - y : ndarray, shape (n_samples) + y : ndarray, shape (n_samples,) Array of labels. alpha : float - Inverse of the cross_validation parameter. + Regularization parameter. alpha is equal to 1 / C. """ n_samples, n_features = X.shape grad = np.empty_like(w) @@ -141,7 +141,6 @@ def _logistic_loss_grad_hess(w, X, y, alpha): z0 = (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w if c is not None: - z0_sum = np.sum(z0) grad[-1] = np.sum(z0) # The mat-vec product of the Hessian @@ -174,8 +173,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha): def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, max_iter=100, tol=1e-4, verbose=0, - solver='liblinear', callback=None, - coef=None, copy=False): + solver='liblinear', coef=None, copy=False): """Compute a Logistic Regression model for a list of regularization parameters. @@ -219,10 +217,6 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, solver : {'lbfgs', 'newton-cg', 'liblinear'} Numerical solver to use. - callback : callable - Function to be called before and after the fit of each regularization - parameter. Must have the signature callback(w, X, y, alpha). - coef: array-lime, shape (n_features,) Initialization value for coefficients of logistic regression. @@ -242,7 +236,7 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, Notes ----- You might get slighly different results with the solver trust-ncg than - with the others since this uses LIBLINEAR penalizes the intercept. + with the others since this uses LIBLINEAR which penalizes the intercept. """ if isinstance(Cs, numbers.Integral): Cs = np.logspace(-4, 4, Cs) @@ -274,8 +268,6 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, coefs = list() for C in Cs: - if callback is not None: - callback(w0, X, y, 1. / C) if solver == 'lbfgs': func = _logistic_loss_and_grad try: @@ -304,8 +296,6 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, else: raise ValueError("solver must be one of {'liblinear', 'lbfgs', " "'newton-cg'}, got '%s' instead" % solver) - if callback is not None: - callback(w0, X, y, 1. / C) coefs.append(w0) return coefs, Cs @@ -335,7 +325,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, 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, integer + Cs: list of floats, ints Each of the values in Cs describes the inverse of regularization strength and must be a positive float. If not provided, then a fixed set of values for Cs are used. @@ -540,14 +530,14 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, _LearntSelectorMixin): - """Logistic Regression (aka logit, MaxEnt) classifier. + """Logistic Regression CV (aka logit, MaxEnt) classifier. This class implements L2 regularized logistic regression using and LBFGS optimizer. Parameters ---------- - Cs: list of floats, integer + Cs: list of floats, ints Each of the values in Cs describes the inverse of regularization strength and must be a positive float. Like in support vector machines, smaller values specify stronger @@ -588,18 +578,17 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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 and C that corresponds to the best score across - each fold is taken and retuned after averaging. + Otherwise the coefs, intercepts and C that correspond to the + best scores across folds are averaged. Attributes ---------- - `coef_` : array, shape (1, n_features) or - (n_classes, n_features) + `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 \ + `coef_` is readonly property derived from `raw_coef_` that follows the internal memory layout of liblinear. `intercept_` : array, shape (1,) or (n_classes,) @@ -626,7 +615,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, an OvA for the corresponding class. Each dict value has shape (n_folds, len(Cs)) - `C_` : array, shape(n_classes,) or (n_classes - 1,) + `C_` : array, shape (n_classes,) or (n_classes - 1,) Array of C that maps to the best scores across every class. See also @@ -683,19 +672,18 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, labels = labels[1:] 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, - method=self.solver, - max_iter=self.max_iter, - tol=self.tol, - verbose=max(0, self.verbose - 1), - scoring=self.scoring, - ) - for label in labels - for train, test in folds - ) + delayed(_log_reg_scoring_path)(X, y, train, test, + pos_class=label, + Cs=self.Cs, + fit_intercept=self.fit_intercept, + method=self.solver, + max_iter=self.max_iter, + tol=self.tol, + verbose=max(0, self.verbose - 1), + scoring=self.scoring) + for label in labels + for train, test in folds + ) coefs_paths, Cs, scores = zip(*fold_coefs_) self.Cs_ = Cs[0] @@ -730,10 +718,10 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, # Take the best scores across every fold and the average of all # coefficients coressponding 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 - ) + 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: diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index d634fdc3831..9f623a17275 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -11,9 +11,11 @@ from sklearn.utils.testing import assert_raises from sklearn.utils.testing import assert_true from sklearn.utils.testing import raises -from sklearn.linear_model.logistic import (LogisticRegression, +from sklearn.linear_model.logistic import ( + LogisticRegression, logistic_regression_path, LogisticRegressionCV, - _logistic_loss_and_grad, _logistic_loss_grad_hess) + _logistic_loss_and_grad, _logistic_loss_grad_hess + ) from sklearn.cross_validation import StratifiedKFold from sklearn.datasets import load_iris, make_classification @@ -197,20 +199,21 @@ def test_logistic_loss_and_grad(): # 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) + 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.) + 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) + 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) @@ -244,10 +247,9 @@ def test_logistic_loss_grad_hess(): 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 - ]) + _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() @@ -256,27 +258,29 @@ def test_logistic_loss_grad_hess(): # 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, 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.) + _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 + """test for LogisticRegressionCV object""" n_samples, n_features = 50, 5 np.random.seed(0) X_ref = np.random.randn(n_samples, n_features) y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) X_ref -= X_ref.mean() X_ref /= X_ref.std() - lr_cv = LogisticRegressionCV(Cs=[1.], fit_intercept=False) + lr_cv = LogisticRegressionCV(Cs=[1.], fit_intercept=False, + solver='liblinear') lr_cv.fit(X_ref, y) - lr = LogisticRegression(C=1, fit_intercept=False) + lr = LogisticRegression(C=1., fit_intercept=False) lr.fit(X_ref, y) - assert_array_almost_equal(lr.coef_, lr_cv.coef_, decimal=2) + assert_array_almost_equal(lr.coef_, lr_cv.coef_) def test_logistic_cv_sparse(): @@ -284,17 +288,15 @@ def test_logistic_cv_sparse(): random_state=np.random.RandomState(0)) X[X < 1.0] = 0.0 csr = sp.csr_matrix(X) - csc = sp.csc_matrix(X) for fit_intercept in [True, False]: clf = LogisticRegressionCV(fit_intercept=fit_intercept) clf.fit(X, y) - for data in [csr, csc]: - clfs = LogisticRegressionCV(fit_intercept=fit_intercept) - clfs.fit(data, y) - assert_array_almost_equal(clfs.coef_, clf.coef_) - assert_array_almost_equal(clfs.intercept_, clf.intercept_) - assert_equal(clfs.C_, clf.C_) + clfs = LogisticRegressionCV(fit_intercept=fit_intercept) + 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(): @@ -329,8 +331,7 @@ def test_intercept_logistic_helper(): assert_almost_equal(hess_interp[-1] + alpha * grad[-1], hess[-1]) -def test_multitask(): - # Right now LogisticRegressionCV cannot handle multiple classes. +def test_multiclass(): X, y = make_classification(n_samples=10, n_features=20, n_informative=10, n_classes=3) clf = LogisticRegressionCV(cv=3) From 7b3bfefaa052f278c5e6a7518e4b5a5ca26b83dd Mon Sep 17 00:00:00 2001 From: MechCoder Date: Fri, 4 Jul 2014 11:53:50 +0200 Subject: [PATCH 36/51] Made the following changes 1. Fixed random state 2. Object dtype support 3. Remove outdated "don't test multiclass" tests --- sklearn/linear_model/logistic.py | 64 ++++++++++++++++----- sklearn/linear_model/tests/test_logistic.py | 25 ++++---- 2 files changed, 64 insertions(+), 25 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 7eaea7aefde..722866052ea 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -18,6 +18,7 @@ from ..preprocessing import LabelEncoder from ..svm.base import BaseLibLinear from ..utils import atleast2d_or_csc, check_arrays from ..utils.extmath import log_logistic, safe_sparse_dot +from ..utils.validation import as_float_array from ..utils.fixes import expit from ..externals.joblib import Parallel, delayed from ..cross_validation import check_cv @@ -28,7 +29,7 @@ from ..metrics import SCORERS # .. some helper functions for logistic_regression_path .. def _intercept_dot(w, X, y): - """Computes y * np.dot(w, X). + """Computes y * np.dot(X, w). It takes into consideration if the intercept should be fit or not. @@ -51,6 +52,7 @@ def _intercept_dot(w, X, y): z = safe_sparse_dot(X, w) if c is not None: z += c + return w, c, y*z @@ -70,6 +72,14 @@ def _logistic_loss_and_grad(w, X, y, alpha): alpha : float Regularization parameter. alpha is equal to 1 / C. + + 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) @@ -104,6 +114,11 @@ def _logistic_loss(w, X, y, alpha): alpha : float Regularization parameter. alpha is equal to 1 / C. + + Returns + ------- + out: float + Logistic loss. """ w, c, yz = _intercept_dot(w, X, y) @@ -128,6 +143,18 @@ def _logistic_loss_grad_hess(w, X, y, alpha): alpha : float Regularization parameter. alpha is equal to 1 / C. + + 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) @@ -217,10 +244,10 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, solver : {'lbfgs', 'newton-cg', 'liblinear'} Numerical solver to use. - coef: array-lime, shape (n_features,) + coef: array-like, shape (n_features,) default None Initialization value for coefficients of logistic regression. - copy: boolean + copy: bool 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 @@ -255,6 +282,9 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, y[mask] = 1 y[~mask] = -1 + # To take care of object dtypes + y = as_float_array(y, copy=False) + if fit_intercept: w0 = np.zeros(X.shape[1] + 1) else: @@ -325,9 +355,10 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, 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, ints + Cs: list of floats | int Each of the values in Cs describes the inverse of - regularization strength and must be a positive float. + 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 @@ -370,6 +401,12 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, y_test[mask] = 1 y_test[~mask] = -1 + # To deal with object dtypes, we need to convert into an array of floats. + X_train = as_float_array(X_train, copy=False) + y_train = as_float_array(y_train, copy=False) + X_test = as_float_array(X_test, copy=False) + 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=method, @@ -532,14 +569,15 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, _LearntSelectorMixin): """Logistic Regression CV (aka logit, MaxEnt) classifier. - This class implements L2 regularized logistic regression using and - LBFGS optimizer. + This class implements L2 regularized logistic regression using liblinear, + newton-cg or LBFGS optimizer. Parameters ---------- - Cs: list of floats, ints + Cs: list of floats | int Each of the values in Cs describes the inverse of regularization - strength and must be a positive float. + 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. @@ -564,14 +602,14 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, tol: float, optional Tolerance for stopping criteria. - max_iter: integer, optional + max_iter: int, optional Maximum number of iterations of the optimization algorithm. - n_jobs : integer, optional + 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 or integer + verbose : bool | int Amount of verbosity. refit : bool @@ -597,7 +635,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, and is of shape(1,) when the problem is binary. `Cs_` : array - Array of C i.e inverse of regularization parameter values used + 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 diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 9f623a17275..e06a0c8ba79 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -165,7 +165,7 @@ def test_consistency_path(): lr = LogisticRegression(C=C, fit_intercept=False, tol=1e-16) lr.fit(X, Y1) lr_coef = lr.coef_.ravel() - assert_array_almost_equal(lr_coef, coefs[i], decimal=1) + assert_array_almost_equal(lr_coef, coefs[i], decimal=3) # test for fit_intercept=True for method in ('lbfgs', 'newton-cg', 'liblinear'): @@ -218,9 +218,10 @@ def test_logistic_loss_and_grad(): def test_logistic_loss_grad_hess(): + rng = np.random.RandomState(0) n_samples, n_features = 50, 5 - X_ref = np.random.randn(n_samples, n_features) - y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) + 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() @@ -270,13 +271,13 @@ def test_logistic_loss_grad_hess(): def test_logistic_cv(): """test for LogisticRegressionCV object""" n_samples, n_features = 50, 5 - np.random.seed(0) - X_ref = np.random.randn(n_samples, n_features) - y = np.sign(X_ref.dot(5 * np.random.randn(n_features))) + 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') + solver='liblinear') lr_cv.fit(X_ref, y) lr = LogisticRegression(C=1., fit_intercept=False) lr.fit(X_ref, y) @@ -285,7 +286,7 @@ def test_logistic_cv(): def test_logistic_cv_sparse(): X, y = make_classification(n_samples=50, n_features=5, - random_state=np.random.RandomState(0)) + random_state=0) X[X < 1.0] = 0.0 csr = sp.csr_matrix(X) @@ -302,7 +303,7 @@ def test_logistic_cv_sparse(): 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=np.random.RandomState(0)) + random_state=0) # Fit intercept case. alpha = 1. @@ -317,14 +318,14 @@ def test_intercept_logistic_helper(): # 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) + 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]) - np.random.seed(0) - grad = np.random.rand(n_features + 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]) From 10a411fe66016a2344670424397e6782e3f166a8 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Tue, 8 Jul 2014 13:23:36 +0200 Subject: [PATCH 37/51] ENH: Logistic Regression now supports newton-cg and lbfgs --- doc/modules/linear_model.rst | 12 ++- sklearn/linear_model/logistic.py | 17 +++- sklearn/linear_model/tests/test_logistic.py | 27 +++++ sklearn/svm/base.py | 106 ++++++++++++++++---- 4 files changed, 133 insertions(+), 29 deletions(-) diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst index 5a1b8580d87..9acd6687822 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 `_. -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,6 +663,12 @@ 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) . +The solvers implemented for Logistic Regression are liblinear (which is a +wrapper around the C++ library, LIBLINEAR), newton-cg and lbfgs. + +The liblinear solver can be used to do L1 or L2 penalized +logistic regression. The lbfgs and newton-cg solvers and do only L1 penalized +regression and are found to converge faster for 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 diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 722866052ea..e06a71e191b 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -230,7 +230,7 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=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 : integer + max_iter : int Maximum number of iterations for the solver. tol : float @@ -446,7 +446,8 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, Parameters ---------- penalty : string, 'l1' or 'l2' - Used to specify the norm used in the penalization. + 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 @@ -479,10 +480,17 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, The 'auto' mode selects weights inversely proportional to class frequencies in the training set. + 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. + solver: {'newton-cg', 'lbfgs', 'liblinear'} + Algorithm to use in the optimization problem. + tol: float, optional Tolerance for stopping criteria. @@ -521,12 +529,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. diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index e06a0c8ba79..49b1e312855 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -382,3 +382,30 @@ def test_ova_iris(): 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_) + + +def test_logreg_newton_lbfgs(): + X, y = make_classification(n_features=50, n_informative=10, 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_logreg_newton_lbfgs_multitask(): + X, y = make_classification(n_features=50, 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=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) diff --git a/sklearn/svm/base.py b/sklearn/svm/base.py index c7934c48381..595e51a5a00 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,11 @@ 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 + self._enc = LabelEncoder() y_ind = self._enc.fit_transform(y) if len(self.classes_) < 2: @@ -686,31 +694,85 @@ 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() + self.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) + + 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) From e4111290b601c05c90455ff557396034c70ae707 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Wed, 9 Jul 2014 14:35:29 +0200 Subject: [PATCH 38/51] ENH: Weighted logistic regression for lbfgs and newton-cg --- .../supervised_learning.rst | 4 +- sklearn/linear_model/logistic.py | 83 +++++++++++++++---- sklearn/svm/base.py | 9 +- 3 files changed, 73 insertions(+), 23 deletions(-) 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/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index e06a71e191b..a3a89404997 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -16,7 +16,7 @@ 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 atleast2d_or_csc, check_arrays +from ..utils import atleast2d_or_csc, check_arrays, compute_class_weight from ..utils.extmath import log_logistic, safe_sparse_dot from ..utils.validation import as_float_array from ..utils.fixes import expit @@ -56,7 +56,7 @@ def _intercept_dot(w, X, y): return w, c, y*z -def _logistic_loss_and_grad(w, X, y, alpha): +def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None): """Computes the logistic loss and gradient. Parameters @@ -86,11 +86,14 @@ def _logistic_loss_and_grad(w, X, y, alpha): 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(log_logistic(yz)) + .5 * alpha * np.dot(w, w) + out = -np.sum(sample_weight * log_logistic(yz)) + .5 * alpha * np.dot(w, w) z = expit(yz) - z0 = (z - 1) * y + z0 = sample_weight * (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w if c is not None: @@ -98,7 +101,7 @@ def _logistic_loss_and_grad(w, X, y, alpha): return out, grad -def _logistic_loss(w, X, y, alpha): +def _logistic_loss(w, X, y, alpha, sample_weight=None): """Computes the logistic loss. Parameters @@ -122,12 +125,15 @@ def _logistic_loss(w, X, y, alpha): """ 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(log_logistic(yz)) + .5 * alpha * np.dot(w, w) + 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): +def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): """Computes the logistic loss, gradient and the Hessian. Parameters @@ -161,17 +167,21 @@ def _logistic_loss_grad_hess(w, X, y, alpha): 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(log_logistic(yz)) + .5 * alpha * np.dot(w, w) + out = -np.sum(sample_weight * log_logistic(yz)) + .5 * alpha * np.dot(w, w) z = expit(yz) - z0 = (z - 1) * y + z0 = sample_weight * (z - 1) * y + grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w if c is not None: grad[-1] = np.sum(z0) # The mat-vec product of the Hessian - d = z * (1 - z) + 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) @@ -200,7 +210,8 @@ def _logistic_loss_grad_hess(w, X, y, alpha): def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, max_iter=100, tol=1e-4, verbose=0, - solver='liblinear', coef=None, copy=False): + solver='liblinear', coef=None, copy=False, + class_weight=None): """Compute a Logistic Regression model for a list of regularization parameters. @@ -253,6 +264,11 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, is called repeatedly with the same data, as y is modified along the path. + class_weight : ndarray, None + Provide a array of weights corresponding to each class, as obtained + from the output of compute_class_weight. It None, then all classes + are assumed to have weight one. + Returns ------- coefs: ndarray, shape (n_cs, n_features) or (n_cs, n_features + 1) @@ -270,9 +286,9 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, X = atleast2d_or_csc(X, dtype=np.float64) X, y = check_arrays(X, y, copy=copy) + n_classes = np.unique(y) if pos_class is None: - n_classes = np.unique(y) if not (n_classes.size == 2): raise ValueError('To fit OvA, use the pos_class argument') # np.unique(y) gives labels in sorted order. @@ -285,6 +301,12 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, # To take care of object dtypes y = as_float_array(y, copy=False) + if class_weight is None: + class_weight = np.ones(len(n_classes)) + + le = LabelEncoder() + sample_weight = class_weight[le.fit_transform(y)] + if fit_intercept: w0 = np.zeros(X.shape[1] + 1) else: @@ -303,21 +325,22 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, try: out = optimize.fmin_l_bfgs_b( func, w0, fprime=None, - args=(X, y, 1. / C), + args=(X, y, 1. / C, sample_weight), iprint=verbose > 0, 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), + args=(X, y, 1. / C, sample_weight), iprint=verbose > 0, pgtol=tol) w0 = out[0] 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), maxiter=max_iter) + w0, args=(X, y, 1./C, sample_weight), maxiter=max_iter) elif solver == 'liblinear': - lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=tol) + lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=tol, + class_weight=class_weight) lr.fit(X, y) if fit_intercept: w0 = np.concatenate([lr.coef_.ravel(), lr.intercept_]) @@ -333,7 +356,7 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, # 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, + max_iter=100, tol=1e-4, class_weight=None, verbose=0, method='liblinear'): """Computes scores across logistic_regression_path @@ -376,6 +399,11 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, tol : float Tolerance for stopping criteria. + class_weight : ndarray, None + Provide a array of weights corresponding to each class, as obtained + from the output of compute_class_weight. It None, then all classes + are assumed to have weight one. + verbose : int Amount of verbosity @@ -411,6 +439,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, fit_intercept=fit_intercept, solver=method, max_iter=max_iter, + class_weight=class_weight, tol=tol, verbose=verbose) scores = list() @@ -594,6 +623,12 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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. @@ -614,6 +649,12 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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. @@ -672,7 +713,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, """ def __init__(self, Cs=10, fit_intercept=True, cv=None, scoring=None, - solver='newton-cg', tol=1e-4, max_iter=100, + solver='newton-cg', tol=1e-4, max_iter=100, class_weight=None, n_jobs=1, verbose=False, refit=True): self.Cs = Cs self.fit_intercept = fit_intercept @@ -680,6 +721,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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 @@ -712,6 +754,9 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, self.classes_ = labels = np.unique(y) n_classes = len(labels) + self.class_weight_ = compute_class_weight( + self.class_weight, self.classes_, y) + if n_classes == 2: # OvA in case of binary problems is as good as fitting # the higher label @@ -726,6 +771,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, method=self.solver, max_iter=self.max_iter, tol=self.tol, + class_weight=self.class_weight_, verbose=max(0, self.verbose - 1), scoring=self.scoring) for label in labels @@ -758,6 +804,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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, copy=True, + class_weight=self.class_weight_, verbose=max(0, self.verbose - 1)) w = w[0] diff --git a/sklearn/svm/base.py b/sklearn/svm/base.py index 595e51a5a00..d84c6c284d2 100644 --- a/sklearn/svm/base.py +++ b/sklearn/svm/base.py @@ -686,8 +686,11 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)): X = check_array(X, accept_sparse='csr', dtype=np.float64, order="C") - self.class_weight_ = compute_class_weight(self.class_weight, - self.classes_, y) + if not isinstance(self.class_weight, np.ndarray): + self.class_weight_ = compute_class_weight(self.class_weight, + self.classes_, y) + else: + self.class_weight_ = self.class_weight if X.shape[0] != y_ind.shape[0]: raise ValueError("X and y have incompatible shapes.\n" @@ -764,7 +767,7 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)): fit_intercept=self.fit_intercept, tol=self.tol, verbose=self.verbose, solver=self.solver, copy=True, - max_iter=self.max_iter) + max_iter=self.max_iter, class_weight=self.class_weight_) coef_ = coef_[0] if self.fit_intercept: From 5859f8b4fe8292fd6db7f00193d16cc4b0f275eb Mon Sep 17 00:00:00 2001 From: MechCoder Date: Thu, 10 Jul 2014 16:33:39 +0200 Subject: [PATCH 39/51] Changed copy default from True to False, updated docstring for sample_weights --- doc/modules/linear_model.rst | 8 +++---- sklearn/linear_model/logistic.py | 41 ++++++++++++++++++++++---------- 2 files changed, 32 insertions(+), 17 deletions(-) diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst index 9acd6687822..cb96cf033da 100644 --- a/doc/modules/linear_model.rst +++ b/doc/modules/linear_model.rst @@ -666,10 +666,10 @@ Similarly, L1 regularized logistic regression solves the following optimization The solvers implemented for Logistic Regression are liblinear (which is a wrapper around the C++ library, LIBLINEAR), newton-cg and lbfgs. -The liblinear solver can be used to do L1 or L2 penalized -logistic regression. The lbfgs and newton-cg solvers and do only L1 penalized -regression and are found to converge faster for high dimensional data. -L1 penalization yields sparse predicting weights. +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. diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index a3a89404997..85c83e070af 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -44,15 +44,12 @@ def _intercept_dot(w, X, y): y : ndarray, shape (n_samples,) Array of labels """ - c = None + c = 0. if w.size == X.shape[1] + 1: c = w[-1] w = w[:-1] - z = safe_sparse_dot(X, w) - if c is not None: - z += c - + z = safe_sparse_dot(X, w) + c return w, c, y*z @@ -73,6 +70,10 @@ def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None): 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 @@ -96,7 +97,9 @@ def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None): z0 = sample_weight * (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w - if c is not None: + + # Case where we fit the intercept. + if grad.shape[0] > n_features: grad[-1] = z0.sum() return out, grad @@ -118,6 +121,10 @@ def _logistic_loss(w, X, y, alpha, sample_weight=None): 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 @@ -150,6 +157,10 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): 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 @@ -177,8 +188,10 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): z0 = sample_weight * (z - 1) * y grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w - if c is not None: - grad[-1] = np.sum(z0) + + # Case where we fit the intercept. + if grad.shape[0] > n_features: + grad[-1] = z0.sum() # The mat-vec product of the Hessian d = sample_weight * z * (1 - z) @@ -199,7 +212,8 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): ret[:n_features] = X.T.dot(dX.dot(s[:n_features])) ret[:n_features] += alpha * s[:n_features] - if c is not None: + # For the fit intercept case. + if s.shape[0] > n_features: ret[:n_features] += s[-1] * dd_intercept ret[-1] = dd_intercept.dot(s[:n_features]) ret[-1] += d.sum() * s[-1] @@ -210,7 +224,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, max_iter=100, tol=1e-4, verbose=0, - solver='liblinear', coef=None, copy=False, + solver='liblinear', coef=None, copy=True, class_weight=None): """Compute a Logistic Regression model for a list of regularization parameters. @@ -255,10 +269,10 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, solver : {'lbfgs', 'newton-cg', 'liblinear'} Numerical solver to use. - coef: array-like, shape (n_features,) default None + coef: array-like, shape (n_features,), default None Initialization value for coefficients of logistic regression. - copy: bool + 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 @@ -440,6 +454,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, solver=method, max_iter=max_iter, class_weight=class_weight, + copy=False, tol=tol, verbose=verbose) scores = list() @@ -803,7 +818,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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, copy=True, + max_iter=self.max_iter, tol=self.tol, class_weight=self.class_weight_, verbose=max(0, self.verbose - 1)) w = w[0] From 4d516e81f45cf7262e207343a038be3901d0cfb9 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Sat, 12 Jul 2014 20:25:59 +0200 Subject: [PATCH 40/51] DOC: Minor changes --- doc/modules/linear_model.rst | 5 +++-- sklearn/linear_model/logistic.py | 34 ++++++++++++++++---------------- 2 files changed, 20 insertions(+), 19 deletions(-) diff --git a/doc/modules/linear_model.rst b/doc/modules/linear_model.rst index cb96cf033da..78a1cd18174 100644 --- a/doc/modules/linear_model.rst +++ b/doc/modules/linear_model.rst @@ -663,8 +663,9 @@ 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) . -The solvers implemented for Logistic Regression are liblinear (which is a -wrapper around the C++ library, LIBLINEAR), newton-cg and lbfgs. +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 diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 85c83e070af..58e8b0ef93f 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -36,13 +36,13 @@ def _intercept_dot(w, X, y): Parameters ---------- w : ndarray, shape (n_features,) or (n_features + 1,) - Coefficient vector + Coefficient vector. X : {array-like, sparse matrix}, shape (n_samples, n_features) - Training data + Training data. y : ndarray, shape (n_samples,) - Array of labels + Array of labels. """ c = 0. if w.size == X.shape[1] + 1: @@ -59,10 +59,10 @@ def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None): Parameters ---------- w : ndarray, shape (n_features,) or (n_features + 1,) - Coefficient vector + Coefficient vector. X : {array-like, sparse matrix}, shape (n_samples, n_features) - Training data + Training data. y : ndarray, shape (n_samples,) Array of labels. @@ -110,10 +110,10 @@ def _logistic_loss(w, X, y, alpha, sample_weight=None): Parameters ---------- w : ndarray, shape (n_features,) or (n_features + 1,) - Coefficient vector + Coefficient vector. X : {array-like, sparse matrix}, shape (n_samples, n_features) - Training data + Training data. y : ndarray, shape (n_samples,) Array of labels. @@ -146,10 +146,10 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): Parameters ---------- w : ndarray, shape (n_features,) or (n_features + 1,) - Coefficient vector + Coefficient vector. X : {array-like, sparse matrix}, shape (n_samples, n_features) - Training data + Training data. y : ndarray, shape (n_samples,) Array of labels. @@ -227,7 +227,7 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, solver='liblinear', coef=None, copy=True, class_weight=None): """Compute a Logistic Regression model for a list of regularization - parameters. + parameters using l2 regularization. This is an implementation that uses the result of the previous model to speed up computations along the set of solutions, making it faster @@ -236,10 +236,10 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, Parameters ---------- X : array-like or sparse matrix, shape (n_samples, n_features) - Input data + Input data. y : array-like, shape (n_samples,) - Input data, target values + Input data, target values. Cs : array-like or integer of shape (n_cs,) List of values for the regularization parameter or integer specifying @@ -380,13 +380,13 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, Training data. y : array-like, shape (n_samples,) or (n_samples, n_targets) - Target labels + Target labels. train : list of indices - The indices of the train set + The indices of the train set. test : list of indices - The indices of the test set + The indices of the test set. pos_class: int, None The class with respect to which we perform a one-vs-all fit. @@ -419,7 +419,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, are assumed to have weight one. verbose : int - Amount of verbosity + Amount of verbosity. method : {'lbfgs', 'newton-cg', 'liblinear'} Decides which solver to use. @@ -516,7 +516,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 From 86dab1ea860d58e079b58c8e77ae0c232001a982 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Mon, 14 Jul 2014 15:07:21 +0200 Subject: [PATCH 41/51] ENH: Added warnings for convergence, added support for l1 penalty if solver is liblinear --- sklearn/linear_model/logistic.py | 107 +++++++++++++++----- sklearn/linear_model/tests/test_logistic.py | 71 ++++++------- sklearn/svm/base.py | 8 ++ sklearn/utils/optimize.py | 44 +++++++- 4 files changed, 159 insertions(+), 71 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 58e8b0ef93f..ff119d70f00 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -8,6 +8,7 @@ Logistic Regression # Manoj Kumar import numbers +import warnings import numpy as np from scipy import optimize, sparse @@ -21,7 +22,7 @@ from ..utils.extmath import log_logistic, safe_sparse_dot from ..utils.validation import as_float_array from ..utils.fixes import expit from ..externals.joblib import Parallel, delayed -from ..cross_validation import check_cv +from ..cross_validation import _check_cv from ..utils.optimize import newton_cg from ..externals import six from ..metrics import SCORERS @@ -225,9 +226,9 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, max_iter=100, tol=1e-4, verbose=0, solver='liblinear', coef=None, copy=True, - class_weight=None): + class_weight=None, dual=False, penalty='l2'): """Compute a Logistic Regression model for a list of regularization - parameters using l2 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 @@ -251,7 +252,7 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, 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 : boolean + 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). @@ -283,13 +284,25 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, from the output of compute_class_weight. It None, then all classes are assumed to have weight one. + 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. + Returns ------- - coefs: ndarray, shape (n_cs, n_features) or (n_cs, n_features + 1) + 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 trust-ncg than @@ -348,13 +361,19 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, args=(X, y, 1. / C, sample_weight), iprint=verbose > 0, 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) + w0, args=(X, y, 1. /C, sample_weight), + maxiter=max_iter, xtol=tol) elif solver == 'liblinear': lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=tol, - class_weight=class_weight) + class_weight=class_weight, dual=dual, + penalty=penalty) lr.fit(X, y) if fit_intercept: w0 = np.concatenate([lr.coef_.ravel(), lr.intercept_]) @@ -364,14 +383,15 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, raise ValueError("solver must be one of {'liblinear', 'lbfgs', " "'newton-cg'}, got '%s' instead" % solver) coefs.append(w0) - return coefs, Cs + 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, method='liblinear'): + verbose=0, method='liblinear', penalty='l2', + dual=False): """Computes scores across logistic_regression_path Parameters @@ -423,6 +443,15 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, method : {'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. """ log_reg = LogisticRegression(fit_intercept=fit_intercept) @@ -455,7 +484,8 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, max_iter=max_iter, class_weight=class_weight, copy=False, - tol=tol, verbose=verbose) + tol=tol, verbose=verbose, + dual=dual, penalty=penalty) scores = list() @@ -482,20 +512,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' + 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) @@ -622,8 +657,10 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, _LearntSelectorMixin): """Logistic Regression CV (aka logit, MaxEnt) classifier. - This class implements L2 regularized logistic regression using liblinear, - newton-cg or LBFGS optimizer. + 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 ---------- @@ -650,6 +687,15 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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`. @@ -727,12 +773,15 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, """ - def __init__(self, Cs=10, fit_intercept=True, cv=None, scoring=None, - solver='newton-cg', tol=1e-4, max_iter=100, class_weight=None, - n_jobs=1, verbose=False, refit=True): + def __init__(self, Cs=10, fit_intercept=True, cv=None, dual=False, + penalty='l2', scoring=None, solver='newton-cg', tol=1e-4, + max_iter=100, class_weight=None, n_jobs=1, verbose=False, + refit=True): 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 @@ -759,11 +808,19 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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 = atleast2d_or_csc(X, dtype=np.float64) X, y = check_arrays(X, y, copy=False) # init cross-validation generator - cv = check_cv(self.cv, X, y, classifier=True) + cv = _check_cv(self.cv, X, y, classifier=True) folds = list(cv) self.classes_ = labels = np.unique(y) @@ -783,6 +840,8 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, pos_class=label, Cs=self.Cs, fit_intercept=self.fit_intercept, + penalty=self.penalty, + dual=self.dual, method=self.solver, max_iter=self.max_iter, tol=self.tol, diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 49b1e312855..d0d0c44df9c 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -283,6 +283,16 @@ def test_logistic_cv(): 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, @@ -290,14 +300,13 @@ def test_logistic_cv_sparse(): X[X < 1.0] = 0.0 csr = sp.csr_matrix(X) - for fit_intercept in [True, False]: - clf = LogisticRegressionCV(fit_intercept=fit_intercept) - clf.fit(X, y) - clfs = LogisticRegressionCV(fit_intercept=fit_intercept) - 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_) + 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(): @@ -332,43 +341,10 @@ def test_intercept_logistic_helper(): assert_almost_equal(hess_interp[-1] + alpha * grad[-1], hess[-1]) -def test_multiclass(): - X, y = make_classification(n_samples=10, n_features=20, n_informative=10, - n_classes=3) - clf = LogisticRegressionCV(cv=3) - clf.fit(X, y) - - assert_array_equal(clf.coef_.shape, (3, 20)) - 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, 20 + 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_shape_attributes_logregcv(): - n_samples, n_features = 10, 20 - clf = LogisticRegressionCV(cv=3) - X, y = make_classification(n_samples=n_samples, n_features=n_features) - clf.fit(X, y) - - assert_array_equal(clf.coef_.shape, (1, n_features)) - assert_array_equal(clf.classes_, [0, 1]) - assert_equal(len(clf.classes_), 2) - - coefs_paths = np.asarray(list(clf.coefs_paths_.values())) - assert_array_equal(coefs_paths.shape, (1, 3, 10, n_features + 1)) - assert_array_equal(clf.Cs_.shape, (10, )) - scores = np.asarray(list(clf.scores_.values())) - assert_array_equal(scores.shape, (1, 3, 10)) - - 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) @@ -383,6 +359,17 @@ def test_ova_iris(): 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_logreg_newton_lbfgs(): X, y = make_classification(n_features=50, n_informative=10, random_state=0) diff --git a/sklearn/svm/base.py b/sklearn/svm/base.py index d84c6c284d2..cb7d11d5aca 100644 --- a/sklearn/svm/base.py +++ b/sklearn/svm/base.py @@ -678,6 +678,14 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)): # 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: diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index 0c43673c0ff..ddc46b53d95 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -17,9 +17,11 @@ 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): """ @@ -47,15 +49,40 @@ def _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval, return ret + def newton_cg(func_grad_hess, func, grad, x0, args=(), xtol=1e-5, eps=1e-4, - maxiter=100, disp=False): + maxiter=100): """ Minimization of scalar function of one or more variables using the Newton-CG algorithm. - func: callable + Parameters + ---------- + func_grad_hess : callable Should return the value of the function, the gradient, and a - callable returning the matvec product of the Hessian + 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. + + xtol : float + Tolerance. + + eps : float, optional + If fhess is approximated, use this value for the step size. + + maxiter : int + Number of iterations. """ avextol = xtol @@ -68,7 +95,14 @@ def newton_cg(func_grad_hess, func, grad, x0, args=(), xtol=1e-5, eps=1e-4, old_old_fval = None # Outer loop: our Newton iteration - while (np.sum(np.abs(update)) > xtol) and (k < maxiter): + while np.sum(np.abs(update)) > xtol: + + # Early stopping + if k > maxiter: + warnings.warn("newton-cg failed to converge. Increase the " + "number of iterations.") + break + # 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) @@ -134,4 +168,4 @@ if __name__ == "__main__": 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) \ No newline at end of file + out = newton_cg(func_grad_hess, func, x0) From 21a14fd82d0f5779bfefcfc1f2f9b1eb70b74019 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Mon, 14 Jul 2014 18:22:06 +0200 Subject: [PATCH 42/51] FIX: Changed tolerance of newton-cg to be compliant with that of lbfgs --- sklearn/linear_model/logistic.py | 19 +++++++++++++--- sklearn/utils/optimize.py | 37 ++++++++++++++++++-------------- 2 files changed, 37 insertions(+), 19 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index ff119d70f00..eb61ed0cac2 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -367,9 +367,9 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, 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, xtol=tol) + 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, @@ -452,6 +452,19 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, Dual or primal formulation. Dual formulation is only implemented for l2 penalty with liblinear solver. Prefer dual=False when n_samples > n_features. + + 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) diff --git a/sklearn/utils/optimize.py b/sklearn/utils/optimize.py index ddc46b53d95..c1f8403438d 100644 --- a/sklearn/utils/optimize.py +++ b/sklearn/utils/optimize.py @@ -50,7 +50,7 @@ def _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval, return ret -def newton_cg(func_grad_hess, func, grad, x0, args=(), xtol=1e-5, eps=1e-4, +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 @@ -75,38 +75,40 @@ def newton_cg(func_grad_hess, func, grad, x0, args=(), xtol=1e-5, eps=1e-4, args: tuple, optional Arguments passed to func_grad_hess, func and grad. - xtol : float - Tolerance. + 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. - """ - avextol = xtol + Returns + ------- + xk : float + Estimated minimum. + """ x0 = np.asarray(x0).flatten() - xtol = len(x0) * avextol - update = [2 * xtol] xk = x0 - k = 0 + k = 1 old_fval = func(x0, *args) old_old_fval = None # Outer loop: our Newton iteration - while np.sum(np.abs(update)) > xtol: - - # Early stopping - if k > maxiter: - warnings.warn("newton-cg failed to converge. Increase the " - "number of iterations.") - break + 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) - maggrad = np.sum(np.abs(fgrad)) + + 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) @@ -151,6 +153,9 @@ def newton_cg(func_grad_hess, func, grad, x0, args=(), xtol=1e-5, eps=1e-4, 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 From 8eb3d30ea89db9ff5550d9954e74b74c1f095387 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Tue, 15 Jul 2014 12:16:27 +0200 Subject: [PATCH 43/51] COSMIT: Utils imports are together --- sklearn/linear_model/logistic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index eb61ed0cac2..5e8fff599df 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -19,11 +19,11 @@ from ..preprocessing import LabelEncoder from ..svm.base import BaseLibLinear from ..utils import atleast2d_or_csc, check_arrays, 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 from ..utils.fixes import expit from ..externals.joblib import Parallel, delayed from ..cross_validation import _check_cv -from ..utils.optimize import newton_cg from ..externals import six from ..metrics import SCORERS From 0634ddf50042642b09047e960af5b366bf3dc46c Mon Sep 17 00:00:00 2001 From: MechCoder Date: Tue, 15 Jul 2014 18:59:58 +0200 Subject: [PATCH 44/51] FIX: Class weights are computed for each OvA --- sklearn/linear_model/logistic.py | 64 +++++++++++++++++--------------- sklearn/svm/base.py | 12 +++--- 2 files changed, 40 insertions(+), 36 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 5e8fff599df..39224b02889 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -279,10 +279,11 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, is called repeatedly with the same data, as y is modified along the path. - class_weight : ndarray, None - Provide a array of weights corresponding to each class, as obtained - from the output of compute_class_weight. It None, then all classes - are assumed to have weight one. + 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 @@ -305,7 +306,7 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, Notes ----- - You might get slighly different results with the solver trust-ncg than + 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): @@ -316,23 +317,29 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, n_classes = np.unique(y) if pos_class is None: - if not (n_classes.size == 2): + 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): + 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 is None: - class_weight = np.ones(len(n_classes)) - - le = LabelEncoder() - sample_weight = class_weight[le.fit_transform(y)] + 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) @@ -391,7 +398,7 @@ 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, method='liblinear', penalty='l2', - dual=False): + dual=False, copy=True): """Computes scores across logistic_regression_path Parameters @@ -433,10 +440,11 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, tol : float Tolerance for stopping criteria. - class_weight : ndarray, None - Provide a array of weights corresponding to each class, as obtained - from the output of compute_class_weight. It None, then all classes - are assumed to have weight one. + 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. @@ -477,18 +485,11 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, y_test = y[test] if pos_class is not None: - # In order to avoid a copy in y, mask test and train separately - mask = (y_train == pos_class) - y_train[mask] = 1 - y_train[~mask] = -1 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. - X_train = as_float_array(X_train, copy=False) - y_train = as_float_array(y_train, copy=False) - X_test = as_float_array(X_test, copy=False) y_test = as_float_array(y_test, copy=False) coefs, Cs = logistic_regression_path(X_train, y_train, Cs=Cs, @@ -496,7 +497,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, solver=method, max_iter=max_iter, class_weight=class_weight, - copy=False, + copy=copy, pos_class=pos_class, tol=tol, verbose=verbose, dual=dual, penalty=penalty) @@ -839,8 +840,8 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, self.classes_ = labels = np.unique(y) n_classes = len(labels) - self.class_weight_ = compute_class_weight( - self.class_weight, self.classes_, y) + 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 @@ -848,6 +849,11 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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, @@ -858,7 +864,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, method=self.solver, max_iter=self.max_iter, tol=self.tol, - class_weight=self.class_weight_, + class_weight=self.class_weight, verbose=max(0, self.verbose - 1), scoring=self.scoring) for label in labels @@ -891,7 +897,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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_, + class_weight=self.class_weight, verbose=max(0, self.verbose - 1)) w = w[0] diff --git a/sklearn/svm/base.py b/sklearn/svm/base.py index cb7d11d5aca..75a6c753c8d 100644 --- a/sklearn/svm/base.py +++ b/sklearn/svm/base.py @@ -694,11 +694,9 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)): X = check_array(X, accept_sparse='csr', dtype=np.float64, order="C") - if not isinstance(self.class_weight, np.ndarray): - self.class_weight_ = compute_class_weight(self.class_weight, - self.classes_, y) - else: - self.class_weight_ = self.class_weight + # Used in the liblinear solver. + self.class_weight_ = compute_class_weight(self.class_weight, + self.classes_, y) if X.shape[0] != y_ind.shape[0]: raise ValueError("X and y have incompatible shapes.\n" @@ -707,7 +705,7 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)): if self.solver not in ['liblinear', 'newton-cg', 'lbfgs']: raise ValueError("Logistic Regression supports only liblinear," - "newton-cg and lbfgs solvers.") + " newton-cg and lbfgs solvers.") if self.solver == 'liblinear': liblinear.set_verbosity_wrap(self.verbose) @@ -775,7 +773,7 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)): 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_) + max_iter=self.max_iter, class_weight=self.class_weight) coef_ = coef_[0] if self.fit_intercept: From f7db8cdf8f0649e66cdd19b82462fb91332b73ac Mon Sep 17 00:00:00 2001 From: MechCoder Date: Thu, 17 Jul 2014 16:51:43 +0200 Subject: [PATCH 45/51] FIX: Liblinear solver for LogisticRegressionCV now works for class_weight==auto --- sklearn/linear_model/logistic.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 39224b02889..7cf00487243 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -71,7 +71,7 @@ def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None): alpha : float Regularization parameter. alpha is equal to 1 / C. - sample_weight: ndarray, shape (n_samples,) optional + 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. @@ -122,7 +122,7 @@ def _logistic_loss(w, X, y, alpha, sample_weight=None): alpha : float Regularization parameter. alpha is equal to 1 / C. - sample_weight: ndarray, shape (n_samples,) optional + 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. @@ -158,7 +158,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): alpha : float Regularization parameter. alpha is equal to 1 / C. - sample_weight: ndarray, shape (n_samples,) optional + 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. @@ -327,9 +327,23 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, # 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): - class_weight = compute_class_weight(class_weight, n_classes, y) - sample_weight = class_weight[le.fit_transform(y)] + 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 @@ -338,8 +352,8 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, # 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)] + 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) From c8540df986d26f46bb980401578031fbc377be41 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Fri, 18 Jul 2014 23:46:32 +0200 Subject: [PATCH 46/51] FIX: Add DataConversionWarning --- sklearn/linear_model/logistic.py | 11 ++++++++++- sklearn/linear_model/tests/test_logistic.py | 3 --- 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 7cf00487243..f6c2564da74 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -20,7 +20,7 @@ from ..svm.base import BaseLibLinear from ..utils import atleast2d_or_csc, check_arrays, 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 +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 @@ -847,6 +847,15 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, X = atleast2d_or_csc(X, dtype=np.float64) X, y = check_arrays(X, y, copy=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) + # init cross-validation generator cv = _check_cv(self.cv, X, y, classifier=True) folds = list(cv) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index d0d0c44df9c..7d1ba738d49 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -132,9 +132,6 @@ 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 = LogisticRegression(random_state=0) clf.fit(X, Y1) clf.coef_[:] = 0 From a52174b22c0ea93f2d83f5f1272b5bdcbd426f30 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Mon, 21 Jul 2014 16:34:43 +0200 Subject: [PATCH 47/51] FIX: Changes due to recent refactoring of the check functions --- sklearn/linear_model/logistic.py | 20 ++++++++++++-------- sklearn/svm/base.py | 18 ------------------ 2 files changed, 12 insertions(+), 26 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index f6c2564da74..e4852c3dc1b 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -17,7 +17,7 @@ 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 atleast2d_or_csc, check_arrays, compute_class_weight +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 @@ -176,6 +176,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): """ 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) @@ -191,7 +192,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): grad[:n_features] = safe_sparse_dot(X.T, z0) + alpha * w # Case where we fit the intercept. - if grad.shape[0] > n_features: + if fit_intercept: grad[-1] = z0.sum() # The mat-vec product of the Hessian @@ -203,7 +204,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): # Precompute as much as possible dX = d[:, np.newaxis] * X - if c is not None: + 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))) @@ -214,7 +215,7 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): ret[:n_features] += alpha * s[:n_features] # For the fit intercept case. - if s.shape[0] > n_features: + if fit_intercept: ret[:n_features] += s[-1] * dd_intercept ret[-1] = dd_intercept.dot(s[:n_features]) ret[-1] += d.sum() * s[-1] @@ -312,8 +313,9 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, if isinstance(Cs, numbers.Integral): Cs = np.logspace(-4, 4, Cs) - X = atleast2d_or_csc(X, dtype=np.float64) - X, y = check_arrays(X, y, copy=copy) + 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: @@ -844,8 +846,8 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, raise ValueError("newton-cg and lbfgs solvers support only " "the primal form.") - X = atleast2d_or_csc(X, dtype=np.float64) - X, y = check_arrays(X, y, copy=False) + 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( @@ -856,6 +858,8 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, ) 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) diff --git a/sklearn/svm/base.py b/sklearn/svm/base.py index 75a6c753c8d..8a92fe62c31 100644 --- a/sklearn/svm/base.py +++ b/sklearn/svm/base.py @@ -710,24 +710,6 @@ class BaseLibLinear(six.with_metaclass(ABCMeta, BaseEstimator)): if self.solver == 'liblinear': liblinear.set_verbosity_wrap(self.verbose) - 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() - self.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='') From 67585f6e6f1ee5912fb520befe241fa74db8c37d Mon Sep 17 00:00:00 2001 From: MechCoder Date: Tue, 22 Jul 2014 16:16:40 +0200 Subject: [PATCH 48/51] MAINT: Improve documentation and coverage --- sklearn/linear_model/logistic.py | 139 +++++++++++++------- sklearn/linear_model/tests/test_logistic.py | 29 +++- sklearn/tests/test_common.py | 7 +- 3 files changed, 123 insertions(+), 52 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index e4852c3dc1b..15b961f7e16 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -51,7 +51,7 @@ def _intercept_dot(w, X, y): w = w[:-1] z = safe_sparse_dot(X, w) + c - return w, c, y*z + return w, c, y * z def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None): @@ -77,10 +77,10 @@ def _logistic_loss_and_grad(w, X, y, alpha, sample_weight=None): Returns ------- - out: float + out : float Logistic loss. - grad: ndarray, shape (n_features,) or (n_features + 1,) + grad : ndarray, shape (n_features,) or (n_features + 1,) Logistic gradient. """ _, n_features = X.shape @@ -128,7 +128,7 @@ def _logistic_loss(w, X, y, alpha, sample_weight=None): Returns ------- - out: float + out : float Logistic loss. """ w, c, yz = _intercept_dot(w, X, y) @@ -164,13 +164,13 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): Returns ------- - out: float + out : float Logistic loss. - grad: ndarray, shape (n_features,) or (n_features + 1,) + grad : ndarray, shape (n_features,) or (n_features + 1,) Logistic gradient. - Hs: callable + Hs : callable Function that takes the gradient as a parameter and returns the matrix product of the Hessian and gradient. """ @@ -226,8 +226,9 @@ def _logistic_loss_grad_hess(w, X, y, alpha, sample_weight=None): def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, max_iter=100, tol=1e-4, verbose=0, - solver='liblinear', coef=None, copy=True, - class_weight=None, dual=False, penalty='l2'): + 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. @@ -243,13 +244,13 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, y : array-like, shape (n_samples,) Input data, target values. - Cs : array-like or integer of shape (n_cs,) + 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 + 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. @@ -261,20 +262,20 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, Maximum number of iterations for the solver. tol : float - Stopping criterion. The iteration will stop when - ``max{|g_i | i = 1, ..., n} <= tol`` + 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 + verbose : int Print convergence message if True. solver : {'lbfgs', 'newton-cg', 'liblinear'} Numerical solver to use. - coef: array-like, shape (n_features,), default None + coef : array-like, shape (n_features,), default None Initialization value for coefficients of logistic regression. - copy: bool, default True + 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 @@ -295,6 +296,18 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, 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) @@ -396,7 +409,8 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, elif solver == 'liblinear': lr = LogisticRegression(C=C, fit_intercept=fit_intercept, tol=tol, class_weight=class_weight, dual=dual, - penalty=penalty) + penalty=penalty, + intercept_scaling=intercept_scaling) lr.fit(X, y) if fit_intercept: w0 = np.concatenate([lr.coef_.ravel(), lr.intercept_]) @@ -413,8 +427,8 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, 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, method='liblinear', penalty='l2', - dual=False, copy=True): + verbose=0, solver='lbfgs', penalty='l2', + dual=False, copy=True, intercept_scaling=1.): """Computes scores across logistic_regression_path Parameters @@ -431,11 +445,11 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, test : list of indices The indices of the test set. - pos_class: int, None + 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 + 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. @@ -451,7 +465,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, term of each coef_ gives us the intercept. max_iter : int - Maximum no. of iterations for the solver. + Maximum number of iterations for the solver. tol : float Tolerance for stopping criteria. @@ -465,7 +479,7 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, verbose : int Amount of verbosity. - method : {'lbfgs', 'newton-cg', 'liblinear'} + solver : {'lbfgs', 'newton-cg', 'liblinear'} Decides which solver to use. penalty : str, 'l1' or 'l2' @@ -477,6 +491,18 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, 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) @@ -510,12 +536,13 @@ def _log_reg_scoring_path(X, y, train, test, pos_class=None, Cs=10, coefs, Cs = logistic_regression_path(X_train, y_train, Cs=Cs, fit_intercept=fit_intercept, - solver=method, + solver=solver, max_iter=max_iter, class_weight=class_weight, copy=copy, pos_class=pos_class, tol=tol, verbose=verbose, - dual=dual, penalty=penalty) + dual=dual, penalty=penalty, + intercept_scaling=intercept_scaling) scores = list() @@ -593,14 +620,14 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, 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) + random_state : int seed, RandomState instance, or None (default) The seed of the pseudo random number generator to use when shuffling the data. - solver: {'newton-cg', 'lbfgs', 'liblinear'} + solver : {'newton-cg', 'lbfgs', 'liblinear'} Algorithm to use in the optimization problem. - tol: float, optional + tol : float, optional Tolerance for stopping criteria. Attributes @@ -608,15 +635,15 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, `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 ----- @@ -683,8 +710,8 @@ class LogisticRegression(BaseLibLinear, LinearClassifierMixin, return np.log(self.predict_proba(X)) -class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, - _LearntSelectorMixin): +class LogisticRegressionCV(LogisticRegression, BaseEstimator, + LinearClassifierMixin, _LearntSelectorMixin): """Logistic Regression CV (aka logit, MaxEnt) classifier. This class implements logistic regression using liblinear, newton-cg or @@ -694,14 +721,14 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, Parameters ---------- - Cs: list of floats | int + 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 + fit_intercept : bool, default: True Specifies if a constant (a.k.a. bias or intercept) should be added the decision function. @@ -726,18 +753,18 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, l2 penalty with liblinear solver. Prefer dual=False when n_samples > n_features. - scoring: callabale + 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'} + solver : {'newton-cg', 'lbfgs', 'liblinear'} Algorithm to use in the optimization problem. - tol: float, optional + tol : float, optional Tolerance for stopping criteria. - max_iter: int, optional + max_iter : int, optional Maximum number of iterations of the optimization algorithm. class_weight : {dict, 'auto'}, optional @@ -760,6 +787,18 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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) @@ -795,7 +834,9 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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. + 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 -------- @@ -804,9 +845,9 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, """ def __init__(self, Cs=10, fit_intercept=True, cv=None, dual=False, - penalty='l2', scoring=None, solver='newton-cg', tol=1e-4, + penalty='l2', scoring=None, solver='lbfgs', tol=1e-4, max_iter=100, class_weight=None, n_jobs=1, verbose=False, - refit=True): + refit=True, intercept_scaling=1.): self.Cs = Cs self.fit_intercept = fit_intercept self.cv = cv @@ -820,6 +861,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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. @@ -864,7 +906,10 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, cv = _check_cv(self.cv, X, y, classifier=True) folds = list(cv) - self.classes_ = labels = np.unique(y) + self._enc = LabelEncoder() + self._enc.fit(y) + + labels = self.classes_ n_classes = len(labels) if n_classes < 2: @@ -888,12 +933,13 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, fit_intercept=self.fit_intercept, penalty=self.penalty, dual=self.dual, - method=self.solver, + 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) + scoring=self.scoring, + intercept_scaling=self.intercept_scaling) for label in labels for train, test in folds ) @@ -930,7 +976,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, else: # Take the best scores across every fold and the average of all - # coefficients coressponding to the best scores. + # coefficients corresponding to the best scores. best_indices = np.argmax(scores, axis=1) w = np.mean([ coefs_paths[i][best_indices[i]] @@ -944,6 +990,7 @@ class LogisticRegressionCV(BaseEstimator, LinearClassifierMixin, 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 7d1ba738d49..5f8fd0546ca 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -339,7 +339,7 @@ def test_intercept_logistic_helper(): def test_ova_iris(): - # Test that our OvA implementation is correct using the iris dataset. + """Test that our OvA implementation is correct using the iris dataset.""" train, target = iris.data, iris.target n_samples, n_features = train.shape @@ -369,7 +369,7 @@ def test_ova_iris(): def test_logreg_newton_lbfgs(): - X, y = make_classification(n_features=50, n_informative=10, random_state=0) + 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) @@ -381,8 +381,8 @@ def test_logreg_newton_lbfgs(): assert_array_almost_equal(clf_n.coef_, clf_lbf.coef_, decimal=3) -def test_logreg_newton_lbfgs_multitask(): - X, y = make_classification(n_features=50, n_informative=10, +def test_logreg_newton_lbfgs_multitask_class_weights(): + 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) @@ -393,3 +393,24 @@ def test_logreg_newton_lbfgs_multitask(): 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_logreg_newton_lbfgs_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) + 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=3) 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 From 2ea804e61f35d87ae4ecb11fcd829ebb8946f42a Mon Sep 17 00:00:00 2001 From: MechCoder Date: Tue, 22 Jul 2014 18:19:04 +0200 Subject: [PATCH 49/51] Update whats new! --- doc/modules/classes.rst | 1 + doc/modules/linear_model.rst | 6 ++++++ doc/whats_new.rst | 6 ++++++ sklearn/linear_model/tests/test_logistic.py | 6 +++--- 4 files changed, 16 insertions(+), 3 deletions(-) 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 78a1cd18174..55a8c5b8e13 100644 --- a/doc/modules/linear_model.rst +++ b/doc/modules/linear_model.rst @@ -692,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/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/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 5f8fd0546ca..7abebd8191f 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -368,7 +368,7 @@ def test_ova_iris(): assert_array_equal(scores.shape, (3, 3, 10)) -def test_logreg_newton_lbfgs(): +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) @@ -381,7 +381,7 @@ def test_logreg_newton_lbfgs(): assert_array_almost_equal(clf_n.coef_, clf_lbf.coef_, decimal=3) -def test_logreg_newton_lbfgs_multitask_class_weights(): +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) @@ -395,7 +395,7 @@ def test_logreg_newton_lbfgs_multitask_class_weights(): assert_array_almost_equal(clf_n.coef_, clf_lbf.coef_, decimal=3) -def test_logreg_newton_lbfgs_class_weights(): +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) From 4de4849e8aae09b09e8e6e469bfbf23ebbea7012 Mon Sep 17 00:00:00 2001 From: MechCoder Date: Wed, 23 Jul 2014 13:05:43 +0200 Subject: [PATCH 50/51] FIX: Fixes for the cross_validation failure --- sklearn/linear_model/logistic.py | 4 ++-- sklearn/linear_model/tests/test_logistic.py | 9 ++++++--- sklearn/utils/estimator_checks.py | 5 +++-- 3 files changed, 11 insertions(+), 7 deletions(-) diff --git a/sklearn/linear_model/logistic.py b/sklearn/linear_model/logistic.py index 15b961f7e16..2162092affa 100644 --- a/sklearn/linear_model/logistic.py +++ b/sklearn/linear_model/logistic.py @@ -389,13 +389,13 @@ def logistic_regression_path(X, y, pos_class=None, Cs=10, fit_intercept=True, out = optimize.fmin_l_bfgs_b( func, w0, fprime=None, args=(X, y, 1. / C, sample_weight), - iprint=verbose > 0, pgtol=tol, maxiter=max_iter) + 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, pgtol=tol) + iprint=(verbose > 0) - 1, pgtol=tol) w0 = out[0] if out[2]["warnflag"] == 1: warnings.warn("lbfgs failed to converge. Increase the number " diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index 7abebd8191f..b7d416cd5c6 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -10,6 +10,7 @@ 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.logistic import ( LogisticRegression, @@ -153,10 +154,11 @@ def test_nan(): def test_consistency_path(): """Test that the path algorithm is consistent""" 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 = logistic_regression_path( + coefs, Cs = f(logistic_regression_path)( X, Y1, 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) @@ -167,7 +169,7 @@ def test_consistency_path(): # test for fit_intercept=True for method in ('lbfgs', 'newton-cg', 'liblinear'): Cs = [1e3] - coefs, Cs = logistic_regression_path( + coefs, Cs = f(logistic_regression_path)( X, Y1, Cs=Cs, fit_intercept=True, tol=1e-16, solver=method) lr = LogisticRegression(C=Cs[0], fit_intercept=True, tol=1e-16) lr.fit(X, Y1) @@ -406,7 +408,8 @@ def test_logistic_regressioncv_class_weights(): 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) + 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) 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) From 2374a09246027c3dc16005df6135dc428a52c75b Mon Sep 17 00:00:00 2001 From: MechCoder Date: Wed, 23 Jul 2014 16:04:30 +0200 Subject: [PATCH 51/51] FIX: Increase testing accuracy --- sklearn/linear_model/tests/test_logistic.py | 27 ++++++++++++--------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/sklearn/linear_model/tests/test_logistic.py b/sklearn/linear_model/tests/test_logistic.py index b7d416cd5c6..f547e6a9394 100644 --- a/sklearn/linear_model/tests/test_logistic.py +++ b/sklearn/linear_model/tests/test_logistic.py @@ -153,28 +153,33 @@ def test_nan(): 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, Y1, Cs=Cs, fit_intercept=False, tol=1e-16, solver=method) + 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, Y1) + lr.fit(X, y) lr_coef = lr.coef_.ravel() - assert_array_almost_equal(lr_coef, coefs[i], decimal=3) + 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, Y1, Cs=Cs, fit_intercept=True, tol=1e-16, solver=method) - lr = LogisticRegression(C=Cs[0], fit_intercept=True, tol=1e-16) - lr.fit(X, Y1) + 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=1) + assert_array_almost_equal(lr_coef, coefs[0], decimal=4) def test_liblinear_random_state(): @@ -392,9 +397,9 @@ def test_logistic_regression_solvers_multiclass(): 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) + 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(): @@ -416,4 +421,4 @@ def test_logistic_regressioncv_class_weights(): 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=3) + assert_array_almost_equal(clf_lib.coef_, clf_lbf.coef_, decimal=4)