code: def geodesic_riemann(A, B, alpha=0.5):
      r"""Affine-invariant Riemannian geodesic between SPD/HPD matrices.
  
      The matrix at position :math:`\alpha` on the affine-invariant Riemannian
      geodesic between two SPD/HPD matrices :math:`\mathbf{A}` and
      :math:`\mathbf{B}` is:
  
      .. math::
          \mathbf{C} = \mathbf{A}^{1/2} \left( \mathbf{A}^{-1/2} \mathbf{B}
                       \mathbf{A}^{-1/2} \right)^\alpha \mathbf{A}^{1/2}
  
      :math:`\mathbf{C}` is equal to :math:`\mathbf{A}` if :math:`\alpha` = 0,
      and :math:`\mathbf{B}` if :math:`\alpha` = 1.
  
      Parameters
      ----------
      A : ndarray, shape (..., n, n)
          First SPD/HPD matrices.
      B : ndarray, shape (..., n, n)
          Second SPD/HPD matrices.
      alpha : float, default=0.5
          The position on the geodesic.
  
      Returns
      -------
      C : ndarray, shape (..., n, n)
          SPD/HPD matrices on the affine-invariant Riemannian geodesic.
      """
      sA, isA = sqrtm(A), invsqrtm(A)
      C = isA @ B @ isA
      D = powm(C, alpha)
      E = sA @ D @ sA
      return E
  
imports: from .base import sqrtm, invsqrtm, powm, logm, expm
label_desc: Use image node_id 3 for calling a global function with example usage: geodesic_riemann(A, B, alpha) and returns: E