# Functional Response Designs via the Analytic Permutation Test

Adam B Kashlak and Sergii Myroshnychenko
Department of Mathematical and Statistical Sciences,
University of Alberta
and
Susanna Spektor
Pilon School of Business,
The Sheridan College Institute of Technology and Advanced Learning  The
authors gratefully acknowledge the Natural Sciences and Engineering Research
Council of Canada and the Pacific Institute for the Mathematical Sciences

###### Abstract

Vast literature on experimental design extends from Fisher and Snedecor to the
modern day. When data lies beyond the assumption of univariate normality,
nonparametric methods including rank based statistics and permutation tests
are enlisted. The permutation test is a versatile exact nonparametric
significance test that requires drastically fewer assumptions than similar
parametric tests. The main downfall of the permutation test is high
computational cost making this approach laborious for complex data and
sophisticated experimental designs and completely infeasible in any
application requiring speedy results such as high throughput streaming data.
We rectify this problem through application of concentration inequalities and
thus propose a computation free permutation test—i.e. a permutation-less
permutation test. This general framework is applied to multivariate, matrix-
valued, and functional data. We improve these concentration bounds via a novel
incomplete beta transform. We extend our theory from 2-sample to $k$-sample
testing through the use of weakly dependent Rademacher chaoses and modified
decoupling inequalities. We test this methodology on classic functional data
sets including the Berkeley growth curves and the phoneme dataset. We further
consider analysis of spoken vowel sound under two experimental designs: the
Latin square and the randomized block design.

Keywords: Concentration of Measure; Functional Data Analysis; Latin Squares
Kahane–Khintchine Inequalities; Randomized Block Design

##  1 Introduction

Exact significance tests date back to the very origins of statistical
hypothesis testing as an alternative to parametric testing. Namely, Fisher's
exact test tests for independence between the rows or columns of a $2\times 2$
contingency table by directly using the hypergeometric distribution instead of
relying on large sample asymptotic statistics such as the chi-squared test. As
a consequence, it obtains the exact p-value of the data without relying on
large sample asymptotics. However, Fisher's exact test is severely limited as
extension to general $r\times c$ tables requires significant amounts of
computational power to enumerate or approximate the entire discrete
distribution (Good, 1956; Agresti, 1992).

Permutation tests comprise a large subclass of such exact significance tests
and have been thoroughly studied (Mielke and Berry, 2007; Basso et al., 2009;
Pesarin and Salmaso, 2010; Brombin and Salmaso, 2013; Good, 2013). Given a
sample $X=\left\\{X_{1},\ldots,X_{n}\right\\}\in\mathcal{X}$ for some measure
space $\mathcal{X}$, a permutation test considers the finite sampling
distribution of a test statistic $T(X)$ over a discrete group where the
distribution of $T$ is invariant for any group action on the observed data
(Kallenberg, 2006)—i.e. for a group $G$, $T(gX)\overset{\text{d}}{=}T(X)$ for
any $g\in G$. A canonical example is one-way ANOVA; see Basso et al. (2009)
Section 5.2 for more details.

The permutation test requires far fewer assumptions than standard parametric
approaches—namely that of exchangability under the null hypothesis—and is thus
robust against deviations from distributional assumptions like normality and
provides guaranteed performance for finite samples. The main limitation is
that of computation. Performing a two sample permutation test for real valued
data is trivial with modern computers. What if we were to perform a $k$ sample
test with $k\choose 2$ post-hoc comparisons taking multiple testing into
account for, say, covariance operators as in Pigoli et al. (2014); Cabassi et
al. (2017) where every permutation requires computation of the singular value
decomposition (SVD) of a large matrix? Furthermore, what if we desire a more
sophisticated experimental design such as a randomized block, Latin square, or
unreplicated factorial design with the addition of multiple testing
corrections? The amount of computation required to get accurate p-values will
be prohibitive. The speech data and design considered in Section 5 would, for
example, require 264 SVDs per permutation and with 66 hypotheses to test at,
say, 2000 permutations each requires nearly 35 million SVDs. For matrices with
dimension $100\times 100$, this would take an estimated 36 hours on a Intel
Core i7-7567U CPU at 3.50GHz. For a $400\times 400$ matrix, it would take 74
days.

In this article, we present a unified methodology for performing computation-
free permutation tests for $k$ sample testing in commutative and non-
commutative $L^{q}$ spaces, which includes multivariate and functional data
and covariance matrices and operators among other data types. Specifically, we
consider the distribution of a test statistic on a discrete space of invariant
group actions. Instead of taking random draws from that space to get a
conditional Monte Carlo estimate (Hemerik and Goeman, 2018) of the p-value, we
apply recent extensions of the Kahane-Khintchine inequality for commutative
and non-commutative Banach spaces (Pisier and Xu, 2003; Garling, 2007;
Spektor, 2016) in order to achieve sub-Gaussian bounds on the tail probability
of our test statistic. Namely, we seek a result like
$\mathrm{P}\left(T(\pi)\geq T^{\star}\right)\leq\exp(-Ct^{2})$ for some
universal constant $C>0$ depending only on the space in which the data lives
irrespective of sample size and dimension. This methodology is presented in
Section 2 for two sample testing within commutative $L^{q}$ spaces—e.g.
univariate, vector valued, and functional data—as well as within non-
commutative $L^{q}$ spaces—e.g. covariance matrices and operators. As such
universal constants are often less than optimal for statistical use, we
introduce an adjustment for these upper bounds based on Talagrand's
concentration inequality (Talagrand, 1996) and the incomplete beta function in
Section 2.4. An extension to testing on $k$-samples is considered in Section 3
making use of Rademacher chaoses and decoupling inequalities (Kwapien, 1987;
De la Pena and Giné, 2012). Section 5 extends these ideas to multi-factor
designs for the analysis of spoken vowel phonemes.

Most previous work on fast or computation-free permutation testing focus on
univariate data in the setting of large scale testing typically applied to
testing for genomics data. The recent work of He et al. (2019) achieves this
goal by using Stolarsky's invariance principle. In Yang et al. (2019), ``very
small'' p-values are approximated via sequential Monte Carlo and the Edgeworth
expansion. In Segal et al. (2018), an asymptotic approximation and a clever
partitioning/resampling scheme is used to approximate small p-values. Density
approximation via Pearson curves (Solomon and Stephens, 1978) has recently
reemerged for p-value approximation in machine learning (Gretton et al., 2012)
and neuroimaging (Winkler et al., 2016) among other areas. While past work is
focused on large scale two-sample testing, this work is motivated by
$k$-sample tests and more sophisticated experimental designs with functional
and operator responses. While permutation tests have been used both for
pointwise and curve-wise analysis of functional data (Cox and Lee, 2008;
Corain et al., 2014; Chakraborty and Chaudhuri, 2015; Pigoli et al., 2014;
Pigoli et al., 2018; Cabassi et al., 2017), approaching statistical hypothesis
testing via analytic estimation of a permutation test p-value in general
Banach spaces has not been deeply explored as of yet.

As a proof of concept for testing within commutative and non-commutative
$L^{q}$ spaces, we consider a variety of simulated and real data sets in
Section 4. In Section 5, our bounds are applied to testing for phonological
differences among twelve spoken vowel sounds performed as a complete
randomized block design on covariance operators with respect to two binary
blocking factors: the speaker's country of origin
$\\{\text{Canada},\text{China}\\}$ and sex $\\{\text{male},\text{female}\\}$.
We also consider a Latin square design for checking the data for within
subject pronunciation changes while running the experiment. Section 5 contains
more detail on the data, experimental design, and its results. Proofs of the
main theorems, the necessary theoretical development, further data
experiments, and a discussion of past results are contained in the
supplementary material.

##  2 Two sample testing

###  2.1 Univariate data

Let $n=m_{1}+m_{2}$ and $X_{1},\ldots,X_{n}\in\mathbb{R}$ be independent
random variables such that $\mathrm{E}X_{i}=\mu_{1}$ for $i\leq m_{1}$ and
$\mathrm{E}X_{i}=\mu_{2}$ for $i\geq m_{1}+1$. We wish to test
$H_{0}:\mu_{1}=\mu_{2}$ versus $H_{1}:\mu_{1}\neq\mu_{2}$. To test these
hypotheses using a permutation test, we treat
$X_{1},\ldots,X_{n}\in\mathbb{R}$ as fixed and consider $\pi\in\mathbb{S}_{n}$
a random permutation uniformly distributed on the symmetric group on $n$
elements. That is, $\pi$ is a bijective map
$\pi:\\{1,\ldots,n\\}\rightarrow\\{1,\ldots,n\\}$. Thus, we can consider the
randomly permuted test statistic

$T(\pi)=\frac{1}{s_{n}}\left[\frac{1}{m_{1}}\sum_{i=1}^{m_{1}}X_{\pi(i)}-\frac{1}{m_{2}}\sum_{i=m_{1}+1}^{n}X_{\pi(i)}\right],$
(2.1)

which is normalized by the sample standard deviation $s_{n}$ for the entire
set $X_{1},\ldots,X_{n}$.111Note that $s_{n}$ is invariant under permutation
and is only included to make the below formulation nicer. The conditional tail
probability is

$\mathrm{P}\left(T(\pi)\geq
t\,|\,X_{1},\ldots,X_{n}\right)=\frac{1}{n!}\sum_{\pi\in\mathbb{S}_{n}}\bm{1}[{T(\pi)\geq
t}].$ (2.2)

Let $T_{0}$ be the test statistic $T(\pi)$ when $\pi$ is the identity—i.e. the
original ordering. Then, the p-value for the above hypothesis test is
$\mathrm{P}\left(T(\pi)\geq T_{0}\right)$, which is often approximated by
randomly generating $N\ll n!$ random permutations from $\mathbb{S}_{n}$
instead of exhaustively enumerating all elements of $\mathbb{S}_{n}$. This
results in an overly conservative test for p-values approaching $1/N$.

To avoid the simulation-based approximation of equation 2.2, we instead prove
a sub-Gaussian bound on the p-value.

######  Theorem 2.1 (Univariate Data).

For $T(\pi)$ from equation 2.1 with $m_{1}=\kappa m_{2}$ for some $\kappa\geq
1$, then $\mathrm{P}\left(T(\pi)\geq
t\right)\leq\exp\left(-nt^{2}/2\lceil\kappa+1\rceil^{3}\right).$

This theorem is extended to more advanced settings, including vectors,
matrices, functional data, and operators, in the following sections. We state
those subsequent theorems for balanced samples, but note that the imbalanced
setting of Theorem 2.1 can also be incorporated with a similar constant
$\kappa$. The proofs are more tedious for imbalanced data, but no additional
innovation is required.

###  2.2 Commutative $L^{q}$ Spaces

To extend our tail bounds beyond the real valued setting, we require some
definitions. Note that both of the following definitions extend to the case of
compact operators on Banach spaces.

######  Definition 2.1 (Matrix Square Root).

Let $A\in\mathbb{R}^{d\times d}$ with $d\geq 2$ be a symmetric positive semi-
definite matrix with eigen-decomposition $A=UD{U}^{\mathrm{T}}$ where
$U=(v_{1}~{}v_{2}~{}\ldots~{}v_{d})$ is the orthonormal matrix of eigenvectors
and $D$ is the diagonal matrix of eigenvalues,
$(\lambda_{1},\ldots,\lambda_{d})$. Then, $A^{1/2}=UD^{1/2}{U}^{\mathrm{T}}$
where $D^{1/2}$ is the diagonal matrix with entries
$(\lambda_{1}^{1/2},\ldots,\lambda_{d}^{1/2})$.

######  Definition 2.2 ($q$-Schatten norm for matrices).

For an arbitrary matrix $A\in\mathbb{R}^{k\times l}$ and $q\in(1,\infty)$, the
$q$-Schatten norm is $\left\lVert
A\right\rVert_{S^{q}}^{q}=\mathrm{tr}[({A}^{\mathrm{T}}A)^{q/2}]=\lVert{\bm{\nu}}\rVert_{\ell^{q}}^{q}=\sum_{i=1}^{\min\\{k,l\\}}\nu_{i}^{q}$
where ${\bm{\nu}}=(\nu_{1},\ldots,\nu_{\min\\{k,l\\}})$ is the vector of
singular values of $A$ and where $\lVert\cdot\rVert_{\ell^{q}}$ is the
standard $\ell^{q}$ norm in $\mathbb{R}^{d}$. In the covariance matrix case
where $A\in\mathbb{R}^{d\times d}$ is symmetric and positive-definite,
$\left\lVert
A\right\rVert_{S^{q}}^{q}=\mathrm{tr}\left(A^{q}\right)=\lVert{\bm{\lambda}}\rVert_{\ell^{q}}^{q}$
where ${\bm{\lambda}}$ is the vector of eigenvalues of $A$.

When $q=\infty$, we have the standard operator norm on
$\ell^{2}(\mathbb{R}^{d})$, $\lVert
A\rVert_{S^{\infty}}=\sup_{v\in\mathbb{R}^{d},\lVert
v\rVert_{\ell^{2}}=1}\lVert
Av\rVert_{\ell^{2}}=\sup_{v\in\mathbb{R}^{d},\lVert
v\rVert_{\ell^{2}}=1}{v}^{\mathrm{T}}Av.$ In the covariance matrix setting,
this coincides with the maximal eigenvalue of $A$.

Let $X_{1},\ldots,X_{n}\in\mathcal{X}$ where
$\\{\mathcal{X},\lVert\cdot\rVert\\}$ is a commutative $L^{q}$ space. The test
statistic of interest is $T_{0}=\lVert\sum_{i\leq
m}X_{i}-\sum_{i>m}X_{i}\rVert$. Then, Theorem 2.1 can be extended to such
settings using a version of the Kahane-Khintchine inequality under a weak
dependency condition from Theorem A.7 proved in the supplementary material.
For simplicity of notation, we assume that the $X_{i}$ are centred about the
sample mean and that the samples are balanced.

######  Theorem 2.2 (Commutative $L^{q}$ Spaces).

Let $m=n/2$, $\lVert\cdot\rVert_{S^{q}}$ be the $q$-Schatten norm for matrices
or operators, and $\varepsilon_{1},\ldots,\varepsilon_{n}$ be Rademacher
random variables such that $\sum_{i=1}^{n}\varepsilon_{i}=0$. Let
$X_{1}(t),\ldots,X_{n}(t)$ be continuous function on a compact interval with
empirical covariance operator
$\hat{\Sigma}(s,t)=(n-1)^{-1}\sum_{i=1}^{n}X_{i}(s)X_{i}(t)$. Let
$q\in[1,\infty)$ with norm $\lVert\cdot\rVert_{L^{q}}$. For
$T(\pi)=\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert_{L^{q}}$.
Then, $\mathrm{P}\left(T(\pi)\geq
t\right)\leq\exp\left(-t^{2}/c\lVert\hat{\Sigma}^{1/2}\rVert_{S^{q}}^{2}\right).$

######  Remark 2.3 (On optimal constants).

The optimal constant $c$ in the above theorem follows from the optimal
constant in the Kahane-Khintchine inequality, which is not currently
known.222It took about 60 years from the advent of the original Khintchine
inequality for optimal constants to be determined. However, it is strongly
conjectured to agree with the optimal constant for the standard Khintchine
inequality. In that case, we would take $c=64$ in the above theorem, which is
$16$ from Theorem 2.1 times $4$ from that fact that $T(\pi)$ is not a
symmetric random variable. For more details, see the proof and discussion in
the supplementary material. We also empirically adjust the p-values in Section
2.4, which is demonstrated to give strong performance in Sections 4 and 5.

###  2.3 Non-Commutative $L^{q}$ Spaces

Following from the previous section, we outline similar tail bounds in non-
commutative $L^{q}$ spaces (Pisier and Xu, 2003). This methodology encompasses
matrix and operator data with emphasis on application to testing for equality
of covariances. Hence, the following theorem is applied to symmetric positive
definite operators in the example below and to the data in Section 5. The test
statistic of interest is still $T_{0}=\lVert\sum_{i\leq
m}X_{i}-\sum_{i>m}X_{i}\rVert$, but with the $X_{i}$ now belonging to a non-
commutative $L^{q}$ space.

######  Theorem 2.3 (Non-Commutative $L^{q}$ Spaces).

Let $\lVert\cdot\rVert_{S^{q}}$ be the $q$-Schatten norm for a matrix or
operator and $\varepsilon_{1},\ldots,\varepsilon_{n}$ be Rademacher random
variables such that $\sum_{i=1}^{n}\varepsilon_{i}=0$. For $d,d^{\prime}>1$,
let $X_{1},\ldots,X_{n}\in\mathbb{R}^{d\times d^{\prime}}$ be a collection of
$n$ fixed matrices (or let $X_{1},\ldots,X_{n}$ be a collection of bounded
linear operators). For
$T(\pi)=\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert_{S^{q}}$,
there exists a universal constant $c>0$ such that
$\mathrm{P}\left(T(\pi)>t\right)\leq\exp\left(-t^{2}/c\mathcal{S}^{2}\right)$
where
$\mathcal{S}=\max\left\\{\lVert(\sum_{i=1}^{n}X_{i}{X}_{i}^{*})^{1/2}\rVert_{S^{q}},\lVert(\sum_{i=1}^{n}{X}_{i}^{*}{X_{i}})^{1/2}\rVert_{S^{q}}\right\\}$
with $X_{i}^{*}$ the adjoint operator.

######  Remark 2.4.

Of particular interest are covariance operators being compact trace-class
self-adjoint operators. Consequently, we have the same bound but with
$\mathcal{S}=\lVert(\sum_{i=1}^{n}X_{i}^{2})^{1/2}\rVert_{S^{q}}.$

###  2.4 Beta and Empirical Beta Adjustment

Inequalities such as the Kahane-Khintchine inequalities are useful tools for
considering the finite sample performance of a statistical method. However,
the biggest impediment to the use of such inequalities, as well as other
concentration inequalities, for statistical inference is the nearly inevitable
loss in power to reject the null due to `universal constants' that are too
large for application. We thus propose a transformation based on the beta
distribution to correct the p-values and recover the lost statistical power.
Proposition 2.5 only applies to univariate data and requires an asymptotic
arguement outside of our finite sample focus. It is included nevertheless to
set the stage for the non-asymptotic beta transform in Theorem 2.4.
Furthermore, this simpler setting yields explicit beta parameters and
demonstrates stellar performance in both the simulated data of Section 4.1 and
the extremely imbalanced small sample setting explored in Kashlak and Yuan
(2020) where, for example, $m_{1}=335$ and $m_{2}=3$.

For a statistical test, if the correct test size is achieved, then a random
null p-value will be distributed as $\mathrm{Uniform}\left[0,1\right]$.
However, our Kahane-Khintchine based null p-values will instead closely follow
a more general $\mathrm{Beta}\left(\alpha,\beta\right)$ distribution. Thus,
identification of the parameters $\alpha$ and $\beta$ will allow us to adjust
the p-values to the null setting to recover lost statistical power. This idea
is spiritually similar to the Pearson curve method (Solomon and Stephens,
1978), but that approach requires estimation of the first 4 central moments
for comparison with the family of generalized Pearson distributions compared
to our more focused use of the beta distribution with Theorem 2.4 proved to
justify such focus. Usage of the Edgeworth expansion (Hall, 2013) is another
method with a long history, but requires some care to note whether or not a
finite number of terms in such an expansion can provide a satisfactory
approximation to the probability density in question (Stuart et al., 1994).

We first consider the univariate case of Section 2.1 before discussing the
more general Banach space setting for the beta transform.

######  Proposition 2.5.

Under the setting of Theorem 2.1 with $n$ sufficiently large,

$\mathrm{P}\left(\exp\left\\{-nT(\pi)^{2}/2\lceil\kappa+1\rceil^{3}\right\\}<u\right)\leq
C_{0}I\left(u;\frac{\lceil\kappa+1\rceil^{3}}{(2+\kappa+\kappa^{-1})},\frac{1}{2}\right)$

where $I(u;\alpha,\beta)$ is the regularized incomplete beta function and
$C_{0}={\left(\frac{\lceil\kappa+1\rceil^{3}}{2+\kappa+\kappa^{-1}}\right)^{1/2}\Gamma\left(\frac{\lceil\kappa+1\rceil^{3}}{2+\kappa+\kappa^{-1}}\right)}{\Gamma\left(\frac{1}{2}+\frac{\lceil\kappa+1\rceil^{3}}{2+\kappa+\kappa^{-1}}\right)^{-1}}.$

Proposition 2.5 allows us to adjust the p-values from Theorem 2.1 so that our
test statistic achieves the desired empirical size. The refined bound is

$\mathrm{P}\left(T(\pi)>t\right)\leq
C_{0}I\left(\mathrm{e}^{-nt^{2}/2\lceil\kappa+1\rceil^{3}};\frac{\lceil\kappa+1\rceil^{3}}{2+\kappa+\kappa^{-1}},\frac{1}{2}\right)$

This adjustment is shown to work in the simulations detailed in Figure 1. For
the more general Banach space setting, we can use Talagrand's concentration
inequality (Talagrand, 1996) to prove the following theorem.

######  Theorem 2.4.

Let $(\mathcal{X},\lVert\cdot\rVert)$ be a Banach space with separable dual
space $\mathcal{X}^{*}$, and let $h:\mathbb{R}\rightarrow\mathbb{R}$ be
monotonically increasing. For any random variable $X$ taking values in
$\mathcal{X}$ such that $\mathrm{E}h(\lVert X\rVert)^{2}<\infty$ and
$\sup_{X\in\mathcal{X}}h(\lVert X\rVert)<U<\infty$ and for $u\in(0,1)$ and
some constants $C,c,\alpha,\beta>0$, $\mathrm{P}\left(\mathrm{e}^{-h(\lVert
X\rVert)/c}<u\right)\leq CI(u;\alpha,\beta)$ where $I(u;\alpha,\beta)$ is the
incomplete beta function for $c$ sufficiently large.

######  Remark 2.6.

Theorem 2.4 requires the Banach space $\mathcal{X}$ to have a separable dual.
This stems from writing the norm as a countable supremum for use within
Talagrand's concentration inequality (Talagrand, 1996). We can directly apply
this result to commutative and non-commutative $L^{q}$ spaces for
$1<q<\infty$. However, $L^{\infty}$ is a standard example of a non-separable
Banach space. For our purposes, we can avoid this issue as it is typical in
functional data analysis to consider the uniform norm on the space of
continuous bounded functions with compact support.

When working in commutative and non-commutative $L^{q}$ spaces, we no longer
have easily defined constants for the righthand bound in Theorem 2.4. Hence,
we instead propose an empirical beta transform, outlined in Algorithm 1, which
estimates these constants. To do this, we must choose a small number $r$ of
permutations to draw uniformly at random from $\mathbb{S}_{n}$. In practice,
we find that $r=10$ or 20 is sufficient to achieve good results on real data.
From these, we compute test statistics sampled from the null setting, which
will yield a collection of $r$ p-values. These p-values can in turn be used to
estimate the parameters for a beta distribution via the method of moments
estimate $\hat{\alpha}$ and $\hat{\beta}$. Lastly, the p-value $p_{0}$
produced by $T_{0}$ can be adjusted by application of the incomplete beta
function: $I(p_{0};\hat{\alpha},\hat{\beta})$. This method was applied to most
of the data examples detailed in Section 4. This transform is also shown to
work well for the construction of wild bootstrap confidence regions for least
squares and ridge regression (Burak and Kashlak, 2021).

Algorithm 1 The Empirical Beta Transform  Compute p-value $p_{0}$ from test
statistic $T_{0}$ using Theorem 2.2 or 2.3.
---
Choose $r>1$, the number of permutations to simulate—e.g. $r=10$.
Draw $\pi_{1},\ldots,\pi_{r}$ from $\mathbb{S}_{n}$ uniformly at random.
Compute p-values $p_{1},\ldots,p_{r}$ from test statistics
$T_{\pi_{1}},\ldots,T_{\pi_{r}}$.
Find the method of moments estimator for $\alpha$ and $\beta$.
Estimate first and second central moments of the $p_{i}$ by $\bar{p}$ and
$s^{2}$.
Estimate $\hat{\alpha}={\bar{p}^{2}(1-\bar{p})}/{s^{2}}-\bar{p}$.
Estimate $\hat{\beta}=[\bar{p}(1-\bar{p})/s^{2}-1][1-\bar{p}]$.
Return the adjusted p-value $I(p_{0};\hat{\alpha},\hat{\beta})$.

##  3 k sample testing

For general one-way ANOVA and more complex experimental designs, we extend the
above two sample tests to $k$ level factors. The two challenges to overcome
are (1) proper multiple testing correction for the $k\choose 2$ pairwise
comparisons and (2) the construction of a global p-value. Classical hypothesis
testing would have us first reject the global hypothesis and follow up with
pairwise post-hoc testing. For permutation tests, we begin with pairwise
testing and combine these tests into a global p-value.

For one-way ANOVA, let $X_{i,j}$ be the $j$th observation from category $i$
for $i=1,\ldots,k$ and $j=1,\ldots,n_{i}$ under the model

$X_{i,j}=\mu+\tau_{i}+\varepsilon_{i,j}$ (3.1)

with global mean $\mu$, $i$th treatment effect $\tau_{i}$ with
$\sum_{i=1}^{k}\tau_{i}=0$, and exchangeable errors $\varepsilon_{i,j}$—i.e.
permutationally invariant (Kallenberg, 2006). We wish to test the following:

Pairwise $\displaystyle H_{0}^{(ij)}:\tau_{i}=\tau_{j}$ $\displaystyle
H_{1}^{(i,j)}:\tau_{i}\neq\tau_{j}$ Global $\displaystyle
H_{0}:\tau_{1}=\ldots=\tau_{k}=0$ $\displaystyle H_{1}:\exists\tau_{i}\neq 0.$

Under the pairwise null $H_{0}^{i,j}$, the difference in category means is
$\bar{X}_{i\cdot}-\bar{X}_{j\cdot}=\bar{\varepsilon}_{i\cdot}-\bar{\varepsilon}_{j\cdot}$.
Thus, the permutation test requires exchangeable errors—i.e. the distribution
of $\bar{\varepsilon}_{i\cdot}-\bar{\varepsilon}_{j\cdot}$ is invariant under
any random permutation. This is weaker than the standard iid setup and, most
critically, does not require normality.

###  3.1 Multiple Pairwise Tests

From Section 2, we can compute test statistics $T_{0}^{(ij)}$ for
$H_{0}^{(ij)}$ and consider the permutation distribution of $T^{(ij)}(\pi)$
for some uniformly distributed $\pi\in\mathbb{S}_{n_{i}+n_{j}}$. For
familywise type I error control, the pairwise statistics come from independent
applications of dependent Rademacher vectors. Hence, we can rely on standard
multiple testing corrections such as the simple Bonferroni correction as
proposed in (Basso et al., 2009, Chapter 5) or the slightly more involved
step-down procedure used in Cabassi et al. (2017). Other methods include
Holm's stepdown method (Holm, 1979), the approach outlined in the canonical
text Westfall and Young (1993), and the more recent Romano and Wolf (2005). In
experimental design, some authors even prefer to forego such corrections and
report raw uncorrected p-values (Wu and Hamada, 2011). The focus of this
article is on computation of the raw p-values and hence, application of one's
favourite multiple testing correction is left to the reader. For the
phonological data analysis in Section 5, we will just consider the raw
p-values and the Bonferroni correction.

###  3.2 Global Test

The k-sample global significance test statistic can be written as a
combination of the pairwise statistics:

$T_{0}=\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}n_{i}n_{j}(T_{0}^{(ij)})^{2}$ (3.2)

To test the significance of $T_{0}$, a permutation framework can be
implemented in one of three ways; see Basso et al. (2009) Chapter 5 for more
details. The first is the pooled method in which the entire data set of
$N=n_{1}+\ldots+n_{k}$ observations is permuted. The second is by aggregation
of the pairwise statistics where each permutation is applied independently to
each pair of samples. The third is the synchronized method which only applies
to balanced designs—i.e. $n_{1}=\ldots=n_{k}$—in which the same permutations
are applied to each category pairing $(i,j)$. This is the preferable approach
when the design is balanced (Basso et al., 2009). As we have already discussed
individual pairwise testing, we focus on the synchronized test in the context
of our Kahane–Khintchine methodology.

######  Remark 3.1.

Beyond univariate data, the above test statistic $T_{0}$ can be considered on
the direct sum of $\kappa={k\choose 2}$ Banach spaces. That is, for a sequence
of Banach spaces $(\mathcal{X}_{i},\lVert\cdot\rVert_{i})$ and elements
$X_{i}\in\mathcal{X}_{i}$, we can define a new Banach space by the $\ell^{2}$
direct sum
$(X_{i})_{i=1}^{n}\in\left(\bigoplus_{i=1}^{\kappa}\mathcal{X}_{i}\right)_{\ell^{2}}$
with norm $\lVert(X_{i})_{i=1}^{n}\rVert=(\sum_{i=1}^{n}\lVert
X_{i}\rVert_{i}^{2})^{1/2}$. See any text on discussing sequences in Banach
spaces such as Diestel et al. (1995) for more details.

The synchronized setting is the preferred approach for balanced designs; see,
for example, Basso et al. (2009); Cabassi et al. (2017). This approach applies
the same permutations to each pairing. Let $X^{(1)},\ldots,X^{(k)}$ be
$m$-long column vectors containing the observations of samples $1,\ldots,k$,
respectively. Then, let $X$ be the $2m\times{k\choose 2}$ matrix with columns

$X=\begin{pmatrix}X^{(1)}&X^{(1)}&\ldots&X^{(k-1)}\\\
X^{(2)}&X^{(3)}&\ldots&X^{(k)}\end{pmatrix}.$

Lastly, let
${\varepsilon}^{\mathrm{T}}=(\varepsilon_{1},\ldots,\varepsilon_{2m})$ such
that $\sum\varepsilon_{i}=0$. The synchronized permuted version of the global
test statistic in Equation 3.2 is then
$\lVert{X}^{\mathrm{T}}\varepsilon\rVert_{\ell^{2}}^{2}=\sum_{i,j=1}^{2m}a_{i,j}\varepsilon_{i}\varepsilon_{j}$
for $a_{i,j}$, the $i,j$th entry in $X{X}^{\mathrm{T}}$. This is a second
order Rademacher chaos (Ledoux and Talagrand, 1991, Section 4.4) except that
the $\varepsilon_{i}$ are not iid. In this case, we still have a sub-Gaussian
bound achievable via a decoupling argument (Kwapien, 1987) with proof in the
supplementary material. See De la Pena and Giné (2012) for more on decoupling
inequalities.

######  Theorem 3.1.

Let $T=\lVert{X}^{\mathrm{T}}\varepsilon\rVert_{\ell^{2}}$ for $X$ the above
$2m\times{k\choose 2}$ matrix and $\varepsilon_{i}$ such that
$\sum_{i=1}^{2m}\varepsilon_{i}=0$. Then, for some universal constant $c$,
$\mathrm{P}\left(T>t\right)\leq\exp\left[-t^{2}/c\mathcal{S}\right]$ where
$\mathcal{S}=\lVert X{X}^{\mathrm{T}}\rVert_{S^{2}}$.

######  Remark 3.2.

Up to constant $c$, this theorem coincides with the result for a two sample
test as for
$X=(x_{1}^{(1)},\ldots,x_{m}^{(1)},x_{1}^{(2)},\ldots,x_{m}^{(2)})^{\text{T}}$,
the term $\mathcal{S}$ equals the sample variance of the $x_{i}^{(j)}$.
However, the constant $c$ emerging from the proof is very large. This
universal constant problem is rectified via the empirical beta transform
presented in Section 2.4.

##  4 Data Examples

###  4.1 Univariate Data

####  4.1.1 Two Sample Test

The performance of Theorem 2.1 on simulated data is displayed in Figure 1 for
balanced and for imbalanced samples averaged over 1000 replications. In the
balanced case, we simulate $m_{1}=m_{2}=100$ Gaussian random variates with
distributions $\mathcal{N}\left(0,1\right)$ and
$\mathcal{N}\left(\mu,1\right)$ for $\mu\in[0,1]$. We compare the classic
student's t-test to the permutation test with 1000 permutations, the bounds
from Theorem 2.1 with $\kappa=1$, and the beta adjusted bound from Proposition
2.5. Notably, the balanced Khintchine bound returns p-values just slightly
larger than the standard t-test while the beta adjusted bound is even tighter.
For the imbalanced case, the sample sizes are now $m_{1}=140,m_{2}=60$ and
$\kappa=2.33$. The imbalanced bound is not as sharp, but the beta adjusted
bound still gives a close approximation to the t-test p-value.

Figure 1:  Univariate two sample test for normal data with balanced sample
sizes $m_{1}=m_{2}=100$ (left) and for imbalanced $m_{1}=140,m_{2}=60$ (right)
comparing the standard t-test (black) to the permutation test (red) with 1000
permutations and to the Khintchine bound, Theorem 2.1, (green) and the beta
adjusted bound, Proposition 2.5, (blue) all across 1000 replications.

####  4.1.2 K Sample Test

The performance of Theorem 3.1 for comparing $k$ samples of size $n$ via a
synchronized permutation test is demonstrated in Figure 2. For this
simulation, $k=4,16$ for the left and right plot, respectively, samples of
size $n=20$ were generated as random Gaussian variates with variance 1 and
with mean 0 for the first $k-1$ sets and with mean $\mu\in[0,2]$ for the $k$th
set. As $\mu$ grows, the p-value for the standard F-test, the synchronized
permutation test, and the beta-adjusted p-value from Theorem 3.1 all decrease
in tandem for $k=4$ with the unadjusted bound above the others. In the $k=16$
case, the beta adjusted bound and the synchronized permutation test return the
same p-values until the lines approach the permutation boundary at
$-\log_{2}(1001)$. More notably, they slightly differ from the classic F-test
as for relatively large $k$ and small $n$ the synchronized permutation test
returns marginally different p-values than the F-test. A total of 1000 random
permutations were generated for the synchronized permutation test, and this
simulation was replicated 1000 times to create these plots.

Figure 2:  Univariate $k$ sample test for normal data with $k=4$ (left) and
$k=16$ (right) balanced samples of size $n=20$. The figure compares the
standard F-test (black) to the synchronized permutation test (red) with 1000
permutations and to the unadjusted (green) and beta adjusted (blue) bounds
from Theorem 3.1.

###  4.2 Berkeley Growth Curves: Functional Means

To demonstrate Theorem 2.2, we apply it to the classic Berkeley growth curve
dataset (Ramsay and Silverman, 2005).333 This data is available in the R
package fda (Ramsay et al., 2018).  This dataset contains measurements of 93
children—39 males and 54 females—taken at 31 time points between the ages of 1
and 18 years. A set of 30 curves was randomly selected from the male curves
and 30 curves from the female curves to test for a difference in the
population mean curves based on those observations. This was repeated 100
times to see the resulting p-values under the $L^{1}$, $L^{2}$, and
$L^{\infty}$ norms. Table 1 displays the percentage of rejections. Applying
Theorem 2.2 results in a reasonable number of rejections under the $L^{1}$
topology. However, differences are not detectable in $L^{2}$ or $L^{\infty}$.
This is rectified via the empirical beta adjustment.

Table 1:  Displayed above are the percentages of rejections by using Theorem 2.2 (left) and the beta adjusted bound (right) at test sizes $5\%,1\%,0.1\%$ for the $L^{1}$, $L^{2}$, and $L^{\infty}$ norms.  | Percentage of Rejections
---|---
| Kahane Bound | Beta Adjusted
Size | 5% | 1% | 0.1% | 5% | 1% | 0.1%
$L^{1}$ | 86% | 42% | 7% | 85% | 55% | 31%
$L^{2}$ | 0% | 0% | 0% | 100% | 88% | 77%
$L^{\infty}$ | 0% | 0% | 0% | 100% | 100% | 98%

###  4.3 Phoneme Curves: Covariance Operators

We apply Theorem 2.3 to the classic phoneme dataset (Ferraty and Vieu, 2006),
which consists of 400 log-periodograms for 5 different phonemes—the vowel from
`dark' aa, the vowel from `water' ao, the plosive d-sound dcl, the fricative
sh-sound sh, the vowel from she iy—sampled at 150 frequencies.444 This data is
available in the R package fds (Shang and Hyndman, 2013).  Using the notation
of the International Phonetic Alphabet (IPA), aa is A, ao is O, dcl is d, sh
is S, and iy is i. To produce covariance operators for testing, we first
randomly permute the order of the 400 curves, then group these curves into
sets of 10 to produce a set of 40 covariance operators for each of the five
phoneme classes. This is replicated 100 times with different random groupings
of curves.

We apply our method after using the empirical beta adjustment from Section 2.4
to each of the 10 pairwise comparisons between phonemes resulting in Table 2.
In the trace norm topology, all $100\times 10$ pairwise tests result in
rejection for a test size of 1%. The Hilbert-Schmidt norm only detects a
significant difference between A and O about 52% of the time whereas the
operator norm fails to detect any significance between those two phonemes. The
difference between phonemes A and O is hardest to identify among the 10
pairings.

Table 2:  The percentage of rejected two sample tests at the 1% level comparing two different phonemes with a sample size of $m_{1}=m_{2}=10$ under the trace, Hilbert-Schmidt, and operator norms.  | Trace Norm |  | Hilbert-Schmidt Norm |  | Operator Norm
---|---|---|---|---|---
| A | O | d | S |  | A | O | d | S |  | A | O | d | S
O | 100 |  |  |  |  | 52 |  |  |  |  | 0 |  |  |
d | 100 | 100 |  |  |  | 93 | 86 |  |  |  | 15 | 23 |  |
S | 100 | 100 | 100 |  |  | 99 | 100 | 100 |  |  | 88 | 97 | 91 |
i | 100 | 100 | 100 | 100 |  | 100 | 100 | 100 | 100 |  | 100 | 100 | 100 | 100

##  5 Phonological differences between vowels

Taking inspiration from the classic phoneme dataset (Hastie et al., 1995;
Ferraty and Vieu, 2006) discussed previously in Section 4.3, we consider a new
data set of log-periodograms for the phonemes of 12 spoken vowels detailed in
Table 3.555Note that this data was collected outside of a proper laboratory
setting to be a proof-of-concept for the proposed methodology as opposed to an
in depth study of language. This data is available at
https://sites.ualberta.ca/~kashlak/kashData.html.

The raw data consists of 12 phonemes recorded 12 times each from 4 different
speakers. The data was recorded on a Tascam DR-05 portable linear PCM audio
recorder as a mono 24-bit wave file sampled at 96 kHz, which is currently
considered high definition audio in contrast to the standard 16-bit 44.1 kHz
audio on compact discs. The primary vowel phoneme was extracted as a 170
millisecond clip corresponding to $16384=2^{14}$ samples. These clips were
transformed into log periodograms via the tuneR package (Ligges et al., 2018)
as displayed in Figure 3 for a single speaker. As is common with functional
data, the raw log-periodograms were first smoothed. In this case, cubic
smoothing splines were used. However, many other smoothing methods can be and
have been applied to functional data.

Figure 3:  The log-periodograms of all 12 vowel phonemes spoken 12 times by
one of the speakers considered over the first 1000 frequencies.

Two experimental designs were employed in the collection of this data and will
be tested in the following subsections. First, the 12 words were vocalized 12
times in a Latin square design. Each row corresponds to a replication of
speaking all of the 12 words, and each column corresponds to the order of the
words within a replication. This was done to test for changes in speech during
the recording period. Secondly, this Latin square design was replicated for 4
different speakers with two binary blocking factors male/female and
Canadian/Chinese. Thus, we have a $12\times 4$ complete randomized block
design with functional responses. The total sample size is $576=12\times
12\times 4$ log-periodogram curves.

Table 3:  The twelve vowel phonemes considered in our dataset along with the 12 spoken words used to produce those vowels.  i | weave | e | waiter | E | wetter | æ | wag
---|---|---|---|---|---|---|---
I | wiggle | 9 | what | u | woohoo | U | woods
3 | world | o | woke | 2 | wonder | 6 | water
Figure 4:  IPA Vowel Chart,
http://www.internationalphoneticassociation.org/content/ipa-chart, available
under a Creative Commons Attribution-Sharealike 3.0 Unported License.
Copyright 2015 International Phonetic Association.

###  5.1 Latin square design for functional means

For an unreplicated Latin square design, we cannot perform a permutation test
for the significance of each factor simultaneously. Exchangability under the
null hypothesis for one factor requires fixing the levels of all other factors
when permuting labels. However, if we fix the Latin square row and column
indices then only a single observation remains leaving nothing to permute. To
rectify this, a stepdown approach as in Basso et al. (2009) chapter 7 for
unreplicated factorial designs can be applied. As a permutation test requires
exchangeable observations under the null hypothesis, test statistics for each
factor are first computed. Beginning with the largest, if that null hypothesis
holds, then this implies that all other null hypotheses hold and hence acts as
the global null allowing for the data to be permuted. If this null is
rejected, then we proceed to test the second largest test statistic while
fixing the levels of the first factor. Once a null is not rejected, this
method stops. Otherwise, all factors can be tested except for the last one as
rejecting all other null hypotheses would leave no room for further
permutations.

For the vowel data, we have 12-level row, column, and vowel factors giving the
model
$y_{ijk}(t)=\mu(t)+\text{row}_{i}(t)+\text{column}_{j}(t)+\text{vowel}_{k}(t)+\xi_{ijk}$
where $y_{ijk}(t)$ is a smoothed log-periodogram, $\mu$ is the global mean,
and the $\xi_{ijk}$ are mean zero exchangeable errors. For all four subjects,
the vowel factor produced a much larger test statistic than the row and column
effects as expected indicating consistency of the speaker during the
experiment. Thus, after rejecting the null of there being no difference among
the spoken vowels, the row or column factor can be considered. For all four
subjects, the row and column effects were not deemed to be statistically
significant—i.e. there were no detectable changes in pronunciation across the
recording session. Pairwise comparison of the vowels for each subject was also
performed. However, of the 66 pairwise hypotheses to test, one subject
rejected 25 nulls, another rejected only 5 nulls, and the last two rejected 0
nulls after taking multiple testing into account. This is in contrast to the
randomized block design discussed in the next section that, making use of the
entire dataset and covariance operators, identifies 60 of the 66 pairings as
significantly different.

Before computing the test statistics and p-values in the randomized block
design discussed next, each log-periodogram was centred by subtracting off the
row and column effects from the Latin square design. This resulted in an
improvement in the reported p-values, which were larger in the case that the
row and column effects were not removed.

###  5.2 Complete Randomized block design for functional data

A complete randomized block design (CRBD) aims to test a treatment effect as
in one-way ANOVA but with the addition of blocking factors to account for
sources of variation unrelated to the treatment of interest. For functional
data, a CRBD can be performed by using the synchronized permutation test for
two-way ANOVA from chapter 6 of Basso et al. (2009) combined with the Kahane-
Khintchine based tail bound. To achieve this, a difference between the
functional means or covariances is computed for each of the $12\choose 2$
vowel pairings while holding the levels of the blocking factors constant. For
each pairing, the test statistics can be summed over the levels of the
blocking factors thus removing any influence from interaction terms even
though they are generally assumed to be negligible in this setting. Theorems
2.2 and 2.3 can be applied for functional means and covariance operators
respectively to bound the pairwise p-values. The computed test statistics can
be aggregated using Theorem 3.1 to get a global p-value. Note that a standard
permutation test would require the computation of $264=66\times 4$ test
statistics via simulation from the symmetric group, which in the case of
covariance operators and Schatten norms implies 264 SVD calculations per
permutation. This is further expanded by, say, performing $132,000=66\times
2000$ permutations to be able to test each hypothesis at the 0.01 level after
correcting for multiple testing. Focusing only on the approximately 35 million
required SVDs, a timing test run on an Intel Core i7-7567U CPU at 3.50GHz
estimates 36 hours of compute time when considering $100$ dimensional matrices
and an estimated run time of 74 days on $400\times 400$ dimensional matrices.

This approach was applied pairwise to the sample covariance operators for each
vowel as past work has emphasized that the covariance structure of speech data
is the best lens to detect phonological differences (Pigoli et al., 2014;
Pigoli et al., 2018). Application of Theorem 2.3 using the trace norm
(1-Schatten norm) and using the empirical beta adjustment from Section 2.4
produced the 66 pairwise p-values displayed in Figure 5. The words are also
grouped by p-value to display which vowel phonemes proved statistically
indistinguishable using our proposed methodology. The use of other Schatten
norms results in lower power—i.e. fewer null hypotheses rejected.

The blocking factors {male,female} and {Canadian, Chinese} can also be
similarly tested without removing the row and column effects from the Latin
square design; otherwise, the mean taken over the entire dataset will be zero.
In trace norm, we get p-values of 0.0002 and 0.00003 for sex and country,
respectively. In Hilbert-Schmidt norm, we get the weaker p-values 0.03 and
0.07.

Figure 5:  $\log_{10}$ p-values for pairwise two sample tests between vowel
pairs under the trace norm.

##  6 Discussion

The p-value has stood for over a century as a pillar of frequentist
statistical methodology. In this article, we approached $k$-sample testing
through application of an analytic approximation to the permutation test
p-value notably without relying on simulation of the permutation distribution
of the test statistic. Experimental design for functional data was the main
motivation for this work as standard simulation-based permutation testing can
be applied but at a high computational cost. Other applications of interest
include online testing where data must be processed, results returned, and
decisions made in real time. The lag resulting from a classic permutation test
is unacceptable in such settings. This methodology is generally applicable to
other complex testing settings including other types of group invariances—e.g.
rotationally invariant test statistics. Furthermore, the duality of hypothesis
testing with confidence sets suggests investigation into using variants of the
Kahane–Khintchine inequality to construct confidence balls for estimators with
finite sample guarantees on the coverage.

## Supplementary Material

Primarily, the supplementary material contains proofs of the main results as
well as auxilary theorems and lemmas. Secondly, it contains further exposition
of this methodology on simulated and real data sets. In particular, these
methods are shown to produce correct p-values in the null setting.

## References

  * Agresti (1992) Agresti, A. (1992).  A survey of exact inference for contingency tables.  Statistical science 7(1), 131–153.
  * Basso et al. (2009) Basso, D., F. Pesarin, L. Salmaso, and A. Solari (2009).  Permutation tests for stochastic ordering and ANOVA: theory and applications with R, Volume 194.  Springer Science & Business Media.
  * Boucheron et al. (2013) Boucheron, S., G. Lugosi, and P. Massart (2013).  Concentration inequalities: A nonasymptotic theory of independence.  Oxford University Press.
  * Bousquet (2003) Bousquet, O. (2003).  Concentration inequalities for sub-additive functions using the entropy method.  In Stochastic inequalities and applications, pp. 213–247. Springer.
  * Brombin and Salmaso (2013) Brombin, C. and L. Salmaso (2013).  Permutation tests in shape analysis, Volume 15.  Springer.
  * Burak and Kashlak (2021) Burak, K. L. and A. B. Kashlak (2021).  Nonparametric confidence regions via the analytic wild bootstrap.  (under review).
  * Cabassi et al. (2017) Cabassi, A., D. Pigoli, P. Secchi, and P. A. Carter (2017).  Permutation tests for the equality of covariance operators of functional data with applications to evolutionary biology.  Electronic Journal of Statistics 11(2), 3815–3840.
  * Chakraborty and Chaudhuri (2015) Chakraborty, A. and P. Chaudhuri (2015).  A wilcoxon–mann–whitney-type test for infinite-dimensional data.  Biometrika 102(1), 239–246.
  * Corain et al. (2014) Corain, L., V. B. Melas, A. Pepelyshev, and L. Salmaso (2014).  New insights on permutation approach for hypothesis testing on functional data.  Advances in Data Analysis and Classification 8(3), 339–356.
  * Cox and Lee (2008) Cox, D. D. and J. S. Lee (2008).  Pointwise testing with functional data using the westfall–young randomization method.  Biometrika 95(3), 621–634.
  * De la Pena and Giné (2012) De la Pena, V. and E. Giné (2012).  Decoupling: from dependence to independence.  Springer Science & Business Media.
  * Diestel et al. (1995) Diestel, J., H. Jarchow, and A. Tonge (1995).  Absolutely summing operators, Volume 43.  Cambridge university press.
  * Ferraty and Vieu (2006) Ferraty, F. and P. Vieu (2006).  Nonparametric functional data analysis: theory and practice.  Springer Science & Business Media.
  * Garling (2007) Garling, D. J. (2007).  Inequalities: a journey into linear analysis.  Cambridge University Press.
  * Giné and Nickl (2016) Giné, E. and R. Nickl (2016).  Mathematical Foundations of Infinite-Dimensional Statistical Models.  Cambridge University Press.
  * Good (1956) Good, I. (1956).  On the estimation of small frequencies in contingency tables.  Journal of the Royal Statistical Society: Series B (Methodological) 18(1), 113–124.
  * Good (2013) Good, P. (2013).  Permutation tests: a practical guide to resampling methods for testing hypotheses.  Springer Science & Business Media.
  * Gretton et al. (2012) Gretton, A., K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola (2012).  A kernel two-sample test.  Journal of Machine Learning Research 13(Mar), 723–773.
  * Haagerup (1981) Haagerup, U. (1981).  The best constants in the Khintchine inequality.  Studia Mathematica 70, 231–283.
  * Hall (2013) Hall, P. (2013).  The bootstrap and Edgeworth expansion.  Springer Science & Business Media.
  * Hastie et al. (1995) Hastie, T., A. Buja, and R. Tibshirani (1995).  Penalized discriminant analysis.  The Annals of Statistics, 73–102.
  * He et al. (2019) He, H. Y., K. Basu, Q. Zhao, and A. B. Owen (2019).  Permutation $p$-value approximation via generalized stolarsky invariance.  The Annals of Statistics 47(1), 583–611.
  * Hemerik and Goeman (2018) Hemerik, J. and J. Goeman (2018).  Exact testing with random permutations.  Test 27(4), 811–825.
  * Holm (1979) Holm, S. (1979).  A simple sequentially rejective multiple test procedure.  Scandinavian journal of statistics, 65–70.
  * Kahane (1964) Kahane, J.-P. (1964).  Sur les sommes vectorielles sigma plus minus un.  Comptes rendus hebdomadaires des seances de l'academie des sciences 259(16), 2577.
  * Kallenberg (2006) Kallenberg, O. (2006).  Probabilistic symmetries and invariance principles.  Springer Science & Business Media.
  * Kashlak and Yuan (2020) Kashlak, A. B. and W. Yuan (2020).  Computation-free nonparametric testing for local and global spatial autocorrelation with application to the canadian electorate.  arXiv preprint arXiv:2012.08647.
  * Klein and Rio (2005) Klein, T. and E. Rio (2005).  Concentration around the mean for maxima of empirical processes.  The Annals of Probability 33(3), 1060–1077.
  * Kwapien (1987) Kwapien, S. (1987).  Decoupling inequalities for polynomial chaos.  The Annals of Probability 15(3), 1062–1071.
  * Latała and Oleszkiewicz (1994) Latała, R. and K. Oleszkiewicz (1994).  On the best constant in the Khintchine–Kahane inequality.  Studia Mathematica 109(1), 101–104.
  * Ledoux and Talagrand (1991) Ledoux, M. and M. Talagrand (1991).  Probability in Banach Spaces: isoperimetry and processes, Volume 23.  Springer.
  * Ligges et al. (2018) Ligges, U., S. Krey, O. Mersmann, and S. Schnackenberg (2018).  tuneR: Analysis of Music and Speech.
  * Mielke and Berry (2007) Mielke, P. W. and K. J. Berry (2007).  Permutation methods: a distance function approach.  Springer Science & Business Media.
  * Pesarin and Salmaso (2010) Pesarin, F. and L. Salmaso (2010).  Permutation tests for complex data: theory, applications and software.  John Wiley & Sons.
  * Pigoli et al. (2014) Pigoli, D., J. A. Aston, I. L. Dryden, and P. Secchi (2014).  Distances and inference for covariance operators.  Biometrika, asu008.
  * Pigoli et al. (2018) Pigoli, D., P. Z. Hadjipantelis, J. S. Coleman, and J. A. Aston (2018).  The statistical analysis of acoustic phonetic data: exploring differences between spoken romance languages.  Journal of the Royal Statistical Society: Series C (Applied Statistics) 67(5), 1103–1145.
  * Pisier and Xu (2003) Pisier, G. and Q. Xu (2003).  Non-commutative lp-spaces.  Handbook of the geometry of Banach spaces 2, 1459–1517.
  * Ramsay and Silverman (2005) Ramsay, J. O. and B. W. Silverman (2005).  Functional data analysis.  New York: Springer.
  * Ramsay et al. (2018) Ramsay, J. O., H. Wickham, S. Graves, and G. Hooker (2018).  fda: Functional Data Analysis.  R package version 2.4.8.
  * Romano and Wolf (2005) Romano, J. P. and M. Wolf (2005).  Exact and approximate stepdown methods for multiple hypothesis testing.  Journal of the American Statistical Association 100(469), 94–108.
  * Segal et al. (2018) Segal, B. D., T. Braun, M. R. Elliott, and H. Jiang (2018).  Fast approximation of small p-values in permutation tests by partitioning the permutations.  Biometrics 74(1), 196–206.
  * Shang and Hyndman (2013) Shang, H. L. and R. J. Hyndman (2013).  fds: Functional data sets.  R package version 1.7.
  * Solomon and Stephens (1978) Solomon, H. and M. A. Stephens (1978).  Approximations to density functions using pearson curves.  Journal of the American Statistical Association 73(361), 153–160.
  * Spektor (2014) Spektor, S. (2014).  Selected Topics in Asymptotic Geometric Analysis and Approximation Theory.  Ph. D. thesis, University of Alberta.
  * Spektor (2016) Spektor, S. (2016).  Restricted Khinchine inequality.  Canadian Mathematical Bulletin 59(1), 204–210.
  * Stuart et al. (1994) Stuart, A., S. Arnold, J. K. Ord, A. O'Hagan, and J. Forster (1994).  Kendall's advanced theory of statistics.  Wiley.
  * Talagrand (1996) Talagrand, M. (1996).  New concentration inequalities in product spaces.  Inventiones mathematicae 126(3), 505–563.
  * Watson (1959) Watson, G. (1959).  A note on gamma functions.  Edinburgh Mathematical Notes 42, 7–9.
  * Westfall and Young (1993) Westfall, P. H. and S. S. Young (1993).  Resampling-based multiple testing: Examples and methods for p-value adjustment, Volume 279.  John Wiley & Sons.
  * Winkler et al. (2016) Winkler, A. M., G. R. Ridgway, G. Douaud, T. E. Nichols, and S. M. Smith (2016).  Faster permutation inference in brain imaging.  Neuroimage 141, 502–516.
  * Wu and Hamada (2011) Wu, C. J. and M. S. Hamada (2011).  Experiments: planning, analysis, and optimization, Volume 552.  John Wiley & Sons.
  * Yang et al. (2019) Yang, J. J., E. M. Trucco, and A. Buu (2019).  A hybrid method of the sequential monte carlo and the edgeworth expansion for computation of very small p-values in permutation tests.  Statistical methods in medical research 28(10-11), 2937–2951.

##  Appendix A Inequalities

###  A.1 Khintchine-type Inequalities

######  Theorem A.1 (Khintchine's Inequality (1923)).

For any $p\in(0,\infty)$, there exist positive finite constants $A_{p}$ and
$B_{p}$ such that for any sequence $x_{1},\ldots,x_{n}\in\mathbb{R}$ (or
$x_{i}\in\mathbb{C}$),

$A_{p}^{p}\lVert
x\rVert_{\ell^{2}}^{p}\leq\mathrm{E}\left\lvert\sum_{i=1}^{n}\varepsilon_{i}x_{i}\right\rvert^{p}\leq
B_{p}^{p}\lVert x\rVert_{\ell^{2}}^{p}$

where $\varepsilon_{1},\ldots,\varepsilon_{n}$ are iid Rademacher random
variables—i.e.
$\mathrm{P}\left(\varepsilon_{i}=1\right)=\mathrm{P}\left(\varepsilon_{i}=-1\right)=1/2$.

For this article, we are only concerned with the upper bound $B_{2p}$ for
$p>2$. In Garling (2007), $B_{2p}=[(2p)!/2^{p}p!]^{1/2p}$ which gives
$B_{2p}<\sqrt{2p}$, but also via Stirling's inequality $B_{p}\sim(p/e)^{1/2}$
as $p\rightarrow\infty$. The expectation in above theorem is with respect to
the $\varepsilon_{i}$ corresponding to a uniform distribution on the $2^{n}$
vertices of the $n$-hypercube. In what follows, we consider expectation over
the uniform distribution on the $n!$ elements of the symmetric group
$\mathbb{S}_{n}$. This will be denoted $\mathrm{E}_{\pi}$ where
$\pi\in\mathbb{S}_{n}$ is treated as a uniform random permutation.

In Spektor (2016), the restricted Khintchine inequality is introduced where it
is required that $\sum_{i=1}^{n}\varepsilon_{i}=0$ introducing a weak
dependency among the $\varepsilon_{i}$. In the proof in Spektor (2016), this
weak dependency doubles the variance by comparing two sets of data. Thus, the
constant becomes $B_{2p}=[(2p)!/p!]^{1/2p}$.

######  Theorem A.2 (Spektor (2016) Theorem 1.1).

For any $p\in[2,\infty)$, there exist positive finite constant $B_{p}$ such
that for any sequence $x_{1},\ldots,x_{n}\in\mathbb{R}$,

$\mathrm{E}\left\lvert\sum_{i=1}^{n}\varepsilon_{i}x_{i}\right\rvert^{p}\leq
B_{p}^{p}\left(\lVert
x\rVert_{\ell^{2}}^{2}-n\bar{x}^{2}\right)^{p/2}=B_{p}^{p}[(n-1)s_{n}^{2}]^{p/2}$
(A.1)

where $\varepsilon_{1},\ldots,\varepsilon_{n}$ are Rademacher random variables
such that $\sum\varepsilon_{i}=0$ and
$s_{n}^{2}=(n-1)^{-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}$ is the sample variance
of $x$.

######  Remark A.3.

In the statistics context, if we divide Inequality A.1 by $m=n/2$, we have
$\mathrm{E}_{\pi}\lvert\bar{x}^{(\pi)}_{1}-\bar{x}^{(\pi)}_{2}\rvert^{p}\leq
B_{p}(2s_{n}^{2}/m)^{p/2}$ where $\bar{x}^{(\pi)}_{1}$ is the average of the
first $m$ of the $x_{\pi(i)}$ for some random permutation $\pi$ and similarly
for $\bar{x}^{(\pi)}_{2}$.

The previous theorem only applies to a balanced two sample setting. In the
following, we extend the ideas in Spektor (2016) to the imbalanced testing
setting. Other such extensions to imbalanced Khintchine inequalities were
considered in Spektor (2014). Note that in the following theorem, the bound on
the right-hand-side is in terms of the smaller of the two sample sizes
$m_{2}<m_{1}$ reducing the power drastically in a highly imbalanced setting.
Nevertheless, it still is seen to be an excellent statistical tool in a
variety of applied settings after the beta correction is applied (Kashlak and
Yuan, 2020).

######  Theorem A.4 (Imbalanced Case).

For $m_{1}>m_{2}>0$, let $n=m_{1}+m_{2}$ and $M=m_{1}-m_{2}$ and let $\kappa
m_{2}=m_{1}$ for some rational $\kappa>1$. Let $\delta_{1},\ldots,\delta_{n}$
be weighted dependent Rademacher random variables such that marginally
$\mathrm{P}\left(\delta_{i}=1/m_{1}\right)=\mathrm{P}\left(\delta_{i}=-1/m_{2}\right)=1/2$
and such that $\sum\delta_{i}=0$—i.e. precisely $m_{1}$ of the $\delta_{i}$
equal $1/m_{1}$ and $m_{2}$ equal $-1/m_{2}$. For any $p\in[2,\infty)$, there
exists a positive finite constant $B_{p}$ such that for any sequence
$x_{1},\ldots,x_{n}\in\mathbb{R}$,666This theorem is also valid for
$x_{i}\in\mathbb{C}$ after standard alterations are made in the proof.

$\mathrm{E}\left\lvert\sum_{i=1}^{n}\delta_{i}x_{i}\right\rvert^{p}\leq
B_{p}\left(\frac{\lceil\kappa+1\rceil^{2}{s}_{n}^{2}}{2m_{2}}\right)^{p/2}$

where $s_{n}^{2}=(n-1)^{-1}\sum_{i=1}^{n}(x_{i}-\bar{x})^{2}$ is the sample
variance of $x$.

######  Lemma A.5.

Let $\xi_{1},\xi_{2}\in[1,\infty]$ be such that $\xi_{1}^{-1}+\xi_{2}^{-1}=1$,
and let $X,Y$ be positive real random variables. Then,

$\min_{\xi_{1},\xi_{2}:\,\xi_{1}^{-1}+\xi_{2}^{-1}=1}\mathrm{E}\left\\{\xi_{1}^{p-1}X^{p}+\xi_{2}^{p-1}Y^{p}\right\\}=\left[(\mathrm{E}X^{p})^{1/p}+(\mathrm{E}Y^{p})^{1/p}\right]^{p}.$

###### Proof.

We note that $\xi_{2}=\xi_{1}/(\xi_{1}-1)$. Then,

$\displaystyle 0$
$\displaystyle=\frac{d}{d\xi_{1}}\left\\{\xi_{1}^{p-1}\mathrm{E}X^{p}+\xi_{2}^{p-1}\mathrm{E}Y^{p}\right\\}$
$\displaystyle=(p-1)\xi_{1}^{p-2}\mathrm{E}X^{p}-(p-1)\xi_{1}^{p-2}\mathrm{E}Y^{p}/(\xi_{1}-1)^{p}$
$\displaystyle\xi_{1}$
$\displaystyle=1+\left(\mathrm{E}Y^{p}/\mathrm{E}X^{p}\right)^{1/p}$
$\displaystyle\xi_{2}$
$\displaystyle=1+\left(\mathrm{E}X^{p}/\mathrm{E}Y^{p}\right)^{1/p}$

Hence,

$\displaystyle\min_{\xi_{1},\xi_{2}:\,\xi_{1}^{-1}+\xi_{2}^{-1}=1}\mathrm{E}\left\\{\xi_{1}^{p-1}X^{p}+\xi_{2}^{p-1}Y^{p}\right\\}$
$\displaystyle=\left[1+\left(\frac{\mathrm{E}Y^{p}}{\mathrm{E}X^{p}}\right)^{1/p}\right]^{p-1}\mathrm{E}X^{p}+\left[1+\left(\frac{\mathrm{E}X^{p}}{\mathrm{E}Y^{p}}\right)^{1/p}\right]^{p-1}\mathrm{E}Y^{p}$
$\displaystyle=\left[(\mathrm{E}Y^{p})^{1/p}+(\mathrm{E}X^{p})^{1/p}\right]^{p-1}(\mathrm{E}X^{p})^{1/p}+\left[(\mathrm{E}X^{p})^{1/p}+(\mathrm{E}Y^{p})^{1/p}\right]^{p-1}(\mathrm{E}Y^{p})^{1/p}$
$\displaystyle=[(\mathrm{E}X^{p})^{1/p}+(\mathrm{E}Y^{p})^{1/p}]^{p}$

∎

###### Proof of Theorem A.4.

We first decompose the weighted Rademacher sum. Without loss of generality,
assume $m_{1}>m_{2}$ and let $n=m_{1}+m_{2}$ and $M=m_{1}-m_{2}$. Also, assume
the $x_{i}$ are centred—i.e. $\sum_{i=1}^{n}x_{i}=0$—and let
$\xi_{1},\xi_{2}>0$ such that $\xi_{1}^{-1}+\xi_{2}^{-1}=1$. Thus, via
convexity, we have

$\displaystyle\mathrm{E}$
$\displaystyle\left\lvert\sum_{i=1}^{n}\delta_{i}x_{i}\right\rvert^{p}$
$\displaystyle=\mathrm{E}_{\pi}\left\lvert\frac{1}{m_{2}}\sum_{i=1}^{m_{2}}x_{\pi(i)}-\frac{1}{m_{1}}\sum_{i=m_{2}+1}^{n}x_{\pi(i)}\right\rvert^{p}$
$\displaystyle=\mathrm{E}_{\pi}\left\lvert\frac{1}{m_{2}}\left\\{\sum_{i=1}^{m_{2}}x_{\pi(i)}-\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\\}-\frac{1}{m_{1}}\sum_{i=2m_{2}+1}^{n}x_{\pi(i)}+\frac{M}{m_{1}m_{2}}\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\rvert^{p}$
$\displaystyle\leq\frac{\xi_{1}^{p-1}}{m_{2}^{p}}\mathrm{E}_{\pi}\left\lvert\sum_{i=1}^{m_{2}}x_{\pi(i)}-\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\rvert^{p}+\frac{\xi_{2}^{p-1}M^{p}}{m_{1}^{p}}\mathrm{E}_{\pi}\left\lvert\frac{1}{M}\sum_{i=2m_{2}+1}^{n}x_{\pi(i)}-\frac{1}{m_{2}}\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\rvert^{p}$
$\displaystyle=\frac{\xi_{1}^{p-1}}{m_{2}^{p}}(\mathrm{I})+\frac{\xi_{2}^{p-1}M^{p}}{m_{1}^{p}}(\mathrm{II}).$

To bound $(\mathrm{I})$, we apply the balanced weakly dependent Khintchine
inequality. Let $I\subset\\{1,\ldots,n\\}$ with cardinality $\lvert
I\rvert=M$. For such an index set $I$, let
$\Pi_{I}=\\{\pi\in\mathbb{S}_{n}\,:\,\pi(\\{2m_{2}+1,\ldots,n\\})=I\\}$. That
is, $\pi\in\Pi_{I}$ maps the final $M$ indices into $I$. Note that
$\lvert\Pi_{I}\rvert={n\choose M}$. As a result,

$\displaystyle\mathrm{E}_{\pi}\left\lvert\sum_{i=1}^{m_{2}}x_{\pi(i)}-\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\rvert^{p}$
$\displaystyle\leq\frac{1}{n!}\sum_{\pi\in\mathbb{S}_{n}}\left\lvert\sum_{i=1}^{m_{2}}x_{\pi(i)}-\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\rvert^{p}$
$\displaystyle\leq\frac{(n-M)!}{n!}\sum_{\lvert
I\rvert=M}\frac{1}{(n-M)!}\sum_{\pi\in\Pi_{I}}\left\lvert\sum_{i=1}^{m_{2}}x_{\pi(i)}-\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\rvert^{p}$
$\displaystyle\leq\frac{M!(n-M)!}{n!}\sum_{\lvert
I\rvert=M}\left[B_{p}(2m_{2}-1)^{p/2}(\sum_{i\notin I}x_{i}^{2})^{p/2}\right]$
$\displaystyle\leq B_{p}(2m_{2}-1)^{p/2}{s}_{n}^{p}.$

As the $x_{i}$ are centred, we have that $(\sum_{i\notin
I}x_{i}^{2})^{p/2}\leq s_{n}^{p}=(\sum x_{i}^{2})^{p/2}$, and hence

$\frac{\xi_{1}^{p-1}}{m_{2}^{p}}\mathrm{E}_{\pi}\left\lvert\sum_{i=1}^{m_{2}}x_{\pi(i)}-\sum_{i=m_{2}+1}^{2m_{2}}x_{\pi(i)}\right\rvert^{p}\leq\xi_{1}^{p-1}B_{p}(2m_{2}^{-1})^{p/2}{s}_{n}^{p}.$

For $(\mathrm{II})$, we first assume that $\kappa$ is a positive integer and
$m_{1}=\kappa m_{2}$ so $M=(\kappa-1)m_{2}$. In this case, we have

$\frac{\xi_{2}^{p-1}M^{p}}{m_{1}^{p}}(\mathrm{II})=\xi_{2}^{p-1}\left(\frac{\kappa-1}{\kappa}\right)^{p}\mathrm{E}_{\pi}\left\lvert\sum_{i=m_{2}+1}^{n}\tilde{\delta}_{i}x_{i}\right\rvert^{p}$

where $\tilde{\delta}_{i}$ are weighted Rademacher random variables with
taking values $1/M$ or $-1/m_{2}$ such that $\sum\tilde{\delta}_{i}=0$.
Applying Lemma A.5 gives

$\mathrm{E}_{\pi}\left\lvert\sum_{i=1}^{n}{\delta}_{i}x_{i}\right\rvert^{p}\leq\left\\{B_{p}^{1/p}\left(\frac{2{s}_{n}^{2}}{m_{2}}\right)^{1/2}+\left(\frac{\kappa-1}{\kappa}\right)\left(\mathrm{E}_{\pi}\left\lvert\sum_{i=m_{2}+1}^{n}\tilde{\delta}_{i}x_{i}\right\rvert^{p}\right)^{1/p}\right\\}^{p}.$

Noting that
$\mathrm{E}_{\pi}\left\lvert\sum_{i=m_{2}+1}^{n}\tilde{\delta}_{i}x_{i}\right\rvert^{p}$
is merely the original term to be bounded but with $m_{1}=\kappa m_{2}$ and
$M=(\kappa-1)m_{2}$ replaced by $(\kappa-1)m_{2}$ and $(\kappa-2)m_{2}$,
respectively, we apply this idea $\kappa-1$ more times to get

$\displaystyle\mathrm{E}_{\pi}\left\lvert\sum_{i=1}^{n}{\delta}_{i}x_{i}\right\rvert^{p}$
$\displaystyle\leq
B_{p}\left(\frac{2{s}_{n}^{2}}{m_{2}}\right)^{p/2}\left\\{1+\left(\frac{\kappa-1}{\kappa}\right)\left(1+\left(\frac{\kappa-2}{\kappa-1}\right)\left(\cdots(1+1/2)\cdots\right)\right)\right\\}^{p}$
$\displaystyle\leq
B_{p}\left(\frac{2{s}_{n}^{2}}{m_{2}}\right)^{p/2}\left\\{1+\frac{\kappa-1}{\kappa}+\frac{\kappa-2}{\kappa}+\ldots+\frac{1}{\kappa}\right\\}^{p}$
$\displaystyle\leq
B_{p}\left(\frac{{s}_{n}^{2}}{m_{2}}\frac{(\kappa+1)^{2}}{2}\right)^{p/2}.$

Noting that $n=m_{1}+m_{2}=(\kappa+1)m_{2}$, we have

$\mathrm{E}_{\pi}\left\lvert\sum_{i=1}^{n}{\delta}_{i}x_{i}\right\rvert^{p}\leq
B_{p}{s}_{n}^{p}\left(\frac{(\kappa+1)^{3}}{2n}\right)^{p/2}.$

Now, consider $\kappa=a+r\in\mathbb{Q}$ with $a\in\mathbb{N}$ and $r\in[0,1)$.
Then,

$\displaystyle B_{p}$
$\displaystyle\left(\frac{2{s}_{n}^{2}}{m_{2}}\right)^{p/2}\left\\{1+\frac{\kappa-1}{\kappa}+\frac{\kappa-2}{\kappa}+\ldots+\frac{r}{\kappa}\right\\}^{p}$
$\displaystyle\leq
B_{p}\left(\frac{2{s}_{n}^{2}}{m_{2}}\right)^{p/2}\left\\{\frac{a}{\kappa}+\frac{a-1}{\kappa}+\frac{a-2}{\kappa}+\ldots+\frac{1}{\kappa}+\frac{(a+1)r}{\kappa}\right\\}^{p}$
$\displaystyle\leq
B_{p}\left(\frac{2{s}_{n}^{2}}{m_{2}}\right)^{p/2}\frac{1}{\kappa^{p}}\left\\{\frac{a(a+1)}{2}+(a+1)r\right\\}^{p}$
$\displaystyle\leq
B_{p}\left(\frac{{s}_{n}^{2}}{m_{2}}\frac{(a+1)^{2}}{2}\right)^{p/2}\left\\{\frac{a+2r}{a+r}\right\\}^{p}$
$\displaystyle=B_{p}\left(\frac{{s}_{n}^{2}}{m_{2}}\frac{(a+1)^{2}}{2}\right)^{p/2}\left\\{1+\frac{r}{\kappa}\right\\}^{p}.$

Noting further that $ra+r<a+r$ so that $1+r/(a+r)<1+1/(a+1)$, we multiply by
$(a+1)$ on each side to get $(a+1)(1+r/(a+r))<a+2$. Hence,

$B_{p}\left(\frac{2{s}_{n}^{2}}{m_{2}}\right)^{p/2}\left\\{1+\frac{\kappa-1}{\kappa}+\frac{\kappa-2}{\kappa}+\ldots+\frac{r}{\kappa}\right\\}^{p}\leq
B_{p}\left(\frac{{s}_{n}^{2}}{m_{2}}\frac{(a+2)^{2}}{2}\right)^{p/2}.$

Hence, for $\kappa\in\mathbb{N}$, we have $\kappa+1=\lceil\kappa+1\rceil$, and
for $\kappa$ a non-integer we have
$a+2=\lfloor\kappa\rfloor+2=\lceil\kappa\rceil+1=\lceil\kappa+1\rceil$. ∎

###  A.2 Kahane-Khintchine-type Inequalities

Kahane extended Khintchine's inequality from the real line to normed spaces
Kahane (1964); Latała and Oleszkiewicz (1994). The optimal value for the
constant $C_{p,p^{\prime}}$ in Theorem A.6 below is not known in the case of
interest for this article, $p>p^{\prime}=2$; however, it has been conjectured
to be the same as in the real case, and as we see from the simulations and
real data experiments, this conjecture seems to hold for our purposes. In what
follows, let $\mathcal{X}$ be a normed space with norm $\lVert\cdot\rVert$.
Those spaces of statistical interest include $\mathbb{R}^{d}$, $L^{2}(0,1)$,
and spaces of matrices and positive definite trace class operators—i.e.
covariance operators.

######  Theorem A.6 (Kahane-Khintchine Inequality (1964)).

For any $p,p^{\prime}\in[1,\infty)$, there exists a universal finite constant
$C_{p,p^{\prime}}>0$ such that for any sequence of
$X_{1},\ldots,X_{n}\in\mathcal{X}$

$\left\\{\mathrm{E}\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert^{p}\right\\}^{1/p}\leq
C_{p,p^{\prime}}\left\\{\mathrm{E}\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert^{p^{\prime}}\right\\}^{1/p^{\prime}}$

where $\varepsilon_{i}$ are iid Rademacher random variables.

In general, we will consider the right hand side with $p^{\prime}=2$, which
bounds the $p$th moments by the second moment. For statistical applications,
we are interested in a few specific setting for this theorem. Namely, if
$\mathcal{X}=\mathbb{R}^{d}$ for $d\geq 2$, then for the $\ell^{q}$ norm with
$q\in[1,\infty]$, we have

$\left\\{\mathrm{E}\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert_{\ell^{q}}^{p}\right\\}^{1/p}\leq
C_{p}\left\lVert\left(\sum_{i=1}^{n}X_{i}{X}^{\mathrm{T}}_{i}\right)^{1/2}\right\rVert_{S^{q}}=C_{p}(n-1)^{1/2}\left\lVert\hat{\Sigma}^{1/2}\right\rVert_{S^{q}}$

where $\lVert\cdot\rVert_{S^{q}}$ is the $q$-Schatten norm and $\hat{\Sigma}$
is the empirical covariance estimator for the $X_{i}$. Similarly, in the
functional data setting, if $X_{i}$ are continuous and in $L^{q}[0,1]$, then
the right hand side becomes
$C_{p}(n-1)^{1/2}\lVert\hat{\Sigma}(s,s^{\prime})^{1/2}\rVert_{S^{q}}$ where
$\hat{\Sigma}:[0,1]^{2}\rightarrow\mathbb{R}$ is the empirical covariance
operator.

For non-commutative Banach spaces (Pisier and Xu, 2003), such as when $X_{i}$
are real valued matrices, we have a slightly different bound. Let
$\mathcal{X}=\mathbb{R}^{d\times d^{\prime}}$. Then, with respect to the
$q$-Schatten norm,

$\left\\{\mathrm{E}\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert_{S^{q}}^{p}\right\\}^{1/p}\leq
C_{p}\max\left\\{\left\lVert\left(\sum_{i=1}^{n}X_{i}{X}^{\mathrm{T}}_{i}\right)^{1/2}\right\rVert_{S^{q}},\left\lVert\left(\sum_{i=1}^{n}{X}^{\mathrm{T}}_{i}{X_{i}}\right)^{1/2}\right\rVert_{S^{q}}\right\\}.$

The above results all have iid $\varepsilon_{i}$. Applying similar methods as
in Spektor (2016) and as in the previous section, we can consider the moment
bounds under weak dependency conditions on the $\varepsilon_{i}$. This theorem
is stated for balanced samples with adjustments for imbalanced samples omitted
as they follow exactly as in the previously discussed real valued setting.

######  Theorem A.7 (Kahane-Khintchine with Weak Dependence).

Let $\varepsilon_{i}$ are Rademacher random variables such that
$\sum\varepsilon_{i}=0$. Furthermore, let $p\in[1,\infty)$.

For commutative Banach spaces there exists a universal finite constant
$C_{p}>0$ such that for any sequence of $X_{1},\ldots,X_{n}\in\mathcal{X}$

$\left\\{\mathrm{E}\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert^{p}\right\\}^{1/p}\leq
C_{p}2^{1/2}\left\lVert\left(\sum_{i=1}^{n}X_{i}X_{i}^{\star}\right)^{1/2}\right\rVert.$

For non-commutative Banach spaces there exists a universal finite constant
$C_{p}>0$ such that for any sequence of $X_{1},\ldots,X_{n}\in\mathcal{X}$

$\left\\{\mathrm{E}\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert^{p}\right\\}^{1/p}\leq
C_{p}2^{1/2}\max\left\\{\left\lVert\left(\sum_{i=1}^{n}X_{i}X_{i}^{\star}\right)^{1/2}\right\rVert,\left\lVert\left(\sum_{i=1}^{n}X_{i}^{\star}X_{i}\right)^{1/2}\right\rVert\right\\}.$

Before proving this theorem, we discuss some preliminary results regarding
Schatten norms. Let $\preceq$ denote positive semi-definite ordering. For
positive semi-definite $q$-Schatten class linear operators $\Gamma$ and
$\Delta$ with $0\preceq\Gamma\preceq\Delta$,

$\lVert\Gamma\rVert_{S^{q}}\leq\lVert\Delta\rVert_{S^{q}},\text{ and
}\lVert(\Gamma\Gamma^{\star})^{1/2}\rVert_{S^{q}}=\lVert\Gamma\Gamma^{\star}\rVert_{S^{q/2}}^{1/2}$

where the square root is well defined as $\Gamma\Gamma^{\star}$ is symmetric
positive semi-definite. Lastly, via direct calculation,

$\displaystyle(\Gamma-\Delta)(\Gamma-\Delta)^{\star}$ $\displaystyle\preceq
2(\Gamma\Gamma^{\star}-\Delta\Delta^{\star})$
$\displaystyle(\Gamma-\Delta)^{\star}(\Gamma-\Delta)$ $\displaystyle\preceq
2(\Gamma^{\star}\Gamma-\Delta^{\star}\Delta).$

###### Proof.

For $\mathbb{S}_{n}$ the symmetric group on $n$ elements, let
$f:\mathbb{S}_{n}\rightarrow\mathbb{R}$ by

$f(\pi):=\left\lVert\sum_{i=1}^{m}X_{\pi(i)}-\sum_{i=m+1}^{2m}X_{\pi(i)}\right\rVert$

For $k=1,\ldots,m$, we define $B_{k,\pi}=X_{\pi(k)}-X_{\pi(k+m)}$ and
$H_{k,\pi}=\sum_{i=k+1}^{m}B_{i,\pi}=\sum_{i=k+1}^{m}X_{\pi(i)}-\sum_{i=m+k+1}^{2m}X_{\pi(i)}$
where $H_{m,\pi}=0$ being an empty sum.

Note that the $B_{k,\pi}$ are symmetric random variables for $\pi$ uniform on
$\mathbb{S}_{n}$. Thus, $\mathrm{E}_{\pi}\lVert
f(\pi)\rVert=\mathrm{E}_{\pi}\lVert
B_{1,\pi}+H_{1,\pi}\rVert=\mathrm{E}_{\pi}\lVert-B_{1,\pi}+H_{1,\pi}\rVert$
and furthermore, letting $\delta_{1},\ldots,\delta_{m}$ be iid Rademacher
random variables,

$\displaystyle\mathrm{E}_{\pi}\lVert f(\pi)\rVert^{p}$
$\displaystyle=\mathrm{E}_{\pi}\mathrm{E}_{\delta_{1}}\lVert\delta_{1}B_{1,\pi}+H_{1,\pi}\rVert$
$\displaystyle=\mathrm{E}_{\pi}\mathrm{E}_{\delta_{1}}\mathrm{E}_{\delta_{2}}\lVert\delta_{1}B_{1,\pi}+\delta_{2}B_{2,\pi}+H_{2,\pi}\rVert$
$\displaystyle=\mathrm{E}_{\pi}\mathrm{E}_{\delta_{1}}\ldots\mathrm{E}_{\delta_{m}}\left\lVert\sum_{i=1}^{m}\delta_{i}B_{i,\pi}\right\rVert$

From here, we consider separately the commutative and non-commutative
settings.

For the commutative setting, we apply the facts about Schatten norms preceding
this proof. Beginning with the classic Kahane-Khintchine inequality from above
with $p^{\prime}=2$, we have

$\mathrm{E}_{\pi}\lVert f(\pi)\rVert_{q}^{p}\leq
C_{p}\left\lVert\left(\sum_{i=1}^{m}B_{i,\pi}{B_{i,\pi}^{\star}}\right)^{1/2}\right\rVert_{S^{q}}^{p}.$

Noting that $B_{i,\pi}{B_{i,\pi}^{\star}}\leq
2(X_{\pi(k)}{X}_{\pi(k)}^{\star}+X_{\pi(k+m)}{X}_{\pi(k+m)}^{\star}),$

$\left\lVert\left(\sum_{i=1}^{m}B_{i,\pi}{B_{i,\pi}^{\star}}\right)^{1/2}\right\rVert_{S^{q}}^{p}=\left\lVert\sum_{i=1}^{m}B_{i,\pi}{B_{i,\pi}^{\star}}\right\rVert_{S^{q/2}}^{p/2}\\\
\leq
2^{p/2}\left\lVert\sum_{i=1}^{2m}X_{i}{X_{i}^{\star}}\right\rVert_{S^{q/2}}^{p/2}=2^{p/2}\left\lVert\left(\sum_{i=1}^{2m}X_{i}{X_{i}^{\star}}\right)^{1/2}\right\rVert_{S^{q}}^{p}$

For the non-commutative setting, we proceed as before using the non-
commutative variant of Kahane-Khintchine and also noting that
$B_{i,\pi}^{\star}{B_{i,\pi}}\leq
2(X_{\pi(k)}^{\star}{X}_{\pi(k)}+X_{\pi(k+m)}^{\star}{X}_{\pi(k+m)}).$

$\displaystyle\mathrm{E}_{\pi}\lVert f(\pi)\rVert_{q}^{p}$ $\displaystyle\leq
C_{p}\max\left\\{\left\lVert\left(\sum_{i=1}^{m}B_{i,\pi}B_{i,\pi}^{\star}\right)^{1/2}\right\rVert_{S^{q}}^{p},\left\lVert\left(\sum_{i=1}^{m}B_{i,\pi}^{\star}B_{i,\pi}\right)^{1/2}\right\rVert_{S^{q}}^{p}\right\\}$
$\displaystyle=C_{p}\max\left\\{\left\lVert\sum_{i=1}^{m}B_{i,\pi}B_{i,\pi}^{\star}\right\rVert_{S^{q}}^{p/2},\left\lVert\sum_{i=1}^{m}B_{i,\pi}^{\star}B_{i,\pi}\right\rVert_{S^{q}}^{p/2}\right\\}$
$\displaystyle\leq
C_{p}2^{p/2}\max\left\\{\left\lVert\sum_{i=1}^{2m}X_{i,\pi}X_{i,\pi}^{\star}\right\rVert_{S^{q}}^{p/2},\left\lVert\sum_{i=1}^{2m}X_{i,\pi}^{\star}X_{i,\pi}\right\rVert_{S^{q}}^{p/2}\right\\}$
$\displaystyle=C_{p}2^{p/2}\max\left\\{\left\lVert\left(\sum_{i=1}^{2m}X_{i,\pi}X_{i,\pi}^{\star}\right)^{1/2}\right\rVert_{S^{q}}^{p},\left\lVert\left(\sum_{i=1}^{2m}X_{i,\pi}^{\star}X_{i,\pi}\right)^{1/2}\right\rVert_{S^{q}}^{p}\right\\}$

∎

####  A.2.1 On Optimal Constants

For the classic Khintchine inequality, the optimal constants due to Haagerup
(1981) coincide with the lower bound imposed by the central limit theorem.
That is, Khintchine's inequality states that $\mathrm{E}\left\lvert\sum
x_{i}\varepsilon_{i}\right\rvert^{p}\leq B_{p}\left\lVert
x\right\rVert_{2}^{p}$ where

$B_{p}=2^{p/2}\Gamma\left\\{(p+1)/2\right\\}/\sqrt{\pi}.$

This coincides precisely with the $p$th absolute moment of a standard normal
random variable—i.e. $\mathrm{E}\lvert Z\rvert^{p}=B_{p}$ for
$Z\sim\mathcal{N}\left(0,1\right)$.

For the Kahane-Khintchine inequality, optimal constants are not currently
known.777For the lower bound, optimal constants are known due to Latała and
Oleszkiewicz (1994). However, it is strongly conjectured that they coincide
with those in the standard Khintchine inequality. Moreover in the multivariate
setting, due again to the central limit theorem, the optimal constant has a
lower bound. Indeed, let $Z\sim\mathcal{N}\left(0,\Sigma\right)$, then

${\mathrm{E}\lVert
Z\rVert_{\ell^{q}}^{q}}\lVert\Sigma^{1/2}\rVert_{S^{q}}^{-1}=2^{p/2}\Gamma\left\\{(p+1)/2\right\\}/\sqrt{\pi}.$

This can be extended into a functional data setting using the fact that the
space of covariance operators arises from the closure of the set of finite
rank operators—i.e. the multivariate setting.

###  A.3 Sub-Gaussian Concentration

Given upper bounded on the $p$th moments of a random permutation statistic, we
want to quantify the concentration behaviour. In particular, we want as sharp
an upper bound as possible to achieve the best statistical power for
hypothesis testing.

We first consider the standard moment bounds to achieve sub-Gaussian
concentration (Boucheron et al., 2013) in Proposition A.8. This is improved if
$X$ is symmetric (Garling, 2007) in Proposition A.9. Lastly, even if the
moment condition is weakened as in Proposition A.11, we still have sub-
Gaussian concentration.

######  Proposition A.8.

For a centred univariate random variable $X\in\mathbb{R}$ such that
$\mathrm{E}{\lvert X\rvert}^{2p}\leq p!C^{p}$ for some constant $C>0$. Then,

$\mathrm{P}\left(X>t\right)\leq\mathrm{e}^{-t^{2}/8C}.$

######  Proposition A.9.

For a centred symmetric univariate random variable $X\in\mathbb{R}$ such that
$\mathrm{E}{\lvert X\rvert}^{2p}\leq p!C^{p}$ for some constant $C>0$. Then,

$\mathrm{P}\left(X>t\right)\leq\mathrm{e}^{-t^{2}/2C}.$

######  Remark A.10.

Note that the difference between the above two propositions is a factor of 4
in the denominator of the exponent. This stems from a standard symmetrization
trick where one considers $X$ and $X^{\prime}$, an iid copy of $X$, so that

$\mathrm{E}\lvert X-X^{\prime}\rvert^{2p}\leq 2^{2p}\mathrm{E}\lvert
X\rvert^{2p}\leq(4C)^{p}p!.$

Thus, the following results can be similarly adjusted for asymmetric random
variables.

######  Proposition A.11.

For a centred symmetric univariate random variable $X\in\mathbb{R}$ such that
$\mathrm{E}{\lvert X\rvert}^{2p}\leq(2p)!C^{p}/p!$ for some constant $C>0$.
Then,

$\mathrm{P}\left(X>t\right)\leq\mathrm{e}^{-t^{2}/4C}.$

###### Proof.

The moment generating function is

$\mathrm{E}\mathrm{e}^{\lambda
Z}=\sum_{p=0}^{\infty}\frac{\lambda^{p}\mathrm{E}Z^{p}}{p!}=\sum_{p=0}^{\infty}\frac{\lambda^{2p}\mathrm{E}Z^{2p}}{(2p)!}\leq\sum_{p=0}^{\infty}\frac{\lambda^{2p}C^{p}}{p!}\leq\mathrm{e}^{\lambda^{2}C}.$

The result follows from Markov's (Chernoff's) Inequality. ∎

##  Appendix B Proofs of main theorems

Now that all of the results from the previous section have been established,
we prove the tail bounds on the test statistics of interest by (1) applying
the appropriate Khintchine-type moment bound and (2) applying the appropriate
sub-Gaussian bound on the moment generating function.

###### Proof of Theorem 2.1.

For the balanced case of $\kappa=1$, let $n=2m$ and
$\varepsilon_{1},\ldots,\varepsilon_{n}$ be Rademacher random variables such
that $\sum_{i=1}^{n}\varepsilon_{i}=0$—i.e. not independent. Then, we can
rewrite $T(\pi)$ from equation 2.1 as

$T(\pi)=\frac{1}{s{m}}\sum_{i=1}^{n}\varepsilon_{i}X_{i}.$

Treating $X_{i}\in\mathbb{R}$ as fixed, we can use Theorem A.2 to bound the
$p$th absolute moment of $T(\pi)$ for $\pi$ uniformly distributed on
$\mathbb{S}_{n}$,

$\mathrm{E}_{\varepsilon}\left\lvert
T(\pi)\right\rvert^{p}=\left(\frac{1}{s{m}}\right)^{p}\mathrm{E}_{\varepsilon}\left\lvert\sum\varepsilon_{i}X_{i}\right\rvert^{p}\leq
B_{p}\left(\frac{\lVert
X\rVert_{2}^{2}-n\bar{X}^{2}}{s^{2}m^{2}}\right)^{p/2}.$

However, the term $\lVert X\rVert_{2}^{2}-n\bar{X}^{2}=(n-1){s}^{2}<2ms^{2}$.
Hence, the result of Spektor (2016) can be equivalently rewritten as

$\mathrm{E}\left\lvert
T(\pi)\right\rvert^{2p}\leq(2/m)^{p}B_{2p}=\frac{2^{p}(2p)!}{m^{p}p!}.$

Applying Proposition A.11 gives the desired result.

For $\kappa>1$—i.e. the imbalanced setting—we apply Theorem A.4 to get moment
bounds

$\mathrm{E}\left\lvert
T(\pi)\right\rvert^{2p}\leq\left(\frac{(\kappa+1)^{2}}{2m_{2}}\right)^{p}\frac{(2p)!}{p!}.$

and Proposition A.11 again to get the desired result. ∎

###### Proof of Theorem 2.2.

As with the previous proof, let $n=2m$ and
$\varepsilon_{1},\ldots,\varepsilon_{n}$ be Rademacher random variables such
that $\sum_{i=1}^{n}\varepsilon_{i}=0$. Our permuted test statistic is
$T(\pi)=\left\lVert\sum_{i=1}^{n}\varepsilon_{i}X_{i}\right\rVert_{q}$. We
apply Theorem A.7, our Kahane-Khintchine variant assuming the above dependency
on $\varepsilon$, in the commutative Banach setting to get

$\mathrm{E}T(\pi)^{p}\leq
C_{p}^{p}2^{p/2}\left\lVert\left(\sum_{i=1}^{n}X_{i}X_{i}^{*}\right)^{1/2}\right\rVert^{p}.$

Note that while the optimal constant is not known, $C_{p}\sim p^{1/2}$ from
the central limit theorem and from the proof in Diestel et al. (1995), Chapter
11. Hence, applying the fact that $((2p)!/p!)^{1/2p}\sim p^{1/2}$ and
Proposition A.11. We have the desired result. ∎

###### Proof of Theorem 2.3.

This proof is identical to that for Theorem 2.2 except we apply the non-
commutative variant of Kahane-Khintchine. ∎

###### Proof of Proposition 2.5.

We note first that $nT_{0}^{2}/(2+\kappa+\kappa^{-1})$ is approximately
$\chi^{2}\left(1\right)$ via the central limit theorem. Hence, for
$Z\sim\chi^{2}\left(1\right)$, some $c>0$, and some $u\in(0,1)$,

$\displaystyle\mathrm{P}\left(\mathrm{e}^{-Z/c}\leq u\right)$
$\displaystyle=\mathrm{P}\left(Z\geq-c\log u\right)$
$\displaystyle=(2\pi)^{-1/2}\int_{-c\log
u}^{\infty}x^{-1/2}\mathrm{e}^{-x/2}dx$
$\displaystyle=\left(\frac{c}{2\pi}\right)^{1/2}\int_{0}^{u}(-\log
y)^{-1/2}y^{c/2-1}dy$
$\displaystyle\leq\left(\frac{c}{2\pi}\right)^{1/2}\int_{0}^{u}(1-y)^{1/2-1}y^{c/2-1}dy$
$\displaystyle=\frac{(c/2)^{1/2}\Gamma(c/2)}{\Gamma((c+1)/2)}I(u;c/2,1/2)$

where we use the inequality $-\log y\geq 1-y$ for $y\in(0,1)$. The coefficient
${(c/2)^{1/2}\Gamma(c/2)}{\Gamma((c+1)/2)}^{-1}\rightarrow 1$ as
$c\rightarrow\infty$. Replacing $c$ with
$2\lceil\kappa+1\rceil^{3}/(2+\kappa+\kappa^{-1})$, we conclude that

$\mathrm{P}\left(\exp\left\\{-nT(\pi)^{2}/2\lceil\kappa+1\rceil^{3}\right\\}<u\right)\leq
C_{0}I\left(u;\frac{\lceil\kappa+1\rceil^{3}}{(2+\kappa+\kappa^{-1})},\frac{1}{2}\right)$

where
$C_{0}={\left(\frac{\lceil\kappa+1\rceil^{3}}{2+\kappa+\kappa^{-1}}\right)^{1/2}\Gamma\left(\frac{\lceil\kappa+1\rceil^{3}}{2+\kappa+\kappa^{-1}}\right)}{\Gamma\left(\frac{1}{2}+\frac{\lceil\kappa+1\rceil^{3}}{2+\kappa+\kappa^{-1}}\right)^{-1}}$.
∎

###### Proof of Theorem 2.4.

Let $Z=h(\lVert X\rVert)\in\mathbb{R}^{+}$, and let $B^{*}$ be a countable
dense subset of the unit ball of the dual space $\mathcal{X}^{*}$, which
consists of bounded linear functionals $\phi$. Then, we can write
$Z=\sup_{\phi\in B^{*}}h(\phi(X))$ being a countable supremum. Via application
of Talagrand's concentration inequality (Talagrand, 1996), we have that

$\mathrm{P}\left(Z\geq\mathrm{E}Z+t\right)\leq\exp\left(\frac{-t^{2}}{a+bt}\right)$

for positive constants $a$ and $b$ depending on $\mathrm{E}h(\lVert
X\rVert)^{2}$ and $\sup_{X\in\mathcal{X}}h(\lVert X\rVert)$.888 Refined values
for such constants can be found in other works (Bousquet, 2003; Klein and Rio,
2005; Giné and Nickl, 2016), but are not pertinent to this discussion.  Noting
that for $t\geq 0$

$\frac{d}{dt}\left\\{\frac{t}{1+bt/a}\right\\}=\frac{1}{(1+bt/a)^{2}}\leq\frac{1}{1+bt/a}=\frac{d}{dt}\left\\{\frac{a}{b}\log(1+bt/a)\right\\},$

we have that

$\exp\left(\frac{-t^{2}}{a+bt}\right)=\exp\left\\{-\frac{1}{b}\left(t-\frac{at}{a+bt}\right)\right\\}\\\
\leq\exp\left\\{-\frac{1}{b}\left(t-\frac{a}{b}\log(1+bt/a)\right)\right\\}=\mathrm{e}^{-t/b}\left(1+\frac{b}{a}t\right)^{a/b^{2}}.$

If $a\notin\mathbb{N}$, then we replace $a$ with $\lceil a\rceil$. Then, we
have that

$\exp\left(\frac{-t^{2}}{a+bt}\right)\leq\mathrm{e}^{-t/b}\left[\left(1+\frac{bt}{a}\right)^{a}\right]^{1/b^{2}}\\\
=\mathrm{e}^{-t/b}\left[\sum_{k=0}^{a}{a\choose
k}\left(\frac{bt}{a}\right)^{k}\right]^{1/b^{2}}\leq\mathrm{e}^{-t/b}\left[\sum_{k=0}^{a}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}.$

If $b=1$, then this is just the distribution function of the Erlang (gamma)
distribution with shape parameter $a$ and scale parameter $1$. More generally,
we have that

$\left[\sum_{k=0}^{a}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}=\mathrm{e}^{t/b}\left[1-\mathrm{e}^{-bt}\sum_{k=a+1}^{\infty}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}$

whose $l$th derivative for $b>1$, denoting the Pochhamer symbol
$(b^{-2})_{j}=\prod_{i=1}^{j}(b^{-2}-i+1)$, can be written as

$\displaystyle\frac{d^{l}}{dt^{l}}$
$\displaystyle\left[\sum_{k=0}^{a}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}$
$\displaystyle=\frac{\mathrm{e}^{t/b}}{b^{l}}\left[1-\mathrm{e}^{-bt}\sum_{k=a+1}^{\infty}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}$
$\displaystyle~{}~{}~{}-\mathrm{e}^{t/b}\sum_{j=1}^{l}b^{j-l}(b^{-2})_{j}\left[1-\mathrm{e}^{-bt}\sum_{k=a+1}^{\infty}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}-j}\frac{d^{j}}{dt^{j}}\left\\{\mathrm{e}^{-bt}\sum_{k=a+1}^{\infty}\frac{1}{k!}\left(bt\right)^{k}\right\\}$
$\displaystyle=\frac{\mathrm{e}^{t/b}}{b^{l}}\left[1-\mathrm{e}^{-bt}\sum_{k=a+1}^{\infty}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}$
$\displaystyle~{}~{}~{}-\mathrm{e}^{t/b}\sum_{j=1}^{l}b^{2j-l}(b^{-2})_{j}\left[1-\mathrm{e}^{-bt}\sum_{k=a+1}^{\infty}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}-j}\mathrm{e}^{-bt}\sum_{i=0}^{j}(-1)^{j-i}\sum_{k=0\vee(a+1-i)}^{\infty}\frac{1}{k!}\left(bt\right)^{k}.$

Thus, for $l\leq a$, we have that

$\left.\frac{d^{l}}{dt^{l}}\left[\sum_{k=0}^{a}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}\right|_{t=0}=\frac{1}{b^{l}}$

as the second term vanishes, and for $l\geq a+1$ and $b>1$, we have that

$\displaystyle\left.\frac{d^{l}}{dt^{l}}\left[\sum_{k=0}^{a}\frac{1}{k!}\left(bt\right)^{k}\right]^{1/b^{2}}\right|_{t=0}$
$\displaystyle=\frac{1}{b^{l}}\left[1-\sum_{j=a+1}^{l}b^{2j}(b^{-2})_{j}\sum_{i=a+1}^{j}(-1)^{j-i}\right]$
$\displaystyle=\frac{1}{b^{l}}\left[1-\sum_{\begin{subarray}{c}j=a+1\\\
j=a+1\text{ mod }2\end{subarray}}^{l}\prod_{i=1}^{j}(1-(i+1)b^{2})\right].$

which is negative for odd $a$. Thus, for $a$ odd—in the case where
$a\in\mathbb{R}^{+}$, we replace $a$ with $2\lfloor a/2\rfloor+1$—we can
finally bound via $a$th order approximation

$\exp\left(\frac{-t^{2}}{a+bt}\right)\leq\mathrm{e}^{-t/b}\left[\left(1+\frac{bt}{a}\right)^{a}\right]^{1/b^{2}}\leq\mathrm{e}^{-t/b}\sum_{k=0}^{a}\frac{1}{k!}\left(\frac{t}{b}\right)^{k}=\int_{t}^{\infty}\frac{x^{a-1}\mathrm{e}^{-x/b}}{b^{a}\Gamma(a)}dx$

being once again the Erlang (gamma) distribution function with shape parameter
$a$ and scale parameter $b$.

As a result, we have for $u\in(0,1)$, $C$ some positive constant, and
$I(u;a,c/b-a)$ the incomplete beta function where $c$ is chosen large enough
so that $c/b-a>0$,

$\displaystyle\mathrm{P}\left(\mathrm{e}^{-(Z-\mathrm{E}{Z})/c}\leq u\right)$
$\displaystyle=\mathrm{P}\left(Z-\mathrm{E}{Z}\geq-c\log u\right)$
$\displaystyle\leq\int_{-c\log
u}^{\infty}\frac{x^{a-1}\mathrm{e}^{-x/b}}{b^{a}\Gamma(a)}dx$
$\displaystyle=\frac{c^{a}}{b^{a}\Gamma(a)}\int_{0}^{u}(-\log
y)^{a-1}{y^{c/b-1}}dy$
$\displaystyle\leq\frac{c^{a}}{b^{a}\Gamma(a)}\int_{0}^{u}(1-y)^{a-1}{y^{c/b-a-1}}dy=CI(u;a,c/b-a)$

where the final inequality comes from $-\log
y=\sum_{k=1}^{\infty}(1-y)^{k}/k\leq(1-y)/y$ for $0<y<1$. ∎

###### Proof of Theorem 3.1.

Let $T=\lVert{X}^{\mathrm{T}}\varepsilon\rVert_{\ell^{2}}$ and note for $n=2m$
that $T^{2}=\sum_{i=1}^{n}a_{i,j}\varepsilon_{i}\varepsilon_{j}$ where
$a_{i,j}$ is the $ij$th entry of $X{X}^{\mathrm{T}}$. This is an homogeneous
Rademacher chaos of order 2.

As in Spektor (2016), we note the following correspondence. Let
$\Omega=\\{\varepsilon\in\\{\pm 1\\}^{n}\,|\,\sum\varepsilon_{i}=0\\},$ then

$\pi\in\mathbb{S}_{n}\longleftrightarrow\left\\{\varepsilon\in\Omega\,|\,\varepsilon_{i}=1\text{
if }\pi(i)\leq m\text{ and }\varepsilon_{i}=-1\text{ if }\pi(i)>m\right\\}.$

Hence, for any $\pi\in\mathbb{S}_{n}$, we can write

$T^{2}(\pi)=\sum_{i\leq m,j\leq m}a_{\pi(i),\pi(j)}-\sum_{i>m,j\leq
m}a_{\pi(i),\pi(j)}-\sum_{i\leq
m,j>m}a_{\pi(i),\pi(j)}+\sum_{i>m,j>m}a_{\pi(i),\pi(j)}$

and consider

$\mathrm{E}_{\varepsilon}\lvert T^{2}\rvert^{p}=\mathrm{E}_{\pi}\lvert
T^{2}(\pi)\rvert^{p}=\mathrm{E}_{\pi}\left\lvert\sum_{i=1}^{m}\left\\{a_{\pi(i),\pi(j)}-a_{\pi(i+m),\pi(j)}-a_{\pi(i),\pi(j+m)}+a_{\pi(i+m),\pi(j+m)}\right\\}\right\rvert^{p}.$

Writing
$b_{k,k^{\prime},\pi}=a_{\pi(k),\pi(k^{\prime})}-a_{\pi(k+m),\pi(k^{\prime})}-a_{\pi(k),\pi(k^{\prime}+m)}+a_{\pi(k+m),\pi(k^{\prime}+m)},$
and $H_{k,\pi}=\sum_{i,j\in I_{k}}b_{i,j,\pi}$ where the sum is over
$I_{k}=\\{1\leq i,j\leq m\,|\,i+j>k+1\\}$ with the empty sum being zero, we
note that

$T^{2}(\pi)=b_{1,1,\pi}+H_{1,\pi}=b_{1,1,\pi}+b_{1,2,\pi}+b_{2,1,\pi}+b_{2,2,\pi}+H_{2,\pi}=\ldots=\sum_{i,j=1}^{m}b_{i,j,\pi}.$

Then,

$\displaystyle\mathrm{E}_{\pi}\lvert T^{2}(\pi)\rvert^{p}$
$\displaystyle=\mathrm{E}_{\pi}\lvert
b_{1,1,\pi}+b_{1,2,\pi}+b_{2,1,\pi}+b_{2,2,\pi}+H_{2,\pi}\rvert^{p}$
$\displaystyle=\mathrm{E}_{\pi}\lvert
b_{1,1,\pi}-b_{1,2,\pi}-b_{2,1,\pi}+b_{2,2,\pi}+H_{2,\pi}\rvert^{p}$
$\displaystyle=\mathrm{E}_{\pi}\mathrm{E}_{\delta}\lvert\delta_{1}\delta_{1}b_{1,1,\pi}+\delta_{1}\delta_{2}b_{1,2,\pi}-\delta_{2}\delta_{1}b_{2,1,\pi}+\delta_{2}\delta_{2}b_{2,2,\pi}+H_{2,\pi}\rvert^{p}$

for $\delta_{1},\delta_{2}$ iid Rademacher random variables. Continuing in
this fashion, we have $\mathrm{E}_{\pi}\lvert
T^{2}(\pi)\rvert^{p}=\mathrm{E}_{\pi}\mathrm{E}_{\delta}\lvert\sum_{i,j=1}^{m}\delta_{i}\delta_{j}b_{i,j,\pi}\rvert^{p}.$
From here we apply Corollary 3 from Kwapien (1987).

First note that as the $\delta_{i}$ are just iid Rademacher random variables,
the standard Khintchine (or Kahane-Khintchine) inequality applies with
coefficient $B_{2p}=((2p)!/2^{p}p!)^{1/2p}$. Then Corollary 3 from Kwapien
(1987) to this degree 2 polynomial chaos implies that

$\left[\mathrm{E}_{\pi}\mathrm{E}_{\delta}\left\lvert\sum_{i,j=1}^{m}\delta_{i}\delta_{j}b_{i,j,\pi}\right\rvert^{p}\right]^{1/p}\leq
B_{p}^{2}C\left[\mathrm{E}_{\pi}\mathrm{E}_{\delta}\left\lvert\sum_{i,j=1}^{m}\delta_{i}\delta_{j}b_{i,j,\pi}\right\rvert^{2}\right]^{1/2}$

where $C$ is a universal constant which for homogeneous degree $d$ polynomial
chaoses is $d^{3d}/d!$ or simply $2^{5}=32$ in our case. The expectation on
the right hand side then becomes

$\left[\mathrm{E}_{\pi}\mathrm{E}_{\delta}\left\lvert\sum_{i,j=1}^{m}\delta_{i}\delta_{j}b_{i,j,\pi}\right\rvert^{2}\right]^{1/2}=\left[\mathrm{E}_{\pi}\sum_{i,j=1}^{m}b_{i,j,\pi}^{2}\right]^{1/2}\leq\left[\sum_{i,j=1}^{m}4a_{i,j}^{2}\right]^{1/2}=2\lVert
X{X}^{\mathrm{T}}\rVert_{S^{2}}.$

Absorbing the $2$ into $C$, we have the moment bounds

$\mathrm{E}_{\pi}\lvert T(\pi)\rvert^{2p}\leq B_{p}^{2p}C^{p}\lVert
X{X}^{\mathrm{T}}\rVert_{S^{2}}^{p}.$

To adapt these moment bounds into a tail bound, we use the standard moment
generating function approach, but in preparation we first recall the Legendre
duplication formula $\Gamma(2p)=2^{2p-1}\Gamma(p)\Gamma(p+1/2)/\sqrt{\pi}$ and
then note the following:

$\frac{2^{p}(\Gamma(p+1))^{2}}{\Gamma(2p+1)(\Gamma(p/2+1))^{2}}=\frac{2^{p}p^{2}(\Gamma(p))^{2}}{2p\Gamma(2p)(p/2)^{2}(\Gamma(p/2))^{2}}=\frac{2^{p+1}(\Gamma(p))^{2}}{p\Gamma(2p)(\Gamma(p/2))^{2}}=\\\
\frac{2^{-p+2}\sqrt{\pi}\Gamma(p)}{p\Gamma(p+1/2)(\Gamma(p/2))^{2}}\leq\frac{\sqrt{\pi}\Gamma(p)}{2^{p-2}p\Gamma(p+1/2)}\frac{p(2\mathrm{e})^{p}}{4\pi
p^{p}}=\frac{\Gamma(p)}{\Gamma(p+1/2)}\frac{\mathrm{e}^{p}}{\sqrt{\pi}p^{p}}\leq\\\
\frac{\sqrt{p+\pi^{-1}}\mathrm{e}^{p}}{\sqrt{\pi}p^{p+1}}\leq\left[\frac{1}{\sqrt{\pi}}+\frac{1}{\pi\sqrt{p}}\right]\frac{\mathrm{e}^{p}}{p^{p+1/2}}\leq\left[\frac{1}{\sqrt{\pi}}+\frac{1}{\pi\sqrt{p}}\right]\frac{\mathrm{e}}{p!}\leq\left[\frac{\mathrm{e}}{\sqrt{\pi}}+\frac{\mathrm{e}}{\pi}\right]\frac{1}{p!},$

because, via Watson's formula (Watson, 1959),

$\frac{\Gamma(p)}{\Gamma(p+1/2)}=\frac{\Gamma(p+1)}{p\Gamma(p+1/2)}\leq\frac{\sqrt{p+\pi^{-1}}}{p}.$

Let $\pi$ and $\pi^{\prime}$ be independent uniform random permutations from
$\mathbb{S}_{n}$. Then, updating $C$ as necessary,

$\displaystyle\mathrm{E}_{\pi}\exp(\lambda T(\pi))$
$\displaystyle\leq\mathrm{E}_{\pi}\exp(\lambda(T(\pi)-T(\pi^{\prime}))$
$\displaystyle\leq\sum_{p=1}^{\infty}\frac{\lambda^{p}}{p!}\mathrm{E}_{\pi}\lvert
T(\pi)-T(\pi^{\prime})\rvert^{p}$
$\displaystyle\leq\sum_{p=1}^{\infty}\frac{\lambda^{2p}2^{2p}C^{p}}{(2p)!}\frac{(p!)^{2}}{2^{p}((p/2)!)^{2}}\lVert
X{X}^{\mathrm{T}}\rVert_{S^{2}}^{p}$
$\displaystyle\leq\sum_{p=1}^{\infty}\frac{\lambda^{2p}C^{p}}{p!}\lVert
X{X}^{\mathrm{T}}\rVert_{S^{2}}^{p}$
$\displaystyle\leq\mathrm{e}^{\lambda^{2}C\lVert
X{X}^{\mathrm{T}}\rVert_{S^{2}}},$

which gives the desired sub-Gaussian concentration as in Proposition A.11. ∎

##  Appendix C Additional Data Experiments

###  C.1 Multivariate Data

We test the performance of the bound in Theorem 2.2 on simulated multivariate
Gaussian data in $\ell^{q}(\mathbb{R}^{12})$ for $q=1,2,\infty$. The sample
size is $m_{1}=m_{2}=50$. Figure 6 displays the result of running such a two-
sample test for each of the three norms compared to the standard permutation
test approximated by sampling 1000 permutations. This was replicated 1000
times and the average $\log_{2}$ p-values are plotted. We see that the Kahane
bound does not achieve as much power as the standard permutation test.
However, after applying the empirical beta adjustment from Section 2.4 with
moments computed via 10 permutations, the computed p-values align perfectly
with the standard permutation test.

Figure 6:  Multivariate two sample test for normal data with balanced sample
sizes $m_{1}=m_{2}=50$ for $\ell^{1}$, $\ell^{2}$, and $\ell^{\infty}$ norms.
The plots compare the permutation test (red) with 1000 permutations to the
Kahane bound from Theorem 2.2 (green) and the beta adjusted Kahane bound
(blue) across 1000 replications.

###  C.2 Berkeley Growth Curves Null Setting

In this section, we repeat the data analysis from Section 4.2. However, we
first remove the sex labels from the Berkeley growth curve dataset. Hence,
when sampling two sets of size 30, each resample will contain both male and
female curves. Thus, there should be on average no statistical difference
between the two sets. Over 100 replications for each of the three norms
$L^{1}$, $L^{2}$, and $L^{\infty}$ as well as the two bounds—unadjusted Kahane
and beta adjusted—we have Figure 7, which plots the empirical p-values against
the theoretical p-values from the uniform distribution on $[0,1]$. We see
large deviations for the unadjusted Kahane bound in the $L^{2}$ and
$L^{\infty}$ norms yielding an overly conservative hypothesis test. Table 4
displays the results of goodness-of-fit testing for the six sets of null
p-values with a similar conclusion.

Kahane Bound Beta Adjusted Bound $L^{1}$ $L^{2}$ $L^{\infty}$ $L^{1}$ $L^{2}$
$L^{\infty}$ KS test 3.6% $<$0.001% $<$0.001% 8.7% 87.9% 9.2% AD test 3.8%
$<$0.001% $<$0.001% 2.2% 79.7% 5.1%

Table 4:  This table contains p-values from both the Kolmogorov-Smirnov and
the Anderson-Darling goodness-of-fit tests for the 100 computed p-values
against the Uniform$[0,1]$ distribution.  Figure 7:  100 simulated null
p-values are plotted against the theoretical values from the uniform
distribution on $[0,1]$. There are massive deviations from uniformity for
unadjusted Kahane bound with the $L^{2}$ and $L^{\infty}$ norms.

###  C.3 Phoneme Curves Null Setting

Similar to Appendix C.2, we aim to test for whether or not the empirical beta
adjusted p-values for the phoneme curves from Section 4.3 follow a
Uniform[0,1] distribution in the null setting. To do this, we repeat the test
from Section 4.3 but remove the label information. Hence, the two samples of
size 40 will comprise operators from both phonemes, and there should be no
significant difference between the two samples.

Table 5 reports p-values from the Kolmogorov-Smirnov and Anderson-Darling
goodness-of-fit tests comparing the empirical distribution of the 100 two
sample test adjusted p-values to a uniform distribution on the unit interval.
Most of these empirical tests yield insignificant p-values especially after
taking multiple testing into account indicating no noticeable deviation from
uniformity. Hence, the empirical beta adjustment is able to account for the
overly conservative nature of the unadjusted Kahane bounds.

Kolmogorov-Smirnov Trace Norm Hilbert-Schmidt Norm Operator Norm A O d S A O d
S A O d S O 10.5 84.8 10.7 d 0.5 43.5 68.8 30.4 29.4 42.2 S 70.5 25.5 31.0
32.7 77.3 30.3 17.5 47.8 55.2 i 16.1 60.3 41.7 71.1 81.0 77.4 0.3 9.1 58.3 6.2
86.2 0.6 Anderson-Darling Trace Norm Hilbert-Schmidt Norm Operator Norm A O d
S A O d S A O d S O 16.1 60.5 1.7 d 0.3 8.0 12.0 2.4 3.1 2.2 S 52.2 19.0 8.7
7.8 2.4 14.0 4.3 17.3 18.0 i 5.6 14.8 11.8 10.4 22.2 15.6 0.8 0.05 65.1 6.5
65.7 1.0

Table 5:  A list of p-values from the Kolmogorov-Smirnov and the Anderson-
Darling tests for the two sample tests comparing two different phonemes with a
sample size of $m_{1}=m_{2}=10$ under the trace, Hilbert-Schmidt, and operator
norms.

###  C.4 Simulated Covariance Operator Data

In this section, we recreate the two-sample simulation study performed in
Pigoli et al. (2014) Section 3 to test our methodology. Let $\Sigma_{M}$ and
$\Sigma_{F}$ be the empirical covariance operators for the male and female
Berkeley growth curves, respectively. For $\gamma\in[0,6]$, we generate two
sets of $n=30$ curves from a Gaussian process with mean zero and with
covariance operator $\Sigma_{M}$ for the first group and
$\Sigma(\gamma)=[\Sigma_{M}^{1/2}+\gamma\Delta][\Sigma_{M}^{1/2}+\gamma\Delta]^{\star}$
where $\Delta=\Sigma_{F}^{1/2}R-\Sigma_{M}^{1/2}$ and $R$ is the operator that
minimizes the Procrustes distance between $\Sigma_{M}$ and $\Sigma_{F}$.
Specifically, $R=UV^{\star}$ where $U$ and $V$ come from the singular value
decomposition of $(\Sigma_{F}^{1/2})^{\star}(\Sigma_{M}^{1/2})=UDV^{\star}$.

For each $\gamma$, we test $H_{0}:\Sigma_{M}=\Sigma(\gamma)$ against
$H_{1}:\Sigma_{M}\neq\Sigma(\gamma)$ via a standard permutation test as in
Pigoli et al. (2014) with $512$ permutations and via our Kahane-Khintchine
bound. This is replicated 1000 times resulting in Figure 8. We see that for
the trace, Hilbert-Schimdt, and operator norms, the power loss for using our
upper bound is not much different from the standard permutation test. At
worst, the p-values are 2-4 times larger than necessary.

Figure 8:  Plotted p-values for a two sample test for equality of covariance
operators coming from Pigoli et al. (2014). From left to right, the trace,
Hilbert-Schmidt, and Operator norms are considered in the three plots.

##  Appendix D Vowel Data

###  D.1 Other Schatten Norms

In this section, we repeat the analysis performed in Section 5 by replacing
the trace norm with both the Hilbert-Schmidt and operator norms. In Figures 9
and 10, we display results analogous to those seen previously for the trace
norm. Most notably, as we consider larger values of $q$ for the $q$-Schatten
norms, the number of null hypotheses that we fail to reject increases
indicating less statistical power to distinguish vowel phonemes. This is in
line with much past work using Schatten norms on functional data (Pigoli et
al., 2014; Pigoli et al., 2018).

Figure 9:  $\log_{10}$ p-values for pairwise two sample tests between vowel
pairs under the Hilbert-Schmidt norm.  Figure 10:  $\log_{10}$ p-values for
pairwise two sample tests between vowel pairs under the operator norm.

###  D.2 Null Setting

To check that our methodology, specifically the empirical beta adjustment from
Section 2.4, achieves the correct empirical size and thus is neither
conservative nor anti-conservative, we first randomize all of the labels
within each of the Latin squares from Section 5. Then, we repeat the same
analysis as before. The 66 p-values produced for each of the 1, 2, and
$\infty$ Schatten norms is displayed in Figure 11. These QQ plots compare our
empirical p-values to the theoretical quantiles of the Uniform[0,1]
distribution. For each of the three norms, we do not see much deviation from
uniformity. Furthermore, for testing goodness-of-fit with the uniform
distribution, the Kolmogorov-Smirnov test returns p-values of 0.434, 0.782,
and 0.290 and the Anderson-Darling test p-values 0.161, 0.511, and 0.241 for
the trace, Hilbert-Schmidt, and operator norms, respectively. None of these
tests are significant indicating no noticeable deviation from uniformity.

Figure 11:  QQ Plots comparing the 66 null p-values from the vowel data after
beta adjustment to the quantiles of the uniform distribution.