MD J Taytard

MD S Valence

MD C Sileo

MD D Rodriguez

MD G Aubertin

MD H Corvol

MD N Beydon

Severe central apnea secondary to cerebellar dysplasia in a child: look past Joubert syndrome

We report the case of a 12 years old girl referred to our pediatric sleep unit with a history of central sleep apnea (CSA), associated to transient episodes of tachypnea on polysomnography recordings.

The patient was otherwise healthy, with no personal or family medical history, and had normal physical and neuropsychological examination. Brain magnetic resonance imaging (MRI) showed signs of cerebellar vermis dysplasia, but without the classical features of the molar tooth sign. The rest of the workup (genetic tests, blood tests, cardiac investigations) was normal, except for an increased peripheral chemosensitivity to carbon dioxide (CO2) and oxygen (O2). The patient was successfully treated with bi-level positive airway pressure (bi-level PAP). This case highlights the importance of performing a brain MRI in case of CSA in order to study the cerebellum, beyond the brainstem area. Cerebellum malformations can be found even in the absence of any other neurological condition.

INTRODUCTION

Central sleep apnea (CSA) in children can be attributed to congenital malformations (Arnold-Chiari malformation or craniofacial abnormalities) or to genetic (congenital central hypoventilation syndrome, Down syndrome, etc.), neurological or endocrine conditions [START_REF] Mclaren | Diagnosis, management and pathophysiology of central sleep apnea in children[END_REF] . Although the ventilatory drive to breathe is mainly under brainstem structures control, cerebellar malformations are sometimes responsible for ventilatory drive impairment, resulting in CSA. In this case, the most frequent diagnosis is Joubert syndrome and related disorders (JSRD) with a sleep pattern of central apnea associated with episodes of tachypnea [START_REF] Rama | Clinical and Molecular Diagnosis of Joubert Syndrome and Related Disorders[END_REF] . Joubert syndrome is characterized by a vermian dysplasia (with typical cleft in the superior vermis) associated with thickening and horizontalization of the superior cerebellar peduncles leading to the so-called "molar tooth" pattern seen on Magnetic Resonance Imaging (MRI) 2 .

In this case report, we present a 12 years old girl with severe CSA and transient episodes of tachypnea associated with cerebellar dysplasia affecting the vermis. In this patient, neither the typical image of "molar tooth" sign nor any of the usual clinical features associated to JRDS were present. She was successfully treated with bi-level positive airway pressure (bi-level PAP) ventilation.

REPORT OF CASE

A 12 years old girl was referred to our pediatric sleep unit with a history of CSA. CSA had been diagnosed in another center at 7 years old, after the patient's parents witnessed apneas during sleep.

There was no significant personal or family medical history, and the patient was not under any medication. She met all developmental milestones and followed a normal scholar education.

At the time of diagnosis, there was no daytime clinical sign, especially no symptoms of sleepiness (Epworth score of 2/24), but apneas witnessed by the parents during sleep and nocturia. Initial polysomnography (PSG) showed a central apnea index (CAI) of 16/h and an obstructive sleep apnea hypopnea index of 6/h, with normal nocturnal gas exchanges. The patient was started on Continuous Positive Airway Pressure (CPAP) without any improvement, and then switched to bilevel PAP with good tolerance. Concurrent brain MRI was considered normal.

Upon her arrival in our department, physical examination was normal (height: 158 cm; weight: 56 kg; body mass index (BMI): 22.4 kg.m -2 , 1.42 Z-score). Awake respiratory rate at rest was of 16 breaths per minute and Oxygen (O2) saturation in room air was 98%. Lung auscultation and neurological examination were both normal.

A general work-up was initiated in our center including the search for PHOX2B mutation which was negative. Cardiac investigations (echocardiography, electrocardiogram with 24-hours holter monitoring, tension holter monitoring) were normal, except for a sinus arrhythmia (with a total of 21 supraventricular extra-systoles/24h) which did not require any specific treatment. Arterialized capillary blood gas sampled at earlobe showed a PaO2 of 93 mmHg and a PaCO2 of 37 mmHg.

Central carbon dioxide (CO2) chemosensitivity was normal but there was an increased peripheral chemosensitivity to CO2 and to O2. The pulmonary function test showed moderate, significantly reversible bronchial obstruction. Blood tests, including thyroid function test, were normal.

Considering the absence of diagnosis, a new brain MRI was performed and showed cerebellar vermis dysplasia and hypoplasia without the characteristic appearance of midbrain "molar tooth" sign on axial images. Neither horizontalization nor thickening of the superior cerebellar peduncles were present (Figure 1). Following the diagnosis of cerebellar malformation, a neuropsychological assessment was found normal and the normal neurological examination was confirmed, without cerebellar or oculomotor apraxia.

New PSG were realized under different conditions (Table 1). All recordings were performed and scored according to the American Association for Sleep Medicine (AASM) recommendations and scoring rules [START_REF] Berry | Rules for scoring respiratory events in sleep: update of the 2007 AASM Manual for the Scoring of Sleep and Associated Events. Deliberations of the Sleep Apnea Definitions Task Force of the American Academy of Sleep Medicine[END_REF] . Acknowledging for the increased peripheral chemosensitivity, PSG were performed with or without bi-level PAP, in room air or with nasal O2 therapy (Figure 2). The latter did not elicit a decrease in central apneas, and showed an increased periodic breathing cycle length compared with the recording in room air (70 vs 50 seconds, respectively). Whatever the ventilatory conditions, nocturnal gas exchanges remained within the normal range.

Bi-level PAP was pursued under S/T mode with dramatic efficacy on CSA (Table 1). Since the clinical tolerance was perfect, the machine parameters were not modified, and set as follows: inspiratory positive airway pressure of 11 cmH2O, expiratory positive airway pressure of 4 cmH2O, and breathing frequency of 14/min.

DISCUSSION

We describe a unique case of CSA associated with cerebellar vermis dysplasia in an otherwise healthy 12 years old girl showing no specific neurological condition.

The AASM defines CSA as the cessation of airflow associated to the absence of chest wall and/or abdominal movements for longer than 20 seconds or lasting more than 2 baseline respiratory cycles if associated to an arousal, an awakening or an O2 desaturation of at least 3% [START_REF] Berry | Rules for scoring respiratory events in sleep: update of the 2007 AASM Manual for the Scoring of Sleep and Associated Events. Deliberations of the Sleep Apnea Definitions Task Force of the American Academy of Sleep Medicine[END_REF] . CAI is considered to be pathological when equal or higher than 5/h 1 . In children, CSA has an estimated prevalence of 1 to 5% 1 .

A brain MRI is systematically recommended to explore CSA, with a special interest for the brainstem area where the respiratory centers are located. However, attention should also be paid to the cerebellum which has been implicated in central chemoreception [START_REF] Calton | Impaired hypercarbic and hypoxic responses from developmental loss of cerebellar Purkinje neurons: implications for sudden infant death syndrome[END_REF] . Unlike the brainstem area, the cerebellum does not contain central pattern generators. It encompasses deep cerebellar nuclei that are partly responsible for the response to stressed breathing and hypoxic or hypercapnic challenges [START_REF] Calton | Impaired hypercarbic and hypoxic responses from developmental loss of cerebellar Purkinje neurons: implications for sudden infant death syndrome[END_REF][START_REF] Harper | The cerebellum and respiratory control[END_REF] . In particular, the fastigial nucleus has been reported to have a potential excitatory effect on breathing via the facilitation of respiratory responses mediated by carotid chemoreceptors [START_REF] Harper | The cerebellum and respiratory control[END_REF] .

Disordered breathing during sleep is an important feature of syndromes including cerebellar malformations, such as the Chiari type II malformation [START_REF] Harper | The cerebellum and respiratory control[END_REF] . Malformations are frequently secondary to genetic diseases and often associated to other neurological dysfunctions of varying severity such as developmental delay, cerebellar ataxia, severe hypotonia, dysarthria or oculomotor apraxia 7,[START_REF] Delrosso | The cerebellum and sleep[END_REF] . JRDS are the most frequent genetic cerebellar diseases involving abnormal respiratory pattern and CSA. A pattern of episodes of tachypnea followed by central apnea is considered as a hallmark of JRDS 8 .

Our patient did not have any clinical and radiological sign of JRDS, therefore, other cerebellar malformations, especially vermian malformations, should be systematically searched for to explain cases of severe CSA syndrome. Although considered typical of JRDS, episodes of transient tachypnea have been described in other rare cases of cerebellar vermis hypoplasia, but usually in association with hypoplastic brainstem in patients presenting other neurological symptoms.

In a recent literature review on CSA in children, McLaren and colleagues pointed out that individuals with high loop gain were prone to ventilatory instability and at risk of periodic breathing [START_REF] Mclaren | Diagnosis, management and pathophysiology of central sleep apnea in children[END_REF] . In our case with increased peripheral chemosensitivity, the O2 therapy, which levels blood O2 content, should have drastically reduced central respiratory events [START_REF] Wellman | Effect of oxygen in obstructive sleep apnea: role of loop gain[END_REF] . The absence of CSA improvement with O2 therapy, alongside with the higher mean periodic breathing cycle length compared to the recording in room air, makes the high loop gain unlikely to be the sole responsible for CSA in our patient. However, as O2 therapy is a possible therapeutic option for patients with periodic breathing, testing central and peripheral chemosensitivity should be considered when investigating CSA. The use of Acetazolamide is another therapeutic option in case of CSA associated with hyperventilation. This treatment has been used in patients with Cheyne-stokes secondary to heart failure or in patients experiencing CSA at high altitude. Acetazolamide is a carboanhydrase inhibitor which interferes with kidney reuptake of bicarbonate, resulting in a metabolic acidosis. It is thought to stimulate chemosensitivity by inducing changes in tissue and systemic pH, and by decreasing peripheral chemoreceptors sensitivity to PaO2, therefore stabilizing breathing control [START_REF] Verhulst | Acetazolamide for severe apnea in Pitt-Hopkins syndrome[END_REF] . However, data on the use of Acetazolamide in children are very scarce, and limited to case reports. It could be considered as a second-line treatment in our patient, if poor tolerance to bi-level PAP therapy appears.

ABBREVIATIONS

  ABBREVIATIONS

Figure 1 .

 1 Figure 1. Brain magnetic resonance imaging. Midline sagittal T2-weighted image showing vermian dysplasia.

Figure 2 .

 2 Figure 2. PSG recordings (A) Recording in room air and without bi-level PAP (B) Recording under oxygen therapy (1L/min) (C) Recording with bi-level PAP

Table 1 . PSG results

 1 Bi-level PAP: bi-level positive airway pressure; TST: total sleep time; AHI: apnea hypopnea index; OAHI: obstructive apnea hypopnea index; CAHI: central apnea hypopnea index; ODI: oxygen desaturation index; SpO2: pulse oximetry; TcPCO2: transcutaneous partial pressure of carbon dioxide

	Sleep architecture	Room air	Oxygen therapy	Bi-level PAP
			(1L/min)
	TST (min)	535	538	483
	Sleep efficiency (%)	96	97	96
	REM sleep (% TST)	21	18	21
	N1 and N2 stages (% TST)	52	50	46
	Micro arousal (h -1 )	17	16	9
	Respiratory characteristics and
	events
	Respiratory rate (min -1 )	16	16	18
	Respiratory cycle length (s)	3.7	3.7	3.3
	Snore index (h -1 )	215	192	130
	AHI (h -1 )	18	17	1.5
	OAHI (h -1 )	8	5	1.5
	CAHI (h -1 )	10	12	0
	ODI (h -1 )	12	9	4.5
	Mean periodic breathing cycle	50	70	-
	length (s)
	Mean apnea length (s)	22	30	0
	Maximum apnea length (s)	46	72	0
	Gas exchanges
	Mean SpO2 (%)	97	98	96
	SpO2 nadir (%)	87	88	91
	Mean TcPCO2 (mmHg)	40	43	45
	Max TcPCO2 (mmHg)	44	47	49
	TcPCO2 > 50 mmHg (% of TST)	0	0	0
	REM: rapid eye movement;


Mokhtar Z Alaya
email: elmokhtar.alaya@utc.fr

Gilles Gasso
email: gilles.gasso@insa-rouen.fr

Maxime Bérar
email: maxime.berar@univ-rouen.fr

Alain Rakotomamonjy

HETEROGENEOUS WASSERSTEIN DISCREPANCY FOR INCOMPARABLE DISTRIBUTIONS

Optimal Transport (OT) metrics allow for defining discrepancies between two probability measures. Wasserstein distance is for longer the celebrated OT-distance frequently-used in the literature, which seeks probability distributions to be supported on the same metric space. Because of its high computational complexity, several approximate Wasserstein distances have been proposed based on entropy regularization or on slicing, and one-dimensional Wassserstein computation. In this paper, we propose a novel extension of Wasserstein distance to compare two incomparable distributions, that hinges on the idea of distributional slicing, embeddings, and on computing the closed-form Wassertein distance between the sliced distributions. We provide a theoretical analysis of this new divergence, called heterogeneous Wasserstein discrepancy (HWD), and we show that it preserves several interesting properties including rotation-invariance. We show that the embeddings involved in HWD can be efficiently learned. Finally, we provide a large set of experiments illustrating the behavior of HWD as a divergence in the context of generative modeling and in query framework.

INTRODUCTION

Optimal Transport-based data analysis has recently found widespread interest in machine learning community, since its significant usefulness to achieve many tasks arising from designing loss functions in supervised learning (Frogner et al., 2015), unsupervised learning (Arjovsky et al., 2017), text classification (Kusner et al., 2015), domain adaptation (Courty et al., 2017), generative models (Arjovsky et al., 2017;Salimans et al., 2018), computer vision (Bonneel et al., 2011;Solomon et al., 2015) among many more applications (Kolouri et al., 2017;Peyré & Cuturi, 2019). Optimal Transport (OT) attempts to match real-world entities through computing distances between distributions, and for that it exploits prior geometric knowledge on the base spaces in which the distributions are valued. Computing OT distance equals to finding the most cost-efficiency way to transport mass from source distribution to target distribution, and it is often referred to as the Monge-Kantorovich or Wasserstein distance (Monge, 1781;Kantorovich, 1942;Villani, 2009).

Matching distributions using Wasserstein distance relies on the assumption that their base spaces must be the same, or that at least a meaningful pairwise distance between the supports of these distributions can be computed. A variant of Wasserstein distance dealing with heterogeneous distributions and overcoming the lack of intrinsic correspondence between their base spaces is Gromov-Wasserstein (GW) distance (Sturm, 2006;Mémoli, 2011). GW distance allows to learn an optimal transport-like plan by measuring how the similarity distances between pairs of supports within each ground space are closed. It is increasingly finding applications for learning problems in shape matching (Mémoli, 2011), graph partitioning and matching (Xu et al., 2019), matching of vocabulary sets between left RIGHT center Figure 1: We measure the discrepancy between two distributions living respectively in R p and R q . Our approach is based on generating random slicing projections distributions in each of the metric spaces R p and R q through the mappings φ and ψ of a random projection vector sampled from an optimal distribution Θ in R d . As each of the projected distribution results in a 1D distribution, we can then compute 1D-Wasserstein distance. It enables us to learn the best projection mappings φ and ψ and to optimize over the distributional part of the generating projection distribution Θ.

different languages (Alvarez-Melis & Jaakkola, 2018), generative models (Bunne et al., 2019), or matching weighted networks (Chowdhury & Mémoli, 2018).

Due to the heterogeneity of the distributions, GW distance uses only the relational aspects in each domain, such as the pairwise relationships to compare the two distributions. As a consequence, the main disadvantage of GW distance is its computational cost as the associated optimization problem is a non-convex quadratic program (Peyré & Cuturi, 2019), and as few as thousand samples can be computationally challenging.

Based on the approach of regularized OT (Cuturi, 2013), in which an entropic penalty is added to the original objective function defining the Wasserstein OT problem, Peyré et al. (2016) propose an entropic version called entropic GW discrepancy, that leads to approximate GW distance. Another approach for scaling up the GW distance is Sliced Gromov-Wasserstein (SGW) discrepancy (Vayer et al., 2019), which leverages on random projections on 1D and on a closed-form solution of the 1D-Gromov-Wasserstein.

In this paper, we take a different approach for measuring the discrepancy between two heterogeneous distributions. Unlike GW distance that compares pairwise distances of elements from each distribution, we consider a method that embeds the metric measure spaces into a one-dimensional space and computes a Wasserstein distance between the two 1D-projected distributions. The key element of our approach is to learn two mappings that transform vectors from the unit-sphere of a latent space to the unit-sphere of the metric space underlying the two distributions of interest, see Figure 1. In a nutshell, we learn to transform a random direction, sampled under an optimal (learned) distribution (optimality being made clear later), from a d-dimensional space to a random direction into the desired spaces. This approach has the benefit of avoiding an ad-hoc padding strategy (completion of 0 of the smaller dimension distributions to fit the high-dimensional one) as in SGW method (Vayer et al., 2019). Another relevant feature of our approach is that the two resulting 1D distributions are now compared through Wasserstein distance. This point, in conjunction, with other key aspect of the method, will lead to a relevant discrepancy between two distributions, called heterogeneous Wasserstein discrepancy (HWD). Although we lose some properties of a distance, we show that HWD is rotation-invariant, that it is robust enough to be considered as a loss for learning generative models between heterogeneous spaces. We also establish that HWD boils down to the recent distributional sliced Wasserstein distance (Nguyen et al., 2020) if the two distributions live in the same space and if some mild constraints are imposed on the mappings.

In summary, our contributions are as follows:

• we propose HWD, a novel slicing-based discrepancy for comparing two distributions living in different spaces. Our chosen formulation is based on comparing 1D random-projected versions of the two distributions using a Wasserstein distance;

• The projection operations are materialized by optimally mapping from one common space to the two spaces of interest. We provide a theoretical analysis of the resulting discrepancy and exhibit its relevant properties;

• Since the discrepancy involves several mappings that need to be optimized, we depict an alternate optimization algorithm for learning them;

• Numerically, we validate the benefits of HWD in terms of comparison between heterogeneous distributions. We show that it can be used as a loss for generative models or shape objects retrieval with better performance and robustness than SGW on those tasks.

BACKGROUND OF OT DISTANCES

For the reader's convenience, we provide here a brief review of the notations and definitions, that will be frequently used throughout the paper. We start by introducing Wasserstein and Gromov-Wasserstein distances with their sliced versions SW and SGW, where we consider these distances in the specific case of Euclidean base spaces (R p , • ) and (R q , • ).

We denote P(X ) and P(Y) the respective sets of probability measures whose supports are contained on compact sets X ⊆ R p and Y ⊆ R q . For r ≥ 1, we denote P r (X ) the subset of measures in P(X ) with finite r-th moment (r ≥ 1), i.e., P r (X ) = η ∈ P(X ) : ´X x r dη(x) < ∞ . For µ ∈ P(X ) and ν ∈ P(Y), we write Π(µ, ν) ⊂ P(X × Y) for the collection of joint probability distributions with marginals µ and ν, known as couplings,

Π(µ, ν) = γ ∈ P(X × Y) : ∀A ⊂ X , B ⊂ Y, γ(A × Y) = µ(A), γ(X × B) = ν(B) .

OT DISTANCES FOR HOMOGENEOUS DOMAINS

We here assume that the distributions µ and ν lie in the same base space, for instance p = q. Taking this into account, we can define the Wasserstein distance and its sliced variant.

Wasserstein distance

The r-th Wasserstein distance is defined on P r (X ) by

W r (µ, ν) = inf γ∈Π(µ,ν) ˆX ×Y x -y r dγ(x, y) 1 r .
(1)

The quantity W r (µ, ν) describes the least amount effort to transform one distribution µ into another one ν. Since the cost distance used between sample supports is the Euclidean one, the infimum in (1) is attained (Villani, 2009), and any probability γ which realizes the minimum is called an optimal transport plan. In a finite discrete setting, Problem (1) can be formulated as a linear program, that is challenging to solve algorithmically as its computational cost is of order O(n 5/2 log n) (Lee & Sidford, 2014), where n is the number of sample supports.

Contrastingly, for the 1D case (i.e. p = 1) of continuous probability measures, the r-th Wasserstein distance has a closed-form solution (Rachev & Rüschendorf, 1998)

, namely, W r (µ, ν) = ( ´1 0 |F -1 µ (u) -F -1 ν (u)| r dt) 1 r
where F -1 µ and F -1 ν are the quantile functions of µ and ν. For empirical distributions, the 1D-Wasserstein distance is simply calculated by sorting the supports of the distributions on the real line, resulting to a complexity of order O(n log n). This nice computational property motivates the use of sliced-Wasserstein (SW) distance (Rabin et al., 2012;Bonneel et al., 2015), where one calculates an (infinity) of 1D-Wasserstein distances between linear projection pushforwards of distributions in question and then computes their average.

To precisely define SW distance, we consider the following notation. Let S p-1 := {u ∈ R p : u = 1} be the unit sphere in p dimension in 2 -norm, and for any vector θ in S p-1 , we define P θ the orthogonal projection onto the real line Rθ = {αθ : α ∈ R}, that is P θ (x) = θ, x , where

•, • stands for the Euclidean inner-product. Let µ θ = P θ #µ the measure on the real line called pushforward of µ by P θ , that is µ θ (A) = µ(P -1 θ (A)) for all Borel set A ⊆ R. We may now define the SW distance.

Sliced Wasserstein distance The r-th order sliced Wasserstein distance between two probability distributions µ, ν ∈ P r (X ) is given by

SW r (µ, ν) = 1 A p ˆSp-1 W r r (µ θ , ν θ )dθ 1 r , (2)
where A p is the area of the surface of S p-1 , i.e., A p = 2π p/2 Γ(p/2) with Γ : R → R, the Gamma function given as Γ(u) = ´∞ 0 t u-1 e -t dt.

Thanks to its computational benefits and its valid metric property (Bonnotte, 2013), the SW distance has recently been used for OT-based deep generative modeling (Kolouri et al., 2019;Deshpande et al., 2019;Wu et al., 2019). Note that the normalized integral in (2) can be seen as the expectation for θ ∼ σ p-1 , the uniform surface measure on

S p-1 , that is SW r (µ, ν) = (E θ∼σ p-1 [W r r (µ θ , ν θ )]) 1 r .
Therefore, the SW distance can be easily approximated via a Monte Carlo sampling scheme by drawing uniform random samples from S p-1 : SW r r (µ, ν)

≈ 1 K K k=1 W r r (µ θ k , ν θ k ) where θ 1 , . . . , θ K i.i.d.
∼ σ p-1 and K is the number of random projections.

OT DISTANCES FOR HETEROGENEOUS DOMAINS

To get benefit from the advantages of OT in many machine learning applications involving heterogeneous and incomparable domains (p = q), the Gromov-Wasserstein distance (Mémoli, 2011) stands for the basic OT distance dealing with this setting.

Gromov-Wasserstein distance The r-th Gromov-Wasserstein distance between two probability distributions µ ∈ P r (X ) and ν ∈ P r (Y) is defined by

GW r (µ, ν) = inf γ∈Π(µ,ν) J r (γ) def. = 1 2 Ẍ 2 ×Y 2 | x -x -y -y | r dγ(x, y)dγ(x , y ) 1 r . (3)
Note that GW r (µ, ν) is a valid metric endowing the collection of all isomorphism classes metric measure spaces of P r (X ) × P r (Y), see Theorem 5 in (Mémoli, 2011). The GW distance learns an optimal transport-like plan which transports samples from a source metric space X into a target metric space Y, by measuring how the similarity distances between pairs of samples within each space are close. Furthermore, GW distance enjoys several geometric properties, particularly translation and rotation invariance. However, its major bottleneck consists in an expensive computational cost, since problem (3) is non-convex and quadratic. A remedy to such a heavy computational burden lies in an entropic regularized GW discrepancy (Peyré et al., 2016), using Sinkhorn iterations algorithm (Cuturi, 2013). This latter needs a large regularization parameter to guarantee a fast computation, which, unfortunately, entails a poor approximation of the true GW distance value.

Another approach to scale up the computation of GW distance is sliced-GW discrepancy (Vayer et al., 2019). The definition of SGW shows 1D-GW distances between projected pushforward of an artifact zero padding of µ or ν distribution. We detail this representation in the following paragraph.

Sliced Gromov-Wasserstein discrepancy Assume that p < q and let ∆ be an artifact zero padding from X onto Y, i.e. ∆(x) = (x 1 , . . . , x p , 0, . . . , 0) ∈ R q . The r-th order sliced Gromov-Wasserstein discrepancy between two probability distributions µ ∈ P r (X ) and ν ∈ P r (Y) is given by

SGW ∆,r (µ, ν) = E θ∼σ q-1 GW r r ((∆#µ) θ , ν θ ) 1 r . ( 4
)
It is worthy to note that SGW ∆,r is depending on the ad-hoc operator ∆, hence the rotation invariance is lost. Vayer et al. (2019) propose a variant of SGW that does not depend on the choice of ∆, called Rotation Invariant SGW (RI-SGW) for p = q, defined as the minimizer of SGW ∆,r over the Stiefel manifold, see (Vayer et al., 2019, Equation 6). In this work, we are interested in calculating an OT-based discrepancy between distributions over distinct domains using the slicing technique. Our approach is different from the SGW one in many points, specifically (and most importantly) we use a 1D-Wasserstein distance between the projected pushforward distributions and not a 1D-GW distance.

In the next section, we detail the setup of our approach.

HETEROGENEOUS WASSERSTEIN DISCREPANCY

Despite the computational benefit of sliced-OT variant discrepancies, they have an unavoidable bottleneck corresponding to an intractable computation of the expectation with respect to uniform distribution of projections. Furthermore, the Monte Carlo sampling scheme can often generate an overwhelming number of irrelevant directions; hence, the larger number of sample projections, the more accurate approximation of sliced-OT values. Recently, Nguyen et al. (2020) have proposed the distributional-SW distance allowing to find an optimal distribution over an expansion area of informative directions. This performs the projection efficiently by choosing an optimal number of important random projections needed to capture the structure of distributions. Our approach for comparing distributions in heterogeneous domains follows a distributional slicing technique combined with OT metric measure embedding (Alaya et al., 2020).

Let us first introduce additional notations. Fix d ≥ 1 and consider two nonlinear mappings φ : S d-1 → S p-1 and ψ : S d-1 → S q-1 . For any constants C φ , C ψ > 0, we define the following probability measure sets:

M C φ = Θ ∈ P(S d-1 ) : E θ,θ ∼Θ [| φ(θ), φ(θ ) |] ≤ C φ and M C ψ = Θ ∈ P(S d-1 ) : E θ,θ ∼Θ [| ψ(θ), ψ(θ ) |] ≤ C ψ . We say that C φ , C ψ are (φ, ψ)- admissible constants if the intersection sets M C φ ∩ M C ψ is not empty.
We hereafter denote µ φ,θ = P φ(θ) #µ and ν ψ,θ = P ψ(θ) #ν the pushforwards of µ and ν by projections over unit sphere P φ(θ) and P ψ(θ) , respectively.

Informal presentation While the distributions µ and ν are valued in different spaces, X ⊂ R p and Y ⊂ R q , any projected distributions will live in real line, enabling the computation of 1D-Wasserstein distance (Figure 1, right). In order to generate random 1D projections in each of the spaces, we map a common random projection distribution from S d-1 into each of the projection spaces S p-1 and S q-1 , through the mappings φ and ψ (see Figure 1, left). Hence, the main components of the heterogeneous Wasserstein discrepancy will be the distribution Θ ∈ P(S d-1 ),

and the two embeddings φ and ψ which will be wisely chosen. The resulting directions φ(θ) and ψ(θ) form the projections P φ(θ) and P ψ(θ) (see Figure 1, center) used to compute several 1D-Wasserstein distances.

DEFINITION AND PROPERTIES

Herein we state the formulation of the proposed discrepancy and exhibit its main theoretical properties.

Definition 1 The heterogeneous Wasserstein discrepancy (HWD) of order r ≥ 1 between µ ∈ P r (X ) and ν ∈ P r (Y) reads as

HWD r (µ, ν) = inf φ,ψ sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r (µ φ,θ , ν ψ,θ ) 1 r . (5)
HWD belongs to a family of projected OT works (Paty & Cuturi, 2019;Rowland et al., 2019;Lin et al., 2021) with a particularity for seeking nonlinear projections minimizing a sliced-OT variant.

HWD further inherits the distributional slicing benefit by finding an optimal probability measure Θ of slices on the unit sphere S d-1 coupled with an optimum couple (φ, ψ) of embeddings. Note that this optimal Θ verifies the double conditions

E θ,θ ∼Θ [| cos(φ(θ), φ(θ ))|] ≤ C φ and E θ,θ ∼Θ [| cos(ψ(θ), ψ(θ ))|] ≤ C ψ . This gives that C φ , C ψ ≤ 1, hence the sets M C φ and M C ψ belong to M 1 = {Θ ∈ P(S d-1
)} the set of all probability measures of the unit sphere S d-1 . It is worthy to note that for small regularizing (φ, ψ)-admissible constants, the measure Θ is forced to distribute more weights to directions that are far from each other in terms of their angles (Nguyen et al., 2020). Now, in order to guarantee the existence of (φ, ψ)-admissible constants, we assume that the couple(φ, ψ)-embeddings are approximately angle preserving.

Assumption 1 (Approximately angle preserving property) For any couple (φ, ψ)-embeddings , assume that there exists two non-negative constants L φ and L ψ , such that the following holds

| φ(θ), φ(θ ) | ≤ L φ | θ, θ | and | ψ(θ), ψ(θ ) | ≤ L ψ | θ, θ |, for all θ, θ ∈ S d-1 .
In Proposition 1, we deliver lower bounds of the regularizing (φ, ψ)-admissible constants C φ and C ψ , depending on the dimension of the latent space d and on the levels (L φ , L ψ ) of approximately angle preserving property. These bounds ensure the non-emptiness of the sets M C φ and M C ψ .

Proposition 1 Let Assumption 1 hold and consider regularizing (φ, ψ)-admissible constants such that

C φ ≥ L φ Γ(d/2) √ πΓ((d+1)/2) and C ψ ≥ L ψ Γ(d/2) √ πΓ((d+1)/2) .
Then the sets M C φ and M C ψ contain the uniform measure σ d-1 and σ = d k=1 1 d δ θ k , where {θ 1 , . . . , θ d } forms any orthonormal basis in R d . Note that by Gautschi's inequality (Gautschi, 1959) for the Gamma function, we have that

C φ ≥ L φ d and C ψ ≥ L ψ d .
Proof of Proposition 1 is presented in Appendix A.2. Together the admissible constants, the levels of angle preserving property, and the dimension d of the latent space form the hyperparameters set of HWD problem. For settings of large d, the admissible constants could take smaller values, that force the measure Θ to focus on far-angle directions. However, for smaller d, we may lose the control on the distributional part, the set M C tends to M 1 the entire set of probability measure on S d-1 , hence it boils down on a standard slicing approach that needs an expensive number of projections to get an accurate approximation. Next, we give a set of interesting theoretical properties characterizing HWD.

Proposition 2 For any r ≥ 1, HWD satisfies the following properties:

(i) HWD r (µ, ν) is finite, that is HWD r (µ, ν) ≤ 2 r-1 r (M r (µ) + M r (ν)) where M r (•) is the r-th moment of the given distribution, i.e. M r (µ) = ´X x r dµ(x) 1 r .
(ii) HWD r (µ, ν) is non-negative, symmetric and verifies HWD r (µ, µ) = 0.

(iii) HWD r (µ, ν) has a discrepancy equivalence given by

1 d 1 r inf φ,ψ max θ∈S d-1 W r (µ φ,θ , ν ψ,θ ) ≤ HWD r (µ, ν) ≤ inf φ,ψ max θ∈S d-1 W r (µ φ,θ , ν ψ,θ ).
(iv) For p = q, HWD is upper bounded by the distributional sliced Wasserstein distance.

(v) HWD r is rotation invariant, namely, HWD r (R#µ, Q#ν) = HWD r (µ, ν), for any

R ∈ O p = {R ∈ R p×p : R R = I p } and Q ∈ O q = {Q ∈ R q×q : Q Q = I q },
the orthogonal group of rotations of order p and q, respectively.

(vi) Let T α and T β be the translations from R p into R p and from R q into R q with vectors α and β, respectively. Then

HWD r (T α #µ, T β #ν) ≤ 2 r-1 HWD r (µ, ν) + α + β ) .
Proof of Proposition 2 is given in Apprendix A.1. From property (i) HWD is finite provided that the distributions in question have a finite r-th moments. Note that the supremum over the probability measure sets M C φ and M C ψ guarantees the property HWD r (µ, µ) = 0. For p = q, if the infimum over the couple (φ, ψ)-embedding in (iii) is realized in the identity mappings, then HWD verifies a metric equivalence with respect to the max-sliced Wasserstein distance (Deshpande et al., 2019).

The property (v) highlights a rotation invariance of HWD, which is well verified by the GW distance.

ALGORITHM

Computing HWD requires a resolution of an optimization problem as given in (5). In what follows, we propose an algorithm for computing an approximation of this discrepancy based on samples X = {x i } n i=1 from µ and samples Y = {y j } m j=1 from ν. At first, let us note that we have min-max optimization to solve; the minimization occuring over the embeddings φ and ψ and the maximization over the distributions on the unit-sphere. This maximization problem is challenging due to both the constraints and because we optimize over distributions. Similarly to Nguyen et al. (2020), we approximate the problem by replacing the constraints with regularization terms and by replacing the optimization over distributions by an optimization over a push-forward of the uniform probability measure σ d-1 by a Borel measurable function f : S d-1 → S d-1 . Hence, assuming that we have drawn from a uniform distribution K directions {θ k } K k=1 , the numerical approximation of HWD is obtained by solving the following problem:

min φ,ψ max f L1 def. = 1 K K k=1 W r r X φ[f (θ k )], Y ψ[f (θ k )] 1/r + min f λC L2 def. = k,k φ[f (θ k )] φ[f (θ k )] + k,k ψ[f (θ k )] ψ[f (θ k )] + min f λa L3 def. = k,,k φ[f (θ k )] φ[f (θ k )] -θ k θ k 2 + k,k ψ[f (θ k )] ψ[f (θ k )] -θ k θ k 2
where the first term in the optimization is related to the sliced Wasserstein, the second term is related to the regularization term associated to

E θ,θ ∼Θ [| φ(θ), φ(θ ) |], and E θ,θ ∼Θ [| ψ(θ), ψ(θ ) |],
and the third term is the angle-preserving regularization term. Note that the min term with respect to f is due to the fact that we want those regularizers to be small. Two hyperparameters λ C and λ a control the impact of these two regularization terms. In practice, φ, ψ and f are parametrized as deep neural networks and the min-max problem is solved by an alternating optimization scheme : (a) optimizing over f with ψ and φ fixed then (b) optimizing over ψ and φ with f fixed. Some details of the algorithms are provided in Algorithm 1.

Regarding computational complexity, if we assume that the mapping f, φ, ψ are already trained, that we have K projections, and that φ[f (θ k )], ψ[f (θ k )] are precomputed, then the computation of HWD (line 25 of Algorithm 1) is in O(K(n log n + np + nq)), where n is the number of samples in X and Y . When taking into account the full optimization process, then the complexity depends on the number of times we compute the full objective function we are optimizing. Each evaluation requires the computation of the sum in L 1 which is O(K(n log n + np + nq)) and the two regularization terms L 2 and L 3 require both O(K 2 (p + q + 2d)). Note that in terms of computational complexity, SGW is O(Kn log n) whereas HWD is O(T N Kn log n), with T × N being the global number of objective function evaluations. Hence, complexity is in favor of SGW. However, one should note that in practice, because we optimize over the distribution of the random projections, we usually need less slices than SGW and thus depending on the problem, T N K can be of the same magnitude than the number of slices involved in SGW (similar findings have been highlighted for Sliced Wasserstein distance (Nguyen et al., 2020)).

NUMERICAL EXPERIMENTS

In this section, we analyze HWD, exhibit its rotation-invariant property, and compare its performance with SGW in a generative model context.

Translation and Rotation

We have used two simple datasets for showing the behavior of HWD with respect to translation and rotation. For translation, we consider two 2D Gaussian distributions one being fixed, the other with varying mean. For rotation, we use two 2D spirals from Scikit-Learn library, one being fixed and the other being rotated from 0 to π/2. For these two cases, we have drawn 500 samples, used 100 random directions for SGW and RI-SGW. For our HWD, we have used only 10 slices and T = 50, N = 5 iterations for each of the alternate optimization. The results we obtain are depicted in Figure 2. From the first panel, we remark that both SGW and RI-SGW are indeed insensitive to translation while HWD captures this translation, which is also verified by property (vi) in Proposition 2. For the spiral problem, as expected HWD and RI-SGW are indeed rotation-invariant while SGW is not.

Generative models For checking whether our distribution discrepancy behaves appropriately, we have used it as a loss function in a generative model. Our task here is to build a model able to generate a distribution defined on a space having a different dimensionality from the target distribution space. As such, we have considered the same toy problems as in Bunne et al. (2019) and investigated two situations: generating 2D distributions from 3D data and the other way around. The 3D target distribution is a Gaussian mixture model with four modes while the 2D ones are 5-mode. Our generative model is composed of a fully-connected neural network with ReLU activation functions.

Algorithm 1 COMPUTING HETEROGENEOUS WASSERSTEIN DISCREPANCY (SEE (5))

1: Input: Source and target samples: (X, µ) and (Y, ν); order r; the set of random direction {θ k } K k=1 ; T number of global iterations; N number of iterations for each of the alternate scheme; 2: Output: HWD

3: function COMPUTE LOSSES(X, Y , φ[f (θ k )] , ψ[f (θ k )]) 4: compute the average of 1D-Wasserstein L 1 between X φ[f (θ k )] and Y ψ[f (θ k )]
for i = 1, • • • , N do 11: compute φ[f (θ k )] and ψ[f (θ k )] 12: L 1 , L 2 , L 3 ← Compute Losses (X, Y , φ[f (θ k )] , ψ[f (θ k )])) 13: L = -L 1 + L 2 + L 3 14: f ← f -γ k ∇L 15:
end for 16: fix f 17: We have considered 3000 samples in the target distributions and batch size of 300. For both HWD and SGW, we have run the algorithm for 30000 iterations with an Adam optimizer, stepsize of 0.001 and default β parameters. For the 4-mode problem, the generator is a MLP with 2 layers while for the 5-mode, as the problem is more complex, it has 3 layers. In each case, we have 256 units on the first layer and then 128. For the hyperparameters, we have set λ C = 1 and λ a = 5 or λ a = 50 depending on the problem. Note that for SGW, we have also added a 2 -norm regularizer on the output of the generator in order to avoid them to drift (see Figures 8 and9 in Appendix C), as the loss is translation-invariant. Examples of generated distributions are depicted in Figure 4. We remark that our HWD is able to produce visually correct distributions whereas SGW struggles in generating the 4 modes and its 3D 5-mode is squeezed on its third dimension. Scalability We consider the non-rigid shape world dataset (Bronstein et al., 2006) which consists of 148 three-dimensional shapes from 12 classes. We draw randomly n ∈ {100, 250, 500, 1000, 1500, 2000} vertices {x i ∈ R 3 } n i=1 on each shape and use them to measure the similarity between a pair of shapes {x i ∈ R 3 } n i=1 and {y j ∈ R 3 } n j=1 . Figure 5 reports the average time to compute on a single core such a similarity for 100 pairs of shapes using respectively GW, SGW and HWD. As expected GW exhibits a slow behavior while the computational burden of HWD is on par with SGW. to different transformations in terms of classification accuracy. To that purpose we employ the Shape Retrieval Contest (SHREC'2010) correspondence dataset, see Bronstein et al. (2010). It includes high resolution (10K-50K) triangular meshes. The shapes are of 3 classes (see Figure 10 in Appendix C) with 9 different transformations and the null shape (no transformation). Each transformation is applied up to five strength levels (weak to strong). Along with the null shape, we consider all strengths of the "isometry", "topology", "scale", "shotnoise" transformations leading to 63 samples. We perform a 1-NN classification.

for i = 1, • • • , N do 18: compute φ[f (θ k )] and ψ[f (θ k )] 19: L 1 , L 2 , L 3 ← Compute Losses (X, Y , φ[f (θ k )] , ψ[f (θ k )]) 20: L = L 1 + L 2 + L 3 21: φ ← φ -γ k ∇L 22: ψ ← ψ -γ k ∇L
Obtained performances over 10 runs are depicted in Figure 6. They highlight the ability of HWD to be robust to perturbations. HWD achieves slightly better mean classification accuracy than SGW with a competitive computation time (see Figure 5). Notice that GW and RISGW are unable to run under reasonable time-budget constraint.

CONCLUSION

We introduce in this paper HWD a novel OT-based discrepancy between distributions lying in different spaces. It takes computational benefits from distributional slicing technique, which amounts to find an optimal number of random projections needed to capture the structure of data distributions. Another feature of this discrepancy consists in projecting the distributions in question through a learning of embeddings enjoying the same latent space. We showed a nice geometrical property verified by the proposed discrepancy, specifically a rotation-invariance. We illustrated through extensive experiments the applicability of this discrepancy on generative modeling and shape objects retrieval. We argue that the implementation part faces the standard deep learning bottleneck of tuning the model's hyperparameters. A future extension line of this work is to deliver theoretical guarantees regarding the regularizing parameters, both of distributional and angle preserving properties.

A PROOFS A.1 PROOF OF PROPOSITION 1

We use the following result:

Lemma 1 [Theorem 3 in (Nguyen et al., 2020)] For uniform measure σ d-1 on the unit sphere S d-1 , we have ˆSd-

1 ×S d-1 | θ, θ |dσ d-1 (θ)dσ d-1 (θ ) = Γ(d/2) √ πΓ((d + 1)/2) . Hence, ˆSd-1 ×S d-1 | φ(θ), φ(θ ) |dσ d-1 (θ)dσ d-1 (θ ) (Assumption 1) ≤ L φ ˆSd-1 ×S d-1 | θ, θ |dσ d-1 (θ)dσ d-1 (θ ) (Lemma 1) ≤ L φ Γ(d/2) √ πΓ((d + 1)/2) .
Therefore, as long as the (φ, ψ)-admissible constants (Gautschi, 1959) for the Gamma function, it yields that

C φ ≥ L φ Γ(d/2) √ πΓ((d+1)/2) and C ψ ≥ L ψ Γ(d/2) √ πΓ((d+1)/2) , we have σ d-1 ∈ M C φ ∩ M C ψ . Now using a Gautschi's inequality
Γ(d/2) √ πΓ((d+1)/2) ≥ 1 √ π(d+1)/2 ≥ 1/d. Let σ = d l=1 1 d δ θ l , where {θ 1 , . . . , θ d } form an orthonormal basis in R d . We then have E θ,θ ∼σ | φ(θ), φ(θ ) | = 1≤k,l≤d 1 d 2 | φ(θ k ), φ(θ l ) | (Assumption 1) ≤ L φ 1≤k,l≤d 1 d 2 | θ k , θ l | = L φ d .
Therefore we get the lower bounds for the (φ, ψ)-admissible constants C φ and C ψ given in Proposition 1, that guarantee

σ d-1 , σ ∈ M C φ ∩ M C ψ .

A.2 PROOF OF PROPOSITION 2

Let us first state the two following lemmas: Lemma 2 writes an integration result using push-forward measures; it relates integrals with respect to a measure η and its push-forward under a measurable map f : X → Y. Lemma 3 proves that the admissible set of couplings between the embedded measures are exactly the embedded of the admissible couplings between the original measures.

Lemma 2 [See Lerner (2014) p. 61] Let f : S → T be a measurable mapping, let η be a measurable measure on S, and let g be a measurable function on T . Then ´T gdf # η = ´S(g • f )dη.

Lemma 3 [Lemma 6 in Paty & Cuturi (2019)] For all φ, ψ and µ ∈ P(X ), ν ∈ P(Y), one has

Π(φ#µ, ψ#ν) = {(φ ⊗ ψ)#γ s.t. γ ∈ Π(µ, ν)},
where φ ⊗ ψ : X × Y → X × Y such that (φ ⊗ ψ(x, y) = (φ(x), ψ(y)) for all x, y ∈ X × Y.

• (i) HWD r (µ, µ) is finite. In one hand, we assume that µ ∈ P r (X ) and ν ∈ P r (Y), hence its r-th

moments are finite, i.e., M r (µ) = ´X x r dµ(x) 1/r < ∞ and M r (ν) = ´Y y r dν(y) 1/r < ∞.
In the other hand, the following holds for all parameter θ ∈ S d-1 and a couple (φ, ψ)-embeddings,

W r r (µ φ,θ , ν ψ,θ ) = inf π∈Π(P φ(θ) #µ,ψ(θ)#ν) ˆR×R |u -u | r dπ(u, u ) (Lemma 3) = inf γ∈Π(µ,ν) ˆX ×Y |φ(θ) x -ψ(θ) y| r dγ(x, y) ≤ 2 r-1 inf γ∈Π(µ,ν) ˆX ×Y |φ(θ) x| r + |ψ(θ) y| r dγ(x, y) = 2 r-1 inf γ∈Π(µ,ν) ˆX |φ(θ) x| r dµ(x) + ˆY |ψ(θ) y| r dν(y) ,
where we use the facts that (s + t) r ≤ 2 r-1 (s r + t r ), ∀s, t ∈ R + and that any γ transport plan has marginals µ on X and ν on Y. By Cauchy-Schwarz inequality, we get

W r r (µ φ,θ , ν ψ,θ ) ≤ 2 r-1 ˆX φ(θ) r x r dµ(x) + ˆY ψ(θ) r y r dν(y) = 2 r-1 M r r (µ) + M r r (ν) .
Then,

E θ∼Θ W r r (µ φ,θ , ν ψ,θ ) 1/r ≤ 2 r-1 r M r r (µ) + M r r (ν) 1/r ≤ 2 r-1 r M r (µ) + M r (ν) . Finally, one has that HWD r (µ, ν) ≤ 2 r-1 r M r (µ) + M r (ν) .
• (ii) Non-negativity and symmetry. Together the non-negativity, symmetry of Wasserstein distance and the decoupling property of iterated infima (or principle of the iterated infima) yield the non-negativity and symmetry of the distributional sliced sub-embedding distance.

• (ii) HWD r (µ, µ) = 0. Let φ and φ two embeddings for projecting the same distribution µ.

Without loss of generality, we suppose that the corresponding (φ, φ )-admissible constants

C φ ≤ C φ , hence M C φ ⊆ M C φ . Using the fact that sup(A ∩ B) ≤ sup A ∧ sup B, (with a ∧ b) = min(a, b)), se have, straightforwardly, HWD r (µ, µ) = inf φ,φ sup Θ∈M C φ ∩M C φ E θ∼Θ W r r (µ φ,θ , µ φ ,θ ) 1 r ≤ inf φ,φ sup Θ∈M C φ E θ∼Θ W r r (µ φ,θ , µ φ ,θ ) 1 r ∧ sup Θ∈M C φ E θ∼Θ W r r (µ φ,θ , µ φ ,θ ) 1 r = inf φ inf φ sup Θ∈M C φ E θ∼Θ W r r (µ φ,θ , µ φ ,θ ) 1 r ∧ sup Θ∈M C φ E θ∼Θ W r r (µ φ,θ , µ φ ,θ ) 1 r ≤ inf φ sup Θ∈M C φ E θ∼Θ W r r (µ φ,θ , µ φ,θ ) 1 r ∧ sup Θ∈M C φ E θ∼Θ W r r (µ φ,θ , µ φ,θ ) 1 r ≤ inf φ sup Θ∈M C φ E θ∼Θ W r r (µ φ,θ , µ φ,θ ) 1 r = 0. • (iii) One has 1 d 1 r inf φ,ψ max θ∈S d-1 W r (µ φ,θ , ν ψ,θ ) ≤ HWD r (µ, ν) ≤ inf φ,ψ max θ∈S d-1 W r (µ φ,θ , ν ψ,θ ). Since M C φ ∩ M C ψ ⊂ M 1 and W r r (µ φ,θ , ν ψ,θ ) ≤ max θ∈S d-1 W r r (µ φ,θ , ν ψ,θ ) we find that sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r (µ φ,θ , ν ψ,θ ) 1 r ≤ sup Θ∈M1 E θ∼Θ W r r (µ φ,θ , ν ψ,θ ) 1 r ≤ max θ∈S d-1 W r r (µ φ,θ , ν ψ,θ ) 1/r ≤ max θ∈S d-1 W r (µ φ,θ , ν ψ,θ ), which entails that HWD r (µ, ν) ≤ inf φ,ψ max θ∈S d-1 W r (µ φ,θ , ν ψ,θ ). Moreover, since the (φ, ψ)- admissible constants C φ and C ψ satisfy C φ ≥ U φ d and C ψ ≥ U ψ d hence σ = d l=1 1 d δ θ l ∈ M C φ ∩ M C ψ , where we set θ 1 = argmax θ∈S d-1 W r (µ φ,θ , ν ψ,θ ). We then obtain HWD r (µ, ν) ≥ inf φ,ψ E θ∼σ W r r (µ φ,θ , ν ψ,θ ) 1 r = inf φ,ψ d l=1 1 d W r r (µ φ,θ l , ν ψ,θ l ) 1 r ≥ 1 d 1/r inf φ,ψ W r (µ φ,θ1 , ν ψ,θ1 ) = 1 d 1/r inf φ,ψ max θ∈S d-1
W r (µ φ,θ , ν ψ,θ ).

• (iv) For p = q, HWD is upper bound by the distributional Wasserstein distance (DSW) . Let us first recall the DSW distance: let C > 0 and set M C = {Θ ∈ P(S d-1 ) :

E θ,θ ∼Θ [| θ, θ |] ≤ C}. DSW r (µ, ν) = sup Θ∈M C E θ∼Θ W r r (µ θ , ν θ ) 1 r .
We have that the case of a identity couple of embeddings, φ = Id, ψ = Id, the probability measure set

M C φ , M C ψ = M C , then it is trivial that HWD r (µ, ν) ≤ DSW r (µ, ν). • (v) Rotation invariance. Note that (R#µ) φ,θ = P φ(θ) #(R#µ) = (P φ(θ) • R)#µ,
and for all x ∈ R p , using the adjoint operator R * , (R * = R -1 ),

(P φ(θ) • R)(x) = φ(θ), R(x) = R * (φ(θ)), x = P R * •φ(θ) (x). Then, (R#µ) φ,θ = (P R * •φ(θ) )#µ. Analogously, one has (Q#ν) ψ,θ = (P Q * •ψ(θ) )#ν. Moreover, M C φ = Θ ∈ P(S d-1 ) : E θ,θ ∼Θ [| φ(θ), φ(θ ) |] = Θ ∈ P(S d-1 ) : E θ,θ ∼Θ [| (R * • φ)(θ), (R * • φ)(θ ) |] = M C R * •φ .
Then we have similarly

M C ψ = M C Q * •ψ . This implies HWD r (R#µ, Q#ν) = inf φ,ψ sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r ((R#µ) φ,θ , (Q#ν) ψ,θ ) 1 r = inf φ,ψ sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r ((P R * •φ(θ) )#µ, P Q * •ψ(θ) )#ν) 1 r = inf φ,ψ sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r (µ R * •φ,θ , µ R * •φ,θ 1 r = inf φ,ψ sup Θ∈M C R * •φ ∩M C Q * •ψ E θ∼Θ W r r (µ R * •φ,θ , ν Q * •φ,θ 1 r = inf φ,ψ sup Θ∈M C R * •φ ∩M C Q * •ψ E θ∼Θ W r r (µ R * •φ,θ , ν Q * •ψ,θ 1 r = inf φ =R * •φ,ψ =Q * •ψ sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r (µ φ ,θ , ν ψ ,θ 1 r = HWD r (µ, ν).
• (vi) Translation quasi-invariance. We have

HWD r (T α #µ, T β #ν) = inf φ,ψ sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r ((T α #µ) φ,θ , (T β #ν) ψ,θ ) 1 r .
By Lemmas 3 and 2 , we have

W r r ((T α #µ) φ,θ , (T β #ν) ψ,θ = inf γ∈Π((Tα#µ) φ,θ ,(T β #ν) ψ,θ )) ˆR2 |u -v| r dγ(u, v) = inf γ∈Π((P φ(θ) •Tα)#µ,(P ψ(θ) •T β )#ν) ˆR2 |u -v| r dγ(u, v) = inf γ∈Π(µ,ν) ˆX ×Y |P φ(θ) • T α (x) -P ψ(θ) • T β )(y)| r dγ(x, y) = inf γ∈Π(µ,ν) ˆX ×Y |(P φ(θ) (x) -P ψ(θ) (y)) + (P φ(θ) (α) -P ψ(θ) (β))| r dγ(x, y) ≤ 2 r-1 inf γ∈Π(µ,ν) ˆX ×Y |(P φ(θ) (x) -P ψ(θ) (y))| r dγ(x, y) + |P φ(θ) (α) -P ψ(θ) (β))| r ≤ 2 r-1 inf γ∈Π(µ,ν) ˆX ×Y |(P φ(θ) (x) -P ψ(θ) (y))| r dγ(x, y) + ( α + β ) r .
Thanks to Minkowski inequality,

sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r ((T α #µ) φ θ , (T β #ν) ψ θ ) 1 r ≤ 2 r-1 sup Θ∈M C φ ∩M C ψ E θ∼Θ inf γ∈Π(µ,ν) ˆX ×Y |(P φ(θ) (x) -P ψ(θ) (y))| r dγ(x, y) 1 r + 2 r-1 ( α + β ) sup Θ∈M C φ ∩M C ψ Θ(S d-1 ) 1 r ≤ 2 r-1 sup Θ∈M C φ ∩M C ψ E θ∼Θ W r r (µ φ,θ , ν ψ,θ ) 1 r + 2 r-1 ( α + β ).
Therefore, we get HWD r (T α #µ, T β #ν) ≤ 2 r-1 HWD r (µ, ν) + 2 r-1 ( α + β ).

B IMPLEMENTATION

This section graphically describes the learning procedure in Algorithm 1. It also provides the training details not exposed in the main body of the paper.

B.1 LEARNING SCHEME

We present in Figure 7 the updated graphics of our approach, highlighting the main components : the distributional part is ensured by a first deep neural network as is each of the mappings. As each of the networks should be learned, we included the part of the loss functions associated with each network (blue fonts correspond to minimization, whereas red fonts correspond to maximization, see Algorithm section).

B.2 TRAINING DETAILS

Our experimental evaluations on shape datasets for scalability contrast GW, SGW and HWD. Regarding or classification under isometry transformations, we additionnally consider RI-SGW. Used hyper-parameters for those experiments are detailed below. Notice that SGW, RI-SGW and HWD rely on K, the number of projections sampled uniformly over the unit sphere. This K may vary from a method to another.

1. SGW: K. 2. RI-SGW: λ RI-SGW , the learning rate and T , the maximal number of iterations for solving 4 over the Stiefeld manifold. 3. HWD: beyond K and the latent space dimension d, it requires the parametrization of φ, ψ and f as deep neural networks and their optimizers. For solving the min-max problem by an alternating optimization scheme we use N inner loops and T number of epochs.

For SGW and RI-SGW we use the code made available by their authors and cite the related reference Vayer et al. (2018) as they require. We use POT toolbox Flamary et al. (2021) to compute GW distance.

Scalability This experiment measures the average running time to compute OT-based distance between two pairs of shapes made of n 3D-vertices. 100 pairs of shapes were considered and n varies in {100, 250, 500, 1000, 1500, 2000}.

We choose K SGW = 1000 (as a default value).

For HWD, the mapping function f is designed as a deep network with 2 dense hidden layers of size 50. Regarding both φ and ψ, they have also the same architecture as f (with adapted input and output layers) but the hidden layers are 10-dimensional. Adam optimizers with default parameters are selected to train them. Finally we consider K HWD = 10, d = 5, T = 50, N = 1 as default values.

Notice also that the regularization parameters λ C and λ a are set to 1.

The used ground cost distance for GW distance is the geodesic distance.

Classification under transformations invariance For this experiment, we consider the same set of hyper-parameters as for Scalability evaluation on shape datasets. Besides, the supplementary competitor RI-SGW was trained by setting K RI-SGW = 1000 = K SGW = 1000, λ RI-SGW = 0.01, T RI-SGW = 500. Notice that due to the high-resolution of the meshes (more than 19K threedimensional vertices), RI-SGW and GW were not able to produce the pairwise-distance matrix used in 1NN classification after several hours. Figure 10: Instances of the shape dataset with null and isometry transformations. The classes are respectively human,dog and horse. For the experiments of Figure 6 we also consider the "topology", "scale", "shotnoise" transformations that respectively amount to deform, to upscale and to add noise to the shapes of each class.

  Figure 2: Examples of distance computation between (left) two-translating Gaussian distributions. (right) two spirals.

Figure 3 :

 3 Figure 3: Examples of target distributions for our generative models (left) 3D 4-mode. (right) 2D 5-mode.

Figure 4 :

 4 Figure 4: Examples of generated distributions across iterations (10, 10000, 20000, and 30000) for two targets.From top to bottom (first-row) HWD for the 4-mode. (second-row) SGW for the 4-mode (third-row) HWD for 3D 5-mode. (fourth-row) SGW for 5-mode. For each row, the last panel shows the evolution of the loss over the 30000 iterations.

Figure 5 :

 5 Figure 5: Computation time with respect to n, the number of vertices on each shape. (Left-panel) Instances of 3D objects. (Right-panel) Running time.

Figure 7 :

 7 Figure 7: The implemented approach. Both the distributional and mappings parts are achieved by deep neural networks. A number K of projections is used to compute 1D-Wasserstein distances.

Figure 8 :

 8 Figure 8: Comparing (top) HWD and (bottom) SGW on generating 2D distributions from 3D target.

Figure 9 :

 9 Figure 9: Comparing (top) HWD and (bottom) SGW on generating 3D distributions from 2D target.

ACKNOWLEDGMENTS

The works of Maxime Bérar, Gilles Gasso and Alain Rakotomamonjy have been supported by the OATMIL ANR-17-CE23-0012 Project of the French National Research Agency (ANR).