A test of weak separability for multi-way functional data, with application to brain connectivity studies

TL;DR

This paper introduces a formal test for weak separability in multi-way functional data, facilitating the analysis of complex brain connectivity patterns with applications to magnetoencephalography signals.

Contribution

It proposes a novel statistical test for weak separability in multi-way functional data, linking tensor methods to brain connectivity analysis.

Findings

The test effectively distinguishes weakly separable structures.

Application to MEG data reveals meaningful brain connectivity insights.

The asymptotic distribution of the test statistic is derived.

Abstract

This paper concerns the modeling of multi-way functional data where double or multiple indices are involved. We introduce a concept of weak separability. The weakly separable structure supports the use of factorization methods that decompose the signal into its spatial and temporal components. The analysis reveals interesting connections to the usual strongly separable covariance structure, and provides insights into tensor methods for multi-way functional data. We propose a formal test for the weak separability hypothesis, where the asymptotic null distribution of the test statistic is a chi-square type mixture. The method is applied to study brain functional connectivity derived from source localized magnetoencephalography signals during motor tasks.

Tables

Table 1. Table 1: Rejection rates for the χ 2 superscript 𝜒 2 \chi^{2} type mixture weak separability test procedure, using V 1 subscript 𝑉 1 V_{1} and choosing ( P n , K n ) subscript 𝑃 𝑛 subscript 𝐾 𝑛 (P_{n},K_{n}) with the fraction of variance explained procedure (FVE) or as ( 2 , 2 ) 2 2 (2,2) , ( 3 , 3 ) 3 3 (3,3) , or ( 4 , 4 ) 4 4 (4,4)

Scenario	Normal				Multivariate $t$
$n=50$	FVE	(2,2)	(3,3)	(4,4)	FVE	(2,2)	(3,3)	(4,4)
$H_0$	0.055	0.020	0.020	0.040	0.025	0.020	0.005	0.045
$\text{cov}({\chi }_{12},{\chi }_{21})=0.065$	0.715	0.785	0.740	0.755	0.440	0.445	0.395	0.370
$\text{cov}({\chi }_{12},{\chi }_{21})=0.13$	1.000	1.000	1.000	1.000	0.935	0.965	0.940	0.940
3 nonzero terms	1.000	1.000	1.000	1.000	0.960	0.990	0.990	0.965
$n=100$
$H_0$	0.035	0.075	0.045	0.050	0.025	0.040	0.020	0.010
$\text{cov}({\chi }_{12},{\chi }_{21})=0.065$	0.985	0.985	0.985	0.985	0.800	0.810	0.785	0.710
$\text{cov}({\chi }_{12},{\chi }_{21})=0.13$	1.000	1.000	1.000	1.000	1.000	0.990	0.990	0.990
3 nonzero terms	1.000	1.000	1.000	1.000	1.000	0.990	1.000	1.000
$n=500$
$H_0$	0.060	0.060	0.040	0.055	0.020	0.045	0.020	0.045
$\text{cov}({\chi }_{12},{\chi }_{21})=0.065$	1.000	1.000	1.000	1.000	0.995	0.995	1.000	0.990
$\text{cov}({\chi }_{12},{\chi }_{21})=0.13$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
3 nonzero terms	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000

Table 2. Table 2: Rejection rates for the χ 2 superscript 𝜒 2 \chi^{2} type mixture weak separability test procedure, using V 2 subscript 𝑉 2 V_{2} and choosing ( P n , K n ) subscript 𝑃 𝑛 subscript 𝐾 𝑛 (P_{n},K_{n}) with the fraction of variance explained procedure (FVE) or as ( 2 , 2 ) 2 2 (2,2) , ( 3 , 3 ) 3 3 (3,3) , or ( 4 , 4 ) 4 4 (4,4)

Scenario	Normal				Multivariate $t$
$n=50$	FVE	(2,2)	(3,3)	(4,4)	FVE	(2,2)	(3,3)	(4,4)
$H_0$	0.030	0.025	0.030	0.025	0.015	0.030	0.015	0.005
$\text{cov}({\chi }_{12},{\chi }_{21})=0.055$	0.515	0.845	0.440	0.465	0.305	0.555	0.225	0.205
$\text{cov}({\chi }_{12},{\chi }_{21})=0.11$	0.995	0.995	0.995	0.990	0.850	0.955	0.825	0.770
3 nonzero terms	1.000	1.000	1.000	1.000	0.965	0.985	0.950	0.970
$n=100$
$H_0$	0.045	0.055	0.035	0.040	0.010	0.050	0.040	0.020
$\text{cov}({\chi }_{12},{\chi }_{21})=0.055$	0.920	0.990	0.930	0.920	0.625	0.900	0.605	0.500
$\text{cov}({\chi }_{12},{\chi }_{21})=0.11$	1.000	1.000	1.000	1.000	0.990	1.000	0.965	0.955
3 nonzero terms	1.000	1.000	1.000	1.000	0.995	1.000	1.000	0.980
$n=500$
$H_0$	0.045	0.065	0.025	0.040	0.025	0.065	0.050	0.035
$\text{cov}({\chi }_{12},{\chi }_{21})=0.055$	1.000	1.000	1.000	1.000	0.970	1.000	0.995	0.995
$\text{cov}({\chi }_{12},{\chi }_{21})=0.11$	1.000	1.000	1.000	1.000	1.000	1.000	0.990	0.995
3 nonzero terms	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000

Table 3. Table 3: Rejection rates for the non-studentized empirical bootstrap weak separability test procedure, using V 1 subscript 𝑉 1 V_{1} and choosing ( P n , K n ) subscript 𝑃 𝑛 subscript 𝐾 𝑛 (P_{n},K_{n}) with the fraction of variance explained procedure (FVE) or as ( 2 , 2 ) 2 2 (2,2) , ( 3 , 3 ) 3 3 (3,3) , or ( 4 , 4 ) 4 4 (4,4)

Scenario	Normal				Multivariate $t$
$n=50$	FVE	(2,2)	(3,3)	(4,4)	FVE	(2,2)	(3,3)	(4,4)
$H_0$	0.035	0.025	0.035	0.025	0.010	0.010	0.005	0.005
$\text{cov}({\chi }_{12},{\chi }_{21})=0.065$	0.670	0.680	0.650	0.640	0.350	0.375	0.330	0.300
$\text{cov}({\chi }_{12},{\chi }_{21})=0.13$	1.000	1.000	1.000	1.000	0.890	0.900	0.885	0.880
3 nonzero terms	1.000	1.000	1.000	1.000	0.920	0.920	0.920	0.915
$n=100$
$H_0$	0.070	0.060	0.060	0.060	0.015	0.010	0.015	0.015
$\text{cov}({\chi }_{12},{\chi }_{21})=0.065$	0.990	0.995	0.985	0.985	0.735	0.785	0.700	0.680
$\text{cov}({\chi }_{12},{\chi }_{21})=0.13$	1.000	1.000	1.000	1.000	0.970	0.970	0.965	0.960
3 nonzero terms	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
$n=500$
$H_0$	0.065	0.050	0.060	0.060	0.055	0.050	0.055	0.055
$\text{cov}({\chi }_{12},{\chi }_{21})=0.065$	1.000	1.000	1.000	1.000	0.995	0.995	0.995	0.995
$\text{cov}({\chi }_{12},{\chi }_{21})=0.13$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
3 nonzero terms	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000

Table 4. Table 4: Rejection rates for the non-studentized empirical bootstrap weak separability test procedure, using V 2 subscript 𝑉 2 V_{2} and choosing ( P n , K n ) subscript 𝑃 𝑛 subscript 𝐾 𝑛 (P_{n},K_{n}) with the fraction of variance explained procedure (FVE) or as ( 2 , 2 ) 2 2 (2,2) , ( 3 , 3 ) 3 3 (3,3) , or ( 4 , 4 ) 4 4 (4,4)

Scenario	Normal				Multivariate $t$
$n=50$	FVE	(2,2)	(3,3)	(4,4)	FVE	(2,2)	(3,3)	(4,4)
$H_0$	0.065	0.045	0.050	0.045	0.015	0.020	0.010	0.010
$\text{cov}({\chi }_{12},{\chi }_{21})=0.055$	0.420	0.785	0.390	0.380	0.220	0.395	0.160	0.150
$\text{cov}({\chi }_{12},{\chi }_{21})=0.11$	0.970	1.000	0.965	0.965	0.730	0.865	0.690	0.675
3 nonzero terms	1.000	1.000	1.000	1.000	0.895	0.925	0.885	0.885
$n=100$
$H_0$	0.025	0.010	0.020	0.020	0.035	0.035	0.030	0.020
$\text{cov}({\chi }_{12},{\chi }_{21})=0.055$	0.950	1.000	0.955	0.955	0.585	0.855	0.575	0.515
$\text{cov}({\chi }_{12},{\chi }_{21})=0.11$	1.000	1.000	1.000	1.000	0.980	0.995	0.980	0.975
3 nonzero terms	1.000	1.000	1.000	1.000	0.985	0.990	0.985	0.985
$n=500$
$H_0$	0.030	0.065	0.030	0.030	0.010	0.035	0.005	0
$\text{cov}({\chi }_{12},{\chi }_{21})=0.055$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
$\text{cov}({\chi }_{12},{\chi }_{21})=0.11$	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000
3 nonzero terms	1.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000

Equations

$$X(s,t)=\mu (s,t)+\sum _{k=1}^{\infty }\sum _{j=1}^{\infty }{\chi }_{jk}{\psi }_j(s){\varphi }_k(t),$$

$$C_{\mathcal{S}}(s,u)={\int }_{\mathcal{T}}C(s,t;u,t)dt,C_{\mathcal{T}}(t,v)={\int }_{\mathcal{S}}C(s,t;s,v)ds.$$

$$X(s,t)=\mu (s,t)+\sum _{j=1}^{\infty }\sum _{k=1}^{\infty }{\tilde{\chi }}_{jk}{f}_j(s)g_k(t),$$

$$C(s,t;u,v)=\sum _{j=1}^{\infty }\sum _{k=1}^{\infty }{\eta }_{jk}{\psi }_j(s){\varphi }_k(t){\psi }_j(u){\varphi }_k(v),$$

$${\mathrm{rank}}_{+}(V)=\min \{\ell :V=V_1+\cdots +V_{\ell };V_i\ge 0,\mathrm{rank}(V_i)=1,\mathrm{for}\mathrm{all}i\},$$

$${\chi }_{i,jk}={\int }_{\mathcal{T}}{\int }_{\mathcal{S}}\{X_i(s,t)-\mu (s,t)\}{\psi }_j(s){\varphi }_k(t)dsdt,$$

$${\hat{\chi }}_{i,jk}={\int }_{\mathcal{T}}{\int }_{\mathcal{S}}\{X_i(s,t)-\bar{X}(s,t)\}{\hat{\psi }}_j(s){\hat{\varphi }}_k(t)dsdt,$$

$$T_n(j,k,j^{\prime },k^{\prime })=n^{-1/2}\sum _{i=1}^n{\hat{\chi }}_{i,jk}{\hat{\chi }}_{i,j^{\prime }{k}^{\prime }}\{(j,k)\ne (j^{\prime },k^{\prime })\},$$

$$C_n=(1/n)\sum _{i=1}^n(X_i-\bar{X})\otimes (X_i-\bar{X}).$$

$$\sum _{j=1,2,\dots }\{\mathrm{E}(⟨X,e_j{⟩}^4){\}}^{1/4}<\infty .$$

$$T_n(j,k,j^{\prime },k^{\prime })=tr\left[\{({\psi }_j\otimes {\psi }_{j^{\prime }})\tilde{\otimes }({\varphi }_k\otimes {\varphi }_{k^{\prime }})\}{\mathcal{Z}}_n\right]+o_p(1),$$

$$T_n(j,k,j,k^{\prime })=$$

$$+$$

$$T_n(j,k,j^{\prime },k)=$$

$$+$$

$$S_n=\sum _{j,j^{\prime }=1,\dots ,P_n;k,k^{\prime }=1,\dots ,K_n;(j,k)<(j^{\prime },k^{\prime })}\{T_n(j,k,j^{\prime },k^{\prime }){\}}^2,$$

$$\text{FVE}(P_n,K_n)=\frac{\frac{1}{n}{\sum }_{i=1}^n{\sum }_{j=1}^{P_n}{\sum }_{k=1}^{K_n}{\hat{\chi }}_{i,jk}^2}{\frac{1}{n}{\sum }_{i=1}^n{\sum }_{j=1}^{\infty }{\sum }_{k=1}^{\infty }{\hat{\chi }}_{i,jk}^2}.$$

$$\frac{\mathrm{E}(\Vert X-\tilde{X}{\Vert }^2)}{\mathrm{E}(\Vert X-\mu {\Vert }^2)}=1-\frac{{\sum }_{j=1}^{P_n}{\sum }_{k=1}^{K_n}{\eta }_{jk}}{{\sum }_{j=1}^{\infty }{\sum }_{k=1}^{\infty }{\eta }_{jk}}.$$

$$\text{FVE}(P_n,K_n)\gtrsim {\text{FVE}}_{\mathcal{S}}(P_n)+{\text{FVE}}_{\mathcal{T}}(K_n)-1,$$

$${\hat{\chi }}_{i,jk}^{\ast }={\int }_{\mathcal{T}}{\int }_{\mathcal{S}}\{X_i^{\ast }(s,t)-{\bar{X}}^{\ast }(s,t)\}{\hat{\psi }}_j^{\ast }(s){\hat{\varphi }}_k^{\ast }(t)dsdt,$$

$$T_n^{\ast }(j,k,j^{\prime },k^{\prime })=n^{-1/2}\sum _{i=1}^n{\hat{\chi }}_{i,jk}^{\ast }{\hat{\chi }}_{i,j^{\prime }{k}^{\prime }}^{\ast }.$$

$$S_n^{\ast }=\sum _{j,j^{\prime }=1,\dots ,P_n;k,k^{\prime }=1,\dots ,K_n;(j,k)<(j^{\prime },k^{\prime })}\{T_n^{\ast }(j,k,j^{\prime },k^{\prime })-T_n(j,k,j^{\prime },k^{\prime }){\}}^2.$$

$$\mathrm{PLV}(\mathrm{s},\mathrm{t})=(1/{\mathrm{n}}_{\mathrm{T}})|\sum _{\mathrm{k}=1}^{{\mathrm{n}}_{\mathrm{T}}}{\mathrm{e}}^{\mathrm{i}\{{\mathrm{B}}_{1,\mathrm{k}}(\mathrm{s},\mathrm{t})-{\mathrm{B}}_{2,\mathrm{k}}(\mathrm{s},\mathrm{t})\}}|.$$

$$C(s,t;u,v)=\sum _{j=1}^{\infty }\sum _{k=1}^{\infty }{\eta }_{jk}{f}_j(s)g_k(t)f_j(u)g_k(v),$$

$$C_{\mathcal{S}}(s,u)=$$

$$=$$

$${\int }_{\mathcal{T}}\sum _{j=1}^{\infty }\sum _{k=1}^{\infty }|{\eta }_{jk}{f}_j(s)g_k(t)f_j(u)g_k(t)|dt$$

$$\le$$

$$=$$

$$\le$$

$$C_{\mathcal{S}}(s,u)={\int }_{\mathcal{T}}C(s,t;u,t)dt=a{C}_1(s,u){\int }_{\mathcal{T}}{C}_2(t,t)dt=a{C}_1(s,u).$$

$$C(s,t;u,v)=$$

Full text

A test of weak separability for multi-way functional data, with application to brain connectivity studies

Brian Lynch, Kehui Chen

University of Pittsburgh

Abstract:

This paper concerns the modeling of multi-way functional data where double or multiple indices are involved. We introduce a concept of weak separability. The weakly separable structure supports the use of factorization methods that decompose the signal into its spatial and temporal components. The analysis reveals interesting connections to the usual strongly separable covariance structure, and provides insights into tensor methods for multi-way functional data. We propose a formal test for the weak separability hypothesis, where the asymptotic null distribution of the test statistic is a chi-square type mixture. The method is applied to study brain functional connectivity derived from source localized magnetoencephalography signals during motor tasks.

Keywords: asymptotics; functional principal component; hypothesis testing; marginal kernel; separable covariance; spatio-temporal data; tensor product

Brian Lynch is a PhD Student, Department of Statistics, University of Pittsburgh, Pittsburgh, PA 15260 (E-mail: bcl28@pitt.edu). Kehui Chen is Assistant Professor, Department of Statistics, University of Pittsburgh, Pittsburgh, PA 15260 (E-mail: khchen@pitt.edu). This work is partially supported by NSF1612458.

1   Introduction

Traditional functional data analysis usually concerns data recorded over a continuum, such as growth curves. Dense and regularly-observed functional data can be recorded in a matrix with dimension $n\times T$, where $n$ is the number of subjects and $T$ is the number of grid points observed for each subject. Multi-way functional data refers to an extension where multiple indices are involved and data can be recorded in a tensor with dimension at least three. Examples include brain imaging data where for each subject $i=1,\ldots,n$, we have observations $X_{i}(s,t)$, with a spatial index $s\in\mathbb{R}^{d}$ and a time index $t\in\mathbb{R}^{1}$. Other examples include repeatedly or longitudinally observed functional data, such as data obtained from tracking apps where subjects’ 24-hour profiles of activities are recorded every day. This type of data can be represented by $X_{i}(s,t)$, where $s$ denotes the day and $t$ denotes the time within a day.

As multi-way functional data become more common with modern techniques, the modeling of this type of data attracts increasing interest. Assume the individual observations $X_{i}(s,t)$ are independent and identically distributed realizations of a random process $X\in L^{2}(\mathcal{S}\times\mathcal{T})$, $s\in\mathcal{S}\subseteq\mathbb{R}^{d_{1}}$, $t\in\mathcal{T}\subseteq\mathbb{R}^{d_{2}}$, with mean $\mu$ and continuous covariance operator $C$. When we can do so without confusion, we use the same symbol for the covariance operator and its kernel function. A well-established tool in functional data analysis is functional principal component analysis. When applied to the multi-way process $X$, functional principal component analysis is based on the Karhunen–Loève representation $X(s,t)=\mu(s,t)+\sum_{l=1}^{\infty}Z_{l}h_{l}(s,t)$, where $Z_{l}\ (l=1,2,\ldots)$ are the (random) uncorrelated coefficients, and $h_{l}(s,t)\ (l=1,2,\ldots)$ are the eigenfunctions of the covariance operator $C$.

To alleviate the difficulties associated with modeling the $(2d_{1}+2d_{2})$-dimensional full covariance function $C(s,t,u,v)$ and characterizing the $(d_{1}+d_{2})$-dimensional eigenfunctions, one usually seeks dimension reduction through factorization of the signal into its spatial and temporal components. Chen et al. (2017) proposed product functional principal component analysis,

[TABLE]

where $\psi_{j}(s)\ (j=1,2,\ldots)$ and $\phi_{k}(t)\ (k=1,2,\ldots)$ are the eigenfunctions of the marginal covariance operators in $L^{2}(\mathcal{S})$ and $L^{2}(\mathcal{T})$, with corresponding marginal kernels

[TABLE]

Here $C_{\mathcal{S}}(s,u)=\sum_{j=1}^{\infty}\lambda_{j}\psi_{j}(s)\psi_{j}(u)$ and $C_{\mathcal{T}}(t,v)=\sum_{k=1}^{\infty}\gamma_{k}\phi_{k}(t)\phi_{k}(v)$, where $\lambda_{1}\geq\lambda_{2}\geq\cdots$ and $\gamma_{1}\geq\gamma_{2}\geq\cdots$ are the eigenvalues. The $\chi_{jk}=\int_{\mathcal{T}}\int_{\mathcal{S}}\{X(s,t)-\mu(s,t)\}\psi_{j}(s)\phi_{k}(t)dsdt$ are the marginal projection scores. To be precise, we should first consider expanding $X(s,t)$ in terms of completed versions of the bases of marginal eigenfunctions, but since it can be shown that the scores $\chi_{jk}$ associated with the extra functions needed to complete the bases are 0, the expansion of $X(s,t)$ in Equation 1 holds.

The above product functional principal component analysis representation is the same as the Karhunen–Loève representation if one makes the separable covariance assumption $C(s,t;u,v)=aC_{1}(s,u)C_{2}(t,v)$, which we call strong separability in contrast to the weak separability that will be proposed in this paper. However, if strong separability is not assumed, the marginal eigenfunctions no longer carry optimal efficiency guarantees (Aston et al., 2012), and can only be proven to have near-optimality under appropriate assumptions (Chen et al., 2017). Moreover, unlike the $Z_{l}$ in multi-way functional principal component analysis, the scores $\chi_{jk}$ are not guaranteed to be mutually uncorrelated.

Factorization of the signal into its spatial ($s$) and temporal ($t$) components, justified using a vague notion of spatial-temporal separability, is a common strategy used in many methods in image analysis and multi-way functional data analysis (Zhang & Zhou, 2005; Lu et al., 2006; Huang et al., 2009; Chen & Müller, 2012; Hung et al., 2012; Allen et al., 2014; Chen et al., 2015, 2017). Despite their empirical success, the rigorous characterization of this separable feature is still mainly restricted to the scope of strong separability, i.e., when the covariance $C(s,t;u,v)$ is separable. There is a large amount of literature on strong separability in related fields (Lu & Zimmerman, 2005; Fuentes, 2006; Srivastava et al., 2009; Hoff et al., 2011; Horváth & Kokoszka, 2012). Tests for strong separability in functional data settings have been proposed recently (Aston et al., 2017; Constantinou et al., 2017).

In this paper, we propose a new concept of weak separability for the process $X$, which can be rigorously tested. We show that under weak separability the eigenfunctions of the full covariance $C(s,t;u,v)$ can be written as tensor products of the marginal eigenfunctions, i.e., $\psi_{j}\otimes\phi_{k}$. This means the Karhunen–Loève representation is the same as the product representation in Equation 1, just as if we had strong separability, and to perform functional principal component analysis we only need to calculate the marginal covariances $C_{\mathcal{S}}(s,u)$ and $C_{\mathcal{T}}(t,v)$ instead of the full covariance $C(s,t;u,v)$. The analysis reveals that if $C(s,t;u,v)$ is separable, then the process $X$ is weakly separable, but the converse is not necessarily true. Indeed, weak separability is a much weaker assumption than separable covariance.

We develop a test for weak separability based on the empirical correlations between the estimated scores $\hat{\chi}_{i,jk}$ and $\hat{\chi}_{i,j^{\prime}k^{\prime}}\ \{i=1,\ldots,n;\ (j,k)\neq(j^{\prime},k^{\prime})\}$. Although the $\hat{\chi}_{i,jk}$ are $n^{1/2}$-consistent estimators of the $\chi_{i,jk}$, the test statistics based on the $\hat{\chi}_{i,jk}$ have different null distributions from their counterparts using the $\chi_{i,jk}$ due to non-negligible estimation errors. The proofs involve expansions of the differences between the estimated marginal eigenfunctions and their true values, i.e., $\hat{\psi}_{j}-\psi_{j}$ and $\hat{\phi}_{k}-\phi_{k}$, as well as multi-way tensor products with indices $(j,k,j^{\prime},k^{\prime})$. A series of careful derivations are carried out to characterize the asymptotic null distribution of the test statistic, which is found to be a $\chi^{2}$ type mixture. No Gaussian assumption on $X$ is imposed. We apply the testing procedure to brain imaging data, where frequency and time-based functional connectivity is constructed from source localized magnetoencephalography signals. The test result supports the use of product functional principal component analysis methods and reveals interesting features about brain connectivity over time and frequency.

2   Weak separability: concepts and properties

For $\mathcal{S}\subseteq\mathbb{R}^{d_{1}}$ and $\mathcal{T}\subseteq\mathbb{R}^{d_{2}}$, we consider the space of square integrable surfaces $L^{2}(\mathcal{S}\times\mathcal{T})$ with the standard inner product $\langle f,g\rangle=\int_{\mathcal{T}}\int_{\mathcal{S}}f(s,t)g(s,t)dsdt$ and the corresponding norm $\|\cdot\|$. The data can be viewed as realizations of a random element $X\in L^{2}(\mathcal{S}\times\mathcal{T})$, which we assume has well defined mean function $\mu$ and covariance operator $C$. We assume the covariance is continuous, and $\mathcal{S}$ and $\mathcal{T}$ are compact. Unless otherwise noted, these assumptions are used in all the lemmas and theorems.

For orthonormal bases $f_{j}\ (j=1,2,\ldots)$ in $L^{2}(\mathcal{S})$ and $g_{k}\ (k=1,2,\ldots)$ in $L^{2}(\mathcal{T})$, the product functions $f_{j}(s)g_{k}(t)\ (j=1,2,\ldots;\ k=1,2,\ldots)$ form an orthonormal basis of $L^{2}(\mathcal{S}\times\mathcal{T})$. We can then have

[TABLE]

where $\tilde{\chi}_{jk}=\int_{\mathcal{T}}\int_{\mathcal{S}}\{X(s,t)-\mu(s,t)\}f_{j}(s)g_{k}(t)dsdt$.

Definition of weak separability: $X(s,t)$* is weakly separable if there exist orthonormal bases $f_{j}\ (j=1,2,\ldots)$ and $g_{k}\ (k=1,2,\ldots)$ such that ${\rm cov}(\tilde{\chi}_{jk},\tilde{\chi}_{j^{\prime}k^{\prime}})=0$ for $j\neq j^{\prime}$ or $k\neq k^{\prime}$, i.e., the scores $\tilde{\chi}_{jk}\ (j=1,2,\ldots;\ k=1,2,\ldots)$ are uncorrelated with each other.*

In the following, we list several important properties of weak separability, which make this concept attractive in many applications. Detailed proofs are given in the appendix.

Lemma 1

If $X$ is weakly separable, the pair of bases $f_{j}\ (j=1,2,\ldots)$ and $g_{k}\ (k=1,2,\ldots)$ that satisfies weak separability is unique up to a sign, and $f_{j}(s)\equiv\psi_{j}(s)$ and $g_{k}(t)\equiv\phi_{k}(t)$, where $\psi_{j}(s)$ and $\phi_{k}(t)$ are the eigenfunctions of the marginal kernels $C_{\mathcal{S}}(s,u)$ and $C_{\mathcal{T}}(t,v)$ as defined in Equation 2. Moreover,

[TABLE]

where $\eta_{jk}={\rm var}(\langle X-\mu,\psi_{j}\otimes\phi_{k}\rangle)$, and the convergence is absolute and uniform.

Lemma 1 shows that natural basis functions for the factorized spatial and temporal effects are eigenfunctions of the marginal kernels. Under weak separability, the eigenfunctions of the covariance $C(s,t;u,v)$ can be written as tensor products of the marginal eigenfunctions, i.e., $\psi_{j}\otimes\phi_{k}$, which could result in a substantial dimension reduction in applications. Lemma 1 also allows us to test the weak separability assumption (see Section 3).

Lemma 2

Strong separability, defined as $C(s,t;u,v)=aC_{1}(s,u)C_{2}(t,v)$ with identifiability constraints $\int_{\mathcal{S}}C_{1}(s,s)ds=1$ and $\int_{\mathcal{T}}C_{2}(t,t)dt=1$, implies weak separability of $X$. And up to a constant scaling, $C_{1}$ and $C_{2}$ are the same as the marginal kernels.

Lemma 2 shows that strong separability is a special case of weak separability, and the following Lemma 3 further illustrates that weak separability is much more flexible than strong separability.

Lemma 3

Define the array $V=(\eta_{jk},\ j=1,2,\ldots;\ k=1,2,\ldots)$. Strong separability is weak separability with the additional assumption that ${\rm rank}(V)=1$. Moreover, under strong separability $V=a\Lambda\Gamma^{T}$, where $\Lambda=(\lambda_{1},\lambda_{2},\ldots)^{T}$ and $\Gamma=(\gamma_{1},\gamma_{2},\ldots)^{T}$ are the eigenvalues of the marginal kernels, and $a=1/\int_{\mathcal{T}}\int_{\mathcal{S}}C(s,t;s,t)dsdt$ is a normalization constant.

When the covariance $C$ is not strongly separable but the process $X$ is weakly separable, we can show that the covariance function is a sum of $L$ separable components, $C(s,t;u,v)=\sum_{l=1}^{L}a^{l}C_{\mathcal{S}}^{l}(s,u)C_{\mathcal{T}}^{l}(t,v),$ where $L\geq{\rm rank}(V)>1$ is the nonnegative rank, defined as

[TABLE]

where $V_{i}\geq 0$ means that $V_{i}$ is entry-wise nonnegative. In applications where one relies on the separable structure of the covariance for ease of computation and interpretation, for example in applications involving the inverse of the covariance, it is not clear whether and how one can modify the concept to work under the weak separability assumption ($L$ additive separable terms). We defer this to future research.

3   Test of weak separability

3.1   Background

Assume we have a sample of independent and identically distributed smooth processes $X_{i}(s,t)\sim X(s,t)$, and the marginal projection scores

[TABLE]

where $\psi_{j}(s)$ and $\phi_{k}(t)$ are the eigenfunctions of the marginal covariances. By the definition of weak separability and Lemma 1, testing weak separability is the same as testing the covariance structure of the marginal projection scores, i.e., $H_{0}:$ ${\rm cov}(\chi_{jk},\chi_{j^{\prime}k^{\prime}})=0$ for $j\neq j^{\prime}$ or $k\neq k^{\prime}$.

The problem of testing covariance structure is a classic problem in multivariate analysis. Suppose we have independent and identically distributed copies of a $p$-variate random variable, with mean $\mu$ and covariance matrix $\Sigma$, and we want to test the null hypothesis that $\Sigma$ is diagonal. Under the traditional multivariate setting where $p$ is fixed and does not increase with $n$, likelihood ratio methods can be used to test the diagonality of $\Sigma$ (Anderson, 1984). The high-dimensional problem has been studied in the context that $p/n\rightarrow\gamma\in(0,\infty)$ or even for $p$ much larger than $n$ (Ledoit & Wolf, 2002; Liu et al., 2008; Cai et al., 2011; Lan et al., 2015). If we were to observe the sample values $\chi_{i,jk}$, a sensible test statistic could be based on the off-diagonal terms of the empirical covariance, i.e., $n^{-1/2}\sum_{i=1}^{n}\chi_{i,jk}\chi_{i,j^{\prime}k^{\prime}}.$ However, unlike in the traditional covariance testing problem, we do not directly observe the sample values $\chi_{i,jk}$. Instead they are estimated from the sample curves $X_{i}(s,t)\ (i=1,\ldots,n)$ as

[TABLE]

where $\bar{X}(s,t)=(1/n)\sum_{i=1}^{n}X_{i}(s,t)$, and $\hat{\psi}_{j}$ and $\hat{\phi}_{k}$ are eigenfunctions of the estimated marginal covariances $\hat{C}_{\mathcal{S}}(s,u)=(1/n)\sum_{i=1}^{n}\int_{\mathcal{T}}\{X_{i}(s,t)-\bar{X}(s,t)\}\{X_{i}(u,t)-\bar{X}(u,t)\}dt$ and $\hat{C}_{\mathcal{T}}(t,v)=(1/n)\sum_{i=1}^{n}\int_{\mathcal{S}}\{X_{i}(s,t)-\bar{X}(s,t)\}\{X_{i}(s,v)-\bar{X}(s,v)\}ds$. In practice, if the data for each subject are observed on arbitrarily dense and equally spaced grid points, and recorded in matrices $X_{i}\ (i=1,\ldots,n)$, the above estimators can be simplified as $\hat{C}_{\mathcal{S}}=(1/n)\sum_{i=1}^{n}(X_{i}-\bar{X})(X_{i}-\bar{X})^{T}$, and $\hat{C}_{\mathcal{T}}=(1/n)\sum_{i=1}^{n}(X_{i}-\bar{X})^{T}(X_{i}-\bar{X})$. The data cannot immediately be written as matrices if the argument $s$ has dimension greater than 1, but as long as the observations are dense in $\mathcal{S}$ one can vectorize them along a certain ordering of $s$, compute the marginal covariances, and reorganize back accordingly.

Although we can prove that the $\hat{\chi}_{i,jk}$ are $n^{1/2}$-consistent estimators of the $\chi_{i,jk}$, test statistics based on $n^{-1/2}\sum_{i=1}^{n}\hat{\chi}_{i,jk}\hat{\chi}_{i,j^{\prime}k^{\prime}}$ have different null distributions from their counterparts using the $\chi_{i,jk}$, and in the following we derive the asymptotic distribution of the former.

3.2   The test statistic and its properties

Let $H$ be a real separable Hilbert space, with inner product $\langle\cdot,\cdot\rangle$. Following standard definitions, we denote the space of bounded linear operators on $H$ as $\mathcal{B}(H)$, the space of Hilbert–Schmidt operators on $H$ as $\mathcal{B}_{HS}(H)$, and the space of trace-class operators on $H$ as $\mathcal{B}_{Tr}(H)$. For any trace-class operator $T\in\mathcal{B}_{Tr}(H)$, we define its trace by $tr(T)=\sum_{i=1,2,\ldots}\langle Te_{i},e_{i}\rangle$, where $e_{i}\ (i=1,2,\ldots)$ is an orthonormal basis of $H$, and it is easy to see that this definition is independent of the choice of basis.

For $H_{1}$ and $H_{2}$ two real separable Hilbert spaces, we use $\otimes$ as the standard tensor product, i.e., for $x_{1}\in H_{1}$ and $x_{2}\in H_{2}$, $(x_{1}\otimes x_{2})$ is the operator from $H_{2}$ to $H_{1}$ defined by $(x_{1}\otimes x_{2})y=\langle x_{2},y\rangle x_{1}$ for any $y\in H_{2}$. With a bit of abuse of notation, we let $H=H_{1}\otimes H_{2}$ denote the tensor product Hilbert space, which contains all finite sums of $x_{1}\otimes x_{2}$, with inner product $\langle x_{1}\otimes x_{2},y_{1}\otimes y_{2}\rangle=\langle x_{1},y_{1}\rangle\langle x_{2},y_{2}\rangle$, for $x_{1},y_{1}\in H_{1}$ and $x_{2},y_{2}\in H_{2}$. For $C_{1}\in\mathcal{B}(H_{1})$ and $C_{2}\in\mathcal{B}(H_{2})$, we let $C_{1}\tilde{\otimes}C_{2}$ denote the unique bounded linear operator on $H_{1}\otimes H_{2}$ satisfying $C_{1}\tilde{\otimes}C_{2}(x_{1}\otimes x_{2})=C_{1}x_{1}\otimes C_{2}x_{2}$ for all $x_{1}\in H_{1},\,x_{2}\in H_{2}.$

We define

[TABLE]

and $\mathcal{Z}_{n}=n^{1/2}(C_{n}-C)$, where the sample covariance operator is defined as

[TABLE]

The following two conditions are needed for the main theorem below and the corollary following it:

Condition 1: For some orthonormal basis $e_{j}\ (j=1,2,\ldots)$ of $L^{2}(\mathcal{S}\times\mathcal{T})$,

[TABLE]

Condition 2: For some integers $P$ and $K$, we have $\delta_{P}=\min_{j=1,\ldots,P}(\lambda_{j}-\lambda_{j+1})>0$ and $\delta_{K}=\min_{k=1,\ldots,K}(\gamma_{k}-\gamma_{k+1})>0$.

Remark: According to Proposition 5 of Mas (2006), Condition 1 implies that $\mathcal{Z}_{n}$ converges to a Gaussian random element in $\mathcal{B}_{Tr}\{L^{2}(\mathcal{S}\times\mathcal{T})\}$.

Theorem 4

Assume Conditions 1 and 2 hold, and that $X$ is weakly separable. For $j,j^{\prime}=1,\ldots,P$ and $k,k^{\prime}=1,\ldots,K$ as defined in Condition 2, we have

(i) for $j\neq j^{\prime}$ and $k\neq k^{\prime}$,

[TABLE]

(ii) for $j=j^{\prime}$ and $k\neq k^{\prime}$,

[TABLE]

(iii) for $j\neq j^{\prime}$ and $k=k^{\prime}$,

[TABLE]

where $Id_{1}$ and $Id_{2}$ are identity operators on $L^{2}(\mathcal{S})$ and $L^{2}(\mathcal{T})$, respectively.

Remark: Since $n^{1/2}tr\left\{(\psi_{j}\otimes\psi_{j^{\prime}})\tilde{\otimes}(\phi_{k}\otimes\phi_{k^{\prime}})C\right\}$ is zero under the null hypothesis, the first term in each case of the above theorem is the same as $n^{1/2}tr\left\{(\psi_{j}\otimes\psi_{j^{\prime}})\tilde{\otimes}(\phi_{k}\otimes\phi_{k^{\prime}})C_{n}\right\}=n^{-1/2}\sum_{i=1}^{n}\chi_{i,jk}\chi_{i,j^{\prime}k^{\prime}}$, i.e., the counterpart of $T_{n}$ as if we had the true marginal projection scores. The second and third terms, if they exist, are non-negligible estimation errors.

Corollary 5

Assume Conditions 1 and 2 hold, and that $X$ is weakly separable. For different sets of $(j,k,j^{\prime},k^{\prime})$, $j,j^{\prime}=1,\ldots,P$; $k,k^{\prime}=1,\ldots,K$, satisfying $(j,k)\neq(j^{\prime},k^{\prime})$, the $T_{n}(j,k,j^{\prime},k^{\prime})$’s are asymptotically jointly Gaussian with mean zero and covariance structure $\Theta$. The formula for $\Theta$ is given in the proof.

3.3   Tests based on $\chi^{2}$ type mixtures

Lemma 6

For $j\neq j^{\prime}$, $\sum_{k=1}^{\infty}{\rm E}(\chi_{jk}\chi_{j^{\prime}k})=0$, and for $k\neq k^{\prime}$, $\sum_{j=1}^{\infty}{\rm E}(\chi_{jk}\chi_{jk^{\prime}})=0$. This also holds in the empirical version such that for $j\neq j^{\prime}$, $\sum_{k=1}^{\infty}T_{n}(j,k,j^{\prime},k)=0$, and for $k\neq k^{\prime}$, $\sum_{j=1}^{\infty}T_{n}(j,k,j,k^{\prime})=0$.

The above lemma does not assume weak separability. Recall the fact that principal component scores in traditional functional principal component analysis are uncorrelated. This lemma is a generalized result for the marginal projection scores.

Due to this linear relationship between the different terms of $T_{n}$, the asymptotic covariance $\Theta$ will be degenerate, and thus the statistic we consider is the sum of squares of the terms of $T_{n}$ without normalizing by the covariance. In practice, for suitably chosen $P_{n}$ and $K_{n}$, we use the statistic defined as

[TABLE]

where $(j,k)<(j^{\prime},k^{\prime})$ means $(j-1)*K_{n}+k<(j^{\prime}-1)*K_{n}+k^{\prime}$.

Take $T_{n}$ to be a long vector of length $m=P_{n}K_{n}(P_{n}K_{n}-1)/2$ created by stacking all of the $T_{n}(j,k,j^{\prime},k^{\prime})\ \{j,j^{\prime}=1,\ldots,P_{n};\ k,k^{\prime}=1,\ldots,K_{n};\ (j,k)<(j^{\prime},k^{\prime})\}$. Then by Corollary 5, $T_{n}\sim N_{m}(0,\Theta)$ under $H_{0}$, where we now take $\Theta$ to be a covariance matrix. Define the spectral decomposition of $\Theta$ as $\Theta=UQU^{T}$, where $Q$ is diagonal with diagonal entries $\sigma_{1},\ldots,\sigma_{m}$, which are the eigenvalues of $\Theta$ ordered from largest to smallest, and $U=[u_{1}\ldots\ u_{m}]$, where the $u_{i}$ are orthonormal column vectors. By Lemma 6, some of the $\sigma_{i}$ are 0. Since $S_{n}=\|T_{n}\|^{2}=\|U^{T}T_{n}\|^{2}$ and $U^{T}T_{n}\sim N_{m}(0,Q)$, we can write $S_{n}=\sum_{i=1}^{m}\sigma_{i}A_{i}$ where the $A_{i}$ are independent and identically distributed $\chi_{1}^{2}$, i.e., the null distribution of $S_{n}$ is a weighted sum of $\chi^{2}$ distributions, which we call a $\chi^{2}$ type mixture.

The Welch–Satterthwaite approximation for a $\chi^{2}$ type mixture (Zhang, 2013) approximates $S_{n}\sim\beta\chi_{d}^{2}$ and determines $\beta$ and $d$ from matching the first 2 cumulants (the mean and the variance). This results in $\beta=tr(\Theta^{2})/tr(\Theta)$ and $d=\{tr(\Theta)\}^{2}/tr(\Theta^{2})$. By using a plug-in estimator of $\Theta$, we can approximate the P-value for our test as an upper tail probability of $\beta\chi_{d}^{2}$. When the first $(P_{n},K_{n})$ terms do not satisfy weak separability, we have $S_{n}{\rightarrow}\infty$ in probability by noticing that for at least one set of $(j,k,j^{\prime},k^{\prime})$, the first term in Equation 5 in the proof of Theorem 4 is on the order of $n^{1/2}$.

The consistent selection of $(P_{n},K_{n})$ for hypothesis testing is a challenging problem. The optimal choice of $(P_{n},K_{n})$ needs to be defined according to the problem at hand and subsequent analysis of interest. Here we focus on the subspace where the subsequent product functional principal component analysis is going to be carried out. A criterion we will use to evaluate a given choice of $(P_{n},K_{n})$ is the fraction of variance explained by the first $P_{n}$ and $K_{n}$ components, defined as

[TABLE]

This definition can be justified by noting its relation to the normalized mean squared $L^{2}$ loss of the truncated process $\tilde{X}(s,t)=\mu(s,t)+\sum_{j=1}^{P_{n}}\sum_{k=1}^{K_{n}}\chi_{jk}\psi_{j}(s)\phi_{k}(t)$. In particular,

[TABLE]

The latter term is approximated by our definition of fraction of variance explained. The above equality only relies on the orthogonality of the eigenfunctions, not the weak separability assumptions. Thus, it still makes sense to consider this definition of fraction of variance explained even when $H_{0}$ is not true.

We also define the marginal fractions of variance explained as $\textrm{FVE}_{S}(P_{n})=\sum_{j=1}^{P_{n}}\hat{\lambda}_{j}/\sum_{j=1}^{\infty}\hat{\lambda}_{j}$ and $\textrm{FVE}_{T}(K_{n})=\sum_{k=1}^{K_{n}}\hat{\gamma}_{k}/\sum_{k=1}^{\infty}\hat{\gamma}_{k}$, where the $\hat{\lambda}_{j}$ are the eigenvalues of $\hat{C}_{\mathcal{S}}$ and the $\hat{\gamma}_{k}$ are the eigenvalues of $\hat{C}_{\mathcal{T}}$. In practice the infinite sums in the denominators of $\textrm{FVE}({P_{n},K_{n}})$, $\textrm{FVE}_{S}(P_{n})$, and $\textrm{FVE}_{T}(K_{n})$ will have to be replaced with the largest number of terms that can reasonably be considered nonzero.

Noting that $\sum_{j=1}^{\infty}{\rm E}(\chi_{jk}^{2})=\gamma_{k}$, $\sum_{k=1}^{\infty}{\rm E}(\chi_{jk}^{2})=\lambda_{j}$ and $\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}{\rm E}(\chi_{jk}^{2})=\sum_{j=1}^{\infty}\lambda_{j}=\sum_{k=1}^{\infty}\gamma_{k}$, we have

[TABLE]

subject to estimation error (to see, for example, that $\sum_{k=1}^{\infty}{\rm E}(\chi_{jk}^{2})=\lambda_{j}$, take $j=j^{\prime}$ in the proof of Lemma 6, with no need to assume weak separability). Therefore, we propose the following fraction of variance explained procedure: First choose $P_{n}$ and $K_{n}$ such that the marginal fractions of variance explained are at least 90%. If this choice results in $\textrm{FVE}(P_{n},K_{n})\geq 90\%$, use these values of $P_{n}$ and $K_{n}$. If not, use the values of $P_{n}$ and $K_{n}$ that have marginal fractions of variance explained at least 95%, in which case $\textrm{FVE}(P_{n},K_{n})$ is expected to be above 90%.

3.4   Bootstrap approximation

As an alternative to asymptotic approximation, we can also consider a bootstrap approach to approximate the distribution of the test statistic. Theorem 4 provides theoretical support for the use of the following empirical bootstrap procedure (Van Der Vaart & Wellner, 1996). Our simulations show that the asymptotic approximation based on the $\chi^{2}$ type mixture has very satisfactory performance. We still present the bootstrap approximation here since it is generally applicable to similar tests where the asymptotic null distributions do not have closed form.

At each step, draw a random sample from the data $X_{1},\ldots,X_{n}$ with replacement. Denote this sample as $X^{*}_{1},\ldots,X^{*}_{n}$. Let

[TABLE]

where $\bar{X}^{*}$ is the sample mean of the $X^{*}_{i}$, and the $\hat{\psi}^{*}_{j}$ and $\hat{\phi}^{*}_{k}$ are the eigenfunctions of the estimated marginal covariances of the $X^{*}_{i}$. The signs of the $\hat{\psi}^{*}_{j}$ and $\hat{\phi}^{*}_{k}$ are chosen to minimize $\|\hat{\psi}_{j}^{*}-\hat{\psi}_{j}\|$ and $\|\hat{\phi}_{k}^{*}-\hat{\phi}_{k}\|$, respectively. Let

[TABLE]

The empirical bootstrap test statistic is calculated as

[TABLE]

This procedure is repeated $B$ times, and the P-value is approximated as the proportion of bootstrap test statistics $S_{n}^{*}$ that are larger than the test statistic $S_{n}$.

Theorem 3.9.13 in Van Der Vaart & Wellner (1996) can be used to prove the validity of the bootstrap procedure, i.e., the conditional random laws (given data) of $S_{n}^{*}$ are asymptotically consistent almost surely for estimating the laws of $S_{n}$, under the null hypothesis. By Theorem 4, we have that under the null hypothesis, $T_{n}$ can be written as $\Phi_{P}^{\prime}\{n^{1/2}(\mathcal{P}_{n}-P)\}+o(1)$ and $T_{n}^{*}-T_{n}$ can be written as $\Phi_{P}^{\prime}\{n^{1/2}(\mathcal{P}_{n}^{*}-\mathcal{P}_{n})\}+o(1)$, where $\Phi_{P}^{\prime}$ is a linear continuous mapping that depends on the three different cases in Theorem 4. Thus, Theorem 3.9.13 applies.

Other than the above non-studentized empirical bootstrap based on $S_{n}$, we have also considered a bootstrap procedure based on a marginally studentized test statistic, in which we divide each term in $S_{n}$ by its corresponding estimated variance $\hat{\theta}(j,k,j^{\prime},k^{\prime})$ (which is the plug-in estimate of ${\theta}(j,k,j^{\prime},k^{\prime})$, the diagonal entry of $\Theta$ corresponding to the asymptotic variance of $T_{n}(j,k,j^{\prime},k^{\prime})$). However, we have found this procedure is much more time consuming, and requires substantially higher sample size to achieve high power, in comparison to the non-studentized empirical bootstrap method. This is not unexpected, since the form of ${\theta}(j,k,j^{\prime},k^{\prime})$ is very complicated and plug-in estimation adds extra variability. Therefore, we do not recommend the marginally studentized empirical bootstrap method.

4   Numerical study

In this section, we conduct numerical experiments to check the size and power of the proposed test for weak separability. Following the notation in Section 3, we are basically testing if the ${\rm cov}(\chi_{jk},\chi_{j^{\prime}k^{\prime}})$ are all zero for $(j,k)\neq(j^{\prime},k^{\prime})$. We consider two different choices for the joint distribution of the ${\chi}_{jk}$. The first is the multivariate normal and the second is the multivariate $t$ distribution. The diagonal values of $\Sigma$, the covariance matrix of the ${\chi}_{jk}$, are determined by the matrix $V=\{{\rm var}(\chi_{jk}),\,j,k=1,\ldots,8\}$. We consider two different choices for $V$, which we denote as $V_{1}$ and $V_{2}$ (specified later). Under $H_{0}$ (when all of the off-diagonal values of $\Sigma$ are 0), $V_{1}$ corresponds to a strongly separable covariance structure, while $V_{2}$ corresponds to a weakly separable structure that is not strongly separable. To study power, for a given choice of $V_{1}$ or $V_{2}$, we take $\textrm{cov}(\chi_{i,12},\chi_{i,21})$ to be the largest positive value such that $\Sigma$ is positive definite, and we also consider half of this value. Alternatively, we let 3 off-diagonal terms, $\textrm{cov}(\chi_{i,12},\chi_{i,21})$, $\textrm{cov}(\chi_{i,11},\chi_{i,22})$, and $\textrm{cov}(\chi_{i,13},\chi_{i,31})$, take their largest positive values such that $\Sigma$ is positive definite.

Empirical rejection rates at the .05 significance level from 200 simulations runs for $n=50,100,500$ are shown in Tables 1 through 4. The fraction of variance explained method described in Section 3.3 ends up with $P_{n}=3$ and $K_{n}=2$ in most trials, and we show results with $(P_{n},K_{n})$ chosen by this procedure, as well as directly setting $(P_{n},K_{n})=(2,2)$, $(3,3)$, or $(4,4)$ for all trials. We see that both the $\chi^{2}$ type mixture approximation and the empirical bootstrap procedure are able to control the type I error under all scenarios and achieve very good power as $n$ or the signal increase, although the empirical bootstrap is slightly less powerful for small $n$. Even when the chosen nonzero off-diagonal covariance terms are set to their maximum values, the other off-diagonal covariance terms of $\Sigma$ are zero, and so the signal is moderate. The test procedures are slightly less powerful in the multivariate $t$ case for small $n$, as the asymptotics likely come into play more quickly for the normal data. The rejection rates are in general stable across different choices of $(P_{n},K_{n})$; although $(P_{n},K_{n})=(2,2)$ seems to have higher power in some cases, the power stabilizes to a reasonable value for larger $(P_{n},K_{n})$.

Details of the simulation settings are as follows: We generate independent samples of data $X_{i}(s,t)=\sum_{j=1}^{8}\sum_{k=1}^{8}\chi_{i,jk}\psi_{j}(s)\phi_{k}(t)\ (i=1,\ldots,n)$, where the scores $\chi_{i,jk}$ are mean 0 random variables that we generate directly. We let $s$ and $t$ take values from 0 to 1 on an evenly spaced grid of 20 points. For the $\psi_{j}$ we use the functions $\psi_{j}(s)=-2^{1/2}\cos\{\pi(n+1)s\}$ for $j$ odd and $\psi_{j}(s)=2^{1/2}\sin(\pi ns)$ for $j$ even. We define the $\phi_{k}$ by taking the first 3 B-spline functions produced by Matlab’s spcol function using order 4 with knots at 0, 0.5, and 1, combining these with the first 5 $\psi_{j}$ as defined above, and orthonormalizing using Gram–Schmidt.

Let ${\chi}_{i}$ be the vector of $\chi_{i,jk}$ for $j,k=1,\ldots,8$. We simulate each ${\chi}_{i}$ independently from either $N(0,\Sigma)$ or the multivariate $t$ distribution. In the latter case, we first simulate a vector $x$ of length 64 from $N(0,\Sigma)$. One standard definition of a multivariate $t$ vector is $x/(u/v)^{1/2}$, where $u$ is a chi-square random variable with $v$ degrees of freedom that is independent of $x$. However, we use $x/\{u/(v-2)\}^{1/2}$ as our multivariate $t$ vector so that its covariance matrix is $\Sigma$. We take $v=6$ in our simulations. For each of 200 trials, we simulate data $X_{i}(s,t)\ (i=1,\ldots,n)$ in the manner described above, estimate the marginal projection scores, calculate the test statistic, and obtain P-values from the test procedures as described in Section 3, using $B=1000$ for the bootstrap procedure.

We choose $V_{1}$ and $V_{2}$ to both give $\lambda_{j}=\exp\{1.2(9-j)\}/\left\{\sum_{j^{\prime}=1}^{8}\exp(1.2j^{\prime})\right\}\ (j=1,\ldots,8)$ and $\gamma_{k}=\exp\{1.6(9-k)\}/\left\{\sum_{k^{\prime}=1}^{8}\exp(1.6k^{\prime})\right\}\ (k=1,\ldots,8)$ as the eigenvalues of the marginal covariances $C_{\mathcal{S}}$ and $C_{\mathcal{T}}$. $V_{1}$ is defined as the rank 1 matrix computed from the outer product of the vectors of $\lambda_{j}$ and $\gamma_{k}$, while $V_{2}$ is a rank 2 matrix with first 2 rows multiples of each other and rows 3 through 8 multiples of each other.

5   Application to brain connectivity studies

Brain imaging analysis is an area where functional data increasingly arise. An important goal in brain imaging studies is to analyze functional connectivity between different regions of the brain. We focus on magnetoencephalography (MEG), which measures neuronal activity by recording magnetic fields generated within the brain. We use MEG data collected by the Human Connectome Project, a study that has compiled a large amount of high quality multi-modal neural data, much of which is freely accessible at https://db.humanconnectome.org (Van Essen et al., 2013; WU-Minn HCP, 2017). We will focus on the motor task data, particularly the trials where subjects moved their right hand. The signal for each trial is recorded from -1.2 to 1.2 seconds in intervals of about 2 ms, where time 0 corresponds to the start of the motion. In the preprocessed sensor-level data (see WU-Minn HCP (2017) for preprocessing details), there are 61 subjects with motor data, and the subjects have an average of 75.38 trials. Our connectivity analysis will focus on two regions of interest: the left primary motor cortex and the right inferior parietal lobule. These regions of interest are spatially separated, likely activated during the task, and potentially functionally connected.

As the MEG sensors are distant from the brain, directly using their signals to represent regions of interest can lead to spurious connectivity measurements. This is due to the volume conduction/field spread problem, in which each sensor picks up the activity of several sources, as well as the common input problem, in which a common source provides input to a pair of signals that do not directly interact (Larson-Prior et al., 2013; Bastos & Schoffelen, 2015). For these reasons, we will use source reconstruction to estimate the signals arising from the cortical surface. Source reconstruction is common in MEG analysis, but it is an inverse problem on which constraints must be placed to obtain a unique solution (Pizzella et al., 2014). The source reconstruction method we use is minimum norm estimation as implemented in the Matlab package FieldTrip (Lin et al., 2004; Oostenveld et al., 2010).

MEG signals are inherently oscillatory, and synchronization at certain frequency ranges of the activity in different regions has been shown to be related to tasks performed by the brain (Pizzella et al., 2014). To study how frequency-based coupling between regions of interest changes over the course of a task, we calculate the time-frequency representations of their signals based on Morlet wavelets, using FieldTrip’s ft_freqanalysis function. The time-frequency representation of a signal is its representation at time $t$ and frequency $s$ as a complex number $A(s,t)e^{iB(s,t)}$, where $A(s,t)$ is the amplitude and $B(s,t)$ is the phase.

Given time-frequency representations $A_{1,k}(s,t)e^{iB_{1,k}(s,t)}$ and $A_{2,k}(s,t)e^{iB_{2,k}(s,t)}$ for two signals recorded in trial $k$, $k=1,\ldots,n_{T}$, we use the phase locking value (Lachaux et al., 1999) to measure their connectivity, calculated as

[TABLE]

The phase locking value takes values from 0 to 1, with 1 indicating complete phase synchrony over trials and 0 indicating no phase synchrony. The phase locking value, like many connectivity measures, is based on an analogue of the cross-correlation function called the coherence, but the phase locking value disregards the amplitudes and considers only the magnitude of the average of the phase differences as unit vectors in the complex plane. The phase locking value has gained popularity due to the belief that phase differences reveal more about functional connectivity than changes in amplitude (Lachaux et al., 1999; Aydore et al., 2013; Bastos & Schoffelen, 2015).

Because we calculate the time-frequency representation using wider time windows for lower frequencies, we are limited in how low of frequencies we can consider, and our preliminary results for power show a lack of activity above 50 Hz. Thus, we calculate the time-frequency representation from 8 to 50 Hz, corresponding to the alpha to gamma low frequency bands. In each trial, the motion usually lasts no longer than about 0.75 seconds. The signal at a time period shortly before time 0 is of interest, as it can represent brain activity when subjects have received the movement cue but have not yet reacted to it. However, the trials are not disjoint, so the signal at times further before 0 overlaps with the signal from the previous trial. Thus, in our analysis we will consider times for each trial on the range of -0.25 to 0.75 seconds. Calculating phase locking value between the two source-reconstructed signals corresponding to our regions of interest, the data we analyze is $\textrm{PLV}_{i}(s,t)$, where $i=1,\ldots,61$; $8\leq s\leq 50$; and $-0.25\leq t\leq 0.75$.

Figure 1 shows the average of the phase locking value matrices over all subjects, as well as slightly smoothed phase locking value matrices for 3 randomly selected subjects. The level of activity seems to vary between subjects. The average phase locking value displays higher synchrony near the beginning of the movement (time 0) in the alpha and beta bands, and the individual subjects’ plots also show higher values near time 0. However, the average has small values overall, which indicates high variability between subjects, and points to the need to study covariance structure and modes of variation.

For the phase locking value data described above, using the fraction of variance explained procedure described in Section 3.3, we choose the number of components to be $P_{n}=7$ and $K_{n}=7$. We apply the weak separability test using both the $\chi^{2}$ type mixture approximation (P-value = 0.5293) and the empirical bootstrap (P-value = 0.9260). The weak separability test does not reject the null hypothesis of weak separability, which supports the use of functional principal component analysis based on products of marginal eigenfunctions. We also apply the strong separability test of Aston et al. (2017) via their R package covsep (Tavakoli, 2016). The resulting P-values are $1.198\times 10^{-4}$ for their chi-square approximation and 0.08 for their non-studentized empirical bootstrap method.

Product functional principal component analysis represents the data $\textrm{PLV}_{i}(s,t)$ with products $\psi_{j}(s)\phi_{k}(t)$ of the marginal eigenfunctions, where $\psi_{j}(s)$ represents the frequency component and $\phi_{k}(t)$ represents the time component. The 3 estimated eigenfunction products that account for the most variance are $\hat{\psi}_{1}(s)\hat{\phi}_{1}(t)$, $\hat{\psi}_{2}(s)\hat{\phi}_{1}(t)$, and $\hat{\psi}_{3}(s)\hat{\phi}_{1}(t)$. The variance explained by $\hat{\psi}_{1}(s)\hat{\phi}_{1}(t)$ is by far the largest. The estimated marginal eigenvalues $\hat{\lambda}_{j}$ and $\hat{\gamma}_{k}$ are plotted in Figure 2. We see the first eigenvalue dominates the others, and there is also a slight drop between the second two $\hat{\lambda}_{j}$ and the rest, reflecting the fact that $\hat{\psi}_{2}(s)\hat{\phi}_{1}(t)$ and $\hat{\psi}_{3}(s)\hat{\phi}_{1}(t)$ are the products that explain the second and third highest amounts of variance, respectively.

The product functions $\hat{\psi}_{1}(s)\hat{\phi}_{1}(t)$, $\hat{\psi}_{2}(s)\hat{\phi}_{1}(t)$, and $\hat{\psi}_{3}(s)\hat{\phi}_{1}(t)$ are plotted in Figure 3. These products capture modes of variation mainly around -0.2 to 0.2 s, from when the subject receives the cue to move to when they just start moving. This variation can be seen more clearly in the first temporal eigenfunction $\hat{\phi}_{1}(t)$ (shown on the bottom right of Figure 4, which plots the individual marginal eigenfunctions), which peaks slightly after 0 s. $\hat{\psi}_{1}(s)\hat{\phi}_{1}(t)$ shows that, within this time range, subjects generally vary in synchrony from the alpha band to the beginning of the gamma low band, peaking within the beta band around 20–30 Hz. $\hat{\psi}_{2}(s)\hat{\phi}_{1}(t)$ shows a contrast between the beta low band and gamma low band. That is, subjects with higher $\chi_{21}$ values have lower synchrony in the beta low band and higher synchrony in the gamma low band. $\hat{\psi}_{3}(s)\hat{\phi}_{1}(t)$ shows a contrast between the alpha band and the beta high band.

6   Discussion

Much of the benefit of using product functional principal component analysis under weak separability is related to ease of interpretation and computation; by representing the eigenfunctions as tensor products of the marginal eigenfunctions, one consumes far fewer degrees of freedom and only needs to compute the marginal covariances instead of the full covariance. When the weak separable assumption does not hold, the product functional principal component analysis scores are correlated, and one expects to have to use more terms in product functional principal component analysis than in conventional functional principal component analysis to explain the same amount of variance. Product functional principal component analysis can still be used under the alternative as a dimension reduction approach, but one needs to be aware of the above issues. We believe that the notion of weak separability will inspire new methodological developments for multi-way functional data analysis, such as multi-way regularization on marginal components (Chen & Lei, 2015).

Although our data example has $s$ and $t$ in $\mathbb{R}$, our test works for scenarios where $d_{1}>1$ or $d_{2}>1$. For example, when modeling brain imaging data $X_{i}(s,t)$ observed on a dense grid of voxels or dipoles over the cortical surface, in which $s\in\mathbb{R}^{2}$ or $s\in\mathbb{R}^{3}$, one can first vectorize along $s$ by ordering the dipoles from 1 to $M$, compute the marginal covariances, perform the hypothesis test of weak separability, and reorganize back to the space $\mathcal{S}$ for interpretation and visualization. However, in scenarios where data are only very sparsely observed on the domain, where individual-subject smoothing is not appropriate, the problem is much more challenging and beyond the scope of this paper. We will possibly pursue it in future work.

Acknowledgement

Data were provided in part by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

7   Appendix: proofs

Proof of Lemma 1

Let $f_{j}\ (j=1,2,\ldots)$ and $g_{k}\ (k=1,2,\ldots)$ be a pair of bases that satisfies weak separability. For $(j,k)\neq(j^{\prime},k^{\prime})$, we have $\langle Cf_{j}\otimes g_{k},f_{j^{\prime}}\otimes g_{k^{\prime}}\rangle=E\left(\langle X-\mu,f_{j}\otimes g_{k}\rangle\langle X-\mu,f_{j^{\prime}}\otimes g_{k^{\prime}}\rangle\right)=0.$ Since the covariance operator $C$ is diagonalized under the orthonormal basis $f_{j}\otimes g_{k}\ (j=1,2,\ldots;\ k=1,2,\ldots)$, by Mercer’s theorem,

[TABLE]

where $\eta_{jk}=\langle Cf_{j}\otimes g_{k},f_{j}\otimes g_{k}\rangle={\rm var}\left(\langle X-\mu,f_{j}\otimes g_{k}\rangle\right)$, and the convergence is absolute and uniform.

The marginal kernel can then be written as

[TABLE]

The exchange of the integral and sums is allowed by the Fubini–Tonelli theorem, by noticing that

[TABLE]

where we use the Cauchy–Schwarz inequality.

Thus, we see that the $f_{j}$ are eigenfunctions of $C_{\mathcal{S}}$ with eigenvalues $\lambda_{j}=\sum_{k=1}^{\infty}\eta_{jk}$. An analogous computation shows that the $g_{k}$ are eigenfunctions of $C_{\mathcal{T}}$ with eigenvalues $\gamma_{k}=\sum_{j=1}^{\infty}\eta_{jk}$.

Proof of Lemma 2

With strong separability, we have $C(s,t;u,v)=a{C}_{1}(s,u){C}_{2}(t,v)$. From the definition of $C_{\mathcal{S}}$, we have

[TABLE]

An analogous argument shows ${C}_{\mathcal{T}}(t,v)=a{C}_{2}(t,v)$. Note that $a=\int_{\mathcal{T}}\int_{\mathcal{S}}C(s,t;s,t)dsdt$. If we use the marginal eigenfunctions $\psi_{j}$ and $\phi_{k}$ as the bases, it is easy to show that when $(j,k)\neq(j^{\prime},k^{\prime})$, ${\rm cov}(\chi_{jk},\chi_{j^{\prime}k^{\prime}})=\int_{\mathcal{T},\mathcal{S},\mathcal{T},\mathcal{S}}C(s,t;u,v)\psi_{j}(s)\phi_{k}(t)\psi_{j^{\prime}}(u)\phi_{k^{\prime}}(v)dsdtdudv=0$. Thus, we have weak separability.

Proof of Lemma 3

When $V$ is of rank 1, $V$ can be written $V=WZ^{T}$, where $W$ and $Z$ are column vectors with entries $(w_{1},w_{2},\ldots)$ and $(z_{1},z_{2},\ldots)$, respectively. Thus, $\eta_{jk}=w_{j}z_{k}$, and under weak separability, Equation 3 can be written

[TABLE]

The above can be normalized to fit the definition of strong separability in Lemma 2.

Under strong separability, from the proof of Lemma 2 we have

[TABLE]

so $\eta_{jk}=\{1/\int_{\mathcal{T}}\int_{\mathcal{S}}C(s,t;s,t)dsdt\}\lambda_{j}\gamma_{k}$, and then $V=\{1/\int_{\mathcal{T}}\int_{\mathcal{S}}C(s,t;s,t)dsdt\}\Lambda\Gamma^{T}$.

Proof of Theorem 4

For $H_{1}$ and $H_{2}$ two real separable Hilbert spaces, we further define the partial trace with respect to $H_{1}$ as the unique bounded linear operator $tr_{1}:\mathcal{B}_{Tr}(H_{1}\otimes H_{2})\rightarrow\mathcal{B}_{Tr}(H_{2})$ satisfying $tr_{1}(C_{1}\tilde{\otimes}C_{2})=tr(C_{1})C_{2}$ for all $C_{1}\in\mathcal{B}_{Tr}(H_{1})$, $C_{2}\in\mathcal{B}_{Tr}(H_{2})$. The partial trace with respect to $H_{2}$ is defined symmetrically and denoted by $tr_{2}$. With the notation of partial trace, we can see that $C_{\mathcal{T}}=tr_{1}(C)$ and $C_{\mathcal{S}}=tr_{2}(C)$. The estimated marginal covariance operators can also be written as $\hat{C}_{\mathcal{S}}=tr_{2}({C}_{n})$ and $\hat{C}_{\mathcal{T}}=tr_{1}({C}_{n})$. We use these equalities in proofs but not in computation. In practice, the estimated marginal covariances are calculated without having to calculate $C_{n}$.

We use similar notation and conditions as used by Aston et al. (2017). However, to derive the asymptotic distribution of their test statistic for strong separability, they focus on deriving the asymptotic distribution of the difference between the sample covariance operator and its strong separable approximation. Then by projecting on the estimated marginal eigenfunctions, they check the requirement for strong separability that $\eta_{jk}=a\lambda_{j}\gamma_{k}$. They do not need further results on the estimation errors of the marginal eigenfunctions and random scores besides that they are consistent. By contrast, our proofs involve the expansion of $\hat{\psi}_{j}-\psi_{j}$ and $\hat{\phi}_{k}-\phi_{k}$, and four-way tensor products with indices $(j,k,j^{\prime},k^{\prime})$.

From Condition 1 in Section 3.2 and the remark following it, $\mathcal{Z}_{n}=n^{1/2}(C_{n}-C)$ converges to a Gaussian random element in $\mathcal{B}_{Tr}\{L^{2}(\mathcal{S}\times\mathcal{T})\}$ with mean 0 and covariance structure $\Sigma_{C}={\rm E}[\{(X-\mu)\otimes(X-\mu)-C\}\tilde{\otimes}\{(X-\mu)\otimes(X-\mu)-C\}].$

For $T_{n}$ as defined in Equation 4,

[TABLE]

Using (5.1.8) in Hsing & Eubank (2015), we have

[TABLE]

where $\mathcal{M}_{j}=\sum_{m\neq j}(\lambda_{j}-\lambda_{m})^{-1}\psi_{m}\otimes\psi_{m}\in\mathcal{B}_{Tr}(\mathcal{S})$ and $\lambda_{j}$ is the $j$th eigenvalue of $C_{\mathcal{S}}$. Analogously,

[TABLE]

where $\mathcal{M^{\prime}}_{k}=\sum_{m\neq k}(\gamma_{k}-\gamma_{m})^{-1}\phi_{m}\otimes\phi_{m}\in\mathcal{B}_{Tr}(\mathcal{T})$ and $\gamma_{k}$ is the $k$th eigenvalue of $C_{\mathcal{T}}$. Here, Condition 2 is used to guarantee that $\mathcal{M}_{j}$ and $\mathcal{M^{\prime}}_{k}$ exist for $j=1,\ldots,P$ and $k=1,\ldots,K$.

Using $\hat{C}_{\mathcal{S}}-C_{\mathcal{S}}=tr_{2}(C_{n}-C)$ and $\hat{C}_{\mathcal{T}}-C_{\mathcal{T}}=tr_{1}(C_{n}-C)$, we can write $T_{n}(j,k,j^{\prime},k^{\prime})$ as

[TABLE]

The first term in the above equation is zero under $H_{0}$, since under $H_{0}$ we have the representation $C(s,t,u,v)=\sum_{j=1}^{\infty}\sum_{k=1}^{\infty}\eta_{jk}\psi_{j}(s)\psi_{j}(u)\phi_{k}(t)\phi_{k}(v),$ where $\eta_{jk}={\rm var}(\chi_{jk})$. Also, by Proposition C.1 in Aston et al. (2017), we have that $tr\{Atr_{1}(T)\}=tr\{(Id_{1}\tilde{\otimes}A)T\}$, where $Id_{1}$ is an identity operator on $\mathcal{S}$, $A\in\mathcal{B}(\mathcal{T})$, and $T\in\mathcal{B}_{Tr}(\mathcal{S}\times\mathcal{T})$. An analogous identity holds for $tr_{2}(T)$. Using these facts, we give a simplified form of $T_{n}(j,k,j^{\prime},k^{\prime})$ under $H_{0}$ for 3 cases:

(Case i) $j\neq j^{\prime}$ and $k\neq k^{\prime}$:

[TABLE]

(Case ii) $j=j^{\prime}$ and $k\neq k^{\prime}$:

[TABLE]

(Case iii) $j\neq j^{\prime}$ and $k=k^{\prime}$:

[TABLE]

In each of the above cases, two or more of the terms in Equation 5 end up being zero due to the orthogonality of the eigenfunctions. The latter 2 cases can be simplified to get the result in the statement of the theorem by noting that $\eta_{jk^{\prime}}(\phi_{k}\otimes\phi_{k^{\prime}})\mathcal{M^{\prime}}_{k}=\eta_{jk^{\prime}}(\gamma_{k}-\gamma_{k^{\prime}})^{-1}\phi_{k}\otimes\phi_{k^{\prime}}$ and $\eta_{jk}(\psi_{j^{\prime}}\otimes\psi_{j})\mathcal{M}_{j^{\prime}}=\eta_{jk}(\lambda_{j^{\prime}}-\lambda_{j})^{-1}\psi_{j^{\prime}}\otimes\psi_{j}$.

Proof of Corollary 5

From Theorem 4, we can see that all the terms of $T_{n}(j,k,j^{\prime},k^{\prime})$ can be written in the form $tr\{(A_{1}\tilde{\otimes}A_{2})\mathcal{Z}_{n}\}$ for some $A_{1}\in\mathcal{B}(\mathcal{S})$ and $A_{2}\in\mathcal{B}(\mathcal{T})$. Since $\mathcal{Z}_{n}$ converges to a Gaussian random element and $tr\{(A_{1}\tilde{\otimes}A_{2})\mathcal{Z}_{n}\}$ is a continuous linear mapping, the $T_{n}(j,k,j^{\prime},k^{\prime})$ are asymptotically jointly Gaussian for different sets of $(j,k,j^{\prime},k^{\prime})$. Let $\Theta$ be the covariance structure of the asymptotic joint distribution of the $T_{n}(j,k,j^{\prime},k^{\prime})$, and define $\mathcal{Z}$ to be a Gaussian random element with the limiting distribution of $\mathcal{Z}_{n}$. By the continuous mapping theorem, $\Theta$ can be calculated from terms of the form

[TABLE]

where $\Sigma_{C}$ is defined as in the proof of Theorem 4.

Recall the Karhunen–Loève expansion of the process

[TABLE]

We define $u_{ij}=\psi_{i}\otimes\psi_{j}\in\mathcal{B}_{HS}(\mathcal{S})$, $v_{ij}=\phi_{i}\otimes\phi_{j}\in\mathcal{B}_{HS}(\mathcal{T})$, $\beta_{i,i^{\prime},j,j^{\prime},k,k^{\prime},l,l^{\prime}}={\rm E}(\chi_{ii^{\prime}}\chi_{jj^{\prime}}\chi_{kk^{\prime}}\chi_{ll^{\prime}})$, and $\eta_{ii^{\prime}}={\rm E}(\chi_{ii^{\prime}}^{2})$. With weak separability, we have

[TABLE]

Each of the trace terms in the above equation can be evaluated using the identities $tr(Id_{1}u_{ij})=I(i=j)$, $tr(Id_{2}v_{i^{\prime}j^{\prime}})=I(i^{\prime}=j^{\prime})$, $tr\{(\psi_{j_{1}}\otimes\psi_{j_{1}^{\prime}})u_{ij}\}=I(i=j_{1})I(j=j_{1}^{\prime})$, and $tr\{(\phi_{k_{1}}\otimes\phi_{k_{1}^{\prime}})v_{i^{\prime}j^{\prime}}\}=I(i^{\prime}=k_{1})I(j^{\prime}=k_{1}^{\prime})$. From these identities and the possible forms of $A_{1}$, $A_{2}$, $B_{1}$, and $B_{2}$ given in Theorem 4, it follows that the second sum is always 0. The first sum can be simplified by considering 9 cases, as follows:

(Case 1) $A_{1}=a_{1}\psi_{j_{1}}\otimes\psi_{j_{1}^{\prime}}$, $A_{2}=a_{2}\phi_{k_{1}}\otimes\phi_{k_{1}^{\prime}}$, $B_{1}=b_{1}\psi_{j_{2}}\otimes\psi_{j_{2}^{\prime}}$, $B_{2}=b_{2}\phi_{k_{2}}\otimes\phi_{k_{2}^{\prime}}$:

[TABLE]

(Case 2) $A_{1}=Id_{1}$, $A_{2}=a_{2}\phi_{k_{1}}\otimes\phi_{k_{1}^{\prime}}$, $B_{1}=b_{1}\psi_{j_{2}}\otimes\psi_{j_{2}^{\prime}}$, $B_{2}=b_{2}\phi_{k_{2}}\otimes\phi_{k_{2}^{\prime}}$:

[TABLE]

(Case 3) $A_{1}=a_{1}\psi_{j_{1}}\otimes\psi_{j_{1}^{\prime}}$, $A_{2}=Id_{2}$, $B_{1}=b_{1}\psi_{j_{2}}\otimes\psi_{j_{2}^{\prime}}$, $B_{2}=b_{2}\phi_{k_{2}}\otimes\phi_{k_{2}^{\prime}}$:

[TABLE]

(Case 4) $A_{1}=a_{1}\psi_{j_{1}}\otimes\psi_{j_{1}^{\prime}}$, $A_{2}=a_{2}\phi_{k_{1}}\otimes\phi_{k_{1}^{\prime}}$, $B_{1}=Id_{1}$, $B_{2}=b_{2}\phi_{k_{2}}\otimes\phi_{k_{2}^{\prime}}$:

[TABLE]

(Case 5) $A_{1}=a_{1}\psi_{j_{1}}\otimes\psi_{j_{1}^{\prime}}$, $A_{2}=a_{2}\phi_{k_{1}}\otimes\phi_{k_{1}^{\prime}}$, $B_{1}=b_{1}\psi_{j_{2}}\otimes\psi_{j_{2}^{\prime}}$, $B_{2}=Id_{2}$:

[TABLE]

(Case 6) $A_{1}=Id_{1}$, $A_{2}=a_{2}\phi_{k_{1}}\otimes\phi_{k_{1}^{\prime}}$, $B_{1}=Id_{1}$, $B_{2}=b_{2}\phi_{k_{2}}\otimes\phi_{k_{2}^{\prime}}$:

[TABLE]

(Case 7) $A_{1}=Id_{1}$, $A_{2}=a_{2}\phi_{k_{1}}\otimes\phi_{k_{1}^{\prime}}$, $B_{1}=b_{1}\psi_{j_{2}}\otimes\psi_{j_{2}^{\prime}}$, $B_{2}=Id_{2}$:

[TABLE]

(Case 8) $A_{1}=a_{1}\psi_{j_{1}}\otimes\psi_{j_{1}^{\prime}}$, $A_{2}=Id_{2}$, $B_{1}=Id_{1}$, $B_{2}=b_{2}\phi_{k_{2}}\otimes\phi_{k_{2}^{\prime}}$:

[TABLE]

(Case 9) $A_{1}=a_{1}\psi_{j_{1}}\otimes\psi_{j_{1}^{\prime}}$, $A_{2}=Id_{2}$, $B_{1}=b_{1}\psi_{j_{2}}\otimes\psi_{j_{2}^{\prime}}$, $B_{2}=Id_{2}$:

[TABLE]

In the above, $a_{1}$, $a_{2}$, $b_{1}$, and $b_{2}$ are scalar constants. Using the above, all the terms in $\Theta$ can be obtained from straightforward but tedious calculations.

To illustrate the calculation of $\Theta(j,k,j^{\prime},k^{\prime},l,m,l^{\prime},m^{\prime})$, the term in $\Theta$ corresponding to the asymptotic covariance of $T_{n}(j,k,j^{\prime},k^{\prime})$ and $T_{n}(l,m,l^{\prime},m^{\prime})$, we consider as an example the case where $j\neq j^{\prime}$, $k\neq k^{\prime}$, $l\neq l^{\prime}$, and $m\neq m^{\prime}$. Here,

[TABLE]

where we have used $A_{1}=\psi_{j}\otimes\psi_{j^{\prime}}$, $A_{2}=\phi_{k}\otimes\phi_{k^{\prime}}$, $B_{1}=\psi_{l}\otimes\psi_{l^{\prime}}$, and $B_{2}=\phi_{m}\otimes\phi_{m^{\prime}}$.

Proof of Lemma 6

Let $X_{N}(s,t)=\mu(s,t)+\sum_{j=1}^{N}\sum_{k=1}^{N}\chi_{jk}\psi_{j}(s)\phi_{k}(t)$, and let $C_{N}$ denote the covariance structure of $X_{N}$. Thus,

[TABLE]

It is easy to show that $C_{N}$ converges to $C$ in Hilbert–Schmidt norm. Let $C_{\mathcal{S},N}=tr_{2}(C_{N})$, which converges to $C_{\mathcal{S}}$ because $tr_{2}$ is continuous and linear. We know that $\langle C_{\mathcal{S}}\psi_{j},\psi_{j^{\prime}}\rangle=0$ for $j\neq j^{\prime}$. Therefore, for any $\epsilon>0$, we can find an $N$ such that $|\langle C_{\mathcal{S},N}\psi_{j},\psi_{j^{\prime}}\rangle|<\epsilon$.

By definition,

[TABLE]

Therefore, $\lim_{N}\sum_{k=1}^{N}{\rm cov}(\chi_{jk},\chi_{j^{\prime}k})=0$, i.e., $\sum_{k=1}^{\infty}{\rm cov}(\chi_{jk},\chi_{j^{\prime}k})=0$ for $j\neq j^{\prime}$.

The same argument holds for the empirical version. Analogous calculations can be done for $k\neq k^{\prime}$ to show that $\sum_{j=1}^{\infty}{\rm cov}({\chi}_{jk},{\chi}_{jk^{\prime}})=0$ and $\sum_{j=1}^{\infty}T_{n}(j,k,j,k^{\prime})=0$.