Paper deep dive
Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics
Aratrika Mustafi, Soumya Mukherjee
Intelligence
Status: succeeded | Model: anthropic/claude-sonnet-4.6 | Prompt: intel-v1 | Confidence: 96%
Last extracted: 3/24/2026, 3:55:20 AM
Summary
This paper proposes a dense associative memory framework for empirical measures (weighted point clouds), where stored patterns and queries are finitely supported probability measures. Retrieval is defined by minimizing a Hopfield-style log-sum-exp energy built from the debiased Sinkhorn divergence. The retrieval dynamics are derived as a Spherical Hellinger Kantorovich (SHK) gradient flow that updates both support locations and weights. The paper proves basin invariance, geometric convergence to local minimizers, exponential storage capacity under a random pattern model, and demonstrates robust recovery from perturbed queries versus a Euclidean Hopfield-type baseline.
Entities (32)
Relation Signals (28)
Aratrika Mustafi → authored → Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics
confidence 99% · Aratrika Mustafi Soumya Mukherjee listed as authors of the paper
Soumya Mukherjee → authored → Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics
confidence 99% · Aratrika Mustafi Soumya Mukherjee listed as authors of the paper
Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics → cites → Hopfield (1982)
confidence 99% · This viewpoint goes back to Hopfield's classical construction... (Hopfield (1982))
Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics → cites → Feydy et al. (2019)
confidence 99% · For more details on the Sinkhorn divergence, we refer our readers to Feydy et al. (2019)
Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics → cites → Tankala and Balasubramanian (2026)
confidence 99% · Tankala and Balasubramanian (2026) develop a distributional dense associative memory
Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics → cites → Ramsauer et al. (2020)
confidence 99% · some modern Hopfield style models (Krotov and Hopfield (2016), Ramsauer et al. (2020))
Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics → cites →
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We propose a dense associative memory for empirical measures (weighted point clouds). Stored patterns and queries are finitely supported probability measures, and retrieval is defined by minimizing a Hopfield-style log-sum-exp energy built from the debiased Sinkhorn divergence. We derive retrieval dynamics as a spherical Hellinger Kantorovich (SHK) gradient flow, which updates both support locations and weights. Discretizing the flow yields a deterministic algorithm that uses Sinkhorn potentials to compute barycentric transport steps and a multiplicative simplex reweighting. Under local separation and PL-type conditions we prove basin invariance, geometric convergence to a local minimizer, and a bound showing the minimizer remains close to the corresponding stored pattern. Under a random pattern model, we further show that these Sinkhorn basins are disjoint with high probability, implying exponential capacity in the ambient dimension. Experiments on synthetic Gaussian point-cloud memories demonstrate robust recovery from perturbed queries versus a Euclidean Hopfield-type baseline.
Tags
Links
- Source: https://arxiv.org/abs/2603.20656v1
- Canonical: https://arxiv.org/abs/2603.20656v1
Trouble viewing inline? Open PDF directly →
Full Text
153,268 characters extracted from source content.
Expand or collapse full text
Sinkhorn Based Associative Memory Retrieval Using Spherical Hellinger Kantorovich Dynamics Aratrika Mustafi Soumya Mukherjee Abstract We propose a dense associative memory for empirical measures (weighted point clouds). Stored patterns and queries are finitely supported probability measures, and retrieval is defined by minimizing a Hopfield-style log-sum-exp energy built from the debiased Sinkhorn divergence. We derive retrieval dynamics as a spherical Hellinger Kantorovich (SHK) gradient flow, which updates both support locations and weights. Discretizing the flow yields a deterministic algorithm that uses Sinkhorn potentials to compute barycentric transport steps and a multiplicative simplex reweighting. Under local separation and PL-type conditions we prove basin invariance, geometric convergence to a local minimizer, and a bound showing the minimizer remains close to the corresponding stored pattern. Under a random pattern model, we further show that these Sinkhorn basins are disjoint with high probability, implying exponential capacity in the ambient dimension. Experiments on synthetic Gaussian point-cloud memories demonstrate robust recovery from perturbed queries versus a Euclidean Hopfield-type baseline. 1 Introduction Associative memories are a class of models designed to store and retrieve information through patterns of association instead of recalling data by specifying its location. These models frame retrieval as a dynamical process : given a partial or corrupted query, the system evolves toward an attractor that represents a stored pattern. The query evolves through a dynamical process that decreases an energy whose local minima encode the memories. This viewpoint goes back to Hopfield’s classical construction, where one designs recurrent interactions that induce a Lyapunov (energy) function whose local minima correspond to memories, yielding content-addressable recall from partial or corrupted cues (Hopfield (1982)). More recently, some modern Hopfield style models (Krotov and Hopfield (2016),Ramsauer et al. (2020)) use the soft-min/ log-sum-exp (LSE) energy to create high-capacity retrieval dynamics and improve efficiency. In the continuous state setting, these retrieval updates can be related to the attention mechanism used in transformer architectures, suggesting a broad role for energy-based retrieval as a reusable computational primitive inside learned systems. A key limitation of much of the associative-memory literature is that patterns are typically vectors. A natural next step is to move beyond vector-valued patterns and treat probability measures as the objects being stored and retrieved. This shift is motivated by the growing role of distributional representations in modern learning, where uncertainty, multimodality, and population-level structure are often more naturally encoded by measures than by single points. In this direction, Tankala and Balasubramanian (2026) develop a distributional dense associative memory for Gaussian distributions in the 2-Wasserstein (Bures-Wasserstein) geometry, defining a log-sum-exp energy over stored distributions and deriving a retrieval dynamics that aggregates optimal transport maps in a Gibbs-weighted manner; their stationary points correspond to self-consistent Wasserstein barycenters. Many practical settings, however, are not well captured by a single Gaussian, instead distributions are often represented nonparametrically through weighted samples - for example, weighted point clouds, histograms, particle approximations of posteriors, or sets of learned feature vectors, suggesting the need for an associative-memory retrieval principle that operates directly on discrete measures. In this paper we pursue this analogous extension of dense-associative-memory retrieval when both stored patterns and queries are empirical measures (weighted point clouds), using a log-sum-exp energy built from the debiased Sinkhorn divergence. Sinkhorn divergences provide a computationally tractable way to compare discrete measures and yield deterministic barycentric projections from entropic optimal-transport couplings. A key issue is that transport only (Wasserstein) dynamics moves support locations while keeping the query weights fixed - an intrinsic limitation when the weights carry information such as saliency, mixture proportions, or coarse discretization effects. To evolve both support points and weights in a principled way while staying within probability measures, we therefore use a spherical Hellinger Kantorovich (SHK) dynamics [see Appendix D, Liero et al. (2018)], which couples Wasserstein (aka Kantorovich) transport with a spherical Hellinger-type reaction component. In the empirical setting, this leads to a fully deterministic retrieval operator that updates support points via Sinkhorn barycentric maps and updates weights via a multiplicative reweighting on the simplex. This combination is attractive because it aligns the retrieval dynamics with the underlying geometry of measures while allowing both where mass is located and how much mass is assigned to each particle, to adapt during recall. With this objective in mind, our goal is to construct a dense associative memory (DAM) capable of storing discrete measures and accurately recovering the corresponding distribution when provided with perturbed/noisy queries. Given a collection of finitely supported discrete distributions (stored patterns) Xii=1N\X_i\_i=1^N and a query finitely supported discrete distribution ξ, we define the LSE functional using the Sinkhorn divergence SεS_ (defined in 2), E(ξ)=−1βlog(∑i=1Nexp(−βSε(ξ,Xi))),E(ξ)\;=\;- 1β \! ( _i=1^N \! (-β\,S_ (ξ,X_i) ) ), (1) with inverse temperature β>0β>0. This is the functional we aim to minimize along the SHK gradient flow. To this end, we develop a practically implementable algorithm (SinkhornAlgo) to carry out the optimization. In parallel, we establish theoretical guarantees ensuring convergence and near-accurate retrieval, and we present empirical evidence to validate these findings. • Local convergence and stability guarantees in SHK geometry - Under standard local assumptions (margin separation, SHK smoothness, bounded gradients, and a local PL inequality), we prove descent and contraction properties for the softmin energy, quantify basin interference through separation margins, and establish geometric convergence of SHK gradient descent iterates to the unique local minimizer within a basin. We also bound the deviation between the basin minimizer and the corresponding stored pattern and provide a local stability bound showing stored patterns are approximate fixed points of one-step retrieval. • Separation results under a random pattern model. We introduce a sampling mechanism for generating general M-atom measures (random supports and weights) and prove that, with high probability, the induced Sinkhorn neighborhoods around stored patterns are pairwise disjoint. This yields an exponential-in-dimension scaling of the number of storable patterns under the stated conditions. • Empirical evidence against a Euclidean Hopfield-type baseline- On synthetic point-cloud patterns sampled from Gaussians, we compare our SinkhornSHK retrieval against a vectorized Euclidean Hopfield-style baseline and observe robust recovery from noisy queries, particularly in regimes where Euclidean similarity is ambiguous but distributional (OT/Sinkhorn) geometry remains discriminative. 2 Setting and notations Let Ω⊂ℝd ^d be open, bounded and convex with diameter D≔supx,y∈Ω‖x−y‖<∞D _x,y∈ \|x-y\|<∞ and boundary ∂Ω∂ , and we equip it with the usual Euclidean geometry induced by the ℓ2 _2 norm. For any ω∈Ωω∈ and any set S⊂ℝdS ^d, define dist(ω,S)=mins∈S‖ω−s‖dist(ω,S)= _s∈ S\|ω-s\|. For any x=(x1,…,xM)∈ΩMx=(x_1,…,x_M)∈ ^M, let d∂(x)≔min1≤m≤Mdist(xm,∂Ω)d_∂(x) _1≤ m≤ Mdist (x_m,∂ ) denote the distance of x from the boundary of Ω and sep(x)≔mini≠j‖xi−xj‖sep(x) _i≠ j \|x_i-x_j \| denote the minimum pairwise separation between any two components of x. Let (x,r)≔y∈Ω:‖y−x‖<rB(x,r) \y∈ :\|y-x\|<r \ and its closure ¯(x,r)≔y∈Ω:‖y−x‖≤r B(x,r) \y∈ :\|y-x\|≤ r \ be the open and closed ball of radius r around x∈Ωx∈ . We denote by ∥⋅∥1\|·\|_1 and ∥⋅∥∞\|·\|_∞ the ℓ1 ^1 and ℓ∞ ^∞ norms, respectively and ∥⋅∥\|·\| will be interpreted in the natural context. Fix any M∈ℕM , 1M>amin>0 1M>a_ >0 and Δmin>0 _ >0. Stored patterns are finitely supported discrete measures Xi=∑m=1Mbi,mδyi,m,bi,m>amin>0,∑m=1Mbi,m=1,X_i= _m=1^Mb_i,m\, _y_i,m,b_i,m>a_ >0, _m=1^Mb_i,m=1, with all the locations yi,my_i,m’s being separated from each other with minimum pairwise distance minm≠k‖yi,m−yi,k‖>Δmin _m≠ k\|y_i,m-y_i,k\|> _ , and the query is also a finitely supported discrete measures ξ=∑m=1Mamδxm,am>amin>0,∑m=1Mam=1.ξ= _m=1^Ma_m\, _x_m,a_m>a_ >0, _m=1^Ma_m=1. with all the locations xmx_m’s being pairwise separated with margin Δmin _ . For any M∈ℕM 1M>amin>0 1M>a_ >0 and Δmin>0 _ >0, define the space of bounded weights ΔM,amin∘≔a∈ℝM:am>amin,∑m=1Mam=1 _M,a_ \a ^M:a_m>a_ , _m=1^Ma_m=1 \ and the space of M pairwise separated locations/atoms LocM,Δmin(Ω)≔x=(x1,…,xM)∈ΩM:minm≠n‖xm−xn‖>Δmin∈ΔM,amin∘×LocM,Δmin(Ω)Loc_M, _ ( ) \x=(x_1,…,x_M)∈ ^M: _m≠ n\|x_m-x_n\|> _ \∈ _M,a_ ×Loc_M, _ ( ). We will often consider the unordered particle parameterization of a finitely supported discrete measure with M atoms (a,x)=((am)m=1M,(xm)m=1M)(a,x)= ( (a_m )_m=1^M, (x_m )_m=1^M ). Let (Ω)P( ) be the collection of probability measures on Ω and M,amin,Δmin(Ω)=ξ∈(Ω):ξ=∑l=1Malδxl for (a1,…,aM)∈ΔM,amin∘,(x1,…,xM)∈LocM,Δmin(Ω).P_M,a_ , _ ( )= \ξ ( ):ξ= _l=1^Ma_l _x_l for (a_1,…,a_M)∈ _M,a_ ,(x_1,…,x_M) _M, _ ( ) \. Throughout the paper, we use the quadratic cost c(x,y)=12‖x−y‖2c(x,y)= 12\|x-y\|^2 for defining the Wasserstein distance, spherical Hellinger Kantorovich distance and the Sinkhorn divergence. For ε>0 >0, define the entropic OT cost OTε(μ,ν)≔minπ∈Π(μ,ν)∫Ω×Ωc(x,y)dπ(x,y)+εKL(π∣μ⊗ν),OT_ (μ,ν) _π∈ (μ,ν) _ × c(x,y)\,dπ(x,y)\;+\; \,KL(π μ ν), where Π(μ,ν) (μ,ν) denotes couplings with marginals μ,νμ,ν. The debiased Sinkhorn divergence is defined as - Sε(μ,ν)≔OTε(μ,ν)−12OTε(μ,μ)−12OTε(ν,ν).S_ (μ,ν) _ (μ,ν)- 12OT_ (μ,μ)- 12OT_ (ν,ν). (2) A standard dual formulation of OTεOT_ uses Schrödinger (entropic OT) potentials (fμ,ν,gμ,ν)(f_μ,ν,g_μ,ν), defined up to an additive constant. We do not require a specific normalization; only gradients ∇f∇ f and differences of potentials will matter. One can define the entropic soft c-transform operator AεA_ via an expression of the form Aε(g,ν)(x)≔−εlog∫Ωexp(gμ,ν(y)−c(x,y)ε)ν(y)(defined up to an additive constant).A_ (g,ν)(x) - _ \! ( g_μ,ν(y)-c(x,y) )\,dν(y) (defined up to an additive constant). Then the optimal potentials (fμ,ν,gμ,ν)(f_μ,ν,g_μ,ν) (Schrödinger potentials) satisfy the Schrödinger system fμ,ν=Aε(gμ,ν,ν)μ-a.e.,gμ,ν=Aε(fμ,ν,μ)ν-a.e..f_μ,ν=A_ (g_μ,ν,ν) μ-a.e., g_μ,ν=A_ (f_μ,ν,μ) ν-a.e.. These potentials are unique up to adding a constant to fμ,νf_μ,ν and subtracting the same constant from gμ,νg_μ,ν (gauge invariance). This does not affect any gradient ∇fμ,ν(x)∇ f_μ,ν(x) or ∇gμ,ν(x)∇ g_μ,ν(x), which is what we ultimately use. The optimal entropic coupling has Gibbs form dπμ,νεd(μ⊗ν)(x,y)=exp(fμ,ν(x)+gμ,ν(y)−c(x,y)ε). dπ _μ,νd(μ ν)(x,y)= \! ( f_μ,ν(x)+g_μ,ν(y)-c(x,y) ). For more details on the Sinkhorn divergence, we refer our readers to Feydy et al. (2019); Hardion and Lavenant (2025). 3 Theoretical Guarantees In this section, we state all our theoretical results. The detailed proofs are available in Section G. Retrieval guarantees largely rely on ensuring sufficient separation in terms of Sinkhorn divergence among the stored patterns. We propose a sampling mechanism that allows the number of patterns N to be exponentially large in the dimension d and Theorem 1 establishes a high-probability separation guarantee for these N randomly generated stored patterns, ensuring that the associated Sinkhorn neighborhoods (basins) are pairwise disjoint. Theorem 1 (Exponential storage capacity and high probability separation of patterns). Fix d≥1,M≥2d≥ 1,M≥ 2, amin>0a_ >0, Δmin>0 _ >0 and let Ω⊂ℝd ^d be open, bounded and convex. Assume there exists c∈Ωc∈ and R>0R>0 such that the closed ball ¯(c,R)⊂Ω B(c,R)⊂ . Fix any σ∈(0,R/4)σ∈(0,R/4), such that σ,Δmin≔(z1,…,zM)∈(0,σ)M:minn≠m‖zn−zm‖>ΔminZ_σ, _ \ (z_1,…,z_M ) (0,σ)^M: _n≠ m \|z_n-z_m \|> _ \ is non-empty, and set R0≔R−2σR_0 R-2σ. Let γ,p∈(0,1)γ,p∈(0,1) and choose N≔⌊2pexp(γ24d)⌋.N 2p ( γ^24d ) . Let X1,…,XNX_1,…,X_N be generated by the sampling mechanism described as SampAlgo (see Sec F). Assume ε>0 >0 is chosen small enough so that εlogM<1−γ16R02, M< 1-γ16R_0^2, and define dmin≔2(1−γ)R0,r≔dmin232−εlogM,Δ≔dmin24.d_ 2(1-γ)R_0, r d_ ^232- M, d_ ^24. Then, with probability atleast 1−p1-p, the following pairwise separation of Sinkhorn neighbourhoods/basins holds true (i.e. Assumption (A1)): For each i, for every ξ∈Bi(r)=ν∈M,amin,Δmin(Ω):Sε(ν,Xi)≤rξ∈ B_i(r)= \ν _M,a_ , _ ( ):S_ (ν,X_i )≤ r \, and every j≠ij≠ i, Sε(ξ,Xj)−Sε(ξ,Xi)≥Δ.S_ (ξ,X_j )-S_ (ξ,X_i )≥ . In particular, Bi(r)∩Bi(j)=∅,∀i≠jB_i(r)∩ B_i(j)= , ∀ i≠ j i.e. the Sinkhorn neighborhoods of XiX_i’s are pairwise disjoint. The Sinkhorn margin separation property provides the foundation required for reliable retrieval, which ensures that if a query is close to a stored pattern, it is separated enough from other patterns to ensure accurate retrieval. In fact, we can prove retrieval guarantees in much greater generality than afforded by the particular sampling algorithm we propose. We can establish these results under some regularity assumptions regarding the energy functional E and the local Sinkhorn energies Fi(⋅)≔Sε(⋅,Xi)F_i(·) S_ (·,X_i), which are presented in Section E. Given the energy functional E along with a query ξ∈M,amin,Δmin(Ω)ξ _M,a_ , _ ( ) and a step-size η, we can define a gradient based evolution of E by equipping the space M,amin,Δmin(Ω)P_M,a_ , _ ( ) of finitely supported M-atom discrete probability measures with a choice of geometry. We choose the spherical Hellinger-Kantorovich (SHK) geometry and the gradient descent based one-step retrieval operator can be defined as Φη(ξ)=Retξ(−ηgradSHKE(ξ)) _η(ξ)=Ret_ξ(- _SHKE(ξ)) (3) where RetξRet_ξ is a retraction map that evolves the query ξ along the negative SHK gradient of E, given by gradSHKE(ξ)grad_SHKE(ξ), for a small step controlled by η and ensures that the resulting object is still a probability distribution (in fact, an element of M,amin,Δmin(Ω)P_M,a_ , _ ( )). We relegate all details to the Appendix. Theorem 2 (Geometric convergence in Sinkhorn divergence to the local minimizer and local basin invariance of gradient descent iterates). Let Assumptions (A2), (A3) and (A4) hold. Define the stored pattern margins wi≔min1≤m≤M(bi,m−amin)>0,di∂≔min1≤m≤Mdist(yi,m,∂Ω)>0w_i _1≤ m≤ M (b_i,m-a_ )>0,d_i^∂ _1≤ m≤ Mdist (y_i,m,∂ )>0 and si≔sep(yi)=minm≠n‖yi,m−yi,n‖2>Δmins_i (y_i )= _m≠ n \|y_i,m-y_i,n \|_2> _ . Let δi,τi>0 _i, _i>0 such that 0<δi<δ¯i≔mindi∂,si−Δmin2>00< _i< δ_i \d_i^∂, s_i- _ 2 \>0 and 0<τi<wi0< _i<w_i. Let r>0r>0 be such that 0<r<riloc(δi,τi)=minaminδi22−εlogM,τi(si−δi)24−εlogM0<r<r^loc_i ( _i, _i )= \ a_ _i^22- M, _i (s_i- _i )^24- M \. Let Ei∗(r)E_i^*(r) be the minimum value of E(ξ)E(ξ) in the local basin Bi(r)B_i(r) and Xi∗(r)X_i^*(r) be a minimizer of E in Bi(r)B_i(r) i.e. E(Xi∗(r))=Ei∗(r)≔infξ∈Bi(r)E(ξ)E(X_i^*(r))=E_i^*(r) _ξ∈ B_i(r)E(ξ). Finally, define ηret,i≔minλ22D2logminmbi,m−τiamin,12Dmindi∂−δi,si−2δi−Δmin _ret,i \ λ^22D^2 _mb_i,m- _ia_ , 12D \d_i^∂- _i,s_i-2 _i- _ \ \. Then for any step size 0<η<min1L,1μ,ηret,i0<η< \ 1L, 1μ, _ret,i\, the following hold true: 1. If ξ(k)∈Bi(r)ξ^(k)∈ B_i(r) for all k∈ℕ∪0k ∪\0\, then the sequence of SHK gradient iterates (ξ(k))k≥0 (ξ^(k) )_k≥ 0 satisfy the explicit geometric bound Sε(ξ(k),Xi∗(r))≤GSη2η(E(ξ(0))−Ei∗(r))1−(1−μη)12⋅(1−ημ)k2,S_ (ξ^(k),X_i^*(r) )≤ GS_η 2η (E (ξ^(0) )-E_i^*(r) )1- (1-μη ) 12·(1-ημ) k2, where Sη≤mineηD2λ2,1aminS_η≤ \e η D^2λ^2, 1 a_ \ and G=D1+D2λ2G=D 1+ D^2λ^2 with D=supx,y∈Ω‖x−y‖D= _x,y∈ \|x-y\| and λ being the relative strength scale of the spherical Hellinger component. 2. If ξ(k)∈Bi(r)ξ^(k)∈ B_i(r) for all k∈ℕ∪0k ∪\0\, the sequence of SHK gradient iterates (ξ(k))k≥0 (ξ^(k) )_k≥ 0 converges weakly to Xi∗(r)X_i^*(r) in (Ω)P( ). 3. If ξ(k)∈Bi(r)ξ^(k)∈ B_i(r) for all k∈ℕ∪0k ∪\0\, then for any δ>0δ>0, the error bound Sε(ξ(k),Xi∗(r))≤δS_ (ξ^(k),X_i^*(r) )≤δ is guaranteed to be achieved once the number of iterations k is greater than or equal to min2ηD2λ2,−logamin+2logG+log(2η(E(ξ(0))−Ei∗(r)))+2log(1δ(1−1−μη))−log(1−μη) \ 2η D^2λ^2,- a_ \+2 G+ (2η (E (ξ^(0) )-E_i^*(r) ))+2 ( 1δ(1- 1-μη))- (1-μη). A simpler sufficient condition on the no. of iterations to achieve the same error bound is k≥min2D2μλ2,−logaminμη+1μηlog(2η(E(ξ(0))−Ei∗(r)))+2μηlog(2Gδμη)k≥ \ 2D^2μλ^2,- a_ μη \+ 1μη (2η (E (ξ^(0) )-E_i^*(r) ))+ 2μη ( 2Gδμη). In addition, for ρi(η,r,ξ(0))≔G2η(E(ξ(0))−Ei∗(r))amin(1−1−μη) _i(η,r,ξ^(0)) G 2η (E (ξ^(0) )-E_i^*(r) ) a_ (1- 1-μη) and α(r,ξ(0))≔r2×μamin2G2(E(ξ(0))−Ei∗(r))α(r,ξ^(0)) r^2× μ a_ 2G^2 (E (ξ^(0) )-E_i^*(r) ), if the step-size η and the initial iterate ξ(0)ξ^(0) satisfies the conditions α(r,ξ(0))>1α(r,ξ^(0))>1, 4α(r,ξ(0))μ(α(r,ξ(0))+1)2≤η<min1/L,1/μ,ηret,i 4α(r,ξ^(0))μ(α(r,ξ^(0))+1)^2≤η< \1/L,1/μ, _ret,i\ and Sε(ξ(0),Xi)≤r−ρi(η,r,ξ(0))S_ (ξ^(0),X_i)≤ r- _i(η,r,ξ^(0)) i.e. ξ(0)∈Bi(r−ρi)ξ^(0)∈ B_i(r- _i), then the sequence of SHK gradient descent iterates (ξ(k))k≥0(ξ^(k))_k≥ 0 all belong to the local basin Bi(r)B_i(r). Consequently, the above 3 properties hold true without the apriori basin invariance assumption ξ(k)∈Bi(r)ξ^(k)∈ B_i(r) for all k∈ℕ∪0k ∪\0\. Theorem 2 provides the core algorithmic guarantee, showing that the SHK gradient descent dynamics of E converges geometrically (in terms of no. of iterations) to the unique local minimizer within that basin and remain invariant inside it. Together, these results connect statistical separation with dynamical stability, yielding rigorous theoretical guarantees for accurate associative-memory retrieval. Theorem 3 (Sinkhorn distance between minimizer in local basin and stored pattern and stability of stored pattern). Let Assumptions (A1) and (A3) hold true. Then, we have that 1. Sε(Xi∗(r),Xi)≤1βlog(1+(N−1)e−βΔ)≤N−1βe−βΔ.S_ (X_i^*(r),X_i)≤ 1β (1+(N-1)e^-β )≤ N-1βe^-β . 2. for η≤ηret,iη≤ _ret,i as defined in Theorem 2, Sε(Xi,Φη(Xi))≤mineηD2/λ2,1aminηG2(N−1)e−βΔ1+(N−1)e−βΔS_ (X_i, _η(X_i))≤ \e^η D^2/λ^2, 1 a_ \η G^2(N-1)e^-β 1+(N-1)e^-β . While Theorem 2 establishes geometric convergence of the SHK gradient descent iterates to the unique local minimizer within a Sinkhorn basin, it does not yet quantify how well this minimizer approximates the original stored pattern. Theorem 3 closes this gap by showing that the basin minimizer remains exponentially close (in terms of the inverse temperature β and margin separation Δ ) to the corresponding stored pattern in Sinkhorn divergence, and that each stored pattern is an approximate fixed point of the retrieval operator. Thus, beyond dynamical convergence, Theorem 3 provides a fidelity guarantee: the attractor reached by the algorithm is not merely stable, but provably close to the intended memory, ensuring accurate and stable associative recall. 4 Conclusion We developed a dense associative memory for empirical measures based on a Sinkhorn log-sum-exp energy and spherical Hellinger Kantorovich gradient dynamics, yielding deterministic transport-reaction retrieval with provable local convergence and separation guarantees. The full retrieval algorithm (pseudo-code and implementation details) along with the numerical experiments is provided in the Appendix (see Sections H and J). Future directions include faster retrieval implementations, adaptive regularization and applying our algorithm on real point clouds. References References L. Chizat (2022) Sparse optimization on measures with over-parameterized gradient descent. Mathematical Programming 194 (1), p. 487–532. Cited by: §D.4. J. Feydy, T. Séjourné, F. Vialard, S. Amari, A. Trouvé, and G. Peyré (2019) Interpolating between optimal transport and mmd using sinkhorn divergences. In The 22nd international conference on artificial intelligence and statistics, p. 2681–2690. Cited by: §G.2, §2. M. Hardion and H. Lavenant (2025) Gradient flows of potential energies in the geometry of sinkhorn divergences. arXiv preprint arXiv:2511.14278. Cited by: Appendix A, §2. J. J. Hopfield (1982) Neural networks and physical systems with emergent collective computational abilities. Proceedings of the national academy of sciences 79 (8), p. 2554–2558. Cited by: §1. D. Krotov and J. J. Hopfield (2016) Dense associative memory for pattern recognition. In Advances in Neural Information Processing Systems, D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Garnett (Eds.), Vol. 29, p. . External Links: Link Cited by: §1. M. Liero, A. Mielke, and G. Savaré (2018) Optimal entropy-transport problems and a new hellinger–kantorovich distance between positive measures. Inventiones mathematicae 211 (3), p. 969–1117. Cited by: §1. Y. Polyanskiy and Y. Wu (2025) Information theory: from coding to learning. Cambridge university press. Cited by: §I.4. H. Ramsauer, B. Schäfl, J. Lehner, P. Seidl, M. Widrich, T. Adler, L. Gruber, M. Holzleitner, M. Pavlović, G. K. Sandve, et al. (2020) Hopfield networks is all you need. arXiv preprint arXiv:2008.02217. Cited by: Appendix J, §1. C. Tankala and K. Balasubramanian (2026) Dense associative memory for gaussian distributions. External Links: 2509.23162, Link Cited by: §1. Appendix : Appendix A First variation of the Sinkhorn divergence Let Fν(μ)≔Sε(μ,ν)F_ν(μ) S_ (μ,ν) with ν fixed. A first variation δFν/δμδ F_ν/δμ is defined (up to an additive constant) by dtFν(μ+tχ)|t=0=∫Ω(δFνδμ(μ)(x))χ(x),∫χ=0. . ddtF_ν(μ+tχ) |_t=0= _ ( δ F_νδμ(μ)(x) )\,dχ(x), dχ=0. It is standard to see that (derivation available in Hardion and Lavenant (2025)) δSε(μ,ν)δμ=fμ,ν−12(fμ,μ+gμ,μ)in C(Ω)/ℝ. δ S_ (μ,ν)δμ=f_μ,ν- 12 (f_μ,μ+g_μ,μ ) C( )/R. (4) Intuitively, δOTε(μ,ν)δμ=fμ,ν _ (μ,ν)δμ=f_μ,ν and δOTε(μ,ν)δν=gμ,ν _ (μ,ν)δν=g_μ,ν. Since μ appears in both marginals of the self term OTε(μ,μ)OT_ (μ,μ), the first variation of 12OTε(μ,μ) 12OT_ (μ,μ) w.r.t. μ is the average potential fμ,μsym≔12(fμ,μ+gμ,μ),f^sym_μ,μ 12 (f_μ,μ+g_μ,μ ), which is invariant under the Sinkhorn gauge transformation (f,g)↦(f+c,g−c)(f,g) (f+c,g-c). When one chooses a symmetric gauge for the self problem (i.e. fμ,μ=gμ,μf_μ,μ=g_μ,μ), which is possible for symmetric costs), (4) reduces to the commonly stated formula fμ,ν−fμ,μf_μ,ν-f_μ,μ in C(Ω)/ℝC( )/R. Appendix B First variation of the log-sum-exp energy Let Z(ξ)=∑i=1Nexp(−βSε(ξ,Xi))Z(ξ)= _i=1^N (-β S_ (ξ,X_i)) and define Gibbs weights wi(ξ)≔exp(−βSε(ξ,Xi))∑j=1Nexp(−βSε(ξ,Xj)).w_i(ξ) (-β S_ (ξ,X_i)) _j=1^N (-β S_ (ξ,X_j)). (5) Differentiating (1) and using (4) gives δEδξ(ξ)=∑i=1Nwi(ξ)(fξ,Xi−12(fξ,ξ+gξ,ξ))=(∑i=1Nwi(ξ)fξ,Xi)−12(fξ,ξ+gξ,ξ),in C(Ω)/ℝ. δ Eδξ(ξ)= _i=1^Nw_i(ξ) (f_ξ,X_i- 12 (f_ξ,ξ+g_ξ,ξ ) )= ( _i=1^Nw_i(ξ)f_ξ,X_i )- 12 (f_ξ,ξ+g_ξ,ξ ), C( )/R. (6) We will denote the centered first variation of the energy functional E(⋅)E(·) as uξ≔δEδξ(ξ)−⟨δEδξ(ξ),ξ⟩=δEδξ(ξ)−∫δEδξ(ξ)ξu_ξ δ Eδξ(ξ)- δ Eδξ(ξ),ξ = δ Eδξ(ξ)- δ Eδξ(ξ)dξ where ⟨⋅,⋅⟩ ·,· is the duality pairing. Appendix C From entropic potentials to deterministic barycentric maps For quadratic cost c(x,y)=12‖x−y‖2c(x,y)= 12\|x-y\|^2, the gradient of a Schrödinger potential has the barycentric form ∇fμ,ν(x)=x−Tμ→νε(x),Tμ→νε(x)≔∫yπμ,νε(dy∣x),∇ f_μ,ν(x)=x-T _μ→ν(x), T _μ→ν(x) y\,π _μ,ν(dy x), (7) where πμ,νε(dy∣x)π _μ,ν(dy x) is the conditional distribution under the optimal entropic coupling. The full derivation is provided in Section I.3. Thus, using (6), the transport velocity becomes v(x)=∑i=1Nwi(ξ)Tξ→Xiε(x)−Tξ→ξε(x).v(x)= _i=1^Nw_i(ξ)\,T _ξ→ X_i(x)-T _ξ→ξ(x). (8) Even though (6) involves the symmetric self potential 12(fξ,ξ+gξ,ξ) 12(f_ξ,ξ+g_ξ,ξ), the transport velocity depends only on gradients. For symmetric costs (such as 12‖x−y‖2 12\|x-y\|^2) and a self-coupling, the optimal entropic plan is symmetric and one has ∇fξ,ξ=∇gξ,ξ∇ f_ξ,ξ=∇ g_ξ,ξ, so ∇12(fξ,ξ+gξ,ξ)=∇fξ,ξ,∇ 12 (f_ξ,ξ+g_ξ,ξ )=∇ f_ξ,ξ, and the self-correction in (8) remains the usual barycentric map Tξ→ξεT _ξ→ξ. This is a deterministic vector field on the query support computed from barycentric projections of Sinkhorn plans. Appendix D Spherical Hellinger-Kantorovich gradient flow (Continuous-time) D.1 Transport++reaction based continuity equation We consider a transport-reaction dynamics of the form ∂tξt+∇⋅(ξtvt)=ξtrt, _t _t+∇·( _tv_t)= _t\,r_t, (9) where vtv_t is a velocity field (transport) and rtr_t is a scalar reaction rate (mass reweighting). A SHK gradient flow of E sets vt(x)=−∇(δEδξ(ξt)(x))=−∇uξt(x),rt(x)=−1λ2(δEδξ(ξt)(x)−⟨δEδξ(ξt),ξt⟩)=−1λ2uξt(x),v_t(x)=-∇\! ( δ Eδξ( _t)(x) )=-∇ u_ _t(x), r_t(x)=- 1λ^2 ( δ Eδξ( _t)(x)- δ Eδξ( _t), _t )=- 1λ^2u_ _t(x), (10) with a scale parameter λ>0λ>0 controlling the relative strength of the spherical Hellinger component. The subtraction of the mean ensures dt∫ξt=0 ddt d _t=0, so probability mass is preserved. D.2 Riemmanian structure induced by Spherical Hellinger-Kantorovich geometry on the space of probability measures Let ξ∈(Ω)ξ ( ) where (Ω)P( ) is the set of all probability measures defined on Ω equipped with the Spherical Hellinger-Kantorovich geometry. A (sufficiently regular) tangent vector at ξ can be represented by a pair (r,v)(r,v) consisting of a vector field v:Ω→ℝdv: ^d (transport velocity) and a scalar field r:Ω→ℝr: (reaction rate), subject to the mass constraint ∫Ωr(x)ξ(x)=0 _ r(x)dξ(x)=0 Given such (v,r)(v,r), the induced infinitesimal change of measure is the distribution ξ˙ ξ defined by the weak form dt∫Ωφξt=∫Ω∇φ(x)⋅v(x)ξt(x)+∫Ωφ(x)r(x)ξt(x) ddt _ d _t= _ ∇ (x)· v(x)d _t(x)+ _ (x)r(x)d _t(x) (11) which corresponds to the PDE ∂tξt+∇⋅(ξtv)=ξtr. _t _t+∇· ( _tv )= _tr. (12) Fix λ>0λ>0. Define the inner product on the tangent space at ξ by ⟨(r,v),(r′,v′)⟩SHK,ξ≔∫Ω(v⋅v′+λ2rr′)ξ (r,v), (r ,v ) _SHK,ξ _ (v· v +λ^2r )dξ for pairs satisfying ∫rξ=∫r′ξ=0 rdξ= r dξ=0 and this induces the metric tensor at ξ to be gξSHK((r,v),(r′,v′))=⟨(r,v),(r′,v′)⟩SHK,ξ.g_ξ^SHK ((r,v), (r ,v ) )= (r,v), (r ,v ) _SHK,ξ. The induced norm on the tangent space at ξ is ‖(r,v)‖SHK,ξ2=∫Ω(‖v‖22+λ2r2)ξ.\|(r,v)\|_SHK,ξ^2= _ (\|v\|_2^2+λ^2r^2 )dξ. D.3 The manifold of finitely supported discrete measures equipped with SHK geometry For any M∈ℕM , 1M>amin>0 1M>a_ >0 and Δmin>0 _ >0, define LocM,Δmin(Ω)≔x=(x1,…,xM)∈ΩM:minm≠n‖xm−xn‖>Δmin∈ΔM,amin∘×LocM,Δmin(Ω)Loc_M, _ ( ) \x=(x_1,…,x_M)∈ ^M: _m≠ n\|x_m-x_n\|> _ \∈ _M,a_ ×Loc_M, _ ( ) as the space of M pairwise separated locations/atoms corresponding to the space of finitely supported discrete probability distributions with exactly M atoms and let ΔM,amin∘≔a∈ℝM:am>amin,∑m=1Mam=1 _M,a_ \a ^M:a_m>a_ , _m=1^Ma_m=1 \ be the space of bounded probability weights associated with the M atoms. Further, define the associated parameter space (ordered particles) to be ℳM≔ΔM,amin∘×LocM,Δmin(Ω).M_M _M,a_ ×Loc_M, _ ( ). For any (a,x)∈ℳM(a,x) _M, we associate to it the discrete measure ∑m=1Mamδxm _m=1^Ma_m _x_m through the parametrization mapping Ξ:ℳM→(Ω)⊂(C1(Ω))∗ :M_M ( )⊂ (C^1( ) )^*, defined as Ξ(a,x)≔∑m=1Mamδxm∈(Ω). (a,x) _m=1^Ma_m _x_m ( ). Here Ck(Ω)C^k( ) is the class of k-times continuously differentiable functions on the domain Ω , and (Ck(Ω))∗ (C^k( ) )^* represents the space of all continuous linear functionals on Ck(Ω)C^k( ) i.e. the dual space of Ck(Ω)C^k( ). Similarly, Cc∞(Ω)C_c^∞ ( ) represents the the class of compactly supported infinitely differentiable functions on the domain Ω , with Cc∞(Ω)C_c^∞ ( ), with its corresponding continuous dual being (Cc∞(Ω))∗ (C_c^∞ ( ) )^*. Since all of our theory is focused on the discrete measures Ξ(a,x) (a,x) and functionals defined using such measures, and all such objects are invariant under permutations of the labels associated with the weight-location pairs, we will consider the quotient space ℳ~M=ℳM/M M_M=M_M/ S_M where M S_M is the symmetric group that acts freely on ℳMM_M by relabeling: σ⋅(a1,…,aM,x1,…,xM)=(aσ−1(1),…,aσ−1(M),xσ−1(1),…,xσ−1(M)).σ· (a_1,…,a_M,x_1,…,x_M )= (a_σ^-1(1),…,a_σ^-1(M),x_σ^-1(1),…,x_σ^-1(M) ). Since the action is free and M S_M is finite, the quotient ℳ~M M_M is a smooth finite-dimensional manifold, and it identifies canonically with the set of probability measures on Ω having exactly M pairwise separated support points and positive weights. For computations and developing the theory without introducing extra notational overhead, it is simplest to work on the ordered cover ℳMM_M. The mathematical objects of interest are permutation invariant, so the theory developed on ℳMM_M descends to the quotient space ℳ~M M_M. Hence, from here on we identify ℳ~M M_M with ℳMM_M itself, which is equivalent to treating (a,x)∈ℳM(a,x) _M as unordered tuples. We also identify (a,x)(a,x) with ξ=Ξ(a,x)ξ= (a,x) when convenient. A tangent vector at (a,x)(a,x) is a pair (δa,δx)∈ℝM×(ℝd)M(δ a,δ x) ^M× (R^d )^M with the simplex constraint ∑m=1Mδam=0. _m=1^Mδ a_m=0. Thus T(a,x)ℳM=(δa,δx):∑mδam=0T_(a,x)M_M= \(δ a,δ x): _mδ a_m=0 \ For a given (a,x)∈ℳM(a,x) _M and (δa,δx)∈T(a,x)ℳM(δ a,δ x)∈ T_(a,x)M_M, we now proceed to compute the differential of Ξ at (a,x)(a,x) in the direction (δa,δx)(δ a,δ x). Let (a(t),x(t))(a(t),x(t)) be any C1C^1 curve in ℳMM_M such that (a(0),x(0))=(a,x) and dt(a(t),x(t))|t=0=(a˙(0),x˙(0))=(δa,δx).(a(0),x(0))=(a,x)\,\ and \,\ ddt(a(t),x(t)) _t=0=( a(0), x(0))=(δ a,δ x). Let us define the induced curve of measures ξt≔Ξ(a(t),x(t))=∑m=1Mam(t)δxm(t). _t (a(t),x(t))= _m=1^Ma_m(t) _x_m(t). For every test function φ∈C1(Ω) ∈ C^1( ), we have that ∫Ωφξt=∑m=1Mam(t)φ(xm(t)). _ d _t= _m=1^Ma_m(t) (x_m(t) ). Differentiating at t=0t=0, we obtain dt|t=0∫Ωφξt= . ddt |_t=0 _ d _t= ∑m=1Ma˙m(0)φ(xm)+∑m=1Mam∇φ(xm)⋅x˙m(0) _m=1^M a_m(0) (x_m )+ _m=1^Ma_m∇ (x_m )· x_m(0) = = ∑m=1Mδamφ(xm)+∑m=1Mam∇φ(xm)⋅δxm. _m=1^Mδ a_m (x_m )+ _m=1^Ma_m∇ (x_m )·δ x_m. Hence dΞ(a,x)(δa,δx)d _(a,x)(δ a,δ x) is the distribution characterized by ⟨dΞ(a,x)(δa,δx),φ⟩=∑m=1Mδamφ(xm)+∑m=1Mam∇φ(xm)⋅δxm. d _(a,x)(δ a,δ x), = _m=1^Mδ a_m (x_m )+ _m=1^Ma_m∇ (x_m )·δ x_m. (13) Equivalently, dΞ(a,x)(δa,δx)=∑m=1Mδamδxm−∇⋅(∑m=1Mamδxmδxm)d _(a,x)(δ a,δ x)= _m=1^Mδ a_m _x_m-∇· ( _m=1^Ma_mδ x_m _x_m ) in the sense of distributions. Comparing Equation 13 with the weak form of the SHK evolution PDE in Equation 11, we can uniquely represent the tangent vector (distribution) dΞ(a,x)(δa,δx)d _(a,x)(δ a,δ x) using any pair (r,v)(r,v) such that r(xm)=δamamr(x_m)= δ a_ma_m and v(xm)=δxmv(x_m)=δ x_m. The mass constraint is automatically satisfied since ∫Ωrξ=∑m=1Mamδamam=∑m=1Mδam=0 _ rdξ= _m=1^Ma_m δ a_ma_m= _m=1^Mδ a_m=0. The uniqueness on the support points x1,…,xMx_1,…,x_M can be established using Lemma 17. Now, we define the pullback metric tensor on ℳMM_M at (a,x)(a,x) induced by gΞ(a,x)SHKg_ (a,x)^SHK as g(a,x)≔Ξ∗gΞ(a,x)SHK,g_(a,x) ^*g_ (a,x)^SHK, which means that, for ζ=(δa,δx)ζ=(δ a,δ x) and η=(δa′,δx′)η= (δ a ,δ x ) belonging to T(a,x)ℳMT_(a,x)M_M, g(a,x)(ζ,η)=gξSHK(dΞ(a,x)ζ,dΞ(a,x)η),ξ=Ξ(a,x).g_(a,x)(ζ,η)=g_ξ^SHK (d _(a,x)ζ,d _(a,x)η ), ξ= (a,x). Using the identification above, we have that r(xm)=δamam,v(xm)=δxm,r′(xm)=δam′am,v′(xm)=δxm′.r (x_m )= δ a_ma_m, v (x_m )=δ x_m, r (x_m )= δ a_m a_m, v (x_m )=δ x_m . Therefore, g(a,x)((δa,δx),(δa′,δx′)) g_(a,x) ((δ a,δ x), (δ a ,δ x ) ) =∑m=1Mam(δxm⋅δxm′+λ2δamamδam′am) = _m=1^Ma_m (δ x_m·δ x_m +λ^2 δ a_ma_m δ a_m a_m ) =∑m=1M(amδxm⋅δxm′+λ2amδamδam′). = _m=1^M (a_mδ x_m·δ x_m + λ^2a_mδ a_mδ a_m ). Consequently, the SHK-induced Riemannian structure on ℳMM_M can be described as follows. Define for (δa,δx),(δa′,δx′)∈T(a,x)ℳM(δ a,δ x), (δ a ,δ x )∈ T_(a,x)M_M : ⟨(δa,δx),(δa′,δx′)⟩(a,x)≔g(a,x)((δa,δx),(δa′,δx′))=∑m=1M(λ2amδamδam′+amδxm⋅δxm′). (δ a,δ x), (δ a ,δ x ) _(a,x) g_(a,x) ((δ a,δ x), (δ a ,δ x ) )= _m=1^M ( λ^2a_mδ a_mδ a_m +a_mδ x_m·δ x_m ). The corresponding norm is ‖(δa,δx)‖(a,x)2=∑m=1M(λ2am(δam)2+am‖δxm‖22).\|(δ a,δ x)\|_(a,x)^2= _m=1^M ( λ^2a_m (δ a_m )^2+a_m \|δ x_m \|_2^2 ). D.4 Defining an appropriate retraction map for ℳMM_M A local retraction is the standard way to define a “first-order accurate exponential map” used in Riemannian gradient descent. Let TℳMTM_M denote the tangent bundle corresponding to the manifold ℳMM_M, given by TℳM≔⊔(a,x)∈ℳMT(a,x)ℳM=((a,x),(δa,δx)):(a,x)∈ℳM,(δa,δx)∈T(a,x)ℳMTM_M _(a,x) _MT_(a,x)M_M= \((a,x),(δ a,δ x)):(a,x) _M,(δ a,δ x)∈ T_(a,x)M_M \, where T(a,x)ℳMT_(a,x)M_M is the tangent space at the parameter (a,x)∈ℳM(a,x) _M. Following the general definition (and consistent with retractions are used in particle evolutions for gradient descent, for e.g. see Chizat (2022)) we use the following definition of a retraction map: A smooth map Ret : TℳM→ℳMTM_M _M is a retraction if for each base point z=(a,x)∈ℳMz=(a,x) _M, its restriction RetzRet_z : TzℳM→ℳMT_zM_M _M satisfies (i) Retz(0)=zRet_z(0)=z, (i) DRetz(0)=IdDRet_z(0)=Id on TzℳMT_zM_M. It need not be well-defined everywhere but there exists an open set U(a,x)⊂T(a,x)ℳMU_(a,x)⊂ T_(a,x)M_M containing the origin on which it is well-defined. For our purpose, we will build Ret as a product of a retraction on ΔM,amin∘ _M,a_ and one on LocM,Δmin(Ω)Loc_M, _ ( ). For the retraction map on LocM,Δmin(Ω)Loc_M, _ ( ), define: Retxpos(δx)≔x+δx, i.e. Ret(x1,…,xM)pos(δx1,…,δxM)=(x1+δx1,…,xM+δxM).Ret_x^pos(δ x) x+δ x, i.e. _ (x_1,…,x_M )^pos (δ x_1,…,δ x_M )= (x_1+δ x_1,…,x_M+δ x_M ). and Retxpos(⋅)Ret_x^pos(·) is well-defined on Ux=δx∈TxLocM(Ω):‖δx‖∞<12min(d∂(x),sep(x)−Δmin)⊂TxLocMU_x= \δ x∈ T_xLoc_M( ):\|δ x\|_∞< 12 (d_∂(x),sep(x)- _ ) \⊂ T_xLoc_M. Clearly Retxpos(0)=xRet_x^pos(0)=x and DRetxpos(0)=IdDRet_x^pos(0)=Id. Therefore, this is a retraction on LocM,Δmin(Ω)Loc_M, _ ( ). Let a∈ΔM,amin∘a∈ _M,a_ . Its tangent space is TaΔM,amin∘=δa∈ℝM:∑mδam=0.T_a _M,a_ = \δ a ^M: _mδ a_m=0 \. Define, for δa∈TaΔM,amin∘δ a∈ T_a _M,a_ , the retraction map Retaw(δa)≔a⊙exp(δa/a)∑j=1Majexp(δaj/aj)Ret_a^w(δ a) a (δ a/a) _j=1^Ma_j (δ a_j/a_j ) (14) where (δa/a)m≔δam/am(δ a/a)_m δ a_m/a_m, exp acts coordinatewise and ⊙ is coordinatewise product. This is well-defined for all δa∈Uaδ a∈ U_a where Ua=δa∈TaΔM,amin∘:‖δa‖∞<12minmlogamamin.U_a= \δ a∈ T_a _M,a_ : \| δ aa \|_∞< 12 _m a_ma_ \.. Further, it can be shown that Retaw(δa)Ret_a^w(δ a) is indeed a retraction map for ΔM,amin∘ _M,a_ . Finally, the product retraction on ℳMM_M is defined as Ret(a,x)(δa,δx)≡RetΞ(a,x)(δa,δx)≔(Retaw(δa),Retxpos(δx))Ret_(a,x)(δ a,δ x) _ (a,x)(δ a,δ x) (Ret_a^w(δ a),Ret_x^pos(δ x) ) which is a retraction on U(a,x)=(δa,δx)∈T(a,x)ℳM:δa∈Ua and δx∈Ux⊂T(a,x)ℳMU_(a,x)= \(δ a,δ x)∈ T_(a,x)M_M:δ a∈ U_a and δ x∈ U_x \⊂ T_(a,x)M_M. D.5 Riemannian gradient of the energy functional E on ℳMM_M Let ξ=Ξ(a,x)ξ= (a,x). For discrete ξ, write um≔uξ(xm)u_m u_ξ (x_m ). Differential of E in particle coordinates : Consider a tangent perturbation (δa,δx)(δ a,δ x). Assume that the functional E defined on (Ω)P( ) admits a C1(Ω)C^1( ) first variation δEδξ(ξ) δ Eδξ(ξ) for distributions of interest ξ∈(Ω)ξ ( ), whose centered version is uξ≔δEδξ(ξ)−⟨δEδξ(ξ),ξ⟩=δEδξ(ξ)−∫ΩδEδξ(ξ)ξu_ξ δ Eδξ(ξ)- δ Eδξ(ξ),ξ = δ Eδξ(ξ)- _ δ Eδξ(ξ)dξ. Now, for a given (a,x)∈ℳM(a,x) _M and (δa,δx)∈T(a,x)ℳM(δ a,δ x)∈ T_(a,x)M_M, (a(t),x(t))(a(t),x(t)) be any C1C^1 curve in ℳMM_M such that (a(0),x(0))=(a,x) and dt(a(t),x(t))|t=0=(a˙(0),x˙(0))=(δa,δx).(a(0),x(0))=(a,x)\,\ and \,\ ddt(a(t),x(t)) _t=0=( a(0), x(0))=(δ a,δ x). and let the induced curve of measures be ξt≔Ξ(a(t),x(t))=∑m=1Mam(t)δxm(t). _t (a(t),x(t))= _m=1^Ma_m(t) _x_m(t). Along the curve of measures ξt _t, by the chain rule, we have that the differential of E∘ΞE is given by d(E∘Ξ)(a,x)[δa,δx]= d (E )_(a,x) [δ a,δ x ]= dt(E∘Ξ)(a(t),x(t))|t=0 ddt (E ) (a(t),x(t) ) _t=0 (15) = = dtE(ξt)|t=0 ddtE ( _t ) _t=0 = = ⟨δEδξ(ξ0),ξ˙0⟩ δ Eδξ( _0), ξ_0 = = ⟨δEδξ(ξ0),dΞ(a,x)(δa,δx)⟩ δ Eδξ( _0),d (a,x)(δ a,δ x) = = ⟨δEδξ(ξ0),∑m=1Mδamδxm−∇⋅(∑m=1Mamδxmδxm)⟩ δ Eδξ( _0), _m=1^Mδ a_m _x_m-∇· ( _m=1^Ma_mδ x_m _x_m ) = = ∑m=1MδEδξ(ξ0)(xm)δam+∑m=1Mam∇δEδξ(ξ0)(xm)⋅δxm.. _m=1^M δ Eδξ( _0) (x_m )δ a_m+ _m=1^Ma_m∇ δ Eδξ( _0) (x_m )·δ x_m.. We note that, since ∑m=1Mδm=0 _m=1^M _m=0, one can replace δEδξ(ξ0) δ Eδξ( _0) by δEδξ(ξ0)+c δ Eδξ( _0)+c in Equation 15 for any constant c without changing the result, which is often termed as gauge-invariance. In particular, d(E∘Ξ)(a,x)[δa,δx]=∑m=1Muξ(xm)δam+∑m=1Mam∇uξ(xm)⋅δxm.d (E )_(a,x) [δ a,δ x ]= _m=1^Mu_ξ (x_m )δ a_m+ _m=1^Ma_m∇ u_ξ (x_m )·δ x_m. (16) We will often abuse notation by using the shorthand ξ≡Ξ[a,x]ξ≡ [a,x] and treating E as a functional directly over ℳMM_M, in which case we will represent the differential of E as dE(ξ)[δa,δx]=∑m=1Muξ(xm)δam+∑m=1Mam∇uξ(xm)⋅δxmdE(ξ)[δ a,δ x]= _m=1^Mu_ξ (x_m )δ a_m+ _m=1^Ma_m∇ u_ξ (x_m )·δ x_m Riemannian gradient of E : The Riemannian gradient gradSHKE(a,x)∈T(a,x)ℳMgrad_SHKE(a,x)∈ T_(a,x)M_M is the unique tangent vector (β,ζ) (β,ζ ) such that for all (δa,δx)(δ a,δ x), g(a,x)((β,ζ),(δa,δx))=⟨(β,ζ),(δa,δx)⟩(a,x)=dE(ξ)[δa,δx].g_(a,x) ( (β,ζ ),(δ a,δ x) )= (β,ζ ),(δ a,δ x) _(a,x)=dE(ξ)[δ a,δ x]. With ξ≡Ξ(a,x)ξ≡ (a,x), define tξ,m=δEδξ(ξ)(xm)+ct_ξ,m= δ Eδξ(ξ)(x_m)+c for any constant c∈ℝc . Then, we have that ∑m(λ2amβmδam+amζm⋅δxm)=∑m(tξ,mδam+am∇tξ,m⋅δxm). _m ( λ^2a_m _mδ a_m+a_m _m·δ x_m )= _m (t_ξ,mδ a_m+a_m∇ t_ξ,m·δ x_m ). Matching the δxmδ x_m terms gives ζm=∇tξ,m=∇uξ(xm). _m=∇ t_ξ,m=∇ u_ξ (x_m ). For weights, note δaδ a is constrained by ∑mδam=0 _mδ a_m=0. The identity ∑m(λ2amβm−tξ,m)δam=0∀δa:∑mδam=0 _m ( λ^2a_m _m-t_ξ,m )δ a_m=0 ∀δ a: _mδ a_m=0 holds if and only if the coefficients km≔λ2amβm−tξ,mk_m λ^2a_m _m-t_ξ,m are all equal to the same constant k, which can be verified by taking δa=ei−ejδ a=e_i-e_j for all i≠ji≠ j. Hence λ2amβm=tξ,m+k. λ^2a_m _m=t_ξ,m+k. Since, we must have ∑m=1Mβm=0 _m=1^M _m=0, we must have that 0=∑m=1Mβm=1λ2∑m=1Mam(tξ,m+k)=1λ2(t¯+k),0= _m=1^M _m= 1λ^2 _m=1^Ma_m (t_ξ,m+k )= 1λ^2( t+k), where t¯≔∑m=1Mamtξ,m. t _m=1^Ma_mt_ξ,m. Therefore k=−t¯k=- t, so βm=amλ2(tξ,m−t¯)=amλ2um=amλ2uξ(xm). _m= a_mλ^2 (t_ξ,m- t )= a_mλ^2u_m= a_mλ^2u_ξ(x_m). So the Riemannian gradient of E is gradSHKE(a,x)=((amλ2um)m=1M,(∇um)m=1M).grad_SHKE(a,x)= ( ( a_mλ^2u_m )_m=1^M, (∇ u_m )_m=1^M ). D.6 SHK gradient descent updates in terms of retraction maps Fix a step size η>0η>0. The retraction-based SHK gradient descent update is as follows. Given ξ(k)=Ξ(a(k),x(k))ξ^(k)= (a^(k),x^(k) ), define (a(k+1),x(k+1))≔Ret(a(k),x(k))(−ηgradSHKE(a(k),x(k)))=Ret(a(k),x(k))(−ηgradSHKE(ξ(k))), (a^(k+1),x^(k+1) ) _ (a^(k),x^(k) ) (- _SHKE (a^(k),x^(k) ) )=Ret_ (a^(k),x^(k) ) (- _SHKE (ξ^(k) ) ), and set ξ(k+1)≔Ξ(a(k+1),x(k+1))ξ^(k+1) (a^(k+1),x^(k+1) ) We can express the SHK gradient descent update using the operator Φη(ξ(k))=ξ(k+1) _η (ξ^(k) )=ξ^(k+1). For a measurable map T:Ω→ΩT: → , the pushforward T#ξT_\#ξ is defined by (T#ξ)(A)=ξ(T−1(A)), equivalently ∫φd(T#ξ)=∫φ∘Tξ. (T_\#ξ )(A)=ξ (T^-1(A) ), equivalently d (T_\#ξ )= Tdξ. For a measurable function ρ:Ω→(0,∞)ρ: →(0,∞), the reweighted measure ρξρξ is defined by (ρξ)(A)≔∫Aρ(x)ξ(x).(ρξ)(A) _Aρ(x)dξ(x). For a given ξ, define the pushforward map Tηξ(x)≔x−η∇uξ(x).T_η^ξ(x) x-η∇ u_ξ(x). and the reweighting map ρηξ(x)≔exp(−ηλ2uξ(x))∫Ωexp(−ηλ2uξ(z))ξ(z). _η^ξ(x) (- ηλ^2u_ξ(x) ) _ (- ηλ^2u_ξ(z) )dξ(z). Then the SHK gradient descent update is Φη(ξ)=(Tηξ)#(ρηξξ). _η(ξ)= (T_η^ξ )_\# ( _η^ξ ). D.7 Particle (empirical-measure) dynamics Let ξt=∑m=1Mam(t)δxm(t) _t= _m=1^Ma_m(t) _x_m(t). Plugging this ansatz into (9) yields the Lagrangian system x˙m(t)=vt(xm(t)),a˙m(t)=am(t)rt(xm(t)),m=1,…,M. x_m(t)=v_t(x_m(t)), a_m(t)=a_m(t)\,r_t(x_m(t)), m=1,…,M. (17) Define the particle-wise first-variation values zm≔(δEδξ(ξ))(xm)=∑i=1Nwi(ξ)(fξ,Xi(xm)−12(fξ,ξ(xm)+gξ,ξ(xm))),z_m\; \; ( δ Eδξ(ξ) )(x_m)= _i=1^Nw_i(ξ) (f_ξ,X_i(x_m)- 12 (f_ξ,ξ(x_m)+g_ξ,ξ(x_m) ) ), (18) and their ξ-average z¯≔∑m=1Mamzm z _m=1^Ma_mz_m. Then (10) gives the weight dynamics a˙m(t)=−1λ2am(t)(zm−z¯),∑mam(t)=1. a_m(t)=- 1λ^2\,a_m(t)\, (z_m- z ), _ma_m(t)=1. (19) Equation (19) is the (mass-preserving) replicator equation and can be viewed as a spherical Hellinger / natural-gradient flow on the simplex. Appendix E Assumptions for analysis of the Sinkhorn and SHK gradient flow based retrieval In this section, we list down some useful assumptions are useful to prove certain properties of the Spherical Hellinger-Kantorovich gradient descent. Define Fi(ξ)=Sε(ξ,Xi)F_i(ξ)=S_ (ξ,X_i) for =1,…,N=1,…,N. Assumption (A1) (Margin separation). The stored patterns X1,…,XNX_1,…,X_N and ε are such that there exists a radius r>0r>0 and a corresponding Δ>0 >0 such that for any ξ∈Bi(r)≔μ∈M,amin,Δmin(Ω):Sε(μ,Xi)≤rξ∈ B_i(r) \μ _M,a_ , _ ( ):S_ (μ,X_i)≤ r \, we have that, for all j≠ij≠ i, Fj(ξ)−Fi(ξ)≥Δ.F_j(ξ)-F_i(ξ)≥ . (20) Assumption (A2) (Local Retraction L-smoothness of energy functional E in SHK geometry). For the same choice of r as in Assumption (A1), there exists L>0L>0 such that for every ξ∈Bi(r)ξ∈ B_i(r), every tangent vector (w,v)∈Tξ(w,v)(w,v)∈ T_ξ(w,v) and every η>0η>0 small enough so that Retξ(η(w,v))Ret_ξ(η(w,v)) is well-defined, we have E(Retξ(η(r,v)))≤E(ξ)+η⟨gradSHKE(ξ),(r,v)⟩SHK,ξ+Lη22‖(r,v)‖SHK,ξ2E (Ret_ξ(η(r,v)) )≤ E(ξ)+η _SHKE(ξ),(r,v) _SHK,ξ+ Lη^22\|(r,v)\|_SHK,ξ^2 (21) where ⟨⋅,⋅⟩SHK,ξ ·,· _SHK,ξ and ∥⋅∥SHK,ξ\|·\|_SHK,ξ is the SHK inner product and norm respectively . Let Ei∗(r)≔infξ∈Bi(r)E(ξ)E_i^*(r) _ξ∈\ B_i(r)E(ξ). We will often denote Ei∗(r)E_i^*(r) using Ei∗E_i^* for convenience. Assumption (A3) (Existence and uniqueness of minimizer of E in local basin). For the same choice of r as in Assumption (A1), there exists a minimizer Xi∗(r)∈Bi(r)X_i^*(r)∈ B_i(r) of E i.e. E(Xi∗(r))=Ei∗(r)≔infξ∈Bi(r)E(ξ)E(X_i^*(r))=E_i^*(r) _ξ∈ B_i(r)E(ξ). Further the minimizer Xi∗(r)X_i^*(r) is unique. Assumption (A4) (PL inequality in local basin). For the same choice of r as in Assumption (A1) and conditional on Assumption (A3) being true, there exists μ>0μ>0 such that for every ξ∈Bi(r)ξ∈ B_i(r), 12‖gradSHKE(ξ)‖SHK,ξ2≥μ(E(ξ)−Ei∗(r)) 12\|grad_SHKE(ξ)\|_SHK,ξ^2≥μ (E(ξ)-E_i^*(r) ) (22) Appendix F Sampling Algorithm to ensure high probability separation of measures in M,amin,Δmin(Ω)P_M,a_ , _ ( ) We construct a sampling model that (i) produces fully general M-atom measures (random weights, random supports) that belong to M,amin,Δmin(Ω)P_M,a_ , _ ( ), and (i) ensures pattern separation with high probability for number of patterns N exponentially large in d. For simplicity and concreteness, we assume that the domain Ω is such that there exists some point c∈Ωc∈ and some R>0R>0 such that the Euclidean ball ¯(c,R)≔x∈Ω:‖x−c‖≤R⊂Ω. B(c,R) \x∈ :\|x-c\|≤ R\⊂ . Sampling algorithm (SampAlgo) : Consider the domain Ω such that (c,R)≔x:‖x−c‖≤R⊆Ω⊂ℝdB(c,R) \x:\|x-c\|≤ R\ ^d. Consider the shape radius σ∈(0,R4)σ∈(0, R4), margin parameter γ∈(0,1)γ∈(0,1) and inradius parameter R0=R−2σR_0=R-2σ. Fix any ε>0 >0 such that εlogM<1−γ16R02 M< 1-γ16R_0^2. For each pattern i=1,…,Ni=1,…,N, 1. Generate random mean using Rademacher random variables : Sample si∈±1ds_i∈\± 1\^d with i.i.d. coordinates, ℙ(si,k=+1)=ℙ(si,k=−1)=1/2P (s_i,k=+1 )=P (s_i,k=-1 )=1/2. Define μi≔c+R0dsi _i c+ R_0 ds_i. 2. Generate random weights : Sample bi=(bi,1,…,bi,M)b_i= (b_i,1,…,b_i,M ) in i.i.d manner using any probability distribution supported on ΔM,amin∘ _M,a_ for some fixed amin>0a_ >0. 3. Sample random point cloud around means : Sample zi=(zi,1,…,zi,M)∈(ℝd)Mz_i=(z_i,1,…,z_i,M)∈(R^d)^M from any probability distribution supported in σ,Δmin≔(z1,…,zM)∈(0,σ)M:minn≠m‖zn−zm‖>ΔminZ_σ, _ \ (z_1,…,z_M ) (0,σ)^M: _n≠ m \|z_n-z_m \|> _ \. 4. Set means using mean correction and define support points: Compute the weighted shape mean z¯i≔∑m=1Mbi,mzi,m z_i _m=1^Mb_i,mz_i,m and define support points yi,m≔μi+(zi,m−z¯i),m=1,…,My_i,m _i+ (z_i,m- z_i ), m=1,…,M. 5. Define the discrete measure : Set the pattern as the discrete Xi≔∑m=1Mbi,mδyi,m∈M,amin,Δmin(Ω)X_i _m=1^Mb_i,m _y_i,m _M,a_ , _ ( ) Appendix G Proof of main theoretical results G.1 Proof of Theorem 1 Proof. We first verify that the probability distributions X1,…,XNX_1,…,X_N generated by the sampling algorithm SampAlgo (see Sec F) belong to M,amin,Δmin(Ω)P_M,a_ , _ ( ). First, by construction of Step 2 in SampAlgo (see Sec F), we have that, for i=1,…,Ni=1,…,N, bi∈ΔM,amin∘b_i∈ _M,a_ . Next, we prove strict pairwise separation of the support points yi,my_i,m. For m≠nm≠ n, yi,n−yi,m=(μi+zi,n−z¯i)−(μi+zi,m−z¯i)=zi,n−zi,m.y_i,n-y_i,m= ( _i+z_i,n- z_i )- ( _i+z_i,m- z_i )=z_i,n-z_i,m. and hence ‖yi,n−yi,m‖=‖zi,n−zi,m‖ \|y_i,n-y_i,m \|= \|z_i,n-z_i,m \| for i=1,…,Ni=1,…,N. Since zi∈σ,Δminz_i _σ, _ , we have that minm≠n‖zi,n−zi,m‖>Δmin _m≠ n \|z_i,n-z_i,m \|> _ . Therefore, we have that minm≠n‖yi,n−yi,m‖>Δmin. _m≠ n \|y_i,n-y_i,m \|> _ . Finally, we prove that every support point yi,my_i,m lies in Ω . Since each zi,m∈¯(0,σ)z_i,m∈ B(0,σ), we have that ‖zi,m‖≤σ \|z_i,m \|≤σ. Further, since z¯i=∑m=1Mbi,mzi,m z_i= _m=1^Mb_i,mz_i,m is a convex combination of the zi,mz_i,m, we have that ‖z¯i‖≤∑m=1Mbi,m‖zi,m‖≤∑m=1Mbi,mσ=σ. \| z_i \|≤ _m=1^Mb_i,m \|z_i,m \|≤ _m=1^Mb_i,mσ=σ. Further, ‖μi−c‖=‖R0dsi‖=R0d‖si‖=R0dd=R0 \| _i-c \|= \| R_0 ds_i \|= R_0 d \|s_i \|= R_0 d d=R_0 Hence, for each m, we have that ‖yi,m−c‖≤‖μi−c‖+‖zi,m‖+‖z¯i‖≤R0+σ+σ=R. \|y_i,m-c \|≤ \| _i-c \|+ \|z_i,m \|+ \| z_i \|≤ R_0+σ+σ=R. Therefore, we have that yi,m∈¯(c,R)⊂Ω.y_i,m∈ B(c,R)⊂ . Therefore, for i=1,…,Ni=1,…,N, the generated pattern XiX_i indeed belongs to M,amin,Δmin(Ω)P_M,a_ , _ ( ). Now, we compute its mean: m(Xi)=∑m=1Mbi,myi,m= m (X_i )= _m=1^Mb_i,my_i,m= ∑m=1Mbi,m(μi+zi,m−z¯i) _m=1^Mb_i,m ( _i+z_i,m- z_i ) = = μi+∑m=1Mbi,mzi,m−∑m=1Mbi,mz¯i _i+ _m=1^Mb_i,mz_i,m- _m=1^Mb_i,m z_i = = μi+z¯i−z¯i=μi. _i+ z_i- z_i= _i. Now, define the event A≔∀1≤i<j≤N,∥μi−μj∥≥dmin=2(1−γ)R0.A \∀ 1≤ i<j≤ N, \| _i- _j \|≥ d_ = 2(1-γ)R_0 \. choosing N≔⌊2pexp(γ24d)⌋N 2p ( γ^24d ) , we have that (N2)exp(−γ2d2)=N(N−1)2exp(−γ2d2)≤p N2 (- γ^2d2 )= N(N-1)2 (- γ^2d2 )≤ p. Therefore, by Lemma 16, we have that ℙ(A)≥1−p.P(A)≥ 1-p. Now , under the event A, the means μi _i are pairwise dmind_ -separated, Lemma 11 applies. Thus, for every i, every ξ∈M,amin,Δmin(Ω)ξ _M,a_ , _ ( ) with Sε(ξ,Xi)≤rS_ (ξ,X_i )≤ r, and every j≠ij≠ i, Sε(ξ,Xj)−Sε(ξ,Xi)≥Δ=dmin24.S_ (ξ,X_j )-S_ (ξ,X_i )≥ = d_ ^24. This proves the margin-separation statement. It remains to prove pairwise disjointness of the basins. Suppose for contradiction that for some i≠ji≠ j there exists ξ∈Bi(r)∩Bj(r).ξ∈ B_i(r)∩ B_j(r). Since ξ∈Bi(r)ξ∈ B_i(r), the margin-separation statement with indices (i,j)(i,j) gives Sε(ξ,Xj)−Sε(ξ,Xi)≥Δ.S_ (ξ,X_j )-S_ (ξ,X_i )≥ . Since ξ∈Bj(r)ξ∈ B_j(r), the same statement with indices (j,i)(j,i) gives Sε(ξ,Xi)−Sε(ξ,Xj)≥Δ.S_ (ξ,X_i )-S_ (ξ,X_j )≥ . Adding these two inequalities yields 0≥2Δ,0≥ 2 , which is impossible since Δ=dmin24>0 = d_ ^24>0. Hence Bi(r)∩Bj(r)=∅ for all i≠j.B_i(r)∩ B_j(r)= for all i≠ j. This completes the proof. ∎ G.2 Proof of Theorem 2 Proof. Assume first that ξ(k)∈Bi(r)ξ^(k)∈ B_i(r) for all k≥0k≥ 0. By Lemma 1, after relabeling the atoms of each iterate, there is a unique ordered representative z(k)=(a(k),x(k))∈Ki(δi,τi)z^(k)= (a^(k),x^(k) )∈ K_i ( _i, _i ) such that ξ(k)=Ξ(z(k)).ξ^(k)= (z^(k) ). Since η<ηret,iη< _ret,i, Lemma 2 implies that the local retraction is well-defined at each iterate. Let ξ(k+1)=Φη(ξ(k))ξ^(k+1)= _η (ξ^(k) ). For each k, define the one-step retraction curve γk(t)≔Retξ(k)(−tηgradSHKE(ξ(k))),t∈[0,1]. _k(t) _ξ^(k) (-t _SHKE (ξ^(k) ) ), t∈[0,1]. Then γk(0)=ξ(k) _k(0)=ξ^(k) and γk(1)=ξ(k+1) _k(1)=ξ^(k+1). Applying Lemma 7 with (δa,δx)=−ηgradSHKE(ξ(k))(δ a,δ x)=- _SHKE (ξ^(k) ) in particle coordinates yields Length(γk)≤e‖δa/a‖∞η‖gradSHKE(ξ(k))‖SHK,ξ(k).Length ( _k )≤ e^\|δ a/a\|_∞η \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k). Now we proceed to bound ‖δa/a‖∞\|δ a/a\|_∞. For the SHK gradient descent direction, the weight update follows δam=−ηamλ2uξ(k)(xm)⟹δamam=−ηλ2uξ(k)(xm)δ a_m=-η a_mλ^2u_ξ^(k) (x_m ) δ a_ma_m=- ηλ^2u_ξ^(k) (x_m ) where uξ≔δEδξ(ξ)−⟨δEδξ(ξ),ξ⟩=δEδξ(ξ)−∫δEδξ(ξ)ξ.u_ξ δ Eδξ(ξ)- δ Eδξ(ξ),ξ = δ Eδξ(ξ)- δ Eδξ(ξ)dξ. Thus ‖δa‖∞≤ηλ2supx∈Ω|uξ(k)(x)|. \| δ aa \|_∞≤ ηλ^2 _x∈ |u_ξ^(k)(x) |. Now, δEδξ(ξ)=∑j=1Nwj(ξ)δSε(ξ,Xj)δξ δ Eδξ(ξ)= _j=1^Nw_j(ξ) δ S_ (ξ,X_j)δξ and hence uξ(x)=∑j=1Nwj(ξ)(ϕi(x)−∫ϕi(y)ξ(y))u_ξ(x)= _j=1^Nw_j(ξ) ( _i(x)- _i(y)dξ(y) ) where ϕi(x)≔(δSε(ξ,Xi)δξ)(x) _i(x) ( δ S_ (ξ,X_i)δξ )(x) From Lemma 2 Part (i), we have that supx∈Ω|uξ(x)|≤D2 _x∈ |u_ξ(x)|≤ D^2. Consequently, we have that ‖δa‖∞≤ηD2λ2. \| δ aa \|_∞≤ η D^2λ^2. Therefore, we have that Length(γk)≤eηD2/λ2η‖gradSHKE(ξ(k))‖SHK,ξ(k).Length ( _k )≤ e^η D^2/λ^2η \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k). Now, using Lemmas 7 and 4, we have that dSHK(ξ(k+1),ξ(k))≤ d_SHK (ξ^(k+1),ξ^(k) )≤ Length(γk)≤mineηD2/λ2,1aminη‖gradSHKE(ξ(k))‖SHK,ξ(k) ( _k )≤ \e^η D^2/λ^2, 1 a_ \η \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k) ≤ ≤ mineηD2/λ2,1aminη⋅2(E(ξ(0))−Ei∗(r))η(1−ημ)k2=Cqk2 \e^η D^2/λ^2, 1 a_ \η· 2 (E (ξ^(0) )-E_i^*(r) )η(1-ημ) k2=Cq k2 where C≔mineηD2/λ2,1amin⋅2η(E(ξ(0))−Ei∗(r))C \e^η D^2/λ^2, 1 a_ \· 2η (E (ξ^(0) )-E_i^*(r) ) and q≔1−μη∈(0,1)q 1-μη∈(0,1). Now, for integers m>km>k, dSHK(ξ(m),ξ(k))≤∑t=km−1dSHK(ξ(t+1),ξ(t))≤∑t=km−1Cqt/2≤C∑t=k∞qt/2=Cqk/21−q.d_SHK (ξ^(m),ξ^(k) )≤ _t=k^m-1d_SHK (ξ^(t+1),ξ^(t) )≤ _t=k^m-1Cq^t/2≤ C _t=k^∞q^t/2=C q^k/21- q. As k→∞k→∞, the RHS converges to zero. Therefore, (ξ(k)) (ξ^(k) ) is Cauchy in dSHKd_SHK. The set Ki(δi,τi)K_i ( _i, _i ) is compact by Lemma 1. Therefore the sequence z(k)z^(k) has a Euclidean-convergent subsequence: z(kj)→z∞=(a∞,x∞)∈Ki(δi,τi).z (k_j )→ z^∞= (a^∞,x^∞ )∈ K_i ( _i, _i ). Define ξ∞≔Ξ(z∞)∈PM(Ω)ξ^∞ (z^∞ )∈ P_M( ) Because Ki(δi,τi)⊂ΔM,amin∘×LocM,Δmin(Ω)K_i ( _i, _i )⊂ _M,a_ ×Loc_M, _ ( ), the limit belongs to the parameter space ΔM,amin∘×LocM,Δmin(Ω) _M,a_ ×Loc_M, _ ( ). Using Lemma 1, we have that dSHK(ξ(kj),ξ∞)≤Li‖z(kj)−z∞‖E→0.d_SHK (ξ (k_j ),ξ^∞ )≤ L_i \|z (k_j )-z^∞ \|_E→ 0. Hence, the subsequence converges to ξ∞ξ^∞ in dSHKd_SHK. Let δ>0δ>0. Since (ξ(k)) (ξ^(k) ) is dSHKd_SHK-Cauchy, there exists N such that dSHK(ξ(m),ξ(n))<δ/2∀m,n≥Nd_SHK (ξ^(m),ξ^(n) )<δ/2 ∀ m,n≥ N. Choose j so large that kj≥Nk_j≥ N and dSHK(ξ(kj),ξ∞)<δ/2.d_SHK (ξ (k_j ),ξ^∞ )<δ/2. Then, for every k≥Nk≥ N, dSHK(ξ(k),ξ∞)≤dSHK(ξ(k),ξ(kj))+dSHK(ξ(kj),ξ∞)<δ.d_SHK (ξ^(k),ξ^∞ )≤ d_SHK (ξ^(k),ξ (k_j ) )+d_SHK (ξ (k_j ),ξ^∞ )<δ. Therefore, ξ(k)ξ^(k) converges to ξ∞ξ^∞ in dSHKd_SHK. Using Lemma 10, each Fi=Sε(⋅,Xi)F_i=S_ (·,X_i) is G-Lipschitz in dSHKd_SHK, and hence, |Fi(ξ(k))−Fi(ξ∞)|≤GdSHK(ξ(k),ξ∞)→0. |F_i (ξ^(k) )-F_i (ξ^∞ ) |≤ Gd_SHK (ξ^(k),ξ^∞ )→ 0. Since Fi(ξ(k))≤rF_i (ξ^(k) )≤ r for all k, Fi(ξ∞)≤rF_i (ξ^∞ )≤ r. Hence, we have that ξ∞∈Bi(r)ξ^∞∈ B_i(r). Again, using Lemma 10, we have that |E(ξ(k))−E(ξ∞)|≤GdSHK(ξ(k),ξ∞)→0. |E (ξ^(k) )-E (ξ^∞ ) |≤ Gd_SHK (ξ^(k),ξ^∞ )→ 0. On the other hand, Lemma 4 gives E(ξ(k))→Ei∗(r)E (ξ^(k) )→ E_i^*(r). Therefore, we have that E(ξ∞)=Ei∗(r)E (ξ^∞ )=E_i^*(r). Therefore, ξ∞ξ^∞ is a minimizer of E on Bi(r)B_i(r). Under Assumption (A3), by uniqueness, we have that ξ∞=Xi∗(r)ξ^∞=X_i^*(r) This proves the convergence of the SHK gradient descent iterates in dSHKd_SHK to Xi∗(r)X_i*(r). Now, since dSHK(ξ(m),ξ(k))≤Cqk/21−qd_SHK (ξ^(m),ξ^(k) )≤ C q^k/21- q and dSHKd_SHK is continuous, we have that dSHK(Xi∗(r),ξ(k))=dSHK(ξ(∞),ξ(k))=limm→∞dSHK(ξ(m),ξ(k))≤Cqk/21−q.d_SHK (X_i^*(r),ξ^(k) )=d_SHK (ξ^(∞),ξ^(k) )= _m→∞d_SHK (ξ^(m),ξ^(k) )≤ C q^k/21- q. Finally, using Lemma 8 and the fact that Sε(Xi∗(r),Xi∗(r))=0S_ (X_i^*(r),X_i^*(r) )=0, we have that Sε(ξ(k),Xi∗(r))≤GdSHK(Xi∗(r),ξ(k))≤GCqk/21−q.S_ (ξ^(k),X_i^*(r) )≤ Gd_SHK (X_i^*(r),ξ^(k) )≤ GC q^k/21- q. Consequently, limk→∞Sε(ξ(k),Xi∗(r))=0 _k→∞S_ (ξ^(k),X_i^*(r))=0 and since the Sinkhorn divergence metrizes weak convergence (Theorem 1 of Feydy et al. (2019)), we have that ξ(k)ξ^(k) converges weakly to Xi∗(r)X_i^*(r). Now, Sε(ξ(k),Xi∗(r))≤GCqk/21−q≤δS_ (ξ^(k),X_i^*(r) )≤ GC q^k/21- q≤δ if k≥ k≥ 2log(mineηD2/λ2,1amin⋅G2η(E(ξ(0))−Ei∗(r))δ(1−1−μη))−log(1−μη) 2 ( \e^η D^2/λ^2, 1 a_ \· G 2η (E (ξ^(0) )-E_i^*(r) )δ(1- 1-μη) )- (1-μη) = = min2ηD2λ2,−logamin+2logG+log(2η(E(ξ(0))−Ei∗(r)))+2log(1δ(1−1−μη))−log(1−μη). \ 2η D^2λ^2,- a_ \+2 G+ (2η (E (ξ^(0) )-E_i^*(r) ))+2 ( 1δ(1- 1-μη))- (1-μη). Using 1+logx≤x1+ x≤ x for 0<x<10<x<1, we have that 1−log(1−μη)≤1μη 1- (1-μη)≤ 1μη. Further, we have that 1−1−μη=μη1+1−μη≥μη2⟹11−1−μη≤2μη1- 1-μη= μη1+ 1-μη≥ μη2 11- 1-μη≤ 2μη. Consequently, we have the sufficient condition k≥min2D2μλ2,−logaminμη+1μηlog(2η(E(ξ(0))−Ei∗(r)))+2μηlog(2Gδμη)k≥ \ 2D^2μλ^2,- a_ μη \+ 1μη (2η (E (ξ^(0) )-E_i^*(r) ))+ 2μη ( 2Gδμη). We will now establish that if Fi(ξ(0))≤r−ρi(η,r,ξ(0)),F_i (ξ^(0) )≤ r- _i(η,r,ξ^(0)), then ξ(k)∈Bi(r)ξ^(k)∈ B_i(r) for every k∈ℕk . We will prove by strong induction that Fi(ξ(k))≤rF_i (ξ^(k) )≤ r for all k. Note that, by assumption Fi(ξ(0))≤r−ρi(η,r,ξ(0))<rF_i (ξ^(0) )≤ r- _i(η,r,ξ^(0))<r. Now, assume that Fi(ξ(t))≤rF_i (ξ^(t) )≤ r for all t=0,1,…,kt=0,1,…,k. Then ξ(t)∈Bi(r)ξ^(t)∈ B_i(r) for those t, so Assumptions (A2) and (A4) apply on each of these iterates, and consequently Lemma 4 applies up to time k, yielding E(ξ(t))−Ei⋆(r)≤(1−μη)t(E(ξ(0))−Ei∗(r)),t=0,…,kE (ξ^(t) )-E_i (r)≤ (1-μη )^t (E (ξ^(0) )-E_i^*(r) ), t=0,…,k We now bound the increment of FiF_i from ξ(t)ξ^(t) to ξ(t+1)ξ^(t+1). To do so, we construct an explicit smooth curve in ℳMM_M connecting them and bound its SHK length. Along any absolutely continuous curve ζ(s)ζ(s) connecting ζ(0)=ξ(t)ζ(0)=ξ^(t) to ζ(1)=ξ(t+1)ζ(1)=ξ^(t+1) restricted to M,amin,Δmin(Ω)P_M,a_ , _ ( ), dsFi(ζ(s))=⟨gradSHKFi(ζ(s)),dsζ(s)⟩SHK,ζ(s)≤‖gradSHKFi(ζ(s))‖SHK,ζ(s)⋅‖dsζ(s)‖SHK,ζ(s). ddsF_i(ζ(s))= _SHKF_i(ζ(s)), ddsζ(s) _SHK,ζ(s)≤ \|grad_SHKF_i(ζ(s)) \|_SHK,ζ(s)· \| ddsζ(s) \|_SHK,ζ(s). Integrating between 0 and 1 and using the global bound ‖gradSHKFi(ξ)‖SHK,ξ≤G \|grad_SHKF_i(ξ) \|_SHK,ξ≤ G from Lemma 8, we obtain Fi(ξ(t+1))−Fi(ξ(t))≤G⋅Length(ζ).F_i (ξ^(t+1) )-F_i (ξ^(t) )≤ G·Length(ζ). Using Lemmas 4 and 7, we have that Fi(ξ(t+1))−Fi(ξ(t))≤ F_i (ξ^(t+1) )-F_i (ξ^(t) )≤ G⋅ηamin⋅2(E(ξ(0))−Ei∗(r))η(1−μη)t2 G· η a_ · 2 (E (ξ^(0) )-E_i^*(r) )η(1-μη) t2 = = G⋅2ηamin(E(ξ(0))−Ei∗(r))(1−μη)t2. G· 2ηa_ (E (ξ^(0) )-E_i^*(r) )(1-μη) t2. Now summing from t=0t=0 to t=kt=k : Fi(ξ(k+1))≤Fi(ξ(0))+G2ηamin(E(ξ(0))−Ei∗(r))∑t=0k(1−μη)t2.F_i (ξ^(k+1) )≤ F_i (ξ^(0) )+G 2ηa_ (E (ξ^(0) )-E_i^*(r) ) _t=0^k(1-μη) t2. Since ∑t=0k(1−μη)t2≤∑t=0∞(1−μη)t2=11−1−μη _t=0^k(1-μη) t2≤ _t=0^∞(1-μη) t2= 11- 1-μη, Fi(ξ(k+1))≤Fi(ξ(0))+G2η(E(ξ(0))−Ei∗(r))amin(1−1−μη)=Fi(ξ(0))+ρi(η,r,ξ(0)).F_i (ξ^(k+1) )≤ F_i (ξ^(0) )+G 2η (E (ξ^(0) )-E_i^*(r) ) a_ (1- 1-μη)=F_i (ξ^(0) )+ _i(η,r,ξ^(0)). By the assumed initialization condition Fi(ξ(0))≤r−ρi(η,r,ξ(0))F_i (ξ^(0) )≤ r- _i(η,r,ξ^(0)), we get Fi(ξ(k+1))≤rF_i (ξ^(k+1) )≤ r. This completes the induction. Thus all iterates remain in Bi(r)B_i(r). The lower bound on η can be derived based on the natural constraint that ρi(η,r,ξ(0))≤r _i(η,r,ξ^(0))≤ r must be satisfied. The, we must have that G2η(E(ξ(0))−Ei∗(r))amin(1−1−μη)≤r⇔2G2η(E(ξ(0))−Ei∗(r))amin(1−1−μη)2≤r2.G 2η (E (ξ^(0) )-E_i^*(r) ) a_ (1- 1-μη)≤ r 2G^2η (E (ξ^(0) )-E_i^*(r) )a_ (1- 1-μη)^2≤ r^2. Define u≔1−μηu 1-μη. Since 0<η<1μ0<η< 1μ, we have that 0<u<10<u<1. Further, η=1−u2μ=(1−u)(1+u)μη= 1-u^2μ= (1-u)(1+u)μ. Therefore, the condition reduces to 1+u1−u≤r2×μamin2G2(E(ξ(0))−Ei∗(r))≕α(r,ξ(0)). 1+u1-u≤ r^2× μ a_ 2G^2 (E (ξ^(0) )-E_i^*(r) ) α(r,ξ^(0)). Because 0<u<10<u<1, we have that 1+u1−u>1 1+u1-u>1. Therefore, a necessary condition is α(r,ξ(0))>1α(r,ξ^(0))>1. Moreover, if α(r,ξ(0))>1α(r,ξ^(0))>1, then the condition is equivalent to 1+u≤α(r,ξ(0))(1−u)⇔u≤α(r,ξ(0))−1α(r,ξ(0))+1⇔1−μη≤(α(r,ξ(0))−1α(r,ξ(0))+1)2.1+u≤α(r,ξ^(0))(1-u) u≤ α(r,ξ^(0))-1α(r,ξ^(0))+1 1-μη≤ ( α(r,ξ^(0))-1α(r,ξ^(0))+1 )^2. Consequently, the condition reduces to α(r,ξ(0))>1α(r,ξ^(0))>1 and η≥1μ×[1−(α(r,ξ(0))−1α(r,ξ(0))+1)2]=4α(r,ξ(0))μ(α(r,ξ(0))+1)2.η≥ 1μ× [1- ( α(r,ξ^(0))-1α(r,ξ^(0))+1 )^2 ]= 4α(r,ξ^(0))μ(α(r,ξ^(0))+1)^2. This completes the proof. ∎ G.2.1 Auxiliary results for proving Theorem 2 For any i=1,…,Ni=1,…,N consider the stored pattern Xi=∑m=1Mbi,mδyi,m∈M,amin,Δmin(Ω)X_i= _m=1^Mb_i,m _y_i,m _M,a_ , _ ( ) whose weight and locations parameters are denoted as bi≔(bi,1,…,bi,M)b_i (b_i,1,…,b_i,M ) and yi≔(yi,1,…,yi,M)y_i (y_i,1,…,y_i,M ), respectively. Let us define the stored pattern margins wi≔min1≤m≤M(bi,m−amin)>0,di∂≔min1≤m≤Mdist(yi,m,∂Ω)>0,si≔sep(yi)=minm≠n‖yi,m−yi,n‖2>Δmin. gatheredw_i _1≤ m≤ M (b_i,m-a_ )>0,\\ d_i^∂ _1≤ m≤ Mdist (y_i,m,∂ )>0,\\ s_i (y_i )= _m≠ n \|y_i,m-y_i,n \|_2> _ . gathered Further, define δ¯i≔mindi∂,si−Δmin2>0 δ_i \d_i^∂, s_i- _ 2 \>0. Then, for any 0<δ<δ¯i0<δ< δ_i and 0<τ<wi0<τ<w_i, let us define riloc(δ,τ)≔minaminδ22−εlogM,τ(si−δ)24−εlogM.r^loc_i(δ,τ) \ a_ δ^22- M, τ (s_i-δ )^24- M \. (23) Lemma 1 (local basin compactness inside parameter space). Fix i∈1,…,Ni∈\1,…,N\, and choose numbers δi,τi _i, _i such that 0<δi<δ¯i,0<τi<wi0< _i< δ_i, 0< _i<w_i. Assume 0<r<riloc(δi,τi)0<r<r_i^loc ( _i, _i ). Then every ξ=∑m=1Mamδxm∈Bi(r)ξ= _m=1^Ma_m _x_m∈ B_i(r) admits, after relabeling of its atoms, a unique ordered representative (a,x)(a,x) satisfying ‖a−bi‖1≤τi,xm∈¯(yi,m,δi) for every m=1,…,M. \|a-b_i \|_1≤ _i, x_m∈ B (y_i,m, _i ) for every m=1,…,M. Consequently, Bi(r)⊂Ξ(Ki(δi,τi)),B_i(r)⊂ (K_i ( _i, _i ) ), where Ki(δi,τi)≔(a,x):∑m=1Mam=1,‖a−bi‖1≤τi,xm∈¯(yi,m,δi)∀m,K_i ( _i, _i ) \(a,x): _m=1^Ma_m=1, \|a-b_i \|_1≤ _i,x_m∈ B (y_i,m, _i )∀ m \, and Ki(δi,τi)K_i ( _i, _i ) is a compact subset of the parameter space ℳM=ΔM,amin∘×LocM,Δmin(Ω)M_M= _M,a_ ×Loc_M, _ ( ). Proof. From Lemma 11, we have that OTε(ξ,ξ)≤εlogMOT_ (ξ,ξ)≤ M and OTε(Xi,Xi)≤εlogMOT_ (X_i,X_i )≤ M. Therefore, Sε(ξ,Xi)≥OTε(ξ,Xi)−εlogM.S_ (ξ,X_i ) _ (ξ,X_i )- M. (24) Now suppose ξ∈Bi(r)ξ∈ B_i(r). Assume for the sake of contradiction that some query atom xm0x_m_0 does not belong to ⋃m=1M¯(yi,m,δi) _m=1^M B (y_i,m, _i ). Then, for every target atom yi,my_i,m, ‖xm0−yi,m‖2>δi \|x_m_0-y_i,m \|_2> _i and hence every unit of mass transported out of row m0m_0 must pay least δi2/2 _i^2/2 in transport cost. Since the row mass equals am0≥amina_m_0≥ a_ , we must have that OTε(ξ,Xi)≥aminδi22.OT_ (ξ,X_i )≥ a_ _i^22. Using (24), Sε(ξ,Xi)≥aminδi22−εlogM.S_ (ξ,X_i )≥ a_ _i^22- M. But r<riloc(δi,τi)≤aminδi22−εlogMr<r^loc_i ( _i, _i )≤ a_ _i^22- M, contradicting Sε(ξ,Xi)≤rS_ (ξ,X_i )≤ r. Hence, each query atom lies in ⋃m¯(yi,m,δi) _m B (y_i,m, _i ). Now, assume for the sake of contradiction that for some n0n_0, no query atom belongs to ¯(yi,n0,δi) B (y_i,n_0, _i ). Then every unit of mass transported into column n0n_0 pays at least δi2/2 _i^2/2. Since the column mass equals bi,n0>amin b_i,n_0>a_min , OTε(ξ,Xi)≥bi,n0δi22≥aminδi22.OT_ (ξ,X_i )≥ b_i,n_0 _i^22≥ a_ _i^22. Again (24) gives Sε(ξ,Xi)≥aminδi22−εlogM>r,S_ (ξ,X_i )≥ a_ _i^22- M>r, which is a contradiction. So every closed ball ¯(yi,n,δi) B (y_i,n, _i ) contains at least one query atom. Since δi<δ¯i≤si−Δmin2<si2 _i< δ_i≤ s_i- _ 2< s_i2, the balls ¯(yi,n,δi) B (y_i,n, _i ) are pairwise disjoint. Based on our arguments above, all M query atoms lie in the union of these M disjoint balls, and each ball contains at least one query atom. Since there are exactly M query atoms, each ball contains exactly one. Therefore, after a unique relabeling, we may assume xm∈¯(yi,m,δi),m=1,…,M.x_m∈ B (y_i,m, _i ), m=1,…,M. (25) Let P=(Pmn)m=1,n=1M,MP= (P_mn )_m=1,n=1^M,M be any coupling matrix between ξ and XiX_i, i.e. P≥0,P=a,P⊤=bi.P≥ 0, P1=a, P 1=b_i. Then, we have that ∑m=1MPmm≤∑m=1Mmin(am,bi,m)=1−12‖a−bi‖1 _m=1^MP_m≤ _m=1^M (a_m,b_i,m )=1- 12 \|a-b_i \|_1 and the off-diagonal mass satisfies ∑m≠nPmn=1−∑m=1MPmm≥12‖a−bi‖1. _m≠ nP_mn=1- _m=1^MP_m≥ 12 \|a-b_i \|_1. (26) Now, fix some m≠nm≠ n. Then, using (25), we have that, ‖xm−yi,n‖2≥‖yi,m−yi,n‖2−‖xm−yi,m‖2≥si−δi \|x_m-y_i,n \|_2≥ \|y_i,m-y_i,n \|_2- \|x_m-y_i,m \|_2≥ s_i- _i. Therefore c(xm,yi,n)=12‖xm−yi,n‖22≥12(si−δi)2.c (x_m,y_i,n )= 12 \|x_m-y_i,n \|_2^2≥ 12 (s_i- _i )^2. So every unit of off-diagonal mass contributes a transport cost of at least 12(si−δi)2 12 (s_i- _i )^2. Using (25), we have that OTε(ξ,Xi)≥12(si−δi)2∑m≠nPmn≥(si−δi)24‖a−bi‖1.OT_ (ξ,X_i )≥ 12 (s_i- _i )^2 _m≠ nP_mn≥ (s_i- _i )^24 \|a-b_i \|_1. Combining with (24), we have that Sε(ξ,Xi)≥(si−δi)24‖a−bi‖1−εlogM.S_ (ξ,X_i )≥ (s_i- _i )^24 \|a-b_i \|_1- M. If ‖a−bi‖1>τi \|a-b_i \|_1> _i, then Sε(ξ,Xi)>(si−δi)24τi−εlogM≥riloc(δi,τi)>rS_ (ξ,X_i )> (s_i- _i )^24 _i- M≥ r_i^loc ( _i, _i )>r which leads to a contradiction. Thus, we must have that ‖a−bi‖1≤τi. \|a-b_i \|_1≤ _i. The defining constraints of Ki(δi,τi)K_i ( _i, _i ) are closed and bounded in the finite-dimensional Euclidean space ℝM×(ℝd)MR^M× (R^d )^M, so Ki(δi,τi)K_i ( _i, _i ) is compact. We only need to verify that Ki(δi,τi)K_i ( _i, _i ) is indeed a subset of the parameter space ΔM,amin∘×LocM,Δmin(Ω) _M,a_ ×Loc_M, _ ( ). If (a,x)∈Ki(δi,τi)(a,x)∈ K_i ( _i, _i ), then, for all m=1,…,Mm=1,…,M, am≥bi,m−τi>amin,a_m≥ b_i,m- _i>a_ , since τi<wi=minm(bi,m−amin) _i<w_i= _m (b_i,m-a_ ). Also, since δi<di∂ _i<d_i^∂, we have that, for every m=1,…,Mm=1,…,M ¯(yi,m,δi)⊂Ω for every m. B (y_i,m, _i )⊂ for every m. Finally, for m≠nm≠ n, we have that ‖xm−xn‖≥‖yi,m−yi,n‖−2δi≥si−2δi>Δmin \|x_m-x_n \|≥ \|y_i,m-y_i,n \|-2 _i≥ s_i-2 _i> _ since δi<si−Δmin2 _i< s_i- _ 2. Therefore, we finally have that Ki(δi,τi)⊂ΔM,amin∘×LocM,Δmin(Ω).K_i ( _i, _i )⊂ _M,a_ ×Loc_M, _ ( ). This completes the proof. ∎ Lemma 2 (uniform retraction domain and local metric upper bound). Assume 0<r<riloc(δi,τi)0<r<r_i^loc ( _i, _i ), and let Ki(δi,τi)K_i ( _i, _i ) be as in Lemma 1. Define ai−≔min1≤m≤Mbi,m−τi>amina_i^- _1≤ m≤ Mb_i,m- _i>a_ and ηret,i≔minλ22D2logai−amin,12Dmindi∂−δi,si−2δi−Δmin. _ret,i \ λ^22D^2 a_i^-a_ , 12D \d_i^∂- _i,s_i-2 _i- _ \ \. Then, the following hold true 1. For every ξ=Ξ(a,x)∈Ξ(Ki(δi,τi))ξ= (a,x)∈ (K_i ( _i, _i ) ), supz∈Ω|uξ(z)|≤D2,supz∈Ω‖∇uξ(z)‖≤D. _z∈ |u_ξ(z) |≤ D^2, _z∈ \|∇ u_ξ(z) \|≤ D. 2. If 0<η<ηret,i0<η< _ret,i, then for every ξ=Ξ(a,x)∈Ξ(Ki(δi,τi))ξ= (a,x)∈ (K_i ( _i, _i ) ), the retraction Retξ(−ηgradSHKE(ξ))Ret_ξ (- _SHKE(ξ) ) is well-defined. 3. If ‖(δa,δx)‖E2≔∑m=1M(δam)2+∑m=1M‖δxm‖2,\|(δ a,δ x)\|_E^2 _m=1^M (δ a_m )^2+ _m=1^M \|δ x_m \|^2, then on Ki(δi,τi)K_i ( _i, _i ), ‖(δa,δx)‖SHK,(a,x)2≤Li2‖(δa,δx)‖E2,Li≔maxλ2ai−,1.\|(δ a,δ x)\|_SHK,(a,x)^2≤ L_i^2\|(δ a,δ x)\|_E^2, L_i \ λ^2a_i^-,1 \. Consequently, for all (a,x),(a′,x′)∈Ki(δi,τi)(a,x), (a ,x )∈ K_i ( _i, _i ), dSHK(Ξ(a,x),Ξ(a′,x′))≤Li‖(a−a′,x−x′)‖E.d_SHK ( (a,x), (a ,x ) )≤ L_i \| (a-a ,x-x ) \|_E. Proof. For j=1,…,Nj=1,…,N, let ϕj,ξ(z)≔δSε(ξ,Xj)δξ(z) _j,ξ(z) δ S_ (ξ,X_j )δξ(z) and let ϕ¯j,ξ(z)≔ϕj,ξ(z)−∫ϕj,ξξ. φ_j,ξ(z) _j,ξ(z)- _j,ξdξ. Then uξ(z)=∑j=1Nwj(ξ)ϕ¯j,ξ(z),∑j=1Nwj(ξ)=1,wj(ξ)≥0.u_ξ(z)= _j=1^Nw_j(ξ) φ_j,ξ(z), _j=1^Nw_j(ξ)=1, w_j(ξ)≥ 0. Following the derivations in the proof of Lemma 8, specifically (31) and (32), for each j=1,…,Nj=1,…,N, we have that supz∈Ω‖∇ϕj,ξ(z)‖≤D,supz∈Ωϕj,ξ(z)−infz∈Ωϕj,ξ(z)≤D2. _z∈ \|∇ _j,ξ(z) \|≤ D, _z∈ _j,ξ(z)- _z∈ _j,ξ(z)≤ D^2. Hence, we obtain supz∈Ω|ϕ¯j,ξ(z)|≤D2 and supz∈Ω‖∇ϕ¯j,ξ(z)‖≤D. _z∈ | φ_j,ξ(z) |≤ D^2\,\ and \,\ _z∈ \|∇ φ_j,ξ(z) \|≤ D. Taking the convex combination with weights wj(ξ)w_j(ξ), we have that, supz∈Ω|uξ(z)|≤D2 and supz∈Ω‖∇uξ(z)‖≤D. _z∈ |u_ξ(z) |≤ D^2\,\ and \,\ _z∈ \|∇ u_ξ(z) \|≤ D. Let gradSHKE(a,x)=((amλ2uξ(xm))m=1M,(∇uξ(xm))m=1M).grad_SHKE(a,x)= ( ( a_mλ^2u_ξ (x_m ) )_m=1^M, (∇ u_ξ (x_m ) )_m=1^M ). The actual retraction step is therefore δam=−ηamλ2uξ(xm),δxm=−η∇uξ(xm).δ a_m=-η a_mλ^2u_ξ (x_m ), δ x_m=-η∇ u_ξ (x_m ). Using the above bound, ‖δa‖∞=maxm|δam|am≤ηD2λ2,‖δx‖∞,2≔maxm‖δxm‖≤ηD. \| δ aa \|_∞= _m |δ a_m |a_m≤η D^2λ^2, \|δ x\|_∞,2 _m \|δ x_m \|≤η D. (27) By Lemma 1, every (a,x)∈Ki(δi,τi)(a,x)∈ K_i ( _i, _i ) satisfies am≥ai−a_m≥ a_i^- for all m=1,…,Mm=1,…,M, d∂(x)≥di∂−δid_∂(x)≥ d_i^∂- _i and sep(x)−Δmin≥si−2δi−Δminsep(x)- _ ≥ s_i-2 _i- _ . Hence, under the local region of validity of the retraction map, as introduced in Section D.4, ‖δa‖∞<12logai−amin \| δ aa \|_∞< 12 a_i^-a_ and ‖δx‖∞,2<12mindi∂−δi,si−2δi−Δmin\|δ x\|_∞,2< 12 \d_i^∂- _i,s_i-2 _i- _ \ are sufficient to make the weight and position retractions well-defined. This is guaranteed by 0<η<ηret,i0<η< _ret,i. On Ki(δi,τi)K_i ( _i, _i ), we have am≥ai−a_m≥ a_i^-and am≤1a_m≤ 1. Therefore, ‖(δa,δx)‖SHK,(a,x)2=∑m=1Mλ2am(δam)2+∑m=1Mam‖δxm‖2≤λ2ai−∑m(δam)2+∑m‖δxm‖2≤Li2‖(δa,δx)‖E2.\|(δ a,δ x)\|_SHK,(a,x)^2= _m=1^M λ^2a_m (δ a_m )^2+ _m=1^Ma_m \|δ x_m \|^2≤ λ^2a_i^- _m (δ a_m )^2+ _m \|δ x_m \|^2≤ L_i^2\|(δ a,δ x)\|_E^2. Now Ki(δi,τi)K_i ( _i, _i ) is convex, since the weight constraints define a convex set, and each position constraint xm∈B¯(yi,m,δi)x_m∈ B (y_i,m, _i ) defines a convex set, and intersection of convex sets is a convex set again. Hence the straight segment joining two points of Ki(δi,τi)K_i ( _i, _i ) stays in Ki(δi,τi)K_i ( _i, _i ). Integrating the above pointwise metric upper bound along that straight segment yields dSHK(Ξ(a,x),Ξ(a′,x′))≤Li‖(a−a′,x−x′)‖E.d_SHK ( (a,x), (a ,x ) )≤ L_i \| (a-a ,x-x ) \|_E. This completes the proof. ∎ G.3 Proof of Theorem 3 Proof. Since Fi(Xi)=Sε(Xi,Xi)=0≤rF_i (X_i )=S_ (X_i,X_i )=0≤ r, we have Xi∈Bi(r)X_i∈ B_i(r). Assumption (A3) guarantees the existence and uniqueness of the minimizer Xi∗(r)∈Bi(r)X_i^*(r)∈ B_i(r). Now, we apply Lemma 5 to ξ=Xi∗(r)ξ=X_i^*(r), giving 0≤Fi(Xi∗(r))−E(Xi∗(r))≤1βlog(1+(N−1)e−βΔ)0≤ F_i (X_i^*(r) )-E (X_i^*(r) )≤ 1β (1+(N-1)e^-β ) Thus Fi(Xi∗(r))≤E(Xi∗(r))+1βlog(1+(N−1)e−βΔ).F_i (X_i^*(r) )≤ E (X_i^*(r) )+ 1β (1+(N-1)e^-β ). But Xi∗(r)X_i^*(r) minimizes E over Bi(r)B_i(r), and Xi∈Bi(r)X_i∈ B_i(r), hence E(Xi∗(r))≤E(Xi).E (X_i^*(r) )≤ E (X_i ). Also E(ξ)≤Fi(ξ)E(ξ)≤ F_i(ξ) for every ξ∈Bi(r)ξ∈ B_i(r), so E(Xi)≤Fi(Xi)=0.E (X_i )≤ F_i (X_i )=0. Combining, we have that Sε(Xi∗(r),Xi)= S_ (X_i^*(r),X_i)= Fi(Xi∗(r)) F_i (X_i^*(r) ) ≤ ≤ E(Xi∗(r))+1βlog(1+(N−1)e−βΔ) E (X_i^*(r) )+ 1β (1+(N-1)e^-β ) ≤ ≤ E(Xi)+1βlog(1+(N−1)e−βΔ) E (X_i )+ 1β (1+(N-1)e^-β ) ≤ ≤ 0+1βlog(1+(N−1)e−βΔ) 0+ 1β (1+(N-1)e^-β ) = = 1βlog(1+(N−1)e−βΔ). 1β (1+(N-1)e^-β ). Now, consider the one-step retraction curve γi(t)≔RetXi(−tηgradSHKE(Xi)),t∈[0,1]. _i(t) _X_i (-t _SHKE (X_i ) ), t∈[0,1]. Then γi(0)=Xi _i(0)=X_i and γi(1)=Φη(Xi) _i(1)= _η(X_i). Note that, as derived in Theorem 2, when η≤ηret,iη≤ _ret,i, the retraction Φη(Xi) _η(X_i) is well defined. Applying Lemma 7 with (δa,δx)=−ηgradSHKE(Xi)(δ a,δ x)=- _SHKE (X_i ) in particle coordinates yields Length(γi)≤mine‖δa‖∞,1aminη‖gradSHKE(Xi)‖SHK,Xi.Length ( _i )≤ \e^\| δ aa\|_∞, 1 a_ \η \|grad_SHKE (X_i ) \|_SHK,X_i. Now we proceed to bound ‖δa‖∞\| δ aa\|_∞. For the SHK gradient descent direction, the weight update follows δam=−ηamλ2uXi(xm)⟹δamam=−ηλ2uXi(xm).δ a_m=-η a_mλ^2u_X_i (x_m ) δ a_ma_m=- ηλ^2u_X_i (x_m ). where uξ≔δEδξ(ξ)−⟨δEδξ(ξ),ξ⟩=δEδξ(ξ)−∫δEδξ(ξ)ξ.u_ξ δ Eδξ(ξ)- δ Eδξ(ξ),ξ = δ Eδξ(ξ)- δ Eδξ(ξ)dξ. Thus ‖δa‖∞≤ηλ2supx∈Ω|uXi(x)|. \| δ aa \|_∞≤ ηλ^2 _x∈ |u_X_i(x) |. Now, δEδξ(ξ)=∑j=1Nwj(ξ)δSε(ξ,Xj)δξ δ Eδξ(ξ)= _j=1^Nw_j(ξ) δ S_ (ξ,X_j)δξ and hence uξ(x)=∑j=1Nwj(ξ)(ϕi(x)−∫ϕi(y)ξ(y))u_ξ(x)= _j=1^Nw_j(ξ) ( _i(x)- _i(y)dξ(y) ) where ϕi(x)≔(δSε(ξ,Xi)δξ)(x) _i(x) ( δ S_ (ξ,X_i)δξ )(x) Following the argument presented in the proof of Lemma 8, specifically (32), we have that, for any x,y∈Ωx,y∈ , |ϕi(x)−ϕi(y)|≤D2| _i(x)- _i(y)|≤ D^2 and hence for any ξ∈M,amin,Δmin(Ω)ξ _M,a_ , _ ( ), we have that supx∈Ω|uξ(x)|≤D2 _x∈ |u_ξ(x)|≤ D^2. Consequently, we have that ‖δa‖∞≤ηD2λ2. \| δ aa \|_∞≤ η D^2λ^2. Using Lemma 8, the fact that Sε(Xi,Xi)=0S_ (X_i,X_i)=0 and following the same argument as in the proof of Theorem 2, we have that Sε(Xi,Φη(Xi))≤ S_ (X_i, _η(X_i))≤ GdSHK(Xi,Φη(Xi)) Gd_SHK(X_i, _η(X_i)) ≤ ≤ GLength(γi) GLength ( _i ) ≤ ≤ mineηD2/λ2,1aminηG‖gradSHKE(Xi)‖SHK,Xi. \e^η D^2/λ^2, 1 a_ \η G \|grad_SHKE (X_i ) \|_SHK,X_i. Now, since Fi(ξ)=Sε(ξ,Xi)F_i(ξ)=S_ (ξ,X_i) has a global minimum at ξ=Xiξ=X_i and FiF_i is differentiable in the SHK sense, we must have that gradSHKFi(Xi)=0grad_SHKF_i(X_i)=0. Now, using Lemma 6, we have that gradSHKE(Xi)=∑j=1Nwj(Xi)gradSHKFj(Xi)=∑j≠iNwj(Xi)gradSHKFj(Xi).grad_SHKE(X_i)= _j=1^Nw_j(X_i)grad_SHKF_j(X_i)= _j≠ i^Nw_j(X_i)grad_SHKF_j(X_i). Consequently, using the gradient bound from Lemma 8, we have that ‖gradSHKE(Xi)‖SHK,Xi≤ \|grad_SHKE(X_i) \|_SHK,X_i≤ ∑j≠iwj(Xi)‖gradSHKFj(Xi)‖SHK,Xi _j≠ iw_j(X_i) \|grad_SHKF_j(X_i) \|_SHK,X_i ≤ ≤ G∑j≠iwj(Xi)≤G(N−1)e−βΔ1+(N−1)e−βΔ G _j≠ iw_j(X_i)≤ G(N-1)e^-β 1+(N-1)e^-β Therefore, we have that Sε(Xi,Φη(Xi))≤ S_ (X_i, _η(X_i))≤ mineηD2/λ2,1aminηG‖gradSHKE(Xi)‖SHK,Xi \e^η D^2/λ^2, 1 a_ \η G \|grad_SHKE (X_i ) \|_SHK,X_i ≤ ≤ mineηD2/λ2,1aminηG2(N−1)e−βΔ1+(N−1)e−βΔ \e^η D^2/λ^2, 1 a_ \η G^2(N-1)e^-β 1+(N-1)e^-β ≤ ≤ mineηD2/λ2,1aminηG2(N−1)e−βΔ \e^η D^2/λ^2, 1 a_ \η G^2(N-1)e^-β where we used u1+u≤u u1+u≤ u for u≥0u≥ 0 in the last inequality. This completes the proof. ∎ Appendix H Discrete retrieval algorithm for empirical measures We now give a fully explicit deterministic retrieval scheme obtained by an explicit Euler discretization of (17). The resulting update alternates a Kantorovich (support) step and a spherical Hellinger (weight) step, both driven by the same Sinkhorn computations. H.1 Discrete Sinkhorn objects Let ξ=∑m=1Mamδxmξ= _m=1^Ma_m _x_m and Xi=∑n=1Mbi,nδyi,nX_i= _n=1^Mb_i,n _y_i,n. An entropic coupling is a matrix Pi∈ℝ+M×MP_i _+^M× M with row/column sums a and bib_i. The barycentric projection map (as derived in Section C) on the query support is Tξ→Xiε(xm)=∑n=1MPi[m,n]amyi,n,Tξ→ξε(xm)=∑ℓ=1MP0[m,ℓ]amxℓ.T _ξ→ X_i(x_m)= _n=1^M P_i[m,n]a_m\,y_i,n, T _ξ→ξ(x_m)= _ =1^M P_0[m, ]a_m\,x_ . Moreover, the dual potentials returned by Sinkhorn (in the first argument) provide the discrete values fξ,Xi(xm)f_ξ,X_i(x_m) and fξ,ξ(xm)f_ξ,ξ(x_m) needed for (18). For the self-coupling OTε(ξ,ξ)OT_ (ξ,ξ), ξ appears in both marginals, so the weight update requires the symmetric combination fξ,ξsym(xm)≔12(fξ,ξ(xm)+gξ,ξ(xm)),f^sym_ξ,ξ(x_m) 12 (f_ξ,ξ(x_m)+g_ξ,ξ(x_m) ), where gξ,ξg_ξ,ξ is the Sinkhorn potential in the second argument. For the self-coupling, the correct (gauge-invariant) contribution to the weight gradient is the symmetric average 12(fξ,ξ+gξ,ξ) 12(f_ξ,ξ+g_ξ,ξ). H.2 Explicit Euler (Kantorovich transport) update Using (8), a step of size η>0η>0 updates the support points by xmk+1=xmk+η(∑i=1NwikTik(xmk)−T0k(xmk)),x_m^k+1=x_m^k+η ( _i=1^Nw_i^k\,T_i^k(x_m^k)-T_0^k(x_m^k) ), (28) where TikT_i^k denotes the barycentric map computed from the Sinkhorn coupling between ξkξ^k and XiX_i, and T0kT_0^k is the barycentric map from the self-coupling between ξkξ^k and itself. H.3 Multiplicative (spherical Hellinger) weight update A first-order discretization of (19) can be implemented as the multiplicative update amk+1∝amkexp(−ηλ2zmk),∑m=1Mamk+1=1,a_m^k+1\; \;a_m^k \! (- ηλ^2\,z_m^k ), _m=1^Ma_m^k+1=1, (29) where zmkz_m^k is defined by (18) at (xk,ak)(x^k,a^k). The normalization removes any additive-constant ambiguity in the potentials and guarantees amk+1>0a_m^k+1>0 whenever amk>0a_m^k>0. H.4 Full deterministic retrieval algorithm Algorithm 1: Entropic DDAM with spherical Hellinger-Kantorovich retrieval (empirical measures) Inputs: stored empirical measures Xii=1N\X_i\_i=1^N; query ξ0=∑mam0δxm0ξ^0= _ma_m^0 _x_m^0; parameters β>0β>0, ε>0 >0, step size η>0η>0, spherical Hellinger scale λ>0λ>0; number of iterations K (or a stopping criterion). For k=0,1,…,K−1k=0,1,…,K-1 do: 1. Sinkhorn couplings and costs. For each i, run Sinkhorn between ξkξ^k and XiX_i to obtain: (a) coupling matrix PikP_i^k, (b) source potential values fik[m]=fξk,Xi(xmk)f_i^k[m]=f_ξ^k,X_i(x_m^k), and (c) the entropic OT cost OTε(ξk,Xi)OT_ (ξ^k,X_i). Compute also the self-coupling between ξkξ^k and itself to obtain P0kP_0^k, both potential vectors f0k[m]=fξk,ξk(xmk)f_0^k[m]=f_ξ^k,ξ^k(x_m^k) and g0k[m]=gξk,ξk(xmk)g_0^k[m]=g_ξ^k,ξ^k(x_m^k) and the entropic self-OT cost OTε(ξk,ξk)OT_ (ξ^k,ξ^k). Finally, for each i, compute the OT cost OTε(Xi,Xi)OT_ (X_i,X_i). 2. Sinkhorn divergences and Gibbs weights. For each i, compute Sε(ξk,Xi)S_ (ξ^k,X_i) using (2), and set wik=exp(−βSε(ξk,Xi))∑j=1Nexp(−βSε(ξk,Xj)).w_i^k= (-β S_ (ξ^k,X_i)) _j=1^N (-β S_ (ξ^k,X_j)). 3. Barycentric maps. For each particle m and each pattern i, compute Tik(xmk)=∑nPik[m,n]amkyi,n,T0k(xmk)=∑ℓP0k[m,ℓ]amkxℓk.T_i^k(x_m^k)= _n P_i^k[m,n]a_m^k\,y_i,n, T_0^k(x_m^k)= _ P_0^k[m, ]a_m^k\,x_ ^k. 4. Support update (transport). Update xmk+1x_m^k+1 using (28). 5. Weight update (reaction). Compute zmk=∑i=1Nwik(fik[m]−12(f0k[m]+g0k[m])),z_m^k= _i=1^Nw_i^k (f_i^k[m]- 12 (f_0^k[m]+g_0^k[m] ) ), then update (amk+1)(a_m^k+1) by (29). Output: retrieved empirical measure ξK=∑mamKδxmKξ^K= _ma_m^K _x_m^K. All steps above are deterministic given the Sinkhorn solver (which itself is deterministic for fixed initialization and tolerance). The support update is a pushforward by a deterministic barycentric map, and the weight update is a deterministic multiplicative reweighting on the simplex. Thus the overall retrieval operator is deterministic. Appendix I Auxiliary results I.1 Dual and optimal potentials A standard dual form of OTεOT_ can be written (up to equivalent normalizations) in terms of potentials f,g∈C(Ω)f,g∈ C( ).One defines the entropic soft c-transform operator AεA_ via an expression of the form Aε(g,ν)(x)≔−εlog∫Ωexp(gμ,ν(y)−c(x,y)ε)ν(y)(defined up to an additive constant).A_ (g,ν)(x) - _ \! ( g_μ,ν(y)-c(x,y) )\,dν(y) (defined up to an additive constant). Then the optimal potentials (fμ,ν,gμ,ν)(f_μ,ν,g_μ,ν) (Schrödinger potentials) satisfy the Schrödinger system fμ,ν=Aε(gμ,ν,ν)μ-a.e.,gμ,ν=Aε(fμ,ν,μ)ν-a.e..f_μ,ν=A_ (g_μ,ν,ν) μ-a.e., g_μ,ν=A_ (f_μ,ν,μ) ν-a.e.. These potentials are unique up to adding a constant to fμ,νf_μ,ν and subtracting the same constant from gμ,νg_μ,ν (gauge invariance). This does not affect any gradient ∇fμ,ν(x)∇ f_μ,ν(x) or ∇gμ,ν(x)∇ g_μ,ν(x), which is what we ultimately use. I.2 Barycentric projection map Let πμ,νεπ _μ,ν be the optimal entropic coupling between μ and ν, which can be disintegrated into a marginal and conditional distribution as follows : πμ,νε(dx,dy)=μ(dx)πμ,νε(dy∣x).π _μ,ν(dx,dy)=μ(dx)\,π _μ,ν(dy x). Define the barycentric projection (conditional mean) map: Tμ→νε(x)≔∫Ωyπμ,νε(dy∣x).T _μ→ν(x) _ y\,π _μ,ν(dy x). This is defined μ-a.e. and takes values in conv(Ω)⊂ℝdconv( ) ^d. I.3 Computing gradient of Schrödinger potentials explicitly for quadratic costs Recall (from the Schrödinger system) that fμ,ν(x)=−εlog∫Ωexp(gμ,ν(y)−c(x,y)ε)ν(y)(up to constant).f_μ,ν(x)=- _ ( g_μ,ν(y)-c(x,y) )\,dν(y) (up to constant). Differentiate with respect to x. Denote Z(x)≔∫Ωexp(gμ.ν(y)−c(x,y)ε)ν(y).Z(x) _ ( g_μ.ν(y)-c(x,y) )\,dν(y). Then fμ,ν(x)=−εlogZ(x)f_μ,ν(x)=- Z(x), so ∇fμ,ν(x)=−ε1Z(x)∇Z(x).∇ f_μ,ν(x)=- 1Z(x)∇ Z(x). Compute ∇Z(x)∇ Z(x): ∇Z(x)=∫Ωexp(gμ,ν(y)−c(x,y)ε)⋅1ε(−∇c(x,y))ν(y).∇ Z(x)= _ ( g_μ,ν(y)-c(x,y) )· 1 (-∇ c(x,y) )\,dν(y). Therefore ∇fμ,ν(x)= ∇ f_μ,ν(x)= −ε1Z(x)∫exp(gμ,ν(y)−c(x,y)ε)⋅1ε(−∇c(x,y))ν(y) - 1Z(x) ( g_μ,ν(y)-c(x,y) )· 1 (-∇ c(x,y) )\,dν(y) = = 1Z(x)∫exp(gμ,ν(y)−c(x,y)ε)∇c(x,y)ν(y). 1Z(x) ( g_μ,ν(y)-c(x,y) )\,∇ c(x,y)\,dν(y). But the conditional distribution πε(dy∣x)π (dy x) has density proportional to exp((gμ,ν(y)−c(x,y))/ε),dν(y) ((g_μ,ν(y)-c(x,y))/ ),dν(y). Hence ∇fμ,ν(x)=∫∇c(x,y)πμ,νε(dy∣x).∇ f_μ,ν(x)= ∇ c(x,y)\,π _μ,ν(dy x). Now plug c(x,y)=12‖x−y‖22c(x,y)= 12\|x-y\|_2^2, so ∇c(x,y)=x−y∇ c(x,y)=x-y. Then ∇fμ,ν(x)=∫(x−y)πμ,νε(dy∣x)=x−∫yπμ,νε(dy∣x)=x−Tμ→νε(x).∇ f_μ,ν(x)= (x-y)\,π _μ,ν(dy x)=x- y\,π _μ,ν(dy x)=x-T _μ→ν(x). Thus we have the fundamental identity: ∇fμ,ν(x)=x−Tμ→νε(x)for c(x,y)=12‖x−y‖22.∇ f_μ,ν(x)=x-T _μ→ν(x) c(x,y)= 12\|x-y\|_2^2. This is precisely why barycentric projections give a transport-map-like representation of entropic OT gradients. Similarly, we have that ∇gμ,ν(y)=y−Tν→με(y)∇ g_μ,ν(y)=y-T_ν→μ^ (y). In particular, we have that 12(∇fξ,ξ(x)+∇gξ,ξ(x))=x−Tξ→ξε(x). 12 (∇ f_ξ,ξ(x)+∇ g_ξ,ξ(x) )=x-T_ξ→ξ^ (x). (30) I.4 Additional lemmas Lemma 3 (One step descent of softmin energy). Let Assumptions (A1) and (A2) hold. Let ξ+=Φη(ξ)=Retξ(−ηgradSHKE(ξ))ξ^+= _η(ξ)=Ret_ξ(- _SHKE(ξ)). Then E(ξ+)≤E(ξ)−η(1−Lη2)‖gradSHKE(ξ)‖SHK,ξ2.E (ξ^+ )≤ E(ξ)-η (1- Lη2 )\|grad_SHKE(ξ)\|_SHK,ξ^2. In particular, if 0<η≤1L0<η≤ 1L, then E(ξ+)≤E(ξ)−η2‖gradSHKE(ξ)‖SHK,ξ2,E (ξ^+ )≤ E(ξ)- η2\|grad_SHKE(ξ)\|_SHK,ξ^2, so E strictly decreases unless gradSHKE(ξ)=(0,0)grad_SHKE(ξ)=(0,0). Proof. Using Equation 21, with (r,v)=−gradSHKE(ξ)(r,v)=-grad_SHKE(ξ), we have that E(ξ+)= E (ξ^+ )= E(Retξ(−gradSHKE(ξ)) E (Ret_ξ(-grad_SHKE(ξ) ) ≤ ≤ E(ξ)−η⟨gradSHKE(ξ),gradSHKE(ξ)⟩SHK,ξ+Lη22‖gradSHKE(ξ)‖SHK,ξ2 E(ξ)-η _SHKE(ξ),grad_SHKE(ξ) _SHK,ξ+ Lη^22\|grad_SHKE(ξ)\|_SHK,ξ^2 = = E(ξ)−η(1−Lη2)‖gradSHKE(ξ)‖SHK,ξ2. E(ξ)-η (1- Lη2 )\|grad_SHKE(ξ)\|_SHK,ξ^2. If η≤1Lη≤ 1L, then 1−Lη2≥121- Lη2≥ 12, giving the stated bound. ∎ Lemma 4 (Energy gap contraction). Let Assumptions (A1), (A2) and (A4) hold and choose 0<η≤min1L,1μ0<η≤ \ 1L, 1μ \. Then for i=1,…,Ni=1,…,N, if ξ(0)∈Bi(r)ξ^(0)∈ B_i(r), ξ(k+1)=Φη(ξ(k))ξ^(k+1)= _η (ξ^(k) ) and ξ(k)∈Bi(r)ξ^(k)∈ B_i(r) for k∈ℕk , we have E(ξ(k+1))−Ei∗(r)≤(1−μη)(E(ξ(k))−Ei∗(r)).E (ξ^(k+1) )-E_i^*(r)≤ (1-μη ) (E (ξ^(k) )-E_i^*(r) ). Hence, for any k∈ℕk , E(ξ(k))−Ei∗(r)≤(1−μη)k(E(ξ(0))−Ei∗(r))E (ξ^(k) )-E_i^*(r)≤ (1-μη )^k (E (ξ^(0) )-E_i^*(r) ) and consequently, limk→∞E(ξ(k))=Ei∗(r) _k→∞E(ξ^(k))=E_i^*(r). Further, for each k∈ℕk , ‖gradSHKE(ξ(k))‖SHK,ξ(k)2≤2η(E(ξ(k))−E(ξ(k+1))). \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)^2≤ 2η (E (ξ^(k) )-E (ξ^(k+1) ) ). Consequently, for each k∈ℕk , ‖gradSHKE(ξ(k))‖SHK,ξ(k)≤2η(E(ξ(k))−Ei∗(r))≤2(E(ξ(0))−Ei∗(r))η(1−ημ)k2. \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)≤ 2η (E (ξ^(k) )-E_i^*(r) )≤ 2 (E (ξ^(0) )-E_i^*(r) )η(1-ημ) k2. Proof. Using Lemma 3, we have that, provided ξ(k)∈Bi(r)ξ^(k)∈ B_i(r), E(ξ(k+1))≤E(ξ(k))−η2‖gradSHKE(ξ(k))‖SHK,ξ(k)2.E (ξ^(k+1) )≤ E (ξ^(k) )- η2 \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)^2. Subtract Ei∗(r)E_i^*(r) from both sides: E(ξ(k+1))−Ei∗(r)≤E(ξ(k))−Ei∗(r)−η2‖gradSHKE(ξ(k))‖SHK,ξ(k)2.E (ξ^(k+1) )-E_i^*(r)≤ E (ξ^(k) )-E_i^*(r)- η2 \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)^2. Using Equation 22, we have that 12‖gradSHKE(ξ(k))‖SHK,ξ(k)2 12 \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)^2 ≥μ(E(ξ(k))−Ei∗(r)) ≥μ (E (ξ^(k) )-E_i^*(r) ) ⟹‖gradSHKE(ξ(k))‖SHK,ξ(k)2 \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)^2 ≥2μ(E(ξ(k))−Ei∗(r)). ≥ 2μ (E (ξ^(k) )-E_i^*(r) ). Thus, provided ξ(k)∈Bi(r)ξ^(k)∈ B_i(r), E(ξ(k+1))−Ei∗(r)≤E(ξ(k))−Ei∗(r)−ημ(E(ξ(k))−Ei∗(r))=(1−μη)(E(ξ(k))−Ei∗(r)).E (ξ^(k+1) )-E_i^*(r)≤ E (ξ^(k) )-E_i^*(r)-ημ (E (ξ^(k) )-E_i^*(r) )= (1-μη ) (E (ξ^(k) )-E_i^*(r) ). Iterating the inequality, we obtain 0≤E(ξ(k))−Ei∗(r)≤(1−μη)k(E(ξ(0))−Ei∗(r))0≤ E (ξ^(k) )-E_i^*(r)≤ (1-μη )^k (E (ξ^(0) )-E_i^*(r) ) and consequently limk→∞E(ξ(k))=Ei∗(r) _k→∞E(ξ^(k))=E_i^*(r). Using Lemma 3, since 0<η≤1L0<η≤ 1L, we have that E(ξ(k+1))≤E(ξ(k))−η2‖gradSHKE(ξ(k))‖SHK,ξ(k)2E (ξ^(k+1) )≤ E (ξ^(k) )- η2 \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)^2 which upon rearrangement gives ‖gradSHKE(ξ(k))‖SHK,ξ(k)2≤2η(E(ξ(k))−E(ξ(k+1))). \|grad_SHKE (ξ^(k) ) \|_SHK,ξ^(k)^2≤ 2η (E (ξ^(k) )-E (ξ^(k+1) ) ). Again, since we assume that ξ(k+1)∈Bi(r)ξ^(k+1)∈ B_i(r), we have that Ei∗(r)≤E(ξ(k+1))E_i^*(r)≤ E (ξ^(k+1) ). Consequently, E(ξ(k))−E(ξ(k+1))≤E(ξ(k))−Ei∗(r).E (ξ^(k) )-E (ξ^(k+1) )≤ E (ξ^(k) )-E_i^*(r). This completes the proof. ∎ Lemma 5 (Control on basin interference and softmin perturbation bounds of energy E inside local basin using separation margin). Let Assumption (A1) hold. Define the Gibbs weights wi(ξ)w_i(ξ) corresponding to any fixed ξ∈Bi(r)ξ∈ B_i(r) as in Equation 5. Then, we have that 1. The weights satisfy wi(ξ)≥11+(N−1)e−βΔ,∑j≠iwj(ξ)≤(N−1)e−βΔ1+(N−1)e−βΔ.w_i(ξ)≥ 11+(N-1)e^-β , _j≠ iw_j(ξ)≤ (N-1)e^-β 1+(N-1)e^-β . 2. The gap between Fi(ξ)F_i(ξ) and E(ξ)E(ξ) satisfies 0≤Fi(ξ)−E(ξ)≤1βlog(1+(N−1)e−βΔ)≤N−1βe−βΔ.0≤ F_i(ξ)-E(ξ)≤ 1β (1+(N-1)e^-β )≤ N-1βe^-β . Proof. Using Equation 20, we have that, for ξ∈Bi(r)ξ∈ B_i(r), ∑j≠ie−βFj(ξ)=e−βFi(ξ)×∑j≠ie−β(Fj(ξ)−Fi(ξ))≤e−βFi(ξ)×∑j≠ie−βΔr≤e−βFi(ξ)×(N−1)e−βΔ. _j≠ ie^-β F_j(ξ)=e^-β F_i(ξ)× _j≠ ie^-β (F_j(ξ)-F_i(ξ) )≤ e^-β F_i(ξ)× _j≠ ie^-β _r≤ e^-β F_i(ξ)×(N-1)e^-β . Therefore, we have that wi(ξ)=e−βFi(ξ)e−βFi(ξ)+∑j≠ie−βFj(ξ)≥11+(N−1)e−βΔw_i(ξ)= e^-β F_i(ξ)e^-β F_i(ξ)+ _j≠ ie^-β F_j(ξ)≥ 11+(N-1)e^-β Therefore, we have that ∑j≠iwj(ξ)=1−wi(ξ)≤(N−1)e−βΔ1+(N−1)e−βΔ. _j≠ iw_j(ξ)=1-w_i(ξ)≤ (N-1)e^-β 1+(N-1)e^-β . Now, we have that E(ξ)=−1βlog(e−βFi(ξ)[1+∑j≠ie−β(Fj(ξ)−Fi(ξ))])=Fi(ξ)−1βlog(1+∑j≠ie−β(Fj(ξ)−Fi(ξ))).E(ξ)=- 1β (e^-β F_i(ξ) [1+ _j≠ ie^-β (F_j(ξ)-F_i(ξ) ) ] )=F_i(ξ)- 1β (1+ _j≠ ie^-β (F_j(ξ)-F_i(ξ) ) ). Since each Fj(ξ)−Fi(ξ)≥ΔF_j(ξ)-F_i(ξ)≥ , 0≤1βlog(1+∑j≠ie−β(Fj−Fi))≤1βlog(1+(N−1)e−βΔ).0≤ 1β (1+ _j≠ ie^-β (F_j-F_i ) )≤ 1β (1+(N-1)e^-β ). Finally, we use log(1+u)≤u (1+u)≤ u to get the last bound. ∎ Lemma 6 (Expression of SHK gradient oof energy functional E in terms of SHK gradients of Sinkhorn divergence). gradSHKE(ξ)=∑j=1Nwj(ξ)gradSHKFj(ξ)grad_SHKE(ξ)= _j=1^Nw_j(ξ)grad_SHKF_j(ξ) Proof. Proof is obvious. ∎ Lemma 7 (Bounding SHK distance in terms of SHK gradient norm along retraction curve). Fix any (a,x)(a,x) and any tangent increment (δa,δx)∈T(a,x)ℳM(δ a,δ x)∈ T_(a,x)M_M. Define the retraction curve γ(t)≔Ret(a,x)(tδa,tδx),t∈[0,1].γ(t) _(a,x)(tδ a,tδ x), t∈[0,1]. Let dSHKd_SHK denote the Riemannian distance induced by the SHK metric on ℳM=ΔM,amin∘×LocM,ΔminM_M= _M,a_ ×Loc_M, _ . Then dSHK((a,x),Ret(a,x)(δa,δx))≤Length(γ)≤mine‖δa‖∞,1amin‖(δa,δx)‖SHK,(a,x).d_SHK ((a,x),Ret_(a,x)(δ a,δ x) ) (γ)≤ \e^\| δ aa\|_∞, 1 a_ \\|(δ a,δ x)\|_SHK,(a,x). where Length(γ)≔∫01‖dtγ(t)‖SHK,γ(t)tLength(γ) _0^1\| ddtγ(t)\|_SHK,γ(t)dt and (δa)m≔δamam( δ aa)_m δ a_ma_m. Proof. By definition of the Riemannian distance as the infimum of lengths over all curves connecting the points, dSHK(γ(0),γ(1))≤Length(γ).d_SHK(γ(0),γ(1)) (γ). So it suffices to bound the length of γ. We have that γ(t)=(a(t),x(t))γ(t)=(a(t),x(t)) with x(t)=Retxpos(tδx)=x+tδx,a(t)=Retaw(tδa).x(t)=Ret_x^pos(tδ x)=x+tδ x, a(t)=Ret_a^w(tδ a). Then, x′(t)≔dtx(t)=δx (t) ddtx(t)=δ x. For the weights, introduce sm≔δamams_m δ a_ma_m Then, s∈ℝMs ^M and ∑mamsm=0 _ma_ms_m=0 since ∑mδam=0 _mδ a_m=0. By the definition of RetawRet_a^w, we have that am(t)=ametsmZ(t),Z(t)≔∑ℓ=1Maℓetsℓ.a_m(t)= a_me^ts_mZ(t), Z(t) _ =1^Ma_ e^ts_ . Differentiating, we have that dtlogam(t)=sm−Z′(t)Z(t),Z′(t)Z(t)=∑ℓ=1Maℓ(t)sℓ=:s¯(t). ddt a_m(t)=s_m- Z (t)Z(t), Z (t)Z(t)= _ =1^Ma_ (t)s_ =: s(t). Thus am′(t)≔dtam(t)=am(t)(sm−s¯(t))a_m (t) ddta_m(t)=a_m(t) (s_m- s(t) ) The SHK gradient norm at γ(t)=(a(t),x(t))γ(t)=(a(t),x(t)) is given by ‖(u,v)‖SHK,(a(t),x(t))2=∑m=1M(λ2am(t)um2+am(t)‖vm‖2).\|(u,v)\|_SHK,(a(t),x(t))^2= _m=1^M ( λ^2a_m(t)u_m^2+a_m(t) \|v_m \|^2 ). Applying it to (u,v)=dtγ(t)≔dtγ(t)=(a′(t),x′(t))(u,v)= ddtγ(t) ddtγ(t)= (a (t),x (t) ), we have that ‖dtγ(t)‖SHK,(a(t),x(t))2= \| ddtγ(t) \|_SHK,(a(t),x(t))^2= ∑m=1Mλ2am(t)(am′(t))2+∑m=1Mam(t)‖xm′(t)‖2 _m=1^M λ^2a_m(t) (a_m (t) )^2+ _m=1^Ma_m(t) \|x_m (t) \|^2 = = λ2∑m=1Mam(t)(sm−s¯(t))2+∑m=1Mam(t)‖δxm‖2 λ^2 _m=1^Ma_m(t) (s_m- s(t) )^2+ _m=1^Ma_m(t) \|δ x_m \|^2 = = λ2[∑m=1Mam(t)sm2−2s¯(t)∑m=1Mam(t)sm+(s¯(t))2]+∑m=1Mam(t)‖δxm‖2 λ^2 [ _m=1^Ma_m(t)s_m^2-2 s(t) _m=1^Ma_m(t)s_m+ ( s(t) )^2 ]+ _m=1^Ma_m(t) \|δ x_m \|^2 = = λ2[∑m=1Mam(t)sm2−(s¯(t))2]+∑m=1Mam(t)‖δxm‖2 λ^2 [ _m=1^Ma_m(t)s_m^2- ( s(t) )^2 ]+ _m=1^Ma_m(t) \|δ x_m \|^2 ≤ ≤ λ2∑m=1Mam(t)sm2+∑m=1Mam(t)‖δxm‖2. λ^2 _m=1^Ma_m(t)s_m^2+ _m=1^Ma_m(t) \|δ x_m \|^2. Let S≔‖s‖∞=‖δa/a‖∞S \|s\|_∞=\|δ a/a\|_∞. Then for all t∈[0,1]t∈[0,1], e−S≤etsm≤eS,e−S≤Z(t)=∑ℓaℓetsℓ≤eS.e^-S≤ e^ts_m≤ e^S, e^-S≤ Z(t)= _ a_ e^ts_ ≤ e^S. Hence ame−SeS≤am(t)=ametsmZ(t)≤ameSe−S, a_me^-Se^S≤ a_m(t)= a_me^ts_mZ(t)≤ a_me^Se^-S, i.e. ame−2S≤am(t)≤ame2S.a_me^-2S≤ a_m(t)≤ a_me^2S. Using am(t)≤ame2Sa_m(t)≤ a_me^2S, ‖dtγ(t)‖SHK,(a(t),x(t))2≤e2S(λ2∑m=1Mamsm2+∑m=1Mam‖δxm‖2)=e2S‖(δa,δx)‖SHK,(a,x)2. \| ddtγ(t) \|_SHK,(a(t),x(t))^2≤ e^2S (λ^2 _m=1^Ma_ms_m^2+ _m=1^Ma_m \|δ x_m \|^2 )=e^2S\|(δ a,δ x)\|_SHK,(a,x)^2. Taking square roots gives, for all t∈[0,1]t∈[0,1], ‖dtγ(t)‖(a(t),x(t))≤eS‖(δa,δx)‖SHK,(a,x). \| ddtγ(t) \|_(a(t),x(t))≤ e^S\|(δ a,δ x)\|_SHK,(a,x). Therefore, Length(γ)=∫01‖dtγ(t)‖SHK,(a(t),x(t))t≤∫01eS‖(δa,δx)‖SHK,(a,x)t=eS‖(δa,δx)‖SHK,(a,x).Length(γ)= _0^1 \| ddtγ(t) \|_SHK,(a(t),x(t))dt≤ _0^1e^S\|(δ a,δ x)\|_SHK,(a,x)dt=e^S\|(δ a,δ x)\|_SHK,(a,x). Now, note that. ∑m=1Mam(t)‖δxm‖2≤maxm‖δxm‖2≤1amin∑m=1Mam‖δxm‖2 _m=1^Ma_m(t) \|δ x_m \|^2≤ _m \|δ x_m \|^2≤ 1a_ _m=1^Ma_m \|δ x_m \|^2 and λ2∑m=1Mam(t)(sm−s¯(t))2≤ λ^2 _m=1^Ma_m(t) (s_m- s(t) )^2≤ λ2∑m=1Mam(t)(sm−s¯(0))2 λ^2 _m=1^Ma_m(t) (s_m- s(0) )^2 ≤ ≤ λ2maxm(sm−s¯(0))2 λ^2 _m (s_m- s(0) )^2 ≤ ≤ λ2amin∑m=1Mam(sm−s¯(0))2=λ2amin∑m=1M(δam)2am. λ^2a_ _m=1^Ma_m (s_m- s(0) )^2= λ^2a_ _m=1^M (δ a_m)^2a_m. Combining, we have that, ‖dtγ(t)‖SHK,(a(t),x(t))≤1amin‖(δa,δx)‖SHK,(a,x) \| ddtγ(t) \|_SHK,(a(t),x(t))≤ 1 a_ \|(δ a,δ x)\|_SHK,(a,x) and hence, Length(γ)≤1amin‖(δa,δx)‖SHK,(a,x).Length(γ)≤ 1 a_ \|(δ a,δ x)\|_SHK,(a,x). ∎ Lemma 8 (Open bounded convex domain implies boundedness of SHK gradient of Sinkhorn divergence). Assume that the domain Ω is open, bounded and convex with diameter D≔supx,y∈Ω‖x−y‖<∞D _x,y∈ \|x-y\|<∞. Then, for ξ,ν∈M,amin,Δmin(Ω)ξ,ν _M,a_ , _ ( ), we have that ‖gradSHKSε(ξ,ν)‖SHK,ξ≤G.\|grad_SHKS_ (ξ,ν)\|_SHK,ξ≤ G. where G≔D1+D2λ2G D 1+ D^2λ^2. In particular, for i=1,…,Ni=1,…,N, we consequently obtain ‖gradSHKFi(ξ)‖SHK,ξ≤G.\|grad_SHKF_i(ξ)\|_SHK,ξ≤ G. Further, for any fixed ξ,ξ′,ν∈M,amin,Δmin(Ω)ξ,ξ ,ν _M,a_ , _ ( ), the following holds |Sε(ξ,ν)−Sε(ξ′,ν)|≤GdSHK(ξ,ξ′) |S_ (ξ,ν)-S_ (ξ ,ν) |≤ Gd_SHK (ξ,ξ ) where dSHKd_SHK is the SHK metric restricted to M,amin,Δmin(Ω)P_M,a_ , _ ( ). Therefore, Sϵ(⋅,ν)S_ε(·,ν) is G-Lipschitz continuous in the SHK metric. Proof. Note that, for any fixed ν∈M,amin,Δmin(Ω)ν _M,a_ , _ ( ) and ξ∈M,amin,Δmin(Ω)ξ _M,a_ , _ ( ), we have that ϕ(x)≔(δSε(ξ,ν)δξ)(x)=fξ,ν(x)−12(fξ,ξ(x)+gξ,ξ(x))φ(x) ( δ S_ (ξ,ν)δξ )(x)=f_ξ,ν(x)- 12 (f_ξ,ξ(x)+g_ξ,ξ(x) ) where wi(ξ)w_i(ξ) are defined as in Equation 5. Note that gradSHKSε(ξ,ν)(x)=(r(x),v(x))grad_SHKS_ (ξ,ν)(x)=(r(x),v(x)) where v(x)=∇ϕ(x)v(x)=∇φ(x) and r(x)=1λ2(ϕ(x)−∫ϕ(y)ξ(y))r(x)= 1λ^2 (φ(x)- φ(y)dξ(y) ). Then, ‖gradSHKSε(ξ,ν)‖SHK,ξ2=∫Ω‖v(x)‖2ξ(x)+λ2∫Ω[r(x)]2ξ(x). \|grad_SHKS_ (ξ,ν) \|_SHK,ξ^2= _ \|v(x)\|^2dξ(x)+λ^2 _ [r(x) ]^2dξ(x). Note that, v(x)=∇(fξ,ν(x)−12(fξ,ξ(x)+gξ,ξ(x)))=−[Tξ→νϵ(x)−Tξ→ξϵ(x)].v(x)=∇ (f_ξ,ν(x)- 12 (f_ξ,ξ(x)+g_ξ,ξ(x) ) )=- [T_ξ→ν^ε(x)-T_ξ→ξ^ε(x) ]. Each barycentric map Tξ→νϵ(x)T_ξ→ν^ε(x) lies in Ω , since it is a conditional expectation defined over the domain Ω , and similarly Tξ→ξϵ(x)∈ΩT_ξ→ξ^ε(x)∈ as well. Hence ‖v(x)‖=‖∇ϕ(x)‖=‖Tξ→νϵ(x)−Tξ→ξϵ(x)‖≤D.\|v(x)\|=\|∇φ(x)\|= \|T_ξ→ν^ε(x)-T_ξ→ξ^ε(x) \|≤ D. (31) Since Ω is convex and open, every line segment joining two points of Ω stays in Ω . By the mean value theorem and the bound ‖∇ϕ(x)‖≤D\|∇φ(x)\|≤ D, for any x,y∈Ωx,y∈ , |ϕ(x)−ϕ(y)|≤supz∈[x,y]‖∇ϕ(z)‖2‖x−y‖2≤D‖x−y‖2≤D2.|φ(x)-φ(y)|≤ _z∈[x,y]\|∇φ(z)\|_2\|x-y\|_2≤ D\|x-y\|_2≤ D^2. (32) Therefore, we have that, |r(x)|=1λ2|∫[ϕ(x)−ϕ(y)]ξ(y)|≤1λ2∫|ϕ(x)−ϕ(y)|ξ(y)≤D2λ2.|r(x)|= 1λ^2 | [φ(x)-φ(y) ]dξ(y) |≤ 1λ^2 |φ(x)-φ(y) |dξ(y)≤ D^2λ^2. Therefore, we have that ‖gradSHKSε(ξ,ν)‖SHK,ξ2≤D2+λ2×D4λ4≤D2+D4λ2. \|grad_SHKS_ (ξ,ν) \|_SHK,ξ^2≤ D^2+λ^2× D^4λ^4≤ D^2+ D^4λ^2. In particular, we can choose ν=Xiν=X_i. For the proof of the Lipschitz continuity of Fν(⋅)≔Sε(⋅,ν)F_ν(·) S_ (·,ν), consider γ:[0,1]→M,amin,Δmin(Ω)γ:[0,1] _M,a_ , _ ( ) to be an absolutely continuous path in M,amin,Δmin(Ω)P_M,a_ , _ ( ) between ξ and ξ′ξ . Then, by the chain rule, dtFν(γ(t))=⟨gradFν(γ(t)),dtγ(t)⟩SHK,γ(t) ddtF_ν(γ(t))= _ν(γ(t)), ddtγ(t) _SHK,γ(t) Apply Cauchy-Schwarz and the gradient bound: |dtFν(γ(t))|≤‖gradSHKFν(γ(t))‖SHK,γ(t)‖dtγ(t)‖SHK,γ(t)≤G‖dtγ(t)‖SHK,γ(t) | ddtF_ν(γ(t)) |≤ \|grad_SHKF_ν(γ(t)) \|_SHK,γ(t) \| ddtγ(t) \|_SHK,γ(t)≤ G \| ddtγ(t) \|_SHK,γ(t) Integrating from 0 to 1, we have that, |Fν(ξ′)−Fν(ξ)|≤G∫01‖dtγ(t)‖SHK,γ(t)t=G Length (γ) |F_ν (ξ )-F_ν(ξ) |≤ G _0^1 \| ddtγ(t) \|_SHK,γ(t)dt=G Length (γ) Taking the infimum over all curves γ from ξ to ξ′ξ gives |Fν(ξ′)−Fν(ξ)|≤GdSHK(ξ,ξ′) |F_ν (ξ )-F_ν(ξ) |≤ Gd_SHK (ξ,ξ ) ∎ Lemma 9 (Softmin function sensitivity). Define softminβ(z1,…,zN)≔−1βlog(∑i=1Nexp(−βzi))softmin_β (z_1,…,z_N ) - 1β ( _i=1^N (-β z_i ) ). Then, we have that |softminβ(z)−softminβ(z′)|≤‖z−z′‖∞∀z,z′∈ℝN. |softmin_β(z)-softmin_β (z ) |≤ \|z-z \|_∞ ∀ z,z ^N. Proof. Let us define Mdiff≔max1≤i≤N(zi−zi′)M_diff _1≤ i≤ N (z_i-z_i ). Then, we have that zi≤zi′+Mdiffz_i≤ z_i +M_diff for i=1,…,Ni=1,…,N. Note that that ∂zisoftminβ(z)=exp(−βzi)∑j=1Nexp(−βzj)>0 ∂ z_isoftmin_β(z)= (-β z_i) _j=1^N (-β z_j)>0. Since every partial derivative is positive, the softminsoftmin function is increasing in each variable ziz_i. Therefore, if two vectors z,w∈ℝNz,w ^N satisfy zi≤wiz_i≤ w_i for i=1,…i=1,…, then softminβ(z)≤softminβ(w).softmin_β(z) _β(w). Further, for any vector w∈ℝNw ^N, by direct computation, we have that softminβ(w+C)=softminβ(w1+C,…,wN+C)=softminβ(w)+Csoftmin_β(w+C1)=softmin_β(w_1+C,…,w_N+C)=softmin_β(w)+C, which show the translation-equivariant nature of softminsoftmin. Therefore, we have that softminβ(z)≤softminβ(Z′+Mdiff)=softminβ(z)+Mdiff.softmin_β(z) _β(Z +M_diff1)=softmin_β(z)+M_diff. Consequently, we have that, ‖softminβ(z)−softminβ(z′)‖∞≤Mdiff≤‖z−z′‖∞.\|softmin_β(z)-softmin_β(z )\|_∞≤ M_diff≤\|z-z \|_∞. ∎ Lemma 10 (Lipschitz continuity of LSE functional with respect to SHK metric). Assume that the domain Ω is open, bounded and convex with diameter D≔supx,y∈Ω‖x−y‖<∞D _x,y∈ \|x-y\|<∞. Then, for stored patterns X1,…,XN∈M,amin,Δmin(Ω)X_1,…,X_N _M,a_ , _ ( ) and Fi≔Sε(⋅,Xi)F_i S_ (·,X_i), we have that |Fi(ξ)−Fi(ξ′)|≤GdSHK(ξ,ξ′) |F_i(ξ)-F_i (ξ ) |≤ Gd_SHK (ξ,ξ ) and |E(ξ)−E(ξ′)|≤GdSHK(ξ,ξ′)∀ξ,ξ′∈PM(Ω). |E(ξ)-E (ξ ) |≤ Gd_SHK (ξ,ξ ) ∀ξ,ξ ∈ P_M( ). where G≔D1+D2λ2G D 1+ D^2λ^2. Therefore,FiF_i and E are both Lipschitz continuous with respect to the SHK metric with Lipschitz constant G, and therefore are continuous as well. Proof. With Fi=Sε(⋅,Xi)F_i=S_ (·,X_i), applying Lemma 9 with zi=Fi(ξ)z_i=F_i(ξ) and zi′=Fi(ξ′)z_i =F_i(ξ ), we have that |E(ξ)−E(ξ′)|≤max1≤i≤N|Fi(ξ)−Fi(ξ′)|. |E(ξ)-E (ξ ) |≤ _1≤ i≤ N |F_i(ξ)-F_i (ξ ) |. Again, by Lemma 8, we have, for i=1,…,Ni=1,…,N, |Fi(ξ)−Fi(ξ′)|≤GdSHK(ξ,ξ′). |F_i(ξ)-F_i (ξ ) |≤ Gd_SHK (ξ,ξ ). Taking maximum over i and combining with the previous inequality, we have that |E(ξ)−E(ξ′)|≤GdSHK(ξ,ξ′). |E(ξ)-E (ξ ) |≤ Gd_SHK (ξ,ξ ). ∎ Lemma 11 (Mean separation of patterns implies margin separation of Sinkhorn divergences). Let X1,…,XN∈M,amin,Δmin(Ω)X_1,…,X_N _M,a_ , _ ( ) and define μi≔m(Xi)∈ℝd _i m (X_i ) ^d. Assume the means are pairwise separated: ‖μi−μj‖≥dmin for all i≠j. \| _i- _j \|≥ d_ for all i≠ j. Fix any ε>0 >0 such that dmin2>32εlogMd_ ^2>32 M and define r≔dmin232−εlogM,Δ≔dmin24r d_ ^232- M, d_ ^24 Then Assumption (A1) holds with the given choice of r and Δ ; i.e. for every i, for every ξ∈Bi(r)ξ∈ B_i(r), and every j≠ij≠ i, Fj(ξ)−Fi(ξ)≥Δ.F_j(ξ)-F_i(ξ)≥ . Proof. Fix an index i and let ξ∈Bi(r)ξ∈ B_i(r), so Sε(ξ,Xi)≤rS_ (ξ,X_i )≤ r. Using Lemma 14 with (μ,ν)=(ξ,Xi)(μ,ν)= (ξ,X_i ), we have that: Sε(ξ,Xi)≥12‖m(ξ)−μi‖2−εlogM.S_ (ξ,X_i )≥ 12 \|m(ξ)- _i \|^2- M. Combining with Sε(ξ,Xi)≤rS_ (ξ,X_i )≤ r, we have that r≥12‖m(ξ)−μi‖2−εlogM⟹12‖m(ξ)−μi‖2≤r+εlogM.r≥ 12 \|m(ξ)- _i \|^2- M 12 \|m(ξ)- _i \|^2≤ r+ M. By the definition of r, we have that r+εlogM=dmin232r+ M= d_ ^232. Hence ‖m(ξ)−μi‖2≤2⋅dmin232=dmin216⟹‖m(ξ)−μi‖≤dmin4. \|m(ξ)- _i \|^2≤ 2· d_ ^232= d_ ^216 \|m(ξ)- _i \|≤ d_ 4. Fix j≠ij≠ i. Then, triangle inequality gives ‖m(ξ)−μj‖≥‖μi−μj‖−‖m(ξ)−μi‖≥dmin−dmin4=3dmin4. \|m(ξ)- _j \|≥ \| _i- _j \|- \|m(ξ)- _i \|≥ d_ - d_ 4= 3d_ 4. Using this along wih Lemma 14, for (ξ,Xj) (ξ,X_j ), we have that Sε(ξ,Xj)≥12‖m(ξ)−μj‖2−εlogM≥12(3dmin4)2−εlogM=9dmin232−εlogM.S_ (ξ,X_j )≥ 12 \|m(ξ)- _j \|^2- M≥ 12 ( 3d_ 4 )^2- M= 9d_ ^232- M. Since ξ∈Bi(r),Sε(ξ,Xi)≤r=dmin232−εlogMξ∈ B_i(r),S_ (ξ,X_i )≤ r= d_ ^232- M. Hence , we have that Sε(ξ,Xj)−Sε(ξ,Xi)≥(9dmin232−εlogM)−(dmin232−εlogM)=8dmin232=dmin24=Δ.S_ (ξ,X_j )-S_ (ξ,X_i )≥ ( 9d_ ^232- M )- ( d_ ^232- M )= 8d_ ^232= d_ ^24= . Thus, Assumption (A1) holds under the given conditions. ∎ Lemma 12 (Entropic OT lower bound in terms of mean differences). For any μ,ν∈(Ω)μ,ν ( ) and any ε>0 >0, OTε(μ,ν)≥12‖m(μ)−m(ν)‖2OT_ (μ,ν)≥ 12\|m(μ)-m(ν)\|^2 where m(μ)≔∫Ωxμ(x)∈conv(Ω)⊂ℝdm(μ) _ xdμ(x) ( ) ^d is the mean of μ. Proof. Given any coupling π∈Π(μ,ν)π∈ (μ,ν) and random pair (X,Y)∼π(X,Y) π, define Z≔X−YZ X-Y. Then, we have that π[Z]=π[X]−π[Y]=m(μ)−m(ν).E_π[Z]=E_π[X]-E_π[Y]=m(μ)-m(ν). By Jensen’s inequality, we have that π[‖Z‖2]≥‖π[Z]‖2=‖m(μ)−m(ν)‖2.E_π [\|Z\|^2 ]≥ \|E_π[Z] \|^2=\|m(μ)-m(ν)\|^2. Consequently, using the fact KL(π∣μ⊗ν)≥0KL(π μ ν)≥ 0, we have that ∫Ω×Ωc(x,y)dπ(x,y)+εKL(π∣μ⊗ν)≥ _ × c(x,y)dπ(x,y)+ (π μ ν)≥ ∫Ω×Ωc(x,y)π(x,y) _ × c(x,y)dπ(x,y) = = 12π[‖X−Y‖2] 12E_π [\|X-Y\|^2 ] ≥ ≥ 12‖m(μ)−m(ν)‖2. 12\|m(μ)-m(ν)\|^2. Since this is true for any π∈Π(μ,ν)π∈ (μ,ν), taking the infimum over π∈Π(μ,ν)π∈ (μ,ν), we have that OTε(μ,ν)≥12‖m(μ)−m(ν)‖2.OT_ (μ,ν)≥ 12\|m(μ)-m(ν)\|^2. ∎ Lemma 13 (Self Entropic OT distance upper bound). If μ=∑m=1Mamδxm∈M,amin,Δmin(Ω)μ= _m=1^Ma_m _x_m _M,a_ , _ ( ) with am>0a_m>0, then OTε(μ,μ)≤εlogM.OT_ (μ,μ)≤ M. Proof. Let πdiag≔∑m=1Mamδ(xm,xm)∈Π(μ,μ) _diag _m=1^Ma_m _ (x_m,x_m )∈ (μ,μ) denote the identity/diagonal coupling and its corresponding transport cost is given by ∫c(x,y)πdiag(x,y)=∑m=1Mam×12‖xm−xm‖2=0. c(x,y)d _diag(x,y)= _m=1^Ma_m× 12 \|x_m-x_m \|^2=0. On the discrete support (xm,xℓ)m,ℓ \ (x_m,x_ ) \_m, , we have that (μ⊗μ)(xm,xℓ)=amaℓandπdiag(xm,xℓ)=am,m=ℓ,0,m≠ℓ.(μ μ) (x_m,x_ )=a_ma_ and _diag (x_m,x_ )= casesa_m,&m= ,\\ 0,&m≠ . cases Consequently, we have that, KL(πdiag∣μ⊗μ)=∑m=1Mamlogamam2=∑m=1Mamlog1am=H(a)KL( _diag μ μ)= _m=1^Ma_m a_ma_m^2= _m=1^Ma_m 1a_m=H(a) where H(a)≔−∑m=1MamlogamH(a) - _m=1^Ma_m a_m is the Shannon entropy of the coupling πdiag _diag corresponding to the probability vector a=(a1,…,am)a=(a_1,…,a_m) or equivalently, that of the discrete distribution ∑m=1Mamδxm _m=1^Ma_m _x_m. Using the strict concavity of x↦logx x and Jensen’s inequality (Polyanskiy and Wu, 2025, Theorem 1.4(b)), we have that H(a)=X∼μlog[1μ(X)]≤log[1μ(X)]=logM.H(a)=E_X μ [ 1μ(X) ]≤ [ 1μ(X) ]= M. Therefore, we have that, OTε(μ,μ)≤∫c(x,y)dπdiag(x,y)+εKL(πdiag∣μ⊗μ)=0+εH(a)≤εlogM.OT_ (μ,μ)≤ c(x,y)d _diag(x,y)+ ( _diag μ μ)=0+ H(a)≤ M. ∎ Lemma 14 (Sinkhorn divergence lower bound in terms of mean differences). For any μ,ν∈M,amin,Δmin(Ω)μ,ν _M,a_ , _ ( ), Sε(μ,ν)≥12‖m(μ)−m(ν)‖2−εlogM.S_ (μ,ν)≥ 12\|m(μ)-m(ν)\|^2- M. Proof. Using Lemma 12 for the first term and Lemma 13 for the last two terms, we have, from the definition of Sε(μ,ν)S_ (μ,ν) : Sε(μ,ν) S_ (μ,ν) = = OTε(μ,ν)−12OTε(μ,μ)−12OTε(ν,ν) _ (μ,ν)- 12OT_ (μ,μ)- 12OT_ (ν,ν) ≥ ≥ 12‖m(μ)−m(ν)‖2−12(εlogM)−12(εlogM) 12\|m(μ)-m(ν)\|^2- 12( M)- 12( M) = = 12‖m(μ)−m(ν)‖2−εlogM. 12\|m(μ)-m(ν)\|^2- M. ∎ Lemma 15 (pairwise separation of random sign vectors). Let si,sj∈±1ds_i,s_j∈\± 1\^d be componentwise independent Rademacher random variables and define μi=c+R0dsi,μj=c+R0dsj. _i=c+ R_0 ds_i, _j=c+ R_0 ds_j. Fix γ∈(0,1)γ∈(0,1). Then ℙ(‖μi−μj‖2<2(1−γ)R02)≤exp(−γ22d).P ( \| _i- _j \|^2<2(1-γ)R_0^2 )≤ (- γ^22d ). Proof. Define the random variable H≔∑k=1dsik≠sjkH _k=1^d1 \s_ik≠ s_jk \. For any k=1,…,dk=1,…,d, since siks_ik and sjks_jk are independent Rademacher random variables, sik≠sjk1 \s_ik≠ s_jk \ is a Bernoulli random variable with probability parameter p=12p= 12. Using the independence across k, we have that H is the sum of d i.i.d Bernoulli random variables with common probability parameter p=12p= 12. Consequently, H∼Binomial(d,12)H (d, 12) and (H)=d2E(H)= d2. Note that |sik−sjk|=0|s_ik-s_jk|=0 if they are equal and 22 is different. Further, μi−μj=R0d(si−sj) _i- _j= R_0 d (s_i-s_j ). Therefore ‖μi−μj‖2=R02d∑k=1d(si,k−sj,k)2=4R02dH \| _i- _j \|^2= R_0^2d _k=1^d (s_i,k-s_j,k )^2= 4R_0^2dH Thus, the event ‖μi−μj‖2<2(1−γ)R02 \| _i- _j \|^2<2(1-γ)R_0^2 is equivalent to 4R02dH<2(1−γ)R02⟺H<1−γ2d=d2−γd2. 4R_0^2dH<2(1-γ)R_0^2 H< 1-γ2d= d2- γ d2. Applying Hoeffding’s inequality to H, we have that, for any t>0t>0, ℙ(H−H≤−t)≤exp(−2t2d) (H-EH≤-t)≤ (- 2t^2d ) ⇔ ℙ(H−d2≤−t)≤exp(−2t2d). (H- d2≤-t)≤ (- 2t^2d ). Choosing t=γd2t= γ d2, we have that ℙ(H−H≤−γd2)≤exp(−γ22d) (H-EH≤- γ d2)≤ (- γ^22d ) ⇔ ℙ(H≤d2−γd2)≤exp(−γ22d) (H≤ d2- γ d2)≤ (- γ^22d ) ⇔ ℙ(H≤1−γ2d)≤exp(−γ22d) (H≤ 1-γ2d)≤ (- γ^22d ) ⇔ ℙ(‖μi−μj‖2<2(1−γ)R02)≤exp(−γ22d). ( \| _i- _j \|^2<2(1-γ)R_0^2 )≤ (- γ^22d ). ∎ Lemma 16 (uniform separation for all pairs of means of stored patterns). Let μ1,…,μN _1,…, _N be defined as in Lemma 15 using the i.i.d. componentwise Rademacher random variables s1,…,sN∈±1ds_1,…,s_N∈\± 1\^d. Fix p∈(0,1)p∈(0,1) and γ∈(0,1)γ∈(0,1). If N≔⌊2pexp(γ24d)⌋.N 2p ( γ^24d ) . then with probability at least 1−p1-p, ‖μi−μj‖≥dmin≔2(1−γ)R0 for all i≠j. \| _i- _j \|≥ d_ 2(1-γ)R_0 for all i≠ j. Proof. Let AijA_ij be the bad event ‖μi−μj‖2<2(1−γ)R02 \ \| _i- _j \|^2<2(1-γ)R_0^2 \. By Lemma 15, ℙ(Aij)≤exp(−γ22d) for each i≠j.P (A_ij )≤ (- γ^22d ) for each i≠ j. By the union bound over (N2)≤N22 N2≤ N^22 pairs, ℙ(∃i<j:Aij)≤N22exp(−γ22d).P (∃\,\ i<j:A_ij )≤ N^22 (- γ^22d ). Under the stated condition on N, the RHS is ≤p≤ p. Therefore, with probability at least 1−p1-p, no bad event occurs, i.e. all pairs satisfy ‖μi−μj‖2≥2(1−γ)R02⟺‖μi−μj‖≥2(1−γ)R0=dmin. \| _i- _j \|^2≥ 2(1-γ)R_0^2 \| _i- _j \|≥ 2(1-γ)R_0=d_ . ∎ Lemma 17 (Linear independence of Dirac distribution and its distributional derivative). Let (x1,…,xM)∈Ω⊂ℝd(x_1,…,x_M)∈ ^d be pairwise distinct location parameters. Suppose ∑i=1Mciδxi+∑i=1Mvi⋅∇δxi=0 in (C∞(ℝd))′, _i=1^Mc_i _x_i+ _i=1^Mv_i·∇ _x_i=0 in (C^∞(R^d) ) , with ci∈ℝc_i and vi∈ℝdv_i ^d. Then ci=0c_i=0 and vi=0v_i=0 for all i=1,…,Mi=1,…,M. Proof. Fix k∈1,…,Mk∈\1,…,M\. Given any x∈Ωx∈ and r≥0r≥ 0, define the Euclidean ball B(x,r)≔y∈Ω:‖y−x‖≤rB(x,r) \y∈ :\|y-x\|≤ r \. Since the xix_i’s are distinct, there exists rk>0r_k>0 such that B(xk,rk)∩xj:j≠k=∅B (x_k,r_k )∩ \x_j:j≠ k \= . Let us choose ψ∈Cc∞(ℝd)⊂C∞(ℝd)ψ∈ C_c^∞ (R^d )⊂ C^∞(R^d) supported in B(xk,rk)B (x_k,r_k ) and equal to 1 in a neighborhood of xkx_k. Then, we have that, ⟨δxi,ψ⟩=ψ(xi)=1,i=kψ(xk)=0,i≠k _x_i,ψ = casesψ(x_i)=1,&i=k\\ ψ(x_k)=0,&i≠ k cases and ⟨∇δxi,ψ⟩=−∇ψ(xi)=0 ∇ _x_i,ψ =-∇ψ (x_i )=0 for i=1,…,Mi=1,…,M, since ψ is constant near xkx_k and vanishes near xix_i for i≠ki≠ k. Consequently, we have that ∑i=1M⟨ciδxi,ψ⟩+∑i=1M⟨vi⋅∇δxi,ψ⟩=⟨0,ψ⟩ _i=1^M c_i _x_i,ψ + _i=1^M v_i·∇ _x_i,ψ = 0,ψ ⇔ ∑i=1Mciψ(xi)−∑i=1Mvi⋅∇ψ(xi)=0 _i=1^Mc_iψ(x_i)- _i=1^Mv_i·∇ψ(x_i)=0 ⇔ ck=0. c_k=0. Now, let us fix any vector u∈x⊂ℝdu∈ x ^d. Then, let us choose ψu∈Cc∞(ℝd)⊂C∞(ℝd) _u∈ C_c^∞ (R^d )⊂ C^∞(R^d) supported in B(xk,rk)B(x_k,r_k) with ψu(xk)=0 _u(x_k)=0 and ∇ψu(xk)=u∇ _u(x_k)=u. Then, we have that ⟨δxi,ψu⟩=0 for i=1,…,M _x_i, _u =0 for i=1,…,M and ⟨∇δxi,ψu⟩=−∇ψu(xi)=−∇ψu(xk)=−u,i=k,−∇ψu(xi)=0,i≠k ∇ _x_i, _u =-∇ _u(x_i)= cases-∇ _u(x_k)=-u,&i=k,\\ -∇ _u(x_i)=0,&i≠ k cases. Consequently, we have that ∑i=1M⟨ciδxi,ψu⟩+∑i=1M⟨vi⋅∇δxi,ψu⟩=⟨0,ψu⟩ _i=1^M c_i _x_i, _u + _i=1^M v_i·∇ _x_i, _u = 0, _u ⇔ ∑i=1Mciψu(xi)−∑i=1Mvi⋅∇ψu(xi)=0 _i=1^Mc_i _u(x_i)- _i=1^Mv_i·∇ _u(x_i)=0 ⇔ −vk⋅u=0. -v_k· u=0. Since this is true for any u∈Ω∈ℝdu∈ ^d, we must have that vk=0v_k=0. Finally, since k was arbitrary, we must have that ci=0c_i=0 and vi=0v_i=0 for all i=1,…,Mi=1,…,M. ∎ Appendix J Numerical Experiments In this section, we demonstrate the empirical performance of our proposed retrieval algorithm, referred to as SinkhornSHK Algo, for finitely supported discrete measures and compare it to a baseline Euclidean geometry based classical Hopfield-type algorithm, which we refer to as Euclidean Algo. The Euclidean algo vectorizes (a,x)∈ℳM(a,x) _M into ξvec=[x11,…,x1d,…,xM1,…,xMd,loga1,…,logam]∈ℝ(d+1)M _vec= [x_11,…,x_1d,…,x_M1,…,x_Md, a_1,…, a_m ] ^(d+1)M and applied the classical Hopfield fixed point algorithm for vector inputs based on Euclidean ℓ2 _2 inner product similarity using Equation 3 of Ramsauer et al. (2020) with the same choice of β as for our proposed SinkhornSHK Algo. We consider a toy experiment where the stored patterns X1,…,XNX_1,…,X_N are uniformly weighted and the support points are sampled from Gaussian distributions. We choose N=5N=5 and data dimension d=2d=2. In Experiment 1, we choose the means of the Gaussian distributions to be (−4.0,−1.0)(-4.0,-1.0),(−2.0,2.2)(-2.0,2.2), (1.0,−6.0)(1.0,-6.0),(4.0,−4.2)(4.0,-4.2) and (4.2,−0.8)(4.2,-0.8), while the covariance matrices were chosen to be [0.600.200.200.90] bmatrix0.60&0.20\\ 0.20&0.90 bmatrix, [0.80−0.15−0.150.55],[0.650.000.000.65],[0.550.100.101.00] bmatrix0.80&-0.15\\ -0.15&0.55 bmatrix, bmatrix0.65&0.00\\ 0.00&0.65 bmatrix, bmatrix0.55&0.10\\ 0.10&1.00 bmatrix and [0.950.000.000.50] bmatrix0.95&0.00\\ 0.00&0.50 bmatrix. M=30M=30 support points were sampled in i.i.d manner from the 5 Gaussian distributions determined by each pair of mean and covariance parameters and the resulting uniformly weighted discrete distributions were set as the patterns to be stored. In Experiment 2, we choose the means of the 5 Gaussian distributions to be all equal to (0,0)(0,0) and the covariance matrices were randomly sampled using random orthogonal matrices coupled with uniformly sampled eigenvalues between 0.15 and 1.75. We sample M=25M=25 support points in i.i.d manner from each of these Gaussian distributions and the resulting uniformly weighted discrete distributions were set as the patterns to be stored. In both the experiments, we first fix a pattern that we want to retrieve, then perturb the support points individually using i.i.d Gaussian noise (sd 0.5 in Experiment 1 and sd 0.2 in Experiment 2) to generate a query distribution that serves as the initial iterate ξ(0)ξ^(0) for both algorithms. We chose β=50β=50, ε=0.05 =0.05 and step-size η=1.3η=1.3 for the SinkhornSHK Algo, and the same β for Euclidean Algo. We do not use the spherical Hellinger update step in these simple experiments since the all discrete measures involved are uniformly weighted. We use a maximum iteration threshold k≤200k≤ 200 for both algorithms, and the Sinkhorn algorithm for computing entropic OT transport plans were capped at 120 iterations. In Experiment 1, we see that both algorithms are able to retrieve the correct discrete distributions when considering the support points that are returned by either algorithm. However, Experiment 2 clearly shows the superiority of SinkhornSHK Algo over Euclidean Algo, since the Sinkhorn Algo is able to converge to the correct pattern even when a noisy query is given. We believe that the ability of Sinkhorn ALgo to leverage the distributional perspective gives it the advantage over Euclidean Algo, since the latter relies on Euclidean inner products and is expected to fail in cases where Euclidean separation between support points is small, but separation in distributional metrics is still feasible. Figure 1: Experiment 1: Sinkhorn Algo and Euclidean Algo are both able to retrieve correct patterns from noisy queries Figure 2: Experiment 2: Sinkhorn Algo succeeds in retrieving correct patterns from noisy queries in all instances, but Euclidean Algo fails in 3 cases.