# Bridging Phylogeny and Taxonomy with Protein-protein Interaction Networks

Long-Huei Chen*, Mohana Prasad Sathya Moorthy*, and Pratyaksh Sharma * All
authors contributed equally.

###### Abstract

The protein-protein interaction (PPI) network provides an overview of the
complex biological reactions vital to an organism’s metabolism and survival.
Even though in the past PPI network were compared across organisms in detail,
there has not been large-scale research on how individual PPI networks reflect
on the species relationships. In this study we aim to increase our
understanding of the tree of life and taxonomy by gleaming information from
the PPI networks. We successful created (1) a predictor of network statistics
based on known traits of existing species in the phylogeny, and (2) a
taxonomic classifier of organism using the known protein network statistics,
whether experimentally determined or predicted de novo. With the knowledge of
protein interactions at its core, our two models effectively connects two
field with widely diverging methodologies–the phylogeny and taxonomy of
species.

##  1 Introduction

###  1.1 Motivation

The wide availability of protein-protein interaction (PPI) networks across
vast number of organisms provide an unique opportunity to increase
understanding of the relationships between species. The taxonomic and
phylogenetic trees describe the similarity in physical characteristics and
evolutionary relationships between the species, respectively. Since the advent
of molecular biology, new techniques emerged that can facilitate the
traditionally classification of species based on morphological and fossil
evidence, such as the molecular clock hypothesis [1], genomic
phylostratigraphy [2] and molecular phylogeny [3].

Past studies indicated that the protein-protein interaction networks of
species is a reflection of the evolution of both individual proteins and the
organisms as a whole. This suggest that the information gleamed from the
networks can be used to classify or predict species in relations with their
evolutionary trajectories and taxonomic classification. This leads to our
observation that the protein-protein interaction network is effectively a
bridge between the phylogeny and taxonomy of species, and the conclusions made
on protein networks can be used to gain insight or make prediction in both
directions between phylogeny and taxonomy. In this study, we would use overall
network statistics (like clustering coefficient) of PPI networks for different
species.

###  1.2 Problem Definition

Briefly our work proceeds as follows: first, we explore how the various
network statistics compare and differ throughout the tree of life and get deep
insights into how protein interactions have evolved over time. Following this,
given the position in the tree of life, we build a predictor for protein
network statistics. Finally we create a taxonomic classifier which takes as
input the protein network statistics (either experimentally determined or
obtained de novo with the predictor) and returns the predicted taxonomy of the
species.

As the construction of PPI networks remains a painstaking task, typically
biologists determine phylogeny and the taxonomic classification of a new
species before the PPI network is created. Constructing PPI networks also
requires huge amounts of financial and research resources. The potential of
our model lies in the calculation of various protein network statistics before
any protein interactions are determined, as long as the phylogeny of an
organism is understood. With the help of our predictor, biologists can predict
statistics about the PPI networks of new species by just knowing its
phylogenetic position. These predicted PPI network statistics can aid in
constructing the PPI networks or finding loopholes/mistakes in the constructed
networks. The predicted values can also be used to further predict interesting
properties of the species such as its taxonomy (a method which we developed in
this work).

We then study the relationship between PPI network characteristics and
taxonomy of the species by building a species classifier that determines the
taxonomy given the PPI network statistics of the species. The classifier can
be used for systematic and automated taxonomy assignment based solely on
protein network statistics and to check the taxonomy determined through other
methods.

The combination of the predictor and the classifier can prove to be another
interesting and high impact tool. Given the phylogenetic position of a species
we can use our predictor to predict its network statistics; then using these
predicted network statistics, our classifier can determine the taxonomy of the
species thereby bridging phylogeny to taxonomy through PPI networks. Even
though we use the PPI network data to build this predictor and classifier, we
can find the taxonomy of a species given its phylogeny without knowing any
information about its PPI networks. Such a tool can explore the
interdependence between the logically distinct fields of Taxonomy and
phylogeny and provide a bridge between the fields.

###  1.3 Protein-protein Interaction (PPI) Networks and its Reflection on the
Evolution of Life

Protein-protein interactions (PPIs) are vital to the operation of all cellular
functions. Proteins are chains of amino acids linked by peptide bonds–the
combinations of which yield a vast array of molecules with countless shapes
and sizes. Proteins act alone or together to perform chemical transformations;
transmit signals; and provide structure/organization, transport, and
recognition among many other cellular roles. However, previous studies have
been done on multiple aspects of the protein network evolution and its
relationship with Information obtained from PPI network databases enables
creation of interaction networks, and allows for systematic investigation of
the complex biological activities within the cell. With the increase in
research on PPIs and advancement in computational tools, interaction networks
have been constructed and analyzed for thousands of species.

###  1.4 Making Sense of the Diversity of Life:
Phylogeny and Taxonomy

In this study we consider two vastly different classification schemes
biologists employ in their study of species relationships: phylogeny and
taxonomy.

A phylogenetic tree describes evolutionary relationships of a group of
species, and a phylogenetic tree depicts the evolutionary relationships among
organisms, while each branch points indicate when new species emerge from a
common ancestor. Therefore, each node in the phylogenetic tree is an actual
species and all non-leaf nodes are ancestors of one or more organisms.
Phylogenetic trees are usually based on morphological or genetic homology and
are build by comparing anatomical traits, genetic difference and by using
molecular clocks.

Taxonomy, on the other hand, is the study of organisms with the goals of
classifying living and extinct organisms according to a set of rules [4]. The
rules specify a hierarchy of groups, which form the non-leaf nodes of the
tree, whereas the organisms are assigned to the tree leaves based on
similarities and dissimilarities of their characteristics. Taxonomy is usually
richly informed by phylogenetics, but remains a methodologically and logically
distinct discipline.

Our understanding of species taxonomy is constantly updated based on the
increasing understanding of evolutionary relationships between species.
Advancements in molecular biology also aid in our study of species
relationships, including the use of molecular clocks (rates of change in RNA
or DNA sequences) and identification of orthologous proteins across the
phylogeny. As a result in this study we are going to consider both the
taxonomy and the phylogeny tree in our quest of a automatic classification of
life, and compare the viability of our model under various goals and settings.
The predictor and the classifier we build can act as a link between them
through PPI networks.

##  2 Related Work

Major research directions related to our theme involves study of protein
networks, phylogenetic tree and taxonomy. But,to our knowledge, this is the
first study attempting to bridge the knowledge of taxonomy and phylogeny using
protein networks.

###  2.1 Emergence of Protein Nodes in Networks

Past research on protein network evolution has focused on single species
network by observing how the protein nodes arise in the evolutionary
timescale. This is done by assigning an evolutionary age to each protein node
using techniques like phylogenetic stratigraphy [2]. Research has shown that
in the course of evolution, the protein nodes are added to the network in
accordance with the preferential-attachment model [5, 6]. Further research
provides details on protein nodes emergence mechanisms like gene duplication
and gene loss [7, 8] or the divergence model [9]. Even though these studies
tend to focus only on a single species and so are of interest primarily in the
biological context, they can in fact supplement and support our observations
during our comparison across the phylogenetic tree using hypotheses formed
from protein emergence models.

###  2.2 Protein Network Structure Evolution

Phylogenetic trees are usually based on morphological or genetic homology.
Research on protein network evolution has also focused on the high level
structure of protein networks. Some studies identify motifs and community
structures that are commonly preserved in the evolution process[10, 11], and
others have looked at the high level structural graph statistics in terms of
the entire graph [12]. The studies suggest that proteins evolved earlier are
more likely to stay on as local clustering structures [13], and that the
overall statistics of the graph is with a power law dynamics with the
evolution [14]. This suggest that our attempt to compare graph statistics is
based on the valid assumption that graph characteristics is related the way
that the interactions links arise, and since the rudimentary species have a
less-evolved network we can expect the network statistics can also reflect the
fact.

###  2.3 Protein Network Alignment Informs the Phylogenetic Tree

Protein network alignment methods [15, 16, 17, 18] can inform our approach and
suggests way to effective summarize the comparison of networks across species.
These methods attempt to assign node-node relationships for homologous
proteins in different species [19, 20], and report a similarity score [21]
that can suggest evolutionary relatedness between species.

###  2.4 Evolution of Protein Network Topology

Our work comparing protein network across species is unique in its focus on
the phylogenetic relationships of species and its relatedness with the protein
networks themselves. The work in previous studies [22, 23] is focused
primarily on the comparison of protein networks, and how it relates to the
evolutionary relationships we’ve come to understand between the species. To
this end the methods have been applied to generate an alternative phylogenetic
tree based solely on the topological characteristics of the protein networks
[24] using tree construction methods such as the average distance algorithm
applied to similarity score between all the PPI networks, without regards to
individual gene/protein evolution. This complements our observation that the
summary statistics of networks can provide information on the species and its
evolution.

(a) Depiction of bacteria, archaea, and eukaryota domains in the phylogenetic
tree

(b) Network diameter

(c) Average degree

(d) Standard deviation in degree distribution

(e) Entropy of degree distribution

(f) Spectral norm

Figure 1: Heatmaps for spatially visualizing the variation of various network
statistics over species in the phylogenetic tree of life. Only selected
figures for network statistics which show interesting patterns are included
for this report.

##  3 Approach

###  3.1 Data collection

#### Protein-protein Interaction Networks

STRING v10 is the data source of choice for studies concerning PPI networks
across multiple species, and version 10 of the database include PPI networks
for more than 2000 species. There are multiple ways to define edges between
the protein nodes, including functional relationships found with experimental
means, known biological pathways, existing studies suggesting functional
codependence, and de novo computational prediction.

We restrict our investigation to the more significant species within the
dataset determined by the number of published work on the species when
searching through PubMed. By limiting our work to species with more than 100
published work, we choose to focus on the 427 most studied organisms since
they are of higher interest to biologists and have more complete protein
networks. The problem remains, however, that the bias in network completeness
(missing nodes in some of the species networks) might introduce artifacts in
the result of our study. We therefore plan to take steps to ensure that the
biases are not interpreted as part of our experimental conclusion.

Further, the data set also contained edges for protein protein interactions
that were anticipated to be existing and not actually found in biological
studies. These edges were removed as we wanted to extract features from the
true network without any assumptions or predictions.

#### Phylogenetic Tree of Life

In recent times, there has been ample work on the problem of constructing
phylogenetic trees from molecular sequence data of species [25, 26, 27]. The
STRING database itself provides a phylogenetic tree of species as accessory
data to the interaction networks. Since the usage of identifiers is consistent
across the interaction networks and the tree of species, it is a natural
choice to work with the STRING tree of species. Figure 1a shows the STRING
tree of species annotated with the three domains (Archaea, Bacteria, and
Eukaryota).

The STRING tree of species, however, does not provide any information on the
similarity of species that are connected by an edge in the tree. To enrich the
STRING tree of life with such information we use the tree of life developed in
[28]. Since the species identifiers used in [28] are not consistent with those
used in STRING, we use heuristics (based on lineage of species) to generate a
mapping between the two sets of species identifiers.

#### Taxonomic Hierarchy

The NCBI has made available the Taxonomy Browser tool [29] which allows
convenient browsing of genome information as well as taxonomic lineage of
species. The species identifiers used in Taxonomy Browser are consistent with
those used in STRING. Thus, to collect information about taxonomic hierarchy,
we

  1. 1.

query the Taxonomy Browser with all the species in the STRING database,

  2. 2.

scrape the lineage of species from the results HTML generated, and

  3. 3.

build a directed tree by adding the lineages as paths in the tree.

The final tree is further processed (as described in section 6), before being
used for construction of our taxonomic lineage classifier.

(a) Average degree

(b) Density

(c) Diameter

(d) Global clustering coefficient

(e) 2-star density

(f) Edge entropy

Figure 2: Variation of select network statistics with distance of species from
the root in the phylogenetic tree of life.

###  3.2 Network Statistics

In the following section, we give a brief summary for each of the PPI network
metrics that form the basis of our study, and are included as features in our
classifier and predictor. With some exploratory work we found that the
following set of features gave the best results.

  1. 1.

Number of nodes: This captures the number of protiens in the species and in
turn the complexity of the species.

  2. 2.

Number of edges: This captures the number of types of protein-protein
interactions, again indicating the complexity of the species.

  3. 3.

Average degree: On an average, how many other proteins does a protein in this
species interact with.

  4. 4.

Maximum degree: How many proteins does the most active protein in this species
interacts with. This refers to the protein with the least specificity.

  5. 5.

Density: This is the ratio of number of edges to the number of possible edges
in the network. This indicates the average specificity of proteins.

  6. 6.

Number of connected components: This captures groups of proteins that donot
interact with one another.

  7. 7.

Fraction of nodes in largest connected component: These capture the complexity
of biochemical reaction in species.

  8. 8.

Full diameter: Diameter is the longest of all shortest paths in the network.

  9. 9.

Global clustering measure: It is ratio of number of closed triads to sum of
closed and open triads.

  10. 10.

Clustering coefficient: The clustering coefficient of a node is defined as the
fraction of existing edges among all the possible edges between the neighbor
nodes.

  11. 11.

2-star density: Ratio of number of 2-stars present to number of 2-stars
possible in the network.

  12. 12.

3-star density: Ratio of number of 3-stars present to the number of 3-stars
possible in the network.

  13. 13.

Entropy and Gini coefficient of degree distribution: Multiple studies focusing
on single species PPI network have shown that proteins emerge in accordance
with the preferential attachment model [5, 6, 30]. Even though in our study of
protein networks we compare degrees between network, not protein nodes, we can
expect the same properties of network growth in preferential attachment. We’re
interested in seeing how the average degree grows with respect to the
evolutionary years and advancement in evolution. To compare degree
distributions across networks, we collect meta statistics over the
distribution like mean, gini coefficient, entropy, etc.

  14. 14.

Assortative mixing : This captures preference for a node to be attached with a
node thats similar to it in terms of degree.

  15. 15.

k-cores features: A k-core of a graph is a maximal connected subgraph of the
graph in which all vertices have degree at least k. We compute the maximum k,
such that k-core exists and use that k, number of nodes in that subgraph,
number of edges in that subgraph as features.

We believe the above features would rightly reflect the complexity of protein-
protein interactions and other biochemical pathways of different species.

##  4 Analyzing the Dataset

We devoted a large amount of effort on the exploratory analyses of the
dataset, and we were able to report some of the important observations we
gathered. We focus on the following analysis in this section:

  1. 1.

Spatial visualization of variation of network features across the phylogenetic
tree: to help us identify macro-level relationships and form hypothesis about
observed patterns.

  2. 2.

Variation of network features at different stages of evolution: we aggregate
the network features for all species at the same number of branch hops from
the root in the phylogenetic tree, and report the variation of these
aggregates with the number of branch hops.

###  4.1 Visualization of network features on the Phylogenetic tree

To show how the various network metrics varies throughout the entire
phylogenetic tree, we generate heat maps of the phylogenetic tree with each
point showing the magnitude of the various network characteristics (see Fig.
1).

We observe that there are obvious differences between species in the Eukaryota
domain and the other prokaryotes (Archaea and Bacteria). This suggest that we
should achieve at least moderate success in our effort to classify organism
based on PPI network statistics.

While we omit a statistical evaluation of the hypothesis that network
properties are different across different domains of species, we outline such
a procedure. A simple null model for comparison would be the phylogenetic tree
with its node labels (and corresponding values of network statistics) shuffled
randomly. The differences in size and distribution of clusters (computed by
binning the network statistic values into small number of bins and using any
of the spectral clustering algorithms) shall give us a meaningful
quantification of our hypothesis.

###  4.2 Variation of network features with degree of evolution

The tree of life with edge weights obtained from [28], gives us a measure of
dissimilarity of all pairs of species that are connected by an edge. Since the
phylogenetic tree encodes evolutionary information such as ancestry and common
descent, we use the tree distance (sum of weights of edges on a path) of a
species from the root of the phylogenetic tree as a proxy for its degree of
evolution. It is observed that Bacteria, in general, are further away from the
root than Eukaryotes and Archaea–which can be reasoned since Bacteria were
among the very first living organisms to appear on Earth.

In search of interesting patterns, we studied the variation of network
statistics of species with the distance from the phylogenetic root of that
species as described above. A selection of interesting results is presented in
Figure 2. It can be seen from the figures that Archaea are species with lowest
distance from the phylogenetic root, and Bacteria are located farther away
from the root along with Eukaryota.

While we did not observe a simple dependence of any network property on the
distance from root: we remark that the average degree, density, 2-star density
tend to decrease as species move farther away from the root, whereas there is
a weak opposite trend for other statistics such as diameter, global clustering
coefficient, and edge entropy. Further, in all the included figures, it can be
seen that there is a rather apparent distinction between different domains
based on various network statistics. For Bacteria and Archaea in particular,
we can see the variation of network statistics within the domain is relatively
small. These observations lead us to conclude with confidence that it is
possible to classify species into taxonomic groups (at least at the level of
domains) using only the network statistics of their PPI networks.

##  5 Predicting Network Statistics from Position of Species in Phylogenetic
Tree

As we identified statistics that are meaningful representation of the protein
networks, our next step is to build a predictor of the same statistics based
on features gleamed based on the phylogeny of species. We construct and test
our model by splitting the 427 species we included into the train set and the
test set with a random 80/20 split. Network statistics of the nodes in the
train set is presumed to be known, whereas statistics of the test set are
those we set out to predict. The prediction is based on the known position of
species in both the train and test sets on the phylogenetic tree of life.
Position refers to the exact location of the species in the tree - the parent
and the children can precisely define it. Lets look at the features that can
capture the position of a node in the tree.

###  5.1 Feature Extraction

####  5.1.1 Sibling Feature

To calculate the sibling features, for every species we find the latest
ancestor that contains more than one descendants in the tree; these
descendants are considered the ”siblings” of the species node in question. We
then take the mean of all the network statistics of these siblings to generate
a set of sibling features that can be used to predict the network statistics
of the particular node.

$feature_{{\mathop{\rm sib}\nolimits}}^{{\mathop{\rm
stats}\nolimits}}=\frac{{\sum\limits_{{\mathop{\rm
siblings}\nolimits}}{stats}}}{{\left|{{\mathop{\rm
siblings}\nolimits}}\right|}}\textrm{ for every }stats$

####  5.1.2 Cousin Feature

Cousins of a node are defined as the nodes with the same level from root in
the tree as the node under consideration. To get the cousin features, we look
at the entire phylogeny to find in the 427 species those that are of the same
hops from root as the species in question. By doing this we can have a
estimation of how the hops from root affects network statistics of a
particular species we set out to predict.

$feature_{{\mathop{\rm cuz}\nolimits}}^{{\mathop{\rm
stats}\nolimits}}=\frac{{\sum\limits_{{\mathop{\rm
cousins}\nolimits}}{stats}}}{{\left|{{\mathop{\rm
cousins}\nolimits}}\right|}}\textrm{ for every }stats$

####  5.1.3 Caveats in Feature Extraction

During feature extraction, our first intuition is to create features based
solely on the species node position in the phylogeny instead of taking the
network statistics into account. This is motivated by the fact that a lot of
the internal nodes in the phylogeny do not include protein networks data, so
it’s impossible to obtain network statistics of these nodes without
estimation. In our attempt to create features that summarize the position of
species nodes in the tree of life, we created a bit vector for every internal
node in the phylogeny indicating which of the 427 species are descendants of
the internal node. We found as a result that the feature vectors are too large
and too sparse (since we include all the internal nodes as possible ancestors)
and perform poorly even with dimensionality reduction using principal
component analysis (PCA).

We therefore conclude that in order to effectively predict the network
statistics of species in the test set we need to tap into the network metrics
instead or relying merely on ancestry in the phylogeny. This is meaningful in
the suggestion that the phylogeny itself is insufficient in accurately
determining the traits of protein networks, and shows that the network nodes
are highly influenced by both the evolutionary relationships and the traits of
existing network statistics of nearby species in the phylogenetic tree.

###  5.2 Predictor Model Evaluation

We first divide the species into a train and test set with a random 80/20
split. From every network statistics, we train a linear support vector machine
(SVM, linear kernel with C = 100) model based on the sibling and cousin
features of the species in the train set. We then predict each of the various
network statistics for species in the test set using the , and then compare
the predicted values with the actual network statistics using the relative
error:

${\textrm{relative error}}=\frac{{|predicted-actual|}}{{actual}}$

The relative errors are then summarized across all species in the test set,
and for each type of network statistics we obtain a mean and standard
deviation of the particular statistics (Table LABEL:table:predictorresults).

TABLE I: The mean and standard deviation of the relative error of the predicted neted values are obtained with SVM models trained using features from the sibling and cousin nodes of species in the train set, and tested using the same features obtained with species in the test set. | Mean of Relative Error | Standard Deviation of Relative Error
---|---|---
Nodes | 0.244 | 0.233
Edges | 0.739 | 0.978
Average Degree | 0.371 | 0.560
Maximum Degree | 0.381 | 0.349
Density | 13.223 | 14.857
Number of Components | 1.771 | 7.575
Maximum Component Size | 0.123 | 0.340
Full Diameter | 0.128 | 0.111
Effective Diameter | 0.106 | 0.181
Global Clustering | 0.405 | 0.395
Clustering Coefficient | 0.136 | 0.103
Star Density 2 | 1.366 | 2.494
Star Density 3 | 14.120 | 41.749
Gini Coefficient of Degree Dist. | 0.089 | 0.061
Edge Entropy | 0.025 | 0.014
Assortative Mixing | 0.186 | 0.235
K-cores Maximum Degree | 0.337 | 0.364
K-cores Maximum Nodes | 0.664 | 0.796
K-cores Maximum Edges | 2.293 | 4.234

###  5.3 Analysis

With results from Table LABEL:table:predictorresults we can conclude that the
predictor works well with most network statistics considered in this study.
The mean error is in the range of 10-20% for most of the network statistics,
and the standard deviation is also considerably low. Interestingly, it seems
that a few of the network statistics are extremely hard to predict, like the
density and star density of the networks.

The reason that this model performs so poorly in a few cases suggests that
some of network statistics are quasi-random and possibly unrelated to the
evolutionary relationships. Our observation was that the star density and the
density of a protein-protein interaction network does not have readily
available biological interpretation; in other words, they are irrelevant in
the context of protein interactions and biochemical reactions in an organism.

The results we obtained can also inform our understanding of the interaction
between evolution and network growth. As we can see in the Table 2, the
diameter of protein network, the maximum component size and a few variance
measures of the degree distribution performs remarkably well in our model;
this combined with the biological interpretation detailed in Section 3.2
suggests that the graph statistics are readily interpretable in the context of
evolutionary biology.

Noting that statistics like number of edges are huge numbers and are generally
orders of magnitude different for different species, we believe predicting the
number of edges with 70% accuracy is still good result considering the
underlying nature of these statistics. Also, we can derive the number of edges
with a better accuracy by using predicted average degree and number of nodes
as they have far lower error rates.

##  6 Taxonomic Lineage Classifier

Figure 3: A partial taxonomic hierarchy shown as a rooted directed tree

One of the glaring observations from the analysis in Sections 4.1-4.2 was that
PPI network statistics differ significantly across species in the Archaea,
Bacteria, and Eukaryota domains. It is apparent that given the PPI network
statistics of a species, it is possible to classify it with reasonable
accuracy into one of the three domains. While in the interest of brevity, we
could not include a similar analysis which compares network statistics across
finer levels of taxonomic hierarchy, we observe that many network statistics
tend to be similar for species in the same taxonomic group and different for
species across different groups. Naturally, we were tempted to ask “Can a
species’ PPI network statistics be used to predict its taxonomic
classification?”

It should be noted that the taxonomic classification of a species and the
construction of its interaction species are two tasks that are fairly
independent of each other. Our work aims to form a bridge between the the
science of taxonomy and molecular interaction biology. Further, if the
taxonomic lineage classifier is coupled with the network statistics predictor
(discussed in Section 5), it will be possible for biologists to discover the
taxonomic lineage of a newly discovered species using only it’s hypothesized
position in the phylogenetic tree.

Figure 4: Performing recursive feature elimination on the classifier for
domains shows us that highest accuracy is achieved when using all of the 19
features.

###  6.1 Approach

The taxonomic hierarchy is a rooted directed tree, where internal nodes
represent taxonomic groups and the leaves represent individual species (see
Figure 3). Each path from the root to a leaf (species) represents the
taxonomic lineage of that species. Thus, the classification task requires us
to classify a species (using its PPI network statistics) into one of such
paths.

We approach the lineage classification task as a hierarchy of classification
tasks, each of which classifies a species into finer taxonomic groups given
that it has been correctly classified to some coarser taxonomic group. In
other words, we construct one classifier for each internal node in the
taxonomic hierarchy (tree), whose role is to classify an example into a
taxonomic group represented by one of the children of this node.

To construct the full taxonomic lineage of an example (network statistics of a
species), we use the classifiers at successive levels as follows:

$c_{i}=\text{classifier}_{c_{i-1}}(x)$ $c_{0}=\text{cellular organisms
(root)}$

where $x$ is the feature vector (of network statistics) to be classified,
$\text{classifier}_{c_{i-1}}$ is the classifier corresponding to node
$c_{i-1}$ in the taxonomic hierarchy, and $c_{i}$ represents a node in the
taxonomic hierarchy at distance $i$ from the root. The lineage of the species
with feature vector $x$ is then $(c_{0},c_{1},\dots,c_{n-1})$ where $n$ is the
height of the taxonomic hierarchy.

Each of the classifiers is a SVM with linear kernel $(C=100)$.
$\text{classifier}_{c_{i}}$ is trained on species for which $i$-th level in
lineage equals $c_{i}$.

The set of features used for classification is the same as those computed by
the predictor in Section 5 (see Table 2). Since the linear kernel SVM allows
us to extract weights corresponding to each feature, we performed recursive
feature elimination to minimize the number of features used. In Figure 4, it
can be seen that it is best to use all of the 19 network statistics as
features.

###  6.2 Caveats

The taxonomic hierarchy described in the previous sections is not exactly
well-structured as described. In particular,

  1. 1.

The length of lineage is different for different species, i.e. the height of
the taxonomic hierarchy is not uniform.

  2. 2.

The levels in the lineage do not necessarily map to the taxonomic levels of
domain, kingdom, phylum, class, order, family, genus, species.

  3. 3.

For many intermediate levels of lineage that do not correspond to any of the
above taxonomic levels, there is just one child node in the taxonomic
hierarchy.

For the purposes of this project, we restrict our dataset to only those
species for which we can extract information for the taxonomic levels of
domain, kingdom, phylum, class, order, family information from the lineage. We
also collapse all the intermediate levels in the lineage which do not
correspond to any of the above taxonomic levels. Since there are not enough
training examples within each genus, we restrict our classification to the
above six taxonomic levels.

Figure 5: Each blue disc marks the cross-validation accuracy of a particular
classifier whose taxonomic level is read from the X-axis. The area of the disc
is proportional to the number of training examples that were used for that
classifier. The weighted average of accuracies is plotted in red. Figure 6: In
contrast to Figure 5, here we show the accuracy for the full lineage
classification. Almost 100% accuracy is achieved at the domain level, whereas
lower than 10% accuracy is achieved for the lineage classification till the
taxonomic level of family.

###  6.3 Results

We perform a five-fold cross validation on the dataset reduced to 1537 species
after filtering. Note that we also include the species with lower publication
counts for this particular problem.

The cross-validation accuracies for each of the individual classifiers (one
for each node in the taxonomic hierarchy) are shown in Figure 5. In the
figure, each blue disc represents the accuracy of a single classifier. The
plot is arranged in order of taxonomic levels and a tree of classifiers can be
observed from the figure. We note that the classifier for discovering the
domain of a species performs exceptionally well. While the classifiers at
lower levels of taxonomy perform relatively poor, but considering that each of
them have much lower amount of training data available and conceivably larger
number of labels to classify into, their performance seems to exceed our
expectations.

The cross-validation accuracy for the full lineage classification is shown in
Figure 6. The plot shows the classification accuracy up till various taxonomic
levels. It can be seen that the classifier (rather, a hierarchy thereof)
performance gets poorer as we seek to discover finer taxonomic classifications
for a species. It is worth mentioning that for a test set of 365 species, our
hierarchy of classifiers was able to predict the full lineage correctly for 15
species, while predicting the domain correctly for 361 species.

###  6.4 Analysis

The results, particularly from Figure 6, lead us to conclude that our
hierarchy of classifiers can be used to predict the lineage (at least
partially) of a species using only the network statistics of its PPI networks.
This is particularly impressive considering the fact that the taxonomy is
constructed not merely with the evolutionary relationships of species, but
rather based on a set of rules on the physical, physiological and other
aspects of species traits.

With our varying degree of success applying the classification model to the
different granularity of the taxonomic hierarchy, we can identify
classification rules that are more relevant to the inherent biological process
of the species. In addition to suggesting the taxonomy of newly discovered
species, the same model can be apply even to organisms that are already
classified, and suggest alternative taxonomy when building a newer model of
classification. The classifier can also be used to make sense of the taxonomic
rules, and can inform proposed modifications to the rules so that they are
more relevant to the biological pathways, and complement existing molecular
biology techniques in terms of species classification.

##  7 Conclusion

With our phylogenetic predictor, we are able to provide an automated model
that predict graph characteristics de novo. The model is unique in the way it
summarizes the statistics based solely on the characteristic of neighboring
nodes on the phylogeny, and can generate with reasonable high accuracy the
graph traits. As the construction of PPI networks for a new species is a
painstaking process and generally occur long after the species is well studies
in both evolutionary relationships and physiological traits, the predictor
model can serve as a useful informant for biologists both during the discovery
of newer species in cases where the PPI network is still incomplete for an
organism.

With the taxonomic classifier we constructed based on characteristics of the
protein networks of species, we are able to suggest taxonomic hierarchy of an
existing species with varying level of success in the different levels. This
is an interesting conclusion considering that the the taxonomy is determined
generally unrelated to the physiology and internal biochemical reactions of
species. The varying level of success, especially on the more granular level,
suggest alternative taxonomy can be in place based on the discoveries enabled
by our model. In addition, when a biologist can determine that network
statistics of an organism with a level of certainly, the classifier can
automatically suggest the taxonomic hierarchy of effectively especially on the
higher level.

With our work on both the phylogeny and taxonomy of biology, we are able to
connect the two fields through the bridge of the protein-protein interaction
networks. Even though the fields of taxonomy and phylogeny are similar in
their goal of making sense of the diversity of life, they employ vastly
different approaches and often draw contrasting conclusions on the
relationships between organisms. The fact that our two model connect the two
fields makes possible a all-encompassing predictor that predicts the taxonomy
of a species based on the known evolutionary relationships, even without prior
knowledge of the protein-protein interaction networks. This conclusion comes
to suggest an inherent relatedness between the fields of study, and we look
forward to exploring the success of our model at different taxonomic levels.
We also plan to engage in exploration in the other direction, going from the
known taxonomy of a species and try to deduce its phylogenetic relationships
towards others. Future work can also include coming up with better features
for representing position in the phylogeny tree, using deep learning models
for learning and investigating methods to improve accuracy of the classifier
at lower levels of taxonomy despite less training data.

## References

  * [1] J. P. Thorpe, “The molecular clock hypothesis: biochemical evolution, genetic differentiation and systematics,” _Annual Review of Ecology and Systematics_ , vol. 13, pp. 139–168, 1982.
  * [2] T. Domazet-Lošo, J. Brajković, and D. Tautz, “A phylostratigraphy approach to uncover the genomic history of major adaptations in metazoan lineages,” _Trends in Genetics_ , vol. 23, no. 11, pp. 533–539, 2007.
  * [3] E. C. Teeling, M. S. Springer, O. Madsen, P. Bates, S. J. O’brien, and W. J. Murphy, “A molecular phylogeny for bats illuminates biogeography and the fossil record,” _Science_ , vol. 307, no. 5709, pp. 580–584, 2005.
  * [4] Black. Taxonomy and phylogeny. [Online]. Available: http://webprojects.oit.ncsu.edu/project/bio181de/Black/
taxonomy/taxonomy.html

  * [5] E. Eisenberg and E. Y. Levanon, “Preferential attachment in the protein network evolution,” _Physical review letters_ , vol. 91, no. 13, p. 138701, 2003.
  * [6] A. Vázquez, “Growing network with local rules: Preferential attachment, clustering hierarchy, and degree correlations,” _Physical Review E_ , vol. 67, no. 5, p. 056104, 2003.
  * [7] J. M. Hancock and M. Simon, “Simple sequence repeats in proteins and their significance for network evolution,” _Gene_ , vol. 345, no. 1, pp. 113–118, 2005.
  * [8] R. Pastor-Satorras, E. Smith, and R. V. Solé, “Evolving protein interaction networks through gene duplication,” _Journal of Theoretical biology_ , vol. 222, no. 2, pp. 199–210, 2003.
  * [9] I. Ispolatov, P. Krapivsky, and A. Yuryev, “Duplication-divergence model of protein interaction network,” _Physical review E_ , vol. 71, no. 6, p. 061911, 2005.
  * [10] S. Wuchty, Z. N. Oltvai, and A.-L. Barabási, “Evolutionary conservation of motif constituents in the yeast protein interaction network,” _Nature genetics_ , vol. 35, no. 2, pp. 176–179, 2003.
  * [11] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon, “Network motifs: simple building blocks of complex networks,” _Science_ , vol. 298, no. 5594, pp. 824–827, 2002.
  * [12] J. Berg, M. Lässig, and A. Wagner, “Structure and evolution of protein interaction networks: a statistical model for link dynamics and gene duplications,” _BMC evolutionary biology_ , vol. 4, no. 1, p. 1, 2004.
  * [13] H. Jeong, S. P. Mason, A.-L. Barabási, and Z. N. Oltvai, “Lethality and centrality in protein networks,” _Nature_ , vol. 411, no. 6833, pp. 41–42, 2001.
  * [14] D. Bu, Y. Zhao, L. Cai, H. Xue, X. Zhu, H. Lu, J. Zhang, S. Sun, L. Ling, N. Zhang _et al._ , “Topological structure analysis of the protein–protein interaction network in budding yeast,” _Nucleic acids research_ , vol. 31, no. 9, pp. 2443–2450, 2003.
  * [15] C.-S. Liao, K. Lu, M. Baym, R. Singh, and B. Berger, “Isorankn: spectral methods for global alignment of multiple protein networks,” _Bioinformatics_ , vol. 25, no. 12, pp. i253–i258, 2009.
  * [16] B. P. Kelley, B. Yuan, F. Lewitter, R. Sharan, B. R. Stockwell, and T. Ideker, “Pathblast: a tool for alignment of protein interaction networks,” _Nucleic acids research_ , vol. 32, no. suppl 2, pp. W83–W88, 2004.
  * [17] M. Kalaev, V. Bafna, and R. Sharan, “Fast and accurate alignment of multiple protein networks,” _Journal of computational biology_ , vol. 16, no. 8, pp. 989–999, 2009.
  * [18] M. Kalaev, M. Smoot, T. Ideker, and R. Sharan, “Networkblast: comparative analysis of protein networks,” _Bioinformatics_ , vol. 24, no. 4, pp. 594–596, 2008.
  * [19] M. Koyutürk, Y. Kim, U. Topkara, S. Subramaniam, W. Szpankowski, and A. Grama, “Pairwise alignment of protein interaction networks,” _Journal of Computational Biology_ , vol. 13, no. 2, pp. 182–199, 2006.
  * [20] R. Singh, J. Xu, and B. Berger, “Pairwise global alignment of protein interaction networks by matching neighborhood topology,” in _Annual International Conference on Research in Computational Molecular Biology_. Springer, 2007, pp. 16–31.
  * [21] Y. Zhang and J. Skolnick, “Tm-align: a protein structure alignment algorithm based on the tm-score,” _Nucleic acids research_ , vol. 33, no. 7, pp. 2302–2309, 2005.
  * [22] R. Milo, S. Itzkovitz, N. Kashtan, R. Levitt, S. Shen-Orr, I. Ayzenshtat, M. Sheffer, and U. Alon, “Superfamilies of evolved and designed networks,” _Science_ , vol. 303, no. 5663, pp. 1538–1542, 2004.
  * [23] S. Maslov and K. Sneppen, “Specificity and stability in topology of protein networks,” _Science_ , vol. 296, no. 5569, pp. 910–913, 2002.
  * [24] O. Kuchaiev, T. Milenković, V. Memišević, W. Hayes, and N. Pržulj, “Topological network alignment uncovers biological function and phylogeny,” _Journal of The Royal Society Interface_ , vol. 7, no. 50, pp. 1341–1354, Sep. 2010.
  * [25] M. Gouy, S. Guindon, and O. Gascuel, “Seaview version 4: a multiplatform graphical user interface for sequence alignment and phylogenetic tree building,” _Molecular biology and evolution_ , vol. 27, no. 2, pp. 221–224, 2010.
  * [26] H. Shimodaira and M. Hasegawa, “Consel: for assessing the confidence of phylogenetic tree selection,” _Bioinformatics_ , vol. 17, no. 12, pp. 1246–1247, 2001.
  * [27] M. Van Oven and M. Kayser, “Updated comprehensive phylogenetic tree of global human mitochondrial dna variation,” _Human mutation_ , vol. 30, no. 2, pp. E386–E394, 2009.
  * [28] L. A. Hug, B. J. Baker, K. Anantharaman, C. T. Brown, A. J. Probst, C. J. Castelle, C. N. Butterfield, A. W. Hernsdorf, Y. Amano, K. Ise _et al._ , “A new view of the tree of life,” _Nature Microbiology_ , vol. 1, p. 16048, 2016.
  * [29] S. Federhen, “The ncbi taxonomy database,” _Nucleic acids research_ , vol. 40, no. D1, pp. D136–D143, 2012.
  * [30] A. Vázquez, A. Flammini, A. Maritan, and A. Vespignani, “Modeling of protein interaction networks,” _Complexus_ , vol. 1, no. 1, pp. 38–44, 2002\.


* [83] [author] van der Vaart, A. W.A. W. and Wellner, J. A.J. A. (1996). Weak Convergence and Empirical Processes: With Applications to Statistics. Springer.
  * [84] [author] van der Vaart, A. W.A. W. (1998). Asymptotic Statistics. Cambridge University Press, Cambridge, UK.
  * [85] [author] van Handel, R.R. (2018). Probability in high dimension. Cambridge University Press.
  * [86] [author] var der Vaart, A. W.A. W. (1996). New Donsker classes. Annals of Probability 24 2128-2124.
  * [87] [author] Villani, C.C. (2008). Optimal Transport: Old and New 338\. Springer Science & Business Media.
  * [88] [author] Rhee, W. S.W. S. and Talagrand, M.M. (1986). Uniform bound in the central limit theorem for Banach space valued dependent random variables. Journal of multivariate analysis 20 303–320.
  * [89] [author] Weed, J.J. and Bach, F.F. (2019). Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance. Bernoulli 25 2620–2648.
  * [90] [author] Yurochkin, M.M., Claici, S.S., Chien, E.E., Mirzazadeh, F.F. and Solomon, J.J. (2019). Hierarchical Optimal Transport for Document Representation. arXiv preprint arXiv:1906.10827.
  * [91] [author] Zhang, P.P., Liu, Q.Q., Zhou, D.D., Xu, T.T. and He, X.X. (2018). On the discrimination-generalization tradeoff in GANs. In Proceedings of the International Conference on Learning Representations (ICLR-2018).
  * [92] [author] Zolotarev, V. M.V. M. (1983). Probability metrics. Teoriya Veroyatnostei i ee Primeneniya 28 264–287.

##  Appendix A Supplementary Proofs for Section 3

###  A.1 Supplementary Proofs for Theorem 1

We start from the following auxiliary lemma.

######  Lemma 16 (Uniform bound on derivatives).

For any $f\in\mathsf{Lip}_{1}$ and any nonzero multi-index
$k=(k_{1},\dots,k_{d})$, we have

$\big{|}D^{k}f_{\sigma}(x)\big{|}\leq\sigma^{-|k|+1}\sqrt{(|k|-1)!},\quad\forall
x\in\mathbb{R}^{d}.$

###### Proof.

Let $H_{m}(z)$ denote the Hermite polynomial of degree $m$ defined by

$H_{m}(z)=(-1)^{m}e^{z^{2}/2}\left[\frac{d^{m}}{dz^{m}}e^{-z^{2}/2}\right],\
m=0,1,\dots.$

Recall that for $Z\sim\mathcal{N}(0,1)$, we have
$\mathbb{E}[H_{m}(Z)^{2}]=m!$. A straightforward computation shows that

$D_{x}^{k}\varphi_{\sigma}(x-y)=\varphi_{\sigma}(x-y)\left[\prod_{j=1}^{d}(-1)^{k_{j}}\sigma^{-k_{j}}H_{k_{j}}\big{(}(x_{j}-y_{j})/\sigma\big{)}\right]$

for any multi-index $k=(k_{1},\dots,k_{d})$, where $D_{x}$ means that the
differential operator is applied to $x$. Hence, we have

$\begin{split}D^{k}f_{\sigma}(x)&=\int
f(y)\varphi_{\sigma}(x-y)\left[\prod_{j=1}^{d}(-1)^{k_{j}}\sigma^{-k_{j}}H_{k_{j}}\big{(}(x_{j}-y_{j})/\sigma\big{)}\right]\,\mathrm{d}y\\\
&=\int f(x-\sigma
y)\varphi_{1}(y)\left[\prod_{j=1}^{d}(-1)^{k_{j}}\sigma^{-k_{j}}H_{k_{j}}(y_{j})\right]\,\mathrm{d}y,\end{split}$

so that, by $1$-Lipschitz continuity of $f$,

$\left|D^{k}f_{\sigma}(x)-D^{k}f_{\sigma}(x^{\prime})\right|\leq\|x-x^{\prime}\|\int\varphi_{1}(y)\left[\prod_{j=1}^{d}\sigma^{-k_{j}}\big{|}H_{k_{j}}(y_{j})\big{|}\right]\,\mathrm{d}y.$

The integral on the RHS equals

$\prod_{j=1}^{d}\sigma^{-k_{j}}\mathbb{E}\Big{[}\big{|}H_{k_{j}}(Z)\big{|}\Big{]}\leq\prod_{j=1}^{d}\sigma^{-k_{j}}\sqrt{\mathbb{E}\left[\big{|}H_{k_{j}}(Z)\big{|}^{2}\right]}=\prod_{j=1}^{d}\sigma^{-k_{j}}\sqrt{k_{j}!}\leq\sigma^{-|k|}\sqrt{|k|!},$

where $Z\sim\mathcal{N}(0,1)$. The conclusion of the lemma follows from
induction on the size of $|k|$. ∎

###### Proof of Lemma 4.

For $\alpha=1$,

$\|f_{\sigma}\|_{\alpha,\mathcal{X}}\leq\sup_{\mathcal{X}}\big{(}\|x\|+\sigma\sqrt{d}\big{)}\vee
1\lesssim_{d}1\vee\sup_{\mathcal{X}}\|x\|.$

For $\alpha>1$,

$\displaystyle\sup_{\overset{|k|=\underline{\alpha}}{x,y\in\mathcal{X}}}\mspace{-4.0mu}\frac{\big{|}D^{k}\mspace{-3.0mu}f_{\sigma}(x)\mspace{-3.0mu}-D^{k}\mspace{-3.0mu}f_{\sigma}(y)\big{|}}{\|x-y\|^{\alpha-\underline{\alpha}}}$
$\displaystyle\leq\mspace{-3.0mu}\left(\sup_{\overset{|k|=\underline{\alpha}}{x,y\in\mathcal{X}}}\mspace{-4.0mu}\frac{\big{|}D^{k}\mspace{-3.0mu}f_{\sigma}(x)\mspace{-2.0mu}-\mspace{-2.0mu}D^{k}\mspace{-3.0mu}f_{\sigma}(y)\big{|}^{\alpha-\underline{\alpha}}}{\|x-y\|^{\alpha-\underline{\alpha}}}\right)\mspace{-3.0mu}\sup_{\overset{|k|=\underline{\alpha}}{x\in\mathcal{X}}}\mspace{-4.0mu}\big{|}D^{k}\mspace{-3.0mu}f_{\sigma}(x)\mspace{-3.0mu}-\mspace{-3.0mu}D^{k}\mspace{-3.0mu}f_{\sigma}(y)\big{|}^{1-\alpha+\underline{\alpha}}$
$\displaystyle\leq
2\left(\sup_{\overset{|k|=\underline{\alpha}+1}{x\in\mathcal{X}}}\big{|}D^{k}\mspace{-3.0mu}f_{\sigma}(x)\big{|}\right)^{\alpha-\underline{\alpha}}\sup_{\overset{|k|=\underline{\alpha}}{x\in\mathcal{X}}}\big{|}D^{k}\mspace{-3.0mu}f_{\sigma}(x)\big{|}^{1-\alpha+\underline{\alpha}}$
$\displaystyle\leq(\sigma^{-\underline{\alpha}}\sqrt{\underline{\alpha}!})^{\alpha-\underline{\alpha}}(\sigma^{-\underline{\alpha}+1}\sqrt{(\underline{\alpha}-1)!})^{1-\alpha+\underline{\alpha}}$
$\displaystyle=\sigma^{-\alpha+1}\big{(}\sqrt{\underline{\alpha}!}\big{)}^{\alpha-\underline{\alpha}}\big{(}\sqrt{(\underline{\alpha}-1)!}\big{)}^{1-\alpha+\underline{\alpha}}.$

Hence,

$\displaystyle\|f_{\sigma}\|_{\alpha,\mathcal{X}}$
$\displaystyle\leq\max_{0\leq|k|\leq\underline{\alpha}}\sup_{x\in\mathcal{X}}\left|D^{k}\mspace{-1.0mu}f_{\sigma}(x)\right|+C_{\alpha}\sigma^{-\alpha+1}$
$\displaystyle\leq\left((\sup_{\mathcal{X}}\|x\|+\sigma\sqrt{d})\vee\sigma^{-\underline{\alpha}+1}\sqrt{(\underline{\alpha}-1)!}\right)+C_{\alpha}\sigma^{-\alpha+1}$
$\displaystyle\lesssim_{\alpha,d}\sigma^{-\alpha+1}\big{(}1\vee\sup_{\mathcal{X}}\|x\|\big{)},$

completing the proof. ∎

###  A.2 Proof of Corollary 1

The following version of Dudley’s entropy integral bound will be used in the
proof of Corollary 1.

######  Lemma 17 (Expectation bound for non-integrable entropy).

For any bounded function class $\mathcal{F}$,

$\mathbb{E}[\|\mathbb{G}_{n}\|_{\mathcal{F}}]\lesssim\mathbb{E}\left[\inf_{\gamma>0}\left\\{\gamma\sqrt{n}+\int_{\gamma}^{\mspace{1.0mu}\sigma_{n}}\sqrt{\log
N(\mathcal{F},\|\cdot\|_{L^{2}(P_{n})},\delta)}\
\,\mathrm{d}\delta\right\\}\right],$

where $\sigma_{n}=\sup_{\mathcal{F}}\|f\|_{L^{2}(P_{n})}$.

This lemma follows by applying [85, Theorem 5.31] after symmetrization and
conditioning on $X_{1},\dots,X_{n}$.

###### Proof of Corollary 1.

We follow the notation used in the proof of Theorem 1. If $\alpha=d/2$, then
by Lemma 17,

$\begin{split}\mathbb{E}\big{[}\|\mathbb{G}_{n}\|_{{\mathcal{F}_{\sigma}}_{j}}\big{]}&\lesssim\mathbb{E}\left[\inf_{\gamma}\left\\{\gamma\sqrt{n}+\int_{\gamma}^{M_{j}^{\prime}P_{n}(I_{j})^{1/2}}\frac{M_{j}^{\prime}P_{n}(I_{j})^{1/2}}{\epsilon}\,\mathrm{d}\epsilon\right\\}\right]\\\
&\leq\mathbb{E}\left[M_{j}^{\prime}P_{n}(I_{j})^{1/2}+M_{j}^{\prime}P_{n}(I_{j})^{1/2}\log\sqrt{n}\right]\\\
&\lesssim M_{j}^{\prime}P(I_{j})^{1/2}\log n\\\
&\lesssim_{d}\sigma^{-d/2+1}M_{j}P(I_{j})^{1/2}\log n,\end{split}$

where the second inequality follows by choosing
$\gamma=M_{j}^{\prime}P_{n}(I_{j})^{1/2}/\sqrt{n}$. The desired moment bound
follows by summing the final expression over $j$. ∎

###  A.3 Proof of Lemma 1 and Corollary 2

###### Proof of Lemma 1.

From the proof of Theorem 1 and the fact that $\mathsf{Lip}_{1}$ is symmetric,
we have $L_{P}^{\mspace{1.0mu}\sigma}=\|G_{P}\|_{{\mathcal{F}_{\sigma}}}$.
Since $G_{P}$ is a tight Gaussian process in
$\ell^{\infty}({\mathcal{F}_{\sigma}})$, ${\mathcal{F}_{\sigma}}$ is totally
bounded for the pseudometric
$d_{P}(f,g)=\sqrt{\mathsf{Var}_{P}(f_{\sigma}-g_{\sigma})}$, and $G_{P}$ is a
Borel measurable map into the space of $d_{P}$-uniformly continuous functions
$\mathcal{C}_{u}({\mathcal{F}_{\sigma}})$ equipped with the uniform norm
$\|\cdot\|_{{\mathcal{F}_{\sigma}}}$. Let $F$ denote the distribution function
of $L_{P}^{\mspace{1.0mu}\sigma}$, and define $r_{0}:=\inf\\{r\geq
0:F(r)>0\\}$.

From [26, Theorem 11.1], $F$ is absolutely continuous on $(r_{0},\infty)$, and
there exists a countable set $\Delta\subset(r_{0},\infty)$ such that
$F^{\prime}$ is positive and continuous on $(r_{0},\infty)\setminus\Delta$.
The theorem however does not exclude the possibility that $F$ has a jump at
$r_{0}$, and we will verify that (i) $r_{0}=0$ and (ii) $F$ has no jump at
$r=0$, which lead to the conclusion. The former follows from p. 57 in [57].
The latter is trivial since

$F(0)-F(0-)=\mathbb{P}\left(L_{P}^{\mspace{1.0mu}\sigma}=0\right)\leq\mathbb{P}\big{(}G_{P}(f)=0\big{)},\quad\forall
f\in{\mathcal{F}_{\sigma}}.$

As $G_{P}$ is Gaussian, we have $\mathbb{P}\big{(}G_{P}(f)=0\big{)}=0$, for
any $f$ such that $\mathsf{Var}_{P}(f)>0$. If $P$ is not a point mass, there
always exists such an $f\in\mathcal{F}_{\sigma}$. ∎

###### Proof of Corollary 2.

From Theorem 3.6.3 in [83] applied to the function class
${\mathcal{F}_{\sigma}}$ together with the continuous mapping theorem, we see
that conditionally on $X_{1},X_{2},\dots$, one has

$\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}\big{(}P_{n}^{B},P_{n}\big{)}=\big{\|}\sqrt{n}(P_{n}^{B}-P_{n})\big{\|}_{\mathcal{F}_{\sigma}}\stackrel{{\scriptstyle
w}}{{\to}}L_{P}^{\mspace{1.0mu}\sigma}$

for almost every realization of $X_{1},X_{2},\dots$ The conclusion follows
from the fact that the distribution function of $L_{P}^{\mspace{1.0mu}\sigma}$
is continuous (cf. Lemma 1) and Polya’s Theorem (cf. Lemma 2.11 in [84]). ∎

###  A.4 Proof of Proposition 1

Case (i) follows from Talagrand’s inequality (see Lemma 1 in [2]). Cases (ii)
and (iii) follow from Theorems 4 and 2 in [3] and [2], respectively, applied
to the function class ${\mathcal{F}_{\sigma}}$ using the envelope function
${F}(x)=\|x\|+\sigma\sqrt{d}$. We omit the details for brevity. ∎

###  A.5 Proof of Corollary 3

If $\mathcal{F}_{\sigma}$ is $P$-Donsker, $Q=P$, and $m/(m+n)\to\lambda$, then
$W_{m,n}$ converges in distribution to the supremum of the tight Gaussian
process $G_{P}$ over $\mathcal{F}_{\sigma}$ [83, Chapter 3.7]. The bootstrap
convergence follows by [86, Theorem 3.7.7]. Lastly, since the distribution of
$L_{H}^{\mspace{1.0mu}\sigma}$ is also Lebesgue absolutely continuous under
the above conditions (Lemma 1) and $H=\lambda P+(1-\lambda)Q=P$ when $P=Q$,
the final statement follows. ∎

##  Appendix B Supplementary Proofs for Section 4

###  B.1 Supplementary Proofs for Theorem 2

###### Proof of Lemma 7.

See Theorem 2.7.4 in [83]. ∎

###### Proof of Lemma 10.

For $I_{j}\subset\big{\\{}x\,:\,\|x\|_{\infty}\in(k-1,k]\big{\\}}=:L_{k}$,
$k\in\mathbb{N}$, we have $\sup_{I_{j}}\|x\|\leq k\sqrt{d}$. Now,

$\sum_{j=1}^{\infty}\left(\sup\nolimits_{I_{j}}\|x\|\right)^{p_{1}}P(I_{j})^{p_{2}}=\sum_{k=1}^{\infty}\sum_{j:I_{j}\subset
L_{k}}\left(\sup\nolimits_{I_{j}}\|x\|\right)^{p_{1}}P(I_{j})^{p_{2}}\lesssim_{d}\sum_{k=1}^{\infty}k^{p_{1}}P_{k},$

where $P_{k}:=\sum_{j:I_{j}\subset L_{k}}P(I_{j})^{p_{2}}$. Setting
$l_{k}:=\operatorname{Card}\\{I_{j}\,:\,I_{j}\subset L_{k}\\}$, we have

$l_{k}=(2k)^{d}-(2k-2)^{d}\lesssim_{d}k^{d-1},$

where the last inequality follows since the expression inside the bracket is a
sequence that converges to $1$, and hence is bounded. By the power mean
inequality

$\left(\sum_{j:I_{j}\subset
L_{k}}\frac{P(I_{j})^{p_{2}}}{l_{k}}\right)^{\frac{1}{p_{2}}}\leq\frac{P(L_{k})}{l_{k}}\quad\implies\quad
P_{k}\leq
l_{k}^{1-p_{2}}P(L_{k})^{p_{2}}\lesssim_{d}k^{(d-1)(1-p_{2})}P(L_{k})^{p_{2}}.$

Hence, a sufficient condition for $\sum_{k=1}^{\infty}k^{p_{1}}P_{k}<\infty$
follows by observing that

$\sum_{k=1}^{\infty}k^{p_{1}+(d-1)(1-p_{2})}P(L_{k})^{p_{2}}\lesssim\int_{1}^{\infty}t^{p_{1}+(d-1)(1-p_{2})}\mathbb{P}(\|X\|_{\infty}>t)^{p_{2}}\
\,\mathrm{d}t.$

If $\mathbb{E}[\|X\|_{\infty}^{l}]=M<\infty$ for any
$l>\frac{(d-1)(1-p_{2})+p_{1}+1}{p_{2}}$, then Markov’s inequality implies
that the above integral is bounded by $M^{p_{2}}$ up to a constant depending
on $l$, $p_{1}$, $p_{2}$, and $d$. ∎

###  B.2 Proof of Theorem 3

Recall the decomposition from (7.3). Term (I) has already been bounded in
(27). We bound terms (II)-(IV) in the following steps. Fix $n\in\mathbb{N}$
and $\sigma,\delta\in(0,1]$.

Step 1 (Discretization error of supremized empirical bootstrap (II)). By [83,
Lemma 3.6.6], we have,

$\mathbb{E}\left[\|\mathbb{G}_{n}^{*}\|_{\mathcal{F}_{\sigma}(\delta)}|\middle|\,X_{1:n}\right]\lesssim\frac{1}{\sqrt{n}}\mathbb{E}\left[\left\|\sum_{i=1}^{n}(N_{i}-N^{\prime}_{i})\delta_{X_{i}}\right\|_{\mathcal{F}_{\sigma}(\delta)}\middle|\,X_{1:n}\right],$

where $N_{i},N^{\prime}_{i}$ are i.i.d. Poisson(1/2) random variables,
independent of $X_{1:n}$. Consider the modified function class:

$\overline{\mathcal{F}_{\sigma}(\delta)}:=\\{(x,w)\mapsto
wf(x):f\in\mathcal{F}_{\sigma}(\delta),w\in\mathbb{Z}\\}.$

With some abuse of notation, let $\mathbb{G}_{n}$ denote also the empirical
process corresponding to the product probability law
$P\otimes\mathsf{Law}(N_{i}-N^{\prime}_{i})$. Then, for the marginal
expectation of the discretization error above, we have

$\mathbb{E}\left[\|\mathbb{G}_{n}^{*}\|_{\mathcal{F}_{\sigma}(\delta)}\right]\lesssim\frac{1}{\sqrt{n}}\mathbb{E}\left[\left\|\sum_{i=1}^{n}(N_{i}-N^{\prime}_{i})\delta_{X_{i}}\right\|_{\mathcal{F}_{\sigma}(\delta)}\right]=\mathbb{E}\left[\|\mathbb{G}_{n}\|_{\overline{\mathcal{F}_{\sigma}(\delta)}}\right].$

To bound this marginal expectation, we will apply local maximal inequality on
restrictions of process to compact partitioning sets (as in Step 2 of the
proof of Theorem 2) and sum up the resulting bounds. Now, let
$\overline{\mathcal{F}_{j,k}(\delta)}:=\big{\\{}(x,w)\mapsto
f(x,w)\mathds{1}_{I_{j}\times\\{k\\}}(x,|w|):j\in\mathbb{N},k\in\mathbb{N}\big{\\}}$,
where $I_{j}$ are an enumeration of unit cubes on integer lattice points in
$\mathbb{R}^{d}$. Note that
$F_{j,k}(x)=kM^{\prime}_{j}\mathds{1}_{I_{j}\times\\{k\\}}(x,|w|)$ is an
envelope function for this class. As in (16), the entropy integral is bounded
by

$\begin{split}J\big{(}\delta,\overline{\mathcal{F}_{j,k}(\delta)},F_{j,k}\big{)}&:=\int_{0}^{\delta}\sup_{Q\in\mathcal{P}_{\mathsf{f}}(\mathbb{R}^{d+1})}\sqrt{\log
N\left(\epsilon\|2F_{j}\|_{L^{2}(Q)},\overline{\mathcal{F}_{j,k}(\delta)},L^{2}(Q)\right)}\,\mathrm{d}\epsilon\\\
&\lesssim_{\alpha,d}\int_{0}^{\delta}\epsilon^{-b}\,\mathrm{d}\epsilon\lesssim_{\alpha,d}\delta^{1-b},\end{split}$
(31)

which, via [22, Theorem 5.2], yields

$\displaystyle\mathbb{E}\big{[}\|\mathbb{G}_{n}^{*}\|_{\mathcal{F}_{\sigma}(\delta)}\big{]}\lesssim\mathbb{E}\big{[}\|\mathbb{G}_{n}\|_{\overline{\mathcal{F}_{\sigma}(\delta)}}\big{]}$
$\displaystyle\lesssim_{\alpha,d}\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}\left[\big{(}kM^{\prime}_{j}p_{k}^{\frac{1}{2}}P(I_{j})^{\frac{1}{2}}\big{)}^{b}\delta^{1-b}\mspace{-3.0mu}+\mspace{-3.0mu}\frac{(kM^{\prime}_{j})^{1+2b}p_{k}^{b}P(I_{j})^{b}}{\delta^{2b}\sqrt{n}}\right]$
$\displaystyle\lesssim_{\alpha,d}\sum_{j=1}^{\infty}\left(M^{\prime}_{j}P(I_{j})^{\frac{1}{2}}\right)^{b}\delta^{1-b}+\frac{1}{\sqrt{n}}(M^{\prime}_{j})^{1+2b}P(I_{j})^{b}\delta^{-2b}$
$\displaystyle\lesssim_{\alpha,d,M,L}\sigma^{-ab}\delta^{1-b}+\frac{\sigma^{-a(1+2b)}\delta^{-2b}}{\sqrt{n}}.$
(32)

Step 2 (Comparison between $\mathbb{G}_{n}^{*}$ and $\mathbb{G}_{n}^{B}$ after
discretization (III)). We now compare the discretized empirical and multiplier
bootstrap processes conditioned on $X_{1:n}$, using Lemma 8. To that end, let
$f_{1},\dots,f_{N_{\sigma,\delta}}$ be an $L^{2}(P)$ $\delta$-net of
$\mathcal{F}_{\sigma}$. Conditioned on $X_{1:n}$, for a bootstrap resample
$X_{1}^{B},\ldots,X_{n}^{B}$, the random vectors
$\big{(}f_{1}(X_{i}^{B})-P_{n}f,\ldots,f_{N_{\sigma,\delta}}(X_{i}^{B})-P_{n}f_{N_{\sigma,\delta}}\big{)}$,
$1\leq i\leq n$, are i.i.d. with conditional mean zero and conditional
covariances given by
$\mathsf{Cov}\big{(}f_{j}(X_{i}^{N}),f_{k}(X_{i}^{B})\,\big{|}\,X_{1:n}\big{)}=\mathsf{Cov}_{P_{n}}(f_{j},f_{k})$.
For the multiplier bootstrap, conditioned on $X_{1:n}$, the distribution of
$\mathbb{G}_{n}^{B}(f),f\in\mathcal{F}_{\sigma}$ is that of a Gaussian process
with covariance function
$\mathsf{Cov}\big{(}\mathbb{G}_{n}^{B}(f),\mathbb{G}_{n}^{B}(g)\big{)}=\mathsf{Cov}_{P_{n}}(f,g)$.
Consequently,
$\big{(}\mathbb{G}_{n}^{B}(f_{1}),\dots,\mathbb{G}_{n}^{B}(f_{N_{\sigma,\delta}})\big{)}$
and $\big{(}f_{1}(X_{i}^{B}),\dots,f_{N_{\sigma,\delta}}(X_{i}^{B})\big{)}$
have the same conditional covariance matrix.

The terms $L_{n}$ and $M_{n,X}(\eta)$ in Lemma 8 can be bounded by
$(2/n)\sum_{i=1}^{n}\big{(}\|X_{i}\|+\sigma\sqrt{d}\big{)}^{3}$. We next bound
$M_{n,Y}(\eta)$. Here, conditioned on $X_{1:n}$,
$Y_{j}:=n^{-1/2}\sum_{i=1}^{n}\eta_{i}\big{(}f_{j}(X_{i})-P_{n}f_{j}\big{)}$
are $\mathcal{N}(0,\sigma_{j}^{2})$, where $\sigma_{j}^{2}\leq
n^{-1}\sum_{i=1}^{n}\big{(}\|X_{i}\|+\sigma\sqrt{d}\big{)}^{2}$, for all
$1\leq j\leq p$. Then, we have the following concentration inequality and
moment bound for Gaussian random variables [41, Example 2.1.19]:

$\displaystyle\mathbb{P}\left(\max_{1\leq j\leq
N_{\sigma,\delta}}|Y_{j}|>\mathbb{E}\left[\max_{1\leq j\leq
N_{\sigma,\delta}}|Y_{j}\,|\middle|X_{1:n}\right]+t\,\middle|X_{1:n}\right)$
$\displaystyle\leq\exp\left(-\frac{t^{2}}{\frac{2}{n}\sum_{i=1}^{n}\big{(}\|X_{i}\|+\sigma\sqrt{d}\big{)}^{2}}\right),$
$\displaystyle\mathbb{E}\left[\max_{1\leq j\leq
N_{\sigma,\delta}}|Y_{j}|\,\middle|X_{1:n}\right]$
$\displaystyle\leq\sqrt{\frac{2}{n}\sum_{i=1}^{n}\big{(}\|X_{i}\|+\sigma\sqrt{d}\big{)}^{2}\log(2N_{\sigma,\delta})}.$

To bound $M_{n,Y}(\eta)$, we integrate the above tail probability bound. To
that end, consider the event

$E:=\left\\{\frac{\eta\sqrt{n}}{\log
N_{\sigma,\delta}}\geq\sqrt{\frac{2}{n}\sum_{i=1}^{n}\big{(}\|X_{i}\|+\sigma\sqrt{d}\big{)}^{2}\log(2N_{\sigma,\delta})}\right\\},$

and note that on $E$, $M_{n,Y}(\eta)\leq
C\big{(}(1/n)\sum_{i=1}^{n}(\|X_{i}\|+\sigma\sqrt{d})^{2}\big{)}^{3/2}$.
Taking $C_{d}=C_{8}(1\vee C)$, this implies

$\begin{split}\mathbb{P}\left(\left\\{\max_{1\leq j\leq
N_{\sigma,\delta}}\mathbb{G}_{n}^{*}(f_{i})\in A^{\eta}\right\\}\cap
E\middle|X_{1:n}\right)&\leq\mathbb{P}\left(\max_{1\leq j\leq
N_{\sigma,\delta}}\mathbb{G}_{n}^{B}(f_{i})\in
A^{(1+C_{7})\eta}\middle|X_{1:n}\right)\\\ &\qquad\qquad+\frac{C_{d}(\log
N_{\sigma,\delta})^{2}}{\eta^{3}\sqrt{n}}\left[\frac{3}{n}\sum_{i=1}^{n}(\|X_{i}\|+\sigma\sqrt{d})^{3}\right].\end{split}$

Adding a term for $E^{c}$, we further obtain

$\begin{split}\mathbb{P}\left(\max_{1\leq j\leq
N_{\sigma,\delta}}\mathbb{G}_{n}^{*}(f_{i})\in
A^{\eta}\middle|X_{1:n}\right)&\leq\mathbb{P}\left(\max_{1\leq j\leq
N_{\sigma,\delta}}\mathbb{G}_{n}^{B}(f_{i})\in
A^{(1+C_{7})\eta}\middle|X_{1:n}\right)\\\ &\qquad+\frac{C_{d}(\log
N_{\sigma,\delta})^{2}}{\eta^{3}\sqrt{n}}\left[\frac{3}{n}\sum_{i=1}^{n}(\|X_{i}\|+\sigma\sqrt{d})^{3}\right]+\mathbb{P}\left(E^{c}\middle|\,X_{1:n}\right).\end{split}$
(33)

Finally, observe that $\mathbb{P}(E^{c}|X_{1:n})=\mathds{1}_{E^{c}}$ since $E$
is $\sigma(X_{1:n})$-measurable, and that Markov’s inequality yields

$\mathbb{P}(E^{c})\leq\mathbb{P}\left(\frac{1}{n}\sum_{i=1}^{n}\left(\|X_{i}\|+\sqrt{d}\right)^{2}>\frac{n\eta^{2}}{2(\log(2N_{\sigma,\delta}))^{3}}\right)\lesssim_{\alpha,d,M,l}\frac{(\log(2N_{\sigma,\delta}))^{3}}{n\eta^{2}}.$
(34)

Step 3 (Comparison between $\mathbb{G}_{n}^{B}$ and $G_{P}$ after
discretization (IV)). Let $Z=\max_{1\leq i\leq
N_{\sigma,\delta}}G_{P}(f_{i})$, and $\widetilde{Z}=\max_{1\leq i\leq
N_{\sigma,\delta}}\mathbb{G}_{n}^{B}(f_{i})$. By Lemma 11, for all $\eta>0$,
and a constant $C_{9}$ depending only on $d$,

$\mathbb{P}\big{(}\widetilde{Z}\in
A^{(1+C_{7})\eta}\big{|}X_{1:n}\big{)}\leq\mathbb{P}\big{(}Z\in
A^{(2+C_{7})\eta}\big{)}+C_{9}\eta^{-1}\sqrt{\Delta(X_{1:n})\log
N_{\sigma,\delta}},$ (35)

where $\Delta(X_{1:n})=\max_{1\leq j,k\leq
N_{\sigma,\delta}}\big{|}\Sigma_{j,k}^{X}(X_{1:n})-\Sigma_{j,k}^{Y}\big{|}$
with $\Sigma^{X}(X_{1:n})$ as the conditional covariance matrix of
$\big{(}\mathbb{G}_{n}^{B}(f_{1}),\dots,\mathbb{G}_{n}^{B}(f_{N_{\sigma,\delta}})\big{)}$,
and $\Sigma_{Y}$ the covariance matrix of
$\big{(}G_{P}(f_{1}),\dots,G_{P}(f_{N_{\sigma,\delta}})\big{)}$. To bound
$\Delta(X_{1:n})$, note that

$\displaystyle\Sigma_{j,k}^{Y}\mspace{-3.0mu}-\mspace{-3.0mu}\Sigma_{j,k}^{X}(X_{1:n})$
$\displaystyle=\mathsf{Cov}_{P}(f_{j},f_{k})\mspace{-3.0mu}-\mspace{-3.0mu}\mathsf{Cov}_{P_{n}}(f_{j},f_{k})$
$\displaystyle=Pf_{j}f_{k}\mspace{-3.0mu}-\mspace{-3.0mu}Pf_{j}Pf_{k}\mspace{-3.0mu}-\mspace{-3.0mu}P_{n}f_{j}f_{k}+P_{n}f_{j}P_{n}f_{k}$
$\displaystyle=\mspace{-3.0mu}-\mspace{-3.0mu}(P_{n}\mspace{-3.0mu}-\mspace{-3.0mu}P)f_{j}f_{k}\mspace{-3.0mu}+\mspace{-3.0mu}(P_{n}\mspace{-3.0mu}-\mspace{-3.0mu}P)f_{j}(P_{n}\mspace{-3.0mu}-\mspace{-3.0mu}P)f_{k}\mspace{-3.0mu}+\mspace{-3.0mu}Pf_{k}(P_{n}\mspace{-3.0mu}-\mspace{-3.0mu}P)f_{j}\mspace{-3.0mu}+\mspace{-3.0mu}Pf_{j}(P_{n}\mspace{-3.0mu}-\mspace{-3.0mu}P)f_{k}.$

Consider the class $(\mathcal{F}_{\sigma})^{\times
2}=\\{fg:f,g\in\mathcal{F}_{\sigma}\\}$, and denote its restriction to the
cube $I_{j}$ by $(\mathcal{F}_{j})^{\times 2}$. With this notation, we have
$\mathbb{P}$-a.s.

$\Delta(X_{1:n})\leq\sum_{j=1}^{\infty}\frac{1}{\sqrt{n}}\left[\|\mathbb{G}_{n}\|_{(\mathcal{F}_{j})^{\times
2}}+\|\mathbb{G}_{n}\|_{\mathcal{F}_{j}}^{2}+2\Big{(}\sup_{\mathcal{F}_{j}}P|f|\Big{)}\|\mathbb{G}_{n}\|_{\mathcal{F}_{j}}\right].$
(36)

Given (36), we bound $\mathbb{E}\big{[}\Delta(X_{1:n})\big{]}$ using empirical
process techniques. We start by controlling the metric entropy of
$(\mathcal{F}_{j})^{\times 2}$. Recall that from Lemma 5 we have

$\log N\big{(}\epsilon
M,C_{M}^{N}(\chi),\|\cdot\|_{\infty}\big{)}\lesssim_{d,N,\text{diam}(\chi)}\epsilon^{-d/N},$

where $C_{M}^{N}(\chi)$ is the Hölder ball of smoothness $N$ and radius $M$.
For $(\mathcal{F}_{\sigma})^{\times 2}$, we take $N=\alpha$ and
$M=2^{\lceil\alpha\rceil}(M^{\prime}_{j})^{2}$, and for any finitely supported
measure $Q$ obtain

$\log
N\big{(}\epsilon(M^{\prime}_{j})^{2}Q(I_{j})^{\frac{1}{2}},(\mathcal{F}_{j})^{\times
2},L^{2}(Q)\big{)}\lesssim_{d,\alpha}\epsilon^{-\frac{d}{\alpha}}.$

Then, using Theorem 2.14.1 in [83], the expectation of the first term in the
RHS of (36) is bounded as

$\mathbb{E}\big{[}\|\mathbb{G}_{n}\|_{(\mathcal{F}_{\sigma})^{\times
2}}\big{]}=\sum_{j=1}^{\infty}\mathbb{E}\big{[}\|\mathbb{G}_{n}\|_{(\mathcal{F}_{j})^{\times
2}}\big{]}\lesssim_{\alpha,d}\sum_{j=1}^{\infty}(M^{\prime}_{j})^{2}P(I_{j})^{\frac{1}{2}}.$

For the second term, (10) combined with [83, Theorem 2.14.5] implies

$\mathbb{E}\left[\|\mathbb{G}_{n}\|_{\mathcal{F}_{j}}^{2}\right]\lesssim_{\alpha,d}\sum_{j=1}^{\infty}(M^{\prime}_{j})^{2}P(I_{j}),$

while the third term can be bounded by (10). Combined, this gives

$\mathbb{E}\big{[}\Delta(X_{1:n})\big{]}\lesssim_{\alpha,d}\frac{1}{\sqrt{n}}\left[\sum_{j=1}^{\infty}(M^{\prime}_{j})^{2}P(I_{j})^{\frac{1}{2}}+\sum_{j=1}^{\infty}(M^{\prime}_{j})^{2}P(I_{j})\right]\lesssim_{\alpha,d,M,l}\sigma^{-2a}n^{-\frac{1}{2}}.$
(37)

Step 4 (Conclusion using Strassens’s theorem). By the previous steps, we have
the following bounds on the terms in RHS of (7.3). First, by the Borel-
Sudakov-Tsirelson inequality [41, Theorem 2.2.7], the Gaussian process
discretization error is bounded by (see (27))

$\text{(I)}\leq\mathbb{P}\Big{(}\|G_{P}\|_{\mathcal{F}_{\sigma}(\delta)}>\eta\Big{)}\leq
e^{-\frac{1}{2\delta^{2}}\big{(}\eta-\mathbb{E}[\|G_{P}\|_{\mathcal{F}_{\sigma}(\delta)}]\big{)}^{2}},$

which holds for every $\eta\geq
C_{\alpha,d,M,l}\sigma^{-ab}\delta^{1-b}\geq\mathbb{E}\big{[}\|G_{P}\|_{\mathcal{F}_{\sigma}(\delta)}\big{]}$
by (19). Second, by Markov’s inequality and (B.2), we have the following bound
for the empirical bootstrap process discretization error:

$\mathbb{E}\big{[}\text{(II)}\big{]}\lesssim_{\alpha,d,M,l}\frac{1}{\eta}\left[\sigma^{-ab}\delta^{1-b}+\frac{\sigma^{-a(1+2b)}\delta^{-2b}}{\sqrt{n}}\right].$

Third, starting from (33), we bound $\log N_{\sigma,\delta}$ and
$\mathbb{P}(E^{c})$ using (14) and (34), respectively. This yield the
following inequality for the comparison between empirical and multiplier
bootstraps:

$\mathbb{E}\big{[}\text{(III)}\big{]}\lesssim_{\alpha,d,M,l}\frac{\sigma^{-4ab}\delta^{-4b}}{\sqrt{n}\eta^{3}}+\frac{\sigma^{-6ab}\delta^{-6b}}{n\eta^{2}}.$

Finally, bounding $\log N_{\sigma,\delta}$ and
$\mathbb{E}\big{[}\Delta(X_{1:n})\big{]}$ in (35) using (14) and (37),
respectively, we obtain

$\mathbb{E}[\text{(IV)}]\lesssim_{\alpha,d,M,l}\frac{\sigma^{-a(1+b)}\delta^{-b}n^{-1/4}}{\eta}.$

As earlier, let $\sigma=\sigma_{n}$. Choosing
$\eta=\eta_{n}=KC_{\alpha,d,M,l}n^{-\frac{1-b}{2(3+b)}}\sigma^{-a}$ for $K>1$
and $\delta=\delta_{n}=n^{-1/(6+2b)}\sigma_{n}^{-a}$ gives

$\mathbb{E}\big{[}\text{(I)}+\text{(II)}+\text{(III)}+\text{(IV)}\big{]}=O(1)/K.$

Choosing $K=K_{n}\to\infty$, Lemma 12 yields the claim (i). For the claim
(ii), we set
$\delta_{n}\mspace{-2.0mu}=\mspace{-2.0mu}n^{-1/(6+2b)}\sigma_{n}^{-a}$ and
$\eta_{n}=2C_{\alpha,d,M,l}n^{-(1-b)/(4(3+b))}\sigma^{-a/2}$, to obtain

$\rho\big{(}\mathsf{Law}\left(\|\mathbb{G}_{n}^{*}\|_{\mathcal{F}_{\sigma}}\middle|\,X_{1:n}\right),\mathsf{Law}\left(\|G_{P}\|_{\mathcal{F}_{\sigma}}\right)\big{)}\leq\eta_{n}\wedge\big{[}\text{(I)+(I)+(III)+(IV)}\big{]},$

where
$\mathbb{E}\left[\text{(I)+(I)+(III)+(IV)}\right]\lesssim_{\alpha,d,M,l}n^{-(1-b)/(4(3+b))}\sigma^{-a/2}$.
Hence,

$\rho\big{(}\mathsf{Law}\left(\|\mathbb{G}_{n}^{*}\|_{\mathcal{F}_{\sigma}}\middle|\,X_{1:n}\right),\mathsf{Law}\left(\|G_{P}\|_{\mathcal{F}_{\sigma}}\right)\big{)}=O_{\mathbb{P}}\big{(}n^{-(1-b)/(4(3+b))}\sigma^{-a/2}\big{)},$

completing the proof. ∎

##  Appendix C Supplementary Proofs for Section 5

###  C.1 Supplementary Proofs for Theorem 4

###### Proof of Lemma 13.

Consider the following bound

$\mathbb{E}[X\mathds{1}_{A}]\leq\mathbb{E}\left[X\mathds{1}_{\big{\\{}X\geq
F_{X}^{-1}(1-u)\big{\\}}}\right]\leq\int_{1-u}^{1}F_{X}^{-1}(t)\,\mathrm{d}t\leq\int_{1-u}^{1}\left(\frac{M}{1-t}\right)^{1/p}\,\mathrm{d}t\leq\frac{M}{1-1/p}u^{1-1/p}.$

∎

###### Proof of Lemma 14.

The proof is identical to that of Theorem 2.7.1 in [83], except for the fact
that the $\epsilon^{1/\alpha}$-net of $B(\epsilon^{1/\alpha})$ is just one
point: namely the origin. The details are omitted for brevity. ∎

###### Proof of Lemma 15.

Without loss of generality let $\|F\|_{L^{2}(P)}=1$. Choose integers
$q_{0}=0,q_{2}$ such that $2^{-q_{2}-2}<\gamma<2^{-q_{2}-1}$. Construct a
nested sequence of partitions $\mathcal{F}=\cup_{i=1}^{N_{q}}\mathcal{F}_{qi}$
such that for $q=q_{0}=0$, we have the entire function class as our partition,
and for each integer $q\geq q_{0}$, $\mathcal{F}_{qi}$ satisfies

$\Big{\|}\sup_{f,g\in\mathcal{F}_{qi}}|f-g|^{*}\Big{\|}_{L^{2}(P)}<2^{-q}\|F\|_{L^{2}(P)},$

where $\mathcal{F}_{01}$ consists of the entire function class. $N_{q}$ can be
chosen to satisfy

$\log N_{q}\leq\sum_{r=1}^{q}\log
N_{[\,]}\big{(}2^{-q}\|F\|_{L^{2}(P)},\mathcal{F},\|\cdot\|_{L^{2}(P)}\big{)}.$

Fix a function $f_{qi}$ in each partition $\mathcal{F}_{qi}$, with $f_{01}=0$,
and let

$\pi_{q}f=f_{qi},\ \ \text{for }f\in\mathcal{F}_{qi},$
$\triangle_{q}f=\sup_{f,g\in\mathcal{F}_{qi}}|f-g|^{*},\ \ \text{for
}f\in\mathcal{F}_{qi}.$

Define the following quantities:

$\displaystyle a_{q}$ $\displaystyle=\frac{2^{-q}}{\sqrt{1+\log N_{q+1}}},$
$\displaystyle A_{q-1}f$
$\displaystyle=\mathds{1}\\{\triangle_{q_{0}}f\leq\sqrt{n}a_{q_{0}},\dots,\triangle_{q-1}f\leq\sqrt{n}a_{q-1}\\},$
$\displaystyle B_{q}f$
$\displaystyle=\mathds{1}\\{\triangle_{q_{0}}f\leq\sqrt{n}a_{q_{0}},\dots,\triangle_{q-1}f\leq\sqrt{n}a_{q-1},\triangle_{q-1}f>\sqrt{n}a_{q-1},\\}$
$\displaystyle B_{q_{0}}f$
$\displaystyle=\mathds{1}\\{\triangle_{q_{0}}f>\sqrt{n}a_{q_{0}}\\}.$

The following decomposition decomposition then holds for each
$f\in\mathcal{F}$ (see [83, Page 241]).

$\begin{split}(f-\pi_{q_{0}}f)&B_{q_{0}}f\ +\
\sum_{q_{0}+1}^{q_{2}}(f-\pi_{q}f)B_{q}f+\sum_{q_{0}+1}^{q_{2}}(\pi_{q}f-\pi_{q-1}f)A_{q-1}f+\
(f-\pi_{q_{2}}f)A_{q_{2}}f.\end{split}$ (38)

In our case $\pi_{q_{0}}f=0$ for each $f\in\mathcal{F}$. Applying
$\mathbb{G}_{n}$ to both sides and taking supremum over the function class, we
get:

$\begin{split}\mathbb{E}\left[\|\mathbb{G}_{n}\|_{\mathcal{F}}\right]\leq\mathbb{E}&\left[\sup_{f\in\mathcal{F}}\
\mathbb{G}_{n}fB_{q_{0}}f\right]\ +\
\sum_{q=1}^{q_{2}}\mathbb{E}\left[\sup_{f\in\mathcal{F}}\mathbb{G}_{n}(f-\pi_{q}f)B_{q}f\right]\\\
&\ +\
\sum_{q=1}^{q_{2}}\mathbb{E}\left[\sup_{f\in\mathcal{F}}\mathbb{G}_{n}(\pi_{q}f-\pi_{q-1}f)A_{q-1}f\right]+\
\mathbb{E}\left[\sup_{f\in\mathcal{F}}\mathbb{G}_{n}(f-\pi_{q_{2}}f)A_{q_{2}}f\right].\end{split}$
(39)

The first term in the display above can simply be bounded as follows:

$\mathbb{E}\left[\sup_{f\in\mathcal{F}}\mathbb{G}_{n}fB_{q_{0}}f\right]\leq
2\sqrt{n}PF\mathds{1}_{\\{2F>\sqrt{n}a_{q_{0}}\\}}.$

Treatment of the second and third terms in (39) is identical to that in Lemma
2.14.3 in [83]. For the final term, note that

$\mathbb{G}_{n}(f-\pi_{q_{2}}f)A_{q_{2}}f\leq\sqrt{n}P_{n}\triangle_{q_{2}}fA_{q_{2}}f+\sqrt{n}P\triangle_{q_{2}}fA_{q_{2}}f\leq
2\sqrt{n}a_{q_{2}}\leq 8\gamma\sqrt{n}.$

∎

###  C.2 Proof of Lemma 2

Suppose that $P$ is supported on the hyperplane $H=\\{x_{0}+v:v\in V_{s}\\}$,
where $V_{s}$ is an $s$-dimensional subspace of $\mathbb{R}^{d}$. Suppose
$K_{\epsilon}$ is large enough so that $P(\|X-x_{0}\|\geq
K_{\epsilon})<\epsilon^{s^{2}/(2(s-2))}$. Then, $B(x_{0},K_{\epsilon})\cap H$
can be covered by
$\big{(}1+K_{\epsilon}/\epsilon\big{)}^{s}\lesssim_{s,d}\big{(}\nicefrac{{K_{\epsilon}}}{{\epsilon}}\big{)}^{s}$
many $\epsilon$-balls.

Let $Z=\|X-x_{0}\|$. Then, $\|Z\|_{\psi_{\beta}}=A<\infty$ by the assumptions
of the lemma. Hence, we have $\mathbb{P}(Z\geq Ax)\leq
C_{\beta}e^{-x^{\beta}}$. Consequently, we may choose
$K_{\epsilon}\lesssim_{s,A,\beta}\log(1/\epsilon)^{1/\beta}$ and obtain

$N_{\epsilon}(P,\epsilon^{s^{2}/(2(s-2))})\lesssim_{s,A,\beta}\epsilon^{-s}\log(1/\epsilon)^{1/\beta}.$

This satisfies condition (29), and also gives

$\limsup_{\epsilon\to 0}\frac{\log
N_{\epsilon}(P,\epsilon^{s^{2}/(2(s-2))})}{-\log\epsilon}\leq s,$

which implies $\dim^{*}_{\mathsf{SW}_{1}}(P)\leq s$.

###  C.3 Proof of Corollary 5

Note that
$\widetilde{s}:=\big{(}\dim_{\mathsf{SW1}}^{*}(P)+s\big{)}/2>\dim_{\mathsf{SW1}}^{*}(P)$.
Since $\widetilde{s}$ depends only on $s$ and $P$, by Theorem 4, we have

$\mathbb{E}\left[\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)\right]\lesssim_{d,s,P}\sigma^{-\widetilde{s}/2+1}(\log
n)^{3/2}.$

Now,

$\displaystyle\mathbb{E}\left[\mathsf{W}_{1}(P_{n},P)\right]$
$\displaystyle\leq\inf_{\sigma>0}\Bigg{\\{}\mathbb{E}\left[\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)\right]+2\sigma\sqrt{d}\Bigg{\\}}$
$\displaystyle\lesssim_{d,s,P}\inf_{\sigma>0}\Bigg{\\{}\sigma^{-\widetilde{s}/2+1}(\log
n)^{3/2}/\sqrt{n}+2\sigma\sqrt{d}\Bigg{\\}}$
$\displaystyle\lesssim_{d,s,P}\left(\frac{(\log
n)^{3}}{n}\right)^{1/\widetilde{s}}$ $\displaystyle\lesssim_{d,s,P}n^{-1/s}.$

##  Appendix D Proofs for Section 6

###  D.1 Preliminaries

The following technical lemmas will be needed.

######  Lemma 18 (Continuity of $\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}$).

The smooth Wasserstein distance $\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}$ is
lower semicontinuous (l.s.c.) relative to the weak convergence on
$\mathcal{P}(\mathbb{R}^{d})$ and continuous in $\mathsf{W}_{1}$. Explicitly,
(i) if $\mu_{k}\stackrel{{\scriptstyle w}}{{\to}}\mu$ and
$\nu_{k}\stackrel{{\scriptstyle w}}{{\to}}\nu$, then

$\liminf_{k\to\infty}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu_{k},\nu_{k})\geq\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu,\nu);$

and (ii) if $\mathsf{W}_{1}(\mu_{k},\mu)\to 0$ and
$\mathsf{W}_{1}(\nu_{k},\nu)\to 0$, then

$\lim_{k\to\infty}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu_{k},\nu_{k})=\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu,\nu).$
(40)

###### Proof of Lemma 18.

Part (i). We first note that if $\mu_{k}\stackrel{{\scriptstyle
w}}{{\to}}\mu$, then $\mu_{k}\ast\mathcal{N}_{\sigma}\stackrel{{\scriptstyle
w}}{{\to}}\mu\ast\mathcal{N}_{\sigma}$. This follows from the facts that weak
convergence is equivalent to pointwise convergence of characteristic
functions, and the Gaussian measure has a nonvanishing characteristic function
$\mathbb{E}_{X\sim\mathcal{N}_{\sigma}}[e^{it\cdot
X}]=e^{-\sigma^{2}\|t\|^{2}/2}\neq 0$ for all $t\in\mathbb{R}^{d}$. Now, if
$\mu_{k}\stackrel{{\scriptstyle w}}{{\to}}\mu$ and
$\nu_{k}\stackrel{{\scriptstyle w}}{{\to}}\nu$, then
$\mu_{k}\ast\mathcal{N}_{\sigma}\stackrel{{\scriptstyle
w}}{{\to}}\mu\ast\mathcal{N}_{\sigma}$ and
$\nu_{k}\ast\mathcal{N}_{\sigma}\stackrel{{\scriptstyle
w}}{{\to}}\nu\ast\mathcal{N}_{\sigma}$. From the lower semicontinuity of
$\mathsf{W}_{1}$ relative to the weak convergence (cf. Remark 6.10 in [87]),
we conclude that
$\liminf_{k\to\infty}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu_{k},\nu_{k})=\liminf_{k\to\infty}\mathsf{W}_{1}(\mu_{k}\ast\mathcal{N}_{\sigma},\nu_{k}\ast\mathcal{N}_{\sigma})\geq\mathsf{W}_{1}(\mu\ast\mathcal{N}_{\sigma},\nu\ast\mathcal{N}_{\sigma})=\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu,\nu)$.

Part (ii). Recall that $\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}$ generates the
same topology as $\mathsf{W}_{1}$, i.e.,

$\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu_{k},\mu)\to 0\ \ \iff\ \
\mathsf{W}_{1}(\mu_{k},\mu)\to 0.$

See Theorem 2 in [43]. So if $\mu_{k}\to\mu$ and $\nu_{k}\to\nu$ in
$\mathsf{W}_{1}$, then
$\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu_{k},\mu)=\mathsf{W}_{1}(\mu_{k}\ast\mathcal{N}_{\sigma},\mu\ast\mathcal{N}_{\sigma})\to
0$ and
$\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\nu_{k},\nu)=\mathsf{W}_{1}(\nu_{k}\ast\mathcal{N}_{\sigma},\nu\ast\mathcal{N}_{\sigma})\to
0$. Thus, by Corollary 6.9 in [87], we have
$\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu_{k},\nu_{k})=\mathsf{W}_{1}(\mu_{k}\ast\mathcal{N}_{\sigma},\nu_{k}\ast\mathcal{N}_{\sigma})\to\mathsf{W}_{1}(\mu_{k}\ast\mathcal{N}_{\sigma},\nu_{k}\ast\mathcal{N}_{\sigma})=\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu,\nu)$.
∎

######  Lemma 19 (Weierstrass criterion for the existence of minimizers).

Let $\mathcal{X}$ be a compact metric space, and let
$f:\mathcal{X}\to\mathbb{R}\cup\\{+\infty\\}$ be l.s.c. (i.e.,
$\liminf_{x\to\overline{x}}f(x)\geq f(\overline{x})$ for any
$\overline{x}\in\mathcal{X}$). Then,
$\mathop{\mathrm{argmin}}_{x\in\mathcal{X}}f(x)$ is nonempty.

###### Proof of Lemma 19.

See, e.g., p. 3 of [77]. ∎

###  D.2 Proof of Theorem 5

By Lemma 19, compactness of $\Theta$, and lower semicontinuity of the map
$\theta\mapsto\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(\omega),Q_{\theta})$
(cf. Lemma 18), we see that
$\mathop{\mathrm{argmin}}_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(\omega),Q_{\theta})$
is nonempty.

To prove the existence of a measurable estimator, we will apply Corollary 1 in
[19]. Consider the empirical distribution as a function on
$\mathcal{X}^{\mathbb{N}}$ with $\mathcal{X}=\mathbb{R}^{d}$, i.e.,
$\mathcal{X}^{\mathbb{N}}\ni x=(x_{1},x_{2},\dots)\mapsto
P_{n}(x)=n^{-1}\sum_{i=1}^{n}\delta_{x_{i}}$. Observe that
$\mathcal{X}^{\mathbb{N}}$ and $\mathbb{R}^{d_{0}}$ are both Polish,
$\mathcal{D}:=\mathcal{X}^{\mathbb{N}}\times\Theta$ is a Borel subset of the
product metric space $\mathcal{X}^{\mathbb{N}}\times\mathbb{R}^{d_{0}}$, the
map $\theta\mapsto\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(x),Q_{\theta})$
is l.s.c. by Lemma 18, and the set
$\mathcal{D}_{x}=\big{\\{}\theta\in\Theta:(x,\theta)\in\mathcal{D}\big{\\}}\subset\mathbb{R}^{d_{0}}$
is $\sigma$-compact (as any subset in $\mathbb{R}^{d_{0}}$ is
$\sigma$-compact). Thus, in view of Corollary 1 of [19], it suffices to verify
that the map
$(x,\theta)\mapsto\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(x),Q_{\theta})$
is jointly measurable.

To this end, we use the following fact: for a real function
$\mathcal{Y}\times\mathcal{Z}\ni(y,z)\mapsto f(y,z)\in\mathbb{R}$ defined on
the product of a separable metric space $\mathcal{Y}$ (endowed with the Borel
$\sigma$-field) and a measurable space $\mathcal{Z}$, if $f(y,z)$ is
continuous in $y$ and measurable in $z$, then $f$ is jointly measurable; see
e.g. Lemma 4.51 in [4]. Equip $\mathcal{P}_{1}(\mathbb{R}^{d})$ with the
metric $\mathsf{W}_{1}$ and the associated Borel $\sigma$-field; the metric
space $(\mathcal{P}_{1}(\mathbb{R}^{d}),\mathsf{W}_{1})$ is separable [87,
Theorem 6.16]. Then, since the map $\mathcal{X}^{\mathbb{N}}\ni x\mapsto
P_{n}(x)\in\mathcal{P}_{1}(\mathbb{R}^{d})$ is continuous (which is not
difficult to verify), the map
$\mathcal{X}^{\mathbb{N}}\times\Theta\ni(x,\theta)\mapsto(P_{n}(x),\theta)\in\mathcal{P}_{1}(\mathbb{R}^{d})\times\Theta$
is continuous and thus measurable. Second, by Lemma 18, the function
$\mathcal{P}_{1}(\mathbb{R}^{d})\times\Theta\ni(\mu,\theta)\mapsto\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu,Q_{\theta})\in[0,\infty)$
is continuous in $\mu$ and l.s.c. (and thus measurable) in $\theta$, from
which we see that the map
$(\mu,\theta)\mapsto\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(\mu,Q_{\theta})$ is
jointly measurable. Conclude that the map
$(x,\theta)\mapsto\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(x),Q_{\theta})$
is jointly measurable. ∎

###  D.3 Proof of Theorem 6

The proof relies on Theorem 7.33 in [74], and is reminiscent of that of
Theorem B.1 in [12]; we present a simpler derivation under our
assumption.555Theorem B.1 in [12] applies Theorem 7.31 in [74]. To that end,
one has to extend the maps
$\theta\mapsto\mathcal{W}_{p}(\widehat{\mu}_{n},\mu_{\theta})$ and
$\theta\mapsto\mathcal{W}_{p}(\mu_{\ast},\mu_{\theta})$ to the entire
Euclidean space $\mathbb{R}^{d_{\theta}}$. The extension was not mentioned in
the proof of [12, Theorem B.1], although this missing step does not affect
their final result. To apply Theorem 7.33 in [74], we extend the map
$\theta\mapsto\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})$ to the
entire Euclidean space $\mathbb{R}^{d_{0}}$ as

$\mathbb{G}_{n}(\theta):=\begin{cases}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})&\text{if
\ $\theta\in\Theta$}\\\ +\infty&\text{if \
$\theta\in\mathbb{R}^{d_{0}}\setminus\Theta$}\end{cases}.$

Likewise, define

$g(\theta):=\begin{cases}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta})&\text{if
\ $\theta\in\Theta$}\\\ +\infty&\text{if \
$\theta\in\mathbb{R}^{d_{0}}\setminus\Theta$}\end{cases}.$

The function $\mathbb{G}_{n}$ is stochastic,
$\mathbb{G}_{n}(\theta)=\mathbb{G}_{n}(\theta,\omega)$, but $g$ is non-
stochastic. By construction, we see that
$\mathop{\mathrm{argmin}}_{\theta\in\mathbb{R}^{d_{0}}}g_{n}(\theta)=\mathop{\mathrm{argmin}}_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})$
and
$\mathop{\mathrm{argmin}}_{\theta\in\mathbb{R}^{d_{0}}}g(\theta)=\mathop{\mathrm{argmin}}_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta})$.
In addition, by Lemma 18, continuity of the map $\theta\mapsto Q_{\theta}$
relative to the weak topology, and closedness of the parameter space $\Theta$,
we see that both $\mathbb{G}_{n}$ and $g$ are l.s.c. (on
$\mathbb{R}^{d_{0}}$). The main step of the proof is to show a.s. epi-
convergence of $\mathbb{G}_{n}$ to $g$. Recall the definition of epi-
convergence (in fact, this is an equivalent characterization; see [74,
Proposition 7.29]):

######  Definition 5 (Epi-convergence).

For extended-real-valued functions $f_{n},f$ on $\mathbb{R}^{d_{0}}$ with $f$
being l.s.c., we say that $f_{n}$ epi-converges to $f$ if the following two
conditions hold:

  1. (i)

$\liminf_{n\to\infty}\inf_{\theta\in\mathcal{K}}f_{n}(\theta)\geq\inf_{\theta\in\mathcal{K}}f(\theta)$
for any compact set $\mathcal{K}\subset\mathbb{R}^{d_{0}}$; and

  2. (ii)

$\limsup_{n\to\infty}\inf_{\theta\in\mathcal{U}}f_{n}(\theta)\leq\inf_{\theta\in\mathcal{U}}f(\theta)$
for any open set $\mathcal{U}\subset\mathbb{R}^{d_{0}}$.

We also need the concept of level-boundedness.

######  Definition 6 (Level-boundedness).

For an extended-real-valued function $f$ on $\mathbb{R}^{d_{0}}$, we say that
$f$ is level-bounded if for any $\alpha\in\mathbb{R}$, the set
$\\{\theta\in\mathbb{R}^{d_{0}}:f(\theta)\leq\alpha\\}$ is bounded (possibly
empty).

We are now in position to prove Theorem 6.

###### Proof of Theorem 6.

By boundedness of the parameter space $\Theta$, both $\mathbb{G}_{n}$ and $g$
are level-bounded by construction as the (lower) level sets are included in
$\Theta$. In addition, by assumption, both $\mathbb{G}_{n}$ and $g$ are proper
(an extended-real-valued function $f$ on $\mathbb{R}^{d_{0}}$ is proper if the
set $\\{\theta\in\mathbb{R}^{d_{0}}:f(\theta)<\infty\\}$ is nonempty). In view
of Theorem 7.33 in [74], it remains to prove that $\mathbb{G}_{n}$ epi-
converges to $g$ a.s. To verify property (i) in the definition of epi-
convergence, recall that $P_{n}\to P$ in $\mathsf{W}_{1}$ (and hence in
$\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}$) a.s. Pick any $\omega\in\Omega$ such
that $P_{n}(\omega)\to P$ in $\mathsf{W}_{1}$. Pick any compact set
$\mathcal{K}\subset\mathbb{R}^{d_{0}}$. Since $\mathbb{G}_{n}(\cdot,\omega)$
is l.s.c., by Lemma 19, there exists $\theta_{n}(\omega)\in\mathcal{K}$ such
that
$\mathbb{G}_{n}(\theta_{n}(\omega),\omega)=\inf_{\theta\in\mathcal{K}}\mathbb{G}_{n}(\theta,\omega)$.
Up to extraction of subsequences, we may assume
$\theta_{n}(\omega)\to\theta^{\ast}(\omega)$ for some
$\theta^{\ast}(\omega)\in\mathcal{K}$. If $\theta^{\ast}(\omega)\notin\Theta$,
then by closedness of $\Theta$, $\theta_{n}(\omega)\notin\Theta$ for all
sufficiently large $n$. Thus, we have

$\liminf_{n\to\infty}\inf_{\theta\in\mathcal{K}}\mathbb{G}_{n}(\theta,\omega)=\liminf_{n\to\infty}\mathbb{G}_{n}(\theta_{n}(\omega),\omega)=+\infty,$

so that
$\liminf_{n\to\infty}\inf_{\theta\in\mathcal{K}}\mathbb{G}_{n}(\theta,\omega)\geq\inf_{\theta\in\mathcal{K}}g(\theta)$.
Next, consider the case where $\theta^{\ast}(\omega)\in\Theta$. In this case,
$\theta_{n}(\omega)\in\Theta$ for all $n$ (otherwise,
$+\infty=\mathbb{G}_{n}(\theta_{n}(\omega),\omega)>\mathbb{G}_{n}(\theta^{\ast}(\omega),\omega)$,
which contradicts the construction of $\theta_{n}(\omega)$). Thus,
$\mathbb{G}_{n}(\theta_{n}(\omega),\omega)=\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(\omega),Q_{\theta_{n}(\omega)})$,
so that

$\displaystyle\liminf_{n\to\infty}\inf_{\theta\in\mathcal{K}}g_{n}(\theta_{n}(\omega),\omega)$
$\displaystyle=\liminf_{n\to\infty}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(\omega),Q_{\theta_{n}(\omega)})$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{\geq}}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta^{\ast}(\omega)})$
$\displaystyle\geq\inf_{\theta\in\mathcal{K}}g(\theta),$ (41)

where (a) follows from Lemma 18.

To verify property (ii) in the definition of epi-convergence, pick any open
set $\mathcal{U}\subset\Theta$. It is enough to consider the case where
$\mathcal{U}\cap\Theta\neq\varnothing$. Let
$\\{\theta_{n}^{\prime}\\}_{n=1}^{\infty}\subset\mathcal{U}$ be a sequence
with
$\lim_{n\to\infty}g(\theta_{n}^{\prime})=\inf_{\theta\in\mathcal{U}}g(\theta)$.
Since $\inf_{\theta\in\mathcal{U}}g(\theta)$ is finite, we may assume that
$\theta_{n}^{\prime}\in\mathcal{U}\cap\Theta$ for all $n$. Thus, we have

$\displaystyle\limsup_{n\to\infty}\inf_{\theta\in\mathcal{U}}g_{n}(\theta,\omega)$
$\displaystyle\leq\limsup_{n\to\infty}g_{n}(\theta_{n}^{\prime},\omega)$
$\displaystyle=\limsup_{n\to\infty}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(\omega),Q_{\theta_{n}^{\prime}})$
$\displaystyle\leq\underbrace{\lim_{n\to\infty}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n}(\omega),P)}_{=0}+\underbrace{\lim_{n\to\infty}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta_{n}^{\prime}})}_{=\inf_{\theta\in\mathcal{U}}g(\theta)}$
$\displaystyle=\inf_{\theta\in\mathcal{U}}g(\theta).$ (42)

Conclude that $\mathbb{G}_{n}$ epi-converges to $g$ a.s. This completes the
proof. ∎

###  D.4 Proof of Theorem 7

Recall that $P=Q_{\theta^{\ast}}$. Condition (ii) implies
$\mathop{\mathrm{argmin}}_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta})=\\{\theta^{\ast}\\}$.
Hence, by Theorem 6, for any neighborhood $N$ of $\theta^{\ast}$,

$\inf_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})=\inf_{\theta\in
N}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})$

with probability approaching one.

Define
$R^{\mspace{1.0mu}\sigma}_{\theta}:=Q^{\mspace{1.0mu}\sigma}_{\theta}-P^{\mspace{1.0mu}\sigma}-\left<\theta-\theta^{\ast},D^{\mspace{1.0mu}\sigma}\right>\in\ell^{\infty}(\mathsf{Lip}_{1,0})$,
and choose $N_{1}$ as a neighborhood of $\theta^{\ast}$ such that

$\big{\|}\big{<}\theta-\theta^{\ast},D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}-\big{\|}R^{\mspace{1.0mu}\sigma}_{\theta}\big{\|}_{\mathsf{Lip}_{1,0}}\geq\frac{1}{2}C,\quad\forall\theta\in
N_{1},$ (43)

for some constant $C>0$. Such $N_{1}$ exists since conditions (iii) and (iv)
ensure the existence of an increasing function $\eta(\delta)=o(1)$ (as
$\delta\to 0$) and a constant $C>0$ such that
$\big{\|}R^{\mspace{1.0mu}\sigma}_{\theta}\big{\|}_{\mathsf{Lip}_{1,0}}\leq\|\theta-\theta^{\ast}\|\eta\big{(}\|\theta-\theta^{\ast}\|\big{)}$
and
$\big{\|}\big{<}t,D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}\geq
C\|t\|$ for all $t\in\mathbb{R}^{d_{0}}$.

For any $\theta\in N_{1}$, the triangle inequality and (43) imply that

$\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})\geq\frac{C}{2}\|\theta-\theta^{\ast}\|-\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P).$
(44)

For
$\xi_{n}:=\frac{4\sqrt{n}}{C}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)$,
consider the (random) set
$N_{2}:=\big{\\{}\theta\in\Theta:\sqrt{n}\|\theta-\theta^{\ast}\|\leq\xi_{n}\big{\\}}$.
Note that $\xi_{n}$ is of order $O_{\mathbb{P}}(1)$ by Theorem 1. By the
definition of $\xi_{n}$, $\inf_{\theta\in
N_{1}}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})$ is unchanged if
$N_{1}$ is replaced with $N_{1}\cap N_{2}$; indeed, if $\theta\in N_{2}^{c}$,
then
$\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})>\frac{C}{2}\frac{\xi_{n}}{\sqrt{n}}-\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)=\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)$,
so that $\inf_{\theta\in
N_{2}^{c}}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})>\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)\geq\inf_{\theta\in
N_{1}}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})$.

Reparametrizing $t:=\sqrt{n}(\theta-\theta^{\ast})$ and setting
$T_{n}:=\big{\\{}t\in\mathbb{R}^{d_{0}}:\|t\|\leq\xi_{n},\
\theta^{\ast}+t/\sqrt{n}\in\Theta\big{\\}}$, we have the following
approximation

$\begin{split}\sup_{t\in
T_{n}}\bigg{|}\sqrt{n}\underbrace{\big{\|}P^{\mspace{1.0mu}\sigma}_{n}-Q^{\mspace{1.0mu}\sigma}_{\theta^{\ast}+t/\sqrt{n}}\big{\|}_{\mathsf{Lip}_{1,0}}}_{=\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta^{\ast}+t/\sqrt{n}})}&-\big{\|}\underbrace{\sqrt{n}\big{(}P^{\mspace{1.0mu}\sigma}_{n}-P^{\mspace{1.0mu}\sigma}\big{)}}_{=G^{\mspace{1.0mu}\sigma}_{n}}-\big{<}t,D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}\bigg{|}\\\
&\leq\sup_{t\in
T_{n}}\sqrt{n}\big{\|}R^{\mspace{1.0mu}\sigma}_{\theta^{\ast}+t/\sqrt{n}}\big{\|}_{\mathsf{Lip}_{1,0}}\\\
&\leq\xi_{n}\eta(\xi_{n}/\sqrt{n})\\\ &=o_{\mathbb{P}}(1).\end{split}$ (45)

Observe that any minimizer $t^{\ast}\in\mathbb{R}^{d_{0}}$ of the function
$h_{n}(t):=\big{\|}G^{\mspace{1.0mu}\sigma}_{n}-\big{<}t,D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}$
satisfies $\|t^{\ast}\|\leq\xi_{n}$; indeed if $\|t^{\ast}\|>\xi_{n}$, then
$h_{n}(t^{\ast})\geq
C\|t^{\ast}\|-\|\mathbb{G}_{n}\|_{\mathsf{Lip}_{1,0}}=C\|t^{\ast}\|-\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)=3\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)=3h_{n}(0)$,
which contradicts the assumption that $t^{\ast}$ is a minimizer of $h_{n}(t)$.
Since by assumption $\theta^{\ast}\in\operatorname{int}(\Theta)$, the set of
minimizers of $h_{n}$ lies inside $T_{n}$. Conclude that

$\inf_{\theta\in\Theta}\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})=\inf_{t\in\mathbb{R}^{d_{0}}}\big{\|}G^{\mspace{1.0mu}\sigma}_{n}-\big{<}t,D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}+o_{\mathbb{P}}(1).$
(46)

Now, from the proof of Theorem 1 and the fact that the map
$G\mapsto(G(f\ast\varphi_{\sigma}))_{f\in\mathsf{Lip}_{1,0}}$ is continuous
(indeed, isometric) from $\ell^{\infty}({\mathcal{F}})$ into
$\ell^{\infty}(\mathsf{Lip}_{1,0})$, we see that
$(G^{\mspace{1.0mu}\sigma}_{n}f)_{f\in\mathsf{Lip}_{1,0}}\to
G_{P}^{\mspace{1.0mu}\sigma}$ weakly in $\ell^{\infty}(\mathsf{Lip}_{1,0})$

Applying the continuous mapping theorem to
$L\mapsto\inf_{t\in\mathbb{R}^{d_{0}}}\big{\|}L-\big{<}t,D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}$
and using the approximation (46), we obtain the conclusion of the theorem.∎

###  D.5 Proof of Corollary 6

The proof relies on the following result on weak convergence of argmin
solutions of convex stochastic functions. The following lemma is a simple
modification of Theorem 1 in [55]. Similar techniques can be found in [71] and
[51].

######  Lemma 20.

Let $H_{n}(t)$ and $H(t)$ be convex stochastic functions on
$\mathbb{R}^{d_{0}}$. Let
$\mathop{\mathrm{argmin}}_{t\in\mathbb{R}^{d_{0}}}H(t)$ be unique a.s., and
suppose that for any finite set of points
$t_{1},\dots,t_{k}\in\mathbb{R}^{d_{0}}$, we have

$(H_{n}(t_{1}),\dots,H_{n}(t_{k}))\stackrel{{\scriptstyle
w}}{{\to}}(H(t_{1}),\dots,H(t_{k})).$

Then, for any sequence $\\{\widehat{t}_{n}\\}_{n\in\mathbb{N}}$ such that
$H_{n}(\widehat{t}_{n})\leq\inf_{t\in\mathbb{R}^{d_{0}}}H_{n}(t)+o_{\mathbb{P}}(1)$,
we have

$\widehat{t}_{n}\stackrel{{\scriptstyle
w}}{{\to}}\mathop{\mathrm{argmin}}_{t\in\mathbb{R}^{d_{0}}}H(t).$

###### Proof of Corollary 6.

By Theorem 6, $\widehat{\theta}_{n}\to\theta^{\ast}$ in probability. From
equation (44) and the definition of $\widehat{\theta}_{n}$, we see that, with
probability approaching one,

$\underbrace{\inf_{\theta\in\Theta}\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\theta})}_{=O_{\mathbb{P}}(1)}+o_{\mathbb{P}}(1)\geq\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\widehat{\theta}_{n}})\geq\frac{C}{2}\sqrt{n}\|\widehat{\theta}_{n}-\theta^{\ast}\|-\underbrace{\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)}_{=O_{\mathbb{P}}(1)},$

which implies that
$\sqrt{n}\|\widehat{\theta}_{n}-\theta^{\ast}\|=O_{\mathbb{P}}(1)$. Let
$H_{n}(t):=\big{\|}G^{\mspace{1.0mu}\sigma}_{n}-\big{<}t,D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}$
and
$H(t):=\big{\|}G_{P}^{\mspace{1.0mu}\sigma}-\big{<}t,D^{\mspace{1.0mu}\sigma}\big{>}\big{\|}_{\mathsf{Lip}_{1,0}}$.
Both $H_{n}(t)$ and $H(t)$ are convex in $t$. Then, from equation (45), for
$\widehat{t}_{n}:=\sqrt{n}(\widehat{\theta}_{n}-\theta^{\ast})=O_{\mathbb{P}}(1)$,
we have

$\sqrt{n}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\widehat{\theta}_{n}})=H_{n}(\widehat{t}_{n})+o_{\mathbb{P}}(1).$

Combining the result (46) and the definition of $\widehat{\theta}_{n}$, we see
that
$H_{n}(\widehat{t}_{n})\leq\inf_{t\in\mathbb{R}^{d_{0}}}H_{n}(t)+o_{\mathbb{P}}(1)$.
Since $\mathbb{G}_{n}^{\mspace{1.0mu}\sigma}$ converges weakly to
$G_{P}^{\mspace{1.0mu}\sigma}$ in $\ell^{\infty}(\mathsf{Lip}_{1,0})$, by the
continuous mapping theorem, we have
$(H_{n}(t_{1}),\dots,H_{n}(t_{k}))\stackrel{{\scriptstyle
w}}{{\to}}(H(t_{1}),\dots,H(t_{k}))$ for any finite number of points
$t_{1},\dots,t_{k}\in\mathbb{R}^{d_{0}}$. By assumption,
$\mathop{\mathrm{argmin}}_{t\in\mathbb{R}^{d_{0}}}H(t)$ is unique a.s. Hence,
by Lemma 20, we have $\widehat{t}_{n}\stackrel{{\scriptstyle
w}}{{\to}}\mathop{\mathrm{argmin}}_{t\in\mathbb{R}^{d_{0}}}H(t)$. ∎

######  Remark 12 (Alternative proofs).

Corollary 6 alternatively follows from the proof of Theorem 7 combined with
the argument given at the end of p. 63 in [70] (plus minor modifications), or
the result of Theorem 8 combined with the argument given at the end of p. 67
in [70]. The proof provided above differs from both these arguments and is
more direct.

###  D.6 Proof of Corollary 7

By triangle inequality,

$\displaystyle\mathbb{P}$
$\displaystyle\left(\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\widehat{\theta}_{n}})-\inf_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta})>\epsilon+t\right)$
$\displaystyle\quad\leq\mathbb{P}\left(\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},Q_{\widehat{\theta}_{n}})+\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)-\inf_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta})>\epsilon+t\right)\leq\mathbb{P}\left(2\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)>t\right).$

Now, for a compactly supported distribution, by (3),
$\mathbb{E}\big{[}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)\big{]}\leq
C_{\sigma,d,P}\frac{1}{\sqrt{n}}.$ Then, for all $t>3C_{\sigma,d,P}/\sqrt{n}$,
by (6),

$\mathbb{P}\big{(}2\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)>t\big{)}\leq\mathbb{P}\big{(}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)-\mathbb{E}[\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P_{n},P)]>t/6\big{)}\leq
e^{-\frac{2nt^{2}}{\operatorname{\mathsf{diam}}(\mathcal{X})^{2}}}.$

Choose $C\geq 1$ such that
$Ce^{-18C_{\sigma,d,P}^{2}/\operatorname{\mathsf{diam}}(\mathcal{X})^{2}}\geq
1$. Then, for all $t$, we have

$\mathbb{P}\left(\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\widehat{\theta}_{n}})-\inf_{\theta\in\Theta}\mathsf{W}_{1}^{\mspace{1.0mu}\sigma}(P,Q_{\theta})>\epsilon+t\right)\leq
Ce^{-\frac{2nt^{2}}{\operatorname{\mathsf{diam}}(\mathcal{X})^{2}}}.$

∎