code: def _geometric_mean(matrices, init=None, max_iter=10, tol=1e-7):
      """Compute the geometric mean of symmetric positive definite matrices.
  
      The geometric mean of n positive definite matrices
      M_1, ..., M_n is the minimizer of the sum of squared distances from an
      arbitrary matrix to each input matrix M_k
  
      gmean(M_1, ..., M_n) = argmin_X sum_{k=1}^N dist(X, M_k)^2
  
      where the used distance is related to matrices logarithm
  
      dist(X, M_k) = ||log(X^{-1/2} M_k X^{-1/2)}||
  
      In case of positive numbers, this mean is the usual geometric mean.
  
      See Algorithm 3 of :footcite:t:`Fletcher2007`.
  
      References
      ----------
      .. footbibliography::
  
      Parameters
      ----------
      matrices : list of numpy.ndarray, all of shape (n_features, n_features)
          List of matrices whose geometric mean to compute. Raise an error if the
          matrices are not all symmetric positive definite of the same shape.
  
      init : numpy.ndarray, shape (n_features, n_features), optional
          Initialization matrix, default to the arithmetic mean of matrices.
          Raise an error if the matrix is not symmetric positive definite of the
          same shape as the elements of matrices.
  
      max_iter : int, default=10
          Maximal number of iterations.
  
      tol : positive float or None, default=1e-7
          The tolerance to declare convergence: if the gradient norm goes below
          this value, the gradient descent is stopped. If None, no  check is
          performed.
  
      Returns
      -------
      gmean : numpy.ndarray, shape (n_features, n_features)
          Geometric mean of the matrices.
  
      """
      # Shape and symmetry positive definiteness checks
      n_features = matrices[0].shape[0]
      for matrix in matrices:
          _check_square(matrix)
          if matrix.shape[0] != n_features:
              raise ValueError("Matrices are not of the same shape.")
          _check_spd(matrix)
  
      # Initialization
      matrices = np.array(matrices)
      if init is None:
          gmean = np.mean(matrices, axis=0)
      else:
          _check_square(init)
          if init.shape[0] != n_features:
              raise ValueError("Initialization has incorrect shape.")
          _check_spd(init)
          gmean = init
  
      norm_old = np.inf
      step = 1.0
  
      # Gradient descent
      for _ in range(max_iter):
          # Computation of the gradient
          vals_gmean, vecs_gmean = linalg.eigh(gmean)
          gmean_inv_sqrt = _form_symmetric(
              np.sqrt, 1.0 / vals_gmean, vecs_gmean
          )
          whitened_matrices = [
              gmean_inv_sqrt.dot(matrix).dot(gmean_inv_sqrt)
              for matrix in matrices
          ]
          logs = [
              _map_eigenvalues(np.log, w_mat)
              for w_mat in whitened_matrices
          ]
          # Covariant derivative is - gmean.dot(logms_mean)
          logs_mean = np.mean(logs, axis=0)
          if np.any(np.isnan(logs_mean)):
              raise FloatingPointError(
                  "Nan value after logarithm operation."
              )
  
          # Norm of the covariant derivative on the tangent space at point gmean
          norm = np.linalg.norm(logs_mean)
  
          # Update of the minimizer
          vals_log, vecs_log = linalg.eigh(logs_mean)
          gmean_sqrt = _form_symmetric(np.sqrt, vals_gmean, vecs_gmean)
          # Move along the geodesic
          gmean = gmean_sqrt.dot(
              _form_symmetric(np.exp, vals_log * step, vecs_log)
          ).dot(gmean_sqrt)
  
          # Update the norm and the step size
          if norm < norm_old:
              norm_old = norm
          elif norm > norm_old:
              step = step / 2.0
              norm = norm_old
          if tol is not None and norm / gmean.size < tol:
              break
      if tol is not None and norm / gmean.size >= tol:
          warnings.warn(
              f"Maximum number of iterations {max_iter} reached without "
              f"getting to the requested tolerance level {tol}."
          )
  
      return gmean
  
imports: import warnings
  from math import floor, sqrt
  import numpy
  from scipy import linalg
  from sklearn.base import BaseEstimator, TransformerMixin, clone
  from sklearn.covariance import LedoitWolf
  from nilearn._utils.docs import fill_doc
  from ..None import signal
  from .._utils.extmath import is_spd
label_desc: Use image node_id 5 for calling a global function with example usage: _geometric_mean(matrices, init, max_iter, tol) and returns: gmean