# Asymptotically Normal Estimation of Local Latent Network Curvature

Steven Wilkins-Reeves<EMAIL_ADDRESS>Department of Statistics, University of
Washington  Tyler McCormick<EMAIL_ADDRESS>Department of Statistics,
University of Washington  Department of Sociology, University of Washington

###### Abstract

Network data, commonly used throughout the physical, social, and biological
sciences, consists of nodes (individuals) and the edges (interactions) between
them. One way to represent network data’s complex, high-dimensional structure
is to embed the graph into a low-dimensional geometric space. The curvature of
this space, in particular, provides insights about the structure in the graph,
such as the propensity to form triangles or present tree-like structures. We
derive an estimating function for curvature based on triangle side lengths and
the length of the midpoint of a side to the opposing corner. We construct an
estimator where the only input is a distance matrix and also establish
asymptotic normality. We next introduce a novel latent distance matrix
estimator for networks and an efficient algorithm to compute the estimate via
solving iterative quadratic programs. We apply this method to the Los Alamos
National Laboratory Unified Network and Host dataset and show how curvature
estimates can be used to detect a red-team attack faster than naive methods,
as well as discover non-constant latent curvature in co-authorship networks in
physics. The code for this paper is available at
https://github.com/SteveJWR/netcurve, and the methods are implemented in the R
package https://github.com/SteveJWR/lolaR.

##  1 Introduction

Social networks consist of a set of actors and the connections between them.
Formally, a set of nodes (or vertices) $V=\\{1,2,\dots,n\\}$ and edges
$E\subset\\{(i,j)|i\not=j,i\in V,j\in V\\}$ define a graph or network
$G=(V,E)$. Beyond just the social network setting where these actors may
represent individuals (e.g. Borgatti et al. (2009)), network models are used
in many areas of the biological and physical sciences where the nodes may
represent objects of interest such as cells (e.g. Bassett et al. (2018)), and
particles (e.g. Papadopoulos et al. (2018)) respectively.

The high-dimensional, complex structure is an inherent feature of network
data. This structure creates challenges in representing and modeling these
rich data. One common approach uses embedding into lower-dimensional geometric
spaces, where both properties of the geometric space and the position of
points in the space reveal insights about the structure in the graph.

In this paper, we focus on the properties of the geometric space, particularly
the notion of curvature. The sectional curvature of a latent space is, broadly
defined, the deviation from a flat (Euclidean) space via the growth of the
circumference of small circles as a function of their radius. Numerous
previous approaches in this area have adopted a “discrete” approach, where the
discreteness arises from the structure of the distance matrix between nodes
(using path length). These curvature definitions include Ollivier-Ricci
Curvature (Ollivier, 2007), Haantjes-Ricci curvature (Saucan et al., 2020) and
Forman-Ricci curvature (Leal et al., 2018).

Despite some of their unknown statistical properties, these ideas have been
useful in a variety of applications such as understanding financial network
instability (Sandhu et al., 2016; Samal et al., 2021), network sampling
(Barkanass et al., 2022), cancer detection in gene regulatory networks (Sandhu
et al., 2015), functional neuroscience (Farooq et al., 2019) and community
detection (Sia et al., 2019; Ni et al., 2019). These definitions are motivated
for networks defined as a static object and do not necessarily measure the
curvature of the underlying latent space. Therefore it is not apparent what
these estimates will converge to (or if they converge at all) when a network
is studied as a random object.

Our approach deviates from these previous approaches to measuring curvature by
developing a method that is explicitly linked to a particular model for
network formation, the latent distance model of Hoff et al. (2002). Defining
curvature based on a specific probabilistic model yields a direct connection
between curvature and the underlying properties of a generative model for the
graph. Specifically, we develop a consistent and asymptotically normal
estimator of the local curvature of the latent distance model, though our
results for curvature estimation apply broadly to settings where distance
matrices can be estimated. We also demonstrate that this estimator provides
valuable insights when used in several existing statistical problems for
network data. Specifically, we develop a hypothesis testing framework to
reveal areas of non-constant curvature, which correspond to notions of
brokerage in the sociology literature (Burt, 1992; Buskens and Van de Rijt,
2008). We also demonstrate that our estimator can be used in changepoint
detection with an application to cybersecurity.

Moving now to discuss the probabilistic framework in more detail, the latent
distance model consists of locations $Z_{i},Z_{j}$, most generally on some
metric space, $\mathfrak{M}$, often taken to be a smooth Riemannian manifold
$\mathcal{M}^{p}$ where the probability of forming an edge $(i,j)\in E$ is
inversely proportional to the distance on the latent manifold,
$d_{\mathcal{M}}(Z_{i},Z_{j})$,

$\widetilde{g}\Biggl{(}P\biggl{(}(i,j)\in E\biggr{)}\Biggr{)}\propto-
d_{\mathcal{M}}(Z_{i},Z_{j})$

where $\widetilde{g}(\cdot)$ is some link function. The latent distance models
and their extensions have been for such problems as modeling social influence
(Sweet and Adhikari, 2020), social media relationships of politicians (Lok et
al., 2021) and neuron connectivity (Aliverti and Durante, 2019) among numerous
others. Where convenient, we condense $d_{\mathcal{M}}(x,y)$ with $d(x,y)$ or
$d_{xy}$ for brevity of notation.

The true latent geometry has implications for the global network properties.
For example, spherical (positively curved) spaces form triangles and clusters
more often than flat spaces, and hyperbolic (negatively curved) spaces tend to
generate tree-like structures (Smith et al., 2017). In other generative models
of networks, the degree distribution and clustering of the network have
implications of epidemic severity in network-based SIR models (Volz et al.,
2011), suggesting links to other properties of interest among network data.

Lubold et al. (2023) constructed a method of identifying the latent curvature.
However, this method requires first picking a sign of the curvature prior to
estimation, relies on constant curvature globally, and does not admit a known
asymptotic result. One study of discrete graph curvature from the physics
literature illustrates consistency to a global manifold in a particular data-
generating case (van der Hoorn et al., 2020). The authors study the
convergence of a modified Ollivier-Ricci curvature to the Ricci curvature of
the underlying space in random geometric graphs under the limit of the
connection radius shrinking to $0$. Smith et al. (2017) also provide a method
for identifying curvature based on simulations that compare the eigenspectrum
of the graph Laplacian to models under spherical, hyperbolic, and Euclidean
geometry.

The remainder of the paper continues as follows. First, in Section 2, we
introduce the statistical model and manifolds of constant curvature. We then
introduce the midpoint curvature equation, which allows for the identification
of curvature given a set of distances. We follow up with some asymptotic
theory when distances are estimated and an exact midpoint is present. We then
establish convergence rates for midpoint set formation. In Section 4, we
consider the practical estimation methods and illustrate consistency of our
estimator.

Next, we discuss downstream statistical tasks such as a testing constant
curvature in Section 5 and detecting changepoints in Section 6. We further
elaborate on these with applications of these methods to co-authorship
networks in physics and an application in cybersecurity.

##  2 Methods

We will begin with an introduction of the models under consideration, then
follow up with an introduction of the estimator and the associated asymptotic
theory. We follow this up with a discussion of the bias associated with a
point that is not precisely the midpoint of two other points and discuss the
rate of convergence of the best approximation to a midpoint of two other
points. Lastly, we discuss the estimation of distance matrices in the
generative network model.

###  2.1 Generative Model

We consider the random symmetric network corresponding to a graph $G=(V,E)$
where $|V|=n$, with adjacency matrix $A\in\\{0,1\\}^{n\times n}$ such that
$A_{ij}=1$ iff $(i,j)\in E$. This random adjacency matrix is generated under a
latent distance model (Hoff et al., 2002) with random effects (gregariousness
parameters)

$P(A_{ij}=1|\nu,Z)=\exp\biggl{(}\nu_{i}+\nu_{j}-d_{\mathcal{M}}(Z_{i},Z_{j})\biggr{)}$
(1)

where

$\displaystyle Z_{i}$ $\displaystyle\sim_{iid}F_{Z}$ $\displaystyle\nu_{i}$
$\displaystyle\sim_{iid}F_{\nu}$

with random effect measure ($F_{\nu}$) with support on the non-positive real
line $\text{supp}(F_{\nu})\subset(\infty,0]$ and where $iid$ refers to
sampling identically and independently from the latent distribution. This
model relates the latent distances of the generative model to the probability
of a connection between two points. Though not discussed in this paper, simple
generalizations can be derived for directed networks. The latent position
measure ($F_{Z}$) has support on some unknown latent manifold
$\mathcal{M}^{p}$ (of dimension $p$). This paper aims to estimate the
curvature of $\mathcal{M}^{p}$. The $\exp(\cdot)$ link function is used as an
example for ease of explanation, however in the supplementary materials we
discuss the extension to the $\operatorname*{\text{expit}}(\cdot)$ model as
well as possible further extensions.

Latent distance models typically assume a class of simply connected Riemannian
manifolds of constant sectional curvature $(\kappa)$, classically, the
Euclidean $\mathbb{E}^{p},(\kappa=0)$ and spherical
$\mathbb{S}^{p}(\kappa),(\kappa>0)$ and more recently the hyperbolic space
$\mathbb{H}^{p}(\kappa)$ ($\kappa<0$). However there are some extensions which
are not restricted to this structure (Fosdick et al., 2016; Lok et al., 2021).
By the Killing-Hopf theorem (Killing, 1891), these are the only manifolds of
constant sectional curvature. We will refer to these latent spaces as the
canonical manifolds. These are selected for two reasons. Firstly there are
techniques for visualization for each of these methods (for $p=2$), thus
highlighting the roles of individuals in the estimated model, and secondly,
closed form expressions are available for the distances, allowing for
estimation of these models. Embedding latent positions in a general Riemannian
manifold would require the solving numerous variational optimization problems
(computing geodesics/distances), preventing their use in practice. For a
broader introduction on Riemannian geometry see Klingenberg (1995).

###  2.2 Models of Canonical Manifolds

Each of these canonical manifolds can be represented using a set of positions
with real-valued vectors and a corresponding distance function. We include
definitions for Euclidean, spherical, and hyperbolic spaces for completeness.

The Euclidean manifold will be equipped with the standard $\ell_{2}$ norm
$\mathbb{E}^{p}$ which describes the set of points for which
$x,y\in\mathbb{R}^{p}$ with distance matrix

$\displaystyle d_{\mathbb{E}^{p}}(x,y)$
$\displaystyle=\sqrt{\sum_{k=1}^{p}(x_{k}-y_{k})^{2}}$

The spherical model with curvature $\kappa>0$ is equivalent to the sphere of
radius $r=\frac{1}{\kappa^{2}}$. A simple way of expressing this is via the
following model. A set of points are embedded in $\mathbb{S}^{p}$ if
$x\in\mathbb{R}^{p+1}$, $\sum_{k=0}^{p}x_{i}^{2}=1$. We add a zero indexed
coordinate for the spherical and hyperbolic models. We compute a signature
$B_{\mathbb{S}^{p}}$ and distances are computed via the signature

$\displaystyle B_{\mathbb{S}^{p}}(x,y)$
$\displaystyle=\sum_{k=0}^{p}x_{i}y_{i}$ $\displaystyle
d_{\mathbb{S}^{p}}(x,y)$
$\displaystyle=\frac{1}{\sqrt{\kappa}}\mathrm{acos}(B_{\mathbb{S}^{p}}(x,y))$

The hyperboloid model with curvature $\kappa<0$ corresponds to a set of points
in $\mathbb{H}^{p}$ if $x\in\mathbb{R}^{p+1}$,
$x_{0}^{2}-\sum_{k=1}^{p}x_{i}^{2}=1,x_{0}>0$. An analogous signature
$B_{\mathbb{H}^{p}}$ exists for the hyperbolic embedding.

$\displaystyle B_{\mathbb{H}^{p}}(x,y)$
$\displaystyle=x_{0}y_{0}-\sum_{k=1}^{p}x_{i}y_{i}$ $\displaystyle
d_{\mathbb{H}^{p}}(x,y)$
$\displaystyle=\frac{1}{\sqrt{-\kappa}}\mathrm{acosh}(B_{\mathbb{H}^{p}}(x,y))$

In our setting, we do not require $\mathcal{M}^{p}$ to have constant
curvature, but we do make some much more mild assumptions on
$\mathcal{M}^{P}$. We consider $\mathfrak{M}=(\mathcal{M}^{p},d)$ the metric
space induced by the Riemannian manifold.

  * •

(A1) (Algebraic Midpoint Property) $\mathfrak{M}$ satisfies the algebraic
midpoint property. For any $x,y\in\mathcal{M}^{p}$ there exists a point $z$
such that
$d_{\mathcal{M}^{p}}(z,x)=d_{\mathcal{M}^{p}}(z,y)=\frac{1}{2}d_{\mathcal{M}^{p}}(y,z)$.

  * •

(A2) (Locally Euclidean) For all $p\in\mathcal{M}^{p}$ there exists some
$\delta>0$ and some $c_{p},C_{p}$ such that for all $\epsilon\leq\delta$:

$c_{p}(\delta)\leq\frac{\text{Vol}(B_{\mathcal{M}^{p}}(\epsilon,q))}{\text{Vol}(B_{\mathbb{E}^{p}}(\epsilon,0))}\leq
C_{p}(\delta)$

and

$\lim_{\epsilon\to
0}\frac{\text{Vol}(B_{\mathcal{M}^{p}}(\epsilon,q))}{\text{Vol}(B_{\mathbb{E}^{p}}(\epsilon,0))}=1$

where $B_{\mathcal{M}^{p}}(\epsilon,q)$ is the $\epsilon$ ball on
$\mathcal{M}^{p}$ centred at a point $q$. Unless otherwise specified, from
here if denoting $B(\epsilon,q)$ we are referring to the ball on the latent
manifold $\mathcal{M}^{p}$.

If $\mathcal{M}^{p}$ is a connected manifold with sectional curvature upper
and lower bounded, then these properties are satisfied. See Section E in the
supplementary materials for details.

###  2.3 Curvature and Estimating Equation

A number of methods exist to verify whether a particular set of distances can
be embedded in a space of constant curvature. These include Schoenberg (1935)
and Begelfor and Werman (2005) as well as Cayley-Menger determinants
(Blumenthal and Gillam, 1943). Lubold et al. (2023) previously used these
methods in order identify whether a set of distances could be globally
embedded in a particular curvature space. In this work, we take a different
approach where we can identify curvature based on a minimal set of points.

We rely on a simple geometric insight. Consider a set of three points that
form a triangle. If we look at the midpoint of any of the sides of this
triangle, we can construct a distance from the midpoint to the opposite
vertex. This length reveals the sectional curvature of the manifold, where a
smaller distance corresponds to a more negatively curved space and a larger
distance corresponding to a more positively curved space. This is visualized
in Figure 1 for Euclidean, spherical and hyperbolic triangles.

Figure 1: Midpoint distances and curvature of the space. The distance $d_{xm}$
is an increasing function of the curvature $\kappa$ for fixed other triangle
side lengths.

This use of midpoints will help us identify the curvature of the latent space.
For example, consider a case where no midpoint is present in the data. We let
$4$ points placed equidistant from each other, we can either embed these in
$\mathbb{E}^{3}$ or $\mathbb{S}^{2}(\kappa)$ for $\kappa>0$. Therefore the
curvature of the latent space cannot be identified from this configuration.
Though conditions like Schoenberg’s, (Schoenberg, 1935) which are used in
Lubold et al. (2023) can determine whether a distance matrix can be embedded
in a space of constant curvature, identifiability concerns remain.

We instead take a much more direct approach at eliciting curvature, through
the use of the existence of a midpoint. This lets us identify a curvature with
only $4$ points in total. The use of a midpoint allows us to work with a sub-
manifold of dimension $2$ which as long as $p\geq 2$ will allow for
identification of the curvature, circumventing the need to identify the
dimensionality of the latent space. We will introduce an estimator for the
curvature of the latent space and derive its asymptotic distribution in the
case of a known midpoint. We follow this up with a discussion regarding the
convergence rate of surrogate (approximate) midpoints.

We now formalize this intuition. For any $3$ points $x,y,z$ which lie in an
unknown Riemannian manifold of dimension $p\geq 2$ of constant sectional
curvature $\kappa$, $(x,y,z)$ can be isometrically embedded in a sub-manifold
of dimension $2$.

######  Definition 1.

A submanifold $\widetilde{\mathcal{M}}\subset\mathcal{M}$ is totally geodesic
if every geodesic in $\widetilde{\mathcal{M}}$ is also a geodesic in
$\mathcal{M}$.

We will use the fact that a totally geodesic submanifold contains all points
along the geodesic, including the midpoint.

######  Lemma 1.

For any $x,y,z\in\mathcal{M}^{p}(\kappa)$ if $m$ is a midpoint lying between
$y$ and $z$. Then $x,y,z,m\in\mathcal{M}^{2}(\kappa)$ where
$\mathcal{M}^{2}(\kappa)$ is a totally geodesic submanifold of dimension 2.

The proof is straightforward and in Section A.1. The main implication here is
that the totally geodesic submanifold allows us to look at distances in a
subspace of dimension $2$ which allow for identification. The three points
$x,y,z$ will fall into one of these sub-manifolds, $m$ which lies on the
geodesic between $y$ and $z$. Since geodesics determine the distance, and
geodesics on the submanifold are the same as geodesics on the total manifold,
then this gives a criteria for identifying the latent curvature through the
fortunate alignment of midpoints.

We now use this fact to derive an equation which will relate the curvature
$\kappa$ to the set of distances between the points $(x,y,z,m)$.

######  Theorem 2 (Midpoint Curvature Equation).

Suppose that points $x,y,z\in\mathcal{M}^{p}(\kappa)$ an unknown Riemannian
manifold of dimension $p$ of constant sectional curvature $\kappa$. Let $m$
denote the midpoint between $y,z$. The following equation holds for
$\kappa\in\mathbb{R}$.

$\displaystyle g(\kappa,d)$
$\displaystyle=\text{Re}\bigg{[}\frac{2\cos(d_{xm}\sqrt{\kappa})}{\kappa}$ (2)
$\displaystyle-\frac{\mathrm{sec}(\frac{d_{yz}}{2}\sqrt{\kappa})(\cos(d_{xy}\sqrt{\kappa})+\cos(d_{xz}\sqrt{\kappa}))}{\kappa}\bigg{]}=0$

Where $d_{jk}$ denoted the distance between points $j,k$ and $\text{Re}[]$
denotes the real part of the equation.

For cases when $\kappa<0$, we take the real part of the above equation, which
is equivalent to replacing the trigonometric functions with their hyperbolic
analogues. Naturally, this can also be solved for $d_{xm}$ to compute the
distance to the midpoint in a space of constant curvature
$d_{xm}(\kappa^{\prime};\boldsymbol{d}(x,y,z))$.

The proof is found in the supplementary materials in Section A.2. In essence,
this equation is related to Toponogov’s theorem (Klingenberg (1995) 2.7.12),
which can be used to define a parameter based on distances $\theta_{T}$ based
on distances of a triangle $x,y,z$.

$\displaystyle\theta_{T}(d)$
$\displaystyle:=\frac{1}{2}d^{2}_{xz}+\frac{1}{2}d^{2}_{xy}+\frac{1}{4}d^{2}_{yz}-d^{2}_{xm}.$

If the sectional curvature of $\mathcal{M}^{p}$ is non-negative, then
$\theta_{T}(d)\geq 0$ and if the curvature is non-positive $\theta_{T}(d)\geq
0$. In Euclidean space, this reduces to the parallelogram law. Toponogov’s
theorem itself does not identify curvature, but if we assume that the points
$x,y,z,m\in\mathcal{M}^{2}(\kappa)\subseteq\mathcal{M}^{p}$ a totally geodesic
submanifold of dimension $2$ which is simply connected and of constant
curvature, then this set of distances can be used to identify the curvature.
The important distinction between this method of identification via Schoenberg
(1935) or through Cayley-Menger determinants (Blumenthal and Gillam, 1943) is
that these do not include the information of a midpoint whereas ours does.
Additionally, these embedding theorems do not necessarily allow for the unique
identification of curvature, but rather if a particular space of constant
curvature can embed a given distance matrix.

To highlight the smoothness of our estimating function we plot a set of
examples. For a unit equilateral triangle, we compute the corresponding
midpoint distance $d_{xm}$ for each curvature space with
$\kappa\in\\{-2,-1,0,1,2\\}$. We see in Figure 2 that our proposed estimating
equation is differentiable around the solution with non-zero derivative,
allowing one to identify the curvature from the $\kappa:g(\kappa,d)=0$.

Figure 2: Example of estimating functions $g$ as a function of $\kappa$.

###  2.4 Curvature Estimation

Now that we have an estimating function which allows us to identify curvature,
we discuss the estimator’s asymptotics.

######  Theorem 3.

Suppose there exist points $x,y,z\in\mathcal{M}^{p}$. Let $m$ denote the
midpoint between $y,z$ where these points are fixed. Let
$\boldsymbol{\widehat{d}}=(\widehat{d}_{xy},\widehat{d}_{xz},\widehat{d}_{yz},\widehat{d}_{xm})$
be the estimated distances and $\boldsymbol{d}=(d_{xy},d_{xz},d_{yz},d_{xm})$
be their true, unknown counterparts. Assume we have a distance estimator such
$\boldsymbol{\widehat{d}}$ and

  * •

(B1) $\sqrt{r(n)}\Big{(}\boldsymbol{\widehat{d}}-\boldsymbol{d}\Big{)}\to
N(0,\Sigma)$

  * •

(B2)
$\kappa<\left(\frac{\pi}{\max\\{d_{xy},d_{xz},d_{yz},d_{xm}\\}}\right)^{2}$

  * •

(B3)
$\frac{d}{d\kappa^{\prime}}d_{xm}(\kappa^{\prime};\boldsymbol{d}(x,y,z))\bigg{|}_{\kappa^{\prime}=\kappa}>0$

Where $r(n)$ is the rate of convergence. Let $\widehat{\kappa}$ to the
solution to of $g(\kappa,\widehat{d})=0$. Then

$\sqrt{r(n)}(\widehat{\kappa}-\kappa)\to_{d}N\left(0,\left(\frac{\partial
g(\kappa,d)}{\partial\kappa}\right)^{-2}\left(\nabla_{d}g(\kappa,d)^{\top}\Sigma\nabla_{d}g(\kappa,d)\right)\right)$
(3)

where $\to_{d}$ refers to convergence in distribution. If $(B1)$ is replaced
by consistency, i.e.
$||{\Big{(}\boldsymbol{\widehat{d}}-\boldsymbol{d}\Big{)}}||_{2}=o_{P}(r(n)^{-1/2})$,
then $\widehat{\kappa}-\kappa=o_{P}(r(n)^{-1/2})$.

The proof is found in the supplementary materials in Section A.3 and is an
application of the implicit function theorem and the delta method. Assumption
(B1) is mild as it only requires asymptotic normality of the distance
estimator. Assumption (B2) is mild as it requires all of the distances to be
embeddable on the sphere of curvature $\kappa$. Assumption (B3) tends to hold
unless the three points $x,y,z$ are co-linear. For a more in depth discussion
of the non-decreasing property of the midpoint, see Lemma 12 in the appendix.
In our applications we will find that $r(n)=\ell^{2}$, where $\ell$ is the
size of the cliques used to construct the distance matrix. However if one has
distance estimates at a faster rate, then these can be plugged in to obtain
faster curvature estimates. Moreover, this is general to many problems where
distances are being estimated. Of course, this requires the knowledge of a
true midpoint as well. In the following sections, we discuss the use of
approximate midpoints, as well as introduce techniques for estimating
distances on the underlying manifold of model (equation (1)).

###  2.5 Midpoint Bias

Until this point we have worked under the assumption that there is a known
midpoint for a particular pair of points. However, in general, this is not
known. Suppose that $m^{\prime}$ is a point which is a surrogate midpoint
($m^{\prime}$) for the midpoint between $y,z$ and we have the analogous
distance $d_{xm^{\prime}}$. By a Taylor series expansion:

$\displaystyle 0$ $\displaystyle=g(\kappa^{\prime},d^{\prime})-g(\kappa,d)$
$\displaystyle=\nabla_{\kappa}g(\kappa,d)(\kappa-\kappa^{\prime})+\nabla_{d_{xm}}g(\kappa,d)(d_{xm}-d_{xm^{\prime}})$
$\displaystyle+o(|\kappa-\kappa^{\prime}|)+o(|d_{xm}-d_{xm^{\prime}}|)$
$\displaystyle\implies|\kappa-\kappa^{\prime}|$
$\displaystyle\approx(\nabla_{\kappa}g(\kappa,d))^{-1}\nabla_{d_{xm}}g(\kappa,d)|d_{xm}-d_{xm^{\prime}}|$
$\displaystyle\leq(\nabla_{\kappa}g(\kappa,d))^{-1}\nabla_{d_{xm}}g(\kappa,d)d_{mm^{\prime}}$

Hence, the bias will scale approximately linearly as a function of
$d_{mm^{\prime}}$ for small values of $d_{mm^{\prime}}$.

We next provide an outline involving how fast we can expect surrogate
midpoints form. We first introduce the definition of geodesic convexity. A
subset $A\subset\mathcal{M}^{p}$ is geodesically convex if the geodesic
between any two points in $A$, is contained within $A$ itself. Here convexity
on a manifold will refer to geodesic convexity. In the following theorem, let
$f_{m,3}(m_{1},m_{2})$ denote the joint density function of a pair of
midpoints with a shared endpoint, and let $f_{m,4}(m_{1},m_{2})$ denote the
joint density of two midpoints without an endpoint shared. These two densities
will be functions of the unknown manifold, $\mathcal{M}$ and the distribution
on the latent positions $Z$.

######  Theorem 4.

Suppose that $Z_{i}\stackrel{{\scriptstyle iid}}{{\sim}}F_{Z}$ are $K$ points
sampled iid from an unknown latent distribution on a simply connected manifold
$\mathcal{M}^{p}$. Denote this set of points
$\\{Z_{i}\\}_{i=1}^{K}=\mathcal{D}_{K}$ Suppose there exists a convex region
$A$ for which

  * •

(C1) $f(z)\geq\alpha>0,\text{ for all }z\in A$

  * •

(C2) $\dim(A)=p\geq 2$

  * •

(C3) $f_{m,3}(m_{1},m_{2}),f_{m,4}(m_{1},m_{2})\leq\alpha_{m}<\infty$

for all $z\in A$ where $f$ is the density function corresponding to $F_{Z}$.

Define the statistic

$\Phi(\mathcal{D}_{K}):=\min_{x,y,z\in\mathcal{D}:x\not=y,x\not=z,y\not=z}d\Bigl{(}m(y,z),x\Bigr{)}$

Then

$\Phi(\mathcal{D}_{K})=\mathcal{O}_{P}(K^{-3/p})$ (4)

Our approach to showing the above result draws on similarities to Cai et al.
(2013) and Brauchart et al. (2015) which discuss convergence of the minimum
distance between any two points sampled uniformly on a hypersphere. The
authors show that $\min_{i,j}d(Z_{i},Z_{j})=\mathcal{O}_{P}(K^{-2/p})$, though
they also derive an exact distribution for $\min_{i,j}d(Z_{i},Z_{j})$ using
the extra assumption of uniformity on the sphere. The authors use a technique
by recursively computing the probability that a point is at least a radius
$\epsilon$ away from from the previous $K$ points. Since at each placement of
a new point, there are $\binom{K}{2}$ midpoints as opposed to $K$ current
points, leading to the faster rate we observe.

Assumption (C1) ensures that there does exist a geodesically convex regions
which midpoints can form. Assumption (C2) ensures that this region is equal to
the dimension of the ambient space. Lastly, assumption (C3) is relatively mild
so long as we have a smooth manifold and a continuous latent density. The
proof is left to the Section A.6 in the supplementary materials but relies on
the following Theorem regarding medians of arbitrarily correlated random
variables which may be of independent interest.

######  Theorem 5.

Consider a set of continuous random variables $X_{i}\in\mathbb{R}$ which are
marginally identical but have an arbitrary dependency. Let
$\mathbb{\widehat{M}}[X]$ denote the sample median of the set of these random
variables. Then

$P(\mathbb{\widehat{M}}[X]\leq t)\leq 2P(X\leq t)$ (5)

We find the proof in the supplementary materials in Section A.5. This bound is
only useful up to $t=\mathbb{M}[X]$ at which point the upper bound is $1$. The
implication here is that we only need most of the midpoints between any two
points to be reasonably separated. We later illustrate how we will find good
midpoints and so we can verify their existence in any observed dataset.

###  2.6 Upper and Lower Bounds on the Curvature

We have shown above that under some mild conditions, we will find good
midpoints eventually. There are, however, some cases, such as some
distributions $F_{Z}$ defined on discrete domains on the latent manifold, for
which good midpoints are not guaranteed to form. Additionally, under any
finite sample, we will not have an exact midpoint alignment. To alleviate
these issues we use a a similar strategy to derive upper and lower bounds for
the curvature as long as the surrogate midpoint exists in the same region of
constant curvature.

######  Theorem 6 (Curvature Bounds).

Let $x,y,z,m^{\prime}\in\mathcal{M}^{p}(\kappa)$ where $m^{\prime}$ need not
be the exact midpoint of $y$ and $z$, but rather a surrogate midpoint. Then
let $d_{jk}$ denote the distance between points
$j,k\in\\{x,y,z,m^{\prime}\\}$. Let $\kappa_{u}$ and $\kappa_{l}$ denote the
solutions

$\displaystyle
g(\kappa_{u},\boldsymbol{d}(x,y,z),d_{xm^{\prime}}+d_{mm^{\prime}}(\kappa_{u},\boldsymbol{d}(m^{\prime},y,z)))=0$
$\displaystyle
g(\kappa_{l},\boldsymbol{d}(x,y,z),d_{xm^{\prime}}-d_{mm^{\prime}}(\kappa_{l},\boldsymbol{d}(m^{\prime},y,z)))=0$

then $\kappa_{l}\leq\kappa\leq\kappa_{u}$.

If $d_{mm^{\prime}}(\kappa,\boldsymbol{d}(m^{\prime},y,z)))=0$ then the upper
and lower bounds converge. Similar to the curvature estimate
$\widehat{\kappa}$, given a noisy estimate of the distances, we can estimate
the upper and lower bounds of the curvature. We will primarily use this for
constant curvature testing as we will see in Section 5.

##  3 Distance Matrix Estimation

In the previous section, we have not required any distances related to the
latent distance model, but rather arbitrary sets of distances. In this section
we introduce a method for estimating a distance matrix under the latent
distance model. If we knew the curvature of the manifold, we could use maximum
likelihood or other statistical techniques to estimate latent positions and,
therefore the latent distances. This approach is commonly used in practice
when the manifold geometry is assumed _a priori_. In our case, however the
goal is to estimate the underlying curvature, so we need a notion of distance
that is consistent across geometries.

We exploit a similar strategy to Lubold et al. (2023) in constructing a
distance matrix. Specifically, we construct distances using connections across
cliques, or fully connected subgraphs. This approach has two advantages.
First, since cliques represent multiple nodes we have multiple opportunities
for there to be connections between two cliques, meaning we can use the
fraction of realized ties between two cliques as an estimate for the
probability of connection between the two cliques. Since our approach relies
on a specific probability model for connections, we can easily convert
connection probabilities into distances using the parametric model. Second, in
general, cliques will correspond to points near each-other in space. This
means that we can consider cliques as “points” on the manifold that can be
used to identify a distance matrix.

We now proceed to describe our approach in detail. In order to estimate the
distance matrix, we first assume we have access to the person-specific effects
and focus only on the distances. We later relax this as it will be negligible
if we can estimate them at a fast enough rate. Under this model, we nodes
within a clique have nearly the same latent position. We will also for
illustration purposes here, first assume that this holds, and later illustrate
the rate at which the diameter of the set of latent positions within a clique
converges. Let $X,Y\in\\{1,2,\dots,n\\}$ denote sets of nodes which denote
non-overlapping cliques. Then the average probability of connection can be
used to identify the latent distance, if we can also estimate the average of
random effects ($\nu$). Let $|W|$ denote the size of $W\in\\{X,Y\\}$. We
define the average probability of connection across cliques $X,Y$ as the
following.

$\displaystyle p_{XY}$ $\displaystyle:=\frac{1}{|X||Y|}\sum_{x\in X}\sum_{y\in
Y}p_{xy}$ $\displaystyle=\frac{1}{|X||Y|}\sum_{x\in X}\sum_{y\in
Y}\exp(\nu_{x}+\nu_{y}-d_{xy})$
$\displaystyle=\exp(-d_{xy})\frac{1}{|X||Y|}\sum_{x\in X}\sum_{y\in
Y}\exp(\nu_{x}+\nu_{y})$ $\displaystyle\implies d_{xy}$
$\displaystyle=-\log(p_{XY})+\gamma_{X}+\gamma_{Y}$ $\displaystyle\gamma_{W}$
$\displaystyle:=\log\left(\frac{1}{|W|}\sum_{w\in
W}\exp(\nu_{w})\right)\quad\text{ for }W\in\\{X,Y\\}$

Since there are $\ell^{2}$ possible connections between a pair of cliques of
size $\ell$, we can construct an asymptotically normal estimator of $p_{XY}$
by just the plug in sample average of connections. However, in order for these
asymptotics to hold, we must estimate $\gamma_{X}$ and $\gamma_{Y}$ at rates
$o_{P}(\frac{1}{\ell})$, so that these are asymptotically negligible nuisance
parameters. A final issue that remains however is the restriction of all
distance estimates to ensure that they preserve a the properties of a metric.
We will illustrate our solutions to these problems in the following sections.

Unless latent positions are nearly (or exactly) in the same location, large
cliques are rare in latent space models. The next lemma illustrates a rate of
convergence of these latent locations relative to the size of a clique.

######  Lemma 7.

Assume the latent distance model as in equation (1)
$\delta=\max_{i,j}d(Z_{i},Z_{j})$ and let $C_{\ell}$ denote the event that
nodes $i=1,\dots,\ell$ form a clique.

Let $\mu_{d}=\mathbb{E}[d(Z_{i},Z_{j})]$ where $Z_{i}$ and $Z_{j}$ are drawn
independently from $F_{Z}$. Then for any $0<\widetilde{\mu}_{d}<\mu_{d}$, if

  * •

(D1) $F_{Z}$ admits a continuous density $f_{Z}$

$\displaystyle\delta|C_{\ell}=o_{P}(\exp(-\widetilde{\mu}_{d}\ell)).$

See the supplementary materials Section A.8 for a proof. The main implication
is that it is reasonable to treat the latent positions as a single point when
nodes are within a clique. The assumption in Lemma 7 ensures that we have some
region in the latent space where positions can be drawn near each other, if
there is a point mass anywhere on the latent space, then this condition is
immediately satisfied, otherwise, if $F_{Z}$ admits a smooth continuous
density ($f_{z}$) then it is also satisfied in any region of positive density
value, since $\mathcal{M}^{p}$ satisfies the locally Euclidean property (A2).

In order to estimate $d_{xy}$ we must estimate the nuisance parameters
$\gamma_{X}$ and $\gamma_{Y}$ at a fast enough rate. If a set of nodes are
generated from a common point, as in a clique, we can use the number of
connections a point has to the rest of the network in order to identify the
relative difference in their random effects. We assume that all points within
a clique share the same latent position. We show in Lemma 7 that this
assumption is reasonable as the latent radius of the points shrinks
exponentially fast. Consider a set of nodes $i\in X$ where $X$ denotes the
indices of a clique. Let $k$ denote the indices of any other node in the
network. Then we can identify the differences of random effect values within a
clique

$\displaystyle P(A_{ik}=1)$
$\displaystyle=\int\exp(\nu_{i}+\nu_{k}+d(z_{i},z_{k}))dF_{Z}(z_{k})dF_{\nu}(\nu_{k})$
$\displaystyle=\exp(\nu_{i})\int\exp(\nu_{k}+d(z_{i},z_{k}))dF_{Z}(z_{k})dF_{\nu}(\nu_{k})$
$\displaystyle\implies\frac{P(A_{ik}=1)}{P(A_{jk}=1)}$
$\displaystyle=\exp(\nu_{i}-\nu_{j}).$

Other than in pathological cases, $\ell\ll n$ where $n$ is number of nodes in
the network. As a result, we are able to measure the differences between the
random effects of members of a clique quite effectively, at a rate of
$O_{P}(\frac{1}{\sqrt{n}})$. Therefore, we can define
$\nu_{i}=\nu_{X}+\Delta\nu_{i}$ which allows us to identify $\Delta\nu_{i}$ if
$\nu_{X}=\min_{i\in X}\nu_{i}$ were known,

$\displaystyle\gamma_{X}$
$\displaystyle=\nu_{X}+\log\left(\frac{1}{|X|}\sum_{i\in
X}\exp\left(\Delta\nu_{i}\right)\right).$

The remaining question is to estimate $\nu_{X}$. We will argue that the
smallest random effect under mild conditions will converge to $0$ at an
exponentially fast rate.

######  Lemma 8.

Let $C_{\ell}$ denote the event that $i\in\\{1,2,\dots,\ell\\}$ form a clique.
Suppose that

  1. 1.

(E1) $F_{\nu}$ admits a continuous density $(f_{\nu})$ on $(-\infty,0]$ and
$f_{\nu}(0)>0$

Let $\mu_{\nu}:=E[\nu]$. Then for any $\mu_{\nu}<\widetilde{\mu}_{\nu}<0$

$\min_{i}\nu_{i}|C_{\ell}=o_{P}\left(\exp\left(-|\widetilde{\mu}_{\nu}|\ell\right)\right)$
(6)

This theorem states that the nodes which we find in a clique, tend to have
near $0$ random effects, and thus it is reasonable to set $\max_{i}\nu_{i}=0$,
when estimating the random effects. Note that this is an exponential rate
which can be arbitrarily close to
$\left(\exp\left(-|\mu_{\nu}|\ell\right)\right)$.

######  Lemma 9.

Let $W\subset\\{1,2,\dots,n\\}$ denote a subset of indices which form a clique
($C_{\ell}$). Let $d_{i}$ denote the degree of node $i\in W$ where $|W|=\ell$
and assume that the points in the clique have a common latent position. Denote
the estimator of $\gamma_{W}$. Let $\mu_{\nu}:=E[\nu]$. Then for any
$\mu_{\nu}<\widetilde{\mu}_{\nu}<0$

$\displaystyle\widehat{\gamma}_{W}$
$\displaystyle=\log\left(\frac{1}{\ell}\sum_{i\in W}\frac{d_{i}}{\max_{j\in
W}d_{j}}\right).$

Then if (E1) in Lemma 8 holds, then define $\widetilde{\mu}_{\nu}$ as in Lemma
8

$\widehat{\gamma}_{W}-\gamma_{W}=\mathcal{O}_{P}\left(\max\left\\{\exp(\widetilde{\mu}_{\nu}\ell),\frac{1}{\sqrt{n}}\right\\}\right)$
(7)

Therefore, estimation of $\gamma_{W}$ within a clique can occur at an
exponential rate, meaning that this estimation will be negligible compared to
the average cross-clique probabilities.

######  Theorem 10.

Let $\widehat{d}$ be a set of clique distances measured via connections
between cliques of size $\ell$ with fixed positions $x$ and $y$ respectively.
Assume that we have indices of cliques $X,Y$ for which $|X|=|Y|=\ell$ and
$X\cap Y=\emptyset$. We also assume the nuisance estimators exist for
$\gamma_{W}$ (condition F1) and that the Lindeberg CLT condition is satisfied
for the counts between cliques

  * •

(F1) $\widehat{\gamma}_{W}-\gamma_{W}=o_{P}(\frac{1}{\ell})$ for $W\in X,Y$

  * •

(F2) $\lim_{\ell\to\infty}\frac{\sum_{x\in X}\sum_{y\in
Y}(A_{xy}-p_{xy})^{2}I(|A_{xy}-p_{xy}|>\epsilon\ell^{2}\sigma_{\ell})}{\sigma_{\ell}}=0$
for all $\epsilon>0$

Where

$\displaystyle p_{XY}$ $\displaystyle=\frac{1}{\ell^{2}}\sum_{x\in
X}\sum_{y\in Y}p_{xy}$ (8) $\displaystyle\sigma_{\ell}$
$\displaystyle=\frac{1}{\ell^{2}}\sum_{x\in X}\sum_{y\in Y}p_{xy}(1-p_{xy})$
(9)

Then

$\sqrt{\ell^{2}\frac{\sigma_{\ell}}{p_{XY}}}\biggl{(}\widehat{d}_{xy}-d_{xy}\biggr{)}\to_{d}N(0,1)$
(10)

The conditions (F1) simply states that we have estimators of the nuisance
parameters, and (F2) simply states that we have the standard Lindeberg
condition for the central limit theorem when random variables are independent,
but not necessarily identical. The proof is found in Section A.4. The main
implication here is that if we can estimate the random effects fast enough
within a clique, then we are able to estimate the random effects at a rate
which is negligibly fast compared to the rate of $\widehat{p}_{XY}$.
Additionally, we note that since in general, the random effects are all
converging to $0$ within a clique, that for large cliques
$p_{XY}\approx\exp(-d_{XY})$, thus allowing for simplified expressions for
$\sigma_{\ell}$ and $p_{XY}$.

For a set of cliques $\mathcal{C}:|\mathcal{C}|=K$, we can compute a distance
matrix $\widehat{D}\in\mathbb{R}^{K\times K}$ such that we best fit the
observed data while maintaining the restriction of the triangle inequalities.
In this characterization and also in Lubold et al. (2023), the distance matrix
estimation does not in fact necessarily respect the restrictions that all
distance matrices must satisfy, namely non-negativity, and the triangle
inequality. For example, we construct an enumeration of the cliques of size
$7$ or greater from the General Relativity co-authorship network of Leskovec
et al. (2007) in Figure 3. Any pair of cliques which does not share an edge
will inherently be estimated to have infinite distance, which prevents
estimation of curvature. However applying the triangle inequality via other
estimated distances in the graph would allow us to construct an upper bound.
This is a challenge observed in Lubold et al. (2023), which required cliques
that are pairwise connected.

Figure 3: Clique subgraph of co-authorship network in ArXiv General
Relativity. Cliques shown by color.

To address this challenge, we posit a similar estimation problem, while
respecting the triangle inequality. The following estimation problem is posed
below. Let $X$ denote a set of indices corresponding to a single clique in a
graph $G$. Let $\mathcal{C}$ denote the set of cliques in an observed graph.
Given a set of fixed effects, we can propose a maximum likelihood optimization
problem for the distance matrix $D\in\mathbb{R}^{K\times K}$ under the
generative model. By the approximation that cliques have a common latent
position, we define the following likelihood function of the clique subset of
the data $\ell_{\mathcal{C}}$

$\displaystyle\ell_{\mathcal{C}}(D;\boldsymbol{\nu})=\sum_{X,Y\in\mathcal{C},i\in
X,j\in Y}A_{ij}\left(\nu_{i}+\nu_{j}-d_{xy}\right)$
$\displaystyle+\sum_{X,Y\in\mathcal{C},i\in X,j\in
Y}(1-A_{ij})\log\left(1-\exp\left(\nu_{i}+\nu_{j}-d_{xy}\right)\right).$

Then we can define the maximum likelihood optimization problem, after
estimating a set of random effects $\boldsymbol{\widehat{\nu}}$

$\displaystyle\widehat{D}$ $\displaystyle=\sup_{D\in\mathbb{R}^{K\times
K}}\ell_{\mathcal{C}}(D;\boldsymbol{\widehat{\nu}})$
$\displaystyle\text{s.t.}d_{xy}$ $\displaystyle\geq
0,\quad\text{Diag}(D)=0,\quad\text{tr}(E_{s}^{\top}D)\geq 0\forall
s\in\mathcal{S}$

where $\mathcal{C}$ is a list of clique indices and $E_{s}$ are matrices full
of zeros other than $3$ indices, $i,j,k$ for which $E[i,j]=E[i,k]=1=-E[j,k]$.
This set of indices allows us to respect the triangle inequality for all
possible triplets of distances. We define $\mathcal{S}$ to be an enumeration
of all such indices $E_{s}$. There are $3\binom{K}{3}$ such restrictions. The
set of feasible distance matrices that satisfy these constraints form a
polyhedron of interior dimension $\binom{K}{2}$. This optimization procedure
is a convex problem under the space of symmetric matrices with many
restrictions. In practice, we use CVXR to solve this system (Fu et al., 2020).
For a greater gain in computational speed, we use the MOSEK solver for the
constrained optimization (Andersen and Andersen, 2000).

Though in its current form, the problem is numerically challenging to solve
due to the number of restrictions. A natural option is to approximate the
likelihood via a second order Taylor expansion and solve this problem
successively. This is exactly Newton’s method. Since the Hessian is diagonal,
$\frac{\partial^{2}}{\partial d_{ij}\partial
d_{kl}}\ell_{\mathcal{C}}(D;\boldsymbol{\nu})=0$ if $d_{ij}\not=d_{kl}$, this
can be made into a more efficient quadratic program which can be solved faster
in CVXR.

Let $g(D;D_{0};\boldsymbol{\nu})$ be the second order Taylor series
approximation to $\ell_{\mathcal{C}}$ about a matrix $D_{0}$. We can
successively solve for $\widehat{D}_{t+1}$ the following constrained
optimization problem

$\displaystyle\widehat{D}_{t+1}$ $\displaystyle=\sup_{D\in\mathbb{R}^{K\times
K}}\widetilde{g}(D,\widehat{D}_{t};\boldsymbol{\widehat{\nu}})$ (11)
$\displaystyle\text{s.t.}d_{xy}$ $\displaystyle\geq
0,\quad\text{Diag}(D)=0,\quad\text{tr}(E_{s}^{\top}D)\geq 0\forall
s\in\mathcal{S}$

which can be iteratively computed until $\ell_{\mathcal{C}}$ increases less
than some threshold. Each iteration is a linear constrained quadratic program
which we implement in CVXR. For further details, see the supplementary
materials in Section B.2.

###  3.1 One-Step Estimation

In practice, running the optimization step for many iterations can be
computationally costly. This is particularly problematic in bootstrap testing
for constant curvature (Section 5) It is well known in maximum likelihood
estimation , one only needs to perform one Newton step for asymptotic
efficiency (Le Cam, 1956). As a result, one can start with any consistent
estimator of the distance matrix $D$ and apply a single Newton step from
equation (11) and obtain the same asymptotic distribution, and therefore in
practice, only one step is needed.

##  4 Practical Considerations and Statistical Tests

In this section we consider some of the practical issues with estimation. In
particular, how to choose a good midpoint, and how to filter points which we
know to have extremely high variance. In general, we would like to find a set
of points for which
$\widehat{d}_{ym}\approx\widehat{d}_{zm}\approx\frac{1}{2}\widehat{d}_{yz}$
with $\widehat{d}_{yz}$ not excessively small.

###  4.1 Choosing Best Midpoint

In practice given a distance matrix we would like search over the space to
find the best midpoint. For an unknown latent geometry, we can exploit the
Fréchet mean of two points

$\displaystyle m^{*}$
$\displaystyle=\operatorname*{arg\,min}_{m\in\mathcal{M}}d^{2}_{ym}+d^{2}_{zm}.$

The midpoint between any two points $y,z$ is also the Fréchet mean. Therefore
we can search over the space of candidate entries of a distance matrix $D$ to
find the best midpoint available. In practice, we also would like to ensure
that the points $y,z$ are not too close to each-other. If $y$ and $z$ are too
close to each-other, then this leads to a very high variance estimator. So a
reasonable option is to normalize this quantity by $d_{yz}$. Given a true
midpoint $m^{*}$, then $d_{ym^{*}}=d_{zm^{*}}=\frac{1}{2}d_{yz}$, then

$\frac{d^{2}_{ym^{*}}+d^{2}_{zm^{*}}}{d^{2}_{yz}}=1/2$ (12)

Therefore in order to compute the best surrogate midpoint set we solve the
following problem

$\widehat{y},\widehat{z},\widehat{m}=\operatorname*{arg\,min}_{y,z,m}\frac{d^{2}_{ym}+d^{2}_{zm}}{d^{2}_{yz}}$
(13)

We can also add a term $\frac{|d_{ym}-d_{zm}|}{d_{yz}}$ which aids in
balancing the lengths of the distance to each point.

###  4.2 Selecting Reference Points

Though all we need are $4$ points in order to identify curvature, we are
restricted by the fact that we have to find pairs with good surrogate
midpoints. However, we will often observe many other potential locations which
we can use for $x$. For example, it is very hard to elucidate curvature when
$x,y,z$ are nearly co-linear, and thus we would like triangles $x,y,z$ which
are closer to equilateral.

We exploit this by considering a scaled version of the triangle inequality.
Let $C_{\Delta}\in[1,2]$ be a constant which determines the flatness of the
allowed triangles $x,y,z$. Then we select only the $x$ such that

$d_{ij}+d_{jk}\geq C_{\Delta}d_{ik}\quad\forall(i,j,k)\in(x,y,z).$ (14)

We can estimate a single curvature my taking the median across the values of
$x$. This tuning parameter allows us to pick the triangles closest to
equilateral, which tend to give the best estimates of the curvature. Letting
$C_{\Delta}=1$ allows for all triangles, no matter how flat they are and
$C_{\Delta}=2$ will only permit exact equilateral triangles. In practice we
find that a value of $C_{\Delta}\in[1.05,1.7]$ is effective and we set the
default to be $1.3$. Setting $C_{\Delta}$ to be too large results in no
triangles found, and setting $C_{\Delta}$ too small will result in using
triangles which are very flat, and often suffer from finite sample bias as
well as high variance in their respective $\widehat{\kappa}$ estimates. In
practice, given a midpoint set $y,z,m$ and a distance matrix $\widehat{D}$ we
let $\boldsymbol{x}=\boldsymbol{x}(\widehat{D},C_{\Delta})$ which are the
points satisfying equation (14), then we can let
$\widehat{\kappa}_{x}=\kappa(\widehat{D},\boldsymbol{d}(x,y,z,m))$, and let
$\widehat{\kappa}_{med}=\text{median}\\{\widehat{\kappa}\\}_{x\in\boldsymbol{x}}$.
Given $\boldsymbol{x}:=\boldsymbol{x}(C_{\Delta},\widehat{D},y,z,m)$ we can
construct the median estimator $\widehat{\kappa}_{med,\boldsymbol{x}}$ of the
curvature in equation (15).

$\widehat{\kappa}_{med,\boldsymbol{x}}=\text{median}_{x\in\boldsymbol{x}}\kappa(\widehat{D},y,z,m,x)$
(15)

We can illustrate this phenomenon using the theoretical variance computed for
the asymptotic results. In Figure 4 the smallest variance reference points in
each of these curvatures are the ones that form nearly equilateral triangles
with the reference points $y,z$. Additionally, the variance of the estimator
tends to be larger in a given reference location as the curvature $\kappa$
decreases.

Hence for a set $(y^{(j)},z^{(j)},m^{(j)})$ we can compute a set of estimates
$\\{\widehat{\kappa}^{(y^{(j)},z^{(j)},m^{(j)})}\\}^{J}_{j=1}$ using all the
$x$ positions such that equation (14) is satisfied. If we wish to compute a
single curvature estimate $\widehat{\kappa}$ then a simple method will be to
return the median over the set of estimates.

Figure 4: Variance as a function of $x$ position. True points $(y,m,z)$ in
black.

###  4.3 Simulations

We construct a number of simulations in order to illustrate the performance of
our curvature estimator, under the full model complexity. This will involve
first sampling the parameters of a latent position cluster model to draw
locations and variances independently, then sampling positions and random
effects from the random latent position cluster model. This is to illustrate
how even in the absence of a true placed, midpoint, that they will still occur
from random draws. We simulate from the following model:

$\displaystyle\mu_{k}$ $\displaystyle\sim
F_{\mu}^{\mathcal{M}^{p}(\kappa)}(\boldsymbol{O},R)\quad$ $\displaystyle
k\in\\{1,\dots,\lceil\sqrt{\rho}50\rceil\\}$ $\displaystyle\sigma^{2}_{k}$
$\displaystyle\sim F_{\sigma^{2}}\quad$ $\displaystyle\boldsymbol{\pi}\sim
Dirichlet(\boldsymbol{2})$ $\displaystyle\rho$
$\displaystyle\in\\{\frac{1}{\sqrt{2}},1,2\\}\quad$
$\displaystyle\kappa\in\\{-2,-1,-0.5,0,0.5,1\\}$

where $F_{\mu}^{\mathcal{M}^{p}(\kappa)}(\boldsymbol{O},R)$ denotes the prior
distribution on the latent positions which is a uniform ball with radius $R$
for $\kappa>0$, and two concentric balls, one with radius $R$ and the other
with $R/2$ with equal probability for $\kappa\leq 0$. This setup facilitates
forming midpoints since the volume of a ball grows exponentially as $\kappa$
decreases. In all simulations, we set $R=2.5$. This process determines a
random latent position cluster model (LPCM). The locations of $\mu_{k}$ as
well as the relative sizes of $\sigma^{2}_{k}$ determine where cliques are
most likely to form in the latent space. The parameter $\rho$ refers to the
scale of the network, allowing for the number of centers to grow with $\rho$,
$\boldsymbol{\pi}$ refers to the vector of cluster probabilities in the
mixture model, $\boldsymbol{\mu}$ and $\boldsymbol{\sigma^{2}}$ refer to the
cluster mean scale parameters and $\kappa$ the true curvature of the space.
The parameters $\boldsymbol{\mu},\boldsymbol{\sigma}^{2},\boldsymbol{\pi}$
determine a latent position cluster model, where the positions, $Z_{i}\sim
F_{Z}:=F_{Z}(\boldsymbol{\mu},\boldsymbol{\sigma}^{2},\boldsymbol{\pi};\kappa)$.
The mixture components correspond to the heat kernels in spherical space (Von-
Mises Fischer distribution), Gaussian distribution in the Euclidean space, and
the wrapped normal distribution in hyperbolic space (Nagano et al., 2019). In
all cases, let $p=3$, so that we are also not forcing the positions to be
drawn on a submanifold of dimension $2$. The randomness the latent position
cluster model incorporates the fact that we are not assigning good midpoints
exactly, but finding them in the data each time we simulate a matrix. We let
$F_{\sigma^{2}}$ denote the gamma distribution with shape parameter $1/4^{2}$.

We then sample a random adjacency matrix $A$ via the sampled LPCM

$\displaystyle n$ $\displaystyle=5000\rho,\quad Z_{i}\sim
F_{Z},\quad\nu_{i}\sim F_{\nu,\rho}$ $\displaystyle p_{ij}$
$\displaystyle=\exp(\nu_{i}+\nu_{j}-d(Z_{i},Z_{j}))$ $\displaystyle A_{ij}$
$\displaystyle\sim\text{Bern}(p_{ij})$

where $F_{Z}=F_{Z}(\boldsymbol{\pi},\boldsymbol{\sigma}^{2},\boldsymbol{\mu})$
is a particular draw of the random latent position cluster model. And
$F_{\nu,\rho}$ is a trimmed normal distribution with mean $-3\rho$ and
standard deviation $3\rho$ trimmed at $0$ and $\log(2/\sqrt{n})$ so that
$\nu_{i}$ remain non-negative and to prevent isolated nodes. We let the
minimum clique size used in the estimator be of size $\ell=(8+4\log(\rho))$,
which tends to generate 30-40 under this model. We set $\ell$ to be the
minimum, but allow for cliques which are larger. We repeat this $500$ times
for each scale and curvature setting. For ease of computation, we find the
cliques using the maximal_cliques function in the igraph R package (Csardi et
al., 2006) clusters of the model. See Section F for additional graph
statistics summaries over the simulations.

Figure 5 presents the results. We see that, as the clique size increases, we
see a reduction in bias and variance. Also, the true curvature impacts the
variance of the estimate at any particular sample size. This is for two
reasons. Firstly, as we saw in Figure 4, the variance of the estimate is
simply larger when $\kappa$ is negative in nearly all regions of the space.
Secondly, due to the vastness of the negatively curved spaces, we tend to have
poorer quality midpoints form as well as fewer reference points $(x)$ which
form good triangles which we can take medians across. In each of the
simulations, we limit the number of cliques used in the estimator at 35 (50
for $\kappa\leq-1)$) for computational convenience for numerous simulations;
though in real data applications, this number can be larger in the practice.
From the distance matrix, we compute the best midpoint set $y,z,m$ according
to equation (12) and compute the median of $\kappa$ estimates for the
$\boldsymbol{x}$ which satisfy equation (14) for $C_{\Delta}=1.5$. We plot the
results in Figure 5.

Figure 5: Consistency of Curvature Estimator. Upper error bars indicate $0.95$
quantile of simulations and lower indicate $0.05$ quantile. Central dots
indicate the mean after trimming outliers larger than $\pm 100$ (0.107% of the
observations)

##  5 Application: Constant Curvature Test

A natural test one may wish to conduct is to verify whether a particular
matrix can be embedded in a space of constant curvature

$H_{0}:\mathcal{M}=\mathcal{M}(\kappa).$

This test could provide a model diagnostic (e.g. testing whether a latent
variable model that assumes a single constant manifold is appropriate) or
provide meaningful insights into heterogeneity in the structure of the graph.

In order to construct an appropriate test, we first consider a method of
subsampling from the distribution of the cliques in order to approximate the
distribution of the distance matrix estimate. This is based on the subsampling
approach of Politis and Romano (1994), which can be used to approximate the
distribution of a random variable via subsampling under conditions weaker than
the bootstrap, similar to the strategy in Lubold et al. (2023). This is
highlighted in Algorithm 1. For simplicity, we illustrate the algorithm with
the subsampling rate $b=\frac{\ell-1}{\ell}$ where $\ell$ is the smallest
clique size.

Let $\mathcal{I}=\\{\textbf{I}_{k}\\}_{k=1}^{K}$ denote a set of indices
corresponding to cliques for which
$\textbf{I}_{k}\cap\textbf{I}_{k^{\prime}}=\emptyset$ when $k\not=k^{\prime}$.
Let $\widehat{\nu}_{\textbf{I}_{k}}$ denote the set of random effects indexed
by $\textbf{I}_{k}$.

Algorithm 1 Subsample Distribution of $\widehat{D}$

1:$G$, $B\geq 1$ and $\\{\textbf{I}_{k}\\}_{k=1}^{K}$

2:for $b\in\\{1,2,\dots,B\\}$ do

3: for $k\in\\{1,2,\dots,K\\}$ do Sample $|\boldsymbol{I}_{k}|-1$ nodes
without replacement from $\boldsymbol{I}_{k}$ and denote this set
$\boldsymbol{I}^{(b)}_{k}$

4: end for

5: Denote the set $\mathcal{I}^{(b)}=\\{\textbf{I}^{(b)}_{k}\\}_{k=1}^{K}$

6: Let $\mathbb{\nu}^{(b)}=\mathbb{\nu}_{\mathcal{I}^{(b)}}$ denote the
corresponding random effects estimates.

7: Let $\widehat{D}^{(b)}_{0}=f_{0}(A,\mathcal{I}^{(b)},\mathbb{\nu}^{(b)})$
Initial Consistent estimate

8: Let
$\widehat{D}^{(b)}_{1}=f_{1}(\widehat{D}^{(b)}_{0};A,\mathcal{I}^{(b)},\mathbb{\nu}^{(b)})$
One-step Estimate

9:end for

10:return
$\\{\widehat{D}^{(b)}\\}_{b=1}^{B}=\\{\widehat{D}^{(b)}_{1}\\}_{b=1}^{B}$

In Algorithm 1 $f_{0}$ denotes the distance estimation by the method of Lubold
et al. (2023) and sparse modification of the distances via the Floyd Warshall
Algorithm so that $\widehat{D}_{0}$ is a metric. This is known as a sparse
metric repair problem (Gilbert and Jain, 2017). For full details, see the
appendix Section B.2. The subsequent step $f_{1}$ simply is a short-hand for
applying the one-step estimation procedure of equation (11).

For a distance matrix $D$, let $\kappa_{u}(D;y,z,m,x)$ and
$\kappa_{l}(D;y,z,m,x)$ denote the upper and lower bound indexed by a set
$y,z,m,x$. Where $y,z,m$ are the endpoints and midpoints, and $x$ is a
reference point. Given a set of samples of a distance matrices
$\\{\widehat{D}^{(b)}\\}_{b=1}^{B}$ and a set of midpoint sets
$y^{(j)},z^{(j)},m^{\prime(j)}$ with corresponding reference points
$\boldsymbol{x}^{(j)}$, we can test whether the curvature is constant between
these regions $j\in\\{1,2,\dots,J\\}$:

$H_{0}:\kappa_{j}=\kappa_{j^{\prime}}\text{ for all
}j,j^{\prime}\in\\{1,2,\dots,J\\}$

using Algorithm 2.

Algorithm 2 Constant Curvature Test

1:$\\{\widehat{D}^{(b)}\\}_{b=1}^{B}$ and
$\\{y^{(j)},z^{(j)},m^{(j)},\boldsymbol{x}^{(j)}\\}_{j=1}^{J}$

2:for $b\in\\{1,2,\dots,B\\}$ do

3: for $j\in\\{1,2,\dots,J\\}$ do

4: for $x\in\\{\boldsymbol{x}^{(j)}\\}$ do

5:
$\widehat{\kappa}_{u,x}^{(b,j)}=\kappa_{u}(\widehat{D}^{(b)};y^{(j)},z^{(j)},m^{(j)},x)$

6:
$\widehat{\kappa}_{l,x}^{(b,j)}=\kappa_{l}(\widehat{D}^{(b)};y^{(j)},z^{(j)},m^{(j)},x)$

7: end for

8: Compute
$\widehat{\kappa}_{u}^{(b,j)}=\text{median}_{x\in\boldsymbol{x}^{(j)}}\widehat{\kappa}_{u,x}^{(b,j)}$

9: Compute
$\widehat{\kappa}_{l}^{(b,j)}=\text{median}_{x\in\boldsymbol{x}^{(j)}}\widehat{\kappa}_{u,x}^{(b,j)}$

10: end for

11: Compute
$\widehat{\kappa}_{u}^{(b)}=\text{min}_{j\in\\{1,2,\dots,J\\}}\widehat{\kappa}_{u}^{(b,j)}$

12: Compute
$\widehat{\kappa}_{l}^{(b)}=\text{max}_{j\in\\{1,2,\dots,J\\}}\widehat{\kappa}_{l}^{(b,j)}$

13:end for

14:Let $\widehat{\kappa}_{u,(m)}$ be the $m^{th}$ order statistic of
$\\{\widehat{\kappa}_{u}^{(b)}\\}_{b=1}^{B}$

15:Let $\widehat{\kappa}_{l,(m)}$ be the $m^{th}$ order statistic of
$\\{\widehat{\kappa}_{l}^{(b)}\\}_{b=1}^{B}$

16:Let $m=\min\\{m:\widehat{\kappa}_{u,(m)}\geq\widehat{\kappa}_{l,(B-m)}\\}$

17:return $p$-value: $\min\\{2m/B,1\\}$

The constant curvature test takes a subsample of the distance matrices
$\\{\widehat{D}\\}_{b=1}^{B}$ generated from Algorithm 1, and a set of
midpoints and reference points
$\\{y^{(j)},z^{(j)},m^{(j)},\boldsymbol{x}^{(j)}\\}_{j=1}^{J}$ to estimate the
upper and lower bounds on the curvature. In practice, each set
$\boldsymbol{x}^{(j)}$ can be selected by the triangle filtering method of
equation (14). To reduce the variance of the estimates, we take the median
across regions at each $x$ and we take the minimum of the upper bounds, (and
maximum of the lower bounds) to bound the corresponding distance matrix. We
repeat this over the set of $b\in\\{1,2,\dots,B\\}$ estimates of the
curvature, bounds and assess the subsampled quantiles of the total upper and
lower bounds. We lastly find the corresponding minmum quantile for which there
is a separation between these sampled bounds and return that as the $p$-value
for the subsampled constant curvature test.

###  5.1 Simulations

Under the same simulation setup as in Section 4.3 we can illustrate the
coverage of the constant curvature test as a function of clique size. Due to
computational complexity, we reduce this to a maximum clique size of $12$. We
see in Figure 6 that this test tends to be overly conservative in small sample
sizes, but returns to nearly nominal coverage when cliques are larger. This is
due to the fact that a poor quality midpoint leads to more conservative bounds
of $\kappa_{l},\kappa_{u}$, however, this can be improved with better aligning
midpoints.

Figure 6: False positive rate for constant curvature test.

###  5.2 Multiplex Network Constant Curvature Tests.

In this next simulation, we consider a model of multiplex (or multiview)
networks. Several methods exist for modelling multiplex networks via
extensions of the latent distance model, however, most assume the same
geometry latent space among views Salter-Townshend and McCormick (2017);
MacDonald et al. (2022). For example, Salter-Townshend and McCormick (2017)
model the multiple relationships between individuals in the Banerjee et al.
(2013) diffusion of microfinance dataset using Euclidean spaces. We illustrate
a simulated example where this is not the case, and how our method can be used
to detect this.

We first construct a latent position model for which multiple views are drawn
from a common set of latent positions, however, these positions are common
coordinates of spheres of curvature $(\kappa_{1},\kappa_{2})=(0.5,1.5)$
respectively. Additional details for the simulation are identical to the
consistency simulations in Section 4.3.

In this case for a set of latent positions, we simulate a multiplex network
for which latent spherical positions are the same for this network set,
however, they are embedded in two spheres of different radii and thus
different curvatures. We simulate $200$ pairs of network and test the
curvature difference and plot the power in Figure 7. In each view we compute
the optimal midpoint set $y^{(r)},z^{(r)},m^{(r)}$ from each view’s distance
matrix $\widehat{D}^{(1)}$ or $\widehat{D}^{(2)}$ via equation 13 and
subsample each distance matrix separately, and concatenate them to conduct the
constant curvature test. In Figure 7 we see that the power of the least grows
with the clique size.

Figure 7: Power of constant curvature of multiview.

###  5.3 Noncanonical manifolds

We next illustrate the ability of our method to detect non-constant curvature
within the latent space. Since many common latent space models are constantly
curved by assumption, it is of interest to verify whether this agrees with the
observed data.

Though the sphere, the hyperbolic and Euclidean spaces all have closed form
solutions for their distances, this is not true of arbitrary manifolds, and we
highlight our example on one which allows us to compute these distances.

We construct a latent space consisting of two adjacent spheres. For any two
points in the same sphere the distance is straightforward to compute. For any
two points $(x,y)$ in opposite spheres, the distance can be computed via the
distance to the origin $(1,0,0)$ in each of the spheres. This is due to the
fact any geodesic must pass through the connecting point, i.e. the origin.

$\displaystyle d_{\mathcal{M}}(x,y)$ $\displaystyle=d_{1}(x,o)+d_{2}(o,y)$

For these simulations, our manifold consists of two adjacent spheres of
curvature $\kappa_{1}=1,\kappa_{2}=1.5$ respectively. This manifold was chosen
as distances were straightforward to compute but also does not contain a
constant curvature. The latent cluster locations are sampled via uniforms
centered on opposite poles. We again simulate $200$ draws from the latent
position cluster model and test for constant curvature by finding the best
three non-overlapping sets minimizing equation (13). We plot the corresponding
power in Figure 8 which again increases as the clique size increases.

Figure 8: Constant Curvature Test Power: Adjacent Spheres.

###  5.4 Geometry of Coauthorship

We consider the sets of co-authorship networks introduced in Leskovec et al.
(2007). These consist of citation networks from High Energy Particle Physics,
General Relativity, Astrophysics, and Condensed Matter Physics, with sizes of
each of the networks seen in Table 1. These networks consist of authors as
nodes where an edge exists whether any pair of authors has co-authored a paper
on ArXiv in any of the specified subject areas between 1993 and 2003. When
these data were collected, these were the top 5 most common subject areas in
Physics.

In previous machine learning applications, hyperbolic network embeddings have
been successful in tasks such as link-prediction or node in citation networks
(Nickel and Kiela, 2017; Chami et al., 2019; Chamberlain et al., 2017). This
is due to the fact that the hierarchical structure (tree-like) structure
generally can be more easily embedded in a negatively curved space. We wish to
answer the question of, under a latent space assumption, is the data generated
from a space of constant curvature?

Physics Sub-field | Number of Nodes (n) | Number of Edges ($|E|$)
---|---|---
Astrophysics | 18771 | 396160
Condensed Matter Physics | 23133 | 186936
General Relativity | 5241 | 28980
High Energy Particle Physics | 12006 | 237010
High Energy Particle Physics (Theory) | 9875 | 51971

Table 1: Physics Co-authorship network sizes.

For each of these networks, we construct an estimate of the distance matrix
with random effects, followed by an estimate of the curvature. For each of
these, we use a triangle filtering constant $C_{\Delta}=1.2$ and a minimum
clique size seen in Table 2. We lastly apply our test to see if the difference
in curvature is present across the networks.

Physics Sub-field |  | Min Clique
---
Size ($\ell$)
| Number of Cliques
---
of Size $\geq\ell$ (K)
| Largest
---
Clique Size
Astrophysics | 19 | 57 | 57
Condensed Matter Physics | 12 | 52 | 26
General Relativity | 7 | 44 | 44
High Energy Particle Physics | 14 | 42 | 239
High Energy Particle Physics (Theory) | 7 | 42 | 32

Table 2: Clique size and number of cliques used to estimate distance matrix.

We estimate the curvature, at the best midpoint set for each of the networks,
along with the following p-values for tests of whether a network has constant
curvature.

Physics Sub-field | Curvature Estimate |  | Constant Curve
---
Test p-value
Astrophysics | -0.01378 | 0.030
Condensed Matter Physics | 0.107 | 1.000
General Relativity | 0.0989 | 1.000
High Energy Particle Physics | -0.986 | 1.000
High Energy Particle Physics (Theory) | 0.1674 | 0.240

Table 3: Curvature estimates from best midpoints and p-values for constant
curvature test $J=3$.

(a) High Energy Particle Physics (Theory) clique curvature labels.

(b) High Energy Particle Physics clique curvature labels.

(c) Astrophysics clique curvature labels.

Figure 9: Cliques Colored to nearest labeled curvature value. Blue, negative
and red positive, with $p$ values of constant curvature decreasing from left
to right.

We see that Astrophysics and General Relativity are estimated to have a large
negative curvature, which however, may not be constant in the case of the
Astrophysics citation network. In the Astrophysics network, at the best 3
midpoint sets, we estimate curvature to be $(-0.01378,-\infty,0.306)$. The
estimate of $-\infty$ comes from the fact there is a minimum distance of
$d_{xm}$ given the other 3, $d_{xy},d_{yz},d_{xz}$, and since the midpoints
are not exact, sometimes this length can be too short to embed in any of the
hyperbolic spaces. If the estimated distance is below this value we call the
estimate $-\infty$, and similarly $\infty$ if it is too large. This highlights
the fact that this network appears to have negatively curved, positively
curved and flat regions, and therefore models which reflect only a single
curvature, may not capture the individual level behavior of the network very
well. In contrast, condensed matter physics and HEP physics seem to have a
slight curvature, though both are nearly flat networks.

##  6 Application: Multiple change point detection

Another natural question one may seek to answer is whether a change in latent
curvature has occurred and being able to distinguish this from changes in
latent positions. Standard Gaussian multiple change point algorithms (for
example that of Harchaoui and Lévy-Leduc (2010)) are not immediately well-
suited to this problem due to the possibility of large outliers.

We apply the method of Fearnhead and Rigaill (2016) which proposes a multiple
changepoint detection algorithm under the presence of outliers. For a
particular network, we may measure a collection of estimates of the curvature
$\\{\widehat{\kappa}_{t}\\}_{t=1}^{T}$. Under this model, we assume that a
network has a constant curvature within a single time point $t$. We construct
an objective function for the changepoint problem

$\displaystyle L(\theta)$
$\displaystyle=\sum_{t=1}^{T}\widetilde{\ell}(\widehat{\kappa}_{t},\theta_{t})$

In our application, we let $\widetilde{\ell}$ be the biweight loss, however,
more general loss functions such as absolute deviation or Huber also available
(see Fearnhead and Rigaill (2016) for a more in depth discussion of robust
change point detection). In multiple change point detection algorithms a
penalty of $\beta$ for the number of segments included is also applied. Let
$J(\theta)$ denote a function of the number of breaks in the sequence
$\theta$. Then the full loss function is

$L(\theta;\beta)=L(\theta)+\beta J(\theta).$ (16)

In our setting, since we tend to have large negative outliers, we also apply
the $c\tanh(\cdot/c)$ function which smoothly truncates the extreme values. In
all simulations and applications where this is applied, we set $c=10$.

We next apply this to a simulation setting where we construct a sequence of
networks with latent positions.

$\displaystyle Z^{(t)}_{i}$
$\displaystyle=\mathcal{F}(Z^{(t-1)}_{i},\epsilon^{(t)}_{i})$

where $Z_{(t-1)i}$ is the location’s previous position, $\epsilon^{(t)}_{i}$
is a noise random variable sampled from the true cluster’s density, and
$\mathcal{F}$ stands for the mean on the sphere (the Fréchet mean). We set
three different curvatures

$\displaystyle\kappa_{1}=1.0\quad\kappa_{2}=0.15\quad\kappa_{3}=1.3$

with time changes at $t_{1}=18,t_{2}=35$ and total time $T=50$.

As illustrated before, we use the robseg package introduced in Fearnhead and
Rigaill (2016) and use the absolute value loss. We see that as $\ell$
increases, we are able to consistently estimate the true curvature function.
We plot the mean of the absolute deveation of the curvature estimate over
$200$ simulations as a function of the clique size using the default tuning
parameter from the Robseg R package. We plot the results in Figure 10 showing
consistency of the true curvature with respect to the mean absolute deviation.

Figure 10: Mean absolute deviation of curvature estimates over time window=.

###  6.1 LANL Netflow and Changepoint detection

For our second application, we apply our estimates of network curvature to Los
Alamos Unified Network and Host Dataset, and illustrate that changes in
curvature can be used to identify a red team attack. A red team attack is a
planned exercise by cybersecurity team that infiltrates a device network as a
test of network security.

This dataset consists of $89$ days of directed communication events between
$n=27436$ devices in the Los Alamos National Laboratory. The dataset consists
of $56$ normal days followed by a red team attack beginning on day $57$ and
continuing to day $89$.

Anomaly detection is of great interest for many cybersecurity applications and
recently latent space models have been shown to be a promising avenue of
detecting changes in node properties (Lee et al., 2019). We take a different
approach and illustrate how changepoint methods can be applied to sequential
estimates of curvature in this dataset, highlighting the role of curvature as
a summary of network behavior.

We consider a dataset preprocessed in the same manner as in Lee et al. (2019).
Edges are defined in this dataset as messages passed between nodes during a
particular time period. In order to maintain enough connections to find
cliques, we consider a connection to be formed if any message was passed in
the previous $4$ days. We then compute the curvature values at each time point
and take the median over the time-point. We scale each estimated value by
$c\tanh(\cdot/c)$ for $c=10$ in order to limit the influence of extreme
negative outliers. We then apply the off the shelf change point algorithm of
Fearnhead and Rigaill (2016). Using a $C_{\Delta}=1.4$ we estimate the
curvature at each time frame. We consider a minimum clique size of $\ell=5$.
Due to the small number of available cliques, and relative sparsity of the
dataset, we restrict the random effects to be $0$ and compute the
corresponding distance matrix.

(a) High regularization LANL changes in curvature estimates. Red team attack
illustrated in purple.

(b) Low regularization LANL changes in curvature estimates.

Figure 11: Two regularization values for LANL Netflow changepoint. Red team
attack illustrated in purple.

(a) Edge density changepoints.

(b) Tri. density changepoints.

Figure 12: Changepoints of daily LANL measurements using simple graph motifs.
Red team attack illustrated in purple.

In Figure 11 we show that the most substantial changepoint in curvature occurs
at the time of the red team attack. In contrast, these changes are much less
substantial in Figure 12 when using simple graph motifs from the daily
averages.

Since the time of detection after the event is most important, we wish to
investigate the time after detection as a function of the number of events
involved (alarm rate). We show that in Figure 13 that our method achieves a
much smaller detection delay given any alarm rate.

Figure 13: Time to first changepoint in days (TimeDelay) as a function of
Alarm Rate.

This suggests that models accounting for changes in network curvature may be a
promising avenue for the development of specialized models to detect anomalous
events in the online setting.

##  7 Discussion

Riemannian sectional curvature is a fundamental property of a manifold, and we
present a novel method to estimate it from a noisy distance matrix. Though our
motivating example involves estimating the distances of a latent distance
model from a random network, the curvature estimates (and constant curvature
tests) in this paper are more general and can be applied whenever one can
either estimate a distance matrix (or bootstrap or sub-sample their distance
matrices).

Though in our method one may worry about selection effects, given that we
choose a midpoint from the data. However, since we do not select for large or
small curvature, therefore selection effects are of minor concern. Another
concern may be the quality of midpoints found in real data. In practice, we
find good midpoints can be found in our real datasets. For example, the
midpoint minimization equation (12) which has a value at the true midpoint of
$0.5$ were found to be $0.5000003$, $0.5039$, $0.50049$, $0.50012$, $0.50043$
in the best midpoints of the 5 physics co-authorship networks respectively,
suggesting we have near midpoints occurring in real data.

Though we develop a test for detecting whether the curvature is constant on
the latent manifold, a natural open question is what should a practitioner do
if they find that their data are not represented well by this model. One might
instead use non-geometric models such as stochastic block models and their
variants, however, this also suggests the development of latent distance
models which are not restricted to constant curvature.

The development of such a model class as something which will scale to large
networks is of further interest. Computing gradients on arbitrary Riemannian
manifolds may prove to be challenging when the positions are in fact latent.
One promising approach we plan to investigate in future work is that of
product spaces for latent space models. This geometry has seen considerable
success in representation learning tasks such as Gu et al. (2018) or Zhang et
al. (2021).

One may consider the interpretation of the curvature parameter $\kappa$ when a
set of points does not lie within a canonical manifold. Our method merely fits
for each $(x,y,z,m)$ the particular constant curvature manifold, which will
embed these 4 distances, and taking the median will estimate the median of the
curvatures of these interpolating spaces. This idea is reminiscent of the
manifold learning methods of Li et al. (2017); Li and Dunson (2019) which use
spherelets to approximate manifolds rather locally linear approximations of
tangent spaces.

When considering the unified network and flow dataset from LANL, we must
transform the original dataset of pairs of individuals and times to one which
considers a set of graphs over time. A more realistic model for the data-
generating mechanism may be one related to a Poisson process (or other point
processes) with intensity parameters $\mu_{ij,t}\propto-d_{ij,t}$.

Future extensions may include methods for fitting models in the product space
geometry (as is done in Gu et al. (2018)) as well as movement towards
statistically consistent node-level definitions of network curvature as well
as more application-focused development of anomaly detection incorporating
curvature into the models.

##  8 Competing interests

No competing interest is declared.

##  9 Acknowledgments

The authors would like to thank Marina Meila, Shane Lubold, Adam Visokay and
Abel Rodriguez for their valuable suggestions. Research reported in this
publication was supported by the National Institute Of Mental Health of the
National Institutes of Health under Award Number DP2MH122405 and the National
Science Foundation under grant NSF SES-2215369. Partial support for this
research came from a Eunice Kennedy Shriver National Institute of Child Health
and Human Development research infrastructure grant, P2C HD042828, to the
Center for Studies in Demography and Ecology at the University of Washington.
Funding was also provided by the National Science Foundation under grant NSF
SES-2215369.

## References

  * Aliverti and Durante [2019] E. Aliverti and D. Durante.  Spatial modeling of brain connectivity data via latent distance models with nodes clustering.  _Statistical Analysis and Data Mining: The ASA Data Science Journal_ , 12(3):185–196, 6 2019.  ISSN 1932-1872.  doi: 10.1002/SAM.11412.  URL https://onlinelibrary.wiley.com/doi/full/10.1002/sam.11412.
  * Andersen [1970] E. B. Andersen.  Asymptotic properties of conditional maximum-likelihood estimators.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 32(2):283–301, 1970.
  * Andersen and Andersen [2000] E. D. Andersen and K. D. Andersen.  The mosek interior point optimizer for linear programming: an implementation of the homogeneous algorithm.  In _High performance optimization_ , pages 197–232. Springer, 2000\.
  * Banerjee et al. [2013] A. Banerjee, A. G. Chandrasekhar, E. Duflo, and M. O. Jackson.  The Diffusion of Microfinance.  _Science_ , 341(6144), 2013.  ISSN 10959203.  doi: 10.1126/SCIENCE.1236498.
  * Barkanass et al. [2022] V. Barkanass, U. Jost, and E. Hancock.  Geometric sampling of networks.  _Journal of Complex Networks_ , 10(4), 6 2022.  ISSN 2051-1310.  doi: 10.1093/COMNET/CNAC014.  URL https://academic.oup.com/comnet/article/10/4/cnac014/6644814.
  * Bassett et al. [2018] D. S. Bassett, P. Zurn, and J. I. Gold.  On the nature and use of models in network neuroscience.  _Nature Reviews Neuroscience 2018 19:9_ , 19(9):566–578, 7 2018.  ISSN 1471-0048.  doi: 10.1038/s41583-018-0038-8.  URL https://www.nature.com/articles/s41583-018-0038-8.
  * Begelfor and Werman [2005] E. Begelfor and M. Werman.  The World is not always Flat or Learning Curved Manifolds.  _School of Engineering and Computer Science, Hebrew University of Jerusalem., Tech. Rep_ , 3(8), 2005.
  * Berger [1962] M. Berger.  An extension of rauch’s metric comparison theorem and some applications.  _Illinois Journal of Mathematics_ , 6(4):700–712, 1962.
  * Blumenthal and Gillam [1943] L. M. Blumenthal and B. E. Gillam.  Distribution of Points in n-Space.  _The American Mathematical Monthly_ , 50(3):181, 3 1943.  ISSN 00029890.  doi: 10.2307/2302400.
  * Borgatti et al. [2009] S. P. Borgatti, A. Mehra, D. J. Brass, and G. Labianca.  Network analysis in the social sciences.  _Science_ , 323(5916):892–895, 2 2009.  ISSN 00368075.  doi: 10.1126/SCIENCE.1165821/ASSET/2421E2F8-2DC7-4CAD-8A84-471B43A3C443/ASSETS/GRAPHIC/323˙892˙F5.JPEG.  URL https://www.science.org/doi/10.1126/science.1165821.
  * Brauchart et al. [2015] J. S. Brauchart, A. B. Reznikov, E. B. Saff, I. H. Sloan, Y. G. Wang, and R. S. Womersley.  Random Point Sets on the Sphere — Hole Radii, Covering, and Separation.  _Experimental Mathematics_ , 27(1):62–81, 12 2015\.  ISSN 1944950X.  doi: 10.1080/10586458.2016.1226209.  URL https://arxiv.org/abs/1512.07470v2.
  * Burt [1992] R. S. Burt.  _Structural holes: The social structure of competition_.  Harvard University Press, Cambridge, MA, 1992.
  * Buskens and Van de Rijt [2008] V. Buskens and A. Van de Rijt.  Dynamics of networks if everyone strives for structural holes.  _American Journal of Sociology_ , 114(2):371–407, 2008.
  * Cai et al. [2013] T. Cai, J. Fan, and T. Jiang.  Distributions of Angles in Random Packing on Spheres.  _Journal of Machine Learning Research_ , 14:1837–1864, 2013\.
  * Chamberlain et al. [2017] B. P. Chamberlain, J. R. Clough, and M. P. Deisenroth.  Neural Embeddings of Graphs in Hyperbolic Space.  In _13th international workshop on mining and learning from graphs held in conjunction with KDD_ , 2017.
  * Chami et al. [2019] I. Chami, R. Ying, C. Ré, and J. Leskovec.  Hyperbolic Graph Convolutional Neural Networks.  _Advances in Neural Information Processing Systems_ , 32, 10 2019\.  ISSN 10495258.  doi: 10.48550/arxiv.1910.12933.  URL https://arxiv.org/abs/1910.12933v1.
  * Csardi et al. [2006] G. Csardi, T. Nepusz, et al.  The igraph software package for complex network research.  _InterJournal, complex systems_ , 1695(5):1–9, 2006.
  * Durante et al. [2013] F. Durante, J. Fernández-Sánchez, and C. Sempi.  A topological proof of Sklar’s theorem.  _Applied Mathematics Letters_ , 26(9):945–948, 9 2013.  ISSN 0893-9659.  doi: 10.1016/J.AML.2013.04.005.
  * Farooq et al. [2019] H. Farooq, Y. Chen, T. T. Georgiou, A. Tannenbaum, and C. Lenglet.  Network curvature as a hallmark of brain structural connectivity.  _Nature Communications 2019 10:1_ , 10(1):1–11, 10 2019.  ISSN 2041-1723.  doi: 10.1038/s41467-019-12915-x.  URL https://www.nature.com/articles/s41467-019-12915-x.
  * Fearnhead and Rigaill [2016] P. Fearnhead and G. Rigaill.  Changepoint Detection in the Presence of Outliers.  _Journal of the American Statistical Association_ , 114(525):169–183, 9 2016.  ISSN 1537274X.  doi: 10.48550/arxiv.1609.07363.  URL https://arxiv.org/abs/1609.07363v2.
  * Fosdick et al. [2016] B. K. Fosdick, T. H. McCormick, T. B. Murphy, T. L. J. Ng, and T. Westling.  Multiresolution network models.  _Journal of Computational and Graphical Statistics_ , 28(1):185–196, 8 2016.  ISSN 15372715.  doi: 10.48550/arxiv.1608.07618.  URL https://arxiv.org/abs/1608.07618v5.
  * Frechet [1957] M. Frechet.  Sur la distance de deux lois de probabilite.  _CR. Acad Sci. Paris_ , 244, 1957.
  * Fu et al. [2020] A. Fu, B. Narasimhan, and S. Boyd.  CVXR: An R package for disciplined convex optimization.  _Journal of Statistical Software_ , 94(14):1–34, 11 2020.  ISSN 15487660.  doi: 10.18637/jss.v094.i14.  URL https://CRAN.R-project.
  * Gallot et al. [2004] S. Gallot, D. Hulin, and J. Lafontaine.  _Riemannian Geometry_.  Universitext. Springer Berlin Heidelberg, Berlin, Heidelberg, 2004.  ISBN 978-3-540-20493-0.  doi: 10.1007/978-3-642-18855-8.  URL http://link.springer.com/10.1007/978-3-642-18855-8.
  * Gilbert and Jain [2017] A. C. Gilbert and L. Jain.  If it ain’t broke, don’t fix it: Sparse metric repair.  In _2017 55th Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pages 612–619. IEEE, 2017.
  * Gromov [2007] M. Gromov, editor.  _Metric Structures for Riemannian and Non-Riemannian Spaces_.  Birkhäuser Boston, 2007.  doi: 10.1007/978-0-8176-4583-0.
  * Gu et al. [2018] A. Gu, F. Sala, B. Gunel, and C. Ré.  Learning mixed-curvature representations in product spaces.  In _International Conference on Learning Representations_ , 2018.
  * Harchaoui and Lévy-Leduc [2010] Z. Harchaoui and C. Lévy-Leduc.  Multiple change-point estimation with a total variation penalty.  _Journal of the American Statistical Association_ , 105(492):1480–1493, 2010.
  * Hoeffding [1940] W. Hoeffding.  Masstabvariarte Korrelationstheorie.  _Schrijl Math. Inst. Univ. Berlin_ , 5(6), 1940.
  * Hoff et al. [2002] P. D. Hoff, A. E. Raftery, and M. S. Handcock.  Latent space approaches to social network analysis.  _Journal of the american Statistical association_ , 97(460):1090–1098, 2002.
  * Killing [1891] W. Killing.  Ueber die Clifford-Klein’schen Raumformen.  _Mathematische Annalen 1891 39:2_ , 39(2):257–278, 6 1891.  ISSN 1432-1807.  doi: 10.1007/BF01206655.  URL https://link.springer.com/article/10.1007/BF01206655.
  * Klingenberg [1995] W. Klingenberg.  _Riemannian Geometry_.  Walter de Gruyter, 1 edition, 12 1995.  doi: 10.1515/9783110905120.
  * Le Cam [1956] L. Le Cam.  On the asymptotic theory of estimation and testing hypotheses.  In _Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Contributions to the Theory of Statistics_ , volume 3, pages 129–157. University of California Press, 1956.
  * Leal et al. [2018] W. Leal, G. Restrepo, P. F. Stadler, and J. Jost.  Forman-Ricci Curvature for Hypergraphs.  _Advances in Complex Systems_ , 24(1), 11 2018.  doi: 10.13140/RG.2.2.27347.84001.  URL http://arxiv.org/abs/1811.07825http://dx.doi.org/10.13140/RG.2.2.27347.84001.
  * Lee et al. [2019] W. Lee, T. H. Mccormick, J. Neil, M. Cole, S. Microsoft, and Y. Cui.  Anomaly Detection in Large Scale Networks with Latent Space Models.  _Technometrics_ , pages 1–23, 11 2019.  ISSN 0040-1706.  doi: 10.1080/00401706.2021.1952900.  URL https://arxiv.org/abs/1911.05522v2.
  * Leskovec et al. [2007] J. Leskovec, J. Kleinberg, and C. Faloutsos.  Graph evolution.  _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , 1(1), 3 2007.  ISSN 15564681.  doi: 10.1145/1217299.1217301.  URL https://dl.acm.org/doi/abs/10.1145/1217299.1217301.
  * Li and Dunson [2019] D. Li and D. B. Dunson.  Geodesic distance estimation with spherelets.  _arXiv preprint arXiv:1907.00296_ , 2019.
  * Li et al. [2017] D. Li, M. Mukhopadhyay, and D. B. Dunson.  Efficient manifold and subspace approximations with spherelets.  _arXiv preprint arXiv:1706.08263_ , 2017.
  * Lok et al. [2021] T. Lok, J. Ng, T. B. Murphy, T. Westling, T. H. Mccormick, and B. Fosdick.  Modeling the social media relationships of Irish politicians using a generalized latent space stochastic blockmodel.  _Annals of Applied Statistics_ , 15(4):1923–1944, 12 2021.  ISSN 1932-6157.  doi: 10.1214/21-AOAS1483.
  * Lubold et al. [2023] S. Lubold, A. G. Chandrasekhar, and T. H. McCormick.  Identifying the latent space geometry of network models through analysis of curvature.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 2023.
  * MacDonald et al. [2022] P. W. MacDonald, E. Levina, and J. Zhu.  Latent space models for multiplex networks with shared structure.  _Biometrika_ , 109(3):683–706, 2022.
  * Nagano et al. [2019] Y. Nagano, S. Yamaguchi, Y. Fujita, and M. Koyama.  A Wrapped Normal Distribution on Hyperbolic Space for Gradient-Based Learning.  _36th International Conference on Machine Learning, ICML 2019_ , 2019-June:8242–8251, 2 2019.  doi: 10.48550/arxiv.1902.02992.  URL https://arxiv.org/abs/1902.02992v2.
  * Ni et al. [2019] C. C. Ni, Y. Y. Lin, F. Luo, and J. Gao.  Community Detection on Networks with Ricci Flow.  _Scientific Reports 2019 9:1_ , 9(1):1–12, 7 2019\.  ISSN 2045-2322.  doi: 10.1038/s41598-019-46380-9.  URL https://www.nature.com/articles/s41598-019-46380-9.
  * Nickel and Kiela [2017] M. Nickel and D. Kiela.  Poincaré Embeddings for Learning Hierarchical Representations.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Ollivier [2007] Y. Ollivier.  Ricci curvature of Markov chains on metric spaces.  _Journal of Functional Analysis_ , 256(3):810–864, 1 2007.  ISSN 00221236.  doi: 10.48550/arxiv.math/0701886.  URL https://arxiv.org/abs/math/0701886v4.
  * Papadopoulos et al. [2018] L. Papadopoulos, M. A. Porter, K. E. Daniels, and D. S. Bassett.  Network analysis of particles and grains.  _Journal of Complex Networks_ , 6(4):485–565, 8 2018.  ISSN 20511329.  doi: 10.1093/COMNET/CNY005.  URL https://academic.oup.com/comnet/article/6/4/485/4959635.
  * Pennec [1999] X. Pennec.  Probabilities and statistics on Riemannian manifolds: Basic tools for geometric measurements.  In _International Workshop on Nonlinear Signal and Image Processing_ , pages 194–198, Antalya, Turkey, 6 1999.  URL http://www-sop.inria.fr/epidaure/personnel/pennec/pennec.html.
  * Politis and Romano [1994] D. N. Politis and J. P. Romano.  Large Sample Confidence Regions Based on Subsamples under Minimal Assumptions.  _https://doi.org/10.1214/aos/1176325770_ , 22(4):2031–2050, 12 1994.  ISSN 0090-5364.  doi: 10.1214/AOS/1176325770.
  * Salter-Townshend and McCormick [2017] M. Salter-Townshend and T. H. McCormick.  Latent space models for multiview network data.  _Annals of Applied Statistics_ , 11(3):1217–1244, 9 2017.  ISSN 1932-6157.  doi: 10.1214/16-AOAS955.  URL https://projecteuclid.org/journals/annals-of-applied-statistics/volume-11/issue-3/Latent-space-models-for-multiview-network-data/10.1214/16-AOAS955.full.
  * Samal et al. [2021] A. Samal, H. K. Pharasi, S. J. Ramaia, H. Kannan, E. Saucan, J. Jost, and A. Chakraborti.  Network geometry and market instability.  _Royal Society Open Science_ , 8(2), 2 2021.  ISSN 20545703.  doi: 10.1098/RSOS.201734.  URL https://royalsocietypublishing.org/doi/10.1098/rsos.201734.
  * Sandhu et al. [2015] R. Sandhu, T. Georgiou, E. Reznik, L. Zhu, I. Kolesov, Y. Senbabaoglu, and A. Tannenbaum.  Graph Curvature for Differentiating Cancer Networks.  _Nature Scientific Reports_ , 5(1):1–13, 7 2015\.  ISSN 2045-2322.  doi: 10.1038/srep12323.  URL https://www.nature.com/articles/srep12323.
  * Sandhu et al. [2016] R. S. Sandhu, T. T. Georgiou, and A. R. Tannenbaum.  Ricci curvature: An economic indicator for market fragility and systemic risk.  _Science Advances_ , 2(5), 5 2016.  ISSN 23752548.  doi: 10.1126/SCIADV.1501495/SUPPL˙FILE/1501495˙SM.PDF.  URL https://www.science.org/doi/10.1126/sciadv.1501495.
  * Saucan et al. [2020] E. Saucan, A. Samal, and J. Jost.  A Simple Differential Geometry for Complex Networks.  _Network Science_ , 9(S1):S106–S133, 4 2020.  ISSN 20501250.  doi: 10.48550/arxiv.2004.11112.  URL https://arxiv.org/abs/2004.11112v2.
  * Schoenberg [1935] I. J. Schoenberg.  Remarks to Maurice Frechet’s Article “Sur La Definition Axiomatique D’Une Classe D’Espace Distances Vectoriellement Applicable Sur L’Espace De Hilbert.  _The Annals of Mathematics_ , 36(3):724, 7 1935\.  ISSN 0003486X.  doi: 10.2307/1968654.
  * Sia et al. [2019] J. Sia, E. Jonckheere, and P. Bogdan.  Ollivier-Ricci Curvature-Based Method to Community Detection in Complex Networks.  _Nature, Scientific Reports 2019 9:1_ , 9(1):1–12, 7 2019.  ISSN 2045-2322.  doi: 10.1038/s41598-019-46079-x.  URL https://www.nature.com/articles/s41598-019-46079-x.
  * Sklar [1959] A. Sklar.  Fonctions de repartition a n dimensions et leurs marges.  _Publications de l’Institut de statistique de l’Université de Paris_ , 8:229–231, 1959.  URL https://ci.nii.ac.jp/naid/10011938360.
  * Smith et al. [2017] A. L. Smith, D. M. Asta, and C. A. Calder.  The Geometry of Continuous Latent Space Models for Network Data.  _Statistical Science_ , 34(3):428–453, 12 2017\.  ISSN 21688745.  doi: 10.48550/arxiv.1712.08641.  URL https://arxiv.org/abs/1712.08641v2.
  * Sweet and Adhikari [2020] T. Sweet and S. Adhikari.  A Latent Space Network Model for Social Influence.  _Psychometrika_ , 85(2):251–274, 6 2020.  ISSN 18600980.  doi: 10.1007/S11336-020-09700-X/FIGURES/11.  URL https://link.springer.com/article/10.1007/s11336-020-09700-x.
  * van der Hoorn et al. [2020] P. van der Hoorn, W. J. Cunningham, G. Lippner, C. Trugenberger, and D. Krioukov.  Ollivier-Ricci curvature convergence in random geometric graphs.  _Physical Review Research_ , 3(1), 8 2020.  doi: 10.1103/PhysRevResearch.3.013211.  URL http://arxiv.org/abs/2008.01209http://dx.doi.org/10.1103/PhysRevResearch.3.013211.
  * van der Vaart [1998] A. van der Vaart.  _Asymptotic Statistics_.  Cambridge University Press, 10 1998.  doi: 10.1017/cbo9780511802256.  URL /core/books/asymptotic-statistics/A3C7DAD3F7E66A1FA60E9C8FE132EE1D.
  * Volz et al. [2011] E. M. Volz, J. C. Miller, A. Galvani, and L. Meyers.  Effects of Heterogeneous and Clustered Contact Patterns on Infectious Disease Dynamics.  _PLOS Computational Biology_ , 7(6):e1002042, 2011\.  ISSN 1553-7358.  doi: 10.1371/JOURNAL.PCBI.1002042.  URL https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002042.
  * Zhang et al. [2021] S. Zhang, Y. Tay, W. Jiang, D.-c. Juan, and C. Zhang.  Switch spaces: learning product spaces with sparse gating.  _arXiv preprint arXiv:2102.08688_ , 2021.

##  Appendix A Proofs of Theorems

###  A.1 Proof of Lemma 1

###### Proof.

Consider first a set of 3 points on $(x,y,z)\in\mathbb{S}^{p}$ which are not
co-linear. We can define an orthogonal matrix (rotation matrix)
$Q\in\mathbb{R}^{p+1\times p+1}$ which we will first construct a rotational
isometry, then use a result which exploits fixed point isometries.

######  Theorem 11 (Theorem 1.10.15 in Klingenberg [1995]).

If $f:(\overline{M},\overline{g})\to(\overline{M},\overline{g})$ is an
isometry on a Riemannian manifold, then the fixed point set of $f$ forms a
totally geodesic submanifold.

WLOG let $x=(1,0,\dots,0)$. We define the orthogonal rotation matrix $Q$ via
the following. Let $q_{i}$ denote the $i^{th}$ un-normalized column vector of
$Q$. We the define the following first $3$ basis functions

$\displaystyle q_{1}$ $\displaystyle=[1,0,\dots,0]^{\top}$ $\displaystyle
q_{2}$ $\displaystyle=y-(y^{\top}q_{1})q_{1}$ $\displaystyle q_{3}$
$\displaystyle=z-\frac{(z^{\top}q_{2})}{||{q_{2}}||_{2}}q_{2}-(z^{\top}q_{1})q_{1}$

Then normalizing each of these vectors to construct the columns of $Q$ and
completing the orthonormal basis with any remaining orthogonal basis vectors
of $\mathbb{R}^{p}$, and we have constructed an orthogonal matrix $Q$ and
$Q^{\top}Q=I=QQ^{\top}$. Thus we have defined an isomorphism
$f_{1}(\cdot):\mathbb{S}^{p}\mapsto\mathbb{S}^{p}$ where

$\displaystyle f_{1}(x)$ $\displaystyle=Q^{\top}x$ $\displaystyle
d(f_{1}(x),f_{1}(y))$
$\displaystyle=\frac{1}{\sqrt{\kappa}}\mathrm{acos}\left(f_{1}(x)^{\top}f_{1}(y)\right)$
$\displaystyle=\frac{1}{\sqrt{\kappa}}\mathrm{acos}\left(x^{\top}QQ^{\top}y\right)$
$\displaystyle=\frac{1}{\sqrt{\kappa}}\mathrm{acos}\left(x^{\top}y\right)$
$\displaystyle=d(f_{1}(x),f_{1}(y))$

hence this rotation is an isomorphism. Then let
$f_{2}(x)=(x_{0},x_{1},x_{2},-x_{3},-x_{4},\cdot,-x_{p})$ denote an additional
isomorphism. Since the composition of isomorphisms is also an isomorphism, we
can define a totally geodesic submanifold of dimension $2$ which consists of
the solution set to $Q^{\top}=[c_{1},c_{2},c_{3},0,0,\dots,0]^{\top}$ for some
$c_{1},c_{2},c_{3}$ by Theorem 11. Then clearly, $x,y,z$ are all fixed points
of this isomorphism by construction, and hence any three points can be
contained within a totally geodesic submanifold of dimension $2$ of the
corresponding curvature.

The proofs for $\mathbb{E}^{p}$ and $\mathbb{H}^{p}$ follows this argument
identically and thus the proof is complete. ∎

###  A.2 Proof of Theorem 2

###### Proof.

Any 3 points $(x,y,z)$ can be embedded isometrically in a totally geodesic
sub-manifold of dimension $2$ when $\mathcal{M}^{p}(\kappa)$ is a canonical
manifold due to Lemma 1.

Case 1: Spherical. Under the spherical model, first place the midpoint at the
origin. Next, given $d_{yz}$ place point $y$ at the point $(y_{0},y_{1},0)$
where $y_{1}=\sqrt{1-y_{0}^{2}}$. Then
$y_{0}=\cos(\sqrt{\kappa}\frac{d_{yz}}{2})$. Next we place point $z$ at
$z=(z_{0},z_{1},0)$ where $z_{0}=y_{0}$ and $z_{1}=-y_{1}$. Third, we solve
for the locations of $x$ such that $d_{xm},d_{xy},d_{xz}$ are satisfied.

We continue with $x=(x_{0},x_{1},x_{2})$. We can set
$x_{0}=\cos(\sqrt{\kappa}\frac{d_{xm}}{2})$. Next, $x_{1}$ can be solved by
letting

$\displaystyle\cos(\sqrt{\kappa}d_{xy})$ $\displaystyle=x_{0}y_{0}+x_{1}y_{1}$
$\displaystyle=\cos(\sqrt{\kappa}d_{xm})\cos(\sqrt{\kappa}d_{yz}/2)+x_{1}\sin(\sqrt{\kappa}d_{yz}/2)$

Similarly

$\displaystyle\cos(\sqrt{\kappa}d_{xz})$ $\displaystyle=x_{0}z_{0}+x_{1}z_{1}$
$\displaystyle=\cos(\sqrt{\kappa}d_{xm})\cos(\sqrt{\kappa}d_{yz}/2)-x_{1}\sin(\sqrt{\kappa}d_{yz}/2)$

solving for $x_{1}$ results in:

$\displaystyle\frac{\cos(\sqrt{\kappa}d_{xy})-\cos(\sqrt{\kappa}d_{xm})\cos(\sqrt{\kappa}d_{yz}/2)}{\cos(\sqrt{\kappa}d_{xz})-\cos(\sqrt{\kappa}d_{xm})\cos(\sqrt{\kappa}d_{yz}/2)}$
$\displaystyle=-1$

which can be further rearranged to

$\displaystyle
2\cos(d_{xm}\sqrt{\kappa})-\text{Sec}(\frac{d_{yz}}{2}\sqrt{\kappa})(\cos(d_{xy}\sqrt{\kappa})+\cos(d_{xz}\sqrt{\kappa}))$
$\displaystyle=0$

We then normalize this by the curvature value $\frac{1}{\kappa}$, which allows
for $g(\kappa,d)$ to be a continuous function of $\kappa$ from the hyperbolic
$\kappa<0$ to spherical space $\kappa>0$.

Case 2: Hyperbolic. The proof of the hyperbolic case follows an identical
method to the spherical case and is omitted for brevity. However this results
in the estimating function:

$\displaystyle\frac{1}{-\kappa}\bigg{(}2\cosh(d_{xm}\sqrt{-\kappa})$
$\displaystyle-\text{Sech}(\frac{d_{yz}}{2}\sqrt{-\kappa})(\cosh(d_{xy}\sqrt{-\kappa})+\cosh(d_{xz}\sqrt{-\kappa}))\bigg{)}=0$

Remark: Under the limit as $\kappa\to 0$ we find that $\lim_{\kappa\to
0}g(\kappa,d)=\frac{1}{2}d_{xy}^{2}+\frac{1}{2}d_{xz}^{2}-\frac{1}{4}d_{yz}^{2}-d_{xm}^{2}$
which gives exactly the parallelogram law in Euclidean space. This also
highlights the necessity that the term $\frac{1}{\kappa}$ plays in maintaining
a smooth equation as a function of $\kappa$ through $0$. ∎

###  A.3 Proof of Theorem 3

###### Proof.

We will use the implicit function theorem to construct an implicit function
for which we can later apply a delta method. We first note that $g(\kappa,d)$
is continuously differentiable when
$\sqrt{\kappa}<\frac{\mathrm{acos}(\pi)}{d_{yz}}$. This boundary, however
corresponds to the maximum distance allowed on a sphere, as given by two anti-
polar points. Therefore to apply the implicit function theorem, what remains
is

$\displaystyle\frac{\partial}{\partial\kappa}g(\kappa,d)$
$\displaystyle\not=0.$

which will hold by a brief application of (B3). As we have derived in
Theorem2, $d_{xm}(\kappa;\boldsymbol{d}(x,y,z))$ is a differentiable,
continuous function.

We note that since by definition, the function

$g(\kappa,d_{xy},d_{xz},d_{yz},d_{xm}(\kappa;\boldsymbol{d}(x,y,z)))=0$

Therefore, it is equivalent for a fixed $d_{xm}$ to write the estimating
function as a function of the equivalent exact midpoint $d_{xm}(\kappa)$111We
drop the $\boldsymbol{d}(x,y,z)$ for brevity

$\displaystyle g(\kappa;d_{xy},d_{xz},d_{yz},d_{xm})$
$\displaystyle=\frac{2\cos(d_{xm}\sqrt{\kappa})}{\kappa}-\frac{2\cos(d_{xm}(\kappa)\sqrt{\kappa})}{\kappa}$

Computing the derivative as a function of $\kappa$

$\displaystyle\frac{d}{d\kappa^{\prime}}g(\kappa^{\prime};d_{xy},d_{xz},d_{yz},d_{xm})\bigg{|}_{\kappa^{\prime}=\kappa}$
$\displaystyle=2\frac{-d_{xm}\sqrt{\kappa}\sin(d_{xm}\sqrt{\kappa})+2\cos(d_{xm}\sqrt{\kappa})}{2\kappa^{2}}$
$\displaystyle+2\frac{(d_{xm}(\kappa)\sqrt{\kappa}+\kappa^{3/2}d^{\prime}_{xm}(\kappa))\sin(d_{xm}\sqrt{\kappa})-2\cos(d_{xm}(\kappa)\sqrt{\kappa})}{2\kappa^{2}}$

Since $g(\kappa,\boldsymbol{d})=0\implies d_{xm}=d_{xm}(\kappa)$,

$\displaystyle\frac{d}{d\kappa^{\prime}}g(\kappa^{\prime};d_{xy},d_{xz},d_{yz},d_{xm})\bigg{|}_{\kappa^{\prime}=\kappa}$
$\displaystyle=2d^{\prime}_{xm}(\kappa)\frac{\sin(d_{xm}(\kappa)\sqrt{\kappa})}{\sqrt{\kappa}}$

Since $\kappa\leq\frac{\pi^{2}}{d_{xm}^{2}}$ then if
$d^{\prime}_{xm}(\kappa)>0$ then $\frac{d}{d\kappa}g(\kappa,\boldsymbol{d})>0$

Next, by the implicit function theorem, there exists an open neighborhood
$\mathcal{D}\in\mathbb{R}^{4}$ and a function $f(d)$ such that $\kappa=f(d)$
and $g(f(d),d)$. Furthermore:

$\displaystyle\nabla_{d}f(d)$ $\displaystyle=-\left(\frac{\partial
g(\kappa,d)}{\partial\kappa}\right)^{-1}\left[\nabla_{d}g(\kappa,d)\right]$

Therefore by the delta method:

$\displaystyle\sqrt{r(n)}(\widehat{\kappa}-\kappa)$
$\displaystyle=\sqrt{r(n)}(f(\widehat{d})-f(d))$
$\displaystyle\to_{W}N(0,\nabla g(d)^{\top}\Sigma\nabla g(d))$

And by the implicit function theorem:

$\displaystyle\nabla g(d)^{\top}\Sigma\nabla g(d)$
$\displaystyle=\left(\frac{\partial
g(\kappa,d)}{\partial\kappa}^{2}\right)^{-1}\left[\nabla_{d}g(\kappa,d)^{\top}\Sigma\nabla_{d}g(\kappa,d)\right]$

If instead $||{\widehat{d}-d}||_{2}=o_{P}(r(n)^{-1/2})$, then we can use the
continuous mapping theorem instead and we have
$|\widehat{\kappa}-\kappa|=o_{P}(r(n)^{-1/2})$. Remark: This form is very
similar to the usual standard asymptotic normality proofs for $Z$ estimators
as in van der Vaart [1998]. However, the main difference is in the fact that
we are not averaging the estimating function, but rather plugging in a
distance estimate which is asymptotically normal.

We find that in practice, condition (B3) holds, however, for some distances
that are plugged in, particularly those which may not satisfy a midpoint
exactly may result in the root of the estimating function be the root of its
partial derivative. In fact, in the Euclidean case, according to the Taylor
series expansion at $\kappa=0$

$\displaystyle g(\kappa;d)$
$\displaystyle=(\frac{d_{xy}^{2}}{2}+\frac{d_{xz}^{2}}{2}-d_{xm}^{2}-\frac{d_{yz}^{2}}{4})+$
$\displaystyle\left(\frac{d_{xm}^{2}}{12}-\frac{5d_{yz}^{2}}{192}+\frac{1}{16}d_{yz}^{2}(d_{xy}^{2}+d_{xz}^{2})-\frac{1}{24}(d_{xy}^{2}+d_{xz}^{2})\right)\kappa$
$\displaystyle+\mathcal{O}(\kappa^{2}).$

Clearly at $\kappa=0$ the solution reduces to the parallelogram law
$(\frac{d_{xy}^{2}}{2}+\frac{d_{xz}^{2}}{2}-d_{xm}^{2}-\frac{d_{yz}^{2}}{4}=0)$.
When we substitute in the corresponding solution at $\kappa=0$.

$\displaystyle\frac{\partial}{\partial\kappa}g(\kappa;d)\Bigg{|}_{\kappa=0}$
$\displaystyle=-\frac{1}{48}\left(d_{yz}^{4}-2d_{yz}^{2}(d_{xz}^{2}+(d_{xy}^{2})+(d_{xz}^{2}-d_{xy}^{2})^{2}\right)$
$\displaystyle=-\frac{1}{48}\bigg{(}(d_{yz}+d_{xz}+d_{xy})(d_{yz}-d_{xz}-d_{xy})$
$\displaystyle\times(d_{yz}-d_{xz}+d_{xy})(d_{yz}+d_{xz}-d_{xy})\bigg{)}$

As long as the triangle inequality is satisfied strictly for $x,y,z$, then
$\frac{\partial}{\partial\kappa}g(\kappa;d)\Bigg{|}_{\kappa=0}>0$. However, if
the triangle inequality is not strict, i.e. $x,y,z$ are co-linear, and
$d_{ij}=d_{ik}+d_{kj}$ for $(i,j.k)$ being some permutation of $(x,y,z)$, then
$\frac{\partial}{\partial\kappa}g(\kappa;d)\Bigg{|}_{\kappa=0}=0$. For
$\kappa\not=0$, we do not have simple closed-form expressions for
$\frac{\partial}{\partial\kappa}g(\kappa;d)\Bigg{|}_{\kappa=0}$ and thus we
leave (B3) as an assumption, however, we believe this trend to hold in the
other canonical manifolds.

∎

######  Lemma 12.

Let $x,y,z\in\mathcal{M}^{p}(\kappa)$ and let $m$ denote the midpoint between
$y$ and $z$. Then let $d_{xm}(\kappa;\boldsymbol{d}(x,y,z))$ denote the
distance to the midpoint $m$ from $x$ as a function of the curvature of the
latent space $\kappa$.

Then $d_{xm}(\kappa;\boldsymbol{d}(x,y,z))$ is an increasing function in
$\kappa$.

###### Proof.

The proof follows from Toponogov’s triangle theorem and the negative curvature
extension [Berger, 1962]. Abbreviated for our purposes in Theorem 13, we
simply let $p=z,x=r,q=m$ and compare two manifolds
$\kappa\leq\kappa^{\prime}$.

Then it immediately follows

$d_{xm}(\kappa;\boldsymbol{d}(x,y,z))\leq
d_{xm}(\kappa^{\prime};\boldsymbol{d}(x,y,z)).$

∎

######  Theorem 13 (Toponogov’s Triangle Comparison Theorem (Theorem 3 of
Berger [1962]).

Let $\mathcal{M}$ be a complete Riemannian manifold, and let
$\mathcal{(}\kappa)$ be the simply connected manifold of constant curvature
$\kappa\in\mathbb{R}$. Let $\Delta_{pqr}$ denote a geodesic triangle with
points $p,q,r\in\mathcal{M}$. Let $\widetilde{\Delta}_{pqr}$ is the geodesic
triangle with side lengths $d_{pq}=\widetilde{d}_{pq}$ and
$d_{pr}=\widetilde{d}_{pr}$. If $\underline{K}(\mathcal{M})$ is the minimum
Riemannian sectional curvature on the manifold, and if
$\kappa\leq\underline{K}(\mathcal{M})$, then

$\widetilde{d}_{qr}\leq d_{qr}.$

Assumption (B3) is therefore very mild and that we will always have a non-
decreasing $d_{xm}$ function of $\kappa$ but that we only have to assume that
the increasingness is strict. In practice, this is only not strict when the
points $x,y,z$ are colinear, and therefore the distance to the midpoint does
not change as a function of the curvature since, effectively, these points lie
along a single geodesic.

###  A.4 Proof of Theorem 10

###### Proof.

The proof here is a straightforward application of the plug in estimator of
$\widehat{p}_{XY}$ and the set of random effects $\widehat{\gamma}_{X/Y}$.
Then by the Lindeberg-Feller CLT:

$\displaystyle\sqrt{\ell^{2}\sigma_{\ell}}(\widehat{p}_{XY}-p_{XY})$
$\displaystyle\to_{d}N\left(0,1\right)$

Plugging this in to our definition of distance

$\displaystyle\widehat{d}_{XY}$
$\displaystyle=-\log(\widehat{p}_{XY})+\widehat{\gamma}_{X}+\widehat{\gamma}_{Y}$
$\displaystyle\sqrt{\ell^{2}}\widehat{d}_{XY}$
$\displaystyle=-\sqrt{\ell^{2}}\log(\widehat{p}_{XY})+\sqrt{\ell^{2}}\widehat{\gamma}_{X}+\sqrt{\ell^{2}}\widehat{\gamma}_{Y}$
$\displaystyle=-\sqrt{\ell^{2}}\log(\widehat{p}_{XY})+o_{P}(1)$

And by the delta method:

$\displaystyle\sqrt{\ell^{2}\frac{\sigma_{\ell}}{p_{XY}}}\left(\widehat{d}_{XY}-d_{XY}\right)$
$\displaystyle\to_{d}N\left(0,1\right)$

∎

###  A.5 Proof of Theorem 5

In order to prove Theorem 5 we first prove a series of lemmas.

Let $F$ denote a marginal distribution (in our case that of of $\Xi(m)$) and
let $X_{i}$ be a random variable with distribution function $F$. By Sklar’s
Theorem [Sklar, 1959, Durante et al., 2013] we can express the joint
distribution ($H$) of a vector of correlated $X$, $H$ as

$\displaystyle H(\boldsymbol{x})$ $\displaystyle=P(X_{1}\leq x_{1},X_{2}\leq
x_{2},\dots,X_{n}\leq x_{n})$
$\displaystyle=\mathbb{C}(F(x_{1}),F(x_{2}),\dots,F(x_{n}))$

where $\mathbb{C}$ is a copula for the joint distribution.

######  Lemma 14.

For a set of random variables $X_{i}$ with equal marginal distribution, with
arbitrary correlation structure denoted via copula $\mathbb{C}$, then
$\mathbb{\widehat{M}}(X)=F^{-1}(\mathbb{\widehat{M}}(U))$ where
$U\sim\mathbb{C}$

###### Proof.

We can express $X=F^{-1}(U)$ where $U$ is sampled from the copula
$\mathbb{C}$. Then if $n$ is odd, $X_{((n+1)/2)}$ is the median, and clearly
$X_{(i)}=F^{-1}(U_{(i)})$ with the subscript $(i)$ denoting the $i^{th}$ order
statistic. Now consider $n$ which is even. If $z$ is a median of $U$ then
$z\in[U_{(n/2)},U_{(n/2+1)}]$ then

$\displaystyle z$ $\displaystyle\in[U_{(n/2)},U_{(n/2+1)}]$ $\displaystyle\iff
F^{-1}(z)$ $\displaystyle\in[F^{-1}(U_{(n/2)}),F^{-1}(U_{(n/2+1)})]$
$\displaystyle\iff F^{-1}(z)$ $\displaystyle\in[X_{(n/2)},X_{(n/2+1)}]$

Hence medians of $X$. Hence medians of the underlying uniform distribution are
mapped to the median of the observed variables. ∎

This lemma suggests that if we take the median of the underlying generative
uniforms of the set of correlated random variables.

######  Lemma 15.

Consider a set $A_{I}=\cap_{i\in I}\\{U_{i}\leq z\\}$ where $I$ is a subset of
$\\{1,2,\dots,n\\}$. Consider a sequence of sets $A_{j}$ such that any pair of
$j,j^{\prime}$ contain at least $1$ overlapping index. Then

$\mathbb{C}(\cup A_{j})\leq\mathbb{W}(\cup A_{j})$ (17)

where $\mathbb{W}$ is the perfectly collinear measure corresponding to the
upper Fréchet Hoeffding Bound [Frechet, 1957, Hoeffding, 1940].

###### Proof.

We prove this by induction. Consider the base case $A_{1},A_{2}$. Suppose
there exists some copula $\mathbb{Q}$ such that $\mathbb{Q}(A_{1}\cup
A_{2})>\mathbb{W}(A_{1}\cup A_{2})$. Firstly under $\mathbb{W}$

$\displaystyle\mathbb{W}(A_{1}\cup A_{2})$
$\displaystyle=\mathbb{W}(A_{1})+\mathbb{W}(A_{2})-\mathbb{W}(A_{1}\cap
A_{2})$ $\displaystyle=\min\\{1,z\\}+\min\\{1,z\\}-\min\\{1,z\\}$
$\displaystyle=z+z-z$ $\displaystyle=z$

Then since $\\{A_{j},\cup A_{j^{\prime}}\\}\subset\\{U_{k}\leq z\\}$ for some
index $k$. Then

$\displaystyle\mathbb{Q}(A_{1}\cup A_{2})$
$\displaystyle\leq\mathbb{Q}(U_{k}\leq z)$

However, then $\mathbb{Q}(U_{k}\leq z)>z$ and thus $\mathbb{Q}$ is not a
copula. Now suppose the following holds for any $\mathbb{Q}$.

$\displaystyle\mathbb{Q}(\cup_{j=1}^{n}A_{j})\leq\mathbb{W}(\cup_{j=1}^{n}A_{j})$

for some $n$. Denote $B=\cup_{j=1}^{n}A_{j}$ then we show:

$\displaystyle\mathbb{Q}(B\cup A_{n+1})\leq\mathbb{W}(B\cup A_{n+1})$

Suppose that there exists $\mathbb{Q}$ such that $\mathbb{Q}(B\cup
A_{n+1})>\mathbb{W}(B\cup A_{n+1})$

$\displaystyle\mathbb{W}(B\cup A_{n+1})$
$\displaystyle=\mathbb{W}(B)+\mathbb{W}(A_{n+1})-\mathbb{W}(B\cap A_{n+1})$

Clearly, $\mathbb{W}(B)=z$, $\mathbb{W}(A_{n+1})=z$ and $\mathbb{W}(B\cap
A_{n+1})=z$. Then $\mathbb{Q}(B\cup
A_{n+1})=\mathbb{Q}(B)+\mathbb{Q}(A_{n+1})-\mathbb{Q}(B\cap A_{n+1})$.
Therefore

$\displaystyle\mathbb{Q}(B\cup A_{n+1})$ $\displaystyle\leq
z+\underbrace{\mathbb{Q}(A_{n+1})-\mathbb{Q}(B\cap A_{n+1})}_{\geq 0}$

which generates a contradiction hence

$\displaystyle\mathbb{Q}(B\cup A_{n+1})\leq\mathbb{W}(B\cup A_{n+1})$

∎

######  Lemma 16.

For any copula, $\mathbb{C}$

$\mathbb{C}(\mathbb{\widehat{M}}(U)\leq z)\leq 2z$ (18)

###### Proof.

Define the events $A_{j}$ as in Lemma 15 of size $\lceil{n/2\rceil}$. Then the
event $\\{\mathbb{\widehat{M}}(U)\leq z\\}$ is equal to the union of all such
$A_{j}$ as the median will be equivalent to the case when at least half of all
uniforms are below $z$. We can partition $A_{j}$ into two sets $S_{1},S_{2}$
for which any two pairs of $A_{j}$ in a set have at least one overlapping
index. This can be done by taking all the sets $A_{j}$ containing
$\\{U_{1}\leq z\\}$ in their index and placing them into set $S_{1}$. Then all
other sets must be placed in set $S_{2}$. Since there are $(n-1)$ possible
remaining indices $\\{2,3,\dots,n\\}$ available for $S_{2}$ then any two
events must have an overlapping index( by the pigeonhole principle). Therefore

$\displaystyle\mathbb{C}(\mathbb{\widehat{M}}(U)\leq z)$
$\displaystyle=\mathbb{C}(\cup A_{i})$ $\displaystyle=\mathbb{C}(S_{1}\cup
S_{2})$ $\displaystyle\leq\mathbb{C}(S_{1})+\mathbb{C}(S_{2})\text{ Union
Bound}$ $\displaystyle\leq\mathbb{W}(S_{1})+\mathbb{W}(S_{2})\text{ Lemma
\ref{lem:union_cop}}$ $\displaystyle=2z$

∎

For our theorem, this is all we need, but a corollary immediately follows.
This relates the concentration of any set of arbitrarily correlated uniforms
to the concentration of a singular component element.

######  Corollary 17.

For any copula $\mathbb{C}$.

$\mathbb{C}(|\mathbb{\widehat{M}}(U)-1/2|>\epsilon)\leq 2P(|U-1/2|>\epsilon)$
(19)

###### Proof.

$P(|U-1/2|>\epsilon)$ is simply the marginal distribution of a uniform
distribution.

$\displaystyle P(|U-1/2|>\epsilon)$ $\displaystyle=\max\\{1-2\epsilon,0\\}$

Next we note that

$\displaystyle\\{|\mathbb{\widehat{M}}(U)-1/2|>\epsilon\\}$
$\displaystyle=\\{\mathbb{\widehat{M}}(U)<1/2-\epsilon\\}\cap\\{\mathbb{\widehat{M}}(U)>1/2+\epsilon\\}$

Note that we can define $1-V=U$ where the measure of $V$,
$\widetilde{\mathbb{C}}$ is also a copula, since this is a joint distribution
of marginally uniform variables.

$\displaystyle\\{\mathbb{\widehat{M}}(U)>1/2+\epsilon\\}$
$\displaystyle=\\{\mathbb{\widehat{M}}(V)<1/2-\epsilon\\}$

hence by a union bound.

$\displaystyle\mathbb{C}(|\mathbb{\widehat{M}}(U)-1/2|>\epsilon)$
$\displaystyle\leq\mathbb{C}(\mathbb{\widehat{M}}(U)<1/2-\epsilon)+\widetilde{\mathbb{C}}(\mathbb{\widehat{M}}(V)<1/2-\epsilon)$
$\displaystyle=2(1/2-\epsilon)+2(1/2-\epsilon)$

then since the copula is non-negative

$\displaystyle\mathbb{C}(|\mathbb{\widehat{M}}(U)-1/2|>\epsilon)$
$\displaystyle\leq 2\max\\{1-2\epsilon,0\\}$
$\displaystyle=2P(|U-1/2|>\epsilon)$

∎

######  Lemma 18.

For any set of marginally identically distributed random variables
$\\{X_{i}\\}_{i=1}^{n}$ with CDF $F$ generated with an arbitrary copula
$\mathbb{C}$, then:

$P(\mathbb{\widehat{M}}(X)\leq t)\leq 2P(X_{1}\leq t)$ (20)

###### Proof.

The proof immediately is a result of Lemma 14 and 16.

$\displaystyle P(\mathbb{\widehat{M}}(X)\leq t)$
$\displaystyle=\mathbb{C}(F^{-1}(\mathbb{\widehat{U}})\leq t)\text{ Lemma
\ref{lem:med_uniform_marginal}}$
$\displaystyle=\mathbb{C}(\mathbb{\widehat{U}}\leq F(t))$ $\displaystyle\leq
2P(U\leq F(t))\text{ Lemma \ref{lem:union_cop2}}$ $\displaystyle=2F(t)$
$\displaystyle=2P(X_{1}\leq t)$

∎

This concludes the proof of Theorem 5

###  A.6 Proof of Theorem 4

###### Proof.

Denote the set of events
$M_{i}(t)=\mathbb{\widehat{M}}[\Xi(Z_{\\{1,\dots,i\\}})]>t$ and let
$G=\cap_{i=1}^{K}M_{i}(t)$. We can upper bound the probability that
$\Phi(\mathcal{D}_{K})\leq t$ by the following:

$\displaystyle P(\Phi(\mathcal{D}_{K})\leq t)$
$\displaystyle=P(\Phi(\mathcal{D}_{K})\leq t|\Phi(\mathcal{D}_{K-1})\leq
t,G)P(\Phi(\mathcal{D}_{K-1})\leq t|G)P(G)$
$\displaystyle+P(\Phi(\mathcal{D}_{K})\leq
t|\Phi(\mathcal{D}_{K-1})>t,G)P(\Phi(\mathcal{D}_{K-1})>t|G)P(G)$
$\displaystyle+P(\Phi(\mathcal{D}_{K})\leq t|G^{c})P(G^{c})$
$\displaystyle\leq\underbrace{P(\Phi(\mathcal{D}_{K})\leq
t|\Phi(\mathcal{D}_{K-1})\leq t,G)}_{(I)}P(\Phi(\mathcal{D}_{K-1})\leq
t|G)+\underbrace{P(G^{c})}_{(II)}$

We first consider $(I)$. Since we assume $G$ at each step, at least half of
the midpoints are at a distance $t$ away. Therefore:

$\displaystyle(I)$ $\displaystyle=P(\Phi(\mathcal{D}_{K})\leq
t|\Phi(\mathcal{D}_{K-1})\leq t,G)$
$\displaystyle\leq\max\\{0,1-\binom{K}{2}\alpha\frac{C_{p}}{2}t^{p}\\}$
$\displaystyle\leq\exp(-\binom{K}{2}\alpha\frac{C_{p}}{2}t^{p})$

By recursion, we can apply this such that:

$\displaystyle\prod_{k=1}^{K}P(\Phi(\mathcal{D}_{K})\leq
t|\Phi(\mathcal{D}_{K-1})>t,G)$
$\displaystyle\leq\exp(-\sum_{k=1}^{K}\binom{k}{2}\alpha\frac{C_{p}}{2}t^{p})$
$\displaystyle=\exp(-\binom{K+1}{3}\alpha\frac{C_{p}}{2}t^{p})$
$\displaystyle\leq\exp(-(K+1)^{3}\alpha\frac{C_{p}}{12}t^{p})$

Secondly let we consider term $(II)$. Then

$\displaystyle P(G^{c})$
$\displaystyle=P(\cup_{k=1}^{K}\\{\mathbb{\widehat{M}}[\Xi(Z_{\\{1,\dots,k\\}})]\leq
t\\})$
$\displaystyle\leq\sum_{k=1}^{K}P(\mathbb{\widehat{M}}[\Xi(Z_{\\{1,\dots,k\\}})]\leq
t)$ $\displaystyle\leq K(P(\mathbb{\widehat{M}}[\Xi(Z_{K})]\leq t))$

Since by definition, the median distance to the closest midpoint will be
decreasing as we add more points. Next, by Lemma 16

$\displaystyle P(\\{\mathbb{\widehat{M}}[\Xi(Z_{K})]\leq t\\})$
$\displaystyle\leq 2P(\Xi(Z_{K})\leq t\\})$
$\displaystyle=2P(\cup_{m^{\prime}\not=m}d(m,m^{\prime})\leq t\\})$
$\displaystyle\leq 2\sum_{m^{\prime}\not=m}P(d(m,m^{\prime})\leq t\\})$

And by assumption in our theorem $P(d(m,m^{\prime})\leq t)\leq C_{p}t^{p}$.
Hence we combine these bounds

$\displaystyle P(\Phi(\mathcal{D}_{K})\leq t)$
$\displaystyle\leq\exp(-(K+1)^{3}\alpha\frac{C_{p}}{12}t^{p})+2K\binom{K}{2}t^{p}$

Letting $t=\frac{C_{2}}{K^{3/q}}$ will allow bounding
$P(\Phi(\mathcal{D}_{K})\leq t)$ by any constant and therefore:

$\Phi(\mathcal{D}_{K})=\mathcal{O}_{P}(K^{-3/q})$

∎

###  A.7 Proof of Theorem 6

###### Proof.

The proof is straightforward, let $\boldsymbol{d}(x,y,z)$ denote the vector of
the distances of the triangle $(x,y,z)$, and let
$\boldsymbol{d}(m^{\prime},y,z)$ be the distances in the triangle
$(m^{\prime},y,z)$. Suppose also that we have the distance $d_{xm^{\prime}}$
as this is the distance to the surrogate point, used as a midpoint. Though we
do not have access to the true distance $d_{xm}$, this is a function of the
curvature, and can be written as $d_{xm}(\kappa;\boldsymbol{d}(x,y,z))$. If we
assume that we have a point $m^{\prime}$, so that
$m^{\prime},x,y,z,\in\mathcal{M}^{p}(\kappa)$, then we can also define
$d_{m^{\prime}m}(\kappa;\boldsymbol{d}(m^{\prime},y,z))$. By the triangle
inequality,

$\displaystyle d_{xm}(\kappa;\boldsymbol{d}(x,y,z))$ $\displaystyle\leq
d_{xm^{\prime}}+d_{m^{\prime}m}(\kappa;\boldsymbol{d}(m^{\prime},y,z))$
$\displaystyle d_{xm}(\kappa;\boldsymbol{d}(x,y,z))$ $\displaystyle\geq
d_{xm^{\prime}}-d_{m^{\prime}m}(\kappa;\boldsymbol{d}(m^{\prime},y,z))$

By Lemma 12 each $d_{xm}$ function is monotone in $\kappa$. Therefore, the
number of solutions in $\kappa$ can either be $\\{0,1,\infty\\}$. If $d_{xm}$
never intersect, then it is $0$, if they intersect once, then it is $1$, or if
they remain equal for a region of $\kappa$ then it is $\infty$. This can be
rearranged so that these solutions satisfy.

$\displaystyle
g(\kappa,\boldsymbol{d}(x,y,z)),d_{xm^{\prime}}+d_{m^{\prime}m}(\kappa;\boldsymbol{d}(m^{\prime},y,z)))$
$\displaystyle=0$ $\displaystyle
g(\kappa,\boldsymbol{d}(x,y,z)),d_{xm^{\prime}}+d_{m^{\prime}m}(\kappa;\boldsymbol{d}(m^{\prime},y,z)))$
$\displaystyle=0$

If there is a single solution, then this becomes a single upper or lower
bound. If there are no solutions, then the upper and lower bounds are
$\pm\infty$, and if there are infinite solutions, we may take the most
conservative of these. ∎

###  A.8 Proof of Lemma 7

The proof draws on a similar structure to that of the consistency result of
Lubold et al. [2023], though we highlight the convergence rate achieved under
mild assumptions.

###### Proof.

Let $\epsilon_{\delta}:=\\{\delta<\delta\\}$. We wish to show that:

$\displaystyle\lim_{\ell\to\infty}\frac{P(\epsilon_{\delta}|C_{\ell})}{P(\epsilon^{c}_{\delta}|C_{\ell})}$
$\displaystyle\to\infty$

for $\delta=O(\exp(-\ell))$.

$\displaystyle P(C_{\ell}|\epsilon_{\delta})$
$\displaystyle=\int\int_{\epsilon_{\delta}}\exp(-\binom{\ell}{2}\overline{d}(z))\exp\left(-\binom{\ell}{2}\overline{\nu}\right)dP_{Z|\epsilon_{\delta}}(z)dP_{\boldsymbol{\nu}}(\nu)$
$\displaystyle=\Psi(P_{\nu})\int_{\epsilon_{\delta}}\exp(-\binom{\ell}{2}\overline{d}(z))dP_{Z|\epsilon_{\delta}}(z)$
where
$\displaystyle\Psi(P_{\nu})=\int\exp\left(-\binom{\ell}{2}\overline{\nu}\right)dP_{\boldsymbol{\nu}}(\nu).$

This can be developed analogously for $P(C_{\ell}|\epsilon_{\delta}^{c})$
allowing us to ignore the random effects in this analysis.

Next, consider the ball of radius $\frac{\delta}{2}$ at a point $q$
$B(\delta,q)$. This describes one set of points for which all points have a
maximal distance of $\delta$. Then $\\{Z_{i}\\}_{i=1}^{\ell}\in
B(\delta,q)\implies\epsilon_{\delta}$ for any $q$. Since $f$ is continuous,
then for some $\widetilde{\delta}$ if $f(q)$ has positive probability $f(q)>0$
for some $q$ then for all $\delta\leq\widetilde{\delta}$, $f(x)\geq
c_{2};\forall x\in B(q,\delta)$ for some constant $c$. Therefore
$P(\epsilon_{\delta})\geq P(\\{Z_{i}\\}_{i=1}^{\ell}\in
B(\delta,q))\geq\left(c_{2}\left(\frac{\delta}{2}\right)^{p}\right)^{\ell}$.

Furthermore from Lemma A.1 in Lubold et al. [2023]

$\displaystyle P(C_{\ell}|\epsilon_{\delta}^{c})$
$\displaystyle\leq\exp\left(-\binom{\ell}{2}\mu_{d}\right)$

Putting this all together, for $P(\epsilon_{\delta})\leq\frac{1}{2}$

$\displaystyle\frac{P(\epsilon_{\delta}|C_{\ell})}{P(\epsilon^{c}_{\delta}|C_{\ell})}$
$\displaystyle=\frac{P(C_{\ell}|\epsilon_{\delta})}{P(C_{\ell}|\epsilon^{c}_{\delta})}\frac{P(\epsilon_{\delta})}{1-P(\epsilon_{\delta})}$
$\displaystyle\geq\frac{P(C_{\ell}|\epsilon_{\delta})}{P(C_{\ell}|\epsilon^{c}_{\delta})}\frac{1}{2}P(\epsilon_{\delta})$
$\displaystyle\geq\frac{1}{2}\frac{P(C_{\ell}|\epsilon_{\delta})}{P(C_{\ell}|\epsilon^{c}_{\delta})}c_{2}^{\ell}\left(\frac{\delta}{2}\right)^{p\ell}$
$\displaystyle\geq\frac{1}{2}\frac{\exp(-\binom{\ell}{2}\delta)}{P(C_{\ell}|\epsilon^{c}_{\delta})}c_{2}^{\ell}\left(\frac{\delta}{2}\right)^{p\ell}$
$\displaystyle\geq\frac{1}{2}\frac{\exp(-\binom{\ell}{2}\delta)}{\exp(-\binom{\ell}{2}\mu_{d})}c_{2}^{\ell}\left(\frac{\delta}{2}\right)^{p\ell}$
$\displaystyle=\frac{1}{2}\exp\left(-\binom{\ell}{2}\delta+\binom{\ell}{2}\mu_{d}+\ell\log(c_{2})+\log(\delta/2)p\ell\right)$

Therefore letting $\delta=2\exp(-\widetilde{\mu}_{d}\frac{\ell-1}{p})$ for
$\widetilde{\mu}_{d}<\mu_{d}$. Then

$\displaystyle\lim_{\ell\to\infty}\frac{P(\epsilon_{\delta}|C_{\ell})}{P(\epsilon^{c}_{\delta}|C_{\ell})}\to\infty$

and therefore the proof is complete. ∎

###  A.9 Proof of Lemma 8

At a high level, the proof structure is nearly identical to Lemma 7, however,
we swap the roles of $\nu$ and $d(Z_{i},Z_{j})$.

###### Proof.

Let
$\varepsilon_{\delta}=\left\\{\min_{i\in\\{1,2,\dots,\ell\\}}\nu_{i}\geq-\delta\right\\}$

If $f_{\nu}(\nu)$ is smooth and has left side derivative around $0$ then for
small enough $\delta,F_{\nu}(\nu>-\delta)\geq c_{3}\delta$. As in the previous
method, we similarly take the ratio and integrate out
$\Psi(P_{\overline{d}})$.

Then consider the set $[-\delta,0]$. Then
$P(\varepsilon_{\delta})\geq\left(c_{3}\left(\delta\right)\right)^{\ell}$.

Next we consider the probability of a clique, conditional on the fact that the
locations did not occur in set $\varepsilon^{c}_{\delta}$.

$\displaystyle P(C_{\ell}|\varepsilon^{c}_{\delta})$
$\displaystyle\leq\int\exp\left(2\binom{\ell}{2}\overline{\nu}\right)dP_{\boldsymbol{\nu}}(\boldsymbol{\nu})$

where $\nu d=\frac{1}{\ell}\sum_{i<j}\nu$ Let $\mu_{\nu}=\mathbb{E}[\nu]$.
This is due to the fact $\varepsilon_{\delta}$ denotes the event that
$\min_{i}\nu<-\delta$ it is reasonable that
$P(C_{\ell}|\epsilon^{c}_{\ell})\geq P(C_{\ell})$ since we are excluding the
events most likely to create a clique. Furthermore by the same argument as in
Lemma A.1 in Lubold et al. [2023]

$\displaystyle P(C_{\ell}|\varepsilon^{c}_{\delta})$
$\displaystyle\leq\exp(2\binom{\ell}{2}\mu_{\nu})$

Returning to the computation of the ratio of probabilities,

$\displaystyle\frac{P(\varepsilon_{\delta}|C_{\ell})}{P(\varepsilon^{c}_{\delta}|C_{\ell})}$
$\displaystyle=\frac{P(C_{\ell}|\varepsilon_{\delta})}{P(C_{\ell}|\varepsilon^{c}_{\delta})}\frac{P(\varepsilon_{\delta})}{1-P(\varepsilon_{\delta})}$
$\displaystyle\geq\frac{1}{2}\frac{P(C_{\ell}|\varepsilon_{\delta})}{P(C_{\ell}|\varepsilon^{c}_{\delta})}P(\varepsilon_{\delta})$
$\displaystyle\geq\frac{1}{2}\frac{P(C_{\ell}|\varepsilon_{\delta})}{P(C_{\ell}|\varepsilon^{c}_{\delta})}c_{3}^{\ell}\left(\delta\right)^{\ell}$
$\displaystyle\geq\frac{1}{2}\frac{\exp(-2\binom{\ell}{2}\delta)}{P(C_{\ell}|\varepsilon^{c}_{\delta})}c_{3}^{\ell}\left(\delta\right)^{\ell}$
$\displaystyle\geq\frac{1}{2}\exp\left(-2\binom{\ell}{2}\delta+2\binom{\ell}{2}|\mu_{\nu}|+\ell\log(c_{3})+\log(\delta)\ell\right)$

therefore we can let $\delta=\exp(-|\widetilde{\mu}_{\nu}|2(\ell-1))$ for any
$|\widetilde{\mu}_{\nu}|<|\mu_{\nu}|$ and therefore the proof is complete.

∎

###  A.10 Proof of Lemma 9

###### Proof.

We note that for a given set of nodes with equal latent positions, conditional
on the node of interest appearing in a clique
$\frac{d_{i}}{d_{j}}-\exp(\nu_{i}-\nu_{j})=\mathcal{O}_{P}(n^{-1/2})$. Then
since $\min_{i}\nu_{i}=\mathcal{O}_{P}(\exp(-|\widetilde{\mu}_{\nu}|\ell))$.
All of the random effects are converging to $0$ and hence
$\min_{i}\nu_{i}=\mathcal{O}_{P}\left(\min\left\\{\exp(-|\widetilde{\mu}_{\nu}|\ell),\frac{1}{\sqrt{n}}\right\\}\right))$.
∎

##  Appendix B Additional Computational Details

###  B.1 Newton Method for $\widehat{\kappa}$

Given a set of distances $\widehat{d}$ we can estimate the curvature using a
newton method. Firstly, we compute the derivative of $g(\kappa,d)$ with
respect to $\kappa$.

$\displaystyle\frac{\partial}{\partial\kappa}g(\kappa,d)$
$\displaystyle=(1/(4\kappa^{2}))\bigg{(}8\cos(d_{xm}\sqrt{\kappa})-4\cos(d_{xz}\sqrt{\kappa})\text{sec}(d_{yz}\frac{\sqrt{\kappa}}{2})$
$\displaystyle+4d_{xm}\sqrt{\kappa}\sin(d_{xm}\sqrt{\kappa})$
$\displaystyle-2d_{xy}\sqrt{\kappa}\text{sec}(d_{yz}\frac{\sqrt{\kappa}}{2})\sin(d_{xy}\sqrt{\kappa})$
$\displaystyle-2d_{xz}\sqrt{\kappa}\text{sec}(d_{yz}\frac{\sqrt{\kappa}}{2})\sin(d_{xz}\sqrt{\kappa})$
$\displaystyle+d_{yz}\sqrt{\kappa}\cos(d_{xz}\sqrt{\kappa})\text{sec}(d_{yz}\frac{\sqrt{\kappa}}{2})\text{tan}(d_{yz}\frac{\sqrt{\kappa}}{2})$
$\displaystyle+\cos(d_{xy}\sqrt{\kappa})\text{sec}((d_{yz}\sqrt{\kappa})/2)(-4+d_{yz}\sqrt{\kappa}\text{tan}(d_{yz}\frac{\sqrt{\kappa}}{2})\bigg{)}$

This allows us to construct a newton method for estimating the root
$\widehat{\kappa}$

$\widehat{\kappa}_{(m+1)}=\widehat{\kappa}_{(m)}-\frac{g(\widehat{\kappa}_{(m)},\widehat{d})}{\frac{\partial}{\partial\kappa}g(\widehat{\kappa}_{(m)},\widehat{d})}$

###  B.2 Distance Estimation

We recall the problem of estimating a distance matrix from a set of cliques
$\mathcal{C}$. Though this problem is convex, due to the $O(K^{3})$
restrictions in the problem, it the problem is often slow to reach a solution
in CVXR. Instead, we solve this problem using a successive second order
approximation.

$\displaystyle f(D):=\sum_{X,Y\in\mathcal{C},i\in X,j\in
Y}\bigg{(}A_{ij}\left(\nu_{i}+\nu_{j}-d_{xy}\right)$
$\displaystyle+(1-A_{ij})\log\left(1-\exp\left(\nu_{i}+\nu_{j}-d_{xy}\right)\right)\bigg{)}$
$\displaystyle\approx\sum_{X,Y\in\mathcal{C},i\in X,j\in
Y}\bigg{(}A_{ij}\left(\nu_{i}+\nu_{j}-D_{0,x,y}\right)$
$\displaystyle+(1-A_{ij})\log\left(1-\exp\left(\nu_{i}+\nu_{j}-D_{0,x,y}\right)\right)$
$\displaystyle+\left((A_{ij}-1)\frac{\exp\left(\nu_{i}+\nu_{j}-D_{0,x,y}\right)}{1-\exp\left(\nu_{i}+\nu_{j}-D_{0,x,y}\right)}-A_{ij}\right)\left(d_{xy}-D_{0,x,y}\right)$
$\displaystyle+(A_{ij}-1)\frac{\exp\left(\nu_{i}+\nu_{j}-D_{0,x,y}\right)}{\left(1-\exp\left(\nu_{i}+\nu_{j}-D_{0,x,y}\right)\right)^{2}}\frac{\left(d_{xy}-D_{0,x,y}\right)^{2}}{2}\bigg{)}$
$\displaystyle:=\widetilde{g}(D,D_{0})$

Hence to compute the global solution $\widehat{D}$, we can iteratively solve
the following optimization problem

$\displaystyle\widehat{D}_{t+1}$
$\displaystyle=\text{argsup}_{D\in\in\mathbb{R}^{K\times
K}}\widetilde{g}(D,\widehat{D}_{t})$ $\displaystyle D_{ij}$ $\displaystyle\geq
0\quad\text{ for all }i,j$ $\displaystyle\text{Diag}(D)$ $\displaystyle=0$
$\displaystyle\text{tr}(E_{s}^{\top}D)$ $\displaystyle\geq 0\forall
s\in\mathcal{S}$

This process can be further sped up by choosing a good initialization matrix.
We can do this by using the unconstrained maximum likelihood estimate
$\widehat{D}_{U}$, which is very fast to compute but does not enforce triangle
inequality restrictions. This can be computed analogously as in theorem 10.
Though many of the distances $\widehat{D}_{U}$ may not satisfy the triangle
inequality, we can trim the distances so that $\widehat{D}_{U}$ form a
distance matrix, and use this as the starting point. The Floyd-Warshall
Algorithm is a possible option for constructing a distance matrix from a noisy
matrix which might not have a distance structure. A natural extension to this