# CovNet: Covariance Networks for Functional Data
on Multidimensional Domains

Soham<EMAIL_ADDRESS>[ Victor M.
<EMAIL_ADDRESS>[ Ecole Polytechnique Fédérale de
Lausanne  Institut de Mathématiques
Ecole Polytechnique Fédérale de Lausanne
e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Covariance estimation is ubiquitous in functional data analysis. Yet, the case
of functional observations over multidimensional domains introduces
computational and statistical challenges, rendering the standard methods
effectively inapplicable. To address this problem, we introduce _Covariance
Networks_ (CovNet) as a modeling and estimation tool. The CovNet model is
_universal_ – it can be used to approximate any covariance up to desired
precision. Moreover, the model can be fitted efficiently to the data and its
neural network architecture allows us to employ modern computational tools in
the implementation. The CovNet model also admits a closed-form
eigendecomposition, which can be computed efficiently, without constructing
the covariance itself. This facilitates easy storage and subsequent
manipulation of a covariance in the context of the CovNet. We establish
consistency of the proposed estimator and derive its rate of convergence. The
usefulness of the proposed method is demonstrated by means of an extensive
simulation study and an application to resting state fMRI data.

62G05, 62M40, 62M45,

15A99, 68T07,

deep learning,

FDA,

neural network,

nonparametric model,

universal approximation,

###### keywords:

[class=AMS]

###### keywords:

and

t1Research supported by a Swiss National Science Foundation grant.

###### Contents

  1. 1 Introduction
  2. 2 Covariance networks
    1. 2.1 Shallow architecture
    2. 2.2 Deep architectures
  3. 3 Practical implementation
  4. 4 Eigendecomposition of the estimated covariance operator
  5. 5 Empirical study
    1. 5.1 Choice of hyperparameters
    2. 5.2 Estimated eigenstructure
    3. 5.3 Application to fMRI data
  6. 6 Asymptotic theory
    1. 6.1 The case of fully observed fields
    2. 6.2 The case of discretely observed fields
  7. 7 Concluding remarks
  8. A Mathematical background
  9. B Bias of the CovNet model: universal approximation and rate of convergence
    1. B.1 Universal approximation
    2. B.2 Rate of convergence of the bias term
  10. C Further details on the implementation of the CovNet models
    1. C.1 Detailed derivation of the loss function
    2. C.2 Simultaneous estimation of the mean
    3. C.3 Computing cost
  11. D Covering numbers
    1. D.1 Covering number of the shallow CovNet class
    2. D.2 Covering number of the deep CovNet class
    3. D.3 Covering number of the deepshared CovNet class
  12. E Proofs of the asymptotic results
    1. E.1 Consistency
    2. E.2 Rates of convergence
    3. E.3 Consistency without boundedness
    4. E.4 The case of discretely observed data
  13. F Additional simulation results

##  1 Introduction

We consider the problem of covariance estimation from a collection of
functional observations defined over a multidimensional domain. To be precise,
let $\mathcal{X}=\\{X(\text{\bf u}):\text{\bf u}\in\mathcal{Q}\\}$ be a
compactly supported random field, i.e., a real-valued second-order stochastic
process on a compact set $\mathcal{Q}\subset\mathbb{R}^{d}$, with covariance
kernel $c(\text{\bf u},\text{\bf v})=\mathrm{Cov}(X(\text{\bf u}),X(\text{\bf
v}))$. We want to estimate $c$ based on an independent and identically
distributed (i.i.d.) sample
$\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\sim\mathcal{X}$. In particular, we
work in the framework of _functional data analysis_ (FDA, see Ramsay and
Silverman, 2002; Hsing and Eubank, 2015), where we assume that $\mathcal{X}$
takes values in $\mathcal{L}_{2}(\mathcal{Q})$, the space of all real-valued
square-integrable functions on $\mathcal{Q}$.

Covariance estimation, along with mean estimation, is a fundamental problem in
functional data analysis and has multifaceted applications, e.g., in
regression, prediction, classification. This problem has been studied
extensively for observations over one-dimensional domains (i.e., $d=1$) or
_curve data_ (see Wang, Chiou and Müller (2016) for a detailed review). The
same is, however, not true for observations over multidimensional domains.
Although, in principal, these two regimes are similar, and most methods for
curve data can be “readily used” for data over multidimensional domains, in
practice, the dimensionality of the problem draws a clear distinction between
the two paradigms. To appreciate this, suppose that we observe the random
fields on a grid of size $K\times\cdots\times K$ in
$\mathcal{Q}\subset\mathbb{R}^{d}$. In this case, the estimation of the
empirical covariance requires $\mathcal{O}(K^{2d})$ computations. The storage
cost for this estimator is also of the order $\mathcal{O}(K^{2d})$ which, for
$d=2$, can be prohibitive even for $K\approx 100$. The problem becomes even
more severe when $d$ is larger ($d\geq 3$), which is increasingly common,
e.g., for observations over spatial volumes or of a spatio-temporal nature.
Moreover, subsequent manipulation, e.g., inversion, needed in applications,
requires computation in the order of $\mathcal{O}(K^{3d})$, leading to a
prohibitive computational burden.

To put things into perspective, consider the $1000$ Functional Connectomes
Project111https://www.nitrc.org/projects/fcon_1000/ which contains functional
magnetic resonance imaging (fMRI) of brains for more than $1200$ individuals.
For each individual, the data consist of 3D brain-scans on a grid of size
$64\times 64\times 33$ taken at $2$ second intervals over $225$ time points.
Covariance estimation is of utmost importance in fMRI studies as it captures
the connectivity patterns in the brain (Aston and Kirch, 2012; Stoehr, Aston
and Kirch, 2021). At the same time, it is extremely difficult to do so
nonparametrically because of the dimensionality of the problem. For instance,
the empirical covariance estimator for the 3D fMRI data would be an object of
size $64\times 64\times 33\times 64\times 64\times 33$ which requires $68$ Gb
memory to compute and store (at $32$-byte precision). This is impossible with
a regular computer which usually has $16$ or $32$ Gb of memory. Also, apart
from looking at the connectivity pattern of the brain as a 3D object, it is
also of importance to check how these patterns evolve over time, i.e., to
consider the full 4D data on a grid of size $64\times 64\times 33\times 225$.
The problem becomes even more severe in this case, where the empirical
covariance would require approximately $3.4\times 10^{6}$ Gb of memory during
computation and for storage (see also Aston and Kirch, 2012; Stoehr, Aston and
Kirch, 2021).

To alleviate this _curse of dimensionality_ , further modeling assumptions are
often made on the underlying covariance, the most popular being that of
_separability_. A _separable model_ assumes that the true covariance over the
multidimensional domain can be factored into several covariances over one-
dimensional domains, i.e., $c(\text{\bf u},\text{\bf
v})=c_{1}(u_{1},v_{1})\times\cdots\times c_{d}(u_{d},v_{d})$ for $\text{\bf
u},\text{\bf v}\in\mathcal{Q}$. This greatly simplifies the problem and
entails enormous computational savings. For instance, in the case of
observations on a grid, a separable model can be estimated with
$\mathcal{O}(dK^{2})$ computations and has the same order of storage
requirements. The gain during application of the model is even better – the
inversion of the model requires $\mathcal{O}(dK^{3})$ computations compared to
the $\mathcal{O}(K^{3d})$ for the empirical covariance. Despite all these
advantages, separability is merely a modeling assumption, which is highly
restrictive and often violated in practice (Aston, Pigoli and Tavakoli, 2017;
Constantinou, Kokoszka and Reimherr, 2017; Rougier, 2017; Bagchi and Dette,
2020). Still, it is often the preferred choice in practice, not because it is
believed to hold, but rather for the savings that it entails (Gneiting, Genton
and Guttorp, 2006; Pigoli et al., 2018). Perhaps it is safe to say that the
popularity of the separable model stems from the non-availability of a better
alternative. It is worth clarifying here that when the data are _sparse_ (each
random field is observed at a few randomly scattered locations) there do exist
methods applicable to multidimensional domains without assuming separability
(for example, Wang, Wong and Zhang (2020) used a penalized method, leading to
the use of tensor-products of splines, and yielding a low-rank approximation
to the covariance). However, such approaches are infeasible in the _dense_
regime, where each random field is measured on the same dense grid (e.g., the
fMRI data), and which is our main interest in this paper. It is the denseness
of the measurements that gives rise to the severe computational challenges
mentioned above (e.g., in a dense regime, the approach of Wang, Wong and Zhang
(2020) is infeasible since it requires computation of the “raw covariances” –
equivalent to the computation of the empirical covariance).

In an effort to deliver a more general yet tractable approach, we propose a
new model for covariance estimation using neural networks. Neural networks
have long been successfully used in nonparametric function estimation, and
recently, they have been shown to successfully overcome the curse of
dimensionality in nonparametric regression (Bauer and Kohler, 2019; Schmidt-
Hieber, 2020). Also, they have been used for mean estimation of functional
data over multidimensional domains (Wang, Cao and Shang, 2021). Motivated by
the success of neural networks, we propose _Covariance Networks_ (CovNet) as a
framework for the estimation of the covariance of multidimensional random
fields. A CovNet is a positive semi-definite function on
$\mathcal{Q}\times\mathcal{Q}$ described by a neural network architecture. In
particular, we define and study three variants: the _shallow_ CovNet model and
the _deep_ CovNet model which differ with respect to the _depth_ of the
network; and the _deepshared_ CovNet model which is a restricted (regularised)
version of the deep CovNet model.

Our framework features several advantages, namely:

  1. 1.

It is genuinely _nonparametric_ – any covariance can be approximated up to
arbitrary precision via a CovNet structure. We establish this so-called
_universal approximation property_ of the CovNet models in Theorems 1 and 2.
Moreover, the proposed model has an explicit functional form. This functional
form has its own advantage in applications such as kriging, where we do not
need to interpolate or smooth the estimated covariance before use.

  2. 2.

Fitting a CovNet to the data is computationally tractable. The models we
introduce can be fitted at the level of the data, without the need to compute
or store any high-order objects. Moreover, the neural network structure allows
us to exploit modern machine learning tools during the estimation. These are
discussed in Section 3.

  3. 3.

The special structure of the CovNet models ensures that the eigendecomposition
of the associated operator can be obtained without the need to explicitly form
the operator itself. Thus, we can access the eigensystem of the fitted CovNet
very easily without ever forming any higher order objects. This allows us to
store and subsequently manipulate the fitted model very easily (see Section
4).

  4. 4.

The CovNet estimators come with theoretical guarantees. In particular, we
establish their consistency and derive their rates of convergence (Section 6).
Our analyses are _fully nonparametric_ – we make no structural assumption on
the underlying covariance $\mathcal{C}$ for our derivations.

The rest of the article is organized as follows. We lay out our methodology in
the next section. In particular, we begin by describing the shallow CovNet
model in Section 2.1 and establish its _universality_. In Section 2.2, we
extend the shallow CovNet model by using _deep_ architectures, leading to the
_deep_ and the _deepshared_ CovNet models. In Section 3, we demonstrate how
the CovNet models can be efficiently estimated in practice. The
eigendecomposition of the CovNet operator is discussed in Section 4. In
Section 5, we demonstrate the usefulness of the proposed CovNet models by
means of a detailed simulation study and an application to the fMRI data. We
establish the theoretical properties of these models in Section 6. Some
concluding remarks are made in Section 7. The proofs of our asymptotic results
are provided in the appendices. The appendices also cover some related
mathematical ideas, as well as some further numerical results.

##  2 Covariance networks

We start with some background concepts, more details can be found in Hsing and
Eubank (2015) and Appendix A. Let $\mathcal{Q}$ be a compact subset of
$\mathbb{R}^{d}$ and let $\mathcal{X}=\\{X(\text{\bf u}):\text{\bf
u}\in\mathcal{Q}\\}$ be a random element of $\mathcal{L}_{2}(\mathcal{Q})$.
For $d=1$, $\mathcal{X}$ is usually referred to as a _random curve_ , whereas
for $d>1$, it is referred to as a _random field_. We assume that $\mathcal{X}$
has finite second moment, i.e., $\mathbb{E}(\|\mathcal{X}\|^{2})<\infty$,
which ensures the existence of its mean $m=\mathbb{E}(\mathcal{X})$ and
covariance $\mathcal{C}=\mathbb{E}\\{(\mathcal{X}-m)\otimes(\mathcal{X}-m)\\}$
(both the expectations are in the Bochner sense, see Hsing and Eubank 2015).
Here, $\|\cdot\|$ is the $\mathcal{L}_{2}$-norm associated with the inner-
product $\langle f,g\rangle=\int_{\mathcal{Q}}f(\text{\bf u})g(\text{\bf
u})\,\mathrm{d}\text{\bf u}$ for $f,g\in\mathcal{L}_{2}(\mathcal{Q})$, and the
tensor product $h\otimes g$ denotes the rank 1 operator $f\mapsto\langle
f,g\rangle h$. The covariance $\mathcal{C}$ is a linear operator from
$\mathcal{L}_{2}(\mathcal{Q})$ onto itself, given by

$\mathcal{C}f(\text{\bf u})=\int_{\mathcal{Q}}c(\text{\bf u},\text{\bf
v})\,f(\text{\bf v})\,\mathrm{d}\text{\bf
v},~{}~{}f\in\mathcal{L}_{2}(\mathcal{Q}).$

Here, $c\in\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})$ defined as
$c(\text{\bf u},\text{\bf v})=\mathrm{Cov}(X(\text{\bf u}),X(\text{\bf v}))$
is the _covariance kernel_ associated with $\mathcal{X}$. We also say that
$\mathcal{C}$ is the _integral operator_ associated with the kernel $c$. The
operator $\mathcal{C}$ is positive semi-definite and the kernel $c$ is non-
negative definite. The Hilbert-Schmidt norm
${|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|}_{2}$
of $\mathcal{C}$ is finite, and
${|\kern-1.07639pt|\kern-1.07639pt|\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}=\|c\|_{\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})}$.
Thus, the covariance operator $\mathcal{C}$ and the covariance kernel $c$ are
linked by an isometric isomorphism. Consequently, we can use $\mathcal{C}$ and
$c$ interchangeably, and the estimation of the covariance $\mathcal{C}$ is
equivalent to the estimation of the kernel $c$. Since the object of interest
is the covariance rather than the mean, we work under the assumption that
$m=\mathbb{E}(\mathcal{X})=0$, unless specifically mentioned.

###  2.1 Shallow architecture

We propose to estimate the covariance kernel $c$ using the following neural
network structure:

$c_{\rm sh}(\text{\bf u},\text{\bf
v})=\sum_{r=1}^{R}\sum_{s=1}^{R}\lambda_{r,s}\,\sigma(\mathbf{w}_{r}^{\top}\text{\bf
u}+b_{r})\,\sigma(\mathbf{w}_{s}^{\top}\text{\bf v}+b_{s}),\qquad\text{\bf
u},\text{\bf v}\in\mathcal{Q},$ (2.1)

where $R\in\mathbb{N}$ is the _width_ , $\sigma:\mathbb{R}\to\mathbb{R}$ is an
activation function, and $\Lambda:=((\lambda_{r,s}))$ is a positive semi-
definite matrix. The parameters $\mathbf{w}_{r}\in\mathbb{R}^{d}$ and
$b_{r}\in\mathbb{R}$ for $r=1,\ldots,R$ are the _weights_ and the _biases_ of
the model (2.1). Positive semi-definiteness of $\Lambda$ readily implies that
$c_{\rm sh}$ is a non-negative definite kernel. For a given activation
function $\sigma$ and width $R\in\mathbb{N}$, we define

$\mathscr{F}^{\rm sh}_{R,\sigma}=\left\\{c_{\rm sh}\text{ of the form
}\eqref{eq:shallow_covnet_kernel}:\Lambda=((\lambda_{r,s}))\succeq
0,\mathbf{w}_{1},\ldots,\mathbf{w}_{R}\in\mathbb{R}^{d},b_{1},\ldots,b_{R}\in\mathbb{R}\right\\},$
(2.2)

to be the class of shallow Covariance Network kernels or _shallow CovNet_
kernels. We also define

$\widetilde{\mathscr{F}}^{\rm
sh}_{R,\sigma}=\left\\{\mathcal{G}:\mathcal{G}\text{ is an integral operator
with kernel }g\in\mathscr{F}^{\rm sh}_{R,\sigma}\right\\},$ (2.3)

to be the class of shallow covariance network operators or shallow CovNet
operators.

We call the structure (2.1) _shallow_ because each of the constituents
$\sigma(\mathbf{w}_{r}^{\top}\cdot+\,b_{r})$ for $r=1,\ldots,R$ of the kernel
(2.1) is a shallow neural network, i.e., a neural network with a single hidden
layer. The special structure of the kernel (2.1) allows us to visualize it as
a neural network with two hidden layers, as depicted in Figure 1. In the first
layer, starting from the inputs u and v, single-layer perceptrons
$\sigma(\mathbf{w}_{r}^{\top}\text{\bf u}+b_{r})$ and
$\sigma(\mathbf{w}_{r}^{\top}\text{\bf v}+b_{r}),r=1,\ldots,R$ are computed.
In the next layer, these outputs are cross-multiplied with the weights
$\lambda_{r,s}$ to produce the final result $c_{\rm sh}(\text{\bf u},\text{\bf
v})$. As one can see, this is a feed-forward neural network (see Anthony and
Bartlett, 1999, Chapter 6), which is completely determined (for fixed $\sigma$
and $R$) by the parameters $\mathbf{w}_{1},\ldots,\mathbf{w}_{R}$,
$b_{1},\ldots,b_{R}$ and $\Lambda=((\lambda_{r,s}))$ (with the added
restriction on $\Lambda$).

Figure 1: A schematic representation of the shallow CovNet structure. The top
layer corresponds to the inputs u and v. In the first hidden layer (second
from the top), the inputs are projected by weights, shifted by bias and
transformed by the activation function $\sigma$ to produce single-layer
perceptrons $\sigma(\mathbf{w}_{r}^{\top}\text{\bf
u}+b_{r}),\sigma(\mathbf{w}_{r}^{\top}\text{\bf v}+b_{r}),r=1,\ldots,R$. The
weights and biases of the two sides (left and right), corresponding to u and
v, are the same. In the second hidden layer (third from the top), the outputs
of the first hidden layer are cross-multiplied. These are then multiplied by
weights $\lambda_{r,s}$ and added to produce the output $c_{\rm sh}(\text{\bf
u},\text{\bf v})$.

As mentioned in the introduction, the shallow CovNet structure (2.1) is a
general model, in the sense that any covariance kernel can be approximated
with arbitrary precision using a shallow CovNet kernel of the form (2.1).
Thus, we do not need to make any assumption on the underlying covariance $c$,
resulting in a completely nonparametric procedure. However, we do need a
particular condition on the activation function $\sigma$ of the network.

######  Definition 1 (Sigmoidal activation).

An activation function $\sigma:\mathbb{R}\to[0,1]$ is said to be sigmoidal if
it is non-decreasing with $\lim_{x\to\infty}\sigma(x)=1\text{ and
}\lim_{x\to-\infty}\sigma(x)=0$.

In probabilistic terms, a sigmoidal function is a cumulative distribution
function. Sigmoidal activations are very common in the literature of neural
networks. One of the most popular activation functions, the _logistic
function_ $\sigma(t)=1/(1+\exp(-t))$ is a sigmoidal function often also
referred to as the _sigmoid_. Many other popularly used activation functions
are also sigmoidal (see Györfi et al., 2002, Chapter 16, for a plethora of
examples). It is worthwhile to note that the definition of sigmoidal functions
is not universal. In this article, whenever we refer to a sigmoidal function,
we mean it in the sense of Definition 1.

If we use a sigmoidal activation function, then any covariance kernel can be
approximated up to arbitrary precision using a shallow CovNet kernel of the
form (2.1). Such a property is often referred to as the _universal
approximation property_ in the computer science literature.

######  Theorem 1 (Shallow CovNet is Universal Approximator).

Let $c:\mathcal{Q}\times\mathcal{Q}\to\mathbb{R}$ be the kernel of the
covariance operator $\mathcal{C}$. Also, assume that the activation function
$\sigma$ is sigmoidal. Then, for every $\epsilon>0$, there exists
$R\in\mathbb{N}$ and $c_{\rm sh}\in\mathscr{F}^{\rm sh}_{R,\sigma}$ such that

$\int_{\mathcal{Q}\times\mathcal{Q}}\big{|}c(\text{\bf u},\text{\bf v})-c_{\rm
sh}(\text{\bf u},\text{\bf v})\big{|}^{2}\,\mathrm{d}\text{\bf
u}\,\mathrm{d}\text{\bf v}\leq\epsilon.$

If in addition $c$ is continuous, then the same conclusion holds uniformly.
That is, for every $\epsilon>0$, we can find $R\in\mathbb{N}$ and $c_{\rm
sh}\in\mathscr{F}^{\rm sh}_{R,\sigma}$ such that

$\sup_{\text{\bf u},\text{\bf v}\in\mathcal{Q}}\big{|}c(\text{\bf u},\text{\bf
v})-c_{\rm sh}(\text{\bf u},\text{\bf v})\big{|}\leq\epsilon.$

######  Remark 1.

The proof of the theorem rests on the universal approximation property of
single hidden layer neural networks on the class of square integrable
functions on $\mathcal{Q}$. The sigmoidal condition on the activation function
ensures this property, but is not necessary (see, e.g., Pinkus, 1999).

The construction of the shallow CovNet in Figure 1 shows an immediate way to
make the structure _deep_ by augmenting more layers. It is well-known that
shallow networks may require a rather large width to approximate a function,
whereas the same precision can be achieved by using a deep and less wide
network (Eldan and Shamir, 2016; Liang and Srikant, 2017; Poggio et al.,
2017). Deep networks can capture more complex structures than the shallow
networks with much fewer parameters. Moreover, during training, shallow
networks are more prone to get stuck at bad local minima, which are usually
avoided by deep networks (Choromanska et al., 2015). In the next section, we
extend the CovNet structure by incorporating deep networks instead of
perceptrons in the construction.

###  2.2 Deep architectures

We start with a brief description of deep neural networks. For an integer
$L>1$, an integer-tuple $\mathbf{p}=(p_{1},\ldots,p_{L})$, matrices
$\mathrm{W}_{1}\in\mathbb{R}^{p_{1}\times d}$,
$\mathrm{W}_{2}\in\mathbb{R}^{p_{2}\times
p_{1}},\ldots,\mathrm{W}_{L}\in\mathbb{R}^{p_{L}\times p_{L-1}}$, vectors
$\mathbf{b}_{1}\in\mathbb{R}^{p_{1}},\ldots,\mathbf{b}_{L}\in\mathbb{R}^{p_{L}},\mathbf{w}_{L+1}\in\mathbb{R}^{p_{L}}$,
and $b_{L+1}\in\mathbb{R}$, define the function
$g:\mathbb{R}^{d}\to\mathbb{R}$ which maps $\text{\bf u}\mapsto g(\text{\bf
u})$ recursively as follows:

$\displaystyle\text{\bf u}_{1}$ $\displaystyle=\sigma(\mathrm{W}_{1}\text{\bf
u}+\mathbf{b}_{1})$ $\displaystyle\text{\bf u}_{l+1}$
$\displaystyle=\sigma(\mathrm{W}_{l+1}\text{\bf
u}_{l}+\mathbf{b}_{l+1})\quad\text{ for }l=1,\ldots,L-1,$ $\displaystyle
g(\text{\bf u})$
$\displaystyle=\sigma(\mathbf{w}_{L+1}^{\top}\mathbf{u}_{L}+b_{L+1}).$ (2.4)

Here, for a vector $\mathbf{z}\in\mathbb{R}^{p}$, $\sigma(\mathbf{z})$
represents the component-wise application of the function $\sigma$. The
function $g$ is a deep neural network, where $L$ is the number of hidden
layers or _depth_ , $p_{1},\ldots,p_{L}$ are the _widths_ of the hidden layers
($p_{\rm max}=\max\\{p_{1},\ldots,p_{L}\\}$ is sometimes referred to as the
width of the network) and
$\mathrm{W}_{1},\ldots,\mathrm{W}_{L},\mathbf{w}_{L+1},\mathbf{b}_{1},\ldots,\mathbf{b}_{L},b_{L+1}$
are the network parameters. A schematic representation of the deep neural
network is shown in Figure 2(a). Starting from the input u, we go to the first
hidden layer by multiplying it with the weight matrix $\mathrm{W}_{1}$, adding
the bias $\mathbf{b}_{1}$ and applying the activation function $\sigma$
component-wise on the resultant. The same structure is repeated for all the
subsequent layers. We define the class

$\displaystyle\mathcal{D}_{L,\mathbf{p}}=\big{\\{}g:\mathbb{R}^{d}\to\mathbb{R}\text{
of the form \eqref{eq:deep_neural_network} with }$
$\displaystyle\mathrm{W}_{1}\in\mathbb{R}^{p_{1}\times
d},\mathbf{b}_{1}\in\mathbb{R}^{p_{1}},\mathrm{W}_{2}\in\mathbb{R}^{p_{2}\times
p_{1}},\mathbf{b}_{2}\in\mathbb{R}^{p_{2}},\ldots,$ $\displaystyle\kern
8.61108pt\mathrm{W}_{L}\in\mathbb{R}^{p_{L}\times
p_{L-1}},\mathbf{b}_{L}\in\mathbb{R}^{p_{L}},\mathbf{w}_{L+1}\in\mathbb{R}^{p_{L}},b_{L+1}\in\mathbb{R}\big{\\}},$
(2.5)

to be the class of all possible deep neural networks with depth $L$ and widths
$p_{1},\ldots,p_{L}$. A network from the class $\mathcal{D}_{L,\mathbf{p}}$
has $\sum_{l=0}^{L}(p_{l}+1)p_{l+1}$ parameters, where $p_{0}=d$ and
$p_{L+1}=1$.

We define the deep CovNet kernel as

$c_{\rm d}(\text{\bf u},\text{\bf
v})=\sum_{r=1}^{R}\sum_{s=1}^{R}\lambda_{r,s}\,g_{r}(\text{\bf
u})\,g_{s}(\text{\bf v}),\qquad\text{\bf u},\text{\bf v}\in\mathcal{Q},$ (2.6)

where $g_{1},\ldots,g_{R}\in\mathcal{D}_{L,\mathbf{p}}$ and
$\Lambda:=((\lambda_{r,s}))$ is positive semi-definite. This is similar to the
shallow CovNet kernel (2.1), except the constituents $g_{r}(\text{\bf u})$ are
deep networks of the form (2.2) instead of the perceptrons
$\sigma(\mathbf{w}_{r}^{\top}\text{\bf u}+b_{r})$. A schematic representation
of the deep CovNet kernel is shown in Figure 2(b). We define the class of deep
CovNet kernels and the corresponding class of operators as

$\displaystyle\mathscr{F}^{\rm d}_{R,L,\mathbf{p},\sigma}$
$\displaystyle=\big{\\{}c_{\rm d}\text{ of the form
\eqref{eq:deep_covnet_kernel}}:\Lambda=((\lambda_{r,s}))\succeq\mathrm{0},g_{1},\ldots,g_{R}\in\mathcal{D}_{L,\mathbf{p}}\big{\\}}$
$\displaystyle\widetilde{\mathscr{F}}^{\rm d}_{R,L,\mathbf{p},\sigma}$
$\displaystyle=\big{\\{}\mathcal{G}:\mathcal{G}\text{ is the integral operator
associated with kernel }g\in\mathscr{F}^{\rm
d}_{R,L,\mathbf{p},\sigma}\big{\\}}.$ (2.7)

######  Remark 2.

As can be seen from the construction of the deep CovNet structure, it is
possible to allow the individual networks $g_{1},\ldots,g_{R}$ to have
different depths and widths, allowing for more flexible models. However, this
complicates the analysis, so we do not pursue this model in this paper.

(a) Schematic representation of a deep neural network
---
(b) Schematic representation of a deep CovNet structure
Figure 2: (a) A schematic representation of the deep neural network. (b) A
schematic representation of the deep CovNet structure. Starting from the
inputs u and v, $R$ individual deep neural networks are fitted to get the
outputs $g_{1}(\text{\bf u}),\ldots,g_{R}(\text{\bf u})$ and $g_{1}(\text{\bf
v}),\ldots,g_{R}(\text{\bf v})$. The networks on left and right (corr. to u
and v) are the same. Cross-products of the outputs of the individual networks
layer are taken in the next layer, which are then multiplied by weights
$\lambda_{r,s}$ and added to produce the output $c_{\rm d}(\text{\bf
u},\text{\bf v})$.

The deep CovNet model is quite rich. Moreover, it retains the universal
approximation property of the shallow CovNet model (Theorem 2). But the number
of parameters of the deep CovNet model can be quite large, making it prone to
overfitting. Thus, some sort of regularization is needed for the deep CovNet
structure to make it more stable. We do this by enforcing _weight sharing_
among the constituents as follows. For an integer $L>1$, integer-tuple
$\mathbf{p}=(p_{1},\ldots,p_{L})$, matrices
$\mathrm{W}_{1}\in\mathbb{R}^{p_{1}\times d}$,
$\mathrm{W}_{2}\in\mathbb{R}^{p_{2}\times
p_{1}},\ldots,\mathrm{W}_{L}\in\mathbb{R}^{p_{L}\times p_{L-1}}$, and vectors
$\mathbf{b}_{1}\in\mathbb{R}^{p_{1}},\ldots,\mathbf{b}_{L}\in\mathbb{R}^{p_{L}}$,
we define the networks $g_{1},\ldots,g_{R}$ jointly as

$\displaystyle\text{\bf u}_{1}$ $\displaystyle=\sigma(W_{1}\text{\bf
u}+\mathbf{b}_{1}),$ $\displaystyle\text{\bf u}_{l+1}$
$\displaystyle=\sigma(W_{l+1}\text{\bf u}_{l}+\mathbf{b}_{l+1})\quad\text{ for
}l=1,\ldots,L-1,$ $\displaystyle g_{r}(\text{\bf u})$
$\displaystyle=\sigma(\bm{\omega}_{r}^{\top}\text{\bf u}_{L}+\beta_{r}),\quad
r=1,\ldots,R,$ (2.8)

where $\bm{\omega}_{r}\in\mathbb{R}^{p_{L}}$ and $\beta_{r}\in\mathbb{R}$ for
$r=1,\ldots,R$. Individually, each of the networks $g_{1},\ldots,g_{R}$ is an
element of the class $\mathcal{D}_{L,\mathbf{p}}$. But collectively, they
share certain patterns among themselves, specifically they share all their
parameters except for the ones in the final layer (see Figure 3). We formally
define the _deepshared_ CovNet kernel as

$c_{\rm ds}(\text{\bf u},\text{\bf
v})=\sum_{r=1}^{R}\sum_{s=1}^{R}\lambda_{r,s}\,g_{r}(\text{\bf
u})\,g_{s}(\text{\bf v}),\qquad\text{\bf u},\text{\bf v}\in\mathcal{Q},$ (2.9)

where $g_{1},\ldots,g_{R}$ are networks with shared structures as defined in
(2.2). A schematic representation of the structure (2.9) is shown in Figure 3.
We also define the class of deepshared CovNet kernels and the corresponding
operators as

$\displaystyle\mathscr{F}^{\rm ds}_{R,L,\mathbf{p},\sigma}$
$\displaystyle=\big{\\{}c_{\rm ds}\text{ of the form
\eqref{eq:deepshared_covnet_kernel}}:\Lambda=((\lambda_{r,s}))\succeq\mathrm{0},g_{1},\ldots,g_{R}\text{
of the form \eqref{eq:deepshared_neural_network}}\big{\\}}$
$\displaystyle\widetilde{\mathscr{F}}^{\rm ds}_{R,L,\mathbf{p},\sigma}$
$\displaystyle=\big{\\{}\mathcal{G}:\mathcal{G}\text{ is the integral operator
associated with kernel }g\in\mathscr{F}^{\rm
ds}_{R,L,\mathbf{p},\sigma}\big{\\}}.$ (2.10)

Figure 3: A schematic representation of the deepshared CovNet structure.
Starting from the inputs u and v, $R$ deep neural networks with shared weights
are fitted to get the outputs $g_{1}(\text{\bf u}),\ldots,g_{R}(\text{\bf u})$
and $g_{1}(\text{\bf v}),\ldots,g_{R}(\text{\bf v})$. These outputs are
combined by cross-multiplication, and finally addition after multiplication by
weights $\lambda_{r,s}$, to produce $c_{\rm ds}(\text{\bf u},\text{\bf v})$.

The shared structure drastically reduces the number of parameters of the
model. A deepshared CovNet kernel from the class $\mathscr{F}^{\rm
ds}_{L,\mathbf{p},R,\sigma}$ requires
$\sum_{l=0}^{L-1}(p_{l}+1)p_{l+1}+R(p_{L}+1)+R(R+1)/2$ parameters, compared to
$R\big{(}\sum_{l=0}^{L}(p_{l}+1)p_{l+1}\big{)}+R(R+1)/2$ for a deep CovNet
kernel from the class $\mathscr{F}^{\rm d}_{L,\mathbf{p},R,\sigma}$. To
appreciate this, assume that each of the hidden layers have the same width
$R$, i.e., $p_{1}=\cdots=p_{L}=R$. In this case, the deepshared kernel
contains $\mathcal{O}(R^{2})$ parameters, compared to $\mathcal{O}(R^{3})$
parameters for the deep kernel.

Similar to the shallow CovNet model, both the deep and the deepshared models
are universal approximators, i.e., they can approximate any covariance kernel
up to any desired accuracy.

######  Theorem 2.

Let $\mathcal{C}$ be a covariance operator on $\mathcal{L}_{2}(\mathcal{Q})$
with kernel $c$. Also, assume that the activation function $\sigma$ is
sigmoidal. Then, for every $\epsilon>0$, we can find a deep CovNet kernel
$c_{\rm d}$ and a deepshared CovNet kernel $c_{\rm ds}$ such that

$\int_{\mathcal{Q}\times\mathcal{Q}}\big{|}c(\text{\bf u},\text{\bf v})-c_{\rm
d}(\text{\bf u},\text{\bf v})\big{|}^{2}\,\mathrm{d}\text{\bf
u}\,\mathrm{d}\text{\bf v}\leq\epsilon\text{ and
}\int_{\mathcal{Q}\times\mathcal{Q}}\big{|}c(\text{\bf u},\text{\bf v})-c_{\rm
ds}(\text{\bf u},\text{\bf v})\big{|}^{2}\,\mathrm{d}\text{\bf
u}\,\mathrm{d}\text{\bf v}\leq\epsilon.$

If in addition $c$ is continuous, then the same conclusion holds uniformly.
That is, for every $\epsilon>0$, we can find $c_{\rm d}$ and $c_{\rm ds}$ such
that

$\sup_{\text{\bf u},\text{\bf v}\in\mathcal{Q}}\big{|}c(\text{\bf u},\text{\bf
v})-c_{\rm d}(\text{\bf u},\text{\bf v})\big{|}\leq\epsilon\text{ and
}\sup_{\text{\bf u},\text{\bf v}\in\mathcal{Q}}\big{|}c(\text{\bf u},\text{\bf
v})-c_{\rm ds}(\text{\bf u},\text{\bf v})\big{|}\leq\epsilon.$

Although both the models are universal approximators, the depth, width and the
number of components $R$ of the approximator for the two models can be
different. Moreover, although this result is similar to the one for the
shallow model, we expect the deep and the deepshared approximator to have a
much smaller number of parameters compared to the shallow approximator (Poggio
et al., 2017).

######  Remark 3.

The proof of the theorem again depends on the universal approximation property
of deep neural networks, and the sigmoidal condition on the activation
function is not necessary. From the proof of the theorem, one can see that the
universal approximation property of the deep (and the deepshared) CovNet
structure is guaranteed whenever the associated class of deep neural networks
has the universal approximation property. In particular, deep and deepshared
CovNet models with the highly popular ReLU activation:
$\sigma(t)=\max\\{t,0\\}$ are also universal approximators.

Theorems 1 and 2 justify the use of covariance networks (shallow, deep or
deepshared) for modeling the covariance kernel $c$. Given a fixed width $R$
and activation function $\sigma$, approximation by a shallow CovNet amounts to
determining a $\mathcal{G}$ in the class $\widetilde{\mathscr{F}}^{\rm
sh}_{R,\sigma}$ that is closest to $\mathcal{C}$ in terms of the Hilbert-
Schmidt norm, i.e.,

$\widehat{\mathcal{C}}^{\rm
sh}_{R,\sigma}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
sh}_{R,\sigma}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{C}-\mathcal{G}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}^{2}.$
(2.11)

Given a sample of random fields
$\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\overset{{\rm
i.i.d.}}{\sim}\mathcal{X}$ in $\mathcal{L}_{2}(\mathcal{Q})$, with covariance
$\mathcal{C}$, we can replace $\mathcal{C}$ in (2.11) by the empirical
covariance operator
$\widehat{\mathcal{C}}_{N}=N^{-1}\sum_{n=1}^{N}\mathcal{X}_{n}\otimes\mathcal{X}_{n}$
to obtain an estimator:

$\widehat{\mathcal{C}}_{R,N}^{\rm
sh}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
sh}_{R,\sigma}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}^{2}.$
(2.12)

We call this the shallow CovNet estimator. Similarly, for given width $R$,
depth $L$ and activation $\sigma$, we can define the estimators based on the
deep and deepshared models as

$\displaystyle\widehat{\mathcal{C}}^{\rm
d}_{L,R,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
d}_{R,L,\mathbf{p},\sigma}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\quad\text{
and }\quad\widehat{\mathcal{C}}^{\rm
ds}_{L,R,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
ds}_{R,L,\mathbf{p},\sigma}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2},$
(2.13)

which we call the deep CovNet and the deepshared CovNet estimators,
respectively. The estimators $\widehat{\mathcal{C}}^{\rm sh}_{R,N}$,
$\widehat{\mathcal{C}}^{\rm d}_{L,R,N}$, and $\widehat{\mathcal{C}}^{\rm
ds}_{L,R,N}$ can be seen as regularized versions of the empirical covariance,
by projection into the corresponding CovNet classes. Nevertheless, it is
crucial to note here that, although the definition of the estimators involve
$\widehat{\mathcal{C}}_{N}$, we never actually need to form the empirical
covariance in order to construct them. The estimators
$\widehat{\mathcal{C}}_{R,N}^{\rm sh}$, $\widehat{\mathcal{C}}^{\rm
d}_{L,R,N}$ and $\widehat{\mathcal{C}}^{\rm ds}_{L,R,N}$ can be computed
directly at the level of the data, without the need to ever store or access
the $2d$-dimensional object $\widehat{\mathcal{C}}_{N}$. We discuss the
implementation details in the next section.

Note that although these estimators depend on the widths $p_{1},\ldots,p_{L}$,
we have suppressed it in the notation for ease of exposition. Also, in our
numerical experiments, we have used $p_{1}=\cdots=p_{L}=R$, which justifies
this notation. This choice is motivated by the empirical evidence that
suggests using the same width for all the hidden layers (see Bengio, 2012,
Section 19.3.2).

######  Remark 4.

Observe the notation in (2.11)–(2.13). In either of the equations, we cannot
guarantee that the minimizer is unique. Here, and throughout the article, by
the notation
$\widehat{\mathcal{G}}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}$,
we mean that $\widehat{\mathcal{G}}$ is an element (out of possibly many) of
the class $\widetilde{\mathscr{F}}$ satisfying
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{G}}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\leq{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}$
for all $\mathcal{G}\in\widetilde{\mathscr{F}}$. Note that this non-uniqueness
does not affect the subsequent development, in particular the asymptotic
theory for the estimator.

##  3 Practical implementation

Note that for all three CovNet models (2.1), (2.6) and (2.9), the covariance
kernel is of the form

$\sum_{r=1}^{R}\sum_{s=1}^{R}\lambda_{r,s}\,g_{r}(\text{\bf
u})\,g_{s}(\text{\bf v}),\qquad\text{\bf u},\text{\bf v}\in\mathcal{Q},$ (3.1)

where $\Lambda:=((\lambda_{r,s}))$ is positive semi-definite, and
$g_{1},\ldots,g_{R}$ are allowed to vary keeping up to the model under
consideration. In particular, $g_{r}(\text{\bf
u})=\sigma(\mathbf{w}_{r}^{\top}\text{\bf u}+b_{r})$ for shallow CovNet,
$g_{r}$’s are the individual deep neural networks from the class
$\mathcal{D}_{L,\mathbf{p}}$ (2.2) for deep CovNet, and $g_{r}$’s are jointly
defined as in (2.2) for deepshared CovNet. We denote the generic class of all
such kernels (with the additional structures on the functions
$g_{1},\ldots,g_{R}$) by $\mathscr{F}_{R}$ and the corresponding class of
operators by $\widetilde{\mathscr{F}}_{R}$. Thus, $\mathscr{F}_{R}$ (resp.,
$\widetilde{\mathscr{F}}_{R}$) can be $\mathscr{F}_{R,\sigma}^{\rm sh}$,
$\mathscr{F}_{R,L,\mathbf{p},\sigma}^{\rm d}$, or
$\mathscr{F}_{R,L,\mathbf{p},\sigma}^{\rm ds}$ (resp.,
$\widetilde{\mathscr{F}}_{R,\sigma}^{\rm sh}$,
$\widetilde{\mathscr{F}}_{R,L,\mathbf{p},\sigma}^{\rm d}$, or
$\widetilde{\mathscr{F}}_{R,L,\mathbf{p},\sigma}^{\rm ds}$) depending on the
situation. Here, and throughout, we suppress the dependence on
$L,p_{1},\ldots,p_{L}$ and $\sigma$ for convenience, unless specifically
mentioned.

For a given $R\in\mathbb{N}$, the CovNet structure (3.1) is completely
determined by the parameters of $g_{1},\ldots,g_{R}$, and the coefficients
$\Lambda=((\lambda_{r,s}))$. In particular, apart from $\Lambda$, these
parameters are $\mathbf{w}_{1},\ldots,\mathbf{w}_{R},b_{1},\ldots,b_{R}$ for
the shallow CovNet (2.1), the weights and the biases of the individual deep
neural networks $g_{1},\ldots,g_{R}$ for the deep CovNet (2.6), and
$W_{1},\ldots,W_{L},\mathbf{b}_{1},\ldots,\mathbf{b}_{L}$ and
$\bm{\omega}_{1},\ldots,\bm{\omega}_{R},\beta_{1},\ldots,\beta_{R}$ for the
deepshared CovNet (2.9). Thus, obtaining the estimators in (2.12) or (2.13) is
equivalent to finding these parameters minimizing the corresponding criterion

$\ell:=\ell(\Theta)={|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}.$

Here, we use $\Theta$ to denote all the estimable parameters (i.e., the
parameters of $g_{1},\ldots,g_{R}$, and $\Lambda$), taking into account the
positive-definiteness of $\Lambda$ (which reduces the number of _free
parameters_).

As already mentioned, we do not need to form the tensor
$\widehat{\mathcal{C}}_{N}$ (or candidate tensor $\mathcal{G}$) to minimize
$\ell$. The trick is to not fit the covariance directly, but to instead fit
the observed fields $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ themselves by
_neural networks with shared structures_. To be precise, consider the fields

$\mathcal{X}^{\rm NN}_{n}(\text{\bf
u})=\sum_{r=1}^{R}\xi_{n,r}\,g_{r}(\text{\bf u}),\qquad n=1,\ldots,N,$ (3.2)

where $\xi_{n,r}\in\mathbb{R}$ for $n=1,\ldots,N$, $r=1,\ldots,R$, and
$g_{1},\ldots,g_{R}$ are the constituents of the CovNet model (3.1) under
consideration. By construction, the fields $\mathcal{X}_{1}^{\rm
NN},\ldots,\mathcal{X}_{N}^{\rm NN}$ are themselves neural networks with
shared components $g_{1},\ldots,g_{R}$ (and hence shared parameters), but
potentially different coefficients $\xi_{n,r}$. Define the operator

$\mathcal{G}_{R,N}^{\rm NN}=\frac{1}{N}\sum_{n=1}^{N}(\mathcal{X}^{\rm
NN}_{n}-\bar{\mathcal{X}}^{\rm NN})\otimes(\mathcal{X}^{\rm
NN}_{n}-\bar{\mathcal{X}}^{\rm NN}),$ (3.3)

which is the empirical covariance based on the neural networks
$\mathcal{X}^{\rm NN}_{1},\ldots,\mathcal{X}^{\rm NN}_{N}$, and let
$\widetilde{\mathscr{F}}^{\rm NN}_{R,N}$ be the class of all such covariance
operators:

$\displaystyle\widetilde{\mathscr{F}}^{\rm NN}_{R,N}=\big{\\{}\text{all
empirical covariance operators of the form
\eqref{eq:covnet_empirical_formulation}}\big{\\}}.$ (3.4)

Because of the shared structure of the networks $\mathcal{X}^{\rm NN}_{n}$,
the kernel of the operator $\mathcal{G}^{\rm NN}_{R,N}$ has the CovNet
structure (3.1):

$g^{\rm NN}_{R,N}(\text{\bf u},\text{\bf
v})=\sum_{r=1}^{R}\sum_{s=1}^{R}\lambda_{r,s}\,g_{r}(\text{\bf
u})\,g_{s}(\text{\bf v}),$

where
$\lambda_{r,s}=N^{-1}\sum_{n=1}^{N}(\xi_{r,n}-\bar{\xi}_{r})(\xi_{s,n}-\bar{\xi}_{s})$.
At the same time, if $N>R$, any CovNet operator from the class
$\widetilde{\mathscr{F}}_{R}$ can be written as an empirical covariance
operator of the form (3.3). Specifically, for every
$\mathcal{G}\in\widetilde{\mathscr{F}}_{R}$, we can find $N$ networks of the
form (3.2) such that $\mathcal{G}$ is the empirical covariance operator of
those $N$ networks. This should be intuitively clear, but we nevertheless
state this formally below, and a detailed construction is shown in Appendix C.

######  Proposition 1.

If $N>R$, then $\widetilde{\mathscr{F}}^{\rm
NN}_{R,N}=\widetilde{\mathscr{F}}_{R}$.

This simple correspondence between the CovNet operator class
$\widetilde{\mathscr{F}}_{R}$ and the class of empirical covariances of neural
networks with shared structure (3.4) is of great consequence in estimating the
CovNet model based on the observed data. Note that the criterion $\ell$ is
non-convex in the parameters, so we cannot find the explicit minimizer.
Instead, we need to rely on some iterative minimization procedure, e.g.,
gradient descent or its variants (Buduma and Locascio, 2017, Chapters 2 and
4). The application of gradient descent requires us to calculate the gradient
of the minimization criterion. But, with modern optimization routines, this
can be done numerically on a computer, without the need to compute the
derivatives analytically. In particular, the special neural network structure
of our method allows us to employ _automatic differentiation_ techniques to
efficiently compute the derivative at machine precision (Baydin et al., 2018).
In essence, the minimizer can be efficiently obtained if we can compute the
criterion efficiently. This is where the empirical covariance formulation
(3.2) and (3.3) come in handy. With this formulation, we compute the criterion
$\ell$ as a function of the parameters of $g_{1},\ldots,g_{R}$, and
coefficients $\xi_{n,r}$ instead of $\lambda_{r,s}$. At each step of gradient
descent, we obtain the fields $\mathcal{X}_{n}^{\rm NN}$ as feed-forward
neural networks. The minimization criterion $\ell$ can be computed by simply
computing inner-products between the observed fields $\mathcal{X}_{n}$ and the
fitted networks $\mathcal{X}^{\rm NN}_{n}$, as shown below. For simplicity, we
assume that the observed fields $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ are
centered, so that the empirical covariance is
$\widehat{\mathcal{C}}_{N}=N^{-1}\sum_{n=1}^{N}\mathcal{X}_{n}\otimes\mathcal{X}_{n}$.
We also assume that the fitted networks $\mathcal{X}^{\rm
NN}_{1},\ldots,\mathcal{X}^{\rm NN}_{N}$ are centered, so that their empirical
covariance is $\mathcal{G}^{\rm NN}_{R,N}=N^{-1}\sum_{n=1}^{N}\mathcal{X}^{\rm
NN}_{n}\otimes\mathcal{X}^{\rm NN}_{n}$. With these, we get the following
formula for the minimization criterion:

$\displaystyle\ell:={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\widehat{\mathcal{C}}_{N}-\mathcal{G}^{\rm
NN}_{R,N}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{2}^{2}=\frac{1}{N^{2}}\sum_{n=1}^{N}\sum_{m=1}^{N}\langle\mathcal{X}_{n},\mathcal{X}_{m}\rangle^{2}+\frac{1}{N^{2}}\sum_{n=1}^{N}\sum_{m=1}^{N}\langle\mathcal{X}^{\rm
NN}_{n},\mathcal{X}^{\rm
NN}_{m}\rangle^{2}-\frac{2}{N^{2}}\sum_{n=1}^{N}\sum_{m=1}^{N}\langle\mathcal{X}_{n},\mathcal{X}^{\rm
NN}_{m}\rangle^{2}.$

The detailed derivations are shown in Appendix C.1.

######  Remark 5.

The alternative formulation also helps us in imposing positive semi-
definiteness on $\Lambda$. After estimating the parameters from the
reformulated problem, we obtain $\lambda_{r,s}$ as
$N^{-1}\sum_{n=1}^{N}(\xi_{n,r}-\bar{\xi}_{r})(\xi_{n,s}-\bar{\xi}_{s})$. By
virtue of this construction, the resulting matrix $\Lambda=(\lambda_{r,s})$ is
automatically positive semi-definite. This is quite useful, as it circumvents
the need to work with a constrained optimisation problem on a cone.

In practice, we observe the data on a grid of size $D=K_{1}\times\cdots\times
K_{d}$, say $\\{\text{\bf u}_{1},\ldots,\text{\bf u}_{D}\\}$. Let us denote
the $i$-th measurement corresponding to the $n$-th field by $X_{ni}$ for
$n=1,\ldots,N$, $i=1,\ldots,D$. So, we can store the observed fields as an
$N\times D$ matrix $\mathbf{X}=((X_{ni}))$. Similarly, the fitted networks
$\mathcal{X}^{\rm NN}_{1},\ldots,\mathcal{X}^{\rm NN}_{N}$ can be evaluated at
the $D$ grid points and all of these can be stored as an $N\times D$ matrix
$\mathbf{X}^{\rm NN}=((X^{\rm NN}_{ni}))$, where $X^{\rm
NN}_{ni}=\mathcal{X}^{\rm NN}_{n}(\text{\bf u}_{i})$. We can approximate
$\langle\mathcal{X}_{n},\mathcal{X}_{m}\rangle$ by the average over the grid
points, i.e., $\langle\mathcal{X}_{n},\mathcal{X}_{m}\rangle\cong
D^{-1}\sum_{i=1}^{D}X_{ni}X_{mi}$. It is easy to see that this is the
$(n,m)$-th element of the $N\times N$ matrix
$D^{-1}\mathbf{X}\mathbf{X}^{\top}$. Similarly, we approximate
$\langle\mathcal{X}^{\rm NN}_{n},\mathcal{X}^{\rm NN}_{m}\rangle$ and
$\langle\mathcal{X}_{n},\mathcal{X}^{\rm NN}_{m}\rangle$ by the corresponding
averages $D^{-1}\sum_{i=1}^{D}X^{\rm NN}_{ni}X^{\rm NN}_{mi}$ and
$D^{-1}\sum_{i=1}^{D}X_{ni}X^{\rm NN}_{mi}$, which are the $(n,m)$-th elements
of $D^{-1}\mathbf{X}^{\rm NN}\mathbf{X}^{\rm NN,\top}$ and
$D^{-1}\mathbf{X}\mathbf{X}^{\rm NN,\top}$, respectively. Thus, apart from the
computation of $\mathbf{X}^{\rm NN}$, the computational cost of $\ell$ is
$\mathcal{O}(N^{2}D)$. Moreover, to store the model, we only need to store the
parameters of $g_{1},\ldots,g_{R}$, and the coefficient matrix $\Lambda$,
which is completely free of the grid size $D$. In particular, this amounts to
a storage cost of $\mathcal{O}(R^{2}+Rd)$ for the shallow CovNet model,
$\mathcal{O}\big{(}R^{2}+Rp_{L}+R\sum_{l=0}^{L-1}(p_{l}+1)p_{l+1}\big{)}$ for
the deep CovNet model, and
$\mathcal{O}\big{(}R^{2}+Rp_{L}+\sum_{l=0}^{L-1}(p_{l}+1)p_{l+1}\big{)}$ for
the deepshared CovNet model ($p_{0}=d$ for the latter two). It is easy to see
the savings relative to the empirical covariance, which requires
$\mathcal{O}(ND^{2})$ computations and $\mathcal{O}(D^{2})$ storage.

In the above discussion, we have not addressed the computational requirements
for $\mathbf{X}^{\rm NN}$. It is not difficult to show that for a fixed set of
parameters, computation of the matrix $\mathbf{X}^{\rm NN}$ needs
$\mathcal{O}(DR(N+d))$ operations for the shallow CovNet,
$\mathcal{O}\big{(}DR(N+\sum_{l=0}^{L-1}p_{l}p_{l+1}+p_{L})\big{)}$ operations
for the deep CovNet, and
$\mathcal{O}\big{(}D(NR+\sum_{l=0}^{L-1}p_{l}p_{l+1}+Rp_{L})\big{)}$
operations for the deepshared CovNet (see Appendix C.3). Thus, for a fixed set
of parameters, the computational cost for the evaluation of $\ell$ remains
linear in the grid size $D$ for all three CovNet models. Of course, we need to
re-evaluate the criterion for each step of the gradient descent algorithm. But
that is also the case for other modern machine learning methods. Moreover, the
computation can be sped up by considering other techniques from machine
learning, such as the stochastic or mini-batch version of gradient descent and
parallel computing (Bengio, 2012; Buduma and Locascio, 2017).

######  Remark 6.

We have not tried to find analytic expression for the derivative of $\ell$ as
a function of the parameters. Instead, we focused more on evaluating the
criterion efficiently, and rely on automatic differentiation to compute the
gradient. There are three reasons for doing this. Firstly, because of the
complex neural network structure, finding analytic expressions for the
gradient is cumbersome. This becomes more relevant for the deep and the
deepshared CovNet structures. Secondly, we have at our disposal modern
optimization routines, which are very efficient in automatic differentiation,
especially with neural network structures such as ours. In our code, we have
used the autograd feature of pytorch (https://pytorch.org/). Finally, even if
we compute the derivatives analytically, when implementing the method on a
computer the accumulation of errors for analytic derivatives may sometimes be
quite large, especially for complex structures such as neural networks.
Automatic differentiation, on the other hand, produces results which are exact
up to machine precision, and thus are preferred to analytic derivatives
(Baydin et al., 2018).

######  Remark 7.

In our derivations, we have assumed that the fields
$\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ as well as $\mathcal{X}_{1}^{\rm
NN},\ldots,\mathcal{X}_{N}^{\rm NN}$ are centered. In practice, we can center
the observed fields by subtracting the mean (empirical or estimated by some
other method), with negligible computational overhead. For the fitted fields
$\mathcal{X}_{1}^{\rm NN},\ldots,\mathcal{X}_{N}^{\rm NN}$, because of their
shared structure, the mean turns out to be $\bar{\mathcal{X}}^{\rm
NN}(\text{\bf u})=\sum_{r=1}^{R}\bar{\xi}_{r}\,g_{r}(\text{\bf u})$. So,
centering the fitted fields boils down to centering the coefficients
$\xi_{n,r}$. We can use another approach, where we do not center the fields
(observed or fitted) beforehand and minimize a slightly different criterion.
In this case, we also get an estimate of the mean as a by-product (see
Appendix C.2 for details).

##  4 Eigendecomposition of the estimated covariance operator

Once we estimate the covariance, it is important to be able to _manipulate_
it, e.g., for regression, prediction, or even for visualization purposes. For
such tasks, typical manipulations involve inverting the covariance operator or
obtaining its eigendecomposition, either of which may be quite demanding in
practice. For instance, for data observed on a grid of size
$D=K_{1}\times\cdots\times K_{d}$, the empirical covariance is stored as a
$D\times D$ matrix. The inversion in this case requires $\mathcal{O}(D^{3})$
operations, which is highly demanding and sometimes even prohibitive. Even if
the inverse is constructed, it is available only at the $D\times D$ pre-
specified locations. To evaluate the inverse (or the covariance itself, for
that matter) at any other location, as required e.g., in kriging, one needs to
apply some sort of interpolation or smoothing on a high-dimensional (in our
case, $\mathbb{R}^{d}\times\mathbb{R}^{d}$) object, which can be even more
demanding than the inversion itself. Finally, the cost of storing the inverse
and/or the eigenfunctions adds another layer of burden.

By contrast, the proposed CovNet estimators enjoy considerable advantage in
this respect. The special form of the CovNet operators allow us to easily
compute their eigendecomposition. Note that our estimated CovNet kernels are
of the form

$\widehat{c}(\text{\bf u},\text{\bf
v})=\sum_{r=1}^{R}\sum_{s=1}^{R}\widehat{\lambda}_{r,s}\,\widehat{g}_{r}(\text{\bf
u})\,\widehat{g}_{s}(\text{\bf v}),\qquad\text{\bf u},\text{\bf
v}\in\mathcal{Q}.$

Thus, for the estimated covariance operator $\widehat{\mathcal{C}}$ and for
any $f\in\mathcal{L}_{2}(\mathcal{Q})$,

$\displaystyle\widehat{\mathcal{C}}f(\text{\bf
u})=\int_{\mathcal{Q}}\widehat{c}(\text{\bf u},\text{\bf v})f(\text{\bf
v})\,\mathrm{d}\text{\bf
v}=\sum_{r=1}^{R}\sum_{s=1}^{R}\widehat{\lambda}_{r,s}\,\widehat{g}_{r}(\text{\bf
u})\int_{\mathcal{Q}}\widehat{g}_{s}(\text{\bf v})\,f(\text{\bf
v})\,\mathrm{d}\text{\bf v}=\sum_{r=1}^{R}a_{r}\,\widehat{g}_{r}(\text{\bf
u}),$

where
$a_{r}=\sum_{s=1}^{R}\widehat{\lambda}_{r,s}\int_{\mathcal{Q}}\widehat{g}_{s}(\text{\bf
v})\,f(\text{\bf v})\,\mathrm{d}\text{\bf v}$. This shows that the
eigenfunctions of $\widehat{\mathcal{C}}$ are of the form $\psi(\text{\bf
u})=\sum_{r=1}^{R}a_{r}\,\widehat{g}_{r}(\text{\bf u})$ for some
$a_{1},\ldots,a_{R}\in\mathbb{R}$. Now, for such a function $\psi$,

$\|\psi\|^{2}=\sum_{r=1}^{R}\sum_{s=1}^{R}a_{r}\,a_{s}\int_{\mathcal{Q}}\widehat{g}_{r}(\text{\bf
u})\,\widehat{g}_{s}(\text{\bf u})\,\mathrm{d}\text{\bf
u}=\sum_{r=1}^{R}\sum_{s=1}^{R}a_{r}\,a_{s}\,\widetilde{g}(r,s)=\mathbf{a}^{\top}\widetilde{\mathrm{G}}\,\mathbf{a},$

where $\widetilde{g}(r,s)=\int_{\mathcal{Q}}\widehat{g}_{r}(\text{\bf
u})\,\widehat{g}_{s}(\text{\bf u})\,\,\mathrm{d}\text{\bf u}$,
$\mathbf{a}=(a_{1},\ldots,a_{R})^{\top}$ and
$\widetilde{\mathrm{G}}=((\widetilde{g}(r,s)))_{1\leq r,s\leq R}$. Also,

$\displaystyle\langle\widehat{\mathcal{C}}\psi,\psi\rangle$
$\displaystyle=\iint_{\mathcal{Q}\times\mathcal{Q}}\widehat{c}(\text{\bf
u},\text{\bf v})\,\psi(\text{\bf u})\,\psi(\text{\bf v})\,\mathrm{d}\text{\bf
u}\,\mathrm{d}\text{\bf v}$
$\displaystyle=\sum_{r=1}^{R}\sum_{s=1}^{R}\widehat{\lambda}_{r,s}\iint_{\mathcal{Q}\times\mathcal{Q}}\widehat{g}_{r}(\text{\bf
u})\,\widehat{g}_{s}(\text{\bf v})\,\psi(\text{\bf u})\,\psi(\text{\bf
v})\,\,\mathrm{d}\text{\bf u}\,\,\mathrm{d}\text{\bf v}$
$\displaystyle=\sum_{r=1}^{R}\sum_{s=1}^{R}\widehat{\lambda}_{r,s}\sum_{i=1}^{R}\sum_{j=1}^{R}a_{i}\,a_{j}\int_{\mathcal{Q}}\widehat{g}_{r}(\text{\bf
u})\,\widehat{g}_{i}(\text{\bf u})\,\mathrm{d}\text{\bf
u}\int_{\mathcal{Q}}\widehat{g}_{s}(\text{\bf v})\,\widehat{g}_{j}(\text{\bf
v})\,\mathrm{d}\text{\bf v}$
$\displaystyle=\sum_{i=1}^{R}\sum_{j=1}^{R}a_{i}\,a_{j}\bigg{(}\sum_{r=1}^{R}\sum_{s=1}^{R}\widehat{\lambda}_{r,s}\,\widetilde{g}(r,i)\,\widetilde{g}(s,j)\bigg{)}=\sum_{i=1}^{R}\sum_{j=1}^{R}a_{i}\,a_{j}(\widetilde{\mathrm{G}}\Lambda\widetilde{\mathrm{G}})_{i,j}=\mathbf{a}^{\top}\widetilde{\mathrm{G}}\Lambda\widetilde{\mathrm{G}}\,\mathbf{a}.$

Thus, finding the leading eigenvalue and eigenfunction of
$\widehat{\mathcal{C}}$ reduces to maximizing
$\mathbf{a}^{\top}\widetilde{\mathrm{G}}\Lambda\widetilde{\mathrm{G}}\,\mathbf{a}$
subject to $\mathbf{a}^{\top}\widetilde{\mathrm{G}}\,\mathbf{a}=1$. This
amounts to solving

$(\widetilde{\mathrm{G}}\Lambda\widetilde{\mathrm{G}}-\eta\widetilde{\mathrm{G}})\,\mathbf{a}=\mathbf{0}.$

Again, if $\psi_{i}(\text{\bf
u})=\sum_{r=1}^{R}a_{i,r}\,\widehat{g}_{r}(\text{\bf u})$, then we can
similarly show that

$\langle\psi_{i},\psi_{j}\rangle=\mathbf{a}_{i}^{\top}\widetilde{\mathrm{G}}\,\mathbf{a}_{j}\text{
and
}\langle\widehat{\mathcal{C}}\psi_{i},\psi_{j}\rangle=\mathbf{a}_{i}^{\top}\widetilde{\mathrm{G}}\Lambda\widetilde{\mathrm{G}}\,\mathbf{a}_{j},$

where $\mathbf{a}_{i}=(a_{i,1},\ldots,a_{i,R})^{\top}$ is the vector of
coefficients of $\psi_{i}$. Thus, finding the subsequent eigenvalues and
eigenfunctions also amounts to solving
$(\widetilde{\mathrm{G}}\Lambda\widetilde{\mathrm{G}}-\eta\widetilde{\mathrm{G}})\,\mathbf{a}=\mathbf{0}$,
with added orthogonality constraints. In summary, finding the eigensystem of
the CovNet operator $\widehat{\mathcal{C}}$ boils down to finding the solution
of the generalized eigenvalue problem (Golub and Van Loan, 2013, Chapter 7)
involving the non-negative definite matrices
$\widetilde{\mathrm{G}}\Lambda\widetilde{\mathrm{G}}$ and
$\widetilde{\mathrm{G}}$. Several optimization routines are available to
obtain the solution. Also, this can be done very efficiently since the
matrices $\Lambda$ and $\widetilde{\mathrm{G}}$ involved in the computations
are of the order $R\times R$, and typical values of $R$ will be much smaller
than $D=K_{1}\times\ldots\times K_{d}$. The matrix $\Lambda$ is obtained
during the estimation procedure. The only bottleneck is the computation of the
matrix $\widetilde{\mathrm{G}}$, which involves the integrals

$\widetilde{g}(r,s)=\int_{\mathcal{Q}}\widehat{g}_{r}(\text{\bf
u})\,\widehat{g}_{s}(\text{\bf u})\,\mathrm{d}\text{\bf u}.$

These are integrals on a compact subset of $\mathbb{R}^{d}$. When $d$ is
moderate, we can approximate the integral using Monte-Carlo methods, while for
large $d$, we can resort to using quasi-Monte-Carlo methods (Dick, Kuo and
Sloan, 2013). In typical FDA applications, $d$ is $2,3$ or $4$ (corresponding
to spatial/spatio-temporal data on $\mathbb{R}^{2}$ and $\mathbb{R}^{3}$), and
it suffices to use Monte-Carlo integration. For this, we generate independent
observations $\text{\bf u}_{1},\ldots,\text{\bf u}_{M}$ distributed uniformly
on $\mathcal{Q}$, and approximate the integral as

$\widetilde{g}(r,s)\simeq\frac{1}{M}\sum_{j=1}^{M}\widehat{g}_{r}(\text{\bf
u}_{j})\,\widehat{g}_{s}(\text{\bf u}_{j}).$

Also, when the functions $\widehat{g}_{1},\ldots,\widehat{g}_{R}$ are bounded
(e.g., when the activation $\sigma$ is sigmoidal), we can control the
approximation error up to any desired accuracy by selecting $M$ large enough.
After generating the observations $\text{\bf u}_{1},\ldots,\text{\bf u}_{M}$,
$\widetilde{\mathrm{G}}$ can be obtained by passing them through
$\widehat{g}_{1},\ldots,\widehat{g}_{R}$ to create an $R\times M$ matrix
$\mathrm{G}$, and then computing the outer-product
$\mathrm{G}{\mathrm{G}}^{\top}$. It can be verified that the overall
computational cost remains linear in $M$ (see Appendix C.3). So, even with a
large value of $M$, the computational time is quite small. Moreover, after
obtaining the eigendecomposition, the complete eigenstructure can be stored
using an $R\times R$ matrix of coefficients
$\widehat{\mathbf{A}}=(\widehat{\mathbf{a}}_{1}^{\top},\ldots,\widehat{\mathbf{a}}_{R}^{\top})^{\top}$
and a vector of eigenvalues
$\widehat{\bm{\eta}}=(\widehat{\eta}_{1},\ldots,\widehat{\eta}_{R})$, in
addition to the already estimated parameters. The usefulness of the
eigendecomposition is shown in Section 5.2.

##  5 Empirical study

We now demonstrate the usefulness of the proposed methods by means of a
variety of numerical examples. We start with some simulated examples, where we
generate the data from a Gaussian process (Adler and Taylor, 2007, Chapter 1)
on $[0,1]^{d}$ with mean $0$ and variance $\mathcal{C}$. We consider the
following five choices for the kernel $c$.

  1. Ex 1

_Brownian sheet_ : $c(\text{\bf u},\text{\bf v})=c_{\rm
bm}(u_{1},v_{1})\times\cdots\times c_{\rm bm}(u_{d},v_{d})$ for $\text{\bf
u},\text{\bf v}\in[0,1]^{d}$, where $c_{\rm bm}(u,v)=\min\\{u,v\\}$ is the
covariance of the standard Brownian motion (Adler and Taylor, 2007, Sec.
1.4.3).

  2. Ex 2

_Rotated Brownian sheet_ : $c(\text{\bf u},\text{\bf
v})=\tilde{c}(\mathrm{O}\text{\bf u},\mathrm{O}\text{\bf v})$, where
$\mathrm{O}$ is a rotation matrix and $\tilde{c}$ is the covariance kernel of
the Brownian sheet from Ex 1.

  3. Ex 3

_Integrated Brownian sheet_ : $c(\text{\bf u},\text{\bf v})=c_{\rm
ibm}(u_{1},v_{1})\times\cdots\times c_{\rm ibm}(u_{d},v_{d})$ for $\text{\bf
u},\text{\bf v}\in[0,1]^{d}$, where $c_{\rm
ibm}(u,v)=(u^{2}/2)\,(v-u/3)\textrm{{1}}\\{u\leq
v\\}+(v^{2}/2)\,(u-v/3)\textrm{{1}}\\{u>v\\}$ is the covariance of the
integrated Brownian motion.

  4. Ex 4

_Rotated integrated Brownian sheet_ : $c(\text{\bf u},\text{\bf
v})=\tilde{c}(\mathrm{O}\text{\bf u},\mathrm{O}\text{\bf v})$, where
$\mathrm{O}$ is a rotation matrix and $\tilde{c}$ is the covariance kernel of
the integrated Brownian sheet from Ex 3.

  5. Ex 5

_Matérn covariance_ : $c_{\nu}(\text{\bf u},\text{\bf
v})=2^{1-\nu}/\Gamma(\nu)\,(\sqrt{2\nu}\,\|\text{\bf u}-\text{\bf
v}\|_{d})^{\nu}\,K_{\nu}(\sqrt{2\nu}\,\|\text{\bf u}-\text{\bf v}\|_{d})$,
where $\Gamma$ is the gamma function, $K_{\nu}$ is the modified Bessel
function of the second kind and $\|\cdot\|_{d}$ is the Euclidean distance on
$\mathbb{R}^{d}$ (Rasmussen and Williams, 2006, Chapter 4). The Matérn
covariance is indexed by the parameter $\nu>0$, which regulates its
smoothness.

Note that the covariance kernels in Ex 1 and 2 are separable. We eliminate the
separability in Ex 3 and 4 by introducing a rotation of the domain. The Matérn
covariance in Ex 5 is stationary and isotropic, but not separable for any
finite $\nu$. On the other hand, none of the other covariances are stationary.
Ex 1 and 3 yield continuous but nowhere differentiable random fields, whereas
Ex 2 and 4 yield continuously differentiable random fields. For Ex 5, the
random fields are $\lceil\nu\rceil-1$ times differentiable in the mean-square
sense.

We carried out our experiments with $d=2$ and $d=3$, which we refer to as 2D
and 3D, respectively. For each experiment, we generated $N$ independent fields
at $K\times\cdots\times K$ regular grid points on $[0,1]^{d}$. Henceforth, we
refer to $K$ as the resolution. We used the three CovNet models (shallow, deep
and deepshared) on the generated data to estimate $\mathcal{C}$. To facilitate
comparison, we also consider the empirical covariance estimator and the best
separable covariance estimator (e.g., Dette, Dierickx and Kutta, 2020). For
each of these estimators, we compute the relative estimation error
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}/{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}$.
Note that

${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}=\iint_{[0,1]^{d}\times[0,1]^{d}}\big{(}\widehat{c}(\text{\bf
u},\text{\bf v})-c(\text{\bf u},\text{\bf v})\big{)}^{2}\,\mathrm{d}\text{\bf
u}\,\mathrm{d}\text{\bf v}~{}~{}\text{ and
}~{}~{}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}=\iint_{[0,1]^{d}\times[0,1]^{d}}c^{2}(\text{\bf
u},\text{\bf v})\,\mathrm{d}\text{\bf u}\,\mathrm{d}\text{\bf v}$

cannot always be computed analytically. So, we use a Monte-Carlo
approximation. We generate $M$ points $(\text{\bf u}_{1},\text{\bf
v}_{1}),\ldots,(\text{\bf u}_{M},\text{\bf v}_{M})$ from the uniform
distribution on $[0,1]^{d}\times[0,1]^{d}$ and approximate

${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\simeq\frac{1}{M}\sum_{i=1}^{M}\big{(}\widehat{c}(\text{\bf
u}_{i},\text{\bf v}_{i})-c(\text{\bf u}_{i},\text{\bf
v}_{i})\big{)}^{2}~{}~{}\text{ and
}~{}~{}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\simeq\frac{1}{M}\sum_{i=1}^{M}c^{2}(\text{\bf
u}_{i},\text{\bf v}_{i}).$

These are then used to approximate the relative errors of the estimators. The
advantage of using Monte-Carlo is that we can control the approximation error
up to any desired accuracy by selecting $M$ large enough. Also, by evaluating
the estimators on a different set of locations than where the data were
generated, we avoid committing an _inverse crime_ (Kaipio and Somersalo,
2005). In particular, we used $M=50000$ in 2D and $M=100000$ in 3D.

(a) Brownian sheet | (b) Rotated Brownian sheet
---|---
|
(c) Integrted Brownian sheet | (d) Rotated integrated Brownian sheet
|
Legend: Empirical  Best separable  Shallow  Deep  DeepShared
Figure 4:  Relative errors of different methods for different examples in 2D.
Results are reported for a fixed sample size of $500$ and varying resolution.
The numbers are averages based on $25$ simulation runs.

We also considered two different setups based on the sample size and the
resolution: (a) fixed resolution $K$ and varying sample size $N$ and (b) fixed
sample size $N$ and varying resolution $K$. Also, for the Matérn example, we
considered different values of $\nu$ with fixed sample size and resolution.
For setup (a), the results are unremarkable – the errors of all the estimators
decrease as $N$ increases. These are reported in Appendix F. The results for
setup (b) are rather interesting and exhibit the superiority of the CovNet
estimators. We show these results for 2D in Figures 4–5, where the reported
numbers are the average relative errors based on $25$ simulation runs. The
results for 3D are qualitatively similar, and we report them in Appendix F.

For the CovNet estimators, the results depend on the choice of hyperparameters
$R$ and $L$. We used $R=5,10,20,40,80$ for the shallow CovNet model, and
$L=2,3,4,R=5,10,20,40$ for the deep and the deepshared CovNet models in our
experiments. In Figures 4–5, we report the best result (i.e., minimum average
estimation error) obtained by each CovNet model. In Section 5.1, we discuss a
practical method to select the hyperparameters and exhibit the corresponding
results. It is seen there that the selection method yields values comparable
to the “best choice”. For all the CovNet models, we used the standard sigmoid
activation function $\sigma(t)=1/(1+\exp(-t))$. For the optimization involved
in fitting these models, we used the ADAM optimizer (Kingma and Ba, 2014)
available in pytorch.

(a) $N=250$, resolution $25\times 25$ | (b) $N=500$, $\nu=0.01$
---|---
|
Legend: Empirical  Best separable  Shallow  Deep  DeepShared
Figure 5:  Relative errors of different methods for the Matérn covariance
model in 2D. In (a), results are for sample size $250$ and resolution
$25\times 25$ with varying smoothness parameter $\nu$. In (b), results are for
sample size $500$ and $\nu=0.01$ with varying resolution. The numbers are
averages based on $25$ simulation runs.

In Figure 4, we show the results for the first four examples (Ex 1–4) in 2D
with fixed sample size $N=500$ and varying resolution $K=5,10,20,40,80$. The
Brownian sheet and the integrated Brownian sheet examples (Ex 1 and 2) are
separable. But even for these examples, the proposed CovNet estimators perform
better than the best separable estimator (Fig. 4(a) and (c)), especially in
low resolutions. This shows the ability of the CovNet model to learn the
underlying pattern, even when we observe the fields at a rather small number
of locations. For the rotated examples, we chose the matrix $\mathrm{O}$ to be
the $45^{\circ}$-rotation along the $x$-axis:

$\mathrm{O}=\begin{pmatrix}1/\sqrt{2}&-1/\sqrt{2}\\\
1/\sqrt{2}&1/\sqrt{2}\end{pmatrix}.$

The absence of separability of $\mathcal{C}$ has dire consequence on the
performance of the best separable estimator. The other estimators are
seemingly unaffected by this, and the CovNet estimators outperform the
empirical estimator. Among the CovNet estimators, the deepshared variant
performed much better than the others. Recall that the integrated Brownian
sheet (both the usual and the rotated) is one order smoother than the
(corresponding version of) Brownian sheet. While this added smoothness
enhances the performance of the CovNet estimators, we see an opposite effect
on the other estimators, especially with small resolutions. The difference
between the performance of the different CovNet models is also lesser in the
smoother examples.

In Figure 5, we show the results for the Matérn covariance model (Ex 5) in 2D.
We consider two different setups. In panel (a), we show the results with
$N=250$ and $K=25$ with varying $\nu$. Here, the empirical covariance performs
very poorly, especially when the surfaces are rougher (i.e., for smaller
values of $\nu$). The CovNet estimators perform better than the best separable
estimator when $\nu$ is small. When $\nu$ is large, i.e., the surfaces are
smoother, the errors of the best separable estimator is almost
indistinguishable from those of the CovNet estimators. However, same relative
error does not mean that the estimators share the same characteristic. In
fact, in this example, the CovNet estimators have an advantage over the other
estimators, which is evident from the eigendecomposition of the estimators
(see Figure 8). Detailed discussion on this is given in Section 5.2. In panel
(b), we report the results for $N=500$ and $\nu=0.01$ with varying resolution
$K$. Here, we again see the superiority of the CovNet estimators, especially
when the resolution is low.

A few words are in order about the cost of storage and manipulation of the
estimators. The storage of the empirical covariance estimator becomes
prohibitive rather quickly. Although we can compute the estimation error of
the empirical covariance relatively easily, it is very costly to manipulate
it, e.g., by inverting, for further applications like kriging. The scenario is
much better for the best separable estimator. But, both the empirical and the
best separable estimators produce a discretized object. Thus, even to evaluate
the estimated covariance at a location outside of the observation grid, one
needs to interpolate or smooth the estimated covariance. Depending on the
smoother used, this can dramatically increase the cost associated with the
estimator. The functional form of the CovNet estimators, on the other hand, do
not suffer from such problems. After estimation, the storage of the model is
quite cheap – one only needs to store the matrices and vectors associated with
the neural network model, which can be done very efficiently. Moreover, using
the eigendecomposition methods discussed in Section 4, we can easily
manipulate the fitted model.

###  5.1 Choice of hyperparameters

The performance of the proposed method depends on the choice of
hyperparameters, namely the number of components $R$ and the depth of the
network $L$ (for deep and deepshared models). Thus, it is important to select
these hyperparameters from the data, which is quite challenging for neural
networks (Bengio, 2012). We can use $V$-fold cross-validation for this
purpose, where we split the data into $V$ parts. One of these $V$ parts is
used as the _validation set_ and the rest are used as the _training set_. For
a particular choice of hyperparameters, the training set is used to fit the
model, and its performance is evaluated on the validation set. This procedure
is repeated for all the $V$ parts to get the average cross-validation score
for a particular set of hyperparameters. Finally, we select the set of
hyperparameters that admit the smallest average cross-validation score.

The alternative formulation of the loss function (cf. Section 3) is again
useful for the cross-validation. Suppose that our data is split as
$\mathcal{X}^{\rm tr}_{1},\ldots,\mathcal{X}^{\rm tr}_{N_{1}}$ and
$\mathcal{X}^{\rm va}_{1},\ldots,\mathcal{X}^{\rm va}_{N_{2}}$, constituting
the training and the validation sets, respectively. The model is fitted on the
training set to produce the estimate $\widehat{\mathcal{G}}^{\rm tr}$. We
evaluate the performance of the model on the validation set by computing the
loss ${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
va}-\widehat{\mathcal{G}}^{\rm tr}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}$,
where $\widehat{\mathcal{C}}^{\rm va}$ is the empirical covariance based on
the validation set. Recall that by construction, both
$\widehat{\mathcal{G}}^{\rm tr}$ and $\widehat{\mathcal{C}}^{\rm va}$ are of
the form

$\widehat{\mathcal{G}}^{\rm
tr}=\frac{1}{N_{1}}\sum_{n=1}^{N_{1}}\mathcal{X}^{\rm
tr,NN}_{n}\otimes\mathcal{X}^{\rm tr,NN}_{n}~{}~{}\text{ and
}~{}~{}\widehat{\mathcal{C}}^{\rm
va}=\frac{1}{N_{2}}\sum_{n=1}^{N_{2}}\mathcal{X}^{\rm
va}_{n}\otimes\mathcal{X}^{\rm va}_{n},$

where $\mathcal{X}^{\rm tr,NN}_{1},\ldots,\mathcal{X}^{\rm tr,NN}_{N_{1}}$ are
the neural networks fitted to the training sample (see Section 3). Here, we
have assumed w.l.o.g. that the observations are centered. Now, it is easy to
see that the loss has the explicit form

${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
va}-\widehat{\mathcal{G}}^{\rm
tr}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}=\frac{1}{N_{1}^{2}}\sum_{n=1}^{N_{1}}\sum_{m=1}^{N_{1}}\langle\mathcal{X}_{n}^{\rm
tr,NN},\mathcal{X}_{m}^{\rm
tr,NN}\rangle^{2}+\frac{1}{N_{2}^{2}}\sum_{n=1}^{N_{2}}\sum_{m=1}^{N_{2}}\langle\mathcal{X}_{n}^{\rm
va},\mathcal{X}_{m}^{\rm
va}\rangle^{2}-\frac{2}{N_{1}N_{2}}\sum_{n=1}^{N_{1}}\sum_{m=1}^{N_{2}}\langle\mathcal{X}_{n}^{\rm
tr,NN},\mathcal{X}_{m}^{\rm va}\rangle^{2},$

which depends only on the inner-products. Thus, we can compute the loss
efficiently, without forming the high-order covariances
$\widehat{\mathcal{C}}^{\rm va}$ or $\widehat{\mathcal{G}}^{\rm tr}$.

Table 1: Relative errors (in %) of the CovNet models with hyperparameters chosen using $5$-fold cross-validation. Difference from the least observed error over the range of hyperparameters is shown in parentheses. Relative errors for the empirical and the best separable estimators are also reported. The reported numbers are based on one simulation run with $N=500$ and $K=25$. Example | Empirical | Best separable | Shallow | Deep | Deepshared
---|---|---|---|---|---
Brownian sheet 2D | $9.58$ | $9.22$ | $9.26\,(0.41)$ | $7.93\,(0.24)$ | $8.72\,(0.35)$
Rotated Brownian sheet 2D | $11.79$ | $65.99$ | $10.03\,(0.64)$ | $10.36\,(0.69)$ | $9.70\,(0.17)$
Integrated Brownian sheet 2D | $7.99$ | $7.98$ | $7.34\,(0.06)$ | $9.58\,(3.37)$ | $7.44\,(0.09)$
Rotated integrated Brownian sheet 2D | $11.02$ | $70.53$ | $6.93\,(0.00)$ | $6.69\,(0.04)$ | $6.15\,(0.36)$
Matern 2D $\nu=0.001$ | $51.74$ | $17.65$ | $12.47\,(0.17)$ | $12.50\,(1.43)$ | $13.11\,(1.74)$
Matern 2D $\nu=0.01$ | $51.52$ | $18.40$ | $14.25\,(0.62)$ | $13.95\,(2.56)$ | $13.05\,(0.00)$
Matern 2D $\nu=0.1$ | $17.00$ | $11.68$ | $12.25\,(0.57)$ | $11.55\,(0.19)$ | $11.38\,(0.37)$
Matern 2D $\nu=1$ | $8.23$ | $8.23$ | $8.14\,(0.00)$ | $8.73\,(0.63)$ | $8.00\,(0.00)$

The results for the proposed cross-validation strategy are shown in Table 1
for the examples in 2D. Note that cross-validation is time consuming, and a
complete simulation study with cross-validation is rather difficult. So, for
each of the examples considered in the previous section, we report relative
errors for the three CovNet models (shallow, deep and deepshared) selected via
cross-validation based on a single simulation run with $500$ observations at a
resolution of $25\times 25$. For each model, we also show the difference from
the least observed relative error over the range of hyperparameters. The
relative errors for the empirical and the best separable estimators are also
reported to facilitate comparison. In all the examples, the average difference
from the best result was less than $1\%$ for all the CovNet models, while the
maximum difference was less than $1.75\%$ for the shallow and the deepshared
models and less than $3.5\%$ for the deep CovNet model. These results clearly
show that the cross-validation method can identify a good set of
hyperparameters in practice.

###  5.2 Estimated eigenstructure

Here, we demonstrate the usefulness of the eigendecomposition of the CovNet
estimators. For this purpose, we consider three examples in 2D, the rotated
Brownian sheet (Ex 2), the rotated integrated Brownian sheet (Ex 4) and the
Matérn covariance (Ex 5) with $\nu=0.01$. For each example, we plot the
eigensurfaces of the CovNet estimators obtained using the method proposed in
Section 4. The reported results are based on $500$ samples. For the first two
examples, we used a resolution of $10\times 10$. The Matérn example with
$\nu=0.01$ is much more rough, and a resolution of $10\times 10$ was too low
for all the methods (see Figure 5(b)). So, for this example, we used a
resolution of $25\times 25$. For comparison, we have also plotted the leading
eigensurfaces of the true covariance $\mathcal{C}$, the empirical estimator
and the best separable estimator. For the CovNet estimators, we selected the
hyperparameters via the cross-validation strategy described in the previous
section.

Figure 6: First six eigensurfaces of the true covariance and different
covariance estimators in the rotated Brownian sheet example with $N=500$ and
resolution $10\times 10$. For the CovNet estimators, the eigenfunctions are
computed using the methods described in Section 4. The plots are heatmaps of
the 2D eigenfunctions.

In Figure 6, we plot the eigensurfaces for the rotated Brownian sheet. Note
that in this case, the eigensurfaces of order seven and beyond explain less
than $1\%$ of the total variation of true covariance. Hence, we report the
first six eigensurfaces for each covariance (the truth and the estimators).
The usefulness of the CovNet estimators is quite evident from these plots.
Here, for each surface, we make observations at $100$ ($10\times 10$)
locations. Clearly, this is not enough for the empirical or the best separable
estimators. In contrast, the shallow and the deepshared CovNet estimators are
able to extract the features of the true covariance.

Figure 7:  First four eigensurfaces of the true covariance and different
covariance estimators in the rotated integrated Brownian sheet example with
$N=500$ and resolution $10\times 10$. For the CovNet estimators, the
eigenfunctions are computed using the methods described in Section 4. The
plots are heatmaps of the 2D eigenfunctions.

The results for the rotated integrated Brownian sheet are shown in Figure 7.
We plot the top four eigensurfaces as the other explain less than $1\%$ of the
total variation. The true covariance is quite smooth in this example, as a
result the empirical covariance does a better job. Even then, the functional
form of the CovNet estimators gives them an edge, which is reflected in the
estimation errors.

Figure 8:  First six eigensurfaces of the true covariance and different
covariance estimators in the Matérn example with $\nu=0.01$, $N=500$, and
resolution $25\times 25$. For the CovNet estimators, the eigenfunctions are
computed using the methods described in Section 4. The plots are heatmaps of
the 2D eigenfunctions.

In Figure 8, we show the results for the Matérn covariance with $\nu=0.01$. In
this case, the underlying process has rough sample paths. This roughness of
the observations has a damaging effect on the performance of the empirical
covariance. And, although the estimation error of the best separable estimator
is relatively low, the estimated eigensurfaces have little resemblance with
the true eigensurfaces. In fact, they are not able to capture the underlying
features, and the roughness of the observations can be clearly seen to affect
the performance. The shallow and the deepshared CovNet estimators do an
excellent job in identifying the salient features, even from the rough
observations.

A few comments are in order for the deep CovNet estimator. In all three
examples, the deep CovNet estimator is seemingly unable to capture the true
eigenstructure. However, the estimation error for this model is rather low.
This is perhaps due to the high complexity of the model, which allows it to
approximate the covariance well enough. But, without any restriction, the
model apparently does not really learn _interesting patterns of variation_
from the data. The added restriction of the deepshared model (in terms of
weight sharing) resolves this problem. The deepshared model is quite rich, but
at the same time it is able to extract interesting traits from the data.

We would also like to point out the favorable computational aspect of the
CovNet estimators in this context. As already mentioned, the
eigendecomposition for the CovNet estimators can be performed without forming
the covariance operators. This is sharply in contrast with the empirical
covariance, for which we need to apply eigendecomposition on a $K^{d}\times
K^{d}$-dimensional object, which can be prohibitive depending on $K$ and $d$.
The best separable estimator is seemingly immune to this problem. But, even
for this estimator, the computation for the eigendecomposition increases with
$K$ in the order of $\mathcal{O}(dK^{3})$ (eigendecomposition of $d$ matrices,
each of the order $K\times K$). The eigendecomposition for the CovNet
estimators, on the other hand, is completely free of $K$ (after estimation of
the model, of course). There is, however, a Monte-Carlo step involved in the
process. But, in all the examples, it took us only a few seconds to obtain the
eigendecomposition, which was much faster than the other two estimators.
Moreover, the functional form of the CovNet models have additional benefits,
as can be seen from the plots.

###  5.3 Application to fMRI data

To explore the usefulness of our methodology in a real setting, we consider
the fMRI data sets from the $1000$ Functional Connectomes Project. These data
sets are available at https://www.nitrc.org/projects/fcon_1000/, and consist
of resting state fMRI scans for more than $1200$ subjects collected at
different locations all over the world. For each subject, the data consist of
3D scans of the brain taken at a resolution of $64\times 64\times 33$ over
$225$ time points. These data sets were previously analyzed by Aston and Kirch
(2012); Stoehr, Aston and Kirch (2021) in the FDA setting, where they checked
for the stationarity of the 3D MRI scans over time for each of $197$
individuals from Beijing, China. For our demonstration, we considered sub69518
from Beijing, which was identified to exhibit stationarity by Aston and Kirch
(2012). This gives us $225$ 3D scans on a grid of size $64\times 64\times 33$,
which we treat as i.i.d. observations. We preprocessed the data set by
removing a polynomial trend of order $3$ from each voxel as suggested by Aston
and Kirch (2012). Further, we scaled the scans to have voxel-wise unit
variance before applying the methodology.

(a) Deepshared CovNet
---
(b) Separable Model
Figure 9: The leading eigenfunction of (a) the deepshared CovNet and (b) the
separable model, fitted to the 3D fMRI data. The hyperparameters of the
deepshared CovNet (depth$=6$, $R=30$) were selected via cross-validation
described in Section 5.1. The plots are heatmaps of 2D slices of the 3D
eigenfunctions, where the slices are taken over the z-axis.

As already mentioned in the introduction, covariance estimation is one of the
most important problems for resting state fMRI data as it enables to
understand the connectivity patterns of the brain. At the same time, the high-
dimensionality of the problem makes it extremely difficult to achieve. In
contrast, the CovNet models can be fitted efficiently to this data. In
particular, we applied the deepshared CovNet model to this data owing to its
superiority over the other CovNet models in the simulations. We selected the
hyperparameters via the cross-validation strategy described in Section 5.1. In
Figure 9(a), we show the leading eigenfunction of the fitted CovNet model. For
comparison, we also fitted a separable covariance model. But, as pointed out
by Aston and Kirch (2012), finding the best separable approximation is also
difficult in 3D. So, we used the separable estimator via marginalization as
proposed by the authors. The leading eigenfunction of the separable estimator
is shown in Figure 9(b).

(a) Second leading eigenfunction
---
(b) Third leading eigenfunction
Figure 10: (a) The second and (b) the third leading eigenfunctions of the
deepshared covariance network fitted to the 3D fMRI data. The second and the
third eigenfunctions explain $26.85\%$ and $19.05\%$ variability,
respectively. The plots are heatmaps of 2D slices of the 3D eigenfunctions,
where the slices are taken over the z-axis.

The plots clearly show that the deepshared CovNet is able to capture a much
richer structure than the separable model. While the leading eigenfunction of
the CovNet model accounts for $42.74\%$ of the variability explained, the same
for the separable model is only $0.0004\%$. In fact, to explain $42.74\%$
variability, we would need $45538$ components in the separable model.
Moreover, the fitted CovNet is able to extract the natural symmetry in the
brain during resting state (Raemaekers et al., 2018), while no such
constraints were imposed a priori.

In Figure 10, we show the second and the third leading eigenfunctions of the
fitted CovNet which, together with the leading eigenfunction, account for
$88.64\%$ of the variability. In contrast, for the separable model, we would
require $84578$ components to explain $88.64\%$ variability. The plots clearly
exhibit the ability of the deepshared CovNet to identify complex structures
from the data.

##  6 Asymptotic theory

We conclude the article by developing asymptotic theory for our estimators. In
particular, we prove that the CovNet estimators are consistent, and derive
their rates of convergence. We will consider two different setups – (a) when
the fields are fully observed and (b) when the fields are observed on a grid
with possible noise contamination.

###  6.1 The case of fully observed fields

Here, our data consist of i.i.d. random fields
$\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ distributed as $\mathcal{X}$. For
convenience, we start by assuming that the random field $\mathcal{X}$ is
bounded, i.e., there exists $\beta_{N}>0$ such that
$\|\mathcal{X}\|^{2}\leq\beta_{N}$ almost surely. This type of boundedness
assumption is quite common in the theoretical analysis of neural networks
(e.g., Györfi et al., 2002; Schmidt-Hieber, 2020). Note that the bound is
allowed to grow with $N$, so this can also be seen as a growing truncation
level. We will eventually remove this boundedness condition.

We also need to impose some condition on the approximating class, for which we
restrict the eigenstructure of the matrix $\Lambda=((\lambda_{r,s}))$ in
(3.1). Specifically, for a constant $\lambda_{N}>0$, we enforce that
$\Lambda\preceq\lambda_{N}\mathrm{I}_{R}$. Note that for any non-negative
definite matrix $\Lambda$, we always have
$\Lambda\preceq\lambda_{\max}\mathrm{I}_{R}$, where $\lambda_{\max}$ is the
largest eigenvalue of $\Lambda$. Thus, our assumption imposes a restriction on
the largest eigenvalue of the matrix $\Lambda$ in (3.1). We write
$\mathscr{F}_{R,\lambda_{N}}$ to denote the restricted class of kernels
corresponding to the generic class $\mathscr{F}_{R}$. In particular,

$\displaystyle\mathscr{F}^{\rm sh}_{R,\lambda_{N}}$
$\displaystyle=\big{\\{}c_{\rm sh}\text{ of the form
\eqref{eq:shallow_covnet_kernel}}:\mathrm{0}\preceq\Lambda=((\lambda_{r,s}))\preceq\lambda_{N}{\rm
I}_{R},\mathbf{w}_{r}\in\mathbb{R}^{d},b_{r}\in\mathbb{R}\big{\\}},$
$\displaystyle\mathscr{F}^{\rm d}_{R,L,\lambda_{N}}$
$\displaystyle=\big{\\{}c_{\rm d}\text{ of the form
\eqref{eq:deep_covnet_kernel}}:\mathrm{0}\preceq\Lambda:=((\lambda_{r,s}))\preceq\lambda_{N}\,\mathrm{I}_{R},g_{1},\ldots,g_{R}\in\mathcal{D}_{L,R}\big{\\}}\text{
and }$ $\displaystyle\mathscr{F}^{\rm ds}_{R,L,\lambda_{N}}$
$\displaystyle=\big{\\{}c_{\rm ds}\text{ of the form
\eqref{eq:deepshared_covnet_kernel}}:\mathrm{0}\preceq\Lambda:=((\lambda_{r,s}))\preceq\lambda_{N}\mathrm{I}_{R},$
$\displaystyle\kern 129.16626ptg_{1},\ldots,g_{R}\text{ of the form
\eqref{eq:deepshared_neural_network} with }p_{1}=\cdots=p_{L}=R\big{\\}},$
(6.1)

denote the restricted shallow, deep, and deepshared CovNet classes. In the
above, we write $\mathcal{D}_{L,R}$ to denote the deep neural network class
$\mathcal{D}_{L,\mathbf{p}}$ when $p_{1}=\cdots=p_{L}=R$. We also denote the
class of operators corresponding to the class of kernels
$\mathscr{F}_{R,\lambda_{N}}$ by $\widetilde{\mathscr{F}}_{R,\lambda_{N}}$.
Recall that our estimators are defined as

$\widehat{\mathcal{C}}^{\rm
sh}_{R,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
sh}_{R,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2},~{}\widehat{\mathcal{C}}^{\rm
d}_{R,L,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
d}_{R,L,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}~{}\text{
and }~{}\widehat{\mathcal{C}}^{\rm
ds}_{R,L,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
ds}_{R,L,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}.$
(6.2)

In the following, we start by proving two different kinds of results. First,
we prove consistency of the estimators under appropriate conditions, and then
derive their rates of convergence.

######  Theorem 3.

Let $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\overset{{\rm
i.i.d.}}{\sim}\mathcal{X}$, where $\mathcal{X}$ takes values in
$\mathcal{L}_{2}(\mathcal{Q})$, and $\mathcal{Q}$ is a compact subset of
$\mathbb{R}^{d}$. Also assume that $\|\mathcal{X}\|^{2}\leq\beta_{N}$ almost
surely, $\mathbb{E}(\mathcal{X})=0$ and
$\mathrm{Cov}(\mathcal{X})=\mathcal{C}$. Let $\widehat{\mathcal{C}}_{R,N}^{\rm
sh}$, $\widehat{\mathcal{C}}_{R,L,N}^{\rm d}$ and
$\widehat{\mathcal{C}}_{R,L,N}^{\rm ds}$ be the shallow, the deep, and the
deepshared CovNet estimators given by (6.2). Suppose that
$R\to\infty,\lambda_{N}\to\infty$ as $N\to\infty$, and define
$\Delta_{N}=\max\\{\beta_{N},|\mathcal{Q}|R\lambda_{N}\\}$.

  1. (A)

If $dR^{2}\Delta_{N}^{4}\log(\Delta_{N})/N\to 0$ as $N\to\infty$, then the
shallow CovNet estimator is weakly consistent for $\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
sh}_{R,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{P}{\to}0$.
Additionally, if $\Delta_{N}^{4}/N^{1-\delta}\to 0$ for some $\delta\in(0,1)$,
then the estimator is strongly consistent for $\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
sh}_{R,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{a.s.}{\to}0$
as $N\to\infty$.

  2. (B)

Let $R>d$. If $L^{4}R^{8}\Delta_{N}^{4}\log^{2}(L\Delta_{N})/N\to 0$ as
$N\to\infty$, then the deep CovNet estimator is weakly consistent for
$\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
d}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{P}{\to}0$.
Additionally, if
$L^{4}R^{8}\Delta_{N}^{4}\log^{2}(L\Delta_{N})/N^{1-\delta}\to 0$ for some
$\delta\in(0,1)$, then the estimator is strongly consistent for $\mathcal{C}$,
i.e., ${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
ds}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{a.s.}{\to}0$
as $N\to\infty$.

  3. (C)

Let $R>d$. If $L^{4}R^{6}\Delta_{N}^{4}\log^{2}(L\Delta_{N})/N\to 0$, then the
deepshared CovNet estimator is weakly consistent for $\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
ds}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{P}{\to}0$.
Additionally, if
$L^{4}R^{6}\Delta_{N}^{4}\log^{2}(L\Delta_{N})/N^{1-\delta}\to 0$ for some
$\delta\in(0,1)$, then the estimator is strongly consistent for $\mathcal{C}$,
i.e., ${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
ds}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{a.s.}{\to}0$
as $N\to\infty$.

The proof of the theorem involves _bias-variance-type decompositions_ for the
estimation error
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}$.
To control the bias term, we need the universal approximation property
(Theorems 1 and 2), but now with the additional restriction on the classes
(6.1). This is ensured by assuming that $R$, $L$ and $\lambda_{N}$ go to
infinity as $N$ diverges (see Remark 11). On the other hand,
$dR^{2}\Delta_{N}^{4}\log(\Delta_{N})/N$,
$L^{4}R^{8}\Delta_{N}^{4}\log^{2}(L\Delta_{N})/N$, and
$L^{4}R^{6}\Delta_{N}^{4}\log^{2}(L\Delta_{N})/N$ are linked to the _variance_
of the estimators. The conditions of the theorem ensures that the variance
also converges to $0$ with the sample size for the different estimators.

######  Remark 8.

There is a small technical ambiguity in the statement of Theorem 3. The
theorem is stated for i.i.d. observations distributed as $\mathcal{X}$ when
$N$ goes to infinity, whereas the bound on $\mathcal{X}$ is also allowed to
evolve with $N$. Thus, the result is to be understood for a triangular
sequence of arrays where, for each $N$, the observations are bounded by a
constant, which in turn is allowed to diverge keeping up with the assumption
of the theorem. However, we avoid stating the theorem in this generality for
ease of exposition. The special case of i.i.d. observations (i.e., when
$\beta_{N}$ is a constant) follows easily from the theorem.

Next, we derive the rate of convergence of the estimators.

######  Theorem 4.

Let $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\overset{{\rm
i.i.d.}}{\sim}\mathcal{X}$, where $\mathcal{X}$ takes values in
$\mathcal{L}_{2}(\mathcal{Q})$, and $\mathcal{Q}$ is a compact subset of
$\mathbb{R}^{d}$. Suppose that $\|\mathcal{X}\|^{2}\leq\beta_{N}$ almost
surely, $\mathbb{E}(\mathcal{X})=0$ and
$\mathrm{Cov}(\mathcal{X})=\mathcal{C}$. Let $\widehat{\mathcal{C}}_{R,N}^{\rm
sh}$, $\widehat{\mathcal{C}}_{R,L,N}^{\rm d}$ and
$\widehat{\mathcal{C}}_{R,L,N}^{\rm ds}$ be the shallow, the deep, and the
deepshared CovNet estimators given by (6.2). Then,

$\displaystyle\mathbb{E}\Big{(}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
sh}_{R,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\Big{)}$
$\displaystyle\leq 2\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
sh}_{R,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}+\mathcal{O}\bigg{(}\frac{dR^{2}\Delta_{N}^{4}\log(N)}{N}\bigg{)},$
$\displaystyle\mathbb{E}\Big{(}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
d}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\Big{)}$
$\displaystyle\leq 2\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
d}_{R,L,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}+\mathcal{O}\bigg{(}\frac{L^{4}R^{8}\Delta_{N}^{4}\log^{2}(N)}{N}\bigg{)}\text{,
and }$
$\displaystyle\mathbb{E}\Big{(}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{\rm
ds}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\Big{)}$
$\displaystyle\leq 2\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
ds}_{R,L,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}+\mathcal{O}\bigg{(}\frac{L^{4}R^{6}\Delta_{N}^{4}\log^{2}(N)}{N}\bigg{)}.$

Here, $\Delta_{N}=\max\\{\beta_{N},|\mathcal{Q}|R\lambda_{N}\\}$ is as defined
in Theorem 3. In the above, for the deep and the deepshared CovNet estimators,
we have assumed that $R>d$.

The theorem clearly shows the bias-variance-type decomposition for the
proposed estimators. To get the exact rates of convergence, we need to
quantify the bias terms, which is an approximation theoretic problem. If, for
example, the bias term is zero for some finite $R,L$ and $\lambda_{N}$, then
the derived rate of convergence is the same as that of the empirical
estimator, except for the logarithmic term. Thus, in this case, our estimator
enjoys a nearly minimax rate of convergence. In general, to get the rate of
convergence of the bias, we need to make further assumptions. There are two
ways of doing this, either by making assumptions on the eigenstructure of the
true covariance or by making assumptions on the smoothness of the underlying
field $\mathcal{X}$ (see Appendix B.2 for details). For this line of
derivations, the rates depend crucially on the approximation error of the
constituents of the CovNet model under consideration. For instance, for the
shallow CovNet, if we assume that the underlying field $\mathcal{X}$ takes
values in $\mathcal{S}^{\alpha}(\mathcal{Q})$, the Sobolev space of functions
of order $\alpha$ on $\mathcal{Q}$ (see Mhaskar, 1996), with almost surely
bounded Sobolev norm (i.e.,
$\|\mathcal{X}\|_{\mathcal{S}^{\alpha}(\mathcal{Q})}^{2}\leq\beta$ a.s.), then
by selecting $\lambda_{N}\asymp R$, we can bound the bias term as
$\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
sh}_{R,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\lesssim
R^{-\alpha/d}$ (see (B.11)). So, for consistency we need
$R^{10}={\scriptstyle{\mathcal{O}}}(N/(d\log(N))$, while the optimal rate is
achieved for $R\asymp(N/d\log(N))^{d/(10d+\alpha)}$. This leads to the rate of
convergence $\mathcal{O}\big{(}(d\log(N)/N)^{\alpha/(10d+\alpha)}\big{)}$.
Similarly, one can use results from Langer (2021); Ohn and Kim (2019) to bound
the bias of the deep and the deepshared CovNet models.

######  Remark 9.

Both Theorems 3 and 4 are derived here for a global minimizer of the loss
function. But, in practice, we are not guaranteed to find a global minima.
Following the proof of Theorem 3, it can be shown that the consistency results
hold as long as the estimator is within ${\scriptstyle{\mathcal{O}}}_{P}(1)$
of the minimizer, i.e., for estimators $\widehat{\mathcal{C}}$ satisfying
${|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\widehat{\mathcal{C}}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\leq\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}+{\scriptstyle{\mathcal{O}}}_{P}(1)$.
Similarly, from the proof of Theorem 4, one can check that for an approximate
minimizer, the rate of convergence gets inflated by the “expected minimization
gap”
$\mathbb{E}\Big{(}{|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}_{N}-\widehat{\mathcal{C}}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}-\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,\lambda_{N}}}{|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\,\big{|}\,\mathscr{X}_{N}\Big{)}$,
where $\mathscr{X}_{N}=\\{\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\\}$ denotes
the observed data, and $\widetilde{C}_{N}$ is distributed identically to
$\widehat{\mathcal{C}}_{N}$ but independently of $\mathscr{X}_{N}$.

Finally, we prove consistency without the boundedness condition on
$\mathcal{X}$. In this case, we need to slightly modify our estimators. To
this extent, let $\mathcal{G}$ be a CovNet operator from the unrestricted
class $\widetilde{\mathscr{F}}_{R}$ with kernel $g(\text{\bf u},\text{\bf
v})=\sum_{r=1}^{R}\sum_{s=1}^{R}\lambda_{r,s}\,g_{r}(\text{\bf
u})\,g_{s}(\text{\bf v})$. For $\lambda_{N}>0$, define
$\mathscr{P}_{\lambda_{N}}\mathcal{G}$ to be the CovNet operator obtained by
thresholding the eigenvalues of $\Lambda:=((\lambda_{r,s}))$ to $\lambda_{N}$.
To be precise, if
$\Lambda=\sum_{i=1}^{R}\eta_{i}\,\mathbf{e}_{i}\,\mathbf{e}_{i}^{\top}$ is the
eigendecomposition of $\Lambda$, then we define
$\Lambda_{\lambda_{N}}=\sum_{i=1}^{R}\min\\{\eta_{i},\lambda_{N}\\}\,\mathbf{e}_{i}\,\mathbf{e}_{i}^{\top}$
to be the $\lambda_{N}$-thresholded version of $\Lambda$. We define
$\mathscr{P}_{\lambda_{N}}\mathcal{G}$ to be the operator with kernel
$g_{\lambda_{N}}(\text{\bf u},\text{\bf
v})=\sum_{r=1}^{R}\sum_{s=1}^{R}\widetilde{\lambda}_{r,s}\,g_{r}(\text{\bf
u})\,g_{s}(\text{\bf v})$, where $\widetilde{\lambda}_{r,s}$ is the $(r,s)$-th
element of the matrix $\Lambda_{\lambda_{N}}$. By construction,
$\mathrm{0}\preceq\Lambda_{\lambda_{N}}\preceq\lambda_{N}\,\mathrm{I}_{R}$ and
consequently, $\mathscr{P}_{\lambda_{N}}\mathcal{G}$ is an element of the
restricted class $\widetilde{\mathscr{F}}_{R,\lambda_{N}}$. We are now ready
to re-define the estimator. Define

$\widehat{\mathcal{C}}^{\rm
sh}_{R,N}\in\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
sh}_{R,\sigma}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2},\quad\widehat{\mathcal{C}}^{\rm
d}_{R,L,N}\in\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
d}_{R,L,\sigma}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\quad\text{
and }\quad\widehat{\mathcal{C}}^{\rm
ds}_{R,L,N}\in\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}^{\rm
ds}_{R,L,\sigma}}{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}_{N}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2},$

to be the shallow, the deep and the deepshared CovNet estimators,
respectively, but without any restriction on the underlying classes. Now, for
a constant $\lambda_{N}>0$, we define our modified estimators as

$\widetilde{\mathcal{C}}^{\rm
sh}_{R,N}=\mathscr{P}_{\lambda_{N}}\widehat{\mathcal{C}}^{\rm
sh}_{R,N},\quad\widetilde{\mathcal{C}}^{\rm
d}_{R,L,N}=\mathscr{P}_{\lambda_{N}}\widehat{\mathcal{C}}^{\rm
d}_{R,L,N}\quad\text{ and }\quad\widetilde{\mathcal{C}}^{\rm
ds}_{R,L,N}=\mathscr{P}_{\lambda_{N}}\widehat{\mathcal{C}}^{\rm ds}_{R,L,N}.$
(6.3)

These modified estimators are consistent, as shown in the following theorem.

######  Theorem 5.

Let $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\overset{{\rm
i.i.d.}}{\sim}\mathcal{X}$, where $\mathcal{X}$ takes values in
$\mathcal{L}_{2}(\mathcal{Q})$ with $\mathbb{E}(\|\mathcal{X}\|^{4})<\infty$.
Also assume that $\mathbb{E}(\mathcal{X})=0$ and
$\mathrm{Cov}(\mathcal{X})=\mathcal{C}$. Let
$\widetilde{\mathcal{C}}_{R,N}^{\rm sh}$,
$\widetilde{\mathcal{C}}_{R,L,N}^{\rm d}$ and
$\widetilde{\mathcal{C}}_{R,L,N}^{\rm ds}$ be the modified shallow, deep, and
deepshared CovNet estimators given by (6.3). Assume that $R>d$, and
$R\to\infty,\lambda_{N}\to\infty$ as $N\to\infty$.

  1. (A)

If $dR^{6}\lambda_{N}^{4}\log(R\lambda_{N})/N\to 0$ as $N\to\infty$, then the
modified shallow CovNet estimator is weakly consistent for $\mathcal{C}$,
i.e., ${|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}^{\rm
sh}_{R,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{P}{\to}0$.
Additionally, if $dR^{6}\lambda_{N}^{4}\log(R\lambda_{N})/N^{1-\delta}\to 0$
for some $\delta\in(0,1)$, then it is strongly consistent for $\mathcal{C}$,
i.e., ${|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}^{\rm
sh}_{R,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{a.s.}{\to}0$
as $N\to\infty$.

  2. (B)

Let $R>d$. If $L^{4}R^{12}\lambda_{N}^{4}\log^{2}(LR\lambda_{N})/N\to 0$ as
$N\to\infty$, then the modified deep CovNet estimator is weakly consistent for
$\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}^{\rm
d}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{P}{\to}0$.
Additionally, if
$L^{4}R^{12}\lambda_{N}^{4}\log^{2}(LR\lambda_{N})/N^{1-\delta}\to 0$ for some
$\delta\in(0,1)$, then it is strongly consistent for $\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}^{\rm
d}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{a.s.}{\to}0$
as $N\to\infty$.

  3. (C)

Let $R>d$. If $L^{4}R^{10}\lambda_{N}^{4}\log^{2}(LR\lambda_{N})/N\to 0$, then
the modified deepshared CovNet estimator is weakly consistent for
$\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}^{\rm
ds}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{P}{\to}0$.
Additionally, if
$L^{4}R^{10}\lambda_{N}^{4}\log^{2}(LR\lambda_{N})/N^{1-\delta}\to 0$ for some
$\delta\in(0,1)$, then it is strongly consistent for $\mathcal{C}$, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}^{\rm
ds}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}\overset{a.s.}{\to}0$
as $N\to\infty$.

######  Remark 10.

Our derived rates are truly nonparametric, with minimal assumptions on the
underlying structure. We assumed $\mathbb{E}(\|\mathcal{X}\|^{4})<\infty$,
which is standard for covariance estimation. Moreover, we made no assumption
on the underlying covariance operator $\mathcal{C}$. As a consequence, our
derived rates are rather slow in terms of the number of parameters of the
models. These can be improved by making further assumptions on the random
field $\mathcal{X}$ or the eigenfunctions of $\mathcal{C}$ (e.g., the ones
used by Bauer and Kohler (2019) or Schmidt-Hieber (2020) in the context of
nonparametric regression). However, such specialized treatments are beyond the
scope of the present article. If the rank of $\mathcal{C}$ is small, which is
very often the case in FDA, then a small $R$ is enough to control the bias
term (see Appendix B). On the other hand, such a small $R$ gives us a
considerable gain in terms of the variance, thus reducing the overall
estimation error. However, one should also note that the derived rates are
only upper bounds, and we do not claim tightness of the bounds.

###  6.2 The case of discretely observed fields

The results derived so far are for fully observed random fields. But in
practice, we observe the fields on a grid, with possible noise contamination.
Here, we develop asymptotic properties of our estimators in this scenario.
W.l.o.g., we assume that $\mathcal{Q}=[0,1]^{d}$, and we observe the data on a
$K_{1}\times\cdots\times K_{d}$ regular grid on $[0,1]^{d}$. To this extent,
let
$\\{T_{1,1}^{K_{1}},\ldots,T_{1,K_{1}}^{K_{1}}\\},\ldots,\\{T_{d,1}^{K_{d}},\ldots,T_{d,K_{d}}^{K_{d}}\\}$
be regular partitions of $[0,1]$ of sizes $K_{1},\ldots,K_{d}$, respectively.
Define $V_{i_{1},\ldots,i_{d}}^{K}=T_{1,i_{1}}^{K_{1}}\times\cdots\times
T_{d,i_{d}}^{K_{d}}$ to be the $(i_{1},\ldots,i_{d})$-th _voxel_ for $1\leq
i_{1}\leq K_{1},\ldots,1\leq i_{d}\leq K_{d}$. The voxels are non-overlapping
(i.e., $V_{i_{1},\ldots,i_{d}}^{K}\cap V_{j_{1},\ldots,j_{d}}^{K}=\emptyset$
for $(i_{1},\ldots,i_{d})\neq(j_{1},\ldots,j_{d})$), and they form a regular
partition of $[0,1]^{d}$. In particular,
$|V_{i_{1},\ldots,i_{d}}^{K}|=(\prod_{i=1}^{d}K_{i})^{-1}$. For each random
field $\mathcal{X}_{n}$, we make a single measurement at each of the voxels.
These measurements are assumed to be of the form

$\widetilde{X}_{n}^{K}[i_{1},\ldots,i_{d}]=X_{n}^{K}[i_{1},\ldots,i_{d}]+E_{n}^{K}[i_{1},\ldots,i_{d}],\qquad
1\leq i_{1}\leq K_{1},\ldots,1\leq i_{d}\leq K_{d},\,n=1,\ldots,N,$ (6.4)

where $X^{K}_{n}[i_{1},\ldots,i_{d}]$ is a discretization of $\mathcal{X}_{n}$
over the $(i_{1},\ldots,i_{d})$-th voxel and $E_{n}^{K}[i_{1},\ldots,i_{d}]$
is the corresponding measurement error or noise. We consider two different
measurement schemes which relate the discrete object
$\mathbf{X}_{n}^{K}=(X_{n}^{K}[i_{1},\ldots,i_{d}])$ to the respective field
$\mathcal{X}_{n}=(X_{n}(\text{\bf u}):\text{\bf u}\in[0,1]^{d})$.

  1. (M1)

Point-wise measurement:

$X_{n}^{K}[i_{1},\ldots,i_{d}]=X_{n}(u_{i_{1}},\ldots,u_{i_{d}}),\qquad 1\leq
i_{1}\leq K_{1},\ldots,1\leq i_{d}\leq K_{d},$

where $(u_{i_{1}},\ldots,u_{i_{d}})\in V_{i_{1},\ldots,i_{d}}^{K}$ is a
location within the $(i_{1},\ldots,i_{d})$-th voxel. For the measurements to
be meaningful, we need to assume that $\mathcal{X}$ has continuous sample
paths (e.g., Hsing and Eubank, 2015).

  2. (M2)

Voxel-wise average:

$X_{n}^{K}[i_{1},\ldots,i_{d}]=\frac{1}{|V_{i_{1},\ldots,i_{d}}^{K}|}\int_{V_{i_{1},\ldots,i_{d}}^{K}}X_{n}(\text{\bf
u})\,\mathrm{d}\text{\bf u},\qquad 1\leq i_{1}\leq K_{1},\ldots,1\leq
i_{d}\leq K_{d}.$

For the measurement errors $E_{n}^{K}[i_{1},\ldots,i_{d}]$, we assume that
they are i.i.d. with mean $0$ and variance $\sigma_{K}^{2}$, and are
uncorrelated with the $X_{n}^{K}[i_{1},\ldots,i_{d}]$’s. In line with our
previous assumptions, we also assume that
$\big{|}E_{n}^{K}[i_{1},\ldots,i_{d}]\big{|}^{2}\leq\beta^{\rm e}_{N,K}$
almost surely.

We denote the measurements corresponding to the $n$-th field $\mathcal{X}_{n}$
by
$\widetilde{\mathbf{X}}_{n}^{K}=(\widetilde{X}_{n}^{K}[i_{1},\ldots,i_{d}])$.
Define
$\widetilde{\mathbf{C}}_{N}^{K}=N^{-1}\sum_{n=1}^{N}\widetilde{\mathbf{X}}_{n}^{K}\otimes\widetilde{\mathbf{X}}_{n}^{K}$
to be the empirical covariance based on the discretely observed data. For a
generic class of CovNet operators $\widetilde{\mathscr{F}}_{R}$, our estimator
is given by

$\widehat{\mathcal{C}}^{K}_{R,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R}}\big{\|}\widetilde{\mathbf{C}}_{N}^{K}-\mathbf{G}^{K}\big{\|}_{\mathrm{F}}^{2},$

where $\|\cdot\|_{\mathrm{F}}$ is the Frobenius norm and $\mathbf{G}^{K}$ is
the discretization of the operator $\mathcal{G}$ over the voxels, defined as

$\mathbf{G}^{K}[i_{1},\ldots,i_{d};j_{1},\ldots,j_{d}]=g(v_{i_{1}},\ldots,v_{i_{d}};v_{j_{1}},\ldots,v_{j_{d}}),\qquad
1\leq i_{1},j_{1}\leq K_{1},\ldots,1\leq i_{d},j_{d}\leq K_{d},$

where $g$ is the kernel corresponding to $\mathcal{G}$ and
$(v_{i_{1}},\ldots,v_{i_{d}})$ is a location in the $(i_{1},\ldots,i_{d})$-th
voxel $V_{i_{1},\ldots,i_{d}}^{K}$. If we define
$\widetilde{\mathcal{C}}_{N}^{K}$ to be the voxel-wise continuation of
$\widetilde{\mathbf{C}}_{N}^{K}$, with kernel

$\widetilde{c}_{N}^{K}(\text{\bf u},\text{\bf
v})=\sum_{i_{1}=1}^{K_{1}}\cdots\sum_{i_{d}=1}^{K_{d}}\sum_{j_{1}=1}^{K_{1}}\cdots\sum_{j_{d}=1}^{K_{d}}\widetilde{C}_{N}^{K}[i_{1},\ldots,i_{d};j_{1},\ldots,j_{d}]\,\textrm{{1}}\\{\text{\bf
u}\in V_{i_{1},\ldots,i_{d}}^{K},\text{\bf v}\in
V_{j_{1},\ldots,j_{d}}^{K}\\},\qquad\text{\bf u},\text{\bf v}\in[0,1]^{d},$

then it is easy to see that

$\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R}}\big{\|}\widetilde{\mathbf{C}}_{N}^{K}-\mathbf{G}^{K}\big{\|}_{\mathrm{F}}^{2}\approx\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R}}{|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}_{N}^{K}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2},$

where the approximation holds when the resolution $K_{1}\times\cdots\times
K_{d}$ is large (see also (E.28) in Appendix E.4). We will derive theoretical
properties for this approximation. In particular, we define

$\displaystyle\widehat{\mathcal{C}}^{{\rm
sh},K}_{R,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,\lambda_{N}}^{\rm
sh}}{|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}_{N}^{K}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2},\,\widehat{\mathcal{C}}^{{\rm
d},K}_{R,L,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,L,\lambda_{N}}^{\rm
d}}{|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}_{N}^{K}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}\text{
and }\widehat{\mathcal{C}}^{{\rm
ds},K}_{R,L,N}\in\operatornamewithlimits{arg\,min}_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,L,\lambda_{N}}^{\rm
ds}}{|\kern-1.07639pt|\kern-1.07639pt|\widetilde{\mathcal{C}}_{N}^{K}-\mathcal{G}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2},$
(6.5)

to be the shallow, deep and deepshared CovNet estimators based on the discrete
measurements. The asymptotic behaviour of these estimators is established in
the following theorem.

######  Theorem 6.

Let $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\overset{{\rm
i.i.d.}}{\sim}\mathcal{X}$, where $\mathcal{X}$ takes values in
$\mathcal{L}_{2}([0,1]^{d})$ with $\mathbb{E}(\mathcal{X})=0$ and
$\mathrm{Cov}(\mathcal{X})=\mathcal{C}$. Let the kernel $c$ of $\mathcal{C}$
is Lipschitz on $[0,1]^{2d}$ with Lipschitz constant $\rho$. Consider the
measurement model (6.4), where the measurement errors
$E_{n}^{K}[i_{1},\ldots,i_{d}]$ are i.i.d. and uncorrelated with
$\mathbf{X}_{1}^{K},\ldots,\mathbf{X}_{N}^{K}$, and satisfy
$\big{|}E_{n}^{K}[i_{1},\ldots,i_{d}]\big{|}^{2}\leq\beta_{N,K}^{\rm e}$
almost surely, $\mathbb{E}(E_{n}^{K}[i_{1},\ldots,i_{d}])=0$,
$\textrm{Var}(E_{n}^{K}[i_{1},\ldots,i_{d}])=\sigma_{K}^{2}$. Suppose that one
of the following two hold.

  1. 1.

$\mathcal{X}$ has continuous sample paths,
$\|\mathcal{X}\|_{\infty}^{2}\leq\beta_{N}$ almost surely, and the
measurements $\mathbf{X}_{1}^{K},\ldots,\mathbf{X}_{N}^{K}$ are obtained from
$\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ via (M1).

  2. 2.

$\|\mathcal{X}\|^{2}\leq\beta_{N}$ almost surely and the measurements
$\mathbf{X}_{1}^{K},\ldots,\mathbf{X}_{N}^{K}$ are obtained from
$\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ via (M2).

Let $\widehat{\mathcal{C}}^{{\rm sh},K}_{R,N}$, $\widehat{\mathcal{C}}^{{\rm
d},K}_{R,L,N}$ and $\widehat{\mathcal{C}}^{{\rm ds},K}_{R,L,N}$ be the
shallow, deep and deepshared CovNet estimators given by (6.5). Then,

$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{{\rm
sh},K}_{R,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}$
$\displaystyle\leq
18\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,\lambda_{N}}^{\rm
sh}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}+\mathcal{O}\bigg{(}\frac{dR^{2}\Delta_{N,K}^{4}\log(N)}{N}\bigg{)}+a_{K},$
$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{{\rm
d},K}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}$
$\displaystyle\leq
18\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,L,\lambda_{N}}^{\rm
d}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}+\mathcal{O}\bigg{(}\frac{L^{4}R^{8}\Delta_{N,K}^{4}\log^{2}(N)}{N}\bigg{)}+a_{K},\quad\text{and}$
$\displaystyle{|\kern-1.07639pt|\kern-1.07639pt|\widehat{\mathcal{C}}^{{\rm
ds},K}_{R,L,N}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}$
$\displaystyle\leq
18\inf_{\mathcal{G}\in\widetilde{\mathscr{F}}_{R,L,\lambda_{N}}^{\rm
ds}}{|\kern-1.07639pt|\kern-1.07639pt|\mathcal{G}-\mathcal{C}|\kern-1.07639pt|\kern-1.07639pt|}_{2}^{2}+\mathcal{O}\bigg{(}\frac{L^{4}R^{6}\Delta_{N,K}^{4}\log^{2}(N)}{N}\bigg{)}+a_{K},$

where $\Delta_{N,K}=\max\\{2(\beta_{N}+\beta_{N,K}^{\rm e}),R\lambda_{N}\\}$
and
$a_{K}=21\rho^{2}\big{(}K_{1}^{-2}+\cdots+K_{d}^{-2}\big{)}+21\sigma_{K}^{4}/(K_{1}\cdots
K_{d})$.

The theorem clearly shows the effect of grid size and noise contamination on
the estimators. The rates are qualitatively the same as in Theorem 4, except
$\Delta_{N}$ is replaced by $\Delta_{N,K}$ and a couple of terms depending on
$(K_{1},\ldots,K_{d})$ are added. The term $\Delta_{N,K}$ can be viewed as a
noise contaminated version of $\Delta_{N}$. If we assume that
$\beta_{N,K}={\scriptstyle{\mathcal{O}}}(\beta_{N})$ (which can be seen as
assuming that the signal-to-noise ratio diverges), then $\Delta_{N,K}$ is
asymptotically equivalent to $\Delta_{N}$. In this case, we can see a clear
separation in the estimation error, one due to the sample size and the other
due to the resolution. Among the remaining two terms,
$\rho^{2}(K_{1}^{-2}+\cdots+K_{d}^{-2})$ arises due to the discretization of
the fields, while $\sigma_{K}^{4}/(K_{1}\cdots K_{d})$ is due to noise
contamination. Thus, for consistency of our estimators, we require
$(K_{1}^{-2}+\cdots+K_{d}^{-2})\to 0$, which is ensured if
$\min\\{K_{1},\ldots,K_{d}\\}\to\infty$. Moreover, the noise level
$\sigma_{K}$ is allowed to diverge, but at a slower rate than $(K_{1}\cdots
K_{d})^{1/4}$. The Lipschitz assumption on $c$ in the theorem is convenient,
but is by no means necessary. It can be verified that consistency of the
estimators holds as long as $c$ is continuous.

##  7 Concluding remarks

We have proposed three new classes of neural network models for covariance
estimation of functional data observed over multidimensional domains. The
advantages of the proposed models include efficient estimation, storage,
manipulation and performance guarantees. Our approach is motivated by the
demonstrated ability of neural networks in solving complex problems. And
indeed, our empirical studies show the superiority of the proposed methods,
especially the deepshared CovNet model. At the same time, our methods will
also be amenable to the shortcomings of neural networks, e.g., lack of
theoretical optimization guarantees for convergence to global minima. But, as
is the case with neural networks, despite these limitations, our experimental
results appear compelling. Any progress in the study of neural networks will,
in principle, translate to a commensurate progress in the understanding of
covariance networks.

Throughout the article, we have used the sigmoidal activation function. But,
most of the results, especially the ones for the deep CovNet models, can be
easily extended to include other activation functions, e.g., the ReLU. In some
preliminary numerical studies, we observed similar performance by the sigmoid
and the ReLU. We prefer the sigmoid because of the smoothness that it
provides, which is often beneficial for functional covariance estimation.

At the level of generality they are derived, our convergence rates are
arguably slow. But, these do not reveal the complete picture and are rather a
reflection of our completely nonparametric treatment of the problem. These
rates can be improved by considering more structured problems, which is now a
topic of interest in theoretical studies of neural networks (Bauer and Kohler,
2019; Schmidt-Hieber, 2020). Such additional structural assumptions may also
allow us to derive approximation errors for the models, which we have not
fully addressed here.

## Acknowledgment

We are grateful to Prof. Sir John A. D. Aston for providing access to the fMRI
data and enlightening us on some of their key aspects.

## Appendices

In these appendices, we give the proofs omitted from the main text, some
further mathematical details and additional simulation results. The
organization is as follows. In the next section, we provide some mathematical
background useful in subsequent developments. In Section B, we discuss the
bias of the CovNet models – we establish the universal approximation property,
and sketch two different ways to derive the rate of convergence of the bias
term. In particular, we derive the rate of convergence of the bias for the
shallow CovNet model. In Section C, we provide some details on the estimation
of the CovNet models and also describe a way to estimate the mean function
from the data using the CovNet models. In Section D, we treat the covering
numbers of certain spaces, as these play a fundamental role in the proofs of
the asymptotic results. In particular, we derive upper bounds on the covering
numbers of the classes of shallow, deep and deepshared CovNet operators. In
Section E, we provide proofs of the asymptotic results. Finally, in Section F,
we provide some additional numerical results which were left out in the main
text.

##  Appendix A Mathematical background

We start by summarizing some definitions and background concepts. More details
can be found in Hsing and Eubank (2015). Let
$\mathcal{Q}\subset\mathbb{R}^{d}$ be a compact set. We denote by
$\mathcal{L}_{2}(\mathcal{Q})$ the space of all real-valued square-integrable
functions on $\mathcal{Q}$. This is a Hilbert space when equipped with the
inner product $\langle f,g\rangle=\int_{\mathcal{Q}}f(\text{\bf u})g(\text{\bf
u})\,\mathrm{d}\text{\bf u}$ for $f,g\in\mathcal{L}_{2}(\mathcal{Q})$. A
linear map $\mathcal{A}$ from $\mathcal{L}_{2}(\mathcal{Q})$ onto itself is
called bounded if there exists a constant $M\geq 0$ such that
$\|\mathcal{A}f\|\leq M\|f\|$ for all $f\in\mathcal{L}_{2}(\mathcal{Q})$,
where $\|\cdot\|$ is the norm induced by the inner product
$\langle\cdot,\cdot\rangle$, i.e., $\|f\|=\sqrt{\langle f,f\rangle}$. A
bounded linear map is referred to as an _operator_. The minimum value of $M$
for which the boundedness condition holds is called the _operator norm_ and is
denoted by
${|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|}_{\infty}$.
An operator $\mathcal{A}$ is _compact_ if there exist orthonormal bases (ONBs)
$(\psi_{j})_{j\geq 1}$ and $(\phi_{j})_{j\geq 1}$ of
$\mathcal{L}_{2}(\mathcal{Q})$ such that $\mathcal{A}f=\sum_{j\geq
1}\lambda_{j}\langle\psi_{j},f\rangle\phi_{j}$. A compact operator
$\mathcal{A}$ is called _Hilbert-Schmidt_ if
${|\kern-1.07639pt|\kern-1.07639pt|\mathcal{A}|\kern-1.07639pt|\kern-1.07639pt|}_{2}=\big{\\{}\sum_{j\geq
1}\|\mathcal{A}\psi_{j}\|^{2}\big{\\}}^{1/2}$ is finite, where
$(\psi_{j})_{j\geq 1}$ is a ONB of $\mathcal{L}_{2}(\mathcal{Q})$. As
indicated by the notation,
${|\kern-1.07639pt|\kern-1.07639pt|\cdot|\kern-1.07639pt|\kern-1.07639pt|}_{2}$
does not depend on the particular choice of the basis, and is called the
_Hilbert-Schmidt norm_. We will use
$\mathcal{B}_{2}(\mathcal{L}_{2}(\mathcal{Q}))$ to denote the class of all
Hilbert-Schmidt operators on $\mathcal{L}_{2}(\mathcal{Q})$. An operator
$\mathcal{A}$ is called positive semi-definite if
$\langle\mathcal{A}f,f\rangle\geq 0$ for all
$f\in\mathcal{L}_{2}(\mathcal{Q})$. A compact, positive semi-definite operator
is called _trace class_ or _nuclear_ if
${|\kern-1.07639pt|\kern-1.07639pt|\mathcal{A}|\kern-1.07639pt|\kern-1.07639pt|}_{1}=\sum_{j\geq
1}\langle\mathcal{A}\psi_{j},\psi_{j}\rangle$ is finite for some ONB
$(\psi_{j})_{j\geq 1}$. Again, the sum is independent of the choice of ONB,
and
${|\kern-1.07639pt|\kern-1.07639pt|\mathcal{A}|\kern-1.07639pt|\kern-1.07639pt|}_{1}$
is called the _trace norm_ of $\mathcal{A}$. One particular type of operators
on $\mathcal{L}_{2}(\mathcal{Q})$, which is of interest to us is the _integral
operator_ defined as $\mathcal{A}f(\text{\bf u})=\int_{\mathcal{Q}}a(\text{\bf
u},\text{\bf v})\,f(\text{\bf v})\,\mathrm{d}\text{\bf v}$ for $\text{\bf
u}\in\mathcal{Q},f\in\mathcal{L}_{2}(\mathcal{Q})$, where
$a\in\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})$ is called the _kernel_ of
the operator $\mathcal{A}$. An integral operator $\mathcal{A}$ is positive
semi-definite if and only if the associated kernel $c$ is non-negative
definite. The operator $\mathcal{A}$ and the kernel $a$ are linked by an
obvious isometry, i.e.,
${|\kern-1.07639pt|\kern-1.07639pt|\mathcal{A}|\kern-1.07639pt|\kern-1.07639pt|}_{2}=\|c\|_{\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})}$.

Now, let $\mathcal{X}=(X(\text{\bf u}):\text{\bf u}\in\mathcal{Q})$ be a
random element in $\mathcal{L}_{2}(\mathcal{Q})$. For $d=1$, $\mathcal{X}$ is
usually referred to as a _random curve_ , whereas for $d>1$, it is referred to
as a _random field_ or _random surface_. We assume that $\mathcal{X}$ has
finite second moment, i.e., $\mathbb{E}(\|\mathcal{X}\|^{2})<\infty$, which
ensures the existence of its mean $m=\mathbb{E}(\mathcal{X})$ and covariance
$\mathcal{C}=\mathbb{E}\\{(\mathcal{X}-m)\otimes(\mathcal{X}-m)\\}$ of
$\mathcal{X}$ (both the expectations are understood in the Bochner sense). The
mean $m$ is an element of $\mathcal{L}_{2}(\mathcal{Q})$. The covariance
$\mathcal{C}$ is the integral operator associated with the _covariance kernel_
$c\in\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})$, where $c(\text{\bf
u},\text{\bf v})=\mathrm{Cov}(X(\text{\bf u}),X(\text{\bf v}))$. Moreover,
$\mathcal{C}$ is positive semi-definite and trace class. In this article, we
are interested in estimating $\mathcal{C}$ based on independent and
identically distributed (i.i.d.) observations
$\mathcal{X}_{1},\ldots,\mathcal{X}_{N}\sim\mathcal{X}$. Because of the
isomorphism linking the integral operator and the associated kernel, the
problem is equivalent to estimating the kernel $c$.

##  Appendix B Bias of the CovNet model: universal approximation and rate of
convergence

Here, we deal with the bias of the CovNet model. We start by proving that all
three CovNet structures can approximate any covariance operator up to
arbitrary precision, a.k.a. the _universal approximation_ property. These
results are instrumental for the consistency of our CovNet estimators. We
start with the shallow CovNet structure and give a detailed proof. The proofs
for the deep and the deepshared structures are similar, and we discuss those
only briefly.

###  B.1 Universal approximation

###### Proof of Theorem 1.

Recall that $c$ is the kernel of the covariance operator $\mathcal{C}$. So, by
the spectral decomposition of $\mathcal{C}$, we get

$c(\text{\bf u},\text{\bf v})=\sum_{i=1}^{\infty}\eta_{i}\,\psi_{i}(\text{\bf
u})\,\psi_{i}(\text{\bf v}),$ (B.1)

where the sum on the right converges in the $\mathcal{L}_{2}$ norm on
$\mathcal{Q}\times\mathcal{Q}$. Here, $\eta_{i}$’s are the eigenvalues of
$\mathcal{C}$ and $\psi_{i}$’s are the corresponding eigenfunctions. Since the
covariance operator $\mathcal{C}$ is trace-class, we get

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{C}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{1}=\sum_{i=1}^{\infty}\eta_{i}=\int_{\mathcal{Q}}c(\text{\bf
u},\text{\bf u})\,\,\mathrm{d}\text{\bf u}<\infty.$

Now, fix $\epsilon>0$ and w.l.o.g. let $\epsilon<1$. Since the sum in (B.1)
converges in the $\mathcal{L}_{2}$ norm, we can find an integer $I$ (depending
on $\epsilon$) such that

$\bigg{\|}c-\sum_{i=1}^{I}\eta_{i}\,\psi_{i}\otimes\,\psi_{i}\bigg{\|}_{\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})}<\frac{\epsilon}{2}.$
(B.2)

Also, since the functions $\psi_{1},\ldots,\psi_{I}$ are in
$\mathcal{L}_{2}(\mathcal{Q})$, we can find a positive constant $M$ (depending
on $\epsilon$) such that $\max_{i=1,\ldots,I}\big{\|}\psi_{i}\big{\|}\leq M$.
Since $\sigma$ is a sigmoidal function, using the density of single hidden
layer neural networks in the class of $\mathcal{L}_{2}$ functions (Györfi et
al., 2002, Theorem 16.2), for each $i=1,\ldots,I$, we can find
$R_{i}\in\mathbb{N}$, coefficients $a_{i,1},\ldots,a_{i,R_{i}}$, weights
$\mathbf{w}_{i,1},\ldots,\mathbf{w}_{i,R_{i}}\in\mathbb{R}^{d}$ and biases
$b_{i,1},\ldots,b_{i,R_{i}}$ such that

$\bigg{\|}\psi_{i}-\sum_{r=1}^{R_{i}}a_{i,r}\,\sigma(\mathbf{w}_{i,r}^{\top}\cdot+\,b_{i,r})\bigg{\|}<\frac{\epsilon}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathcal{C}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{1}(4M+2)}.$
(B.3)

Define $\widehat{\psi}_{i}(\text{\bf
u})=\sum_{r=1}^{R_{i}}a_{i,r}\,\sigma(\mathbf{w}_{i,r}^{\top}\text{\bf
u}+b_{i,r})$ for $i=1,\ldots,I$, and

$\displaystyle\widehat{c}_{R}(\text{\bf u},\text{\bf
v})=\sum_{i=1}^{I}\eta_{i}\,\widehat{\psi}_{i}(\text{\bf
u})\,\widehat{\psi}_{i}(\text{\bf v})$
$\displaystyle=\sum_{i=1}^{I}\eta_{i}\sum_{r=1}^{R_{i}}\sum_{s=1}^{R_{i}}a_{r,i}\,a_{s,i}\,\sigma(\mathbf{w}_{r,i}^{\top}\text{\bf
u}+b_{r,i})\,\sigma(\mathbf{w}_{s,i}^{\top}\text{\bf v}+b_{s,i})$
$\displaystyle=\sum_{r=1}^{R}\sum_{s=1}^{R}\lambda_{r,s}\,\sigma(\mathbf{w}_{r}^{\top}\text{\bf
u}+b_{r})\,\sigma(\mathbf{w}_{s}^{\top}\text{\bf v}+b_{s}),$ (B.4)

where $R=\sum_{i=1}^{I}R_{i}$. Here,
$\\{\mathbf{w}_{1},\ldots,\mathbf{w}_{R}\\}$ is the collection of all the
weights of the $I$ neural networks
$\widehat{\psi}_{1},\ldots,\widehat{\psi}_{I}$, and $\\{b_{1},\ldots,b_{R}\\}$
is the collection of all the biases. The associated matrix
$\Lambda=((\lambda_{r,s}))$ is block-diagonal with blocks
$\Lambda_{i}=\eta_{i}\,(a_{r,i}\,a_{s,i})_{1\leq r,s\leq R_{i}}$. Since
$\eta_{i}>0$, each of the $\Lambda_{i}$’s are positive semi-definite, which in
turn shows that $\Lambda$ is positive semi-definite. Define $c_{I}(\text{\bf
u},\text{\bf v})=\sum_{i=1}^{I}\eta_{i}\,\psi_{i}(\text{\bf
u})\,\psi_{i}(\text{\bf v})$. Using

$\|f\otimes f-g\otimes g\|_{\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})}\leq
2\|f\|\,\|f-g\|+\|f-g\|^{2},$ (B.5)

we get

$\displaystyle\|c_{I}-\widehat{c}_{R}\|_{\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})}$
$\displaystyle=\bigg{\|}\sum_{i=1}^{I}\eta_{i}\big{\\{}\psi_{i}\otimes\psi_{i}-\widehat{\psi}_{i}\otimes\widehat{\psi}_{i}\big{\\}}\bigg{\|}_{\mathcal{L}_{2}(\mathcal{Q}\times\mathcal{Q})}~{}~{}~{}~{}\text{(by
\eqref{eq:universal_approximation_proof_c_hat_L2})}$