Title: THE SHAPE OF DATA: INTRINSIC DISTANCE FOR DATA DISTRIBUTIONS

Abstract: The ability to represent and compare machine learning models is crucial in order to quantify subtle model changes, evaluate generative models, and gather insights on neural network architectures. Existing techniques for comparing data distributions focus on global data properties such as mean and covariance; in that sense, they are extrinsic and uni-scale. We develop a first-of-its-kind intrinsic and multi-scale method for characterizing and comparing data manifolds, using a lower-bound of the spectral Gromov-Wasserstein inter-manifold distance, which compares all data moments. In a thorough experimental study, we demonstrate that our method effectively discerns the structure of data manifolds even on unaligned data of different dimensionality, and showcase its efficacy in evaluating the quality of generative models.

1. INTRODUCTION:

The geometric properties of neural networks provide insights about their internals (Morcos et al., 2018; Wang et al., 2018) and help researchers in the design of more robust models (Arjovsky et al., 2017; Bińkowski et al., 2018) . Generative models are a natural example of the need for geometric comparison of distributions. As generative models aim to reproduce the true data distribution P d by means of the model distribution P g (z; Θ), more delicate evaluation procedures are needed. Oftentimes, we wish to compare data lying in entirely different spaces, for example to track model evolution or compare models having different representation space. In order to evaluate the performance of generative models, past research has proposed several extrinsic evaluation measures, most notably the Fréchet (Heusel et al., 2017) and Kernel (Bińkowski et al., 2018) Inception Distances (FID and KID). Such measures only reflect the first two or three moments of distributions, meaning they can be insensitive to global structural problems. We showcase this inadvertence in Figure 1 : here FID and KID are insensitive to the global structure of the data distribution. Besides, as FID and KID are based only on extrinsic properties they are unable to compare unaligned data manifolds. In this paper, we start out from the observation that models capturing the multi-scale nature of the data manifold by utilizing higher distribution moment matching, such as MMD-GAN (Li et al., 2017) and Sphere-GAN (Park & Kwon, 2019), perform consistently better than their single-scale counterparts. On the other hand, using extrinsic information can be misleading, as it is dependent on factors external to the data, such as representation. To address this drawback, we propose IMD, an Intrinsic Multi-scale Distance, that is able to compare distributions using only intrinsic information about the data, and provide an efficient approximation thereof that renders computational complexity nearly linear. We demonstrate that IMD effectively quantifies difference in data distributions in three distinct application scenarios: comparing word vectors in languages with unaligned vocabularies, tracking dynamics of intermediate neural network representations, and evaluating generative models.

3. MULTI-SCALE INTRINSIC DISTANCE:

At the core of deep learning lies the manifold hypothesis, which states that high-dimensional data, such as images or text, lie on a low-dimensional manifold (Narayanan & Mitter, 2010; Belkin & Niyogi, 2002; 2007) . We aim to provide a theoretically motivated comparison of data manifolds based on rich intrinsic information. Our target measure should have the following properties: intrinsic -it is invariant to isometric transformations of the manifold, e.g. translations or rotations.

multi-scale -it captures both local and global information. We expose our method starting out with heat kernels, which admit a notion of manifold metric and can be used to lower-bound the distance between manifolds.

3.1. HEAT KERNELS ON MANIFOLDS AND GRAPHS:

Based on the heat equation, the heat kernel captures all the information about a manifold's intrinsic geometry (Sun et al., 2009) . Given the Laplace-Beltrami operator (LBO) ∆ X on a manifold X , the heat equation is ∂u ∂t = ∆ X u for u : R + ×X → R + . A smooth function u is a fundamental solution of the heat equation at point x ∈ X if u satisfies both the heat equation and the Dirac condition u(t, x ) → δ(xx) as t → 0 + . We assume the Dirichlet boundary condition u(t, x) = 0 for all t and x ∈ ∂X . The heat kernel k X : X ×X × R + → R + 0 is the unique solution of the heat equation; while heat kernels can be defined on hyperbolic spaces and other exotic geometries, we restrict our exposition to Euclidean spaces X = R d , on which the heat kernel is defined as:

EQUATION

For a compact X including submanifolds of R d , the heat kernel admits the expansion k X (x, x , t) = ∞ i=0 e -λit φ i (x)φ i (x ), where λ i and φ i are the i-th eigenvalue and eigenvector of ∆ X . For t 0 + , according to Varadhan's lemma, the heat kernel approximates geodesic distances. Importantly for our purposes, the Heat kernel is multi-scale: for a local domain D with Dirichlet condition, the localized heat kernel k D (x, x , t) is a good approximation of k X (x, x , t) if either (i) D is arbitrarily small and t is small enough, or (ii) t is for arbitrarily large and D is big enough. Formally, Definition 1 Multi-scale property (Grigor'yan, 2006; Sun et al., 2009 ) (i) For any smooth and relatively compact domain D ⊆ X , lim t→0 k D (x, x , t) = k X (x, x , t) (ii) For any t ∈ R + and any

x, x ∈ D 1 localized heat kernel k D1 (x, x , t) ≤ k D2 (x, x , t) if D 1 ⊆ D 2 . Moreover, if {D n } is an expanding and exhausting sequence ∞ i=1 D i = X and D i-1 ⊆ D i , then lim i→∞ k Di (x, x , t) = k X (x, x , t) for any t.

Heat kernels are also defined for graphs in terms of their Laplacian matrices. An undirected graph is a pair

G = (V, E), where V = (v 1 , . . . , v n ), n = |V |, is the set of vertices and E ⊆ (V × V ) the set of edges. The adjacency matrix of G is a n × n matrix A having A ij = 1 if (i, j) ∈ E and A ij = 0 otherwise. The normalized graph Laplacian is the matrix L = I-D -1 2 AD -1 2

, where D is the diagonal matrix in which entry D ii holds the degree of node i, i.e, D ii = n j=1 A ij . Since the Laplacian matrix is symmetric, its eigenvectors φ 1 , . . . , φ n , are real and orthonormal. Thus, it is factorized as L = ΦΛΦ , where Λ is a diagonal matrix with the sorted eigenvalues λ 1 ≤ . . . ≤ λ n , and Φ is the orthonormal matrix Φ = (φ 1 , . . . , φ n ) having the eigenvectors of L as its columns. The heat kernel on a graph is also given by the solution to the heat equation on a graph, which requires an eigendecomposition of its Laplacian:

H t = e -tL = Φe -tΛ Φ = i e -tλi φ i φ i .

A useful invariant of the heat kernel is the heat kernel trace hkt X : X × R + 0 → R + 0 , defined by a diagonal restriction as hkt X (t) = X k X (x, x, t)dx = ∞ i=0 e -λit or, in the discrete case, hkt L (t) = Tr(H t ) = i e -tλi . Heat kernels traces (HKTs) have been successfully applied to the analysis of 3D shapes (Sun et al., 2009) and graphs (Tsitsulin et al., 2018) . The HKT contains all the information in the graph's spectrum, both local and global, as the eigenvalues λ i can be inferred therefrom (Mémoli, 2011, Remark 4.8) . For example, if there are c connected components in the graph, then lim t→∞ hkt L (t) = c.

3.2. CONVERGENCE TO THE LAPLACE-BELTRAMI OPERATOR:

An important property of graph Laplacians is that it is possible to construct a graph among points sampled from a manifold X such that the spectral properties of its Laplacian resemble those of the Laplace-Beltrami operator on X . Belkin and Niyogi (Belkin & Niyogi, 2002) proposed such a construction, the point cloud Laplacian, which is used for dimensionality reduction in a technique called Laplacian eigenmaps. Convergence to the LBO has been proven for various definitions of the graph Laplacian, including the one we use (Belkin & Niyogi, 2007; Hein et 

3.3. SPECTRAL GROMOV-WASSERSTEIN DISTANCE:

Even while it is a multi-scale metric on manifolds, the heat kernel can be spectrally approximated by finite graphs constructed from points sampled from these manifolds. In order to construct a metric between manifolds, Mémoli (2011) suggests an optimal-transport-theory-based "meta-distance": a spectral definition of the Gromov-Wasserstein distance between Riemannian manifolds based on matching the heat kernels at all scales. The cost of matching a pair of points (x, x ) on manifold M to a pair of points (y, y ) on manifold N at scale t is given by their heat kernels k M , k N :

Γ(x, y, x , y , t) = |k M (x, x , t) -k N (y, y , t)| .

The distance between the manifolds is then defined in terms of the infimal measure coupling

d GW (M, N ) = inf µ sup t>0 e -2(t+t -1 ) Γ L 2 (µ×µ) ,

where the infimum is sought over all measures µ on M × N marginalizing to the standard measures on M and N . For finite spaces, µ is a doubly-stochastic matrix. This distance is lower-bounded (Mémoli, 2011) in terms of the respective heat kernel traces as:

EQUATION

This lower bound is the scaled L ∞ distance between the heat trace signatures hkt M and hkt N .

The scaling factor e -2(t+t -1 ) favors medium-scale differences, meaning that this lower bound is not sensitive to local perturbations. The maximum of the scaling factor occurs at t = 1, and more than 1 -10 -8 of the function mass lies between t = 0.1 and t = 10.

3.4. HEAT TRACE ESTIMATION:

Calculating the heat trace signature efficiently and accurately is a challenge on a large graph as it involves computing a trace of a large matrix exponential, i.e. Tr(e -tL ). A naive approach would be to use an eigendecomposition exp(-tL) = Φ exp(-tΛ)Φ , which is infeasible for large n.

Recent work (Tsitsulin et al., 2018) 

EQUATION

where the function of interest f (•) = exp(•) and v i are n v random vectors drawn from a distribution p(v) with zero mean and unit variance. A typical choice for p(v) is Rademacher or a standard normal distribution. In practice, there is little difference, although in theory Rademacher has less variance, but Gaussian requires less random vectors (Avron & Toledo, 2011) .

To estimate the quadratic form v i f (L)v i in (3) with a symmetric real-valued matrix L and a smooth function f , we plug in the eigendecomposition L = ΦΛΦ , rewrite the outcome as a Riemann-Stieltjes integral and apply the m-point Gauss quadrature rule (Golub & Welsch, 1969):

EQUATION

where µ j = [Φ v i ] j and µ(t) is a piecewise constant function defined as follows Given the symmetric matrix L and an arbitrary starting unit-vector q 0 , the m-step Lanczos algorithm computes an n × m matrix Q = [q 0 , q 1 , . . . , q m-1 ] with orthogonal columns and an m × m tridiagonal symmetric matrix T, such that Q LQ = T. The columns of Q constitute an orthonormal basis for the Krylov subspace K that spans vectors {q 0 , Lq 0 , . . . , L m-1 q 0 }; each q i vector is given as a polynomial in L applied to the initial vector q 0 : q i = p i (L)q 0 . These Lanczos polynomials are orthogonal with respect to the integral measure µ(t). As orthogonal polynomials satisfy the three term recurrence relation, we obtain p k+1 as a combination of p k and p k-1 . The tridiagonal matrix storing the coefficients of such combinations, called the Jacobi matrix J, is exactly the tridiagonal symmetric matrix T. A classic result tells us that the nodes θ k and the weights ω k of the Gauss quadrature are the eigenvalues of T, λ k , and the squared first components of its normalized eigenvectors, τ 2 k , respectively (see Golub & Welsch (1969) ; Wilf (1962); Golub & Meurant (2009)). Thereby, setting q 0 = v i , the estimate for the quadratic form becomes:

µ(t) =      0, if t < a = λ n i j=1 µ 2 j , if λ i ≤ t < λ i-1 n j=1 µ 2 j , if b = λ 1 ≤ t and

EQUATION

Applying ( 5) over n v random vectors in the Hutchinson trace estimator (3) yields the SLQ estimate:

EQUATION

We derive error bounds for the estimator based on the Lanczos approximation of the matrix exponential, and show that even a few Lanczos steps, i.e., m = 10, are sufficient for an accurate approximation of the quadratic form. However, the trace estimation error is theoretically dominated by the error of the Hutchinson estimator, e.g. for Gaussian p(v) the bound on the number of samples to guarantee that the probability of the relative error exceeding is at most δ is 8 (b) distances from the simple English Wikipedia visualized for IMD, FID, and KID. We consider 16 languages: Polish, Russian, Greek, Hungarian, Turkish, Arabic, Hebrew, English, Simple English, Swedish, German, Spanish, Dutch, Portugese, Vietnamese, and Waray-Waray.

3.5. PUTTING IMD TOGETHER:

We employ the heretofore described advances in differential geometry and numerical linear algebra to create IMD (Multi-Scale Intrinsic Distance), a fast, intrinsic method to lower-bound the spectral Gromov-Wasserstein distance between manifolds.

Algorithm 1 IMD algorithm. function IMDESC(X) G ← kNN(X) L ← Laplacian(G) return Γ = slq(L, s, n v ) function IMDIST(X, Y ) hkt X ← IMDist(X) hkt Y ← IMDist(Y ) return sup e -2(t+t -1 ) |hkt X -hkt Y |

We describe the overall computation of IMD in Algorithm 1. Given data samples in R d , we build a kNN graph G by OR-construction such that its Laplacian spectrum approximates the one of the Laplace-Beltrami operator of the underlying manifold (Ting et al., 2010) , and then compute hkt G (t) = i e -λit ≈ Γ. We compare heat traces in the spirit of Equation ( 2 

4. EXPERIMENTS:

We evaluate IMD on the ability to compare intermediate representations of machine learning models. For instance, in a recommender system we could detect whether a problem is related to the representation or the the classifier in the end of a pipeline. In this section, we show the effectiveness of our intrinsic measure on multiple tasks and show how our intrinsic distance can provide insights beyond previously proposed extrinsic measures.

Summary of experiments.:

We examine the ability of IMD1 to measure several aspects of difference among data manifolds. We first consider a task from unsupervised machine translation with unaligned word embeddings and show that IMD captures correlations among language kinship (affinity or genealogical relationships). Second, we showcase how IMD handles data coming from data sources of unequal dimensionalities. Third, we study how IMD highlights differences among image data representations across initializations and through training process of neural networks.

4.1. COMPARING UNALIGNED LANGUAGE MANIFOLDS:

The problem of unaligned representations is particularly severe in the domain of natural language processing as the vocabulary is rarely comparable across different languages or even different documents. We employ IMD to measure the relative closeness of pairs of languages based on the word embeddings with different vocabularies. Figure 2 (a) shows a heatmap of pairwise IMD scores. IMD detects similar languages (Slavic, Semitic, Romanic, etc.) despite the lack of ground truth vocabulary alignment. On the other hand, Figure 11 in Appendix C shows that FID and KID, are not able to distinguish the intrinsic language-specific structure in word embeddings. Detailed description and setting of the experiment can be found in Appendix C. To this end, we measure IMD distance between English word embeddings of varying dimensions. Figure 3 shows the heatmap of the scores between sets of word vectors of different dimensionalities. Closer dimensionalities have lower distance scores for both metrics. However, IMD better highlights gradual change of the size of word vectors, e.g. word vectors of size 4 and 8 are clearly closer to each other than embeddings of size 4 and 16 in terms of IMD, which is not true for PIP.

4.3. TRACKING THE EVOLUTION OF IMAGE MANIFOLDS:

Next, we employ IMD to inspect the internal dynamics of neural networks. We investigate the stability of output layer manifolds across random initializations. We train 10 instances of the VGG-16 (Simonyan & Zisserman, 2015) network using different weight initializations on the CIFAR-10 and CIFAR-100 datasets. We compare the average IMD scores across representations in each network layer relative to the last layer. As Figure 4 (left) shows, for both CIFAR-10 and CIFAR-100, the convolutional layers exhibit similar behavior; IMD shows that consequent layers do not monotonically contribute to the separation of image representations, but start to do so after initial feature extraction stage comprised of 4 convolutional blocks. A low variance across the 10 networks trained from different random initializations indicates stability in the network structure. We now move on to apply IMD to evaluation of generative models. First, we evaluate the sensitivity of IMD, FID, and KID to simple image transformations as a proxy to more intricate artifacts of modern generative models. We progressively blur images from the CIFAR-10 training set, and measure the distance to the original data manifold, averaging outcomes over 100 subsamples of 10k images each. To enable comparison across methods, we normalize each distance measure such that the distance between CIFAR-10 and MNIST is 1. Figure 5 reports the results at different levels σ of Gaussian blur. We additionally report the normalized distance to the CIFAR-100 training set (dashed lines

4.4. EVALUATING GENERATIVE MODELS:

). FID and KID quickly drift away from the original distribution and match MNIST, a dataset of a completely different nature. Contrariwise, IMD is more robust to noise and follows the datasets structure, as the relationships between objects remain mostly unaffected on low blur levels. Moreover, with both FID and KID, low noise (σ = 1) applied to CIFAR-10 suffices to exceed the distance of CIFAR-100, which is similar to CIFAR-10. IMD is much more robust, exceeding that distance only with σ = 2. Next, we turn our attention to the sample-based evaluation of generative models.

We then train the WGAN (Arjovsky et al., 2017) and WGAN-GP (Gulrajani et al., 2017) models on four datasets: MNIST, FashionMNIST, CIFAR10 and CelebA. We sample 10k samples, Y, from each GAN. We then uniformly subsample 10k images from the corresponding original dataset, X, and compute the IMD, KID and FID scores between X and Y. Table 1 reports the average measure and its 99% confidence interval across 100 runs. IMD, as well as both FID and KID, reflect the fact that WGAN-GP is a more expressive model. We provide details on architecture, training, and generated samples in Appendix C. Additionally, in Appendix C we demonstrate superiority of IMD on synthetic data. 

4.5. INTERPRETING IMD:

To understand how IMD operates, we investigate the behavior of heat kernel traces of different datasets that are normalized by a null model. Tsitsulin et al. (2018) proposed a normalization by the heat kernel trace of an empty graph, which amounts to taking the average, rather than the sum, of the original heat kernel diagonal. However, this normalization is not an appropriate null model as it ignores graph connectivity. We propose a heat kernel normalization by the expected heat kernel of an Erdős-Rényi graph (further details in the Appendix D).

Figure 6 depicts the obtained normalized hkt g for all datasets we work with. We average results over 100 subsamples of 10k images each. For t = 10, i.e., at a medium scale, CelebA is most different from the random graph, while for large-scale t values, which capture global community structure,

dhktg(t) dt

reflects the approximate number of clusters in the data. Surprisingly, CIFAR-100 comes close to CIFAR-10 for large t values; we have found that this is due to the fact that the pre-trained Inception network does not separate the CIFAR-100 data classes well enough. We conclude that the heat kernel trace is interpretable if we normalize it with an appropriate null model.

4.6. VERIFYING STABILITY AND SCALABILITY OF IMD:

In addition to the complexity analysis in Section 3.5, we assess the scaling and sample stability of IMD. Since IMD, like FID, is a lower bound to an optimal transport-based metric, we cannot hope for an unbiased estimator. However, we empirically verify, in Figure 7 (left), that IMD does not diverge too much with increased sample size. Most remarkably, we observe that IMD with approximate kNN (Dong et al., 2011) does not induce additional variance, while it diverges slightly further than the exact version as the number of samples grows.

In terms of scalability, Figure 7 (right) shows that the theoretical complexity is supported in practice. Using approximate kNN, we break the O(n 2 ) performance of KID. FID's time complexity appears constant, as its runtime is dominated by the O(d 3 ) matrix square root operation. Geometry score (GS) fails to perform scalably, as its runtime grows exponentially. Due to this prohibitive computational cost, we eschew other comparison with GS. Furthermore, as IMD distance is computed through a low-dimensional heat trace representation of the manifold, we can store HKT for future comparisons, thereby enhancing performance in the case of many-to-many comparisons.

5. DISCUSSION AND FUTURE WORK:

We introduced IMD, a geometry-grounded, first-of-its-kind intrinsic multi-scale method for comparing unaligned manifolds, which we approximate efficiently with guarantees, utilizing the Stochastic Lanczos Quadrature. We have shown the expressiveness of IMD in quantifying the change of data representations in NLP and image processing, evaluating generative models, and in the study of neural network representations. Since IMD allows comparing diverse manifolds, its applicability is not limited to the tasks we have evaluated, while it paves the way to the development of even more expressive techniques founded on geometric insights.

APPENDIX A TRACE ESTIMATION ERROR BOUNDS:

We will use the error of the Lanczos approximation of the action of the matrix exponential f (L)v = exp -tL v to estimate the error of the trace. We first rewrite quadratic form under summation in the trace approximation to a convenient form,

EQUATION

Because the Krylov subspace K m (L, v) is built on top of vector v with Q as an orthogonal basis of K m (L, v), i.e. q 0 = v and v ⊥ q i for i ∈ (1, . . . , m -1), the following holds

EQUATION

Thus, the error in quadratic form estimate v f (L)v is exactly the error of Lanczos approximation f (L)v ≈ Qf (T)e 1 . To obtain the error bounds, we use the Theorem 2 in Hochbruck & Lubich (1997), which we recite below.

Theorem 2 Let L be a real symmetric positive semi-definite matrix with eigenvalues in the interval [0, 4ρ]. Then the error in the m-step Lanczos approximation of exp -tL v, i.e. m = exp -tL v -Q m exp -tTm e 1 , is bounded in the following ways:

EQUATION

Since v is a unit vector, thanks to Cauchy-Bunyakovsky-Schwarz inequality, we can upper-bound the error of the quadratic form approximation by the error of the exp -tL v approximation, i. We now derive the absolute error between the Hutchinson estimate of Equation ( 3) and the SLQ of Equation ( 6):

Tr nv (f (L)) -Γ = n n v nv i=1 v i f (L)v i - nv i=1 e 1 f (T (i) )e 1 ≤ n n v nv i=1 v i f (L)v i -e 1 f (T (i) )e 1 ≤ n n v nv i=1 m = n m ,

where T (i) is the tridiagonal matrix obtained with Lanczos algorithm with starting vector v i .

EQUATION

Finally, we formulate SLQ as an ( , δ) estimator,

1 -δ ≤ Pr Tr(f (L)) -Tr nv (f (L)) ≤ 2 Tr(f (L)) ≤ Pr Tr(f (L)) -Tr nv (f (L)) + Tr nv (f (L)) -Γ ≤ 2 Tr(f (L)) + 2 Tr(f (L)) ≤ Pr Tr(f (L)) -Γ ≤ Tr(f (L)) ,

For the normalized Laplacian L, the minimum eigenvalue is 0 and f min (0) = exp(0) = 1, hence m ≤ 2 , and the eigenvalue interval has ρ = 0.5. We can thus derive the appropriate number of Lanczos steps m to achieve error , We can see that for most of the temperatures t, very few Lanczos steps m are sufficient, i.e. we can set m = 10. However, the error from the Hutchinson estimator dominates the overall error. Figure 9 shows the error of trace estimation does not change with m and for t = 0.1 is around 10 -3 . In case of a Rademacher p(v), the bound on the number of random samples is n v ≥ 6

≤    20e -m 2 /(2.5t) , √ 2t ≤ m ≤ t (11a) 40t -1 e -0.5t 0.5et m m , m ≥ t (11b)

2 log(2/δ) (Roosta-Khorasani & Ascher, 2015). Employing 10k vectors results in the error bound of roughly 10 -2 . In practice, we observe the performance much better than given by the bound, see Figure 9 .

One particular benefit of small m value is that we do not have to worry about the orthogonality loss in the Lanczos algorithm which often undermines its convergence. Since we do only a few Lanczos iterations, the rounding errors hardly accumulate causing little burden in terms of orthogonality loss between the basis vectors of the Krylov subspace. 

B VARIANCE REDUCTION:

We reduce variance of the randomized estimator through control variates. The idea is to use Taylor expansion to substitute a part of the trace estimate with its easily computed precise value,

EQUATION

= slq exp(-tL) + slq tLslq

EQUATION

where we use the fact that L F = Tr(L L) and that the trace of normalized Laplacian is equal to n. It does reduce the variance of the trace estimate for smaller temperatures t ≤ 1.

To obtain this advantage over the whole range of t, we utilize the following variance reduction form:

EQUATION

where there exists an alpha that is optimal for every t, namely setting α = 1/ exp(t). We can see the variance reduction that comes from this procedure in the Figure 12 .

C EXPERIMENTS DISCUSSION:

Here we include additional results that did not find their way to the main paper body.

C.1 FID AND KID FAIL TO FIND STRUCTURE IN UNALIGNED CORPORA:

Figure 11 shows the matrix of distances for FID and KID aligned and colored in the same way as Figure 2 (a). FID and KID can not find meaningful structure in the data in the same way as IMD as they rely on extrinsic data properties.

C. We use gensim ( Řehůřek & Sojka, 2010) to learn word vectors on the latest Wikipedia corpus snapshot on 16 languages: Polish, Russian, Greek, Hungarian, Turkish, Arabic, Hebrew, English, Simple English, Swedish, German, Spanish, Dutch, Portugese, Vietnamese, and Waray-Waray. We then compute FID, KID and IMD scores on all the pairs, we average 100 runs for the heatmap figures 2. For the different dimensionality experiment, we learn vectors on the English Wikipedia of sizes equal to the powers of 2 from 4 to 512. After that we compute IMD and covariance error, i.e. normalized PIP loss, between the pairs of sizes to generate the heatmap figure 3 .

C.3 VANILLA GAN ON TORUS:

We provide an additional experiment clearly showing the case where IMD is superior to its main competitors, FID and KID.

We train two vanilla GANs on the points of a 3D torus. The bad GAN fails to learn the topology of the dataset it tries to mimic, yet previous metrics cannot detect this fact. IMD, on the contrary, can tell the difference. Figure 10 shows the points sampled from the GAN with some of the points inside the hole. KID and FID confidence intervals overlap for good and bad GANs, meanwhile IMD scores are clearly distinct from each other.

C.4 NORMALIZATION DETAILS:

For the purpose of normalizing IMD, we need to approximate that graph's eigenvalues. Coja-Oghlan (2007) proved that λ 1 ≤ 1 -c d-1 /2 ≤ λ 2 ≤ λ n ≤ 1 + c d-1 /2 for the core of the graph for some We train all our models on a single server with NVIDIA V100 GPU with 16Gb memory and 2 × 20 core Intel E5-2698 v4 CPU. For the experiment summarized in Table 1 in the Section 4.1 we train WGAN and WGAN-GP models on 4 datasets: MNIST, FashionMNIST, CIFAR10 and CelebA and sample 10k samples, Y, from each of the GANs. We uniformly subsample 10k images from the original datasets, X, and compute the IMD, KID and FID scores between X and Y. We report the mean as well as the 99% confidence interval across 100 runs.

Below we report the architectures, hyperparameters and generated samples of the models used for the experiments. We train each of the GANs for 200 epochs on MNIST, FMNIST and CIFAR-10, and for 50 epochs on CelebA dataset. For WGAN we use RMSprop optimizer with learning rate of 5×10 -5 . For WGAN-GP we use Adam optimizer with learning rate of 10 -4 , β 1 = 0.9, β 2 = 0.999. 

E GRAPH EXAMPLE:



:

Our code is available open-source: https://github.com/xgfs/imd.

Figures and Tables:

Figure 1: Two distributions having the same first 3 moments, meaning FID and KID scores are close to 0.

θ k are the quadrature's nodes and ω k are the corresponding weights. We obtain ω k and θ k with the m-step Lanczos algorithm (Golub & Meurant, 2009), which we describe succinctly.

Figure2: (a) IMD distances between language pairs for unaligned Wikipedia word embeddings and (b) distances from the simple English Wikipedia visualized for IMD, FID, and KID. We consider 16 languages: Polish, Russian, Greek, Hungarian, Turkish, Arabic, Hebrew, English, Simple English, Swedish, German, Spanish, Dutch, Portugese, Vietnamese, and Waray-Waray.

), i.e., |hkt G1 (t) -hkt G2 (t)| for t ∈ (0.1, 10) sampled from a logarithmically spaced grid.Constructing exact kNN graphs is an O(dn 2 ) operation; however, approximation algorithms take near-linear time O(dn 1+ω )(Dong et al., 2011;Aumüller et al., 2019). In practice, with approximate kNN graph construction(Dong et al., 2011), computational time is low while result variance is similar to the exact case. The m-step Lanczos algorithm on a sparse n × n kNN Laplacian L with one starting vector has O(knm) complexity, where kn is the number of nonzero elements in L. The symmetric tridiagonal matrix eigendecomposition incurs an additional O(m log m)(Coakley & Rokhlin, 2013). We apply this algorithm over n v starting vectors, yielding a complexity of O(n v (m log m + kmn)), with constant k = 5 and m = 10 by default. In effect, IMD's time complexity stands between those of two common GAN evaluation methods: KID, which is O(dn 2 ) and FID, which is O(d 3 + dn). The time complexity of Geometry Score is unspecified in Khrulkov & Oseledets (2018), yet in Section 4.6 we show that its runtime grows exponentially in sample size.

Figure 3: Comparison of IMD and PIP loss on word embeddings of different dimension. IMD detects subtle changes in the dimensionality. Comparing data having different dimensionality is cumbersome, even when representations are aligned. We juxtapose IMD by PIP loss (Yin & Shen, 2018) which allows the comparison of aligned representations for word embeddings.To this end, we measure IMD distance between English word embeddings of varying dimensions. Figure3shows the heatmap of the scores between sets of word vectors of different dimensionalities. Closer dimensionalities have lower distance scores for both metrics. However, IMD better highlights gradual change of the size of word vectors, e.g. word vectors of size 4 and 8 are clearly closer to each other than embeddings of size 4 and 16 in terms of IMD, which is not true for PIP.

Figure 4: (left) IMD score across convolutional layers of the VGG-16 network on CIFAR-10 and CIFAR-100 datasets; (right) training progression in terms of accuracy (dotted) and IMD (solid) on CIFAR-10 and CIFAR-100 datasets for VGG-16 and ResNet-20, with respect to VGG-16.

Figure 5: FID, KID and IMD on the CIFAR-10 dataset with Gaussian blur.

Figure 6: Plotting the normalized heat trace allows interpretation of medium-and global-scale structure of datasets. Best viewed in color.

Figure 7: Stability and scalability experiment: (left) stability of FID, KID and IMD wrt. sample size on CIFAR-10 and CIFAR-100 dataset; (right) scalability of FID, KID and IMD wrt. sample size on synthetic datasets.

e. |v f (L)ve 1 Uf (Λ)U e 1 | ≤ exp -tL v -Q m exp -tTm e 1 = m . Following the argumentation in Ubaru et al. (2017), we obtain a condition on the number of Lanczos steps m by setting m ≤ 2 f min (λ), where f min (λ) is the minimum value of f on [λ min , λ max ].

Figure 8:Errors (solid) and error bounds (dotted) for the approximation of matrix exponential action with varying temperature t.

Figure8shows the tightness of the bound for the approximation of the matrix exponential action on vector v, m = exp(-tL) -Q m exp(-tT m )e 1 . We can see that for most of the temperatures t, very few Lanczos steps m are sufficient, i.e. we can set m = 10. However, the error from the Hutchinson estimator dominates the overall error. Figure9shows the error of trace estimation does not change with m and for t = 0.1 is around 10 -3 . In case of a Rademacher p(v), the bound on the number of random samples is n v ≥ 62 log(2/δ) (Roosta-Khorasani & Ascher, 2015). Employing 10k vectors results in the error bound of roughly 10 -2 . In practice, we observe the performance much better than given by the bound, see Figure9.

Figure 9: Trace estimation errors (solid) and error bounds (dotted) for: (left) the number of Lanczos steps m with fixed number of random vectors n v = 100; (right) the number of random vectors n v in Hutchinson estimator with fixed number of Lanczos steps m = 10. Lines correspond to varying temperatures t.

Figure 10: Bad GAN produces samples inside the torus hole (red). FID and KID cannot detect such behaviour.

Figure11: FID and KID are not able to capture language affinity from unaligned word2vec embeddings.

Figure 12: Variance of the trace estimate.

Figure 13: CIFAR-10 graph colored with true class labels.

Figure 13 provides visual proof that the 5NN graph reflects the underlying manifold structure of the CIFAR-10 dataset. Clusters in the graph exactly correspond to CIFAR-10 classes.

Figure 14: MNIST samples (left: WGAN, right: WGAN-GP)

be the i th eigenvalue and eigenvector, respectively, of the point cloud Laplacian L tn ; let λ i and φ i be the i th eigenvalue and eigenvector of the LBO ∆. Then, there exists t n → 0 such that Still, the point cloud Laplacian involves the creation of an O(n 2 ) matrix; for the sake of scalability, we use the k-nearest-neighbours (kNN) graph by OR-construction (i.e., based on bidirectional kNN relationships among points), whose Laplacian converges to the LBO for data with sufficiently high intrinsic dimension(Ting et al., 2010). As for the choice of k, a random geometric kNN graph is connected when k ≥ log n /log 7 ≈ 0.5139 log n(Balister et al., 2005); k = 5 yields connected graphs for all sample sizes we tested.

Table content: <table><tr><td colspan="2">al., 2007; Coifman &amp; Lafon, 2006; Ting et al., 2010). We recite the convergence results for the point cloud Laplacian from Belkin &amp; Niyogi (2007):</td></tr><tr><td>Theorem 1 Let λ tn n,i and φ tn n,i lim n→∞</td><td>λ tn n,i = λ i</td></tr><tr><td>lim n→∞</td><td>φ tn n,i -φ i 2 = 0</td></tr></table>

, we use the Stochastic Lanczos Quadrature (SLQ)(Ubaru et al., 2017;Golub & Meurant, 2009). This method combines the Hutchinson trace estimator (Hutchinson, 1989; Adams et al., 2018) and the Lanczos algorithm for eigenvalues. We aim to estimate the trace of a matrix function with a Hutchinson estimator:

Table content: <table/>

IMD agrees with KID and FID across varying datasets for GAN evaluation.We now examine the last network layers during training with different initializations. Figure4(right) plots the VGG16 validation errors and IMD scores relative to the final layer representations of two pretrained networks, VGG16 itself with last layer dimension d = 512 and ResNet-20 with d = 64 and ∼50 times less parameters. We observe that even in such unaligned spaces, IMD correctly identifies the convergence point of the networks. Surprisingly, we find that, in terms of IMD, VGG-16 representations progress towards not only the VGG-16 final layer, but the ResNet-20 final layer representation as well; this result suggests that these networks of distinct architectures share similar final structures.

Table content: <table><tr><td/><td colspan="2">MNIST</td><td colspan="2">FashionMNIST</td><td colspan="2">CIFAR10</td><td colspan="2">CelebA</td></tr><tr><td>Metric IMD</td><td colspan="8">WGAN WGAN-GP WGAN WGAN-GP WGAN WGAN-GP WGAN WGAN-GP 57.74 ± 0.47 10.77 ± 0.42 118.14 ± 0.52 13.45 ± 0.54 18.10 ± 0.36 10.84 ± 0.42 10.11 ± 0.33 2.84 ± 0.31</td></tr><tr><td colspan="2">KID ×10 3 47.26 ± 0.07 FID 31.75 ± 0.07</td><td>5.53 ± 0.03 8.95 ± 0.03</td><td>119.93 ± 0.14 152.44 ± 0.12</td><td>25.49 ± 0.07 35.31 ± 0.07</td><td>93.89 ± 0.09 101.43 ± 0.09</td><td>59.59 ± 0.09 80.65 ± 0.09</td><td>217.28 ± 0.14 205.63 ± 0.09</td><td>92.71 ± 0.08 85.55 ± 0.08</td></tr></table>

2 WORD EMBEDDING EXPERIMENT DETAILS.

Table content: <table><tr><td>metric FID KID IMD</td><td>good GAN 0.00529 ± 0.00070 0.00172 ± 0.00073 9.02059 ± 1.5195</td><td>bad GAN 0.00627 ± 0.00076 0.00259 ± 0.00077 14.0732 ± 2.1706</td></tr></table>

Peer review (iclr2020):
1. Define/explain data manifold: in the abstract, ``data moments* suggests the data is treated as random variables; the authors could explain more on how they make the connection between random variables and manifolds 2. introduction: include an example of ``models having different representation space* 3. introduction: elaborate more on ``unaligned data manifolds* 4. introduction: the distinction between ``intrinsic* and ``extrinsic* can be made more explicit; as currently written, the difference between ``extrinsic*/``multi-scale*/``utilizing higher moments* isn*t clear 5. the authors mentioned the hypothesis that high-dimensional data lies on a low-dimensional manifold, but in (3.1) they only considered without justifying this restriction 6. (3.1) `` * change the dot to a variable 7. (3.1) `` and for all * change to ``for all and 8. equation 1, write ``we restrict our exposition to Euclidean spaces * and remove `` * 9. (3.1) the phrase ``a compact X including * is confusing as is not compact under Euclidean topology 10. (3.1) ``for a local domain with Dirichlet condition * change to ``for a local domain with Dirichlet condition* 11. (3.1) ``Definition 1* is more of a proposition/theorem. 12. (3.1) notation: in the body paragraph the authors use but in ``definition 1* they used 13. (3.1) `` disconnected components* to `` connected components* 14. (3.1) the authors could elaborate more about how HKT contains all the information in the graph*s spectrum. Currently, they only stated the limit equal to the number of connected components. 15. Github code is helpful but it would be better if they can include an algorithm section 16. A figure of the architecture could be included under the ``putting IMD together* section