# Authors: The scikit-learn developers # SPDX-License-Identifier: BSD-3-Clause """ Newton solver for Generalized Linear Models """ import warnings from abc import ABC, abstractmethod import numpy as np import scipy.linalg import scipy.optimize from ..._loss.loss import HalfSquaredError from ...exceptions import ConvergenceWarning from ...utils.optimize import _check_optimize_result from .._linear_loss import LinearModelLoss class NewtonSolver(ABC): """Newton solver for GLMs. This class implements Newton/2nd-order optimization routines for GLMs. Each Newton iteration aims at finding the Newton step which is done by the inner solver. With Hessian H, gradient g and coefficients coef, one step solves: H @ coef_newton = -g For our GLM / LinearModelLoss, we have gradient g and Hessian H: g = X.T @ loss.gradient + l2_reg_strength * coef H = X.T @ diag(loss.hessian) @ X + l2_reg_strength * identity Backtracking line search updates coef = coef_old + t * coef_newton for some t in (0, 1]. This is a base class, actual implementations (child classes) may deviate from the above pattern and use structure specific tricks. Usage pattern: - initialize solver: sol = NewtonSolver(...) - solve the problem: sol.solve(X, y, sample_weight) References ---------- - Jorge Nocedal, Stephen J. Wright. (2006) "Numerical Optimization" 2nd edition https://doi.org/10.1007/978-0-387-40065-5 - Stephen P. Boyd, Lieven Vandenberghe. (2004) "Convex Optimization." Cambridge University Press, 2004. https://web.stanford.edu/~boyd/cvxbook/bv_cvxbook.pdf Parameters ---------- coef : ndarray of shape (n_dof,), (n_classes, n_dof) or (n_classes * n_dof,) Initial coefficients of a linear model. If shape (n_classes * n_dof,), the classes of one feature are contiguous, i.e. one reconstructs the 2d-array via coef.reshape((n_classes, -1), order="F"). linear_loss : LinearModelLoss The loss to be minimized. l2_reg_strength : float, default=0.0 L2 regularization strength. tol : float, default=1e-4 The optimization problem is solved when each of the following condition is fulfilled: 1. maximum |gradient| <= tol 2. Newton decrement d: 1/2 * d^2 <= tol max_iter : int, default=100 Maximum number of Newton steps allowed. n_threads : int, default=1 Number of OpenMP threads to use for the computation of the Hessian and gradient of the loss function. Attributes ---------- coef_old : ndarray of shape coef.shape Coefficient of previous iteration. coef_newton : ndarray of shape coef.shape Newton step. gradient : ndarray of shape coef.shape Gradient of the loss w.r.t. the coefficients. gradient_old : ndarray of shape coef.shape Gradient of previous iteration. loss_value : float Value of objective function = loss + penalty. loss_value_old : float Value of objective function of previous itertion. raw_prediction : ndarray of shape (n_samples,) or (n_samples, n_classes) converged : bool Indicator for convergence of the solver. iteration : int Number of Newton steps, i.e. calls to inner_solve use_fallback_lbfgs_solve : bool If set to True, the solver will resort to call LBFGS to finish the optimisation procedure in case of convergence issues. gradient_times_newton : float gradient @ coef_newton, set in inner_solve and used by line_search. If the Newton step is a descent direction, this is negative. """ def __init__( self, *, coef, linear_loss=LinearModelLoss(base_loss=HalfSquaredError(), fit_intercept=True), l2_reg_strength=0.0, tol=1e-4, max_iter=100, n_threads=1, verbose=0, ): self.coef = coef self.linear_loss = linear_loss self.l2_reg_strength = l2_reg_strength self.tol = tol self.max_iter = max_iter self.n_threads = n_threads self.verbose = verbose def setup(self, X, y, sample_weight): """Precomputations If None, initializes: - self.coef Sets: - self.raw_prediction - self.loss_value """ _, _, self.raw_prediction = self.linear_loss.weight_intercept_raw(self.coef, X) self.loss_value = self.linear_loss.loss( coef=self.coef, X=X, y=y, sample_weight=sample_weight, l2_reg_strength=self.l2_reg_strength, n_threads=self.n_threads, raw_prediction=self.raw_prediction, ) @abstractmethod def update_gradient_hessian(self, X, y, sample_weight): """Update gradient and Hessian.""" @abstractmethod def inner_solve(self, X, y, sample_weight): """Compute Newton step. Sets: - self.coef_newton - self.gradient_times_newton """ def fallback_lbfgs_solve(self, X, y, sample_weight): """Fallback solver in case of emergency. If a solver detects convergence problems, it may fall back to this methods in the hope to exit with success instead of raising an error. Sets: - self.coef - self.converged """ opt_res = scipy.optimize.minimize( self.linear_loss.loss_gradient, self.coef, method="L-BFGS-B", jac=True, options={ "maxiter": self.max_iter - self.iteration, "maxls": 50, # default is 20 "iprint": self.verbose - 1, "gtol": self.tol, "ftol": 64 * np.finfo(np.float64).eps, }, args=(X, y, sample_weight, self.l2_reg_strength, self.n_threads), ) self.iteration += _check_optimize_result("lbfgs", opt_res) self.coef = opt_res.x self.converged = opt_res.status == 0 def line_search(self, X, y, sample_weight): """Backtracking line search. Sets: - self.coef_old - self.coef - self.loss_value_old - self.loss_value - self.gradient_old - self.gradient - self.raw_prediction """ # line search parameters beta, sigma = 0.5, 0.00048828125 # 1/2, 1/2**11 eps = 16 * np.finfo(self.loss_value.dtype).eps t = 1 # step size # gradient_times_newton = self.gradient @ self.coef_newton # was computed in inner_solve. armijo_term = sigma * self.gradient_times_newton _, _, raw_prediction_newton = self.linear_loss.weight_intercept_raw( self.coef_newton, X ) self.coef_old = self.coef self.loss_value_old = self.loss_value self.gradient_old = self.gradient # np.sum(np.abs(self.gradient_old)) sum_abs_grad_old = -1 is_verbose = self.verbose >= 2 if is_verbose: print(" Backtracking Line Search") print(f" eps=16 * finfo.eps={eps}") for i in range(21): # until and including t = beta**20 ~ 1e-6 self.coef = self.coef_old + t * self.coef_newton raw = self.raw_prediction + t * raw_prediction_newton self.loss_value, self.gradient = self.linear_loss.loss_gradient( coef=self.coef, X=X, y=y, sample_weight=sample_weight, l2_reg_strength=self.l2_reg_strength, n_threads=self.n_threads, raw_prediction=raw, ) # Note: If coef_newton is too large, loss_gradient may produce inf values, # potentially accompanied by a RuntimeWarning. # This case will be captured by the Armijo condition. # 1. Check Armijo / sufficient decrease condition. # The smaller (more negative) the better. loss_improvement = self.loss_value - self.loss_value_old check = loss_improvement <= t * armijo_term if is_verbose: print( f" line search iteration={i+1}, step size={t}\n" f" check loss improvement <= armijo term: {loss_improvement} " f"<= {t * armijo_term} {check}" ) if check: break # 2. Deal with relative loss differences around machine precision. tiny_loss = np.abs(self.loss_value_old * eps) check = np.abs(loss_improvement) <= tiny_loss if is_verbose: print( " check loss |improvement| <= eps * |loss_old|:" f" {np.abs(loss_improvement)} <= {tiny_loss} {check}" ) if check: if sum_abs_grad_old < 0: sum_abs_grad_old = scipy.linalg.norm(self.gradient_old, ord=1) # 2.1 Check sum of absolute gradients as alternative condition. sum_abs_grad = scipy.linalg.norm(self.gradient, ord=1) check = sum_abs_grad < sum_abs_grad_old if is_verbose: print( " check sum(|gradient|) < sum(|gradient_old|): " f"{sum_abs_grad} < {sum_abs_grad_old} {check}" ) if check: break t *= beta else: warnings.warn( ( f"Line search of Newton solver {self.__class__.__name__} at" f" iteration #{self.iteration} did no converge after 21 line search" " refinement iterations. It will now resort to lbfgs instead." ), ConvergenceWarning, ) if self.verbose: print(" Line search did not converge and resorts to lbfgs instead.") self.use_fallback_lbfgs_solve = True return self.raw_prediction = raw if is_verbose: print( f" line search successful after {i+1} iterations with " f"loss={self.loss_value}." ) def check_convergence(self, X, y, sample_weight): """Check for convergence. Sets self.converged. """ if self.verbose: print(" Check Convergence") # Note: Checking maximum relative change of coefficient <= tol is a bad # convergence criterion because even a large step could have brought us close # to the true minimum. # coef_step = self.coef - self.coef_old # change = np.max(np.abs(coef_step) / np.maximum(1, np.abs(self.coef_old))) # check = change <= tol # 1. Criterion: maximum |gradient| <= tol # The gradient was already updated in line_search() g_max_abs = np.max(np.abs(self.gradient)) check = g_max_abs <= self.tol if self.verbose: print(f" 1. max |gradient| {g_max_abs} <= {self.tol} {check}") if not check: return # 2. Criterion: For Newton decrement d, check 1/2 * d^2 <= tol # d = sqrt(grad @ hessian^-1 @ grad) # = sqrt(coef_newton @ hessian @ coef_newton) # See Boyd, Vanderberghe (2009) "Convex Optimization" Chapter 9.5.1. d2 = self.coef_newton @ self.hessian @ self.coef_newton check = 0.5 * d2 <= self.tol if self.verbose: print(f" 2. Newton decrement {0.5 * d2} <= {self.tol} {check}") if not check: return if self.verbose: loss_value = self.linear_loss.loss( coef=self.coef, X=X, y=y, sample_weight=sample_weight, l2_reg_strength=self.l2_reg_strength, n_threads=self.n_threads, ) print(f" Solver did converge at loss = {loss_value}.") self.converged = True def finalize(self, X, y, sample_weight): """Finalize the solvers results. Some solvers may need this, others not. """ pass def solve(self, X, y, sample_weight): """Solve the optimization problem. This is the main routine. Order of calls: self.setup() while iteration: self.update_gradient_hessian() self.inner_solve() self.line_search() self.check_convergence() self.finalize() Returns ------- coef : ndarray of shape (n_dof,), (n_classes, n_dof) or (n_classes * n_dof,) Solution of the optimization problem. """ # setup usually: # - initializes self.coef if needed # - initializes and calculates self.raw_predictions, self.loss_value self.setup(X=X, y=y, sample_weight=sample_weight) self.iteration = 1 self.converged = False self.use_fallback_lbfgs_solve = False while self.iteration <= self.max_iter and not self.converged: if self.verbose: print(f"Newton iter={self.iteration}") self.use_fallback_lbfgs_solve = False # Fallback solver. # 1. Update Hessian and gradient self.update_gradient_hessian(X=X, y=y, sample_weight=sample_weight) # TODO: # if iteration == 1: # We might stop early, e.g. we already are close to the optimum, # usually detected by zero gradients at this stage. # 2. Inner solver # Calculate Newton step/direction # This usually sets self.coef_newton and self.gradient_times_newton. self.inner_solve(X=X, y=y, sample_weight=sample_weight) if self.use_fallback_lbfgs_solve: break # 3. Backtracking line search # This usually sets self.coef_old, self.coef, self.loss_value_old # self.loss_value, self.gradient_old, self.gradient, # self.raw_prediction. self.line_search(X=X, y=y, sample_weight=sample_weight) if self.use_fallback_lbfgs_solve: break # 4. Check convergence # Sets self.converged. self.check_convergence(X=X, y=y, sample_weight=sample_weight) # 5. Next iteration self.iteration += 1 if not self.converged: if self.use_fallback_lbfgs_solve: # Note: The fallback solver circumvents check_convergence and relies on # the convergence checks of lbfgs instead. Enough warnings have been # raised on the way. self.fallback_lbfgs_solve(X=X, y=y, sample_weight=sample_weight) else: warnings.warn( ( f"Newton solver did not converge after {self.iteration - 1} " "iterations." ), ConvergenceWarning, ) self.iteration -= 1 self.finalize(X=X, y=y, sample_weight=sample_weight) return self.coef class NewtonCholeskySolver(NewtonSolver): """Cholesky based Newton solver. Inner solver for finding the Newton step H w_newton = -g uses Cholesky based linear solver. """ def setup(self, X, y, sample_weight): super().setup(X=X, y=y, sample_weight=sample_weight) if self.linear_loss.base_loss.is_multiclass: # Easier with ravelled arrays, e.g., for scipy.linalg.solve. # As with LinearModelLoss, we always are contiguous in n_classes. self.coef = self.coef.ravel(order="F") # Note that the computation of gradient in LinearModelLoss follows the shape of # coef. self.gradient = np.empty_like(self.coef) # But the hessian is always 2d. n = self.coef.size self.hessian = np.empty_like(self.coef, shape=(n, n)) # To help case distinctions. self.is_multinomial_with_intercept = ( self.linear_loss.base_loss.is_multiclass and self.linear_loss.fit_intercept ) self.is_multinomial_no_penalty = ( self.linear_loss.base_loss.is_multiclass and self.l2_reg_strength == 0 ) def update_gradient_hessian(self, X, y, sample_weight): _, _, self.hessian_warning = self.linear_loss.gradient_hessian( coef=self.coef, X=X, y=y, sample_weight=sample_weight, l2_reg_strength=self.l2_reg_strength, n_threads=self.n_threads, gradient_out=self.gradient, hessian_out=self.hessian, raw_prediction=self.raw_prediction, # this was updated in line_search ) def inner_solve(self, X, y, sample_weight): if self.hessian_warning: warnings.warn( ( f"The inner solver of {self.__class__.__name__} detected a " "pointwise hessian with many negative values at iteration " f"#{self.iteration}. It will now resort to lbfgs instead." ), ConvergenceWarning, ) if self.verbose: print( " The inner solver detected a pointwise Hessian with many " "negative values and resorts to lbfgs instead." ) self.use_fallback_lbfgs_solve = True return # Note: The following case distinction could also be shifted to the # implementation of HalfMultinomialLoss instead of here within the solver. if self.is_multinomial_no_penalty: # The multinomial loss is overparametrized for each unpenalized feature, so # at least the intercepts. This can be seen by noting that predicted # probabilities are invariant under shifting all coefficients of a single # feature j for all classes by the same amount c: # coef[k, :] -> coef[k, :] + c => proba stays the same # where we have assumned coef.shape = (n_classes, n_features). # Therefore, also the loss (-log-likelihood), gradient and hessian stay the # same, see # Noah Simon and Jerome Friedman and Trevor Hastie. (2013) "A Blockwise # Descent Algorithm for Group-penalized Multiresponse and Multinomial # Regression". https://doi.org/10.48550/arXiv.1311.6529 # # We choose the standard approach and set all the coefficients of the last # class to zero, for all features including the intercept. n_classes = self.linear_loss.base_loss.n_classes n_dof = self.coef.size // n_classes # degree of freedom per class n = self.coef.size - n_dof # effective size self.coef[n_classes - 1 :: n_classes] = 0 self.gradient[n_classes - 1 :: n_classes] = 0 self.hessian[n_classes - 1 :: n_classes, :] = 0 self.hessian[:, n_classes - 1 :: n_classes] = 0 # We also need the reduced variants of gradient and hessian where the # entries set to zero are removed. For 2 features and 3 classes with # arbitrary values, "x" means removed: # gradient = [0, 1, x, 3, 4, x] # # hessian = [0, 1, x, 3, 4, x] # [1, 7, x, 9, 10, x] # [x, x, x, x, x, x] # [3, 9, x, 21, 22, x] # [4, 10, x, 22, 28, x] # [x, x, x, x, x, x] # The following slicing triggers copies of gradient and hessian. gradient = self.gradient.reshape(-1, n_classes)[:, :-1].flatten() hessian = self.hessian.reshape(n_dof, n_classes, n_dof, n_classes)[ :, :-1, :, :-1 ].reshape(n, n) elif self.is_multinomial_with_intercept: # Here, only intercepts are unpenalized. We again choose the last class and # set its intercept to zero. self.coef[-1] = 0 self.gradient[-1] = 0 self.hessian[-1, :] = 0 self.hessian[:, -1] = 0 gradient, hessian = self.gradient[:-1], self.hessian[:-1, :-1] else: gradient, hessian = self.gradient, self.hessian try: with warnings.catch_warnings(): warnings.simplefilter("error", scipy.linalg.LinAlgWarning) self.coef_newton = scipy.linalg.solve( hessian, -gradient, check_finite=False, assume_a="sym" ) if self.is_multinomial_no_penalty: self.coef_newton = np.c_[ self.coef_newton.reshape(n_dof, n_classes - 1), np.zeros(n_dof) ].reshape(-1) assert self.coef_newton.flags.f_contiguous elif self.is_multinomial_with_intercept: self.coef_newton = np.r_[self.coef_newton, 0] self.gradient_times_newton = self.gradient @ self.coef_newton if self.gradient_times_newton > 0: if self.verbose: print( " The inner solver found a Newton step that is not a " "descent direction and resorts to LBFGS steps instead." ) self.use_fallback_lbfgs_solve = True return except (np.linalg.LinAlgError, scipy.linalg.LinAlgWarning) as e: warnings.warn( f"The inner solver of {self.__class__.__name__} stumbled upon a " "singular or very ill-conditioned Hessian matrix at iteration " f"{self.iteration}. It will now resort to lbfgs instead.\n" "Further options are to use another solver or to avoid such situation " "in the first place. Possible remedies are removing collinear features" " of X or increasing the penalization strengths.\n" "The original Linear Algebra message was:\n" + str(e), scipy.linalg.LinAlgWarning, ) # Possible causes: # 1. hess_pointwise is negative. But this is already taken care in # LinearModelLoss.gradient_hessian. # 2. X is singular or ill-conditioned # This might be the most probable cause. # # There are many possible ways to deal with this situation. Most of them # add, explicitly or implicitly, a matrix to the hessian to make it # positive definite, confer to Chapter 3.4 of Nocedal & Wright 2nd ed. # Instead, we resort to lbfgs. if self.verbose: print( " The inner solver stumbled upon an singular or ill-conditioned " "Hessian matrix and resorts to LBFGS instead." ) self.use_fallback_lbfgs_solve = True return def finalize(self, X, y, sample_weight): if self.is_multinomial_no_penalty: # Our convention is usually the symmetric parametrization where # sum(coef[classes, features], axis=0) = 0. # We convert now to this convention. Note that it does not change # the predicted probabilities. n_classes = self.linear_loss.base_loss.n_classes self.coef = self.coef.reshape(n_classes, -1, order="F") self.coef -= np.mean(self.coef, axis=0) elif self.is_multinomial_with_intercept: # Only the intercept needs an update to the symmetric parametrization. n_classes = self.linear_loss.base_loss.n_classes self.coef[-n_classes:] -= np.mean(self.coef[-n_classes:])