# On the Identifiability and Interpretability of Gaussian Process Models

Jiawen Chen
Department of Biostatistics
University of North Carolina at Chapel Hill
<EMAIL_ADDRESS>
&Wancen Mu
Department of Biostatistics
University of North Carolina at Chapel Hill
<EMAIL_ADDRESS>
Yun Li
Department of Biostatistics
University of North Carolina at Chapel Hill
<EMAIL_ADDRESS>
&Didong Li
Department of Biostatistics
University of North Carolina at Chapel Hill
<EMAIL_ADDRESS>

###### Abstract

In this paper, we critically examine the prevalent practice of using additive
mixtures of Matérn kernels in single-output Gaussian process (GP) models and
explore the properties of multiplicative mixtures of Matérn kernels for multi-
output GP models. For the single-output case, we derive a series of
theoretical results showing that the smoothness of a mixture of Matérn kernels
is determined by the least smooth component and that a GP with such a kernel
is effectively equivalent to the least smooth kernel component. Furthermore,
we demonstrate that none of the mixing weights or parameters within individual
kernel components are identifiable. We then turn our attention to multi-output
GP models and analyze the identifiability of the covariance matrix $A$ in the
multiplicative kernel $K(x,y)=AK_{0}(x,y)$, where $K_{0}$ is a standard single
output kernel such as Matérn. We show that $A$ is identifiable up to a
multiplicative constant, suggesting that multiplicative mixtures are well
suited for multi-output tasks. Our findings are supported by extensive
simulations and real applications for both single- and multi-output settings.
This work provides insight into kernel selection and interpretation for GP
models, emphasizing the importance of choosing appropriate kernel structures
for different tasks.

##  1 Introduction

Gaussian processes (GPs) have emerged as a powerful and popular tool in
machine learning, spatial statistics, functional data analysis, etc., due to
their versatility, flexibility, and interpretability as a nonparametric method
(Rasmussen and Williams, 2006; Banerjee et al., 2014). GPs provide an
intuitive means of modeling uncertainty, allowing the generation of predictive
distributions for unseen data points without requiring explicit model
specification. Consequently, GPs have found applications in a variety of
contexts.

The key to harnessing the power of GPs lies in the choice of kernel functions,
also known as covariance functions. Rather than specifying a model, the GP
framework revolves around selecting an appropriate kernel, transforming the
problem into one of parameter estimation. Over time, researchers have
developed various kernels, each tailored to specific scenarios, such as
spatial data, time series data, and others (Genton, 2001). For instance, for
spatiotemporal, a variety of kernels have been developed to capture unique
characteristics and patterns (Gneiting, 2002; Stein, 2005).

The Radial Basis Function (RBF) kernel and Matérn kernels serve as prime
examples of the diversity within the kernel family (Cressie and Wikle, 2015).
The RBF kernel, renowned for its infinite differentiability, yields smooth
functions, while Matérn kernels control the degree of smoothness by a kernel
parameter, thereby accommodating both smooth and nonsmooth functions (Stein
(1999)). With each kernel bearing its own set of advantages and disadvantages,
kernel selection, as a nontrivial task, often demands extensive domain
knowledge.

Inspired by the potential of harnessing the strengths of multiple kernels
concurrently, researchers have ventured into methods that involve kernel
combinations, such as spectral mixtures (Wilson and Adams, 2013; Samo and
Roberts, 2015; Remes et al., 2017), addition and/or multiplication of kernels
(Duvenaud et al., 2011) such as mixture of RBF (Duvenaud et al., 2013),
mixture of RBF, periodic (Per), linear (Lin), and rational quadratic (RQ)
(Kronberger and Kommenda, 2013), mixture of RBF, RQ, Matérn and Per (Cheng et
al., 2019; Verma and Engelhardt, 2020), mixture of RBF, Matérn 1/2, 3/2, 5/2,
as well as more sophisticated methods like Neural Kernel Networks (NKN,Sun et
al. (2018) and Automatic Bayesian Covariance Discovery (ABCD, Lloyd et al.
(2014). For instance, within the realm of spectral mixtures, Remes et al.
(2017) proposed the Generalized Spectral Mixture (GSM) kernel, which is a
product of three components all parametrized by GPs, representing
frequencies,length-scales and mixture weights. In the category of summation
and/or multiplication, Verma and Engelhardt (2020) utilized the sum of the RBF
kernel and Matérn kernel with varying smoothness as the final kernel in their
t-distribution Gaussian process latent variable model (tGPLVM), an extension
to GPLVM (Lawrence, 2003; Lalchand et al., 2022), to characterize the latent
features in single-cell RNA sequencing data.

In addition, methodological advances have facilitated the efficient discovery
of optimal kernel mixtures. For example, NKN is reported to be more efficient
than the Automatic Statistician, a gradient-based method. Simpson et al.
(2021) utilized a transformer-based framework to generate mixture kernel
recommendations. These studies underscore the advantages of employing mixed
kernels in GP modeling, showcasing improved model fitting and more accurate
predictions. By combining multiple kernels, the strengths of individual
kernels can be leveraged to capture complex data patterns and relationships,
potentially exceeding the capabilities of single kernels. Furthermore, the
automatic selection of kernels optimizes model performance, reducing the
reliance on extensive prior knowledge.

Another essential benefit of using mixture kernels in GP modeling lies in the
improved interpretability. In contrast to complex, data-driven kernels,
mixture kernels allow the decomposition of intricate patterns into simpler,
distinct base kernels that are more readily interpretable. This technique,
often termed decomposition, enables a more comprehensive understanding of the
underlying data structure by simplifying complex patterns into their
constitutive components. An illustrative example of this approach is the work
of Duvenaud et al. (2013), who applied this decomposition technique for
structure discovery in time series data. Their proposed mixture kernel
comprised several base kernels, including RBF, periodic, linear, and rational
quadratic kernels, enabling the dissection of time-series data patterns into
components such as long-term trends, annual periodicity, and medium-term
deviations. In a similar vein, the ABCD approach also employed decomposition,
albeit with a different set of base kernels.

The identifiability of parameters within a single Matérn kernel has previously
been explored, with the microergodic parameter uniquely identified as the only
identifiable parameter, as outlined in Stein (1999). Tang et al. (2022)
examined the identifiablity of parameters within a single Matérn kernel with
nuggets. Despite the prevalent use and interpretation of mixture kernels,
their theoretical properties including identifiability and interpretability of
parameters within kernel components, to the best of our knowledge, remain
underexplored. A related critical question pertains to the common practice of
using kernels with varying degrees of smoothness for enhanced flexibility. In
this work, we turn our attention to the additive kernel in univariate GPs and
multiplicative (separable) kernel in multivariate GPs, specifically focusing
on the mixture of the widely-used Matérn kernel. We highlight the following
novel findings:

  * •

The smoothness of an additive mixture kernel is completely determined by the
least smooth component.

  * •

For additive mixture of Matérn kernels, the identifiability is confined only
to a single parameter, also known as the microergodic parameter that is
associated with the least smooth kernel component.

  * •

For multivariate GPs with multiplicative separable kernels, the multiplicative
matrix that controls the correlation structure among the response variables is
identifiable up to a multiplicative constant.

  * •

Our conclusions extend beyond the specific case of Matérn kernels,
demonstrating applicability to a wider range of mixture kernels.

Our study aims to deepen the understanding of mixture kernel identifiability
and interpretability, as well as to clarify the practical benefits of using
kernels with varying degrees of smoothness. Our theoretical assertions are
supported by both simulations and real-world applications. Details regarding
the proofs of our theories and numerical experiments can be found in the
Appendix.

##  2 Gaussian process, kernels and parameter inference

This section serves to define key terms and outline the parameter inference
algorithm integral to the forthcoming theory and simulation sections.

A GP is a random function where any finite set of its realizations follows a
multivariate normal distribution, characterized by a mean function $\mu$ and a
covariance function $K$ (Rasmussen and Williams, 2006).

######  Definition 1.

$f$ is said to follow a Gaussian process in domain $\Omega$ with a mean
function $\mu:\Omega\rightarrow\mathbb{R}$ and a covariance function
$K:\Omega\times\Omega\rightarrow\mathbb{R}$ if for any
$x_{1},\dots,x_{n}\in\Omega$,

$[f(x_{1}),\dots,f(x_{n})]^{\top}\sim N(v,\Sigma),\leavevmode\nobreak\
\text{where }v=[\mu(x_{1}),\dots,\mu(x_{n})]^{\top},\leavevmode\nobreak\
\Sigma_{ij}=K(x_{i},x_{j}).$

For the purpose of this study, we assume $\mu=0$, adhering to common practice
and for the sake of simplicity. If the mean is not zero, a data transformation
can be applied to achieve this. This approach is often adopted in machine
learning and statistical analysis to simplify calculations and comply with
certain algorithmic requirements (Rasmussen and Williams, 2006; Murphy, 2012).

As a consequence, the behavior of a GP is primarily determined by its kernel
function. In this study, we assume the domain $\Omega=\mathbb{R}^{p}$ and our
focus is on the selection of kernels that are widely recognized and commonly
used in machine learning.

######  Definition 2.

The RBF kernel, also known as the squared exponential kernel of Gaussian
kernel is given by:

$K(x,x^{\prime})=\sigma^{2}\exp(-\alpha\left\|x-x^{\prime}\right\|^{2}).$

The parameter $\sigma^{2}$ is called the spatial variance or partial sill that
controls the point-wise variance, while $\alpha$ is called the (length) scale,
range, or decay that controls the spatial dependency. A key characteristic of
RBF is its smoothness, defined as follows:

######  Definition 3.

A Gaussian process $f$ is said to be mean-square continuous (MSC) if
$\mathbb{E}(f(x+h)-f(x))^{2}\xrightarrow{h\to 0}0$ for any
$x\in\mathbb{R}^{p}$. $f$ is said to be mean-square differentiable (MSD) if
$\lim_{h\to 0}\frac{f(x+h)-f(x)}{h}$ exists in the mean-square topology, and
the limiting process is called the derivative process of $f$, denoted by
$f^{\prime}$. Similarly, the $d$-times mean-square differentiable ($d$-MSD)
GPs can be defined inductively.

In particular, the above-defined RBF kernel is infinitely differentiable.
However, oversmoothing can be problematic in prediction (Stein, 1999), leading
to the popularity of the following flexible family of kernels that allow for
varying smoothness:

######  Definition 4.

The Matérn kernel is given by

$K(x,x^{\prime})=\frac{\sigma^{2}2\left(\frac{\alpha}{2}\left\|x-x^{\prime}\right\|\right)^{\nu}}{\Gamma(\nu)}K_{\nu}(\alpha\left\|x-x^{\prime}\right\|)$

where $K_{\nu}$ is the modified Bessel function of the second kind. When
$\nu=1/2$, the Matérn kernel becomes the exponential kernel:

$K(x,x^{\prime})=\sigma^{2}\exp(-\alpha\left\|x-x^{\prime}\right\|)$

The additional parameter $\nu$ is called the smoothness, since the smoothness
of a Matérn GP is exactly $\lceil\nu\rceil-1$. Inferring $\nu$ is known to be
a particularly challenging problem, both from the theoretical and empirical
perspectives (Zhang, 2004). In practice, $\nu$ is usually set to
$1/2,3/2,5/2,\cdots$ due to the simplified analytic form of the modified
Bessel function.

The most straightforward estimator of the parameters $\sigma^{2}$ and $\alpha$
is the maximum likelihood estimator (MLE). In practice, optimizers such as
Adam (Kingma and Ba, 2015) or stochastic gradient descent (SGD) are used to
maximize the log-likelihood.

##  3 Univariate GP: additive kernels

In this section, we investigate the smoothness of mixture kernels and the
identifiability of parameters in the additive mixture of Matérn kernels in the
context of a univariate response variable. We note here that all the theorems
presented are framed in an asymptotic context.

Consider $K_{1},\cdots,K_{L}$ as $L$ kernels with $d_{l}$-MSD. Define
$K=\sum_{l=1}^{L}w_{l}K_{l}$ where $\sum_{l=1}^{L}w_{l}=1$ and $w_{l}\geq 0$.

######  Theorem 1.

$K$ is $d\coloneqq\min_{l}\\{d_{l}\\}$-times MSD, but not $d+1$-times MSD, so
the smoothness of $K$ is determined by the least smooth component.

Beyond smoothness, we discuss the identifiability of the parameters in the
mixture kernel, which relies on the notion of equivalence of measures.

######  Definition 5.

Two measures $\mathbb{P}_{1}$ and $\mathbb{P}_{2}$ are said to be equivalent,
i.e., $\mathbb{P}_{1}\equiv\mathbb{P}_{2}$, if they are absolutely continuous
with respect to each other, i.e.,
$\mathbb{P}_{1}(B)=0\Longleftrightarrow\mathbb{P}_{2}(B)=0$. Two GPs are
equivalent if the corresponding Gaussian random measures are equivalent.

As a consequence, two equivalent GPs cannot be distinguished by any finite
number of realizations (Stein, 1999). Specifically, given a family of GPs
parametrized by $\theta$, if
$\mathbb{P}_{\theta_{1}}\Longleftrightarrow\mathbb{P}_{\theta_{2}}$ with
$\theta_{1}\neq\theta_{2}$, then $\theta$ is not identifiable since we cannot
distinguish between $\theta_{1}$ and $\theta_{2}$. As a corollary, there does
not exist any consistent estimator for $\theta$.

Now we can study the identifiability of the mixture Matérn kernel
$K=\sum_{l=1}^{L}w_{l}K_{l}$, where $K_{l}$ is the Matérn kernel with
parameters $(\sigma^{2}_{l},\alpha_{l},\nu_{l})$, given known $\nu_{l}$’s in
ascending order $\nu_{1}<\nu_{2}\cdots<\nu_{L}$ with $\nu_{l}-\nu_{l-1}\geq 1$
for different smoothness.

######  Theorem 2.

Let $K=\sum_{l}w_{l}K_{l}$ and
$\widetilde{K}=\sum_{l}\widetilde{w}_{l}\widetilde{K}_{l}$ be two kernels
represented as linear combinations of Matérn kernels $K_{l}$’s and
$\widetilde{K}_{l}$ with parameters $(\sigma^{2}_{l},\alpha_{l},\nu_{l})$ and
$(\widetilde{\sigma}^{2}_{l},\widetilde{\alpha}_{l},\nu_{l})$.

  1. (i)

When $p=1,2,3$, Then $K\equiv\widetilde{K}$ if
$w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}=\widetilde{w}_{1}\widetilde{\sigma}_{1}^{2}\widetilde{\alpha}_{1}^{2\nu_{1}}.$

  2. (ii)

When $p\geq 5$, Then $K\equiv\widetilde{K}$ if

$w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}=\widetilde{w}_{1}\widetilde{\sigma}_{1}^{2}\widetilde{\alpha}_{1}^{2\nu_{1}},\leavevmode\nobreak\
w_{2}\sigma^{2}_{2}\alpha_{2}^{2\nu_{2}}-\nu_{1}w_{1}\sigma^{2}_{1}\alpha_{1}^{2(\nu_{1}+1)}=\widetilde{w}_{2}\widetilde{\sigma}^{2}_{2}\widetilde{\alpha}_{2}^{2\nu_{2}}-\nu_{1}\widetilde{w}_{1}\widetilde{\sigma}^{2}_{1}\widetilde{\alpha}_{1}^{2(\nu_{1}+1)}.$

Consequently, no single parameter is identifiable. The only identifiable
parameter when $p\leq 3$ is $w_{1}\sigma^{2}_{1}\alpha_{1}^{2\nu_{1}}$, known
as the microergodic parameter (Stein, 1999).

Interestingly, the case of $p=4$ is missing here. We hypothesize that it
aligns with the case of $p\geq 5$. However, providing a proof remains
challenging. It is also worth noting that the issue of identifiability
persists even for a single Matérn kernel, as highlighted in the existing
literature (Zhang, 2004; Anderes, 2010; Li et al., 2023).

######  Corollary 1.

The mixture kernel $K$ is equivalent to $K_{1}$, that is, the mixture kernel
is equivalent to its least smooth component. Furthermore, the mean squared
error (MSE) of the mixture kernel $K$ is asymptotically equal to the MSE of
$K_{1}$.

Theorem 2 implies that interpreting $w_{l}$ as the weight of each component is
often misleading due to its nonidentifiability. Furthermore, Corollary 1
suggests the use of a single Matérn kernel in lieu of the linear mixture
kernel. Both assertions are supported by the simulation study in Section 5 and
the real-world application in Section 6.

An alternative strategy is to assume the same smoothness across mixing
components, i.e., $\nu_{1}=\nu_{2},\cdots,\nu_{L}$. In this case, we prove
that no single parameter is identifiable. The complete theorem, proof, and
simulation are in the Appendix.

When considering real-world data, it is often the case that the observed
outcomes are noisy, which is frequently modeled as an i.i.d. Gaussian with
variance denoted by $\tau^{2}$, commonly referred as the "nugget". The
inclusion of such a noise term is essential to capture the inherent
variability in the data. Equivalently, the model can be formulated as
noiseless by adjusting the kernel to be
$K(x,x^{\prime})+\tau^{2}\mathbf{1}_{\\{x=x^{\prime}\\}}$. The following
corollary shows that the noise, or nugget, does not impact the identifiability
and interpretability of the mixture of Matérn kernels as shown in Theorem 2.

######  Corollary 2.

Let $\tau^{2}$ and $\widetilde{\tau}^{2}$ be the noise variance of $K$ and
$\widetilde{K}$. If $\tau^{2}\neq\widetilde{\tau}^{2}$, then
$K\not\equiv\widetilde{K}$. If $\tau^{2}=\widetilde{\tau}^{2}$, the previous
results in Theorem 2 hold.

In essence, while the presence of noise or "nugget" captures real-world data
variability, it does not obscure or alter the fundamental characteristics of
the mixture of Matérn kernels as established in our theorems.

##  4 Multivariate GP: multiplicative kernels

In this section, we shift our focus to multivariate GPs, also known as multi-
output or multi-task GPs, which are applicable to datasets with multiple
response variables.

######  Definition 6.

$f:\Omega\to\mathbb{R}^{m}$ is said to follow a $m$-variate GP in domain
$\Omega$ with zero mean and cross-covariance function
$K:\Omega\times\Omega\to\mathrm{PD}(m)$, where $\mathrm{PD}(m)$ is the space
of all $m$ by $m$ positive definite matrices, if for any $x_{1},\cdots,x_{n}$:

$[f(x_{1}),\cdots,f(x_{n})]^{\top}\sim MN(0,\Sigma),$

where $MN$ represents matrix Gaussian distribution and $\Sigma$ consists of
$m\times m$ blocks $\Sigma_{ij}=K(x_{i},x_{j})$.

The construction of valid cross-covariance kernels is a pivotal yet
challenging aspect of multivariate GPs (Gneiting, 2002; Apanasovich and
Genton, 2010; Genton and Kleiber, 2015). A popular kernel that is widely used
due to its simplicity and seemingly interpretability is the multiplicative
kernel, also known as a separable kernel. This kernel admits the form
$K(x,y)=AK_{0}(x,y)$ where $K_{0}$ is a standard kernel for univariate GP like
Matérn while $A\in\mathrm{PD}(m)$ is a positive definite matrix reflecting the
correlation between response variables. The following theorem investigates the
identifiability of the multiplicative cross-covariance kernel.

######  Theorem 3.

  1. (i)

When $K_{0}$ is Matérn with parameters $(\sigma^{2},\alpha,\nu)$ where $\nu$
is given, then $K\equiv\widetilde{K}$ if
$\sigma^{2}\alpha^{2\nu}A=\widetilde{\sigma}^{2}\widetilde{\alpha}^{2\nu}\widetilde{A}$.
That is, the identifiable parameter is $\sigma^{2}\alpha^{2\nu}A$. Hence $A$
is identifiable up to a multiplicative constant and the correlation structure
is identifiable.

  2. (ii)

Let $K_{0}$ be an arbitrary kernel with spectral density $\rho_{0}$
satisfying:

$\exists\gamma\in\mathcal{W}_{[-b,b]^{m}},\leavevmode\nobreak\
0<b<\infty,\leavevmode\nobreak\ c_{1},c_{2}>0,\leavevmode\nobreak\
s.t.,\leavevmode\nobreak\ c_{1}\gamma(\omega)^{2}\leq\rho_{0}(\omega)\leq
c_{2}\gamma(\omega)^{2},$

where $\mathcal{W}_{[-b,b]^{d}}$ is the space of Fourier transforms of
$L^{2}(\mathbb{R}^{p})$ functions with compact support $[-b,b]^{d}$ (see the
Appendix for more details). If $\theta$ is a mircroergodic parameter of
$K_{0}$, such that
$\int_{\mathbb{R}^{p}}\frac{1}{\gamma^{4}(\omega)}\left(\frac{\rho(\omega)}{\theta}-\frac{\widetilde{\rho}(\omega)}{\widetilde{\theta}}\right)^{2}\mathrm{d}\omega<0$
if $\theta=\widetilde{\theta}$, then $K\equiv\widetilde{K}$ if
$\sigma^{2}\alpha^{2\nu}A=\widetilde{\sigma}^{2}\widetilde{\alpha}^{2\nu}\widetilde{A}$.
That is, the identifiable parameter is $\theta A$, hence $A$ is identifiable
up to a multiplicative constant and the correlation structure is identifiable.

Note that (i) is a special case of (ii), as the Matérn kernel satisfies the
additional assumptions in (ii). This theorem positively supports the
utilization of the multiplicative kernel and the interpretation that $A_{ij}$
measures the correlation between the $i$-th response variable and the $j$-th
response variable.

##  5 Simulation

In this section, we will validate the proposed theories by conducting three
simulation studies. The first simulation demonstrates Theorem 1 using a
mixture of Matérn kernels, highlighting that the smoothness of the sample path
(e.g., a realization of the GP) is determined by the least smooth component.
The second simulation supports Theorem 2 using a mixture kernel of Matérn
kernels with varying smoothness. It shows that none of the single parameters
are identifiable, and only the proposed microergodic parameter can be
identified. Finally, the third simulation showcases Theorem 3, suggesting that
the multiplicative matrix is identifiable up to a multiplicative constant.

###  5.1 Simulation 1 - Univariate GP: Smoothness

In the first simulation, we show that the smoothness of the sample path is
determined by the smoothness of the least smooth kernel. We generate 200
samples as an equidistant sequence ranging from $0-10$. We consider the
following mixture kernel: $K=\sum_{l=1}^{3}w_{l}K_{l}$, where $K_{l}$ is the
Matérn kernel with smoothness parameter $\nu_{l}$ being $\frac{2l-1}{2}$ and
$\sum_{l=1}^{L}w_{l}=1,\leavevmode\nobreak\ w_{l}\geq 0$. We randomly sample
$Y$ from multivariate normal distribution with kernel function as the mixture
kernel, Matérn with $\nu=1/2,3/2,5/2$. The resulting sample paths are
summarized in Figure 1.

Figure 1: Smoothness of Matérn kernels. Sample path of (A) Mixture kernel (B)
Matérn $1/2$ (C) Matérn $3/2$ (D) Matérn $5/2$. (E) Numerical examination of
Theorem 1.

Our results clearly indicate that the Matérn kernel with $\nu=1/2$ exhibits
the least degree of smoothness (continuous but not differentiable), and the
smoothness increases with $\nu$. The smoothness of the mixture kernel is
predominantly influenced by the Matérn kernel with $\nu=1/2$. When comparing
the sample path of the mixture kernel to those of the Matérn kernels with
$\nu=3/2$ and $\nu=5/2$, it is evident that the mixture kernel demonstrates a
degree of smoothness similar to the Matérn kernels with $\nu=1/2$.

To better demonstrate the smoothness of the mixture kernel and its least
component, we further examine the continuity and MSD of the mixture kernel and
its mixing components empirically. For a fixed $x_{0}=0$, let $x_{i}=1/i$ for
$i=1,2,\dots$, and we generate $y_{i}^{l}$ from the GP, where $l=1,\dots,T$
denotes the index of replicates. This allows us to approximate $\lim_{x\to
0}\mathbb{E}(f(x)-f(0))^{2}$ by
$\beta_{i}=\frac{1}{T}\sum_{l=1}^{T}(y_{i}^{l}-y_{0}^{l})^{2}$. As per
Definition 3, the GP is continuous if and only if $\beta_{i}\to 0$. Similarly,
we can approximate $\lim_{x\to
0}\mathbb{E}\left(\frac{f(x)-f(0)}{x}\right)^{2}$ by
$\gamma_{i}=\frac{1}{T}\sum_{l=1}^{T}\left(\frac{y_{i}^{l}-y_{0}^{l}}{x_{i}}\right)^{2}$.
The GP is mean-square differentiable if and only if
$\lim_{i\to\infty}\gamma_{i}$ exists. Specifically, for the mixture of Matérn
$1/2$ and $3/2$ and Matérn $1/2$, while $\beta_{i}\to 0$, $\gamma_{i}$ does
not converge (Figure 1 first two columns), which indicates that both the
mixture kernel and Matérn $1/2$ are continuous but not differentiable.
However, for Matérn $3/2$ (Figure 1 third column), $\beta_{i}\to 0$ and
$\gamma_{i}$ converges, implying that Matérn $3/2$ is continuous and
differentiable.

These empirical observations strongly support the claims made in Theorem 1,
which suggests that the inclusion of smoother kernel components in the mixture
does not inherently enhance the overall smoothness of the mixture kernel.

###  5.2 Simulation 2 - Univariate GP: Components with different smoothness

In the second simulation, our aim is to assess parameter identifiability in a
GP with a mixture kernel consisting of Matérn kernels with distinct
smoothness. The mixture kernel is denoted as $K=\sum_{l=1}^{3}w_{l}K_{l}$,
where $K_{l}$ is a Matérn kernel with parameters
$(\sigma_{l}^{2},\alpha_{l},\nu_{l})$. For this simulation, we set
$\nu_{1}=1/2,\nu_{2}=3/2,\nu_{3}=5/2$. Theorem 2 indicates that the only
identifiable parameter is $w_{l}\sigma_{l}^{2}\alpha_{l}$, while all other
parameters remain unidentifiable. In this scenario, the identifiable
parameter, also known as the microergodic parameter, is associated with the
least smooth component, i.e., the Matérn $1/2$ kernel. We will assess
parameter inference by comparing estimated values to their true values across
different training sample sizes. For this experiment, we generate $n$ samples,
ranging from ${20,50,100,500}$. For each sample size, we replicate the
simulation 100 times.

Figure 2: Parameter estimation. (A) Only the microergodic parameter
$w_{1}\sigma_{1}^{2}\alpha_{1}$ is identifiable. (B-F) other parameters are
not identifiable. The red line indicates the true value.

The results provide persuasive support for Theorem 2 that, in the context of
additive mixture of Matérn kernels with distinct smoothness, only the MLE of
the microergodic parameter $w_{1}\sigma_{1}^{2}\alpha_{1}$ for the least
smooth component converges to the true value. The MLEs for all other
parameters do not converge to their respective true values, highlighting the
unique identifiability of the parameters in the least smooth component within
the additive mixture kernel framework. Such results are robust in terms of
optimizer choice and consistent as the sample size increases to $3000$
(details in the Appendix).

###  5.3 Simulation 3 - Multivariate GP

Next, we turn our attention to multivariate GPs with a separable kernel. We
define the separable kernel as $K=AK_{0}$. Here, we select Matérn $1/2$ as
$K_{0}$, characterized by the parameters $(\sigma^{2},\alpha)$. Our Theorem 3
suggests that only $A\alpha^{2\nu}\sigma^{2}$ is identifiable. To verify this,
we adopt a bivariate setup where $A$ is a $2\times 2$ positive definite
matrix. In this simulation, the sample size $n$ varies from
${50,100,200,400}$. For each sample size, we replicate the simulation 100
times.

Figure 3: Parameter estimation. (A-C) Only the microergodic parameter
$A\sigma^{2}\alpha^{2\nu}$ is identifiable (D-E) All other parameters are not
identifiable. The red dotted line indicates the true value.

The parameter estimates are summarized in Figure 3. Evidently, only the MLE of
the microergodic parameter $A\sigma^{2}\alpha^{2\nu}$ converges to its true
value. These findings reinforce the identifiability and interpretability of
the separable kernel. This understanding is critical for the interpretation
and application of separable kernels in real-world situations.

##  6 Application

Motivated by the findings in our theoretical and simulation studies, we
proceed to compare the performance of the mixture kernel with varying
smoothness and the least smooth component within that mixture only, in real-
world data applications. In this section, we will revisit several applications
in previous studies (Rasmussen and Williams (2006), Wilson et al. (2014), Sun
et al. (2018)), illustrating that the mixture kernel achieves prediction
performance comparable to the least smooth kernel component within the
mixture.

###  6.1 Application 1 - Image prediction

For our first real-data application, we delve into the domain of image
analysis. GP has been widely used to predict missing image (Wilson et al.
(2014), Sun et al. (2018)). This task of image prediction is an exemplary
context for assessing kernel performance in discerning local correlation
patterns, providing a visually intuitive framework for comparative analysis.

We employ a handwritten zero image from the MNIST dataset (LeCun et al.
(1998)) for this analysis; we extract a $20\times 20$ section from the center
of the original $100\times 100$ image to serve as our test image. This results
in a training dataset of $9600$ pixels and a testing dataset with $400$ pixels
(Figure 4A). The input for this analysis comprises the 2-D pixel locations,
with the corresponding pixel intensities serving as the response. To
reconstruct the missing area, we train GP regression models with both the
mixture kernel with varying smoothness and a single Matérn $1/2$ kernel on the
training dataset, and subsequently make predictions on the testing dataset.
Specifically, we employ a mixture of three Matérn kernels of smoothness
$1/2,\leavevmode\nobreak\ 3/2,\leavevmode\nobreak\ 5/2$. Our analysis seeks to
demonstrate the similarities in performance between the mixture kernel and the
least smooth component within that mixture, i.e., Matérn $1/2$, in terms of
prediction in the missing image area.

Figure 4: Texture exploration in diamond plate. (A) Training image. (B)
Prediction of the mixture kernel. The area in the pink box represents the
testing image area. (C) Predication of the Matérn $1/2$ kernel.

The predictions generated by both the mixture kernel and the Matérn 1/2 kernel
demonstrate a remarkable similarity, as shown in Figure 4. This observation
suggests that the mixture kernel and a single Matérn $1/2$ kernel have
comparable performance. The close alignment between their predictions implies
that using the mixture kernel does not yield significant advantages over the
least smooth kernel component in this specific application, further
reinforcing our theoretical findings (Theorem 2) regarding the dominance of
the least smooth kernel component in the mixture kernel.

###  6.2 Application 2 - Manua Loa CO2

In this section, we extend our comparative analysis to the widely recognized
Moana Loa CO2 dataset from the Global Monitoring Laboratory’s Repository as
our benchmark for regression analysis (Tans and Keeling (2023)). Rasmussen and
Williams (2006) employed this dataset to demonstrate that constructing a
complex covariance function by combining several different types of simple
covariance functions could yield an excellent fit to complex data.
Consequently, we consider it to be a suitable choice for testing our theorem.

Specifically, we used the decimal date from 1960 to 2020 as our input and the
monthly average CO2 as our response. In this context, we compare the
performance of Matérn mixture $1/2+3/2,\leavevmode\nobreak\
1/2+3/2+5/2,\leavevmode\nobreak\ 3/2+5/2$ and three single Matérn kernels
($1/2,\leavevmode\nobreak\ 3/2,\leavevmode\nobreak\ 5/2$). The dataset
consists of 720 records. We undertook an extensive assessment, varying the
training sample size from $5\%$ to $95\%$ and keeping the test sample size as
the remaining data. For each training sample size, we conducted 10 simulations
using different random splits of the dataset as replications. This strategy
allows for a robust and thorough evaluation of the prediction accuracy of both
the mixture kernel and the single kernel component in the context of a
regression problem with varying training sample sizes. We use the MSE to
quantify the prediction accuracy (Figure 5).

Figure 5: Manua Loa CO2 prediction MSE in test dataset with 20%, 50% and 75%
training dataset. Boxplot are colored by different Matérn kernels.

The Matérn mixture kernel $1/2+3/2+5/2$ aligned closely with $1/2+3/2$ across
all experiments. Furthermore, these three Matérn mixture kernels
($1/2+3/2+5/2$, $1/2+3/2$, $1/2$) converge when the training sample size ratio
reaches $50\%$. This empirical finding supports the theoretical equivalence of
these three GPs as posited in Theorem 2. Similarly, the Matérn mixture
$3/2+5/2$ kernel and the standalone Matérn $3/2$ kernel converge from a $65\%$
sample size, lending empirical weight to their theoretical equivalence and
their equivalence in MSE as stated in Corollary 1.

Matérn $3/2$ kernel has exceptional performance across all experiments.
However, its influential performance is diluted when coupled with the less
smooth component, Matérn $1/2$. On the contrary, when mixed with the smoother
Matérn $5/2$ component, the worse performance of the latter is overshadowed by
Matérn $3/2$. These behaviors lend empirical weight to our Theorem 2,
suggesting that in the asymptotic sense, the mixture kernel is dominated by
its least smooth component.

A noteworthy trend is the inequivalance of the Matérn mixture kernel and its
least component with a limited training sample (training sample size $<50\%$).
We consider such scenario as the "finite" sample scenario in contrast to our
focus on the infinite sample scenario, i.e. asymptotic theory. Theoretically,
while the Matérn mixture kernels $1/2+3/2+5/2$, $1/2+3/2$, and the standalone
Matérn $1/2$ kernel are deemed equivalent, and the Matérn mixture kernel
$3/2+5/2$ is akin to the Matérn $3/2$ kernel, their practical agreement
points, i.e. $50\%$ and $65\%$, differed with limited training samples.
Although securing a conclusive theoretical support for the finite sample
regime is challenging, we managed to provide some insights to interpret our
findings. To do so, we have to delve further into the proof of Stein (1993)’s
Theorem 1, which guides us to Theorem 3.1 in Stein (1990). The proof suggests
that the relative difference between the MSEs rests on the tail of the series.
This difference is influenced by both the sample size (denoted as $N$ in Stein
(1990)) and by $b_{jk}$ and $\mu_{j}$, as defined on the same page. For a
fixed sample size, smaller values of $(b_{jk}+\mu_{j}\mu_{k})^{2}$ result in a
smaller relative difference between MSEs. By the definition of $b_{jk}$ and
$\mu_{j}$, the more "different" the two kernels are, in other words, the more
"significant" the additional smoother components become — the greater the
difference between the MSEs will be.

##  7 Discussion and future work

In this study, we delved into the identifiability and interpretability of
parameters in GPs with various mixture kernels in the asymptotic scenario and
fixed domain, including the additive kernel and separable kernels. We
formulated a series of theorems clarifying the identifiable parameters in
these kernel structures, and further corroborated our theorems through
multiple simulations and real-data applications. Our simulation results
convincingly demonstrated that in GPs with mixture kernels, the only
identifiable parameter, known as the microergodic parameter, is associated
with the least smooth kernel component. This discovery has profound
implications for parameter interpretation in the context of GPs with mixture
kernels. Empirical evidence from image data analysis and Manua Loa CO2
regression studies further reinforce these discoveries. Despite the inclusion
of kernels with varying smoothness in the mixture kernel, the performance of
the mixture kernel closely paralleled that of the least smooth kernel
component within it when the sample size is large enough, regardless of its
performance. This observation suggests that the inclusion of kernels with
different smoothness does not necessarily improve the prediction accuracy. In
fact, due to various real world factors including limited training samples,
optimization and more, determining a clear winner in performance between
Matérn mixture kernel and single Matérn kernel proves to be a challenging
task. Lastly, in the case of multivariate GPs with separable kernels, our
theoretical and simulation results show that the correlation structure is
identifiable up to a multiplicative constant. This result underlines that the
interpretability of separable kernels mainly resides in the relative
correlation structure, rather than individual parameters.

Although our study has provided substantial insight on the identifiability and
interpretability of parameters in GPs with various mixture kernel types, it
also opens several exciting avenues for future research. First, our analysis
primarily focuses on Matérn kernels, and it would be intriguing to extend this
framework to other families of kernels, such as the periodic kernel. Such an
extension would provide a more comprehensive understanding of parameter
identifiability and predictive performance across a broader spectrum of kernel
types. Second, our work has so far considered the cases where $p\leq 3$ or
$p\geq 5$. However, extending this to $p=4$ presents a substantial challenge
due to the lack of mathematical tools to determine whether two Gaussian random
measures are equivalent or not. Third, our observations in the Mauna Loa CO2
example underscore the need for further exploration in the finite sample
scenario. Potential avenues for future research include investigating the
convergence rates of both the Matérn mixture kernel and the single Matérn
kernel. Lastly, while we have identified the microergodic parameters that are
theoretically identifiable in mixture kernels, an important direction for
future work involves finding consistent estimators for these kernel
parameters. This would entail developing novel estimation techniques or
adapting existing ones to reliably estimate the parameters in practice,
thereby enhancing the practical utility of our theoretical findings. These
endeavors will not only extend the theoretical foundations of GPs with mixture
kernels, but will also broaden their applicability across various real-world
scenarios.

## Acknowledgments and Disclosure of Funding

We express our sincere gratitude to anonymous reviewers and the Area Chair for
their insightful comments and constructive feedback, which significantly
enhanced the quality and clarity of this work. DL was supported by the
NIH/NCATS award UL1 TR002489, NIH/NHLBI award R01 HL149683 and NIH/NIEHS award
P30 ES010126. YL was partially supported by NIH awards R01AG079291 and
U01HG011720. All the codes for replicating the analysis and parameter
estimates are available at: https://github.com/JiawenChenn/GP_mixture_kernel.

## References

  * Anderes (2010) Anderes, E. (2010).  On the consistent separation of scale and variance for Gaussian random fields.  The Annals of Statistics.
  * Apanasovich and Genton (2010) Apanasovich, T. V. and M. G. Genton (2010).  Cross-covariance functions for multivariate random fields based on latent dimensions.  Biometrika 97(1), 15–30.
  * Bachoc et al. (2022) Bachoc, F., E. Porcu, M. Bevilacqua, R. Furrer, and T. Faouzi (2022).  Asymptotically equivalent prediction in multivariate geostatistics.  Bernoulli 28(4), 2518–2545.
  * Banerjee et al. (2014) Banerjee, S., B. P. Carlin, and A. E. Gelfand (2014).  Hierarchical modeling and analysis for spatial data.  CRC press.
  * Cheng et al. (2019) Cheng, L., S. Ramchandran, T. Vatanen, N. Lietzén, R. Lahesmaa, A. Vehtari, and H. Lähdesmäki (2019).  An additive Gaussian process regression model for interpretable non-parametric analysis of longitudinal data.  Nature Communications 10(1), 1798.
  * Cressie and Wikle (2015) Cressie, N. and C. K. Wikle (2015).  Statistics for spatio-temporal data.  John Wiley & Sons.
  * Duvenaud et al. (2013) Duvenaud, D., J. Lloyd, R. Grosse, J. Tenenbaum, and G. Zoubin (2013, 17–19 Jun).  Structure discovery in nonparametric regression through compositional kernel search.  In S. Dasgupta and D. McAllester (Eds.), Proceedings of the 30th International Conference on Machine Learning, Volume 28 of Proceedings of Machine Learning Research, Atlanta, Georgia, USA, pp. 1166–1174. PMLR.
  * Duvenaud et al. (2011) Duvenaud, D. K., H. Nickisch, and C. Rasmussen (2011).  Additive Gaussian processes.  In J. Shawe-Taylor, R. Zemel, P. Bartlett, F. Pereira, and K. Weinberger (Eds.), Advances in Neural Information Processing Systems, Volume 24. Curran Associates, Inc.
  * Genton (2001) Genton, M. G. (2001).  Classes of kernels for machine learning: a statistics perspective.  Journal of machine learning research 2(Dec), 299–312.
  * Genton and Kleiber (2015) Genton, M. G. and W. Kleiber (2015).  Cross-Covariance Functions for Multivariate Geostatistics.  Statistical Science 30(2), 147 – 163.
  * Gneiting (2002) Gneiting, T. (2002).  Nonseparable, stationary covariance functions for space–time data.  Journal of the American Statistical Association 97(458), 590–600.
  * Kingma and Ba (2015) Kingma, D. P. and J. Ba (2015).  Adam: A method for stochastic optimization.  In Y. Bengio and Y. LeCun (Eds.), 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings.
  * Kronberger and Kommenda (2013) Kronberger, G. and M. Kommenda (2013).  Evolution of covariance functions for Gaussian process regression using genetic programming.
  * Lalchand et al. (2022) Lalchand, V., A. Ravuri, and N. D. Lawrence (2022).  Generalised gplvm with stochastic variational inference.  In International Conference on Artificial Intelligence and Statistics, pp. 7841–7864. PMLR.
  * Lawrence (2003) Lawrence, N. (2003).  Gaussian process latent variable models for visualisation of high dimensional data.  In S. Thrun, L. Saul, and B. Schölkopf (Eds.), Advances in Neural Information Processing Systems, Volume 16. MIT Press.
  * LeCun et al. (1998) LeCun, Y., C. Cortes, and C. J. Burges (1998).  The mnist database of handwritten digits.  http://yann.lecun.com/exdb/mnist/.  Accessed: 2023-08-01.
  * Li et al. (2023) Li, D., W. Tang, and S. Banerjee (2023).  Inference for gaussian processes with matérn covariogram on compact riemannian manifolds.  Journal of Machine Learning Research 24(101), 1–26.
  * Lloyd et al. (2014) Lloyd, J. R., D. Duvenaud, R. Grosse, J. B. Tenenbaum, and Z. Ghahramani (2014).  Automatic construction and natural-language description of nonparametric regression models.
  * Murphy (2012) Murphy, K. P. (2012).  Machine learning: a probabilistic perspective.  MIT press.
  * Pedregosa et al. (2011) Pedregosa, F., G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay (2011).  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research 12, 2825–2830.
  * Rasmussen and Williams (2006) Rasmussen, C. E. and C. K. I. Williams (2006).  Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press.
  * Remes et al. (2017) Remes, S., M. Heinonen, and S. Kaski (2017).  Non-stationary spectral kernels.  In I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Eds.), Advances in Neural Information Processing Systems, Volume 30. Curran Associates, Inc.
  * Samo and Roberts (2015) Samo, Y.-L. K. and S. Roberts (2015).  Generalized spectral kernels.
  * Simpson et al. (2021) Simpson, F., I. Davies, V. Lalchand, A. Vullo, N. Durrande, and C. E. Rasmussen (2021).  Kernel identification through transformers.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P. Liang, and J. W. Vaughan (Eds.), Advances in Neural Information Processing Systems, Volume 34, pp. 10483–10495. Curran Associates, Inc.
  * Stein (1990) Stein, M. (1990).  Uniform asymptotic optimality of linear predictions of a random field using an incorrect second-order structure.  The Annals of Statistics 18(2), 850–872.
  * Stein (1993) Stein, M. L. (1993).  A simple condition for asymptotic optimality of linear predictions of random fields.  Statistics & Probability Letters 17(5), 399–404.
  * Stein (1999) Stein, M. L. (1999).  Interpolation of spatial data: some theory for kriging.  Springer Science & Business Media.
  * Stein (2005) Stein, M. L. (2005).  Space–time covariance functions.  Journal of the American Statistical Association 100(469), 310–321.
  * Sun et al. (2018) Sun, S., G. Zhang, C. Wang, W. Zeng, J. Li, and R. Grosse (2018, 10–15 Jul).  Differentiable compositional kernel learning for Gaussian processes.  In J. Dy and A. Krause (Eds.), Proceedings of the 35th International Conference on Machine Learning, Volume 80 of Proceedings of Machine Learning Research, pp. 4828–4837. PMLR.
  * Tang et al. (2022) Tang, W., L. Zhang, and S. Banerjee (2022).  On identifiability and consistency of the nugget in gaussian spatial process models.  Journal of the Royal Statistical Society Series B: Statistical Methodology 83(5), 1044–1070.
  * Tans and Keeling (2023) Tans, P. and R. Keeling (2023).  Trends in atmospheric carbon dioxide.  https://gml.noaa.gov/ccgg/trends/data.html.  Accessed: 2023-08-01.
  * Verma and Engelhardt (2020) Verma, A. and B. E. Engelhardt (2020).  A robust nonlinear low-dimensional manifold for single cell RNA-seq data.  BMC Bioinformatics 21(1), 324.
  * Wilson and Adams (2013) Wilson, A. and R. Adams (2013, 17–19 Jun).  Gaussian process kernels for pattern discovery and extrapolation.  In S. Dasgupta and D. McAllester (Eds.), Proceedings of the 30th International Conference on Machine Learning, Volume 28 of Proceedings of Machine Learning Research, Atlanta, Georgia, USA, pp. 1067–1075. PMLR.
  * Wilson et al. (2014) Wilson, A. G., E. Gilboa, A. Nehorai, and J. P. Cunningham (2014).  Fast kernel learning for multidimensional pattern extrapolation.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K. Weinberger (Eds.), Advances in Neural Information Processing Systems, Volume 27. Curran Associates, Inc.
  * Yadrenko and Balakrishnan (1983) Yadrenko, M. I. and A. V. Balakrishnan (1983).  Spectral theory of random fields.
  * Zhang (2004) Zhang, H. (2004).  Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics.  Journal of the American Statistical Association 99(465), 250–261.

##  Appendix A Proof of Theorem 1

We first state a Lemma for determining the smoothness for a stationary GP.

######  Lemma 1 (Stein (1999)).

Let $\rho$ be the spectral density of a GP with kernel $K$, then $K$ is
d-times mean squared differentiable if and only if

$\int\|\omega\|^{2d}\rho(\omega)\mathrm{d}\omega<\infty.$

###### Proof of Theorem 1.

Let $\rho_{l}$ be the spectral density of $K_{l}$, then the spectral density
of $K$ is given by $\rho=\sum_{l=1}^{L}\rho_{l}$ by the linearity of Fourier
transform. Then by the smoothness of $K_{l}$ and Lemma 1, we have
$\int\|\omega\|^{2d_{l}}\rho_{l}(\omega)\mathrm{d}\omega<\infty.$ As a result,

$\displaystyle\int\|\omega\|^{2d}\rho(\omega)\mathrm{d}\omega$
$\displaystyle=\sum_{l=1}^{L}\int\|\omega\|^{2d}\rho_{l}(\omega)\mathrm{d}\omega=\sum_{l=1}^{L}\int_{\|\omega\|\leq
1}\|\omega\|^{2d}\rho_{l}(\omega)\mathrm{d}\omega+\sum_{l=1}^{L}\int_{\|\omega\|\geq
1}\|\omega\|^{2d}\rho_{l}(\omega)\mathrm{d}\omega$
$\displaystyle\leq\sum_{l=1}^{L}\int_{\|\omega\|\leq
1}\rho_{l}(\omega)\mathrm{d}\omega+\sum_{l=1}^{L}\int\|\omega\|^{2d_{l}}\rho_{l}(\omega)\mathrm{d}\omega<\infty,$

that is, $K$ is $d$-times MSD. Similarly, $K$ is not $d+1$-times MSD since

$\displaystyle\int\|\omega\|^{2(d+1)}\rho(\omega)\mathrm{d}\omega$
$\displaystyle=\sum_{l=1}^{L}\int\|\omega\|^{2(d+1)}\rho_{l}(\omega)\mathrm{d}\omega\geq\int\|\omega\|^{2(d+1)}\rho_{1}(\omega)\mathrm{d}\omega=\infty.$

∎

##  Appendix B Proof of Theorem 2

We first state a Lemma, also known as the integral test to determine whether
two GPs are equivalent or not.

######  Lemma 2 (Yadrenko and Balakrishnan (1983), Zhang (2004)).

Let $K_{1}$ and $K_{2}$ be the kernels of two GPs $P_{1}$ and $P_{2}$ with
spectral densities $\rho_{1}$ and $\rho_{2}$. Then $P_{1}\equiv P_{2}$ if

  1. 1.

There exists $r>0$ such that $\|\omega\|^{r}\rho_{1}(\omega)$ is bounded away
from both zero and infinity as $\omega\to\infty$.

  2. 2.

There exists $\delta>0$ such that
$\int_{\|\omega\|>\delta}\left(\frac{\rho_{2}(\omega)-\rho_{1}(\omega)}{\rho_{1}(\omega)}\right)^{2}\mathrm{d}\omega<\infty$.

###### Proof of Theorem 2.

Recall that the spectral densities of $K_{l}$ are

$\rho_{l}(\omega)=\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}},\leavevmode\nobreak\
\rho(\omega)=\sum_{l=1}^{L}w_{l}\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}},\leavevmode\nobreak\
\widetilde{\rho}(\omega)=\sum_{l=1}^{L}\widetilde{\omega}_{l}\frac{\widetilde{\sigma}_{l}^{2}\widetilde{\alpha}_{l}^{2\nu_{l}}}{(\widetilde{\alpha}_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}.$

To apply the integral test, first, let $r=2(\nu_{1}+\frac{p}{2})$, then

$\rho(\omega)\|\omega\|^{r}=\sum_{l=1}^{L}w_{l}\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}\|\omega\|^{r}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}\xrightarrow[]{\omega\to\infty}w_{1}\sigma_{1}^{2}\alpha_{l}^{2\nu_{1}},$

that is, bounded away from both $0$ and $\infty$ as $\omega\to\infty$. Then we
check the relative difference between the spectral densities $\rho$ and
$\widetilde{\rho}$. Observe that when,

$\displaystyle\rho(\omega)$
$\displaystyle=\sum_{l=1}^{L}w_{l}\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}=\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}\Pi_{j\neq
l}(\alpha_{j}^{2}+\|\omega\|^{2})^{\nu_{j}+\frac{p}{2}}}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}\eqqcolon\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}p_{l}(\omega)}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}},$

where $p_{l}\coloneqq\Pi_{j\neq
l}(\alpha_{j}^{2}+\|\omega\|^{2})^{\nu_{j}+\frac{p}{2}}$. Similarly,

$\widetilde{\rho}(\omega)=\frac{\sum_{l=1}^{L}\widetilde{w}_{l}\widetilde{\sigma}_{l}^{2}\widetilde{\alpha}_{l}^{2\nu}\widetilde{p}_{l}(\omega)}{\Pi_{l=1}^{L}(\widetilde{\alpha}_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}},$

where $\widetilde{p}_{l}\coloneqq\Pi_{j\neq
l}(\widetilde{\alpha}_{j}^{2}+\|\omega\|^{2})^{\nu_{j}+\frac{p}{2}}$. Observe
that when $\|\omega\|\to\infty$,

$\displaystyle\frac{\Pi_{l=1}^{L}(\widetilde{\alpha}_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}=\frac{\Pi_{l=1}^{L}\left(1+\frac{\widetilde{\alpha}_{l}^{2}}{\|\omega\|^{2}}\right)^{\nu_{l}+\frac{p}{2}}}{\Pi_{l=1}^{L}\left(1+\frac{\alpha_{l}^{2}}{\|\omega\|^{2}}\right)^{\nu_{l}+\frac{p}{2}}}=\frac{1+O(\|\omega\|^{-2})}{1+O(\|\omega\|^{-2})}=1+O(\|\omega\|^{-2}).$

Furthermore, let $\gamma_{l}=\sum_{j\neq l}(\nu_{j}+\frac{p}{2})$, then
$p_{l}(\omega)=O(\|\omega\|^{2\gamma_{l}})$ where
$\gamma_{1}>\gamma_{2}\cdots>\gamma_{L}$ and $\gamma_{l}-\gamma_{l-1}\geq 1$.
As a result, the leading term that dominates
$\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}p_{l}(\omega)$ is
$w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}p_{1}(\omega)$, as
$\|\omega\|\to\infty$. Moreover,
$\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}p_{l}(\omega)=\|\omega\|^{2\gamma_{1}}\left(w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}+O(\|\omega\|^{-2})\right)$
since $2\gamma_{l}-2\gamma_{1}\geq 2,\leavevmode\nobreak\ \forall l>1$.

Now we can analyze the relative difference:

$\displaystyle\left|\frac{\rho(\omega)-\widetilde{\rho}(\omega)}{\rho(\omega)}\right|$
$\displaystyle=\left|\frac{\widetilde{\rho}(\omega)}{\rho(\omega)}-1\right|=\left|\frac{\sum_{l=1}^{L}\widetilde{w}_{l}\widetilde{\sigma}_{l}^{2}\widetilde{\alpha}_{l}^{2\nu_{l}}\widetilde{p}_{l}(\omega)}{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}p_{l}(\omega)}\frac{\Pi_{l=1}^{L}(\widetilde{\alpha}_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu_{l}+\frac{p}{2}}}-1\right|$
$\displaystyle=\left|\frac{\sum_{l=1}^{L}\widetilde{w}_{l}\widetilde{\sigma}_{l}^{2}\widetilde{\alpha}_{l}^{2\nu_{l}}\widetilde{p}_{l}(\omega)}{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu_{l}}p_{l}(\omega)}(1+O(\|\omega\|^{-2}))-1\right|$
$\displaystyle=\left|\frac{\|\omega\|^{2\gamma_{1}}(\widetilde{w}_{1}\widetilde{\sigma}_{1}^{2}\widetilde{\alpha}_{1}^{2\nu_{1}}+O(\|\omega\|^{-2})}{\|\omega\|^{2\gamma_{1}}({w}_{1}{\sigma}_{1}^{2}{\alpha}_{1}^{2\nu_{1}}+O(\|\omega\|^{-2})}(1+O(\|\omega\|^{-2}))-1\right|$
$\displaystyle=\left|\frac{\widetilde{w}_{1}\widetilde{\sigma}_{1}^{2}\widetilde{\alpha}_{1}^{2\nu_{1}}}{w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}}(1+O(\|\omega\|^{-2}))(1+O(\|\omega\|^{-2}))-1\right|$
$\displaystyle=\left|\frac{\widetilde{w}_{1}\widetilde{\sigma}_{1}^{2}\widetilde{\alpha}_{1}^{2\nu_{1}}}{w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}}-1+O(\|\omega\|^{-2})\right|.$

As a result, if
$w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}=\widetilde{w}_{1}\widetilde{\sigma}_{1}^{2}\widetilde{\alpha}_{1}^{2\nu_{1}}$,

$\int_{\omega\in\mathbb{R}^{p}:\|\omega\|>1}\left|\frac{\rho(\omega)-\widetilde{\rho}(\omega)}{\rho(\omega)}\right|^{2}\mathrm{d}\omega\approx\int_{\omega\in\mathbb{R}^{p}:\|\omega\|>1}\frac{1}{\|\omega\|^{4}}\mathrm{d}\omega<\infty,$

when $p=1,2,3$, so $K\equiv\widetilde{K}$. That is, none of the parameters are
identifiable, while the parameter that might be identifiable is
$w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}$, also known as the microergodic
parameter.

Now we turn to the case for $p\geq 5$. By Anderes (2010), it suffices to
anaylize the principal irregular terms for $K$. For the spectral density of
each individual kernel component, the principal irregular term is

$G_{\nu_{l}}(\omega)=\frac{-\pi}{2^{2\nu_{l}}\sin(\nu_{l}\pi)\Gamma(\nu_{l})\Gamma(\nu_{l}+1)}\omega^{2\nu_{l}}\eqqcolon
C_{l}\omega^{2\nu_{l}}.$

Furthermore,

$K_{l}(x+h,x^{\prime})=\sigma_{l}^{2}G_{\nu_{l}}(|\alpha_{l}h|)-\nu_{l}\sigma_{l}^{2}G_{\nu_{l}+1}(|\alpha_{l}h|)+\iota_{l}(|\alpha_{l}h|),$

where $\iota_{l}(t)=p(|h|)+o(G_{\nu_{l}+1}(|h|))$ as $|h|\to 0$ and $p$ is a
polynomial with even degree. By linearity, we have

$\displaystyle K(x+h,x^{\prime})$
$\displaystyle=\sum_{l=1}^{L}w_{l}K_{l}(x+h,x^{\prime})$
$\displaystyle=\sum_{l=1}^{L}w_{l}\left(\sigma_{l}^{2}G_{\nu_{l}}(|\alpha_{l}h|)-\nu_{l}\sigma^{2}_{l}G_{\nu_{l}+1}(|\alpha_{l}h|)+\iota_{l}(t)\right)$
$\displaystyle=w_{1}\sigma_{1}^{2}G_{\nu_{1}}(|\alpha_{1}h|)-w_{1}\nu_{1}\sigma_{1}^{2}G_{\nu_{1}+1}(|\alpha_{1}h|)+w_{2}\sigma_{2}^{2}G_{\nu_{2}}(|\alpha_{2}h|)+\iota(t),$

where $\iota(t)=p(|h|)+o(G_{\nu_{2}+1}(|h|))$ as $|h|\to 0$. They by Theorem 4
of Anderes (2010), there exists consistent estimators for
$w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}$ and
$w_{2}\sigma^{2}_{2}\alpha_{2}^{2\nu_{2}}-\nu_{1}w_{1}\sigma_{1}^{2}\alpha_{1}^{2(\nu_{1}+1)}$
when $0<2(2\nu_{2}-2\nu_{1})<p$, that is $0<4<p$. This completes the proof. ∎

##  Appendix C Proof of Corollary 1

We first state a Lemma for comparing MSE of two best linear predictor with two
spectral density.

######  Lemma 3 (Stein (1993)).

Suppose $Z$ is a zero mean stationary process in $\mathbb{R}^{d}$ and
$x_{1},x_{2},\dots$ is a dense sequence of points in a bounded subset of
$\mathbb{R}^{d}$. $x_{0}$ is a point in the bounded set but not the sequence.
Let $\hat{Z}(x_{0},n,\rho)$ be the best linear predictor of $Z(x_{0})$ based
on $Z(x_{1}),\dots,Z(x_{n})$ assuming $\rho$ is the spectral density for $Z$.
$e(x_{0},n,\rho)\coloneqq Z(x_{0})-\hat{Z}(x_{0},n,\rho)$. If there exists
c>0, such that

$\lim_{|\omega|\to\infty}\frac{\rho_{1}(\omega)}{\rho_{0}(\omega)}=c$

and $\rho_{0}$ satisfies condition 1 in Lemma 2, then

$\lim_{n\to\infty}\frac{\mathbb{E}_{\rho_{0}}e(x_{0},n,\rho_{1})^{2}}{\mathbb{E}_{\rho_{0}}e(x_{0},n,\rho_{0})^{2}}=1$

###### Proof of Corollary 1.

For $K_{1}=\mathrm{Mat}(w_{1}\sigma_{1}^{2},\alpha_{1},\nu_{1})$, by Theorem
2, $K_{1}\equiv K$ since the microergodic parameters match.

Regarding MSE, let $\rho$ be the spectral density of the mixture kernel $K$,
and $\rho_{1}$ be the spectral density of the first mixing component $K_{1}$.
We assume $\rho_{1}$ as the true spectral density. Then by the proof of
Theorem 2,

$\lim_{|\omega|\to\infty}\frac{\rho(\omega)}{\rho_{1}(\omega)}=\frac{{w}_{1}{\sigma}_{1}^{2}{\alpha}_{1}^{2\nu_{1}}}{w_{1}\sigma_{1}^{2}\alpha_{1}^{2\nu_{1}}}=1,$

then the MSE of GP with $K_{1}$ is asymptocially equal to the MSE of GP with
$K$. ∎

##  Appendix D Proof of Corollary 2

We first state a Lemma regarding the equivalence of Gaussian measure with
nuggets.

######  Lemma 4 (Tang et al. (2022)).

Let ${S}$ be a closed set, $G_{S}(m,K)$ denotes the Gaussian measure of the
random field on ${S}$ with mean function $m$ and covariance function $K$. Let
$w(s)$ be a mean square continuous process on $S$ under $G_{S}(m_{1},K_{1})$
and $\chi$ be a dense sequence of points in ${S}$. Then

  1. 1.

if $\tau^{2}_{1}\neq\tau_{2}^{2}$, then
$G_{\chi}(m_{1},K_{1},\tau_{1}^{2})\perp G_{\chi}(m_{2},K_{2},\tau_{2}^{2})$.

  2. 2.

if $\tau^{2}_{1}=\tau_{2}^{2}$, then $G_{\chi}(m_{1},K_{1},\tau_{1}^{2})\equiv
G_{\chi}(m_{2},K_{2},\tau_{2}^{2})$ if and only if $G_{S}(m_{1},K_{1})\equiv
G_{S}(m_{2},K_{2})$.

###### Proof of Corollary 2.

We denote the adjusted kernel with noise of $K$ and $\widetilde{K}$ as
$K_{\tau}$ and $\widetilde{K}_{\widetilde{\tau}}$.

From Lemma 4(1), where we set $m_{1}=m_{2}=0$, if
$\tau^{2}\neq\widetilde{\tau}^{2}$, then
$K_{\tau}\not\equiv\widetilde{K}_{\widetilde{\tau}}$. If
$\tau^{2}=\widetilde{\tau}^{2}$, then
$K_{\tau}\equiv\widetilde{K}_{\widetilde{\tau}}$ if and only if
$K\equiv\widetilde{K}$. Then the previous results in Theorem 2 holds. ∎

##  Appendix E Additional theorem - mixture of Matérn kernels with the same
smoothness

Here we introduce a special scenario for Theorem 4 where all component in the
mixture kernel have same smoothness. The simulation study for the following
theorems are included in Section G: simulation 4.

######  Theorem 4.

For $K=\sum_{l=1}^{L}K_{l}$ where
$K_{l}=\mathrm{Mat}(\sigma^{2}_{l},\alpha_{l},\nu)$,then

  1. (i)

When $p=1,2,3$, the only identifiable parameter is
$\sum_{l=1}^{L}w_{l}\sigma^{2}_{l}\alpha_{l}^{2\nu}$.

  2. (ii)

When $p>4$, the only identifiable parameters are
$\sum_{l=1}^{L}w_{l}\sigma^{2}_{l}\alpha_{l}^{2\nu}$ and
$\sum_{l=1}^{L}w_{l}\sigma^{2}_{l}\alpha_{l}^{2\nu+2}$.

As a consequence, none of the parameters are identifiable.

###### Proof.

The proof is similar to the proof of Theorem 2. Recall that the spectral
densities of $K_{l}$ are

$\rho_{l}(\omega)=\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}},\leavevmode\nobreak\
\rho(\omega)=\sum_{l=1}^{L}w_{l}\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}},\leavevmode\nobreak\
\widetilde{\rho}(\omega)=\sum_{l=1}^{L}\widetilde{\omega}_{l}\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}.$

To apply the integral test, first, let $r=2(\nu+\frac{p}{2})$, then

$\rho(\omega)\|\omega\|^{r}=\sum_{l=1}^{L}w_{l}\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu}\|\omega\|^{r}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}\xrightarrow[]{\omega\to\infty}\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu},$

that is, bounded away from both $0$ and $\infty$ as $\omega\to\infty$. Then we
check the relative difference between the spectral densities $\rho$ and
$\widetilde{\rho}$. Observe that

$\displaystyle\rho(\omega)$
$\displaystyle=\sum_{l=1}^{L}w_{l}\frac{\sigma_{l}^{2}\alpha_{l}^{2\nu}}{(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}=\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}\Pi_{j\neq
l}(\alpha_{j}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}\eqqcolon\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}},$

where $p_{l}\coloneqq\Pi_{j\neq
l}(\alpha_{j}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}$. Observe that

$\displaystyle|p_{l}(\omega)|$
$\displaystyle\geq\|\omega\|^{2(L-1)(\nu+\frac{p}{2})}$ (1)
$\displaystyle\frac{p_{l}(\omega)}{\|\omega\|^{2(L-1)(\nu+\frac{p}{2})}}$
$\displaystyle=\Pi_{j\neq
l}\left(1+\frac{\alpha_{j}^{2}}{\|\omega\|^{2}}\right)^{\nu+\frac{p}{2}}=1+O(\|\omega\|^{-2}).$
(2)

Now we can analyze the relative difference:

$\displaystyle\left|\frac{\rho(\omega)-\widetilde{\rho}(\omega)}{\rho(\omega)}\right|$
$\displaystyle=\left|\frac{\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}-\frac{\sum_{l=1}^{L}\widetilde{w}_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}}{\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)}{\Pi_{l=1}^{L}(\alpha_{l}^{2}+\|\omega\|^{2})^{\nu+\frac{p}{2}}}}\right|$
$\displaystyle=\left|\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)-\sum_{l=1}^{L}\widetilde{w}_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)}{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)}\right|$
$\displaystyle\leq\frac{1}{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}}\left|\frac{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)-\sum_{l=1}^{L}\widetilde{w}_{l}\sigma_{l}^{2}\alpha^{2\nu}p_{l}(\omega)}{\|\omega\|^{2(L-1)(\nu+\frac{p}{2})}}\right|$
$\displaystyle=\frac{1}{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}}\left|\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}(1+O(\|\omega\|^{2})-\sum_{l=1}^{L}\widetilde{w}_{l}\sigma_{l}^{2}\alpha^{2\nu}(1+O(\|\omega\|^{2})\right|$
$\displaystyle=\frac{1}{\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}}\left|\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}-\sum_{l=1}^{L}\widetilde{w}_{l}\sigma_{l}^{2}\alpha^{2\nu}+O(\|\omega\|^{-2})\right|.$

As a result, if
$\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha^{2\nu}=\sum_{l=1}^{L}\widetilde{w}_{l}\sigma_{l}^{2}\alpha^{2\nu}$,

$\int_{\omega\in\mathbb{R}^{p}:\|\omega\|>1}\left|\frac{\rho(\omega)-\widetilde{\rho}(\omega)}{\rho(\omega)}\right|^{2}\mathrm{d}\omega\approx\int_{\omega\in\mathbb{R}^{p}:\|\omega\|>1}\frac{1}{\|\omega^{4}\|}\mathrm{d}\omega<\infty,$

when $p=1,2,3$, so $K\equiv\widetilde{K}$. That is, none of the parameters are
identifiable, while the parameter that might be identifiable is
$\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu}$, also known as the
microergodic parameter.

Now we turn to the case for $p\geq 5$. The proof is similar to the proof of
Theorem 2: it suffices to analyze the principal irregular terms for $K$. For
the spectral density of each individual kernel component, the principal
irregular term is

$G_{\nu}(\omega)=\frac{-\pi}{2^{2\nu}\sin(\nu\pi)\Gamma(\nu)\Gamma(\nu+1)}\omega^{2\nu}\eqqcolon
C_{l}\omega^{2\nu}.$

Furthermore,

$K_{l}(x+h,x^{\prime})=\sigma_{l}^{2}G_{\nu}(|\alpha_{l}h|)-\nu\sigma^{2}G_{\nu+1}(|\alpha_{l}h|)+\iota(|\alpha_{l}h|),$

where $\iota(t)=p(|h|)+o(G_{\nu+1}(|h|))$ as $|h|\to 0$ and $p$ is a
polynomial with even degree. By linearity, we have

$\displaystyle K(x+h,x^{\prime})$
$\displaystyle=\sum_{l=1}^{L}w_{l}K_{l}(x+h,x^{\prime})$
$\displaystyle=\sum_{l=1}^{L}w_{l}\left(\sigma_{l}^{2}G_{\nu}(|\alpha_{l}h|)-\nu\sigma^{2}_{l}G_{\nu+1}(|\alpha_{l}h|)+\iota(t)\right)$
$\displaystyle=\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}G_{\nu}(|\alpha_{l}h|)-\sum_{l=1}^{L}\nu
w_{l}\sigma_{l}^{2}G_{\nu+1}(|\alpha_{l}h|)+\iota(t),$

where $\iota(t)=p(|h|)+o(G_{\nu+1}(|h|))$ as $|h|\to 0$. They by Theorem 4 of
Anderes (2010), there exists consistent estimators for
$\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu}$ and
$\sum_{l=1}^{L}w_{l}\sigma^{2}_{l}\alpha_{l}^{2(\nu+1)}$ when
$0<2(2(\nu+1)-2\nu)<p$, that is $0<4<p$. This completes the proof. ∎

##  Appendix F Proof of Theorem 3

(i) is a direct corollary of Theorem 3 of Bachoc et al. (2022), where
$\sigma_{ii}\sigma_{jj}\rho_{ij}=A_{ij}\sigma^{2}$ and
$\alpha_{ij}=\frac{1}{\alpha}$ for any $i,j=1,\cdots,m$.

Now we prove (ii). Let $\theta$ be the microergodic parameter of $K_{0}$, that
is,
$K_{0}\equiv\widetilde{K}_{0}\Longleftrightarrow\theta=\widetilde{\theta}$,
where $\widetilde{K}_{0}$ is characterized by $\widetilde{\theta}$. Let
$\rho_{0}$ and $\widetilde{\rho_{0}}$ be the spectral densities of $K$ and
$\widetilde{K}_{0}$, then the matrix spectral densities of $K$ and
$\widetilde{K}$ are $P\coloneqq A\rho_{0}$ and
$\widetilde{P}\coloneqq\widetilde{A}\widetilde{\rho}_{0}$. By the assumption,
that is, there exists constants $b,\leavevmode\nobreak\
c_{1},\leavevmode\nobreak\ c_{2}>0$ and $\gamma\in\mathcal{W}_{[-b,b]^{m}}$
such that $c_{1}\gamma^{2}(\omega)\leq\rho_{0}(\omega)\leq
c_{2}\gamma^{2}(\omega)$, we claim that there exists
$q_{1},\leavevmode\nobreak\ q_{2}>0$ such that

$q_{1}\gamma^{2}(\omega)\mathrm{I}_{m}\leq P(\omega)\leq
q_{2}\gamma^{2}(\omega)\mathrm{I}_{m},\leavevmode\nobreak\
\forall\omega\in\mathbb{R}^{p}.$ (3)

From the assumption,

$\displaystyle P(\omega)=A\rho_{0}\leq c_{2}\gamma^{2}(\omega)A\leq
c_{2}\operatorname{eig}_{1}(A)\gamma^{2}(\omega),$

where $\operatorname{eig}_{1}(A)$ is the largest eigenvalue of $A$. Similarly,
$c_{1}\operatorname{eig}_{m}(A)\gamma^{2}(\omega)\mathrm{I}_{m}\leq P(\omega)$
where $\operatorname{eig}_{m}(A)>0$ is the smallest eigenvalue of $A$, which
is positive by the positive definiteness of $A$. Then claim (3) holds where
$q_{1}\coloneqq c_{1}\operatorname{eig}_{m}(A)$ and $q_{1}\coloneqq
c_{2}\operatorname{eig}_{1}(A)$. Now we can prove (ii). If $\theta
A=\widetilde{\theta}\widetilde{A}$, then for any $i,j=1,\cdots,m$,

$\displaystyle\int_{\mathbb{R}^{p}}\frac{1}{\gamma^{4}(\omega)}\left(P_{ij}(\omega)-\widetilde{P}_{ij}(\omega)\right)^{2}\mathrm{d}\omega$
$\displaystyle=\int_{\mathbb{R}^{p}}\frac{1}{\gamma^{4}(\omega)}\left(A_{ij}\rho_{0}(\omega)-\widetilde{A}_{ij}\widetilde{\rho}_{0}(\omega)\right)^{2}\mathrm{d}\omega$
$\displaystyle=\int_{\mathbb{R}^{p}}\frac{1}{\gamma^{4}(\omega)}\left(\theta
A_{ij}\frac{\rho_{0}(\omega)}{\theta}-\widetilde{\theta}\widetilde{A}_{ij}\frac{\widetilde{\rho}_{0}(\omega)}{\widetilde{\theta}}\right)^{2}\mathrm{d}\omega$
$\displaystyle=\int_{\mathbb{R}^{p}}\frac{1}{\gamma^{4}(\omega)}\left(\theta
A_{ij}\frac{\rho_{0}(\omega)}{\theta}-{\theta}{A}_{ij}\frac{\widetilde{\rho}_{0}(\omega)}{\widetilde{\theta}}\right)^{2}\mathrm{d}\omega$
$\displaystyle=(\theta
A_{ij})^{2}\int_{\mathbb{R}^{p}}\frac{1}{\gamma^{4}(\omega)}\left(\frac{\rho_{0}(\omega)}{\theta}-\frac{\widetilde{\rho}_{0}(\omega)}{\widetilde{\theta}}\right)^{2}\mathrm{d}\omega$
$\displaystyle<\infty,$

where the last inequality holds from the assumption that $\theta$ is the
microergodic parameter of $K_{0}$. As a result, by Theorem 2 of Bachoc et al.
(2022), $\theta A$ is the microergodic parameter of $K$.

##  Appendix G Experiment details

###  G.1 Simulation 1

In the first simulation, we examined a mixture kernel defined as follows:
$K=\sum_{l=1}^{3}w_{l}K_{l}$, wherein $K_{l}$ represents the Matérn kernel
with parameters $(\sigma^{2}_{l},\alpha_{l},\nu_{l})$. For this simulation, we
assigned the values $\nu_{1}=1/2,\nu_{2}=3/2,\nu_{3}=5/2$ for the smoothness
parameter.

The weights for each kernel component, represented as $(w_{1},w_{2},w_{3})$,
were chosen as $(0.03,0.33,0.63)$. This set-up presents the case where despite
$w_{1}$ being the smallest, which is associated with the kernel that exhibits
the least smoothness, our result further substantiates our claim about the
dominance of the kernel with the lowest smoothness over the influence of
weights. In terms of the scale parameters, for
$(\sigma^{2}_{1},\sigma^{2}_{2},\sigma^{2}_{3})$, we picked $(3,3,3)$. We
selected $(\alpha_{1},\alpha_{2},\alpha_{3})$ to be $(1,1,1)$. This uniform
selection across the components aids in isolating the effect of the smoothness
and weights in our analysis.

###  G.2 Simulation 2

In the second simulation, we consider the following mixture kernel:
$K=\sum_{l=1}^{3}w_{l}K_{l}$, where $K_{l}$ is the Matérn kernel with
parameter $(\sigma^{2}_{l},\alpha_{l},\nu_{l})$. For this simulation, the
values for the smoothness parameter, $\nu$, were set as
$\nu_{1}=1/2,\nu_{2}=3/2,\nu_{3}=5/2$. Regarding the choice of the weight
parameters $(w_{1},w_{2},w_{3})$, we selected $(0.1,0.3,0.6)$ to reflect the
difference in contribution of each kernel to the overall mixture. Although the
weight of the Matérn $1/2$ kernel ($w_{1}$) is small, it still remains
dominant due to its lesser smoothness. This setup allowed us to test our
hypothesis that smoothness plays a more significant role than weight in
determining the identifiability of parameters.
$(\sigma^{2}_{1},\sigma^{2}_{2},\sigma^{2}_{3})$ are set to $(16,4,1)$, and
$(\alpha_{1},\alpha_{2},\alpha_{3})$ are set to $(4,2,1)$. This variation
further facilitated the examination of our proposition that the parameters
associated with the least smooth kernel converge to their true values, while
others do not.

In this scenario, we generate $X\in\mathbb{R}^{2}$ from an equal-spaced
sequence ranging from $-10$ to $10$, with a random variation sampled from
$\mathrm{unif}(-\frac{1}{5n},\frac{1}{5n})$. The sample size, $n$, takes on
values from the set ${20,50,100,500}$. For each sample size, we replicate the
simulation 100 times. Subsequently, $Y$ is simulated from $N(0,K+\epsilon)$.
We consider $\epsilon$ as a fixed value ($\epsilon=0.1$) and include it for
numerical robustness. This $\epsilon$ is also added in the training process.
For parameter initialization, we initialize $w=(0.2,0.3,0.5)$,
$\sigma^{2}=(5.0067,10,15)$ and $\alpha=(0.7615,0.4702,0.3280)$. These initial
values were close to the true parameters, yet sufficiently distinct to
illustrate the efficacy of the learning process. Here we use the SGD optimizer
with $0.005$ learning rate and $1000$ epochs. All parameter estimates are
summarised in the Figure S1.

Figure S1: All parameter estimation in simulation 2.

We further examined the simulation in a larger sample size and also used a
different optimizer L-BFGS (learning rate 0.5 for $n=50,500$, 1.0 for
$n=100$). The results are consistent with the previous results.

Figure S2: Parameter estimation in simulation 2 with larger sample size.
Figure S3: Parameter estimation in simulation 2 with L-BFGS optimizer.

###  G.3 Simulation 3

In the third simulation, we consider the following mixture kernel: $K=AK_{0}$,
where $K_{0}$ is the Matérn kernel with parameter $(\alpha,\sigma^{2},\nu)$.
For this simulation, we set $\nu=1/2$. The ground truth for the parameters is
assigned as follows: the matrix $A$ is set as $[[5,1][1,5]]$, which serves as
a symmetric, positive-definite structure to facilitate the properties of the
covariance matrix; $\sigma^{2}$ is set to $10$, and $\alpha$ is set to $1$.

In this simulation, we generate $X\in\mathbb{R}^{2}$ from an equal-spaced
sequence ranging from $-10$ to $10$, with a random variation sampled from
$\mathrm{unif}(-\frac{1}{5n},\frac{1}{5n})$. The sample size, $n$, takes on
values from the set ${50,100,200,400}$. For each sample size, we replicate the
simulation 100 times. Subsequently, $Y$ is simulated from $N(0,K+\epsilon)$.
We consider $\epsilon$ as a fixed value ($\epsilon=0.5$). This $\epsilon$ is
also added in the training process. In terms of parameter initialization, we
elect to start with $A$ as an identity matrix, $\sigma^{2}=1$ and $\alpha=10$.
The initialization of $A$ as an identity matrix ensures a simple, non-
informative starting point, while the initial $\sigma^{2}$ and $\alpha$ are
chosen to be significantly distinct from the ground truth to assess the
robustness of the learning process. Here we use the SGD optimizer with $0.001$
learning rate and $2000$ epochs. The parameter estimates of $A$ are summarised
in the Figure S4.

Figure S4: Parameter estimation in simulation 3.

###  G.4 Additional simulation (Simulation 4) for mixture of Matérn kernels
with the same smoothness

In this additional simulation 4, we aim to evaluate the parameter
identifiability in a GP with a mixture kernel consisting of Matérn kernels
with the same smoothness. We denote the mixture as
$K=\sum_{l=1}^{L}w_{l}K_{l}$, where $K_{l}$ is a Matérn kernel with parameters
$(\sigma_{l}^{2},\alpha_{l},\nu)$. Our theorem suggests that only the
microergodic parameter $\sum_{l=1}^{L}w_{l}\sigma^{2}_{l}\alpha_{l}$ might be
identifiable, while other parameters are not identifiable. Consequently, we
will evaluate the parameter estimation by comparing it with the true value for
different training sample sizes.

Here we use mixture of three Matérn $1/2$ kernels. We generate
$X\in\mathbb{R}^{2}$ from an equal-spaced sequence with values ranging from
$-10$ to $10$, with a random variation sampled from
$\mathrm{unif}(-\frac{1}{10n},\frac{1}{10n})$. The sample size $n$ takes
values from ${20,50,100,500}$. For each sample size, we replicate the
simulation 100 times. Subsequently, we simulate $Y$ from a normal distribution
with a mean of 0 and covariance matrix $K+\epsilon$. We treat $\epsilon$ as a
fixed value ($\epsilon=0.1$) and include it for numerical robustness. This
$\epsilon$ is also added in the training process. The ground truth parameters
are given as $w=(0.2,0.3,0.5)$, $\sigma^{2}=(16,4,1)$ and $\alpha=(2,1,4)$.
Here the weight and sigma are set to true value,
$(\alpha_{1},\alpha_{2},\alpha_{3})$ are set to $(4,2,1)$. The motivation
behind our specific parameter initialization strategy is to intentionally
introduce some degree of initial discrepancy. This approach allows us to
critically observe the convergence behavior of the parameters. Here we use the
Adam optimizer with $0.01$ learning rate and $1000$ epochs.

Figure S5: Parameter estimation in simulation 4. (A) Only the microergodic
parameter $\sum_{l=1}^{L}w_{l}\sigma^{2}_{l}\alpha_{l}$ might be identifiable
(B-D) All other parameters $(\sigma_{l}^{2},\alpha_{l},w_{l})_{l=1}^{L}$ are
not identifiable. Figure S6: All parameter estimation in simulation 4.

Our results offer compelling evidence that, in the case of a mixture kernel
consisting of three Matérn 1/2 kernels, only the MLE of the microergodic
parameter $\sum_{l=1}^{L}w_{l}\sigma_{l}^{2}\alpha_{l}^{2\nu}$ converges to
the true value (Figure S5A). The MLEs of all other parameters do not
demonstrate convergence to their corresponding true values (Figure S5B-D).
This result underscores the importance of understanding the identifiability of
parameters in such mixture kernels and highlights the need for careful
consideration when interpreting the estimated parameters in Gaussian process
models. When treating mixture kernel consist of same smoothness, we could not
interprete every single model. All parameter estimates are summarised in
Figure S6.

###  G.5 Application 1

In application 1, we focus on image analysis. We employ a hand-written zero
from the MNIST dataset (LeCun et al. (1998)). In this application, we use
mixture kernel with Matérn kernels of smoothness $1/2$, $3/2$ and $5/2$ and
compare its performance with Matérn $1/2$ kernel. For both kernels, we used
the SGD optimizer with learning rate $1e^{-05}$. The training epochs are
$40000$. During both the training and prediction stages, we introduced a term,
$\epsilon$, into the covariance computation to guarantee its positive-
definiteness. For both kernels, we set $\epsilon$ to $0.01$. For mixture
kernel, the parameters are initialized as $(w_{1},w_{2},w_{3})=(1/3,1/3,1/3)$,
$(\sigma^{2}_{1},\sigma^{2}_{2},\sigma^{2}_{3})=(1,1,1)$ and
$(\alpha_{1},\alpha_{2},\alpha_{3})=(1,1,1)$.

Here the estimated parameters for mixture kernel are
$(w_{1},w_{2},w_{3})=(9.9954e^{-01},2.3194e^{-04},2.3194e^{-04})$,
$(\sigma^{2}_{1},\sigma^{2}_{2},\sigma^{2}_{3})=(1.2148e^{+02},4.9960e^{+01},4.9256e^{+01})$,
$(\alpha_{1},\alpha_{2},\alpha_{3})=(1.3941e^{-01},1.1203e^{-02},1.1020e^{-02})$.
The estimated parameters for Matérn 1/2 kernel are
$\alpha=3.0639e^{-02}$,$\sigma^{2}=2.7927e^{+02}$.

###  G.6 Application 2

In application 2, we focus on the Moana Loa CO2 dataset (Tans and Keeling
(2023)). In this application, we compare the performance of Matérn mixture
$1/2+3/2,\leavevmode\nobreak\ 1/2+3/2+5/2,\leavevmode\nobreak\ 3/2+5/2$ and
three single Matérn kernels ($1/2,\leavevmode\nobreak\
3/2,\leavevmode\nobreak\ 5/2$). We use sklearn package (Pedregosa et al.
(2011)) for this analysis. The the parameters are initialized as
$(\sigma^{2}_{1},\sigma^{2}_{2},\sigma^{2}_{3})=(10,500,500)$ and
$(\alpha_{1},\alpha_{2},\alpha_{3})=(4,0.1,0.1)$ for both single kernel and
mixture kernel. Other parameters remain default in sklearn.


# Partial sums of Horadam sequences: sum-free representations via generating
functions

Helmut Prodinger  Helmut Prodinger, Mathematics Department, Stellenbosch
University, 7602 Stellenbosch, South Africa<EMAIL_ADDRESS>
###### Abstract.

Horadam sequences and their partial sums are computed via generating
functions. The results are as simple as possible.

##  1\. Introduction

Horadam sequences [2, 1] $W_{n}=W_{n}(a,b;p,q)$ are defined via

$W_{n}=pW_{n-1}+qW_{n-2},\quad n\geq 2,\quad W_{0}=a,\ W_{1}=b.$

These numbers are of course a generalization of Fibonacci numbers, Lucas
numbers and many others. The characteristic equation

$X^{2}-pX-q=0$

is essential, and the two roots are

$\lambda=\frac{p+\sqrt{p^{2}+4q}}{2},\quad\mu=\frac{p-\sqrt{p^{2}+4q}}{2}.$

We define

$F_{n}=\frac{\lambda^{n}-\mu^{n}}{\lambda-\mu}\quad\text{and}\quad
L_{n}=\lambda^{n}+\mu^{n},$

as these sequences resemble Fibonacci resp. Lucas numbers, and each solution
of the recursion may be expressed as a linear combination of these two
sequences. We have $F_{0}=0$, $F_{1}=1$, $L_{0}=2$, $L_{1}=p$. For instance

$W_{n}=\bigl{(}b-\frac{ap}{2}\bigr{)}F_{n}+\frac{a}{2}L_{n}.$

The paper [1] concentrates on finding expressions for

$\sum_{n\leq k\leq n+m}W_{k}=\sum_{0\leq k\leq n+m}W_{k}-\sum_{0\leq k\leq
n-1}W_{k}.$

In the rest of this short paper, we will find simple expressions for

$S_{n}:=\sum_{0\leq k\leq n}W_{k}$

using generating functions. The results do not contain summations, and can be
expressed with the sequences $F_{n}$ and $L_{n}$.

##  2\. Generating functions

Standard computations produce

$W(z)=\sum_{k\geq 0}W_{k}z^{k}=\frac{a+z(b-pa)}{1-pz-qz^{2}};$

furthermore

$F(z)=\sum_{k\geq 0}F_{k}z^{k}=\frac{z}{1-pz-qz^{2}}\quad\text{and}\quad
L(z)=\sum_{k\geq 0}L_{k}z^{k}=\frac{2-pz}{1-pz-qz^{2}}.$

By general principles,

$\displaystyle S(z)$ $\displaystyle=\sum_{k\geq
0}S_{k}z^{k}=\frac{1}{1-z}\frac{a+z(b-pa)}{1-pz-qz^{2}}$
$\displaystyle=\frac{a-pa+b}{1-p-q}\frac{1}{1-z}+\frac{-b-qa+qz(pa-
a-b)}{1-p-q}\frac{1}{1-pz-qz^{2}}$
$\displaystyle=\frac{a-pa+b}{1-p-q}\frac{1}{1-z}$ $\displaystyle-\frac{2qa-
pqa+2qb+pb}{2(1-p-q)}\frac{z}{1-pz-
qz^{2}}-\frac{qa+b}{2(1-p-q)}\frac{2-pz}{1-pz-qz^{2}}.$

Reading off the coefficient of $z^{n}$ on both sides leads to

$\displaystyle S_{n}=\frac{a-pa+b}{1-p-q}-\frac{2qa-
pqa+2qb+pb}{2(1-p-q)}F_{n}-\frac{qa+b}{2(1-p-q)}L_{n}.$

This answers the finite sum problem addressed in [1] completely, since the
answer is

$S_{n+m}-S_{n-1}=-\frac{2qa-
pqa+2qb+pb}{2(1-p-q)}(F_{n+m}-F_{n-1})-\frac{qa+b}{2(1-p-q)}(L_{n+m}-L_{n-1}).$

## References

  * [1] C. Cooper.  Finite sums of consecutive terms of a second order linear recurrence relation.  Integers, 21, #A114, 2021.
  * [2] A. F. Horadam.  Basic properties of a certain generalized sequence of numbers.  Fibonacci Quart., 3, 161–176, 1965.


* [165] ANITA, P. W. Gorham et al., “Constraints on the ultrahigh-energy cosmic neutrino flux from the fourth flight of ANITA,” Phys. Rev. D 99 (2019) no. 12, 122001, arXiv:1902.04005.
  * [166] ANITA, P. W. Gorham et al., “Characteristics of Four Upward-pointing Cosmic-ray-like Events Observed with ANITA,” Phys. Rev. Lett. 117 (2016) no. 7, 071101, arXiv:1603.05218.
  * [167] Pierre Auger, A. Aab et al., “Probing the origin of ultra-high-energy cosmic rays with neutrinos in the EeV energy range using the Pierre Auger Observatory,” JCAP 10 (2019) 022, arXiv:1906.07422.
  * [168] IceCube, M. G. Aartsen et al., “Differential limit on the extremely-high-energy cosmic neutrino flux in the presence of astrophysical background from nine years of IceCube data,” Phys. Rev. D 98 (2018) no. 6, 062003, arXiv:1807.01820.
  * [169] M. H. Reno, J. F. Krizmanic, and T. M. Venters, “Cosmic tau neutrino detection via Cherenkov signals from air showers from Earth-emerging taus,” Phys. Rev. D 100 (2019) no. 6, 063010, arXiv:1902.11287.
  * [170] L. A. Anchordoqui et al., “Performance and science reach of the Probe of Extreme Multimessenger Astrophysics for ultrahigh-energy particles,” Phys. Rev. D 101 (2020) no. 2, 023012, arXiv:1907.03694.
  * [171] A. Neronov, D. V. Semikoz, L. A. Anchordoqui, J. Adams, and A. V. Olinto, “Sensitivity of a proposed space-based Cherenkov astrophysical-neutrino telescope,” Phys. Rev. D 95 (2017) no. 2, 023004, arXiv:1606.03629.
  * [172] M. Beneke, A. P. Chapovsky, A. Signer, and G. Zanderighi, “Effective theory calculation of resonant high-energy scattering,” Nucl. Phys. B 686 (2004) 205–247, arXiv:hep-ph/0401002.
  * [173] G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, “Theory and phenomenology of two-Higgs-doublet models,” Phys. Rept. 516 (2012) 1–102, arXiv:1106.0034.
  * [174] G. Senjanovic and R. N. Mohapatra, “Exact Left-Right Symmetry and Spontaneous Violation of Parity,” Phys. Rev. D 12 (1975) 1502.
  * [175] M. Magg and C. Wetterich, “Neutrino Mass Problem and Gauge Hierarchy,” Phys. Lett. B 94 (1980) 61–64.
  * [176] J. Schechter and J. W. F. Valle, “Neutrino Masses in SU(2) x U(1) Theories,” Phys. Rev. D 22 (1980) 2227.
  * [177] T. P. Cheng and L.-F. Li, “Neutrino Masses, Mixings and Oscillations in SU(2) x U(1) Models of Electroweak Interactions,” Phys. Rev. D 22 (1980) 2860.
  * [178] G. Lazarides, Q. Shafi, and C. Wetterich, “Proton Lifetime and Fermion Masses in an SO(10) Model,” Nucl. Phys. B 181 (1981) 287–300.
  * [179] R. N. Mohapatra and G. Senjanovic, “Neutrino Masses and Mixings in Gauge Models with Spontaneous Parity Violation,” Phys. Rev. D 23 (1981) 165.
  * [180] K. S. Babu, P. S. B. Dev, S. Jana, and A. Thapa, “Non-Standard Interactions in Radiative Neutrino Mass Models,” arXiv:1907.09498.
  * [181] Y. Cai, J. Herrero-García, M. A. Schmidt, A. Vicente, and R. R. Volkas, “From the trees to the forest: a review of radiative neutrino mass models,” Front.in Phys. 5 (2017) 63, arXiv:1706.08524.
  * [182] C. Klein, M. Lindner, and S. Ohmer, “Minimal Radiative Neutrino Masses,” JHEP 03 (2019) 018, arXiv:1901.03225.
  * [183] S. Jana, P. K. Vishnu, and S. Saad, “Minimal realizations of Dirac neutrino mass from generic one-loop and two-loop topologies at $d=5$,” JCAP 04 (2020) 018, arXiv:1910.09537.
  * [184] A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31 (1980) 260.
  * [185] M. Dine, W. Fischler, and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B 104 (1981) 199–202.
  * [186] D. Abercrombie et al., “Dark Matter Benchmark Models for Early LHC Run-2 Searches: Report of the ATLAS/CMS Dark Matter Forum,” arXiv:1507.00966.
  * [187] A. Zee, “A Theory of Lepton Number Violation, Neutrino Majorana Mass, and Oscillation,” Phys. Lett. 93B (1980) 389. [Erratum: Phys. Lett.95B,461(1980)].
  * [188] J. C. Pati and A. Salam, “Unified Lepton-Hadron Symmetry and a Gauge Theory of the Basic Interactions,” Phys. Rev. D8 (1973) 1240–1251.
  * [189] J. C. Pati and A. Salam, “Lepton Number as the Fourth Color,” Phys. Rev. D10 (1974) 275–289. [Erratum: Phys. Rev.D11,703(1975)].
  * [190] H. Georgi and S. L. Glashow, “Unity of All Elementary Particle Forces,” Phys. Rev. Lett. 32 (1974) 438–441.
  * [191] H. Georgi, H. R. Quinn, and S. Weinberg, “Hierarchy of Interactions in Unified Gauge Theories,” Phys. Rev. Lett. 33 (1974) 451–454.
  * [192] H. Georgi, “The State of the Art-Gauge Theories,” AIP Conf. Proc. 23 (1975) 575–582.
  * [193] H. Fritzsch and P. Minkowski, “Unified Interactions of Leptons and Hadrons,” Annals Phys. 93 (1975) 193–266.
  * [194] K. D. Lane, “An Introduction to technicolor,” in Theoretical Advanced Study Institute (TASI 93) in Elementary Particle Physics: The Building Blocks of Creation - From Microfermius to Megaparsecs.  6, 1993.  arXiv:hep-ph/9401324.
  * [195] R. Barbier et al., “R-parity violating supersymmetry,” Phys. Rept. 420 (2005) 1–202, arXiv:hep-ph/0406039.
  * [196] Belle, A. Abdesselam et al., “Measurement of $\mathcal{R}(D)$ and $\mathcal{R}(D^{\ast})$ with a semileptonic tagging method,” arXiv:1904.08794.
  * [197] I. Bigaran and R. R. Volkas, “Getting chirality right: Single scalar leptoquark solutions to the $(g-2)_{e,\mu}$ puzzle,” Phys. Rev. D 102 (2020) no. 7, 075037, arXiv:2002.12544.
  * [198] K. S. Babu, P. S. B. Dev, S. Jana, and A. Thapa, “Unified framework for $B$-anomalies, muon $g-2$ and neutrino masses,” JHEP 03 (2021) 179, arXiv:2009.01771.
  * [199] I. Doršner, S. Fajfer, and S. Saad, “$\mu\to e\gamma$ selecting scalar leptoquark solutions for the $(g-2)_{e,\mu}$ puzzles,” Phys. Rev. D 102 (2020) no. 7, 075007, arXiv:2006.11624.
  * [200] I. Bischer and W. Rodejohann, “General neutrino interactions from an effective field theory perspective,” Nucl. Phys. B 947 (2019) 114746, arXiv:1905.08699.
  * [201] ATLAS, G. Aad et al., “Search for a heavy charged boson in events with a charged lepton and missing transverse momentum from $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” Phys. Rev. D 100 (2019) no. 5, 052013, arXiv:1906.05609.
  * [202] K. S. Babu, B. Dutta, and R. N. Mohapatra, “A theory of R(D∗, D) anomaly with right-handed currents,” JHEP 01 (2019) 168, arXiv:1811.04496.
  * [203] X. He, G. C. Joshi, H. Lew, and R. Volkas, “NEW Z-prime PHENOMENOLOGY,” Phys. Rev. D 43 (1991) 22–24.
  * [204] R. Foot, “New Physics From Electric Charge Quantization?,” Mod. Phys. Lett. A 6 (1991) 527–530.
  * [205] X.-G. He, G. C. Joshi, H. Lew, and R. R. Volkas, “Simplest Z-prime model,” Phys. Rev. D 44 (1991) 2118–2132.
  * [206] S. Baek, N. G. Deshpande, X. G. He, and P. Ko, “Muon anomalous g-2 and gauged $L_{\mu}-L_{\tau}$ models,” Phys. Rev. D 64 (2001) 055006, arXiv:hep-ph/0104141.
  * [207] A. Davidson, “$B-L$ as the fourth color within an $SU(2)_{L}\times U(1)_{R}\times U(1)$ model,” Phys. Rev. D20 (1979) 776.
  * [208] R. E. Marshak and R. N. Mohapatra, “Quark - Lepton Symmetry and B-L as the U(1) Generator of the Electroweak Symmetry Group,” Phys. Lett. 91B (1980) 222–224.
  * [209] F. Wilczek and A. Zee, “Conservation or Violation of $B^{-}$l in Proton Decay,” Phys. Lett. B 88 (1979) 311–314.
  * [210] R. N. Mohapatra and R. E. Marshak, “Local B-L Symmetry of Electroweak Interactions, Majorana Neutrinos and Neutron Oscillations,” Phys. Rev. Lett. 44 (1980) 1316–1319. [Erratum: Phys. Rev. Lett.44,1643(1980)].
  * [211] E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal, “Neutrino Masses and Mixings Dynamically Generated by a Light Dark Sector,” Phys. Lett. B791 (2019) 210–214, arXiv:1808.02500.
  * [212] E. Bertuzzo, S. Jana, P. A. N. Machado, and R. Zukanovich Funchal, “Dark Neutrino Portal to Explain MiniBooNE excess,” Phys. Rev. Lett. 121 (2018) no. 24, 241801, arXiv:1807.09877.
  * [213] M. Berbig, S. Jana, and A. Trautner, “The Hubble tension and a renormalizable model of gauged neutrino self-interactions,” Phys. Rev. D 102 (2020) no. 11, 115008, arXiv:2004.13039.
  * [214] ATLAS, G. Aad et al., “Search for high-mass dilepton resonances using 139 fb-1 of $pp$ collision data collected at $\sqrt{s}=$13 TeV with the ATLAS detector,” Phys. Lett. B 796 (2019) 68–87, arXiv:1903.06248.
  * [215] K. S. Babu, A. Friedland, P. A. N. Machado, and I. Mocioiu, “Flavor Gauge Models Below the Fermi Scale,” JHEP 12 (2017) 096, arXiv:1705.01822.
  * [216] F. Elahi and A. Martin, “LHC constraints on a $(B-L)_{3}$ gauge boson,” Phys. Rev. D 100 (2019) no. 3, 035016, arXiv:1905.10106.
  * [217] K. S. Babu, D. Gonçalves, S. Jana, and P. A. N. Machado, “Neutrino Non-Standard Interactions: Complementarity Between LHC and Oscillation Experiments,” Phys. Lett. B 815 (2021) 136131, arXiv:2003.03383.
  * [218] K. S. Babu and S. Jana, “Enhanced Di-Higgs Production in the Two Higgs Doublet Model,” JHEP 02 (2019) 193, arXiv:1812.11943.
  * [219] CMS, A. M. Sirunyan et al., “Search for dark matter produced in association with heavy-flavor quark pairs in proton-proton collisions at $\sqrt{s}=13$ TeV,” Eur. Phys. J. C 77 (2017) no. 12, 845, arXiv:1706.02581.
  * [220] ALEPH, DELPHI, L3, OPAL, LEP, G. Abbiendi et al., “Search for Charged Higgs bosons: Combined Results Using LEP Data,” Eur. Phys. J. C 73 (2013) 2463, arXiv:1301.6065.
  * [221] CMS, A. M. Sirunyan et al., “Search for supersymmetry in events with a $\tau$ lepton pair and missing transverse momentum in proton-proton collisions at $\sqrt{s}=$ 13 TeV,” JHEP 11 (2018) 151, arXiv:1807.02048.
  * [222] ATLAS, G. Aad et al., “Search for charged Higgs bosons decaying into a top quark and a bottom quark at $\sqrt{\mathrm{s}}$ = 13 TeV with the ATLAS detector,” JHEP 06 (2021) 145, arXiv:2102.10076.
  * [223] ATLAS, M. Aaboud et al., “Search for charged Higgs bosons decaying via $H^{\pm}\to\tau^{\pm}\nu_{\tau}$ in the $\tau$+jets and $\tau$+lepton final states with 36 fb-1 of $pp$ collision data recorded at $\sqrt{s}=13$ TeV with the ATLAS experiment,” JHEP 09 (2018) 139, arXiv:1807.07915.
  * [224] CMS, V. Khachatryan et al., “Search for a light charged Higgs boson decaying to $\mathrm{c}\overline{\mathrm{s}}$ in pp collisions at $\sqrt{s}=8$ TeV,” JHEP 12 (2015) 178, arXiv:1510.04252.
  * [225] J. Fuentes-Martin, A. Greljo, J. Martin Camalich, and J. D. Ruiz-Alvarez, “Charm physics confronts high-pT lepton tails,” JHEP 11 (2020) 080, arXiv:2003.12421.
  * [226] ATLAS, G. Aad et al., “Search for pairs of scalar leptoquarks decaying into quarks and electrons or muons in $\sqrt{s}$ = 13 TeV $pp$ collisions with the ATLAS detector,” JHEP 10 (2020) 112, arXiv:2006.05872.
  * [227] ATLAS, M. Aaboud et al., “Searches for third-generation scalar leptoquarks in $\sqrt{s}$ = 13 TeV pp collisions with the ATLAS detector,” JHEP 06 (2019) 144, arXiv:1902.08103.
  * [228] L. Buonocore, U. Haisch, P. Nason, F. Tramontano, and G. Zanderighi, “Lepton-Quark Collisions at the Large Hadron Collider,” Phys. Rev. Lett. 125 (2020) no. 23, 231804, arXiv:2005.06475.
  * [229] ATLAS, G. Aad et al., “Search for new non-resonant phenomena in high-mass dilepton final states with the ATLAS detector,” JHEP 11 (2020) 005, arXiv:2006.12946. [Erratum: JHEP 04, 142 (2021)].
  * [230] A. Angelescu, D. Bečirević, D. A. Faroughy, and O. Sumensari, “Closing the window on single leptoquark solutions to the $B$-physics anomalies,” JHEP 10 (2018) 183, arXiv:1808.08179.
  * [231] CMS, A. M. Sirunyan et al., “Constraints on models of scalar and vector leptoquarks decaying to a quark and a neutrino at $\sqrt{s}=$ 13 TeV,” Phys. Rev. D 98 (2018) no. 3, 032005, arXiv:1805.10228.
  * [232] CMS, S. Chatrchyan et al., “Identification of b-Quark Jets with the CMS Experiment,” JINST 8 (2013) P04013, arXiv:1211.4462.
  * [233] L. A. Anchordoqui, C. A. Garcia Canal, H. Goldberg, D. Gomez Dumm, and F. Halzen, “Probing leptoquark production at IceCube,” Phys. Rev. D 74 (2006) 125021, arXiv:hep-ph/0609214.
  * [234] U. K. Dey and S. Mohanty, “Constraints on Leptoquark Models from IceCube Data,” JHEP 04 (2016) 187, arXiv:1505.01037.
  * [235] V. Barger and W.-Y. Keung, “Superheavy Particle Origin of IceCube PeV Neutrino Events,” Phys. Lett. B 727 (2013) 190–193, arXiv:1305.6907.
  * [236] B. Dutta, Y. Gao, T. Li, C. Rott, and L. E. Strigari, “Leptoquark implication from the CMS and IceCube experiments,” Phys. Rev. D 91 (2015) 125015, arXiv:1505.00028.
  * [237] N. Mileo, A. de la Puente, and A. Szynkman, “Implications of a Electroweak Triplet Scalar Leptoquark on the Ultra-High Energy Neutrino Events at IceCube,” JHEP 11 (2016) 124, arXiv:1608.02529.
  * [238] P. S. B. Dev, D. K. Ghosh, and W. Rodejohann, “R-parity Violating Supersymmetry at IceCube,” Phys. Lett. B 762 (2016) 116–123, arXiv:1605.09743.
  * [239] U. K. Dey, D. Kar, M. Mitra, M. Spannowsky, and A. C. Vincent, “Searching for Leptoquarks at IceCube and the LHC,” Phys. Rev. D 98 (2018) no. 3, 035014, arXiv:1709.02009.
  * [240] B. Chauhan, B. Kindra, and A. Narang, “Discrepancies in simultaneous explanation of flavor anomalies and IceCube PeV events using leptoquarks,” Phys. Rev. D 97 (2018) no. 9, 095007, arXiv:1706.04598.
  * [241] D. Bečirević, B. Panes, O. Sumensari, and R. Zukanovich Funchal, “Seeking leptoquarks in IceCube,” JHEP 06 (2018) 032, arXiv:1803.10112.
  * [242] K. S. Babu, P. S. Dev, S. Jana, and Y. Sui, “Zee-Burst: A New Probe of Neutrino Nonstandard Interactions at IceCube,” Phys. Rev. Lett. 124 (2020) no. 4, 041805, arXiv:1908.02779.
  * [243] U. K. Dey, N. Nath, and S. Sadhukhan, “Charged Higgs effects in IceCube: PeV events and NSIs,” JHEP 09 (2021) 113, arXiv:2010.05797.
  * [244] (IceCube Collaboration)*, IceCube, R. Abbasi et al., “All-flavor constraints on nonstandard neutrino interactions and generalized matter potential with three years of IceCube DeepCore data,” Phys. Rev. D 104 (2021) no. 7, 072006, arXiv:2106.07755.
  * [245] D. B. Clark, E. Godat, and F. I. Olness, “ManeParse : A Mathematica reader for Parton Distribution Functions,” Comput. Phys. Commun. 216 (2017) 126–137, arXiv:1605.08012.
  * [246] J. Ellis, K. Sakurai, and M. Spannowsky, “Search for Sphalerons: IceCube vs. LHC,” JHEP 05 (2016) 085, arXiv:1603.06573.
  * [247] K. J. Golec-Biernat and M. Wusthoff, “Saturation effects in deep inelastic scattering at low Q**2 and its implications on diffraction,” Phys. Rev. D 59 (1998) 014017, arXiv:hep-ph/9807513.
  * [248] J. Bartels, K. J. Golec-Biernat, and H. Kowalski, “A modification of the saturation model: DGLAP evolution,” Phys. Rev. D 66 (2002) 014001, arXiv:hep-ph/0203258.
  * [249] E. Iancu, K. Itakura, and S. Munier, “Saturation and BFKL dynamics in the HERA data at small x,” Phys. Lett. B 590 (2004) 199–208, arXiv:hep-ph/0310338.
  * [250] N. Armesto, C. Merino, G. Parente, and E. Zas, “Charged current neutrino cross-section and tau energy loss at ultra-high energies,” Phys. Rev. D 77 (2008) 013001, arXiv:0709.4461.
  * [251] A. Y. Illarionov, B. A. Kniehl, and A. V. Kotikov, “Ultrahigh-energy neutrino-nucleon deep-inelastic scattering and the Froissart bound,” Phys. Rev. Lett. 106 (2011) 231802, arXiv:1105.2829.
  * [252] M. M. Block, L. Durand, and P. Ha, “Connection of the virtual $\gamma^{*}p$ cross section of ep deep inelastic scattering to real $\gamma$p scattering, and the implications for $\nu$N and ep total cross sections,” Phys. Rev. D 89 (2014) no. 9, 094027, arXiv:1404.4530.
  * [253] C. A. Argüelles, F. Halzen, L. Wille, M. Kroll, and M. H. Reno, “High-energy behavior of photon, neutrino, and proton cross sections,” Phys. Rev. D 92 (2015) no. 7, 074040, arXiv:1504.06639.
  * [254] D. Binosi and L. Theussl, “JaxoDraw: A Graphical user interface for drawing Feynman diagrams,” Comput. Phys. Commun. 161 (2004) 76–86, arXiv:hep-ph/0309015.
  * [255] J. Edsjö, “Calculation of neutrino cross sections and the nusigma neutrino-nucleon scattering Monte Carlo,”. http://wimpsim.astroparticle.se/code/nucross3.pdf.
  * [256] V. Bertone, R. Gauld, and J. Rojo, “Neutrino Telescopes as QCD Microscopes,” JHEP 01 (2019) 217, arXiv:1808.02034.
  * [257] R. M. Barnett, “Evidence in Neutrino Scattering for Righthanded Currents Associated with Heavy Quarks,” Phys. Rev. D 14 (1976) 70.