code: def _multinomial_grad_hess(w, X, Y, alpha, sample_weight):
      """
      Computes the gradient and the Hessian, in the case of a multinomial loss.
  
      Parameters
      ----------
      w : ndarray, shape (n_classes * n_features,) or
          (n_classes * (n_features + 1),)
          Coefficient vector.
  
      X : {array-like, sparse matrix}, shape (n_samples, n_features)
          Training data.
  
      Y : ndarray, shape (n_samples, n_classes)
          Transformed labels according to the output of LabelBinarizer.
  
      alpha : float
          Regularization parameter. alpha is equal to 1 / C.
  
      sample_weight : array-like, shape (n_samples,) optional
          Array of weights that are assigned to individual samples.
  
      Returns
      -------
      grad : array, shape (n_classes * n_features,) or
          (n_classes * (n_features + 1),)
          Ravelled gradient of the multinomial loss.
  
      hessp : callable
          Function that takes in a vector input of shape (n_classes * n_features)
          or (n_classes * (n_features + 1)) and returns matrix-vector product
          with hessian.
  
      References
      ----------
      Barak A. Pearlmutter (1993). Fast Exact Multiplication by the Hessian.
          http://www.bcl.hamilton.ie/~barak/papers/nc-hessian.pdf
      """
      n_features = X.shape[1]
      n_classes = Y.shape[1]
      fit_intercept = w.size == (n_classes * (n_features + 1))
  
      # `loss` is unused. Refactoring to avoid computing it does not
      # significantly speed up the computation and decreases readability
      loss, grad, p = _multinomial_loss_grad(
          w, X, Y, alpha, sample_weight
      )
      sample_weight = sample_weight[:, np.newaxis]
  
      # Hessian-vector product derived by applying the R-operator on the gradient
      # of the multinomial loss function.
      def hessp(v):
          v = v.reshape(n_classes, -1)
          if fit_intercept:
              inter_terms = v[:, -1]
              v = v[:, :-1]
          else:
              inter_terms = 0
          # r_yhat holds the result of applying the R-operator on the multinomial
          # estimator.
          r_yhat = safe_sparse_dot(X, v.T)
          r_yhat += inter_terms
          r_yhat += (-p * r_yhat).sum(axis=1)[:, np.newaxis]
          r_yhat *= p
          r_yhat *= sample_weight
          hessProd = np.zeros(
              (n_classes, n_features + bool(fit_intercept))
          )
          hessProd[:, :n_features] = safe_sparse_dot(r_yhat.T, X)
          hessProd[:, :n_features] += v * alpha
          if fit_intercept:
              hessProd[:, -1] = r_yhat.sum(axis=0)
          return hessProd.ravel()
  
      return grad, hessp
  
imports: import numbers
  import warnings
  import numpy
  from scipy import optimize, sparse
  from .base import LinearClassifierMixin, SparseCoefMixin, BaseEstimator
  from .sag import sag_solver
  from ..feature_selection.from_model import _LearntSelectorMixin
  from ..preprocessing import LabelEncoder, LabelBinarizer
  from ..svm.base import _fit_liblinear
  from ..utils import check_array, check_consistent_length, compute_class_weight
  from ..utils import check_random_state
  from ..utils.extmath import logsumexp, log_logistic, safe_sparse_dot, softmax, squared_norm
  from ..utils.extmath import row_norms
  from ..utils.optimize import newton_cg
  from ..utils.validation import check_X_y
  from ..exceptions import NotFittedError
  from ..utils.fixes import expit
  from ..utils.multiclass import check_classification_targets
  from ..externals.joblib import Parallel, delayed
  from ..model_selection import check_cv
  from ..externals import six
  from ..metrics import SCORERS
label_desc: Use image node_id 7 for calling a global function with example usage: _multinomial_grad_hess(w, X, Y, alpha, sample_weight) and returns: grad,  hessp, hessProd