# A bipartite ranking approach to the two-sample problem

Stéphan Clémençon, Myrto Limnios, Nicolas Vayatis

► To cite this version:

Stéphan Clémençon, Myrto Limnios, Nicolas Vayatis. A bipartite ranking approach to the two-sample problem. Electronic Journal of Statistics, 2025, 19 (1), pp.2733-2779. 10.1214/25-EJS2392. hal-05409025

HAL Id: hal-05409025

https://telecom-paris.hal.science/hal-05409025v1

Submitted on 16 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A bipartite ranking approach to the two-sample problem

Stephan Clémençon<fnref n="1" />, Myrto Limnios*<fnref n="2" /> and Nicolas Vayatis<fnref n="3" />

**Abstract:** The two-sample problem consists in testing whether two independent samples are drawn from the same (unknown) distribution. Its study in high-dimension is the subject of much attention, especially because the information acquisition processes at work in the Big Data era often involve various poorly controlled sources, leading to datasets possibly exhibiting strong sampling bias. While the efficiency of classic methods relying on computing a discrepancy measure between the empirical distributions of each sample, is negatively impacted by increasing dimensionality, we develop a two-step approach based on statistical learning and an extension of *rank tests*. By dividing the initial samples in two, a bipartite ranking algorithm first learns a real-valued scoring function inducing a preorder on the multivariate space. Then, a rank statistic based on the scores of the remaining observations, tests for differences in distribution. Because the ranking algorithm learns how to map the data onto the real line as the likelihood ratio between the original multivariate distributions, the approach resists to large dimensions (ignoring ranking model bias issues) and preserves the advantages of univariate rank tests. We prove nonasymptotic error bounds based on recent results for two-sample linear rank-processes, and experimentally show how the promoted approach surpasses state-of-the-art methods.

**MSC2020 subject classifications:** Primary 62G10, 68Q32; secondary 60E15, 62C12.

**Keywords and phrases:** Bipartite ranking, nonparametric statistical hypothesis testing, two-sample linear rank statistic/process, two-sample problem.

Received January 2024.

# 1. Introduction

The (nonparametric) statistical hypothesis testing problem referred to as the *two-sample problem* is of central importance in statistics and machine-learning. Based on two independent i.i.d. samples {$X_1, \dots, X_n$} and {$Y_1, \dots, Y_m$} of random variables defined on the same probability space ($\Omega, \mathcal{F}, \mathcal{P}$), valued in the same space $\mathcal{Z}$, usually $\mathbb{R}^d$ with $d \ge 1$, and drawn from (unknown) probability distributions $G$ and $H$ respectively, with $n, m \ge 1$, it consists in testing whether the distributions $G$ and $H$ are the same or not, i.e., in testing the

<footnote n="1">¹Telecom Paris, LTCI, Institut Polytechnique de Paris, 19 place Marguerite Perey, Palaiseau, 91120, France, e-mail: stephan.clemencon@telecom-paris.fr</footnote>

<footnote n="2">²Ecole Polytechnique Fédérale de Lausanne, Institute of Mathematics, Station 8, EPFL, Lausanne, 1015, Switzerland, e-mail: myrto.limnios@epfl.ch</footnote>

<footnote n="3">³Université Paris Saclay, ENS Paris Saclay, CNRS, Centre Borelli, 91190, Gif-sur-Yvette, France, e-mail: nicolas.vayatis@ens-paris-saclay.fr</footnote>
null (composite) hypothesis $\mathcal{H}_0 : G = H$ against the alternative $\mathcal{H}_1 : G \neq H$. Easy to formulate, this generic problem is ubiquitous. It finds applications in many areas, clinical trials in particular, in order to determine whether the fluctuations of a collection of measurements performed over two statistical populations subject to different treatments are simply due to the sampling phenomenon or else to the effect of the treatment. It may also be used in the context of multimodal data fusion, to decide whether two datasets can be pooled or not. Indeed, selection and sampling bias issues are also a major concern in machine learning now. As recently highlighted by theoretical and empirical works (see e.g. Bertail et al. (2021) or Wang et al. (2019)), a poor control on the acquisition process of training data, even massive, can significantly jeopardize the generalization ability of the learned predictive rules. Various dedicated approaches have been introduced in the statistical literature, see e.g. section 6.9 in Lehmann and Romano (2005) or Chapter 3.7 in van der Vaart and Wellner (1996). Many of them consist in computing first nonparametric estimators $\hat{G}_n$ and $\hat{H}_m$ of the underlying distributions (possibly smoothed versions of the empirical distributions, typically), evaluating next a distance or information-theoretic pseudo-metric $\mathcal{D}(\hat{G}_n, \hat{H}_m)$ between the latter in order to measure their dissimilarity (e.g. the two-sample Kolmogorov-Smirnov statistic) and rejecting $\mathcal{H}_0$ for ‘large’ values of the statistic $\mathcal{D}(\hat{G}_n, \hat{H}_m)$, see Biau and Gyorfi (2005) or Gretton et al. (2012) for instance. Beyond computational difficulties and the necessity of identifying a proper standardization in order to make the statistic asymptotically pivotal, i.e., its limit distribution is parameter-free (this generally requires in practice the use of resampling bootstrap techniques), the major issue one faces when trying to implement such plug-in procedures is related to the curse of dimensionality, see Ramdas et al. (2015); Huang et al. (2023) in particular. Indeed, such procedures involve the consistent estimation of distributions on a feature space of possibly very large dimension $d \in \mathbb{N}^*$.

The approach promoted and analyzed in this article is very different and is inspired from *rank tests* in the univariate case, see, e.g., section 6.8 in Lehmann and Romano (2005). Referring to $G$ and $\{X_1, \dots, X_n\}$ as the ‘positive’ distribution and the ‘positive’ sample and to $H$ and $\{Y_1, \dots, Y_m\}$ as the ‘negative’ distribution and the ‘negative’ sample, such methods consist in computing a functional of the ‘positive ranks’, i.e., the ranks of the positive instances among the pooled sample of size $N = n + m$. The rationale behind the use of such statistics, called *two-sample rank statistics*, for univariate two-sample problems naturally lies in the fact that they are pivotal under the null assumption (i.e. ignoring ties, the positive ranks are uniformly distributed on $\{1, \dots, N\}$ under $\mathcal{H}_0$), providing tests with zero bias. The notion of two-sample *linear* rank statistics, namely the (possibly weighted) sum of the images of the positive ranks divided by $(N+1)$ by a score generating function $u \in (0,1) \mapsto \phi(u)$, has been extensively studied in statistics. For an appropriate choice of $\phi$, the related rank test is known to have asymptotic power with optimal properties, describing its capacity to detect alternatives close to the null assumption, see e.g. Hájek and Sidák (1967) or Chapter 15 in van der Vaart (1998). For instance, the popular Mann-Whitney-Wilcoxon ‘ranksum’ statistic, widely used to test whether a distribution is stochastically larger than another one, corresponds to the case $\phi(u) = u$ and is optimal to detect asymptotically small shifts. The methodology we propose here contrasts with straightforward extensions of such techniques to the multivariate framework based on preliminary given *depth functions* (see e.g. Liu and Singh (1993)). We start from the observation that, in the univariate case, two-sample rank statistics are summaries of the statistical version of the ROC curve relative to the pair of
probability measures $(H, G)$. Under $H_0$, the ROC curve coincides with the diagonal of the unit square $[0, 1]^2$. In particular, up to an affine transform, the ranksum statistic is nothing else than the area under the empirical ROC curve (AUC in abbreviated form). The method uses the 2-split trick and is implemented in two steps. It consists in learning first from half of the samples how to rank the multivariate observations in $Z$ by means of a scoring function $s: Z \to \mathbb{R}$ transporting the natural order on $\mathbb{R}$ onto the feature space $Z$ like (any increasing transform of) the likelihood ratio $dG/dH(z)$ would do (i.e. the larger the score $s(Z)$ of an observation $Z$ drawn either from $H$ or $G$, the likelier $Z$ is drawn from $G$ ideally). Next, one applies a rank test to the (univariate) scores of the other half of the samples. The first step of the method thus boils down to solving a bipartite ranking problem, which can be classically formulated as a ROC/AUC optimization problem, see Clémençon and Vayatis (2009a) or Clémençon and Vayatis (2010) for instance. The second step consists in testing whether the obtained empirical ROC curve, significantly deviates from the diagonal, where the score generating function $\phi(u)$ used in the rank test determines the nature of the deviation.

By means of concentration results for two-sample linear rank processes, i.e., collections of two-sample linear rank statistics, recently established in Clémençon, Limnios and Vayatis (2021), the (unbiased) two-stage testing procedure is analyzed from a nonasymptotic perspective. For nonparametric classes of alternative hypotheses, we prove bounds for the type-I and II errors of tests based on a scoring function $s(z)$ maximizing a statistical counterpart of a bipartite ranking performance criterion (summarizing the ROC curve) taking the form of a two-sample linear rank statistic. The capacity of the two-stage method proposed to detect ‘small’ deviations from the null assumption, preserved even in very high dimension to a certain extent, is also thoroughly investigated from an empirical angle. An extensive experimental study is presented, covering a wide variety of two-sample problems and comparing the performance of the ranking-based tests to that of alternative nonparametric methods documented in the literature. It should be highlighted that, if the concept of ROC curve is widely used to evaluate the merits of any univariate test statistic (and find an appropriate trade-off between the two types of error), the approach developed in this article is the first attempt to use this notion to devise testing procedures in a general multivariate framework. Notice finally that a very preliminary version of the two-stage testing method, limited to AUC optimization and ‘ranksum’ test statistics, has been previously outlined in the conference paper Clémençon, Depecker and Vayatis (2009). The present article presents ranking-based tests and analyzes their performance in a much more general framework: building on recent results established in Clémençon, Limnios and Vayatis (2021) a wider class of ranking-based tests are considered, a finite-sample analysis of their power is carried out here and more complete experimental results are presented.

The paper is organized as follows. In section 2, the main notations are set out, the statistical framework of the two-sample problem is recalled at length and traditional methods, rank tests in the univariate case in particular, are briefly reviewed. Insights into the connection between bipartite ranking and the two-sample problem are also provided, together with an explanation of the rationale behind the approach we promote. The novel testing procedure we propose is described, and theoretical guarantees are established in section 3. Experimental results are presented and discussed in section 4, while several concluding remarks are collected in section 5 Technical proofs and details are deferred to the Appendix section.
## 2. Background and preliminaries

In this section, the main notations are introduced and the two-sample problem is formulated in a nonparametric framework. The existing methods to solve it are briefly reviewed, with a particular attention to *rank tests* in the 1 to *d*-dimensional case when $\mathcal{Z} \subset \mathbb{R}^d$, and their interpretation through ROC analysis. Concepts and results related to the bipartite ranking task, viewed as the problem of optimizing (summaries of) the ROC curve are also recalled, insofar as the methodology proposed and analyzed in the subsequent section is based on the latter. Here and throughout, the indicator function of any event $\mathcal{E}$ is denoted by $\mathbb{I}\{\mathcal{E}\}$, the Dirac mass at any point $x$ by $\delta_x$, the generalized inverse of any cumulative distribution function $W(t)$ on $\mathbb{R} \cup \{+\infty\}$ by $W^{-1}(u) = \inf\{t \in (-\infty, +\infty] : W(t) \ge u\}, u \in [0, 1]$. We denote the floor and ceiling functions by $u \in \mathbb{R} \mapsto \lfloor u \rfloor$ and by $u \in \mathbb{R} \mapsto \lceil u \rceil$ respectively. For any bounded function $\psi : (0, 1) \to \mathbb{R}$, we also set $||\psi||_\infty = \sup_{u \in (0,1)} |\psi(u)|$. Throughout the paper, bold symbols refer to multivariate variables, e.g., we write $X$ and $Y$ when considering univariate random variables, while $\mathbf{X}$ and $\mathbf{Y}$ are used when the space which they take their values in the multivariate measurable space $\mathcal{Z}$, usually subset of $\mathbb{R}^d$ with $d \ge 2$ possibly.

### 2.1. The two-sample problem - nonparametric formulation

Let $\{\mathbf{X}_1, \dots, \mathbf{X}_n\}$ and $\{\mathbf{Y}_1, \dots, \mathbf{Y}_m\}$ be independent i.i.d. samples drawn from probability distributions $H$ and $G$ on a measurable space $\mathcal{Z} \subset \mathbb{R}^d$. In the most general version of the two-sample problem, one makes no assumption about the distributions $H$ and $G$ and the goal pursued is to test the composite hypothesis:

$$ \mathcal{H}_0: G = H \quad \text{against the alternative} \quad \mathcal{H}_1: G \neq H, \qquad (1) $$

based on the two samples. A classic approach is to consider a probability (pseudo-) metric $\mathcal{D}$ on the space of probability distributions on $\mathcal{Z}$. Based on the simple observation that $\mathcal{D}(G, H) = 0$ under the null hypothesis, a natural testing procedure consists in computing estimates of the distributions $G$ and $H$, typically the empirical distributions:

$$ \hat{G}_n = \frac{1}{n} \sum_{i=1}^{n} \delta_{\mathbf{X}_i} \quad \text{and} \quad \hat{H}_m = \frac{1}{m} \sum_{j=1}^{m} \delta_{\mathbf{Y}_j}, $$

and rejecting $\mathcal{H}_0$ for ‘large’ values of the statistic $\mathcal{D}(\hat{G}_n, \hat{H}_m)$, see Biau and Gyorfi (2005) for instance. Various metrics or pseudo-metrics can be considered for measuring dissimilarity between two probability distributions. We refer to Rachev (1991) for an excellent account of metrics in spaces of probability measures and their applications. Typical examples include the chi-square distance, the Kullback-Leibler divergence, the Hellinger distance, the Kolmogorov-Smirnov distance and its generalizations of the following type:

$$ \text{MMD}(G, H) = \sup_{f \in \mathcal{F}} \left| \int_{X \in \mathcal{Z}} f(x)G(dx) - \int_{X \in \mathcal{Z}} f(x)H(dx) \right| \qquad (2) $$

where $\mathcal{F}$ denotes a supposedly rich enough class of functions $f: \mathcal{Z} \to \mathbb{R}$, so that $\text{MMD}(G, H) = 0$ if and only if $G = H$, see Th. 5 in Gretton et al. (2012). In the univariate case, when $\mathcal{F}$ is the set of half lines, the quantity (2), referred to as *Maximum Mean
Discrepancy if $\mathcal{F}$ is a particular Reproducing Kernel Hilbert Space (RKHS), naturally reduces to the usual Kolmogorov-Smirnov statistic. Beyond computational difficulties and the necessity of identifying a proper standardization in order to make the test statistic (2) asymptotically pivotal, i.e., its limit distribution is parameter-free, and determining an appropriate critical threshold (this generally requires in practice the use of bootstrap techniques), the major issue one faces when trying to implement such plug-in procedures is related to the curse of dimensionality. Indeed, such procedures involve the consistent estimation of distributions on a feature space of possibly very large dimension $d \in \mathbb{N}^*$. This difficulty can however be circumvented to a certain extent when a unit ball of a RKHS $\mathcal{H}$ is chosen for $\mathcal{F}$ in order to allow for efficient computation of the supremum (2), see Gretton et al. (2007, 2012) and Bach, Harchaoui and Moulines (2008), with recent advances in Schrab et al. (2023). In a broader scope, the two-sample problem formulates also the problem of testing the independence between two random samples, for which numerous methods have also been developed, see e.g. covariance-based distances Székely and Rizzo (2007); Székely and Rizzo (2013), later generalized to metric spaces by Lyons (2013); Jakobsen (2017). In particular, an important extension of MMD-type of distances defined as Hilbert-Schmidt Independence Criterion (HSIC) was introduced and analyzed by Gretton et al. (2005a,b, 2007), and later showed to encompass covariance-based distances for measuring independence in Sejdinovic et al. (2013). Ramdas et al. (2015) provided empirical evidence on the negative impact of the increase of dimension $d$ on the empirical power of kernel-based tests however, recently investigated and theoretically confirmed by Huang et al. (2023). Other measures for testing independence have been proposed that go beyond the purpose of this paper. Lastly, a related problem in computer science literature refers to the two-sample problem as property testing, see for instance Rubinfeld and Sudan (1996); Goldreich, Goldwasser and Ron (1998). The methodology promoted in the present paper for testing (1) is very different in nature and is inspired from traditional techniques based on the notion of rank statistics in the particular one-dimensional case recalled below for clarity.

## 2.2. The univariate case - rank tests and ROC analysis

A classic approach to the two-sample problem in the one-dimensional setup consists in ranking the observed data using the natural order on $\mathbb{R}$, and taking the decision depending on the ranks of the positive instances among the pooled sample:

$$ \forall i \in \{1, \dots, n\}, \quad \text{Rank}(X_i) = N \hat{F}_N(X_i), $$

where $\hat{F}_N(t) = (1/N)(\sum_{i=1}^n \mathbb{I}\{X_i \le t\} + \sum_{j=1}^m \mathbb{I}\{Y_j \le t\})$ and $N = n + m$. Assuming that the distributions $G$ and $H$ are continuous for simplicity (the ties occur with probability zero), the idea underlying rank tests lies in the simple fact that, under the null hypothesis $H_0$, the ranks of positive instances are distribution-free, uniformly distributed over $\{1, \dots, N\}$. A popular choice is to consider the sum of ‘positive ranks’, leading to the well-known rank-sum Mann-Whitney-Wilcoxon statistic, see Wilcoxon (1945),

$$ \hat{W}_{n,m} = \sum_{i=1}^{n} \text{Rank}(X_i). \qquad (3) $$
It is widely used to test $H_0$ against the alternative stipulating that one of the two distributions is stochastically larger than the other one. In the situation where $G(dt)$ is stochastically larger<fnref n="1" /> than $H(dt)$, i.e., when $H(t) \ge G(t)$ for all $t \in \mathbb{R}$, the test statistic (3) is expected to take ‘large’ values. Tables for the distribution of the statistics $\hat{W}_{n,m}$ under $H_0$ being available (even in the case where some observations are tied, by assigning the mean rank to ties, see Cheung and Klotz (1997)), no asymptotic approximation result is thus needed for building a test at an appropriate level. In the case where the two c.d.f. are linked by the relationship $G(t) = H(t - \theta)$ with $\theta \ge 0$ (i.e. when the treatment effect is modeled as additive such that one tests $H_0: \theta = 0$), the test statistic (3) is asymptotically uniformly most powerful in the limit experiment when $\theta$ tends to 0 at rate $O(1/\sqrt{N})$ as $n, m$ simultaneously tend to infinity s.t. $n/N \to p$ for some $0 < p < 1$, see e.g. Corollary 15.14 in section 15.5 of van der Vaart (1998). Other functionals of the ‘positive ranks’ can be used as test statistics for the two-sample problem. In particular, the class of two-sample linear rank statistics defined below forms a rich collection of functionals.

**Definition 2.1.** Let $\phi: [0, 1] \to [0, 1]$ be a non-decreasing function. The two-sample linear rank statistic with score-generating function $\phi(u)$ based on the random samples $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$ is given by

$$ \hat{W}_{n,m}^{\phi} = \sum_{i=1}^{n} \phi \left( \frac{\text{Rank}(X_i)}{N+1} \right). \qquad (4) $$

For $\phi(u) = u$, the statistic (4) coincides with $\hat{W}_{n,m}/(N+1)$. Under $H_0$, the statistics (4) are all distribution-free, which makes them particularly useful to detect differences between the distributions $G$ and $H$. They can be used to design unbiased tests at chosen levels $\alpha \in (0, 1)$ by tabulating their distribution under the null assumption. The choice of the score-generating function $\phi$ can be guided by the type of difference between the two distributions (e.g. in scale, in location) one possibly expects, and may then lead to locally most powerful testing procedures, capable of detecting certain types of ‘small’ deviations from $H_0$. Refer to Chapter 9 in Serfling (1980) or to Chapter 13 in van der Vaart (1998) for an account of the asymptotic theory of rank statistics. For instance, choosing $\phi(u) = u \times \mathbb{I}\{u \ge u_0\}$ for $u_0 \in (0, 1)$ or $\phi(u) = u^q$ with $q > 1$ enhances the role played by the highest ranks, see Clémençon and Vayatis (2007) or Rudin (2006). We also emphasize that concentration properties of two-sample linear rank processes, i.e., collections of two-sample linear rank statistics, have recently been studied in Clémençon, Limnios and Vayatis (2021), motivated by the interpretation of (4) as a scalar statistical summary of the ROC curve relative to the pair $(H, G)$. Based on these results, an exponential tail probability bound for (4) is proved in the Appendix section (see Theorem B.2 therein), which is used in the theoretical analysis carried out in section 3.

<footnote n="1">¹Given two distribution functions $H(dt)$ and $G(dt)$ on $\mathbb{R} \cup \{+\infty\}$, it is said that $G(dt)$ is stochastically larger than $H(dt)$ if and only if for any $t \in \mathbb{R}$, we have $G(t) \le H(t)$. We then write: $H \le_{sto} G$. Classically, a necessary and sufficient condition for $G$ to be stochastically larger than $H$ is the existence of a coupling $(X, Y)$ of $(G, H)$, i.e., a pair of random variables defined on the same probability space with first and second marginals equal to $H$ and $G$ respectively, such that $X \le Y$ with probability one.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig 1. Examples of pairs of distributions and their related ROC curves. The 'negative' distribution H is represented in blue and three typical examples of 'positive' distributions with their associated ROC curves are in red (location model), brown (location and scale model) and green (mixture model).</figcaption></figure>

**ROC analysis** The ROC curve is a gold standard tool to describe the dissimilarity between two univariate probability distributions G and H. This criterion of functional nature, ROC<sub>H,G</sub>, can be defined as the Probability-Probability plot: $t \in \mathbb{R} \mapsto (1 - H(t), 1 - G(t))$, where possible jumps are connected by line segments, ensuring that the resulting curve is continuous. With this convention, one may then see the ROC curve related to the pair of d.f. (H, G) as the graph of a càd-làg (i.e. right-continuous and left-limited) non-decreasing mapping valued in [0, 1], defined by $\alpha \in (0, 1) \mapsto 1 - G \circ H^{-1}(1 - \alpha)$ at points $\alpha$ such that $G \circ H^{-1}(1 - \alpha) = 1 - \alpha$. Denoting by $\mathcal{Z}_H$ and $\mathcal{Z}_G$ the supports of H and G respectively, observe that it connects the point $(0, 1 - G(\mathcal{Z}_H))$ to $(H(\mathcal{Z}_G), 1)$ in the unit square $[0, 1]^2$ and that, in absence of plateau (which we assume here for simplicity, rather than restricting the feature space to G's support), the curve $\alpha \in (0, 1) \mapsto \text{ROC}_{G,H}(\alpha)$ is the image of $\alpha \in (0, 1) \mapsto \text{ROC}_{H,G}(\alpha)$ by the reflection with the main diagonal of the Euclidean plane (i.e. the line of equation $\beta = \alpha$) as axis. Notice that the curve ROC<sub>H,G</sub> coincides with the main diagonal of $[0, 1]^2$ iff. $H = G$, i.e., $\mathcal{H}_0$ is true. Hence, the concept of ROC curve offers a visual tool to examine the differences between two distributions in a pivotal manner, see Fig. 1. For instance, the univariate distribution $G(dt)$ is stochastically larger than $H(dt)$ iff. the curve ROC<sub>H,G</sub> is everywhere above the main diagonal and ROC<sub>H,G</sub> coincides with the left upper corner of the unit square iff. the essential supremum of the distribution H is smaller than the essential infimum of the distribution G.

Based on the samples $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$, a statistical version of ROC<sub>H,G</sub> is obtained by computing the empirical c.d.f. $\hat{H}_m(t) = (1/m) \sum_{j=1}^m I\{Y_j \le t\}$ and $\hat{G}_n(t) = (1/n) \sum_{i=1}^n I\{X_i \le t\}$ for $t \in \mathbb{R}$, and plotting the empirical ROC curve:

$$ \widehat{\text{ROC}}_{H,G} = \text{ROC}_{\hat{H}_m, \hat{G}_n} \qquad (5) $$

Observe that the ROC curve (5) is fully determined by the set of ranks occupied by the positive instances within the pooled sample $\{\text{Rank}(X_i) : i = 1, \dots, n\}$. Breakpoints of the piecewise linear curve (5) necessarily belong to the set of gridpoints:

$$ \{(j/m, i/n) : j \in \{1, \dots, m-1\}, i \in \{1, \dots, n-1\} \} . $$
Denote by $X_{(i)}$ the order statistics related to the sample $\{X_1, \dots, X_n\}$, i.e., $\text{Rank}(X_{(n)}) > \dots > \text{Rank}(X_{(1)})$, the ROC curve (5) is the continuous piecewise linear curve connecting the jump points of the step curve:

$$ \alpha \in [0, 1] \mapsto \sum_{j=1}^{m} \hat{\gamma}_j \times \mathbb{I}\{\alpha \in [(j-1)/m, j/m)\}, \quad (6) $$

where, for all $j \in \{1, \dots, m\}$, we set

$$ \hat{\gamma}_j = \frac{1}{n} \sum_{i=1}^{n} \mathbb{I}\{j \ge N - \text{Rank}(X_{(n-i+1)}) - i + 2\}. $$

The empirical ROC curve (5) is expressed as a function of the Rank($X_i$)'s only. As a consequence, any of its summary is a two-sample rank statistic, that is a measurable function of the 'positive ranks'. In particular, choosing the empirical counterpart of the Area Under the ROC Curve (AUC in short) is the popular ROC curve summary, given by

$$ \text{AUC}_{H,G} := \int_{0}^{1} \text{ROC}_{H,G}(\alpha) d\alpha = \mathbb{P}\{Y < X\} + \frac{1}{2} \mathbb{P}\{X = Y\}. \quad (7) $$

Indeed, the AUC of the empirical ROC curve (5), also termed the rate of concording pairs, can be easily shown to coincide with the *ranksum* Mann-Whitney-Wilcoxon statistic (3) up to an affine transformation

$$ \hat{W}_{n,m} = nm\text{AUC}_{\hat{H}_m, \hat{G}_n} + \frac{n(n+1)}{2}. $$

More generally, two-sample linear rank statistics (4) related to score generating functions different from $\phi(u) = u$ provide alternative summaries of the empirical ROC curve and measure different ways of deviating from the main diagonal of the unit square, which coincides with $\text{ROC}_{H,H}$ under the null assumption. Notice incidentally that, under $\mathcal{H}_0$, we have

$$ \text{AUC}_{H,G} - \frac{1}{2} = \int_{\alpha=0}^{1} \{\text{ROC}_{H,G}(\alpha) - \alpha\} d\alpha = 0. \quad (8) $$

Like (3), the statistics (4) are pivotal and in order to quantify their fluctuations as the full sample sizes $n$, $m$ increase, the fraction of 'positive'/ 'negative' observations in the pooled dataset must be controlled. Let $p \in (0, 1)$ be the 'theoretical' fraction of positive instances. For $N \ge 1/p$, we suppose that $n = \lfloor pN \rfloor$ and $m = \lceil (1-p)N \rceil = N - n$. Define the mixture probability distribution $F = pG + (1-p)H$. As $N$ tends to infinity, the asymptotic mean of $\hat{W}_{n,m}/n$ is:

$$ \begin{aligned} W_{\phi}(H, G) &= \mathbb{E}[\phi \circ F(X)] \\ &= \frac{1}{p} \int_{0}^{1} \phi(u) du - \frac{1-p}{p} \int_{0}^{1} \phi \left( p(1 - \text{ROC}_{H,G}(\alpha)) + (1-p)(1-\alpha) \right) d\alpha \\ &= \frac{1}{p} \int_{0}^{1} \phi(u) du - \frac{1-p}{p} \int_{0}^{1} \phi \left( p(\alpha - \text{ROC}_{H,G}(\alpha)) + 1 - \alpha \right) d\alpha, \end{aligned} \quad (9) $$

see section 3 in Clémençon, Limnios and Vayatis (2021). In absence of any ambiguity about the pair $(H, G)$ of univariate distributions considered, we write $W_{\phi}$ for the sake of simplicity. Observe that, under $\mathcal{H}_0$, we have

$$ W_{\phi} = \int_{0}^{1} \phi(u) du. $$
**Remark 1.** (On $W_\phi$-criteria and distances in the ROC space) Notice that Eq. (9) shows that the quantity $(W_\phi(H, G) - \int_0^1 \phi(u)du)/(1-p)$ can be seen as a dissimilarity measure between the curve $ROC_{H,G}$ and the diagonal, which is not a distance in the ROC space in general for an arbitrary choice of the score-generating function $\phi(u)$. However, when $G$ is stochastically larger than $H$, it coincides with the classic $L_1$-distance w.r.t. Lebesgue measure on $(0, 1)$ between the ROC curve and the diagonal in the specific case where $\phi(u) = u$, as shown by (8). In addition, when $\phi(u) = u^q$ for $q \ge 1$, it can also be related to the (shifted) $L_q$-distance. In this case, we have:

$$W_{\phi}(H, G) - \frac{1}{q+1} = \frac{1-p}{p} \left( \frac{1}{q+1} - \|1 - \mathrm{Id}_{(0,1)} + p(\mathrm{Id}_{(0,1)} - \mathrm{ROC}_{H,G})\|_{q}^{q} \right),$$

where $|| \cdot ||_q$ denotes the $L_q$-norm w.r.t. Lebesgue measure on $(0, 1)$ and $\mathrm{Id}_{(0,1)}$ the identity function on $(0, 1)$. Observe finally that, by virtue of the finite increment theorem, we have:

$$\left| W_{\phi}(H, G) - \int_{0}^{1} \phi(u) du \right| \le (1-p) ||\phi'|_{\infty, (0,1)} ||\mathrm{Id}_{(0,1)} - \mathrm{ROC}_{H,G}||_1,$$

as soon as the (supposedly non-decreasing) score-generating function $\phi$ is of class $C^1$ on $(0, 1)$.

**Building a two-sample test in the ROC space** Under $H_0$, the theoretical ROC curve coincides with the main diagonal of $[0, 1]^2$: $ROC_{H,H}(\alpha) = \alpha$, for all $\alpha \in (0, 1)$ and any distribution $H$ on $\mathbb{R}$. In addition, since the empirical ROC curve is itself a function of the ranks, it is also a (functional) pivotal statistic under the null assumption. When the probability distribution $H$ is continuous, all the possible empirical ROC curves are then equiprobable under $H_0$. Hence, in the situation where $\binom{N}{n}$ is not too large, since the ensemble $C_{n,m}$ of all possible empirical ROC curves based on positive and negative samples of respective sizes $n$ and $m$ is of cardinality $\binom{N}{n}$, all broken lines included in it can be enumerated (see Fig. 2) and for any $\alpha \in \{i/\binom{N}{n} : i = 1, \dots, \binom{N}{n}\}$, one can build a tolerance/prediction region $R_\alpha \subset C_{n,m}$ of level $\alpha$, i.e., a subset $R_\alpha \subset C_{n,m}$ of cardinality $\alpha\binom{N}{n}$. Then, a test rejecting $H_0$ when the empirical ROC curve $\widehat{ROC}_{H,G}$ falls outside $R_\alpha$ can be considered. A natural way of building a critical region in the space $D([0, 1])$ of càd-làg mappings $[0, 1] \to [0, 1]$, defining a test of hypothesis $H_0$ at level $\alpha$, is to: fix a pseudo-distance $D$ on $D([0, 1])$, sort the $\binom{N}{n}$ curves in $C_{n,m}$ by increasing distance to the first diagonal, and keep the subset $R_\alpha$ formed by the $\lceil \binom{N}{n} (1-\alpha) \rceil$ curves closest to the diagonal in the sense of the chosen distance $D$. When choosing the distance defined by $L_1$-norm, one naturally recovers the Mann-Whitney-Wilcoxon test. However, many functional distances can be considered for this purpose.

**Remark 2.** (The AUC as a probability metric) One may easily show that

$$\mathrm{AUC}_{H,G} = \frac{1}{2} + \int_{-\infty}^{\infty} (H(t) - G(t)) \mathrm{d}H(t), \qquad (10)$$

see the Appendix section for further details. Hence, when $H$ is stochastically smaller than $G$, the quantity $\mathrm{AUC}_{H,G} - 1/2$ is equal to the $L_1(H)$-distance between the c.d.f. $H(t)$ and $G(t)$.
<figure><img src="image_3.png" /><figcaption>Fig 2. Examples of empirical ROC curves simulations under the null hypothesis.</figcaption></figure>

**Extensions of two-sample rank tests to the multivariate framework** Given the absence of any 'natural order' on $\mathbb{R}^d$ as soon as $d \ge 2$, existing methods explored different concepts of multivariate ranks statistics. Beyond the approaches based on component-wise ranks and copula (see e.g. Puri and Sen (1971) or Lung-Yut-Fong, Lévy-Leduc and Cappé (2015)), provably valid under strong assumptions on the probabilistic model, the concept of statistical depth aims at defining an ordering of the points in the support of a multivariate distribution $P(dx)$ on $\mathbb{R}^d$ in order to emulate ranks, see Mosler (2013). Precisely, a depth function relative to $P$ is a bounded non-negative Borel-measurable mapping $D_P: \mathbb{R}^d \to \mathbb{R}_+$ that defines a preorder for multivariate points in $\mathbb{R}^d$ and hopefully determines the centrality of any point $x \in \mathbb{R}^d$ with respect to the probability measure $P$. Thus, points $x \in \mathbb{R}^d$ near the 'center' of the mass are the deepest, i.e., $D_P(x)$ is among the highest values taken by the depth function. Originally introduced in the seminal contribution Tukey (1975), the half-space depth of $x$ in $\mathbb{R}^d$ relative to $P$, is the minimum of the mass $P(S)$ taken over all closed half-spaces $S \subset \mathbb{R}^d$ such that $x \in S$. Many alternatives have been developed, refer to Liu (1990); Koshevoy and Mosler (1997); Chaudhuri (1996); Oja (1983); Vardi and Zhang (2000); Chernozhukov et al. (2017); Beirlant et al. (2020) for instance.

In Zuo and Serfling (2000), an axiomatic nomenclature of statistical depths has been devised, providing a systematic way of comparing their merits and drawbacks. As proposed by Liu and Singh (1993), assuming that a certain notion of depth has been selected, a natural way of extending two-sample rank tests to the multivariate framework consists in considering the largest of the two samples available, the **X** sample say, using next part of it to compute the sampling version $D_{\hat{S}}$ of the depth relative to $S$. Then, apply a univariate two-sample rank test to the (univariate) sample formed by the depth values of the $X_i$'s, that have not been involved in the depth estimation step and that formed by the depth values of the $Y_j$'s. The main limitation of depth-based rank tests naturally lies in the impact of the notion of depth considered: different choices highlight different features of the distributions, leading to possibly different decisions. We further refer to Hallin et al. (2021), section 1.2 therein, for a review on the vast literature dedicated to the concept of spatial ranks (and signs), see e.g. Möttönen and Oja (1995); Möttönen, Oja and Tienari (1997); Möttönen, Oja and Serfling (2005). Such models are usually held for testing classic parametric alternatives to the homogeneity hypothesis $H_0$ (e.g. location, scale). We refer to Oja (2010) for a comprehensive review of these approaches and to Chakraborty and Chaudhuri (2015, 2017) when considering high-dimensional settings.
A new extension of rank statistics to multivariate settings for hypothesis testing under nonparametric alternatives, is based on measure transport theoretic formulations. Briefly, by learning the empirical mapping minimizing the optimal transport of any absolute continuous distribution to: the spherical uniform distribution on the *d*-dimensional unit ball defines the center-outward ranks (Hallin et al. (2021)), or to the *d*-dimensional unit square defines the rank maps (Deb and Sen (2021)). Either of the latter multivariate rank statistics, are then plugged into classical multivariate extensions of typically distance-based test statistics. Refer to Hallin et al. (2021), subsection 2.1 and Appendix A.2, for a comprehensive understanding of homogeneity tests based on center-outward distributions. These methods are the main modern nonparametric rank-based methods for the two-sample problem. Importantly, they are the only nonparametric approaches being distribution-free under alternatives so far, in addition to important (asymptotic) properties, see Hallin et al. (2021).

Formulations of the two-sample problem for testing statistical independence between two random samples resulted in important contributions. We specifically refer to Shi et al. (2022) defining *generalized symmetric (test) statistics*, corresponding to plugging into classic statistics a center-outward generalization of rank statistics. It encompasses rank-based works such as Deb and Sen (2021); Leung and Drton (2018) but also HSIC-based tests such as Sejdinovic et al. (2013); Gretton et al. (2007), see Remark 2.1 in particular. However, and as already pointed out in section 2.1, the power of the proposed methods decreases when the dimension *d* of the feature space increases. The related test statistics take the form of *U*-statistics with kernel naturally depending on *d*. Huang et al. (2023) proved that the asymptotic distribution has different regimes, defined by a ratio that depends on the underlying distributions (unknown in our testing setting under the alternative *H*₁). The asymptotic distribution thus leads to erroneous rejection threshold. For the sake of completeness, we finally mention the works devoted to (semi-) parametric ranks based on the Mahalanobis statistical distance, applicable to the class of elliptical distributions only, see Um and Randles (1998); Hallin and Paindaveine (2002a,b, 2008) for instance.

As shall be seen, the approach sketched in the next subsection, and analyzed at length in section 3, shares some similarities with the method relying on statistical depth, except that the mapping used to ‘project’ the multivariate observations onto the real line is specifically learned from the data in order to detect at best the deviations in distribution between the two samples. The following subsection precisely explains how any scoring function $s: \mathcal{Z} \to \mathbb{R}$, solution of the bipartite ranking problem related to the pair $(H, G)$, permits to extend the use of ROC analysis and two-sample rank statistics to the two-sample problem in the multivariate setup.

## 2.3. Bipartite ranking - the rationale behind our approach

The goal of bipartite ranking is to learn, based on the ‘positive’ and ‘negative’ samples {$X_1, \dots, X_n$} and {$Y_1, \dots, Y_m$}, how to score any new observations $Z_1, \dots, Z_k$, being each either ‘positive’ or else ‘negative’, that is to say drawn either from $G$ or else from $H$, without prior knowledge, so that positive instances are mostly at the top of the resulting list with large probability. A natural way of defining a total preorder<fnref n="2" /> on $\mathcal{Z}$ is to map it

<footnote n="2">²A preorder $\ll$ on a set $\mathcal{Z}$ is a reflexive and transitive binary relation on $\mathcal{Z}$. It is said to be *total*, when either $z \ll z'$ or else $z' \ll z$ holds true, for all $(z, z') \in \mathcal{Z}^2$.</footnote>
with the natural order on $\mathbb{R} \cup \{+\infty\}$ by means of a scoring rule, i.e., a measurable mapping $s: \mathcal{Z} \to (-\infty, \infty]$. A preorder $\ll_s$ on $\mathcal{Z}$ is then defined by: for all $(x, x') \in \mathcal{Z}$, $x \ll_s x'$ iff. $s(x) \le s(x')$. We denote by $\mathcal{S}$ the set of all scoring functions. The capacity of a candidate $s(z)$ in $\mathcal{S}$ to discriminate between the positive and negative statistical populations is generally evaluated by means of the ROC curve $\text{ROC}_{H_s, G_s}(\alpha) = \text{ROC}(s, \alpha)$, where $H_s$ and $G_s$ the pushforward probability distributions of $H$ and $G$ by the mapping $s(z)$. It offers a visual tool for assessing ranking performance: the closer to the left upper corner of the unit square the curve $\text{ROC}(s, \cdot)$, the better the scoring rule $s$. Therefore, the ROC curve conveys a partial preorder on the set of all scoring functions: for all pairs of scoring functions $s_1$ and $s_2$, one says that $s_2$ is more accurate than $s_1$ when $\text{ROC}(s_1, \alpha) \le \text{ROC}(s_2, \alpha)$ for all $\alpha \in [0, 1]$. It follows from a standard Neyman-Pearson argument that the most accurate scoring rules are increasing transforms of the likelihood ratio $\Psi(z) = dG/dH(z)$. Precisely, Clémençon and Vayatis (2009a), Proposition 4 therein, proved that the optimal set of elements is

$$ \mathcal{S}^* = \{s \in \mathcal{S}, \forall (z, z') \in \mathcal{Z}^2, \Psi(z) < \Psi(z') \Rightarrow s^*(z) < s^*(z')\}. \quad (11) $$

And, for all $(s, \alpha) \in \mathcal{S} \times (0, 1)$:

$$ \text{ROC}(s, \alpha) \le \text{ROC}^*(\alpha), $$

where $\text{ROC}^*(\cdot) = \text{ROC}(\Psi, \cdot) = \text{ROC}(s^*, \cdot)$ for any $s^* \in \mathcal{S}$. Recall that this optimal curve is non-decreasing and concave, thus always above the main diagonal of the unit square (more properties are recalled in Appendix A). Now, the bipartite ranking task can be reformulated in a more quantitative manner: the objective pursued is to build a scoring function $s(z)$, based on the training random samples $\{\mathbf{X}_1, \dots, \mathbf{X}_n\}$ and $\{\mathbf{Y}_1, \dots, \mathbf{Y}_m\}$, with a ROC curve as close as possible to $\text{ROC}^*$. A typical measure for the deviation between the two curves is to consider the distance in sup norm:

$$ d_{\infty}(s, s^*) = \sup_{\alpha \in (0,1)} |\text{ROC}(s, \alpha) - \text{ROC}^*(\alpha)|. \quad (12) $$

This quantity is a distance between ROC curves (or between the related equivalence classes of scoring functions, the ROC curve of any scoring function being invariant by strictly increasing transform) not between the scoring functions themselves. Since the curve $\text{ROC}^*$ is unknown in practice, the major difficulty is that no straightforward statistical counterpart of the (functional) loss (12) is available. Clémençon and Vayatis (2009a) (see also Clémençon and Vayatis (2010)) proved that bipartite ranking can be viewed as a superposition of cost-sensitive classification problems and ‘discretized’ in an adaptive manner, thus applying empirical risk minimization with statistical guarantees in the $d_{\infty}$-sense. The price of that procedure is an additional bias term inherent to the approximation step. Alternatively, the performance of a candidate scoring rule $s$ can be measured by means of the $L_1$-norm in the ROC space. Observing that, in this case, the loss can be decomposed as follows:

$$ d_1(s, s^*) = \int_0^1 |\text{ROC}(s, \alpha) - \text{ROC}^*(\alpha)| d\alpha = \int_0^1 \text{ROC}^*(\alpha) d\alpha - \int_0^1 \text{ROC}(s, \alpha) d\alpha, \quad (13) $$

minimizing the $L_1$-distance to the optimal ROC curve boils down to maximizing the area under the curve $\text{ROC}(s, \cdot)$:

$$ \text{AUC}(s) := \text{AUC}_{H_s, G_s} = \mathbb{P}\{s(\mathbf{Y}) < s(\mathbf{X})\} + \frac{1}{2}\mathbb{P}\{s(\mathbf{Y}) = s(\mathbf{X})\}, $$
where X and Y are random variables defined on the same probability space, independent, with respective distributions G and H, denoting by $G_s$ and $H_s$ the distributions of $s(\mathbf{X})$ and $s(\mathbf{Y})$ respectively. The scalar performance criterion AUC(s) defines a total preorder on $\mathcal{S}$ and its maximal value is denoted by AUC$^* = \text{AUC}(s^*)$, with $s^* \in \mathcal{S}^*$. Bipartite ranking through maximization of empirical versions of the AUC criterion has been studied in several articles, including Agarwal et al. (2005) or Clémençon, Lugosi and Vayatis (2008). Considering a class $\mathcal{S}_0 \subset \mathcal{S}$, upper confidence bounds for the AUC deficit (i.e. the difference between the maximum AUC$^*$ and the AUC) of scoring rules obtained by solving the problem:

$$ \max_{s \in \mathcal{S}_0} \text{AUC}_{\hat{H}_{s,m}, \hat{G}_{s,n}}, \qquad (14) $$

where $\hat{H}_{s,m} = (1/m) \sum_{j=1}^{m} \delta_{s(\mathbf{Y}_j)}$ and $\hat{G}_{s,n} = (1/n) \sum_{i=1}^{n} \delta_{s(\mathbf{X}_i)}$, have been established in particular. Notice that maximizing the empirical version of AUC(s) over a class $\mathcal{S}_0$ of scoring rule candidates boils down to maximizing the rank-sum criterion:

$$ \hat{W}_{n,m}(s) := \sum_{i=1}^{n} \text{Rank}(s(\mathbf{X}_i)), $$

where $\text{Rank}(s(\mathbf{X}_i)) = N\hat{F}_{s,N}(s(\mathbf{X}_i))$ for $i \in \{1, \dots, n\}$ and $\hat{F}_{s,N}(t) = (1/N) \sum_{i=1}^{n} \mathbb{I}\{s(\mathbf{X}_i) \le t\} + (1/N) \sum_{j=1}^{m} \mathbb{I}\{s(\mathbf{Y}_j) \le t\}$ for $t \in \mathbb{R}$. In Clémençon, Limnios and Vayatis (2021) (see also Clémençon and Vayatis (2009b)), the performance of scoring rules $s(z)$ maximizing alternative scalar summaries of the ROC curve, of the form of two-sample linear rank statistics based on the univariate samples $\{s(\mathbf{X}_1), \dots, s(\mathbf{X}_n)\}$ and $\{s(\mathbf{Y}_1), \dots, s(\mathbf{Y}_m)\}$ as well, has been investigated. Precisely, considering a non-decreasing score-generating function $\phi(u)$, the $W_\phi$-ranking performance criterion is defined as:

$$ \begin{aligned} W_{\phi}(s) &= \mathbb{E}[(\phi \circ F_s)(s(\mathbf{X}))] \\ &= \frac{1}{p} \int_{0}^{1} \phi(u) du - \frac{1-p}{p} \int_{0}^{1} \phi(p(1-\text{ROC}(s, \alpha)) + (1-p)(1-\alpha)) d\alpha \quad (15) \end{aligned} $$

where $F_s = pG_s + (1-p)H_s$ for any $s \in \mathcal{S}$. Equipped with this notation, we have $W_\phi(s) \le W_\phi^*(H, G) := W_\phi(s^*)$ for any $s^* \in \mathcal{S}^*$, and when $\phi$ is strictly increasing, the set of maximizers of the criterion $W_\phi$ coincides with $\mathcal{S}^*$, see Proposition 6 in Clémençon, Limnios and Vayatis (2021). For simplicity, we write $W_\phi^*(H, G) = W_\phi^*$, when there is no ambiguity about the pair $(H, G)$ of probability distributions considered. Whereas maximizing the quantity (15) boils down to performing AUC maximization when $\phi(u) = u$, some specific patterns of the preorder induced by a scoring function $s(z)$ can be more or less enhanced depending on the score-generating function $\phi$ chosen. For $\phi(u) = u^q$ with $q > 1$, or any other score generating function $\phi$ that rapidly vanishes near 0 and takes much higher values near 1, for instance, the value of (15) is essentially determined by the behavior of the curve $\text{ROC}(s, \cdot)$ near 0 (i.e. the probability that $s(\mathbf{X})$ takes high values), as discussed in subsection 2.3 of Clémençon, Limnios and Vayatis (2021).

**Remark 3.** (On ranking bipartite vs. binary classification) ROC analysis can be carried out in a situation where H and G have disjoint support: the curve then has vertical or horizontal segments. We point out that the bipartite ranking problem is formulated in a framework that allows the likelihood ratio to be zero (over all or part of the support of G) and/or infinite
(over all or part of the support of $H$). In such a case, the optimal curve ROC\* has horizontal and/or vertical segments. That said, if the distributions have disjoint supports, the scoring problem is really just a classification problem: the (binary) Bayes classifier is then an optimal scoring function. This is why, in the literature, analysis of the performance of scoring rules is often limited to situations where the bipartite ranking problem is quite distinct from binary classification, i.e., situations where $H$ and $G$ are absolutely continuous with respect to each other.

**Ranking-based two-sample rank tests** The two-sample test procedures rely on the observation that deviations of the curve ROC\* from the main diagonal of $[0, 1]^2$, as well as those of $W_{\phi}^*$ from $\int_0^1 \phi(u)du$ for appropriate score generating functions $\phi$, provide a natural way of measuring the dissimilarity between $G$ and $H$ in theory. As revealed by the proposition below, such deviations are equal to zero as soon as the null assumption is fulfilled.

**Proposition 2.2.** The following assertions are equivalent.

(i) The assumption ‘$H_0: H = G$’ holds true.

(ii) The optimal ROC curve relative to the bipartite ranking problem defined by the pair $(H, G)$ coincides with the diagonal of $[0, 1]^2$

$$ \forall \alpha \in (0, 1), \text{ROC}^*(\alpha) = \alpha. $$

(iii) For any score-generating function $\phi(u)$, we have

$$ W_{\phi}^{*} = \int_{0}^{1} \phi(u) du . $$

(iv) There exists a strictly increasing score-generating function $\phi(u)$, such that:

$$ W_{\phi}^{*} = \int_{0}^{1} \phi(u) du . $$

(iv) We have $\text{AUC}^* = 1/2$.

In addition, we have:

$$ \text{AUC}^* - \frac{1}{2} = \mathbb{E}[|\Psi(\mathbf{Y}) - 1|] . \qquad (16) $$

We also recall that the optimal ROC curve related to the pair of distributions $(H, G)$ is the same as that related to the pair of univariate distributions $(H_{s^*}, G_{s^*})$ and that $dG/dH(z) = dG_{s^*}/dH_{s^*}(s^*(z))$ for any $s^* \in S^*$, see Corollary 5 in Clémençon and Vayatis (2009a). Hence, the optimal curve ROC\* is a very natural and exhaustive way of measuring the dissimilarity between two multivariate distributions, extending the basic ROC analysis for univariate distributions defined on $\mathbb{R}$ as recalled in subsection 2.2 and illustrated by the example below.

**Example 1.** (Multivariate Gaussian populations) Consider two Gaussian distributions $H$ and $G$ on $\mathbb{R}^d$ with same positive definite covariance matrix $\Gamma$ and respective means $\theta_-$ and $\theta_+$ in $\mathbb{R}^d$, supposed to be distinct. As an increasing transform of the loglikelihood ratio, the scoring function

$$ z \in \mathbb{R}^d \mapsto s(z) = \langle z, \Gamma^{-1}(\theta_+ - \theta_-) \rangle $$
is optimal, denoting by $\langle \cdot, \cdot \rangle$ the usual Euclidean inner product on $\mathbb{R}^d$. Since it is linear, the pushforward measures $H_s$ and $G_s$ are both univariate Gaussian distributions. Denoting $\Delta(t) = (1/\sqrt{2\pi}) \int_{-\infty}^t \exp(-u^2/2)du$, $t \in \mathbb{R}$, the c.d.f. of the centered standard univariate Gaussian distribution, one may immediately check that the optimal ROC curve is given by, for all $\alpha \in (0, 1)$,

$$ \mathrm{ROC}^*(\alpha) = 1 - \Delta \left( \Delta^{-1}(1-\alpha) - \sqrt{\langle \theta_+ - \theta_-, \Gamma^{-1}(\theta_+ - \theta_-) \rangle} \right). \quad (17) $$

In addition, we have

$$ \mathrm{AUC}^* = 1 - \mathbb{E} \left[ \Delta \left( V - \sqrt{\langle \theta_+ - \theta_-, \Gamma^{-1}(\theta_+ - \theta_-) \rangle} \right) \right], $$

with $V$ a univariate standard Gaussian random variable.

Hence, Proposition 2.2 permits to reformulate the nonparametric test problem (1) in various ways, as follows for instance:

$$ \mathcal{H}_0: \mathrm{AUC}^* = \frac{1}{2} \quad \text{vs.} \quad \mathcal{H}_1: \mathrm{AUC}^* > \frac{1}{2}, \quad (18) $$

or, equivalently, as:

$$ \mathcal{H}_0: W_{\phi}^* = \int_0^1 \phi(u)du \quad \text{vs.} \quad \mathcal{H}_1: W_{\phi}^* > \int_0^1 \phi(u)du, \quad (19) $$

for any given strictly increasing score generating function $\phi(u)$. It is noteworthy that the formulations above are unilateral, in contrast with the classic rank-sum test in the univariate case: indeed, $G^*$ is always stochastically larger than $H^*$, whereas, given two arbitrary univariate probability distributions, one is not necessary stochastically larger than the other.

As the optimal ROC curve and its summaries, such as the quantities $\mathrm{AUC}^*$ or $W_{\phi}^*$, are unknown in practice, we propose an approach for solving the two-sample problem implemented in two steps. By splitting the samples $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$ into two halves: 1) first, solve the bipartite ranking problem based on the first halves of the ‘positive’ and ‘negative’ samples producing a scoring function $\hat{s}(z)$, as described in the preceding subsection; 2) then, perform a univariate rank-based test on the remaining data samples mapped by the obtained scoring function $\hat{s}(z)$. The subsequent sections provide both theoretical and empirical evidence that, beyond the fact that they are nearly unbiased (in the sense that the fixed level $\alpha$ always coincides approximately with their type-I errors), such testing procedures detect very small deviations from the null assumption.

### 3. Ranking-based rank tests for the two-sample problem

In this section, we describe at length the two-sample methodology foreshadowed by the observations made in the preceding section and discuss its possible practical implementations. Its theoretical properties (level and power) are next analyzed from a nonasymptotic angle under specific assumptions.
## Ranking-based Two-sample Rank Tests

**Input.** Two independent and i.i.d. samples $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$ of sizes $n, m \ge 2$ and valued in $\mathcal{Z}$; subsample sizes $n' < n$ and $m' < m$; bipartite ranking $\mathcal{A}$ algorithm operating on the class $\mathcal{S}_0$ of scoring functions defined on $\mathcal{Z}$; univariate two-sample rank test $\Phi_{\alpha}^{\phi}$ of level $\alpha \in (0, 1)$.

**2-split trick.** Divide each of the original samples into two subsamples

$$ \{X_1, \dots, X_{n'}\} \cup \{X_{1+n'}, \dots, X_n\} \text{ and } \{Y_1, \dots, Y_{m'}\} \cup \{Y_{1+m'}, \dots, Y_m\} $$

1. **Bipartite Ranking.** Run the bipartite ranking algorithm $\mathcal{A}$ based on the training data $\mathcal{D}_{n',m'} = \{X_1, \dots, X_{n'}\} \cup \{Y_1, \dots, Y_{m'}\}$, producing the scoring function

$$ \hat{s} = \mathcal{A}(\mathcal{D}_{n',m'}). \qquad (20) $$

2. **Univariate Rank Test.** Form the univariate samples

$$ \{\hat{s}(X_{1+n'}), \dots, \hat{s}(X_n)\} \text{ and } \{\hat{s}(Y_{1+m'}), \dots, \hat{s}(Y_m)\}, $$

the outcome of the test being finally determined by computing the binary quantity

$$ \Phi_{\alpha}(\text{Rank}(\hat{s}(X_{1+n'})), \dots, \text{Rank}(\hat{s}(X_n))), \qquad (21) $$

where $\text{Rank}(t) = \sum_{i=1+n'}^{n} \mathbb{I}\{\hat{s}(X_i) \le t\} + \sum_{j=1+m'}^{m} \mathbb{I}\{\hat{s}(Y_j) \le t\}$.

<caption>Fig 3. Ranking-based two-sample rank test procedure.</caption>

### 3.1. Method and implementations

We explain how the general idea sketched in section 2.3 can be applied effectively, based on the observation of two independent i.i.d. samples $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$ with $n, m \ge 1$. Let $\alpha \in (0, 1)$ be the target level, i.e., the desired type-I error. As previously discussed, two ingredients are essentially involved in the testing procedure:

1. A bipartite ranking algorithm $\mathcal{A}: \mathcal{Z}^{n+m} \to \mathcal{S}_0$ operating on a class $\mathcal{S}_0 \subset \mathcal{S}$ and assigning to any set of training observations $\mathcal{D}_{n,m} = \{x_1, \dots, x_n\} \cup \{y_1, \dots, y_m\}$ a scoring function $\mathcal{A}(\mathcal{D}_{n,m})$ in $\mathcal{S}_0$;

2. A two-sample rank test $\Phi_{\alpha}^{\phi}: (-\infty, +\infty]^{n+m} \to \{0, 1\}$ of level $\alpha \in (0, 1)$ with outcome depending on $\{\text{Rank}(x_1), \dots, \text{Rank}(x_n)\}$ for any pooled univariate dataset $\mathcal{D}_{n,m} = \{x_1, \dots, x_n\} \cup \{y_1, \dots, y_m\}$ and score-generating function $\phi$.

Equipped with these two components, the methodology, relying on the computation of a Ranking-based two-sample rank tests, is implemented in two main steps, as summarized in Fig. 3. Before discussing at length the various forms the general principle described above may take, a few remarks are in order.

**Remark 4. (Bipartite ranking algorithms)** As mentioned in section 2.3, the vast majority of bipartite ranking algorithms documented in the statistical learning literature solve M-estimation problems over specific classes $\mathcal{S}_0$ of scoring functions. The criterion one seeks to maximize is the AUC or a (smoothed / concavified / penalized) variant, such as (15), whose set of optimal elements coincide with a subset of $\mathcal{S}^*$. See for instance Freund et al. (2003), Rakotomamonjy (2004), Rudin et al. (2005), Rudin (2006) or Burges, Ragno and Le (2007). Generalization
results in the form of confidence upper bounds for the deficit of empirical maximizers have been established under various complexity assumptions for $S_0$ in Cl\'emençon, Lugosi and Vayatis (2008), Agarwal et al. (2005), Cl\'emençon and Vayatis (2007), Cl\'emençon, Depecker and Vayatis (2011), Menon and Williamson (2016) and Cl\'emençon, Limnios and Vayatis (2021). Stronger theoretical guarantees (i.e. bounds for the sup-norm deviation (12)) have also been established for alternative approaches, considering ROC optimization as a continuum of cost-sensitive binary classification problems and combining M-estimation with nonlinear approximation methods, see Cl\'emençon and Vayatis (2009a), Cl\'emençon and Vayatis (2010), Cl\'emençon and Vayatis (2009c) or Cl\'emençon, Depecker and Vayatis (2013a).

**Remark 5. (2-split trick)** As recalled above, nearly optimal scoring functions are generally learned by means of M-estimation techniques. Consequently, their dependence on the training observations may be complex and can hardly be explicit in general. For this reason, a 2-split trick is used in order to make the analysis of the fluctuations of the quantity (21) tractable. Hence, conditioned upon the subsamples used in the bipartite ranking step of the procedure, the functional (21) is a two-sample rank statistic. We also underline that the issue of choosing appropriately the subsample dedicated to bipartite ranking and the one used for the rank test is as crucial from a practical perspective as difficult, insofar as it is hard to know in advance the complexity of the ranking task.

We now propose several ways of implementing the methodology summarized in Fig. 3, which will be next studied theoretically in specific situations and whose performance will be empirically investigated at length in section 4.

**Ranking-based two-sample linear rank tests** The simplest implementation consists in considering a test based on a univariate two-sample linear rank statistic (4), characterized by a given score-generating function $\phi$, see Definition 2.1. As recalled in section 2.2, such a statistic is pivotal under the homogeneity assumption $H_0$ in the univariate case. Its null probability distribution can easily be tabulated, even in the case where $(n, m)$ takes very large values given the computing power now at disposal. For all $n, m \ge 1$ and any $\alpha \in (0, 1)$, one may thus determine the quantile:

$$q_{n, m}^{\phi}(\alpha) = \inf_{t \ge 0} \left\{ \mathbb{P}_{H_0} \left\{ \frac{1}{n} \hat{W}_{n,m}^{\phi} - \int_0^1 \phi(u) du \le t \right\} \ge 1 - \alpha \right\}, \quad (22)$$

as well as the critical region:

$$\left\{ \frac{1}{n} \hat{W}_{n,m}^{\phi} > \int_{0}^{1} \phi(u) du + q_{n,m}^{\phi}(\alpha) \right\} \quad (23)$$

occuring with probability less than $\alpha$ under $H_0$ in the univariate case and defining the test at level $\alpha$:

$$\Phi_{\alpha}^{\phi}(\mathcal{D}_{n,m}) = \mathbb{I} \left\{ \frac{1}{n} \hat{W}_{n,m}^{\phi} > \int_{0}^{1} \phi(u) du + q_{n,m}^{\phi}(\alpha) \right\}, \quad (24)$$

based on the univariate samples $\mathcal{D}_{n,m} = \{X_1, \dots, X_n\} \cup \{Y_1, \dots, Y_m\}$. As discussed in section 2.3, only a unilateral test $\Phi_{\alpha}^{\phi}$ is relevant in the multivariate case, given the reformulations (18) or (19) of the two-sample testing problem. This contrasts with the univariate situation for
which no bipartite ranking step is required. Clémençon, Limnios and Vayatis (2021) investigated a natural bipartite ranking approach, consisting in maximizing a statistical version of the performance criterion (15) based on the (multivariate) training data $\mathcal{D}_{n',m'}$ over the class $\mathcal{S}_0$, i.e., in solving the optimization problem:

$$ \max_{s \in \mathcal{S}_0} \hat{W}_{n',m'}^{\phi}(s), \qquad (25) $$

where we set for any scoring function $s(z)$:

$$ \hat{W}_{n',m'}^{\phi}(s) = \sum_{i=1}^{n'} \phi \left( \frac{\mathrm{Rank}(s(\mathbf{X}_i))}{N' + 1} \right), \qquad (26) $$

with $N' = n' + m'$, the quantity Rank($s(\mathbf{X}_i)$) / ($N' + 1$) being a natural empirical counterpart of $F_s(s(\mathbf{X}_i))$ for $i = 1, \dots, n'$. Clémençon, Limnios and Vayatis (2021) studied the generalization capacity of solutions of the problem (25) and (gradient ascent based) optimization strategies for approximately solving (25). Hence, by considering a solution $\hat{s}$ of (25) obtained at Step 1, the test built at Step 2 based on the scored two samples:

$$ \mathcal{D}_{n'',m''}(\hat{s}) = \{\hat{s}(\mathbf{X}_{1+n'}), \dots, \hat{s}(\mathbf{X}_n)\} \cup \{\hat{s}(\mathbf{Y}_{1+m'}), \dots, \hat{s}(\mathbf{Y}_m)\}, \qquad (27) $$

with $n'' = n - n'$ and $m'' = m - m'$, writes

$$ \Phi_{\alpha}^{\phi}(\mathcal{D}_{n'',m''}(\hat{s})) = \mathbb{I} \left\{ \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) > \int_{0}^{1} \phi(u) du + q_{n'',m''}^{\phi}(\alpha) \right\}. \qquad (28) $$

Under specific assumptions, in particular related to the class $\mathcal{S}_0$ and the score-generating function $\phi$, the nonasymptotic properties of the test (28) are investigated in section 3.2.

**Remark 6.** (On numerical aspects of empirical $W_\phi$-maximization) As will be proved in the subsequent theoretical analysis and next empirically shown in Section 4, the approach proposed above is successful (and significantly outperforms its competitors) as soon as the scoring function learned $\hat{s}$ is not too far from optimality. The performance gain in terms of power depends on the success of step 1, which of course has a computational cost. Complexity of empirical $W_\phi$-maximization and related numerical optimization issues have been discussed at length in Clémençon, Limnios and Vayatis (2021), see subsection 4.2 and section 5 therein. More generally, the computational complexity of most bipartite ranking algorithms mentioned in Remark 4 is documented in the literature, see e.g. Clémençon, Depecker and Vayatis (2013b).

**Remark 7.** (Combining multiple ranking-based two-sample linear rank tests) As highlighted in Clémençon, Limnios and Vayatis (2021), depending on the score-generating function $\phi$ chosen, the quantity $W_\phi^*$ summarizes ROC* in a certain fashion. As one hardly knows in advance which $\phi$ may capture best the way ROC* possibly deviates from the diagonal in practice (through the scalar quantity $W_\phi^* - \int_0^1 \phi(u) du \ge 0$), a natural strategy could consist in performing simultaneously several ranking-based two-sample linear rank tests, implementing popular principles in multiple hypothesis testing and ensemble learning. For instance, $K \ge 1$ score generating functions $\phi_1, \dots, \phi_K$ can be considered, together with levels $\alpha_1, \dots, \alpha_K$ to form the ensemble of tests $\{\Phi_{\alpha_k}^{\phi_k} : k = 1, \dots, K\}$ and the combination

$$ \sup_{1 \le k \le K} \Phi_{\alpha_k}^{\phi_k}(\mathcal{D}_{n'',m''}(\hat{s})). \qquad (29) $$
Although one may guarantee that the type I error of (29) is less than $\alpha$ (by choosing the $\alpha_k$'s so that $\sum_{k=1}^K \alpha_k = \alpha$), one faces significant difficulties when investigating its properties, due to the dependence of the tests combined and of the different $\phi_k$. The study of such approaches are thus left for further research.

**Ranking-based tests in the ROC space** As underlined in section 2.2, the empirical ROC curve is itself a two-sample rank statistic (and consequently pivotal) under $H_0$. Hence, the test involved at Step 2 can be based on a confidence region $R_{n'',m''}(\alpha) \subset C_{n'',m''}$ at level $1-\alpha$ for the empirical ROC curve based on univariate samples of sizes $n''$ and $m''$ under $H_0$. Using the scoring function $\hat{s}$ produced at Step 1, one plots the empirical ROC curve related to the univariate distributions:

$$ \hat{H}_{\hat{s},n''} = \frac{1}{n''} \sum_{i=1+n'}^{n} \delta_{\hat{s}(\mathbf{X}_i)} \quad \text{and} \quad \hat{G}_{\hat{s},m''} = \frac{1}{m''} \sum_{j=1+m'}^{n} \delta_{\hat{s}(\mathbf{Y}_j)} \quad (30) $$

and the critical region then writes

$$ \{\text{ROC}_{\hat{H}_{\hat{s},m''}, \hat{G}_{\hat{s},n''} \notin R_{n'',m''}(\alpha)}\}. \quad (31) $$

As pointed out in section 2.2, given a (pseudo-)metric in the ROC space (e.g. the sup norm), the confidence region $R_{n'',m''}(\alpha)$ could naturally correspond to the set of piecewise linear curves in $C_{n'',m''}$ at a distance smaller than a specific threshold $t_\alpha$ from the diagonal (and possibly above the diagonal, just like ROC*). In the case of the $L_1$-distance, this implementation coincides with the previous one when choosing $\phi(u) = u$.

## 3.2. Theoretical guarantees - nonasymptotic error bounds

The properties of the ranking-based two-sample linear rank tests described in section 3.1 are now analyzed from a nonasymptotic perspective. Let $\alpha \in (0, 1)$. The test (28) is of level $\alpha$ by construction as proved by Theorem 3.1 below.

**Theorem 3.1.** (Type-I error bound) Let $\phi(u)$ be a score-generating function and $n, m \ge 2$. Fix $\alpha \in (0, 1)$. Under the null hypothesis $H_0$, the type-I error of the test (28) is less than $\alpha$

$$ \mathbb{P}_{H_0} \left\{ \Phi_{\alpha}^{\phi} (\mathcal{D}_{n'',m''}(\hat{s})) = 1 \right\} \le \alpha, \quad (32) $$

for all $1 \le n'' < n$ and $1 \le m'' < m$.

The proof directly follows from a conditioning argument detailed in Appendix section B. The distribution of (4) can be easily tabulated and the quantile (22) numerically computed for any $n, m \ge 1$, see section 3.1. Consider the following assumption related to the smoothness of the score-generating function $\phi$.

**Assumption 1.** The score-generating function $\phi: [0, 1] \to \mathbb{R}$, is non-decreasing and twice continuously differentiable.

Under Assumption 1, the following result provides an upperbound for (22) that decays to 0, as simultaneously $n, m$ tend to infinity (so that $n/(n+m) \to p$ at the rate $1/(n+m)$ for $p \in (0, 1)$).
**Proposition 3.2.** Suppose Assumption 1 fulfilled. Let $p \in (0, 1)$ and $N \ge 1/p$. Set $n = \lfloor pN \rfloor$ and $m = \lceil (1-p)N \rceil = N - n$. Then, for any $\alpha \in (0, 1)$, we have:

$$q_{n,m}^{\phi}(\alpha) \le \sqrt{\frac{\log(18/\alpha)}{CN}}, \quad (33)$$

where $C = (1/8) \min(p/\|\phi\|_{\infty}^2, 1/(p\|\phi'\|_{\infty}^2), 1/((1-p)\|\phi'\|_{\infty}^2))$.

Refer to Appendix B.6 for the proof. It straightforwardly results from the tail probability bound for two-sample linear rank statistics established in Appendix B, Theorem B.2 under Assumption 1, the bound being of the expected order given the CLT satisfied by (4), refer to e.g. Theorem 13.25 in van der Vaart (1998). Notice that, as pointed out in Remark 2 in Clémençon, Limnios and Vayatis (2021), the validity of the bound can be extended to situations where $\phi$ and its derivatives present a countable number of possible discontinuities on $[0, 1]$ at the cost of much more technicalities (see e.g. Clémençon and Vayatis (2007)). Observe also that the constant $C$ involved in the upper bound in Proposition 3.2 decreases (and the upper bound for the quantile increases) as $\|\phi'\|_{\infty}$ increases. This is to be expected, since the dispersion of the two-sample rank statistic $\hat{W}_{n,m}^{\phi}$ naturally increases with the variability of the score-generating function $\phi$ on which it is based.

In addition, the type-II error can also be controlled using nonasymptotic results proved in Clémençon, Limnios and Vayatis (2021). They provided probability inequalities for the maximal deviations between the empirical criterion (26) and the theoretical one (15), over a class $S_0$ of scoring functions of controlled complexity, and upper confidence bounds for the deficit of $W_\phi$-ranking performance of solutions of problem (25). The latter results rely on linearization techniques applied to the statistic (26) and concentration bounds for two-sample U-processes. As will be shown below and given $\hat{s} \in S_0$, these results establish tail bounds for the quantity

$$\frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) - \int_{0}^{1} \phi(u) du = \left\{ \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) - W_{\phi}(\hat{s}) \right\} + \left\{ W_{\phi}(\hat{s}) - W_{\phi}^{*} \right\} + \left\{ W_{\phi}^{*} - \int_{0}^{1} \phi(u) du \right\}. \quad (34)$$

The (non-negative) third term on the right hand side of Eq. (34) quantifies the deviation from the homogeneity hypothesis $H_0$, while the second one corresponds to the error inherent in the bipartite ranking step (Step 1 in Fig. 3). The following additional assumptions are required to apply those results.

**Assumption 2.** For all $s \in S_0$, the random variables $s(X)$ and $s(Y)$ are absolutely continuous w.r.t the Lebesgue measure.

**Assumption 3.** The class of scoring functions $S_0$ is a VC-class of finite VC-dimension $V < +\infty$.

We refer to section 2.6.2 in van der Vaart and Wellner (1996) for the definition of VC-classes of functions. In Clémençon, Limnios and Vayatis (2021) it is precisely proved, under the assumptions above and when $n'/(n' + m') \sim n''/(n'' + m'') \sim p$, that the deficit of $W_\phi$-ranking performance of solutions of empirical maximizers $\hat{s}$, i.e., the second term on the right hand side of (34), is of order $O_P(1/\sqrt{n' \wedge m'})$ (when neglecting the possible bias model,
arising from the fact that $S_0 \cap S^*$ might be empty, see Corollary 7 therein), and that the first term is of order $O_{\mathbb{P}}(1/\sqrt{n'' \wedge m''})$ (see the argument of Theorem 5 therein). Whereas Assumption 3 is a classical complexity assumption, which it is not surprising to see assumed to obtain concentration results, we point out that Assumption 2 permits to control the variance of the U-processes indexed by $s$ in $S_0$ involved in the decomposition of the deviation process $\{\hat{W}_{n',m'}(s)/n' - W_{\phi}(s)\}_{s \in S_0}$ and to apply concentration inequalities for U-processes (see Lemmas 15 and 16 in Clémençon, Limnios and Vayatis (2021)). Incidentally, we also point out that Assumption 2 is fulfilled as soon as $X$ and $Y$ are continuous random vectors and the scoring functions $s$ in $S_0$ are sufficiently smooth. This includes in particular all the situations considered in Section 4 with linear or quadratic scoring functions.

Considering the quantity $W_{\phi}^* - \int_0^1 \phi(u)du$ to describe the departure from the null assumption $H_0$ (see Proposition 2.2) and the bias model $W_{\phi}^* - \sup_{s \in S_0} W_{\phi}(s)$ inherent in the bipartite ranking step (when formulated as empirical $W_{\phi}$-ranking performance maximization), we introduce the two (nonparametric) classes of pairs of probability distributions on $\mathcal{Z}$.

**Definition 3.3.** Let $\varepsilon > 0$ and $\phi(u)$ be a score-generating function. We denote by $H_1(\varepsilon)$ the set of alternative hypotheses corresponding to all pairs $(H, G)$ of probability distributions on $\mathcal{Z}$ such that

$$W_{\phi}^* - \int_0^1 \phi(u)du \geq \varepsilon, \qquad (35)$$

where we recall $W_{\phi}^* = W_{\phi}(s^*) = W_{\phi}^*(H, G)$ for any $s^* \in S^*$.

**Definition 3.4.** Let $\delta \geq 0$, $\phi(u)$ be a score-generating function and $S_0$ be a class of scoring functions. We denote by $\mathcal{B}(\delta)$ the set of all pairs $(H, G)$ of probability distributions on $\mathcal{Z}$ such that

$$W_{\phi}^* - \sup_{s \in S_0} W_{\phi}(s) \leq \delta. \qquad (36)$$

We point out that the class $H_1(\varepsilon)$ defined above is a large nonparametric collection of pairs of distributions $(H, G)$: reformulating (35), when $\phi(u) = u$ and $\mathcal{Z} = \mathbb{R}^d$ for instance, it includes in particular all pairs of $d$-dimensional Gaussian distributions $(H, G)$ with same positive definite covariance $\Gamma$ with means $\theta_-$ and $\theta_+$ sufficiently far away from each other, i.e., such that $\langle \theta_+ - \theta_-, \Gamma^{-1}(\theta_+ - \theta_-) \rangle \geq \kappa^{-1}(\varepsilon)$, where $\kappa^{-1}$ means the inverse of the function $\kappa: \sigma \in \mathbb{R}_+ \mapsto 1/2 - \mathbb{E}[\Delta(V - \sigma)]$ strictly increasing from 0 to $1/2$, see Example 1. Incidentally, observe that, for such specific pairs of distributions, we also have $W_{\phi}^* = \sup_{s \in S_0} W_{\phi}(s)$ (in other terms, (36) holds with $\delta = 0$), when $S_0$ encompasses the class of all linear scoring functions on $\mathbb{R}^d$.

**Remark 8.** (On departure condition (35)) As discussed in Remark 1, the quantity $W_{\phi}^* - \int_0^1 \phi(u)du$ can be interpreted as a distance between $ROC^* (= ROC_{H,G}^*)$ and the diagonal of the ROC space for very specific choices of the score generating function $\phi$ only. In particular, when $\phi(u) = u$, condition (35) involved in Definition 3.3 precisely means that $ROC^*$ is outside the open $L_1$-ball centered at the diagonal with radius $\varepsilon/(1-p)$. We further refer to Appendix section C.1 for additional analysis when considering other score-generating functions.
The theorem below provides a rate bound for the type-II error of the ranking-based rank test (28) of size $\alpha$, depending on the sizes $N'$ and $N'' = N - N'$ of the two pooled samples involved in the procedure, i.e., that are respectively used for bipartite ranking and for performing the rank test based on the scoring function learned.

**Theorem 3.5.** (Type-II error bound) Let $\phi(u)$ be a score-generating function and $\varepsilon > \delta > 0$. Fix $\alpha \in (0, 1)$. Suppose that Assumptions 1-3 are fulfilled. Let $p \in (0, 1)$ such that $N' \wedge N'' \ge 1/p$. Set $n' = \lfloor pN' \rfloor$ and $m' = \lceil (1-p)N' \rceil = N' - n'$, as well as $n'' = \lfloor pN'' \rfloor$ and $m'' = \lceil (1-p)N'' \rceil = N'' - n''$. Then, there exist constants $C_1$ and $C_2 \ge 24$, depending on $(\phi, \mathcal{V})$, such that the type-II error of the test (28) is uniformly bounded as follows:

$$ \sup_{(H,G) \in \mathcal{H}_1(\varepsilon) \cap \mathcal{B}(\delta)} \mathbb{P}_{H,G} \left\{ \Phi_{\alpha}^{\phi}(\mathcal{D}_{n'',m''}(\hat{s})) = 0 \right\} \le 18 \exp \left\{ -\frac{CN''}{16} (\varepsilon - \delta)^2 \right\} \\ + C_2 \exp \left\{ -\frac{N'}{8C_2} p (p \wedge (1-p)) (\varepsilon - \delta) \log \left( 1 + \frac{\varepsilon - \delta}{32C_1(p \wedge (1-p))} \right) \right\}, \quad (37) $$

as soon as $N'' \ge 4 \log(18/\alpha)/(C(\varepsilon - \delta)^2)$ and $N' \ge 16C_1^2/(p(\varepsilon - \delta)^2)$, the constant C being that involved in Proposition 3.2, the $C_j$'s those involved in Theorem B.5.

Refer to Appendix B.7 for the detailed proof. In view of Definitions 3.3 and 3.4, by requiring $\varepsilon > \delta$, we consider here pairs $(H, G)$ for which the class $S_0$ induces a model bias in the bipartite ranking problem that is sufficiently small, namely smaller than the minimum departure $\varepsilon$ from the null assumption. The bound stated above reveals that, when $\varepsilon - \delta > 0$ is held fixed, the type-II error uniformly vanishes exponentially fast as $N'$ and $N''$ simultaneously increase to infinity: the second term on the right-hand side of the bound (37) is inherent in the learning stage, while the first term corresponds to a bound for the type-II error of a univariate rank test. Hence, for $N'$ large enough, the type-II error of the ranking-based rank test is of the same order of magnitude as that of a rank test based on the univariate samples $\{\Psi(\mathbf{X}_{1+n'})\}, \dots, \Psi(\mathbf{X}_n)\}$ and $\{\Psi(\mathbf{Y}_{1+m'})\}, \dots, \Psi(\mathbf{Y}_m)\}$.

The two bounds exhibit a different behavior as the departure level $\varepsilon$ from the null assumption tends to 0 or to $+\infty$: when the pooled sample sizes $N'$ and $N''$ are held fixed, the bound related to the type-II error of the univariate rank test becomes negligible compared to the learning bound as $\varepsilon$ increases to infinity, while it deteriorates faster than it when $\varepsilon$ decays to 0. Observe also that in the result stated above, the fraction $p$ of positive instances within the pooled sample involved in the bipartite ranking step (Step 1) is the same as that within the pooled sample used in the rank test (Step 2), the same rank statistic $\hat{W}^{\phi}$ being used for both steps in the present analysis. However, if the choice of $N''$ can be guided by the value of the size $\alpha$ considered in view of the discrete pivotal conditional distribution of the test statistic used, the bound established in Theorem 3.5 does not permit tuning in practice the amount of the $N$ training observations that should be used for each Step 1 and Step 2.

As will be discussed in section 4, simply dividing the samples into two parts, with a larger part of the pooled dataset kept for the learning stage (e.g. $N' \sim 4N/5$) may permit to obtaining satisfactory results. The bound (37) explicitly shows that the rank-based testing step is independent on the underlying dimension $d$, importantly differing from existing methods that do suffer from increasing choices of the dimension $d$, see Section 2. The second term, related to the bipartite ranking procedure, is however dependent on the dimension through the constants induced by Assumption 3. Specifically, only constants $C_1$ and $C_2$ depend on $d$ when
assuming the class $S_0$ to be of VC-type. We finally point out that similar results to Theorem 3.5 can be established for the ranking-based test corresponding to the critical region (31) in the ROC space when the latter is based on the sup-norm, exactly in the same way except that bounds (in sup norm) established in Clémençon and Vayatis (2010) or in Clémençon and Vayatis (2009a) must be used instead of those in Clémençon, Limnios and Vayatis (2021).

**Remark 9. (Relation to minimax testing)** The proposed definition for the alternative assumption $H_1$ in the sense of the $W_\phi$-criterion is novel compared to the nonparametric minimax testing literature. The latter formulates $H_1$ by measuring the dissimilarity of the underlying distributions using a distance w.r.t. a particular metric, see Lam-Weil, Carpentier and Sriperumbudur (2022) for local minimax separation rate defined by $L_1$-norm for discrete distributions, Carpentier et al. (2018) for $L_2$-norm in sparse linear regression for instance. Classical two-sample testing literature based on rank statistics formulates lower-bounds for specific measures of efficiencies and choices of underlying distributions. For instance, Deb and Sen (2021) studies lower bounds on the Asymptotic Relative Efficiency (ARE) for multivariate Gaussian with independent coordinates (see Theorem 3.4 therein), and for symmetric elliptical distributions (see Theorem 3.4 therein). In the univariate setting, asymptotic results on the ARE and Pitman efficiency are proved with particular choices of distributions for the location problem, refer to Chernoff and Savage (1958) and to Hodges and Lehmann (1956) in particular. As shown by these results and discussed in Deb and Sen (2021), it is not possible to derive lower bounds for the ARE, both in the univariate (Hodges and Lehmann (1956)) and multivariate settings (Deb and Sen (2021), Proposition 3.5 therein applied to the contamination model with comparison to the Hotelling test). In fact, Definition 3.3 offers a similar minimax interpretation by considering the $W_\phi$-criterion in the ROC space as measure of departure for the null. It is particularly highlighted when choosing $\phi(u) = u$, see Eq. (16) and Remark 5 in Clémençon, Limnios and Vayatis (2021), that is explicitely showed by Example 1 discussed below Definition 3.4. Thus, one could formulate minimax properties of the type-II error of the proposed statistic, see the dedicated discussion in Limnios (2022), Chap. 6.4.2 therein.

**Remark 10. (Limiting the bias model effect)** As quantified by the parameter $\delta$ in the analysis above, the power of the learning-based testing procedure proposed is impacted by the choice of the class $S_0$. In order to limit the effect of model bias (i.e. of the choice of class $S_0$), classic (resampling-based) model selection techniques in bipartite ranking can be used in Step 1, see section 4 in Clémençon, Depecker and Vayatis (2011) for instance. However, as empirically illustrated in the experimental section 4, even if the model bias is significant, the learning-based approach promoted in this article outperforms its competitors in terms of power, especially in high-dimensional situations. From a practical perspective, we would recommend using the tree-based ranking algorithm RANKING FOREST in Step 1 (see Clémençon, Depecker and Vayatis (2013a)), whose hyperparameters can be easily tuned, that generally yields the best results in Section 4.

# 4. Illustrative numerical experiments

In this section, we illustrate the methodology previously described, and analyzed in the case of $W_\phi$-ranking performance optimization, by displaying the results of various numerical experiments, based on synthetic datasets. Beyond the empirical evidence of its capacity to detect (small) departures from the homogeneity assumption successfully, compared to certain
popular two-sample tests in the multivariate setup, they aim at showing the impact of the various ingredients involved, the score generating function $\phi$ chosen and the bipartite ranking algorithm used especially. We first detail the variants of the ranking-based approach to the two-sample problem examined in the experiments. Next the two-sample problems considered are specified, and the empirical results are finally summarized and discussed. Additional experimental results are postponed to the Appendix section C.2. All the Python codes used to carry out these experiments are available online<fnref n="3" /> for reproducibility purpose.

**4.1. Practical implementations - hyperparameters**

The two-stage testing procedure introduced in subsection 3.1 can be implemented in various ways, depending on the bipartite ranking algorithm $\mathcal{A}$ used and the score-generating function chosen.

**Bipartite ranking algorithms (step 1)** As recalled in subsection 2.3 (see also Remark 5), most techniques documented in the machine learning literature formulate bipartite ranking as a pairwise classification problem, see Clémençon, Lugosi and Vayatis (2008) and Menon and Williamson (2016). Hence, popular classification algorithms (e.g. SVM, neural networks, boosting) can be applied to pairs $(Z, Z')$ to which label $-1$ is assigned when $Z$ and $Z'$ are drawn from the same multivariate distribution ($H$ or $G$) and label $+1$ otherwise. In particular, for the present experiments we implemented the linear version of RankSVM with $L_2$ loss (rSVM2, see Joachims (2002)), RankNN (rNN, see Burges et al. (2005)) and RankBoost (rBoost, see Freund et al. (2003)). Beyond the pairwise approach, recursive partitioning methods based on oriented binary trees have been proposed to optimize an adaptively discretized version of the ROC curve (with remarkable generalization guarantees in sup norm), see Clémençon and Vayatis (2009a) and Clémençon, Depecker and Vayatis (2011) (see also Clémençon and Vayatis (2010) and Clémençon and Vayatis (2009c)). Here we have implemented the stabler ensemble learning version referred to as Ranking Forest (rForest, see Clémençon, Depecker and Vayatis (2013a)), whose practical performance is assessed in Clémençon, Depecker and Vayatis (2013b) in particular. In addition, as explained in Clémençon, Limnios and Vayatis (2021) (see section 4 therein), a regularized version of the empirical $W_\phi$-ranking performance criterion can be obtained by a kernel smoothing procedure and maximized through a (stochastic) gradient ascent algorithm, akin to that theoretically analyzed in subsection 3.2. We further refer to section 5 therein, for the dedicated empirical study of the ability of algorithms to approximate the true scoring function $s^*$, by the optimal $\hat{s}$ maximizing versions of $W_\phi$, to recover ROC$^*$.

**Univariate two-sample rank tests (step 2)** Let $\hat{s}$ be the outcome of the first step. To perform Step 2, we considered the test statistic $\hat{W}_{n,m}^\phi(\hat{s})$ with the following score generating functions: $\phi_{MWW}(u) = u$, which corresponds to the classic Mann-Whitney-Wilcoxon statistic (MWW, Wilcoxon (1945)), and $\phi_{RTB}(u) = u \mathbb{I}\{u \ge u_0\}$ for $u_0 \in (0, 1)$ (RTB, Clémençon and Vayatis (2007)) focusing on top ranks and the behavior of the ROC curve near the origin.

For a given level $\alpha$, the type-I error is controlled by taking as critical threshold the quantile $q_{n'',m''}^\phi(\alpha)$, while the type-II error is estimated by Monte Carlo for the pairs of instrumental

<footnote n="3">³https://github.com/MyrtoLimnios/twosampleranktest</footnote>
distributions described below. As discussed in subsection 2.2, the rank statistics are pivotal and the quantiles $q_{n'',m''}^{\phi}(\alpha)$ are tabulated for classic choices of $\phi$ (easily accessible by means of the SciPy open-access library available in Python for instance) or can be easily calculated, even for very large values of $n''$, $m''$, using modern computing frameworks.

We first consider both probabilities when the *dissimilarity/discrepancy* parameter $\epsilon > 0$ varies with fixed design, then we fix the dissimilarity parameter and let the dimensionality $d$ increase.

**Two-split trick** Since the most challenging step is undoubtedly the first one, consisting in learning to rank observations nearly in the same order as that induced by the likelihood ratio $\Psi(z)$, a much larger fraction of the data is allocated to bipartite ranking. Namely, we take $N'/N = 4/5$ (so that $N'' = N/5$).

**Benchmark tests** The performance of the ranking-based methodology proposed is compared to that of four state-of-the-art multivariate and nonparametric two-sample tests: namely the unbiased (quadratic) Maximum Mean Discrepancy (MMD) test with Gaussian kernels (see Gretton et al. (2007, 2012)), the graph-based Wald-Wolfowitz runs test (FR) generalized to the multivariate setting as proposed in Friedman and Rafsky (1979), the metric-based Energy test (Energy) (see Székely and Rizzo (2013)) and the depth-based procedure (Tukey) documented in Liu and Singh (1993) to extend univariate two-sample rank tests with the Tukey statistical depth (Tukey (1975)), see subsection 2.2 for further details. Notice that both MMD and Energy are not exactly distribution-free tests under $H_0$, therefore they are calibrated by means of a permutation technique. We also point out that Tukey is implemented by using the same tabulations of rank statistics as those used to apply the ranking-based methodology.

**Evaluation criteria** The frequencies of type-I error ($\epsilon = 0$) and of type-II error ($\epsilon > 0$) of each testing procedure at all levels $\alpha \in (0, 1)$ have been computed over $B \ge 1$ Monte-Carlo replications for each experiment and the distributions of $p$-values obtained are also reported. The impact of an increase of the dimension $d$ on the power, for fixed $\epsilon$ and sample size $N$, is investigated as well.

**Experimental parameters** For all the experiments, the pooled sample is of size $N = 2000$ and balanced, i.e., $n = m = N/2$ (corresponding to $p = 1/2$). The number of Monte Carlo replications is $B = 100$. The parameter $u_0$ for the RTB score-generating function varies in the set $\{0.7, 0.8, 0.9\}$. For the benchmark tests, the null distribution is estimated over $B_{perm} = 1000$ permutations and the related hyperparameter is optimized over the range $\{10^{-3}, 10^{-2}, 0.10, 1, 5, 10, 15, 20, 25, 30, 10^2, 10^3\}$. The figures also depict the pointwise confidence intervals at level 95%.

## 4.2. Synthetic datasets

We illustrate the performance of the family of rank-based tests proposed through the following two-sample problems. Various location and scale Gaussian models are considered in order to compare it to that of the (optimal) likelihood ratio statistic, see Example 1. We also test the homogeneity samples drawn from mixture models and from less classic (heavy-tailed) statistical distributions. See the Appendix section C.2 for additional information on
the parameters chosen for the distributions and the plots of the true ROC curves (including ROC"). For $d \ge 1$, denote by $0_d \in \mathbb{R}^d$ the null vector, by $\mathbf{1}_d \in \mathbb{R}^d$ the unit vector (i.e. with all
coordinates equal to 1), by $\mathbb{I}_d \in \mathbb{R}^{d \times d}$ the identity matrix and by $S_d^+(\mathbb{R})$ the cone of positive
semi-definite $d \times d$ matrices with real entries.

**Location Gaussian models** The two samples $\mathbf{X} \sim \mathcal{N}_d(\mu_X, \Sigma)$ and $\mathbf{Y} \sim \mathcal{N}_d(\mu_Y, \Sigma)$ are drawn independently, with $\varepsilon \in \{0.0, 0.02, 0.05, 0.08, 0.1\}, \Sigma \in S_d^+(\mathbb{R})$, as follows.

(L1) We set $\mu_Y = 0_d$ and $\mu_X = (\varepsilon/\sqrt{d}) \times \mathbf{1}_d$. Two models for the covariance matrix are considered: the first coordinate is negatively correlated with all the others and the other $d-1$ coordinates are: mutually independent (L1-); positively correlated (L1+).

For such models, the class $\mathcal{S}^*$ can be explicit, as detailed in Example 1.

**Scale Gaussian models** The two samples $\mathbf{X} \sim \mathcal{N}_d(0_d, \Sigma_X)$ and $\mathbf{Y} \sim \mathcal{N}_d(0_d, \Sigma_Y)$ are drawn independently with $\Sigma_X, \Sigma_Y$ in $S_d^+(\mathbb{R})$, defined as follows.

(S1) Decreasing correlation. Set $\Sigma_{X,i,j} = \alpha^{|i-j|}$ and $\Sigma_{Y,i,j} = \beta^{|i-j|}$ for $i, j \le d$, with $d=20$, $\beta = 0.2$ and $\alpha = \beta + \varepsilon, \varepsilon \in \{0, 0.1, 0.2, 0.3\}$.

(S2) Equi-correlated samples. Set $\Sigma_X = (1-\alpha)\mathbb{I}_d + \alpha\mathbf{1}_d\mathbf{1}_d^T$ and $\Sigma_Y = (1-\beta)\mathbb{I}_d + \beta\mathbf{1}_d\mathbf{1}_d^T$, with $d \in \{20, 30, 60, 100\}$, $\beta = 0.3$ and $\alpha = \beta + \varepsilon, \varepsilon \in \{0, 0.05, 0.1, 0.15\}$.

As for the Gaussian location model, an optimal scoring function can easily be explicit in this
case, the quadratic function $s_{\theta^*}(\cdot) = \langle \cdot, \theta^* \cdot \rangle$, where $\theta^* = \Sigma_X^{-1} - \Sigma_Y^{-1}$ for instance.

**Non Gaussian models** We also consider generative models with heavier tails to build other
two-sample problems.

(T1) *Cauchy distribution.* We generate $\mathbf{X} = (X_1, X_2, X_3)$ and $\mathbf{Y} = (Y_1, Y_2, Y_3)$ independently, such that $X_1, X_2 \stackrel{i.i.d}{\sim}$ Cauchy($\varepsilon, 1$) and $X_3, Y_1, Y_2, Y_3 \stackrel{i.i.d}{\sim}$ Cauchy(0, 1) with $\varepsilon \in \{0, 0.05, 0.15, 0.20\}, d=3$.

(T2) Lognormal location. We generate $\mathbf{X}$ and $\mathbf{Y}$, so that $\log(\mathbf{X})$ and $\log(\mathbf{Y})$ are drawn from the (L1) model with $\varepsilon \in \{0, 0.1, 0.2, 0.3\}, d = 4$.

(T3) Lognormal scale. We generate $\mathbf{X}$ and $\mathbf{Y}$, so that $\log \mathbf{X}$ and $\log \mathbf{Y}$ are drawn from the (S1) model with $\varepsilon \in \{0, 0.1, 0.2, 0.3\}, d = 20$.

The model (L1) is used in Clémençon, Limnios and Vayatis (2021) in the bipartite ranking
context, the models (S1), (S2) and (T1) are considered in Deb and Sen (2021), whereas the
models (T2) and (T3) are widely used in the univariate setting. Additional numerical experi-
ments, corroborating the empirical results described in this section, can be found in Appendix
section C.2.

**4.3. Results and discussion**

We now discuss the numerical results obtained by the approach for the two-sample problems
we have proposed and previously described. We compare its performance, depending on the
bipartite ranking algorithm chosen for Step 1 and the function $\phi$, to some state-of-the-art
(SoA) nonparametric tests on the basis of the type-I/II error counts. Focus is naturally on
<caption>Table 1</caption>

Rejection rates of the homogeneity assumption under $H_0$ and $H_1$ at significance level $\alpha = 0.05$ for the location models (L1-) and (L1+) and for the scale models (S1) and (S2), with $n = m = 1000$, $\pm$ their standard deviation at 95%. For ranking methods, only the results associated to MWW test are presented. Methods in bold minimize among all algorithms the rejection rate under $H_0$ while maximizing under $H_1$.

<table><thead><tr><th rowspan="2">Rejection rate of the null<br>Method</th><th colspan="4">Model (L1-), d = 6</th><th colspan="4">Model (L1+), d = 6</th></tr><tr><th>Under H<sub>0</sub><br>ε = 0.0</th><th>ε = 0.02</th><th>Under H<sub>1</sub><br>ε = 0.05</th><th>ε = 0.08</th><th>Under H<sub>0</sub><br>ε = 0.0</th><th>ε = 0.05</th><th>Under H<sub>1</sub><br>ε = 0.08</th><th>ε = 0.1</th></tr></thead><tbody><tr><td>rSVM2</td><td>0.11 ± 0.031</td><td><strong>0.22 ± 0.040</strong></td><td>0.60 ± 0.048</td><td>0.96 ± 0.019</td><td>0.10 ± 0.030</td><td><strong>0.23 ± 0.041</strong></td><td>0.37 ± 0.047</td><td><strong>0.66 ± 0.046</strong></td></tr><tr><td>rForest</td><td><strong>0.05 ± 0.021</strong></td><td>0.15 ± 0.035</td><td><strong>0.71 ± 0.044</strong></td><td><strong>0.98 ± 0.014</strong></td><td><strong>0.05 ± 0.021</strong></td><td>0.19 ± 0.038</td><td><strong>0.48 ± 0.049</strong></td><td>0.62 ± 0.047</td></tr><tr><td>rBoost</td><td>0.09 ± 0.028</td><td>0.07 ± 0.025</td><td>0.08 ± 0.026</td><td>0.11 ± 0.031</td><td>0.06 ± 0.023</td><td>0.10 ± 0.029</td><td>0.07 ± 0.025</td><td>0.06 ± 0.023</td></tr><tr><td>MMD</td><td><strong>0.05 ± 0.021</strong></td><td>0.04 ± 0.019</td><td>0.03 ± 0.017</td><td>0.04 ± 0.019</td><td>0.08 ± 0.026</td><td>0.06 ± 0.023</td><td>0.06 ± 0.023</td><td>0.06 ± 0.023</td></tr><tr><td>Energy</td><td>0.06 ± 0.023</td><td>0.05 ± 0.021</td><td>0.03 ± 0.017</td><td>0.06 ± 0.023</td><td>0.06 ± 0.023</td><td>0.05 ± 0.021</td><td>0.04 ± 0.019</td><td>0.02 ± 0.014</td></tr><tr><td>FR</td><td>0.09 ± 0.028</td><td>0.09 ± 0.028</td><td>0.09 ± 0.028</td><td>0.06 ± 0.023</td><td>0.08 ± 0.026</td><td>0.10 ± 0.029</td><td>0.09 ± 0.028</td><td>0.07 ± 0.025</td></tr><tr><td>Tukey</td><td>0.06 ± 0.023</td><td>0.05 ± 0.021</td><td>0.08 ± 0.026</td><td>0.07 ± 0.025</td><td><strong>0.05 ± 0.021</strong></td><td>0.08 ± 0.026</td><td>0.07 ± 0.025</td><td>0.09 ± 0.028</td></tr></tbody><thead><tr><th rowspan="2">Rejection rate of the null<br>Method</th><th colspan="4">Model (S1), d = 20</th><th colspan="4">Model (S2), d = 30</th></tr><tr><th>Under H<sub>0</sub><br>ε = 0.0</th><th>ε = 0.1</th><th>Under H<sub>1</sub><br>ε = 0.2</th><th>ε = 0.3</th><th>Under H<sub>0</sub><br>ε = 0.0</th><th>ε = 0.05</th><th>Under H<sub>1</sub><br>ε = 0.1</th><th>ε = 0.15</th></tr></thead><tbody><tr><td>rSVM2</td><td>0.06 ± 0.023</td><td>0.03 ± 0.017</td><td>0.06 ± 0.023</td><td>0.09 ± 0.028</td><td>0.04 ± 0.019</td><td>0.04 ± 0.019</td><td>0.04 ± 0.019</td><td>0.04 ± 0.019</td></tr><tr><td>rForest</td><td>0.05 ± 0.021</td><td><strong>0.30 ± 0.045</strong></td><td>0.95 ± 0.021</td><td><strong>1.00</strong></td><td><strong>0.03 ± 0.017</strong></td><td>0.19 ± 0.039</td><td>0.72 ± 0.044</td><td><strong>1.00</strong></td></tr><tr><td>rBoost</td><td>0.04 ± 0.021</td><td>0.02 ± 0.014</td><td>0.08 ± 0.026</td><td>0.10 ± 0.029</td><td>0.09 ± 0.028</td><td>0.04 ± 0.019</td><td>0.05 ± 0.021</td><td>0.03 ± 0.017</td></tr><tr><td>rNN</td><td><strong>0.03 ± 0.017</strong></td><td>0.27 ± 0.043</td><td><strong>0.98 ± 0.014</strong></td><td><strong>1.00</strong></td><td>0.04 ± 0.019</td><td><strong>0.20 ± 0.039</strong></td><td>0.54 ± 0.049</td><td><strong>1.00</strong></td></tr><tr><td>MMD</td><td>0.04 ± 0.019</td><td>0.08 ± 0.026</td><td>0.93 ± 0.025</td><td><strong>1.00</strong></td><td>0.04 ± 0.019</td><td>0.13 ± 0.033</td><td><strong>0.79 ± 0.040</strong></td><td><strong>1.00</strong></td></tr><tr><td>Energy</td><td>0.06 ± 0.023</td><td>0.07 ± 0.025</td><td>0.39 ± 0.048</td><td><strong>1.00</strong></td><td>0.04 ± 0.019</td><td>0.11 ± 0.030</td><td>0.18 ± 0.037</td><td>0.93 ± 0.025</td></tr><tr><td>FR</td><td>0.10 ± 0.029</td><td>0.19 ± 0.038</td><td>0.82 ± 0.037</td><td><strong>1.00</strong></td><td>0.06 ± 0.023</td><td>0.06 ± 0.023</td><td>0.28 ± 0.044</td><td>0.78 ± 0.040</td></tr><tr><td>Tukey</td><td>0.08 ± 0.026</td><td>0.00</td><td>0.02 ± 0.014</td><td>0.04 ± 0.019</td><td>0.04 ± 0.019</td><td>0.04 ± 0.019</td><td>0.04 ± 0.019</td><td>0.05 ± 0.021</td></tr></tbody></table>

the ability of each method to reject $H_0$ for small departures from it (the type of departure considered varying as well across the experiments), i.e., when $\epsilon \to 0$, while controlling the type-I error. Precisely, we discuss for all methods: 1) their ability to control the type-I error for the range of levels $\alpha \in (0, 1)$ (see the graphs), 2) the distribution of their $p$-values depending on $\epsilon$ (see the boxplots), and 3) their ability to reject the null hypothesis for the range of levels $\alpha \in (0, 1)$ as the dimension of the feature space $d$ increases (see the graphs).

The various tests are tuned to control the type-I error while maximizing the power at fixed level $\alpha = 0.05$. They are compared for the most difficult two-sample problems considered, i.e., the smallest value of $\epsilon$ (see the tables). Here, the results are displayed for the score generating function $\phi(u) = u$, which corresponds to the MWW test statistic. In the Appendix section, some results are also given for other score-generating functions in order to show how the choice of $\phi$ may impact the performance of the two-sample test (Step 2), in particular that of the RTB version. In addition to the figures and tables, we successively review the results obtained for each of the two-sample problems considered.

We first analyze the results for the Gaussian location models. For such two-sample problems, as evidenced by Table 1, ranking-based tests clearly manifest a greater capacity to reject the null hypothesis for small values of $\epsilon$ than SoA tests, while achieving a better control of the type-I error, see Figures 4 (a) and 5 (a), in particular rForest. The distributions of the $p$-values exhibit lower variance and higher power for very small $\epsilon = 0.02$ (L1-) for ranking-based methods compared to SoA tests, see Fig. 4 (b-d), as well as Fig. 5 (b-d) for (L1+), gathered in Table 1. Notice finally, that ranking-based algorithm show higher performance when the covariance matrix is sparse, e.g., model (L1-) with $d = 6$, compared to (L1+) with $d = 6$.

For the Gaussian scale models, we investigate the performance of the methods in higher dimensions. The results clearly reveal that all ranking-based tests control better the type-
<figure><img src="image_4.png" /><figcaption>Fig 4. Rejection rates of the homogeneity assumption under $H_0$ against significance level $\alpha \in (0, 1)$ in Fig. (a). Boxplots of p-values under $H_1$ in Fig. (b-d) for (L1-) with $\phi(u) = u$, $n = m = 1000$, $d = 6$, $B = 100$. The ranking algorithms are rSVM2, rForest, rBoost compared to MMD, Energy, FR and Tukey.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig 5. Rejection rates of the homogeneity assumption under $H_0$ against significance level $\alpha \in (0, 1)$ in Fig. (a). Boxplots of p-values under $H_1$ in Fig. (b-d) for (L1+) with $\phi(u) = u$, $n = m = 1000$, $d = 6$, $B = 100$. The ranking algorithms are rSVM2, rForest, rBoost compared to MMD, Energy, FR and Tukey.</figcaption></figure>

I error than SoA procedures, see Fig. 6 (a). In lower dimensions, the power of SoA tests is competitive, see Table 2, except for Tukey that exhibits a high p-value variance and a low power. In particular, both rForest and rNN have the highest rejection rates for the smaller $\epsilon$ (see columns (S2) $\epsilon = 0.1$ and (S3) $\epsilon = 0.05$) compared to the very low rates of SoA methods. Additionally, when the dimension $d$ of the feature space increases, the empirical power of ranking-based test rNN is always greater than that of any other method, whatever the test level $\alpha \in (0, 1)$, see Fig. 11 for (S2) and $d > 30$. However, results show that SoA methods are sensitive to the sparsity structure of the underlying distributions. In spite of higher dimension $d$ compared to the location models, MMD, FR and Energy have higher rejection rates for dense covariance matrices such as (S1). We refer to Limnios and
<figure><img src="image_6.png" /><figcaption>Fig 6. Rejection rates of the homogeneity assumption under $H_0$ against significance level $\alpha \in (0, 1)$ in Fig. (a). Boxplots of p-values under $H_1$ in Fig. (b-d) for (S1) with $\phi(u) = u$, $n = m = 1000$, $d = 20$, $B = 100$. The ranking algorithms are rSVM2, rForest, rBoost, rNN compared to MMD, Energy, FR and Tukey.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig 7. Rejection rates of the homogeneity assumption under $H_0$ against significance level $\alpha \in (0, 1)$ in Fig. (a). Boxplots of p-values under $H_1$ in Fig. (b-d) for (S2) with $\phi(u) = u$, $n = m = 1000$, $d = 30$, $B = 100$. The ranking algorithms are rSVM2, rForest, rBoost, rNN compared to MMD, Energy, FR and Tukey.</figcaption></figure>

Clémençon (2024) for further discussion.

For the third class of models, we analyze the distributions of the p-values in Fig. 8 for (T1), Fig. 9 for (T2) with (L1+), and Fig. 10 for (T3) with (S1). rBoost and MMD perform similarly for (T1), i.e., multivariate Cauchy distribution. rForest and MMD have comparable results for both (T2) and (T3). Overall the three models and under the alternative, rForest shows higher empirical rejection rates for the smallest $\epsilon$. Lastly, the rejection rate under the null is better controlled by ranking-based algorithms see Fig. 9 (a) and Fig. 10 (a).

Finally, we tested the impact of the increase of $d$ on the rejection rate under alternatives for all methods w.r.t. the range of levels $\alpha \in (0, 1)$ for the (S2) two-sample problem. Figure 11 clearly shows that for the smallest $\epsilon = 0.05$, the empirical power of the ranking-based method
<caption>Table 2</caption>

Rejection rates of the homogeneity assumption under $H_0$ and $H_1$ at significance level $\alpha = 0.05$ for models (T1-3), with $n = m = 1000$, $\pm$ their standard deviation at 95%. For ranking methods, only the results associated to MWW test are presented. Methods in bold minimize among all algorithms the rejection rate under $H_0$ while maximizing under $H_1$.

<table><thead><tr><th rowspan="3">Rejection rate of the null Method</th><th colspan="4">Model (T1), d = 3</th></tr><tr><th>Under H<sub>0</sub></th><th colspan="3">Under H<sub>1</sub></th></tr><tr><th>ε = 0.0</th><th>ε = 0.05</th><th>ε = 0.15</th><th>ε = 0.20</th></tr></thead><tbody><tr><td>rSVM2</td><td><strong>0.05 ± 0.021</strong></td><td>0.06 ± 0.023</td><td>0.06 ± 0.023</td><td>0.14 ± 0.034</td></tr><tr><td>rForest</td><td>0.06 ± 0.023</td><td><strong>0.08 ± 0.026</strong></td><td>0.14 ± 0.034</td><td>0.17 ± 0.037</td></tr><tr><td>rBoost</td><td><strong>0.05 ± 0.021</strong></td><td>0.06 ± 0.023</td><td>0.26 ± 0.043</td><td>0.41 ± 0.048</td></tr><tr><td><strong>MMD</strong></td><td>0.06 ± 0.023</td><td>0.04 ± 0.019</td><td><strong>0.35 ± 0.047</strong></td><td><strong>0.57 ± 0.048</strong></td></tr><tr><td>Energy</td><td><strong>0.05 ± 0.021</strong></td><td>0.05 ± 0.021</td><td>0.09 ± 0.028</td><td>0.15 ± 0.035</td></tr><tr><td>FR</td><td><strong>0.05 ± 0.021</strong></td><td>0.05 ± 0.021</td><td>0.08 ± 0.026</td><td>0.06 ± 0.023</td></tr><tr><td>Tukey</td><td>0.06 ± 0.023</td><td>0.07 ± 0.025</td><td>0.04 ± 0.019</td><td>0.07 ± 0.025</td></tr></tbody><thead><tr><th rowspan="3">Rejection rate of the null Method</th><th colspan="4">Model (T2), d = 4</th></tr><tr><th>Under H<sub>0</sub></th><th colspan="3">Under H<sub>1</sub></th></tr><tr><th>ε = 0.0</th><th>ε = 0.10</th><th>ε = 0.20</th><th>ε = 0.30</th></tr></thead><tbody><tr><td>rSVM2</td><td>0.03 ± 0.017</td><td>0.06 ± 0.023</td><td>0.04 ± 0.019</td><td>0.03 ± 0.017</td></tr><tr><td>rForest</td><td>0.05 ± 0.021</td><td>0.12 ± 0.032</td><td>0.43 ± 0.048</td><td>0.80 ± 0.039</td></tr><tr><td>rBoost</td><td>0.06 ± 0.023</td><td>0.09 ± 0.028</td><td>0.18 ± 0.037</td><td>0.34 ± 0.046</td></tr><tr><td><strong>MMD</strong></td><td>0.05 ± 0.021</td><td><strong>0.18 ± 0.038</strong></td><td><strong>0.56 ± 0.048</strong></td><td><strong>0.92 ± 0.026</strong></td></tr><tr><td>Energy</td><td>0.06 ± 0.023</td><td>0.08 ± 0.026</td><td>0.16 ± 0.038</td><td>0.32 ± 0.045</td></tr><tr><td>FR</td><td><strong>0.01 ± 0.010</strong></td><td>0.02 ± 0.014</td><td>0.08 ± 0.026</td><td>0.28 ± 0.044</td></tr><tr><td>Tukey</td><td>0.04 ± 0.019</td><td>0.12 ± 0.032</td><td>0.27 ± 0.043</td><td>0.46 ± 0.049</td></tr></tbody><thead><tr><th rowspan="3">Rejection rate of the null Method</th><th colspan="4">Model (T3), d = 20</th></tr><tr><th>Under H<sub>0</sub></th><th colspan="3">Under H<sub>1</sub></th></tr><tr><th>ε = 0.0</th><th>ε = 0.1</th><th>ε = 0.2</th><th>ε = 0.3</th></tr></thead><tbody><tr><td>rSVM2</td><td>0.05 ± 0.021</td><td>0.06 ± 0.023</td><td>0.05 ± 0.021</td><td>0.13 ± 0.033</td></tr><tr><td><strong>rForest</strong></td><td>0.04 ± 0.019</td><td><strong>0.23 ± 0.041</strong></td><td><strong>0.99 ± 0.010</strong></td><td><strong>1.00</strong></td></tr><tr><td>rBoost</td><td>0.05 ± 0.021</td><td>0.07 ± 0.025</td><td>0.11 ± 0.031</td><td>0.24 ± 0.042</td></tr><tr><td><strong>MMD</strong></td><td>0.07 ± 0.025</td><td>0.15 ± 0.035</td><td>0.83 ± 0.037</td><td>1.00</td></tr><tr><td>Energy</td><td>0.06 ± 0.023</td><td>0.09 ± 0.028</td><td>0.40 ± 0.048</td><td>0.98 ± 0.014</td></tr><tr><td>FR</td><td>0.07 ± 0.025</td><td>0.09 ± 0.028</td><td>0.63 ± 0.047</td><td>1.00</td></tr><tr><td>Tukey</td><td><strong>0.03 ± 0.017</strong></td><td>0.05 ± 0.021</td><td>0.04 ± 0.019</td><td>0.00</td></tr></tbody></table>

rNN is always above the SoA procedures, with a notable difference for the highest dimension considered, namely $d = 100$. When $\epsilon$ increases, then rNN and MMD have similar performance. Table 3 gathers the estimated power for $\alpha = 0.05$ and illustrates the non-decreasing power of the ranking-based method with the dimension, where SoA methods overall have constant/decreasing power w.r.t. the dimension, see Energy $\epsilon = 0.1$.

To conclude, we empirically illustrated the competitiveness of ranking-based rank tests, especially for the location and pathological statistical problems, as it overall shows a clear control of the empirical type-I error for a large range of levels $\alpha$, while resulting to similar or higher rejection rates under alternatives for very small deviations $\epsilon \to 0$. From these experiments, the high performer ranking-based algorithm for Step 1 shows to be rForest, while the lower is rBoost. The comparative SoA tests overall control the type-I error for the majority of models but at the price of algorithmic complexity due to the high number of permutations. Additionally, their empirical power for very small $\epsilon$ is lower than this of the ranking-based methods. The depth-based test (Tukey) particularly under performs. Lastly,
<figure><img src="image_8.png" /><figcaption>Fig 8. Rejection rates of the homogeneity assumption under $\mathcal{H}_0$ against significance level $\alpha \in (0, 1)$ in Fig. (a). Boxplots of p-values under $\mathcal{H}_1$ in Fig. (b-d) for (T1) with $\phi(u) = u$, $n = m = 1000$, $d = 3$, $B = 100$. The ranking algorithms are rSVM2, rForest, rBoost compared to MMD, Energy, FR and Tukey.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig 9. Rejection rates of the homogeneity assumption under $\mathcal{H}_0$ against significance level $\alpha \in (0, 1)$ in Fig. (a). Boxplots of p-values under $\mathcal{H}_1$ in Fig. (b-d) for (T2) with (LI-) with $\phi(u) = u$, $n = m = 1000$, $d = 4$, $B = 100$. The ranking algorithms are rSVM2, rForest, rBoost compared to MMD, Energy, FR and Tukey.</figcaption></figure>

the distributions of the p-values show large variance for decreasing $\varepsilon$ and are especially larger for the SoA methods, independently of the underlying probabilistic model.

## 5. Conclusion

In this paper, we have developed a fully novel and flexible approach to the two-sample problem, ubiquitous in statistical and machine-learning applications. It relies on the observation that the null/homogeneity assumption $\mathcal{H}_0$ ceases to hold true as soon as the optimal ROC curve for the bipartite ranking problem related to the probability distributions H and G of the two samples observed does not coincide with the diagonal (and is then necessarily above
<figure><img src="image_10.png" /><figcaption>Fig 10. Rejection rates of the homogeneity assumption under $H_0$ against significance level $\alpha \in (0, 1)$ in Fig. (a). Boxplots of p-values under $H_1$ in Fig. (b-d) for (T3) with (S1) with $\phi(u) = u$, $n = m = 1000$, $d = 20$, $B = 100$. The ranking algorithms are rSVM2, rForest, rBoost compared to MMD, Energy, FR and Tukey.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig 11. Rejection rates of the homogeneity assumption under $H_1$ against significance level $\alpha \in (0, 1)$ for increasing dimension $d \in \{30, 60, 100\}$ of the feature space for (S2) with $\phi(u) = u$, $n = m = 1000$, $B = 100$ and $\varepsilon = 0.05$ (a-c), $\varepsilon = 0.1$ (e-g). The ranking algorithm is rNN compared to MMD, Energy, FR and Tukey.</figcaption></figure>

it at some points of the ROC space). As the curve $ROC_{H,G}^*$ is not directly observable, the methodology proposed involves a statistical learning step to solve the bipartite ranking problem based on a first fraction of the samples. It next consists in evaluating the departure from $H_0$ using an empirical version of (a scalar summary of) $ROC_{H,G}^*$. In practice, the second
<caption>Table 3</caption>

Rejection rates of the homogeneity assumption under $H_1$ at significance level $\alpha = 0.05$ for the scale model (S2), with $n = m = 1000$, $\pm$ their standard deviation at 95%. For $\tau$NN, only the results associated to MWW test are presented. Methods in bold maximize among all algorithms the rejection rate under $H_1$.

<table><thead><tr><th rowspan="2">Rejection rate of the null<br>Method</th><th colspan="3">Under H<sub>1</sub>, (ε = 0.05)</th><th colspan="3">Under H<sub>1</sub>, (ε = 0.1)</th></tr><tr><th>d = 30</th><th>d = 60</th><th>d = 100</th><th>d = 30</th><th>d = 60</th><th>d = 100</th></tr></thead><tbody><tr><td>rNN</td><td><b>0.20 ± 0.039</b></td><td><b>0.23 ± 0.041</b></td><td><b>0.42 ± 0.048</b></td><td>0.54 ± 0.049</td><td><b>0.84 ± 0.036</b></td><td>0.93 ± 0.025</td></tr><tr><td>MMD</td><td>0.13 ± 0.033</td><td>0.09 ± 0.028</td><td>0.18 ± 0.037</td><td><b>0.89 ± 0.031</b></td><td>0.76 ± 0.042</td><td><b>0.99 ± 0.010</b></td></tr><tr><td>Energy</td><td>0.04 ± 0.019</td><td>0.05 ± 0.021</td><td>0.11 ± 0.031</td><td>0.27 ± 0.43</td><td>0.22 ± 0.040</td><td>0.23 ± 0.041</td></tr><tr><td>FR</td><td>0.06 ± 0.023</td><td>0.11 ± 0.030</td><td>0.13 ± 0.033</td><td>0.39 ± 0.048</td><td>0.43 ± 0.048</td><td>0.38 ± 0.047</td></tr><tr><td>Tukey</td><td>0.08 ± 0.026</td><td>0.09 ± 0.028</td><td>0.11 ± 0.030</td><td>0.11 ± 0.030</td><td>0.07 ± 0.025</td><td>0.11 ± 0.030</td></tr></tbody></table>

step simply boils down to performing a classic two-sample rank test based on the ranks of the fraction of the samples not used in the learning stage, obtained by means of the ranking rule constructed in the first step.

A sound theoretical analysis has been carried out, showing that the two types of errors made by the method can be controlled in a nonasymptotic fashion. It is complemented by an experimental study providing strong empirical evidence of the merits of the ranking-based method compared to some popular alternatives: in general, it exhibits a performance that better withstands high-dimensionality, better resists as one gets closer and closer to the null hypothesis and adapts to a wide variety of ways of departing from homogeneity. However, several issues remain to be investigated so that the general approach introduced here can systematically perform well in practice.

One may state two open problems in particular. Understanding how to dispatch the data at disposal, for the statistical learning purpose and for performing the rank test, as well as developing a method to combine efficiently several ranking-based rank tests to provably improve performance (insofar as one does not know in advance which score generating function $\phi$ may capture best the possible type of departure from homogeneity exhibited), requires further research.

## Appendix A: Properties of the ROC curves

For the sake of completeness, basic properties of ROC curves are recalled. See e.g. Proposition 17 in Clémençon and Vayatis (2009a).

**Proposition A.1.** (Properties of the ROC curve) For any probability distributions H and G and any scoring function s : $Z \to \mathbb{R}$, the following properties hold:

1. **Limit values.** We have:

$$ROC(s, 0) = 0 \text{ and } ROC(s, 1) = 1.$$

2. **Invariance.** For any strictly increasing function $T : \mathbb{R} \to \mathbb{R}$, we have, for all $\alpha \in (0, 1)$: $ROC(T \circ s, \alpha) = ROC(s, \alpha)$.

3. **Concavity.** If the likelihood ratio $dG_s/dH_s$ is a monotone function then the ROC curve is concave. In particular, the optimal curve ROC* is concave.

4. **Linear parts.** If the likelihood ratio $dG_s/dH_s$ is constant on some interval in the range of the scoring function s then the ROC curve will present a linear part on the corre-
sponding domain. Furthermore, ROC* is linear on $[\alpha_1, \alpha_2]$ iff. $\Psi(z) = (dG/dH)(z)$ is constant on the subset

$$ \{z \in \mathbb{Z} / Q^*(\alpha_2) \leq \Psi(z) \leq Q^*(\alpha_1)\}, $$

where $Q^*(\alpha)$ means the quantile at level $1 - \alpha \in (0, 1)$ of the r.v. $\Psi(\mathbf{Y})$.

5. **Differentiability.** Assume that the distributions $H$ and $G$ are continuous. Then, the ROC curve of a scoring function $s$ is differentiable iff. the distributions $H_s$ and $G_s$ are continuous.

# Appendix B: Technical proofs

## B.1. Proof of formula 10

Given two real-valued, independent and continuous random variables $X \sim G$ and $Y \sim H$, the proof is straightforward by noticing that:

$$ \int_{-\infty}^{\infty} \{H(t) - G(t)\}dH(t) = \frac{1}{2} - \int_{-\infty}^{\infty} G(t)dH(t). $$

Hence, we have:

$$ \begin{aligned} \mathrm{AUC}_{H,G} &= \mathbb{P}\{Y \le X\} = \mathbb{E}[\mathbb{E}[\mathbb{I}\{Y \le X\}] | Y] = \mathbb{E}_{Y \sim H}[1 - G(Y)] \\ &= \frac{1}{2} + \int_{-\infty}^{\infty} \{H(t) - G(t)\}dH(t). \end{aligned} $$

## B.2. Proof of Proposition 2.2

The equivalence between assertions (i) and (ii) can be directly deduced from the following result, proved in Clémençon and Vayatis (2009a) (see Corollary 5 and Proposition 6's proof therein), recalled here for the sake of clarity.

**Lemma B.1.** (Clémençon and Vayatis (2009a)) It holds with probability one:

$$ \Psi(\mathbf{Z}) = \frac{dG_{\Psi}}{dH_{\Psi}}(\Psi(\mathbf{Z})), $$

$\mathbf{Z}$ denoting either $\mathbf{X}$ or $\mathbf{Y}$.

The equivalence of assertion (ii) with the other assertions is immediate using formula (15).

## B.3. An exponential inequality for two-sample linear rank statistics

The result stated below provides a tail bound of exponential type for a two-sample linear rank statistic (see Definition 2.1), when recentered by its asymptotic mean (9).

**Theorem B.2.** Let $p \in (0, 1)$ and $N \ge 1/p$. Set $n = \lfloor pN \rfloor$ and $m = \lceil (1 - p)N \rceil = N - n$. Let $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$ be two independent i.i.d. random samples, drawn
from univariate probability distributions G and H respectively. Suppose that Assumption 1 is fulfilled. Then, for all $t > 0$, we have:

$$ \mathbb{P} \left\{ \frac{1}{n} \widehat{W}_{n,m}^{\phi} - W_{\phi} > t \right\} \le 18 \exp \{-CNt^2\}, \quad (38) $$

$$ \text{where } C = 8^{-1} \min (p/\|\phi\|_{\infty}^2, (p/\|\phi'\|_{\infty}^2)^{-1}, ((1-p)/\|\phi'\|_{\infty}^2)^{-1}). $$

*Proof.* The proof is based on the linearization technique used to establish the concentration inequality for two-sample linear rank processes in Clémençon, Limnios and Vayatis (2021), see Theorem 5 therein. Consider $\phi$ fulfilling Assumption 1. Summing $\phi$'s Taylor expansions of order 2, evaluated at $N\widehat{F}_N(X_i)/(N+1)$ around $F(X_i)$ for $1 \le i \le n$ leads to an a.s. decomposition of the statistic $\widehat{W}_{n,m}^{\phi}$, see Eq. (B.3,4) in Clémençon, Limnios and Vayatis (2021). After the application of Hoeffding's decomposition to its first order term, which takes the form of a two-sample U-statistic, we write the decomposition as follows:

$$ \frac{1}{n} \widehat{W}_{n,m}^{\phi} - W_{\phi} = \widehat{W}_{\phi} - W_{\phi} + \frac{1}{n} \left( \widehat{V}_{n}^{X} - \mathbb{E} \left[ \widehat{V}_{n}^{X} \right] \right) + \frac{1}{n} \left( \widehat{V}_{m}^{Y} - \mathbb{E} \left[ \widehat{V}_{m}^{Y} \right] \right) + \frac{1}{n} \mathcal{R}_{n,m}, \quad (39) $$

where

$$ \widehat{W}_{\phi} = \frac{1}{n} \sum_{i=1}^{n} (\phi \circ F)(X_i), $$

$$ \widehat{V}_{n}^{X} = \frac{n-1}{N+1} \sum_{i=1}^{n} \int_{X_i}^{+\infty} (\phi' \circ F)(u) dG(u), $$

$$ \widehat{V}_{m}^{Y} = \frac{n}{N+1} \sum_{j=1}^{m} \int_{Y_j}^{+\infty} (\phi' \circ F)(u) dG(u). $$

The quantity $\mathcal{R}_{n,m}$ corresponds to the Taylor-Lagrange residual term plus additional statistics of order $O_p(N^{-1})$ inherited from Hoeffding's decomposition. The lemma below provides a tail bound for the latter.

**Lemma B.3.** Suppose that the assumptions of Theorem B.2 are satisfied. Then, for all $t > 0$ and $N \ge 2$, we have, if $Nt \ge 128 \|\phi'\|_{\infty}^2 / (p \|\phi''\|_{\infty})$:

$$ \mathbb{P} \{|\mathcal{R}_{n,m}| > t\} \le 12 \exp \left\{ -\frac{Nt}{12\kappa_p \|\phi''\|_{\infty}^2} \right\}, \quad (40) $$

and otherwise

$$ \mathbb{P} \{|\mathcal{R}_{n,m}| > t\} \le 12 \exp \left\{ -\frac{\alpha_p N^2 t^2}{512 \|\phi'\|_{\infty}^2} \right\}, \quad (41) $$

$$ \text{where } \alpha_p = \min(p/(1-p), 1), \quad \kappa_p = \max(p, 1-p). $$

For any $t > 0$, a straightforward application of the classic exponential tail inequality, proved in Hoeffding (1963), to each of the first three terms on the right hand side of (39) yields for any $t > 0$

$$ \mathbb{P} \left\{ |\widehat{W}_{\phi} - W_{\phi}| > t \right\} \le 2 \exp \left\{ - \frac{2pNt^2}{\|\phi\|_{\infty}^2} \right\}, $$
$$
\begin{align*}
\mathbb{P}\left\{\frac{1}{n}\left|\hat{V}_n^X - \mathbb{E}\left[\left|\hat{V}_n^X\right|\right]\right| > t\right\} &\le 2 \exp\left\{-\frac{2Nt^2}{p\|\phi'\|_{\infty}^2}\right\}, \\
\mathbb{P}\left\{\frac{1}{n}\left|\hat{V}_m^Y - \mathbb{E}\left[\left|\hat{V}_m^Y\right|\right]\right| > t\right\} &\le 2 \exp\left\{-\frac{2Nt^2}{(1-p)\|\phi'\|_{\infty}^2}\right\}.
\end{align*}
$$

Then, the application of the four tail bounds above with threshold $t/4$, combined with expansion (39) and the union bound, concludes the proof. □

B.4. *Proof of Lemma B.3*

Consider the nonsymmetric bounded kernels on $\mathbb{R}^2$:

$$
k(x, x') = \mathbb{I}\{x' \le x\} (\phi' \circ F)(x),
$$

$$
\ell(x, y) = \mathbb{I}\{y \le x\} (\phi' \circ F)(x) .
$$

The remainder process is decomposed as follows

$$
\mathcal{R}_{n,m} = \hat{R}_{n,m} + \frac{n(n-1)}{N+1} U_n(k) + \frac{nm}{N+1} U_{n,m}(\ell) + \hat{T}_{n,m},
$$

where

$$
\begin{align*}
\hat{R}_{n,m} = \left(\frac{n}{N+1} - p\right) \sum_{i=1}^{n} \left(\mathbb{E}[k(X_i, X) | X_i] - \mathbb{E}[k(X, X')]\right) \\
\qquad + \left(\frac{m}{N+1} - 1 + p\right) \sum_{i=1}^{n} \left(\mathbb{E}[\ell(X_i, Y) | X_i] - \mathbb{E}[\ell(X, Y)]\right)
\end{align*}
$$

$(X, X')$ and $Y$ denoting two independent r.v.'s with distributions $G \otimes G$ and $H$ respectively,
independent from the $X_i$'s, $U_n(k)$ is the one-sample degenerate $U$-statistic of order 2 based
on $\{X_1, \dots, X_n\}$ with kernel $k$, $U_{n,m}(\ell)$ is the two-sample degenerate $U$-statistic of degree
(1, 1) based on the two samples $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$ with kernel $\ell$, $\hat{T}_{n,m}$ is the
Taylor-Lagrange integral remainder of the expansion of $\phi$ at order 2, namely

$$
\hat{T}_{n,m} = \sum_{i=1}^{n} \int_{F(X_i)}^{N \hat{F}_N(X_i)/(N+1)} \left( \frac{N \hat{F}_N(X_i)}{N+1} - u \right) \phi''(u) du. \quad (42)
$$

We clearly have

$$
|\mathcal{R}_{n,m}| \le |\hat{R}_{n,m}| + p^2 N |U_n| + p(1-p)N |U_{n,m}| + |\hat{T}_{n,m}|.
$$

We now exhibit a tail bound for each of the terms involved in the sum on the right hand side of the bound above. First, observing that the mappings $x \mapsto \mathbb{E}[k(x, X)] = G(x)\phi' \circ F(x)$ and $x \mapsto \mathbb{E}[k(x, Y)] = H(x)\phi' \circ F(x)$ are bounded by $\|\phi'\|_{\infty}$, that $|n/(N+1) - p| \le 1/N$ and $|m/(N+1) - (1-p)| \le 1/N$, Hoeffding's inequality applied twice and combined with the union bound directly yields, for any $t > 0$,

$$
\mathbb{P}\left\{\frac{1}{n}\left|\hat{R}_{n,m}\right| > t/4\right\} \le 4 \exp \left\{-\frac{pN^3t^2}{32\|\phi'\|_{\infty}^2}\right\}. \quad (43)
$$

Tail bounds for the two degenerate *U*-processes are provided by Lemma 27 in Clémençon, Limnios and Vayatis (2021), recalled below for clarity.
**Lemma B.4.** Let P and Q be two probability distributions on measurable spaces $\mathcal{X}$ and $\mathcal{Y}$ respectively. Consider the degenerate two-sample U-statistic of degree (1, 1) with a bounded kernel $\theta: \mathcal{X} \times \mathcal{Y} \to \mathbb{R}$ based on the independent i.i.d. random samples $\{X_1, \dots, X_n\}$ and $\{Y_1, \dots, Y_m\}$, drawn from P and Q respectively. For all $t > 0$, we then have:

$$ \mathbb{P}\{U_{n,m}(\theta) \geq t\} \leq \exp\left\{-\frac{nmt^2}{32c_{\theta}^2}\right\}, $$

where $c_{\theta} = \sup_{(x,y) \in \mathcal{X} \times \mathcal{Y}} |\theta(x, y)| < +\infty$.

A direct application of Lemma B.4 to $U_{n,m}$ yields, for all $t > 0$,

$$ \mathbb{P}\left\{\frac{p(1-p)N}{n} |U_{n,m}| > \frac{t}{4}\right\} \leq 2 \exp\left\{-\frac{pN^2t^2}{512(1-p)\|\phi'\|_{\infty}^2}\right\}, \quad (44) $$

and similarly, by virtue of Lemma 3 in Nolan and Pollard (1987), we have

$$ \mathbb{P}\left\{\frac{p^2 N}{n} |U_n| > \frac{t}{4}\right\} \leq 2 \exp\left\{-\frac{N^2 t^2}{512 \|\phi'\|_{\infty}^2}\right\}. \quad (45) $$

Finally, from Eq. (42) we derive

$$ \frac{1}{n} |\hat{T}_{n,m}| \le \|\phi''\|_{\infty} \left( \sup_{t \in \mathbb{R}} \left( \hat{F}_N(t) - F(t) \right)^2 + \frac{1}{(N+1)^2} \right) \\ \le 3p^2 \|\phi''\|_{\infty} \sup_{t \in \mathbb{R}} \left( \hat{G}_n(t) - G(t) \right)^2 + 3(1-p)^2 \|\phi''\|_{\infty} \sup_{t \in \mathbb{R}} \left( \hat{H}_m(t) - H(t) \right)^2 + \frac{13 \|\phi''\|_{\infty}}{N^2}. $$

Considering the terms involved on the right hand side of the bound above, by applying the classic Dvoretzky-Kiefer-Wolfowitz inequality twice and noticing that the third term is negligible compared to the first two terms, we obtain:

$$ \mathbb{P}\left\{\frac{1}{n} |\hat{T}_{n,m}| > \frac{t}{4}\right\} \leq 4 \exp\left\{-\frac{Nt}{12\kappa_p \|\phi''\|_{\infty}}\right\}, \quad (46) $$

where $\kappa_p = \max(p, 1-p)$. Combining the union bound with equations (43), (44), (45) and (46) concludes the proof.

## B.5. Proof of Theorem 3.1

It suffices to observe that we almost-surely have

$$ \begin{aligned} \mathbb{P}_{H_0} \left\{ \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) > \int_0^1 \phi(u) du + q_{n'',m''}^{\phi}(\alpha) \mid \mathbf{X}_i, \mathbf{Y}_j, \forall i \le n', \forall j \le m' \right\} \\ = \mathbb{P}_{H_0} \left\{ \frac{1}{n''} \hat{W}_{n'',m''} > \int_0^1 \phi(u) du + q_{n'',m''}^{\phi}(\alpha) \right\} \le \alpha \end{aligned} $$

and take the expectation in the equation above w.r.t. $\{\mathbf{X}_1, \dots, \mathbf{X}_{n'}, \mathbf{Y}_1, \dots, \mathbf{Y}_{m'}\}$.
## B.6. Proof of Proposition 3.2

Notice that, by virtue of Theorem B.2, we have

$$ \mathbb{P}_{H_0} \left\{ \frac{1}{n} \hat{W}_{n,m}^{\phi} - \int_{0}^{1} \phi(u) du > \sqrt{\frac{\log(18/\alpha)}{CN}} \right\} \le \alpha, $$

which yields the desired bound.

## B.7. Proof of Theorem 3.5

Let $\alpha \in (0, 1)$ and $\hat{s} \in S_0$ be a solution of Step 1 based on the two samples $\mathcal{D}_{n',m'} = \{\mathbf{X}_1, \dots, \mathbf{X}_{n'}\} \cup \{\mathbf{Y}_1, \dots, \mathbf{Y}_{m'}\}$. First, observe that, for all $(H, G) \in \mathcal{H}_1(\varepsilon)$, we have:

$$ \begin{aligned} \mathbb{P}_{H,G} \{ \Phi_{\alpha}^{\phi}(\mathcal{D}_{n'',m''}(\hat{s})) = 0 \} &= \mathbb{P}_{H,G} \left\{ \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) - \int_{0}^{1} \phi(u) du \le q_{n'',m''}^{\phi}(\alpha) \right\} \\ &\le \mathbb{P}_{H,G} \left\{ 2 \sup_{s \in S_0} \left| \frac{1}{n'} \hat{W}_{n',m'}(s) - W_{\phi}(s) \right| + \left| \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) - W_{\phi}(\hat{s}) \right| \ge \varepsilon - \delta - \sqrt{\frac{\log(18/\alpha)}{CN''}} \right\} \\ &\le \mathbb{P}_{H,G} \left\{ 2 \sup_{s \in S_0} \left| \frac{1}{n'} \hat{W}_{n',m'}(s) - W_{\phi}(s) \right| \ge \frac{\varepsilon - \delta}{4} \right\} + \mathbb{P}_{H,G} \left\{ \left| \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) - W_{\phi}(\hat{s}) \right| \ge \frac{\varepsilon - \delta}{4} \right\}, \end{aligned} \quad (47) $$

where $C = 8^{-1} \min(p/\|\phi\|_{\infty}^2, (p/\|\phi'\|_{\infty}^2)^{-1}, ((1-p)/\|\phi'\|_{\infty}^2)^{-1})$, using decomposition (34) above and the bound

$$ W_{\phi}^{*} - W_{\phi}(\hat{s}) \le 2 \sup_{s \in S_0} \left| \frac{1}{n'} \hat{W}_{n',m'}(s) - W_{\phi}(s) \right| + \delta, $$

as well as Proposition 3.2, Definition 3.3, the union bound and $N'' \ge 4 \log(18/\alpha)/(C(\varepsilon - \delta)^2)$. By applying Theorem B.2 to the univariate samples:

$$ \{\hat{s}(\mathbf{X}_{1+n'}), \dots, \hat{s}(\mathbf{X}_n)\} \text{ and } \{\hat{s}(\mathbf{Y}_{1+m'}), \dots, \hat{s}(\mathbf{Y}_m)\}, $$

when conditioning on $\mathcal{D}_{n',m'} = \{\mathbf{X}_1, \dots, \mathbf{X}_{n'}\} \cup \{\mathbf{Y}_1, \dots, \mathbf{Y}_{m'}\}$ (and consequently the scoring function $\hat{s}$ produced at Step 1), we obtain that we almost-surely have

$$ \mathbb{P}_{H,G} \left\{ \left| \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) - W_{\phi}(\hat{s}) \right| \ge \frac{\varepsilon - \delta}{4} \mid \mathcal{D}_{n',m'} \right\} \le 18 \exp \left\{ -\frac{CN'' (\varepsilon - \delta)^2}{16} \right\} $$

as soon as $N'' \ge 4 \log(18/\alpha)/(C(\varepsilon - \delta)^2)$. By taking the expectation, the bound above yields

$$ \mathbb{P}_{H,G} \left\{ \left| \frac{1}{n''} \hat{W}_{n'',m''}^{\phi}(\hat{s}) - W_{\phi}(\hat{s}) \right| \ge \frac{\varepsilon - \delta}{4} \right\} \le 18 \exp \left\{ -\frac{CN'' (\varepsilon - \delta)^2}{16} \right\}. \quad (48) $$

We also recall the following result providing a tail bound for the maximal deviations between the $W_{\phi}$-ranking performance and its empirical version, where the detailed constants are in the corresponding proof.
**Theorem B.5.** (Clémençon, Limnios and Vayatis (2021), Theorem 5) Suppose that the assumptions of Theorem 3.5 are fulfilled. Then, there exist constants $C_1$, $C_2 \ge 24$, depending on $(\phi, \mathcal{V})$, such that for all $C_4 \ge C_1$ depending on $\phi$, and for all $t > 0$:

$$ \mathbb{P} \left\{ \sup_{s \in S_0} \left| \frac{1}{n'} \hat{W}_{n',m'}^{\phi}(s) - W_{\phi}(s) \right| > t \right\} \le C_2 \exp\{-C_3 p N' t^2\}, $$

provided that $C_1/\sqrt{pN'} \le t \le C_4(p \wedge (1-p))$, where $C_3 = \log(1+C_4/(4C_1))/(C_2C_4) > 0$ depends on $\phi, \mathcal{V}$.

Applying the theorem above with $t = (\varepsilon - \delta)/8 = C_4(p \wedge (1-p))$, we get:

$$ \mathbb{P}_{H,G} \left\{ 2 \sup_{s \in S_0} \left| \frac{1}{n'} \hat{W}_{n',m'}(s) - W_{\phi}(s) \right| \ge \frac{\varepsilon - \delta}{4} \right\} \\ \le C_2 \exp \left\{ -\frac{N' p (p \wedge (1-p))}{8C_2} (\varepsilon - \delta) \log \left( 1 + \frac{\varepsilon - \delta}{32C_1 (p \wedge (1-p))} \right) \right\}, \quad (49) $$

as soon as $N' \ge 16C_1^2/(p(\varepsilon - \delta)^2)$. Combining (49) with (48) and (47), we obtain the desired bound

$$ \mathbb{P}_{H,G} \left\{ \Phi_{\alpha}^{\phi}(\mathcal{D}_{n'',m''}(\hat{s})) = 0 \right\} \le 18 \exp \left\{ \frac{-CN'' (\varepsilon - \delta)^2}{16} \right\} \\ + C_2 \exp \left\{ -\frac{N' p (p \wedge (1-p))}{8C_2} (\varepsilon - \delta) \log \left( 1 + \frac{\varepsilon - \delta}{32C_1 (p \wedge (1-p))} \right) \right\}. $$

# Appendix C: Additional material

## C.1. Example 1 - continued

We illustrate the departure condition formulated in Definition 3.3 in the particular simple model described in Example 1, and following Remark 8. We show that the quantity $W_{\phi}^* - \int \phi(u)du$ may permit to quantify (even very small) departures from the null assumption in a very rich fashion. Consider no model bias, i.e., the parameter $\delta$ to be set to 0. Then, by plugging in the explicit formulation of ROC* of Eq. (17), in the formula (15), we can easily compute numerically the exact value of the departure from the null in the sense of the $W_{\phi}$-criterion. Precisely, consider specific choices of the score-generating function $\phi$, namely for $\phi(u) = u$, $\phi(u) = u\mathbb{I}\{u \ge u_0\}$ with $u_0 \in (0, 1)$ and $\phi_q(u) = u^q$ with $q \in \mathbb{N}^*$. The intuition is to assign more weight to the X's having higher rank in the pooled sample, which may correspond to observations being more informative to distinguish the two samples. As explained in section 2 and in Clémençon, Limnios and Vayatis (2021) (see section 5 therein), the last two choices of score-generating functions give greater importance to parts of ROC curves with steeper slopes, i.e., near the origin $\alpha = 0$. The numerical results are gathered in Figure 12 below. Figure 12 shows that for a fixed departure from the null assumption, that is, for a fixed value of $\|\theta\|_2$, with $\theta := \theta_+ - \theta_-$, in the present example, the deviation of $W_{\phi}^*$ scales differently depending on the choice of $\phi$. In particular, when fixing $\|\theta\|_2$ (or equivalently $\varepsilon$), it shows that $\phi(u) = u$ (MWW) yields a more conservative test, while for very small departures
<figure><img src="image_12.png" /><figcaption>Fig 12. Graph of the deviation $W_{\phi}^* - \int_{0}^{1} \phi(u)du$ against the departure from the null quantified by $\|\theta\|_2$ with $p = 1/2$. The score-generating functions are $\phi(u) = u$ (MWW), $\phi(u) = u\mathbb{I}\{u \ge u_0\}$ with $u_0 \in \{0.80, 0.90, 0.99\}$, and $\phi_q(u) = u^q$ with $q \in \{3, 5\}$.</figcaption></figure>

$\|\theta\|_2$, choosing $\phi(u) = u\mathbb{I}\{u \ge u_0\}$ or $\phi_q(u) = u^q$ with large $u_0$ and $q$, would empirically result in a test with higher power. As expected, when only considering the most informative observations $\mathbf{X}$, i.e., those related to the $(1-u_0)$-quantile for small $u_0$, the resulting criterion is less sensitive to the difficulty of the problem, i.e., to small values of $\|\theta\|_2$.

## C.2. Numerical experiments

This section gathers additional material to section 4, and in particular numerical illustration of ROC* for the models considered, followed by additional numerical experiments in the last paragraph.

### Exact parameters for the synthetic datasets

* (L1-) For $d = 4$, the diagonals of $\Sigma$ are (2, 6, 1, 5), (-1, 0, 0), (-1, 0), (-1); $d = 6$. For $d = 6$, the matrix diagonals are equal to (2, 6, 1, 5, 4, 3), (-1, 0, 0, 0, 0), (-1, 0, 0), (-1, 0, 0), (-1, 0), (-1).

* (L1+) For $d = 4$, the diagonals of $\Sigma$ are (6, 4, 5, 3), (-2, 4, 2), (-3, 0), (-2); $d = 6$. For $d = 6$, the matrix diagonals are equal to (6, 5, 5, 3, 2, 3), (-2, 4, 2, 1, 0), (-3, 0, 0, 1), (-2, 1, 1), (-3, 2), (-2).

### ROC* curves for the location and the scale models

**Additional numerical experiments** We present additional numerical experiments in this paragraph following the experimental set up of section 4.

(T4) Location model. We generate $\mathbf{X}$ from a multivariate Student distribution of location vector $\mu_X = (\varepsilon/\sqrt{d}) \times \mathbf{1}_d$ and scale matrix that of (L1-), and degrees of freedom equal to $n$. $\mathbf{Y}$ is drawn from a centered multivariate Gaussian $\mu_Y = 0_d$ of equal variance matrix, with $\varepsilon \in \{0.1, 0.3, 0.5\}$ and $d = 6$. Results are showed in Figure 14.
<figure><img src="image_13.png" /><figcaption>Fig 13. True ROC curves (ROC*) for the location and the scale models depending on the discrepancy parameter $\epsilon \in \{0.2, 0.3, 0.6\}$ for (L1+) and (L1-), $\epsilon \in \{0.2, 0.4, 0.6\}$ for (S1) and $\epsilon \in \{0.05, 0.1, 0.2\}$ for (S2).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig 14. Rejection rates of the homogeneity assumption under alternatives against significance level $\alpha \in (0, 1)$ for the model (T4) with $\phi(u) = u$. Parameters: $n = m = 1000$, $d = 6$, $B = 100$, $\epsilon \in \{0.1, 0.3, 0.5\}$. The ranking algorithms are rSVM2, rForest, rBoost (left) compared to MMD, Energy and Tukey (right).</figcaption></figure>

(L1) We generate **X** and **Y** from models (L1+) and (L1-), and where the score-generating function implemented for the *Step 2* is RTB: $\phi: u \mapsto uI\{u \ge u_0\}$, with $u_0 \in \{0.7, 0.8, 0.9\}$. We used rForest for *Step 1* as bipartite ranking algorithm. We show the graphs of the empirical rejection rate of the homogeneity assumption under the alternatives w.r.t. the level of the test $\alpha \in (0, 1)$. Figure 15 shows the results.

## Acknowledgments

The authors would like to thank both referees, the Associate Editor and the Editor Grace Y. Yi for constructive comments that improved the paper.
<caption>Table 4</caption>

Rejection rates of the homogeneity assumption under $H_1$ at significance level $\alpha = 0.05$ for the location model (T4), with $n = m = 1000$, $\pm$ their standard deviation at 95%. Methods in bold maximize among all algorithms the rejection rate under $H_1$.

<table><thead><tr><th rowspan="2">Rejection rate of the null<br>Method</th><th colspan="3">Under H<sub>1</sub>, (d = 6)</th></tr><tr><th>&epsilon; = 0.1</th><th>&epsilon; = 0.3</th><th>&epsilon; = 0.5</th></tr></thead><tbody><tr><td>rSVM</td><td><b>0.99 &pm; 0.010</b></td><td><b>1.00 &pm; 0.000</b></td><td><b>1.00 &pm; 0.000</b></td></tr><tr><td>rForest</td><td>0.14 &pm; 0.038</td><td>1.00 &pm; 0.000</td><td>1.00 &pm; 0.000</td></tr><tr><td>rBoost</td><td>0.17 &pm; 0.038</td><td>0.78 &pm; 0.042</td><td><b>1.00 &pm; 0.00</b></td></tr><tr><td>MMD</td><td>0.06 &pm; 0.024</td><td>0.06 &pm; 0.024</td><td>0.06 &pm; 0.024</td></tr><tr><td>Energy</td><td>0.03 &pm; 0.017</td><td>0.03 &pm; 0.017</td><td>0.01 &pm; 0.010</td></tr><tr><td>Tukey</td><td>0.04 &pm; 0.020</td><td>0.14 &pm; 0.035</td><td>0.18 &pm; 0.039</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Fig 15. Rejection rates of the homogeneity assumption under alternatives against significance level $\alpha \in (0, 1)$ for the models (L1-) (a) and (L1+) (b) under $H_1$. Step 1 with rForest, Step 2 with score-generating function RTB $\phi: u \mapsto u\mathbb{I}\{u \ge u_0\}$ with $u_0 \in \{0.7, 0.8, 0.9\}$ resp. corresponding to RTB7, RTB8, RTB9. Comparison of the curves with MWW.</figcaption></figure>

## Funding

Stephan Clémençon was supported by PEPR IA Foundry and ANR project Far-See. Myrto Limnios was supported by Novo Nordisk Foundation Grant NNF20OC0062897.
References

AGARWAL, S., GRAEPEL, T., HERBRICH, R., HAR-PELED, S. and ROTH, D. (2005). Generalization Bounds for the Area Under the ROC Curve. *Journal of Machine Learning Research* **6** 393–425. MR2249826

BACH, F., HARCHAOUI, Z. and MOULINES, E. (2008). Testing for homogeneity with kernel Fischer discriminant analysis. In *Advances in Neural Information Processing Systems* 2008. MIT Press, Cambridge, MA.

BEIRLANT, J., BUITENDAG, S., DEL BARRIO, E., HALLIN, M. and KAMPER, F. (2020). Center-outward quantiles and the measurement of multivariate risk. *Insurance: Mathematics and Economics* **95** 79–100. MR4154446

BERTAIL, P., CLÉMENÇON, S., GUYONVARCH, Y. and NOIRY, N. (2021). Learning from Biased Data: A Semi-Parametric Approach. In *Proceedings of the 38th International Conference on Machine Learning. Proceedings of Machine Learning Research* **139** 803–812. PMLR.

BIAU, G. and GYORFI, L. (2005). On the asymptotic properties of a nonparametric $l_1$-test statistic of homogeneity. *IEEE Transactions on Information Theory* **51** 3965–3973. MR2239012

BURGES, C., RAGNO, R. and LE, Q. (2007). Learning to Rank with Nonsmooth Cost Functions. In *Advances in Neural Information Processing Systems* **19**. MIT Press.

BURGES, C., SHAKED, T., RENSHAW, E., LAZIER, A., DEEDS, M., HAMILTON, N. and HULLENDER, G. (2005). Learning to Rank Using Gradient Descent. In *Proceedings of the 22nd International Conference on Machine Learning* 89–96.

CARPENTIER, A., COLLIER, O., COMMINGES, L., TSYBAKOV, A. and WANG, Y. (2018). Minimax Rate of Testing in Sparse Linear Regression. *Automation and Remote Control* **80**. MR4032408

CHAKRABORTY, A. and CHAUDHURI, P. (2015). A Wilcoxon-Mann-Whitney-type test for infinite-dimensional data. *Biometrika* **102** 239–246. MR3335109

CHAKRABORTY, A. and CHAUDHURI, P. (2017). Tests for high-dimensional data based on means, spatial signs and spatial ranks. *The Annals of Statistics* **45** 771–799. MR3650400

CHAUDHURI, P. (1996). On a geometric notion of quantiles for multivariate data. *Journal of the American Statistical Association* **91** 862–872. MR1395753

CHERNOFF, J. and SAVAGE (1958). Asymptotic normality and efficiency of certain non-parametric test statistics. *Ann. Math. Stat.* **29** 972–994. MR0100322

CHERNOZHUKOV, V., GALICHON, A., HALLIN, M. and HENRY, M. (2017). Monge–Kantorovich depth, quantiles, ranks and signs. *The Annals of Statistics* **45** 223–256. MR3611491

CHEUNG, Y. K. and KLOTZ, J. H. (1997). The Mann Whitney Wilcoxon distribution using linked list. *Statistica Sinica* **7** 805–813.

CLÉMENÇON, S., DEPECKER, M. and VAYATIS, N. (2009). AUC maximization and the two-sample problem. In *Advances in Neural Information Processing Systems. Lecture Notes in Computer Science* **3559** 1–15. Springer.

CLÉMENÇON, S., DEPECKER, M. and VAYATIS, N. (2011). Adaptive partitioning schemes for bipartite ranking. *Machine Learning* **43** 31–69. MR3108203

CLÉMENÇON, S., DEPECKER, M. and VAYATIS, N. (2013a). Ranking Forests. *Journal of Machine Learning Research* **14** 39–73. MR3033325
CLÉMENÇON, S., DEPECKER, M. and VAYATIS, N. (2013b). An Empirical Comparison of Learning Algorithms for Nonparametric Scoring – The TreeRank Algorithm and Other Methods. *Pattern Analysis and its Applications* **16** 475–496. MR3117191

CLÉMENÇON, S., LIMNIOS, M. and VAYATIS, N. (2021). Concentration inequalities for two-sample rank processes with application to bipartite ranking. *Electronic Journal of Statistics* **15** 4659–4717. MR4317718

CLÉMENÇON, S., LUGOSI, G. and VAYATIS, N. (2008). Ranking and empirical risk minimization of U-statistics. *The Annals of Statistics* **36** 844–874. MR2396817

CLÉMENÇON, S. and VAYATIS, N. (2007). Ranking the best instances. *Journal of Machine Learning Research* **8** 2671–2699. MR2383571

CLÉMENÇON, S. and VAYATIS, N. (2009a). Tree-based ranking methods. *IEEE Transactions on Information Theory* **55** 4316–4336. MR2582885

CLÉMENÇON, S. and VAYATIS, N. (2009b). Empirical performance maximization based on linear rank statistics. In *Advances in Neural Information Processing Systems. Lecture Notes in Computer Science* **3559** 1–15. Springer.

CLÉMENÇON, S. and VAYATIS, N. (2009c). Adaptive estimation of the optimal ROC curve and a bipartite ranking algorithm. In *Proceedings of ALT’09*. MR2564229

CLÉMENÇON, S. and VAYATIS, N. (2010). Overlaying Classifiers: a practical approach to optimal scoring. *Constructive Approximation* **32** 619–648. MR2726448

DEB, N. and SEN, B. (2021). Multivariate Rank-based Distribution-free Nonparametric Testing using Measure Transportation. *Journal of the American Statistical Association* 1–16. MR4571116

FREUND, Y., IYER, R. D., SCHAPIRE, R. E. and SINGER, Y. (2003). An efficient boosting algorithm for combining preferences. *Journal of Machine Learning Research* **4** 933–969. MR2125342

FRIEDMAN, J. H. and RAFSKY, L. C. (1979). Multivariate Generalizations of the Wald-Wolfowitz and Smirnov Two-Sample Tests. *The Annals of Statistics* **7** 697–717. MR0532236

GOLDREICH, O., GOLDWASSER, S. and RON, D. (1998). Property Testing and Its Connection to Learning and Approximation. *Journal of the ACM* **45** 653–750. MR1675099

GRETTON, A., BOUSQUET, O., SMOLA, A. and SCHÖLKOPF, B. (2005a). Measuring Statistical Dependence with Hilbert-Schmidt Norms. In *Algorithmic Learning Theory: 16th International Conference, ALT 2005* **3734** 63–78. MR2255909

GRETTON, A., SMOLA, A., BOUSQUET, O., HERBRICH, R., BELITSKI, A., AUGATH, M., MURAYAMA, Y., PAULS, J., SCHÖLKOPF, B. and LOGOTHETIS, N. (2005b). Kernel Constrained Covariance for Dependence Measurement. In *Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics* **R5** 112–119. Proceedings of Machine Learning Research.

GRETTON, A., BORGWARDT, K. M., RASCH, M. J., SCHOLKOPF, B. and SMOLA, A. (2007). A kernel method for the two-sample problem. In *Advances in Neural Information Processing Systems 19*. MIT Press, Cambridge, MA.

GRETTON, A., BORGWARDT, K. M., RASCH, M. J., SCHOLKOPF, B. and SMOLA, A. (2012). A kernel two-sample problem. *Journal of Machine Learning Research* **13** 723–773. MR2913716

HÁJEK, J. and SIDÁK, Z. (1967). *Theory of Rank Tests*. Academic Press. MR0229351
HALLIN, M. and PAINDAVEINE, D. (2002a). Optimal tests for multivariate location based on interdirections and pseudo-Mahalanobis ranks. *The Annals of Statistics* **30** 1103–1133.
MR1926170

HALLIN, M. and PAINDAVEINE, D. (2002b). Optimal procedures based on interdirections and pseudo-Mahalanobis ranks for testing multivariate elliptic white noise against ARMA dependence. *Bernoulli* **8** 787–815. MR1963662

HALLIN, M. and PAINDAVEINE, D. (2008). Optimal rank-based tests for homogeneity of scatter. *The Annals of Statistics* **36** 1261–1298. MR2418657

HALLIN, M., DEL BARRIO, E., CUESTA-ALBERTOS, J. and MATRÁN, C. (2021). Distribution and quantile functions, ranks and signs in dimension d: A measure transportation approach. *The Annals of Statistics* **49** 1139–1165. MR4255122

HODGES, J. L. and LEHMANN, E. L. (1956). The Efficiency of Some Nonparametric Competitors of the t-Test. *The Annals of Mathematical Statistics* **27** 324–335. MR0079383

HOEFFDING, W. (1963). Probability Inequalities for Sums of Bounded Random Variables. *Journal of the American Statistical Association* **58** 13–30. MR0144363

HUANG, K. H., LIU, X., DUNCAN, A. B. and GANDY, A. (2023). A High-dimensional Convergence Theorem for U-statistics with Applications to Kernel-based Testing. In *Proceedings of Thirty Sixth Conference on Learning Theory. Proceedings of Machine Learning Research* **195** 3827–3918. *Proceedings of Machine Learning Research.*

JAKOBSEN, M. E. (2017). Distance Covariance in Metric Spaces: Non-Parametric Independence Testing in Metric Spaces (Master’s thesis).

JOACHIMS, T. (2002). Optimizing Search Engines Using Clickthrough Data. In *Proceedings of the Eighth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. KDD '02* 133–142. Association for Computing Machinery.

KOSHEVOY, G. and MOSLER, K. (1997). Zonoid trimming for multivariate distributions. *The Annals of Statistics* **25** 1998–2017. MR1474078

LAM-WEIL, J., CARPENTIER, A. and SRIPERUMBUDUR, B. K. (2022). Local minimax rates for closeness testing of discrete distributions. *Bernoulli* **28** 1179–1197. MR4388934

LEHMANN, E. L. and ROMANO, J. P. (2005). *Testing Statistical Hypotheses*. Springer.
MR2135927

LEUNG, D. and DRTON, M. (2018). Testing independence in high dimensions with sums of rank correlations. *The Annals of Statistics* **46** 280–307. MR3766953

LIMNIOS, M. (2022). Rank processes and statistical applications in high dimension, Theses, Université Paris-Saclay.

LIMNIOS, M. and CLÉMENÇON, S. (2024). On Ranking-based Tests of Independence. In *Proceedings of The 27th International Conference on Artificial Intelligence and Statistics. Proceedings of Machine Learning Research* **238** 577–585. PMLR.

LIU (1990). On a notion of data depth based upon random simplices. *The Annals of Statistics* **18** 405–414. MR1041400

LIU, R. Y. and SINGH, K. (1993). A quality index based on data depth and multivariate rank tests. *Journal of the American Statistical Association* **88** 252–260. MR1212489

LUNG-YUT-FONG, A., LÉVY-LEDUC, C. and CAPPÉ, O. (2015). Homogeneity and change-point detection tests for multivariate data using rank statistics. *Journal de la société française de statistique* **156** 133–162. MR3436651

LYONS, R. (2013). Distance covariance in metric spaces. *The Annals of Probability* **41** 3284–3305. MR3127883
MENON, A. K. and WILLIAMSON, R. C. (2016). Bipartite Ranking: A Risk Theoretic Perspective. *Journal of Machine Learning Research* **7** 1–102. MR3580348

MOSLER, K. (2013). Depth statistics. In *Robustness and Complex Data Structures: Festschrift in Honour of Ursula Gather* 17–34. Springer.

MÖTTÖNEN, J. and OJA, H. (1995). Multivariate spatial sign and rank methods. *Journal of Nonparametric Statistics* **5** 201–213. MR1346895

MÖTTÖNEN, J., OJA, H. and TIENARI, J. (1997). On the Efficiency of Multivariate Spatial Sign and Rank Tests. *The Annals of Statistics* **25** 542–552. MR1439313

MÖTTÖNEN, J., OJA, H. and SERFLING, R. (2005). Multivariate generalized spatial signed-rank methods. *Journal of Statistical Research* **39** 19–42. MR2195206

NOLAN, D. and POLLARD, D. (1987). *U-Processes: Rates of Convergence*. *The Annals of Statistics* **15** 780–799. MR0888439

OJA, H. (1983). Descriptive statistics for multivariate distributions. *Statistics and Probability Letters* **1** 327–332. MR0721446

OJA, H. (2010). *Multivariate Nonparametric Methods with R: An approach based on spatial signs and ranks*. *Lecture Notes in Statistics* **199**. Springer-Verlag New York. MR2598854

PURI, M. L. and SEN, P. K. (1971). *Nonparametric methods in multivariate analysis*. John Wiley and Sons. MR0298844

RACHEV, A. (1991). *Probability Metrics and the Stability of Stochastic Models*. Wiley. MR1105086

RAKOTOMAMONJY, A. (2004). Optimizing Area Under ROC Curve with SVMs. In *Proceedings of the First Workshop on ROC Analysis in AI*.

RAMDAS, A., REDDI, S. J., PÓCZOS, B., SINGH, A. and WASSERMAN, L. (2015). On the Decreasing Power of Kernel and Distance Based Nonparametric Hypothesis Tests in High Dimensions. In *Proceedings of the Twenty-Ninth AAAI Conference on Artificial Intelligence*. AAAI’15 3571–3577. AAAI Press.

RUBINFELD, R. and SUDAN, M. (1996). Robust Characterizations of Polynomials with Applications to Program Testing. *SIAM Journal on Computing* **25** 252–271. MR1379300

RUDIN, C. (2006). Ranking with a P-Norm Push. In *Proceedings of COLT 2006. Lecture Notes in Computer Science* **4005** 589–604.

RUDIN, C., CORTES, C., MOHRI, M. and SCHAPIRE, R. E. (2005). Margin-Based Ranking and Boosting Meet in the Middle. In *Proceedings of COLT 2005. Lecture Notes in Computer Science* **3559** 63–78. Springer.

SCHRAB, A., KIM, I., ALBERT, M., LAURENT, B., GUEDJ, B. and GRETTON, A. (2023). MMD Aggregated Two-Sample Test. *Journal of Machine Learning Research* **24** 1–81. MR4633583

SEJDINOVIC, D., SRIPERUMBUDUR, B., GRETTON, A. and FUKUMIZU, K. (2013). Equivalence of distance-based and RKHS-based statistics in hypothesis testing. *The Annals of Statistics* **41** 2263–2291. MR3127866

SERFLING, R. J. (1980). *Approximation theorems of mathematical statistics*. John Wiley and Sons. MR0595165

SHI, H., HALLIN, M., DRTON, M. and HAN, F. (2022). On universally consistent and fully distribution-free rank tests of vector independence. *The Annals of Statistics* **50** 1933–1959. MR4474478

SZÉKELY, G. J. and RIZZO, M. L. (2007). Measuring and testing dependence by correlation of distances. *The Annals of Statistics* **35** 2769–2794. MR2382665
SZÉKELY, G. J. and RIZZO, M. L. (2013). Energy statistics: A class of statistics based on distances. *Journal of Statistical Planning and Inference* **143** 1249–1272. MR3055745

TUKEY, J. W. (1975). Mathematics and the picturing of data. In *Proceedings of the International Congress of Mathematicians* **2** 523–531. MR0426989

UM, Y. and RANDLES, R. H. (1998). Nonparametric Tests for the Multivariate Multi-sample Location Problem. *Statistica Sinica* **8** 801–812. MR1651509

VAN DER VAART, A. (1998). *Asymptotic Statistics*. Cambridge University Press. MR1652247

VAN DER VAART, A. and WELLNER, J. (1996). *Weak Convergence and Empirical Processes*. Springer-Verlag New York. MR1385671

VARDI, Y. and ZHANG, C. H. (2000). The multivariate $L_1$-median and associated data depth. *Proceedings of the National Academy of Sciences* **97** 1423–1426. MR1740461

WANG, M., DENG, W., HU, J., TAO, X. and HUANG, Y. (2019). Racial Faces in the Wild: Reducing Racial Bias by Information Maximization Adaptation Network. In *2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019* 692–702. IEEE.

WILCOXON, F. (1945). Individual comparisons by ranking methods. *Biometrics* **1** 80–83. MR0025133

ZUO, Y. and SERFLING, R. (2000). General notions of statistical depth function. *The Annals of Statistics* **28** 461–482. MR1790005