Matrix factorisation and the interpretation of geodesic distance
Year: 2022

Abstract:
Given a graph or similarity matrix, we consider the problem of recovering a notion of true distance between the nodes, and so their true positions. We show that this can be accomplished in two steps: matrix factorisation, followed by nonlinear dimension reduction. This combination is effective because the point cloud obtained in the first step lives close to a manifold in which latent distance is encoded as geodesic distance. Hence, a nonlinear dimension reduction tool, approximating geodesic distance, can recover the latent positions, up to a simple transformation. We give a detailed account of the case where spectral embedding is used, followed by Isomap, and provide encouraging experimental evidence for other combinations of techniques.

heading: 1 Introduction
text: Assume we observe, or are given as the result of a computational procedure, data in the form of a symmetric matrix A 2 Rn⇥n which we relate to unobserved vectors Z1, . . . , Zn 2 Z ⇢ Rd, where d ⌧ n, by
A ij = f(Zi, Zj) +E ij , 1  i  j  n, (1)
for some real-valued symmetric function f which will be called a kernel, and a matrix of unobserved perturbations, E 2 Rn⇥n. As illustrative examples, A could be the observed adjacency matrix of a graph with n vertices, a similarity matrix associated with n items, or some matrix-valued estimator of covariance or correlation between n variates.
Broadly stated, our goal is to recover Z1, . . . , Zn given A, up to identifiability constraints. We seek a method to perform this recovery which can be implemented in practice without any knowledge of f , nor of any explicit statistical model for A, and to which we can plug in various matrix-types, be it adjacency, similarity, covariance, or other.
At face value, this may seem rather a lot to ask. Without assumptions, neither f , the Zi’s, nor even d are identifiable [24]. However, under quite general assumptions, we will show that a practical solution is provided by the key combination of two procedures: i) matrix factorisation of A, such as spectral embedding, followed by ii) nonlinear dimension reduction, such as Isomap [57].
This combination is effective because of the following facts, only the first of which is already known. The matrix factorisation step approximates high-dimensional images of the Zi’s, living near a ddimensional manifold [45]. Under regularity and non-degeneracy assumptions on the kernel, which for our analysis is taken to be positive definite, geodesic distance in this manifold equals Euclidean
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
geodesic distance on Z , up to simple transformations of the coordinates, for example scaling. Thus a dimension reduction technique which approximates in-manifold geodesic distances can be expected to recover the Zi’s, up to a simple transformation. If those assumptions fail, as they must with real data, we may still find that those approximate geodesic distances are useful because they reflect a geometry implicit in the kernel.
There are many works which address recovery of the Zi’s when one does consider a particular matrix-type, an explicit statistical model, and/or specific form for f . When A is an adjacency matrix and f(Zi, Zj) is the probability of an edge between nodes i and j, the construct (1) is a latent position model [22] and various scenarios have been studied: Z is a sphere and f(x, y) is a function of hx, yi [10]; estimation by graph distance when f(x, y) is a function of kx yk [11]; spectral estimation when f(x, y) is a possibly indefinite inner-product [53, 55, 12, 46]; inference under a logistic model via MCMC [22] or variational approximation [48]. The case where Z is a discrete set of points corresponds to the very widely studied stochastic block model [23], see [5, 33, 44, 53, 8, 14] and references therein; by contrast the methods we propose are directed at the case where Z is continuous. The opposite problem of estimating f , with Zi unknown but assumed uniform on [0, 1], is known as graphon estimation [19].
Concerning the case when A is a covariance matrix, Latent Variable Gaussian Process models [31, 30, 59] use f(Zi, Zj) to define the population covariance between the ith and jth variates under a hierarchical model. When f(x, y) = hx, yi, this reduces to a latent variable optimisation counterpart of Probabilistic PCA [58] and the maximum likelihood estimate of the Zi’s is obtained from the eigendecomposition of the empirical covariance matrix, which in our setting could be A. When f is nonlinear, but fixed e.g. to a Radial Basis Function kernel, gradient and variational techniques are available [31, 30, 59]. See the same references for onward connections to kernel PCA [50].
Our contribution and precursors. The general practice of using dimensionality reduction and geodesic distance to extract latent features from data is common in data science and machine-learning. Examples arise in genomics [38], neuroscience [63], speech analysis [28, 21] and cyber-security [9]. What sets our work apart from these contributions is that we establish a new rigorous basis for this practice.
It has been understood for several years that the high-dimensional embedding obtained by matrix factorisation can be related via a feature map to the Zi of model (1) [26, 15, 34, 56, 62, 32]. However, it was only recently observed that the embedding must therefore concentrate about a low dimensional set, in the Hausdorff sense [45]. Our key mathematical contribution is to describe the topology and geometry of this set, proving it is a topological manifold and establishing how in-manifold geodesic distance is related to geodesic distance in Z . Riemannian geometry underlying kernels was sketched in [7] but without consideration of geodesic distances or rigorous proofs, and not in the context of latent position estimation. The work [13, 61] was an inspiration to us, suggesting Isomap as a tool for analysing spectral embeddings, under a latent structure model in which Z is a one-dimensional curve in Rd and f is the inner product. A key feature of our problem setup is that f is unknown, in which case the manifold is not available in closed form and typically lives in an infinite-dimensional space. To our knowledge, we are the first to show why spectral embedding followed by Isomap might recover the true Zi’s, or a useful transformation thereof, in general.
We complement our mathematical results with experimental evidence, obtained from both simulated and real data, suggesting that alternative combinations of matrix-factorisation and dimensionreduction techniques work too. For the former, we consider the popular node2vec [20] algorithm, a likelihood-based approach for graphs said to perform matrix factorisation implicitly [43, 64]. For the latter, we consider the popular t-SNE [37] and UMAP [39] algorithms. As predicted by the theory, a direct low-dimensional matrix factorisation, whether using spectral embedding or node2vec, is less successful.

heading: 2 Proposed methods and their rationale
text: 2.1 Spectral embedding, as estimating (Zi)
In our mathematical setup (precise details come later) the kernel f will be nonnegative definite and will be the associated Mercer feature map. It is well known that each point (Zi) then lives in an
infinite-dimensional Hilbert space, which will be denoted `2, and the inner-product h (Zi), (Zj)i2 in this space equals f(Zi, Zj).
The spectral embedding procedure. For p  n, we define the p-dimensional spectral embedding of A to be X̂ = [X̂1, . . . , X̂n]> = Û|Ŝ|1/2 2 Rn⇥p, where |Ŝ| 2 Rp⇥p is a diagonal matrix containing the absolute values of the p largest eigenvalues of A, by magnitude, and Û 2 Rn⇥p is a matrix containing corresponding orthonormal eigenvectors, in the same order. The R packages irlba and RSpectra provide fast solutions which can exploit sparse inputs.
One should think of X̂i as approximating the vector of first p components of (Zi), denoted p(Zi), up to orthogonal transformation, and this can be formalised to a greater or lesser extent depending on what assumptions are made. There are several situations, e.g. f any polynomial [45], the cosine kernel used in Section 4.1, the degree-corrected [29] or mixed-membership [6] stochastic block model, in which only the first p0 (say) components of (·) are nonzero, where typically p0 d. If, after n reaches p0, we embed into p = p0 dimensions, then with k · k denoting the Euclidean norm, we have [18]:
max i2{1,...,n}
kQX̂i p(Zi)k = OP ⇣ (logn)c
n1/2
⌘ , (2)
for a universal constant c 1, orthogonal matrix Q, under regularity assumptions on the Zi’s, f and E (that Zi are i.i.d., f(Zi, Zj) has finite expectation, and the perturbations Eij are independent and centered with exponential tails). This encompasses the case where A is binary, for example a graph adjacency matrix [36, 46]. Similar results are available in the cases where A is a Laplacian [46, 40] or covariance matrix [17]. The methods of this paper are based, in practice, on the distances kX̂i X̂jk, which are invariant to orthogonal transformation and so for the purposes of validating kX̂i X̂jk ⇡ k p(Zi) p(Zj)k the presence of Q in (2) is immaterial.
For X̂i to converge to (Zi) more generally, we must let its dimension p grow with n and, at least given the present state of literature, accept weaker consistency results, for example, convergence in Wasserstein distance between QX̂1, . . .QX̂n and p(Z1), . . . , p(Zn) [32]. Uniform consistency results, in the style of (2), are also available for indefinite [46], bipartite, and directed graphs [25]. These are left for future work because of the complications of Q no longer being orthogonal.
Rank selection. In real data, where n is typically fixed, there is no ‘best’ way of selecting p, as discussed for example in [42]. The method of [65], based on profile-likelihood, provides a popular, practical choice, taking as input the spectrum of A, and is implemented in the R package igraph.

heading: 2.2 Isomap, as estimating Zi
text: We propose to recover Z1, . . . , Zn through the following procedure.
Algorithm 1 Isomap procedure
input p-dimensional points X̂1, . . . , X̂n 1: Compute the neighbourhood graph of radius ✏: a weighted graph connecting i and j, with weight
kX̂i X̂jk, if kX̂i X̂jk  ✏ 2: Compute the matrix of shortest paths on the neighbourhood graph, D̂M 2 Rn⇥n 3: Apply classical multidimensional scaling (CMDS) to D̂M into Rd
return d-dimensional points Ẑ1, . . . , Ẑn
Theorem 1 below establishes, under general assumptions, that M := (Z) is a topological manifold of Hausdorff dimension exactly that of Z (as opposed to an upper bound [45]). It also proposes a ‘change of metric’, often trivial, under which a path on Z and its image on M have the same length.
This result explains how Ẑ1, . . . , Ẑn can estimate Z1, . . . , Zn. First, think of a path from X̂i to X̂j on the neighbourhood graph as a noisy, discrete version of some corresponding continuous path from (Zi) to (Zj) on M. The length of the first (the sum of the weights of the edges) is approximately equal to the length of the second if measured in the standard metric: an infinitesmal step from x to x + dx has length hdx, dxi1/22 . By inversion of we can trace a third path, taking us from Zi to Zj on Z . But to make its length agree with the first two, we must pick a non-standard metric: an infinitesmal step from z to z + dz must be regarded to have length hdz,Hzdzi1/2 = p dz>Hzdz, where
H ij z :=
@ 2 f
@x(i)@y(j)
(z,z)
, i, j 2 {1, . . . , d}.
The exciting news for practical purposes is that, through elementary calculus, we might establish that Hz is constant in z (e.g. if f is translation-invariant) or even proportional to the identity (e.g. if f is just a function of Euclidean distance). In the latter case, if Z is convex, it is not too difficult to see that the length of the shortest path, from (Zi) to (Zj) on M, must be proportional to the Euclidean distance between Zi and Zj . Thus, the matrix of shortest paths obtained by Isomap (Step 2) approximates the matrix of Euclidean distances between the Zi (up to scaling), from which the Zi themselves can be recovered (Step 3), up to scaling, rotation, and translation. Of course, we have taken a few liberties in this argument, such as assuming to be invertible and dz>Hzdz to be positive. We will address these rigorously in the next section.
Dimension and radius selection. For visualisation in our examples we will pick d = 1 or d = 2. For other applications, d can be estimated as the approximate rank of D̂M after double-centering [60], e.g. via [65], again. In practice, we suggest picking ✏ just large enough for the neighbourhood graph to be connected or, if the data have outliers, a fixed quantile (such as 5%) of D̂M, removing nodes outside the largest connected component. The same recommendations apply if the k-nearest neighbour graph is used instead of the ✏-neighbourhood graph.

heading: 3 Theory

heading: 3.1 Setup and assumptions
text: The usual inner-product and Euclidean norm on Rd are denoted h·, ·i and k · k. `2 is the set of x = [x(1) x(2) · · · ]> 2 RN such that kxk2 := (
P1 k=1 |x(k)|2)1/2 < 1. For x, y 2 `2 , hx, yi2 :=P1
k=1 x (k)
y (k). With d 1, let Z be a compact subset of Rd, and let eZ be a closed ball in
Rd, centered at the origin, such that Z ⇢ eZ . Let f : eZ ⇥ eZ ! R be a symmetric, continuous, nonnegative-definite function.
By Mercer’s Theorem, e.g., [51, Thm 4.49], there exist nonnegative real numbers ( k)k 1 and functions (uk)k 1, with each uk : eZ ! R, which are orthonormal with respect to the inner-product (uj , uk) 7! R Z̃ uj(x)uk(x)dx and such that
f(x, y) = 1X
k=1
kuk(x)uk(y), x, y 2 eZ, (3)
where the convergence is absolute and uniform. For x 2 Z let (x) := [ 1/21 u1(x) 1/2 2 u2(x) · · · ]> 2 RN. The image of Z by is denoted M. Observe from (3) that for any x, y 2 Z , f(x, y) = h (x), (y)i2 , k (x)k22 = f(x, x), and the latter is finite for any x 2 Z since f is continuous and Z is compact, hence M ⇢ `2. The following definitions are standard in metric geometry [16]. For any a, b 2 `2, a path in M with end-points a, b is a continuous function : [0, 1] ! M such that 0 = a and 1 = b. With n 1, a non-decreasing sequence t0, . . . , tn such that t0 = 0 and tn = 1, is called a partition. Given a path and a partition P = (t0, . . . , tn), define ( ,P) := Pn k=1 k tk tk 1k2. The length of
(w.r.t. k · k2) is l( ) := supP ( ,P), where the supremum is over all possible partitions. The geodesic distance in M between a and b is defined to be the infimum of l( ) over all paths in M with end-points a, b. DijM denotes geodesic distance in M between (Zi) and (Zj). Similarly a path ⌘ in Z with end-points x, y is a continuous function ⌘ : [0, 1] ! Z such that ⌘0 = x, ⌘1 = y, and with (⌘,P) := Pn k=1 k⌘tk ⌘tk 1k the length of ⌘ (w.r.t. k · k) is l(⌘) := supP (⌘,P). D ij Z denotes geodesic distance in Z between Zi and Zj . If Z is convex, DijZ = kZi Zjk. Assumption 1. For all x, y 2 Z with x 6= y, there exists a 2 Z such that f(x, a) 6= f(y, a). Assumption 2. f is C2 on eZ and for every z 2 Z , the matrix Hz is positive definite.
Assumption 1 is used to show is injective. Assumption 2 has various implications, loosely speaking it ensures a non-degenerate relationship between path-length in M and path-length in Z . Concerning our general setup, if f is not required to be nonnegative definite, but is still symmetric, then a representation formula like (3) is available under e.g. trace-class assumptions [45]. It is of interest to generalise our results to that scenario but technical complications are involved, so we leave it for future research.

heading: 3.2 Path-lengths in M and in Z
text: Theorem 1. Let assumptions 1 and 2 hold. Then is a bi-Lipschitz homeomorphism between Z and M. Let a, b be any two points in M such that there exists a path in M with end-points a, b of finite length. Let be any such path and define ⌘ : [0, 1] ! Z by ⌘t := 1( t). Then ⌘ is a path in Z with l(⌘) < 1. For any ✏ > 0 there exists a partition P✏ such that for any partition P = (t0, . . . , tn) satisfying P✏ ✓ P ,
l( ) nX
k=1
D ⌘tk ⌘tk 1 ,H⌘tk 1 (⌘tk ⌘tk 1) E1/2  ✏. (4)
If and ⌘ are continuously differentiable, the following two equalities hold:
l( ) =
Z 1
0 k̇tk2dt =
Z 1
0 h⌘̇t,H⌘t ⌘̇ti
1/2 dt. (5)
The proof is in the appendix. In the terminology of differential geometry, the collection of innerproducts (x, y) 7! hx,Hzyi, indexed by z 2 Z , constitute a Riemannian metric on Z , see e.g. [41, Ch.1 and p.121] for background, and the right hand side of (5) is the length of the path ⌘ in Z with respect to this Riemannian metric. The theorem tells us that finding geodesic distance in M, i.e. minimising l( ) with respect to for fixed end-points say a = (Zi), b = (Zj), is equivalent to minimizing path-length in Z between end-points 1(a) = Zi, 1(b) = Zj under the Riemannian metric. In section 3.3 we will show how, for various classes of kernels, this minimisation can be performed in closed form leading to explicit relationships between DijM and D ij Z . In Section 3.3 we focus on (5) rather than (4) only for ease of exposition, it is shown in the appendix that exactly the same conclusions can be derived directly from (4). To avoid repetition Assumption 1 will be taken to hold throughout section 3.3 without further mention.

heading: 3.3 Geodesic distance in M and in Z
text: Translation invariant kernels. Suppose Z ⇢ Rd is compact and convex, and f(x, y) = g(x y) where g : Rd ! R is C2. In this situation Hz is constant in z, and equal to the Hessian of g evaluated at the origin. The positive-definite part of Assumption 2 is therefore satisfied if and only if
g has a local maximum at the origin. We now use Theorem 1 to find geodesic distance in M between generic points a, b 2 M for this class of translation-invariant kernels. To do this we obtain a lower bound on l( ) over all paths in M with end-points a, b, then show there exists such a path whose length achieves this lower bound.
Let G = V⇤1/2 where V⇤V> is the eigendecomposition of the Hessian of g evaluated at the origin, so Hz = GG> for all z. Then from (5), for a generic path in M with end-points a, b,
l( ) =
Z 1
0 h⌘̇t,H⌘t ⌘̇ti
1/2 dt =
Z 1
0
⌦ ⌘̇t,GG > ⌘̇t ↵1/2 dt =
Z 1
0 kG>⌘̇tkdt (6)
Z 1
0 G
> ⌘̇tdt = G
> Z 1
0 ⌘̇tdt
= kG >(⌘1 ⌘0)k = kG>[ 1(b) 1(a)]k, (7)
where the inequality is due to the triangle inequality for the k · k norm. The right-most term in (7) is independent of , other than through the end-points a, b. To see there exists a path whose length equals this lower bound, take
⌘̃t := 1(a) + t[ 1(b) 1(a)], t 2 [0, 1], (8)
which is well-defined as a path in Z , since for this class of kernels Z is assumed convex. With ̃t := (⌘̃t), ̃ is clearly a path in M. Differentiating ⌘̃t w.r.t. t and substituting into (6) shows l(̃) = kG>[ 1(b) 1(a)]k, as required. The following proposition summarises our conclusions in the case that a = (Zi) and b = (Zj), for any Zi, Zj .
Proposition 1. If Z is compact and convex, and f(x, y) = g(x y) where g : Rd ! R is C2 and has a local maximum at the origin, then DijM is equal to Euclidean distance between Zi and Zj , up to their linear transformation by G>. In the particular case where g(x y) = h(kx yk2) with h 0(0) < 0, we have DijM / D ij Z = kZi Zjk.
Inner-product kernels. Suppose Z = {x 2 Rd : kxk = 1} and f(x, y) = g(hx, yi), where g : R ! R is C2 and such that g0(1) > 0. In this case Hz = g0(1)I + g00(1)zz>, and by a result of [27], any kernel of the form f(x, y) = g(hx, yi) on the sphere is positive-definite iif g(x) = P1 n=0 anx
n for some nonnegative (an)n 0, implying g00(1) 0. Assumption 2 then holds. To derive the geodesic distance in M, first write out (5) for a generic path in M with end-points a, b:
l( ) =
Z 1
0 h⌘̇t,H⌘t ⌘̇ti
1/2 dt =
Z 1
0
⇣ g 0(1)k⌘̇tk2 + g00(1) |h⌘t, ⌘̇ti|2 ⌘1/2 dt. (9)
Since Z is a radius-1 sphere centered at the origin we must have k⌘tk = 1 for all t, so 0 = 12 d dtk⌘tk 2 = 12 d dt h⌘t, ⌘ti = h⌘t, ⌘̇ti. Therefore the r.h.s. of (9) is in fact equal to g 0(1)1/2 R 1 0 k⌘̇tkdt. This is minimised over all possible paths ⌘ in Z with end-points
1(a), 1(b) when ⌘ is a shortest (with respect to Euclidean distance) circular arc in Z , in which caseR 1 0 k⌘̇tkdt = arccos ⌦ 1(a), 1(b) ↵ . Thus we have:
Proposition 2. If Z = {x 2 Rd : kxk = 1} and f(x, y) = g(hx, yi) where g : R ! R is C2 and such that g0(1) > 0, then DijM = g 0(1)1/2DijZ = g 0(1)1/2 arccos hZi, Zji.
Additive kernels. Suppose Z is the Cartesian product of intervals Zi = [c i , c + i ], i = 1, . . . , d and f(x, y) = Pd
i=1 ↵ifi(x (i) , y (i)), x = [x(1) · · · x(d)]>, where each ↵i > 0 and each fi is positivedefinite and C2. For this class of kernels Hz is diagonal with Hiiz = ↵i @2fi @x(i)@y(i)
(z(i),z(i))
where
z (i) is the ith element of the vector z. The positive-definite part of Assumption 2 thus holds if these diagonal elements are strictly positive for all z 2 Z . Let us introduce the vector of monotonically increasing transformations:
(z) := [ 1(z (1)) · · · d(z(d))]>, i(z(i)) := ↵1/2i
Z z(i)
c i
@ 2 fi
@x(i)@y(i)
1/2
(⇠,⇠)
d⇠, (10)
and write the vector of corresponding inverse transformations 1. For a generic path in M with end-points a, b and as usual ⌘t := 1( t), define ⇣t := (⌘t). Then:
l( ) =
Z 1
0 h⌘̇t,H⌘t ⌘̇ti
1/2 dt =
Z 1
0 k⇣̇tkdt (11) Z 1
0 ⇣̇t dt
= k⇣1 ⇣0k = k 1(b) 1(a)k, (12)
where the second and last equalities hold due to the definition of ⇣t, and denotes elementwise composition. The quantity k 1(b) 1(a)k is thus a lower bound on path-length l( ) for any path in M with end-points a, b. To see that there exists a path whose length achieves this lower bound, and hence that it is the geodesic distance in M between a, b, define ̃t := 1(⇣̃t) where ⇣̃t := 1(a) + t[ 1(b) 1(a)]. Differentiating ⇣̃t w.r.t. t and substituting into (11) yields l(̃) = k 1(b) 1(a)k, as required. Our conclusions in the case that a = (Zi) and b = (Zj) are summarised by: Proposition 3. If Z is the Cartesian product of d intervals Zi ⇢ R, i = 1, . . . , d, and f(x, y) =Pd
i=1 ↵ifi(x (i) , y (i)) where for each i, ↵i > 0, fi is positive definite and C2, and @ 2fi @x(i)@y(i)
(z,z)
>
0 for all z 2 Zi, then DijM is equal to the Euclidean distance between Zi and Zj up to their coordinatewise-monotone transformation by . If ↵i @ 2fi @x(i)@y(i)
(z,z)
is a positive constant over all
z 2 Zi, and over all i = 1, . . . , d, then DijM / D ij Z = kZi Zjk.

heading: 4 Experiments

heading: 4.1 Simulated data: latent position network model
text: This section shows the theory at work in a pedagogical example, where A is an adjacency matrix. Consider an undirected random graph following a latent position model with kernel f(x, y) = ⇢n{cos(x(1) y(1)) + cos(x(2) y(2)) + 2}/4, operating on R2 ⇥ R2. Here, the sequence ⇢n is a sparsity factor which will either be constant, ⇢n = 1, or shrinking to zero sufficiently slowly,
reflecting dense (degree grows linearly) and sparse (degree grows sublinearly) regimes respectively. The kernel is clearly translation-invariant, satisfies Assumptions 1 and 2 and has finite rank, p = 5. The true latent positions Zi 2 R2 are deterministic and equally spaced points on a grid over the region Z = [ ⇡ + 0.25,⇡ 0.25]⇥ [ ⇡ + 0.25,⇡ 0.25], this range chosen to give valid probabilities and an interesting bottleneck in the 2-dimensional manifold M. From Proposition 1, the geodesic distance between (Zi) and (Zj) on M is equal to the Euclidean distance between Zi and Zj , up to scaling, specifically DijM = ⇢nD ij Z/2.
Focusing first on the dense regime, for each n = 100, 400, 1600, 6400, we simulate a graph, and its spectral embedding X̂1, . . . , X̂n into p = 5 dimensions. The first three dimensions of this embedding are shown in Figure 2a), and we see that the points gather about a two-dimensional, curved, manifold. To approximate geodesic distance, we compute a neighbourhood graph of radius ✏ from X̂1, . . . , X̂n, choosing ✏ as small as possible subject to maintaining a connected graph. Figure 2b) shows that approximated geodesic distances roughly agree with the true Euclidean distances in latent space Z (up to the predicted scaling of 1/2), whereas there is significant distortion between ambient and latent distance (kX̂i X̂jk versus kZi Zjk). Finally, we recover the estimated latent positions in R2 using Isomap, which we align with the original by Procrustes (orthogonal transformation with scaling), in Figure 2c). As the theory predicts, the recovery error vanishes as n increases.
In the appendix, we perform the same experiment in a sparse regime, n⇢n = !{log4 n}, chosen to ensure the spectral embedding is still consistent, and the recovery error still shrinks, but more slowly. We also implement other related approaches: UMAP, t-SNE applied to spectral embedding, node2vec directly into two dimensions and node2vec in five dimensions followed by Isomap (Figure 6). We use default configurations for all other hyperparameters. Together, the results support the central recommendation of this paper: to use matrix factorisation (e.g. spectral embedding, node2vec), followed by nonlinear dimension reduction (e.g. Isomap, UMAP, t-SNE).

heading: 4.2 Real data: the global flight network
text: For this example, we use a publically available, clean version [52] (with license details therein) of the crowdsourced OpenSky dataset. Six undirected graphs are extracted, for the months November 2019 to April 2020, connecting any two airports in the world for which at least a single flight was recorded (n ⇡ 10, 000 airports in each). For each graph the adjacency matrix A is spectrally embedded into p = 10 dimensions, degree-corrected by spherical projection [35, 33, 49], after which we apply Isomap with ✏ chosen as the 5% distance quantile.
The dimension p was chosen as a loose upper bound on the dimension estimate returned by the method of [65] on any of the individual graphs, as to facilitate comparison we prefer to use the same dimension throughout (e.g. to avoid artificial differences in variance) and so follow the recommendation of [12] to err on the large side (experiments with different choices of p are in the appendix). In the degreecorrection step we aim to remove network effects related to “airport popularity” (e.g. population, economic and airport size), after which geographic distance might be expected to be the principal factor deciding whether there should be a flight between two airports. Therefore, after applying Isomap, we might hope that Ẑi could estimate the geographical location of airport i. The experiments take a few hours on a standard desktop (and the same is a loose upper bound on compute time for the other experiments in the paper).
The embedding for January is shown in Figures 3)a-b), before and after Isomap. Both recover real geographic features, notably the relative positioning of London, Paris, Madrid, Lisbon, Rome. However, the embedding of the US is warped, and only after using Isomap do we see New York and Los Angeles correctly positioned at opposite extremes of the country. As further confirmation, Figures 3)c-d) show the results restricted to North America, with the points coloured by latitude and longitude. For this continent, the accuracy of longitude recovery is particularly striking.
In line with the recommendation to overestimate p [12], one obtains a similar figure using p = 20 before applying Isomap (Figure 10, appendix), whereas the corresponding figure with p = 2 is too poor to show.
In the appendix, the embeddings of all six months are shown, aligned using two-dimensional Procrustes, showing an important structural change in April. A statistical analysis of inter-continental
geodesic distances suggests the change reflects severe flight restrictions in South America and Africa, at the beginning COVID-19 pandemic (the 11th March 2020 according to the WHO).

heading: 4.3 Real data: correlations between daily temperatures around the globe
text: In this example A is a correlation matrix and we demonstrate a simple model-checking diagnostic informed by our theory. The raw data consist of average temperatures over each day, for several thousand days, recorded in cities around the globe. These data are open source, originate from Berkeley Earth [1] and the particular data set analyzed is from [3]. We used open source latitude and longitude data from [4]. See those references for license details. Removing cities and dates for which data are missing yields a temperature time-series of 1450 common days for each of the n = 2211 cities. A is the matrix of Pearson correlation coefficients between the n time-series.
Figure 4b) shows the spectral embedding of A, with p = 2 and points coloured by the latitude of the corresponding cities. Two visual features are striking: the concentration of points around a curved manifold, and a correspondence between latitude and location on the manifold. Figure 4c) shows latitude against estimated latent positions from Isomap using the k-nearest neighbour graph with k = 200, which is roughly 10% of n, and with d = 1. A clear monotone relationship appears. Our theoretical results can explain this phenomenon. Notice that when d = 1, the additive structure of the kernel in Proposition 3 disappears. Hence Proposition 3 shows that in the case d = 1, for any kernel
(meeting the basic requirements of section 3.1, including Assumptions 1 and 2 of course) DijM is equal to Euclidean distance between Zi and Zj up to their monotone transformation by defined in (10). In turn, this implies that if A did actually follow the model (1) for some f and some “true” latent positions Z1, . . . , Zn to which we had access, we should observe a monotone relationship between those true latent positions and the estimated positions Ẑ1, . . . , Ẑn from Isomap.
The empirical monotonicity in Figure 4c) thus can be interpreted as indicating latitudes of the cities are plausible as “true” latent positions underlying the correlation matrix A, without us having to specify anything in practice about f . In further analysis (appendix) no such empirical monotone relationship is found between longitude and estimated latent position; this indicates longitude does not influence the correlations captured in A. A possible explanation is the daily averaging of temperatures in the underlying data: correlations between these average temperatures may be insensitive to longitude due to the rotation of the earth.

heading: 5 Conclusion
text: Our research shows how matrix-factorisation and dimension-reduction can work together to reveal the true positions of objects based only on pairwise similarities such as connections or correlations. For the sake of exposition and reproducibility, we have used public datasets which can be interpreted without specialist domain knowledge, but the methods are potentially relevant to any scientific application involving large similarity matrices. Thinking about societal impact, our results highlight the depth of information in principle available about individuals, given network metadata, and we hope to raise awareness of these potential issues of privacy.
Concerning the limitations of the methods we have discussed, the bound in (2) indicates that n should be “large” in order for QX̂i to approximate p(Zi) well. n also has an impact on the performance of Isomap: heuristically one needs a high density and a large number of points on or near the manifold to get a good estimate of the geodesic distance. Thus the methods we propose are likely to perform poorly when n is small, corresponding e.g. to a graph with a small number of vertices.
On the other hand, in applications involving large networks or matrices (e.g. cyber-security, recommender systems), the data encountered are often sparse. This is good news for computational feasibility but bad news for statistical accuracy. In particular, for a graph, (2) can only be expected to hold under logarithmically growing degree, the information-theoretic limit below which no algorithm can obtain asymptotically exact embeddings [5]. What this means in practice (as illustrated in our numerical results for the simulated data example in the appendix) is that the manifold may be very hard to distinguish, even for a large graph. Missing data may have a similarly negative impact on discerning manifold structure. There could be substantial estimation gains in better integrating the factorisation and manifold estimation steps to overcome these difficulties.
Another limitation is that our theory restricts attention to the case of positive-definite kernels. When A is say a correlation or covariance matrix the positive-definite assumption on the kernel is of course natural, but when A is say an adjacency matrix, it is a less natural assumption. The implications of removing the positive-definite assumption are the subject of ongoing research.
Funding Transparency Statement
Nick Whiteley and Patrick Rubin-Delanchy’s research was supported by Turing Fellowships from the Alan Turing Institute. Annie Gray’s research was supported by a studentship from Compass, the EPSRC Centre for Doctoral Training in Computational Statistics and Data Science.

Review summaries:
This paper considers the problem of sensing a notion of the true distance between nodes of a  network. This allows to estimate the true position of these nodes up to pairwise scaling. The  cradle of this approach is essentially organized into two steps: a) compute a matrix factorization  of the matrix associated with the graph, i.e., spectral embedding, and b) perform a nonlinear  dimensionality reduction step. The main observation, subject to necessary conditions, is that  Euclidean distances of the spectral embeddings agree with geodesic distances between the  graph nodes. Numerical experiments provide evidence that the proposed approach works. 

This paper considers the problem of finding the latent variables $Z$ given noisy observations of $f(Z_i, Z_j)$ with unknown $f$, and suggests applying matrix factorization (e.g. spectral embedding) + nonlinear dimension reduction (e.g. isomap). It shows mathematically that under some assumptions on the kernel $f$, finding geodesic distance in the manifold (which the matrix factorization step approximates) is equivalent to finding the geodesic distance in the original latent space. Experiments on simulated and real-world data are provided.

The paper considers the problem of recovering a notion of distance between the nodes of a graph. The authors write their main problem very abstractly and provide an excellent description of particular A's, f's, and Z_i's that are of interest in a wide range of applications. The main contribution is to prove that a low-dimensional embedding is related to a topological manifold.  In this regard, it is a follow-up to results by Rubin-Delanchy in [21] from 2020.  The paper seems to be well written, though I was not able to follow all the technical details in the supplementary material as it is pretty far from my expertise.