{'fraction_non_alphanumeric': 0.03580346355215465, 'fraction_numerical': 0.010833668948852194, 'mean_word_length': 4.433369803063457, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wassmap: Wasserstein Isometric Mapping for Image Manifold Learning

In this paper, we propose Wasserstein Isometric Mapping (Wassmap), a nonlinear dimensionality reduction technique that provides solutions to some drawbacks in existing global nonlinear dimensionality reduction algorithms in imaging applications. Wassmap represents images via probability measures in Wasserstein space, then uses pairwise Wasserstein distances between the associated measures to produce a low-dimensional, approximately isometric embedding. We show that the algorithm is able to exactly recover parameters of some image manifolds including those generated by translations or dilations of a fixed generating measure. Additionally, we show that a discrete version of the algorithm retrieves parameters from manifolds generated from discrete measures by providing a theoretical bridge to transfer recovery results from functional data to discrete data. Testing of the proposed algorithms on various image data manifolds show that Wassmap yields good embeddings compared with other global and local techniques.

thought of as prima facie infinite-dimensional. That is, one obtains where x is the (discrete) image, H : X → R n is an imaging (or discretization) operator mapping a Banach space X (often L p (R m ) for some p and m, or more commonly L p (Ω) for some compact Ω ⊂ R m ) and η is some noise (often treated as stochastic and is based on the imaging operator as well as other external factors). The noise η can come from multiple sources including background noise, e.g., randomly occurring features such as non-diseased tissue that are not of primary interest for the task [49]; such noise typically has much higher intrinsic dimension than the signal of interest. Image data can also be corrupted by electronic and quantum noise, which is particularly prevalent in scientific and clinical medical imaging where, for instance, radiation dose may prevent the usage of strong light sources [4].
Many dimensionality reduction methods operate by forming an ε-neighborhood of k-nearest neighbor graph over sample points and embedding this graph (or one derived from it) into some much smaller dimensional space (or sometimes an infinite dimensional space). Examples which use variations of this procedure are Isomap [55], Local Linear Embedding [50], Laplacian Eigenmaps [6], UMAP [43], and Diffusion Maps [14]. Each step above has been an avenue of substantial research. While these methods have enjoyed great success in many areas, there are a few drawbacks, especially for imaging applications. First, the graph formation step is typically done in a heuristic fashion and can be problematic in that it is very sensitive to parameter tuning, i.e., choosing ε or k and how to weight the edges.
Second, the most common framework above assumes that {x i } comes from a Riemannian manifold embedded in Euclidean space. Under this assumption, variants of the above procedure are designed so that (hopefully) the graph-theoretic geodesics closely approximate the manifold geodesics between data points. Bernstein et al. [7] prove that if the sampling density of the points {x i } is sufficiently small with respect to the minimum radius of curvature of the ambient manifold and a prescribed tolerance, the graph geodesics of an ε-neighborhood graph can approximate the manifold geodesics between all pairs (x i , x j ) within the prescribed tolerance. Their results show that success typically requires dense sampling of the manifold, which is often unrealistic in image applications due to sparsity of sampling. Additionally, images of the same objects can have different dimensions (n) in the imaging domain under different imaging systems, which may lead to quite different results in the dimensionality reduction procedure. Many algorithms will downsample images to alleviate this issue, but this can lead to information loss and is not necessary in our proposed framework.
Finally, such models tacitly assume that Euclidean distances between data vectors are semantically meaningful. However, this assumption may be invalid in many applications. Indeed, small variation in pixel intensities results in large Euclidean distances, but the images may be semantically the same. For instance, in object recognition, one would expect a model to understand that two images of a car are the same object even if the car is translated in the frame of one of the images. These two images can have large Euclidean distance, even though they are semantically identical.
1.2. Functional image manifolds in Wasserstein space. We propose the following paradigm shift from the previous discussion. In contrast to many imaging techniques which assume that imaged data is on a manifold without reference to the function space underlying them, we assume a functional manifold hypothesis that {x i } ⊂ R n is obtained from imaging a functional manifold M. A natural question arises: what function space naturally represents image data? Our analysis below assumes that images correspond to probability measures with finite p-th moment; i.e., that x = H(µ) as in (1.1) where µ ∈ W p (R m ), the p-Wasserstein space of probability measures with finite p-th moment M p (µ) := R m |x| p dµ(x) < ∞. The p-Wasserstein space is equipped with the Wasserstein metric arising from Optimal Transport Theory [57]. Given two measures µ, ν ∈ W p (R m ), define the set of couplings Γ(µ, ν) := {γ ∈ P(R 2m ) : π 1 γ = µ, π 2 γ = ν} where P(R 2m ) is the set of all probability measures on R 2m , π 1 is the projection onto the first m coordinates, and π 2 onto the last m coordinates. Thus Γ(µ, ν) is the set of all joint probability measures on R 2m whose marginals are µ and ν [51]. Then we may define Our initial assumption will be that measures are absolutely continuous; however, we will provide a bridge to transfer results from these to arbitrary measures in section 4. Treating images as probability measures and considering their p-Wasserstein distances mitigates issues of geodesic blowup inherent in assuming L 2 as the ambient function space (see [17]) as W p is a length space, meaning it contains all geodesics, which are necessarily finite [51,54]. Additionally, the Wasserstein distance contains more semantic meaning in that the Wasserstein distance between images captures how much energy is required to morph one image into another, and the displacement interpolant (the line from one measure to another in W p ) provides a more natural nonlinear path between measures (see [32], for instance). Our initial theoretical and experimental results indicate that Wasserstein distances have significant advantages over other choices in terms of recovering image manifold parametrizations and providing good low-dimensional embeddings of image manifolds (section 5).
Additionally, use of Wasserstein distance and optimal transport theory provides one with a significant and powerful theoretical framework due to the substantial work on optimal transport related to PDEs and other fields (e.g., [48,51,56,57]). In addition to the bulk of theory developed, the use of optimal transport in the past few years has yielded a plethora of advances to the state-of-the-art in many subfields of Machine Learning (ML). For example, use of Wasserstein distances in training of Generative Adversarial Networks (GANs) leads to substantial improvement and stability of such networks [3,35], and in image processing, use of optimal transport ideas has enabled linearization of nonlinear classification problems [31,32].
Note that the Wasserstein manifold assumption is general, and subsumes a setting which is natural in imaging applications. In many cases, we may readily assume that objects being imaged are compactly supported, nonnegative, and integrable (e.g., we may consider a car to be an element of L ≥0 1 (Ω) for a compact set Ω ⊂ R 3 , where the function value at a given point is the density of the car in that location). Thus, one could consider the data {x i } ⊂ R n to be obtained from imaging a functional manifold M ⊂ L ≥0 1 (Ω) for some compact set Ω ⊂ R m . Assuming the images have unit L 1 -norm implies that we may view this manifold as a subset of P(Ω), the set of probability measures supported on Ω via the mapping f i → µ i such that dµ i = f i dx (with dx being the Lebesgue measure on R m ). These measures have finite p-th moment since Therefore, each µ i is an element of the p-Wasserstein space W p (R m ).

Main results.
We briefly summarize our main results here. For some definitions and background on optimal transport and Wasserstein distance, see section 2. We propose a variant of Isomap called Wassmap which uses Wasserstein distances instead of Euclidean distances. For theory in this work, we treat the case when the graph geodesic computation is excluded, which is a Wasserstein distance based Multidimensional Scaling algorithm; however in experiments we illustrate the behavior of our algorithm with and without the graph geodesic step. This algorithm and the assumption that images correspond to elements of W p are used to explore settings in which image manifold parametrizations can be exactly recovered up to rigid transformation by our algorithm. The main results of this paper are the following (here we state things informally, and rigorous theorems follow in later sections).
• Given discrete samples from a smooth submanifold of W p that is isometric to Euclidean space, the Functional Wassmap Algorithm (Algorithm 3.1) yields an isometric embedding and recovers the parameter set governing the manifold up to rigid transformation. • For manifolds generated by translates or dilations of a fixed absolutely continuous measure, Functional Wassmap (Algorithm 3.1) recovers the translation set or a scaled version of the dilation set up to rigid transformation. This result holds for general p for translations, p = 2 for anisotropic dilations, and general p for isotropic dilations. • For submanifolds of W p generated by "nice" parametrized sets of diffeomorphisms acting on a generating measure, the Wasserstein distances between pairs of measures are the same whether the measure is discrete or absolutely continuous. The final result provides a bridge to transfer results for absolutely continuous measures to discrete measures which arise in practice in imaging applications, e.g., as obtained via (1.1). More specifically, it implies that in certain cases, if Functional Wassmap (Algorithm 3.1) recovers an image manifold parametrization for absolutely continuous measures, then it recovers the parametrization for arbitrary measures, and Discrete Wassmap (Algorithm 4.1) recovers the parametrization for discrete measures.
After discussing prior art in the next subsection, the rest of the paper is organized as follows: section 2 describes in brief the background of Wasserstein distances and other details from optimal transport theory, section 3 describes the Functional Wassmap algorithm and contains the main results and proofs related to it, section 4 describes the Discrete Wassmap algorithm and the theorem transferring Wasserstein computations from the continuous to the discrete measure case. Section 5 contains experiments, section 6 discusses computational aspects of the algorithm, and we end with a brief conclusion section.
1.4. Prior art. Most nonlinear dimensionality reduction methods assume data on or near a low-dimensional manifold in Euclidean space and utilize Euclidean distances between points to estimate manifold geodesics. Isomap, described by Tenenbaum et al. [55] is one of the most classical of these algorithms and is the inspiration of this work. Bernstein et al. [7] showed that dense sampling is required to well-approximate geodesics in the Isomap procedure, still under the assumption of Euclidean manifolds. Zha and Zhang [63] proposed continuum Isomap, assuming continuous sampling of the manifold, and utilizing an integral operator formulation of Isomap and Multi-dimensional Scaling (MDS) [42]. Continuum Isomap therein illuminates the theory of classical Isomap maintaining the Euclidean manifold assumption, but is not a practical algorithm. Donoho and Grimes [17] utilized the functional manifold hypothesis that data lives in a submanifold of L 2 (R m ) as a theoretical tool to study the performance of Isomap. Due to the fact that geodesics formed by the metric induced by the L 2 metric can blow up, e.g., in the case of translates of indicator functions, the authors require convolution of the input measures by Gaussians. They also utilize normalization with respect to a reference geodesic, which our framework does not require.
The works of Kolouri et al. [30,31,32] and others [1,28,45] consider absolutely continuous measures, but instead of working with the Wasserstein distances directly, they work with L 2 distances between the optimal transport maps. This approach can speed up computations [28] compared to our approach, but the theory does not transfer to arbitrary measures as is done here. As with the Donoho and Grimes framework, these works also require a reference image to define the transport distance, whereas our method of utilizing Wasserstein distances directly avoids this. Additionally, exact recovery results for image manifold parametrizations appear to be easier to obtain with our framework than these approaches.
Recently, Kileel et al. [29] study the problem of manifold learning with arbitrary norms. Their assumption remains that the data manifold is embedded in Euclidean space, but they construct an analogue of the graph Laplacian which utilizes an arbitrary norm as opposed to the standard Euclidean norm. In similarity with our work, Kileel et al. are motivated by the fact that Euclidean distances may lack semantic meaning or may lead to inflated computational load compared with other norms. Additionally, in their experiments they employ an approximation of the W 1 distance using wavelet expansions for sparse representations of image data. This is computationally faster than W p approximations; however, the results presented herein are primarily aimed at understanding exact recovery of image manifold parametrizations up to rigid motion, and uniqueness of transport maps holds in p ∈ (1, ∞) but not in the case p = 1, so we focus on this range of p; though our method could be applied to the p = 1 case, but we leave this to future work.
Wang et al. [58] utilized opimal transport metrics to compare images of nuclear chromatin. Part of the method they applied is the same as our Algorithm 4.1 without the optional graph geodesic step. They used W 2 distance to quantitatively measure the difference of nuclei, then combined MDS and Fisher Discriminant Analysis to study the distributions of nuclei. They apply the method for classification, which is different from our end goal of retrieving parameters generating the image manifolds.
2. Basics of Wasserstein Distances and Optimal Transport. Given a measure space X, we denote the space of all finite measures on X by M(X). Given a measure space Y , and a continuous map T : X → Y , the pushforward of a measure µ ∈ M(X) via the map T , denoted T # µ ∈ M(Y ), is the measure which satisfies Suppose X and Y are measure spaces, c : X × Y → R + is a cost function, and µ ∈ M(X) and ν ∈ M(Y ); then the Monge Problem is to find the optimal transport map T : X → Y which minimizes The Monge problem does not always admit a solution, even in seemingly innocuous cases such as discrete measures. To get around this difficulty, Kantorovich proposed the following relaxation, which we call the Kantorovich Problem: where Π(µ, ν) is the class of transport plans, or couplings: Here, π x , π y are the projections of X × Y onto X and Y (i.e., the marginals on X and Y ), respectively. Of use to us will also be the Dual Problem to the Kantorovich Problem: Finally, intimately tied to all of these problems is the Wasserstein Distance. Here, we specialize to the concrete case X = Y = R n and utilize the 2 -norm (quadratic) cost function, i.e., c(x, y) := |x − y| 2 (here and throughout, | · | denotes the Euclidean norm on R m where m may be determined from context). For µ, ν ∈ P(R m ) with finite 2-nd moment, the 2-Wasserstein distance is defined by (1.2).
, and at least one of which is absolutely continuous. Then there exists an optimal transport map T from µ to ν and a unique optimal transport plan π ∈ Π(µ, ν). Additionally, For p ∈ (1, ∞), the optimal transport map is unique and can be characterized by the following theorem of Gangbo and McCann [20]. A function ψ : If µ is absolutely continuous with respect to Lebesgue measure then a map T , which solves the Monge problem, pushing µ forward to ν is uniquely determined µ-almost everywhere by the requirement that it be of the form 3. Functional Wassmap: Algorithm and Theory. In this section, we consider the problem of when an image manifold treated as a submanifold of the Wasserstein space is isometric to Euclidean space. We state the Functional Wassmap algorithm below in full generality, but for treatment of theoretical guarantees we will restrict the class of manifolds we consider. First, we consider the case when the geodesics on the manifold M ⊂ W p (R m ) are given by the W p distance between measures. Below, Step 4 is the optional step of forming a graph whose nodes are the measures and there is some rule for determining neighborhoods of nodes and setting edge weights (for instance, ε-neighborhood or k-nearest neighbors). We use APSP to stand for All-pairs shortest path, i.e., computing graph-theoretic distances between all nodes (for example via Dijkstra's algorithm). Curved manifolds in Wasserstein space will require this graph geodesic step, whereas classical results regarding Multidimensional Scaling imply that flat manifolds do not require this step. For our theoretical results, we restrict to cases when the metric space (M, W p ) is isometric up to a constant to a subset of Euclidean space (Θ, | · | R d ), and suppress the graph geodesic step in Algorithm 3.1.
Steps 5-7 of Algorithm 3.1 are the classical Multi-dimensional Scaling Algorithm, or MDS. An important result for MDS is the following.
A distance matrix is Euclidean provided there exists a point set  Proof. The isometry condition implies existence of a global constant c > 0 such that and uniqueness up to rigid transformation is given by Theorem 3.2.
In subsequent results, we will compute the global scaling factor c for some image manifolds, in which case we obtain recovery of {cθ i } up to rigid transformation. [17] developed a theoretical framework for understanding the behavior of Isomap on image manifolds. They studied whether a normalized version of geodesic distance is equivalent to the Euclidean distance, in which case Isomap recovers the underlying parametrization of image manifolds. They show several positive cases including translation, pivoting and morphing boundaries of black objects on white backgrounds. They also show that Isomap may fail when the parameter space is not convex or the image manifold is not flat (for example, the dilation manifold of rectangles or ellipses).

Comparison to other techniques. Donoho and Grimes
In comparison, Wassmap does not require normalization when computing distances between images. For translation manifolds, Wassmap retreives the underlying parameters without requiring the parameter space to be convex. Wassmap also recovers translation and dilation manifolds generated by a base measure which has nonsmooth pdf, like the indicator function of a domain, whereas Isomap fails in this case due to geodesic blowup.

Translation manifolds.
Given a fixed generating measure µ 0 ∈ W p (R m ) and translation set Θ ⊂ R m , define This simple translation manifold satisfies dim(M) = dim(span(Θ)). We show the following: The crux of the proof of this theorem is the following lemma. This lemma is known [48,Remark 2.19], but for completeness we present the full proof.
) and φ is c-concave. By Theorem 2.2, T is the optimal transport map, whence Proof of Theorem 3.4. Combine Lemma 3.5 and Corollary 3.3.

Dilation manifolds.
Here we will consider dilation manifolds. We begin with the general case of anisotropic dilations along coordinate axes, but our results here are only valid for the case p = 2. Given a dilation set Θ ⊂ R m + (θ ∈ Θ that has strictly positive entries ϑ 1 , . . . , ϑ m ), we define the corresponding manifold with a fixed absolutely continuous generator where we use the slight abuse of notation and consider dµ 0 (D θ ·) = f (D θ x)dx and the dilation matrix is defined by Recall that M p (µ) := R m |x| p dµ(x) is the p-th moment of a measure µ ∈ P(R m ), and let P i µ denote the i-th marginal of µ defined by In the results below, the p-th moment of the i-th marginal is thus M 2 (P i µ) := R m |x i | 2 dµ(x). This choice of notation rather than µ i is to avoid confusion, as subscripts 0 and θ will be used frequently in the sequel.
This theorem can be derived from the following lemma.
and similarly for dµ θ . Then Proof. The proof proceeds by using (MP) to find an upper bound for the Wasserstein distance in question, and (DP) to find a lower bound. These being the same, we use Theorem 2.1 to conclude the result. To show the upper bound, we use the fact that (MP) has a solution, and note that the map Hence, T is the pushforward from µ θ to µ θ . By (MP) and Theorem 2.1, The penultimate equality follows from substituting x → D θ x. Now we use (DP) to find a lower bound for the Wasserstein distance by setting These are easily seen to be in L 1 (µ θ ) and L 1 (µ θ ), respectively. Additionally, hence φ and ψ are feasible solutions to (DP). Finally, by (DP) and Theorem 2.1, and the lemma is proved.
Proof of Theorem 3.6. Lemma 3.7 implies that W 2 (µ θ , µ θ ) = |Sθ − Sθ |, so then the Wasserstein distance matrix arising from {µ θi } is a Euclidean distance matrix with point configuration The uniqueness of recovery of this set up to rigid transformation is given by Theorem 3.2.
Note that the matrix S is determined by the generator µ 0 of the manifold. Thus, in order to retrieve the parameters {θ i } N i=1 , information on µ 0 is required. In certain conditions, Algorithm 3.1 recovers {θ i } m i=1 up to a constant.
Remark 3.8. Lemma 3.7 may not generalize to any p > 1 when the dilations are not isotropic (i.e., do not satisfy ϑ i = ϑ for all i). By Theorem 2.2, φ and T are linked by should be the gradient of some function φ. This is not true when the dilations are anisotropic.
Remark 3.10. Note that if the dilations occur only along certain coordinates, i.e., Θ is supported on a d-dimensional coordinate plane for some 1 ≤ d < n, then one can specify the embedding dimension in Corollary 3.9 to be d rather than m. In this case, one recovers the isometric projection of Θ into R d which ignores the undilated coordinates. For example, if Θ only has elements other than 1 in coordinates {i 1 , . . . , i d }, then Functional Wassmap will recover (up to rigid transformation) For isotropic dilations, the result holds for arbitrary p ∈ (1, ∞).
p θ i } is equivalent up to rigid transformation to {Sθ i } where S is as in Theorem 3.6, so the conclusion of Corollary 3.11 is not contradictory.
We end this subsection by giving some concrete examples. The first is for the simple case when the density function of µ 0 is symmetric, giving a concrete example of Corollary 3.9.
The proof of this proposition follows from direct calculation of the moments in Corollary 3.9 and so is omitted. The converse of Proposition 3.12 is not necessarily true. That is, the condition for some c ∈ R does not imply that µ 0 is symmetric across x 1 = x 2 . Indeed, consider the following example: suppose dµ 0 = 1 For further illustration, the following corollary, easily obtained by computing the relevant second moments from Theorem 3.6, shows what one recovers for a dilation manifold when the generating measure is the indicator function of a domain suitably normalized.
Corollary 3.13. Let A be a rectangle in R m with endpoints a i,1 , a i,2 on the i-th coordinate axis, and let dµ 0 = 1 where S A is the diagonal matrix whose diagonal entries are defined as Note that if the parameter set is a lattice in R m , i.e., Θ = α 1 Z × · · · × α m Z, then Functional Wassmap will recover the set α 1 M 3.5. Rotation manifolds. We will show in subsequent experiments that the discrete Wassmap algorithm is capable of recovering the underlying circle governing a rotational manifold. However, at present, the authors do not have a proof analogous to the above results for this case. Let a rotation manifold be defined as follows: Consider the following: Theorem 3.14 ([51, Theorem 1.22]). Suppose µ, ν ∈ W 2 (R m ) and µ gives no mass to (d − 1)-surfaces of class C 2 . Then there exists a unique optimal transport map T from µ to ν, and it is of the form T = ∇u for a convex function u.
A direct consequence of Theorem 3.14 is that a rotation matrix R θ is not the optimal transport map from a measure to its pushforward under R θ as this is not the gradient of a convex function. Consequently, exactly computing W 2 (µ 0 (R θ ·), µ 0 (R θ ·)) is nontrivial.
On the other hand, we can give an upper bound for the Wasserstein distance of a rotated version of a fixed measure with itself. Restrict to clockwise rotation in R 2 by angle θ ∈ (0, 2π), and let R θ be the resulting rotation matrix. One can verify that

Discrete Wassmap: Algorithm and Theory.
In imaging practice, one obtains discrete vectors rather than continuous distributions, so a practical version of Algorithm 3.1 must take this into account. To do this, one must consider how to form a probability measure from a given image. Given a two-dimensional (planar) or multidimensional (e.g., volumetric) image in pixel/voxel representation, that is g = [g 1 , . . . , g D ] ∈ R D where D is the total number of pixels or voxels, we will assign a discrete measure P (g) ∈ W 2 (R m ) by selecting a set of D locations x n ∈ R m , and assigning mass g n > 0 to a corresponding physical location x n and normalizing: where δ xn is a Dirac mass at location x n . The locations x n are most conveniently assumed, at least initially, to lie on a regular grid in the ambient space R m . Given two images (with common ambient dimension m but not necessarily the same D), the problem of computing the Wasserstein distance between µ i = P (g i ) and µ j = P (g j ) reduces to a discrete optimization problem for which many algorithms exist [21,48]. Below we summarize the Discrete Wassmap algorithm which mimics the procedure described above. Note that we state the algorithm for image input, but one could equally well state it simply for discrete probability measures input in which case one simply skips the measure construction step.  For instance, if one translates a discrete measure, is the Wasserstein distance the same as in the absolutely continuous case (the magnitude of the translation)? Here we show that this is the case for a wide variety of discrete measures and transformations of them. We will state our results in terms of the pushforward operators defining the transformation of a base measure.
The following theorem provides a bridge which allows one to transfer results on recovery of Wasserstein image manifold parametrizations from manifolds generated by absolutely continuous measures to those generated by arbitrary measures. Note that there is no requirement on the generating measure µ 0 aside from the fact that it lies in W p ; it may have a mix of continuous and discrete spectra, and need not have compact support. Theorem 4.1 shows that if Wasserstein distances between absolutely continuous measures generated by a given parameter set depend only on the parameters, then the Wasserstein distances between arbitrary measures likewise depend only on the parameters. Thus, if Functional Wassmap recovers a parameter set for absolutely continuous generating measure µ 0 , then Discrete Wassmap recovers the manifold generated analogously from a discrete measure µ 0 . Additionally, since Theorem 4.1 holds for arbitrary measures, we find that if Functional Wassmap recovers a parameter set for absolutely continuous generating measure, then in fact it also recovers the parameter set for arbitrary generating measure. Below, we let be the multivariate Gaussian kernel on R m , and * represents convolution of measures.  Let p ∈ [1, ∞). Suppose that for all absolutely continuous µ 0 ∈ W p (R m ), M(µ 0 , Θ) = {T θ µ 0 : θ ∈ Θ} is a smooth submanifold of W p (R m ), that T θ is Lipschitz for all θ, and that for all θ, θ ∈ Θ and all absolutely continuous µ 0 , W p (T θ µ 0 , T θ µ 0 ) = f (θ, θ , µ 0 ) for some function f dependent only upon θ, θ , and µ 0 , for which lim σ→0 f (θ, θ , µ 0 * g σ ) = f (θ, θ , µ 0 ). Then for any ν 0 ∈ W p (R m ), W p (T θ ν 0 , T θ ν 0 ) = f (θ, θ , ν 0 ).
The crux of the proof of this theorem is the lemma below. We take ν σ ν to mean weak convergence of measures, which by the Portmanteau Theorem is equivalent to the statement ν σ (A) → ν(A) for all continuity sets A of ν (i.e., ν(∂A) = 0).
Proof. By [51, Lemma 5.2], W p (T θ µ, T θ µ) = lim σ→0 W p (T θ µ * g σ , T θ µ * g σ ). While it is not true in general that W p (T θ µ * g σ , T θ µ * g σ ) = W p (T θ (µ * g σ ), T θ (µ * g σ )), we will show that the limits of these two expressions is the same. Considering their difference and utilizing the reverse triangle inequality then triangle inequality below, we have We claim that lim σ→0 W p (T θ µ, T θ (µ * g σ )) = lim σ→0 W p (T θ µ, T θ µ * g σ ) = 0 for any θ. By [51,Lemma 5.11], it is sufficient to prove the following: Here and subsequently all integrals are over R m . For 3) and 4), we will first prove the statements for integer p and show that it can be generalized to any p ∈ [1, ∞).
Proof of 1) Suppose A is a measurable set of µ, then by definition T −1 (A) is a continuity set of T θ µ. Then, where convergence follows from the fact that µ * g σ µ for any µ. Item 2 is well-known and follows from a simple computation, so is omitted. Proof of 3) First, we note that by direct computation, |T (x + y)| p g σ (y)dydµ(x) and where the second equality follows from g σ being a probability density function.
With these observations in hand, if p is an integer, we have that the difference | |x| p d(T θ (µ * g σ ) − T θ µ)| is bounded as follows: The second inequality follows from utilizing the fact that T is Lipschitz. The final integral can be bounded by a sum of integrals involving products of |x| p−1 and |y| α for various integer α ∈ [1, p]. That this quantity goes to 0 as σ → 0 follows from the facts that µ has finite p − 1 moment and that any Gaussian moment tends to 0. Indeed, by substitution, The integral above is a constant depending only upon p and m, so the conclusion follows by application of the Dominated Convergence Theorem. Proof of 4) By similar argument, but noting that we see that for integer p, The finite moment of g σ tends to 0 as before. For non-integer p, the conclusion follows from the fact that |a p − b p | ≤ |a p − b p | + |a p − b p | and µ ∈ W p has finite p -th moment.
With this Lemma we are now in a position to finish the proof of the main theorem of this section.
Proof of Theorem 4.1. Let g σ be the multivariate Gaussian with variance σ as before. Then by Lemma 4.2, we have By combining the above results, we readily see that Discrete Wassmap recovers the translation set for translation of discrete measures, and the scaled dilation set for dilation image manifolds.

Experiments.
To demonstrate our theoretical results, we provide several experiments 1 using both synthetically generated two-dimensional image data and the standard MNIST digits dataset [33]. For each synthetic experiment, a fixed absolutely continuous base measure µ 0 ∈ W 2 (R 2 ) with density f 0 (x) is selected, then a manifold M(µ 0 , Θ) is sampled by applying the parametric transformation T θ to µ 0 for a finite number of θ values {θ 1 , . . . , θ N } ⊂ Θ, resulting in the measures µ θi and corresponding densities f θi . These (continuum) images are subsequently discretized by performing a spatial sampling, selecting {x 1 , . . . , x D } ⊂ R 2 , evaluating each density f θi (x) at these points, then forming the discrete measure (4.1).
Comparisons are shown to Euclidean MDS, Isomap, and Diffusion Maps. MDS and Isomap are the nearest, most faithful comparison to our method as they are also global algorithms, but the other methods mentioned are local methods, which we employ for a more comprehensive comparison. For Isomap, the kNN graph is used for geodesic estimation, with k tuned to give the best visual result. For Diffusion Map embeddings, we employ the standard Gaussian kernel with the automatic epsilon selection algorithm described in [8]. Our experiments use the Wassmap variant without the graph geodesic computation step with the exception of the MNIST experiment in subsection 5.4 ( Figure 5). For further illustration of Wassmap with and without use of geodesics, see the supplemental material. We use Euclidean distances in all methods except for Wassmap. Note that methods other than Wassmap assume a 'pixel' representation of images; that is, each image is treated as an element of R D for some fixed D. One can obtain such a representation by following the steps outlined above but keeping points of zero density (for more discussion, see the supplemental material and code). In all figures, points in the original parameter set are color coded, and the corresponding points in the embedding depend on the initial point color in a one-to-one correspondence. Thus one can see where initial parameter points end up in each embedding.

Translation manifold.
In this set of experiments, we take the base measure µ 0 to be the indicator function of a disc of radius 1, that is dµ 0 = 1 Both Figures 1 and 2 show that Wassmap recovers the underlying translation grid up to rigid motion as predicted by Theorem 3.4; in Figure 1 a rotation appears, but the side-lengths of the embedded grid are 2 as in the original parameter set Θ 1 .
In both experiments, the translated discs overlap; consequently, the other methods produces embeddings that appear coherent despite failing to recover the parameter set exactly. In particular, the other embeddings of Θ 1 appear to be morphed grids, and the scale is dilated in a way that the Wassmap embedding is not. Two advantages of Wassmap in this case are that it is not subject to careful parameter tuning, and the size of the disk and its relation to the parameter grid does not matter, whereas for other methods parameter tuning may play an important role, and the translates of the disk must overlap significantly for the pairwise Euclidean distances to be meaningful. Figure 2 illustrates that Wassmap is capable of recovering nonconvex translation parameter sets in contrast to both discrete and continuum Isomap [17]. The Isomap embedding is largely incoherent in this case, while the Euclidean MDS and Diffusion Map embeddings exhibit significant skewing as well as overlapping points in the embedding.

Dilation manifold.
To illustrate the case of dilations, we consider the same base measure as in the translation case (disc support function centered at (0, 0)), but  now apply the anisotropic dilation transformation D θ as discussed in subsection 3.4, where the parameters ϑ 1 , ϑ 2 come from a regular 4x4 subgrid of [0.5, 2] × [0. 5,4]. The dilation parameter grid and result of the different embeddings are shown in Figure 3. The Euclidean MDS embedding is the next best compared to Wassmap, the latter of which recovers the structure of the parameter set faithfully, but all other embeddings are relatively poor. Note that the dilation grid has size 1.5 × 3.5, and the Wassmap embedding has size approximately 1.75 × .75. One can compute the projected second moment of the base measure dµ 0 = 1 π 1 D dx as (M 2 (P i µ 0 )) 1 2 = 1 2 . Corollary 3.9 states that Wassmap recovers the dilation grid up to this factor and a global rotation. Thus we see that the Wassmap embedding does recover the original dilation grid multiplied by 0.5 (the moment term) rotated by π/2.

Rotation manifold.
To illustrate the case of rotational manifolds, we consider the base measure µ 0 as the indicator of an ellipse with major radius 1 and minor radius 0.5, centered at (x, y) = (0, 1) This measure is rotated about the origin to obtain the sampled manifold at uniform angles θ i ∈ [0, 2π]; embedding results are shown in Figure 4. We see that all methods approximately or exactly recover a 0 . 5 0 . 8 7 5 1 . circular manifold. This experiment provides evidence that Wassmap is capable of recovering rotational manifolds, though at present we are not able to prove this as discussed in subsection 3.5.

Embedding MNIST.
Here we show the effect of Wassmap on embedding MNIST handwritten digits [33]. We randomly sample 250 handwritten 1s, 2s, 7s and 9s from MNIST and compute the 2-dimensional Wassmap embeddings corresponding to two different choices of k when forming the kNN graph. We also computed the Isomap embeddings based on a kNN graph with the same k; recall that for both Wassmap and Isomap, the large k limit corresponds to skipping the Geodesic computation step and using the raw pairwise distances. Figure 5 shows the resulting embeddings, with Euclidean MDS and Diffusion Maps also shown for comparison. Note that Wassmap produces an embedding for which, for any value of k, the classes are relatively easily separated by a kernel SVM or nearest neighbor classifier, whereas the Isomap embedding appears sensitive to the choice of k and results in nontrivial class overlaps for any k. Both the Euclidean MDS and Diffusion Map embeddings also struggle to separate classes, particularly 9s and 7s. A full analysis of classification performance and its dependence on the embedding dimension d is a subject of future work.  6. Computational aspects. From a practical perspective, the biggest drawback of the proposed approach is that of computational cost. Here, we discuss several relaxations and approximations that can be done that speed up Wassmap.
6.1. Approximations of Wasserstein distances. The experiments done here used the exact linear program solver to compute Wasserstein distances, but one can trade off accuracy of the embedding with speed of computation by utilizing approximation algorithms that approximate each of the Wasserstein distances. In general, an exact Wasserstein distance computation between discrete measures with n points of mass carries complexity Ω(n 3 log n) without enforcing additional constraints on the measures [47]. However, one can utilize any Wasserstein distance approximation in the Wassmap framework; some examples are entropic regularization and Sinkhorn distances [2,9,15,36,53], multiscale methods [21,23,52], and linearization techniques [24]. Sinkhorn distances can be computed in O(n 2 ) time, but many approximate Wasserstein distance methods do not have concrete computational complexities, which makes it difficult to write down an overall complexity for Wassmap.  Figure 3. The Wasserstein distance matrix is approximated by the Sinkhorn distance with different regularization parameters λ. The right embedding took about a third as long as the standard Wassmap embedding.
6.2. Approximating the MDS kernel matrix. One note regarding Algorithm 3.1 is that the Wasserstein distance computations can be done in parallel, which can greatly speedup the computation time, although computing O(N 2 ) Wasserstein distances is still prohibitive for large N .
One could use the Linear Optimal Transport (LOT) framework developed by Wang et al. [59] and furthered by several others [1,28,32,45,46] as an alternative. Cloninger and Moosmüller show that for some types of object manifolds, the LOT distance is equal to the Wasserstein distance between manipulated images, hence our theoretical recovery results could be obtained for LOT distances in these cases as well.
Here, the LOT distance is defined as the L 2 -norm between the transport (Monge) maps from each image to a fixed reference measure (e.g., a Gaussian). In this way, one could compute O(N ) LOT distances rather than O(N 2 ) Wasserstein distances prior to MDS. One potential drawback of this approach is that for more curved manifolds in Wasserstein space, or for diffeomorphisms that do not satisfy the compatibility condition of [45], the LOT approximation is not exact, and the embedding may suffer from this.
An alternative would be to use the Nyström method [22,60] to approximate the squared distance matrix (B) of MDS directly. A Nyström approximation of a symmetric matrix is of the form B ≈ CW † r C T where C = W (:, I) is a column submatrix of B, W = B(I, I) is the intersection of C and C T , and W † r is the Moore-Penrose pseudoinverse of W r , which is the truncated SVD of W of rank r. It is known that in cases where a kernel matrix is incoherent and approximately low-rank, a rank r Nyström approximation of an N × N matrix can yield a good approximation by only computing O(r log N ) columns of the kernel matrix. For the Wasserstein distance matrix, this would allow one to perform only O(rN log N ) Wasserstein distance computations, which is slightly more than using LOT but considerably less than computing the full distance matrix. As an illustration, we rerun the dilation experiment above but with a 12 × 12 grid in the parameter space (Figure 7), and we sample 4 out of 144 (log N ) columns to form the Nyström approximation to the Wasserstein distance matrix, which we then embed via MDS. For each Wasserstein computation we use the Sinkhorn algorithm of Cuturi. One can see that the shape of the embedding remains, though there are some errors.
Note that this Nyström approximation is of the squared Wasserstein distance matrix in Algorithm 4.1, which is different than the method of Altschuler et al. [2], which uses the Nyström method to approximate the cost matrix for a single Wasserstein distance computation between a pair of measures. These methods can be combined to yield a faster, more scalable approximation of our Wassmap algorithm.

Updating an MDS embedding.
In practice, if one already has a lowdimensional embedding via MDS, one might want to update the embedding when a new image is obtained. This can easily be done in the following way: first compute the Wasserstein distance from the new image to the other images (this requires N W 2 computations), and perform a rank-one update of the SVD of the matrix B in Algorithm 3.1 (or Algorithm 4.1), which can be done in O(N d + d 3 ) flops [11], which is much faster than recomputing a new SVD of B at a cost of O(N 2 d) (here d is the embedding dimension, which is the rank of the embedding matrix in MDS).
6.4. Choosing the embedding dimension. In any method requiring choice of embedding dimension, there are several methods one can employ. For example, a scree plot of the singular values of the distance matrix can give an idea by looking for an elbow, or bend in the graph. More sophisticated techniques involve estimating the dimension of the manifold by local PCA or something similar [37,38,39].

Conclusion and Future
Outlook. This paper proposed the use of Wasserstein distances in the Isomap algorithm (and its precursor MDS) as a more suitable measure of distance between images. The resulting Wassmap algorithm and its variants were shown to recover (up to rigid transformation) several parametrizations of image manifolds, including translation and dilation sets. We provided a bridge which transfers functional manifold recovery results to discrete recovery, which illustrate that the Discrete Wassmap algorithm recovers parametrizations of image manifolds generated by discrete measures. The practical experiments illustrate the effectiveness of the proposed framework on various synthetic and benchmark data. There is more to be explored regarding Wassmap, including its potential to recover rotation manifolds, those generated by composition of different operations (e.g., translation plus dilation or rotation), and manifolds generated by some class of parametrized diffeomorphisms acting on one or multiple generators. It also remains to explore the effects of additive noise (η in (1.1)) and the structure of the imaging operator H.
Future work will also explore the use of Wasserstein distances in other manifold learning paradigms, including local methods such as LLE and tSNE, as well as use of W p for other p ∈ [1, ∞) (for example, Kileel et al. [29] approximate the classic Earthmover's Distance W 1 in the Laplacian eigenmap setting). One could easily replace Euclidean distances with Wasserstein distances in a naïve way in any manifold learning algorithm such as Laplacian eigenmaps or Diffusion Maps. We have done this in a simple example in the Supplemental Material here, but the results are not competitive with other methods for this particular task of parametrization recovery. Better understanding how to use these algorithms for data coming from submanifolds of Wasserstein space would be an interesting avenue of future study.
Additional studies will be done on choosing ε adaptively for the neighborhood graph step [44], and combination of Wasserstein distance based algorithms with task performance such as classification or clustering (see [40] for recent work in this direction). This experiment shows what is already well-understood regarding Isomap and MDS: flat manifolds can be recovered via MDS embeddings, while curved manifolds require the graph geodesic step in Isomap.
Appendix C. Second rotation experiment. Here we show a second rotation experiment where we take µ 0 to be the indicator function of an ellipse with major radius 1 and minor radius 0.5, but whose initial center is at the origin (x, y) = (0, 0). We rotate the ellipse by N = 21 uniformly sampled angles in [0, 2π), and the results of the embeddings are shown in Figure 9. Due to rotational symmetry of the base figure, we would not expect any method to recover the 'correct' manifold in this case. Appendix D. On low-rankness of Wasserstein distance matrices. In section 6.2, we discuss the use of the Nyström method for approximating the Wasserstein distance matrix W of the Wassmap algorithm. Use of the Nyström method requires that the kernel matrix in question be approximately low-rank. Here we show scree plots the spectrum of W from various experiments in the main paper.  One can see that the translation distance matrix is well-approximated by a lowrank matrix. The MNIST distance matrix can be approximated by an approximately rank 25 matrix. The dilation matrix has relatively flat spectrum after the initial drop; because of the scale of the singular values (O(1)), this matrix is somewhat harder to approximate, but one can see that one gains no approximation power from taking more than about 10 singular values until one takes almost all of them.
Note that these scree plots can also be used to estimate the dimensionality of the embedding as discussed in section 6.4. tioned above, a naïve extension of our idea to algorithms like Laplacian eigenmaps or Diffusion Maps may not be the correct generalization of those methods. For illustration, we repeat the dilation experiment here and in the bottom center use a simple variant of the Diffusion map in which W 2 distances are used instead of Euclidean. One can see in Figure 11 that the embedding does not do any better than the diffusion map with Euclidean distances, and we attribute this to the fact that a more sophisticated understanding of how to approximate submanifolds of Wasserstein space via spectral graph theory may need to be developed. Appendix F. Using W 1 . We have mainly illustrated the Wassmap algorithm for quadratic Wasserstein distances. However, since W 1 distances are also commonly used, it makes sense to consider these. Here, we repeat the dilation example for Wassmap using W 1 distances and we see in Figure 12 that the embedding fails to recover the grid exactly up to scaling; however, the embedding does still maintain the general characteristics of the dilation set. It is perhaps unsurprising that the qualitative behavior of the W 1 distance based embedding is similar to that using W 2 distances. Indeed, all W p distances are based on smoothly distorting one image to the other. This experiment incidates that further study of the case p = 1 is warranted.   on (a, b). A subset of 64 such zeros in shown in Figure 13. Figure 14 shows the two-dimensional Wassmap and Isomap embeddings of {f θi } as defined above, along with the scree plots for each. The Wassmap embedding shows a more distinctly clustered embedding. Isomap demonstrates similar structures but does not seem able to separate subtle morphological variations as easily. The scree plots demonstrate the improved embedding efficiency of Wassmap versus Isomap. This evidence suggests that Wassmap may be capable of finding the structure of manifolds generated by a restricted family of diffeomorphisms, but future exploration is needed in this case.  Comparison of the two-dimensional embeddings and scree plots for Wassmap and Isomap for the grid deformation family. The Wassmap technique produces a much smoother embedding with clear geometric structure, while the Isomap method produces a less coherent embedding with diminshed cluster separation. The scree plots show that Wassmap is more dimensionally efficient.