Delaunay Weighted Two-sample Test for High-dimensional Data by Incorporating Geometric Information

Considering the simplest case $\mathcal{M}=\mathcal{R}^{d}$, for each internal point $\textbf{p}\in\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})=\mathcal{S}_{1}\cup\cdots\cup\mathcal{S}_{m}$, there exists a Delaunay simplex $\mathcal{S}(\textbf{p};\mathbb{Z})\in\mathcal{DT}_{\mathcal{R}^{d}}(\mathbb{Z})$ such that $\textbf{p}\in\mathcal{S}(\textbf{p};\mathbb{Z})$. We define the Delaunay weight vector for $\textbf{p}\in\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})$ with respect to $\mathbb{Z}$ as a vector $\textbf{w}(\textbf{p};\mathbb{Z})=(w_{1}(\textbf{p};\mathbb{Z}),\ldots,w_{n}(\textbf{p};\mathbb{Z}))^{\mathsf{T}}$ such that

• (a)

  $w_{i}(\textbf{p};\mathbb{Z})\in[0,1]$ if $\textbf{z}_{i}\in\mathbb{V}(\mathcal{S}(\textbf{p};\mathbb{Z})),$ and $w_{i}(\textbf{p};\mathbb{Z})=0$ if $\textbf{z}_{i}\notin\mathbb{V}(\mathcal{S}(\textbf{p};\mathbb{Z}))$, where $\mathbb{V}(\mathcal{S}(\textbf{p};\mathbb{Z}))$ denotes the vertex set of the simplex $\mathcal{S}(\textbf{p};\mathbb{Z})$;

• (b)

  $\sum_{i=1}^{n}w_{i}(\textbf{p};\mathbb{Z})=1$;

• (c)

  $\sum_{i=1}^{n}w_{i}(\textbf{p};\mathbb{Z})\cdot\textbf{z}_{i}=\textbf{p}$.

According to Abadie and L’Hour, (2021), the Delaunay weight vector $\textbf{w}(\textbf{p};\mathbb{Z})$ for $\textbf{p}\in\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})$ corresponds to the convex combination of $\mathbb{Z}$ with the smallest compound discrepancy $\sum_{i=1}^{n}w_{i}\|\textbf{p}-\textbf{z}_{i}\|^{2}$ among all convex combinations satisfying $\textbf{p}=\sum_{i=1}^{n}w_{i}\textbf{z}_{i}$. That is, $\textbf{w}(\textbf{p};\mathbb{Z})$ selects at most $d+1$ points in $\mathbb{Z}$ (those with $w_{i}>0$) that are convexly combined to form p and closest to p in the sense of the minimal compound discrepancy.

Conditions (a)–(c) only define the Delaunay weight vector for an internal point $\textbf{p}\in\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})$. For generalization to an external point $\textbf{p}\notin\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})$, we utilize the projection idea from Section 2.3 of Abadie and L’Hour, (2021). Let $\textbf{w}^{(\lambda)}(\textbf{p};\mathbb{Z})=(w^{(\lambda)}_{1}(\textbf{p};\mathbb{Z}),\ldots,w^{(\lambda)}_{n}(\textbf{p};\mathbb{Z}))^{\mathsf{T}}$ be the solution to the optimization problem,

subject to $\textbf{w}=(w_{1},\ldots,w_{n})^{\mathsf{T}}\in[0,1]^{n}$ and $\sum_{i=1}^{n}w_{i}=1$, where $\|\cdot\|$ denotes the Euclidean norm. The Delaunay weight vector $\textbf{w}(\textbf{p};\mathbb{Z})$ for $\textbf{p}\in\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})$ can be rewritten as $\textbf{w}(\textbf{p};\mathbb{Z})=\lim\limits_{\lambda\to 0}\textbf{w}^{(\lambda)}(\textbf{p};\mathbb{Z})$. For $\textbf{p}\notin\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})$, Abadie and L’Hour, (2021) showed that $\lim\limits_{\lambda\to 0}\sum_{i=1}^{n}{w}_{i}^{(\lambda)}(\textbf{p};\mathbb{Z})\textbf{z}_{i}=\textbf{p}^{*},$ where $\textbf{p}^{*}$ is the projection of p on the convex hull $\mathcal{H}_{\mathcal{R}^{d}}(\mathbb{Z})$. Therefore, we have

Based on (1), we can generalize the definition of the Delaunay weight for all $\textbf{p}\in\mathcal{R}^{d}$ by substituting condition (c) with

• (c^†)

  $\sum_{i=1}^{n}w_{i}(\textbf{p};\mathbb{Z})\cdot\textbf{z}_{i}=\textbf{p}^{*}.$

However, when $\mathcal{M}$ is not Euclidean, condition (c^†) is no longer applicable because it relies on the linear operation of Cartesian coordinates, which is not well-defined on $\mathcal{M}$. To generalize the Delaunay weight vector to nonlinear manifolds, we make use of the logarithmic map.

For the special case that $\mathcal{M}=\mathcal{R}^{d}$, (c^†) and (c^‡) are equivalent because

where $\overrightarrow{\textbf{p}^{*}\textbf{z}_{i}}$ denotes the vector from $\textbf{p}^{*}$ to $\textbf{z}_{i}$. Therefore, condition (c^‡) is indeed a generalization of (c^†) on a nonlinear manifold. Since the Delaunay triangulation $\mathcal{DT}_{\mathcal{M}}(\mathbb{Z})$ is unique for generic $\mathbb{Z}$ (Leibon and Letscher,, 2000), the Delaunay weight vector $\textbf{w}(\textbf{p};\mathbb{Z})$ is unique for generic $\mathbb{Z}$ and any point $\textbf{p}\in\mathcal{M}$, as shown in Section S1 of the supplementary material.

Based on Definition 2, the Delaunay weight matrix of $\mathbb{Z}$ is defined as $\boldsymbol{\Gamma}_{\mathbb{Z}}=(\gamma_{ij})_{n\times n}$, with

where $\gamma_{ij}$ measures the geometric proximity of $\textbf{z}_{j}$ to $\textbf{z}_{i}$ and $\mathbb{Z}\setminus\{\textbf{z}_{i}\}$ represents the remaining set of $\mathbb{Z}$ excluding $\textbf{z}_{i}$. Because $\gamma_{ij}$ and $\gamma_{ji}$ are usually not the same, the Delaunay weight matrix $\boldsymbol{\Gamma}_{\mathbb{Z}}$ is a directed weighted adjacency matrix of $\mathbb{Z}$. Unlike the kernel, local regression and $k$-NN based approaches which only consider local distances to $\textbf{z}_{i}$, the Delaunay weight $\textbf{w}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ evaluates the geometric proximity of the vertices of the Delaunay simplex $\mathcal{S}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ to $\textbf{z}_{i}$ using both the distance and direction information. In particular, $\gamma_{ij}$ and $\gamma_{ik}$ can be different even if $\|\log_{\textbf{z}_{i}^{*}}(\textbf{z}_{j})\|=\|\log_{\textbf{z}_{i}^{*}}(\textbf{z}_{k})\|$. Indeed, we see from Definition 2 that $\textbf{w}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ depends on the vectors $\log_{\textbf{z}_{i}^{*}}(\textbf{z}_{j})$ rather than merely the distances $\|\log_{\textbf{z}_{i}^{*}}(\textbf{z}_{j})\|$, for $j\neq i$.

To gain more insights on the Delaunay weight, we see that, when $\textbf{z}_{i}$ is an inner point of the Delaunay simplex $\mathcal{S}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$, the non-zero elements of $\textbf{w}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ correspond to the points in $\mathbb{Z}\setminus\{\textbf{z}_{i}\}$ from all directions surrounding $\textbf{z}_{i}$. In contrast, under sparse or nonuniform sampling, nearest neighbors of $\textbf{z}_{i}$ cannot cover all directions from $\textbf{z}_{i}$ because they tend to cluster in certain directions. We illustrate this point using an example with a sample $\mathbb{Z}=\{\textbf{z}_{1},\ldots,\textbf{z}_{8}\}$ in $\mathcal{M}=\mathcal{R}^{2}$. Figure 3 shows the Delaunay simplex $\mathcal{S}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ that contains the projection of $\textbf{z}_{i}$ on the convex hull $\mathcal{H}_{\mathcal{M}}(\mathbb{Z}\setminus\{\textbf{z}_{i}\})$. Based on $\mathcal{S}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ ($i=1,\ldots,8$), the weight vectors $\textbf{w}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ are concatenated to form the Delaunay weight matrix $\boldsymbol{\Gamma}_{\mathbb{Z}}$, as shown in Figure 4. For $\textbf{z}_{2}\in\mathcal{H}_{\mathcal{M}}(\mathbb{Z}\setminus\{\textbf{z}_{2}\})$, the three nearest neighbors of $\textbf{z}_{2}$ are $\textbf{z}_{3},\textbf{z}_{7}$ and $\textbf{z}_{8}$, which are all from the top-right direction of $\textbf{z}_{2}$ and the information from the bottom-left direction is not captured. In contrast, the Delaunay simplex $\mathcal{S}(\textbf{z}_{2};\mathbb{Z}\setminus\{\textbf{z}_{2}\})$ is formed by $\textbf{z}_{5},\textbf{z}_{7}$ and $\textbf{z}_{8}$, which covers all directions from $\textbf{z}_{2}$, and the relative directions of $\overrightarrow{\textbf{z}_{2}\textbf{z}_{k}}$ $(k=5,7,8)$ are used in constructing $\textbf{w}(\textbf{z}_{2};\mathbb{Z}\setminus\{\textbf{z}_{2}\})$. Therefore, the Delaunay weight utilizes the local geometric information in a more comprehensive way than the $k$-NN.

For the pooled sample $\mathbb{Z}=\{\textbf{z}_{1},\ldots,\textbf{z}_{n}\}$ with the group indicator $\boldsymbol{\delta}=\{\delta_{1},\ldots,\delta_{n}\}$, the $k$-NN test statistic in Schilling, (1986) is

which has inspired a family of nonparametric two-sample tests for detecting distribution differences via a proximity graph of the pooled sample (Rosenbaum,, 2005; Chen and Friedman,, 2017; Chen et al.,, 2018). Based on the Delaunay weight $\gamma_{ij}$, we construct the test statistic

where the new geometric proximity measure $\gamma_{ij}$ accounts for the manifold structure and direction information.

In practice, the underlying manifold $\mathcal{M}$ is typically unknown, so that exact computation of the Delaunay weight matrix $\boldsymbol{\Gamma}_{\mathbb{Z}}$ via Definition 2 (and thus $T_{\text{DW}}$) is infeasible. Because $\boldsymbol{\Gamma}_{\mathbb{Z}}$ is constructed by Delaunay weight vectors $\textbf{w}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$ which only measure the geometric proximity between $\textbf{z}_{i}$ and the points locally surrounding $\textbf{z}_{i}$ (i.e., the vertices of $\mathcal{S}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\})$), we can utilize manifold learning techniques to approximate $\boldsymbol{\Gamma}_{\mathbb{Z}}$ by $\boldsymbol{\Gamma}_{\widetilde{\mathbb{Z}}}$, where $\widetilde{\mathbb{Z}}=\{\tilde{\textbf{z}}_{1},\ldots,\tilde{\textbf{z}}_{n}\}\subset\mathcal{R}^{d}$ is a low-dimensional Euclidean representation of $\mathbb{Z}$ whose Delaunay weight matrix can be computed. Generally speaking, any nonlinear dimension reduction technique that aims to preserve the local geometric structure of $\mathbb{Z}$ can be applied to obtain $\widetilde{\mathbb{Z}}$. Often, a $k$-NN nonlinear dimension reduction technique may yield a proximity graph of $\widetilde{\mathbb{Z}}$ with more than one connected components. Under such cases, elements in $\boldsymbol{\Gamma}_{\widetilde{\mathbb{Z}}}$ corresponding to all pairs of points from different components must be zero, which incurs information loss and power deterioration. To circumvent this issue, we suggest adopting a procedure that can return a connected low-dimensional representation $\widetilde{\mathbb{Z}}$. Specifically, we modify the robust manifold learning algorithm in Dimeglio et al., (2014), as detailed in Section S2 of the supplementary material, to estimate the geodesic distances ${d}_{\mathcal{M}}(\textbf{z}_{i},\textbf{z}_{j})$ and apply the classical MDS (Kruskal,, 1964) to obtain $\widetilde{\mathbb{Z}}\subset\mathcal{R}^{{d}}$. For the unknown intrinsic dimension $d$, any existing method can be applied to estimate $d$, e.g., Facco et al., (2017) as summarized in Section S2 of the supplementary material. Our numerical experiments show that our method is robust to estimation of $d$.

To compute $\boldsymbol{\Gamma}_{\widetilde{\mathbb{Z}}}$, we note the equivalence of conditions (c^†) and (c^‡) in (2) for the low-dimensional Euclidean representation $\widetilde{\mathbb{Z}}$. Once we know the Delaunay simplex $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$ that contains $\tilde{\textbf{z}}_{i}^{*}$, which is the projection of $\tilde{\textbf{z}}_{i}$ on the convex hull $\mathcal{H}_{\mathcal{R}^{d}}(\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$, we can easily compute $\textbf{w}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})=(\tilde{\gamma}_{i1},\ldots,\tilde{\gamma}_{i,i-1},\tilde{\gamma}_{i,i+1},\ldots,\tilde{\gamma}_{in})^{\mathsf{T}}$ by solving the linear equations,

The remaining difficulty in computing $\boldsymbol{\Gamma}_{\widetilde{\mathbb{Z}}}$ lies in finding the Delaunay simplex $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$, for $i=1,\ldots,n$.

Although the DELAUNAYSPARSE algorithm by Chang et al., (2020) can efficiently find $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$ for $\tilde{\textbf{z}}_{i}\in\mathcal{H}_{\mathcal{R}^{d}}(\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$, it does not work for those $\tilde{\textbf{z}}_{i}\notin\mathcal{H}_{\mathcal{R}^{d}}(\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$. To tackle this difficulty, we incorporate the inverse stereographic projection to the DELAUNAYSPARSE algorithm and develop the stereographic projected DELAUNAYSPARSE algorithm to find $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$, regardless of whether $\tilde{\textbf{z}}_{i}\in\mathcal{H}_{\mathcal{R}^{d}}(\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$ or not. Without loss of generality, we assume the Euclidean representation $\widetilde{\mathbb{Z}}$ is centered at $\textbf{0}^{d}$. The stereographic projected DELAUNAYSPARSE algorithm is given as follows.

• 1.

  We first compute $r_{\text{max}}=\max_{i=1,\ldots,n}\|\tilde{\textbf{z}}_{i}\|$, set a scaling parameter $\eta>0$ and apply the inverse stereographic projection,

  $$\phi_{\eta r_{\text{max}}}^{-1}:\mathcal{R}^{d}\to\overline{\mathcal{B}}_{\eta r_{\text{max}}}\setminus\{(\mathbf{0}_{d},\eta r_{\text{max}})^{\mathsf{T}}\},\quad\tilde{\textbf{z}}\mapsto\Bigg{(}\frac{\eta^{2}r_{\text{max}}^{2}}{\eta^{2}r_{\text{max}}^{2}+\|\tilde{\textbf{z}}\|^{2}}\cdot\tilde{\textbf{z}},\frac{\eta r_{\text{max}}\|\tilde{\textbf{z}}\|^{2}}{\eta^{2}r_{\text{max}}^{2}+\|\tilde{\textbf{z}}\|^{2}}\Bigg{)}^{\mathsf{T}},$$

  (6)

  to project $\widetilde{\mathbb{Z}}$ to a $d$-sphere $\overline{\mathcal{B}}_{\eta r_{\text{max}}}=\{\textbf{z}\in\mathcal{R}^{d+1};\|\textbf{z}-(\mathbf{0}_{d},\eta r_{\text{max}}/2)^{\top}\|=\eta r_{\text{max}}/2\}$. Let $\check{\textbf{z}}_{i}=\phi_{\eta r_{\text{max}}}^{-1}(\tilde{\textbf{z}}_{i})$, for $i=1,\ldots,n$, and let the reference point $\check{\textbf{z}}_{0}$ be $(\mathbf{0}_{d},\eta r_{\text{max}})$. We define the augmented point set as $\check{\mathbb{Z}}=\{\check{\textbf{z}}_{0},\check{\textbf{z}}_{1},\ldots,\check{\textbf{z}}_{n}\}\subset\overline{\mathcal{B}}_{\eta r_{\text{max}}}$. A graphical illustration of obtaining $\check{\mathbb{Z}}$ from ${\widetilde{\mathbb{Z}}}$ with $d=2$ is exhibited in Figure 5 (a).

• 2.

  As shown in Figure 5 (b), for any $i=1,\ldots,n$, the union of simplicies in the Delaunay triangulation $\mathcal{DT}_{\overline{\mathcal{B}}_{\eta r_{\text{max}}}}(\check{\mathbb{Z}}\setminus\{\check{\textbf{z}}_{i}\})$ covers the whole sphere $\overline{\mathcal{B}}_{\eta r_{\text{max}}}$. By projecting $\mathcal{DT}_{\overline{\mathcal{B}}_{\eta r_{\text{max}}}}(\check{\mathbb{Z}}\setminus\{\check{\textbf{z}}_{i}\})$ back to $\mathcal{R}^{d}$ under the convention that the geodesic path $\ell_{\overline{\mathcal{B}}_{\eta r_{\text{max}}}}(\check{\textbf{z}}_{j},\check{\textbf{z}}_{0})$ (for any $j\neq i$) is mapped to the line in $\mathcal{R}^{d}$ from $\tilde{\textbf{z}}_{j}$ in the opposite direction of $\mathbf{0}_{d}$, an approximation of $\mathcal{DT}_{\mathcal{R}^{d}}({\widetilde{\mathbb{Z}}}\setminus\{\tilde{\textbf{z}}_{i}\})$ is obtained as shown in Figure 5 (c). Figure 5 (d) shows that such an approximation converges to $\mathcal{DT}_{\mathcal{R}^{d}}({\widetilde{\mathbb{Z}}}\setminus\{\tilde{\textbf{z}}_{i}\})$ as the scaling parameter $\eta$ tends to infinity.

• 3.

  Because $\mathcal{DT}_{\overline{\mathcal{B}}_{\eta r_{\text{max}}}}(\check{\mathbb{Z}}\setminus\{\check{\textbf{z}}_{i}\})$ covers $\overline{\mathcal{B}}_{\eta r_{\text{max}}}$, for any $i=1,\ldots,n$, there must exist a Delaunay simplex $\mathcal{S}^{*}\in\mathcal{DT}_{\overline{\mathcal{B}}_{\eta r_{\text{max}}}}(\check{\mathbb{Z}}\setminus\{\check{\textbf{z}}_{i}\})$ that contains $\check{\textbf{z}}_{i}$ on $\overline{\mathcal{B}}_{\eta r_{\text{max}}}$. Thus, we use the breadth first search procedure similar to the DELAUNAYSPARSE algorithm (Chang et al.,, 2020) to compute the Delaunay simplex $\mathcal{S}^{*}$. Finally, if $\check{\textbf{z}}_{0}\notin\mathbb{V}(\mathcal{S}^{*})$, $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$ on $\mathcal{R}^{d}$ is obtained by the simplex formed by $\{\tilde{\textbf{z}}_{j}:\check{\textbf{z}}_{j}\in\mathbb{V}(\mathcal{S}^{*})\}$; if $\check{\textbf{z}}_{0}\in\mathbb{V}(\mathcal{S}^{*})$, $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$ is obtained by the simplex formed by $\{\tilde{\textbf{z}}_{j}:\check{\textbf{z}}_{j}\in\mathbb{V}(\mathcal{S}^{\dagger})\}$, where $\mathcal{S}^{\dagger}$ is the neighbor simplex of $\mathcal{S}^{*}$ opposite to $\check{\textbf{z}}_{0}$.

The stereographic projected DELAUNAYSPARSE algorithm is detailed in Section S3 of the supplementary material. Note that we only need to generate the local (but not full) Delaunay triangulation to find the Delaunay simplex $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$, and thus the computation is not intensive. Specifically, Section S3 of the supplementary material shows that the computational complexity of the stereographic projected DELAUNAYSPARSE algorithm is $O(n^{2+1/d}d^{2})$. As the computational complexity of solving (5) via Gaussian elimination is $O(nd^{3})$ and the sample size $n$ is usually much larger than the intrinsic dimension $d$, we conclude that the overall computational complexity for computing the approximation $\boldsymbol{\Gamma}_{\widetilde{\mathbb{Z}}}$ is $O(n^{2+1/d}d^{2})$. As shown in Figure 5 (d), the difference between $\mathcal{DT}_{\mathcal{R}^{d}}({\widetilde{\mathbb{Z}}}\setminus\{\tilde{\textbf{z}}_{i}\})$ and its approximation with $\eta>5$ is negligible. Because computation of the Delaunay weight is unrelated to $\eta$ provided that the simplex $\mathcal{S}(\tilde{\textbf{z}}_{i};\widetilde{\mathbb{Z}}\setminus\{\tilde{\textbf{z}}_{i}\})$ is correctly identified, we suggest $\eta=10$ in practice.

The Delaunay weighted test statistic in (4), ${T}_{\text{DW}}=\sum_{j=1}^{n}\sum_{i\neq j}\gamma_{ij}\times I(\delta_{i}=\delta_{j})$, aims to detect differences between the underlying distributions $F$ and $G$ via the deviation of the sum of within-group weights from its expectation under the null. Figure 4 gives an illustrative example of the Delaunay weight matrix $\boldsymbol{\Gamma}_{\mathbb{Z}}$ with group indicators $\delta_{i}=1$, for $i=1,\ldots,4$, and $\delta_{i}=0$ otherwise. If we partition $\Gamma_{\mathbb{Z}}$ at the 4-th row and column, ${T}_{\text{DW}}$ equals the sum of all Delaunay weights in the two diagonal blocks of $\boldsymbol{\Gamma}_{\mathbb{Z}}$. Generally, off-diagonal elements in $\boldsymbol{\Gamma}_{\mathbb{Z}}$ are identically distributed under the null. In contrast, under $H_{1}:F\neq G$, more $\textbf{z}_{i}$’s in $\mathbb{X}$ are located in the region where the density of $F$ is higher than $G$. For these $\textbf{z}_{i}$’s in $\mathbb{X}$, the expected proportion of data points from $\mathbb{X}$ in $\mathbb{V}(\mathcal{S}(\textbf{z}_{i};\mathbb{Z}\setminus\{\textbf{z}_{i}\}))$ is larger than that in $\mathbb{Z}\setminus\{\textbf{z}_{i}\}$. Similar arguments also hold for $\mathbb{Y}$. Thus, ${T}_{\text{DW}}$ is stochastically larger under the alternative than under the null.

The exact distribution of $T_{\text{DW}}$ under the null is complex and difficult to derive, and thus we implement a permutation procedure to compute the $p$-value.

1. 1.

   Compute the observed value of test statistic $T_{\text{DW}}$ via (4).

2. 2.

   For $b=1,\ldots,B$, we randomly generate a permutation of group indicators $\boldsymbol{\delta}^{(b)}=\{\delta_{1}^{(b)},\ldots,\delta_{n}^{(b)}\}$, where $\sum_{i=1}^{n}\delta_{i}^{(b)}=n_{1}$ and $\sum_{i=1}^{n}(1-\delta_{i}^{(b)})=n_{0}$, and compute $T^{(b)}_{\text{DW}}=\sum_{j=1}^{n}\sum_{i\neq j}\gamma_{ij}\times I(\delta_{i}^{(b)}=\delta_{j}^{(b)})$ under $\boldsymbol{\delta}^{(b)}$.

3. 3.

   Compute $p{\textrm{-value}}=\{\sum_{b=1}^{B}I(T_{\text{DW}}\leq T^{(b)}_{\text{DW}})+1\}/(B+1),$ where $1$ is added in both the numerator and denominator to calibrate that under $H_{0}$, $\text{Pr}(p\text{-value}\leq\alpha)\leq\alpha$ for all $\alpha\in[0,1]$. Indeed, under $H_{0}$, $T_{\text{DW}},T_{\text{DW}}^{(1)},\ldots,,T_{\text{DW}}^{(B)}$ are i.i.d. given the pooled sample $\mathbb{Z}$ and sample sizes $(n_{1},n_{0})$. Therefore, the $p\text{-value}$ is uniformly distributed on $\{1/(B+1),2/(B+1),\ldots,1\}$, which converges to $\text{Unif}[0,1]$ as $B\to\infty$.

Under the unknown manifold setting where $\Gamma_{\mathbb{Z}}$ cannot be computed exactly, we replace it by $\Gamma_{\widetilde{\mathbb{Z}}}$ for approximate inference. Unlike most of the existing tests based on pairwise distances only, the Delaunay weight takes the relative directions among data points into account, and thus our test gains more power in detecting the direction difference between $F$ and $G$.