Relative-Translation Invariant Wasserstein Distance
Binshuai Wang Qiwei Di Ming Yin Mengdi Wang Quanquan Gu Peng Wei

#1mw:#1

[
    [
    September 9, 2024
=====================




§ ABSTRACT 
We introduce a new family of distances, relative-translation invariant Wasserstein distances (RW_p), for measuring the similarity of two probability distributions under distribution shift. Generalizing it from the classical optimal transport model, we show that RW_p distances are also real distance metrics defined on the quotient set 𝒫_p(ℝ^n)/∼ and invariant to distribution translations. When p=2, the RW_2 distance enjoys more exciting properties, including decomposability of the optimal transport model, translation-invariance of the RW_2 distance, and a Pythagorean relationship between RW_2 and the classical quadratic Wasserstein distance (W_2). Based on these properties, we show that a distribution shift, measured by W_2 distance, can be explained in the bias-variance perspective. In addition, we propose a variant of the Sinkhorn algorithm, named RW_2 Sinkhorn algorithm, for efficiently calculating RW_2 distance, coupling solutions, as well as W_2 distance. We also provide the analysis of numerical stability and time complexity for the proposed algorithm. Finally, we validate the RW_2 distance metric and the algorithm performance with three experiments. We conduct one numerical validation for the RW_2 Sinkhorn algorithm and show two real-world applications demonstrating the effectiveness of using RW_2 under distribution shift: digits recognition and similar thunderstorm detection. The experimental results report that our proposed algorithm significantly improves the computational efficiency of Sinkhorn in certain practical applications, and the RW_2 distance is robust to distribution translations compared with baselines.

















§ INTRODUCTION


















Optimal transport (OT) theory and Wasserstein distance <cit.> provide a rigorous measurement of similarity between two probability distributions. Numerous state-of-the-art machine learning applications are developed based on the OT formulation and Wasserstein distances, including domain adaptation, score-based generative model, Wasserstein generative adversarial networks, Fréchet inception distance (FID) score, Wasserstein auto-encoders, distributionally robust Markov decision processes, distributionally robust regressions, graph neural networks based objects tracking, etc <cit.>. However, the classical Wasserstein distance has major limitations in certain machine learning and computer vision applications. For example, a meteorologist often focuses on identifying similar weather patterns in a large-scale geographical region <cit.>, where he/she cares more about the “shapes” of weather events rather than their exact locations. The weather events are represented as images or point clouds from the radar reflectivity map. Here the classical Wasserstein distance is not useful since the relative location difference or relative translation between two very similar weather patterns will add to the Wasserstein distance value. Another example is the inevitable distribution shift in real-world datasets. A distribution shift may be introduced by sensor calibration error, environment changes between train and test datasets, simulation to real-world (sim2real) deployment, etc. Motivated by these practical use cases and the limitations of Wasserstein distances, we ask the following research question:

Can we find a new distance metric and a corresponding efficient algorithm to measure the similarity between probability distributions (and their supports) regardless of their relative translation?











To answer this research question, we introduce the relative translation optimal transport (ROT) problem and the corresponding relative-translation invariant Wasserstein distance RW_p. We then focus on the quadratic case (p=2) and identify three exciting properties of the RW_2 distance. We leverage these properties to design a variant of the Sinkhorn algorithm to compute RW_2 distance, coupling solutions, as well as W_2 distance. In addition, we provide analysis and numerical experiment results to demonstrate the effectiveness of the new RW_2 distance against translation shifts. Finally, we show the scalability and practical usage of the RW_2 in a real-world meteorological application.






  
Contributions. The main contributions of this paper are highlighted as follows: (a) we introduce a family of new similarity metrics, relative-translation invariant Wasserstein (RW_p) distances, which are real distance metrics like the Wasserstein distance and invariant to the relative translation of two distributions; (b) we identify three useful properties of the quadratic case RW_2 to support our algorithm design: decomposability of the ROT problem, translation-invariance of both the ROT problem solution and the resulted RW_2, and Pythagorean relationship between RW_2 and the classical W_2; (c) we show that the RW_2 can be used to analyze and explain a general distribution shift (measured by W_2) in the bias-variance perspective; and (d) we propose an efficient variant of Sinkhorn algorithm, named the RW_2 Sinkhorn, for calculating RW_2 distance, coupling solutions as well as W_2 distance with significantly reduced computational complexity and enhanced numerical stability under relative translations. Empirically, we report promising performance from the proposed RW_2 distance when the relative translation is large, and the RW_2 Sinkhorn algorithm in illustrative numerical examples and a large-scale real-world task for similar weather detection. Figure <ref> shows our major findings in this work.













  
Notations.
Assume that 𝒫_p(ℝ^n) is the set of all probability distributions with finite moments of order p defined on the space ℝ^n and ℳ(ℝ^n) is the set of all probability distributions with finite supports defined on the space ℝ^n. ℝ_*^m_1 × m_2 represents the set of all m_1× m_2 matrices with non-negative entries. [μ] represents the equivalence class (orbit) of μ under the shift equivalence relation in 𝒫_p(ℝ^n). μ̅ represents the mean of probability distribution μ (data points). m denotes a vector in ℝ^m where all elements are ones. ./ represents the component-wise division.



























  
Related work.
Optimal transport theory is a classical area of mathematics with strong connections to probability theory, diffusion processes and PDEs. Due to the vast literature, we refer readers to <cit.> for comprehensive reviews. Computational OT methods have been widely explored, including Greenkhorn algorithm <cit.>, Network Simplex method <cit.>, Wasserstein gradient flow <cit.>, neural network approximation <cit.>. Significant research has also been conducted on Wasserstein distances, such as the sliced Wasserstein distance <cit.>, Gromov-Wasserstein distance <cit.>, etc. Other important topics include Wasserstein barycenter <cit.> and unbalanced optimal transport <cit.>.



§ PRELIMINARIES






Before delving into the details of our proposed method, it is essential to focus on the groundwork with an introduction to key aspects of classical optimal transport theory and formulations. This foundation will support the subsequent derivations and proofs presented in Section 3.



 §.§ Optimal Transport Theory


The optimal transport theory focuses on finding the minimal-cost transport plans for moving one probability distribution to another probability distribution in a metric space. The core of this theory involves a cost function, denoted as c(x,y), alongside two probability distributions, μ(x) and ν(y). The optimal transport problem is to find the transport plans (coupling solutions) that minimize the cost of moving the distribution μ(x) to ν(y), under the cost function c(x,y). Although the cost function can take any non-negative form, our focus will be on those derived from the p-norm, expressed as x-y_p^p for p ≥ 1, because the optimal transport problem is well-defined <cit.>.


Assuming μ(x) as the source distribution and ν(y) as the target distribution, μ, ν∈𝒫_p(ℝ^n), we can formulate the optimal transport problem as a functional optimization problem, detailed below:
[p-norm optimal transport problem <cit.>]
    
    OT (μ, ν, p) =γ∈Γ∈ℝ^n min∫_ℝ^nx-y_p^p d γ(x,y),

with Γ(μ, ν) = {γ∈𝒫_p(ℝ^n×ℝ^n)|∫_ℝ^nγ(x,y) dx = ν(y), ∫_ℝ^nγ(x,y) dy  = μ(x),  γ(x,y) ≥ 0}.

Here γ(x,y) represents the transport plan (or the coupling solution), indicating the amount of probability mass transported from source support x to target support y. The objective function is to minimize the total transport cost, which is the integrated product cost of distance and transported mass across all source-target pairs (x, y). 




After the foundational optimal transport problem is outlined, we can introduce a family of real metrics, the Wasserstein distances, for measuring the distance between probability distributions on the set 𝒫_p(ℝ^n). These distances are defined based on the optimal transport problem. 
[Wasserstein distances <cit.>]
    The Wasserstein distance between μ and ν is the pth root of the minimal total transport cost from μ to ν, denoted as W_p, (p≥ 1):
    
    W_p(μ, ν) = OT(μ, ν, p)^1/p.

The Wasserstein distance is a powerful tool for assessing the similarity between probability distributions. It is a real metric admitting the properties of indiscernibility, non-negativity, symmetry, and triangle inequality <cit.>. Meanwhile, it is well-defined for any probability distribution pairs, including discrete-discrete, discrete-continuous, and continuous-continuous.

For practical machine learning applications, the functional optimization described in Equation (<ref>) can be adapted into a discrete optimization framework. This adaptation involves considering the distributions, μ and ν, as comprised of finite supports, {x_i}_i=1^m_1 and {y_j}_j=1^m_2, with corresponding probability masses {a_i}_i=1^m_1 and {b_j}_j=1^m_2, respectively, where m_1 and m_2 are the number of supports (data points). Since all m_1 and m_2 are finite numbers, we can use an m_1 × m_2 matrix C to represent the cost between supports, where each entry represents the transporting cost from x_i to y_j, i.e., C_ij = x_i - y_j_p^p. This discrete version of the optimal transport problem can then be expressed as a linear programming problem, denoted as OT(μ, ν, p):

    OT(μ, ν, p) = P_ij∈Π(μ, ν)min∑_i=1^m_1∑_j=1^m_2 C_ij P_ij,

with Π(μ, ν) = {P ∈ℝ_*^m_1 × m_2 | P m_1 = a, P^⊤m_2 = b}, 
where Π(μ, ν) is the feasible set of this problem, vectors a and b are the probability masses of μ and ν, respectively. coupling solutions P_ij indicates the amount of probability mass transported from the source point x_i to the target point y_j. This linear programming approach provides a scalable and efficient way for solving discrete optimal transport problems in various data-driven applications.



 §.§ Sinkhorn Algorithm


Equation (<ref>) formulates a linear programming problem, which is commonly solved by simplex methods or interior-point methods <cit.>. Because of the special structure of the feasible set Π(μ, ν), another approach for solving this problem is to transform it into a matrix scaling problem by adding an entropy regularization in the objective function <cit.>. The matrix scaling problem can be solved by the Sinkhorn algorithm, which is an iterative algorithm that enjoys both efficiency and scalability. In detail, the Sinkhorn algorithm will initially assign u^(0) and v^(0) with vector m_1 and m_2, then the vector u^(k) and v^(k) (k≥ 1) are updated alternatively by the following equations:

    u^(k+1)← a./K v^(k),    v^(k+1)← b./K^⊤ u^(k+1),


where K_ij = e^-C_ij/λ (λ is the coefficient of the entropy regularized term) and the division is component-wise. When the convergence precision is satisfied, the coupling solution P will be calculated by the matrix diag(u) K diag(v). It has been proved the solution calculated by the Sinkhorn algorithm can converge to the exact coupling solution of the linear programming model, as λ goes to zero <cit.>. One caveat of this calculation is the exponent operation, which may cause “division by zero", we will show how we can improve the numerical stability in Section <ref>.




§ RELATIVE TRANSLATION OPTIMAL TRANSPORT AND RW_P DISTANCES


Here we present the relative translation optimal transport model and the RW_p distances. We first introduce the general case of the relative translation optimal transport problem and the RW_p distances. We will then focus on the quadratic RW_2 case and its properties. 




 §.§ Relative Translation Optimal Transport (ROT) Formulation and RW_p Distances


As discussed in Section 1, the classical optimal transport (OT) problem is not useful when there is a relative translation between two distributions (or the two datasets known as their supports). We introduce the relative translation optimal transport problem, ROT(μ, ν, p), which is formulated to find the minimal total transport cost under any translation. 
[Relative translation optimal transport problem]
    Continuing with the previous notations,
    
    ROT(μ, ν, p) = s ∈ℝ^n infP ∈Π(μ, ν)min∑_i=1^m_1∑_j=1^m_2x_i + s - y_j _p^p P_ij,

    where variable s represents the translation of source distribution μ, and variables P_ij represent the coupling solution between the support x_i and the support y_j.



The ROT problem can be viewed as a generalized form of the classical OT in Equation (<ref>). Different from the one-stage optimization in the classical OT problem, there are two stages in this optimization. The inner optimization is exactly the classical OT, whereas the outer optimization finds the best relative translation for the source distribution to minimize the total transport cost.






    For Equation (<ref>), the domain of the variable s can be confined to a compact set {s∈ℝ^n |  s_p ≤ 2 max_ijx_i-y_j_p}. Thus, we have
    
    ROT(μ, ν, p) = s ∈ℝ^n minP ∈Π(μ, ν)min∑_i=1^m_1∑_j=1^m_2x_i + s - y_j _p^p P_ij,

    where the minimum can be achieved.

The proof of Theorem <ref> is provided in Appendix A.


From the perspective of group theory, we could have a better view of which space the ROT problem is defined on. Assume that ∼ is the translation relation on the set 𝒫_p(ℝ^n). When distribution μ can be translated to distribution μ', we denote it by μ∼μ'. Because the translation is an equivalence relation defined on the set 𝒫_p(ℝ^n), we may divide set 𝒫_p(ℝ^n) by the translation relation, which leads to a quotient set, 𝒫_p(ℝ^n)/∼. 𝒫_p(ℝ^n)/∼ consists of the equivalence class of distributions, and each equivalence class, denoted by [μ], contains all mutually translatable probability distributions. Therefore, the ROT problem can also be regarded as an OT problem defined on the quotient set, 𝒫_p(ℝ^n)/∼, which tries to find the minimal total transport cost between [μ] and [ν]. Figure <ref> illustrates this idea.
We can see that the value of the ROT problem is invariant to translations. This is because the ROT problem is only associated with the equivalence classes of probability distributions.


Building upon the ROT model, we introduce a new family of Wasserstein distances to measure the minimal total transport cost between different equivalence classes of probability distributions. As mentioned above, the value of the ROT problem is invariant to any relative translations, thus, we name the corresponding Wasserstein distances as relative-translation invariant Wasserstein distances, denoted by RW_p, (p ≥ 1):
[Relative-translation invariant Wasserstein distances]
    
    RW_p(μ, ν) = ROT(μ, ν, p)^1/p.



Similar to the situation where W_p is a real metric on 𝒫_p(ℝ^n), we can obtain the following theorem.

    RW_p is a real metric on the quotient set 𝒫_p(ℝ^n)/∼.

The proof of Theorem <ref> is provided in Appendix A. The RW_p can also be regarded as a pseudo-metric defined on the 𝒫_p(ℝ^n), where the indiscernibility is no longer admitted.

It should be noted that we can not incorporate rotation in our equation <ref> since rotation will violate the metric properties and convexity.







 §.§ Quadratic ROT and Properties of the RW_2 Distance



We show three useful properties in the quadratic case of ROT and the resulted RW_2 distance: decomposability of the ROT optimization model (Theorem <ref>), translation-invariance of coupling solutions of the ROT problem (Corollary <ref>), and Pythagorean relationship of RW_2 and W_2 (Corollary <ref>). Decomposability will motivate our algorithm design for efficiently solving an ROT problem, translation-invariance provides RW_2 the robustness against distribution shift, and the Pythagorean relationship helps us better understand a distribution shift in the perspective of bias-variance.



    Continuing with previous notations, the two-stage optimization problem in quadratic ROT can be decomposed into two independent single-stage optimization problems: 
    
    ROT(μ, ν, 2) = 
            s ∈ℝ^nminP ∈Π(μ, ν)min H(P, s) = P ∈Π(μ, ν)min (E(P)) + s ∈ℝ^n min (V(s))

    where horizontal function H(P, s) is described by H(P, s) =  ∑_i=1^m_1∑_j=1^m_2x_i + s - y_j _2^2 P_ij,
    the overall function E(P) is described by,
    E(P) = ∑_i=1^m_1∑_j=1^m_2x_i- y_j _2^2 P_ij,
    and vertical function V(s) is described by, 
    V(s) = s_2^2 + 2 s · (μ̅ - ν̅).
    Moreover, the variables P_ij are fully determined by P ∈Π(μ, ν)min E(P), which is the classical optimal transport problem between μ and ν; and the variable s is fully determined by s ∈ℝ^n min V(s), which is a minimization problem of a quadratic function in ℝ^n and the minimum is achieved when s = ν̅ - μ̅.

Function E(P), V(P,s) and V(s) are illustrated in Figure <ref>. Due to the page limit, the proof of Theorem <ref> is provided in Appendix A.


Theorem <ref> is the core idea for our algorithm design in Section 4. It indicates that the coupling solutions P_ij to the OT problem are always the same as its ROT version, and verse versa. In other words, 

    The coupling solutions to the quadratic ROT problem are invariant to any translation of distributions.


Corollary <ref> not only guarantees the robustness of RW_2 against translational shifts but also suggests that the coupling solution of a ROT problem (including the classical OT problem) can be calculated by finding a more “proper” cost matrix of a translated probability distribution pair. This helps us improve the numerical stability and reduce the time complexity in many practical conditions. We provide a detailed analysis in Section <ref> and demonstrate it in Section <ref>.




Let s be the minimizer ν̅ - μ̅, it follows that,
    
    W_2^2(μ, ν)= μ̅ - ν̅_2^2 + RW_2^2(μ,ν).






Corollary <ref> indicates that there exists a Pythagorean relationship among three types of distances, W_2, RW_2, and L_2, as illustrated in Figure <ref>.


Corollary <ref> provides a refinement to understand a distribution shift (measured by W_2) in the perspective of bias and variance decomposition. The L_2 Euclidean distance between the expectations of two distributions corresponds to the “bias” between two distributions, and the value of RW_2 corresponds to the difference of “variances” or the “shapes” of two distributions. When the source distribution is a Dirac distribution, the situation degenerates to a classical bias-variance decomposition.










§ RW_2 SINKHORN ALGORITHM

Based on the first two properties of RW_2, we design a variant of the Sinkhorn algorithm, named RW_2 Sinkhorn, to compute RW_2 under translational shift. Furthermore, we show the new algorithm can be used as a technique to compute W_2 distance, and it provides enhanced numerical stability and a reduction of time complexity when the difference between the means of two distributions is large.
r0.5
0.5















 §.§ RW_2 Algorithm Design

Based on Theorem <ref> and Corollary <ref>, we propose the RW_2 Sinkhorn algorithm for computing RW_2 distance and coupling solution P_ij, which is described in Algorithm <ref>. The key idea of this algorithm involves precomputing the difference between the means of two distributions, as shown in Line 3. Subsequently, it addresses a specific instance of the optimal transport problem where the means of the two distributions are identical by a regular Sinkhorn algorithm. It is important to note that alternative algorithms, such as the network-simplex algorithm or the auction algorithm <cit.>, can also be employed to complete the specific instance procedure.



 §.§ RW_2 Technique for W_2 Computation

With the observation of Corollary <ref>, we can also propose a new approach to compute the W_2 distance from the right side of the Equation (<ref>). When ν̅ - μ̅_2 is large enough, this new RW_2 based technique has advantages than the regular Sinkhorn in terms of numerical stability and time complexity. We provide an analysis of this new approach in the rest of this section. In addition, experiment in Section 5.1 validates the analysis of our proposed RW_2 Sinkhorn algorithm for computing W_2 distance.



 §.§ Numerical Stability


The division by zero is a known numerical issue of the Sinkhorn algorithm, as shown in Equation (<ref>), and infinitesimal value often occurs in the exponential process of the (negative) cost matrix, that is, K← e^-C/λ. The results of the Corollary <ref> suggest that it is possible to switch to another “mutually translated" cost matrix under a relative translation s to increase the numerical stability while preserving the same coupling solutions.


To measure the numerical stability of the matrix, we introduce g(K), defined by the product of all entries K_ij. Since all entries K_ij are in the exponential format and lie in the range of (0,1], g(K) is also in the range of (0,1]. As g(K) increases, most entries K_ij deviate from zero, which means numerical computation will be more stable. Because of the following transformation, 

    g(K)=∏_i=1^m_1∏_j=1^m_2 K_ij = ∏_i=1^m_1∏_j=1^m_2exp(-C_ij/λ) = exp(-∑_i=1^m_1∑_j=1^m_2x_i +s - y_j_2^2/λ),

one can choose the relative translation s=y̅ - x̅, so that g(K) can achieve its maximum. When we assume the probability mass on each data is equal, y̅ - x̅ is the same as  ν̅ - μ̅.



 §.§ Complexity Analysis

<cit.> provided a proof for the time complexity of the optimal transport model by Sinkhorn algorithm with τ approximation, which is O(m^2C_∞^3(log m) τ^-3), where C_∞ = max_ij C_ij, with assuming m = m_1 = m_2 for the sake of simplicity. When the translated cost matrix has a smaller infinity norm C_∞, we can see that the time complexity will be reduced. When two distributions are sub-Gaussian and ν̅ - μ̅_2 is large enough, we can have a more rigorous analysis for the complexity analysis, which is provided in Appendix C. Our results from the first experiment in Section <ref> also empirically validate this idea.





§ EXPERIMENTS

To comprehensively evaluate and validate our proposed method, we conducted three experiments: numerical validation, digit recognition, and similar weather pattern detection. The three experiments are motivated differently: the first experiment is used to validate the numerical stability and computational time of the RW_2 Sinkhorn algorithm; the second experiment serves as validation of the robustness of RW_2 distance to the distribution translation in a real-world application; and the third experiment is to show the RW_2 can scale up to large-scale datasets to identify the similar weather patterns. All tests are executed on a workstation with a 2.60 GHz Intel Core i7 processor and 16GB of RAM using a single computation thread.



 §.§ Numerical Validation

r0.35



    < g r a p h i c s >



Schematic illustration of the first experiment. Assume that the distributions μ and ν are the same type of distribution. We compare the performance of Algorithm 1 and the classical Sinkhorn by translating the distribution μ along vector s, i.e., s =ν̅ - μ̅.



First we demonstrate the benefits of applying the RW_2 Sinkhorn algorithm for computing W_2 distance on specially designed examples. Let us consider two sets of data points randomly sampled from two identical distributions μ and ν, each with 1,000 samples. To compare the performance of Algorithm 1 and the classical Sinkhorn, we translate the distribution μ with a vector s, where the length of translation s ranges within [0, 3]. The idea is illustrated in Figure <ref>.





  
Experimental setting
We compare the two versions of Sinkhorn algorithms in W_2 error and running time. We repeat each experiment and data sampling 10 times. We test a pair of Gaussian distributions in ℝ (Figure 3(a) and 3(d)), a pair of Gaussians in ℝ^10 (Figure 3(b) and 3(e)), and a pair of uniform distributions in ℝ^10 (Figure 3(c) and 3(e)). We set λ = 0.1 and ϵ = 0.01 for both algorithms.





  
Experiment results
Figure <ref> shows that when the length of the translation increases, the proposed RW_2 Sinkhorn algorithm significantly outperforms the regular Sinkhorn algorithm in both the running time and computational errors. We also test the performance of the RW_2 Sinkhorn algorithm on other distributions. See more results in Appendix B.





 §.§ Digit Recognition

We test the robustness of the RW_2 distance on the MNIST dataset, perturbed by random translations. Each digital image is a 28 x 28 grid and can be converted into a discrete probability distribution by normalizing the intensity on each pixel by the sum of all intensities. We consider a subset of N data points randomly sampled from the MNIST dataset, where N = {100, 1000}. We embed each image in a larger 84 x 84 grid and translate each image with a random vector s (varying lengths and directions). We test performance on the length of s ranging within {0, 4, 8, 12, 16, 20, 24, 28} pixels.





  
Experimental setting
We report the mean and standard deviation of classification error using a 1/4 test ratio scheme, repeated 10 times. We use the nearest neighbor search as a classifier with RW_2 distance, as well as other distances (L_1, L_2, W_1, W_2) as the baselines. We set the λ = 0.1 and ϵ = 0.1 for the RW_2 Sinkhorn algorithm.





  
Experiment results
Figure <ref> shows that RW_2 distance significantly outperforms the other distances when the translation s is large in both sample sizes N = 100 and N = 1,000.





 §.§ Similar Thunderstorm Pattern Detection

We apply RW_2 distance on the real-world thunderstorm dataset, to show that RW_2 can be used for identifying similar weather patterns and focus more on shape similarity compared with W_2 distance. Our data are radar images from MULTI-RADAR/MULTI-SENSOR SYSTEM (MRMS) <cit.> with a 150 km × 150 km rectangular area centered at the Dallas Fort Worth International Airport (DFW), where each pixel represents a 3 km × 3 km area. The data is assimilated every 10 minutes tracking a long time range of history from 2014 to 2022, with 205,848 images in total. Vertically Integrated Liquid Density (VIL density) and reflectivity are two common measurements for assessing thunderstorm intensity, with threshold values of 3 kg· m^-3 and 35 dBZ, respectively <cit.>. For the sake of simplicity, we will focus on reflectivity as our main indicator for thunderstorm intensity. 

We consider two types of thunderstorm events: thunderstorm snapshots and thunderstorm sequences. More comparison results can be found in the Appendix B.





  
Experimental setting We apply the RW_2 Sinkhorn algorithm to compute both the W_2 distance and the RW_2 distance to identify the most similar thunderstorm event given a reference event, with λ = 0.1 and ϵ = 0.01. For the thunderstorm sequence identification, we consider a one-hour thunderstorm event as our reference, consisting of 6 images in time order.





  
Thunderstorm snapshot results
Given the same reference thunderstorm snapshot, Figure <ref> shows that the similar events identified by RW_2 focus more on shape similarity compared with W_2.







  
Thunderstorm sequence results
Given a reference thunderstorm sequence, Figure <ref> shows that the similar sequences identified by RW_2 focus more on shape similarity compared with W_2.








§ CONCLUSIONS

In this paper, we introduce a new family of distances, relative-translation invariant Wasserstein (RW_p) distances, for measuring the pattern similarity between two probability distributions (and their data supports). Generalizing from the classical optimal transport model, we show that the proposed RW_p distances are real distance metrics defined on the quotient set 𝒫_p(ℝ^n)/∼ and invariant to the translations of distributions. When p=2, this distance enjoys more useful properties, including decomposability of the reformulated optimal transport model, translation-invariance of coupling solutions and RW_2, and Pythagorean relationship of RW_2 and W_2 distances. Based on these properties, we show a distribution shift, measured by W_2 distance, which can be explained from the perspective of bias-variance. In addition, we propose a variant of the Sinkhorn algorithm, named RW_2 Sinkhorn algorithm, for efficiently calculating RW_2 distance, coupling solutions, as well as W_2 distance. We provide the analysis of numerical stability and time complexity for the proposed algorithm. Finally, we validate the RW_2 distance metric and the algorithm performance with illustrative and real-world experiments. The experimental results report that our proposed algorithm significantly improves the computational efficiency of Sinkhorn in practical applications with large translations, and the RW_2 distance is robust to distribution translations compared with baselines.






plainnat
 







§ MAIN PROOFS






 §.§ Theorem <ref>




It is clear to verify that when s_p ≥ 2 max_ijx_i-y_j_p, for any i,j, (1≤ i ≤ m_1, 1 ≤ j ≤ m_2), it follows that x_i + s-y_j_p ≥s _p - x_i - y_j_p ≥ 2max_ijx_i-y_j_p - x_i - y_j_p ≥x_i - y_j_p. In other words, when s_p ≥ 2max_ijx_i-y_j_p, the relative distance between each pair of support x_i and y_j are always greater than or equal to the non-translated distance, which implies the total transport cost for the translated case will also greater than or equal to the cost for the non-translated distance. Since we are trying to find the minimal value, we can only focus on the compact set {s∈ℝ^n | s_p ≤ 2 max_ijx_i-y_j_p}.




 §.§ Theorem <ref>


With the previous notations, firstly, we will show that the translation relation ∼ is an equivalence relation on set 𝒫_p(ℝ^n).

Equivalence relation requires reflexivity, symmetry, and transitivity, and the following observations show translation relation is indeed an equivalence relation.

    
  * Reflexivity, (x∼ x).
    
    For any distribution μ∈𝒫_p(ℝ^n), it can translate to itself with zero vector.
    
    
  * Symmetry, (x∼ y  y ∼ x).
    
    For any distribution μ, ν∈𝒫_p(ℝ^n), if μ can be translated to ν, then ν can also be translated to μ.

    
  * Transitivity, (x∼ y  and  y ∼ z  x ∼ z).
    
    For any distribution μ, ν, η∈𝒫_p(ℝ^n), if μ can be translated to ν, and ν can be translated to η, then μ can also be translated to η.


Based on the property of equivalence relation, it is clear that set 𝒫_p(ℝ^n)/∼ is well-defined. Let [μ] be an element in set M/∼, where μ is a representative of  [μ], i.e. [μ] is the set of distributions that can be mutually translated from μ. Noticing that W_p(·, ·) is a real distance metric on 𝒫_p(ℝ^n) <cit.>, it implies that W_p(·, ·) satisfied with identity, positivity, symmetry, and the triangle inequality. Based on  W_p(·, ·), we show RW_p(·, ·) satisfies with identity, positivity, symmetry, and triangle inequality with respect to elements in 𝒫_p(ℝ^n)/∼.

For any μ, ν, η∈𝒫_p(ℝ^n)/∼,

    
  * Identity,
        
    RW_p([μ], [μ])  =  μ∈ [μ],μ∈ [μ]min. [W_p(μ, μ)] = 0.

        
    
  * Positivity,
        
    RW_p([μ], [ν])  =  μ∈ [μ],ν∈ [ν]min. [W_p(μ, ν)] ≥ 0.


    
    
    
    
    
    
    
    
    
    
  * Symmetry,
    
    RW_p(μ, ν) = μ∈ [μ], ν∈ [ν]min.[W_p(μ, ν)] =ν∈ [ν], μ∈ [μ]min. [W_p(ν, μ )]= RW_p(ν, μ).

    
  * Triangle inequality,
    
    RW_p(μ, ν)  
    
                =   μ∈ [μ], ν∈ [ν]min.[W_p(μ, ν)] 
    ≤   η, η' ∈ [η]min.μ∈ [μ], ν∈ [ν]min.[W_p(μ, η
    ) + W_p(η, η') + W_p(η', ν)] 
    
                =   μ∈ [μ], ν∈ [ν], η, η' ∈ [η]min.[W_p(μ, η
    ) + 0 + W_p(η', ν)] 
    
                =    μ∈ [μ], η∈ [η]min.[W_p(μ, η)]  + ν∈ [ν], η' ∈ [η]min.[  W_p(η', ν)] 
    
                =    μ∈ [μ], η∈ [η]min.[W_p(μ, η)]  + ν∈ [ν], η∈ [η]min.[  W_p(η, ν)] 
    
                =     RW_p(μ, η) + RW_p(η, ν).






 §.§ Theorem <ref>


With the previous notations, firstly, we show the two-stage optimization problem, s ∈ℝ^nminP ∈Π(μ, ν)min H(P, s), can be decomposed into two independent one-stage optimization problems, P ∈Π(μ, ν)min (E(P)) and s ∈ℝ^n min (V(s)).

For the objective function H(P, s), we expand it with respect to s,

    H(P, s) 
    
        =    ∑_i=1^m_1∑_j=1^m_2x_i + s -y_j _2^2 P_ij
    
        =    ∑_i=1^m_1∑_j=1^m_2(x_i-y_j_2^2+ s_2^2 + 2 s · (x_i-y_j)) P_ij
    
        =    ∑_i=1^m_1∑_j=1^m_2x_i-y_j_2^2 P_ij
         +∑_i=1^m_1∑_j=1^m_2s_2^2 P_ij + 2∑_i=1^m_1∑_j=1^m_2 s · (x_i-y_j)  P_ij.


We can rewrite the second and the third terms in Equation (<ref>) under the condition P ∈Π(μ, ν), which implies that,

    ∑_i=1^m_1∑_j=1^m_2 P_ij =1, ∑_j=1^m_2 P_ij = a_i, ∑_i=1^m_1 P_ij = b_j, 1≤ i ≤ m_1, 1 ≤ j ≤ m_2.


For the second term, it follows that

    ∑_i=1^m_1∑_j=1^m_2s_2^2 P_ij =  s_2^2 · (∑_i=1^m_1∑_j=1^m_2 P_ij) = s_2^2 · 1= s_2^2.


For the third term, it follows that

    2∑_i=1^m_1∑_j=1^m_2 s · (x_i-y_j)  P_ij
       = 2 s ·∑_i=1^m_1∑_j=1^m_2 (x_i-y_j)  P_ij
       = 2 s · (∑_i=1^m_1∑_j=1^m_2 x_i · P_ij - ∑_i=1^m_1∑_j=1^m_2 y_j ·  P_ij) 
       = 2 s · (∑_i=1^m_1 x_i · (∑_j=1^m_2 P_ij) - ∑_j=1^m_2 y_j · (∑_i=1^m_1P_ij)) 
        = 2 s · (∑_i=1^m_1 x_i · a_i - ∑_j=1^m_2 y_j · b_j) 
        = 2 s · (μ̅ - ν̅).


Thus, we have the following transformation,

    s ∈ℝ^nminP ∈Π(μ, ν)min H(P, s) 
        =  s ∈ℝ^nminP ∈Π(μ, ν)min (∑_i=1^m_1∑_j=1^m_2x_i-y_j_2^2 P_ij
         +∑_i=1^m_1∑_j=1^m_2s_2^2 P_ij + 2∑_i=1^m_1∑_j=1^m_2 s · (x_i-y_j)  P_ij)  
        =  s ∈ℝ^nminP ∈Π(μ, ν)min∑_i=1^m_1∑_j=1^m_2x_i-y_j_2^2 P_ij
         + s ∈ℝ^nminP ∈Π(μ, ν)min (∑_i=1^m_1∑_j=1^m_2s_2^2 P_ij + 2∑_i=1^m_1∑_j=1^m_2 s · (x_i-y_j)  P_ij)  
        =  s ∈ℝ^nminP ∈Π(μ, ν)min∑_i=1^m_1∑_j=1^m_2x_i-y_j_2^2 P_ij
         + s ∈ℝ^nmin  ( s_2^2+ 2 s · (μ̅ - ν̅))  
        =  P ∈Π(μ, ν)min∑_i=1^m_1∑_j=1^m_2x_i-y_j_2^2 P_ij
         + s ∈ℝ^nmin  ( s_2^2 + 2 s · (μ̅ - ν̅))  
        =  P ∈Π(μ, ν)min (E(P)) + s ∈ℝ^n min (V(s))


Since V(s) is a high-dimensional quadratic function of variable s, it is easy to follow that the minimum is achieved when s = ν̅ - μ̅.






















§ ADDITIONAL EXPERIMENT RESULTS




 §.§ Additional experiment results for Section 5.1 - Numerical Validation






 §.§ Additional experiment results for Section 5.3 - Similar Thunderstorm Pattern Detection













§ COMPLEXITY ANALYSIS FOR ALGORITHM 1 UNDER SUB-GAUSSIAN DISTRIBUTIONS








This section is organized as follows. In section <ref>, we state and prove the theorem regarding the time complexity of Algorithm 1. We leave the definitions and theorems used in the proof to section <ref>.


 §.§ Theoretical Results of Time Complexity


Assuming the two distributions are sub-Gaussian and ν̅ - μ̅_2^2 is large enough, we can prove that with high probability, the translated cost matrix has a smaller norm C_∞, thus the time complexity can be reduced.

Let μ, ν be two high-dimensional sub-Gaussian distributions in ^n. (X_1, X_2, …, X_m_1), (Y_1, Y_2, …, Y_m_2) are i.i.d data sampled from μ and ν separately. Let μ̅= μ, ν̅= ν, X̅ = ∑_i=1^m_1 X_i / m_1, Y̅ = ∑_i=1^m_2 Y_i / m_2. Assume μ-μ̅_ψ_2 < ∞, ν - ν̅_ψ_2 < ∞. l = μ̅- ν̅_2. If 

    l    ≥ C√(n)[1 +  μ-μ̅_ψ_2 + ν-ν̅_ψ_2] 
            + C [√(log(4m_1/δ))·μ - μ̅_ψ_2 + √(log(4m_2/δ)))·ν - ν̅_ψ_2],

where C is an absolute constant,
then with probability at least 1-δ, we have

    max_i,jX_i - X̅ - Y_j + Y̅_2 ≤max_i,jX_i - Y_j_2.

 

    Sub-Gaussian distributions represent a broad class of distributions that encompass many common types, including multivariate normal distribution, multivariate symmetric Bernoulli, and uniform distribution on the sphere. Theorem <ref> demonstrates that when the translation is significant, the maximum absolute entry of the cost matrix C_∞ = max_ij |C_ij| tends to decrease. Consequently, our RW_2 method achieves better time complexity compared to W_2. This theoretical finding is consistent with our experimental results, as shown in Figure <ref>.


For i=1,2,…, m_1, X_i - μ̅ is a sub-Gaussian random vector. Using Theorem <ref> and taking a union bound over all the random vectors, we have for all X_i with probability at least 1- δ/4, the following inequality holds

    X_i - μ̅_2 ≤ c(√(n) + √(log(4m_1/δ)))·μ - μ̅_ψ_2.

Similarly, we have for all Y_j, with probability at least 1-δ, the following inequality holds

    Y_j - ν̅_2 ≤ c(√(n) + √(log(4m_2/δ)))·ν - ν̅_ψ_2.

Using Theorem <ref>, ∑_i=1^m (X_i - μ̅) is a sub-Gaussian random vector, with ∑_i=1^m_1 (X_i - μ̅)_ψ_2≤√(C ∑_i=1^m_1X_i - μ̅_ψ_2^2) Then using Theorem <ref>, with probability at least 1-δ/4, we have

    ∑_i=1^m_1 X_i - m_1μ̅_2 
           ≤ c(√(n) + √(log(1/δ)))·∑_i=1^m_1 (X_i - μ̅)_ψ_2
        = c(√(n) + √(log(1/δ)))·√(C ∑_i=1^m_1X_i - μ̅_ψ_2^2)
        = c'(√(n) + √(log(1/δ)))·√(m_1)μ-μ̅_ψ_2,

where c' is an absolute constant. Similarly, with probability at least 1-δ, we have

    ∑_j=1^m_2 Y_j - m_2ν̅_2 ≤ c'(√(n) + √(log(1/δ)))·√(m_2)ν-ν̅_ψ_2,

where c' is an absolute constant. In the following proof, we consider the union bound of all the high-probability events above, such that (<ref>), (<ref>), (<ref>) and (<ref>) hold. It occurs with probability at least 1-δ.

First, for max_i,jX_i - Y_j_2, we have

    max_i,jX_i - Y_j_2    ≥max_i,jμ̅-ν̅_2 - X_i - μ̅_2 - Y_j-ν̅_2
       
       ≥ l - [c(√(n) + √(log(4m_1/δ)))·μ - μ̅_ψ_2] 
            - [c(√(n) + √(log(4m_2/δ)))·ν - ν̅_ψ_2]
       = l - 2c√(n) - c√(log(4m_1/δ)))·μ - μ̅_ψ_2
            - c√(log(4m_2/δ)))·ν - ν̅_ψ_2,

where the first inequality holds due to the triangle inequality. The second inequality holds due to (<ref>) and (<ref>).

For max_i,jX_i - Y_j - X̅_i + Y̅_j_2, we have

    max_i,jX_i - Y_j - X̅_i + Y̅_j_2 ≤max_i,jX_i-μ̅_2 + Y_j-ν̅_2
            + X̅_i - μ̅_2 + Y̅_j - ν̅_2
       ≤[c(√(n) + √(log(4m_1/δ)))·μ - μ̅_ψ_2] + [c(√(n) + √(log(4m_2/δ)))·ν - ν̅_ψ_2]
            + [c'√(n) + √(log(1/δ))/√(m_1)μ-μ̅_ψ_2] + [c'√(d) + √(log(1/δ))/√(m_2)ν-ν̅_ψ_2]
       ≤ C√(n)[1 +  μ - μ̅_ψ_2/√(m_1) + ν - ν̅_ψ_2/√(m_2)]
           + C [√(log(4m_1/δ))·μ - μ̅_ψ_2 + √(log(4m_2/δ)))·ν - ν̅_ψ_2],

where the first inequality holds due to (<ref>), (<ref>), (<ref>) and (<ref>).
Therefore, we have the following conclusion:
As long as 

    l    ≥ C√(n)[1 +  μ-μ̅_ψ_2 + ν-ν̅_ψ_2] 
            + C [√(log(4m_1/δ))·μ - μ̅_ψ_2 + √(log(4m_2/δ)))·ν - ν̅_ψ_2],

where C is an absolute constant, we can conclude that

    max_i,jX_i - Y_j - X̅_i + Y̅_j_2 ≤max_i,jX_i - Y_j_2.

This completes the proof of Theorem <ref>.



 §.§ High Dimensional Probability Basics

In this section, we introduce some basic knowledge we will use in the proof of Theorem.
The results mostly come from <cit.>.


We first introduce a broad and widely used distribution class.
[Sub-Gaussian]
    A random variable X that
satisfies one of the following equivalent properties is called a subgaussian random variable. 

    
  (a) There exists K_1 > 0 such that the tails of X satisfy


    {|X| ≥ t}≤ 2 exp(-t^2/K^2_1)  for all  t ≥ 0.


  (b) There exists K_2 > 0 such that the moments of X satisfy


    X_L^p = ( |X|^p)^1/p≤ K_2 √(p) for all  p ≥ 1.


  (c)
 There exists K_3 > 0 such that the moment-generating function (MGF) of X^2 satisfies
 
    exp(λ^2X^2) ≤exp(K^2_3λ^2)  for all λ such that  |λ| ≤1/K_3.


  (d) There exists K_4 > 0 such that the MGF of X^2
is bounded at some point,
namely,

    exp(X^2/K^2_4) ≤ 2.


  (e) Moreover, if X = 0, the following property is also equivalent.
There exists K_5 > 0 such that the MGF of X satisfies

    exp(λ X) ≤exp(K^2_5λ^2)  for all λ∈.


The parameters K_i > 0 appearing in
these properties differ from each other by at most an absolute constant factor.

The sub-gaussian norm of X, denoted X_ψ_2, is defined to be

    X_ψ_2 = inf{
    t > 0 : exp(X^2/t^2) ≤ 2}.



    A random vector X ∈^d
is sub-Gaussian
if for any vector ∈^d
the inner product ⟨ X, ⟩ is a sub-Gaussian random variable. And the
corresponding ψ_2 norm of X is defined as

    X_ψ_2 = sup__2=1⟨ X, ⟩_ψ_2.































































    




















    Let X_1,…,X_N ∈^d be
independent, mean zero, sub-Gaussian random vectors. Then
∑_i=1^N X_i is also a
sub-Gaussian random vector, and


    ∑_i=1
    ^N X_i^2_ψ_2≤ C ∑_i=1^N X_i^2_ψ_2.

where C is an absolute constant.

    For any vector ∈, _2 = 1, consider ⟨∑_i=1
^N X_i, ⟩. Using independence, we have for all λ,

    exp(λ∑_i=1
    ^N ⟨ X_i, ⟩)    = ∏_i=1^N exp(λ⟨ X_i, ⟩)
       ≤∏_i=1^N exp(C ⟨ X_i, ⟩_ψ_2^2 λ^2)
        = exp(Cλ^2 ∑_i=1^N ⟨ X_i, ⟩_ψ_2^2),

where C is an absolute constant and the first inequality holds due to property (e) of the sub-Gaussian variables. Taking supreme over , we prove that ∑_i=1^N X_i is also a
sub-Gaussian random vector. Moreover,

    ∑_i=1
    ^N X_i^2_ψ_2≤ C ∑_i=1^N X_i^2_ψ_2.

where C is an absolute constant.



    Let X ∈^d be a sub-Gaussian random vector. Then with probability at least 1-δ, 

    X_2 ≤ c(√(d) + √(log(1/δ)))·X_ψ_2.



    Let B_d be the d-dimensional unit ball, N
be a 1/2-covering of B_d in 2-norm with covering number = N(B_d, ·_2, 1/2). Therefore,

    ∀∈ B_d, ∃∈ N,  s.t.  - ≤ 1/2.

Using Lemma <ref>, we have

    N ≤  5^d.

Using the fact _2 = max__2 ≤ 1⟨, ⟩, we have

    X_2    = max_∈ B_d⟨, X ⟩
       ≤max_∈ N⟨, X⟩ + max_∈ (1/2) B_d⟨
    , X ⟩
        = max_∈ N⟨, X⟩ + 1/2max_∈  B_d⟨
    , X ⟩.

Therefore, we have

    X_2 ≤ 2 max_∈ N⟨, X⟩.

Then we can provide a high probability upper bound for the Euclidean norm of the random vector X by considering the probability (X_2 ≥ t).

    (X_2 ≥ t)    ≤(max_∈ N⟨, X⟩≥t/2)
       ≤(∃∈ N, 
    ⟨, X⟩≥t/2)
       ≤∑_∈ N(
    ⟨, X⟩≥t/2)
       ≤ N exp(- ct^2/X_ψ_2^2)
       ≤ 5^d exp(- ct^2/X_ψ_2^2),

where c is an absolute constant. Here the first inequality holds due to (<ref>). The second inequality holds due to {max_∈ N⟨, X⟩≥ 
t/2}⊆{∃∈ N, 
⟨, X⟩≥ 
t/2}. The third inequality holds due to the union bound. The fourth inequality holds due to the definition of the sub-Gaussian vector and the property 
(a) of a sub-Gaussian variable. The last inequality holds due to (<ref>).

Finally, let t = √([d log 5 + log(1/δ)]/c)·X_ψ_2. We have with probability at least 1-δ,

    X_2 ≥ t.

Finally, using √(a+b)≤√(a) + √(b), we complete the proof of Theorem <ref>.

[ϵ-covering]
Let (V, ·) be a normed space, and Θ⊂ V. V_1,…, V_N is an
ϵ-covering of Θ if Θ⊆∪_i=1^N V_i, or equivalently, ∀θ∈Θ, ∃ i such that θ - V_i≤ϵ.

[Covering number]
The covering number is defined by

    N(Θ, ·, ϵ) := min{n : ∃ϵ-covering over Θ of size  n}.



    Let B_d be the d-dimensional Euclidean unit ball. Consider N(B_d, ·_2, ϵ). When ϵ≥ 1, N(B_d, ·_2, ϵ) =
1. When ϵ < 1, we have

    (1/ϵ)^d≤ N(B_d, ·_2, ϵ) ≤(1+2/ϵ)^d.
































 




























































































































    

























    
















































§ PROPOSED OPTIMAL TRANSPORT FRAMEWORK WITH RW DISTANCE FOR SPATIAL-TEMPORAL DATASETS

After the above definitions and theorems, we have introduced a new metric called RW distance to define the similarity between two distributions. We will introduce how to use RW distance to identify similar spatial-temporal data patterns, by treating them as discrete probability distributions.

Generally speaking, each spatial-temporal data is described as “a series or sequence of vectors”. Here, we use a set of spatial-temporal vectors to represent spatial-temporal data d. i.e., 

    d = { (t_i, x_i)| 1 ≤ i ≤ m},

where t_i is a scale value representing the time information of this vector, and x_i is a multi-dimensional vector that stands for the space information, and assumes there are m spatial-temporal vectors in the data d.

By treating each spatial-temporal vector as a support point, we can treat the spatial-temporal data d as a probability distribution, where the probability mass of each support point can be distributed uniformly or customized based on different objectives. Here, we assume the probability mass of spatial-temporal vectors is distributed uniformly.

Since time is independent of space information, it is unwise to directly treat raw data as distributions. Thus, we introduce a time weight parameter, w, to customize different demands of spatial and temporal similarity. Formally, the time axis dilation operation is defined as follows:


  
Definition 5 The time axis dilation operation ω(·).

Assume that there is spatial-temporal data d and a time weight w, i.e., 

    d = {(t_1, x_1), ..., (t_m, x_m) }, w ∈ [0,+∞).

Then the time axis dilation operation:

    ω(d): = {(w· t_1, x_1), ..., (w· t_m, x_m) }.


The advantage of this time weight is that we can adjust the spatial and temporal similarity by changing the value of w. Specifically, 


    
  * when w goes to 0, which leads to w · t_i (1≤ i ≤ m) also going to 0, the cost of matching along the time axis will be cheaper. Thus, the weighted data will reflect more of the spatial information;
    
  * when w goes to infinity, which leads to w · t_i (1≤ i ≤ m) also going to infinity, the cost of matching along the time axis will be more expensive. Thus, the weighted data will reflect more of the time series information.


It is important to consider whether different orders of move-to-origin shift and time axis dilation to distribution will lead to different results. In fact, the following theorem shows that the two operations are independent and the order does not matter.



  
Theorem 6 Given any data d, the time axis dilation operation ω(d) and move-to-origin shift τ (d) are commutative, i.e.,

    ω∘τ (d) =τ∘ω  (d).


The proof of Theorem 6 can be found in the supplementary material.

The independence of the above two operations provides a flexible framework to unify different definitions of similarity, as shown in Table 1.












the shift in terms of their shift-invariant difference with toleration of the shifts of distributions.










The histogram method <cit.> is another common method to use for measuring the similarity of high-dimensional temporal sequences. The basic idea is to extract a few essential features and use those features to build feature vectors for the temporal sequences. Finally, we compute the distance between feature vectors to compute the similarities. The advantage is that it is a flexible framework for choosing all kinds of global or local features, such as motion, time length, level of clustering, shape, etc. However, it relies on human knowledge, and often difficult to find the proper features from the dataset when it is lacking understanding and human experience. In a real-world application, it would be too complicated to tune this heavy multi-parameter model properly to make it work.

Dynamic time warping is a traditional method of measuring the similarity between two temporal sequences in time series analysis, which may allow different time lengths. The main idea of this method is that it allows to have one-to-many and many-to-one matches for different time lengths. The goal of this method is to  minimize the sum of the distances of mapping between the two sequences. A common strategy to find the best matches is via dynamic programming, which has the time complexity of O(nm) <cit.>, where n and m are the lengths of two input sequences, separately. This method has many different variants in real applications, including SparseDTW <cit.>, FastDTW <cit.>, etc. However, because of the curse of dimensionality, it is not an easy way to provide a one-to-many or many-to-one optimal matching in high-dimensional space. This method is usually applied in low-dimensional cases and is hard to extend to a high-dimension space. In addition, the distance between two temporal sequences computed by DTW is not a true metric and is sensitive to the shift or disturbance of absolute coordinates of two temporal sequences.


Due to the space limit, we leave the remaining results of k = 1 and k = 5 to the supplementary material. Table 2 shows the results when k =3.



We choose the Modified National Institute of Standards and Technology database (MNIST) <cit.> as the dataset for this task. MNIST is a large dataset of handwritten digits containing 60,000 images of 10 digits, where each image is a greyscale image with 28x28 pixels. The k-nearest neighbors algorithm (k-NN) is used as the basic framework with different distance metrics, including Euclidean distance (L_2), the Wasserstein distances (W_1 and W_2) <cit.>, and the shift-invariant Wasserstein distance (RW_2). The difference between W_2 and RW_2 is shown in Figure 4.


We collect N images randomly from the MNIST dataset and then randomly split the collected images into train/test set by a ratio of 3:1. Using the k-NN framework as a classifier, (k= 1, 3, 5), where the similarity measures include L_2, W_1, W_2, RW_2, we can compute the accuracy of collected data on the test set. Finally, we repeat the above procedures many times to compute the average accuracy and deviations.


    
  * Purely spatial similarity. The results are shown in Figure <ref> and Figure <ref>.
    
  * Spatial-temporal similarity. Assume that the time weight is 0.1. The results are shown in Figures <ref>, <ref> <ref>.












For the airline operational consideration, W_2 and RW may satisfy different demands. Since most of the flights are scheduled routinely at the same time for each day and W_2 will be a good option since it sticks to the specific time and locations. For the purposes of clustering thunderstorms in a local area, RW_2 might be a better metric, since it can get rid of the impact of coordinate shifts.













    Note that the existence of RW_p distances is inherited from the underlying geometry structure of Euclidean space. (If we change to another geometry structure, such as a spherical surface, those properties may no longer be admitted.)




