Source code for geomstats.geometry.invariant_metric

"""Left- and right- invariant metrics that exist on Lie groups.

Lead authors: Nicolas Guigui and Nina Miolane.
"""

import geomstats.backend as gs
from geomstats.geometry.matrices import Matrices
from geomstats.geometry.riemannian_metric import RiemannianMetric
from geomstats.integrator import integrate
from geomstats.numerics.bvp import ScipySolveBVP
from geomstats.numerics.geodesic import ExpODESolver, LogODESolver
from geomstats.numerics.ivp import ScipySolveIVP
from geomstats.vectorization import repeat_out

EPSILON = 1e-6


[docs] class InvariantMetricMatrixExpODESolver(ExpODESolver): """An exp solver adapted to _InvariantMetricMatrix."""
[docs] def exp(self, tangent_vec, base_point): r"""Compute Riemannian exponential of tan. vector wrt to base point. If :math:`\gamma` is a geodesic, then it satisfies the Euler-Poincare equation [Kolev]_: .. math:: \dot{\gamma}(t) = (dL_{\gamma(t)}) X(t) \dot{X}(t) = ad^*_{X(t)}X(t) where :math:`ad^*` is the dual adjoint map with respect to the metric. For a right-invariant metric, :math:`dR` is used instead of :math:`dL` and :math:`ad^*` is replaced by :math:`-ad^*`. The exponential map is approximated by numerical integration of this equation, with initial conditions :math:`\dot{\gamma}(0)` given by the argument `tangent_vec` and :math:`\gamma(0)` by `base_point`. Parameters ---------- tangent_vec : array-like, shape=[..., n, n] Tangent vector at a base point. base_point : array-like, shape=[..., n, n] Point in the group. Optional, defaults to identity if None. Returns ------- exp : array-like, shape=[..., n, n] Point in the group equal to the Riemannian exponential of tangent_vec at the base point. References ---------- https://en.wikipedia.org/wiki/Runge–Kutta_methods .. [Kolev] Kolev, Boris. “Lie Groups and Mechanics: An Introduction.” Journal of Nonlinear Mathematical Physics 11, no. 4, 2004: 480–98. https://doi.org/10.2991/jnmp.2004.11.4.5. """ base_point, left_angular_vel = self._space.metric._pre_exp( tangent_vec, base_point ) return super().exp(left_angular_vel, base_point)
[docs] def geodesic_ivp(self, tangent_vec, base_point): """Geodesic curve for initial value problem.""" base_point, left_angular_vel = self._space.metric._pre_exp( tangent_vec, base_point ) return super().geodesic_ivp(left_angular_vel, base_point)
[docs] class InvariantMetricMatrixLogODESolver(LogODESolver): """A log solver adapted to _InvariantMetricMatrix."""
[docs] def log(self, point, base_point): """Logarithm map.""" left_angular_vel = super().log(point, base_point) return self._space.to_tangent( self._space.tangent_translation_map( base_point, left=self._space.metric.left )(left_angular_vel), base_point, )
class _InvariantMetricMatrix(RiemannianMetric): """Class for invariant metrics on Matrix Lie groups. This class supports both left and right invariant metrics which exist on Lie groups. Parameters ---------- space : LieGroup Group to equip with the invariant metric. metric_mat_at_identity : array-like, shape=[dim, dim] Matrix that defines the metric at identity. Optional, defaults to identity matrix if None. left : bool Whether to use a left or right invariant metric. Optional, default: True. References ---------- .. [Milnor] Milnor, John. “Curvatures of Left Invariant Metrics on Lie Groups.” Advances in Mathematics 21, no. 3, 1976: 293–329. https://doi.org/10.1016/S0001-8708(76)80002-3. .. [Kolev] Kolev, Boris. “Lie Groups and Mechanics: An Introduction.” Journal of Nonlinear Mathematical Physics 11, no. 4, 2004: 480–98. https://doi.org/10.2991/jnmp.2004.11.4.5. .. [Gallier] Gallier, Jean, and Jocelyn Quaintance. Differential Geometry and Lie Groups: A Computational Perspective. Geonger International Publishing, 2020. https://doi.org/10.1007/978-3-030-46040-2. """ def __init__(self, space, metric_mat_at_identity=None, left=True): super().__init__(space=space) if metric_mat_at_identity is None: metric_mat_at_identity = gs.eye(space.dim) self.metric_mat_at_identity = metric_mat_at_identity self.left = left self._instantiate_solvers() def _instantiate_solvers(self): self.log_solver = InvariantMetricMatrixLogODESolver( self._space, n_nodes=100, use_jac=False, integrator=ScipySolveBVP(tol=1e-8) ) self.exp_solver = InvariantMetricMatrixExpODESolver( self._space, integrator=ScipySolveIVP(atol=1e-8, point_ndim=2) ) @property def reshaped_metric_matrix(self): """Diagonal metric matrix reshaped to a symmetric matrix of size n. Reshape a diagonal metric matrix of size `dim x dim` into a symmetric matrix of size `n x n` where :math:`dim= n (n -1) / 2` is the dimension of the space of skew symmetric matrices. The non-diagonal coefficients in the output matrix correspond to the basis matrices of this space. The diagonal is filled with ones. This useful to compute a matrix inner product. Returns ------- symmetric_matrix : array-like, shape=[n, n] Symmetric matrix. """ if Matrices.is_diagonal(self.metric_mat_at_identity): metric_coeffs = gs.diagonal(self.metric_mat_at_identity) metric_mat = gs.abs( self._space.lie_algebra.matrix_representation(metric_coeffs) ) return metric_mat raise ValueError("This is only possible for a diagonal matrix") def inner_product_at_identity(self, tangent_vec_a, tangent_vec_b): """Compute inner product at tangent space at identity. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] First tangent vector at identity. tangent_vec_b : array-like, shape=[..., n, n] Second tangent vector at identity. Returns ------- inner_prod : array-like, shape=[...] Inner-product of the two tangent vectors. """ tan_b = tangent_vec_b metric_mat = self.metric_mat_at_identity if Matrices.is_diagonal(metric_mat) and self._space.lie_algebra is not None: tan_b = tangent_vec_b * self.reshaped_metric_matrix return Matrices.frobenius_product(tangent_vec_a, tan_b) def inner_product(self, tangent_vec_a, tangent_vec_b, base_point=None): """Compute inner product of two vectors in tangent space at base point. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] First tangent vector at base_point. tangent_vec_b : array-like, shape=[..., n, n] Second tangent vector at base_point. base_point : array-like, shape=[..., n, n] Point in the group. Optional, defaults to identity if None. Returns ------- inner_prod : array-like, shape=[...,] Inner-product of the two tangent vectors. """ if base_point is None: return self.inner_product_at_identity(tangent_vec_a, tangent_vec_b) tangent_translation = self._space.tangent_translation_map( base_point, left=self.left, inverse=True ) tangent_vec_a_at_id = tangent_translation(tangent_vec_a) tangent_vec_b_at_id = tangent_translation(tangent_vec_b) return self.inner_product_at_identity(tangent_vec_a_at_id, tangent_vec_b_at_id) def structure_constant(self, tangent_vec_a, tangent_vec_b, tangent_vec_c): r"""Compute the structure constant of the metric. For three tangent vectors :math:`x, y, z` at identity, compute :math:`<[x,y], z>`. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_c : array-like, shape=[..., n, n] Tangent vector at identity. Returns ------- structure_constant : array-like, shape=[...,] """ bracket = Matrices.bracket(tangent_vec_a, tangent_vec_b) return self.inner_product_at_identity(bracket, tangent_vec_c) def dual_adjoint(self, tangent_vec_a, tangent_vec_b): r"""Compute the metric dual adjoint map. For two tangent vectors at identity :math:`x,y`, this corresponds to the vector :math:`a` such that :math:`\forall z, <[x,z], y > = <a, z>`. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at identity. Returns ------- ad_star : array-like, shape=[..., n, n] Tangent vector at identity corresponding to :math:`ad_x^*(y)`. References ---------- .. [Kolev] Kolev, Boris. “Lie Groups and Mechanics: An Introduction.” Journal of Nonlinear Mathematical Physics 11, no. 4, 2004: 480–98. https://doi.org/10.2991/jnmp.2004.11.4.5. .. [Gallier] Gallier, Jean, and Jocelyn Quaintance. Differential Geometry and Lie Groups: A Computational Perspective. Geonger International Publishing, 2020. https://doi.org/10.1007/978-3-030-46040-2. """ basis = self.normal_basis(self._space.lie_algebra.basis) return -gs.einsum( "i...,ijk->...jk", gs.array( [ self.structure_constant(tan, tangent_vec_a, tangent_vec_b) for tan in basis ] ), gs.array(basis), ) def connection_at_identity(self, tangent_vec_a, tangent_vec_b): r"""Compute the Levi-Civita connection at identity. For two tangent vectors at identity :math:`x,y`, one can associate left (respectively right) invariant vector fields :math:`\tilde{x}, \tilde{y}`. Then the vector :math:`(\nabla_\tilde{x}(\tilde{x}))_{ Id}` is computed using the lie bracket and the dual adjoint map. This is a bilinear map that characterizes the connection [Gallier]_. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at identity. Returns ------- nabla : array-like, shape=[..., n, n] Tangent vector at identity. References ---------- .. [Gallier] Gallier, Jean, and Jocelyn Quaintance. Differential Geometry and Lie Groups: A Computational Perspective. Geonger International Publishing, 2020. https://doi.org/10.1007/978-3-030-46040-2. """ sign = 1.0 if self.left else -1.0 return ( sign / 2 * ( Matrices.bracket(tangent_vec_a, tangent_vec_b) - self.dual_adjoint(tangent_vec_a, tangent_vec_b) - self.dual_adjoint(tangent_vec_b, tangent_vec_a) ) ) def connection(self, tangent_vec_a, tangent_vec_b, base_point=None): r"""Compute the Levi-Civita connection of invariant vector fields. For two tangent vectors at a base point :math:`p, x,y`, one can associate left (respectively right) invariant vector fields :math: `\tilde{x}, \tilde{y}`. Then the vector :math:`(\nabla_\tilde{x}( \tilde{x}))_{p}` is computed using the invariance of the connection and its value at identity [Gallier]_. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at `base_point`. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at `base_point`. base_point : array-like, shape=[..., n, n] Point in the group. Returns ------- nabla : array-like, shape=[..., n, n] Tangent vector at `base_point`. References ---------- .. [Gallier] Gallier, Jean, and Jocelyn Quaintance. Differential Geometry and Lie Groups: A Computational Perspective. Geonger International Publishing, 2020. https://doi.org/10.1007/978-3-030-46040-2. """ if base_point is None: return self.connection_at_identity(tangent_vec_a, tangent_vec_b) translation_map = self._space.tangent_translation_map( base_point, left=self.left, inverse=True ) tan_a_at_id = translation_map(tangent_vec_a) tan_b_at_id = translation_map(tangent_vec_b) translation_map = self._space.tangent_translation_map( base_point, left=self.left, inverse=False ) value_at_id = self.connection_at_identity(tan_a_at_id, tan_b_at_id) return translation_map(value_at_id) def curvature_at_identity(self, tangent_vec_a, tangent_vec_b, tangent_vec_c): r"""Compute the curvature at identity. For three tangent vectors at identity :math:`x,y,z`, the curvature is defined by :math:`R(x, y)z = \nabla_{[x,y]}z - \nabla_x\nabla_y z + \nabla_y\nabla_x z`. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_c : array-like, shape=[..., n, n] Tangent vector at identity. Returns ------- curvature : array-like, shape=[..., n, n] Tangent vector at identity. """ bracket = Matrices.bracket(tangent_vec_a, tangent_vec_b) bracket_term = self.connection_at_identity(bracket, tangent_vec_c) left_term = self.connection_at_identity( tangent_vec_a, self.connection_at_identity(tangent_vec_b, tangent_vec_c) ) right_term = self.connection_at_identity( tangent_vec_b, self.connection_at_identity(tangent_vec_a, tangent_vec_c) ) return bracket_term - left_term + right_term def curvature(self, tangent_vec_a, tangent_vec_b, tangent_vec_c, base_point=None): r"""Compute the curvature. For three tangent vectors at a base point :math:`x,y,z`, the curvature is defined by :math:`R(x, y)z = \nabla_{[x,y]}z - \nabla_x\nabla_y z + \nabla_y\nabla_x z`. It is computed using the invariance of the connection and its value at identity. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at `base_point`. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at `base_point`. tangent_vec_c : array-like, shape=[..., n, n] Tangent vector at `base_point`. base_point : array-like, shape=[..., n, n] Point on the group. Optional, default is the identity. Returns ------- curvature : array-like, shape=[..., n, n] Tangent vector at `base_point`. """ if base_point is None: return self.curvature_at_identity( tangent_vec_a, tangent_vec_b, tangent_vec_c ) translation_map = self._space.tangent_translation_map( base_point, left=self.left, inverse=True ) tan_a_at_id = translation_map(tangent_vec_a) tan_b_at_id = translation_map(tangent_vec_b) tan_c_at_id = translation_map(tangent_vec_c) translation_map = self._space.tangent_translation_map( base_point, left=self.left, inverse=False ) value_at_id = self.curvature_at_identity(tan_a_at_id, tan_b_at_id, tan_c_at_id) return translation_map(value_at_id) def sectional_curvature_at_identity(self, tangent_vec_a, tangent_vec_b): """Compute the sectional curvature at identity. For two orthonormal tangent vectors at identity :math:`x,y`, the sectional curvature is defined by :math:`< R(x, y)x, y>`. Non-orthonormal vectors can be given. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at identity. Returns ------- sectional_curvature : array-like, shape=[...,] Sectional curvature at identity. References ---------- https://en.wikipedia.org/wiki/Sectional_curvature .. [Milnor] Milnor, John. “Curvatures of Left Invariant Metrics on Lie Groups.” Advances in Mathematics 21, no. 3, 1976: 293–329. https://doi.org/10.1016/S0001-8708(76)80002-3. """ curvature = self.curvature_at_identity( tangent_vec_a, tangent_vec_b, tangent_vec_a ) num = self.inner_product(tangent_vec_b, curvature) denom = ( self.squared_norm(tangent_vec_a) * self.squared_norm(tangent_vec_a) - self.inner_product(tangent_vec_a, tangent_vec_b) ** 2 ) condition = gs.isclose(denom, 0.0) denom = gs.where(condition, EPSILON, denom) return gs.where(~condition, num / denom, 0.0) def sectional_curvature(self, tangent_vec_a, tangent_vec_b, base_point=None): """Compute the sectional curvature. For two orthonormal tangent vectors at a base point :math:`x,y`, the sectional curvature is defined by :math:`<R(x, y)x, y>`. Non-orthonormal vectors can be given. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at `base_point`. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at `base_point`. base_point : array-like, shape=[..., n, n] Point in the group. Optional, default is the identity Returns ------- sectional_curvature : array-like, shape=[...,] Sectional curvature at `base_point`. References ---------- https://en.wikipedia.org/wiki/Sectional_curvature .. [Milnor] Milnor, John. “Curvatures of Left Invariant Metrics on Lie Groups.” Advances in Mathematics 21, no. 3, 1976: 293–329. https://doi.org/10.1016/S0001-8708(76)80002-3. """ if base_point is None: return self.sectional_curvature_at_identity(tangent_vec_a, tangent_vec_b) translation_map = self._space.tangent_translation_map( base_point, inverse=True, left=self.left ) tan_a_at_id = translation_map(tangent_vec_a) tan_b_at_id = translation_map(tangent_vec_b) return self.sectional_curvature_at_identity(tan_a_at_id, tan_b_at_id) def curvature_derivative_at_identity( self, tangent_vec_a, tangent_vec_b, tangent_vec_c, tangent_vec_d ): r"""Compute the covariant derivative of the curvature at identity. For four tangent vectors at identity :math:`x, y, z, t`, the covariant derivative of the curvature :math:`(\nabla_x R)(y, z)t` is computed using Leibniz formula. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_c : array-like, shape=[..., n, n] Tangent vector at identity. tangent_vec_d : array-like, shape=[..., n, n] Tangent vector at identity. Returns ------- curvature_derivative : array-like, shape=[..., n, n] Tangent vector at identity. """ first_term = self.connection_at_identity( tangent_vec_a, self.curvature_at_identity(tangent_vec_b, tangent_vec_c, tangent_vec_d), ) second_term = self.curvature_at_identity( self.connection_at_identity(tangent_vec_a, tangent_vec_b), tangent_vec_c, tangent_vec_d, ) third_term = self.curvature_at_identity( tangent_vec_b, self.connection_at_identity(tangent_vec_a, tangent_vec_c), tangent_vec_d, ) fourth_term = self.curvature_at_identity( tangent_vec_b, tangent_vec_c, self.connection_at_identity(tangent_vec_a, tangent_vec_d), ) return first_term - second_term - third_term - fourth_term def curvature_derivative( self, tangent_vec_a, tangent_vec_b, tangent_vec_c, tangent_vec_d, base_point=None, ): r"""Compute the covariant derivative of the curvature. For four tangent vectors at a base point :math:`x, y, z, t`, the covariant derivative of the curvature :math:`(\nabla_x R)(y, z)t` is computed at the base point using Leibniz formula. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] Tangent vector at `base_point`. tangent_vec_b : array-like, shape=[..., n, n] Tangent vector at `base_point`. tangent_vec_c : array-like, shape=[..., n, n] Tangent vector at `base_point`. tangent_vec_d : array-like, shape=[..., n, n] Tangent vector at `base_point`. base_point : array-like, shape=[..., n, n] Point on the group. Returns ------- curvature_derivative : array-like, shape=[..., n, n] Tangent vector at identity. """ if base_point is None: return self.curvature_derivative_at_identity( tangent_vec_a, tangent_vec_b, tangent_vec_c, tangent_vec_d ) translation_map = self._space.tangent_translation_map( base_point, inverse=True, left=self.left ) tan_a_at_id = translation_map(tangent_vec_a) tan_b_at_id = translation_map(tangent_vec_b) tan_c_at_id = translation_map(tangent_vec_c) tan_d_at_id = translation_map(tangent_vec_d) value_at_id = self.curvature_derivative_at_identity( tan_a_at_id, tan_b_at_id, tan_c_at_id, tan_d_at_id ) translation_map = self._space.tangent_translation_map( base_point, inverse=False, left=self.left ) return translation_map(value_at_id) def _pre_exp(self, tangent_vec, base_point): if base_point is None: base_point = self._space.identity left_angular_vel = tangent_vec else: left_angular_vel = self._space.tangent_translation_map( base_point, left=self.left, inverse=True )(tangent_vec) return base_point, self._space.to_tangent(left_angular_vel) def log(self, point, base_point): r"""Compute Riemannian logarithm of a point from a base point. The log is computed by solving an optimization problem. The cost function to be optimized is defined by: .. math:: L(v) = \Vert exp_x(v) - y \Vert^2 where :math:`x,y` are respectively `base_point` and `point`, an extrinsic 2-norm is used, and exp is computed by integration of the Euler-Poincare equation [Kolev]_. Parameters ---------- point : array-like, shape=[..., n, n] Point in the group. base_point : array-like, shape=[..., n, n] Point in the group, from which to compute the log. Optional, default: identity. Returns ------- log : array-like, shape=[..., n, n] Tangent vector at the base point equal to the Riemannian logarithm of point at the base point. References ---------- .. [Kolev] Kolev, Boris. “Lie Groups and Mechanics: An Introduction.” Journal of Nonlinear Mathematical Physics 11, no. 4, 2004: 480–98. https://doi.org/10.2991/jnmp.2004.11.4.5. """ if hasattr(self._space, "are_antipodals") and not gs.all( ~self._space.are_antipodals(point, base_point) ): raise ValueError( "The Logarithm map is not well-defined for" f" antipodal matrices: {point} and {base_point}." ) return self.log_solver.log(point, base_point) def parallel_transport( self, tangent_vec, base_point, direction=None, end_point=None, n_steps=10, step="rk4", return_endpoint=False, ): r"""Compute the parallel transport of a tangent vec along a geodesic. Approximate solution for the parallel transport of a tangent vector a along the geodesic between two points `base_point` and `end_point` or alternatively defined by :math:`t \mapsto exp_{(base\_point)}( t*direction)`. The parallel transport equation is written entirely in the Lie algebra and solved with an integration scheme. Parameters ---------- tangent_vec : array-like, shape=[..., n, n] Tangent vector at base point to be transported. base_point : array-like, shape=[..., n, n] Point on the manifold. direction : array-like, shape=[..., n, n] Tangent vector at base point, along which the parallel transport is computed. Optional, default: None end_point : array-like, shape=[..., n, n] Point on the manifold. Point to transport to. Unused if `tangent_vec_b` is given Optional, default: None n_steps : int Number of integration steps to take. Optional, default : 10. step : str, {'euler', 'rk2', 'rk4'} Scheme to use for the approximation of the solution of the ODE Optional, default : rk4 return_endpoint : bool Whether the end-point of the geodesic should be returned. Optional, default : False. Returns ------- transported_tangent_vec: array-like, shape=[..., n, n] Transported tangent vector at `exp_(base_point)(tangent_vec_b)`. end_point : array-like, shape=[..., n, n] `exp_(base_point)(tangent_vec_b)`, only returned if `return_endpoint` is set to `True`. See Also -------- geomstats.integrator References ---------- .. [GP21] Guigui, Nicolas, and Xavier Pennec. “A Reduced Parallel Transport Equation on Lie Groups with a Left-Invariant Metric.” 5th conference on Geometric Science of Information, Paris 2021. Springer. Lecture Notes in Computer Science. https://hal.inria.fr/hal-03154318. """ if direction is None: if end_point is not None: tangent_vec_b_ = self.log(end_point, base_point) else: raise ValueError( "Either an end_point or a tangent_vec_b must be given to define the" " geodesic along which to transport." ) else: tangent_vec_b_ = direction translation_map = self._space.tangent_translation_map( base_point, left=self.left, inverse=True ) left_angular_vel_a = self._space.to_tangent(translation_map(tangent_vec)) left_angular_vel_b = self._space.to_tangent(translation_map(tangent_vec_b_)) def acceleration(state): """Compute the right-hand-side of the parallel transport eq.""" omega = state[..., 1, :, :] zeta = state[..., 2, :, :] new_state = self.geodesic_equation(state[..., :2, :, :]) gam_dot = new_state[..., 0, :, :] omega_dot = new_state[..., 1, :, :] zeta_dot = -self.connection_at_identity(omega, zeta) return gs.stack([gam_dot, omega_dot, zeta_dot], axis=-3) base_point, left_angular_vel_a, left_angular_vel_b = gs.broadcast_arrays( base_point, left_angular_vel_a, left_angular_vel_b ) initial_state = gs.stack( [base_point, left_angular_vel_b, left_angular_vel_a], axis=-3 ) force = lambda state, _: acceleration(state) flow = integrate(force, initial_state, n_steps=n_steps, step=step) gamma = flow[-1][..., 0, :, :] zeta_t = flow[-1][..., 2, :, :] transported = self._space.tangent_translation_map( gamma, left=self.left, inverse=False )(zeta_t) return (transported, gamma) if return_endpoint else transported def geodesic_equation(self, state): r"""Compute the geodesic ODE associated with the invariant metric. This is a reduced geodesic equation written entirely in the Lie algebra. It is known as Euler-Poincare equation [Kolev]_. .. math:: \dot{\gamma}(t) = (dL_{\gamma(t)}) X(t) \dot{X}(t) = ad^*_{X(t)}X(t) Parameters ---------- state : array-like, shape=[..., dim] Tangent vector at the position. Returns ------- geodesic_ode : array-like, shape=[..., dim] Value of the vector field to be integrated at position. References ---------- .. [Kolev] Kolev, Boris. “Lie Groups and Mechanics: An Introduction.” Journal of Nonlinear Mathematical Physics 11, no. 4, 2004: 480–98. https://doi.org/10.2991/jnmp.2004.11.4.5. """ sign = 1.0 if self.left else -1.0 basis = self.normal_basis(self._space.lie_algebra.basis) point = state[..., 0, :, :] vector = state[..., 1, :, :] velocity = self._space.tangent_translation_map(point, left=self.left)(vector) coefficients = gs.stack( [ self.structure_constant(vector, basis_vector, vector) for basis_vector in basis ], axis=-1, ) acceleration = gs.einsum("...i,ijk->...jk", coefficients, basis) return gs.stack([velocity, sign * acceleration], axis=-3) class _InvariantMetricVector(RiemannianMetric): """Class for invariant metrics on Lie groups represented by vectors. Parameters ---------- space : LieGroup Group to equip with the invariant metric left : bool Whether to use a left or right invariant metric. Optional, default: True. """ def __init__(self, space, left=True): super().__init__(space=space) self.metric_mat_at_identity = gs.eye(space.dim) self.left = left @staticmethod def inner_product_at_identity(tangent_vec_a, tangent_vec_b): """Compute inner product at tangent space at identity. Parameters ---------- tangent_vec_a : array-like, shape=[..., dim] First tangent vector at identity. tangent_vec_b : array-like, shape=[..., dim] Second tangent vector at identity. Returns ------- inner_prod : array-like, shape=[..., dim] Inner-product of the two tangent vectors. """ return gs.dot(tangent_vec_a, tangent_vec_b) def metric_matrix(self, base_point=None): """Compute inner product matrix at the tangent space at a base point. Parameters ---------- base_point : array-like, shape=[..., dim], optional Point in the group (the default is identity). Returns ------- metric_mat : array-like, shape=[..., dim, dim] Metric matrix at base_point. """ if base_point is None: return self.metric_mat_at_identity base_point = self._space.regularize(base_point) jacobian = self._space.jacobian_translation(point=base_point, left=self.left) inv_jacobian = gs.linalg.inv(jacobian) inv_jacobian_transposed = Matrices.transpose(inv_jacobian) return Matrices.mul( inv_jacobian_transposed, self.metric_mat_at_identity, inv_jacobian ) def left_exp_from_identity(self, tangent_vec): """Compute the exponential from identity with the left-invariant metric. Compute Riemannian exponential of a tangent vector at the identity associated to the left-invariant metric. If the method is called by a right-invariant metric, it uses the left-invariant metric associated to the same inner-product matrix at the identity. Parameters ---------- tangent_vec : array-like, shape=[..., dim] Tangent vector at identity. Returns ------- exp : array-like, shape=[..., dim] Point in the group. """ return self._space.regularize_tangent_vec_at_identity(tangent_vec=tangent_vec) def exp_from_identity(self, tangent_vec): """Compute Riemannian exponential of tangent vector from the identity. Parameters ---------- tangent_vec : array-like, shape=[..., dim] Tangent vector at identity. Returns ------- exp : array-like, shape=[..., dim] Point in the group. """ if self.left: exp = self.left_exp_from_identity(tangent_vec) else: opp_left_exp = self.left_exp_from_identity(-tangent_vec) exp = self._space.inverse(opp_left_exp) return self._space.regularize(exp) def exp(self, tangent_vec, base_point=None): """Compute Riemannian exponential of tan. vector wrt to base point. Parameters ---------- tangent_vec : array-like, shape=[..., dim] Tangent vector at a base point. base_point : array-like, shape=[..., dim] Point in the group. Optional, defaults to identity if None. Returns ------- exp : array-like, shape=[..., dim] Point in the group equal to the Riemannian exponential of tangent_vec at the base point. """ identity = self._space.identity if base_point is None: base_point = identity else: base_point = self._space.regularize(base_point) if gs.allclose(base_point, identity): return self.exp_from_identity(tangent_vec) tangent_vec_at_id = self._space.tangent_translation_map( point=base_point, left=self.left, inverse=True )(tangent_vec) exp_from_id = self.exp_from_identity(tangent_vec_at_id) if self.left: exp = self._space.compose(base_point, exp_from_id) else: exp = self._space.compose(exp_from_id, base_point) return self._space.regularize(exp) def left_log_from_identity(self, point): """Compute Riemannian log of a point wrt. id of left-invar. metric. Compute Riemannian logarithm of a point wrt the identity associated to the left-invariant metric. If the method is called by a right-invariant metric, it uses the left-invariant metric associated to the same inner-product matrix at the identity. Parameters ---------- point : array-like, shape=[..., dim] Point in the group. Returns ------- log : array-like, shape=[..., dim] Tangent vector at the identity equal to the Riemannian logarithm of point at the identity. """ point = self._space.regularize(point) return self._space.regularize_tangent_vec_at_identity(tangent_vec=point) def log_from_identity(self, point): """Compute Riemannian logarithm of a point wrt the identity. Parameters ---------- point : array-like, shape=[..., dim] Point in the group. Returns ------- log : array-like, shape=[..., dim] Tangent vector at the identity equal to the Riemannian logarithm of point at the identity. """ point = self._space.regularize(point) if self.left: return self.left_log_from_identity(point) inv_point = self._space.inverse(point) left_log = self.left_log_from_identity(inv_point) return -left_log def log(self, point, base_point=None): """Compute Riemannian logarithm of a point from a base point. Parameters ---------- point : array-like, shape=[..., dim] Point in the group. base_point : array-like, shape=[..., dim], optional Point in the group, from which to compute the log, (the default is identity). Returns ------- log : array-like, shape=[..., dim] Tangent vector at the base point equal to the Riemannian logarithm of point at the base point. """ identity = self._space.identity if base_point is None: base_point = identity else: base_point = self._space.regularize(base_point) if gs.allclose(base_point, identity): return self.log_from_identity(point) point = self._space.regularize(point) if self.left: point_near_id = self._space.compose(self._space.inverse(base_point), point) else: point_near_id = self._space.compose(point, self._space.inverse(base_point)) log_from_id = self.log_from_identity(point_near_id) return self._space.tangent_translation_map(base_point, left=self.left)( log_from_id ) def inner_product(self, tangent_vec_a, tangent_vec_b, base_point=None): """Compute inner product of two vectors in tangent space at base point. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] First tangent vector at base_point. tangent_vec_b : array-like, shape=[..., n, n] Second tangent vector at base_point. base_point : array-like, shape=[..., n, n] Point in the group. Optional, defaults to identity if None. Returns ------- inner_prod : array-like, shape=[...,] Inner-product of the two tangent vectors. """ if base_point is None: return self.inner_product_at_identity(tangent_vec_a, tangent_vec_b) tangent_translation = self._space.tangent_translation_map( base_point, left=self.left, inverse=True ) tangent_vec_a_at_id = tangent_translation(tangent_vec_a) tangent_vec_b_at_id = tangent_translation(tangent_vec_b) return self.inner_product_at_identity(tangent_vec_a_at_id, tangent_vec_b_at_id)
[docs] class InvariantMetric: """Class for invariant metrics on Lie groups. This class supports both left and right invariant metrics which exist on Lie groups. If `point_type == 'vector'`, points are parameterized by the Riemannian logarithm for the canonical left-invariant metric. Parameters ---------- space : LieGroup Group to equip with the invariant metric metric_mat_at_identity : array-like, shape=[dim, dim] Matrix that defines the metric at identity. Optional, defaults to identity matrix if None. left : bool Whether to use a left or right invariant metric. Optional, default: True. """ def __new__(cls, space, metric_mat_at_identity=None, left=True): """Instantiate a special euclidean group. Select the object to instantiate depending on the point_type. """ if space.point_ndim == 1: return _InvariantMetricVector(space, left=left) return _InvariantMetricMatrix( space, left=left, metric_mat_at_identity=metric_mat_at_identity, )
[docs] class BiInvariantMetric(RiemannianMetric): """Class for bi-invariant metrics on compact Lie groups. Compact Lie groups and direct products of compact Lie groups with vector spaces admit bi-invariant metrics [Gallier]_. Products Lie groups are not implemented. Other groups such as SE(3) admit bi-invariant pseudo-metrics. Parameters ---------- space : LieGroup The group to equip with the bi-invariant metric References ---------- .. [Gallier] Gallier, Jean, and Jocelyn Quaintance. Differential Geometry and Lie Groups: A Computational Perspective. Geonger International Publishing, 2020. https://doi.org/10.1007/978-3-030-46040-2. """ def __init__(self, space): self._check_implemented(space) super().__init__(space=space) if self._space.point_ndim == 1: # keeps behavior before removing inheritance self.left = True def _check_implemented(self, space): # TODO (nguigs): implement it for SE(3) if not ("SpecialOrthogonal" in space.__str__() or "SO" in space.__str__()): raise ValueError("The bi-invariant metric is only implemented for SO(n)")
[docs] def exp(self, tangent_vec, base_point=None): """Compute Riemannian exponential of tangent vector from the identity. For a bi-invariant metric, this corresponds to the group exponential. Parameters ---------- tangent_vec : Tangent vector at identity. base_point : array-like, shape=[..., {dim, [n, n]}] Point in the group. Optional, default : identity. Returns ------- exp : array-like, shape=[..., {dim, [n, n]}] Point in the group. References ---------- .. [Gallier] Gallier, Jean, and Jocelyn Quaintance. Differential Geometry and Lie Groups: A Computational Perspective. Geonger International Publishing, 2020. https://doi.org/10.1007/978-3-030-46040-2. """ return self._space.exp(tangent_vec, base_point)
[docs] def log(self, point, base_point=None): """Compute Riemannian logarithm of a point wrt the identity. For a bi-invariant metric this corresponds to the group logarithm. Parameters ---------- point : array-like, shape=[..., {dim, [n, n]}] Point in the group. base_point : array-like, shape=[..., {dim, [n, n]}] Point in the group. Optional, default : identity. Returns ------- log : array-like, shape=[..., {dim, [n, n]}] Tangent vector at the identity equal to the Riemannian logarithm of point at the identity. References ---------- .. [Gallier] Gallier, Jean, and Jocelyn Quaintance. Differential Geometry and Lie Groups: A Computational Perspective. Geonger International Publishing, 2020. https://doi.org/10.1007/978-3-030-46040-2. """ log = self._space.log(point, base_point) return self._space.to_tangent(log, base_point)
[docs] def inner_product_at_identity(self, tangent_vec_a, tangent_vec_b): """Compute inner product at tangent space at identity. Parameters ---------- tangent_vec_a : array-like, shape=[..., {dim, [n, n]}] First tangent vector at identity. tangent_vec_b : array-like, shape=[..., {dim, [n, n]}] Second tangent vector at identity. Returns ------- inner_prod : array-like, shape=[...] Inner-product of the two tangent vectors. """ if self._space.point_ndim == 1: return gs.dot(tangent_vec_a, tangent_vec_b) return Matrices.frobenius_product(tangent_vec_a, tangent_vec_b) / 2
[docs] def inner_product(self, tangent_vec_a, tangent_vec_b, base_point=None): """Compute inner product of two vectors in tangent space at base point. Parameters ---------- tangent_vec_a : array-like, shape=[..., n, n] First tangent vector at base_point. tangent_vec_b : array-like, shape=[..., n, n] Second tangent vector at base_point. base_point : array-like, shape=[..., n, n] Point in the group. Optional, defaults to identity if None. Returns ------- inner_prod : array-like, shape=[...,] Inner-product of the two tangent vectors. """ if base_point is None or self._space.point_ndim == 2: inner_prod = self.inner_product_at_identity(tangent_vec_a, tangent_vec_b) if base_point is not None: return repeat_out( self._space.point_ndim, inner_prod, base_point, tangent_vec_a, tangent_vec_b, ) return inner_prod tangent_translation = self._space.tangent_translation_map( base_point, left=True, inverse=True ) tangent_vec_a_at_id = tangent_translation(tangent_vec_a) tangent_vec_b_at_id = tangent_translation(tangent_vec_b) return self.inner_product_at_identity(tangent_vec_a_at_id, tangent_vec_b_at_id)
[docs] def parallel_transport( self, tangent_vec, base_point, direction=None, end_point=None ): r"""Compute the parallel transport of a tangent vec along a geodesic. Closed-form solution for the parallel transport of a tangent vector a along the geodesic between the base point and an end point, or alternatively defined by :math:`t \mapsto exp_{(base\_point)}(t*direction)`. As a compact Lie group endowed with its canonical bi-invariant metric is a symmetric space, parallel transport is achieved by a geodesic symmetry, or equivalently, one step of the pole ladder scheme. Parameters ---------- tangent_vec : array-like, shape=[..., n, n] Tangent vector at base point to be transported. base_point : array-like, shape=[..., n, n] Point on the manifold. direction : array-like, shape=[..., n, n] Tangent vector at base point, along which the parallel transport is computed. Optional, default: None. end_point : array-like, shape=[..., n, n] Point on the manifold. Point to transport to. Optional, default: None. Returns ------- transported_tangent_vec: array-like, shape=[..., n, n] Transported tangent vector at `end_point=exp_(base_point)(tangent_vec_b)`. """ if direction is None: if end_point is not None: direction = self.log(end_point, base_point) else: raise ValueError( "Either an end_point or a tangent_vec_b must be given to define the" " geodesic along which to transport." ) midpoint = self.exp(1.0 / 2.0 * direction, base_point) transposed = Matrices.transpose(tangent_vec) transported_vec = Matrices.mul(midpoint, transposed, midpoint) return (-1.0) * transported_vec
[docs] def injectivity_radius(self, base_point=None): """Compute the radius of the injectivity domain. This is is the supremum of radii r for which the exponential map is a diffeomorphism from the open ball of radius r centered at the base point onto its image. In the case of a bi-invariant metric, it does not depend on the base point. Parameters ---------- base_point : array-like, shape=[..., n, n] Point on the manifold. Returns ------- radius : array-like, shape=[...,] Injectivity radius. """ radius = gs.array(gs.pi * self._space.dim**0.5) return repeat_out(self._space.point_ndim, radius, base_point)