{'fraction_non_alphanumeric': 0.04079813767874958, 'fraction_numerical': 0.013382108413701363, 'mean_word_length': 4.7333740085417935, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transport-based analysis, modeling, and learning from signal and data distributions

Transport-based techniques for signal and data analysis have received increased attention recently. Given their abilities to provide accurate generative models for signal intensities and other data distributions, they have been used in a variety of applications including content-based retrieval, cancer detection, image super-resolution, and statistical machine learning, to name a few, and shown to produce state of the art in several applications. Moreover, the geometric characteristics of transport-related metrics have inspired new kinds of algorithms for interpreting the meaning of data distributions. Here we provide an overview of the mathematical underpinnings of mass transport-related methods, including numerical implementation, as well as a review, with demonstrations, of several applications.

Software accompanying this tutorial is available at [134].

B. Why transport?
In recent years numerous techniques for signal and image analysis have been developed to address important learning and estimation problems. Researchers working to find solutions to these problems have found it necessary to develop techniques to compare signal intensities across different signal/image coordinates. A common problem in medical imaging, for example, is the analysis of magnetic resonance images with the goal of learning brain morphology differences between healthy and diseased populations. Decades of research in this area have culminated with techniques such as voxel and deformation-based morphology [6], [7] which make use of nonlinear registration methods to understand differences in tissue density and locations [66], [149]. Likewise, the development of dynamic time warping techniques was necessary to enable the comparison of time series data more meaningfully, without confounds from commonly encountered variations in time [80]. Finally, researchers desiring to create realistic models of facial appearance have long understood that appearance models for eyes, lips, nose, etc. are significantly different and must thus be dependent on position relative to a fixed anatomy [30]. The pervasive success of these, as well as other techniques such as optical flow [10], level-set methods [33], deep neural networks [142], for example, have thus taught us that 1) nonlinearity and 2) modeling the location of pixel intensities are essential concepts to keep in mind when solving modern regression problems related to estimation and classification.
We note that the methodology developed above for modeling appearance and learning morphology, time series analysis and predictive modeling, deep neural networks for classification of sensor data, etc., is algorithmic in nature. The transportrelated techniques reviewed below are nonlinear methods that, unlike linear methods such as Fourier, wavelets, and dictionary models, for example, explicitly model jointly signal intensities as well as their locations. Furthermore, they are often based on the theory of optimal mass transport from which fundamental principles can be put to use. Thus they hold the promise to ultimately play a significant role in the development of a theoretical foundation for certain subclasses of modern learning and estimation problems.

C. Overview and outline
As detailed below in section II the optimal mass transport problem first arose due to Monge [109]. It was later expanded by Kantarovich [74], [75], [77], [76] and found applications in operations research and economics. Section III provides an overview of the mathematical principles and formulation of optimal transport-related metrics, their geometric interpretation, and related embedding methods and signal transforms. We also explain Brenier's theorem [22], which helped pave the way for several practical numerical implementation algorithms, which are then explained in detail in section IV. Finally, in section V we review and demonstrate the application of transport-based techniques to numerous problems including: image retrieval, registration and morphing, color and texture analysis, image denoising and restoration, morphometry, super resolution, and machine learning. As mentioned above, software implementing the examples shown can be downloaded from [134].

II. A BRIEF HISTORICAL NOTE
The optimal mass transport problem seeks the most efficient way of transforming one distribution of mass to another, relative to a given cost function. The problem was initially studied by the French Mathematician Gaspard Monge in his seminal work "Mémoire sur la théorie des déblais et des remblais" [109] in 1781. Between 1781 and 1942, several mathematicians made valuable contributions to the mathematics of the optimal mass transport problem among which are Arthur Cayley [27], Charles Duplin [44], and Paul Appell [5]. In 1942, Leonid V. Kantorovich, who at that time was unaware of Monge's work, proposed a general formulation of the problem by considering optimal mass transport plans, which as opposed to Monge's formulation allows for mass splitting [74]. Six years later and in 1948, Kantorovich became aware of Monge's work and made the connection between his work and that of the French mathematician in his short paper, "A Problem of Monge" [75]. Kantorovich shared the 1975 Nobel Prize in Economic Sciences with Tjalling Koopmans for his work in the optimal allocation of scarce resources. Kantorovich's contribution is considered as "the birth of the modern formulation of optimal transport" [163] and it made the optimal mass transport problem an active field of research in the following years [138], [96], [85], [128]. Among the published works in this era is the prominent work of Sudakhov [151] in 1979 on the existence of the optimal transport map. Later on, a gap was found in Sudakhov's work [151] which was eventually filled by Ambrosio in 2003 [1].
A significant portion of the theory of the optimal mass transport problem was developed in the Nineties. Starting with Brenier's seminal work on characterization, existence, and uniqueness of optimal transport maps [22], followed by Cafarrelli's work on regularity conditions of such mappings [25], Gangbo and McCann's work on geometric interpretation of the problem [62], [63], and Evans and Gangbo's work on formulating the problem through differential equations and specifically the p-Laplacian equation [50], [51]. A more thorough history and background on the optimal mass transport problem can be found in Rachev and Rüschendorf's book "Mass Transportation Problems" [129], Villani's book "Optimal Transport: Old and New" [163], and Santambrogio's book " Optimal transport for applied mathematicians" [139], or in shorter surveys such as that of Bogachev and Kolesnikov's [18] and Vershik's [162].
The significant contributions in mathematical foundations of the optimal transport problem together with recent advancements in numerical methods [36], [15], [14], [114], [160] have spurred the recent development of numerous data analysis techniques for modern estimation and detection (e.g. classification) problems. Figure 1 plots the number of papers related to the topic of optimal transport that can be found in the public domain per year demonstrating significant growth in these techniques.

III. FORMULATION OF THE PROBLEM AND METHODOLOGY
In this section we first review the two classical formulations of the optimal transport problem (i.e. Monge's and Kantorovich's formulations). Next, we review the geometrical characteristics of the problem, and finally review the transport based signal/image embeddings.
A. Optimal Transport: Formulation 1) Monge's formulation: The Monge optimal mass transportation problem is formulated as follows. Consider two probability measures µ and ν defined on measure spaces X and Y . In most applications X and Y are subsets of R d and µ and ν have density functions which we denote by I 0 and I 1 , dµ(x) = I 0 (x)dx and dν(x) = I 1 (x)dx, (originally representing the height of a pile of soil/sand and the depth of an excavation). Monge's optimal transportation problem is to find a measurable map f : X → Y that pushes µ onto ν and minimizes the following objective function, where c : X×Y → R + is the cost functional, and M P := {f : X → Y | f # µ = ν} where f # µ represents the pushforward of measure µ and is characterized as, If µ and ν have densities and f is smooth and one to one then the equation above can be written in a differential form as almost everywhere, where Df is the Jacobian of f (See Figure 2, top panel). Such measurable maps f ∈ M P are sometimes called 'transport maps' or 'mass preserving maps'. Simply put, the Monge formulation of the problem seeks the best pushforward map that rearranges measure µ into measure ν while minimizing a specific cost function. Monge considered the Euclidean distance as the cost function in his original formulation, c(x, f (x)) = |x − f (x)|. Note that both the objective function and the constraint in Equation (1) are nonlinear with respect to f (x). Hence, for over a century the answers to questions regarding existence and characterization of the Monge's problem remained unknown. It should be mentioned that, for certain measures the Monge's formulation of the optimal transport problem is illposed; in the sense that there is no transport map to rearrange one measure to another. For instance, consider the case where µ is a Dirac mass while ν is not. Kantorovich's formulation alleviates this problem by finding the optimal transport plan as opposed to the transport map.
2) Kantorovich's formulation: Kantorovich formulated the transportation problem by optimizing over transportation plans, where a transport plan is a probability measure γ ∈ P (X × Y ) with marginals µ and ν. One can think of γ as the joint distribution of I 0 and I 1 describing how much 'mass' is being moved to different coordinates. That is let A be a measurable subset of X and similarly B ⊆ Y . The quantity γ(A, B) tells us how much 'mass' in set A is being moved to set B. Now let γ ∈ P (X × Y ) be a plan with marginals µ and ν, i.e., where Γ X : X × Y → X and Γ Y : X × Y → Y are the canonical projections. Let Γ(µ, ν) be the set of all such plans. Kantorovich's formulation can then be written as, The minimizer of the optimization problem above, γ * , is called the optimal transport plan. Note that unlike the Monge problem, in Kantorovich's formulation the objective function and the constraints are linear with respect to γ(x, y). Moreover, Kantorovich's formulation is in the form of a convex optimization problem. Note that the Monge problem is more restrictive than the Kantorovich problem. That is, in Monge's version, mass from a single location in one domain is being sent to a single location in the domain of the other measure. Kantorovich's formulation considers transport plans which can deal with arbitrary measurable sets and has the ability to distribute mass from the one location in one density to multiple locations in another (See Figure 2, bottom panel). For any transport map f : X → Y there is an associated transport plan, given by Furthermore when an optimal transport map f * exists, with the optimal transport then γ * given by (4) is an optimal transportation plan [163]. Note that the converse does not always hold. Finally, the following relationship holds between Equations (1) and (3) [18], [163], The Kantorovich problem is especially interesting in a discrete setting, that is for probability measures of the form µ = M i=1 p i δ xi and ν = N j=1 q j δ yj , where δ xi is a dirac measure centered at x i , the Kantorovich problem can be written as, where γ ij identifies how much of the mass particle m i at x i needs to be moved to y j (See Figure 2, bottom panel). Note that the optimization above has a linear objective function and linear constraints, therefore it is a linear programming problem. We note that the problem is convex, but not strictly so, and the constraint provides a polyhedral set of M × N matrices, Γ(µ, ν).
In practice, a non-discrete measure is often approximated by a discrete measure and the Kantorovich problem is solved through the linear programming optimization expressed in equation (7). An important result that allows this discretization is the following. If µ h and ν h are sequences of measures that weakly converge to µ and ν respectively and γ h is the optimal transport plan between µ h and ν h then, up to subsequences, γ h converges weakly to the optimal transport plan between µ and ν, γ * (see for example [163,Thm 5.20]). This result enables the application of discrete OT methods to more general measures.
Here we note that the Kantorovich problem can be extended to multi-marginal problems where instead of working in X ×Y one defines χ := X 1 × X 2 × ... × X N , the cost c : χ → R + , and marginals are defined as (Γ Xi ) # γ = µ i . Multi-marginal transport problems and Martingale optimal transport problems [119], [139], [29] are two variations of the transport problem, which have recently attracted attention in economics [119], [72].
3) Dual formulation: The focus of the Kantorovich problem is to minimize the cost of transportation. Villani [163] describes the problem by an intuitive example using the transportation of bread between bakeries and cafés. Consider a company that owns bakeries throughout a city, which produce loaves that should be transported to cafés. If this company were to take delivery into its own hands, then it would like to minimize the transportation costs. Minimizing the cost of the transportation between bakeries and cafés is the focus of Kantorovich's problem. Now imagine that a second company comes into play, which offers to take care of all transportation (maybe they have a few tricks to do the transportation in a very efficient manner) offering competitive prices for buying bread at bakeries and selling them at cafés. Let φ(x) be the cost of bread bought at bakery x, and ψ(y) be its selling price at café y. Before the second company showed up, the first company's selling price was equal to φ(x) + c(x, y), which is the price of the bread plus the delivery fee (transportation cost). Now if the second company wants to be competitive, it needs to sell bread to cafés with a price which is lower or equal to that which the first company offers. Meaning that, At the same time, the second company of course wants to maximize its profit, which can be formulated as, Note that maximizing the profit is the dual problem [76], [129], [163], [139], [18] of the Kantorovich formulation. In addition, since the primal problem is convex we have, It is apparent from this scenario that the optimal prices for the dual problem must be tight meaning that, the second company cannot increase the selling price, ψ(y), or decrease the buying price, φ(x). Formally, The equations above imply that the pair ψ and φ are related through a c-transform (see [163] for the definition and more detail), which coincides with the Legendre-Fenchel transform [20] when c(x, y) = − x, y . Finally, for the quadratic transport cost, c(x, y) = 1 2 |x − y| 2 , and when an optimal transport map exists, the optimal arguments of the dual problem and the optimal transport map are related through, f * (x) = ∇( 1 2 |x| 2 − φ * (x)). 4) Basic properties: The existence of a solution for the Kantorovich problem follows from Ulham's [158] and Prokhorov's [122] theorems. Ulham's theorem states that Γ(µ, ν) is a tight set of probability measures defined on X ×Y . On the other hand, Prokhorov's theorem [122] states that Γ(µ, ν) is relatively compact, and hence the limit of a sequence of transport plans (γ n ) ⊂ Γ(µ, ν) where γ n → γ is also in Γ(µ, ν). Using these theorems one can show the following theorem [2], [163]: For a lower semicontinuous and bounded from below cost function, c(x, y) there exists a minimizer to the Kantorovich problem.
In engineering applications the common cost functions almost always satisfy the existence of a transport plan. A further important question is regarding the existence of an optimal transport map instead of a plan. Brenier [22] addressed this problem for the special case where c(x, y) = |x − y| 2 . Bernier's results was later relaxed to more general cases by Gangbo and McCann [63], which led to the following theorem: Theorem Let µ and ν be two Borel probability measures on compact measurable supports X and Y , respectively. When c(x, y) = h(x−y) for some strictly convex function h and µ is absolutely continuous with respect to the Lebesgue measure, then there exists a unique optimal transportation map f * : In addition, the optimal transport plan is unique, and thus equal to γ(x, y) = (Id × f * ) # µ (See Eq. 4). Moreover, f * is characterized by the gradient of a c-concave function φ : X → R as follows, For a proof see [2], [63], [163]. Note that, h(x) = 1 2 |x| 2 is rather a special case, because the gradient of h is equal to identity,∇h = Id, and the optimal transport map is simply characterized as f * (x) = x − ∇φ(x).
B. Optimal Mass Transport: Geometric properties 1) Wasserstein metric: Let Ω be a subset of R d on which the measures we consider are defined. In most applications Ω is the domain where the signal is defined and thus bounded. Let P p (Ω) be the set of Borel probability measures on Ω, with finite p'th moment, that is the set of probability measures µ on R d such that Ω |x| p dµ(x) < ∞. The p-Wasserstein metric, W p , for p ≥ 1 on P p (Ω) is then defined as using the optimal transportation problem (3) with the cost function c(x, y) = |x − y| p . For µ and ν in P p (Ω), For any p ≥ 1, W p is a metric on P p (Ω). The metric space (P p (Ω), W p ) is referred to as the p-Wasserstein space.
If Ω is bounded then for any p ≥ 1, W p metrizes the weak convergence of measures on P (Ω). That is the convergence with respect to W p is equivalent to weak convergence of measures.
Note that the p-Wasserstein metric can equivalently be defined using the dual Kantorovich problem, For the specific case of p = 1 the p-Wasserstein metric is also known as the Monge-Rubinstein [163] metric, or the earth mover distance [136], which can also be expressed in the following form where Lip(φ) is the Lipschitz constant for φ. The p-Wasserstein metric for one-dimensional probability measures is specifically interesting due to its simple and unique characterization. For one-dimensional probability measures µ and ν on R the optimal transport map has a closed form solution. Let F µ be the cumulative distribution function of a measure µ ∈ P p (R) Note that this is a nondecreasing function going from 0 to 1, which is continuous if µ has a density. We define the pseudoinverse of F µ as follows: for z ∈ (0, 1), If µ has positive density then F µ is increasing (and thus invertible) and the inverse of the function F µ is equal to the pseudoinverse we just defined. In other words the pseudoinverse is a generalization of the notion of the inverse of a function. The pseudoinverse (i.e. the inverse if the densities of µ and ν are positive) provides a closed form solution for the p-Wasserstein distance: The closed-form solution of the p-Wasserstein distance in one dimension is an attractive property, as it alleviates the need for optimization. This property was employed in the Sliced Wasserstein metrics as defined below. 2) Sliced-Wasserstein Metric: The idea behind the Sliced Wasserstein metric is to first obtain a family of onedimensional representations for a higher-dimensional probability distribution through projections (slicing the measure), and then calculate the distance between two input distributions as a functional on the Wasserstein distance of their onedimensional representations. In this sense, the distance is obtained by solving several one-dimensional optimal transport problems, which have closed-form solutions.
Slicing a measure is closely related to the well known Radon transform in the imaging and image processing community [19], [87]. The d-dimensional Radon transform R maps a function |dx ≤ ∞} into the set of its integrals over the hyperplanes of R n and is defined as, (20) here θ ⊥ is the subspace orthogonal to θ, and S d−1 is the unit sphere in R d . Note that R : . One can slice a probability measure µ defined on R × S d−1 into its conditional measures with respect to the uniform measure on S d−1 to obtain a measure µ θ , which satisfies, . The Sliced Wasserstein metric, for continuous probability measures µ and ν on R d is then defined as, where p ≥ 1, and W p is the Wasserstein metric, which for one dimensional measures µ θ and ν θ has a closed form solution (See Figure 3). For more details and definitions of the Sliced Wasserstein metric we refer the reader to [127], [93], [19], [87]. 3) Wasserstein spaces, geodesics, and Riemannian structure: In this section we assume that Ω is convex. Here we highlight that the p-Wasserstein space (P p (Ω), W p ) is not just a metric space, but has additional geometric structure. In particular for any p ≥ 1 and any µ 0 , µ 1 ∈ P p (Ω) there exists a continuous path between µ 0 and µ 1 whose length is the distance between µ 0 and µ 1 .
Furthermore the space with p = 2 is special as it possesses a structure of an formal, infinite dimensional, Riemannian manifold. That structure was first noted by Otto [117] who developed the formal calculations for using this structure. The weak setting in which the notions of tangent space, Riemannian metric were made rigorous was developed by Ambrosio, Gigli and Savaré [3]. A slightly less general, but more approachable introduction is available in [2] and also in [139]. The notion of curvature was developed in [64], [101].
Here we review the two main notions, which have a wide use. Namely we characterize the geodesics in (P p (Ω), W p ) and in the case p = 2 describe what is the local, Riemannian metric of (P 2 (Ω), W 2 ). Finally we state the seminal result of Benamou and Brenier [14] who provided a characterization of geodesics via action minimization which is useful in computations and also gives an intuitive explanation of the Wasserstein metric.
We first recall the definition of the length of a curve in a metric space. Let (X, d) be a metric space and µ : [a, b] → X. Then the length of µ, denoted by L(µ) is A metric space (X, d) is a geodesic space if for any µ 0 and µ 1 there exists a curve µ : In particular the length of µ is equal to the distance from µ 0 to µ 1 . Such curve µ is called a geodesic. The existence of geodesic is useful as it allows one to define the average of µ 0 and µ 1 as the midpoint of the geodesic connection them. An important property of (P p (Ω), W p ) is that it is a geodesic space and that geodesics are easy to characterize. Namely they are given by the displacement interpolation (a.k.a. McCann interpolation). We first describe them for µ 0 in P p (Ω) which has a density: dµ 0 = I 0 dx and arbitrary µ 1 ∈ P p (Ω). Then there exists a unique transportation map f * µ 0 = µ 1 which minimizes the transportation cost for the given p ≥ 1. The geodesic is obtained by moving the mass at constant speed from x to f * (x). More precisely, let t ∈ [0, 1] and x ∈ Ω let be the position at time t of the mass initially at x. Note that f * 0 is identity mapping and f * 1 = f * . Pushing forward the mass by f * t : provides the desired geodesic from µ 0 to µ 1 . We remark that the velocity of each particle ∂ t f * t = f * (x)−x which is nothing but the displacement of the optimal transportation map.
If µ 0 does not have a density, we need to use the optimal transportation plan γ ∈ Γ(µ 0 , µ 1 ) which minimizes the ptransportation cost. For t ∈ [0, 1] let F t : Ω × Ω → Ω be the convex interpolation: F t (x, y) = (1 − t)x + ty. Then µ * (t) = F t γ is the desired geodesic. The intuitive explanation is that mass from x which ends up at y is being transported at constant speed, along the line segment from x to y. Figure 4 conceptualizes the geodesic between two measures in P 2 (Ω), and visualizes it for three different pairs of measures.
An important fact regarding the 2-Wasserstein space is Otto's presentation of a formal Riemannian metric for this space [117]. It involves shifting to Lagrangian point of view. To explain, consider path µ(t) in P 2 (Ω) with smooth densities I(x, t). Then s(x, t) = ∂I ∂t (x, t), the perturbation of µ(t), can be thought as a tangent vector. Instead of thinking of increasing/decreasing the density this perturbation can be viewed as resulting form moving the mass by a vector field. In other words consider vector fields v(x, t) such that There are many such vector fields. Otto defined the size of s( · , t) as the (square root of) the minimal kinetic energy of the vector field that produces the perturbation to density s. That is s, s = min v satisfies (23) Utilizing the Riemmanian manifold structure of P 2 (Ω) together with the inner product presented in Equation (24) the 2-Wasserstein metric can be reformulated into finding the minimizer of the following action among all curves in P 2 (Ω) connecting µ and ν [14], where the first constraint is the continuity equation.

C. Optimal Transport: Embeddings and Transforms
The optimal transport problem and specifically the 2-Wasserstein metric and the Sliced 2-Wasserstein metric have been recently used to define bijective nonlinear transforms for signals and images [165], [133], [87], [89]. In contrast to commonly used linear signal transformation frameworks (e.g. Fourier and Wavelet transforms) which only employ signal intensities at fixed coordinate points, thus adopting an 'Eulerian' point of view, the idea behind the transportbased transforms is to consider the intensity variations together with the locations of the intensity variations in the signal. Therefore, such transforms adopt a 'Lagrangian' point of view for analyzing signals. Here we briefly describe these transforms and some of their prominent properties.
1) The linear optimal transportation framework: The linear optimal transportation (LOT) framework was proposed by Wang et al. [165]. The framework was successfully used in [11], [155], [118], [156] for pattern recognition in biomedical images and specifically histopathology and cytology images. Later, it was extended in [89] as a generic framework for pattern recognition and was used in [88] for single-frame super-resolution reconstruction of face images. The framework was modified by Seguy and Cuturi [145], where they studied the principal geodesics of probability measure. The LOT framework provides an invertible Lagrangian transform for images. It was initially proposed as a method to simultaneously amend the computationally expensive requirement of calculating pairwise 2-Wasserstein distances between N probability measures for pattern recognition purposes, and to allow for the construction of generative models for images involving Fig. 5. Graphical representation of the LOT framework. The framework embeds the probability measures ν i in the tangent space of P 2 (Ω) at a fixed probability measure µ with density I 0 . As a consequence, the Euclidean distance between the embedded functionsĨ 1 andĨ 2 provides an approximation for the 2-Wasserstein distance, W 2 (ν 1 , ν 2 ).
textures and shapes. For a given set of probability measures ν i ∈ P 2 (Ω), for i = 1, ..., N , and a fixed template probability measure µ, the transform projects the probability measures to the tangent space at µ, T µ . The projections are acquired by finding the optimal velocity fields corresponding to the optimal transport plans between µ and each probability measure in the set.
Formally, the framework provides a linear embedding for P 2 (Ω) with respect to a fixed measure µ ∈ P 2 (Ω). Meaning that, the Euclidean distance between the embedded measures, ν i , and the fixed measure, µ, is equal to W 2 (µ, ν i ) and the Euclidean distance between two embedded measures is, generally speaking, an approximation of their 2-Wasserstein distance. The geometric interpretation of the LOT framework is presented in Figure 5. The linear embedding then facilitates the application of linear techniques such as principal component analysis (PCA) and linear discriminant analysis (LDA) to probability measures.
2) The cumulative distribution transform: Park et al. [133] considered the LOT framework for one-dimensional probability measures, and since in dimension one the transport maps are explicit, they were able to characterize the properties of the transformed densities in the tangent space. Here, we briefly review their results. Similar to the LOT framework, let ν i for i = 1, ..., N , and µ be absolutely continuous probability measures defined on R, with corresponding positive probability densities I i and I 0 . The framework first calculates the optimal transport maps between I i and I 0 using Then the forward and inverse transport-based transform, denoted as the cumulative distribution transform (CDT) by Park et al. [133], for these density functions with respect to the fixed template µ is defined  as, where ). Note that the L 2 norm of the transformed signals,Ĩ i corresponds to the 2-Wasserstein distance between µ and ν i . In contrast to the higher-dimensional LOT, the Euclidean distance between two transformed (embedded) signalsĨ i andĨ j , however, is the exact 2-Wasserstein distance between ν i and ν j (See [133] for a proof) and not just an approximation. Hence, the transformation is isometric with respect to the 2-Wasserstein metric. This isometric nature of the CDT was utilized in [90] to provide positive definite kernels for n-dimensional measures.
From a signal processing point of view, the CDT is a nonlinear signal transformation that captures certain nonlinear variations in signals including translation and scaling. Specifically, it gives rise to the transformation pairs presented in Table I. From Table I one can observe that although I(t−τ ) is nonlinear in τ , its CDT representationĨ(t)+τ I 0 (t) becomes affine in τ (similar effect is observed for scaling). In effect, the Lagrangian transformations (compositions) in original signal space are rendered into Eulerian perturbations in transform space, borrowing from the PDE parlance. Furthermore, Park et al. [133] demonstrated that the CDT facilitates certain pattern recognition problems. More precisely, the transformation turns certain not linearly separable and disjoint classes of signals into linearly separable ones. Formally, let C be a set of measurable maps and let P, Q ⊂ P 2 (Ω) be sets of positive probability density functions born from two positive probability density functions p 0 , q 0 ∈ P 2 (Ω) (mother density functions) as follows, The sets P and Q are disjoint but not necessarily linearly separable in the signal space. A main result of [133] states that the signal classes P and Q are guaranteed to be linearly separable in the transform space (regardless of the choice of the reference signal I 0 ) if C satisfies the following conditions, Fig. 7. The cumulative percentage of the face dataset in Figure 6 in the image space, the Radon transform space, the Ridgelet transform space, and the Radon-CDT transform space.
The set of translations C = {f |f (x) = x + τ, τ ∈ R}, and scaling C = {f |f (x) = ax, a ∈ R + }, for instance, satisfy above conditions. We refer the reader to [133] for further reading. The top panel in Figure 6 demonstrates the linear separation property of the CDT and demonstrate the capability of such nonlinear transformation. The signal classes P and Q are chosen to be the set of all translations of a single Gaussian and a Gaussian mixture including two Gaussian functions with a fixed mean difference, respectively. The discriminant subspace is calculated for these classes and it is shown that while the signal classes are not linearly separable in the signal domain, they become linearly separable in the transform domain.
3) The Radon cumulative distribution transform: The CDT framework was extended to higher dimensional density functions through the sliced-Wasserstein distance in [87], and was denoted as the Radon-CDT. It is shown in [87] that similar characteristics of the CDT, including the linear separation property, also hold for the Radon-CDT. Figure 6 clarifies the linear separation property of the Radon-CDT and demonstrate the capability of such transformations. Particularly, Figure 6 shows a facial expression dataset with two classes (i.e. neutral and smiling expressions) and its corresponding representations in the LDA discriminant subspace calculated from the images (bottom left panel), the Radon-CDT of the dataset and the corresponding representation of the transformed data in the LDA discriminant subspace (bottom right panel). It is clear that the image classes become more linearly separable in the transform space. In addition, the cumulative percentage variation of the dataset in the image space, the Radon transform space, the Ridgelet transform space, and the Radon-CDT space are shown in Figure 7. It can be seen that the variations in the dataset could be explained with fewer components in the nonlinear transform spaces as opposed to the linear ones.
IV. NUMERICAL METHODS There exists a variety of fundamentally different approaches to finding optimal transportation maps and plans. Below we present the notable approaches. In the table II we summarize the expected computational complexity of the given algorithm.  Monge-Ampere solver N Quadratic cost. One in [16] is proved to be convergent. Accuracy is an issue due to wide stencil used. A. Linear programming Leonid Kantorovich's extensive monologue on 'Mathematical Methods in the Organization and Planning of Production' in 1939 [76], George B Dantzig's pioneering work on the simplex algorithm in 1947 [38], John von Neumann's duality theory in 1947 [112], and Tjalling Koopmans' work on 'Activity analysis of production and allocation' in 1951 [55] were the founding pillars of the linear programming as we know it today. The term 'linear programming', was coined by Koopmans in [55]. The linear programming problem, is an optimization problem with a linear objective function and linear equality and inequality constraints. The set of feasible points to a linear programming problem forms a possibly unbounded convex set. The optimal solution can be thought of as the intersection of the hyper plane corresponding to the linear objective function with an extreme point (i.e. a corner point) of the feasible convex set of solutions.
Several numerical methods exist for solving linear programming problems [161], among which are the simplex method [38] and its variations [35], [39], [13] and the interior-point methods [78], [102]. The computational complexity of the mentioned numerical methods, however, scales at best cubically in the size of the domain. Hence, assuming the density measures have N particles the number of unknowns γ ij s is N 2 and the computational complexities of the solvers are at best O(N 3 logN ) [136], [36]. The computational complexity of the linear programming methods is a very important limiting factor for the applications of the Kantorovich problem.
We point out that in the special case where M = N and the mass is equidistributed, the optimal transport problem simplifies to a one to one assignment problem, which can be solved by the auction algorithm of Bertsekas [17] or similar methods in O (N 2 logN ). In addition, Pele and Werman showed that a thresholded version of the Kantorovich problem can be solved through the min-cost-flow algorithm in O(N 2 ) [121]. Several multiscale approaches and sparse approximation approaches have recently been introduced to improve the computational performance of the linear programming solvers, including Schmitzer [143], ]. The work by Oberman and Ruan [114] provides an algorithm which is claimed to be of linear complexity for the linear programming problem. Here we provide a brief description of this method.
1) Multi-Scale algorithm: For a coarse grid of size h, the space Ω ⊂ R d is uniformly discretized into hypercubes with centers x i presented with the set Ω h = {x i } N i=1 and where the distance between neighboring points is h. Measures µ and ν are discretized by One then performs the linear program (7) with inputs p i = p h i , q i = q h i , x i = x h i and y i = y h i to calculate the optimal map γ h between µ h and ν h .
In the next step, the discretization is refined to h 2 . Since γ h is expected to converge to the optimal transport plan between µ and ν, γ * , then γ h 2 should be close to γ h . Let S h be the support of γ h on the grid Ω h × Ω h andS h be the set of neighbors (in space) of S h . The Oberman and Ruan method assumes that is the projection onto the refined grid Ω 2 h 2 (where one point in Ω h is projected to 2 d points in Ω h 2 ). One then solves the linear program on the refined grid. For better resolution repeat the refinement procedure. The reason that the multi-scale procedure is advantageous is that the linear programming problem solved at the finest level is expected to have O(N ) variables as opposed to O(N 2 ) variables that the original full linear program has.

B. Entropy regularized solution
Cuturi's work [36] provides a fast and easy to implement variation of the Kantorovich problem, which has attracted ample attention recently. Cuturi [36] proposed a variation of the Kantorovich problem by considering the transportation problem from a maximum-entropy perspective. He suggested to regularize the Wasserstein metric by the entropy of the transport plan. This modification simplifies the problem and enables much faster numerical schemes with complexity O(N 2 ) [36] or O (N log(N )) using the convolutional Wasserstein distance Fig. 8. The geometric interpretation of the entropy regularization of the Kantorovich problem. For λ = 0 the optimal transport plan is located at an extreme point (corner point) of the feasible set Γ(µ, ν). Increasing λ corresponds to constraining the space to the convex subset of Γ(µ, ν), for which the KL-divergence between the transport plans and γ * ∞ = µ × ν is smaller than 1 λ . The entropy regularized optimal transport plan, γ * λ is shown for a one-dimensional example with different values of λ.
presented in [147] (compared to O(N 3 ) of the linear programming methods), where N is the number of delta masses in each of the measures. The price one pays is that it is difficult to obtain high accuracy approximations of the optimal transport plan. Formally, the entropy regularized p-Wasserstein distance, otherwise known as the Sinkhorn distance, between probability measures µ and ν defined on the metric space (Ω, d) is defined as, W p p,λ (µ, ν) = inf γ∈Γ(µ,ν) Ω×Ω d p (x, y)γ(x, y)dxdy − λh(γ) (28) where h(γ) is the entropy of the plan and is defined as, h(γ) = − Ω×Ω γ(x, y) ln(γ(x, y))dxdy.
We note that this is not a true metric since W p p,λ (µ, µ) > 0. Since the entropy term is strictly concave, the overall optimization problem in (28) becomes strictly convex. It is straightforward to show that ( see [36]) the entropy regularized p-Wasserstein distance in Equation (28) can be reformulated as, where K λ (x, y) = exp(− d p (x,y) λ ) and KL(γ|K λ ) is the Kullback-Leibler (KL) divergence between γ and K λ and is defined as, Figure 8 provides a geometric interpretation of the entropy regularization as discussed in [36]. It can be seen that, the regularizer enforces the plan to be within 1 λ radius in the KLdivergence sense from the transport plan γ * ∞ = µ × ν. In the discrete setting, where µ = N i=1 p i δ xi and ν = N j=1 q j δ yj , γ and K λ are N × N matrices and hence the number of unknowns is N 2 . However, using Equation (31) it can be shown that the optimal transport plan γ is of the form D v K λ D w where D v and D w are diagonal matrices with diagonal entries v, w ∈ R N [36]. Therefore, due to the new formulation the number of unknowns decreases from N 2 to 2N . The new problem can then be solved through the iterative proportional fitting procedure (IPFP) [43], [137], otherwise known as the RAS algorithm [9], or alternatively through the Sinkhorn-Knopp algorithm [146], [84]. The Sinkhorn-Knopp algorithm is specifically interesting as it is computationally efficient and very simple to implement.
The entropy regularization of the transportation problem has recently attracted ample attention from both application and numerical analysis point of view. One can refer to Benamou et al. [15] and Solomon et al. [147] for further numerical analysis, and Cuturi and Doucet [37], Rabin et al. [125], and Rabin and Papadakis [124] and others for applications of the entropy regularized transportation problem. It should be mentioned, however, that the number of iterations for the entropy regularized solvers increases as the regularization parameter goes to zero [114]. In addition, these method become numerically unstable for small regularization parameters [15], [114]. This is due to the fact that the elements of matrix K λ go to zero as λ goes to zero.

C. Fluid dynamics solution
Benamou and Brenier [14] presented a fluid dynamic formulation of the 2-Wassestein metric using the continuity equation as discussed in Section III-B3 and Equation (25). They provided a numerical scheme for optimizing the problem in (25), using the augmented Lagrangian method [56]. A brief overview of the method and its derivation is as follows. Let φ(x, t) be the space-time dependent Lagrange multiplier for constraints in Equation (25). Using φ(x, t) the Lagrangian is given by, where m = Iv and we used integration by parts together with the equality constraints in Equation (25) to obtain the second line. Note that, Using the Legendre transform [20] on |m| 2 /2I one can write, Define ψ := {I, m}, q := {a, b}, and their corresponding inner product to be, then it is straightforward to show that Equation (33) can be rewritten as, where ∇ t,x φ = {∂ t φ, ∇ x φ}. In the formulation above, ψ can be treated as the Lagrange multiplier for a new constraint, namely ∇ t,x φ = q. Thus the augmented Lagrangian can be written as, and finally the corresponding saddle point problem is, Benamou and Brenier [14] used a variation of the Uzawa algorithm [47] to solve the problem above. We note that recent methods based on Bregman iteration such as [166], [49] could also be used for solving the saddle point problem in (37). Due to the space-time nature of the fluid dynamic formulation, the solvers require a space-time discretization (spatialtemporal grid), which increases the computational complexity of such solvers. However, the fluid dynamics solution enables one to handle situations for which there exist barriers in the domain. Take for instance transportation of a distribution through a maze, in which the mass can not be transported over the maze walls.

D. Flow minimization (AHT)
Angenent, Haker, and Tannenbaum (AHT) [4] proposed a flow minimization scheme to obtain the optimal transport map from the Monge problem. The method was used in several image registration applications [71], [167], pattern recognition [89], [165], and computer vision [88]. The method was polished in several follow up papers [71], [159]. The idea behind the method, is to first obtain an initial mass preserving transport map using the Knothe-Rosenblatt coupling [135], [163] and then update the initial map to obtain a curl free mass preserving transport map that minimizes the transport cost. A brief review of the method is provided here.
Let µ and ν be continuous probability measures defined on convex domains X, Y ⊆ R d with corresponding positive densities I 0 and I 1 . In order to find the optimal transport map, f * , AHT starts with an initial transport map, f 0 : X → Y calculated from the Knothe-Rosenblatt coupling [135], [163], for which (f 0 ) # µ = ν. Then it updates f 0 to minimize the transport cost. The goal, however, is to update f 0 in a way that it remains a transport map from µ to ν. AHT defines s(x, t), where for a fixed time, t 0 , s(x, t 0 ) : X → X is a transport map from µ to itself. The initial transport map is then updated through s(x, t), starting from s(x, 0) = x, such that f 0 (s −1 (x, t)) minimizes the transport cost. Following simple calculations, one can show (see [71]) that for s(x, t) to be a MP mapping from I 0 to itself, ∂s ∂t should have the following form, for some vector field ξ on X with div(ξ) = 0 and ξ, n = 0 on the boundary, where n is the vector normal to the boundary. From (38), it follows that the time derivative of f (x, t) = f 0 (s −1 (x, t)) satisfies, AHT then differentiate the Monge objective function, with respect to t, which after rearranging will lead to, Substituting (38) in the above equation we get, Writing the Helmholtz decomposition of f as f = ∇φ + χ and using the divergence theorem we get, thus, ξ = χ decreases M the most. In order to find χ, AHT first find φ and then subtract its gradient from f . Given that div(χ) = 0 and χ, n = 0 and taking the divergence of the Helmholtz decomposition of f leads to, ∇φ, n = f, n on the boundary (44) where ∆ is the Laplace operator. Therefore, is the solution to (44). Substituting ξ back into (39) we have, which for f ∈ R 2 can be even simplified further, where ⊥ indicates rotation by 90 degrees. Finally, the transport map f is updated with a gradient descent scheme, where is the gradient descent step size. AHT show that for infinitesimal step size, , f (x, t) converges to the optimal transport map. Fig. 9. Visualization of the iterative update of the transport potential and correspondingly the transport displacement map through CWVP iterations.

E. Gradient descent on the dual problem
Chartrand, Wohlberg, Vixie, and Bollt (CWVB) [28] provide a straight-forward algorithm for computing the optimal transport. It should be mentioned, however, that the method has not been studied in depth. The method reformulates the Kantorovich's dual formulation as described in Eq. 9 into minimization of a convex, continuous functional. More importantly, the derivative of the continuous functional can be explicitly found, and thus a gradient descent approach is developed for computing the optimal transport. A brief review of the method and its derivation is presented here.
For the strictly convex cost function, c(x, y) = 1 2 |x − y| 2 , the dual problem can be written as, where η c (y) := max x∈X (x·y −η(x)) is the Legendre-Fenchel transform of η(x). Note the relationship between η and φ in Section III-A3 is φ(x) = 1 2 x 2 − η(x). The method utilizes the following property of the Legendre-Fenchel transform [20], and derives the functional derivative of M as, where Hη cc is the Hessian matrix of η cc , I is the identity matrix, and η cc (x) = max y∈Y (y · x − η c (y)). Note that if η is a concave function (respectively if φ is a convex function) then η cc = η (and φ cc = φ). CWVB initializes η 0 (x) = 1 2 |x| 2 and updates the potential field η through gradient descent, where k is the step size and is calculated with a line search method. We note that when µ and ν are absolutely continuous probability measures the optimal transport map is obtained from η * through, f * (x) = ∇η * (x). Figure 9 visualizes the iterations of the CWVB method for two face images taken from YaleB face database (normalized to integrate to 1) [94]. We note that the number of iterations required for convergence of the CWVB method depends on the smoothness of the input densities. In [28] the authors proposed a multi-scale scheme (i.e. using Gaussian pyramids) to update η through different levels of smoothness and hence decrease the number of iterations needed for convergence.

F. Monge-Ampere equation
The Monge-Ampére partial differential equation (PDE) is defined as, for some functional h, and where Hφ is the Hessian matrix of φ. The Monge-Ampére PDE is closely related to the Monge problem. According to Bernier's theorem (discussed in Section III-A4) when µ and ν are absolutely continuous probability measures on sets X, Y ⊂ R n , the optimal transport map than minimizes the 2-Wasserstein metric is uniquely characterized as the gradient of a convex function φ : X → Y . Moreover, we showed that the mass preserving constraint of the Monge problem can be written as det(Df )I 1 (f ) = I 0 . Combining these results one can write, where D∇φ = Hφ, and therefore the equation shown above is in the form of the Monge-Ampére PDE. Now if φ is a convex function on X satisfying ∇φ(X) = Y and solving the Equation (53) then f * = ∇φ is the optimal transportation map from µ to ν. The geometrical constraint on this problem is rather unusual in PDEs and is often referred to as the optimal transport boundary conditions. The equation (53) is a nonlinear elliptic equation and thus one hopes that methods achieving complexity O(N ) (as do the solvers for the Laplace equation) can be designed. Several authors have proposed numerical methods to obtain the optimal transport map through solving the Monge-Ampére PDE in Equation (53) [100], [157], [26], [140], [152], [16]. In particular the scheme in [16] is monotone, has compexity O(N ) (up to logarithms) and is provably convergent. We conclude by remarking that several regularity results on the optimal transport maps, including the pioneering work of Caffarreli [25], and the more recent results by De Philippis and Figalli [40] were established through the Monge-Ampére equation.

G. Other approaches
It is also worth pointing out the approach to optimal transport with quadratic cost using the semi-discrete approximation. Namely several works [8], [106], [83], [97] have considered the situation in which one of the measures considered has density dµ = I 0 dx, while the other is a sum of delta masses, µ = q i δ yi . It turns out that there exists weights w i such that the optimal transport map x → y can be described via a power diagram. More precisely the set of x mapping to y i is the following cell of the power diagram: The main observation (originating from [8] and stated in [106]) is that the weights w i are minimizers of the following unconstrained convex functional Works of Kitagawa, Mérigot, and Thibert [83], Merigot [106], and Levy [97] use Newton based schemes and multiscale approaches to minimize the functional. The need to integrate over the power diagram makes the implementation somewhat geometrically delicate. Nevertheless recent implementation by Lévy [97] gives impressive results in terms of speed. We also note that this approach provides the transportation mapping (not just the approximation of a plan).
Finally, we mention a recent work of Aude, Cuturi, Peyré and Bach who investigated the possibility of using stochastic gradient descent on several formulations of OT problem with quadratic cost on large data sets.

V. APPLICATIONS
In this section we review some recent applications of the optimal transport problem in signal and image processing, computer vision, and machine learning.

A. Image retrieval
One of the earliest applications of the optimal transport problem was in image retrieval. Rubner, Tomasi, and Guibas [136] employed the discrete Wasserstein metric, which they denoted as the Earth Mover's Distance (EMD), to measure the dissimilarity between image signatures. In image retrieval applications, it is common practice to first extract features (i.e. color features, texture feature, shape features, etc.) and then generate high dimensional histograms or signatures (histograms with dynamic/adaptive binning), to represent images. The retrieval task then simplifies to finding images with similar representations (i.e. small distance between their histograms/signatures). The Wasserstein metric is specifically suitable for such applications as it can compare histograms/signatures of different sizes (histograms with different binning). This unique capability turns the Wasserstein metric into an attractive candidate in image retrieval applications [136], [121], [98]. In [136], the Wasserstein metric was compared with common metrics such as the Jeffrey's divergence, the χ 2 statistics, the L 1 distance, and the L 2 distance in an image retrieval task; and it was shown that the Wasserstein metric achieves the highest precision/recall performance amongst all.
Speed of computation is an important practical consideration in image retrieval applications. For almost a decade, the high computational cost of the optimal transport problem overshadowed its practicality in large scale image retrieval applications. Recent advancements in numerical methods including the work of Pele and Werman [121], Merigot [106], and Cuturi [36], among many others, have reinvigorated optimal transport-based distances as a feasible and appealing candidate for large scale image retrieval problems.

B. Registration and Morphing
Image registration deals with finding a common geometric reference frame between two or more images and it plays an important role in analyzing images obtained at different times or using different imaging modalities. Image registration and more specifically biomedical image registration is an active research area. Registration methods find a transformation f that maximizes the similarity between two or more image representations (e.g. image intensities, image features, etc.). Comprehensive recent surveys on image registration can be found in [149], [34], [115]. Among the plethora of registration methods, nonrigid registration methods are especially important given their numerous applications in biomedical problems. They can be used to quantify the morphology of different organs, correct for physiological motion, and allow for comparison of image intensities in a fixed coordinate space (atlas). Generally speaking, nonrigid registration is a nonconvex and non-symmetric problem, with no guarantee on existence of a globally optimal transformation.
Haker et al. [70] proposed to use the Monge problem for image warping and elastic registration. Utilizing the Monge problem in an image warping/registration setting has a number of advantages. First, the existence and uniqueness of the global transformation (the optimal transport map) is known. Second, the problem is symmetric, meaning that the optimal transport map for warping I 0 to I 1 is the inverse of the optimal transport map for warping I 1 to I 0 . Lastly, it provides a landmark-free and parameter-free registration scheme with a built-in mass preservation constraint. These advantages motivated several follow-up work to investigate the application of the Monge problem in image registration and warping [69], [71], [167], [159], [111], [68]. Figure 10 shows an example of morphing a Magnetic Resonance Imaging (MRI) image of a brain into another. The axial slices of three dimensional images are shown, as well as slices of the determinant of the Jacobian of the optimal transport map and the deformation field.
In addition to images, the optimal mass transport problem has also been used in point cloud and mesh registration [103], [54], [93], [150], which have various applications in shape analysis and graphics. In these applications, shape images (2D or 3D binary images) are first represented whether with sets of weighted points (i.e. point clouds), using clustering techniques such as K-means or Fuzzy C-means, or with meshes. Then, a regularized variation of the optimal transport problem is solved to match such representations. The regularization on the transportation problem is often imposed to enforce the neighboring points (or vertices) to remain near to each other after the transformation [54]. As another example, Su et al. [150] proposed a composition of conformal maps and optimal transport maps and introduced the Conformal Wasserstein Shape Space, which enables efficient shape matching and comparison.

C. Color transfer and texture synthesis
Texture mixing and color transfer are among appealing applications of the optimal transport framework in image analysis, graphics, and computer vision. Here we briefly discuss these applications.
1) Color transfer: The purpose of color transfer [132] is to change the color palette of an image to impose the feel and look of another image. The color transfer is generally performed through finding a map, which morphs the color distribution of the first image into the second one. For grayscale images, the color transfer problem simplifies to a histogram matching problem, which is solved through the one-dimensional optimal transport formulation [42]. In fact, the classic problem of histogram equalization is indeed a one-dimensional transport problem [42]. The color transfer problem, on the other hand, is concerned with pushing the three-dimensional color distribution of the first image into the second one. This problem can also be formulated as an optimal transport problem as demonstrated in [126], [123], [52].
A complication that occurs in the color transfer on real images, however, is that a perfect match between color distributions of the images is often not satisfying. This is due to the fact that a color transfer map may not transfer the colors of neighboring pixels in a coherent manner, and may introduce artifacts in the color transferred image. Therefore, the color transfer map is often regularized to make the transfer map spatially coherent [123]. Figure 11 shows a simple example of grayvalue and color transfer via optimal transport framework. It can be seen that the cumulative distribution of the grayvalue and color transferred images are similar to that of the input image.
2) Texture synthesis and mixing: Texture synthesis is the problem of synthesizing a texture image that is visually similar to an exemplar input texture image, and has various applications in computer graphics and image processing [45], [61]. Many methods have been proposed for texture synthesis, among which are synthesis by recopy and synthesis by statistical modeling. Texture mixing, on the other hand, considers the problem of synthesizing a texture image from a collection of input texture images in a way that the synthesized texture provides a meaningful integration of the colors and textures of the input texture images. Metamorphosis is one of the successful approaches in texture mixing, which performs the mixing via identifying correspondences between elementary features (i.e. textons) among input textures and progressively morphing between the shapes of elements. In other approaches, texture images are first parametrized through a tight frame (often steerable wavelets) and statistical modeling is performed on the parameters. Rabin et al. [127], for instance, used the optimal transport framework, specifically the sliced Wasserstein barycenters, on the coefficients of a steerable wavelet tight frame for texture mixing.
Other successful approaches include the random phase and spot noise texture modeling [53], which model textures as stationary Gaussian random fields [61]. Briefly, these methods are based on the assumption that the visual texture perception is based on the spectral magnitude of the texture image. Therefore, utilizing the spectral magnitude of an input image and randomizing its phase will lead to a new synthetic texture image which is visually similar to the input image. Ferradans et al. [53] utilized this assumption together with the Wasserstein geodesics to interpolate between spectral magnitude of texture images, and provide synthetic mixed texture images. Figure 12 shows an example of texture missing via the Wasserstein geodesic between the spectral magnitudes of the input texture images. The in-between images are synthetically generated using the random phase technique. Fig. 12. An example of texture mixing via optimal transport using the method presented in Ferradans et al. [53] D. Image denoising and restoration The optimal transport problem has also been used in several image denoising and restoration problems [95], [126], [153]. The goal in these applications is to restore or reconstruct an image from noisy or incomplete observation. Rabin and Peyré [126], for instance, introduced an optimal transport-based framework to enforce statistical priors for constrained functional minimization. More precisely, they utilized the Sliced Wasserstein metric to measure the adequacy of the statistics of the reconstructed signal/image with the given constraints. In a similar approach, Swoboda and Schörr proposed a variational optimization problem with the Total Variation (TV) as the regularizer and the Wasserstein distance (dual formulation) to enforce prior statistics on reconstructed signal/image. They relaxed the optimization into a convex/concave saddle point problem and solved it via a proximal gradient descent.
In another approach, Lellmann et al. [95] utilized the Kantorovich-Rubinsten discrepancy term together with a Total Variation term in the context of image denoising. They called their method Kantorovich-Rubinstein-TV (KR-TV) denoising. It should be noted that, the Kantorovich-Rubinstein metric is closely related to the 1-Wasserstein metric (for one dimensional signals they are equivalent). The KR term in their proposed functional provides a fidelity term for denoising while the TV term enforces a piecewise constant reconstruction.

E. Transport based morphometry
Given their suitability for comparing mass distributions, transport-based approaches for performing pattern recognition of morphometry encoded in image intensity values have also recently emerged. Recently described approaches for transport-based morphometry (TBM) [165], [11], [89] work by computing transport maps or plans between a set of images and a reference or template image. The transport plans/maps are then utilized as an invertible feature/transform onto which pattern recognition algorithms such as principal component analysis (PCA) or linear discriminant analysis (LDA) can be applied. In effect, it utilizes the LOT framework described earlier in Section III-C1. These techniques has been recently employed to decode differences in cell and nuclear morphology for drug screening [11], and cancer detection histopathology [164], [118] and cytology [156], amongst other applications including the analysis of galaxy morphologies [89], for example.
We note the strong similarity between deformation-based methods which have long been used to analyzed radiological images [66], [73], for example. The difference being that TBM allows for numerically exact, uniquely defined solutions for the transport plans or maps used. That is, images can be matched with little perceptible error. The same is not true in methods that rely on registration via the computation of deformations, given the significant topology differences commonly found in medical images. Moreover, TBM allows for comparison of the entire intensity information present in the images (shapes and textures), while deformation-based methods are usually employed to deal with shape differences. Figure 13 shows a schematic of the TBM steps applied to a cell nuclei dataset. It can be seen that the TBM is capable of modeling the variation in the dataset. In addition, it enables one to visualize the classifier, which discriminates between image classes (in this case malignant versus benign).

F. Super-Resolution
Super-resolution is the process of reconstructing a highresolution image from one or several corresponding lowresolution images. Super-resolution algorithms can be broadly categorized into two major classes namely "multi-frame" super resolution and "single-frame" super resolution, based on the number of low-resolution images they require to reconstruct the corresponding high-resolution image. The transport-based morphometry approach was used for single frame super resolution in [88] to reconstruct high-resolution faces from very low resolution input face images. The authors utilized the transport-based morphometry in combination with subspace learning techniques to learn a nonlinear model for the highresolution face images in the training set.
In short, the method consists of a training and a testing phase. In the training phase, it uses high resolution face images and morphs them to a template high-resolution face through optimal transport maps. Next, it learns a subspace for the calculated optimal transport maps. A transport map in this subspace can then be applied to the template image to synthesize a high-resolution face image. In the testing phase, the goal is to reconstruct a high-resolution image from the low-resolution input image. The method searches for a synthetic high-resolution face image (generated from the transport subspace) that provides a corresponding lowresolution image which is similar to the input low-resolution image. Figure 14 shows the steps used in this method and demonstrates reconstruction results.

G. Machine-Learning and Statistics
The optimal transport framework has recently attracted ample attention from the machine learning and statistics communities [148], [31], [60], [87], [91], [110]. Some applications of the optimal transport in these arenas include various Fig. 13. The schematic of the TBM framework. The optimal transport maps between input images I 1 , ..., I N and a template image I 0 is calculated. Next, linear statistical modeling such as principal component analysis (PCA), linear discriminant analysis (LDA), and canonical correlation analysis (CCA) is performed on the optimal transport maps. The resulting transport maps obtained from the statistical modeling step are then applied to the template image to visualize the results of the analysis in the image space.
1) Learning: Transport-based distances have been recently used in several works as a loss function for regression, classification, etc. Montavon, Müller, and Cuturi [110] for instance utilized the dual formulation of the entropy regularized Wasserstein distance to train restricted Boltzmann machines (RBMs). Boltzmann machines are probabilistic graphical models (Markov random fields) that can be categorized as stochastic neural networks and are capable of extracting hierarchical features at multiple scales. RBMs are bipartite graphs which are special cases of Boltzmann machines, which define parameterized probability distributions over a set of d-binary input variables (observations) whose states are represented by h binary output variables (hidden variables). RBMs' parameters are often learned through information theoretic divergences such as KL-Divergence. Montavon et al. [110] proposed an alternative approach through a scalable entropy regularized Wasserstein distance estimator for RBMs, and showed the practical advantages of this distance over the commonly used information divergence-based loss functions.
In another approach, Frogner et al. [60] used the entropy regularized Wasserstein loss for multi-label classification. They proposed a relaxation of the transport problem to deal with unnormalized measures by replacing the equality constraints in Equation (7) with soft penalties with respect to KL-divergence. In addition, Frogner et al. [60] provided statistical bounds on the expected semantic distance between the prediction and the groundtruth. In yet another approach, Kolouri et al. [91] utilized the sliced Wasserstein metric and provided a family of positive definite kernels, denoted as Sliced-Wasserstein Kernels, and showed the advantages of learning with such kernels. The Sliced-Wasserstein Kernels were shown to be effective in various machine learning tasks including classification, clustering, and regression.
Solomon et al. [148] considered the problem of graph-based semi-supervised learning, in which graph nodes are partially labeled and the task is to propagate the labels throughout the nodes. Specifically, they considered a problem in which the labels are histograms. This problem arises for example in traffic density prediction, in which the traffic density is observed for few stop lights over 24 hours in a city and the city is interested in predicting the traffic density in the un-observed stop lights. They pose the problem as an optimization of a Dirichlet energy for distribution-valued maps based on the 2-Wasserstein distance, and present a Wasserstein propagation scheme for semi-supervised distribution propagation along graphs.
2) Domain adaptation : Domain adaptation is one of the fundamental problems in machine learning which has gained proper attention from the machine learning research community in the past decade [120]. Domain adaptation is the task of transferring knowledge from classifiers trained on available labeled data to unlabeled test domains with data distributions that differ from that of the training data. The optimal transport framework is recently presented as a potential major player in domain adaptation problems [54], [31], [32]. Courty, Flamary, and Davis [31], for instance, assumed that there exists a non-rigid transformation between the source and target distributions and find this transformation using an entropy regularized optimal transport problem. They also proposed a label-aware version of the problem in which the transport plan is regularized so a given target point (testing exemplar) is only Fig. 14. In the training phase, optimal transport maps that morph the template image to high-resolution training face images are calculated (i). Principal component analysis (PCA) is used to learn a linear subspace for transport maps for which a linear combination of obtained eigenmaps can be applied to the template image to obtain synthetic face images (ii). A geometric interpretation of the problem is depicted in panel (iii), and reconstruction results are shown in panel (iv). associated with source points (training exemplars) belonging to the same class. Courty et al. [31] showed that domain adaptation via regularized optimal transport outperform the state-of-the-art results in several challenging domain adaptation problems.
3) Bayesian inference: Another interesting and emerging application of the optimal transport problem is in Bayesian inference [46], [81], [131], [116]. In Bayesian inference, one critical step is the evaluation of expectations with respect to a posterior probability function, which leads to complex multidimensional integrals. These integrals are commonly solved through the Monte Carlo numerical integration, which requires independent sampling from the posterior distribution. In practice, sampling from a general posterior distribution might be difficult, therefore the sampling is performed via a Markov Chain which converges to the posterior probability after certain number of steps. This leads to the celebrated Markov Chain Monte Carlo (MCMC) method. The downside of MCMC is that the samples are not independent and hence the convergence of the empirical expectation is slow. El Moselhy and Marzouk [46] proposed a transport-based method that evades the need for Markov chain simulation by allowing direct sampling from the posterior distribution. The core idea in their work is to find a transport map (via a regularized Monge formulation), which pushes forward the prior measure to the posterior measure. Then, sampling the prior distribution and applying the transport map to the samples, will lead to a sampling scheme from the posterior distribution. Similar ideas were used in [81], [116]. Figure 15 shows the basic idea behind these methods. Left panel shows the prior distribution p and the posterior distribution q and the corresponding transport map f that pushes p into q. One million samples, x i , were generated from distribution p and the top-right panel shows the empirical distribution of these samples denoted asp. The bottom-right panel shows the empirical distribution of transformed samples, y i = f (x i ), denoted asq.

4)
Hypothesis testing: Wasserstein distance is used for goodness of fit testing in [41] and for two sample testing in [130]. Ramdas et al. [130] presented connections between the entropy regularized Wasserstein distance, multivariate Energy distance, and the kernel maximum mean discrepancy (MMD) [67], and provided a "distribution free" univariate Wasserstein test statistic. These and other applications of transport-related concepts show the promise of the mathematical modeling technique in the design of statistical data analysis methods to tackle modern learning problems.

VI. SUMMARY AND CONCLUSIONS
Transport-related methods and applications have come a long way. While earlier applications focused primarily in civil engineering and economics problems, they have recently begun to be employed in a wide variety of problems related to signal and image analysis, and pattern recognition. In this tutorial, seven main areas of application were reviewed: image retrieval V-A, registration and morphing V-B, color transfer and texture analysis V-C, image restoration V-D, transportbased morphometry V-E, image super-resolution V-F, and machine learning and statistics V-G. As evidenced by the number of papers published (see Fig 1) per year, transport and related techniques have received increased attention in recent years. Overall, researchers have found that the application of transport-related concepts can be helpful to solve problems in diverse applications. Given recent trends, it seems safe to expect that the the number of different application areas will continue to grow.
In its most general form, the transport-related techniques reviewed in this tutorial can be thought as mathematical models for signals, images, and in general data distributions. Transport-related metrics involve calculating differences not only of pixel or distribution intensities, but also "where" they are located in the corresponding coordinate space (a pixel coordinate in an image, or a particular axis in some arbitrary feature space). As such, the geometry (e.g. geodesics) induced by such metrics can give rise to dramatically different algorithms and data interpretation results. The interesting performance improvements recently obtained could motivate the search for a more rigorous mathematical understanding of transport-related metrics and applications.
We note that the emergence of numerically precise and efficient ways of computing transport-related metrics and geodesics, presented in section IV also serves as an enabling mechanism. Coupled with the fact that several mathematical properties of transport-based metrics have been extensively studied, we believe that the ground is set of their increased use as foundational tools or building blocks based on which complex computational systems can be built. The confluence of these emerging ideas may spur a significant amount of innovation in a world where sensor and other data is becoming abundant, and computational intelligence to analyze these is in high demand. We believe transport-based models while become an important component of the ever expanding tool set available to modern signal processing and data science experts.