Estimating and Inferring the Maximum Degree of Stimulus-Locked Time-Varying Brain Connectivity Networks

TL;DR

This paper introduces a method to estimate stimulus-locked brain connectivity networks from natural viewing fMRI data, addressing challenges of separating stimulus effects from intrinsic signals, and tests the maximum degree of these networks.

Contribution

It proposes a novel inferential approach to test the maximum degree in stimulus-locked brain networks, controlling type I error and achieving high power asymptotically.

Findings

Method effectively estimates stimulus-locked networks.

Type I error is controlled in the testing procedure.

High power achieved asymptotically.

Abstract

Neuroscientists have enjoyed much success in understanding brain functions by constructing brain connectivity networks using data collected under highly controlled experimental settings. However, these experimental settings bear little resemblance to our real-life experience in day-to-day interactions with the surroundings. To address this issue, neuroscientists have been measuring brain activity under natural viewing experiments in which the subjects are given continuous stimuli, such as watching a movie or listening to a story. The main challenge with this approach is that the measured signal consists of both the stimulus-induced signal, as well as intrinsic-neural and non-neuronal signals. By exploiting the experimental design, we propose to estimate stimulus-locked brain network by treating non-stimulus-induced signals as nuisance parameters. In many neuroscience applications, it is often important to identify brain regions that are connected to many other brain regions during cognitive process. We propose an inferential method to test whether the maximum degree of the estimated network is larger than a pre-specific number. We prove that the type I error can be controlled and that the power increases to one asymptotically. Simulation studies are conducted to assess the performance of our method. Finally, we analyze a functional magnetic resonance imaging dataset obtained under the Sherlock Holmes movie stimuli.

Tables

Table 1. Table 1: The type I error and power for testing the maximum degree of the graph at the 0.05 significance level are calculated as the proportion of falsely rejected and correctly rejected null hypotheses, respectively, over 500 data sets. Simulation results with d = 50 𝑑 50 d=50 and k = { 5 , 6 } 𝑘 5 6 k=\{5,6\} , over a range of n 𝑛 n are shown.

		$n=400$	$n=600$	$n=800$	$n=1000$	$n=1500$
$k$=5	Type I error	0.014	0.024	0.030	0.034	0.028
	Power	0.068	0.182	0.690	0.976	1
$k$=6	Type I error	0.032	0.040	0.034	0.028	0.018
	Power	0.050	0.142	0.446	0.898	1

Equations

$$\mathbf{X}(z)\mid Z=z\sim N_d\{0,{\mathbf{\Sigma }}_X(z)\},$$

$$\mathbf{X}(z)=\mathbf{S}(z)+\mathbf{E}(z),\mathbf{S}(z)\mid Z=z\sim N_d\{0,\mathbf{\Sigma }(z)\},\mathbf{E}(z)\mid Z=z\sim N_d\{0,{\mathbf{L}}_X(z)\},$$

$$\begin{array}{rl}\mathbf{X}(z)=\mathbf{S}(z)+{\mathbf{E}}_X(z),\mathbf{Y}(z)=\mathbf{S}(z)+& {\mathbf{E}}_Y(z),\mathbf{S}(z)|Z=z\sim N_d\{0,\mathbf{\Sigma }(z)\},\\ {\mathbf{E}}_X(z)|Z=z\sim N_d\{0,{\mathbf{L}}_X(z)\},& {\mathbf{E}}_Y(z)|Z=z\sim N_d\{0,{\mathbf{L}}_Y(z)\},\end{array}$$

$$\mathbb{E}[\mathbf{X}(z)\{\mathbf{Y}(z){\}}^T\mid Z=z]=\mathbb{E}[\mathbf{S}(z)\{\mathbf{S}(z){\}}^T\mid Z=z]+\mathbb{E}[{\mathbf{E}}_X(z)\{{\mathbf{E}}_Y(z){\}}^T\mid Z=z]=\mathbf{\Sigma }(z).$$

$$\hat{\mathbf{\Sigma }}(z)=\frac{{\sum }_{i\in [n]}{K}_h(Z_i-z){\mathbf{X}}_i{\mathbf{Y}}_i^T}{{\sum }_{i\in [n]}{K}_h(Z_i-z)},$$

$$K(u)=0.75\cdot \left(1-u^2\right)\cdot {\text{mathds}1}_{\{|u|\le 1\}},$$

$${\hat{\mathbf{\Theta }}}_j(z)=\underset{\mathbf{\theta }\in {\mathbb{R}}^d}{\mathrm{argmin}}\Vert \mathbf{\theta }{\Vert }_1\mathrm{subject}\mathrm{to}{\Vert \hat{\mathbf{\Sigma }}(z)\cdot \mathbf{\theta }-{\mathbf{e}}_j\Vert }_{\infty }\le \lambda ,$$

$${\mathrm{cv}}_{\lambda }=\frac{1}{L}\sum _{\ell =1}^L\sum _{i\in C_{\ell }}\Vert \hat{\mathbf{\Sigma }}(z_i,\lambda {)}^{(\ell )}\hat{\mathbf{\Theta }}{(z_i,\lambda )}^{(-\ell )}-{\mathbf{I}}_d{\Vert }_{\max },$$

$$\begin{array}{rl}& H_0:\mathrm{for}\mathrm{all}z\in [0,1],\mathrm{the}\mathrm{maximum}\mathrm{degree}\mathrm{of}\mathrm{the}\mathrm{graph}\mathrm{is}\mathrm{not}\mathrm{greater}\mathrm{than}k,\\ & H_1:\mathrm{there}\mathrm{exists}\mathrm{a}{z}_0\in [0,1]\mathrm{such}\mathrm{that}\mathrm{the}\mathrm{maximum}\mathrm{degree}\mathrm{of}\mathrm{the}\mathrm{graph}\mathrm{is}\mathrm{greater}\mathrm{than}k.\end{array}$$

$${\hat{\mathbf{\Theta }}}_{jk}^{\mathrm{de}}(z)={\hat{\mathbf{\Theta }}}_{jk}(z)-\frac{{\left\{{\hat{\mathbf{\Theta }}}_j(z)\right\}}^T\left\{\hat{\mathbf{\Sigma }}(z){\hat{\mathbf{\Theta }}}_k(z)-{\mathbf{e}}_k\right\}}{{\left\{{\hat{\mathbf{\Theta }}}_j(z)\right\}}^T{\hat{\mathbf{\Sigma }}}_j(z)},$$

$$T_E=\underset{z\in [0,1]}{\sup }\underset{(j,k)\in E(z)}{\max }\sqrt{nh}\cdot \left|{\hat{\mathbf{\Theta }}}_{jk}^{\mathrm{de}}(z)-{\mathbf{\Theta }}_{jk}(z)\right|\cdot \left\{\frac{1}{n}\sum _{i\in [n]}{K}_h(Z_i-z)\right\}.$$

$$T_E^B=\underset{z\in [0,1]}{\sup }\underset{(j,k)\in E(z)}{\max }\sqrt{nh}\cdot \left|\frac{{\sum }_{i\in [n]}{\left\{{\hat{\mathbf{\Theta }}}_j(z)\right\}}^T{K}_h(Z_i-z)\left\{{\mathbf{X}}_i{\mathbf{Y}}_i^T{\hat{\mathbf{\Theta }}}_k(z)-{\mathbf{e}}_k\right\}{\xi }_i/n}{{\left\{{\hat{\mathbf{\Theta }}}_j(z)\right\}}^T{\hat{\mathbf{\Sigma }}}_j(z)}\right|.$$

$$c(1-\alpha ,E)=\inf \left(t\in \mathbb{R}|P\left[T_E^B\le t|\{(Z_i,{\mathbf{X}}_i,{\mathbf{Y}}_i){\}}_{i\in [n]}\right]\ge 1-\alpha \right).$$

$$c(1-\alpha ,E)=\inf \left[t\in \mathbb{R}\mid P(T_E^B)\le t\mid \{(Z_i,{\mathbf{X}}_i,{\mathbf{Y}}_i){\}}_{i\in [n]}\ge 1-\alpha \right],$$

$$\mathcal{R}(z)=\left\{e\in E(z)\mid \sqrt{nh}\cdot |{\hat{\mathbf{\Theta }}}_e^d(z)|\cdot \sum _{i\in [n]}{K}_h(Z_i-z)/n>c(1-\alpha ,E)\right\}.$$

$$\mathbf{X}(Z_i)=\mathbf{S}(Z_i)+{\mathbf{E}}_X(Z_i)\mathrm{and}\mathbf{Y}(Z_i)=\mathbf{S}(Z_i)+{\mathbf{E}}_Y(Z_i).$$

$$\begin{array}{rl}& H_0:\mathrm{for}\mathrm{all}z\in [z_{\min },z_{\max }],\mathrm{the}\mathrm{maximum}\mathrm{degree}\mathrm{of}\mathrm{the}\mathrm{graph}\mathrm{is}\mathrm{no}\mathrm{greater}\mathrm{than}k,\\ & H_1:\mathrm{there}\mathrm{exists}\mathrm{a}{z}_0\in [z_{\min },z_{\max }]\mathrm{such}\mathrm{that}\mathrm{the}\mathrm{maximum}\mathrm{degree}\mathrm{of}\mathrm{the}\mathrm{graph}\mathrm{is}\mathrm{greater}\mathrm{than}k.\end{array}$$

$$\begin{array}{rl}& H_0:\mathrm{for}\mathrm{all}z\in [0,1],\mathrm{the}\mathrm{maximum}\mathrm{degree}\mathrm{of}\mathrm{the}\mathrm{graph}\mathrm{is}\mathrm{no}\mathrm{greater}\mathrm{than}15,\\ & H_1:\mathrm{there}\mathrm{exists}\mathrm{a}{z}_0\in [0,1]\mathrm{such}\mathrm{that}\mathrm{the}\mathrm{maximum}\mathrm{degree}\mathrm{of}\mathrm{the}\mathrm{graph}\mathrm{is}\mathrm{greater}\mathrm{than}15.\end{array}$$

$$\int K(u)du=1,\int uK(u)du=0,\int u^{l}K(u)du<\infty ,\int K^l(u)du<\infty .$$

$$\underset{z\in [0,1]}{\sup }\underset{j,k\in [d]}{\max }\max \left\{|{\mathbf{\Sigma }}_{jk}(z)|,|{\dot{\mathbf{\Sigma }}}_{jk}(z)|,|{\ddot{\mathbf{\Sigma }}}_{jk}(z)|\right\}\le M_{\sigma }.$$

$$\underset{z\in [0,1]}{\sup }{\Vert \hat{\mathbf{\Sigma }}(z)-\mathbf{\Sigma }(z)\Vert }_{\max }={\mathcal{O}}_P\left\{h^2+\sqrt{\frac{\log (d/h)}{nh}}\right\}.$$

$${\mathcal{U}}_{s,M}=\left\{\mathbf{\Theta }\in {\mathbb{R}}^{d\times d}\mid \mathbf{\Theta }\succ 0,\Vert \mathbf{\Theta }{\Vert }_2\le \rho ,\underset{j\in [d]}{\max }\Vert {\mathbf{\Theta }}_j{\Vert }_0\le s,\underset{j\in [d]}{\max }\Vert {\mathbf{\Theta }}_j{\Vert }_1\le M\right\}.$$

$$\underset{z\in [0,1]}{\sup }{\Vert \hat{\mathbf{\Theta }}(z)-\mathbf{\Theta }(z)\Vert }_{\max }={\mathcal{O}}_P\left\{h^2+\sqrt{\frac{\log (d/h)}{nh}}\right\};$$

$$\underset{z\in [0,1]}{\sup }\underset{j\in [d]}{\max }{\Vert {\hat{\mathbf{\Theta }}}_j(z)-{\mathbf{\Theta }}_j(z)\Vert }_1={\mathcal{O}}_P\left[s\cdot \left\{h^2+\sqrt{\frac{\log (d/h)}{nh}}\right\}\right];$$

$$\underset{z\in [0,1]}{\sup }\underset{j\in [d]}{\max }{\Vert {\left\{{\hat{\mathbf{\Theta }}}_j(z)\right\}}^T\hat{\mathbf{\Sigma }}(z)-{\mathbf{e}}_j\Vert }_{\infty }={\mathcal{O}}_P\left\{h^2+\sqrt{\frac{\log (d/h)}{nh}}\right\}.$$

$$\underset{n\to \infty }{\lim }\underset{\mathbf{\Theta }(\cdot )\in {\mathcal{U}}_{s,M}}{\sup }{P}_{\mathbf{\Theta }(\cdot )}\left\{T_E\ge c(1-\alpha ,E)\right\}\le \alpha .$$

$$\underset{n\to \infty }{\lim }{P}_{\mathrm{null}}(\mathrm{Algorithm}\text{ref}\mathrm{al}:\mathrm{stepdown}\mathrm{rejects}\mathrm{the}\mathrm{null}\mathrm{hypothesis})\le \alpha .$$

$${\text{Sig}}_{\mathrm{deg}}(\mathbf{\Theta }):=\underset{E^{\prime }\subseteq E(\mathbf{\Theta }),\text{Deg}(E)>k}{\max }\underset{e\in E^{\prime }}{\min }|{\mathbf{\Theta }}_e|,$$

$$\mathcal{G}_{1}(\theta)=\Big{[}\bm{\Theta}(\cdot)\in\mathcal{U}_{s,M}\,\Big{|}\,\text{Sig}_{\mathrm{deg}}\{\bm{\Theta}(z_{0})\}\geq\theta\text{ for some $z_{0}\in[0,1]$}\Big{]}.$$

$$\underset{n\to \infty }{\lim }\underset{\mathbf{\Theta }\in {\mathcal{G}}_1(\theta )}{\inf }{\mathbb{P}}_{\mathbf{\Theta }}(\mathrm{Algorithm}\text{ref}\mathrm{al}:\mathrm{stepdown}\mathrm{rejects}\mathrm{the}\mathrm{null}\mathrm{hypothesis})=1,$$

Full text

Keywords: Gaussian multiplier bootstrap; Hypothesis testing; Inter-subject; Latent variables; Maximum degree; Subject specific effects.

1 Introduction

In the past few decades, much effort has been put into understanding task-based brain connectivity networks. For instance, in a typical visual mapping experiment, subjects are presented with a simple static visual stimulus and are asked to maintain fixation at the visual stimulus, while their brain activities are measured. Under such highly controlled experimental settings, numerous studies have shown that there are substantial similarities across brain connectivity networks constructed for different subjects (Press et al. 2001, Hasson et al. 2003). However, such experimental settings bear little resemblance to our real-life experience in several aspects: natural viewing consists of a continuous stream of perceptual stimuli; subjects can freely move their eyes; there are interactions among viewing, context, and emotion (Hasson et al. 2004). To address this issue, neuroscientists have started measuring brain activity under continuous natural stimuli, such as watching a movie or listening to a story (Hasson et al. 2004, Simony et al. 2016, Chen et al. 2017). The main scientific question is to understand the dynamics of the brain connectivity network that are specific to the continuous natural stimuli.

In the neuroscience literature, a typical approach for constructing a brain connectivity network is to calculate a sample covariance matrix for each subject: the covariance matrix encodes marginal relationships for each pair of brain regions within each subject. More recently, graphical models have been used in modeling brain connectivity networks: graphical models encode conditional dependence relationships between each pair of brain regions, given the others (Rubinov & Sporns 2010). A graph consists of $d$ nodes, each representing a random variable, as well as a set of edges joining pairs of nodes corresponding to conditionally dependent variables. There is a vast literature on learning the structure of static undirected graphical models, and we refer the reader to Drton & Maathuis (2017) for a detailed review.

Under natural continuous stimuli, it is often of interest to estimate a dynamic brain connectivity network, i.e., a graph that changes over time. A natural candidate for this purpose is the time-varying Gaussian graphical model (Zhou et al. 2010, Kolar et al. 2010). The time-varying Gaussian graphical model assumes

[TABLE]

where $\bm{\Sigma}_{X}(z)$ is the covariance matrix of $\bm{X}(z)$ given $Z=z$, and $Z\in[0,1]$ has a continuous density. The inverse covariance matrix $\{\bm{\Sigma}_{X}(z)\}^{-1}$ encodes conditional dependence relationships between pairs of random variables at time $Z=z$: $\{\bm{\Sigma}_{X}(z)\}^{-1}_{jk}=0$ if and only if the $j$th and $k$th variables are conditionally independent given the other variables at time $Z=z$.

In natural viewing experiments, the main goal is to construct a brain connectivity network that is locked to the processing of external stimuli, referred to as stimulus-locked network (Simony et al. 2016, Chen et al. 2017, Regev et al. 2018). Constructing a stimulus-locked network can better characterize the dynamic changes of brain patterns across the continuous stimulus (Simony et al. 2016). The main challenge in constructing stimulus-locked network is the lack of highly controlled experiments that remove spontaneous and individual variations. The measured blood-oxygen-level dependent (BOLD) signal consists of not only signal that is specific to the stimulus, but also intrinsic neural signal (random fluctuations) and non-neuronal signal (physiological noise) that are specific to each subject. The intrinsic neural signal and non-neuronal signal can be interpreted as measurement error or latent variables that confound the stimuli-specific signal. We refer to non-stimulus-induced signals as subject specific effects throughout the manuscript. Thus, directly fitting (1) using the measured data will yield a time-varying graph that primarily reflects intrinsic BOLD fluctuations within each brain rather than BOLD fluctuations due to the natural continuous stimulus.

In this paper, we exploit the experimental design aspect of natural viewing experiments and propose to estimate a dynamic stimulus-locked brain connectivity network by treating the intrinsic neural signal and non-neuronal signal as nuisance parameters. Our proposal exploits the fact that the same stimulus will be given to multiple independent subjects, and that the intrinsic neural signal and non-neuronal signal for different subjects are independent. Thus, this motivates us to estimate a brain connectivity network across two brains rather than within each brain. In fact, this approach has been considered in Simony et al. (2016) and Chen et al. (2017) where they estimated brain connectivity networks by calculating covariance for brain regions between two brains.

After estimating the stimulus-locked brain connectivity network, the next important question is to infer whether there are any regions of interest that are connected to many other regions of interest during cognitive process (Hagmann et al. 2008). These highly connected brain regions are referred to as hub nodes, and the number of connections for each brain region is referred to as degree. Identifying hub brain regions that are specific to the given natural continuous stimulus will lead to a better understanding of the cognitive processes in the brain, and may shed light on various cognitive disorders. In the existing literature, several authors have proposed statistical methodologies to estimate networks with hubs (see, for instance, Tan et al. 2014). In this paper, we instead focus on developing a novel inferential framework to test the hypothesis whether there exists at least one time point such that the maximum degree of the graph is greater than $k$.

Our proposed inferential framework is motivated by two major components: (1) the Gaussian multiplier bootstrap for approximating the distribution of supreme of empirical processes (Chernozhukov et al. 2013, 2014b), and (2) the step-down method for multiple hypothesis testing problems (Romano & Wolf 2005). In a concurrent work, Neykov et al. (2019) proposed a framework for testing general graph structure on a static graph. In Appendix A, we will show that our proposed method can be extended to testing a large family of graph structures similar to that of Neykov et al. (2019).

2 Stimulus-Locked Time-Varying Brain Connectivity Networks

2.1 A Statistical Model

Let $\bm{X}(z)$, $\bm{S}(z)$, $\bm{E}(z)$ be the observed data, stimulus-induced signal, and subject specific effects at time $Z=z$, respectively. Assume that $Z$ is a continuous random variable with a continuous density. For a given $Z=z$, we model the observed data as the summation of stimulus-induced signal and the subject specific effects:

[TABLE]

where $\bm{\Sigma}(z)$ is the covariance matrix of the stimulus-induced signal, and $\mathbf{L}_{X}(z)$ is the covariance matrix of the subject specific effects. We assume that $\bm{S}(z)$ and $\bm{E}(z)$ are independent for all $z$. Thus, estimating the stimulus-locked brain connectivity network amounts to estimating $\{\bm{\Sigma}(z)\}^{-1}$. Fitting the model in (1) using the observed data will yield an estimate of $\{\bm{\Sigma}(z)+\mathbf{L}_{X}(z)\}^{-1}$, and thus, (1) fails to estimate the stimulus-locked brain connectivity network $\{\bm{\Sigma}(z)\}^{-1}$.

To address this issue, we exploit the experimental design aspect of natural viewing experiments. In many studies, neuroscientists often measure brain activity for multiple subjects under the same continuous natural stimulus (Chen et al. 2017, Simony et al. 2016). Let $\bm{X}(z)$ and $\bm{Y}(z)$ be measured data for two subjects at time point $Z=z$. Since the same natural stimulus is given to both subjects, this motivates the following statistical model:

[TABLE]

where $\bm{S}(z)$ is the stimulus-induced signal, and $\bm{E}_{X}(z)$ and $\bm{E}_{Y}(z)$ are the subject specific effects at $Z=z$. Model (3) motivates the calculation of inter-subject covariance between two subjects rather than the within-subject covariance. For a given time point $Z=z$, we have

[TABLE]

That is, we estimate $\bm{\Sigma}(z)$ via the inter-subject covariance by treating $\mathbf{L}_{X}(z)$ and $\mathbf{L}_{Y}(z)$ as nuisance parameters. In the neuroscience literature, several authors have been calculating inter-subject covariance matrix to estimate marginal dependencies among brain regions that are stimulus-locked (Chen et al. 2017, Simony et al. 2016). They have found that calculating the inter-subject covariance is able to better capture the stimulus-locked marginal relationships for pairs of brain regions.

For simplicity, throughout the paper, we focus on two subjects. When there are multiple subjects, we can split the subjects into two groups, and obtain an average of each group to estimate the stimulus-locked brain network. We also discuss a $U$-statistic type estimator for the case when there are multiple subjects in Appendix B.

2.2 Inter-Subject Time-Varying Gaussian Graphical Models

We now propose inter-subject time-varying Gaussian graphical models for estimating stimulus-locked time-varying brain networks. Let $(Z_{1},\bm{X}_{1},\bm{Y}_{1}),\ldots,(Z_{n},\bm{X}_{n},\bm{Y}_{n})$ be $n$ independent realizations of the triplets $(Z,\bm{X},\bm{Y})$. Both subjects share the same $Z_{1},\ldots,Z_{n}$ since they are given the same continuous stimulus. Let $K:\mathbb{R}\rightarrow\mathbb{R}$ be a symmetric kernel function. To obtain an estimate for $\bm{\Sigma}(z)$, we propose the inter-subject kernel smoothed covariance estimator

[TABLE]

where $K_{h}(Z_{i}-z)=K\{(Z_{i}-z)/h\}/h$, $h>0$ is the bandwidth parameter, and $[n]=\{1,\ldots,n\}$. For simplicity, we use the Epanechnikov kernel

[TABLE]

where $\operatorname{\mathds{1}}_{\{|u|\leq 1\}}$ is an indicator function that takes value one if $|u|\leq 1$ and zero otherwise. The choice of kernel is not essential as long as it satisfies regularity conditions in Section 5.1.

Let $\bm{\Theta}(z)=\{\bm{\Sigma}(z)\}^{-1}$. Given the kernel smoothed inter-subject covariance estimator in (27), there are multiple approaches to obtain an estimate of the inverse covariance matrix $\bm{\Theta}(z)$. We consider the CLIME estimator proposed by Cai et al. (2011). Let $\mathbf{e}_{j}$ be the $j$th canonical basis in $\mathbb{R}^{d}$. For a vector $\mathbf{v}\in\mathbb{R}^{d}$, let $\|\mathbf{v}\|_{1}=\sum_{j=1}^{d}|v_{j}|$ and let $\|\mathbf{v}\|_{\infty}=\max_{j}|v_{j}|$. For each $j\in[d]$, the CLIME estimator takes the form

[TABLE]

where $\lambda>0$ is a tuning parameter that controls the sparsity of $\widehat{\bm{\Theta}}_{j}(z)$. We construct an estimator for the stimulus-locked brain network as $\widehat{\bm{\Theta}}(z)=[\{\widehat{\bm{\Theta}}_{1}(z)\}^{T},\ldots,\{\widehat{\bm{\Theta}}_{d}(z)\}^{T}]$.

There are two tuning parameters in our proposed method: a bandwidth parameter $h$ that controls the smoothness of the estimated covariance matrix, and a tuning parameter $\lambda$ that controls the sparsity of the estimated network. The bandwidth parameter $h$ can be selected according to the scientific context. For instance, in many neuroscience applications that involve continuous natural stimuli, we select $h$ such that there are always at least 30% of the time points that have non-zero kernel weights. In the following, we propose a $L$-fold cross-validation type procedure to select $\lambda$. We first partition the $n$ time points into $L$ folds. Let $C_{\ell}$ be an index set containing time points for the $\ell$th fold. Let $\bm{\Theta}(z)^{(-\ell)}$ be the estimated inverse covariance matrix using data excluding the $\ell$th fold, and let $\bm{\Sigma}(z)^{(\ell)}$ be the estimated kernel smoothed covariance estimated using data only from the $\ell$th fold. We calculate the following quantity for various values of $\lambda$ :

[TABLE]

where $\|\cdot\|_{\max}$ is the element-wise max norm for matrix. From performing extensive numerical studies, we find that picking $\lambda$ that minimizes the above quantity tend to be too conservative. We instead propose to pick the smallest $\lambda$ with $\mathrm{cv}_{\lambda}$ smaller than the minimum plus two standard deviation.

2.3 Inference on Maximum Degree

We consider testing the hypothesis:

[TABLE]

In the existing literature, many authors have proposed to test whether there is an edge between two nodes in a graph (see, Neykov et al. 2018, and the references therein). Due to the $\ell_{1}$ penalty used to encourage a sparse graph, classical test statistics are no longer asymptotically normal. We employ the de-biased test statistic

[TABLE]

where $\widehat{\bm{\Theta}}_{j}(z)$ is the $j$th column of $\widehat{\bm{\Theta}}(z)$. The subtrahend in (9) is the bias introduced by imposing an $\ell_{1}$ penalty during the estimation procedure.

We use (9) to construct a test statistic for testing the maximum degree of a time-varying graph. Let $G(z)=\{V,E(z)\}$ be an undirected graph, where $V=\{1,\ldots,d\}$ is a set of $d$ nodes and $E(z)\subseteq V\times V$ is a set of edges connecting pairs of nodes. Let

[TABLE]

The edge set $E(z)$ is defined based on the hypothesis testing problem. In the context of testing maximum degree of a time-varying graph as in (8), $E(z)=V\times V$, and therefore the maximum is taken over all possible edges between pairs of nodes. Throughout the manuscript, we will use the notation $E(z)$ to indicate some predefined known edge set. This general edge set will be different for testing different graph structures, and we refer the reader to Appendix A for details.

Since the test statistic (10) involves taking the supreme over $z$ and the maximum over all edges in $E(z)$, it is challenging to evaluate its asymptotic distribution. To this end, we generalize the Gaussian multiplier bootstrap proposed in Chernozhukov et al. (2013) and Chernozhukov et al. (2014b) to approximate the distribution of the test statistic $T_{E}$. Let $\xi_{1},\ldots,\xi_{n}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}N(0,1)$. We construct the bootstrap statistic as

[TABLE]

We denote the conditional $(1-\alpha)$-quantile of $T^{B}_{E}$ given $\{(Z_{i},\bm{X}_{i},\bm{Y}_{i})\}_{i\in[n]}$ as

[TABLE]

The quantity $c(1-\alpha,E)$ can be calculated numerically using Monte-Carlo. In Section 5.2, we show that the quantile of $T_{E}$ in (10) can be estimated accurately by the conditional $(1-\alpha)$-quantile of the bootstrap statistic.

We now propose an inference framework for testing hypothesis problem of the form (8). Our proposed method is motivated by the step-down method in Romano & Wolf (2005) for multiple hypothesis tests. The details are summarized in Algorithm 1. Algorithm 1 involves evaluating all values of $z\in[0,1]$. In practice, we implement the proposed method by discretizing values of $z\in[0,1]$ into a large number of time points. We note that there will be approximation error by taking the maximum over the discretized time points instead of the supremum of the continuous trajectory. The approximation error could be reduced to arbitrarily small if we increase the density of discretization.

In Section 5.2, we will show that Algorithm 1 is able to control the type I error at a pre-specified value $\alpha$. Moreover, the power of the proposed inferential method increases to one as we increase the number of time points $n$. In fact, the proposed inferential method can be generalized to testing a wide variety of structures that satisfy the monotone graph property. Some examples of monotone graph property are that the graph is connected, the graph has no more than $k$ connected components, the maximum degree of the graph is larger than $k$, the graph has no more than $k$ isolated nodes, and the graph contains a clique of size larger than $k$. This generalization will be presented in Appendix A.

3 Simulation Studies

We perform numerical studies to evaluate the performance of our proposal using the inter-subject covariance relative to the typical time-varying Gaussian graphical model using within-subject covariance. To this end, we define the true positive rate as the proportion of correctly identified non-zeros in the true inverse covariance matrix, and the false positive rate as the proportion of zeros that are incorrectly identified to be non-zeros. To evaluate our testing procedure, we calculate the type I error rate and power as the proportion of falsely rejected $H_{0}$ and correctly rejected $H_{0}$, respectively, over a large number of data sets.

To generate the data, we first construct the inverse covariance matrix $\bm{\Theta}(z)$ for $z=\{0,0.2,0.5\}$. At $z=0$, we set $(d-2)/4$ off-diagonal elements of $\bm{\Theta}(0)$ to equal 0.3 randomly with equal probability. At $z=0.2$, we set an additional $(d-2)/4$ off-diagonal elements of $\bm{\Theta}(0)$ to equal 0.3. At $z=0.5$, we randomly select two columns of $\bm{\Theta}(0.2)$ and add $k+1$ edges to each of the two columns. This guarantees that the maximum degree of the graph is greater than $k$. To ensure that the inverse covariance matrix is smooth, for $z\in[0,0.2]$, we construct $\bm{\Theta}(z)$ by taking linear interpolations between the elements of $\bm{\Theta}(0)$ and $\bm{\Theta}(0.2)$. For $z\in[0.2,0.5]$, we construct $\bm{\Theta}(z)$ in a similar fashion based on $\bm{\Theta}(0.2)$ and $\bm{\Theta}(0.5)$. The construction is illustrated in Figure 1.

To ensure that the inverse covariance matrix is positive definite, we set $\bm{\Theta}_{jj}(z)=|\Lambda_{\min}\{\bm{\Theta}(z)\}|+0.1$, where $\Lambda_{\min}\{\bm{\Theta}(z)\}$ is the minimum eigenvalue of $\bm{\Theta}(z)$. We then rescale the matrix such that the diagonal elements of $\bm{\Theta}(z)$ equal one. The covariance $\bm{\Sigma}(z)$ can be obtained by taking the inverse of $\bm{\Theta}(z)$ for each value of $z$. Model (3) involves the subject specific covariance matrix $\mathbf{L}_{X}(z)$ and $\mathbf{L}_{Y}(z)$. For simplicity, we assume that these covariance matrices stay constant over time. We generate $\mathbf{L}_{X}$ by setting the diagonal elements to be one and the off-diagonal elements to be 0.3. Then, we add random perturbations $\epsilon_{k}\epsilon_{k}^{T}$ to $\mathbf{L}_{X}$ for $k=1,\ldots,10$, where $\epsilon_{k}\sim N_{d}(\mathbf{0},\mathbf{I}_{d})$. The matrix $\mathbf{L}_{Y}$ is generated similarly.

To generate the data according to (3), we first generate $Z_{i}\sim\mathrm{Unif}(0,1)$. Given $Z_{1},\ldots,Z_{n}$, we generate $\bm{S}(Z_{i})\mid Z=Z_{i}\sim N_{d}\{\mathbf{0},\bm{\Sigma}(Z_{i})\}$. We then simulate $\bm{E}_{X}(Z_{i})\mid Z=Z_{i}\sim N_{d}(\mathbf{0},\mathbf{L}_{X})$ and $\bm{E}_{Y}(Z_{i})\mid Z=Z_{i}\sim N_{d}(\mathbf{0},\mathbf{L}_{Y})$. Finally, for each value of $Z$, we generate

[TABLE]

Note that both $\bm{X}(Z_{i})$ and $\bm{Y}(Z_{i})$ share the same generated $\bm{S}(Z_{i})$ since both subjects will be given the same natural continuous stimulus. In the following sections, we will assess the performance of our proposal relative to that of typical approach for time-varying Gaussian graphical models using the. within-subject covariance matrix as input. We then evaluate the proposed inferential procedure in Section 2.3 by calculating its type I error and power.

3.1 Estimation

To mimic the data application we consider, we generate the data with $n=945$, $d=172$, and $k=10$. Given the data $(Z_{1},\bm{X}_{1},\bm{Y}_{1}),\ldots,(Z_{n},\bm{X}_{n},\bm{Y}_{n})$, we estimate the covariance matrix at $Z=z$ using the inter-subject kernel smoothed covariance estimator as defined in (27). To obtain estimates of the inverse covariance matrix $\widehat{\bm{\Theta}}(Z_{1}),\ldots,\widehat{\bm{\Theta}}(Z_{n})$, we use the CLIME estimator as described in (6), implemented using the R package clime. There are two tuning parameters $h$ and $\lambda$: we set $h=1.2\cdot n^{-1/5}$ and vary the tuning parameter $\lambda$ to obtain the ROC curve in Figure 2. The smoothing parameter $h$ is selected such that there are always at least 30% of the time points that have non-zero kernel weights. We compare our proposal to time-varying Gaussian graphical models with the kernel smoothed within-subject covariance matrix. The true and false positive rates, averaged over 100 data sets, are in Figure 2.

From Figure 2, we see that our proposed method outperforms the typical approach for time-varying Gaussian graphical models by calculating the within-subject covariance matrix. This is because the typical approach is not estimating the parameter of interest, as discussed in Section 2.2. Our proposed method treats the subject specific effects as nuisance parameters and is able to estimate the stimulus-locked graph accurately.

3.2 Testing the Maximum Degree of a Time-Varying Graph

In this section, we evaluate the proposed inferential method in Algorithm 1 by calculating its type I error and power. In all of our simulation studies, we consider $d=50$ and $B=500$ bootstrap samples, across a range of samples $n$. Similarly, we select the smoothing parameter to be $h=1.2\cdot n^{-1/5}$. The tuning parameter $\lambda$ is then selected using the cross-validation criterion defined in (7). The tuning parameter $\lambda=0.9\cdot[h^{2}+\sqrt{\{\log(d/h)\}/(nh)}]$ is selected for one of the simulated data set. For computational purposes, we use this value of tuning parameter across all replications.

We construct the test statistic $T_{E}$ and the Gaussian multiplier bootstrap statistic $T^{B}_{E}$ as defined in (10) and (11), respectively. Both the statistics $T_{E}$ and $T^{B}_{E}$ involve evaluating the supreme over $z\in[0,1]$. In our simulation studies, we approximate the supreme by taking the maximum of the statistics over 50 evenly spaced grid $z\in[z_{\min},z_{\max}]$, where $z_{\min}=\min\left\{Z_{i}\right\}_{i\in[n]}$ and $z_{\max}=\max\left\{Z_{i}\right\}_{i\in[n]}$.

Our testing procedure tests the hypothesis

[TABLE]

For power analysis, we construct $\bm{\Theta}(z)$ according to Figure 1 by randomly selecting two columns of $\bm{\Theta}(0.2)$ and adding $k+1$ edges to each of the two columns. This ensure that the maximum degree of the graph is greater than $k$. To evaluate the type I error under $H_{0}$, instead of adding $k+1$ edges to the two columns, we instead add sufficient edges such that the maximum degree of the graph is no greater than $k$. For the purpose of illustrating the type I error and power in the finite sample setting, we increase the signal-to-noise ratio of the data by reducing the effect of the nuisance parameters in the data generating mechanism described in Section 3. The type I error and power for $k=\{5,6\}$, averaged over $500$ data sets, are reported in Table 1. We see that the type I error is controlled and that the power increases to one as we increase the number of time points $n$.

4 Sherlock Holmes Data

We analyze a brain imaging data set studied in Chen et al. (2017). This data set consists of fMRI measurements of 17 subjects while watching audio-visual movie stimuli in an fMRI scanner. More specifically, the subjects were asked to watch a 23-minute segment of BBC television series Sherlock, taken from the beginning of the first episode of the series. The fMRI measurements were taken every 1.5 seconds of the movie, yielding $n=945$ brain images for each subject. To understand the dynamics of the brain connectivity network under natural continuous stimuli, we partition the movie into 26 scenes (Chen et al. 2017). The data were pre-processed for slice time correction, motion correction, linear detrending, high-pass filtering, and coregistration to a template brain (Chen et al. 2017). Furthermore, for each subject, we attempt to mitigate issues caused by non-neuronal signal sources by regressing out the average white matter signal.

There are measurements for 271,633 voxels in this data set. For interpretation purposes, we reduce the dimension from 271,633 voxels to $d=172$ regions of interest (ROIs) as described in Baldassano et al. (2015). We map the $n=945$ brain images taken across the 23 minutes into the interval $[0,1]$ chronologically. We then standardize each of the 172 ROIs to have mean zero and standard deviation one.

We first estimate the stimulus-locked time-varying brain connectivity network. To this end, we construct the inter-subject kernel smoothed covariance matrix $\widehat{\bm{\Sigma}}(z)$ as defined in (27). Since there are 17 subjects, we randomly split the 17 subjects into two groups, and use the averaged data to construct (27). Note that we could also construct a brain connectivity network for each pair of subjects separately. We then obtain estimates of the inverse covariance matrix using the CLIME estimator as in (6). We set the smoothing parameter $h=1.2\cdot n^{-1/5}$ so that at least 30% of the kernel weights are non-zero across all time points $Z$. For the sparsity tuning parameter, our theoretical results suggest picking $\lambda=C\cdot\{h^{2}+\sqrt{\log(d/h)/nh}\}$ to guarantee a consistent estimator. We select the constant $C$ by considering a sequence of numbers using a $5$-fold cross-validation procedure described in (7), and this yields $\lambda=1.4\cdot\{h^{2}+\sqrt{\log(d/h)/(nh)}\}$. Heatmaps of the estimated stimulus-locked brain connectivity networks for three different scenes in Sherlock are in Figure 3.

From Figure 3, we see that there are quite a number of connections between brain regions that remain the same across different scenes in the movie. It is also evident that the graph structure changes across different scenes. We see that most brain regions are very sparsely connected, with the exception of a few ROIs. This raises the question of identifying whether there are hub ROIs that are connected to many other ROIs under audio-visual stimuli.

To answer this question, we perform a hypothesis test to test whether there are hub nodes that are connected to many other nodes in the graph across the 26 scenes. If there are such hub nodes, which ROIs do they correspond to. More formally, we test the hypothesis

[TABLE]

The number 15 is chosen since we are interested in testing whether there is any brain region that is connected to more than 10% of the total number of brain regions. We apply Algorithm 1 with 26 values of $z$ corresponding to the middle of the 26 scenes. Figure 4 shows the ROIs that have more than 12 rejected edges across the 26 scenes based on Algorithm 1. Since the maximum degree of the rejected nodes in some scenes are larger than 15, we reject the null hypothesis that the maximum degree of the graph is no greater than 15. In Figure 5, we plot the sagittal snapshot of the brain connectivity network, visualizing the rejected edges from Algorithm 1 and the identified hubs ROIs.

From Figure 4, we see that the rejected hub nodes (nodes that have more than 15 rejected edges) correspond to the frontal pole (7), temporal fusiform cortex (16, 100), lingual gyrus (17), and precuneus (102) regions of the brain. Many studies have suggested that the frontal pole plays significant roles in higher order cognitive operations such as decision making and moral reasoning (among others, Okuda et al. 2003). The fusiform cortex is linked to face and body recognition (see, Iaria et al. 2008, and the references therein). In addition, the lingual gyrus is known for its involvement in processing of visual information about parts of human faces (McCarthy et al. 1999). Thus, it is not surprising that both of these ROIs have more than 15 rejected edges since the brain images are collected while the subjects are exposed to an audio-visual movie stimulus.

Compared to the lingual gyrus, temporal fusiform cortex, and the frontal pole, the precuneus is the least well-understood brain literature in the current literature. We see from Figure 4 that the precuneus is the most connected ROI across many scenes. This is supported by the observation in Hagmann et al. (2008) where the precuneus serves as a hub region that is connected to many other parts of the brain. In recent years, Lerner et al. (2011) and Ames et al. (2015) conducted experiments where subjects were asked to listen to a story under an fMRI scanner. Their results suggest that the precuneus represents high-level concepts in the story, integrating feature information arriving from many different ROIs of the brain. Interestingly, we find that the precuneus has the highest number of rejected edges during the first half of the movie and that the number of rejected edges decreases significantly during the second half of the movie. Our results correspond well to the findings of Lerner et al. (2011) and Ames et al. (2015) in which the precuneus is active when the subjects comprehend the story. However, it also raises an interesting scientific question for future study: is the precuneus active only when the subjects are trying to comprehend the story, that is, once the story is understood, the precuneus is less active.

5 Theoretical Results

We establish uniform rates of convergence for the proposed estimators, and show that the testing procedure in Algorithm 1 is a uniformly valid test. We study the asymptotic regime in which $n$, $d$, and $s$ are allowed to increase. In the context of the Sherlock Holmes data set, $n$ is the total number of brain images obtained under the continuous stimulus, $d$ is the number of brain regions, and $s$ is the maximum number of connections for each brain region in the true stimulus-locked brain connectivity network. The current theoretical results assume that $Z$ is a random variable with continuous density. Our theoretical results can be easily generalized to the case when $\{Z_{i}\}_{i\in[n]}$ are fixed.

5.1 Theoretical Results on Parameter Estimation

Our proposed estimator involves a kernel function $K(\cdot)$: we require $K(\cdot)$ to be symmetric, bounded, unimodal, and compactly supported. More formally, for $l=1,2,3,4$,

[TABLE]

In addition, we require the total variation of $K(\cdot)$ to be bounded, i.e., $\|K\|_{\mathrm{TV}}<\infty$, where $\|K\|_{\mathrm{TV}}=\int|\dot{K}|$. In other words, we require the kernel function to be a smooth function. The Epanechnikov kernel we consider in (5) for analyzing Sherlock data satisfies (13). A unimodal kernel function is extremely plausible in our setting: for instance, to estimate brain network in the “police press conference scene”, we expect the brain images within that scene to play a larger role than brain images that are far away from the scene. One practical limitation of the conditions on the kernel function is the symmetric kernel condition. When we are estimating a stimulus-locked brain network for a particular time point, the ideal case is to weight the previous images more heavily than the future brain images. The scientific reasoning is that there may be some time lag for information processing. In order to capture this effect, a more carefully designed kernel function is needed and it is out of the scope of this paper.

Next, we impose regularity conditions on the marginal density $f_{Z}(\cdot)$.

Assumption 1**.**

There exists a constant $\underline{f}_{Z}$ such that $\inf_{z\in[0,1]}\;f_{Z}(z)\geq\underline{f}_{Z}>0.$ Furthermore, $f_{Z}$ is twice continuously differentiable and that there exists a constant $\bar{f}_{Z}<\infty$ such that $\max\;\{\|f_{Z}\|_{\infty},\|\dot{f}_{Z}\|_{\infty},\|\ddot{f}_{Z}\|_{\infty}\}\leq\bar{f}_{Z}$.

Next, we impose smoothness assumptions on the inter-subject covariance matrix $\bm{\Sigma}(\cdot)$. Our theoretical results hold for any positive definite subject-specific covariance matrices $\mathbf{L}_{X}(z)$ and $\mathbf{L}_{Y}(z)$, since these matrices are treated as nuisance parameters.

Assumption 2**.**

There exists a constant $M_{\sigma}$ such that

[TABLE]

In other words, we assume that the inter-subject covariance matrices are smooth and do not change too rapidly in neighboring time points. This assumption clearly holds in a dynamic brain network where we expect the brain network to change smoothly over time. Assumptions 1 and 2 on $f(z)$ and $\bm{\Sigma}(z)$ are standard assumptions in the nonparametric statistics literature (see, for instance, Chapter 2 of Pagan & Ullah 1999).

The following theorem establishes the uniform rates of convergence for $\widehat{\bm{\Sigma}}(z)$.

Theorem 1**.**

Assume that $h=o(1)$ and that $\log^{2}(d/h)/(nh)=o(1)$. Under Assumptions 1-2, we have

[TABLE]

Theorem 1 guarantees that our estimator always converges to the population parameter under the max norm, if the smoothing parameter $h$ goes to zero asymptotically. For instance, this will satisfy if $h=C\cdot n^{-1/5}$ for some constant $C>0$. The quantity ${\sup}_{{z\in[0,1]}}\|\widehat{\bm{\Sigma}}(z)-\bm{\Sigma}(z)\|_{\max}$ can be upper bounded by the summation of two terms: ${\sup}_{{z\in[0,1]}}\|\mathbb{E}[\widehat{\bm{\Sigma}}(z)]-\bm{\Sigma}(z)\|_{\max}$ and ${\sup}_{{z\in[0,1]}}\|\widehat{\bm{\Sigma}}(z)-\mathbb{E}[\widehat{\bm{\Sigma}}(z)]\|_{\max}$, which are known as the bias and variance terms, respectively, in the kernel smoothing literature (see, for instance, Chapter 2 of Pagan & Ullah 1999). The term $h^{2}$ on the upper bound corresponds to the bias term and the term $\sqrt{\log(d/h)/(nh)}$ corresponds to the variance term.

Next, we establish theoretical results for $\widehat{\bm{\Theta}}(z)$. Recall that the stimulus-locked brain connectivity network is encoded by the support of the inverse covariance matrix $\bm{\Theta}(z)$: $\bm{\Theta}_{jk}(z)=0$ if and only if the $j$th and $k$th brain regions are conditionally independent given all of the other brain regions. We consider the class of inverse covariance matrices:

[TABLE]

Here, $\|\bm{\Theta}\|_{2}$ is the largest singular value of $\bm{\Theta}$ and $\|\bm{\Theta}_{j}\|_{0}$ is the number of non-zeros in $\bm{\Theta}_{j}$.

Brain connectivity networks are usually densely connected due to the intrinsic-neural and non-neuronal signals, such as background processing. Instead of assuming an overall sparse brain network, we assume that the stimulus-locked brain network $\bm{\Theta}(z)$ is sparse, and allow the intrinsic brain network unrelated to the stimulus to be dense. The sparsity assumption on stimulus-locked brain network is plausible in this setting since it characterizes brain activities that are specific to the stimulus. For instance, we may believe that only certain brain regions are active under cognitive process. The other conditions are satisfied since $\bm{\Theta}(z)$ is the inverse of a positive definite covariance matrix.

Given Theorem 1, the following corollary establishes the uniform rates of convergence for $\widehat{\bm{\Theta}}(z)$ using the CLIME estimator as defined in (6). It follows directly from the proof of Theorem 6 in Cai et al. (2011).

Corollary 1**.**

Assume that $\bm{\Theta}(z)\in\mathcal{U}_{s,M}$ for all $z\in[0,1]$. Let $\lambda\geq C\cdot\{h^{2}+\sqrt{\log(d/h)/(nh)}\}$ for $C>0$. Under the same conditions in Theorem 1,

[TABLE]

[TABLE]

[TABLE]

In the real data analysis, Corollary 1 is helpful in terms of selecting the sparsity tuning parameter $\lambda$: it motivates a sparsity tuning parameter of the form $\lambda\geq C\cdot\{h^{2}+\sqrt{\log(d/h)/(nh)}\}$ to guarantee statistically consistent estimated stimulus-locked brain network. To select the constant $C$, we consider a sequence of number and select the appropriate $C$ using a data-driven cross-validation procedure in (7).

5.2 Theoretical Results on Topological Inference

In this section, we first show that the distribution of the test statistic $T_{E}$ can be approximated by the conditional $(1-\alpha)$-quantile of the bootstrap statistic $T^{B}_{E}$. Next, we show that the proposed testing method in Algorithm 1 is valid in the sense that the type I error can be controlled at a pre-specified level $\alpha$.

Recall from (12) the definition of $c(1-\alpha,E)$. The following theorem shows that the Gaussian multiplier bootstrap is valid for approximating the quantile of the test statistic $T_{E}$ in (10). Our results are based on the series of work on Gaussian approximation on multiplier bootstrap in high dimensions (see, e.g., Chernozhukov et al. 2013, 2014b). We see from (10) that $T_{E}$ involves taking the supremum over $z\in[0,1]$ and a dynamic edge set $E(z)$. Due to the dynamic edge set $E(z)$, existing theoretical results for the Gaussian multiplier bootstrap methods cannot be directly applied. We construct a novel Gaussian approximation result for the supreme of empirical processes of $T_{E}$ by carefully characterizing the capacity of the dynamic edge set $E(z)$.

Theorem 2**.**

Assume that $\sqrt{nh^{5}}+s\cdot\sqrt{nh^{9}}=o(1)$. In addition, assume that $s\sqrt{\log^{4}(d/h)/(nh^{2})}+\log^{22}(s)\cdot\log^{8}(d/h)/(nh)=o(1)$. Under the same conditions in Corollary 1, we have

[TABLE]

Some of the scaling conditions are standard conditions in nonparametric estimation (Tsybakov 2009). The most notable scaling conditions are $s\sqrt{\log^{4}(d/h)/(nh^{2})}=o(1)$ and $\log^{22}(s)\cdot\log^{8}(d/h)/(nh)=o(1)$: these conditions arise from Gaussian approximation on multiplier bootstrap (Chernozhukov et al. 2013). These scaling conditions will hold asymptotically as long as the number of brain images $n$ is much larger than the maximum degree in the graph $s$. This corresponds well with the real data analysis where we expect only certain ROIs are active during information processing.

Recall the hypothesis testing problem in (8). We now show that the type I error of the proposed inferential method for testing the maximum degree of a time-varying graph can be controlled at a pre-specified level $\alpha$.

Theorem 3**.**

Assume that the same conditions in Theorem 2 hold. Under the null hypothesis in (8), we have

[TABLE]

To study the power analysis of the proposed method, we define the signal strength of a precision matrix $\bm{\Theta}$ as

[TABLE]

where $\text{Deg}(E)$ is the maximum degree of graph $G=(V,E)$. Under the alternative hypothesis in (8), $\mathrm{there\;exists\;a\;}z_{0}\in[0,1]\mathrm{\;such\;that\;}$ the maximum degree of the graph is greater than $k$. We define the parameter space under the alternative:

[TABLE]

The following theorem presents the power analysis of Algorithm 1.

Theorem 4**.**

Assume that the same conditions in Theorem 2 hold and select the smoothing parameter such that $h=o(n^{-1/5})$. Under the alternative hypothesis in (8) and the assumption that $\theta\geq C\sqrt{\log(d/h)/nh}$, where $C$ is a fixed large constant, we have

[TABLE]

for any fixed $\alpha\in(0,1)$.

The signal strength condition defined in (24) is weaker than the typical minimal signal strength condition required on testing a single edge on a conditional independent graph, $\min_{e\in E(\bm{\Theta})}|\bm{\Theta}_{e}|$. The condition in (24) requires only that there exists a subgraph whose maximum degree is larger than $k$ and the minimal signal strength on that subgraph is above certain level. In our real data analysis, this requires only the edges for brain regions that are highly connected to many other brain regions to be strong, which is plausible since these regions should have high brain activity.

6 Discussion

We consider estimating stimulus-locked brain connectivity networks from data obtained under natural continuous stimuli. Due to lack of highly controlled experiments that remove all spontaneous and individual variations, the measured brain signal consists of not only stimulus-induced signal, but also intrinsic neural signal and non-neuronal signals that are subjects specific. Typical approach for estimating time-varying Gaussian graphical models will fail to estimate the stimulus-locked brain connectivity network accurately due to the presence of subject specific effects. By exploiting the experimental design aspect of the problem, we propose a simple approach to estimating stimulus-locked brain connectivity network. In particular, rather than calculating within-subject smoothed covariance matrix as in the typical approach for modeling time-varying Gaussian graphical models, we propose to construct the inter-subject smoothed covariance matrix instead, treating the subject specific effects as nuisance parameters.

To answer the scientific question on whether there are any brain regions that are connected to many other brain regions during the given stimulus, we propose an inferential method for testing the maximum degree of a stimulus-locked time-varying graph. In our analysis, we found that several interesting brain regions such as the fusiform cortex, lingual gyrus, and precuneus are highly connected. From the neuroscience literature, these brain regions are mainly responsible for high order cognitive operations, face and body recognition, and serve as control region that integrates information from other brain regions. We have also extended the proposed inferential framework to testing various topological graph structures. These are detailed in Appendix A.

The practical limitation of our proposed method is on the Gaussian assumption on the data. While we focus on the time-varying Gaussian graphical model in this paper, our framework can be extended to other types of time-varying graphical models such as the time-varying discrete graphical model or the time-varying nonparanormal graphical model (Kolar et al. 2010, Lu et al. 2015b). Another limitation is the independence assumption on the data across time points. All of our theoretical results can be generalized to the case when the data across time points are correlated, and we leave such generalization for future work.

Acknowledgement

We thank Hasson’s lab at the Princeton Neuroscience Institute for providing us the fMRI data set under audio-visual movie stimulus. We thank Janice Chen for very helpful conversations on preprocessing the fMRI data set and interpreting the results of our analysis.

Appendix A Inference on Topological Structure of Time-Varying Graph

In this section, we generalize Algorithm 1 in the main manuscript to testing various graph structures that satisfy the monotone graph property. In A.1, we briefly introduce some concepts on graph theory. These include the notion of isomorphism, graph property, monotone graph property, and critical edge set. In A.2, we provide a test statistic and an estimate of the quantile of the proposed test statistic using the Gaussian multiplier bootstrap. We then develop an algorithm to test the dynamic topological structure of a time-varying graph which satisfies the monotone graph property.

A.1 Graph Theory

Let $G=(V,E)$ be an undirected graph where $V=\{1,\ldots,d\}$ is a set of nodes and $E\subseteq V\times V$ is a set of edges connecting pairs of nodes. Let $\mathcal{G}$ be the set of all graphs with the same number of nodes. For any two graphs $G=(V,E)$ and $G^{\prime}=(V,E^{\prime})$, we write $G\subseteq G^{\prime}$ if $G$ is a subgraph of $G^{\prime}$, that is, if $E\subseteq E^{\prime}$. We start with introducing some concepts on graph theory (see, for instance, Chapter 4 of Lovász 2012).

Definition 1**.**

Two graphs $G=(V,E)$ and $G^{\prime}=(V,E^{\prime})$ are said to be isomorphic if there exists permutations $\pi:V\rightarrow V$ such that $(j,k)\in E$ if and only if $\{\pi(j),\pi(k)\}\in E^{\prime}$.

The notion of isomorphism is used in the graph theory literature to quantify whether two graphs have the same topological structure, up to any permutation of the vertices (see Chapter 1.2 of Bondy & Murty 1976). We provide two concrete examples on the notion of isomorphism in Figure 6.

Next, we introduce the notion of graph property. A graph property is a property of graphs that depends only on the structure of the graphs, that is, a graph property is invariant under permutation of vertices. A formal definition is given as follows.

Definition 2**.**

For two graphs $G$ and $G^{\prime}$ that are isomorphic, a graph property is a function $\mathcal{P}:\mathcal{G}\rightarrow\{0,1\}$ such that $\mathcal{P}(G)=\mathcal{P}(G^{\prime})$. A graph $G$ satisfies the graph property $\mathcal{P}$ if $\mathcal{P}(G)=1$.

Some examples of graph property are that the graph is connected, the graph has no more than $k$ connected components, the maximum degree of the graph is larger than $k$, the graph has no more than $k$ isolated nodes, the graph contains a clique of size larger than $k$, and the graph contains a triangle. For instance, the two graphs in Figures 6(i) and 6(ii) are isomorphic and satisfy the graph property of being connected.

Definition 3**.**

For two graphs $G\subseteq G^{\prime}$, a graph property $\mathcal{P}$ is monotone if $\mathcal{P}(G)=1$ implies that $\mathcal{P}(G^{\prime})=1$.

In other words, we say that a graph property is monotone if the graph property is preserved under the addition of new edges. Many graph property that are of interest such as those given in the paragraph immediately after Definition 2 are monotone. In Figure 7, we present several examples of graph property that are monotone by showing that adding additional edges to the graph does not change the graph property. For instance, we see from Figure 7(a) that the existing graph with gray edges are connected. Adding the red edges to the existing graph, the graph remains connected and therefore the graph property is monotone. Another example is the graph with maximum degree at least three as in Figure 7(c). We see that adding the red dash edges to the graph preserves the graph property of having maximum degree at least three.

For a given graph $G=(V,E)$, we define the class of edge sets satisfying the graph property $\mathcal{P}$ as

[TABLE]

Finally, we introduce the notion of critical edge set in the following definition.

Definition 4**.**

Given any edge set $E\subseteq V\times V$, we define the critical edge set of $E$ for a given monotone graph property $\mathcal{P}$ as

[TABLE]

For a given monotone graph property $\mathcal{P}$, the critical edge set is the set of edges that will change the graph property of the graph once added to the existing graph. We provide two examples in Figure 8. Suppose that $\mathcal{P}$ is the graph property of being connected. In Figure 8(a), we see that the graph is not connected, and thus $\mathcal{P}(G)=0$. Adding any of the red dash edges in Figure 8(b) changes $\mathcal{P}(G)=0$ to $\mathcal{P}(G)=1$.

A.2 An Algorithm for Topological Inference

Throughout the rest of the paper, we denote $G(z)=\{V,E(z)\}$ as the graph at $Z=z$. We consider hypothesis testing problem of the form

[TABLE]

where $G(\cdot)$ is the true underlying graph and $\mathcal{P}$ is a given monotone graph property as defined in Definition 3. We provide two concrete examples of the hypothesis testing problem in (23).

Example 1**.**

Number of connected components:

[TABLE]

Example 2**.**

Maximum degree of the graph:

[TABLE]

We now propose an algorithm to test the topological structure of a time-varying graph. The proposed algorithm is very general and is able to test the hypothesis problem of the form in (23). Our proposed algorithm is motivated by the step-down algorithm in Romano & Wolf (2005) for testing multiple hypothesis simultaneously. The main crux of our algorithm is as follows. By Definition 4, the critical edge set $\mathcal{C}\{E_{t-1}(z),\mathcal{P}\}$ contains edges that may change the graph property from $\mathcal{P}\{G(z)\}=0$ to $\mathcal{P}\{G(z)\}=1$. Thus, at the $t$-th iteration of the proposed algorithm, it suffices to test whether the edges on the critical edge set $\mathcal{C}\{E_{t-1}(z),\mathcal{P}\}$ are rejected. Let $E_{t}(z)=E_{t-1}(z)\;\cup\;\mathcal{R}(z)$, where $\mathcal{R}(z)$ is the rejected edge set from the critical edge set $\mathcal{C}\{E_{t-1}(z),\mathcal{P}\}$. Since $\mathcal{P}$ is a monotone graph property, if there exists a $z_{0}\in[0,1]$ such that $E_{t}(z_{0})\in\mathscr{P}$, we directly reject the null hypothesis $H_{0}:\mathcal{P}\{G(z)\}=0$ for all $z$. This is due to the definition of monotone graph property that adding more edges does not change the graph property. If $E_{t}(z_{0})\notin\mathscr{P}$, we repeat this process until the null hypothesis is rejected or no more edges in the critical edge set are rejected. We summarize the procedure in Algorithm 2.

Finally, we generalize the theoretical results in Theorems 3 and 4 to the general testing procedure in Algorithm 2. Given a monotone graph property $\mathcal{P}$, let

[TABLE]

We now show that the type I error of the proposed inferential method in Algorithm 2 can be controlled at a pre-specified level $\alpha$.

Theorem 5**.**

Under the same conditions in Theorem 2, we have

[TABLE]

In order to study the power analysis for testing graph structure that satisfies the monotone graph property, we define signal strength of a precision matrix $\bm{\Theta}$ as

[TABLE]

Under $H_{1}:\mathrm{there\;exists\;a\;}z_{0}\in[0,1]\mathrm{\;such\;that\;}\mathcal{P}\{G(z_{0})\}=1$, we define the parameter space

[TABLE]

Again, we emphasize that the signal strength defined in (24) is weaker than the typical minimal signal strength for testing a single edge in a graph $\min_{e\in E(\bm{\Theta})}|\bm{\Theta}_{e}|$. $\text{Sig}(\bm{\Theta})$ only requires that there exists a subgraph satisfying the property of interest such that the minimal signal strength on that subgraph is above certain level. For example, for $\mathcal{P}(G)=1$ if and only if $G$ is connected, it suffices for $\bm{\Theta}$ belongs to $\mathcal{G}_{1}(\theta;\mathcal{P})$ if the minimal signal strength on a spanning tree is larger than $\theta$. The following theorem presents the power analysis of our test.

Theorem 6**.**

Assume that the same conditions in Theorem 2 hold and select the smoothing parameter $h=o(1/n^{-1/5})$. Assume that $\theta\geq C\sqrt{\log(dn)}/n^{2/5}$ for some sufficiently large constant $C$. Under the alternative hypothesis $H_{1}:\mathcal{P}(G)=1$ in (23), we have

[TABLE]

for any fixed $\alpha\in(0,1)$.

Thus, we have shown in Theorem 6 that the power of the proposed inferential method increases to one asymptotically.

Appendix B A $U$-Statistic Type Estimator

The main manuscript primarily concerns the case when there are two subjects. In this section, we present a $U$-statistic type inter-subject covariance to accommodate the case when there are more than two subjects. First, we note that the same natural stimuli is given to all subjects. This motivates the following statistical model for each $Z=z$:

[TABLE]

where $\bm{X}^{(\ell)}$, $\bm{E}^{(\ell)}$, and $\bm{L}^{(\ell)}(z)$ are the data, subject specific effect, and the covariance matrix for the subject specific effect for the $\ell$th subject, respectively. Suppose that there $N$ subjects. Then, the following $U$-statistic type inter-subject covariance matrix can be constructed to estimate $\bm{\Sigma}(z)$:

[TABLE]

We leave the theoretical analysis of the above estimator for future work.

Appendix C Preliminaries

In this section, we define some notation that will be used throughout the Appendix. Let $[n]$ denote the set $\{1,\ldots,n\}$ and let $[d]$ denote the set $\{1,\ldots,d\}$. For two scalars $a,b$, we define $a\vee b=\max(a,b)$. We denote the $\ell_{q}$-norm for the vector $\mathbf{v}$ as $\|\mathbf{v}\|_{q}=(\sum_{j\in[d]}|v_{j}|^{q})^{1/q}$ for $1\leq q<\infty$. In addition, we let $\mathrm{supp}(\mathbf{v})=\{j:v_{j}\neq 0\}$, $\|\mathbf{v}\|_{0}=|\mathrm{supp}(\mathbf{v})|$, and $\|\mathbf{v}\|_{\infty}={\max}_{j\in[d]}\;|v_{j}|$, where $|\mathrm{supp}(\mathbf{v})|$ is the number of non-zero elements in $\mathbf{v}$. For a matrix $\mathbf{A}\in\mathbb{R}^{n_{1}\times n_{2}}$, we denote the $j$th column as $\mathbf{A}_{j}$. We denote the Frobenius norm of $\mathbf{A}$ by $\|\mathbf{A}\|_{F}^{2}=\sum_{i\in[n_{1}]}\sum_{j\in[n_{2}]}A_{ij}^{2}$, the max norm $\|\mathbf{A}\|_{\max}=\max_{i\in[n_{1}],j\in[n_{2}]}|A_{ij}|$, and the operator norm $\|\mathbf{A}\|_{2}=\sup_{\|\mathbf{v}\|_{2}=1}\;\|\mathbf{A}\mathbf{v}\|_{2}$. Given a function $f$, let $\dot{f}$ and $\ddot{f}$ be the first and second-order derivatives, respectively. For $1\leq p<\infty$, let $\|f\|_{p}=(\int f^{p})^{1/p}$ denote the $L_{p}$ norm of $f$ and let $\|f\|_{\infty}=\sup_{x}|f(x)|$. The total variation of $f$ is defined as $\|f\|_{\mathrm{TV}}=\int|\dot{f}|$. We use the Landau symbol $a_{n}=\mathcal{O}(b_{n})$ to indicate the existence of a constant $C>0$ such that $a_{n}\leq C\cdot b_{n}$ for two sequences $a_{n}$ and $b_{n}$. We write $a_{n}=o(b_{n})$ if ${\lim}_{n\rightarrow\infty}\;a_{n}/b_{n}\rightarrow 0$. Let $C,C_{1},C_{2},\ldots$ be generic constants whose values may vary from line to line.

Let

[TABLE]

For notational convenience, for fixed $j,k\in[d]$, let

[TABLE]

[TABLE]

and let

[TABLE]

Recall from 5 that $K(\cdot)$ can be any symmetric kernel function that satisfies (13) and that $K_{h}(Z_{i}-z)=K\{(Z_{i}-z)/h\}/h$. By the definition of $\widehat{\bm{\Sigma}}(z)$ in (27), we have

[TABLE]

In addition, let

[TABLE]

[TABLE]

[TABLE]

and let

[TABLE]

For two functions $f$ and $g$, we define its convolution as

[TABLE]

In our proofs, we will use the following property of the derivative of a convolution

[TABLE]

Finally, our proofs use the following inequality

[TABLE]

where the first inequality holds by an application of Jensen’s inequality.

Appendix D Proof of Results in 5.1

In this section, we establish the uniform rate of convergence for $\widehat{\bm{\Sigma}}(z)$ and $\widehat{\bm{\Theta}}(z)$ over $z\in[0,1]$. To prove Theorem 1, we first observe that

[TABLE]

The first term is known as the variance term and the second term is known as the bias term in the kernel smoothing literature (see, for instance, Chapter 2 of Pagan & Ullah 1999). Both the variance and bias terms involve evaluating the quantity $\mathbb{E}\{\widehat{\bm{\Sigma}}_{jk}(z)\}$. From (32), we see that $\widehat{\bm{\Sigma}}_{jk}(z)$ involves the quotient of two averages and it is not straightforward to evaluate its expectation. The following lemma quantifies $\mathbb{E}\{\widehat{\bm{\Sigma}}_{jk}(z)\}$ in terms of the expectations of its numerator and its denominator.

Lemma 1**.**

Under the following conditions

[TABLE]

we have

[TABLE]

We note that (42) only holds under the two conditions in (41). In the proof of Theorem 1, we will show that the two conditions in (41) hold for $n$ sufficiently large. To obtain upper bounds for the bias and variance terms in (40), we use the following intermediate lemmas.

Lemma 2**.**

Assume that $h=o(1)$. Under Assumptions 1-2, we have

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Lemma 3**.**

Assume that $h=o(1)$ and $\log^{2}(d/h)/(nh)=o(1)$. Under Assumptions 1-2, there exists a universal constant $C>0$ such that

[TABLE]

with probability at least $1-3/d$.

The proofs of Lemmas 1-3 are deferred to Sections 1-3, respectively. We now provide a proof of Theorem 1.

D.1 Proof of Theorem 1

Recall from (40) that

[TABLE]

It suffices to obtain upper bounds for $I_{1}$ and $I_{2}$.

We first verify that the two conditions in (41) hold. By Lemma 2, we have

[TABLE]

where the last inequality follows from Assumption 1. Moreover,

[TABLE]

for sufficiently large $n$, where the first inequality is obtained by an application of Lemma 2, the second inequality is obtained by an application of Lemma 3, and the last inequality is obtained by the scaling assumptions $h=o(1)$ and $\log(d/h)/(nh)=o(1)$.

Upper bound for $I_{1}$: By (55) in the proof of Lemma 1, we have

[TABLE]

Thus, by Lemma 1, we have

[TABLE]

where $I_{13}=\mathcal{O}[\{\mathbb{G}_{n}(w_{z})\mathbb{G}_{n}(g_{j,zk})\}+\mathbb{G}_{n}^{2}(g_{z,jk})+\mathbb{E}\{\mathbb{G}_{n}(w_{z})\cdot\mathbb{G}_{n}(g_{j,zk})\}+\mathbb{E}\{\mathbb{G}_{n}^{2}(g_{z,jk})\}]/n$.

We now provide upper bounds for $I_{11}$, $I_{12}$, and $I_{13}$. By an application of Lemmas 2 and 3, we obtain

[TABLE]

Similarly, we have

[TABLE]

For $I_{13}$, we have

[TABLE]

where the first and second inequalities follow from Lemmas 2 and 3, respectively. Combining (49), (50), and (51), we have

[TABLE]

with probability at least $1-3/d$.

Upper bound for $I_{2}$: By Lemmas 1 and 2, we have

[TABLE]

where the first inequality follows from (43) and (44), the second equality follows from (45) and (46), and the last inequality follows from the assumption that $h=o(1)$.

Combining the upper bounds (52) and (53), we obtain

[TABLE]

with probability at least $1-3/d$.

Appendix E Proof of Technical Lemmas in Appendix D

In this section, we provide the proofs of Lemmas 1-3.

E.1 Proof of Lemma 1

The proof of the lemma uses the following fact

[TABLE]

From (32), we have

[TABLE]

Under the conditions (41) and by applying (54), we have

[TABLE]

Note that $\mathbb{E}\{\mathbb{G}_{n}(f)\}=0$ by the definition of $\mathbb{G}_{n}(f)$ in (28). Taking expectation on both sides of (55), we obtain

[TABLE]

as desired.

E.2 Proof of Lemma 2

To prove Lemma 2, we write the expectation as an integral and apply Taylor expansion to the density function and the covariance function. We will show that the higher-order terms of the Taylor expansion can be bounded by $\mathcal{O}(h^{2})$. We start by proving (43).

Proof of (43): Recall from (29) the definition of $g_{z,jk}(Z_{i},X_{ij},Y_{ik})=K_{h}(Z_{i}-z)X_{ij}Y_{ik}$. Thus, we have

[TABLE]

where the third equality hold using the fact that the subject-specific effects are independent between two subjects, and the last equality holds by a change of variable, $u=(Z-z)/h$. Applying Taylor expansions to $\bm{\Sigma}_{jk}(uh+z)$ and $f_{Z}(uh+z)$, we have

[TABLE]

and

[TABLE]

where $z^{\prime}$ and $z^{\prime\prime}$ are between $z$ and $uh+z$. Substituting (57) and (58) into the last expression of (56), we have

[TABLE]

By (13), we have $\int uK(u)du=0$ and $\int u^{l}K(u)du<\infty$ for $l=1,2,3,4$. By Assumptions 1 and 2, we have

[TABLE]

Substituting (60) into (59) and bounding the other higher-order terms by $\mathcal{O}(h^{2})$, we obtain

[TABLE]

for all $z\in[0,1]$ and $j,k\in[d]$. This implies that

[TABLE]

The proof of (44) follows from the same set of argument.

Proof of (45): Recall from (29) the definition of $w_{z}(Z_{i})=K_{h}(Z_{i}-z)$. Thus, we have

[TABLE]

where the second to the last equality follows from the fact that $Z_{i}$ and $Z_{i^{\prime}}$ are independent.

We now obtain an upper bound for $I_{1}$. By (13) and Assumptions 1-2, we have

[TABLE]

where the last equality holds by a change of variable. Moreover, by (43) and (44), we have

[TABLE]

Substituting (62) and (63) into (61), and taking the supreme over $z\in[0,1]$ and $j,k\in[d]$ on both sides of the equation, we obtain

[TABLE]

where the last equality holds by the scaling assumption of $h=o(1)$. The proof of (46) follows from the same set of argument.

E.3 Proof of Lemma 3

The proof of Lemma 3 involves obtaining upper bounds for the supreme of the empirical processes $\mathbb{G}_{n}(w_{z})$ and $\mathbb{G}_{n}(g_{z,jk})$. To this end, we apply the Talagrand’s inequality in Lemma 20. Let $\mathcal{F}$ be a function class. In order to apply Talagrand’s inequality, we need to evaluate the quantities $\eta$ and $\tau^{2}$ such that

[TABLE]

Talagrand’s inequality in Lemma 20 provides an upper bound for the supreme of an empirical process in terms of its expectation. By Lemma 21, the expectation can then be upper bounded as a function of the covering number of the function class $\mathcal{F}$, denoted as $N\{\mathcal{F},L_{2}(Q),\epsilon\}$. The following lemmas provide upper bounds for the supreme of the empirical processes $\mathbb{G}_{n}(w_{z})$ and $\mathbb{G}_{n}(g_{z,jk})$, respectively. The proofs are deferred to Sections E.3.1 and E.3.2, respectively.

Lemma 4**.**

Assume that $h=o(1)$ and $\log(d/h)/(nh)=o(1)$. Under Assumptions 1-2, for sufficiently large $n$, there exists a universal constant $C>0$ such that

[TABLE]

with probability at least $1-1/d$.

Lemma 5**.**

Assume that $h=o(1)$ and $\log^{2}(d/h)/(nh)=o(1)$. Under Assumptions 1-2, for sufficiently large $n$, there exists a universal constant $C>0$ such that

[TABLE]

with probability at least $1-2/d$.

Applying Lemmas 4 and 5, we obtain

[TABLE]

with probability at least $1-3/d$, as desired.

E.3.1 Proof of Lemma 4

The proof of Lemma 4 uses the set of arguments as detailed in the beginning of E.3. Recall from (29) and (31) the definition of $w_{z}(Z_{i})=K_{h}(Z_{i}-z)$ and $k_{z}(Z_{i})=w_{z}(Z_{i})-\mathbb{E}\{w_{z}(Z)\}$, respectively. We consider the class of function

[TABLE]

First, note that

[TABLE]

where the first inequality holds by (13) and Lemma 2, and the last inequality holds by the scaling assumption $h=o(1)$ for sufficiently large $n$.

Next, we obtain an upper bound for the variance of $k_{z}(Z_{i})$. Note that

[TABLE]

where we apply the inequality $(x-y)^{2}\leq 2x^{2}+2y^{2}$ for two scalars $x,y$. By Lemma 2, we have $I_{2}\leq 2\{\bar{f}_{Z}+\mathcal{O}(h^{2})\}^{2}$. Also, by a change of variable and second-order Taylor expansion on the marginal density $f_{Z}(\cdot)$, we have

[TABLE]

Thus, for sufficiently large $n$ and the assumption that $h=o(1)$, we have

[TABLE]

By Lemma 16, the covering number for the function class $\mathcal{K}$ satisfies

[TABLE]

We are now ready to obtain an upper bound for the supreme of the empirical process, $\sup_{z\in[0,1]}\left|\mathbb{G}_{n}(w_{z})\right|$. By Lemma 21 with $A=2\cdot\|K\|_{\mathrm{TV}}\cdot C_{K}^{4/5}\cdot\bar{f}_{Z}^{1/5}/\|K\|_{\infty}$, $\|F\|_{L_{2}(\mathbb{P}_{n})}=2\cdot\|K\|_{\infty}/h$, $V=5$, $\sigma^{2}_{P}=3\cdot\bar{f}_{Z}\cdot\|K\|_{2}^{2}/h$, for sufficiently large $n$, we obtain

[TABLE]

where $C>0$ is some sufficiently large constant. By Lemma 20 with $\tau^{2}=3\bar{f}_{Z}\cdot\|K\|_{2}^{2}/h$, $\eta=2\cdot\|K\|_{\infty}/h$, $\mathbb{E}[Y]\leq C\cdot\sqrt{\log(1/h)/(nh)}$, and picking $t=\sqrt{\log(d)/n}$, for sufficiently large $n$, we have

[TABLE]

with probability $1-1/d$, where the last expression holds by the assumption that $\log(d/h)/(nh)=o(1)$ and $h=o(1)$. Multiplying both sides of the above equation by $\sqrt{n}$ completes the proof of Lemma 4.

E.3.2 Proof of Lemma 5

The proof of Lemma 5 uses the set of arguments as detailed in the beginning of E.3. For convenience, we prove Lemma 5 by conditioning on the event

[TABLE]

Since $X_{ij}$ and $Y_{ij}$ conditioned on $Z$ are Gaussian random variables, the event $\mathcal{A}$ occurs with probability at least $1-1/d$ for sufficiently large constant $M_{X}>0$.

Recall from (29) and (30) the definition of $g_{z,jk}(Z_{i},X_{ij},Y_{ik})=K_{h}(Z_{i}-z)X_{ij}Y_{ik}$ and $q_{z,jk}(Z_{i},X_{ij},Y_{ik})=g_{z,jk}(Z_{i},X_{ij},Y_{ik})-\mathbb{E}\{g_{z,jk}(Z,X_{j},Y_{k})\}$, respectively. We consider the function class

[TABLE]

We first obtain an upper bound for the function class

[TABLE]

where the second inequality holds by Assumptions 1-2 and Lemma 2, the third inequality holds by (13) and by conditioning on the event $\mathcal{A}$, and the last inequality holds by the scaling assumption $h=o(1)$ for sufficiently large $n$.

Next, we obtain an upper bound for the variance of $q_{z,jk}(Z_{i},X_{ij},Y_{ik})$. Note that

[TABLE]

where we apply the inequality $(x-y)^{2}\leq 2x^{2}+2y^{2}$ for two scalars $x,y$. By Lemma 2, we have $I_{2}\leq 2\left\{\bar{f}_{Z}\cdot M_{\sigma}+\mathcal{O}(h^{2})\right\}^{2}$. Also, by a change of variable and second-order Taylor expansion on the marginal density $f_{Z}(\cdot)$ as in (68), we have

[TABLE]

where the first inequality follows from the fact that $|\mathbb{E}(X_{j}^{2}Y_{k}^{2}\mid Z)|\leq\kappa$ for some $\kappa<\infty$ since these are Gaussian random variables, and the second inequality follows from (68). Thus, for sufficiently large $n$ and the assumption that $h=o(1)$, we have

[TABLE]

By Lemma 17, the covering number for the function class $\mathcal{Q}$ satisfies

[TABLE]

We now obtain an upper bound for the supreme of the empirical process, $\underset{z\in[0,1]}{\sup}\;\underset{j,k\in[d]}{\max}\;\left|\mathbb{G}_{n}(g_{z,jk})\right|$. By Lemma 21 with $A=2\cdot\|K\|_{\mathrm{TV}}\cdot C_{K}^{4/5}\cdot\bar{f}_{Z}^{1/5}\cdot M_{\sigma}^{1/5}\cdot d^{1/10}/\|K\|_{\infty}$, $\|F\|_{L_{2}(\mathbb{P}_{n})}=2\cdot\|K\|_{\infty}\cdot M^{2}_{X}\cdot\log d/h$, $V=5$, $\sigma^{2}_{P}=(3\kappa/h)\cdot\bar{f}_{Z}\cdot\|K\|_{2}^{2}$, for sufficiently large $n$, we obtain

[TABLE]

where the last inequality holds by the assumption $\log(d/h)/nh=o(1)$. By Lemma 20 with $\tau^{2}=3\cdot\kappa\cdot\bar{f}_{Z}\cdot\|K\|_{2}^{2}/h$, $\eta=2\cdot\|K\|_{\infty}\cdot M^{2}_{X}\cdot\log d/h$, $\mathbb{E}[Y]\leq C\cdot\sqrt{\log(d/h)/(nh)}$, and picking $t=\sqrt{\log d/n}$, for sufficiently large $n$, we have

[TABLE]

with probability at least $1-2/d$. The second inequality holds by the assumption that $\log^{2}(d/h)/(nh)=o(1)$. Multiplying both sides of the equation by $\sqrt{n}$, we completed the proof of Lemma 5.

Appendix F Proof of Theorem 2

In this section, we provide the proof of Theorem 2. To prove Theorem 2, we use a similar set of arguments in the series of work on Gaussian multiplier bootstrap of the supreme of empirical process (see, for instance, Chernozhukov et al. 2013, 2014a, 2014b). Recall from (10) and (11) that

[TABLE]

and

[TABLE]

respectively, where $\xi_{i}\sim N(0,1)$. Note that for notational convenience, we drop the subscript $E$ from $T_{E}$ and $T_{E}^{B}$ throughout the proof.

We aim to show that $T^{B}$ is a good approximation of $T$. However, $T$ and $T^{B}$ are not exact averages. To apply the results in Chernozhukov et al. 2014a, we define four intermediate processes:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\xi_{i}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}N(0,1)$.

To prove Theorem 2, we show that $T_{00}$ is a good approximation of $T$ and that $T_{00}^{B}$ is a good approximation of $T^{B}$. We then show that there exists a Gaussian process $W$ such that both $T_{00}^{B}$ and $T_{00}$ can be accurately approximated by $W$. This is done by applications of Theorems A.1 and A.2 in Chernozhukov et al. (2014a). The following summarizes the chain of empirical and Gaussian processes that we are going to study

[TABLE]

The following lemma provides an approximation error between the statistic $T$ and the intermediate empirical process $T_{00}$.

Lemma 6**.**

Assume that $h^{2}+\sqrt{\log(d/h)/nh}=o(1)$. Under Assumptions 1-2, for sufficiently large $n$, there exists a universal constant $C>0$ such that

[TABLE]

with probability at least $1-1/d$.

Proof.

The proof is deferred to F.2. ∎

We now apply Theorems A.1 and A.2 in Chernozhukov et al. (2014a) to show that there exists a Gaussian process $W$ such that the quantities $|T_{00}-W|$ and $|T_{00}^{B}-W|$ can be controlled, respectively. The results are stated in the following lemmas.

Lemma 7**.**

Assume that $\log^{6}s\cdot\log^{4}(d/h)/(nh)=o(1)$. Under Assumptions 1-2, for sufficiently large $n$, there exists universal constants $C,C^{\prime}>0$ such that

[TABLE]

Proof.

The proof is deferred to F.3. ∎

Lemma 8**.**

Assume that $\log^{4}(s)\cdot\log^{3}(d/h)/(nh)=o(1)$. Under Assumptions 1-2, for sufficiently large $n$, there exists universal constants $C,C^{\prime\prime}>0$ such that

[TABLE]

with probability at least $1-3/n$.

Proof.

The proof is deferred to F.4. ∎

Finally, the following lemma provides an upper bound on the difference between $T^{B}$ and $T_{00}^{B}$, conditioned on the data $\left\{(Z_{i},\bm{X}_{i},\bm{Y}_{i})\right\}_{i\in[n]}$.

Lemma 9**.**

Assume that $s\cdot\sqrt{h^{3}\log^{3}(d/h)}+s\cdot\sqrt{\log^{4}(d/h)/nh^{2}}+\sqrt{h^{5}\log n}=o(1)$. Under Assumptions 1-2, for sufficiently large $n$, there exists universal constants $C,C^{\prime\prime}>0$ such that, with probability at least $1-1/d$,

[TABLE]

Proof.

The proof is deferred to F.5. ∎

With Lemmas 6-9, we are now ready to prove Theorem 2.

F.1 Proof of Theorem 2

Recall that for notational convenience, we drop the subscript $E$ from $T_{E}$ and $T_{E}^{B}$ throughout the proof. In this section, we show that $T$ can be well-approximated by the $(1-\alpha)$-conditional quantile of $T^{B}$, i.e., $P\{T\geq c(1-\alpha)\}\leq\alpha$. For notational convenience, we let $r=r_{1}+r_{2}+r_{3}+r_{4}$, where

[TABLE]

These are the scaling that appears in Lemmas 6-9. By Lemmas 6 and 7, it can be shown that

[TABLE]

since $r_{2}\geq r_{1}$ and $r_{2}\geq 1/d$. With some abuse of notation, throughout the proof, we write $P_{\xi}(T^{B}\geq t)$ to indicate $P[T^{B}\geq t\mid\{(Z_{i},\bm{X}_{i},\bm{Y}_{i})\}_{i\in[n]}]$. By Lemmas 8 and 9, we have

[TABLE]

since $r_{2}\geq r_{3}$ and $r_{2}\geq 2/d+1/n$. Define the event

[TABLE]

and note that $P(\mathcal{E})\geq 1-2/d-4/n$ by Lemmas 8 and 9. Throughout the proof, we condition on the event $\mathcal{E}$.

By the triangle inequality, we obtain

[TABLE]

where the last inequality follows from (84). By a similar argument and by (85), we have

[TABLE]

where the last inequality follows from the fact that $P(X\leq t-\epsilon)-P(X\leq t)\geq-P(|X-t|\leq\epsilon)$ for any $\epsilon>0$. Thus, combining (86) and (87), we obtain

[TABLE]

It remains to show that the quantity $P_{\xi}\{|T^{B}-c(1-\alpha)|\leq r_{2}\}$ converges to zero as we increase $n$.

By the definition of $T_{00}$ and from (35), we have

[TABLE]

Let $\widehat{\sigma}^{2}_{z,jk}=\sum_{i=1}^{n}J^{2}_{z,jk}(Z_{i},\bm{X}_{i},\bm{Y}_{i})/n$ be the conditional variance, and let $\underline{\sigma}=\inf_{z,jk}\widehat{\sigma}_{z,jk}$ and $\bar{\sigma}=\sup_{z,jk}\widehat{\sigma}_{z,jk}$. By Lemma A.1 of Chernozhukov et al. (2014b) and Theorem 3 of Chernozhukov et al. (2013), we obtain

[TABLE]

We first calculate the quantity $\bar{\sigma}$. By (110), we have

[TABLE]

Moreover, by (110), we have

[TABLE]

Define the function class $\mathcal{J}^{\prime}=\{J^{2}_{z,jk}(\cdot)\mid z\in[0,1],j,k\in[d]\}$. By Lemmas 15, 18 and 19, we have

[TABLE]

Thus, applying Lemma 21 with $\sigma^{2}_{P}=C\cdot\log^{4}s/h^{2}$ and $\|F\|_{L_{2}(\mathbb{P}_{n})}\leq C\cdot d^{2}\cdot(d^{17/24}\cdot\log^{3/4}d/h^{11/12})^{24}$, we have

[TABLE]

By an application of the Markov’s inequality, we obtain

[TABLE]

Thus, we have with probability at least $1-C\cdot\left\{\log^{5}(d/h)/(nh^{2})\right\}^{1/4}$,

[TABLE]

where the last inequality follows from (115) for sufficiently large $n$. By Lemma 10, we have $\inf_{z,j,k}\mathbb{E}\{J^{2}_{z,jk}(Z,\bm{X},\bm{Y})\}\geq c>0.$ Therefore, we have

[TABLE]

with probability at least $1-C\cdot\left\{\log^{5}(d/h)/(nh^{2})\right\}^{1/4}$.

Next, we calculate the quantity $\mathbb{E}[T^{B}_{00}\mid\{(Z_{i},\bm{X}_{i},\bm{Y}_{i})\}_{i\in[n]}]$. By Dudley’s inequality (see, e.g., Corollary 2.2.8 in Van Der Vaart & Wellner 1996) and (116), we obtain

[TABLE]

Moreover, by Lemma 9, we have

[TABLE]

with probability at least $1-2/d-1/n$. Substituting (94), (95), and (96) into (89), we obtain

[TABLE]

Thus, substituting (97) into (88), we have

[TABLE]

By the scaling assumptions, $r_{2}=o(1)$ and $\log^{22}s\cdot\log^{8}(d/h)/(nh)=o(1)$. Thus, this implies that

[TABLE]

which implies that

[TABLE]

as desired.

F.2 Proof of Lemma 6

In this section, we show that $|T-T_{00}|$ is upper bounded by the quantity

[TABLE]

with high probability for sufficiently large constant $C>0$. By the triangle inequality, we have $|T-T_{00}|\leq|T-T_{0}|+|T_{0}-T_{00}|$. Thus, is suffices to obtain upper bounds for the terms $|T-T_{0}|$ and $|T_{0}-T_{00}|$.

Upper Bound for $|T-T_{0}|$: Let $\widetilde{\bm{\Theta}}_{k}=\left(\widehat{\bm{\Theta}}_{1k},\ldots,\widehat{\bm{\Theta}}_{(j-1)k},\bm{\Theta}_{jk},\widehat{\bm{\Theta}}_{(j+1)k},\ldots,\widehat{\bm{\Theta}}_{dk}\right)^{T}\in\mathbb{R}^{d}$. Then, the statistics $T$ can be rewritten as

[TABLE]

To obtain an upper bound on the difference between $T$ and $T_{0}$, we make use of the following inequality:

[TABLE]

Recall from (80) that

[TABLE]

Applying (99) with $x=\{\widehat{\bm{\Theta}}_{j}(z)\}^{T}\{\widehat{\bm{\Sigma}}(z)\widetilde{\bm{\Theta}}_{k}-\mathbf{e}_{k}\}$, $\delta=\{\widehat{\bm{\Theta}}_{j}(z)\}^{T}\widehat{\bm{\Sigma}}_{j}(z)-1$, and $y=\{\bm{\Theta}_{j}(z)\}^{T}\{\widehat{\bm{\Sigma}}(z)-\bm{\Sigma}(z)\}\bm{\Theta}_{k}(z)$, and by the triangle inequality, we have

[TABLE]

It remains to obtain upper bounds for $I_{1}$ and $I_{2}$ in (100).

Upper bound for $I_{1}$: By Corollary 1, we have

[TABLE]

Moreover, by Lemmas 4 and 2, we have

[TABLE]

with probability at least $1-1/d$. Thus, by Holder’s inequality, we have

[TABLE]

with probability greater than $1-4/d$, where the third inequality holds by Theorem 1, (101), and (102).

Upper bound for $I_{2}$: To obtain an upper bound for $I_{2}$, we first decompose the quantity $\sqrt{nh}\cdot\{\widehat{\bm{\Theta}}_{j}(z)\}^{T}\{\widehat{\bm{\Sigma}}(z)\widetilde{\bm{\Theta}}_{k}-\mathbf{e}_{k}\}$ into the following

[TABLE]

Next, we show that $I_{21}$ converges to zero and that the difference between $I_{22}$ and the term $\sqrt{nh}\cdot\{\bm{\Theta}_{j}(z)\}^{T}\{\widehat{\bm{\Sigma}}(z)-\bm{\Sigma}(z)\}\bm{\Theta}_{k}(z)$ is small.

Upper bound for $I_{21}$: By Holder’s inequality and Corollary 1, we have

[TABLE]

with probability at least $1-1/d$.

Decomposition of $I_{22}$: By adding and subtracting terms, we have

[TABLE]

Similar to (104), we have

[TABLE]

where the second inequality holds by Holder’s inequality, Corollary 1, and the fact that $\bm{\Theta}(z)\in\mathcal{U}_{s,M}$.

Combining the results (104)-(106), we have

[TABLE]

where the third inequality follows from (102).

Combining the upper bounds for $I_{1}$ in (103) and $I_{2}$ in (107), we have

[TABLE]

with probability at least $1-1/d$.

Upper bound for $|T_{0}-T_{00}|$: Recall from (81) the definition of $T_{00}$

[TABLE]

Using the triangle inequality $||x|-|y||\leq|x-y|$, we obtain

[TABLE]

where the second inequality follows from an application of Holder’s inequality, the third inequality follows from the fact that $\bm{\Theta}(z)\in\mathcal{U}_{s,M}$, the first equality follows by an application of Lemma 2, and the last inequality follows from Assumption 2 and that $h^{2}=o(1)$.

Thus, combining (108) and (109), there exists a constant $C>0$ such that

[TABLE]

with probability at least $1-1/d$.

F.3 Proof of Lemma 7

Recall from (81) the definition

[TABLE]

Recall from (35) that $J_{z,jk}(Z_{i},\bm{X}_{i},\bm{Y}_{i})=J_{z,jk}^{(1)}(Z_{i},\bm{X}_{i},\bm{Y}_{i})-J_{z,jk}^{(2)}(Z_{i})$, where $J_{z,jk}^{(1)}(Z_{i},\bm{X}_{i},\bm{Y}_{i})$ and $J_{z,jk}^{(2)}(Z_{i})$ are as defined in (33) and (34), respectively. Let $\mathcal{J}=\left\{J_{z,jk}\mid z\in[0,1],\;j,k\in[d]\right\}$. Then the intermediate empirical average $T_{00}$ can be written as

[TABLE]

In this section, we show that there exists a Gaussian process $W$ such that

[TABLE]

with high probability. To this end, we apply Theorem A.1 in Chernozhukov et al. (2014a), which involves the following quantities

• •

upper bound for $\underset{z\in[0,1]}{\sup}\;\underset{j,k\in[d]}{\max}\;\|J_{z,jk}(Z_{i},\bm{X}_{i},\bm{Y}_{i})\|_{\infty}$;

• •

upper bound for $\underset{z\in[0,1]}{\sup}\;\underset{j,k\in[d]}{\max}\;\mathbb{E}\left\{J_{z,jk}^{2}(Z,\bm{X},\bm{Y})\right\}$;

• •

covering number for the function class $\mathcal{J}$.

Let ${\mathcal{S}}_{j}(z)$ and ${\mathcal{S}}_{k}(z)$ to be the support of $\bm{\Theta}_{j}(z)$ and $\bm{\Theta}_{k}(z)$, respectively. Note that the cardinality for both sets are less than $s$. We now obtain the above quantities.

Upper bound for $\underset{z\in[0,1]}{\sup}\;\underset{j,k\in[d]}{\max}\;\|J_{z,jk}(Z_{i},\bm{X}_{i},\bm{Y}_{i})\|_{\infty}$: We have with probability at least $1-1/(2s)$,

[TABLE]

where the first inequality follows by Holder’s inequality and the definition of $q_{z,jk}$ and $k_{z}$ and the second inequality follows from (67) and (74). Note that since we are only taking max over the set ${\mathcal{S}}_{j}(z)$ and ${\mathcal{S}}_{k}(z)$, instead of a $\log d$ factor from (74), we obtain a $\log(2s)$ factor.

Upper bound for $\underset{z\in[0,1]}{\sup}\;\underset{j,k\in[d]}{\max}\;\mathbb{E}\{J_{z,jk}^{2}(Z,\bm{X},\bm{Y})\}$: By an application of the inequality $(x-y)^{2}\leq 2x^{2}+2y^{2}$, we have

[TABLE]

To obtain an upper bound for $I_{1}$, we need an upper bound for $\underset{z\in[0,1]}{\sup}\;\underset{j,k\in[d]}{\max}\;\mathbb{E}\{\underset{j\in{\mathcal{S}}_{j}(z),k\in{\mathcal{S}}_{k}(z)}{\max}\;q_{z,jk}^{2}\}$. Recall from (29) the definition of $g_{z,jk}(Z_{i},X_{ij},Y_{ik})=K_{h}(Z_{i}-z)X_{ij}Y_{ik}$ and that $q_{z,jk}(Z_{i},X_{ij},Y_{ik})=g_{z,jk}(Z_{i},X_{ij},Y_{ik})-\mathbb{E}\{g_{z,jk}(Z,X_{j},Y_{k})\}$. Thus, we have

[TABLE]

where we apply the fact that $(x-y)^{2}\leq 2x^{2}+2y^{2}$ to obtain the last inequality. By Lemma 2, we have $2\underset{z\in[0,1]}{\sup}\;\underset{j,k\in[d]}{\max}\;\mathbb{E}^{2}(g_{z,jk})\leq 2\left\{\bar{f}_{Z}\cdot M_{\sigma}+\mathcal{O}(h^{2})\right\}^{2}$. Moreover, we have

[TABLE]

with probability at least $1-1/(2s)$, where the second inequality follows from an application of Lemma 2.

Thus, by Holder’s inequality, we have

[TABLE]

where the second inequality holds by the fact that $\bm{\Theta}(z)\in\mathcal{U}_{s,M}$.

Similarly, to obtain an upper bound for $I_{2}$, we use the fact from (69) that

[TABLE]

By Holder’s inequality, we have

[TABLE]

where the second inequality holds by Assumption 2 and by the fact that $\bm{\Theta}(z)\in\mathcal{U}_{s,M}$, and the last inequality holds by (113).

Combining the upper bounds for $I_{1}$ (112) and $I_{2}$ (114), we have

[TABLE]

for sufficiently large $C>0$.

Covering number of the function class $\mathcal{J}$: First, we note that the function class $\mathcal{J}$ is generated from the addition of two function classes

[TABLE]

Thus, to obtain the covering number of $\mathcal{J}$, we first obtain the covering number for the function classes $\mathcal{J}^{(1)}_{jk}$ and $\mathcal{J}^{(2)}_{jk}$. Then, we apply Lemma 15 to obtain the covering number of the function class $\mathcal{J}$. From Lemma 18, we have with probability at least $1-1/d$,

[TABLE]

Moreover, from Lemma 19, we have

[TABLE]

Applying Lemma 15 with $a_{1}=d^{5/4}\cdot\log^{3/2}d/h^{1/2}$, $v_{1}=6$, $a_{2}=d^{1/6}/h^{4/3}$, and $v_{2}=6$, we have

[TABLE]

where we multiply $d^{2}$ on the right hand side since the function class $\mathcal{J}$ is taken over all $j,k\in[d]$.

Application of Theorem A.1 in Chernozhukov et al. (2014a): Applying Theorem A.1 in Chernozhukov et al. (2014a) with $a=d^{65/24}\cdot\log^{7/4}d/h^{17/12}$, $b=C\cdot\log s/\sqrt{h}$, $\sigma_{J}=C\cdot\log s$, and

[TABLE]

for sufficiently large constant $A,C>0$, there exists a random process $W$ such that for any $\gamma\in(0,1)$,

[TABLE]

for some absolute constant $C^{\prime}$. Picking $\gamma=\left\{\log^{6}s\cdot\log^{4}(d/h)/(nh)\right\}^{1/8}$, we have

[TABLE]

as desired.

F.4 Proof of Lemma 8

Recall from the proof of Lemma 7 that

[TABLE]

We note that

[TABLE]

where $\xi_{i}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}N(0,1)$. To show that the term $|W-T_{00}^{B}|$ can be controlled, we apply Theorem A.2 in Chernozhukov et al. (2014a).

Let

[TABLE]

as defined in Theorem A.2 in Chernozhukov et al. (2014a). From the proof of Lemma 7, we have $b=C\cdot\log s/\sqrt{h}$, $K_{n}=C\cdot\log(d/h)$, and $\sigma_{J}=C\cdot\log s$. Since $b^{2}K_{n}=C\cdot\log^{2}s\cdot\log(d/h)/h\leq n\cdot\log^{2}s$ for sufficiently large $n$, there exists a constant $C^{\prime\prime}>0$ such that

[TABLE]

with probability at least $1-3/n$. Choosing $\delta=\left\{\log^{4}(s)\cdot\log^{3}(d/h)/(nh)\right\}^{1/8}$, we have

[TABLE]

with probability at least $1-3/n$.

F.5 Proof of Lemma 9

In this section, we show that $|T^{B}-T^{B}_{00}|$ is upper bounded by the quantity

[TABLE]

with high probability for sufficiently large constant $C>0$. Throughout the proof of this lemma, we conditioned on the data $\{(Z_{i},\bm{X}_{i},\bm{Y}_{i})\}_{i\in[n]}$. By the triangle inequality, we have $|T^{B}-T_{00}^{B}|\leq|T^{B}-T_{0}^{B}|+|T_{0}^{B}-T_{00}^{B}|$. Thus, it suffices to obtain upper bounds for the terms $|T^{B}-T_{0}^{B}|$ and $|T_{0}^{B}-T_{00}^{B}|$.

Upper bound for $|T^{B}-T_{0}^{B}|$: Recall from (79) and (82) that

[TABLE]

and that

[TABLE]

respectively. Using the triangle inequality, we have

[TABLE]

where the second inequality holds by another application of the triangle inequality and inequality in (99). We now obtain upper bounds for $I_{1}$, $I_{2}$, and $I_{3}$.

Upper bound for $I_{1}$: By an application of Holder’s inequality, we have

[TABLE]

where the last inequality follows from the fact that $\bm{\Theta}(z)\in\mathcal{U}_{s,M}$ and by an application of Corollary 1. For notational convenience, we use the notation as defined in (36)

[TABLE]

Then, we have

[TABLE]

We note that conditioned on the data $\{(Z_{i},\bm{X}_{i},\bm{Y}_{i})\}_{i\in[n]}$, the above expression is a Gaussian process. It remains to bound the supreme of the Gaussian process

[TABLE]

in probability.

To this end, we apply the Dudley’s inequality (see, e.g., Corollary 2.2.8 in Van Der Vaart & Wellner 1996) and the Borell’s inequality (see, e.g., Proposition A.2.1 in Van Der Vaart & Wellner 1996), which involves the following quantities:

• •

upper bound on the conditional variance $\sum_{i\in[n]}W_{z,jk}^{2}(Z_{i},X_{ij},Y_{ik})/n;$

• •

the covering number of the function class

[TABLE]

under the $L_{2}$ norm on the empirical measure.

**Upper bound for the conditional variance $\sum_{i=1}^{n}W_{z,jk}^{2}(Z_{i},X_{ij},Y_{ik})/n:$ ** By the definition of $W_{z,jk}(Z_{i},X_{ij},Y_{ik})$ in (119), we have

[TABLE]

with probability at least $1-1/d$. Note that the second inequality holds by the fact that $(x-y)^{2}\leq 2x^{2}+2y^{2}$, and the third inequality holds by (13) and Assumption 2, and the fact that $\max(X_{ij},Y_{ij})\leq M_{X}\cdot\sqrt{\log d}$ with probability at least $1-1/d$.

Covering number of the function class $\mathcal{W}$: To obtain the covering number of the function class $\mathcal{W}$ under the $L_{2}$ norm on the empirical measure, it suffices to obtain the covering number $\underset{Q}{\sup}\;N\{\mathcal{W},L_{2}(Q),\epsilon\}$. First, we note that $W_{z,jk}=\sqrt{h}\cdot\left\{g_{z,jk}-w_{z}\cdot\bm{\Sigma}_{jk}(z)\right\}$. From Lemma 16, we have $\mathcal{K}_{1}=\{w_{z}(\cdot)\mid z\in[0,1]\}$ and that

[TABLE]

Also, From Lemma 17, we have $\mathcal{G}_{1,jk}=\{g_{z,jk}(\cdot)\mid z\in[0,1]\}$ and that

[TABLE]

Moreover, by Assumption 2, $\bm{\Sigma}_{jk}(z)$ is $M_{\sigma}$-Lipschitz. Thus, applying Lemmas 14 and 15, we obtain

[TABLE]

where the term $d^{2}$ appear on the right hand side because the function class $\mathcal{W}$ is over $j,k\in[d]$.

Applying Dudley’s inequality and Borell’s inequality: Applying Dudley’s inequality (see Corollary 2.2.8 in Van Der Vaart & Wellner 1996) with (120) and (121), we have

[TABLE]

Applying (39) with $b_{1}=C\cdot\sqrt{\log^{2}d/h}$ and $b_{2}=d^{2/9}\cdot\log^{4/9}d/h^{17/18}$, we have

[TABLE]

for some sufficiently large $C>0$.

By Borell’s inequality (see Proposition A.2.1 in Van Der Vaart & Wellner 1996), for $\lambda>0$, we have

[TABLE]

where $\sigma^{2}_{X}$ is the upper bound on the conditional variance. Picking $\lambda=C\cdot\sqrt{\frac{\log^{3}(d/h)}{h}}$, we have

[TABLE]

Thus, substituting (123) into (118), we have

[TABLE]

with probability $1-1/d$.

Upper bound for $I_{2}$: By an application of Holder’s inequality, we have

[TABLE]

where the second inequality holds by triangle inequality and Corollary 1, and the last inequality holds by another application of Corollary 1 and the assumption that $h^{2}+\sqrt{\log(d/h)/(nh)}=o(1)$.

Recall the definition of $g_{z,jk}(Z_{i},X_{ij},Y_{ik})=K_{h}(Z_{i}-z)X_{ij}Y_{ik}$. Conditioned on the data $\{(Z_{i},\bm{X}_{i},\bm{Y}_{i})\}_{i\in[n]}$, we note that

[TABLE]

Similar to the upper bound for $I_{1}$, we apply Dudley’s inequality and Borell’s inequality to bound the supreme of the Gaussian process in the last expression.

To this end, we need to obtain an upper bound for the conditional covariance. By (74), we have

[TABLE]

with probability at least $1-1/d$. In addition, by an application of Lemma 17, the covering number for the class of function $\{\sqrt{h}\cdot g_{z,jk}(\cdot)\;|\;z\in[0,1],\;j,k\in[d]\}$ is

[TABLE]

By an application of Dudley’s inequality, we have

[TABLE]

Applying (39) with $b_{1}=\sqrt{M_{X}^{4}\cdot\|K\|_{\infty}^{4}\cdot\log^{2}d/h}$ and $b_{2}=d^{1/2}\cdot\log d/h^{1/2}$, we have

[TABLE]

By Borell’s inequality (see Proposition A.2.1 in Van Der Vaart & Wellner 1996), we have

[TABLE]

Picking $\lambda=C\cdot\sqrt{\frac{\log^{3}(d/h)}{h}}$, we have

[TABLE]

Thus, by (125) and (129), we have

[TABLE]

with probability at least $1-1/d$.

Upper bound for $I_{3}$: By an application of Holder’s inequality, we have

[TABLE]

where the second inequality holds by the fact that $\bm{\Theta}(z)\in\mathcal{U}_{s,M}$ and by an application of Corollary 1, and the third inequality holds by (123).

Thus, combining (124), (130), and (131), we have

[TABLE]

with probability at least $1-3/d$.

Upper bound for $|T_{0}^{B}-T_{00}^{B}|$: Recall from (83) that

[TABLE]

By the triangle inequality, we have

[TABLE]

where the last inequality holds by applying Holder’s inequality and Lemma 2. Since $\xi_{i}\stackrel{{\scriptstyle\mathrm{i.i.d.}}}{{\sim}}N(0,1)$, by the Gaussian tail inequality, we have

[TABLE]

Thus, substituting the above expression into (133), we obtain

[TABLE]

with probability at least $1-1/n$.

Combining the upper bounds: Combining the upper bounds (132) and (134), and applying the union bound, we have

[TABLE]

as desired.

F.6 Lower Bound of the Variance

We aim to show that the variance of $J_{z,jk}$ defined in (35) is bounded from below.

Lemma 10**.**

Under the same conditions of Theorem 2, there exists a constant $c>0$ such that $\inf_{z}\min_{j,k}\operatorname{{\rm Var}}(J_{z,jk})\geq c>0$.

Proof.

In this proof, we will apply Isserlis’ theorem (Isserlis 1918). Given $\mathbf{T}\sim N(0,\bm{\Sigma})$, Isserlis’ theorem implies that for any vectors $\mathbf{u},\mathbf{v}\in\mathbb{R}^{d}$,

[TABLE]

According to the definition of $J_{z,jk}$ in (35), it can be decomposed into $J_{z,jk}(Z_{i},\bm{X}_{i},\bm{Y}_{i})=J_{z,jk}^{(1)}(Z,\bm{X}_{i},\bm{Y}_{i})-J_{z,jk}^{(2)}(Z_{i})$. Recall that

[TABLE]

and

[TABLE]

We will calculate $\operatorname{{\rm Var}}\{J_{z,jk}^{(2)}(Z)\}$, $\operatorname{{\rm Var}}\{J_{z,jk}^{(1)}(Z,\bm{X},\bm{Y})\}$, and $\operatorname{\rm Cov}\{J_{z,jk}^{(1)}(Z,\bm{X},\bm{Y}),J_{z,jk}^{(2)}(Z)\}$ separately.

We first calculate $\operatorname{{\rm Var}}\{J_{z,jk}^{(2)}(Z)\}$. Following a similar method as the proof of Lemma 2, we have $\mathbb{E}\{K_{h}(Z-z)\}=f_{Z}(z)+\mathcal{O}(h^{2})$ and $\mathbb{E}\{K^{2}_{h}(Z-z)\}=h^{-1}f_{Z}(z)\int K^{2}(u)du+\mathcal{O}(1).$ This implies that

[TABLE]

Next, we proceed to calculate the variance of $J_{z,jk}^{(1)}(Z)$. By a change of variable and Taylor’s expansion, we obtain

[TABLE]

Note that each term in the integrant that involves $\int uK(u)du$ is equal to zero since $\int uK(u)du=0$ by assumption. For terms with $\bm{\Sigma}(z)$, we have

[TABLE]

For terms that involve $\dot{\bm{\Sigma}}(z)$ and $\ddot{\bm{\Sigma}}(z^{\prime})$, we have

[TABLE]

since the maximum eigenvalue of $\bm{\Theta}(z)$ is bounded by $\rho$ by assumption. Thus, combining the above into (137), we have

[TABLE]

Next, we bound the second moment. By the Isserlis’ theorem in (135), and by taking the conditional expectation, we have

[TABLE]

Following a similar argument as in (138), we can derive

[TABLE]

Thus, we have

[TABLE]

Now we begin to bound the $\operatorname{\rm Cov}\{J_{z,jk}^{(1)}(Z),J_{z,jk}^{(2)}(Z)\}$. By using a similar argument as (138), we have

[TABLE]

Combining with (142) and (138), and using the covariance formula, we have that

[TABLE]

Using (136), (141) and (143), we have

[TABLE]

where the last inequality is because $\rho$ is smaller than the minimum eigenvalue of $\bm{\Sigma}(z)$ for any $z\in[0,1]$ and $\inf_{z\in[0,1]}\;f_{Z}(z)\geq\underline{f}_{Z}>0$ by Assumption 1. Since the lower bound above is uniformly true over $z,j,k$, the lemma is proven.

∎

Appendix G Proof of Theorem 5

In this section, we show that the proposed procedure in Algorithm 1 is able to control the type I error below a pre-specified level $\alpha$. We first define some notation that will be used throughout the proof of Theorem 3. Let $E^{*}(z)$ be the true edge set at $Z=z$. That is, $E^{*}(z)$ is the set of edges induced by the true inverse covariance matrix $\bm{\Theta}(z)$. Recall from Definition 4 that the critical edge set is defined as

[TABLE]

where $\mathscr{P}=\{E\subseteq V\times V\mid\mathcal{P}(G)=1\}$ is the class of edge sets satisfying the graph property $\mathcal{P}$.

Suppose that Algorithm 1 rejects the null hypothesis at the $T$th iteration. That is, there exists $z_{0}\in[0,1]$ such that $E_{T}(z_{0})\in\mathscr{P}$ but $E_{T-1}(z_{0})\notin\mathscr{P}$. To prove Theorem 3, we state the following two lemmas on the properties of critical edge set.

Lemma 11**.**

Let $E_{T}(z_{0})\in\mathscr{P}$ for some $z_{0}\in[0,1]$. Then, at least one rejected edge in $E_{T}(z_{0})$ is in the critical edge set $\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$.

Lemma 12**.**

Let $\bar{e}\in\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$ be the first rejected edge in the critical edge set $\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$. Suppose that $\bar{e}$ is rejected at the $l$th step of Algorithm 1. Then, $\mathcal{C}\{E^{*}(z),\mathcal{P}\}\subseteq\mathcal{C}\{E_{l-1}(z),\mathcal{P}\}$ for all $z\in[0,1]$.

The proofs of Lemmas 11 and 12 are deferred to Sections G.2 and G.3, respectively. We now provide the proof of Theorem 3.

G.1 Proof of Theorem 3

Suppose that Algorithm 1 rejects the null hypothesis at the $T$th iteration. That is, $E_{T}(z_{0})\in\mathscr{P}$ and $E_{T-1}(z_{0})\notin\mathscr{P}$. By Lemma 11, there is at least one edge in $E_{T}(z_{0})$ that is also in the critical edge set $\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$. We denote the first rejected edge in the critical edge set as $\bar{e}$, i.e., $\bar{e}\in\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$ and suppose that $\bar{e}$ is rejected at the $l$th iteration of Algorithm 1. We note that $l$ is not necessarily $T$. Thus, we have

[TABLE]

where the first inequality follows by Lemma 11, the second inequality follows from the $l$th step of Algorithm 1, and the last inequality follows directly from Lemma 12.

Under the null hypothesis, $\bm{\Theta}_{e}(z)=0$ for any $e\in\mathcal{C}\{E^{*}(z),\mathcal{P}\}$. By Theorem 2, we have

[TABLE]

as desired.

G.2 Proof of Lemma 11

To prove Lemma 11, it suffices to show that the intersection between the two sets $E_{T}(z_{0})$ and $\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$ is not an empty set, i.e., $E_{T}(z_{0})\cap\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}\neq\emptyset$. To this end, we let $F=E_{T}(z_{0})\cup E^{*}(z_{0})$ and let $E_{T}(z_{0})\setminus E^{*}(z_{0})=\{e_{1},e_{2},\ldots,e_{k}\}$. We note that the set $E_{T}(z_{0})\setminus E^{*}(z_{0})$ is not an empty set since $E_{T}(z_{0})\in\mathscr{P}$ but $E^{*}(z_{0})\notin\mathscr{P}$.

Using the fact that $\mathcal{P}$ is monotone and that $E_{T}(z_{0})\in\mathscr{P}$, we have $F\in\mathscr{P}$ since adding additional edges to $E_{T}(z_{0})$ does not change the graph property of $E_{T}(z_{0})$. Then, we have

[TABLE]

Since $E^{*}(z_{0})\notin\mathscr{P}$ and $F\in\mathscr{P}$, there must exists an edge set $\{e_{1},\ldots,e_{k_{0}}\}$ for $k_{0}\leq k$ that changes the graph property of $E^{*}(z_{0})$ from $E^{*}(z_{0})\notin\mathscr{P}$ to $E^{*}(z_{0})\cup\{e_{1},\ldots,e_{k_{0}}\}\in\mathscr{P}$.

Thus, there must exists at least an edge $\bar{e}\in\{e_{1},\ldots,e_{k_{0}}\}$ such that $\bar{e}\in\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$ since adding the set of edges $\{e_{1},\ldots,e_{k_{0}}\}$ changes the graph property of $E^{*}(z_{0})$. Also, $\bar{e}\in E_{T}(z_{0})$ by construction. Thus, we conclude that $E_{T}(z_{0})\cap\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}\neq\emptyset$.

G.3 Proof of Lemma 12

Let $\bar{e}\in\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$ be the first rejected edge in the critical edge set $\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$ for some $z_{0}\in[0,1]$. Suppose that $\bar{e}$ is rejected at the $l$th step of Algorithm 1. We want to show that $\mathcal{C}\{E^{*}(z),\mathcal{P}\}\subseteq\mathcal{C}\{E_{l-1}(z),\mathcal{P}\}$ for all $z\in[0,1]$. It suffices to show that $\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}\subseteq\mathcal{C}\{E_{l-1}(z_{0}),\mathcal{P}\}$. In other words, we want to prove that for any $e^{\prime}\in\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$, $e^{\prime}\in\mathcal{C}\{E_{l-1}(z_{0}),\mathcal{P}\}$. We first note the following fact

[TABLE]

By the definition of the critical edge set (144), we construct a set $E^{\prime}$ such that $E^{*}(z_{0})\supseteq E^{\prime}$, $E^{\prime}\in\mathscr{P}$, and $E^{\prime}\setminus\{e^{\prime}\}\notin\mathscr{P}$, for any $e^{\prime}\in\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$. By the definition of monotone property, we have $E^{\prime}\cup E_{l-1}(z_{0})\in\mathscr{P}$. Since $\mathcal{C}\{E^{\prime}\cup E_{l-1}(z_{0}),\mathcal{P}\}\subseteq\mathcal{C}\{E_{l-1}(z_{0}),\mathcal{P}\}$, to show that $e^{\prime}\in\mathcal{C}\{E_{l-1}(z_{0}),\mathcal{P}\}$, it is equivalent to showing $e^{\prime}\in\mathcal{C}\{E^{\prime}\cup E_{l-1}(z_{0}),\mathcal{P}\}$. That is, we want to show

[TABLE]

This is equivalent to showing

[TABLE]

There are two cases: (1) $E_{l-1}(z_{0})\setminus E^{\prime}=\emptyset$ and (2) $E_{l-1}(z_{0})\setminus E^{\prime}\neq\emptyset$. For the first case, (146) is true by the construction of $E^{\prime}$. For the second case, we prove by contradiction.

Suppose that $(E^{\prime}\setminus e^{\prime})\cup\{E_{l-1}(z_{0})\setminus E^{\prime}\}\in\mathscr{P}$. Let $E_{l-1}(z_{0})\setminus E^{\prime}=\{e^{\prime}_{1},\ldots,^{\prime}_{k}\}$. By the definition of monotone property, we have

[TABLE]

Since $E^{\prime}\setminus\{e^{\prime}\}\notin\mathscr{P}$ by construction, and that $(E^{\prime}\setminus\{e^{\prime}\})\cup(E_{l-1}(z_{0})\setminus E^{\prime})\in\mathscr{P}$, there must exists an edge set $\{e_{1},\ldots,e_{k_{0}}\}$ for $k_{0}\leq k$ that changes the graph property of $E^{\prime}\setminus\{e^{\prime}\}\notin\mathscr{P}$ to $(E^{\prime}\setminus\{e^{\prime}\})\cup\{e_{1}^{\prime},\ldots,e_{k_{0}}^{\prime}\}\in\mathscr{P}$.

Since $e^{\prime}_{k_{0}}\in E_{l-1}(z_{0})\setminus E^{\prime}$ and that $E^{*}(z_{0})\subseteq E^{\prime}$ by construction, we have $e^{\prime}_{k_{0}}\notin E^{*}(z_{0})$. Thus, $e^{\prime}_{k_{0}}\in\mathcal{C}\{E^{*}(z_{0}),\mathcal{P}\}$. This contradicts the fact that

[TABLE]

Appendix H Proof of Theorem 6

By the definition in (25), if $\bm{\Theta}(\cdot)\in\mathcal{G}_{1}(\theta;\mathcal{P})$, there exists an edge set $E^{\prime}_{0}$ and $z_{0}\in[0,1]$ satisfying

[TABLE]

and we will determine the magnitude tf constant $C$ later. We aim to show that $\mathcal{P}\{E_{0}^{\prime}\cap\mathcal{C}(\emptyset,\mathcal{P})\}=\mathcal{P}(E_{0}^{\prime})=1$. First, there exists a subgraph $E_{0}^{\prime\prime}\subset E_{0}^{\prime}$ such that $\mathcal{P}(E_{0}^{\prime\prime})=\mathcal{P}(E_{0}^{\prime})=1$ and for any $\widetilde{E}\subset E_{0}^{\prime\prime}$, $\mathcal{P}(\widetilde{E})=0$. We can construct such $E_{0}^{\prime\prime}$ by deleting edges from $E_{0}^{\prime}$ until it is impossible to further deleting any edge such that the property $\mathcal{P}$ is still true. By Definition 4, $E_{0}^{\prime\prime}\subseteq\mathcal{C}(\emptyset,\mathcal{P})$ and therefore $E_{0}^{\prime\prime}\subseteq E_{0}^{\prime}\cap\mathcal{C}(\emptyset,\mathcal{P})$. By monotone property, we have $\mathcal{P}\{E_{0}^{\prime}\cap\mathcal{C}(\emptyset,\mathcal{P})\}=\mathcal{P}(E_{0}^{\prime})=1$ since $\mathcal{P}(E^{\prime\prime}_{0})=\mathcal{P}(E^{\prime}_{0})=1$. Consider the following event

[TABLE]

According to Algorithm 1, the rejected set in the first iteration at $z_{0}$ is

[TABLE]

Under the event $\mathcal{E}_{1}$, we have $E_{0}^{\prime}\cap\mathcal{C}(\emptyset,\mathcal{P})\subseteq E_{1}(z_{0})$ and since $\mathcal{P}\{E_{0}^{\prime}\cap\mathcal{C}(\emptyset,\mathcal{P})\}=\mathcal{P}(E_{0}^{\prime})$, we have $\mathcal{P}\{E_{1}(z_{0})\}=\mathcal{P}(E_{0}^{\prime})=1$. Therefore,

[TABLE]

It suffices to bound $\mathbb{P}(\mathcal{E}_{1})$ then. We consider two events

[TABLE]

We have $\mathbb{P}(\mathcal{E}_{1})\geq\mathbb{P}(\mathcal{E}_{2}\cup\mathcal{E}_{3}).$ By Lemmas 2 and 3, we have

[TABLE]

Combining with (15) in Corollary 1, we have with probability at least $1-6/d$,

[TABLE]

For any fixed $\alpha\in(0,1)$ and sufficiently large $d,n$, as $\mathcal{C}(\emptyset,\mathcal{P})\subseteq V\times V$, we have

[TABLE]

Thus $\mathbb{P}(\mathcal{E}_{3})>1-6/d$. Similarly, we also have $\mathbb{P}(\mathcal{E}_{2})>1-3/d$. By (148), we have

[TABLE]

Therefore, we complete the proof of the theorem.

Appendix I Technical Lemmas on Covering Number

In this section, we present some technical lemmas on the covering number of some function classes. Lemma 13 provides an upper bound on the covering number for the class of function of bounded variation. Lemma 14 provides an upper bound on the covering number of a class of Lipschitz function. Lemma 15 provides an upper bound on the covering numbers for function classes generated from the product and addition of two function classes.

Lemma 13**.**

(Lemma 3 in Giné & Nickl 2009) Let $K:\mathbb{R}\rightarrow\mathbb{R}$ be a function of bounded variation. Define the function class $\mathcal{F}_{h}=\left[K\left\{(t-\cdot)/h)\right\}\mid t\in\mathbb{R}\right]$. Then, there exists $C_{K}<\infty$ independent of $h$ and of $K$ such that for all $0<\epsilon<1$,

[TABLE]

where $\|K\|_{\mathrm{TV}}$ is the total variation norm of the function $K$.

Lemma 14**.**

Let $f(l)$ be a Lipschitz function defined on $[0,1]$ such that $|f(l)-f(l^{\prime})|\leq L_{f}\cdot|l-l^{\prime}|$ for any $l,l^{\prime}\in[0,1]$. We define the constant function class $\mathcal{F}=\{g_{l}:=f(l)\mid l\in[0,1]\}$. For any probability measure $Q$, the covering number of the function class $\mathcal{F}$ satisfies

[TABLE]

where $\epsilon\in(0,1)$.

Proof.

Let $\mathcal{N}=\left\{\frac{i\epsilon}{L_{f}}\;|\;i=1,\ldots,\frac{L_{f}}{\epsilon}\right\}$. By definition of $\mathcal{N}$, for any $l\in[0,1]$, there exists an $l^{\prime}\in\mathcal{N}$ such that $|l-l^{\prime}|\leq\epsilon/L_{f}$. Thus, we have

[TABLE]

This implies that $\{g_{l}\;|\;l\in\mathcal{N}\}$ is an $\epsilon$-cover of the function class $\mathcal{F}$. To complete the proof, we note that the cardinality of the set $|\mathcal{N}|\leq L_{f}/\epsilon$. ∎

Lemma 15**.**

Let $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ be two function classes satisfying

[TABLE]

for some $C_{1},C_{2},a_{1},a_{2},v_{1},v_{2}>0$ and any $0<\epsilon<1$. Define $\|F_{\ell}\|_{\infty}=\underset{f\in\mathcal{F}_{\ell}}{\sup}\;\|f\|_{\infty}$ for $\ell=1,2$ and $U=\|\mathcal{F}_{1}\|_{\infty}\vee\|\mathcal{F}_{2}\|_{\infty}$. For the function classes $\mathcal{F}_{\times}=\{f_{1}f_{2}\;|\;f_{1}\in\mathcal{F}_{1},\;f_{2}\in\mathcal{F}_{2}\}$ and $\mathcal{F}_{+}=\{f_{1}+f_{2}\;|\;f_{1}\in\mathcal{F}_{1},\;f_{2}\in\mathcal{F}_{2}\}$, we have for any $\epsilon\in(0,1)$,

[TABLE]

and

[TABLE]

Lemma 16**.**

Let $w_{z}(u)=K_{h}(u-z)$. We define the function classes

[TABLE]

Given Assumptions 1-2, we have for any $\epsilon\in(0,1)$,

[TABLE]

and

[TABLE]

Moreover, let $k_{z}(u)=w_{z}(u)-\mathbb{E}\{w_{z}(Z)\}$ and let $\mathcal{K}=\{k_{z}(\cdot)\mid z\in[0,1]\}$. We have

[TABLE]

Proof.

The covering number for the function class $\mathcal{K}_{1}$ is obtained by an application of Lemma 13. To obtain the covering number for $\mathcal{K}_{2}$, we show that the constant function $\mathbb{E}\{w_{z}(Z)\}$ is Lipschitz. The covering number is obtained by applying Lemma 14. Finally, we note that the function class $\mathcal{K}$ is generated from the addition of the two function classes $\mathcal{K}_{1}$ and $\mathcal{K}_{2}$. The covering number can be obtained by an application of Lemma 15. The details are deferred to I.1. ∎

Lemma 17**.**

Let $g_{z,jk}(u,X_{ij},Y_{ik})=K_{h}(u-z)X_{ij}Y_{ik}$. We define the function classes

[TABLE]

Given Assumptions 1-2, for all $\epsilon\in(0,1)$,

[TABLE]

and

[TABLE]

with probability at least $1-1/d$. Moreover, let $q_{z,jk}(u,X_{ij},Y_{ik})=g_{z,jk}(u,X_{ij},Y_{ik})-\mathbb{E}\{g_{z,jk}(Z,X_{j},Y_{k})\}$ and let $\mathcal{G}_{jk}=\left\{q_{z,jk}(\cdot)\mid z\in[0,1]\right\}$. We have

[TABLE]

with probability at least $1-1/d$.

Proof.

The proof uses the same set of argument as in the proof of Lemma 16. The probability statement comes from the fact that we upper bound the random variable $X_{j}$ by $M_{X}\cdot\sqrt{\log d}$ for some constant $M_{X}>0$. The details are deferred to I.2. ∎

Lemma 18**.**

Let $J_{z,jk}^{(1)}(u,\bm{X}_{i},\bm{Y}_{i})=\sqrt{h}\cdot\left\{\bm{\Theta}_{j}(z)\right\}^{T}\cdot\left[K_{h}(u-z)\bm{X}_{i}\bm{Y}_{i}^{T}-\mathbb{E}\left\{K_{h}(Z-z)\bm{X}\bm{Y}^{T}\right\}\right]\cdot\bm{\Theta}_{k}(z)$ and let $\mathcal{J}^{(1)}_{jk}=\{J_{z,jk}^{(1)}\;|\;z\in[0,1]\}$. Given Assumptions 1-2, for all $\epsilon\in(0,1)$

[TABLE]

with probability at least $1-1/d$, where $C>0$ is a generic constant that does not depend on $d$, $h$, and $n$.

Proof.

The proof is deferred to I.3. ∎

Lemma 19**.**

Let $J_{z,jk}^{(2)}(u)=\sqrt{h}\cdot\left\{\bm{\Theta}_{j}(z)\right\}^{T}\cdot\left[K_{h}(u-z)-\mathbb{E}\left\{K_{h}(Z-z)\right\}\right]\cdot\bm{\Sigma}(z)\cdot\bm{\Theta}_{k}(z)$ and let $\mathcal{J}^{(2)}_{jk}=\{J_{z,jk}^{(2)}\mid z\in[0,1]\}.$ Given Assumptions 1-2, for all probability measures $Q$ on $\mathbb{R}$ and all $0<\epsilon<1$,

[TABLE]

where $C>0$ is a generic constant that does not depend on $d$, $h$, and $n$.

Proof.

We first note that $\mathcal{J}^{(2)}_{jk}$ is a function class generated from the product of two function classes $\mathcal{K}$ as in Lemma 16 and $\varTheta_{jk}=\{\bm{\Theta}_{jk}(z)\;|\;z\in[0,1]\}$. To obtain the covering number of $\varTheta_{jk}$, we show that the constant function $\bm{\Theta}_{jk}(z)$ is Lipschitz and apply Lemma 14. We then apply Lemma 15 to obtain the covering number of $\mathcal{J}^{(2)}_{jk}$. The details are deferred to I.4. ∎

I.1 Proof of Lemma 16

Let $w_{z}(u)=K_{h}(u-z)$ and that $k_{z}(u)=w_{z}(u)-\mathbb{E}\{w_{z}(Z)\}$. We first obtain the covering number for the function classes

[TABLE]

Then, we apply Lemma 15 to obtain the covering number of the function class

[TABLE]

Covering number for $\mathcal{K}_{1}$: By an application of Lemma 13, the covering number for $\mathcal{K}_{1}$ is

[TABLE]

Covering number for $\mathcal{K}_{2}$: First, note that $\mathbb{E}\{w_{z}(Z)\}=\int K_{h}(z-Z)f_{Z}(Z)dZ=(K_{h}*f_{Z})(z)$ is a function of $z$ generated by the convolution $(K_{h}*f_{Z})(z)$. By the property of the derivative of a convolution as in (38), we have

[TABLE]

where the last expression is obtained by an application of Young’s inequality. The expression in (151) depends on the quantity $\|\dot{K}_{h}\|_{1}$, which is equal to the following expression

[TABLE]

where the second inequality holds by a change of variable, and $\|K\|_{\mathrm{TV}}$ is the total variation of the function $K(\cdot)$. Substituting (152) into (151) and by Assumption 1, we have

[TABLE]

Thus, for any $z_{1},z_{2}\in[0,1]$, we have

[TABLE]

implying that $\mathbb{E}\{w_{z}(Z)\}$ is a Lipschitz continuous function with Lipschitz constant $h^{-1}\cdot\|K\|_{\mathrm{TV}}\cdot\bar{f}_{Z}$. By Lemma 14, an upper bound for the covering number of $\mathcal{K}_{2}$ is

[TABLE]

Covering number of the function class $\mathcal{K}$: The function class $\mathcal{K}$ can be written as

[TABLE]

By an application of Lemma 15 with $C_{1}=(2\cdot C_{K}\cdot\|K\|_{\mathrm{TV}})^{4}$, $C_{2}=2\cdot\bar{f}_{Z}\cdot\|K\|_{\mathrm{TV}}$, $a_{1}=a_{2}=h^{-1}$, $v_{1}=4$, and $v_{2}=1$, along with (150) and (154), we obtain

[TABLE]

I.2 Proof of Lemma 17

Throughout the proof, we condition on the event

[TABLE]

Since $\bm{X}$ and $\bm{Y}$ conditioned on $Z$ are Gaussian random variables, the event $\mathcal{A}$ occurs with probability at least $1-1/d$ for sufficiently large constant $M_{X}>0$.

Recall that $g_{z,jk}(u,X_{ij},Y_{ik})=K_{h}(u-z)X_{ij}Y_{ik}$ and that $q_{z,jk}(u,X_{ij},Y_{ik})=K_{h}(u-z)X_{ij}Y_{ik}-\mathbb{E}\{K_{h}(Z-z)X_{j}Y_{k}\}$. We first obtain the covering number of the function classes

[TABLE]

and

[TABLE]

Then, we apply Lemma 15 to obtain the covering number of the function class

[TABLE]

Covering number for $\mathcal{G}_{1,jk}$: Conditioned on the event $\mathcal{A}$ in (155), we have

[TABLE]

By an application of Lemma 13, the covering number for $\mathcal{G}_{1,jk}$ is

[TABLE]

Covering number for $\mathcal{G}_{2,jk}$: We now obtain the covering number for $\mathcal{G}_{2,jk}$ by showing that the function $\mathbb{E}\{g_{z,jk}(Z,X_{j},Y_{k})\}$ is Lipschitz. First, note that

[TABLE]

where $\varphi_{jk}(Z)=f_{Z}(Z)\cdot\bm{\Sigma}_{jk}(Z)$ and $K_{h}*\varphi_{jk}$ is the convolution between $K_{h}$ and $\varphi_{jk}$. Similar to (151)-(153), we have

[TABLE]

where the first inequality is obtained by an application of Young’s inequality, and the last expression is obtained by (152) and Assumptions 1-2.

Equation 157 implies that for any $z_{1},z_{2}\in[0,1]$,

[TABLE]

implying that $\mathbb{E}\{g_{z,jk}(Z,X_{j},Y_{k})\}$ is a Lipschitz continuous function with Lipschitz constant $h^{-1}\cdot\|K\|_{\mathrm{TV}}\cdot\bar{f}_{Z}\cdot M_{\sigma}$. By an application of Lemma 14, we have

[TABLE]

Covering number of the function class $\mathcal{G}_{jk}$: The function class $\mathcal{G}_{jk}$ can be written as

[TABLE]

By an application of Lemma 15 with $C_{1}=(2\cdot C_{K}\cdot\|K\|_{\mathrm{TV}}\cdot M^{2}_{X})^{4}$, $C_{2}=2\cdot\bar{f}_{Z}\cdot\|K\|_{\mathrm{TV}}\cdot M_{\sigma}$, $a_{1}=h^{-1}\cdot\log d$, $a_{2}=h^{-1}$, $v_{1}=4$, and $v_{2}=1$, along with (156) and (158), we obtain

[TABLE]

as desired.

I.3 Proof of Lemma 18

Similar to the proof of Lemma 17, we condition on the event

[TABLE]

The event $\mathcal{A}$ holds with probability at least $1-1/d$.

Recall that $J_{z,jk}^{(1)}(u,\bm{X}_{i},\bm{Y}_{i})=\sqrt{h}\cdot\left\{\bm{\Theta}_{j}(z)\right\}^{T}\cdot\left[K_{h}(u-z)\bm{X}_{i}\bm{Y}_{i}^{T}-\mathbb{E}\left\{K_{h}(Z-z)\bm{X}\bm{Y}^{T}\right\}\right]\cdot\bm{\Theta}_{k}(z)$ and let $\mathcal{J}^{(1)}_{jk}=\left\{J_{z,jk}^{(1)}\;|\;z\in[0,1]\right\}$. To obtain the covering number of the function class $\mathcal{J}_{jk}^{(1)}$, we consider bounding the covering number of a larger class of function. To this end, we define $\bm{\Phi}_{\omega}^{(1)}(u,\bm{X}_{i},\bm{Y}_{i})=\sqrt{h}\cdot\left[K_{h}(u-z)\bm{X}_{i}\bm{Y}_{i}^{T}-\mathbb{E}\left\{K_{h}(Z-z)\bm{X}\bm{Y}^{T}\right\}\right]$ to be a $d\times d$ matrix. We denote the $(j,k)$th element of $\Phi_{\omega}^{(1)}(u,\bm{X}_{i},\bm{Y}_{i})$ as $\Phi_{\omega,jk}^{(1)}(u,X_{ij},Y_{ik})=\sqrt{h}\cdot q_{\omega,jk}(u,X_{ij},Y_{ik})$, where $q_{\omega,jk}(u,X_{ij},Y_{ik})=K_{h}(u-\omega)X_{ij}Y_{ik}-\mathbb{E}\{K_{h}(Z-\omega)X_{j}Y_{k}\}$. We aim to obtain an $\epsilon$-cover $\mathcal{N}^{(1^{\prime})}$ for the following function class

[TABLE]

In other words, we show that for any $(\omega_{1},z_{1})\in[0,1]^{2}$, there exists $(\omega_{2},z_{2})\in\mathcal{N}^{(1^{\prime})}$ such that

[TABLE]

Given any $j,k\in[d]$, $\omega,\omega^{\prime},z,z^{\prime}\in[0,1]$, by the triangle inequality, we have

[TABLE]

We now obtain the upper bounds for $I_{1},I_{2}$, and $I_{3}$.

Upper bound for $I_{1}$ and $I_{3}$: First, we note that by Holder’s inequality, we have

[TABLE]

Since $\bm{\Theta}(z)\in\mathcal{U}(s,M,\rho)$, we have

[TABLE]

Moreover, for any $z_{1},z_{2}\in[0,1]$, we have

[TABLE]

where the second to the last inequality follows from the fact that $\bm{\Theta}(z)\in\mathcal{U}(s,M,\rho)$ and the last inequality follows from Assumption 2. Finally, from (74) and the definition of $\Phi_{\omega_{1},jk}^{(1)}(\cdot)=\sqrt{h}\cdot q_{\omega_{1},jk}(\cdot)$, we have

[TABLE]

Combining (161)-(163), we have

[TABLE]

We note that $I_{3}$ can be upper bounded the same way as $I_{1}$.

Upper bound for $I_{2}$: Recall from (160) that

[TABLE]

where the inequality holds by Holder’s inequality and the definition of $\bm{\Phi}_{\omega}^{(1)}(u,\bm{X}_{i},\bm{Y}_{i})$. Let $\varPhi_{jk}^{(1)}=\left\{\sqrt{h}\cdot q_{\omega,jk}(\cdot)\;|\;\omega\in[0,1]\right\}$ and recall from Lemma 17 that we constructed an $\epsilon$-cover $\mathcal{N}^{(1^{\prime\prime})}\subset[0,1]$ for the function class $\varPhi_{jk}^{(1)}$ with cardinality $\left|\mathcal{N}^{(1^{\prime\prime})}\right|=\left(\frac{4\cdot\|K\|_{\mathrm{TV}}\cdot C_{K}^{4/5}\cdot\bar{f}_{Z}^{1/5}\cdot M_{\sigma}^{1/5}\cdot M^{8/5}_{X}\cdot\log^{4/5}d}{\sqrt{h}\cdot\epsilon}\right)^{5}$. Since the construction of the $\epsilon$-cover in Lemma 17 is independent of the indices $j$ and $k$, we have that for any $j,k\in[d]$ and $\omega_{1}\in[0,1]$, there exists a $\omega_{2}\in\mathcal{N}^{(1^{\prime\prime})}$ such that

[TABLE]

Thus, by (161) and (165), we have

[TABLE]

Covering number of the function class $\mathcal{J}_{jk}^{(1)}$: Since $\mathcal{J}_{jk}^{(1)}\subset\mathcal{J}_{jk}^{(1^{\prime})}$, the covering number of $\mathcal{J}_{jk}^{(1)}$ is upper bounded by the covering number of $\mathcal{J}_{jk}^{(1^{\prime})}$. It suffices to construct an $\epsilon$-cover of the function class $\mathcal{J}_{jk}^{(1^{\prime})}$. In the following, we show that $\mathcal{N}^{(1^{\prime})}=\mathcal{N}^{(1^{\prime\prime})}\times\left\{i\cdot\epsilon\cdot\sqrt{h}\;|\;i=1,\ldots,\frac{1}{\epsilon\cdot\sqrt{h}}\right\}$ is an $\epsilon$-cover of $\mathcal{J}_{jk}^{(1^{\prime})}$. For any $(\omega_{1},z_{1})\in[0,1]^{2}$, there exists $(\omega_{2},z_{2})\in\mathcal{N}^{(1^{\prime})}$ such that (165) holds and that $|z_{1}-z_{2}|\leq\sqrt{h}\cdot\epsilon$. Thus, combining (164) and (166), we have

[TABLE]

where $C>0$ is a generic constant that does not depend on $d$, $h$, and $n$.

Thus, we have

[TABLE]

Since $\mathcal{J}_{jk}^{(1)}\subset\mathcal{J}_{jk}^{(1^{\prime})}$, the above expression implies that

[TABLE]

as desired.

I.4 Proof of Lemma 19

First, we note that

[TABLE]

where $k_{z}(u)=K_{h}(u-z)-\mathbb{E}\{K_{h}(Z-z)\}$. Let $\mathcal{J}^{(2)}_{jk}=\left\{J_{z,jk}^{(2)}\;|\;z\in[0,1]\right\}$. Furthermore, recall that $\mathcal{K}=\left\{k_{z}(\cdot)\;|\;z\in[0,1]\right\}$ and let $\varTheta_{jk}=\left\{\bm{\Theta}_{jk}(z)\;|\;z\in[0,1]\right\}$. The function class $\mathcal{J}^{(2)}_{jk}$ can be written as

[TABLE]

It suffices to obtain the covering number for $\mathcal{K}$ and $\varOmega_{jk}$, and apply Lemma 15.

Covering number of the function class $\mathcal{K}$: By Lemma 16, we have

[TABLE]

Covering number of the function class $\varTheta_{jk}$: We show that $\bm{\Theta}_{jk}(z)$ is Lipschitz, and apply Lemma 14 to obtain the covering number for $\varTheta_{jk}$. Similar to (162), for any $z_{1},z_{2}\in[0,1]$, we have

[TABLE]

where the last inequality follows from Assumption 2. Since $\bm{\Theta}_{jk}(z)$ is $\rho^{2}\cdot d\cdot M_{\sigma}$-Lipschitz, by Lemma 14, we have

[TABLE]

Covering number of the function class $\mathcal{J}_{jk}^{(2)}$: We now apply Lemma 15 to obtain the covering number of $\mathcal{J}_{jk}^{(2)}$. Applying Lemma 15 with $a_{1}=d$, $v_{1}=1$, $C_{1}=M_{\sigma}\cdot\rho^{2}$, $a_{2}=h^{-1}$, $v_{2}=5$, $C_{2}=\left(4\cdot\|K\|_{\mathrm{TV}}\cdot C_{K}^{4/5}\cdot\bar{f}_{Z}^{1/5}\right)^{5}$, and $U=\frac{2}{h}\cdot\|K\|_{\infty}$, along with (169) and (170), we have

[TABLE]

where $C>0$ is a generic constant that does not depend on $n,d$, and $h$. This implies that

[TABLE]

as desired.

Appendix J Technical Lemmas on Empirical Process

In this section, we present some existing tools on empirical process. The following lemma states that the supreme of any empirical process is concentrated near its mean. It follows directly from Theorem 2.3 in Bousquet (2002).

Lemma 20**.**

(Theorem A.1 in Van de Geer 2008) Let $X_{1},\ldots,X_{n}$ be independent random variables and let $\mathcal{F}$ be a function class such that there exists $\eta$ and $\tau^{2}$ satisfying

[TABLE]

Define

[TABLE]

Then, for any $t>0$,

[TABLE]

The above inequality involves evaluating the expectation of the supreme of the empirical process. The following lemma follows directly from Theorem 3.12 in Koltchinskii (2011). It provides an upper bound on the expectation of the supreme of the empirical process as a function of its covering number.

Lemma 21**.**

(Lemma F.1 in Lu et al. 2015a) Assume that the functions in $\mathcal{F}$ defined on $\mathcal{X}$ are uniformly bounded by a constant $U$ and $F(\cdot)$ is the envelope of $\mathcal{F}$ such that $|f(x)|\leq F(x)$ for all $x\in\mathcal{X}$ and $f\in\mathcal{F}$. Let $\sigma^{2}_{P}=\underset{f\in\mathcal{F}}{\sup}\;\mathbb{E}(f^{2})$. Let $X_{1},\ldots,X_{n}$ be $\mathrm{i.i.d.}$ copies of the random variables $X$. We denote the empirical measure as $\mathbb{P}_{n}=\frac{1}{n}\sum_{i\in[n]}\delta_{X_{i}}$. If for some $A,V>0$ and for all $\epsilon>0$ and $n\geq 1$, the covering entropy satisfies

[TABLE]

then for any $\mathrm{i.i.d.}$ sub-gaussian mean zero random variables $\xi_{1},\ldots,\xi_{n}$, there exists a universal constant $C$ such that

[TABLE]

Furthermore, we have

[TABLE]