API Reference

Aggregation

groupoid.aggregation.TransportGroupoidAggregator dataclass

Federated aggregator using explicit invertible point actions.

The current point-valued aggregation path is scientifically supported for caller-supplied matrices that act invertibly on the chosen point representation and preserve the manifold domain used by the Karcher mean. The preregistered S^2 benchmark exercises this contract with explicit SO(3) rotations. Arbitrary square matrices are not thereby validated as geometric transport morphisms.

consistency_threshold is a threshold on the basis-dependent holonomy defect. A finite threshold decision is not invariant under general non-orthogonal changes of frame; it must not be interpreted as a canonical cohomological verdict.

Source code in groupoid/aggregation.py

@dataclass
class TransportGroupoidAggregator:
    """Federated aggregator using explicit invertible point actions.

    The current point-valued aggregation path is scientifically supported for
    caller-supplied matrices that act invertibly on the chosen point
    representation and preserve the manifold domain used by the Karcher mean.
    The preregistered S^2 benchmark exercises this contract with explicit
    SO(3) rotations.  Arbitrary square matrices are not thereby validated as
    geometric transport morphisms.

    ``consistency_threshold`` is a threshold on the basis-dependent holonomy
    defect.  A finite threshold decision is not invariant under general
    non-orthogonal changes of frame; it must not be interpreted as a canonical
    cohomological verdict.
    """

    manifold: object
    graph: nx.DiGraph
    base_node: str
    consistency_threshold: float = 1e-6
    track_divergence: bool = False
    morphisms: dict[tuple[str, str], Morphism] = field(default_factory=dict)
    _round_idx: int = field(default=0, init=False)
    _prev_divergence: _persistence.PersistenceSummary | None = field(
        default=None, init=False, repr=False
    )

    def register_transport(
        self, source: str, target: str, matrix: npt.NDArray[np.float64]
    ) -> None:
        """Register an invertible candidate point action between two clients.

        Registration establishes only the algebraic prerequisites that can be
        checked without seeing a point: a finite square matrix with a finite
        inverse. If the opposite orientation is already registered, the two
        matrices must also satisfy the groupoid inverse law numerically. During
        aggregation the actual forward and return actions are required to map
        the transported points back onto ``self.manifold``. Passing these checks
        validates the exercised point actions, not every possible manifold point.
        """
        candidate = np.asarray(matrix, dtype=float)
        if candidate.ndim != 2 or candidate.shape[0] != candidate.shape[1]:
            raise InvalidPointTransportError(
                f"transport {source}->{target} must be a square matrix; "
                f"got shape {candidate.shape}"
            )
        if not np.all(np.isfinite(candidate)):
            raise InvalidPointTransportError(
                f"transport {source}->{target} contains non-finite values"
            )
        try:
            with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
                candidate_inverse = np.linalg.inv(candidate)
        except np.linalg.LinAlgError as exc:
            raise InvalidPointTransportError(
                f"transport {source}->{target} is singular and cannot define "
                "the inverse point action required by aggregation"
            ) from exc
        if not np.all(np.isfinite(candidate_inverse)):
            raise InvalidPointTransportError(
                f"transport {source}->{target} has a non-finite numerical inverse "
                "and cannot define the return point action required by aggregation"
            )

        if source != target:
            # For a self-loop the ``(target, source)`` key is this same arrow,
            # so an existing entry is a replacement of the same orientation
            # rather than an independently supplied reverse arrow.
            reverse = self.morphisms.get((target, source))
            if reverse is not None:
                validate_reciprocal_transports(
                    candidate,
                    reverse.transport_map,
                    source=source,
                    target=target,
                )

        self.morphisms[(source, target)] = Morphism(
            source=source,
            target=target,
            transport_map=candidate,
        )
        logger.debug("Registered candidate point action {} -> {}", source, target)

    def register_transport_from_points(
        self,
        source: str,
        target: str,
        source_point: npt.NDArray[np.float64],
        target_point: npt.NDArray[np.float64],
        method: str = "pole",
        n_rungs: int = 2,
    ) -> npt.NDArray[np.float64]:
        """Deprecated compatibility stub; tangent transport is not a point action.

        Earlier releases assembled a square ambient array from transported
        tangent basis vectors and silently registered it as an invertible
        point-valued morphism.  On an embedded manifold such as S^2, the exact
        projector extension of tangent parallel transport is rank-deficient in
        the ambient representation and sends the base point's normal direction
        to zero.  It therefore cannot satisfy the point-action and inverse
        contract used by :meth:`aggregate`.

        Use :meth:`register_transport` with an explicitly justified point action
        (for example, the SO(3) rotations used by the S^2 benchmark).  Tangent-
        vector utilities remain available in :mod:`groupoid.transport`.
        """
        raise UnsupportedTransportRepresentationError(
            "register_transport_from_points() is disabled because tangent-vector "
            "parallel transport does not by itself define the invertible point "
            "action required by this aggregator. Register an explicit, "
            "representation-correct point action instead."
        )

    def _require_manifold_point(
        self,
        point: npt.NDArray[np.float64],
        *,
        context: str,
    ) -> None:
        """Fail closed unless ``point`` belongs to the configured manifold."""
        belongs = getattr(self.manifold, "belongs", None)
        if not callable(belongs):
            raise InvalidPointTransportError(
                "the point-valued aggregation contract requires a manifold.belongs() check"
            )
        if not bool(np.all(np.asarray(belongs(point)))):
            raise InvalidPointTransportError(f"{context} is outside the configured manifold")

    def _apply_point_action(
        self,
        matrix: npt.NDArray[np.float64],
        point: npt.NDArray[np.float64],
        *,
        context: str,
    ) -> npt.NDArray[np.float64]:
        """Apply a registered linear point action and verify its exercised image."""
        try:
            with np.errstate(over="ignore", invalid="ignore"):
                mapped = matrix @ point
        except ValueError as exc:
            raise InvalidPointTransportError(
                f"{context} has incompatible matrix/point dimensions"
            ) from exc
        result = np.asarray(mapped, dtype=float)
        if not np.all(np.isfinite(result)):
            raise InvalidPointTransportError(f"{context} produced non-finite coordinates")
        self._require_manifold_point(result, context=context)
        return result

    def _get_transport_to_base(self, node: str) -> npt.NDArray[np.float64] | None:
        """Compute the composite registered point action from ``node`` to base."""
        if node == self.base_node:
            return None

        try:
            path = nx.shortest_path(self.graph.to_undirected(), node, self.base_node)
        except nx.NetworkXNoPath as exc:
            raise DisconnectedClientGraphError(
                f"client graph is disconnected: no transport path from "
                f"{node} to base {self.base_node}"
            ) from exc
        composite: Morphism | None = None

        for i in range(len(path) - 1):
            src, tgt = path[i], path[i + 1]
            if (src, tgt) in self.morphisms:
                morphism = self.morphisms[(src, tgt)]
            elif (tgt, src) in self.morphisms:
                morphism = inverse(self.morphisms[(tgt, src)])
            else:
                raise ValueError(f"No transport map for edge ({src}, {tgt})")

            composite = morphism if composite is None else compose(composite, morphism)

        return composite.transport_map if composite is not None else None

    def check_consistency(self, client_params: dict[str, npt.NDArray[np.float64]]) -> float:
        """Return the current cycle-basis holonomy defect.

        ``client_params`` is retained in the signature for API compatibility;
        the defect depends only on the graph and registered matrices.  A value
        near zero is not, by itself, a proof that the graph is connected, that
        bridge transports are present, or that the matrices define valid point
        actions.
        """
        transport_maps = {(m.source, m.target): m.transport_map for m in self.morphisms.values()}
        defect = cycle_basis_holonomy_defect(self.graph, transport_maps)
        logger.info("Cycle-basis holonomy defect = {:.2e}", defect)
        return defect

    def aggregate(
        self,
        client_params: dict[str, npt.NDArray[np.float64]],
        weights: dict[str, float] | None = None,
    ) -> FederatedRound:
        """Run one point-valued aggregation round under the explicit transport contract."""
        self._round_idx += 1
        logger.info("Starting aggregation round {}", self._round_idx)

        for node, params in client_params.items():
            self._require_manifold_point(params, context=f"client {node} input")

        defect = self.check_consistency(client_params)
        passes_threshold = defect < self.consistency_threshold

        if not passes_threshold:
            logger.warning(
                "Cycle-basis holonomy defect {:.2e} exceeds configured threshold {:.2e}; "
                "this is a representation-dependent diagnostic, not a canonical verdict",
                defect,
                self.consistency_threshold,
            )

        transported: dict[str, npt.NDArray[np.float64]] = {}
        orthogonality_residuals: dict[str, float] = {}
        for node, params in client_params.items():
            if node == self.base_node:
                transported[node] = params
                orthogonality_residuals[node] = 0.0
            else:
                # _get_transport_to_base returns None only for node ==
                # base_node, which this else branch excludes; a disconnected
                # graph raises DisconnectedClientGraphError instead. This guard
                # is therefore defensively unreachable.
                transform = self._get_transport_to_base(node)
                if transform is None:  # pragma: no cover - unreachable defensive guard (see above)
                    raise ValueError(f"No transport path from {node} to {self.base_node}")
                transported[node] = self._apply_point_action(
                    transform,
                    params,
                    context=f"forward point action {node}->{self.base_node}",
                )
                orthogonality_residuals[node] = float(
                    np.linalg.norm(transform @ transform.T - np.eye(transform.shape[0]), "fro")
                )

        nodes = sorted(transported.keys())
        param_stack = np.stack([transported[node] for node in nodes])

        divergence: _persistence.PersistenceSummary | None = None
        if self.track_divergence:
            divergence = _persistence.track_divergence(
                param_stack, previous_summary=self._prev_divergence
            )
            self._prev_divergence = divergence

        if weights is not None:
            normalized_weights = np.array([weights.get(node, 1.0) for node in nodes])
            normalized_weights = normalized_weights / normalized_weights.sum()
        else:
            normalized_weights = None

        global_params = karcher_mean(self.manifold, param_stack, weights=normalized_weights)
        self._require_manifold_point(global_params, context="Karcher mean output")

        local_updates: dict[str, npt.NDArray[np.float64]] = {}
        for node in client_params:
            if node == self.base_node:
                local_updates[node] = global_params
            else:
                # Same defensively-unreachable guard as the forward loop above:
                # None is returned only for the base node, excluded here.
                transform = self._get_transport_to_base(node)
                if transform is None:  # pragma: no cover - unreachable defensive guard (see above)
                    raise ValueError(f"No transport path from {node} to {self.base_node}")
                try:
                    with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
                        inverse_transform = np.linalg.inv(transform)
                except np.linalg.LinAlgError as exc:
                    raise InvalidPointTransportError(
                        f"composite transport for {node}->{self.base_node} became singular"
                    ) from exc
                if not np.all(np.isfinite(inverse_transform)):
                    raise InvalidPointTransportError(
                        f"return point action {self.base_node}->{node} has a non-finite "
                        "numerical inverse"
                    )
                local_updates[node] = self._apply_point_action(
                    inverse_transform,
                    global_params,
                    context=f"return point action {self.base_node}->{node}",
                )

        result = FederatedRound(
            global_params=global_params,
            local_updates=local_updates,
            h1_norm=defect,
            is_consistent=passes_threshold,
            transport_residuals=orthogonality_residuals,
            round_idx=self._round_idx,
            divergence=divergence,
        )

        logger.info(
            "Round {} complete: cycle-basis defect={:.2e}, passes_threshold={}",
            self._round_idx,
            defect,
            passes_threshold,
        )
        return result

aggregate(client_params, weights=None)

Run one point-valued aggregation round under the explicit transport contract.

Source code in groupoid/aggregation.py

def aggregate(
    self,
    client_params: dict[str, npt.NDArray[np.float64]],
    weights: dict[str, float] | None = None,
) -> FederatedRound:
    """Run one point-valued aggregation round under the explicit transport contract."""
    self._round_idx += 1
    logger.info("Starting aggregation round {}", self._round_idx)

    for node, params in client_params.items():
        self._require_manifold_point(params, context=f"client {node} input")

    defect = self.check_consistency(client_params)
    passes_threshold = defect < self.consistency_threshold

    if not passes_threshold:
        logger.warning(
            "Cycle-basis holonomy defect {:.2e} exceeds configured threshold {:.2e}; "
            "this is a representation-dependent diagnostic, not a canonical verdict",
            defect,
            self.consistency_threshold,
        )

    transported: dict[str, npt.NDArray[np.float64]] = {}
    orthogonality_residuals: dict[str, float] = {}
    for node, params in client_params.items():
        if node == self.base_node:
            transported[node] = params
            orthogonality_residuals[node] = 0.0
        else:
            # _get_transport_to_base returns None only for node ==
            # base_node, which this else branch excludes; a disconnected
            # graph raises DisconnectedClientGraphError instead. This guard
            # is therefore defensively unreachable.
            transform = self._get_transport_to_base(node)
            if transform is None:  # pragma: no cover - unreachable defensive guard (see above)
                raise ValueError(f"No transport path from {node} to {self.base_node}")
            transported[node] = self._apply_point_action(
                transform,
                params,
                context=f"forward point action {node}->{self.base_node}",
            )
            orthogonality_residuals[node] = float(
                np.linalg.norm(transform @ transform.T - np.eye(transform.shape[0]), "fro")
            )

    nodes = sorted(transported.keys())
    param_stack = np.stack([transported[node] for node in nodes])

    divergence: _persistence.PersistenceSummary | None = None
    if self.track_divergence:
        divergence = _persistence.track_divergence(
            param_stack, previous_summary=self._prev_divergence
        )
        self._prev_divergence = divergence

    if weights is not None:
        normalized_weights = np.array([weights.get(node, 1.0) for node in nodes])
        normalized_weights = normalized_weights / normalized_weights.sum()
    else:
        normalized_weights = None

    global_params = karcher_mean(self.manifold, param_stack, weights=normalized_weights)
    self._require_manifold_point(global_params, context="Karcher mean output")

    local_updates: dict[str, npt.NDArray[np.float64]] = {}
    for node in client_params:
        if node == self.base_node:
            local_updates[node] = global_params
        else:
            # Same defensively-unreachable guard as the forward loop above:
            # None is returned only for the base node, excluded here.
            transform = self._get_transport_to_base(node)
            if transform is None:  # pragma: no cover - unreachable defensive guard (see above)
                raise ValueError(f"No transport path from {node} to {self.base_node}")
            try:
                with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
                    inverse_transform = np.linalg.inv(transform)
            except np.linalg.LinAlgError as exc:
                raise InvalidPointTransportError(
                    f"composite transport for {node}->{self.base_node} became singular"
                ) from exc
            if not np.all(np.isfinite(inverse_transform)):
                raise InvalidPointTransportError(
                    f"return point action {self.base_node}->{node} has a non-finite "
                    "numerical inverse"
                )
            local_updates[node] = self._apply_point_action(
                inverse_transform,
                global_params,
                context=f"return point action {self.base_node}->{node}",
            )

    result = FederatedRound(
        global_params=global_params,
        local_updates=local_updates,
        h1_norm=defect,
        is_consistent=passes_threshold,
        transport_residuals=orthogonality_residuals,
        round_idx=self._round_idx,
        divergence=divergence,
    )

    logger.info(
        "Round {} complete: cycle-basis defect={:.2e}, passes_threshold={}",
        self._round_idx,
        defect,
        passes_threshold,
    )
    return result

check_consistency(client_params)

Return the current cycle-basis holonomy defect.

client_params is retained in the signature for API compatibility; the defect depends only on the graph and registered matrices. A value near zero is not, by itself, a proof that the graph is connected, that bridge transports are present, or that the matrices define valid point actions.

Source code in groupoid/aggregation.py

def check_consistency(self, client_params: dict[str, npt.NDArray[np.float64]]) -> float:
    """Return the current cycle-basis holonomy defect.

    ``client_params`` is retained in the signature for API compatibility;
    the defect depends only on the graph and registered matrices.  A value
    near zero is not, by itself, a proof that the graph is connected, that
    bridge transports are present, or that the matrices define valid point
    actions.
    """
    transport_maps = {(m.source, m.target): m.transport_map for m in self.morphisms.values()}
    defect = cycle_basis_holonomy_defect(self.graph, transport_maps)
    logger.info("Cycle-basis holonomy defect = {:.2e}", defect)
    return defect

register_transport(source, target, matrix)

Register an invertible candidate point action between two clients.

Registration establishes only the algebraic prerequisites that can be checked without seeing a point: a finite square matrix with a finite inverse. If the opposite orientation is already registered, the two matrices must also satisfy the groupoid inverse law numerically. During aggregation the actual forward and return actions are required to map the transported points back onto self.manifold. Passing these checks validates the exercised point actions, not every possible manifold point.

Source code in groupoid/aggregation.py

def register_transport(
    self, source: str, target: str, matrix: npt.NDArray[np.float64]
) -> None:
    """Register an invertible candidate point action between two clients.

    Registration establishes only the algebraic prerequisites that can be
    checked without seeing a point: a finite square matrix with a finite
    inverse. If the opposite orientation is already registered, the two
    matrices must also satisfy the groupoid inverse law numerically. During
    aggregation the actual forward and return actions are required to map
    the transported points back onto ``self.manifold``. Passing these checks
    validates the exercised point actions, not every possible manifold point.
    """
    candidate = np.asarray(matrix, dtype=float)
    if candidate.ndim != 2 or candidate.shape[0] != candidate.shape[1]:
        raise InvalidPointTransportError(
            f"transport {source}->{target} must be a square matrix; "
            f"got shape {candidate.shape}"
        )
    if not np.all(np.isfinite(candidate)):
        raise InvalidPointTransportError(
            f"transport {source}->{target} contains non-finite values"
        )
    try:
        with np.errstate(over="ignore", invalid="ignore", divide="ignore"):
            candidate_inverse = np.linalg.inv(candidate)
    except np.linalg.LinAlgError as exc:
        raise InvalidPointTransportError(
            f"transport {source}->{target} is singular and cannot define "
            "the inverse point action required by aggregation"
        ) from exc
    if not np.all(np.isfinite(candidate_inverse)):
        raise InvalidPointTransportError(
            f"transport {source}->{target} has a non-finite numerical inverse "
            "and cannot define the return point action required by aggregation"
        )

    if source != target:
        # For a self-loop the ``(target, source)`` key is this same arrow,
        # so an existing entry is a replacement of the same orientation
        # rather than an independently supplied reverse arrow.
        reverse = self.morphisms.get((target, source))
        if reverse is not None:
            validate_reciprocal_transports(
                candidate,
                reverse.transport_map,
                source=source,
                target=target,
            )

    self.morphisms[(source, target)] = Morphism(
        source=source,
        target=target,
        transport_map=candidate,
    )
    logger.debug("Registered candidate point action {} -> {}", source, target)

register_transport_from_points(source, target, source_point, target_point, method='pole', n_rungs=2)

Deprecated compatibility stub; tangent transport is not a point action.

Earlier releases assembled a square ambient array from transported tangent basis vectors and silently registered it as an invertible point-valued morphism. On an embedded manifold such as S^2, the exact projector extension of tangent parallel transport is rank-deficient in the ambient representation and sends the base point's normal direction to zero. It therefore cannot satisfy the point-action and inverse contract used by :meth:aggregate.

Use :meth:register_transport with an explicitly justified point action (for example, the SO(3) rotations used by the S^2 benchmark). Tangent- vector utilities remain available in :mod:groupoid.transport.

Source code in groupoid/aggregation.py

def register_transport_from_points(
    self,
    source: str,
    target: str,
    source_point: npt.NDArray[np.float64],
    target_point: npt.NDArray[np.float64],
    method: str = "pole",
    n_rungs: int = 2,
) -> npt.NDArray[np.float64]:
    """Deprecated compatibility stub; tangent transport is not a point action.

    Earlier releases assembled a square ambient array from transported
    tangent basis vectors and silently registered it as an invertible
    point-valued morphism.  On an embedded manifold such as S^2, the exact
    projector extension of tangent parallel transport is rank-deficient in
    the ambient representation and sends the base point's normal direction
    to zero.  It therefore cannot satisfy the point-action and inverse
    contract used by :meth:`aggregate`.

    Use :meth:`register_transport` with an explicitly justified point action
    (for example, the SO(3) rotations used by the S^2 benchmark).  Tangent-
    vector utilities remain available in :mod:`groupoid.transport`.
    """
    raise UnsupportedTransportRepresentationError(
        "register_transport_from_points() is disabled because tangent-vector "
        "parallel transport does not by itself define the invertible point "
        "action required by this aggregator. Register an explicit, "
        "representation-correct point action instead."
    )

groupoid.aggregation.FederatedRound dataclass

Result of a single federated aggregation round.

h1_norm, is_consistent, and transport_residuals are retained as compatibility field names. h1_norm stores the cycle-basis holonomy Frobenius defect, not a canonical H^1 norm. is_consistent means only that this representation-dependent defect is below the configured numerical threshold. transport_residuals stores ||T T^T - I||_F for the composite forward maps, so its precise meaning is an orthogonality defect.

Source code in groupoid/aggregation.py

@dataclass
class FederatedRound:
    """Result of a single federated aggregation round.

    ``h1_norm``, ``is_consistent``, and ``transport_residuals`` are retained
    as compatibility field names. ``h1_norm`` stores the cycle-basis holonomy
    Frobenius defect, not a canonical H^1 norm. ``is_consistent`` means only
    that this representation-dependent defect is below the configured numerical
    threshold. ``transport_residuals`` stores ``||T T^T - I||_F`` for the
    composite forward maps, so its precise meaning is an orthogonality defect.
    """

    global_params: npt.NDArray[np.float64]
    local_updates: dict[str, npt.NDArray[np.float64]]
    h1_norm: float
    is_consistent: bool
    transport_residuals: dict[str, float]
    round_idx: int = 0
    divergence: _persistence.PersistenceSummary | None = None

    @property
    def cycle_basis_holonomy_defect(self) -> float:
        """Primary name for the scalar stored in the legacy ``h1_norm`` field."""
        return self.h1_norm

    @property
    def passes_consistency_threshold(self) -> bool:
        """Whether the defect is below the configured representation-specific threshold."""
        return self.is_consistent

    @property
    def orthogonality_residuals(self) -> dict[str, float]:
        """Primary name for the legacy ``transport_residuals`` values."""
        return self.transport_residuals

cycle_basis_holonomy_defect property

Primary name for the scalar stored in the legacy h1_norm field.

orthogonality_residuals property

Primary name for the legacy transport_residuals values.

passes_consistency_threshold property

Whether the defect is below the configured representation-specific threshold.

groupoid.aggregation.InvalidPointTransportError

Bases: ValueError

Raised when a registered matrix cannot serve as the required point action.

Source code in groupoid/aggregation.py

class InvalidPointTransportError(ValueError):
    """Raised when a registered matrix cannot serve as the required point action."""

groupoid.aggregation.UnsupportedTransportRepresentationError

Bases: RuntimeError

Raised when tangent transport is requested as a point-valued morphism.

Source code in groupoid/aggregation.py

class UnsupportedTransportRepresentationError(RuntimeError):
    """Raised when tangent transport is requested as a point-valued morphism."""

Manifold Operations

groupoid.manifold.karcher_mean(manifold, points, weights=None, max_iter=100, tol=1e-06)

Compute the Karcher (Frechet) mean on a Riemannian manifold.

Parameters

manifold : geomstats manifold A geomstats manifold instance with a metric. points : np.ndarray Array of shape (n_points, *point_shape) on the manifold. weights : np.ndarray or None Optional weights for the mean computation. max_iter : int Maximum iterations for gradient descent. tol : float Convergence tolerance.

Returns

np.ndarray The Karcher mean point on the manifold.

Source code in groupoid/manifold.py

def karcher_mean(
    manifold: LevelSet,
    points: npt.NDArray[np.float64],
    weights: npt.NDArray[np.float64] | None = None,
    max_iter: int = 100,
    tol: float = 1e-6,
) -> npt.NDArray[np.float64]:
    """Compute the Karcher (Frechet) mean on a Riemannian manifold.

    Parameters
    ----------
    manifold : geomstats manifold
        A geomstats manifold instance with a metric.
    points : np.ndarray
        Array of shape (n_points, *point_shape) on the manifold.
    weights : np.ndarray or None
        Optional weights for the mean computation.
    max_iter : int
        Maximum iterations for gradient descent.
    tol : float
        Convergence tolerance.

    Returns
    -------
    np.ndarray
        The Karcher mean point on the manifold.
    """
    from geomstats.learning.frechet_mean import FrechetMean

    n_points = points.shape[0]
    logger.debug(
        "Computing Karcher mean of {} points on {}",
        n_points,
        type(manifold).__name__,
    )

    estimator = FrechetMean(manifold)

    # Forward convergence controls to the underlying gradient-descent
    # optimizer when the installed geomstats version exposes it. Older
    # versions without an `optimizer` attribute fall back to their defaults.
    optimizer = getattr(estimator, "optimizer", None)
    if optimizer is not None:
        if hasattr(optimizer, "max_iter"):
            optimizer.max_iter = max_iter
        if hasattr(optimizer, "epsilon"):
            optimizer.epsilon = tol

    if weights is not None:
        estimator.fit(points, weights=weights)
    else:
        estimator.fit(points)

    mean: npt.NDArray[np.float64] = estimator.estimate_
    logger.debug("Karcher mean converged")
    return mean

Groupoid

groupoid.groupoid.Morphism

Bases: BaseModel

A matrix-labelled arrow between two nodes.

This container implements algebraic composition and inversion. Its mere construction does not certify that transport_map is a geometrically valid action on a particular manifold representation. The point-valued aggregation pipeline imposes that stronger contract when matrices are registered and exercised.

Source code in groupoid/groupoid.py

class Morphism(BaseModel):
    """A matrix-labelled arrow between two nodes.

    This container implements algebraic composition and inversion.  Its mere
    construction does not certify that ``transport_map`` is a geometrically
    valid action on a particular manifold representation.  The point-valued
    aggregation pipeline imposes that stronger contract when matrices are
    registered and exercised.
    """

    model_config = ConfigDict(arbitrary_types_allowed=True)

    source: str
    target: str
    transport_map: npt.NDArray[np.float64]

    def __repr__(self) -> str:
        return f"Morphism({self.source} -> {self.target})"

    __str__ = __repr__

groupoid.groupoid.compose(f, g)

Compose two morphisms f then g.

Source code in groupoid/groupoid.py

def compose(f: Morphism, g: Morphism) -> Morphism:
    """Compose two morphisms ``f`` then ``g``."""
    if f.target != g.source:
        raise CompositionError(f"Cannot compose: {f.target} != {g.source}")

    logger.debug("Composing {} with {}", f, g)
    composed_map = g.transport_map @ f.transport_map
    return Morphism(
        source=f.source,
        target=g.target,
        transport_map=composed_map,
    )

groupoid.groupoid.inverse(f)

Return the matrix inverse of a morphism.

numpy.linalg.LinAlgError is raised if the stored matrix is singular. Geometric validity of the inverse as a manifold point action is a separate contract enforced by the aggregation layer when such an action is used.

Source code in groupoid/groupoid.py

def inverse(f: Morphism) -> Morphism:
    """Return the matrix inverse of a morphism.

    ``numpy.linalg.LinAlgError`` is raised if the stored matrix is singular.
    Geometric validity of the inverse as a manifold point action is a separate
    contract enforced by the aggregation layer when such an action is used.
    """
    logger.debug("Inverting {}", f)
    inv_map = np.linalg.inv(f.transport_map)
    return Morphism(
        source=f.target,
        target=f.source,
        transport_map=inv_map,
    )

groupoid.groupoid.NonReciprocalTransportError

Bases: ValueError

Raised when opposite registered arrows violate the groupoid inverse law.

Source code in groupoid/groupoid.py

class NonReciprocalTransportError(ValueError):
    """Raised when opposite registered arrows violate the groupoid inverse law."""

Cycle Holonomy Diagnostic

groupoid.cohomology.cycle_basis_holonomy_defect(graph, transport_maps)

Return the maximum cycle-basis holonomy Frobenius defect.

For the cycle basis B = nx.cycle_basis(graph.to_undirected()), this function computes

D_B(T) = max_{gamma in B} ||Hol_T(gamma) - I||_F.

This is a representation-dependent diagnostic, not a canonical norm on first cohomology. Its magnitude depends in general on the selected cycle basis and is not invariant under arbitrary invertible changes of frame. Orthogonal conjugation preserves the Frobenius magnitude; cycle reversal also preserves it when the holonomy itself is orthogonal.

Exact zero has a stronger meaning than the magnitude. On a connected graph, when every underlying undirected edge carries one invertible connection map (represented by one orientation or by a reciprocal pair) and every cycle emitted by :func:networkx.cycle_basis is evaluated completely, D_B(T) == 0 is equivalent to flat transport: every closed-loop holonomy is identity.

That justification is specific to the NetworkX Paton implementation this module exercises. The emitted list is not in general the fundamental-cycle basis of one fixed spanning tree; the argument instead uses the emission order, in which every emitted cycle contributes exactly one chord not present in any earlier emitted cycle, so the induced constraint system is triangular. It is not claimed for arbitrary graph-theoretic cycle bases or for future NetworkX implementations whose emission order may differ.

The equivalence does not certify graph connectedness, bridge completeness, point-action validity, or any finite numerical threshold.

Parameters

graph Directed client/transport graph. Cycle selection is performed on its undirected projection. transport_maps Maps (source, target) edge tuples to square transport matrices. Reverse traversal uses the matrix inverse. If both orientations of an underlying edge are supplied, they must be mutual numerical inverses.

Returns

float Maximum Frobenius distance ||Hol(gamma) - I||_F over the selected basis cycles. An acyclic graph returns 0.0 because it has no cycle holonomy to test.

Raises

IncompleteCocycleError If a selected basis cycle contains an edge with no transport map in either direction. NonReciprocalTransportError If both orientations of an underlying edge are supplied but do not satisfy the groupoid inverse law. numpy.linalg.LinAlgError If a reverse-oriented edge must be traversed but its registered matrix is singular.

Source code in groupoid/cohomology.py

def cycle_basis_holonomy_defect(
    graph: nx.DiGraph,
    transport_maps: dict[tuple[str, str], npt.NDArray[np.float64]],
) -> float:
    """Return the maximum cycle-basis holonomy Frobenius defect.

    For the cycle basis ``B = nx.cycle_basis(graph.to_undirected())``, this
    function computes

        D_B(T) = max_{gamma in B} ||Hol_T(gamma) - I||_F.

    This is a representation-dependent diagnostic, not a canonical norm on
    first cohomology.  Its magnitude depends in general on the selected cycle
    basis and is not invariant under arbitrary invertible changes of frame.
    Orthogonal conjugation preserves the Frobenius magnitude; cycle
    reversal also preserves it when the holonomy itself is orthogonal.

    Exact zero has a stronger meaning than the magnitude. On a connected
    graph, when every underlying undirected edge carries one invertible
    connection map (represented by one orientation or by a reciprocal pair)
    and every cycle emitted by :func:`networkx.cycle_basis` is evaluated
    completely, ``D_B(T) == 0`` is equivalent to flat transport: every
    closed-loop holonomy is identity.

    That justification is specific to the NetworkX Paton implementation this
    module exercises. The emitted list is not in general the fundamental-cycle
    basis of one fixed spanning tree; the argument instead uses the emission
    order, in which every emitted cycle contributes exactly one chord not
    present in any earlier emitted cycle, so the induced constraint system is
    triangular. It is not claimed for arbitrary graph-theoretic cycle bases or
    for future NetworkX implementations whose emission order may differ.

    The equivalence does not certify graph connectedness, bridge completeness,
    point-action validity, or any finite numerical threshold.

    Parameters
    ----------
    graph
        Directed client/transport graph.  Cycle selection is performed on its
        undirected projection.
    transport_maps
        Maps ``(source, target)`` edge tuples to square transport matrices.
        Reverse traversal uses the matrix inverse. If both orientations of an
        underlying edge are supplied, they must be mutual numerical inverses.

    Returns
    -------
    float
        Maximum Frobenius distance ``||Hol(gamma) - I||_F`` over the selected
        basis cycles.  An acyclic graph returns ``0.0`` because it has no cycle
        holonomy to test.

    Raises
    ------
    IncompleteCocycleError
        If a selected basis cycle contains an edge with no transport map in
        either direction.
    NonReciprocalTransportError
        If both orientations of an underlying edge are supplied but do not
        satisfy the groupoid inverse law.
    numpy.linalg.LinAlgError
        If a reverse-oriented edge must be traversed but its registered matrix
        is singular.
    """
    undirected = graph.to_undirected()
    for u, v in undirected.edges():
        if u == v:
            # A self-loop is one directed transport, not two opposite arrows.
            # ``(u, v)`` and ``(v, u)`` are the same key here, so applying the
            # reciprocal-pair test would spuriously demand an involution.
            continue
        if (u, v) in transport_maps and (v, u) in transport_maps:
            validate_reciprocal_transports(
                transport_maps[(u, v)],
                transport_maps[(v, u)],
                source=u,
                target=v,
            )

    cycles = nx.cycle_basis(undirected)

    if not cycles:
        logger.debug("No cycles in graph; cycle-basis holonomy defect = 0 trivially")
        return 0.0

    max_holonomy_defect = 0.0

    for cycle in cycles:
        n = len(cycle)
        edge_maps: list[npt.NDArray[np.float64]] = []
        for i in range(n):
            u = cycle[i]
            v = cycle[(i + 1) % n]

            if (u, v) in transport_maps:
                edge_maps.append(transport_maps[(u, v)])
            elif (v, u) in transport_maps:
                edge_maps.append(np.linalg.inv(transport_maps[(v, u)]))
            else:
                raise IncompleteCocycleError(
                    f"Incomplete cycle transport: no transport map for edge ({u}, {v}) "
                    f"on cycle {cycle}; holonomy is undefined. Supply the edge "
                    "map in either direction before computing the cycle-basis defect."
                )

        holonomy = edge_maps[0]
        for transport in edge_maps[1:]:
            holonomy = transport @ holonomy

        dim = holonomy.shape[0]
        deviation = float(np.linalg.norm(holonomy - np.eye(dim), ord="fro"))
        max_holonomy_defect = max(max_holonomy_defect, deviation)

    logger.debug("Cycle-basis holonomy defect = {:.6e}", max_holonomy_defect)
    return max_holonomy_defect

Deprecated compatibility alias

groupoid.cohomology.compute_h1(graph, transport_maps)

Deprecated compatibility alias for :func:cycle_basis_holonomy_defect.

Earlier GROUPOID releases called the returned scalar an H^1 or H^1 norm. The numerical value is preserved for compatibility, but that mathematical interpretation is superseded: the value is the basis-dependent cycle-holonomy defect defined by :func:cycle_basis_holonomy_defect.

Source code in groupoid/cohomology.py

def compute_h1(
    graph: nx.DiGraph,
    transport_maps: dict[tuple[str, str], npt.NDArray[np.float64]],
) -> float:
    """Deprecated compatibility alias for :func:`cycle_basis_holonomy_defect`.

    Earlier GROUPOID releases called the returned scalar an ``H^1`` or
    ``H^1 norm``.  The numerical value is preserved for compatibility, but
    that mathematical interpretation is superseded: the value is the
    basis-dependent cycle-holonomy defect defined by
    :func:`cycle_basis_holonomy_defect`.
    """
    warnings.warn(
        "compute_h1() is deprecated: it returns a cycle-basis holonomy "
        "Frobenius defect, not a canonical H^1 norm. Use "
        "cycle_basis_holonomy_defect().",
        DeprecationWarning,
        stacklevel=2,
    )
    return cycle_basis_holonomy_defect(graph, transport_maps)

Sheaf

groupoid.sheaf.Sheaf

A cellular sheaf on a graph.

Assigns vector spaces to nodes and linear restriction maps to edges.

Source code in groupoid/sheaf.py

class Sheaf:
    """A cellular sheaf on a graph.

    Assigns vector spaces to nodes and linear restriction maps to edges.
    """

    def __init__(self, graph: nx.DiGraph) -> None:
        self.graph = graph
        self._restriction_maps: dict[tuple[str, str], npt.NDArray[np.float64]] = {}
        self._sections: dict[str, npt.NDArray[np.float64]] = {}

    def set_restriction_map(
        self, source: str, target: str, matrix: npt.NDArray[np.float64]
    ) -> None:
        """Set the restriction map for an edge."""
        self._restriction_maps[(source, target)] = matrix
        logger.debug("Set restriction map {} -> {}", source, target)

    def get_restriction_map(self, source: str, target: str) -> npt.NDArray[np.float64]:
        """Get the restriction map for an edge."""
        return self._restriction_maps[(source, target)]

    def set_section(self, node: str, value: npt.NDArray[np.float64]) -> None:
        """Set a section value at a node."""
        self._sections[node] = value

    def get_section(self, node: str) -> npt.NDArray[np.float64]:
        """Get the section value at a node."""
        return self._sections[node]

    def restrict(
        self, section: npt.NDArray[np.float64], source: str, target: str
    ) -> npt.NDArray[np.float64]:
        """Apply the restriction map to a section.

        Parameters
        ----------
        section : np.ndarray
            The section value at the source node.
        source : str
            Source node identifier.
        target : str
            Target node identifier.

        Returns
        -------
        np.ndarray
            The restricted section at the target node.
        """
        R = self._restriction_maps[(source, target)]
        result: npt.NDArray[np.float64] = R @ section
        return result

    def restrict_along_path(
        self, section: npt.NDArray[np.float64], path: list[str]
    ) -> npt.NDArray[np.float64]:
        """Restrict a section along a path of nodes.

        Parameters
        ----------
        section : np.ndarray
            The section value at path[0].
        path : list[str]
            Ordered list of nodes forming a path in the graph.

        Returns
        -------
        np.ndarray
            The section restricted to path[-1].
        """
        result = section
        for i in range(len(path) - 1):
            result = self.restrict(result, path[i], path[i + 1])
        return result

get_restriction_map(source, target)

Get the restriction map for an edge.

Source code in groupoid/sheaf.py

def get_restriction_map(self, source: str, target: str) -> npt.NDArray[np.float64]:
    """Get the restriction map for an edge."""
    return self._restriction_maps[(source, target)]

get_section(node)

Get the section value at a node.

Source code in groupoid/sheaf.py

def get_section(self, node: str) -> npt.NDArray[np.float64]:
    """Get the section value at a node."""
    return self._sections[node]

restrict(section, source, target)

Apply the restriction map to a section.

Parameters

section : np.ndarray The section value at the source node. source : str Source node identifier. target : str Target node identifier.

Returns

np.ndarray The restricted section at the target node.

Source code in groupoid/sheaf.py

def restrict(
    self, section: npt.NDArray[np.float64], source: str, target: str
) -> npt.NDArray[np.float64]:
    """Apply the restriction map to a section.

    Parameters
    ----------
    section : np.ndarray
        The section value at the source node.
    source : str
        Source node identifier.
    target : str
        Target node identifier.

    Returns
    -------
    np.ndarray
        The restricted section at the target node.
    """
    R = self._restriction_maps[(source, target)]
    result: npt.NDArray[np.float64] = R @ section
    return result

restrict_along_path(section, path)

Restrict a section along a path of nodes.

Parameters

section : np.ndarray The section value at path[0]. path : list[str] Ordered list of nodes forming a path in the graph.

Returns

np.ndarray The section restricted to path[-1].

Source code in groupoid/sheaf.py

def restrict_along_path(
    self, section: npt.NDArray[np.float64], path: list[str]
) -> npt.NDArray[np.float64]:
    """Restrict a section along a path of nodes.

    Parameters
    ----------
    section : np.ndarray
        The section value at path[0].
    path : list[str]
        Ordered list of nodes forming a path in the graph.

    Returns
    -------
    np.ndarray
        The section restricted to path[-1].
    """
    result = section
    for i in range(len(path) - 1):
        result = self.restrict(result, path[i], path[i + 1])
    return result

set_restriction_map(source, target, matrix)

Set the restriction map for an edge.

Source code in groupoid/sheaf.py

def set_restriction_map(
    self, source: str, target: str, matrix: npt.NDArray[np.float64]
) -> None:
    """Set the restriction map for an edge."""
    self._restriction_maps[(source, target)] = matrix
    logger.debug("Set restriction map {} -> {}", source, target)

set_section(node, value)

Set a section value at a node.

Source code in groupoid/sheaf.py

def set_section(self, node: str, value: npt.NDArray[np.float64]) -> None:
    """Set a section value at a node."""
    self._sections[node] = value

Sheaf Laplacian

groupoid.laplacian.build_sheaf_laplacian(sheaf, stalk_dim)

Build the sheaf Laplacian matrix.

For a sheaf F on graph G with n nodes and stalk dimension d, the sheaf Laplacian is a (nd) x (nd) block matrix defined as:

L_F = delta^T @ delta

where delta is the connection coboundary. For each edge (u, v) with restriction (transport) map R = R_{uv}: stalk(u) -> stalk(v), the coboundary acts as (delta x){(u,v)} = x_v - R{uv} x_u, so L = delta^T @ delta has blocks (summed over incident edges):

L[u,u] += R_{uv}^T @ R_{uv}    (source diagonal)
L[v,v] += I                    (target diagonal)
L[u,v] += -R_{uv}^T            (off-diagonal)
L[v,u] += -R_{uv}              (off-diagonal)

L is symmetric positive semi-definite for ANY restriction maps (it is delta^T delta); its kernel is the space of transport-consistent global sections (x_v = R_{uv} x_u on every edge).

Parameters

sheaf A Sheaf instance with restriction maps set. stalk_dim Dimension of each stalk (vector space at each node).

Returns

np.ndarray The sheaf Laplacian matrix of shape (nd, nd).

Source code in groupoid/laplacian.py

def build_sheaf_laplacian(sheaf: Sheaf, stalk_dim: int) -> npt.NDArray[np.float64]:
    """Build the sheaf Laplacian matrix.

    For a sheaf F on graph G with n nodes and stalk dimension d,
    the sheaf Laplacian is a (n*d) x (n*d) block matrix defined as:

        L_F = delta^T @ delta

    where delta is the connection coboundary. For each edge (u, v) with
    restriction (transport) map R = R_{uv}: stalk(u) -> stalk(v), the
    coboundary acts as (delta x)_{(u,v)} = x_v - R_{uv} x_u, so
    L = delta^T @ delta has blocks (summed over incident edges):

        L[u,u] += R_{uv}^T @ R_{uv}    (source diagonal)
        L[v,v] += I                    (target diagonal)
        L[u,v] += -R_{uv}^T            (off-diagonal)
        L[v,u] += -R_{uv}              (off-diagonal)

    L is symmetric positive semi-definite for ANY restriction maps (it is
    delta^T delta); its kernel is the space of transport-consistent global
    sections (x_v = R_{uv} x_u on every edge).

    Parameters
    ----------
    sheaf
        A Sheaf instance with restriction maps set.
    stalk_dim
        Dimension of each stalk (vector space at each node).

    Returns
    -------
    np.ndarray
        The sheaf Laplacian matrix of shape (n*d, n*d).
    """
    nodes = sorted(sheaf.graph.nodes())
    n = len(nodes)
    node_idx = {node: i for i, node in enumerate(nodes)}
    N = n * stalk_dim

    L = np.zeros((N, N))

    for u, v in sheaf.graph.edges():
        i, j = node_idx[u], node_idx[v]
        R = sheaf.get_restriction_map(u, v)
        i_slice = slice(i * stalk_dim, (i + 1) * stalk_dim)
        j_slice = slice(j * stalk_dim, (j + 1) * stalk_dim)

        # L = delta^T delta for coboundary (delta x)_(u,v) = x_v - R_uv x_u:
        L[i_slice, i_slice] += R.T @ R  # source diagonal: R^T R
        L[j_slice, j_slice] += np.eye(stalk_dim)  # target diagonal: I
        L[i_slice, j_slice] += -R.T  # off-diagonal: -R^T
        L[j_slice, i_slice] += -R  # off-diagonal: -R

    logger.debug("Built sheaf Laplacian: {}x{} ({} nodes, stalk_dim={})", N, N, n, stalk_dim)
    return L

groupoid.laplacian.spectral_analysis(sheaf, stalk_dim, tol=1e-10)

Compute spectral decomposition of the sheaf Laplacian.

Parameters

sheaf A Sheaf instance with restriction maps. stalk_dim Dimension of each stalk. tol Tolerance for identifying zero eigenvalues.

Returns

SpectralSummary Full spectral summary including connectivity and consensus rate.

Source code in groupoid/laplacian.py

def spectral_analysis(
    sheaf: Sheaf,
    stalk_dim: int,
    tol: float = 1e-10,
) -> SpectralSummary:
    """Compute spectral decomposition of the sheaf Laplacian.

    Parameters
    ----------
    sheaf
        A Sheaf instance with restriction maps.
    stalk_dim
        Dimension of each stalk.
    tol
        Tolerance for identifying zero eigenvalues.

    Returns
    -------
    SpectralSummary
        Full spectral summary including connectivity and consensus rate.
    """
    L = build_sheaf_laplacian(sheaf, stalk_dim)
    eigenvalues, eigenvectors = np.linalg.eigh(L)

    # Sort by magnitude
    idx = np.argsort(eigenvalues)
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # Kernel dimension (number of zero eigenvalues)
    kernel_dim = int(np.sum(np.abs(eigenvalues) < tol))

    # Spectral gap and algebraic connectivity
    nonzero_eigs = eigenvalues[np.abs(eigenvalues) >= tol]
    if len(nonzero_eigs) > 0:
        algebraic_connectivity = float(nonzero_eigs[0])
        spectral_gap = float(nonzero_eigs[0])
    else:
        algebraic_connectivity = 0.0
        spectral_gap = 0.0

    # Consensus rate: exponential convergence rate of sheaf diffusion
    # x(t+1) = (I - epsilon * L) @ x(t), converges as exp(-lambda_1 * t)
    consensus_rate = algebraic_connectivity

    logger.info(
        "Spectral analysis: kernel_dim={}, spectral_gap={:.4f}, connectivity={:.4f}",
        kernel_dim,
        spectral_gap,
        algebraic_connectivity,
    )

    return SpectralSummary(
        eigenvalues=eigenvalues,
        eigenvectors=eigenvectors,
        spectral_gap=spectral_gap,
        algebraic_connectivity=algebraic_connectivity,
        kernel_dimension=kernel_dim,
        consensus_rate=consensus_rate,
    )

groupoid.laplacian.sheaf_diffusion_step(sheaf, sections, stalk_dim, step_size=0.1)

One step of sheaf diffusion (Laplacian smoothing).

Drives local sections toward global consistency by flowing along the negative gradient of the sheaf Laplacian energy:

E(x) = x^T L_F x = sum_{(i,j)} ||R_{ij} x_i - x_j||^2

Parameters

sheaf Sheaf with restriction maps. sections Current section values at each node. stalk_dim Dimension of each stalk. step_size Diffusion step size (must be < 1/lambda_max for stability).

Returns

dict[str, npt.NDArray[np.float64]] Updated section values after one diffusion step.

Source code in groupoid/laplacian.py

def sheaf_diffusion_step(
    sheaf: Sheaf,
    sections: dict[str, npt.NDArray[np.float64]],
    stalk_dim: int,
    step_size: float = 0.1,
) -> dict[str, npt.NDArray[np.float64]]:
    """One step of sheaf diffusion (Laplacian smoothing).

    Drives local sections toward global consistency by flowing along
    the negative gradient of the sheaf Laplacian energy:

        E(x) = x^T L_F x = sum_{(i,j)} ||R_{ij} x_i - x_j||^2

    Parameters
    ----------
    sheaf
        Sheaf with restriction maps.
    sections
        Current section values at each node.
    stalk_dim
        Dimension of each stalk.
    step_size
        Diffusion step size (must be < 1/lambda_max for stability).

    Returns
    -------
    dict[str, npt.NDArray[np.float64]]
        Updated section values after one diffusion step.
    """
    L = build_sheaf_laplacian(sheaf, stalk_dim)
    nodes = sorted(sheaf.graph.nodes())

    # Stack sections into vector
    x = np.concatenate([sections[n] for n in nodes])

    # Diffusion step: x' = x - step_size * L @ x
    x_new = x - step_size * L @ x

    # Unstack
    result = {}
    for idx, node in enumerate(nodes):
        result[node] = x_new[idx * stalk_dim : (idx + 1) * stalk_dim]

    return result

Tangent-Vector Parallel Transport

groupoid.transport.schild_ladder(manifold, tangent_vec, base_point, end_point, n_rungs=1)

Parallel transport a tangent vector via Schild's ladder.

This is a discrete tangent-vector approximation. On the currently tested S^2 configuration its direction is substantially coarser than pole ladder and does not converge to the analytic value as n_rungs increases; see LIMITATIONS.md for the measured behavior.

Source code in groupoid/transport.py

def schild_ladder(
    manifold: Any,
    tangent_vec: npt.NDArray[np.float64],
    base_point: npt.NDArray[np.float64],
    end_point: npt.NDArray[np.float64],
    n_rungs: int = 1,
) -> npt.NDArray[np.float64]:
    """Parallel transport a tangent vector via Schild's ladder.

    This is a discrete tangent-vector approximation.  On the currently tested
    S^2 configuration its direction is substantially coarser than pole ladder
    and does not converge to the analytic value as ``n_rungs`` increases; see
    ``LIMITATIONS.md`` for the measured behavior.
    """
    metric = manifold.metric
    direction = metric.log(end_point, base_point)

    current_base = base_point
    current_vec = tangent_vec

    for _k in range(n_rungs):
        step = direction * (1.0 / n_rungs)
        next_base = metric.exp(step, current_base)

        u = metric.exp(-current_vec, current_base)
        midpoint = metric.exp(step / 2.0, current_base)
        log_mid_u = metric.log(u, midpoint)
        u_prime = metric.exp(-log_mid_u, midpoint)
        current_vec = metric.log(u_prime, next_base)

        current_base = next_base
        direction = metric.log(end_point, current_base)

    result: npt.NDArray[np.float64] = current_vec
    return result

groupoid.transport.pole_ladder(manifold, tangent_vec, base_point, end_point, n_rungs=1)

Parallel transport a tangent vector via pole ladder.

This routine is validated only as a tangent-vector transport approximation. In the repository's S^2 validation case it closely matches analytic Levi-Civita parallel transport in direction and magnitude, with the documented small off-tangent approximation residual.

Source code in groupoid/transport.py

def pole_ladder(
    manifold: Any,
    tangent_vec: npt.NDArray[np.float64],
    base_point: npt.NDArray[np.float64],
    end_point: npt.NDArray[np.float64],
    n_rungs: int = 1,
) -> npt.NDArray[np.float64]:
    """Parallel transport a tangent vector via pole ladder.

    This routine is validated only as a tangent-vector transport approximation.
    In the repository's S^2 validation case it closely matches analytic
    Levi-Civita parallel transport in direction and magnitude, with the
    documented small off-tangent approximation residual.
    """
    metric = manifold.metric
    direction = metric.log(end_point, base_point)

    current_base = base_point
    current_vec = tangent_vec

    for _k in range(n_rungs):
        step = direction * (1.0 / n_rungs)
        next_base = metric.exp(step, current_base)

        pole = metric.exp(-current_vec, current_base)
        mid_of_geodesic = metric.exp(step / 2.0, current_base)
        log_pole_from_mid = metric.log(pole, mid_of_geodesic)
        reflected = metric.exp(-log_pole_from_mid, mid_of_geodesic)
        current_vec = metric.log(reflected, next_base)

        current_base = next_base

    return current_vec

groupoid.transport.compute_tangent_transport_matrix(manifold, base_point, end_point, method='pole', n_rungs=2)

Assemble an ambient-coordinate operator for tangent-vector transport.

This helper currently supports only vector-shaped point representations: both base_point and end_point must be one-dimensional coordinate arrays of the same shape. Each ambient coordinate vector is projected into the tangent space at base_point and transported to end_point. The transported vectors are stored as columns of a square ambient array. For matrix-valued or otherwise structured manifold points, use the ladder functions directly on tangent objects with the native point shape instead of this matrix helper.

Only the action of this array on tangent vectors is geometrically supported. It is not a generic point action and must not be registered as an invertible groupoid morphism for point-valued aggregation. On an embedded d-dimensional manifold represented in an m-dimensional ambient space with d < m, the exact projector extension has rank at most d and is therefore singular. Numerical ladder error can perturb that rank; such accidental ambient invertibility has no geometric significance.

Returns

np.ndarray Square ambient-coordinate array whose supported interpretation is the tangent-vector operator described above.

Raises

ValueError If the points are not one-dimensional coordinate arrays of the same shape.

Source code in groupoid/transport.py

def compute_tangent_transport_matrix(
    manifold: Any,
    base_point: npt.NDArray[np.float64],
    end_point: npt.NDArray[np.float64],
    method: str = "pole",
    n_rungs: int = 2,
) -> npt.NDArray[np.float64]:
    """Assemble an ambient-coordinate operator for tangent-vector transport.

    This helper currently supports only vector-shaped point representations: both
    ``base_point`` and ``end_point`` must be one-dimensional coordinate arrays of
    the same shape. Each ambient coordinate vector is projected into the tangent
    space at ``base_point`` and transported to ``end_point``. The transported
    vectors are stored as columns of a square ambient array. For matrix-valued or
    otherwise structured manifold points, use the ladder functions directly on
    tangent objects with the native point shape instead of this matrix helper.

    Only the action of this array on tangent vectors is geometrically supported.
    It is **not** a generic point action and must not be registered as an
    invertible groupoid morphism for point-valued aggregation.  On an embedded
    d-dimensional manifold represented in an m-dimensional ambient space with
    d < m, the exact projector extension has rank at most d and is therefore
    singular.  Numerical ladder error can perturb that rank; such accidental
    ambient invertibility has no geometric significance.

    Returns
    -------
    np.ndarray
        Square ambient-coordinate array whose supported interpretation is the
        tangent-vector operator described above.

    Raises
    ------
    ValueError
        If the points are not one-dimensional coordinate arrays of the same
        shape.
    """
    if base_point.ndim != 1 or end_point.ndim != 1 or base_point.shape != end_point.shape:
        raise ValueError(
            "compute_tangent_transport_matrix() supports only vector-shaped "
            "point representations: base_point and end_point must be 1D arrays "
            "with the same shape. Use the ladder functions directly for "
            "structured point representations."
        )

    transport_fn = pole_ladder if method == "pole" else schild_ladder
    dim = base_point.shape[0]
    operator = np.zeros((dim, dim))

    for i in range(dim):
        e_i = np.zeros(dim)
        e_i[i] = 1.0
        tangent = manifold.to_tangent(e_i, base_point)
        transported = transport_fn(
            manifold,
            tangent,
            base_point,
            end_point,
            n_rungs=n_rungs,
        )
        operator[:, i] = transported

    logger.debug(
        "Tangent transport ambient operator computed ({} method, {} rungs)",
        method,
        n_rungs,
    )
    return operator

Deprecated compatibility alias

groupoid.transport.compute_transport_matrix(manifold, base_point, end_point, method='pole', n_rungs=2)

Deprecated alias for :func:compute_tangent_transport_matrix.

The historical name suggested a generic invertible transport matrix. The returned array is only validated as an ambient representation of a tangent-vector transport operator.

Source code in groupoid/transport.py

def compute_transport_matrix(
    manifold: Any,
    base_point: npt.NDArray[np.float64],
    end_point: npt.NDArray[np.float64],
    method: str = "pole",
    n_rungs: int = 2,
) -> npt.NDArray[np.float64]:
    """Deprecated alias for :func:`compute_tangent_transport_matrix`.

    The historical name suggested a generic invertible transport matrix.  The
    returned array is only validated as an ambient representation of a
    tangent-vector transport operator.
    """
    warnings.warn(
        "compute_transport_matrix() is deprecated as a generic transport "
        "matrix. Use compute_tangent_transport_matrix() for tangent-vector "
        "semantics.",
        DeprecationWarning,
        stacklevel=2,
    )
    return compute_tangent_transport_matrix(
        manifold,
        base_point,
        end_point,
        method=method,
        n_rungs=n_rungs,
    )

Riemannian Optimizers

groupoid.optimizer.RiemannianSGD dataclass

Riemannian stochastic gradient descent.

Updates parameters by computing the Riemannian gradient (projection of Euclidean gradient onto tangent space) and retracting back to the manifold via the exponential map. With momentum, the velocity is parallel-transported into each new iterate's tangent space (see the module docstring), so it accumulates geometry-consistently across steps.

Parameters

manifold A geomstats manifold instance. lr Learning rate. momentum Momentum coefficient (0 = no momentum).

Source code in groupoid/optimizer.py

@dataclass
class RiemannianSGD:
    """Riemannian stochastic gradient descent.

    Updates parameters by computing the Riemannian gradient (projection
    of Euclidean gradient onto tangent space) and retracting back to
    the manifold via the exponential map. With momentum, the velocity is
    parallel-transported into each new iterate's tangent space (see the
    module docstring), so it accumulates geometry-consistently across
    steps.

    Parameters
    ----------
    manifold
        A geomstats manifold instance.
    lr
        Learning rate.
    momentum
        Momentum coefficient (0 = no momentum).
    """

    manifold: Any
    lr: float = 0.01
    momentum: float = 0.0
    _velocity: npt.NDArray[np.float64] | None = field(default=None, init=False, repr=False)

    def step(
        self, point: npt.NDArray[np.float64], euclidean_grad: npt.NDArray[np.float64]
    ) -> npt.NDArray[np.float64]:
        """Perform one optimization step.

        Parameters
        ----------
        point
            Current point on the manifold.
        euclidean_grad
            Euclidean gradient (will be projected to tangent space).

        Returns
        -------
        np.ndarray
            Updated point on the manifold.
        """
        # Project gradient onto tangent space
        riemannian_grad = self.manifold.to_tangent(euclidean_grad, point)

        # Apply momentum. _velocity was parallel-transported into `point`'s
        # tangent space at the end of the previous step; the to_tangent wrap
        # keeps the combination numerically tangent (and tolerates callers
        # stepping from a point other than the previous iterate).
        vel: npt.NDArray[np.float64]
        if self.momentum > 0:
            if self._velocity is None:
                vel = riemannian_grad
            else:
                vel = self.manifold.to_tangent(
                    self.momentum * self._velocity + riemannian_grad, point
                )
            update = -self.lr * vel
        else:
            update = -self.lr * riemannian_grad

        # Retract to manifold via exponential map
        new_point: npt.NDArray[np.float64] = self.manifold.metric.exp(update, point)

        # Carry the velocity to the new iterate's tangent space.
        if self.momentum > 0:
            self._velocity = _transport_moment(self.manifold, vel, point, new_point)

        return new_point

step(point, euclidean_grad)

Perform one optimization step.

Parameters

point Current point on the manifold. euclidean_grad Euclidean gradient (will be projected to tangent space).

Returns

np.ndarray Updated point on the manifold.

Source code in groupoid/optimizer.py

def step(
    self, point: npt.NDArray[np.float64], euclidean_grad: npt.NDArray[np.float64]
) -> npt.NDArray[np.float64]:
    """Perform one optimization step.

    Parameters
    ----------
    point
        Current point on the manifold.
    euclidean_grad
        Euclidean gradient (will be projected to tangent space).

    Returns
    -------
    np.ndarray
        Updated point on the manifold.
    """
    # Project gradient onto tangent space
    riemannian_grad = self.manifold.to_tangent(euclidean_grad, point)

    # Apply momentum. _velocity was parallel-transported into `point`'s
    # tangent space at the end of the previous step; the to_tangent wrap
    # keeps the combination numerically tangent (and tolerates callers
    # stepping from a point other than the previous iterate).
    vel: npt.NDArray[np.float64]
    if self.momentum > 0:
        if self._velocity is None:
            vel = riemannian_grad
        else:
            vel = self.manifold.to_tangent(
                self.momentum * self._velocity + riemannian_grad, point
            )
        update = -self.lr * vel
    else:
        update = -self.lr * riemannian_grad

    # Retract to manifold via exponential map
    new_point: npt.NDArray[np.float64] = self.manifold.metric.exp(update, point)

    # Carry the velocity to the new iterate's tangent space.
    if self.momentum > 0:
        self._velocity = _transport_moment(self.manifold, vel, point, new_point)

    return new_point

groupoid.optimizer.RiemannianAdam dataclass

Riemannian Adam optimizer.

Adapts the Adam optimizer to Riemannian manifolds by maintaining exponential moving averages of the Riemannian gradient and its squared norm, with updates via the exponential map. The first moment is parallel-transported into each new iterate's tangent space (see the module docstring); the second moment is a scalar gradient-norm average and needs no transport because parallel transport is an isometry (norms are invariant).

Parameters

manifold A geomstats manifold instance. lr Learning rate. beta1 Exponential decay rate for first moment. beta2 Exponential decay rate for second moment. eps Small constant for numerical stability.

Source code in groupoid/optimizer.py

@dataclass
class RiemannianAdam:
    """Riemannian Adam optimizer.

    Adapts the Adam optimizer to Riemannian manifolds by maintaining
    exponential moving averages of the Riemannian gradient and its
    squared norm, with updates via the exponential map. The first moment
    is parallel-transported into each new iterate's tangent space (see
    the module docstring); the second moment is a scalar gradient-norm
    average and needs no transport because parallel transport is an
    isometry (norms are invariant).

    Parameters
    ----------
    manifold
        A geomstats manifold instance.
    lr
        Learning rate.
    beta1
        Exponential decay rate for first moment.
    beta2
        Exponential decay rate for second moment.
    eps
        Small constant for numerical stability.
    """

    manifold: Any
    lr: float = 0.001
    beta1: float = 0.9
    beta2: float = 0.999
    eps: float = 1e-8
    _m: npt.NDArray[np.float64] | None = field(default=None, init=False, repr=False)
    _v: float = field(default=0.0, init=False, repr=False)
    _t: int = field(default=0, init=False, repr=False)

    def step(
        self, point: npt.NDArray[np.float64], euclidean_grad: npt.NDArray[np.float64]
    ) -> npt.NDArray[np.float64]:
        """Perform one optimization step.

        Parameters
        ----------
        point
            Current point on the manifold.
        euclidean_grad
            Euclidean gradient.

        Returns
        -------
        np.ndarray
            Updated point on the manifold.
        """
        self._t += 1

        # Riemannian gradient
        grad = self.manifold.to_tangent(euclidean_grad, point)
        grad_norm_sq = float(np.sum(grad**2))

        # Update biased first moment (tangent vector). Seed from zero like the
        # second moment below: m_1 = (1 - beta1) * grad, so the bias correction
        # m_hat = m_1 / (1 - beta1) recovers `grad` on the first step. Seeding
        # m_1 = grad directly would leave the 1/(1-beta1) factor uncancelled and
        # inflate the first update by ~10x (beta1=0.9).
        # _m was parallel-transported into `point`'s tangent space at the end
        # of the previous step; the to_tangent wrap keeps the combination
        # numerically tangent (and tolerates callers stepping from a point
        # other than the previous iterate).
        first_moment: npt.NDArray[np.float64]
        if self._m is None:
            first_moment = (1 - self.beta1) * grad
        else:
            first_moment = self.manifold.to_tangent(
                self.beta1 * self._m + (1 - self.beta1) * grad, point
            )

        # Update biased second moment (scalar, norm-based)
        self._v = self.beta2 * self._v + (1 - self.beta2) * grad_norm_sq

        # Bias correction
        m_hat = first_moment / (1 - self.beta1**self._t)
        v_hat = self._v / (1 - self.beta2**self._t)

        # Adaptive update
        update = -self.lr * m_hat / (np.sqrt(v_hat) + self.eps)
        update = self.manifold.to_tangent(update, point)

        new_point: npt.NDArray[np.float64] = self.manifold.metric.exp(update, point)

        # Carry the first moment to the new iterate's tangent space.
        self._m = _transport_moment(self.manifold, first_moment, point, new_point)

        return new_point

step(point, euclidean_grad)

Perform one optimization step.

Parameters

point Current point on the manifold. euclidean_grad Euclidean gradient.

Returns

np.ndarray Updated point on the manifold.

Source code in groupoid/optimizer.py

def step(
    self, point: npt.NDArray[np.float64], euclidean_grad: npt.NDArray[np.float64]
) -> npt.NDArray[np.float64]:
    """Perform one optimization step.

    Parameters
    ----------
    point
        Current point on the manifold.
    euclidean_grad
        Euclidean gradient.

    Returns
    -------
    np.ndarray
        Updated point on the manifold.
    """
    self._t += 1

    # Riemannian gradient
    grad = self.manifold.to_tangent(euclidean_grad, point)
    grad_norm_sq = float(np.sum(grad**2))

    # Update biased first moment (tangent vector). Seed from zero like the
    # second moment below: m_1 = (1 - beta1) * grad, so the bias correction
    # m_hat = m_1 / (1 - beta1) recovers `grad` on the first step. Seeding
    # m_1 = grad directly would leave the 1/(1-beta1) factor uncancelled and
    # inflate the first update by ~10x (beta1=0.9).
    # _m was parallel-transported into `point`'s tangent space at the end
    # of the previous step; the to_tangent wrap keeps the combination
    # numerically tangent (and tolerates callers stepping from a point
    # other than the previous iterate).
    first_moment: npt.NDArray[np.float64]
    if self._m is None:
        first_moment = (1 - self.beta1) * grad
    else:
        first_moment = self.manifold.to_tangent(
            self.beta1 * self._m + (1 - self.beta1) * grad, point
        )

    # Update biased second moment (scalar, norm-based)
    self._v = self.beta2 * self._v + (1 - self.beta2) * grad_norm_sq

    # Bias correction
    m_hat = first_moment / (1 - self.beta1**self._t)
    v_hat = self._v / (1 - self.beta2**self._t)

    # Adaptive update
    update = -self.lr * m_hat / (np.sqrt(v_hat) + self.eps)
    update = self.manifold.to_tangent(update, point)

    new_point: npt.NDArray[np.float64] = self.manifold.metric.exp(update, point)

    # Carry the first moment to the new iterate's tangent space.
    self._m = _transport_moment(self.manifold, first_moment, point, new_point)

    return new_point

groupoid.optimizer.curvature_adaptive_lr(manifold, point, base_lr, tangent_vec)

Adapt learning rate based on local sectional curvature.

In regions of high positive curvature, geodesics converge and we should take smaller steps, so the base rate is damped by 1 / (1 + kappa). In flat or negatively curved regions the base rate is returned unchanged; no enlargement is applied.

Parameters

manifold A geomstats manifold with a curvature method. point Current point on the manifold. base_lr Base learning rate to adapt. tangent_vec Direction of the update.

Returns

float Adapted learning rate.

Source code in groupoid/optimizer.py

def curvature_adaptive_lr(
    manifold: Any,  # geomstats manifold; no upstream type stubs
    point: npt.NDArray[np.float64],
    base_lr: float,
    tangent_vec: npt.NDArray[np.float64],
) -> float:
    """Adapt learning rate based on local sectional curvature.

    In regions of high positive curvature, geodesics converge and we
    should take smaller steps, so the base rate is damped by
    ``1 / (1 + kappa)``. In flat or negatively curved regions the base
    rate is returned unchanged; no enlargement is applied.

    Parameters
    ----------
    manifold
        A geomstats manifold with a curvature method.
    point
        Current point on the manifold.
    base_lr
        Base learning rate to adapt.
    tangent_vec
        Direction of the update.

    Returns
    -------
    float
        Adapted learning rate.
    """
    try:
        if hasattr(manifold.metric, "sectional_curvature"):
            # Use a random orthogonal vector for the plane
            random_vec = manifold.to_tangent(np.random.randn(*tangent_vec.shape), point)
            kappa = manifold.metric.sectional_curvature(tangent_vec, random_vec, point)
            kappa = float(np.mean(kappa)) if hasattr(kappa, "__len__") else float(kappa)

            # Scale: lr / (1 + max(kappa, 0)) damps in positive curvature
            adapted: float = base_lr / (1.0 + max(kappa, 0.0))
            logger.debug("Curvature-adapted LR: {:.6f} (kappa={:.4f})", adapted, kappa)
            return adapted
    except (AttributeError, NotImplementedError):
        pass

    return base_lr

Persistent Homology

groupoid.persistence.compute_persistence(points, max_dim=1, max_edge_length=np.inf)

Compute persistent homology of a point cloud.

Parameters

points Array of shape (n_points, n_features) representing model parameters or their embeddings. max_dim Maximum homological dimension to compute. max_edge_length Maximum edge length for the Rips filtration.

Returns

PersistenceSummary Topological summary including Betti numbers and persistence.

Source code in groupoid/persistence.py

def compute_persistence(
    points: npt.NDArray[np.float64],
    max_dim: int = 1,
    max_edge_length: float = np.inf,
) -> PersistenceSummary:
    """Compute persistent homology of a point cloud.

    Parameters
    ----------
    points
        Array of shape (n_points, n_features) representing model
        parameters or their embeddings.
    max_dim
        Maximum homological dimension to compute.
    max_edge_length
        Maximum edge length for the Rips filtration.

    Returns
    -------
    PersistenceSummary
        Topological summary including Betti numbers and persistence.
    """
    from ripser import ripser

    logger.debug("Computing persistence: {} points, max_dim={}", points.shape[0], max_dim)

    result = ripser(points, maxdim=max_dim, thresh=max_edge_length)
    diagrams = result["dgms"]

    # Compute Betti numbers (count features that persist at infinity)
    betti_0 = int(np.sum(diagrams[0][:, 1] == np.inf)) if len(diagrams) > 0 else 0
    betti_1 = int(np.sum(diagrams[1][:, 1] == np.inf)) if len(diagrams) > 1 else 0

    # Compute total and max persistence (excluding infinite features)
    all_finite = []
    for dgm in diagrams:
        finite_mask = dgm[:, 1] < np.inf
        if np.any(finite_mask):
            lifetimes = dgm[finite_mask, 1] - dgm[finite_mask, 0]
            all_finite.extend(lifetimes.tolist())

    total_persistence = float(sum(all_finite)) if all_finite else 0.0
    max_persistence = float(max(all_finite)) if all_finite else 0.0

    # Concatenate all diagrams for storage, RETAINING the homology
    # dimension as a third column. ripser returns a per-dimension list
    # (diagrams[d] holds dimension d); a bare np.vstack would discard
    # that label and make H0 and H1 bars indistinguishable. We append a
    # dim column so downstream consumers (e.g. track_divergence) can
    # compare like dimensions instead of pooling them.
    labelled = []
    for dim, dgm in enumerate(diagrams):
        if dgm.shape[0] == 0:
            continue
        dim_col = np.full((dgm.shape[0], 1), float(dim))
        labelled.append(np.hstack([dgm, dim_col]))
    full_diagram = np.vstack(labelled) if labelled else np.empty((0, 3))

    logger.debug(
        "Persistence: beta_0={}, beta_1={}, total={:.4f}",
        betti_0,
        betti_1,
        total_persistence,
    )

    return PersistenceSummary(
        betti_0=betti_0,
        betti_1=betti_1,
        total_persistence=total_persistence,
        max_persistence=max_persistence,
        diagram=full_diagram,
    )

groupoid.persistence.track_divergence(current_params, previous_summary=None, max_dim=1)

Track federation divergence across rounds.

Computes persistence of current parameter distribution and, if a previous summary exists, computes the bottleneck distance to measure how much the topological structure has changed.

Parameters

current_params Array of shape (n_clients, n_features). previous_summary PersistenceSummary from the previous round, if available. max_dim Maximum homological dimension.

Returns

PersistenceSummary Updated summary with bottleneck distance to previous round.

Source code in groupoid/persistence.py

def track_divergence(
    current_params: npt.NDArray[np.float64],
    previous_summary: PersistenceSummary | None = None,
    max_dim: int = 1,
) -> PersistenceSummary:
    """Track federation divergence across rounds.

    Computes persistence of current parameter distribution and, if a
    previous summary exists, computes the bottleneck distance to
    measure how much the topological structure has changed.

    Parameters
    ----------
    current_params
        Array of shape (n_clients, n_features).
    previous_summary
        PersistenceSummary from the previous round, if available.
    max_dim
        Maximum homological dimension.

    Returns
    -------
    PersistenceSummary
        Updated summary with bottleneck distance to previous round.
    """
    summary = compute_persistence(current_params, max_dim=max_dim)

    if previous_summary is not None:
        from persim import bottleneck

        # Compare H0 diagrams (connected components) ONLY. We select bars
        # whose homology dimension is 0 and which are finite (the single
        # infinite H0 bar -- the whole space's surviving component -- is
        # excluded so persim never has to match infinities). Selecting on
        # the dimension label is essential: without it, H1 (loop) bars
        # leak into this pool and silently contaminate the "H0 divergence"
        # with loop-structure changes. See PersistenceSummary.diagram_for_dim.
        current_h0 = summary.diagram_for_dim(0, finite_only=True)
        prev_h0 = previous_summary.diagram_for_dim(0, finite_only=True)

        if len(current_h0) > 0 and len(prev_h0) > 0:
            dist = bottleneck(current_h0, prev_h0)
            summary.bottleneck_to_previous = float(dist)
            logger.info("Bottleneck distance to previous round: {:.4f}", dist)
        else:
            summary.bottleneck_to_previous = 0.0

    return summary