BEAUTY Powered BEAST

Our study of the BET framework is inspired by the celebrated Euler’s formula,

which is often regarded as one of the most beautiful equations in mathematics. In particular, when $x=\pi$, one has Euler’s identity, $e^{i\pi}+1=0$, which connects the five most important numbers in mathematics $0,1,i,e,\pi$ in one simple yet deep equation. Beside the beauty of this equation, how is it useful for statisticians? To see that, consider any binary variable $A$ (not necessarily symmetric) which takes values $-1$ or $1$. Through the parity of the sine and cosine functions, one can easily show the following binary Euler’s equation. Since we were not aware of any reference of this equation in literature, we formally state it below.

Lemma 1.2 generalizes Euler’s formula with additional randomness from a binary variable and reduces its complex exponentiation to its linear polynomial. To the best of our knowledge, no other random variables enjoy the same remarkable attribute. Moreover, note that the random variable $e^{iAx}$ in (1.5) is closely related to characteristic functions, particularly when it is combined with the binary expansion in Theorem 1.1. For example, for $U=\sum_{d=1}^{\infty}2^{-d}{A_{d}}\sim\text{Unif}[-1,1]$ and for any $t\in\mathbb{R},$ we have

The complex exponent of $U$ can be approximated by a polynomial of the binary variables in its binary expansion! Moreover, we show in Section 2 that this approximation is universal for any $p$-dimensional vector supported within $[-1,1]^{p}$. We refer this universal binary interaction approximation of the complex exponent and the characteristic function as the Binary Expansion Approximation of UniformiTY (BEAUTY) in Theorem 2.2.

Based on the BEAUTY, in this paper we make the following three main novel contributions to the problem of nonparametric tests of independence:

1. A unification of important nonparamatric tests of independence. In Section 3, we show that many important tests of independence in literature can be approximated by some quadratic forms of symmetry statistics, which are shown to be complete sufficient statistics for dependence in Zhang (2019). In particular, each of these test statistics corresponds to a different deterministic weight matrix in the quadratic form, which in turn dictates the power properties of the test. Therefore, this deterministic weight in existing test statistics creates the key issue on uniformity and robustness of the test, as it may favor certain alternatives but cause a substantial loss of power for other alternatives. Following this observation, we consider a test statistic that has data-adaptive weights to make automatic adjustments under different situations so as to achieve a robust power. We refer this test as the Binary Expansion Adaptive Symmetry Test (BEAST), as described in Section 4.

2. A benchmark of feasible power from the BEAST with oracle. We begin by considering the test of uniformity. By utilizing the properties of the binary expansion filtration, we show in a heuristic asymptotic study of the BEAST a surprising fact that for any given alternative, the Neyman-Pearson test for testing uniformity can be approximated by a weighted sum of symmetry statistics. We thus develop the BEAST through an oracle approach over this Neyman-Pearson one-dimensional projection of symmetry statistics, which quantifies a boundary of feasible power performance. Numerical studies in Section 5 show that the BEAST with oracle leads a wide range of prevailing tests by a surprisingly huge margin under all alternatives we considered. This enormous margin thus provides helpful information about the potential of substantial power improvement for each alternative. To the best of our knowledge, there is no other type of similar approach or results to study the potential performance of a test of uniformity or independence. Therefore, the BEAST with oracle sets a novel and useful benchmark for the feasible power under any alternative. Moreover, it provides guidance for choosing suitable weights to boost the power of the test.

3. A powerful and robust BEAST from a regularized resampling approximation of the oracle. Motivated by the form of the BEAST with oracle, we construct the practical BEAST to approximate the optimal power by approximating the oracle weights in testing uniformity. The proposed BEAST combines the ideas of resampling and regularization to obtain data-adaptive weights that adjusts the statistic towards the oracle under each alternative. Here resampling helps the approximation of the sampling distribution of the oracle test statistic, and regularization screens the noise in the estimation of optimal weights. This test is applied to the problem to the test of independence. Simulation studies in Section 5 demonstrate that the BEAST improves the power of many existing tests of univariate or multivariate independence against many common forms of non-uniformity, particularly multimodal and nonlinear ones. Besides its robust power, the BEAST provides clear and meaningful interpretations of statistical significance, which we demonstrate in Section 6. We conclude our paper with discussions in Section 7. Details of notation, theoretical proofs and additional numerical results are deferred to Supplementary materials.

As a robust version of the Pearson correlation, the Spearman’s $\rho$ statistic leads to a test with high asymptotic relative efficiency compared to the optimal test with Pearson correlation under bivariate normal distribution (Lehmann and Romano, 2006). We show below it can be approximated by a quadratic form of symmetry statistics.

When $U_{1}$ and $U_{2}$ are marginally uniformly distributed over $[-1,1]$, Spearman’s $\rho$ can be written as the correlation between $U_{1}$ and $U_{2}$, i.e.,

where $\mathcal{L}_{2,D,\text{spe}}=\{\Lambda=\Lambda_{1}\raisebox{0.5pt}{\raisebox{-.9pt} {r}⃝ }\Lambda_{2}:\Lambda_{1},\Lambda_{2}\in\mathbf{B}^{1\times D},\text{ where }\Lambda_{1}{\bm{1}}=1\text{ and }\Lambda_{2}{\bm{1}}=1\}$ consists of $2\times D$ matrices whose rows are both binary vectors with only one unique $1$, and the $D^{2}$-dimensional vector $\bm{r}_{D}$ has entry $2^{-(d_{1}+d_{2})}$ corresponding to $\mathbf{E}[A_{1,d_{1}}A_{2,d_{2}}].$ The test based on Spearman’s $\rho$ rejects the null when the estimate of $\rho$ has a large absolute value. This test statistic can be approximated with

which is a quadratic form with a rank-one weight matrix $\mathbf{W}_{\rho,D}=\bm{r}_{D}\bm{r}_{D}^{T}.$

Although the test based on Spearman’s $\rho$ has a higher power against the linear form of dependency particularly present in bivariate normal distributions, we see from $\mathcal{L}_{2,D,\text{spe}}$ and $\mathbf{W}_{\rho,D}$ that this test only considers $D^{2}$ out of $(2^{D}-1)^{2}$ cross interactions of binary variables in $\mathcal{L}_{2,D,\text{cross}}$. Thus this test is not capable of detecting complex nonlinear forms of dependency.

When $U_{1}$ and $U_{2}$ are $\text{Unif}[-1,1]$ distributed, the binary expansion up to depth $D$ effectively leads to a discretization of $[-1,1]^{2}$ into a $2^{D}\times 2^{D}$ contingency table. Classical tests for contingency tables such as $\chi^{2}$-test can thus be applied. Similar tests include Fisher’s exact test and its extensions (Ma and Mao, 2019). Multivariate extensions of these methods include Gorsky and Ma (2018); Lee et al. (2023).

In Zhang (2019), it is shown that the $\chi^{2}$-statistic at depth $D$ can be written as the sum of squares of symmetry statistics for cross interactions. Thus,

where $\mathcal{L}_{2,D,\text{cross}}$ is the collection of all cross interactions. The weight matrix for $Q_{\chi^{2}}$ is thus the identity matrix $\mathbf{I}_{(2^{D}-1)\times(2^{D}-1)}$.

The Max BET proposed in Zhang (2019) can be approximated by a quadratic form with another diagonal weight matrix, which we explain in the Supplementary Materials. These tests with diagonal weights can detect signals among the squared symmetry statistics, but might be powerless for signals from their cross products.

To study the dependency between a $p_{1}$-dimensional vector $\bm{U}_{1}$ and a $p_{2}$-dimensional vector $\bm{U}_{2}$, in Székely et al. (2007), a class of measures of dependence is defined as

where $\phi_{(\bm{U}_{1},\bm{U}_{2})}(\bm{t}_{1},\bm{t}_{2})$ is the characteristic function of the joint distribution of $(\bm{U}_{1},\bm{U}_{2}),$ $w(\bm{t}_{1},\bm{t}_{2})$ is a suitable weight function, and $\phi_{\bm{U}_{k}}(\bm{t}_{k})$ is the characteristic function of $\bm{U}_{k},k=~{}1,2$. Note that $\mathcal{V}^{2}(\bm{U}_{1},\bm{U}_{2})=0$ if and only if $\bm{U}_{1}$ and $\bm{U}_{2}$ are independent. The distance correlation is then defined through $\mathcal{V}^{2}(\bm{U}_{1},\bm{U}_{2})$ and admits some desirable properties such as universal consistency against alternatives with finite expectation.

When $\bm{U}_{1}\sim\text{Unif}[-1,1]^{p_{1}}$ and $\bm{U}_{2}\sim\text{Unif}[-1,1]^{p_{2}}$, by Theorem 2.2, the term corresponding to $\Lambda={\bm{0}}$ cancels with $\phi_{\bm{U}_{1}}(\bm{t}_{1})\phi_{\bm{U}_{2}}(\bm{t}_{2})$, and we can write (3.2) as

where $\mathcal{L}_{p_{1}+p_{2},D,\text{unif}}=\{\Lambda\in\mathbf{B}^{(p_{1}+p_{2})\times D}:\Lambda\neq{\bm{0}}^{p\times D}\}$, and the weight matrix $\mathbf{W}_{\mathcal{V}^{2},p_{1},p_{2},D}$ consists of constants $w_{\Lambda_{1},\Lambda_{2}}$’s from the integration over $\bm{t}_{1}$ and $\bm{t}_{2}$. The test is significant when the empirical quadratic form $Q_{\mathcal{V}^{2},p1,p2,D}$ is large, where

Note that $\mathbf{W}_{\mathcal{V}^{2},p_{1},p_{2},D}$ here depends only on the weight function $w(\bm{t}_{1},\bm{t}_{2})$ and is deterministic. Hence, the test based on $Q_{\mathcal{V}^{2},p1,p2,D}$ will have a high power when the vector of $\mathbf{E}[A_{\Lambda}]$’s from the alternative distribution lies in the subspace spanned by eigenvectors of $\mathbf{W}_{\mathcal{V}^{2},p_{1},p_{2},D}$ corresponding to its largest eigenvalues. On the other hand, if instead the signals lie in the subspace spanned by eigenvectors of $\mathbf{W}_{\mathcal{V}^{2},p_{1},p_{2},D}$ corresponding to its lowest eigenvalues, then the power of the test could be considerably compromised. Therefore, a deterministic weight over symmetry statistics becomes a general uniformity issue of existing test statistics. In the next section, we study data-adaptive weights with the aim to improve the power by setting proper weights both among diagonal and off-diagonal entries in the matrix.

The unification in Section 3 inspires us to consider a class of nonparametric statistics for the goodness-of-fit test as a weighted sum of symmetry statistics. Since the properties of this form of statistics are closely related to the first two moments of the binary interaction variables in the filtration, we consider the collection of all nontrivial binary interactions $\mathcal{L}=\mathcal{L}_{p,D,\text{unif}}=\{\Lambda\in\mathbf{B}^{p\times D}:\Lambda\neq{\bm{0}}^{p\times D}\}$ and study the moment properties of the corresponding binary random vector $\bm{A}_{\mathcal{L}}.$

We begin by studying the connection between the $(2^{pD}-1)\times 1$ vector $\bm{A}_{\mathcal{L}}$ and the multinomial distribution from the corresponding discretization with $2^{pD}$ categories. We order the indices $\Lambda$’s in $\bm{A}_{\mathcal{L}}$ by the integer corresponding to the binary vector representation $vec(\Lambda^{T}),$ where $vec(\cdot)$ is the vectorization function. For example, the last (i.e. the $(2^{pD}-1)$th) entry in $\bm{A}_{\mathcal{L}}$ corresponds to the $\Lambda={\bm{1}}^{p\times D}$. We also denote the $2^{pD}\times 1$ vector of cell probabilities in the multinomial distribution by $\bm{p}_{c}.$ Label the entries in $\bm{p}_{c}$ by binary matrices $\Lambda\in\mathbf{B}^{p\times D}$ through $\Lambda=\Lambda_{1}\raisebox{0.5pt}{\raisebox{-.9pt} {r}⃝ }\ldots\raisebox{0.5pt}{\raisebox{-.9pt} {r}⃝ }\Lambda_{p},$ where each realization of the $2^{D}\times 1$ vector $\Lambda_{j}$ labels one of the $2^{D}$ intervals for dimension $j$ from low to high according to 1 plus the integer corresponding to the binary representation of $\Lambda_{j}^{T}$. We define a $2^{pD}\times 1$ random vector $\bm{Z}=(Z_{\Lambda})\sim Multinomial(1,\bm{p}_{c})$ to denote one draw from the $2^{pD}$ intervals from the discretization. With the above notation, we develop the general binary interaction design (BID) equation, which extends the two-dimensional case in Zhang (2019).

The Hadamard matrix $\mathbf{H}$ is also referred to as the Walsh matrix in engineering, where the linear transformation with $\mathbf{H}$ is referred to as the Hadamard transform (Lynn, 1973; Golubov et al., 2012; Harmuth, 2013). The earliest referral to the Hadamard matrix we found in the statistical literature is Pearl (1971), and it is also closely related to the orthogonal full factorial design (Cox and Reid, 2000; Box et al., 2005). In our context of testing independence, the BID equation can be regarded as a transformation from the physical domain to the frequency domain, which turns the focus to global forms of non-uniformity instead of local ones. In developing statistics, this transformation facilitates regularizations through thresholding, as $\bm{\mu}_{\mathcal{L}}={\bm{0}}$ is equivalent to uniformity $\bm{p}_{c}=1/2^{pD}{\bm{1}}.$ This transformation also enables clear interpretations of statistical significance with the form of dependency, as shown in Zhang (2019).

To study the power of the test of uniformity, we further study the properties of the first two moments of $\bm{A}_{\mathcal{L}}.$ Let $\bm{\mu}_{\mathcal{L}}=\mathbf{E}[\bm{A}_{\mathcal{L}}]$ and $\bm{\Sigma}_{\bm{\mu}_{\mathcal{L}}}=\mathbf{E}[\bm{A}_{\mathcal{L}}\bm{A}_{\mathcal{L}}^{T}]$ denote the vector of expectations and the matrix of second moments of $\bm{A}_{\mathcal{L}}$ respectively. We summarize some properties of $\bm{\mu}_{\mathcal{L}}$ and $\bm{\Sigma}_{\bm{\mu}_{\mathcal{L}}}$ in the following theorem.

To the best of our knowledge, the results in Theorem 4.2, despite their simplicity, have not been documented in literature. These simple results unveil interesting insights of the first two moments of binary variables in the filtration. The quadratic form in (4.4) characterizes the functional relationship between $\bm{\mu}_{c}$ and $\bm{\Sigma}_{\bm{\mu}_{c}}$. The two equations in (4.5) show that for binary variables, the Hotelling $T^{2}$ quadratic form is a monotone function of the harmonic mean of the cell probabilities in the corresponding multinomial distribution. The inequalities in (4.6) reveal the eigen structure of $\bm{\Sigma}_{\bm{\mu}_{\mathcal{L}}}$ when the signal $\bm{\mu}_{\mathcal{L}}$ is weak. The Taylor expansion in (4.7) provides the asymptotic behavior of $(\bm{\Sigma}_{\bm{\mu}_{\mathcal{L}}}-\bm{\mu}_{\mathcal{L}}\bm{\mu}_{\mathcal{L}}^{T})^{-1}\bm{\mu}_{\mathcal{L}}$ when the joint distribution is close to the uniform distribution. These insights shed important lights on how we can develop a powerful goodness-of-fit test, as we explain in Sections 4.2 and 4.3.

In this section, we study how to construct a powerful robust nonparametric test of uniformity based on what we learned in Sections 3 and 4.1. As discussed in Section 3, the deterministic weights of symmetry statistics in existing tests create an issue on the uniformity and robustness: They make the test powerful for some alternatives but not for others. Therefore, we construct a test statistic with data-adaptive weights, which allow the test to adjust itself towards the alternative to improve the power. We refer this class of statistics as the binary expansion adaptive symmetry test (BEAST).

We construct our test through an oracle approach. Suppose we know from an oracle $\bm{\mu}_{\mathcal{L}}$ and thus $\bm{\Sigma}_{\bm{\mu}_{\mathcal{L}}}$ as shown in Theorem 4.1. Note that by Theorem 4.1, 4.2 and 4.3 in Zhang (2019), under $H_{0,D}$ of (1.4), the distribution of $\bar{\bm{S}}_{\mathcal{L}}$ is either pairwise independent centered binomial or uncorrelated centered hypergeometric, depending on whether the marginal distributions are known. Therefore, for fixed $p$ and $D$, with a large $n$ and the central limit theorem on $\bar{\bm{S}}_{\mathcal{L}}=\bm{S}_{\mathcal{L}}/n$, we approximately have a simple-versus-simple hypothesis testing problem:

According to the fundamental Neyman-Pearson Lemma (Neyman and Pearson, 1933), the corresponding most powerful (MP) test is the likelihood ratio test. We thus consider the data-relevant part of the log-likelihood ratio of the above two distributions,

For a large $n$, the dominating term in $f_{\bar{\bm{S}}_{\mathcal{L}}}(\bm{\mu}_{\mathcal{L}})$ is $\bm{\mu}_{\mathcal{L}}^{T}(\bm{\Sigma}_{\bm{\mu}_{\mathcal{L}}}-\bm{\mu}_{\mathcal{L}}\bm{\mu}_{\mathcal{L}}^{T})^{-1}\bm{S}_{\mathcal{L}}.$ By (4.7) in Theorem 4.2, the first order Taylor expansion of this term is precisely $\bm{\mu}_{\mathcal{L}}^{T}\bm{S}_{\mathcal{L}}$! This implies that the MP test rejects when $\bar{\bm{S}}_{\mathcal{L}}$ is colinear with $\bm{\mu}_{\mathcal{L}}.$ The above heuristics thus suggests that we consider the oracle test statistic $B_{\text{oracle}}=\bm{\mu}_{\mathcal{L}}^{T}\bar{\bm{S}}_{\mathcal{L}}/\|\bm{\mu}_{\mathcal{L}}\|_{2}$.

In our simulation studies in Section 5, since we know the form of the alternative distribution, we can estimate $\bm{\mu}_{\mathcal{L}}$ with high accuracy through an independent simulation. That is, with the known alternative distribution $\mathbf{P}_{\bm{U}}$, for a large $K$ we simulate $\bm{V}_{1},\ldots,\bm{V}_{K}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathbf{P}_{\bm{U}}.$ From the binary expansion of $\bm{V}_{1},\ldots,\bm{V}_{K},$ we obtain the vector of symmetry statistics $\tilde{\bm{S}}_{\mathcal{L}}$ and an estimate of $\bm{\mu}_{\mathcal{L}}$ denoted by $\tilde{\bm{\mu}}_{\mathcal{L}}=\tilde{\bm{S}}_{\mathcal{L}}/n.$ The oracle test statistic from simulations is then $\widetilde{B}_{\text{oracle}}=\tilde{\bm{\mu}}_{\mathcal{L}}^{T}\bar{\bm{S}}_{\mathcal{L}}/\|\tilde{\bm{\mu}}_{\mathcal{L}}\|.$

We show in simulations that even when $D$ is as small as $3$, $\widetilde{B}_{\text{oracle}}$ is powerful, and numerically it outperforms all existing competitors under consideration across a wide spectrum of alternatives and noise levels. For example, for the cases when the joint distributions are Gaussian with linear dependency, the power curves of $\widetilde{B}_{\text{oracle}}$ dominate those of the distance correlation when $p=2$ and the F-test when $p=3,$ which are known to be optimal. Compared to existing tests, the gain of the BEAST with oracle in power suggests that suitably chosen deterministic weights for the alternative provide a unified yet simple solution to improve the power. To the best of our knowledge, this is the first time that such a benchmark on the feasible power performance is available for the problem of testing uniformity.

Besides the useful insight about the feasible limit of power, the oracle also provides insights on the optimal weights under each alternative. For example, in simulations we find high colinearity between the approximate oracle weight vector $\tilde{\bm{\mu}}_{\mathcal{L}}$ and that of the Spearman’s $\rho$, $\bm{r}_{D}$, as found in Section 3.1. This weight vector makes the one-sided test with $\widetilde{B}_{\text{oracle}}$ more powerful than the two-sided test with Spearman’s $\rho.$

Although the optimal weight $\bm{\mu}_{\mathcal{L}}$ or $\tilde{\bm{\mu}}_{\mathcal{L}}$ is unknown in practice, an unbiased and asymptotically efficient estimate of $\bm{\mu}_{\mathcal{L}}$ is $\bar{\bm{S}}_{\mathcal{L}}.$ This motivates us to develop an approximation of $\widetilde{B}_{\text{oracle}}$ through resampling and regularization, which we discuss next.

In practice, we are agnostic about $\bm{\mu}_{\mathcal{L}}.$ Blindly replacing $\tilde{\bm{\mu}}_{\mathcal{L}}$ in $\widetilde{B}_{\text{oracle}}$ with $\bar{\bm{S}}_{\mathcal{L}}$ will result in colinearity with itself and the statistic reduces to the classical $\chi^{2}$-test statistic. Traditionally, the data-splitting strategy has often been employed for this type of situations to facilitate data-driven decision (Hartigan, 1969; Cox, 1975), i.e., half of the data is used to calibrate the statistical procedure such as screening the null features (Wasserman and Roeder, 2009; Barber and Candès, 2019), determining the proper weights for individual hypotheses (Ignatiadis et al., 2016), recovering the optimal projection for dimension reducion (Huang, 2015), and estimating the latent loading for factor models (Fan et al., 2019), while a statistical decision is implemented using the remaining half. However, the single data-splitting procedure only uses half of the data for decision making, which inevitably bears undesirable randomness and therefore leads to power loss for hypothesis testing. Some recent efforts have shown that this shortcoming can be lessened by using multiple splittings (Romano and DiCiccio, 2019; Liu et al., 2019; Dai et al., 2020).

Motivated by the principle of multiple splitting, we propose to approximate $B_{\text{oracle}}$ through resampling: We replace $\bm{\mu}_{\mathcal{L}}$ in $B_{\text{oracle}}$ with $\bar{\bm{S}}_{\mathcal{L}}$, and we replace $\bar{\bm{S}}_{\mathcal{L}}$ in $B_{\text{oracle}}$ with its resampling version $\bar{\bm{S}}_{\mathcal{L}}^{*}$. Important resampling methods include bootstrap (Efron and Tibshirani, 1994) and subsampling (Politis et al., 1999). Bootstrap and subampling are known to have similar performance in approximating the sampling distribution of the target statistic. In this paper, we use the subsampling method to facilitate the calculation of the empirical copula distribution when the marginal distributions are unknown Nelsen (2007). In addition to the above consideration, one intuition behind this resampling approach is to help distinguish the alternative distribution from the null: Under the null, since $\bm{\mu}_{\mathcal{L}}={\bm{0}},$ we expect the magnitude of $\bar{\bm{S}}_{\mathcal{L}}$ and $\bar{\bm{S}}_{\mathcal{L}}^{*}$ to be small and not very colinear after regularization. On the other hand, under the alternative, since $\bm{\mu}_{\mathcal{L}}\neq{\bm{0}},$ we expect the two estimations of $\bm{\mu}_{\mathcal{L}}$ to be both colinear with $\bm{\mu}_{\mathcal{L}}$ and thus to be highly colinear themselves. Therefore, the magnitude of the test statistic could be different to help distinguish the alternative distribution from the null.

In addition, we apply regularization to accommodate sparsity, i.e., the non-uniformity can be explained by a few binary interactions $\Lambda$’s with $\mathbf{E}[A_{\Lambda}]\neq 0.$ This sparsity assumption is often reasonable in the BET framework, since $\mathbf{E}[A_{\Lambda}]=0$ is equivalent to the symmetry of distribution according to the interaction $\Lambda$. Thus, sparsity over $\mathbf{E}[A_{\Lambda}]$’s is equivalent to a highly symmetric distribution. For example, if a multivariate distribution is symmetric in every direction, then each one-dimensional projection of this distribution has a real characteristic function. By Theorem 2.2, we have $\mathbf{E}[A_{\Lambda}]=0$ for all $\Lambda$ involving an even number of binary variables (${\bm{1}}_{p}^{T}\Lambda{\bm{1}}_{D}$ is even). Many global forms of dependency also correspond to sparse structures in $\bm{\mu}_{\mathcal{L}}.$

The estimation of $\bm{\mu}_{\mathcal{L}}$ under the sparsity assumption is closely related to the normal mean problem, where many good regularization based methods are readily available (Wasserman, 2006). For example, in Donoho and Johnstone (1994), it is shown that estimation with soft thresholding is nearly optimal. We denote the vector-valued soft thresholding function by ${\mathcal{T}}(\bm{x},\lambda)$ for $q\times 1$ vector $\bm{x}$ and threshold $\lambda>0,$ so that ${\mathcal{T}}_{\ell}(\bm{x},\lambda)=\text{sign}(x_{\ell})(|x_{\ell}|-\lambda)_{+},~{}\ell=1,\ldots,q.$ In construction of our test statistic, we choose to use soft-thresholding as a regularization step to screen the small observations in $\bar{\bm{S}}_{\mathcal{L}}$ and $\bar{\bm{S}}_{\mathcal{L}}^{*}$ due to the null distribution or due to the sparsity $\mathbf{E}[A_{\Lambda}]=0$ for certain interaction $\Lambda$’s under the alternative, thus improves the power of the test statistic.

In summary, we consider the approximation of $B_{\text{oracle}}$ through subsampling, while using regularization to obtain a good estimate of the optimal weight vector ${\mathcal{T}}(\bar{\bm{S}}_{\mathcal{L}},\lambda)/\|{\mathcal{T}}(\bar{\bm{S}}_{\mathcal{L}},\lambda)\|_{2}.$ The detailed steps are listed below.

1. Step 1:

   From $n$ observations of $\bm{U}_{1},,\ldots,\bm{U}_{n},$ obtain $m$ subsamples of size $r$: $\bm{U}_{1,k}^{*},\ldots,\bm{U}_{r,k}^{*}$, $k=1,\ldots,m.$ For each subsample $k,$ base on the binary expansions of $\bm{U}_{1,k}^{*},\ldots,\bm{U}_{r,k}^{*}$, find the vector of average symmetry statistics $\bar{\bm{S}}_{\mathcal{L},k}^{*}$.

2. Step 2:

   Take the average over $m$ subsamples to obtain $\bar{\bm{S}}_{\mathcal{L}}^{*}={m}^{-1}\sum_{k=1}^{m}\bar{\bm{S}}_{\mathcal{L},k}^{*}.$ Apply the soft-thresholding function to get an estimate of $\bm{\mu}_{\mathcal{L}}$ as ${\mathcal{T}}(\bar{\bm{S}}_{\mathcal{L}}^{*},\lambda).$

3. Step 3:

   The BEAST statistic $B_{\lambda}$ is obtained as

   $$B_{\lambda}={\mathcal{T}}(\bar{\bm{S}}_{\mathcal{L}},\lambda)^{T}{\mathcal{T}}(\bar{\bm{S}}^{*}_{\mathcal{L}},\lambda)/\|{\mathcal{T}}(\bar{\bm{S}}_{\mathcal{L}},\lambda)\|_{2}.$$

   (4.8)