On the Finite-Sample Performance of
Measure Transportation-Based
Multivariate Rank Tests
Multivariate rank tests

Marc Hallin and Gilles Mordant

Marc Hallin ECARES and Département de Mathématique, Université libre de Bruxelles,
Avenue F.D. Roosevelt 50, 1050
Brussels, Belgium
Gilles Mordant IMS, Universität Göttingen,
Goldschmidtstraße 7, 37077
Göttingen, Germany

Extending to dimension 2 and higher the dual univariate concepts of ranks and quantiles has remained an open problem for more than half a century. Based on measure transportation results, a  solution has been proposed recently under the name center-outward ranks and quantiles which, contrary to previous proposals, enjoys all the properties that make univariate ranks a successful tool for statistical inference. Just as their univariate counterparts (to which they reduce in dimension one), center-outward ranks allow for the construction of distribution-free and asymptotically efficient tests for a variety of problems where the density of some noise or innovation remains unspecified. The actual implementation of these tests involves the somewhat arbitrary choice of a grid. While the asymptotic impact of that choice is nil, its  finite-sample consequences are not. In this note, we investigate the  finite-sample impact of that choice in the typical context of the multivariate two-sample location problem.

§ INTRODUCTION

§.§ David Tyler, beyond affine equivariance and elliptical symmetry

The closely related concepts of affine equivariance and elliptical symmetry  played a central role in the development of robust multivariate statistics over the past 60 years.[[Tukey, 1960], [Huber, 1964], and [Hampel, 1968] generally are considered as   laying the foundations of modern robust statistics; see [Ronchetti, 2006] for a historical perspective and [Stigler, 1973] for an account of the pre-Tukey era.] A critical  attitude towards this dominant role of elliptical densities constitutes a red thread running through all of David's contributions to   multivariate anlysis[Significantly, “Robust Multivariate Statistics: Beyond Ellipticity and Affine Equivariance” is the title of one of David's NSF grants.]—an attitude that actually  takes  place in a broader debate on the ordering of the real space in dimension $d\geq 2$. Such  ordering  is an essential issue if the univariate concepts of distribution and quantile functions, ranks, and signs,   all heavily depending on the canonical ordering of the real line,  are to be extended to a multivariate context.

§.§ Ordering the real space in dimension $d\geq 2$

The problem of ordering $\mathbb R ^d$ for $d\geq 2$, hence ranking multivariate observations,  has a long history in statistics. Many attempts have been made to define adequate  multivariate concepts of ranks.

The notion of  ranks, however, is not an isolated one, as it is inseparable from that of empirical quantiles, quantile regions (collections of points with ranks less than or equal to  some given value), and quantile contours (collections of points with ranks  equal to some given value). A sound definition thus should include all these concepts, along with their population versions—the population distribution and quantile functions $F$ and $Q$—and their mutual relations (a  quantile function is the inverse of a distribution function $(Q=F^{-1})$; a population distribution function and its empirical version are asymptotically related via a Glivenko-Cantelli result, etc.).  Among the key properties of any successful concept are the distribution-freeness (within the class of absolutely continuous distributions $\rm P$, say) of the ranks and the distribution-freeness of the push-forward[We adopt here the convenient terminology and notation of measure transportation: the push-forward $F\#{\rm P}$ of $\rm P$ by $F$ is the distribution of $F(Z)$ where $Z\sim{\rm P}$, i.e.,  $F(Z)\sim F\#{\rm P}$ if $Z\sim\rm P$.]    of a distribution $\rm P$ by its distribution function $F$. Without these distribution-freeness properties,  the level of a quantile $Q(\tau)=F^{-1}(\tau)$  depends on the distribution $\rm P$ characterized by $F$ and can be anything  larger or smaller than $\tau$: as a quantile of order $\tau$, thus, $Q(\tau)$ is totally meaningless.

Appealing as they are, none of the attempts that had been made until recently—marginal ranks, spatial ranks, elliptical (or Mahalanobis) ranks, ... —is satisfying the desired properties; actually, none of them is  even enjoying distribution-freeness![The Mahalanobis ranks and signs [Hallin & Paindaveine, 2002] are enjoying distribution-freeness over the class of elliptically symmetric distributions only.] Nor do the various depth  concepts: the probability content of a depth contour of given depth strongly depends on the underlying $\rm P$, which hinders its interpretation as a quantile contour.

Based on measure transportation results (mainly, a theorem by [McCann, 1995]) [Chernozhukov et al., 2017], [Hallin, 2017], and [Hallin et al., 2021] recently introduced the  concepts  of center-outward ranks and signs, distribution and quantile functions which, for the first time, satisfy all the desired properties (see   [Hallin et al., 2021]  and the review by [Hallin, 2022] for details) and   further triggered the development of several  appealing multivariate, distribution-free statistical procedures, among which   Faugeras & Rüschendorf, 2017, Shi et al., 2021, Shi et al., 2021, Ghosal & Sen, 2019, Deb & Sen, 2021, Deb et al., 2020, Deb et al., 2021, Masud & Aeron, 2021, Hallin et al., 2020, Hallin et al., 2020, Hallin et al., 2020.
These are the concepts we are considering in this paper and now describe under their various versions.

§ CENTER-OUTWARD RANKS AND SIGNS

For the simplicity of exposition, we throughout consider  distributions $\rm P$ on $\mathbb R ^d$ in the family $\mathcal P ^d$ of Lebesgue-absolutely continuous distributions with nonvanishing densities,  that is,  with a density $f$ such that for all $B>0$ there exist $0<m^-_B\leq~\!m^+_B<~\!\infty$ such that $m^-_B\leq f(\zb)\leq m^+_B$ for all $\zb$ such that $\Vert\zb\Vert\leq B$. That assumption can be relaxed, though: see [del Barrio et al., 2020].

§.§ A   measure transportation-based concept of distribution and quantile functions

The basic idea behind the definition of the center-outward distribution and quantile functions of a probability measure $\rm P\in\mathcal P ^d$ is quite simple. For $d=1$, the distribution function $F$ of $\rm P$ is the unique monotone increasing function pushing $\rm P$ forward to the uniform ${\rm U}_{[0,1]}$ over $[0,1]$—namely, $F\#{\rm P}={\rm U}_{[0,1]}$. Rather than $F$, however, which is based on a left-to-right ordering of $\mathbb R$ that does not extend to $\mathbb R ^d$ for $d\geq 2$, we consider the center-outward distribution function $F_\pm:=2F-1$, which contains the same information as $F$ and is the unique monotone increasing function pushing $\rm P$ forward to the uniform ${\rm U}_{[-1,1]}$ over $[-1,1]$. A monotone increasing function is the gradient (the derivative) of a convex function: the center-outward distribution function $F_\pm$ this is the unique gradient of a convex function such that $F_\pm\#{\rm P}={\rm U}_{[-1,1]}$.
The interval $[-1,1]$
is, for $d=1$, the closed unit ball $\bar{\mathbb S}_d$, where ${\mathbb S}_d:=\big\{\ub \vert \Vert\ub\Vert \leq 1\big\}$ and, denoting by ${\rm U}_d$ the spherical uniform[The spherical uniform ${\rm U}_d$ over $\bar{\mathbb S}_d$ is the spherical distribution with center $\mathbf 0$ and  radial density the uniform over $[0,1]$: it is thus the product of a uniform over $[0,1]$ for the distances to the origin and a uniform over the unit (hyper)sphere for the directions.] over  $\bar{\mathbb S}_d$, the spherical uniform over $\bar{\mathbb S}_1$ coincides with the uniform over $[-1,1]$:  namely, ${\rm U}_1={\rm U}_{[-1,1]}$.

A celebrated theorem by  [McCann, 1995] tells us that, for arbitrary dimension $d\in~\!\mathbb N$ and arbitrary $\rm P\in\mathcal P ^d$, there exist a ($\rm P$-a.s., here Lebesgue-a.e.) unique gradient  of a convex function $\Fb_\pm$ such that $\Fb_\pm\#{\rm P}={\rm U}_d$. Obviously, for $d=1$, $\Fb_\pm$ coincides with the univariate $F_\pm$, whence the notation. Call $\Fb_\pm$ the center-outward distribution function of $\rm P$. It follows from [Figalli, 2018] that—except perhaps at $\Fb_\pm^{-1}({\mathbf 0})$ (a set of points with Lebesgue measure zero)—$\Fb_\pm\#{\rm P}$ is a homeomorphism, hence admits a continuous (except perhaps at $\Fb_\pm^{-1}({\mathbf 0})$) inverse $\Qb_\pm:=\Fb_\pm^{-1}$: call $\Qb_\pm$ the center-outward quantile function of $\rm P$. Clearly, $\Qb_\pm\#{\rm U}_d={\rm P}$.

This, with the spherical uniform ${\rm U}_d$ as a reference distribution (extending ${\rm U}_{[-1,1]}$)    is the concept proposed in  [Hallin, 2017] and [Hallin et al., 2021], where  we refer to for further properties of $\Fb_\pm$ and $\Qb_\pm$ justifying their qualification as distribution and quantile functions.

Other   choices are possible for the reference $\rm U$, though.  Replacing ${\rm U}_d$ with an arbitrary compactly supported  absolutely continuous reference distribution $\rm U$, [Chernozhukov et al., 2017], in a very general  approach, propose, under the name of Monge-Kantorovich  vector rank and  Monge-Kantorovich  quantile functions,  measure transportation-based definitions of a broad class of  analogues, $\Fb_{\text{\rm MK}}$ and $\Qb_{\text{\rm MK}}$, say, of $\Fb_\pm$ and $\Qb_\pm$. For nonspherical $\rm U$'s, however, the Monge-Kantorovich  quantile functions do not enjoy all the features expected from a quantile function;[On this point, see Section 3.4 in [Hallin, 2022].]  [Chernozhukov et al., 2017] therefore also introduce a concept of Monge-Kantorovich  depth $D_{\text{\rm MK}}$—a transformation-retransformation version (based on the Monge-Kantorovich  vector rank function) of  classical Tukey depth $D_{\text{\rm Tukey}}$.  For spherical $\rm U$'s, Monge-Kantorovich  depth and quantile contours coincide. More precisely, defining  $\delta (\tau):=D_{\text{\rm Tukey}}(\ub_\tau)$ where $\ub_\tau$ is such that ${\rm U}\big( \big\{ \ub\,\big\vert\, \Vert\ub\Vert\leq\Vert\ub_\tau\Vert\big\}\big) = \tau$, one has $\Big\{\zb\, \big\vert\,  \Vert\Fb_{\text{\rm MK}}\Vert =~\!\tau\Big\}=~\!D_{\text{\rm MK}}^{-1}(\delta(\tau))$. Recurring to depth in order to construct quantile regions and contours, thus, is not necessary in the case of  a spherical reference $\rm U$ which, in that respect, offers a better conceptual coherence between the resulting notions of vector ranks and quantiles. As far as rank tests are concerned, however, this can be considered a minor concern.

The choice for $\rm U$ of the nonspherical  Lebesgue uniform ${\rm U}_{[0,1]^d}$ over the unit  (in the canonical basis) hypercube $[0,1]^d$—call it the cubic uniform—yields a vector rank function $\Fb_{\text{\rm MK}}$  that  reduces, for $d=1$, to the classical   distribution   function $F$ just as $\Fb_\pm$ reduces to $F_\pm$. Despite poor equivariance properties[Contrary to $\Fb_\pm$, which is nicely equivariant, the rank vector function $\Fb_{\text{\rm MK}}$ associated with  the cubic uniform ${\rm U}_{[0,1]^d}$ is highly non-equivariant under orthogonal transformations.] its use has been advocated by several authors: see, e.g., [Faugeras & Rüschendorf, 2017], [Carlier et al., 2016],
[Deb et al., 2020, Deb et al., 2021], [Deb & Sen, 2021].

§.§ Multivariate ranks and signs

Denote by $\Zb\n:=(\Zb\n_1,\ldots,\Zb\n_n)$ an i.i.d. sample with distribution $\rm P \in \mathcal P ^d$. The empirical counterpart $\Fb_\pm\n$ of $\Fb_\pm$ is obtained as the solution of an optimal pairing problem between the sample values $\Zb\n_1,\ldots,\Zb\n_n$ and a “regular” grid ${\mathfrak G}\n$ with gridpoints ${\scriptstyle{\mathfrak{G}}}\n_1,\ldots,{\scriptstyle{\mathfrak{G}}}\n_n$. Precisely, $\left(\Fb_\pm\n(\Zb\n_1),\ldots, \Fb_\pm\n(\Zb\n_n)\right)$ is defined as the minimizer  $\left(
{\scriptstyle{\mathfrak{G}}}\n_{\pi^*(1)},\ldots,{\scriptstyle{\mathfrak{G}}}\n_{\pi^*(n)}\right)$, over the $n!$ possible permutations $\pi\in\Pi_n$ of the integers $\{1,\ldots,n\}$, of
%\sum_{i=1}^n\big\Vert {\bf Z}\n_i - {\bf F} _{\!{\pm}}\n({\bf Z}\n_i)
% \min_{\pi\in\Pi_n}
\sum_{i=1}^n\big\Vert {\bf Z}\n_{i} - {\scriptstyle{\mathfrak{G}}}\n_{\pi(i)})
\big\Vert ^2

The choice of the grid ${\mathfrak G}\n$, of course, depends on the reference distribution $\rm U$ adopted in the definitions of Section <ref>: in particular, the uniform discrete distribution over the $n$ gridpoints ${\scriptstyle{\mathfrak{G}}}\n_1,\ldots,{\scriptstyle{\mathfrak{G}}}\n_n$ should converge weakly to $\rm U$ as $n\to\infty$. Our objective is to investigate the finite-sample performance of two-sample location tests based on

(Ti) $\:\,\,$the empirical center-outward distribution function ${\bf F}\n_\pm$ associated with the spherical uniform reference distribution $\rm U =\rm U _d$;

(Tii) $\:$the empirical Monge-Kantorovich vector ranks ${\bf F}\n_{\scriptscriptstyle{\square}}$ associated  with the cubic uniform  reference distribution $\rm U =\rm U _{[0,1]^d}$;

(Tiii) $\:$the empirical Monge-Kantorovich vector ranks ${\bf F}\n_{\pm{\cal N}}$ associated  with the Gaussian ${\mathcal N}({\mathbf 0}, {\mathbf I}_d)$ reference considered as a spherical distribution;

(Tiv) $\:$the empirical Monge-Kantorovich vector ranks ${\bf F}\n_{\scriptscriptstyle{\square}{\cal N}}$ associated  with the Gaussian ${\mathcal N}({\mathbf 0}, {\mathbf I}_d)$ reference considered as  a product of univariate standard normal distributions.

The grids we are using for these four cases are constructed as follows (see Section <ref> for details on Halton sequences and the choice of $n_R$ and $n_S$):

($\mathfrak G$i) $\:\, \rm U = \rm U _d$: (a) factorize $n$ into $n=n_Rn_S + n_0$ with $n_0<\min(n_R,n_S)$; (b) generate a Halton sequence ${\mathfrak S}_{(n_S)}:=({\bf u}_1,\ldots {\bf u}_{n_S})$ over the unit (hyper)sphere ${\mathcal S _{d-1}}$; (c) the grid ${\mathfrak G}\n$ consists of the intersections of these $n_S$ unit vectors with   the $n_R$ hyperspheres centered at $\mathbf 0$, with radii ${j}/{(n_R  +1)}$,
$j=1,\ldots ,n_R$—along with $n_0$ copies of the origin;

($\mathfrak G$ii) $\:\!\! \rm U = \rm  U _{[0,1]^d}$: the grid ${\mathfrak G}\n$ is a Halton sequence over $[0,1]^d$;

($\mathfrak G$iii) $\:\!\!\!\hspace{-.2mm} \rm U ={\mathcal N}({\mathbf 0}, {\mathbf I}_d)$, spherical grid: the grid ${\mathfrak G}\n$ is the image, by  the radial transformation $\zb\mapsto \sqrt{F_{\chi^2_d}^{-1}(\Vert\zb\Vert)}\zb$,  of the spherical grid constructed in (i), where $F_{\chi^2_d}$ denotes the chi-square distribution function with $d$ degrees of freedom;

($\mathfrak G$iv) $\:\!\!\! \rm U ={\mathcal N}({\mathbf 0}, {\mathbf I}_d)$, cubic grid: the grid ${\mathfrak G}\n$ is the image, by componentwise application of the standard normal quantile transformation $z_i\mapsto \Phi^{-1}(z_i)$,[As usual, we denote by $\Phi$ the standard normal distribution function, by $\Phi^{-1}$ the standard normal quantile function.]
of a Halton sequence over  $[0,1]^d$.

The grid ${\mathfrak G}\n$ in ($\mathfrak G$i) reduces, for $d=1$, to

\[
\left\{\pm 1/(\lceil n/2\rceil +1), \ldots, \pm \lceil n/2\rceil/(\lceil n/2\rceil +1)\right\}
\]

along with the origin in case $n$ is odd; that grid is of the form

\[
\left\{2\left(1/(n+1)
\right)-1,\ldots, 2\left(n/(n+1)
\right)-1
\right\}
\]

where $ \left\{\left(1/(n+1)
\right),\ldots, \left(n/(n+1)
\right)
\right\}
$ is the grid producing traditional univariate ranks to which the grid ${\mathfrak G}\n$ in ($\mathfrak G$ii) also reduces for $d=1$.

In ($\mathfrak G$i) and ($\mathfrak G$iii), the grid ${\mathfrak G}\n$ is  spherical; as a consequence, ${\bf F}\n_\pm(\Zb\n_i)$ and ${\bf F}\n_{\pm {\cal N}}(\Zb\n_i)$ in (Ti) and (Tiii) naturally factorize as

{\bf F}\n_\pm(\Zb\n_i)=\Vert{\bf F}\n_\pm(\Zb\n_i) \Vert\frac{{\bf F}\n_\pm(\Zb\n_i)}{\Vert{\bf F}\n_\pm(\Zb\n_i) \Vert}=: \frac{R\n_{i\pm}}{n_R+1}{\Sb\n_{i\pm}}
\vspace{-1mm} $$

where  $R\n_{i\pm}= (n_R+1)\Vert{\bf F}\n_\pm(\Zb\n_i)\Vert$, ranging from 0 or 1 (according as $n_0\neq 0$ or $n_0=0$) to $n_R$, is the center-outward rank of $\Zb\n_i$ and $\Sb\n_{i\pm}$ (a unit vector) has the interpretation of a (multivariate) center-outward sign and

\begin{equation}\label{vdWDef}
{\bf F}\n_{\pm {\cal N}}(\Zb\n_i)=\Vert{\bf F}\n_{\pm {\cal N}}(\Zb\n_i) \Vert
\frac{{\bf F}\n_{\pm {\cal N}}(\Zb\n_i)}{
\Vert{\bf F}\n_{\pm {\cal N}}(\Zb\n_i) \Vert
}=: J_{\text{\rm vdW}}\left(\frac{R\n_{i\pm{\cal N}}}{n_R+1}\right){\Sb\n_{i\pm{\cal N}}}
\vspace{-1mm} \end{equation}

where $J_{\text{\rm vdW}}= \sqrt{F_{\chi^2_d}^{-1}}$ is the univariate normal or van der Waerden score function, ${R\n_{i\pm{\cal N}}}$ the rank of $\Vert{\bf F}\n_{\pm {\cal N}}(\Zb\n_i) \Vert$ among the $n_R$ distinct values of $\Vert{\bf F}\n_{\pm {\cal N}}(\Zb\n_i) \Vert$ for $i=1,\ldots,n$ and ${\Sb\n_{i\pm{\cal N}}}$ similarly has the interpretation of a multivariate   sign. Being based on different transport maps, however, neither $R\n_{i\pm}$ and ${R\n_{i\pm{\cal N}}}$ nor $\Sb\n_{i\pm}$ and ${\Sb\n_{i\pm{\cal N}}}$ need coincide.

No similar factorization into ranks and signs occurs with the vector ranks ${\bf F}\n_{\scriptscriptstyle{\square}}$ and ${\bf F}\n_{\scriptscriptstyle{\square}{\cal N}}$ in (Tii) and (Tiv).

§.§ Distribution-free tests based on center-outward ranks and signs

[Hallin et al., 2020] propose, for multiple-output regression models with unspecified noise distribution $\rm P\in \mathcal P ^d$,  fully distribution-free yet, for adequate choice of   scores, parametrically efficient center-outward rank tests of the null hypothesis of no-treatment effect based on the empirical center-outward distribution functions $\Fb\n_\pm$ (hence, the center-outward ranks and signs).

The particular case of  two-sample location is treated by [Deb et al., 2021] who also consider tests based on the empirical Monge-Kantorovich vector ranks   ${\bf F}\n_{\text{\rm MK}}$ associated with various reference distributions.

§.§.§ Score functions

In line with the classical theory developed, e.g., by [Hájek & Šidák, 1967], rank-based statistics, irrespective of the reference distribution, involve score functions or scores.  Depending on the context, a   score function is  a   mapping $\Jb$ from the unit ball ${\mathbb S}_d$ or the unit cube $[0,1]^d$

to ${\mathbb R}^d$ satisfying some mild regularity assumptions    (continuity, square integrability, etc.: see, e.g. [Hallin et al., 2020], Assumption 3.1). The only score functions we are considering here are the Wilcoxon, the spherical  van der Waerden, and the marginal van der Waerden  score functions

$$\Jb_{\text{\rm W}} (\ub)\!:=\!\ub,\quad\!\! \Jb^\pm_{\text{\rm vdW}}(\ub)\!:=\!\sqrt{F_{\chi^2_d}^{-1}}(\Vert\ub\Vert)\frac{\ub}{\Vert\ub\Vert},\text{ and } \Jb^{\scriptscriptstyle\square}_{\text{\rm vdW}}(\ub)\!:=\!\left(\Phi^{-1}(u_1),\ldots,\Phi^{-1}(u_d)
\right),

where $F_{\chi^2_d}$ and $\Phi$  stand for the (univariate) chi-square ($d$ degrees of freedom) and standard normal distribution functions, respectively.

§.§.§ Test statistics

For simplicity, our investigation here
is limited to that particular case of   two-sample location models, where   the $n$ observations   are i.i.d. under the null and  consist of two samples, $\Zb\n_1,\ldots,\Zb\n_{n_1}$ and $\Zb\n_{n_1+1},\ldots,\Zb\n_{n_1+n_2}$, with $n_1+n_2=n$. The classical procedure for this problem is Hotelling's test, based on  a quadratic statistic of the form

$$\left(T _{\text{\rm Hot}}\n\right)^2= {\Deltab}_{\text{\rm Hot}}^{(n)\prime}\left({\boldsymbol\Sigma}\n_{\text{\rm Hot}}\right)^{-1} {\Deltab}\n_{\text{\rm Hot}}

where ${\boldsymbol\Sigma}\n_{\text{\rm Hot}}$ is the estimated (under the null) covariance matrix of

\[
{\Deltab}\n_{\text{\rm Hot}}:=\frac{1}{n_1}\sum_{i=1}^{n_1}\Zb\n_i - \frac{1}{n_2}\sum_{i=n_1+1}^{n}\Zb\n_i .
\vspace{-1mm}  \]

The Hotelling test is parametrically efficient under Gaussian assumptions; it remains asymptotically valid,[Asymptotically valid, here, means pointwise (with respect to the actual density of the observations) asymptotically correct nominal probability  levels, not uniformly  asymptotically correct nominal probability  levels.] however, under mild moment assumptions and therefore qualifies as a pseudo-Gaussian procedure.

For score functions $\Jb$, the center-outward rank-based test statistics in Section 5.3.1 of [Hallin et al., 2020]

are of the form

\begin{equation}%\[
\label{2samplestatpm}
\left(\tenq{T}^{(n)}_{\Jb\pm}\right)^2=\tenq{\Deltab}_{\Jb\pm}^{(n)\prime}\left({\boldsymbol\Sigma}_{\scriptscriptstyle\tenq{\scriptstyle\Deltab}_{\Jb\pm}}\right)^{-1}\tenq{\Deltab}\n_{\Jb\pm}%\frac{nd}{n_1n_2\int_0^1 J^2(u)du }
%\frac{nd}{n_1n_2\!\int_0^1 \! J^2(u){\rm d}u}\left(
%\sum_{i=1}^{n_1}  J\Big(\frac{R\n_{i;{\pm}}}{n_R+1}\Big) {{\bf S}\n_{i;{\pm}}}
%-\frac{n_1}{n}\sum_{i=1}^n J\Big(\frac{R\n_{i;{\pm}}}{n_R+1}\Big){\bf S}\n_{i;{\pm}}\right)
\end{equation}

where ${\boldsymbol\Sigma}_{\scriptscriptstyle\tenq{\scriptstyle\Deltab}_\Jb}$ is the exact or asymptotic covariance of

\begin{equation}
\label{2sampledelta}
\tenq{\Deltab}\n_{\Jb\pm}:=\frac{1}{n_1}\sum_{i=1}^{n_1} \Jb(\Fb\n_\pm(\Zb\n_i))
-\frac{1}{n_2}\sum_{i=n_1+1}^{n} \Jb(\Fb\n_\pm(\Zb\n_i)) .
\vspace{-1mm}
\end{equation}

Since the quadratic form (<ref>) is invariant under affine  transformations of $\tenq{\Deltab}\n_{\Jb\pm}$ and since the sum $\sum _{i=1}^n\Jb(\Fb\n_\pm(\Zb\n_i))$ is a deterministic constant that only depends on $\Jb$ and the grid used in the definition of $\Fb\n_\pm$, the same test statistic can be based on

$$\tenq{\Deltab}^{(n)}_\Jb= \frac{1}{n_1} \sum_{i=1}^{n_1} \Jb(\Fb\n_\pm(\Zb\n_i)) -\frac{1}{n}\sum_{i=1}^n  \Jb(\Fb\n_\pm(\Zb\n_i)),
\vspace{-1mm} $$

yielding the test statistic described in Section 5.3.1 of [Hallin et al., 2020] which, in the particular case of  Wilcoxon and van der Waerden scores $\Jb_{\text{\rm W}} $ and  $\Jb^\pm_{\text{\rm vdW}} $,  we denote as $\left(\tenq{T}^{(n)}_{\text{\rm W}\pm}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm vdW}\pm}\right)^2$, respectively.

For the same testing problem,  [Deb et al., 2021] consider statistics of the form (<ref>), but also

(a) based on the empirical Monge-Kantorovich vector ranks ${\bf F}\n_{\scriptscriptstyle\square}
$ associated  with the cubic uniform  reference $\rm U ={\rm U}_{[0,1]^d}$, statistics $\tenq{T}^{(n)}_{\Jb{\scriptscriptstyle\square}}$ and $\tenq{\Deltab}\n_{\Jb{\scriptscriptstyle\square}}$ of the same form as $\tenq{T}^{(n)}_{\Jb\pm}$ and $\tenq{\Deltab}\n_{\Jb\pm}$ in (<ref>) and (<ref>) but  with $\Jb(\Fb\n_{\scriptscriptstyle\square}(\Zb\n_i)) $ substituting $\Jb(\Fb\n_\pm(\Zb\n_i))$;

denote by $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}}\right)^2$  the particular cases of  the Wilcoxon and cubic van der Waerden statistics obtained for the scores $\Jb_{\text{\rm W}} $ and  $\Jb^{\scriptscriptstyle\square}_{\text{\rm vdW}}$, respectively;

(b) based on the empirical Monge-Kantorovich vector ranks ${\bf F}\n_{\pm{\mathcal N}}
$ and ${\bf F}\n_{\scriptscriptstyle\square{\mathcal N}}$ associated with the spherical Gaussian reference ${\mathcal N}({\mathbf 0}, {\mathbf I}_d)$ considered as spherical or as a product of independent uniforms,  statistics $\tenq{T}^{(n)}_{\Jb{\pm}{\mathcal N}}$ and $\tenq{T}^{(n)}_{\Jb{\scriptscriptstyle\square}{\mathcal N}}$   of the same form as $\tenq{T}^{(n)}_{\Jb\pm}$ and $\tenq{\Deltab}\n_{\Jb\pm}$ in (<ref>) and (<ref>) but  with $\Jb(\Fb\n_{\pm{\mathcal N}}(\Zb\n_i)) $ and $\Jb(\Fb\n_{\scriptscriptstyle\square{\mathcal N}}(\Zb\n_i)) $,  respectively,  substituting $\Jb(\Fb\n_\pm(\Zb\n_i))$; this, for Wilcoxon scores $\Jb(\ub)=\ub\vspace{-1mm}$, yields the van der Waerden statistics $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}{\mathcal N}}\right)^2$.

Although they are based on Wilcoxon (identity) scores, the terminology “van der Waerden statistic” for $\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}\vspace{1mm}$ and $\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}{\mathcal N}}$ seems more appropriate than the terminology “Wilcoxon statistic”  used by [Deb et al., 2021], and is   in line with the traditional terminology of rank-based inference. Both $\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}}\vspace{1mm}$ and $\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}$ indeed result from a transport from the sample values to a grid of Gaussian quantiles of the form ($\mathfrak G$iii). For   $\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}}$, the transport is $\Jb\circ \Fb\n_\pm$, which, as a rule, is not an optimal one (not the gradient of a convex function) while, for  $\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}$, the transport is the optimal one ${\bf F}\n_{\scriptscriptstyle{\square}{\cal N}}$;  the difference between $\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}}\vspace{1mm}$ and $\tenq{T}^{(n)}_{\text{\rm vdW}\pm{\mathcal N}}$    thus essentially consists in the way the transport to the spherical Gaussian grid is performed.  A similar remark can be made for $\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}}\vspace{1mm}$ and $\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}{\mathcal N}}$.

§ FINITE-SAMPLE PERFORMANCE: TWO-SAMPLE LOCATION SIMULATIONS

It clearly appears that  choices are to be made before performing a rank test based on the concepts  of multivariate ranks developed on the previous sections: center-outward ranks?  vector ranks? which ones? with which scores? The analysis of asymptotic performance does not help much, as the same local powers can be achieved by several  alternatives. The objective of this paper is to determine whether finite-sample performance can help us with these choices. We restrict ourselves to the two-sample location problem, Wilcoxon and van der Waerden scores, but the conclusions are quite likely to hold for the general case of multiple-output linear models considered in [Hallin et al., 2020].

Before explaining how  simulations were conducted, let us provide some details on the way the grids described in Section <ref> were constructed. Recall that the aim of these grids is to provide a discrete approximation of the chosen continuous reference distribution.

§.§ Halton sequences on the cube and the sphere (($\mathfrak G$ii) and ($\mathfrak G$iv)  grids)

The grid constructions ($\mathfrak G$ii) and ($\mathfrak G$iv) involve Halton sequences on the hypercube $[0,1]^d$.  Halton sequences are pseudo-random numbers with low discrepancy of classical use in methods such as Monte Carlo simulations.  We used the implementation available in the package
by McDonald & McDonald, 2020.
The grid construction in ($\mathfrak G$i), hence also in ($\mathfrak G$iii), requires a $n_S$-point “Halton sequence” over the hypersphere ${\mathcal S}^{d-1}$. To obtain such a grid, we first generate a $n_S$-point Halton sequence over $[0,1]^{d-1}$ then componentwise perform the standard normal quantile transformation $u_j\mapsto z_j:=\Phi^{-1}(u_j)$. This yields $n_S$ points $\zb_1,\ldots,\zb_{n_S}$, with

$$\zb_j:= (\Phi^{-1}(u_{j1}),\ldots, \Phi^{-1}(u_{jd})).$$

The   resulting unit vectors $\zb_j/\Vert\zb_j \Vert$, $j=1,\ldots,n_S$ constitute the desired sequence over ${\mathcal S}^{d-1}$.
[The justification is the fact that if the distribution of $\Zb$ is a product of independent univariate standard normal marginals, then $\Zb$ is spherical Gaussian ${\mathcal N}({\boldsymbol 0}, {\bf I})$ and hence $\Zb/\Vert\Zb\Vert$ is uniform over ${\mathcal S}^{d-1}$.

§.§ Factorization of $n$ (($\mathfrak G$i) and ($\mathfrak G$iii) grids)

As for the grid constructions ($\mathfrak G$i) and ($\mathfrak G$iii), they require a factorization of $n$ into $n_Rn_S +~\!n_0$ with $n_0 <\min(n_R, n_S)$. Intuition suggests choosing $n_R$ and $n_S$ of  order $n^{1/d}$ and $n^{(d-1)/d}$, respectively. This, however, is of little help for finite $n$. Since the grid is supposed to provide an approximation of the spherical uniform, we rather proceed by minimizing the Wasserstein distance to the spherical uniform as proposed in [Mordant, 2021].

More precisely, considering the grid with $n_R$ radial points described in ($\mathfrak G$i), denote by  ${\rm G}\n_{n_R}$ the   discrete measure placing a probability mass $1/n$ on each of the $n$ gridpoints    except for the origin which receives probability mass  $n_0/n$.  We adopt here the strategy suggested in [Mordant, 2021]  by selecting the grid with $n_R^*$ radial points, where

\begin{equation}\label{minfact}
n_R^*:=\argmin_{ 1\le n_R\le n} \ W_2( {\rm G}\n_{n_R},\mathrm{U}_d)\vspace{-1mm}
\end{equation}

($W_2$, as usual,  stands for the Wasserstein distance of order two). For $d\geq 3$, that distance $W_2( {\rm G}\n_{n_R},\mathrm{U}_d)$ does not only depend on $n_R$ (hence on $n_S$) but also on the $n_S$ points chosen (as explained in Section <ref>)   on the hypersphere ${\mathcal S}^{d-1}$.   The minimization in (<ref>) is feasible, as $n$, $n_R$, $n_S$, and $n_0$ all are integers. A similar strategy is adopted for the construction of the spherical Gaussian grids ($\mathfrak G$iii).

Table <ref> provides, for dimensions $d=2$ and $d=5$, various sample sizes,

and  reference distributions  the spherical uniform (($\mathfrak G$i) grids) and the spherical Gaussian (($\mathfrak G$iii) grids),  the “optimal values"  obtained via (<ref>) for $n_R$, $n_S$, and $n_0$. These values are in line with the intuition that the “optimal” $n_R$ behaves like $n^{1/d}$ while the role of distances to the center rapidly decreases as the dimension $d$ increases.

Reference distribution       $d$       $\quad n=50\quad $     $\quad n=100\quad $      $\quad n=200\quad $      $\quad n=300\quad $      $\quad n=400\quad $

$n_R=4$     $n_R=6$     $n_R=9$    $n_R=11$    $n_R=12$

$\mathrm{U}_d$, ($\mathfrak G$i) grid    2     $n_S=12$   $n_S=16$     $n_S=22$    $n_S=27$    $n_S=33$

$n_0=2$   $n_0=4$     $n_0=2$    $n_0=3$    $n_0=4$

$n_R=4$    $n_R=7$     $n_R=11$     $n_R=14$     $n_R=18$

${\mathcal N}({\mathbf 0}, {\mathbf I}_d)$,  ($\mathfrak G$iii) grid     2     $n_S=12$    $n_S=14$     $n_S=18$    $n_S=21$    $n_S=22$

$n_0=2$    $n_0=2$     $n_0=2$    $n_0=6$    $n_0=4$

$n^{1/d}$          $n^{1/2}\!=7.071$    $n^{1/2}\!=10.000$
$n^{1/2}\!=14.142$    $n^{1/2}\!=17.321$    $n^{1/2}\!=20.000$

$n_R=2$     $n_R=2$      $n_R=2$      $n_R=3$      $n_R=3$

$\mathrm{U}_d$, ($\mathfrak G$i) grid      5      $n_S=25$     $n_S=50$      $n_S=100$      $n_S=100$      $n_S=133$

$n_0=0$     $n_0=0$      $n_0=0$      $n_0=0$      $n_0=1$

$n_R=1$      $n_R=1$       $n_R=1$       $n_R=2$       $n_R=2$

${\mathcal N}({\mathbf 0}, {\mathbf I}_d)$, ($\mathfrak G$iii) grid    5       $n_S=50$      $n_S=100$       $n_S=200$       $n_S=150$       $n_S=200$

$n_0=0$      $n_0=0$       $n_0=0$       $n_0=0$       $n_0=0$

$n^{1/d}$          $n^{1/5}\!=2.187$    $n^{1/5}\!=2.512$
$n^{1/5}\!=2.885$    $n^{1/5}\!=3.129$    $n^{1/5}\!=3.314$

Optimal (in the sense of (<ref>)) values of  $n_R$, $n_S$, and $n_0$ as  functions of the sample size $n$, the dimension $d$, and the reference distributions (($\mathfrak G$i) or ($\mathfrak G$iii) grids).

§.§ Simulations

Based on the  grids obtained along the lines described in Sections <ref> and <ref>,   the distribution-free critical values of the various rank tests under study were computed from 40 000 replications. Throughout, we chose $n_1=n_2=n/2$. The optimal maps between the sample and the grids were  obtained via an exact solver relying on the so-called Hungarian method that is implemented in the -package  by Hornik, 2021. We now  turn to the empirical evaluation of the performance of the various rank-based Wilcoxon and van der Waerden  tests for the two-sample location problem.

The objective of our simulations is, essentially,   empirical answers to the following two questions:

(a) should we use spherical grids (($\mathfrak G$i) or ($\mathfrak G$iii)) or cubic (($\mathfrak G$ii) or ($\mathfrak G$iv)) ones?

(b) should we, in line with the Hájek tradition,  privilege transports to the uniform combined  with scores  (as in $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}}\right)^2$), or, as recommended by [Deb et al., 2021], should we rather consider “direct transports”  to the “scored distribution,” that is, choose as reference distribution  the push-forward of the uniform by the score as in $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}{\mathcal N}}\right)^2$?

§ WILCOXON-TYPE  TESTS

The Wilcoxon tests are based on the identity score function $\Jb (\ub)=\ub$ and uniform (either spherical or cubic) reference distributions, yielding (see Section <ref>) the test statistics $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2$.

§.§ The bivariate case

In this section, we evaluate the performance of the  Wilcoxon tests based on $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2\vspace{1mm}$ for samples of size $n_1=n_2=n/2$ with $n=100$, 200, and 400. The first sample is drawn from a centered distribution, the second one  from a shifted version with shift $(\eta, \eta)\pr$, $\eta>0$ of the same.  The number of replications is $N=500$.

§.§.§ Spherical Gaussian samples

The first sample is drawn from $\mathcal{N} ( (0, 0)\pr ,  {\bf I}_2)$, the second one from $\mathcal{N} ( (\eta, \eta)^\top ,  {\bf I}_2)$ with $\eta>0$.   Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in Figure <ref>. All three tests display, essentially, the same performance: although Wilcoxon, in principle, is strictly less powerful than Hotelling (which in this case is finite-sample optimal), no significant loss of efficiency is detected.

§.§.§ Spherical Student samples

The first sample is drawn from a centered spherical Student with 2.1 degrees of freedom $ t_{2.1} ( (0,0)\pr, {\bf I}_2)$,   the second one from the shifted version $ t_{2.1} ((\eta, \eta)\pr , {\bf I}_2)$ of the same distribution.   Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in Figure <ref>. The Wilcoxon tests substantially outperform Hotelling, an advantage that does not disappear with increasing $n$. Although the underlying distribution is spherical,   very slight superiority of $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2\vspace{1mm}$ over $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$.

§.§.§ “Banana-shaped”  samples

The first sample is drawn from a centered “banana-shaped”   mixture

\begin{align*}
\nonumber
0.3 %(1-2p)%
\, \normal_2\left(
\begin{pmatrix} \phantom{-}0\phantom{.7} \\ -0.7 \end{pmatrix},
\begin{pmatrix}
0.35^2&0 \\
0& 0.35^2
\end{pmatrix}
\right)&+
0.35 %p
\, \normal_2\left(
\begin{pmatrix} -0.9 \\ \phantom{-}0.3 \end{pmatrix},
\begin{pmatrix}
\phantom{-}0.358&-0.55 \\
-0.55\phantom{0}& \phantom{-}1.02
\end{pmatrix}
\right)\\
\label{eq:banana}
0.35 %p
\, \normal_2\left(
\begin{pmatrix} \phantom{-}0.9 \\ \phantom{-}0.3 \end{pmatrix},
\begin{pmatrix}
\phantom{-}0.358& \phantom{-}0.55 \\
\phantom{-}0.55\phantom{0}& \phantom{-}1.02
\end{pmatrix}
\right).
\end{align*}

of three Gaussian components.
The second sample is drawn from a shifted version (shift $(\eta, \eta)\pr$, $\eta>0$) of the same mixture.   Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in Figure <ref>. The conclusions are the same  as in the previous case, except that the (slight) advantage now belongs to $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$, despite the fact that the actual distribution is highly nonspherical.

§.§.§ Samples   with independent Cauchy marginals

The first sample is drawn from a product of two independent Cauchy,   the second one from the shifted version of the same distribution.   Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>. With a rejection probability uniformly less than the nominal $5\%$ level, Hotelling, as expected, performs  miserably. In this independent component situation,  $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2\vspace{-.8mm}$ does outperform $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$.

§.§.§ Nonspherical Gaussian samples

The first sample is drawn from a $ \mathcal{N}( (0,0)\pr, \Sigmab)$ distribution,    the second sample is drawn from a $\mathcal{N}( (\eta,\eta)\pr, \Sigmab)$ one;    $\text{vech}(\Sigmab) = (1,0.8,1)\pr$.   Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in Figure <ref>.  The results are essentially the same as in the spherical case (Section <ref>). Note the loss of power  in the three tests under study, due to the non-specification of the population covariance matrix; that loss, however, is uniform over the three tests.

§.§.§ Spherical Cauchy samples

The first sample is drawn from a centered spherical Student with one degree of freedom $ t_{1} ( (0,0)\pr, {\bf I}_2)$ (spherical Cauchy),   the second one from the shifted version $ t_{1} ((\eta, \eta)\pr , {\bf I}_2)$ of the same distribution.   Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in Figure <ref>.  The performance of Hotelling, again, is a disaster; although the actual distribution is spherical, $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2\vspace{-.8mm}$ still outperforms $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$.

§.§ Wilcoxon-type statistics in dimension $d=5$

We essentially adopted the same simulation settings as before, with $n_1=n_2=n/2$. A sample size of $n=100$ in dimension $d=5$ is very small, though, and we considered sample   sizes $n=200$, 400, and 800.

§.§.§ Spherical Gaussian samples

Here, the first sample is drawn from the $ \mathcal{N} ( {\bf 0}, {\bf I}_5)$ distribution, the second one   from  the $ \mathcal{N} ( {\eta}{\mathbf 1} , {\bf I}_5)$ distribution, where ${\mathbf 1}$ denotes a 5-variate vector of ones.   Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>.  In the “small sample” case $(n=200)$, the optimality of  Hotelling over Wilcoxon is perceptible (more so than in dimension $d=2$); this superiority, however, fades away with growing $n$: again, under Gaussian assumptions, abandoning the parametrically optimal Hotelling test in favor of the rank-based Wilcoxon one has no visible cost in terms of power.

§.§.§ Spherical Student samples

The first sample is drawn from the centered spherical Student with 2.1 degrees of freedom $ t_{2.1} ({\boldsymbol 0}, {\bf I}_5)$, the second one   from  the shifted $ t_{2.1} ({\eta}{\mathbf 1} , {\bf I}_5)$ distribution, where ${\mathbf 1}$ denotes a 5-variate vector of ones.   Rejection frequencies over $N=500$ replications

shown (as  functions of $\eta$) in  Figure <ref>.  The conclusions are quite similar to those

in dimension $d=2$: although the moments of order two still are finite, the power of Hotelling deteriorates with respect to the Gaussian case. Despite the spherical nature of the actual distribution, slight advantage of $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2\vspace{-.8mm}$ over $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$.

§.§.§ Nonspherical Gaussian samples

The first sample is drawn from the $ \mathcal{N} ( {\bf 0}, \Sigmab)$ distribution, the second one   from  the $ \mathcal{N} ( {\eta}{\mathbf 1} , \Sigmab)$ distribution, where ${\mathbf 1}$ denotes a 5-variate vector of ones and  $\Sigmab$ is a correlation matrix with all off-diagonal entries equal to 0.5.    Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>.  Here again, the slight advantage of Hotelling over Wilcoxon very rapidly fades away with growing $n$, and the three tests yield very similar performances; in particular, no significant difference can be detected between   $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2\vspace{-.8mm}$.

§.§.§ Samples with independent Cauchy marginals

The first sample is drawn from a product of five independent Cauchy distributions,   the second one from the shifted version of the same.   Rejection frequencies over $500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>.  The performance of Hotelling, as in dimension $d=2$, is terrible. The advantage (which is in line with the independent component nature of the distribution) of $\left(\tenq{T}^{(n)}_{\text{\rm W}{\scriptscriptstyle\square}}\right)^2\vspace{-.8mm}$ over  $\left(\tenq{T}^{(n)}_{\text{\rm W}{\pm}}\right)^2\vspace{1mm}$ is even more significant than in dimension $d=2$.

§ VAN DER WAERDEN-TYPE  TESTS

Below we are considering four distinct tests of the van der Waerden type, based (see Section <ref>) on $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}}\right)^2\vspace{-1mm}$ (spherical uniform reference density, ($\mathfrak G$i) grid), $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}}\right)^2$ (cubic uniform reference density, ($\mathfrak G$ii) grid), $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}\right)^2$ (spherical Gaussian reference density, spherical grid ($\mathfrak G$iii)), and $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}{\mathcal N}}\right)^2\vspace{.6mm}$  (spherical Gaussian reference density, cubic grid ($\mathfrak G$iv)).

§.§ Bivariate case

§.§.§ Spherical Gaussian samples

Same Gaussian samples as in Section <ref>.
Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>. Indistinctiveness between the performances of Hotelling and the various rank-based tests is even more pronounced than for the Wilcoxon tests: definitely, performing rank-based van der Waerden tests does not imply any loss of efficiency in the Gaussian case.

§.§.§ Spherical Student samples

Same Student samples as in Section <ref>.
Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>.
As in dimension $d=2$, powers are dropping (compared with the Gaussian case). The power of Hotelling, however,  deteriorates much more than that   of the various versions of van der Waerden  tests. The latter all yield very similar performance.

§.§.§ “Banana-shaped” samples

Same “banana-shaped” mixtures as in Section <ref>.
Rejection frequencies
over $500$
replications are
shown (as  functions of $\eta$) in  Figure <ref>.  The conclusions are quite
similar as in the previous case: the empirical power curves of the various van der Waerden tests are essentially indistinguishable, while significantly outperforming the Hotelling ones.

§.§.§ Samples   with independent Cauchy marginals

Same  samples with independent Cauchy marginals as in Section <ref>.
Rejection frequencies over $N=500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>.
Again, all powers are much less than in the Gaussian case, but Hotelling is totally inefficient.

§.§.§ Nonspherical Gaussian samples

Same  correlated Gaussian samples  as in Section <ref>.
Rejection frequencies over $N=~\!500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>.
The non-specification of the covariance matrix apparently has no impact on the comparative performance of Hotelling and its rank-based van der Waerden competitors, which all coincide.

§.§.§ Spherical Cauchy samples

Same  spherical Cauchy samples   as in Section <ref>.
Rejection frequencies over $500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>. All tests perform similarly except, of course, for Hotelling, which fails completely.

§.§ van der Waerden-type statistics in dimension $d=5$

§.§.§ Spherical Gaussian samples

Same  spherical Gaussian samples   as in Section <ref>.
Rejection frequencies over $N=~\!500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>. A “small sample” superiority of  Hotelling for $n=200$ rapidly disappears as $n$ increases; all van der Waerden tests yield the same performance.

§.§.§ Spherical Student samples

Same  spherical Student samples   as in Section <ref>.
Rejection frequencies over $500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>. All van der Waerden tests roughly yield the same performance, with a slight adavantage in favor of $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}\right)^2$ in “small samples.”

§.§.§ Nonspherical Gaussian samples

Same  spherical Gaussian samples   as in Section <ref>.
Rejection frequencies over $N=~\!500$ replications are
shown (as  functions of $\eta$) in  Figure <ref>. The slight advantage for $n=200$ of Hotelling under spherical Gaussian has almost disappeared. All van der Waerden tests yield similar performance.

§.§.§ Samples with independent Cauchy marginals

Same  Cauchy samples   as in Section <ref>.
Rejection frequencies over $N=~\!500$
replications are
shown (as  functions of $\eta$) in  Figure <ref>. The conclusions drawn for $d=2$ still hold, with a very slight superiority of the “direct transportation” test $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}{\mathcal N}}\right)^2$ over  the Gaussian score ones $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\pm}}\right)^2$ and $\left(\tenq{T}^{(n)}_{\text{\rm vdW}{\scriptscriptstyle\square}}\right)^2$.

§ CONCLUSIONS

While confirming the advantages and excellent performance of rank tests over their daily practice pseudo-Gaussian counterparts, the simulations of the previous sections provide  empirical answers to several  questions of great practical importance.

The  choice of the grid (whether spherical ($\mathfrak G$i), cubic ($\mathfrak G$ii), or Gaussian (($\mathfrak G$iii) and ($\mathfrak G$iv)) seems to have relatively  little impact on the performance of the corresponding Wilcoxon tests (the only significant case being the Cauchy one), and no impact at all on the performance of van der Waerden tests. In particular, there is no evidence that Wilcoxon tests based on spherical grids ($\mathfrak G$i) are preferable under spherical distributions while Wilcoxon tests based on  cubic grids ($\mathfrak G$ii)  are preferable under distributions with  independent components: see, e.g., the Cauchy case   (Sections <ref> and <ref>).

[Carlier et al., 2016]
Carlier, G., Chernozhukov, V., & Galichon, A. (2016).
Vector quantile regression: an optimal transport approach.
Ann. Statist., 44, 1165–92.

[Chernozhukov et al., 2017]
Chernozhukov, V., Galichon, A., Hallin, M., & Henry, M. (2017).
Monge-Kantorovich depth, quantiles, ranks and signs.
Ann. Statist., 45, 223–256.

[Deb et al., 2021]
Deb, N., Bhattacharya, B. B., & Sen, B. (2021).
Efficiency lower bounds for distribution-free Hotelling-type
two-sample tests based on optimal transport.

[Deb et al., 2020]
Deb, N., Ghosal, P., & Sen, B. (2020).
Measuring association on topological spaces using kernels and
geometric graphs.

[Deb & Sen, 2021]
Deb, N., & Sen, B. (2021).
Multivariate rank-based distribution-free nonparametric testing using
measure transportation.
Journal of the American Statistical Association, 0(0), 1–16.

[del Barrio et al., 2020]
del Barrio, E., González-Sanz, A., & Hallin, M. (2020).
A note on the regularity of optimal-transport-based center-outward
distribution and quantile functions.
J. Multivariate Anal., 180, 104671, 13.

[Faugeras & Rüschendorf, 2017]
Faugeras, O., & Rüschendorf, L. (2017).
Markov morphisms: a combined copula and mass transportation approach
to multivariate quantiles.
Mathematica Applicanda, 45, 21–63.

[Figalli, 2018]
Figalli, A. (2018).
On the continuity of center-outward distribution and quantile
Nonlinear Anal., 177 , 413–21.

[Ghosal & Sen, 2019]
Ghosal, P., & Sen, B. (2019).
Multivariate ranks and quantiles using optimal transportation and
applications to goodness-of-fit testing.

[Hájek & Šidák, 1967]
Hájek, J., & Šidák, Z. (1967).
Theory of Rank Tests.
Academic Press, New York.

[Hallin, 2017]
Hallin, M. (2017).
On distribution and quantile functions, ranks, and signs in
$\mathbb{R}^d$: a measure transportation approach.

[Hallin, 2022]
Hallin, M. (2022).
Measure transportation and statistical decision theory.
Annual Review of Statistics and its Applications, in
press, 9.

[Hallin et al., 2021]
Hallin, M., Del Barrio, E., Cuesta-Albertos, J., & 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.

[Hallin et al., 2020]
Hallin, M., Hlubinka, D., & Hudecová, Š. (2020a).
Fully distribution-free center-outward rank tests for multiple-output
regression and MANOVA.

[Hallin et al., 2020]
Hallin, M., La Vecchia, D., & Liu, H. (2020b).
Center-outward R-estimation for semiparametric VARMA models.
Journal of the American Statistical Association, in

[Hallin et al., 2020]
Hallin, M., La Vecchia, D., & Liu, H. (2020c).
Rank-based testing for semiparametric var models: a measure
transportation approach.

[Hallin & Paindaveine, 2002]
Hallin, M., & Paindaveine, D. (2002).
Optimal tests for multivariate location based on interdirections and
pseudo-Mahalanobis ranks.
Ann. Statist., 30, 1103–33.

[Hampel, 1968]
Hampel, F. R. (1968).
Contributions to the Theory of Robust Estimation.
Ph.D. thesis, University of California, Berkeley.

[Hornik, 2021]
Hornik, K. (2021).
clue: Cluster ensembles.
R package version 0.3-60.

[Huber, 1964]
Huber, P. (1964).
Robust estimation of a location parameter.
Annals of Mathematical Statistics, 35, 73–101.

[Masud & Aeron, 2021]
Masud, S. B., & Aeron, S. (2021).
Soft and subspace robust multivariate rank tests based on entropy
regularized optimal transport.

[McCann, 1995]
McCann, R. J. (1995).
Existence and uniqueness of monotone measure-preserving maps.
Duke Math. J., 80, 309–324.

[McDonald & McDonald, 2020]
McDonald, T., & McDonald, A. (2020).
SDraw: Spatially balanced samples of spatial objects.
R package version 2.1.13.

[Mordant, 2021]
Mordant, G. (2021).
Transporting Probability Measures: some contributions to
statistical inference.
Ph.D. thesis, Université catholique de Louvain.

[Ronchetti, 2006]
Ronchetti, E. (2006).
The historical development of robust statistics.
In A. Rossman, & B. Chance (Eds.) ICOTS-7 Proceedings.

[Shi et al., 2021]
Shi, H., Drton, M., & Han, F. (2021a).
Distribution-free consistent independence tests via center-outward
ranks and signs.
Journal of the American Statistical Association, in

[Shi et al., 2021]
Shi, H., Hallin, M., Drton, M., & Han, F. (2021b).
On universally consistent and fully distribution-free rank tests of
vector independence.
Annals of Statistics, to appear.

[Stigler, 1973]
Stigler, S. M. (1973).
Simon Newcomb, Percy Daniell, and the history of robust
estimation 1885-1920.
Journal of the American Statistical Association, 68,

[Tukey, 1960]
Tukey, J. W. (1960).
A survey of sampling from contaminated distributions.
In I. Olkin (Ed.) Contributions to Probability and
Statistics, (pp. 448–485). Palo Alto: Stanford University Press.


# Photon-number-resolving segmented avalanche-photodiode detectors

Rajveer Nehra  Chun-Hung Chang  Department of Physics, University of Virginia,
382 McCormick Rd, Charlottesville, VA 22903  Andreas Beling  Department of
Electrical and Computer Engineering, University of Virginia, 351 McCormick Rd,
Charlottesville, VA 22903  Olivier Pfister  Department of Physics, University
of Virginia, 382 McCormick Rd, Charlottesville, VA 22903

(August 29, 2017)

###### Abstract

We investigate the feasibility and performance of photon-number-resolved
photodetection employing avalanche photodiodes (APDs) with low dark counts.
The main idea is to split $\bf n$ photons over $\bf m$ modes such that every
mode has no more than one photon, which is detected alongside propagation by
an APD. We characterize performance by evaluating the purities of positive-
operator-valued measurements (POVMs), in terms of APD number and photon loss.

##  I Introduction

Quantum measurements are essential to quantum technology. Photon-Number-
Resolution (PNR) has become a need in various fields including linear optical
quantum computing Knill  _et al._ (2001), quantum metrology and
sensingZwinkels  _et al._ (2010), quantum cryptography O’Brien  _et al._
(2009), quantum imagingG.  _et al._ (2010), and quantum communication
Sangouard  _et al._ (2007) and conditional state preparation Kok  _et al._
(2002).
A PNR detector produces a signal proportional to the number of incident
photons. Currently, PNR detectors can be realized with silicon
photomultipliers Ramilli  _et al._ (2010), superconducting nanowires Divochiy
_et al._ (2008); Natarajan  _et al._ (2013); Mattioli  _et al._ (2015), linear
mode avalanche photodiodes (APDs), quantum-dot field-effect transistors
Kardynał  _et al._ (2007); Shields  _et al._ (2000), and transtion-edge
sensors (TES) Rosenberg  _et al._ (2005); Lita  _et al._ (2008). Several
methods based on spatial- and time-multiplexing have been proposed for PNR
measurements using APDs Achilles  _et al._ (2003); Fitch  _et al._ (2003);
Piacentini  _et al._ (2015); Mattioli  _et al._ (2015).
In this paper we propose a way for PNR using segmented APDs with low dark
currents. In the figure shown below we have APDs on waveguide. The set up
whittles down n photons, one at a time.

\begin{overpic}[width=212.47617pt]{waveguide.png} \put(73.0,13.0){$\alpha$}
\put(48.0,13.0){$\alpha$} \put(23.0,13.0){$\alpha$} \put(7.0,9.0){$|n\rangle$}
\end{overpic} Figure 1: Segmented detector. Guided optics are used to detect
photons alongside propagation by APDs. The design goal is to keep all
undetected photons in the waveguide, for further detection.

Our segmented photodetector features multiple APDs fed by the same waveguide,
as shown in Fig. 1. Photons that are not absorbed in the first APD couple back
into the waveguide and will be absorbed in one of the following detectors. The
crucial advantage of this configuration is that nonideal quantum efficiency of
the APDs doesn’t amount to photon loss, unlike butt-coupled PNR detectors in
which temporally or spatially split photons impinge on APDs on the end of
their path Fitch  _et al._ (2003); Achilles  _et al._ (2003); Piacentini  _et
al._ (2015). Hence, this linear array is essentially a long detector divided
into $m$ short detectors, each with an individual read-out. We envision that
such a segmented photodetector will become feasible in large-scale integrated
photonic platforms using either monolithic or heterogeneous integration of
APDs on low-loss waveguides.

We will model the segmented detector using beam splitters as shown in Fig.2
where waveguide coupling with APDs can be modeled by the reflectivitys of the
beam splitters, and the transmissivity per APD by $\eta$.

Figure 2: Model of a PNR segmented photodetector.
$R_{j}+T_{j}=r_{j}^{2}+t_{j}^{2}=1,\forall j\in[1,m]$.

Before we address the model in earnest, we need to reiterate that $\eta$ won’t
include quantum efficiency of the APD. This can be seen by a detail model of
each APD coupling, pictured in Fig.3. The APD’s field transmissivity is
$\alpha$ (such that the APD’s quantum efficiency is $1-\alpha^{2}$). In fact,
$\alpha$ cannot be too large because APDs are not PNR detectors and we want
the quantum efficiency $1-\alpha^{2}\ll 1/n$, for $n$ incident photons, which
will ensure that photons are detected no more than one at a time.

\begin{overpic}[width=104.07117pt]{coupling.png}
\put(99.0,21.0){$tt^{\prime}+rr^{\prime}\sqrt{1-\alpha^{2}}$}
\put(60.0,-3.0){$tr^{\prime}-rt^{\prime}\sqrt{1-\alpha^{2}}$}
\put(3.0,33.0){$|n\rangle$} \put(28.0,16.0){$(r,t)$}
\put(54.0,16.0){$(r^{\prime},t^{\prime})$} \put(50.0,50.0){$\alpha$}
\end{overpic} Figure 3: Model for detection alongside propagation. If this is
the $\rm j^{th}$ APD, then $t_{j}=tt^{\prime}+rr^{\prime}\sqrt{1-\alpha^{2}}$,
where $t_{j}=T_{j}^{1/2}$ in Fig.2.

In a design reminiscent of interaction-free measurements proposed by Elitzur
and Vaidman Elitzur and Vaidman (1993), we require that the bottom output of
the exit beam splitter by nulled by destructive interference. The condition
can be achieved by choosing parameters $(r,t,r^{\prime},t^{\prime})$ of the
beam splitters, and absorption coefficient $\alpha$, such that

$tr^{\prime}-rt^{\prime}\sqrt{1-\alpha^{2}}=0.$ (1)

If this is the case, then the detection process truly takes place alongside
propagation and finite quantum efficiency — necessary here to attain PNR
performance — won’t contribute to photon loss.

In the next section, we recall basic properties of the quantum formalism
needed to evaluate PNR performance. We then return to the segmented detector
model, present our theoretical characterization results, and conclude.

##  II DETECTOR CHARACTERIZATION USING POSITIVE-OPERATOR-VALUED-MEASUREMENTS
(POVMs)

The most general measurements in quantum physics are known as POVMs van Enk
(2017a, b). Each measurement outcome k is given by an Hermitian operator, POVM
element $\hat{\Pi}_{k}$, with nonnegative eigenvalues such that the
probability of an outcome for a quantum state $\hat{\rho}$ is given by

$p_{k}=Tr(\rho\hat{\Pi}_{k}).$ (2)

These operators statisfy the completeness property, $\sum_{k}\hat{\Pi}_{k}=I$
which makes the sum of probabilities $\sum_{k}p_{k}=1$ for different outcomes
to be unity. Therefore, these measurements completely describe all possible
outcomes for any quantum measurement. For a phase insensitive detector the
POVM element for $k$ clicks is given as

$\Pi_{k}=\sum_{n=0}P(k|n)|n\rangle\langle n|,$ (3)

where $P(k|n)$ represents the conditional probability of getting $k$ clicks
given n-photon input state. The POVMs are more general from the projective
measurements in couple of ways. First, unlike the projective measurements the
POVMs are not orthogonal measurements. For outcomes $k$ and $k^{\prime}$, the
pairwise POVMs need not satisfy,

$\hat{\Pi}_{k}\hat{\Pi}_{k}^{\prime}\neq\delta_{kk^{\prime}}\Pi_{k}.$ (4)

Note that the orthogonal POVMs are essentially the projective measurements,
therefore we can define the purity of the POVM for outcome $k$ as

$Purity(\Pi_{k})=\frac{\big{[}Tr(\Pi_{k}^{2})\big{]}}{{\big{[}Tr(\Pi_{k})\big{]}}^{2}}$
(5)

The purity satisfies $0\leq Purity(\Pi_{k})\leq 1$, where a value of 1 denotes
a pure, i.e., projective POVM.

##  III Model

In order to study the PNR performance of the segmented detector we consider 2
cases: first, the beam splitters all have the same reflectivity and, second,
the reflectivity of the $\rm j^{th}$ beam splitter is given by
$R_{j}=\frac{1}{m-j+1}$ such that the last beam splitter in Fig.2 in balanced.
It can be shown that the field splitting generated by this is identical to
that of a symmetric beam splitter tree. We study each case in turn in the
lossless case before investigating the symmetric case in the presence of
losses.

###  III.1 SEGMENTED DETECTOR WITH IDENTICAL BEAM SPLITTERS

In Fig.2 the quantum input mode has annihilation operator $a_{1}$ and is in
the input Fock state $|n\rangle$, and the other $m-1$ input modes, of $a_{2}$,
$a_{3}$, …, $a_{m}$, are vacuum ones. We consider $m-1$ identical beam
splitters (T, R) and $\eta=1$ (no losses) for all modes. The input quantum
state is

$|n\rangle=\frac{{a^{\dagger}_{1}}^{n}}{\sqrt{n!}}|0\rangle^{\otimes^{m}}.$
(6)

In the Heisenberg picture, after m-1 beam splitters $a^{\dagger}_{1}$ evolves
to

$\displaystyle\underbrace{(U_{m-1}...U_{2}U_{1}}_{U}$ $\displaystyle
a^{\dagger}_{1}\underbrace{U^{\dagger}_{1}U^{\dagger}_{2}...U^{\dagger}_{m-1})}_{U\dagger}$
(7)
$\displaystyle=(t^{m-1}a^{\dagger}_{1}+rt^{m-2}a^{\dagger}_{2}+\dots+ra^{\dagger}_{m}).$

It is easy to see that the output quantum state is

$\displaystyle|\psi\rangle_{out}$
$\displaystyle=\frac{1}{\sqrt{n!}}(Ua^{\dagger}_{1}U^{\dagger})^{n}|0\rangle^{\otimes^{m}}$
(8)
$\displaystyle=\frac{1}{\sqrt{n!}}{(t^{m-1}a^{\dagger}_{1}+rt^{m-2}a^{\dagger}_{2}+\dots+ra^{\dagger}_{m}})^{n}|0\rangle^{\otimes^{m}}.$
(9)

Using multinomial expansion we get

$\displaystyle|\psi\rangle_{out}=$
$\displaystyle\underbrace{\sum_{n_{1}=0}^{n}\cdots\sum_{n_{m}=0}^{n}}_{\sum_{j=1}^{m}n_{j}=n}\left\\{\frac{\sqrt{n!}}{n_{1}!n_{2}!...n_{m}!}\right.$
(10)
$\displaystyle\left.[t^{\sum_{j=1}^{m}2(m-j)n_{m-j+1}}r^{(n-n_{m})}]\prod_{i=1}^{m}(a^{\dagger}_{i})^{n_{i}}\right\\}|0\rangle^{\otimes^{m}},$

where each $n_{i}$ can take any value from 0 to n. In order to have PNR each
$n_{i}$ must have at most one photon. The probability of splitting $n$ photons
over $m$ modes is

$P(\\{0,1\\}^{m},n,K)=n!t^{\sum_{j=1}^{m}2(m-j)n_{m-j+1}}r^{2(n-n_{m})},$ (11)

where $K=(n_{1},n_{2},\cdots,n_{m})^{T}$ is an $m$-dimensional vector. The
total probability of photon-number-resolution is

$\displaystyle P(\\{0,1\\}^{m},n)$
$\displaystyle=n!\underbrace{\sum_{n_{1}=0}^{1}\cdots\sum_{n_{m}=0}^{1}}_{\sum_{j=1}^{m}n_{j}=n}t^{2(\sum_{j=1}^{m}(m-j)n_{m-j+1})}r^{2(n-n_{m})}$
(12)
$\displaystyle=n!\underbrace{\sum_{n_{1}=0}^{1}\cdots\sum_{n_{m}=0}^{1}}_{\sum_{j=1}^{m}n_{j}=n}T^{\sum_{j=1}^{m}(m-j)n_{m-j+1}}R^{(n-n_{m})}.$
(13)

Note that our goal is to split the $n$ input photons over $m$ modes such that
each more have either zero or one photon, no matter what modes get the
photons. When all beam splitters are identical, the most constrained is the
first one as its reflectivity must be much less than $1/n$. However, this
approach isn’t optimal since the subsequent modes will gradually see fewer
photons and can therefore afford larger reflectivities without running the
risk of detecting more than one photon. Also, at the end of the device, the
last beam splitter should clearly be balanced since the constraint is
symmetric for both its output ports. Bearing all this in mind, a symmetrized
device appears to be the optimal choice. We investigate it next

###  III.2 SYMMETRIC SEGMENTED DETECTOR

####  III.2.1 Lossless case

Now we consider the case where the beam splitters have reflectivity
$R_{j}=\frac{1}{m-j+1}$, still in the lossless case ($\eta=1$). It is easy to
see that the probability of getting $k$ APD clicks (where each click may
result from one or several simultaneous photons) from $m$ modes, given $n$
input photons, $P_{m}(k|n)$ is

$P_{m}(k|n)={m\choose
k}\underbrace{\sum_{n_{1}=1}^{n}\cdots\sum_{n_{m}=1}^{n}}_{\sum_{j=1}^{k}n_{j}=n}\frac{n!}{\prod_{i=1}^{k}n_{i}!}X,$
(14)

where X is given by

$\displaystyle X=[(R_{1})^{n_{1}}$
$\displaystyle(T_{1})^{n-n_{1}}][(R_{2})^{n_{2}}(T_{2})^{n-n_{1}-n_{2}}]\cdots$
(15)
$\displaystyle[(R_{m-1})^{n_{m-1}}(T_{m-1})^{n-{n_{1}}-{n_{2}}\cdots-{n_{m-1}}}].$

Using $R_{j}=\frac{1}{(m-j+1)}$, and ${R_{j}+T_{j}=1}$ yields the
simplification

$X=\frac{1}{m^{n}}.$ (16)

Considering the number of possible configurations, we get

$P_{m}(k|n)=\frac{n!}{m^{n}}{m\choose
k}\underbrace{\sum_{n_{1}=1}^{n}\cdots\sum_{n_{m}=1}^{n}}_{\sum_{j=1}^{k}n_{j}=n}\frac{1}{\prod_{i=1}^{k}n_{i}!}.$
(17)

Figure 4: Conditional probabilities $P_{m}(k|n)$ versus $n$, for $m=100$ APDs
and $\eta=1$.

These conditional probabilities are plotted in Fig.4 for $m=100$.
Unsurprisingly, they increase for larger photon numbers as $m$ increases (see
supplemental material, for $m$ up to 2000). This translates directly into the
POVM purities, plotted in Fig.5.

Figure 5: POVM purity, versus click number $k$, for different numbers $m$ of
APDs.

As can be seen in that figure, reasonably good PNR performance can be reached
for $n\sim 10$ with $m\sim 1000$.

####  III.2.2 Lossy case

We now consider the effect of photon losses in each detection mode, i.e.,
$\eta<1$. Again, $\eta$ should not be misconstrued to be the APD quantum
efficiency, which plays no role in photon losses. (An instance of photon loss
would be if a photon went undetected and exited the waveguide, losing its
chance for further detection.) We assume that the parameter $\eta$ is
independent of the photon number. The probability to get zero clicks in one
mode is

$P_{1}(0|n,\eta)=(1-\eta)^{n}.$ (18)

Likewise, the probability to get one click in one mode is Scully and Lamb
(1969)

$P_{1}(1|n,\eta)=\sum_{k=1}^{n}{{n}\choose{k}}{\eta}^{k}{(1-\eta)^{n-k}}=1-(1-\eta)^{n}.$
(19)

It is important to note that in the sum over $k$ starts with 1 here because we
neglected dark counts. Therefore, Eq. (17) can be generalized to

$\displaystyle
P_{m}(k|n,\eta)=n!\left(\frac{1-\eta}{m}\right)^{n}{{m}\choose{k}}$
$\displaystyle\underbrace{\sum_{n_{1}=0}^{n}\cdots\sum_{n_{m}=0}^{n}}_{\sum_{j=1}^{k}n_{j}=n}\left\\{\frac{1}{\prod_{j=1}^{k}n_{j}!}\right.$
(20)
$\displaystyle\left.\prod_{l=n_{1}}^{n_{k}}\left[\left(\frac{1}{1-\eta}\right)^{l}-1\right]\right\\}.$

It is worth noting that the probability of getting zero clicks is still the
same as the case of one detector, i.e., Eq. 18. For $\eta<1$ the computer
simulations run extremely slowly for higher values of $m$ (reminiscent of the
boson sampling problem). We calculated the conditional probabilities at $m=50$
for $\eta=0.9,0.99,0.999$, and plot the highest losses for illustration in
Fig.6. Plots for other values of $\eta$’s are included in the supplementary
material.

Figure 6: Conditional probabilities $P_{m}(k|n)$ versus $n$, for $m=50$ APDs
and $\eta=0.9$.

The degradation of the count probability with photon loss is evident, compared
to Fig.4. Also recall that $\eta=0.9$ means 10% loss per mode which is a very
poor performance. The purity calculation, displayed in Fig.7, is particularly
illuminating.

Figure 7: POVM purity, versus click number $k$, for several values of $\eta$
at $m=50$

Indeed, it is clear that, as $\eta$ increases beyond the low $eta=0.9$ level,
the photon losses have decreasing to negligible ($\eta=0.999$) effect on
purity, which is essentially limited by $m$, as per Fig.5. This is an
interesting result. It is likely that the same level of photon loss may have a
more detrimental effect as $m$ increases, however, the exact scaling of this
effect is not yet known, due to the long computation times for the nonideal
case.

##  IV Conclusion

We carried out the theoretical evaluation of the photon-count POVM for the
segmented detector. Results show that PNR detection is indeed achievable in
the ideal case. This opens a new path to PNR devices that operate at room
temperature and can be manufactured with available integrated photonic
technology. The number of integrated APDs appears to be, to a large extent,
the dominant factor toward high-quality PNR detection. While photon losses
must also be taken in to account, of course, it is important to note that they
do not include the quantum efficiency of the APDs, by design of the segmented
detector. The reduction of photon losses will therefore only involve passive
optical design considerations, a notable difference with butt-coupled tree-
splitting detectors.

##  V acknowledgements

We are grateful to Joe Campbell, Seth Bank, Aye L. Win, Rafael Alexander, Ben
Godek, Sharon S. Philip, Bargav Jayaraman, and Oshin Jakhete for stimulating
discussions. This work was supported by the U.S. Defense Advanced Research
Projects Agency (DARPA).

## References

  * Knill  _et al._ (2001) E. Knill, R. Laflamme, and G. J. Milburn, Nature 409, 46 (2001).
  * Zwinkels  _et al._ (2010) J. C. Zwinkels, E. Ikonen, N. P. Fox, G. Ulm, and M. L. Rastello, Metrologia 47, R15 (2010).
  * O’Brien  _et al._ (2009) J. L. O’Brien, A. Furusawa, and J. Vuckovic, Nat Photon 3, 687 (2009).
  * G.  _et al._ (2010) B. G., G. M., and R. BercheraI., Nat Photon 4, 227 (2010).
  * Sangouard  _et al._ (2007) N. Sangouard, C. Simon, J. c. v. Minář, H. Zbinden, H. de Riedmatten, and N. Gisin, Phys. Rev. A 76, 050301 (2007).
  * Kok  _et al._ (2002) P. Kok, H. Lee, and J. P. Dowling, Phys. Rev. A 65, 052104 (2002).
  * Ramilli  _et al._ (2010) M. Ramilli, A. Allevi, V. Chmill, M. Bondani, M. Caccia, and A. Andreoni, J. Opt. Soc. Am. B 27, 852 (2010).
  * Divochiy  _et al._ (2008) A. Divochiy, F. Marsili, D. Bitauld, A. Gaggero, R. Leoni, F. Mattioli, A. Korneev, V. Seleznev, N. Kaurova, O. Minaeva, G. Gol’tsman, K. G. Lagoudakis, M. Benkhaoul, F. Levy, and A. Fiore, Nat Photon 2, 302 (2008).
  * Natarajan  _et al._ (2013) C. M. Natarajan, L. Zhang, H. Coldenstrodt-Ronge, G. Donati, S. N. Dorenbos, V. Zwiller, I. A. Walmsley, and R. H. Hadfield, Opt. Express 21, 893 (2013).
  * Mattioli  _et al._ (2015) F. Mattioli, Z. Zhou, A. Gaggero, R. Gaudio, S. Jahanmirinejad, D. Sahin, F. Marsili, R. Leoni, and A. Fiore, Superconductor Science and Technology 28, 104001 (2015).
  * Kardynał  _et al._ (2007) B. E. Kardynał, S. S. Hees, A. J. Shields, C. Nicoll, I. Farrer, and D. A. Ritchie, Applied Physics Letters 90, 181114 (2007), http://dx.doi.org/10.1063/1.2735281 .
  * Shields  _et al._ (2000) A. J. Shields, M. P. O’Sullivan, I. Farrer, D. A. Ritchie, R. A. Hogg, M. L. Leadbeater, C. E. Norman, and M. Pepper, Applied Physics Letters 76, 3673 (2000), http://dx.doi.org/10.1063/1.126745 .
  * Rosenberg  _et al._ (2005) D. Rosenberg, A. E. Lita, A. J. Miller, and S. W. Nam, Phys. Rev. A 71, 061803 (2005).
  * Lita  _et al._ (2008) A. E. Lita, A. J. Miller, and S. W. Nam, Opt. Express 16, 3032 (2008).
  * Achilles  _et al._ (2003) D. Achilles, C. Silberhorn, C. Śliwa, K. Banaszek, and I. A. Walmsley, Opt. Lett. 28, 2387 (2003).
  * Fitch  _et al._ (2003) M. J. Fitch, B. C. Jacobs, T. B. Pittman, and J. D. Franson, Phys. Rev. A 68, 043814 (2003).
  * Piacentini  _et al._ (2015) F. Piacentini, M. P. Levi, A. Avella, M. López, S. Kück, S. V. Polyakov, I. P. Degiovanni, G. Brida, and M. Genovese, Opt. Lett. 40, 1548 (2015).
  * Elitzur and Vaidman (1993) A. C. Elitzur and L. Vaidman, Foundations of Physics 23, 987 (1993).
  * van Enk (2017a) S. J. van Enk, arXiv:1705.09640 [quant-ph] (2017a), arXiv:1705.09640 [quant-ph] .
  * van Enk (2017b) S. J. van Enk, arXiv:1705.09033 [quant-ph] (2017b), arXiv:1705.09033 [quant-ph] .
  * Scully and Lamb (1969) M. O. Scully and W. E. Lamb, Jr, Phys. Rev. 179, 368 (1969).