class_docstr: 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.
docstr: 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.
docstr_tok: ["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","."]
code: 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
  
imports: import warnings
  from abc import ABC, abstractmethod
  import numpy
  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
label_desc: Use image node_id 9 for calling the NewtonSolver obj's underlying member method code with example usage: obj.solve(X, y, sample_weight) and returns: self