# Artificial proto-modelling with simplified-model results from the LHC

Sabine Kraml${}^{1}\,$111Speaker André Lessa2 Wolfgang Waltenberger3,4 1
Laboratoire de Physique Subatomique et de Cosmologie, Université Grenoble-
Alpes,
CNRS/IN2P3, 53 Avenue des Martyrs, F-38026 Grenoble, France
2 Centro de Ciências Naturais e Humanas, Universidade Federal do ABC, Santo
André, SP - Brasil
3 Institut für Hochenergiephysik, Österreichische Akademie der Wissenschaften,
Nikolsdorfer Gasse 18, A-1050 Wien, Austria
4 University of Vienna, Faculty of Physics, Boltzmanngasse 5, A-1090 Wien,
Austria

###### Abstract

We present a novel approach to identify potential dispersed signals of new
physics in the slew of published LHC results. It employs a random walk
algorithm to introduce sets of new particles, dubbed “proto-models”, which are
tested against simplified-model results from ATLAS and CMS searches for new
physics by exploiting the SModelS software framework. A combinatorial
algorithm identifies the set of analyses and/or signal regions that maximally
violates the Standard Model hypothesis, while remaining compatible with the
entirety of LHC constraints in our database. Crucial to the method is the
ability to construct a reliable likelihood in proto-model space; we explain
the various approximations which are needed depending on the information
available from the experiments, and how they impact the whole procedure.

##  1 Introduction

Searches for new physics at the LHC are usually pursued on a channel-by-
channel basis. A plethora of experimental analyses are thus performed in many
different final states, and the outcomes interpreted as limits on new
particles in the context of specific Beyond the Standard Model (BSM)
scenarios, again for each search channel separately. The disadvantage of this
kind of hypothesis testing—which is also done on the theory side by comparing
signal predictions of specific BSM incarnations to the experimental limits—is
that only a small part of the overall available data is used. Small effects of
dispersed signals, which will show up simultaneously in different final states
and/or signal regions, might easily be missed, or disregarded as statistical
fluctuations.

A more global exploration of the LHC data is in order. As an attempt in this
direction, we presented in [1] the prototype of a statistical learning
algorithm that identifies potential dispersed signals in the slew of published
LHC analyses, building candidate “proto-models ” from them, while remaining
compatible with the entirety of LHC results. Such proto-models may then be
scrutinised further in dedicated analyses and, in case of a discovery, help to
eventually unravel the concrete underlying BSM theory. The ultimate goal is a
data-driven bottom-up approach to the quest of new physics with minimal
theoretical bias. Being easily extendable, it should eventually also allow one
to fold in additional information from other (future) experiments to
continuously improve the picture of what data is telling us about BSM physics.

##  2 The walker algorithm

At present, our statistical learning algorithm is based on the concept of
simplified models, exploiting the SModelS [2, 3, 4, 5] software framework and
its large database of experimental results. It employs a Markov Chain Monte
Carlo (MCMC)-type random walk through proto-model space, adding or removing
new particles, and randomly changing their cross sections and branching
ratios. This is coupled to a combinatorial algorithm, which identifies the set
of analyses and signal regions that maximally violates the Standard Model (SM)
hypothesis.

The algorithm to comb through proto-model parameter space in order to identify
the models that best fit the data is dubbed the walker. It is composed of
several building blocks, or “machines” that interact with each other in a
well-defined fashion:

  1. 1.

Starting with the Standard Model, the builder creates proto-models, randomly
adding or removing new particles and changing any of the proto-model
parameters.

  2. 2.

The proto-model is then passed on to the critic, which checks the model
against the database of simplified-model results to determine an upper bound
on an overall signal strength ($\mu_{\rm max}$).

  3. 3.

The combiner identifies all possible combinations of results and constructs a
combined likelihood for each subset.

An essential aspect in our procedure is that that neither the number nor the
kind of the new particles is fixed. Instead, the BSM particle content, the
particle masses, production cross sections and decay branching ratios (BRs)
are taken as free parameters of the proto-models. This implies a parameter
space of varying dimensionality for the MCMC-type random walks!

##  3 Proto-model construction rules

Proto-models are not intended to be fully consistent theoretical models;
therefore their properties are not bound by higher-level theoretical
assumptions on the underlying BSM theory. Nonetheless, for practical purposes,
we have to make a number of assumptions.

Concretely, in this work we assume that all particles either decay promptly or
are fully stable at detector scales. Moreover, since we make extensive use of
simplified-model results from searches for supersymmetry (SUSY), we impose
that all BSM particles are odd under a $\mathcal{Z}_{2}$-type symmetry, so
they are always pair produced and always cascade decay to the lightest state.
The Lightest BSM Particle (LBP) is taken to be stable and electrically and
color neutral, and hence is a dark matter candidate. Finally, only particles
with masses within LHC reach are considered part of a specific proto-model.

In the current version of the algorithm, we allow proto-models to consist of
up to 20 BSM particles: light quark partners $X_{q}$ ($q=u,d,c,s$); heavy
quark partners $X_{b}^{i}$, $X_{t}^{i}$ ($i=1,2$); charged lepton and neutrino
partners $X_{\ell}$, $X_{\nu_{\ell}}$ ($\ell=e,\mu,\tau$); a color-octet gluon
partner $X_{g}$; and electroweak partners $X_{W}^{i}$, $X_{Z}^{j}$ ($i=1,2$;
$j=1,2,3$).222Here positively and negatively charged $X_{W}^{i}$ are counted
as one instance; they are taken as mass-degenerate and have the same decay
branching ratios, but their production cross-sections are free parameters. A
priori they may be same-spin or opposite-spin partners to the SM particles.

##  4 Simplified-model LHC results

In order to confront individual proto-models with a large number of LHC
results, we make use of SModelS. Not employing any MC event simulation, this
is computationally cheap, making it feasible to test hundreds of thousands of
proto-models within reasonable time. The experimental results stored in the
SModelS database fall into two main categories:

  * •

Upper Limit (UL) results: these are the 95% confidence level limits on the
production cross sections (times BR) for simplified model topologies as a
function of the BSM masses obtained by the experimental collaborations.

  * •

Efficiency Map (EM) results: these correspond to signal efficiencies (more
precisely, acceptance $\times$ efficiency, ${\cal{A}}\times\epsilon$, values)
for simplified topologies as a function of the BSM masses for the signal
regions considered by the corresponding experimental analysis.

The SModelS database v1.2.4 used here includes results from 40 ATLAS and 46
CMS experimental searches at $\sqrt{s}=8$ and 13 TeV, corresponding to about
250 upper limit maps and 1,700 individual efficiency maps; see [1] for
details.

##  5 Construction of a global likelihood

Key to our procedure is the construction of an approximate global likelihood
for the signal,

$\mathrm{L}_{\mathbf{BSM}}(\mu|D)=P\left(D|\mu s+b+\theta\right)p(\theta)\,,$
(1)

which describes the plausibility of the signal strength $\mu$, given the data
$D$. Here, $\theta$ denotes the nuisance parameters describing systematic
uncertainties in the signal ($s$) and background ($b$), while $p(\theta)$
corresponds to their probability distribution function.333Note that $s$, $b$
and $\theta$ in Eq. (1) are multi-dimensional quantities. The principle which
we follow for combining individual likelihoods to a global one is illustrated
in Figure 1.

Figure 1: Left: an illustration of the signal likelihood as a function of the
signal strength ($\mu$) for two results (A and B) and the corresponding
combined likelihood. The value of $\hat{\mu}$ given by the maximum of the
combined likelihood (green curve) is also shown as well as the upper limit on
$\mu$ ($\mu_{max}$) given by Result B alone. Right: a similar example, but
where it is not possible to combine both likelihoods, or the likelihood for
Result B is not available. In this case the likelihood for Result A is
truncated at $\mu_{max}$ (obtained from Result B). All the curves are
normalized to the corresponding SM likelihood value ($\mu=0$).

###  5.1 Likelihoods for individual analyses

The extent to which we can compute likelihoods for individual analyses
crucially depends on the information available from the experimental
collaborations. Simplified-model EMs allow us to determine the expected number
of events for the hypothesised proto-model. Together with the number of
observed events, the expected backgrounds, and the uncertainties thereon, we
can then construct a simplified likelihood by assuming a Gaussian distribution
for the uncertainties and a Poissonian for the data; the nuisances are
profiled over. This is a priori done per signal region. For some CMS analyses,
a covariance matrix is available, which allows one to combine different signal
regions, still in a simplified likelihood approach [6]. ATLAS has recently
started to provide the full statistical models [7] for some analyses, with
which one can compute the likelihood at nearly 444Essentially up to small
differences resulting from uncertainties in the EMs. the same fidelity as in
the experiment.

If only ULs are available for an analysis, things look less good. If the
expected ULs are available in addition to the observed ones, we can describe
the likelihood as a truncated Gaussian, although this is often a crude
approximation.[1] However, if only observed ULs are available, it is simply
not possible to construct a reasonable likelihood function; such results can
only be used to determine $\mu_{\rm max}$ for the critic.

###  5.2 Combining analyses

In principle, to combine the likelihoods from individual analyses, we would
need to know their correlations. Lacking this information, we treat cross-
analyses correlations in an approximate, binary way—a given pair of analyses
is either considered to be approximately uncorrelated, in which case it may be
combined (by multiplying the respective likelihoods), or it is not at all
considered for combination. We always assume results from different LHC runs
and/or from different experiments to be approximately uncorrelated. Moreover,
we also treat results with clearly different final states (e.g., fully
hadronic final states vs. final states with leptons) to be uncorrelated.
Figure 2 illustrates this for the Run 2 analyses considered in this work.

Figure 2: Binary correlation matrices for the 23 ATLAS (left) and 28 CMS
(right) Run 2 analyses considered in this work, showing which ones are taken
as approximately uncorrelated (green bins) and which are not (red bins).
Analyses from different LHC runs or different experiments are always taken as
uncorrelated. White bins denote analyses, for which no likelihood can be
computed because only observed ULs are available.

With these assumptions, we construct an approximate combined likelihood for
subsets of LHC results,
$\mathrm{L}_{\mathbf{BSM}}(\mu)=\prod_{i=1}^{n}\mathrm{L}_{i}(\mu)$, where the
product is over all $n$ uncorrelated analyses and $\mu$ is the global signal
strength. We refer to such subsets as combinations of results. Two important
rules apply: 1. any pair of results in the subset must be considered as
uncorrelated; and 2. any result which is allowed to be added to the
combination, must be added. Information from all the other analyses, which are
not included in the combination, is accounted for as a constraint on the
global signal strength $\mu$, i.e. $\mu_{\rm max}$.

A comment is in order at this point. Indeed, while a good number of analyses
can be combined in the way described above, we also see from the white bins in
Fig. 2, that about half of the Run 2 analyses in the SModelS v1.2.4 database
consist of observed ULs only, and thus no proper likelihood can be computed
for them. For Run 1, this concerns about one third of all analyses. Clearly,
more EM-type results would be welcome to improve this situation.

##  6 Test statistics

Last but not least, in order to guide the MCMC-type walk and identify the
proto-models that best fit the data, we need a test statistic $K$. While $K$
should increase for models which better satisfy all the constraints (which
includes better fitting potential dispersed signals), it is also desirable to
enforce the law of parsimony by reducing the test statistic of models with too
many degrees of freedom. We define the test statistic as

$K:=\max\limits_{\forall c\in C}K^{c}\qquad{\rm with}\qquad
K^{c}=2\ln\frac{\mathrm{L}^{c}_{\mathbf{BSM}}(\hat{\mu})\cdot\pi(\mathbf{BSM})}{\mathrm{L}^{c}_{\mathbf{SM}}\cdot\pi(\mathbf{SM})}\,.$
(2)

Here $\mathrm{L}^{c}_{\mathbf{BSM}}$ is the likelihood for a combination $c$
of experimental results given the proto-model, evaluated at the signal
strength value $\hat{\mu}$, which maximizes the likelihood and satisfies
$0\leq\hat{\mu}<\mu_{max}$. $\mathrm{L}^{c}_{\mathbf{SM}}$ is the
corresponding SM likelihood, given by $\mathrm{L}^{c}_{\mathbf{BSM}}(\mu=0)$.
Finally, $\pi(\mathbf{SM})$ and $\pi(\mathbf{BSM})$ denote respectively the
priors for the SM and the proto-model. The total set of combinations of
results, $C$, is determined as explained in Section 5.2.

The proto-model prior $\pi(\mathbf{BSM})$ should penalize the test statistic
for newly introduced particles, branching ratios, or signal strength
multipliers, while $\pi(\mathbf{SM})\equiv 1$. Here, we choose

$\pi(\mathbf{M})=\exp\left[-\left(\frac{n_{\mathrm{particles}}}{2}+\frac{n_{\mathrm{BRs}}}{4}+\frac{n_{\mathrm{SSMs}}}{8}\right)\right]\,,$
(3)

where $n_{\mathrm{particles}}$ is the number of new particles present in the
proto-model, $n_{\mathrm{BRs}}$ is the number of non-trivial branching
ratios,555This means, for only one decay mode with 100% BR,
$n_{\mathrm{BRs}}=0$. and $n_{\mathrm{SSMs}}$ the number of signal strength
multipliers. This way, one particle with one non-trivial decay and two
production modes is equivalent to one free parameter in the Akaike Information
Criterion. For the SM,
$n_{\mathrm{particles}}=n_{\mathrm{BRs}}=n_{\mathrm{ssm}}=0$, and
$\pi(\mathbf{SM})\equiv 1$. The test statistic thus roughly corresponds to a
$\Delta\chi^{2}$ of the proto-model with respect to the SM, with a penalty for
the new degrees of freedom.

##  7 Results

For first physics results, we performed 10 runs of the walker algorithm, each
employing $50$ walkers and $1,000$ steps/walker. Figure 3 displays the mass
spectra of the proto-models with the highest $K$ value from each run. Besides
the $X_{Z}^{1}$ as the LBP, all models include one top-partner, $X_{t}^{1}$,
and one light-flavor quark partner, $X_{d,c}$, and their test statistics are
at $K=6.76\pm 0.08$ thus showing the stability of the algorithm.

Figure 3: Particle content and masses for the proto-models with highest test-
statistic $K$ obtained in each of the 10 runs performed over the real
database. The corresponding $K$ values are shown at the top. Note that for all
practical purposes $X_{d}$ and $X_{c}$ are indistinguishable.

For concreteness, let us take a closer look at the proto-model with $K=6.9$
from run 9.666Note that the scenarios in Fig. 3 are driven by the same data
and are statistically all equally plausible. The $X_{d}$ with mass $\approx
700$ GeV is introduced in order to fit small excesses in the
multijet+$E_{T}^{\rm miss}$ analyses from ATLAS at $\sqrt{s}=8$ and 13 TeV [8,
9] and from CMS at $\sqrt{s}=8$ TeV [10]. The $X_{t}^{1}$ with a mass of about
$1.2$ TeV is introduced to fit the $\approx 1.5\,\sigma$ and $2\,\sigma$
excesses observed in the 13 TeV CMS and ATLAS stop searches [11, 12]; indeed
this is what mostly drives the $K$ value. Despite corresponding to small
excesses, identifying the presence of such potential dispersed signals is one
of the main goals of the algorithm presented here. The fact that these
excesses appear in distinct ATLAS and CMS analyses and can be explained by the
introduction of a single top partner (actually with SUSY-like cross sections)
is another interesting outcome.

Finally, we can also compute a global $p$-value for the SM hypothesis by
running over synthetic “background-only” data. In the present setup, we find
that $p({\rm SM})\approx 0.19$; see [1] for details.

##  8 Conclusions

In view of the null results (so far) in the numerous channel-by-channel
searches for new particles, it becomes increasingly relevant to change
perspective and attempt a more global approach to find out where BSM physics
may hide. To this end, we presented a novel statistical learning algorithm
that is capable of identifying potential dispersed signals in the slew of
published LHC analyses. The task of the algorithm is to build candidate proto-
models from small excesses in the data, while at the same time remaining
consistent with all other constraints. At present, this is based on the
concept of simplified models, exploiting the SModelS software framework and
its large database of simplified-model results from ATLAS and CMS searches for
new physics.

Acknowledgments: S.K. was supported in part by the IN2P3 project “Théorie –
BSMGA”.

## References

## References

  * [1] W. Waltenberger, A. Lessa, S. Kraml, _Artificial Proto-Modelling: Building Precursors of a Next Standard Model from Simplified Model Results_ , _JHEP_ 03 (2021) 207 [2012.12246].
  * [2] S. Kraml et al., _SModelS: a tool for interpreting simplified-model results from the LHC and its application to supersymmetry_ , _Eur. Phys. J. C_ 74 (2014) 2868 [1312.4175].
  * [3] F. Ambrogi et al., _SModelS v1.1 user manual: Improving simplified model constraints with efficiency maps_ , _Comput. Phys. Commun._ 227 (2018) 72 [1701.06586].
  * [4] F. Ambrogi et al., _SModelS v1.2: long-lived particles, combination of signal regions, and other novelties_ , _Comput. Phys. Commun._ 251 (2020) 106848 [1811.10624].
  * [5] G. Alguero et al., _New developments in SModelS_ , _PoS_ TOOLS2020, 022 [2012.08192].
  * [6] CMS collaboration, _Simplified likelihood for the re-interpretation of public CMS results_ , Tech. Rep. CMS-NOTE-2017-001.
  * [7] ATLAS Collaboration, _Reproducing searches for new physics with the ATLAS experiment through publication of full statistical likelihoods_ , Tech. Rep. ATL-PHYS-PUB-2019-029.
  * [8] ATLAS collaboration, _Search for squarks and gluinos with the ATLAS detector in final states with jets and missing transverse momentum using $\sqrt{s}=8$ TeV proton–proton collision data_, _JHEP_ 09 (2014) 176 [1405.7875].
  * [9] ATLAS collaboration, _Search for squarks and gluinos in final states with jets and missing transverse momentum using 36 fb -1 of $\sqrt{s}=13$ TeV $pp$ collision data with the ATLAS detector_, _Phys. Rev. D_ 97 (2018) 112001 [1712.02332].
  * [10] CMS collab., _Search for new physics in the multijet and missing transverse momentum final state in proton-proton collisions at $\sqrt{s}=8$ TeV_, _JHEP_ 06 (2014) 055 [1402.4770].
  * [11] CMS collaboration, _Search for supersymmetry in proton-proton collisions at 13 TeV using identified top quarks_ , _Phys. Rev. D_ 97 (2018) 012007 [1710.11188].
  * [12] ATLAS collaboration, _Search for top-squark pair production in final states with one lepton, jets, and missing transverse momentum using 36 fb -1 of $\sqrt{s}=13$ TeV $pp$ collision data with the ATLAS detector_, _JHEP_ 06 (2018) 108 [1711.11520].


]https://polly.phys.msu.ru/en/labs/Tribelsky/

# Fall of Quantum Particle to the Center: Exact solution

Michael I. Tribelsky  E-mail: mitribel_at_gmail.com [ M. V. Lomonosov Moscow
State University, Moscow, 119991, Russia  National Research Nuclear University
MEPhI (Moscow Engineering Physics Institute), Moscow, 115409, Russia

###### Abstract

A fall of a particle to the center of a singular potential is one of a few
fundamental problems of quantum mechanics. Nonetheless, its solution is not
complete yet. The known results just indicate that if the singularity of the
potential is strong enough, the spectrum of the Schrodinger equation is not
bounded from below. However, the wave functions of the problem do not admit
the limiting transition to the ground state. Therefore, the unboundedness of
the spectrum is only a necessary condition. To prove that a quantum particle
indeed can fall to the center, a wave function describing the fall should be
obtained explicitly. This is done in the present paper. Specifically, an exact
solution of the time-dependent Schrodinger equation corresponding to the fall
is obtained and analyzed. A law for the collapse of the region of the wave
function localization to a single point is obtained explicitly. It is shown
that the known necessary conditions for the particle to fall simultaneously
are sufficient.

Introduction. It is impossible to overrate the importance of the Schrödinger
equation (SE). To say nothing that it is the keystone of entire non-
relativistic quantum mechanics, it arises in many other quite classical
problems, for example, in the stability analysis of solutions of various
nonlinear problems in combustion [1], optical damage of glass initiated by
small inclusions [2], Gunn domains in semiconductors [3, 4], structure
stability of the interphase boundary in liquid crystals [5], etc. Regarding
the analogy between the quantum scattering of a particle and the light
scattering by an obstacle, this is a well-known fact.

Among the vast diversity of problems related to SE and its applications, there
are several revealing its fundamental properties. They have paramount
importance and enter into all main textbooks on quantum mechanics. The problem
of a “fall” of a particle to the origin of a singular potential, also known as
quantum collapse, is the one of such a kind, see, e.g., Ref. [6]. Nowadays, in
this problem, the interest is shifted to its various generalizations [7, 8, 9,
10, 11, 12, 13, 14]. However, the critical question of whether the collapse
indeed may happen actually remains open even in the simplest case of a single
quantum particle in a centrally symmetric potential; see below. The goal of
this paper is to fill the gap by obtaining an exact solution to this problem.

Preliminary. In classical mechanics, the problem of a fall to the origin of a
particle moving in a centrally symmetric potential $U(r)$ is trivial since the
governing equations are exactly integrable at any $U(r)$ [15]. The integration
shows that the particle may pass through the origin, provided

$[r^{2}U(r)]_{r\rightarrow 0}<-{M^{2}}/(2m),$ (1)

where $M$ is the angular momentum and $m$ stands for the particle mass [15].

In quantum mechanics, the case is more tricky. Firstly, just a few potentials
make SE exactly integrable. Secondly, a particle’s localization in a given
space point must not contradict the uncertainty relations. The latter implies
a divergence of the uncertainty of the particle momentum if the localization
occurs.

Estimates of the spectrum of SE equation with

$U(r)=-\beta/r^{s};\;\;(\beta>0)$ (2)

reveal that at $s>2$, there are bound states with energy $E<0$, which is not
limited from below. This indicates a possibility of the collapse. At $s<2$,
the spectrum is bounded from below, and the collapse cannot happen. The case
$s=2$ requires more accurate consideration [6].

It should be stressed that potential (2) with $s\geq 2$ does exist in nature.
For example, the potential singular as $1/r^{2}$ describes the interaction of
a point electric charge with a particle with zero total electric charge but a
fixed finite dipole moment [8], to say nothing about various applications of
SE in non-quantum problems. Thus, in addition to purely academic interest, the
quantum collapse may have a pretty practical meaning.

Note now that while the spectrum bounded from below does mean that the fall to
the center is not possible, the unboundedness of the spectrum from below is
just the necessary condition for the fall to happen. To make the fall sure,
one has to build the corresponding wave function explicitly. The latter is
hardly possible with the help of the stationary-state eigenfunctions since for
the given problem they do not have any limit at $r\rightarrow 0$. This
behavior prevents one from implementing a direct limiting transition to the
ground state with $E\rightarrow-\infty$.

On top of that, even if the expression for the stationary-state eigenfunction
with $E=-\infty$ were obtained, it would not mean yet that the fall can indeed
occur. To prove the occurrence of the fall, one has to show that there is at
least a single initial state whose wave function has a non-zero projection to
this eigenfunction, i.e., that the state corresponding to the localization of
the particle at $r=0$ may be indeed excited. Last but not least is the
question about the law describing the contraction with the time of the region
of the wave function spatial localization if the fall does happen. Thus,
strictly speaking, the question about a quantum particle fall to the origin of
a centrally symmetric potential remains open.

The difficulties mentioned above have arisen owing to the attempt to describe
an essentially time-dependent process of the fall with the help of the
stationary-state eigenfunctions. In a sense, it looks like trying to force a
square peg into a round hole. To save the day, one has to consider a
spatiotemporal evolution of a wave package, and this consideration should not
be based on the projection of the corresponding wave function onto the
stationary-state eigenmodes.

Here I should confess that the idea is not mine. When I was a junior
scientist, my adviser Dr. Sergei I. Anisimov [16] from Landau Institute, told
me about this approach to the problem. Moreover, he said that together with
Dr. Igor E. Dzyaloshinskii [17], they had found a solution describing the
collapse of the wave package.

Many years later, in connection with light scattering by nanoparticles, I came
across the problem mathematically analogous to the quantum collapse. I
remembered this conversation with Dr. Anisimov and asked him for details and
references. The reply was that these results had never been published and
details he did not remember. Now, when both Dr. Anisimov and Dr.
Dzyalashinskii passed away, I decide to reobtain their results and make them
available to a broad readership as a small token of my great respect to these
distinguished scholars.

Problem formulation. Thus, we have to start from the full time-dependent
version of SE:

$i\hbar\frac{\partial\Psi}{\partial
t}=\hat{H}\Psi;\;\;\hat{H}\equiv-\frac{\hbar^{2}}{2m}\Delta+U(r),$ (3)

where $\Delta$ stands for the Laplasian. Eq. (3) should be supplemented by the
initial condition $\Psi(\mathbf{r},t)=\Psi_{0}(\mathbf{r})$ and the standard
boundary conditions stipulating convergence of $\int|\Psi|^{2}d^{3}r$ for the
states corresponding to the bounded in space motion.

Suppose that the fall does take place. It means the collapse of the region of
the wave function localization to a point. Then, in the course of the
collapse, the wave function “forgets” the characteristic spatial scale of
$\Psi_{0}(\mathbf{r})$, if any. In other words, at the final stage of the
collapse, the dynamic of the wave function may become self-similar and hence
universal, i.e., independent of the initial condition, provided the latter
gives rise to the collapse.

Then, we have to check whether Eq. (3) has self-similar solutions and, if so,
to find them explicitly. To prove the possibility of the fall, obtaining just
a single solution of such a kind is enough.

Self-similar problem. We count the time down from the moment of the complete
collapse so that the collapse corresponds to $t=0$, and its dynamics is
described by $-\infty<t<0$. Once again, we consider the potential (2). Let us
look for a self-similar solution in the form

$\Psi=R(\xi)Y_{\ell}^{m}(\theta,\varphi).$ (4)

Here $Y_{\ell}^{m}(\theta,\varphi)$ is the spherical harmonic function and
$\xi=r/(-\chi t)^{\nu}$, where $\chi$ and $\nu$ are unknown yet constants.
Substitution of Eq. (4) into Eq. (3) indicates that the resulting equation for
$R(\xi)$ is reduced to an ODE depending on the single self-similar variable
$\xi$ solely if and only if

$s=2,\;\nu=1/2,\;\chi=const\cdot\hbar/m.$ (5)

The dimensionless constant in the expression for $\chi$ may be any. In what
follows, it is convenient to suppose it equals unity. Then, the corresponding
governing equation for $R(\xi)$ reads:

$R^{\prime\prime}+\left(\frac{2}{\xi}+i\xi\right)R^{\prime}+\frac{\gamma
R}{\xi^{2}}=0,$ (6)

where $\gamma\equiv\beta-\ell(\ell+1)$ and prime denotes $d/d\xi$.

Figure 1: $R(\xi)$ in several characteristic ranges of $\xi$; $\gamma=1/2$
(1), $\gamma=1$ (2), and $\gamma=2$ (3).

Eq. (4) is exactly integrable. Its solution satisfying the condition
$R(\xi)\rightarrow 0$ at $\xi\rightarrow\infty$ is

$\displaystyle R(\xi)$ $\displaystyle=$
$\displaystyle\frac{C}{\sqrt{\xi}}\left[{\xi}^{-\frac{i\alpha}{2}}{2}^{\frac{i\alpha}{2}}\Gamma\left(1+\frac{i\alpha}{2}\right)\Gamma\left(\frac{5-i\alpha}{4}\right)\\!\,_{1}F_{1}\left(-\frac{1+i\alpha}{4};\frac{2-i\alpha}{2};-\frac{i\xi^{2}}{2}\right)\right.$
$\displaystyle-$
$\displaystyle\left.\xi^{\frac{i\alpha}{2}}\exp\left(-\frac{\pi\alpha}{4}\right)\Gamma\left(1-\frac{i\alpha}{2}\right)\Gamma\left(\frac{5+i\alpha}{4}\right)\\!\,_{1}F_{1}\left(-\frac{1-i\alpha}{4};\frac{2+i\alpha}{2};-\frac{i\xi^{2}}{2}\right)\right].$

Here $C$ is a constant of integration, $\Gamma(z)$ stands for the Euler gamma
function, ${}_{1}F_{1}(a,b,z)$ designate the Kummer confluent hypergeometric
function of the first kind [18], and $\alpha\equiv\sqrt{4\gamma-1}>0$. The
positiveness of $\alpha$ follows from the fact that for $U(r)=\beta/r^{2}$ the
necessary conditions for the fall to occur derived from the unboundedness of
the spectrum from below stipulate that $\gamma>1/4$ [6].

Eq. (Fall of Quantum Particle to the Center: Exact solution) requires a
comment: It includes real positive quantities 2 and $\xi$ in imaginary power
$\pm i\alpha/2$. These expressions may be transformed as follows:

$\displaystyle z^{\pm i\frac{\alpha}{2}}$ $\displaystyle=$
$\displaystyle(e^{i2\pi n+\log z})^{\pm i\frac{\alpha}{2}}=e^{\mp\pi\alpha
n}e^{i\frac{\alpha}{2}\log z}$ $\displaystyle=$ $\displaystyle e^{\mp\pi\alpha
n}\left[\cos\left(\frac{\alpha}{2}\log z\right)\mp
i\sin\left(\frac{\alpha}{2}\log z\right)\right],$

where $n$ is arbitrary integer and $z>0$.

Eq. (Fall of Quantum Particle to the Center: Exact solution) has an infinite
number of branches corresponding to different values of $n$. At any $n$, all
these branches have the same phase equal to $({\alpha}/{2})\log z$. Therefore,
the real and imaginary parts of any branch have a number of zeros which in the
same manner increases without limit as $z\rightarrow 0$. However, the modulus
of these branches varies for different $n$ as $\exp({\mp\pi\alpha n})$.

For the sake of simplicity, anywhere in Eq. (Fall of Quantum Particle to the
Center: Exact solution) and in what follows, only the single branch with $n=0$
is selected. Just for this branch Eq. (Fall of Quantum Particle to the Center:
Exact solution) satisfies the condition $R(\xi)\rightarrow 0$ at
$\xi\rightarrow\infty$. For other branches to fulfill this condition, the
prefactor at the Kummer functions in Eq. (Fall of Quantum Particle to the
Center: Exact solution) should be changed accordingly.

Eq. (Fall of Quantum Particle to the Center: Exact solution) is a solution to
the problem, provided the obtained wave function may be normalized, i.e.,
$\int_{0}^{\infty}|R(\xi)|^{2}\xi^{2}d\xi$ converges. Based on the explicit
form of $R(\xi)$ given by Eq. (Fall of Quantum Particle to the Center: Exact
solution) one may select just two “dangerous” areas, which might result in
divergence of the integral, namely the asymptotics at $\xi\rightarrow\infty$
and at $\xi\rightarrow 0$.

Analysis of the behavior of the integrand at $\xi\rightarrow\infty$ is
relatively trivial. Employing the known asymptotical expression for
${}_{1}F_{1}(a,b,z)$ [18] one readily obtains that

$R(\xi)=C_{\infty}\frac{\exp(i\xi^{2}/2)}{\xi^{3}}+\ldots,$ (9)

where $C_{\infty}$ is a constant and ellipsis indicates dropped higher order
in $1/\xi$ terms. Thus, at the upper limit
$\int_{0}^{\infty}|R(\xi)|^{2}\xi^{2}d\xi$ converges.

The case $\xi\rightarrow 0$ is more tricky owing to the singularity described
by Eq. (Fall of Quantum Particle to the Center: Exact solution) . The
asymptotical behavior of Eq. (6) at $\xi\rightarrow 0$ is

$\displaystyle R(\xi)$ $\displaystyle=$
$\displaystyle\frac{C}{\sqrt{\xi}}\left[{\xi}^{-\frac{i\alpha}{2}}{2}^{\frac{i\alpha}{2}}\Gamma\left(1+\frac{i\alpha}{2}\right)\Gamma\left(\frac{5-i\alpha}{4}\right)\right.$
$\displaystyle-$
$\displaystyle\left.\xi^{\frac{i\alpha}{2}}\exp\left(-\frac{\pi\alpha}{4}\right)\Gamma\left(1-\frac{i\alpha}{2}\right)\Gamma\left(\frac{5+i\alpha}{4}\right)\right].$

Importantly, however, that according to Eqs. (Fall of Quantum Particle to the
Center: Exact solution), (Fall of Quantum Particle to the Center: Exact
solution) $|R(\xi)|^{2}$ in the vicinity of $\xi=0$ is majorized by
$C_{0}/\xi$, where $C_{0}$ is a constant. Therefore, the singularity at
$\xi=0$ is integrable, and $\int_{0}^{\infty}|R(\xi)|^{2}\xi^{2}d\xi$
converges at the lower limit too. In other words, at any finite $t<0$ the
quantity $\langle\Psi|\Psi\rangle$ remains finite.

As an example, plots $|R(\xi)|^{2}$ at several values of $\gamma$ are
presented in Fig. 1. It is interesting to note that despite Re$\,R(\xi)$ and
Im$\,R(\xi)$ both exhibit an infinite number of oscillations as $\xi$ tends to
zero, this is not the case for $|R(\xi)|^{2}$, which is a monotonically
decaying function of $\xi$.

Discussion. The obtained $\Psi$ makes it possible to calculate the average
value of any operator $\hat{A}$:
$\langle\hat{A}\rangle=\langle\Psi|\hat{A}|\Psi\rangle/\langle\Psi|\Psi\rangle$,
provided the integral $\langle\Psi|\hat{A}|\Psi\rangle$ converges, i.e.,
$\langle\hat{A}\rangle$ exists. However, it should be remembered that
$\Psi=R(\xi)Y_{\ell}^{m}(\theta,\varphi)$, while in scalar products the
integration is taken over $\mathbf{r}$-space, which requires the corresponding
replacement of the variable. Because of that the scalar products become time-
dependent. Then, one readily obtains that

$\langle\hat{r}\rangle=C_{r}\sqrt{-\chi t},\;\;\langle\hat{p}\rangle=\hbar
C_{p}/\sqrt{-\chi t},$ (11)

(remember that $t<0$). Here $\hat{p}$ stands for modulus of the momentum and
$C_{r,p}$ are dimensionless constants of the order of unity. Thus, the size of
the spatial region of the particle localization contracts as $\sqrt{-t}$ while
its momentum diverges as $1/\sqrt{-t}$ so that $\Delta p\Delta r$ remains a
constant of the order of $\hbar$, in agreement with the uncertainty relations.

Lets us now calculate the average energy of the falling particle 111The
discussed quantum state is not stationary; its energy is not an eigenvalue of
the Hamiltonian. Therefore for this state only average energy is meaningful.,
namely $E=\langle\hat{H}\rangle$. According to Eq. (3)

$E=\langle\Psi|\hat{H}|\Psi\rangle/\langle\Psi|\Psi\rangle=i\hbar\langle\Psi|\partial\Psi/\partial
t\rangle/\langle\Psi|\Psi\rangle.$ (12)

Bearing in mind that for the self-similar $\Psi$, given by Eq. (4)

$\frac{\partial}{\partial
t}=\frac{\partial}{\partial\xi}\frac{\partial\xi}{\partial
t}=-\frac{\xi}{2t}\frac{\partial}{\partial\xi}$

and employing the asymptotics (9), (Fall of Quantum Particle to the Center:
Exact solution), it is easy to see that at any $t<0$ integral
$\langle\Psi|\partial\Psi/\partial t\rangle$ converges. Seemingly, this gives
rise to a paradoxical conclusion that $E$ is not conserved and varies in time
as $1/t$.

However, this is not the case. The point is that for the given quantum state
$E=0$. Owing to the very complicated structure of Eq. (6) it is rather
difficult to prove the vanishing of $E$ by the direct calculation of the
corresponding integrals. Much more simple proof may be done based on the
general properties of SE. Indeed, suppose in a certain arbitrary moment of
time $E$ has a finite negative value. Since $U(r)\rightarrow 0$ at
$r\rightarrow\infty$ it means that at $r>r_{0}$, where $r_{0}$ satisfies the
condition $U(r_{0})=E$, the classically inaccessible region lies. In this
case, at $r\gg r_{0}$, the wave function must exponentially decay. However,
Eq. (9) does not exhibit an exponentially decaying asymptotic at
$\xi\rightarrow\infty$. Similarly, if $E$ is a finite positive quantity, the
asymptotic at $\xi\rightarrow\infty$ should correspond to a traveling plane
wave, which also disagrees with Eq. (9). Thus, the only remaining option is
$E=0$.

Conclusions. The performed study gives rise to the following outcome. The fall
of a quantum particle to the center does take place. In the potential
$U(r)=-\beta/r^{2}$, it occurs under the necessary conditions formulated in
Ref. [6]. It means that these conditions simultaneously are sufficient. The
fall is described by the self-similar solution of SE given by Eqs. (4), (Fall
of Quantum Particle to the Center: Exact solution). In this case, the radius
of the wave function localization region collapses as $\sqrt{-\chi t}$, while
the characteristic value of the momentum increases as $1/\sqrt{-\chi t}$. Here
$t<0$ is the time counted down from the moment of the full collapse. Finally,
note that since SE is invariant for time reversal accompanied with complex
conjugation, these transformations applied to Eqs. (4), (Fall of Quantum
Particle to the Center: Exact solution) generate the wave function describing
the escape of the particle from the center at $t>0$.

###### Acknowledgements.

Acknowledgements. This work is supported by the Russian Foundation for Basic
Research (Projects No. 20-02-00086) for the analytical study, the Russian
Science Foundation (Project No. 21-12-00151) for the symbolic computer
calculations, the Moscow Engineering Physics Institute Academic Excellence
Project (agreement with the Ministry of Education and Science of the Russian
Federation of 27 August 2013, Project No. 02.a03.21.0005) for the package of
computer graphic and the Russian Science Foundation (Project No. 19-72-30012)
for the provision of user facilities.

## References

  * Zeldovich and Barenblatt [1959] Y. B. Zeldovich and G. I. Barenblatt, Combustion and flame 3, 61 (1959).
  * Anisimov and Makshantsev [1973] S. Anisimov and B. Makshantsev, Sov. Phys. Solid State 15, 743 (1973).
  * Knight and Peterson [1967] B. Knight and G. Peterson, Physical Review 155, 393 (1967).
  * Tribelsky [2021] M. I. Tribelsky,  (2021),
10.20944/preprints202107.0347.v2.

  * Kats and Tribel’skii [1990] E. Kats and M. Tribel’skii, Sov. Phys. JETP 71, 108 (1990).
  * Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory, §18, §35_ (Elsevier, 2013).
  * Camblong  _et al._ [2000] H. E. Camblong, L. N. Epele, H. Fanchiotti, and C. A. García Canal, Phys. Rev. Lett. 85, 1590 (2000).
  * Sakaguchi and Malomed [2011a] H. Sakaguchi and B. A. Malomed, Physical Review A 83, 013607 (2011a).
  * Sakaguchi and Malomed [2011b] H. Sakaguchi and B. A. Malomed, Physical Review A 84, 033616 (2011b).
  * Sakaguchi and Malomed [2013] H. Sakaguchi and B. A. Malomed, Physical Review A 88, 043638 (2013).
  * Astrakharchik and Malomed [2015] G. Astrakharchik and B. A. Malomed, Physical Review A 92, 043632 (2015).
  * Malomed [2018] B. Malomed, Condensed Matter 3, 15 (2018).
  * Shamriz  _et al._ [2020a] E. Shamriz, Z. Chen, B. A. Malomed, and H. Sakaguchi, Condensed Matter 5, 20 (2020a).
  * Shamriz  _et al._ [2020b] E. Shamriz, Z. Chen, and B. A. Malomed, Physical Review A 101, 063628 (2020b).
  * Landau and Lifshitz [1976] L. D. Landau and E. M. Lifshitz,  _Mechanics, §14_ (Pergamon Press, Oxford New York, 1976).
  * Dzyaloshinskii  _et al._ [2019] I. E. Dzyaloshinskii, V. E. Zakharov, N. A. Inogamov, I. V. Kolokolov, E. A. Kuznetsov, V. V. Lebedev, S. P. Novikov, R. Z. Sagdeev, A. A. Starobinskii, S. M. Stishov, V. E. Fortov, and I. M. Khalatnikov, Physics-Uspekhi 62, 1282 (2019).
  * Andreev  _et al._ [2021] A. F. Andreev, S. A. Brazovskii, G. E. Volovik, E. I. Kats, P. S. Kondratenko, N. M. Kreines, V. P. Mineev, S. P. Novikov, S. P. Obukhov, L. P. Pitaevskii, D. E. Khmel’nitskii, and G. M. Eliashberg, Physics-Uspekhi 64, 214 (2021).
  * Buchholz [1969] H. Buchholz,  _The Confluent Hypergeometric Function_ (Springer Berlin Heidelberg, 1969).
  * Note [1] The discussed quantum state is not stationary; its energy is not an eigenvalue of the Hamiltonian. Therefore for this state only average energy is meaningful.


remarkRemark hypothesisHypothesis claimClaim MuyGPsA. Muyskens, B. Priest, I.
Goumiri, and M. Schneider ex_supplement

# MuyGPs: Scalable Gaussian Process Hyperparameter Estimation Using Local
Cross-Validation††thanks: Submitted to the editors April 29, 2021. This work
was performed under the auspices of the U.S. Department of Energy by Lawrence
Livermore National Laboratory under Contract DE-AC52-07NA27344 with IM release
number LLNL-JRNL-822013. Funding for this work was provided by LLNL Laboratory
Directed Research and Development grant 19-SI-004.

Amanda Muyskens Engineering Division, Lawrence Livermore National Laboratory
Livermore, CA 94550, USA<EMAIL_ADDRESS>Benjamin Priest Center for
Applied Scientific Computing, Lawrence Livermore National Laboratory
Livermore, CA 94550, USA Imène Goumiri Physics Division, Lawrence Livermore
National Laboratory
Livermore, CA 94550, USA Michael Schneider Physics Division, Lawrence
Livermore National Laboratory
Livermore, CA 94550, USA

###### Abstract

Gaussian processes (GPs) are non-linear probabilistic models popular in many
applications. However, naïve GP realizations require quadratic memory to store
the covariance matrix and cubic computation to perform inference or evaluate
the likelihood function. These bottlenecks have driven much investment in the
development of approximate GP alternatives that scale to the large data sizes
common in modern data-driven applications. We present in this manuscript
MuyGPs, a novel efficient GP hyperparameter estimation method. MuyGPs builds
upon prior methods that take advantage of the nearest neighbors structure of
the data, and uses leave-one-out cross-validation to optimize covariance
(kernel) hyperparameters without realizing a possibly expensive likelihood. We
describe our model and methods in detail, and compare our implementations
against the state-of-the-art competitors in a benchmark spatial statistics
problem. We show that our method outperforms all known competitors both in
terms of time-to-solution and the root mean squared error of the predictions.

###### keywords:

computational efficient, gaussian process, interpolation, land surface
temperature, machine learning, spatial statistics

##  1 Introduction

Gaussian process (GP) models are widely popular nonlinear models ubiquitous in
spatial statistics [8, 29] and commonly utilized in machine learning
applications [27, 22]. Gaussian process regression (GPR) tractably accounts
for the correlation among all observed data, making them favored models in the
interpolation of highly-nonlinear responses. GPR is popular in practice due to
its sample efficiency and closed-form posterior distribution inference.
Consequently, GPs are widely employed in applications where observations are
expensive to record or simulate and where overconfident prediction or
extrapolation has a high cost.

A process is said to follow a GP if any finite set of $n$ realizations of the
process follows a multivariate normal distribution. We conventionally assume
that the covariance has a parametric form depending on a pairwise _kernel
function_ $K_{\theta}(\cdot,\cdot)$ depending upon hyperparameters $\theta$.
Unfortunately, learning $\theta$ is typically difficult and computationally
expensive. The computing and storing a covariance has $O(n^{2})$ cost, while
realizations of GPR and evaluating the likelihood required in the conventional
training of $\theta$ have $O(n^{3})$ cost. Consequently, training a full GP
model is prohibitively expensive for large $n$ using common hardware.

The design of computationally efficient methods of Gaussian process estimation
is an active area of research. Many approximate methods for GP inference
attempt to sparsify the correlation matrix. Some approaches partition the
domain into spatially contiguous partition blocks and estimate independent
models within those blocks [14, 5, 28]. Alternately, others randomly sample
data partitions across the domain for independent Bayesian stationary GP
analyses and subsequently use the geometric median of the subset posteriors as
a means of model averaging [15]. Covariance tapering [6, 7] alters the
covariance matrix so that near-zero entries are considered independent and
therefore set to zero. Finally, the locally approximate GP (laGP) [13, 12]
uses a similar sparsification - although it does not realize an overall
covariance matrix - by fitting independent models within a local window of
each prediction location and assuming that all other observations.

Other methods assume instead that the precision matrix, the inverse of the
covariance matrix, is sparse. One approach is to assume that data observed on
a regular grid is spatially autogressive, which yields a sparse precision
matrix, and to employ multiple resolutions of radial basis functions for
efficient prediction [26]. Next, [21] rely on the equivalence between
stochastic partial differential equations (PDEs) and the Matérn covariance
fields. This method induces sparsity in the precision matrix by fitting
piecewise linear basis functions on a triangularization of the domain. Other
approaches still rewrite the multivariate normal likelihood in terms of
chained conditional distributions, inducing sparsity by modeling only a subset
of conditioning sets, such as those induced by the nearest neighbors structure
of the data [32, 2]. This is mathematically equivalent to a sparse precision
matrix defined through its Cholesky decomposition.

Spectral methods form effective models on gridded data. Such methods forgo
modeling the correlation function in favor of the spectral density, which is
the Fourier transformation of the covariance matrix. This allows estimating
$\theta$ via the Whittle Likelihood [34] or stochastic score approximation
[30]. Although the Whittle likelihood requires a full grid, [16] use an
iterative imputation scheme to allow for missing values. Further, [24]
demonstrate a fine grid can be used to approximately apply spectral methods to
irregularly-spaced observations.

Variational low-rank approximations to the covariance are particularly popular
in the machine learning literature [20, 31, 10]. The spatial statistics
literature refers to such models as predictive processes, where a small number
of knot locations across the domain induce a low-rank Nyström approximation to
the dense covariance matrix [1]. A related approach fits compactly supported
basis functions on recursively partitioned subregions of the domain, using
computations that can be parallelized [19, 18]. This basis function selection
uses a low-rank approach similar to that of the predictive process [1].
Another method that relies on multiple resolutions for estimation uses a
spatial process assumed to be the sum of a set of resolutions of Gaussian
basis functions, which they refer to as basic areal units (BAUs) [35, 36].

Other methods do not rely on Gaussian processes, or any statistical
distribution at all. For example, the Gapfill method relies on quantile
regression within a local window of the observation [9] . Since this method is
originally described for space-time data, but the data is often over space
only, the ”time” dimension is approximated by shifting the original image
[17].

Heaton et al. [17] directly compared many of these methods on a benchmark land
surface temperature dataset. This survey is particular notable because the
authors coordinated a blind competition, where each team implemented the
method in which they invented or otherwise owned expertise. Therefore, there
was no chance of a method underperforming due to misunderstanding or misuse of
the model or method. The central data problem on which the paper focuses is
the prediction of missing measurements from a large land surface temperature
gridded dataset sourced from the MODIS satellite. Instead of sampling the
testing set randomly across the domain, the authors use realistic cloud
coverage from an image taken on another day to form a realistic and
challengingly irregularly-shaped region of missing data. Ultimately, the
authors conclude that there is no one best method for computationally
efficient GP estimation as some methods produce estimates quickly while others
are more accurate in terms of root mean squared error (RMSE). Additionally,
[4] demonstrate their new methods on the same dataset. These methods are
improved in terms of time of estimation from most of the original methods, but
are not the top competitors in terms of RMSE.

Figure 1: Full (left) and training (right) land surface temperature datasets
used in [17] as a benchmark competition and sourced from the MODIS satellite
on August 4, 2016.

We have developed a GP estimation method designed to remain both accurate and
fast as data size grows. We invoke sparsity through locality neighborhoods
inspired by methods like [13]. However, instead of assuming independent models
within those local neighborhoods, we assume an overall stationary GP model
across the domain. Further, we entirely avoid maximum likelihood evaluation in
favor of leave-one-out cross-validation, where we produce predictions using
only the local neighborhoods. This amounts to the assumption that the weights
that are applied to the data vector in order to obtain the predictions
(kriging weights) are sparse, rather than the covariance matrix or precision
matrix as in previous approximation methods. We formally optimize the
hyperparameter values against the leave-one-out cross-validation objective
function. We limit ourselves in this work to mean squared error loss for our
objective functions, but other loss functions are possible. We then use these
estimated hyperparameters to realize responses at the prediction locations,
which are similarly formulated with only their nearest neighbors.

Leave-one-out cross-validation is typically a computationally expensive
process when predictions are obtained using the entire dataset. In fact, it
can be more expensive than the original maximum likelihood estimation problem
($O(n^{4})$), since a large matrix ($(n-1)\times(n-1)$)must be inverted for
each prediction. However, by limiting predictions to be based on only their
$k$ nearest neighbors, computation of each prediction reduces to many small,
parallelizable solves instead of one prohibitively large solve. We can further
improve our time-to-solution by relying upon a batch of $b$ training examples
for the leave-one-out training procedure, reducing the overall time complexity
to $O(bk^{3})$. In addition to making the time complexity formally independent
of $n$, this dramatically improves cross validation speed when $b\ll n$ and
$k\ll n$. Additionally, obtaining the $k$ nearest neighbors is more efficient
to compute than likelihood evaluation, and can be further accelerated by
approximate KNN algorithms. We find in our experiments that this deceptively
simple method is effective, fast, and leaves open the opportunity for further
future acceleration by way of distributed computation.

In this manuscript, we develop a novel method of approximate Gaussian
processes hyperparameter estimation that is efficient, scalable and accurate.
In section 2 we review a general stationary GP model. In section 3 we describe
our novel hyperparameter estimation method. In section 4 we demonstrate the
performance of our method in the benchmark dataset from [17], and show it
performs favorably in accuracy and computation time in comparison to all of
the existing state-of-the-art competitors. Finally in 5, we summarize our work
and outline future advancements that could be made in order to improve on our
methods.

##  2 Background: Gaussian Processes

We will consider throughout a univariate response
$Y:\mathcal{X}\rightarrow\mathbb{R}$, where
$\mathcal{X}\subseteq\mathbb{R}^{p}$ is the observation space. For notational
convenience and by convention we assume that $Y$ is de-trended and therefore
has zero mean. Extensions to non-zero and multivariate processes are trivial,
so we will avoid them for the sake of clarity. We say that $Y$ follows a
Gaussian process if the response at any finite set of $n$ points
$X=(\mathbf{x}_{1},\dots,\mathbf{x}_{n})\in\mathcal{X}^{n}$ follows a
multivariate normal distribution. That is,

(1)
$Y(X)=(Y(\mathbf{x}_{1}),\dots,Y(\mathbf{x}_{n}))^{T}\sim\mathcal{N}\left(\widetilde{0},K_{\theta}(X,X)\right),$

where $\mathcal{N}$ is the multivariate Gaussian distribution, $\widetilde{0}$
is the $n$-dimensional zero vector, and $K_{\theta}(X,X)$ is an $n\times n$
positive definite, symmetric covariance matrix between the elements of $X$
that is controlled non-linearly through kernel function
$K_{\theta}(\cdot,\cdot)$ and hyperparameters $\theta$.

Similarly, any finite set of $n^{*}$ unobserved data
$X^{*}=(\mathbf{x}^{*}_{1},\dots,\mathbf{x}^{*}_{n^{*}})\in\mathcal{X}^{n^{*}}$
is also jointly normal with observed data $X$ by the GP assumption. Thus, the
conditional distribution for the response at the new locations $X^{*}$ given
responses observed at $X$ is also multivariate normal with mean and variance

(2) $\displaystyle\widehat{Y}(X^{*}\mid X)$
$\displaystyle=K_{\theta}(X^{*},X)K_{\theta}(X,X)^{-1}Y(X),\text{ and}$ (3)
$\displaystyle\text{Var}(\widehat{Y}(X^{*}\mid X))$
$\displaystyle=K_{\theta}(X^{*},X^{*})-K_{\theta}(X^{*},X)K_{\theta}(X,X)^{-1}K_{\theta}(X,X^{*}),$

where $K_{\theta}(X^{*},X)=K_{\theta}(X,X^{*})^{T}$ is the cross-covariance
matrix between the elements of $X^{*}$ and $X$, and $K_{\theta}(X^{*},X^{*})$
is the covariance matrix between the elements of $X^{*}$, similar to
$K_{\theta}(X,X)$.

Note that the conditional mean in Equation 2 is the best linear unbiased
predictor (BLUP) for $Y(X^{*}|X)$, the conditional distribution of the
response at $X^{*}$ given data at $X$, even when the normality assumption of
the Gaussian process is violated. The $n$-dimensional row vectors of
$K_{\theta}(X^{*},X)K_{\theta}(X,X)^{-1}$ are referred to as the _kriging
weights_ , and these are the vectors our method assumes are sparse for
computational efficiency. Taking the posterior mean prediction given in
Equation 2 consists of computing the inner product of the kriging weights with
the de-trended data vector $Y(X)$. Most covariances exhibit no general closed
form solution to directly compute these weights without forming and inverting
the matrix $K_{\theta}(X,X)$, which can be prohibitively expensive in large
training data.

By the properties of the GP assumption, the log-likelihood of the observations
$Y(X)$ follow a jointly multivariate normal distribution:

(4)
$log(L(\theta,Y(X)))=-\frac{p}{2}log(2\pi)-\frac{1}{2}log(|K_{\theta}(X,X)|)-\frac{1}{2}Y(X)^{T}K_{\theta}(X,X)^{-1}Y(X).$

Often $K_{\theta}(X,X)$ is assumed to be stationary and isotropic. That is,
the covariance function is assumed to be a function of only the distances
between the elements of $X$. Stationarity implies that for locations
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}$,

(5)
$K_{\theta}(\mathbf{x}_{i},\mathbf{x}_{j})=\phi_{\theta}(||\mathbf{x}_{i}-\mathbf{x}_{j}||_{2}),$

where $\phi_{\theta}$ is a functional form with parameters $\theta$. One of
the most common covariance forms $\phi_{\theta}$ is the Matérn covariance
function. With covariance hyperparameters
$\theta=\\{\sigma^{2},\rho,\nu,\tau^{2}\\}$, and where $d$ is the distance
between two locations in $X$,

(6)
$\phi_{\theta}(d)=\sigma^{2}\left[\frac{2^{1-\nu}}{\Gamma(\nu)}\Bigg{(}\sqrt{2\nu}\frac{d}{\rho}\Bigg{)}^{\nu}K_{\nu}\Bigg{(}\sqrt{2\nu}\frac{d}{\rho}\Bigg{)}+\tau^{2}\mathbb{I}\\{d=0\\}\right],$

where $K_{\nu}$ is the modified Bessel function of the second kind. As
$\nu\to\infty$, this form converges pointwise to the well-known Gaussian (RBF)
covariance function.

Conventional approaches to training GP models estimate the covariance
parameters using the log-likelihood in Equation 4 via maximum likelihood
estimation, Bayesian analysis using Markov Chain Monte Carlo (MCMC), or by way
of grid search cross validation. However, these estimation methods are too
expensive to compute in large data since they require at least $O(n^{3})$
computation and $O(n^{2})$ memory. Investigators have developed scalable
approximate methods and models like those aforementioned in the introduction
in order to perform such estimation in large datasets. However, trends thus
far in the literature indicate a trade-off between speed and accuracy: those
methods with the fastest-time-to-solution on benchmark data yield posterior
mean predictions that are not competitive in terms of RMSE with more expensive
but accurate approximate methods [17]. We aim to obtain the best of both
worlds: accurate stationary Gaussian process predictions that maintain best-
in-class efficiency and scalability, while maintaining the accurate
uncertainty quantification.

##  3 MuyGPs

We describe a novel approximate method for training stationary Gaussian
process hyperparameters using leave-one-out cross-validation restricted to
local predictions. Our methodology derives from the union of two insights:
optimization by way of leave-one-out cross-validation allows us to avoid
evaluating the likelihood in Equation 4, and restriction to the $k$ nearest
neighbors of a prediction location limits the cost of computing the kriging
weights $K_{\theta}(X^{*},X)K_{\theta}(X,X)^{-1}$ in Equation 2 to $O(k^{3})$.
While other investigators may have exploited both of these observations in
different ways, our MuyGPs estimation method is the first to take advantage of
both simultaneously to accelerate kernel hyperparameter estimation by
enforcing sparsity in the kriging weights.

Leave-one-out cross-validation seeks hyperparameter values that minimize the
sum of an out-of-sample loss associated with computing the posterior
distribution of each training observation conditioned on all of the others. To
formally describe the procedure, let $\theta$ denote the hyperparameters that
require estimation and $K_{\theta}(\cdot,\cdot)$ a GP kernel function of
interest. In the Matérn kernel, define
$\theta=(\sigma^{2},\nu,\ell,\tau^{2})^{T}$, but this method is applicable to
other kernel forms. Here, we fix $\sigma^{2}=1$ in estimation of the other
parameters here since posterior mean predictions do not depend on overall
scale parameter $\sigma^{2}$. Hence, $\sigma^{2}$ is not estimable via any
cross-validation method for all kernel forms. We will introduce a different
efficient estimation protocol for $\sigma^{2}$ at the end of this section
after its fellow parameters have been estimated.

We must formulate the leave-one-out prediction of $Y(\mathbf{x}_{i})$ given
the set of all training points excluding $\mathbf{x}_{i}$. Define
$X_{-i}=(\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{i-1},\mathbf{x}_{i+1},\dots,\mathbf{x}_{n})$
and
$Y(X_{-i})=(Y(\mathbf{x}_{1}),Y(\mathbf{x}_{2}),\dots,Y(\mathbf{x}_{i-1}),Y(\mathbf{x}_{i+1}),\dots,Y(\mathbf{x}_{n}))$
to be the training locations and observed responses excluding the $i$th
training observation. Then we modify Equation 2 to obtain the mean GP
prediction by regressing the training labels on the corresponding inputs,

(7) $\widehat{Y}(\mathbf{x}_{i}\mid
X_{-i})=K_{\theta}(\mathbf{x}_{i},X_{-i})K_{\theta}(X_{-i},X_{-i})^{-1}Y(X_{-i}).$

Here $K_{\theta}(\mathbf{x}_{i},X_{-i})$ is the cross-covariance between
$\mathbf{x}_{i}$ and $X_{-i}$, while $K_{\theta}(X_{-i},X_{-i})$ is the
covariance among points in $X_{-i}$, both in terms of
$K_{\theta}(\cdot,\cdot)$. Although many criterion for predictions accuracy
are possible, we select the mean squared error criterion

(8)
$Q(\theta)=\frac{1}{n}\sum_{i=1}^{n}\left(Y(\mathbf{x}_{i})-\widehat{Y}(\mathbf{x}_{i}\mid
X_{-i})\right)^{2}.$

The hyperparameters are estimated to minimize the leave-one-out cross-
validation loss

(9) $\widehat{\theta}=\min_{\theta}Q(\theta).$

In our experiments we utilize the L-BFGS-B algorithm [37] to minimize Equation
9.

When there are a large number of observations, realizing an instance of the
cross-validation loss in Equation 8 so described is even more computationally
expensive than the loglikelihood in Equation 4; the procedures have complexity
$O(n^{4})$ and $O(n^{3})$, respectively. This suggests that the optimization
given by Equation 9 is even less efficient than traditional maximum likelihood
estimation. In order to achieve our scalability objectives, we replace the
full kriging of Equation 7 with _local kriging_ , using the $k$ nearest
neighbor locations of $\mathbf{x}_{i}$ instead of all of the $n-1$ locations
denoted by $X_{-i}$. This modifies the complexity of leave-one-out cross-
validation (computing $Q(\theta)$) to $O(nk^{3})$, which is much more scalable
than likelihood-based approaches when $k\ll n$.

We now illustrate our method precisely. Let $N_{i}$ be the set of $k$ indices
in $\\{1,\dots,i-1,i+1,\dots,n\\}$ indicating those elements of $X_{-i}$ that
are nearest to $\mathbf{x}_{i}$ in terms of distance. Similarly, define
$X_{N_{i}}$ as the set of training observations nearest to $\mathbf{x}_{i}$,
and let $Y(X_{N_{i}})$ be their corresponding responses. This allows us to
modify our leave-one-out prediction in Equation 7 to

(10) $\widehat{Y}_{NN}(\mathbf{x}_{i}\mid
X_{N_{i}})=K_{\theta}(\mathbf{x}_{i},X_{N_{i}})K_{\theta}(X_{N_{i}},X_{N_{i}})^{-1}Y(X_{N_{i}}),$

where the kriging weights and observed responses are defined in terms of
$N_{i}$.

Modifying Equation 8 in terms of the nearest neighbor leave-one-out
predictions from Equation 10, the resulting objective function requires
$O(nk^{3})$ operations to evaluate. We can further reduce our training
complexity overhead of $Q(\theta)$ by utilizing batching, a common technique
in machine learning. Let $B$ be a subset of $b$ indices randomly sampled
without replacement from $\\{1,\dots,n\\}$. We can then modify Equation 8 by
summing only over the nearest neighbor leave-one-out squared error of
locations $\mathbf{x}_{i}$ such that $i\in B$, whereby we obtain the modified
loss function

(11) $Q_{B}(\theta)=\frac{1}{b}\sum_{i\in
B}\left(Y(\mathbf{x}_{i})-\widehat{Y}_{NN}(\mathbf{x}_{i}\mid
X_{N_{i}})\right)^{2}.$

This approximation introduces some additional variability into the
optimization problem since the batch indices are randomly selected, but the
computational savings can be significant if $b\ll n$, bringing the cost of
evaluating the loss function to $O(bk^{3})$, which importantly is independent
of $n$. It is also important to note that the nearest neighbor index sets
$N_{i}$ utilized in the evaluations of $\widehat{Y}_{NN}(\mathbf{x}_{i})$ for
$i\in B$ in Equation 11 still range over the full $n$ data points $X$, rather
than only those data indicated by $B$. If nearest neighbor candidates were
restricted to batched points only, then the nearest neighbors of each batched
observation would be artificially distant from one another. Hence, more data
is ultimately used to estimate the hyperparameters than is included directly
in the batch.

Although $\sigma^{2}$ is not needed for the mean predictions, the prediction
errors depend linearly on its estimation since for batched observations it is
defined as

(12) $\text{Var}(\widehat{Y}_{NN}(\mathbf{x}_{i}\mid
X_{N_{i}}))=K_{\theta}(\mathbf{x}_{i},\mathbf{x}_{i})-K_{\theta}(\mathbf{x}_{i},X_{N_{i}})K_{\theta}(X_{N_{i}},X_{N_{i}})^{-1}K_{\theta}(X_{N_{i}},\mathbf{x}_{i}).$

After all other parameters are estimated via minimization of (11), we obtain
$\widehat{\sigma^{2}}$. Define
$\Omega_{\theta}=\frac{K_{\theta}}{\sigma^{2}}$. In general, there is a closed
fom equation for the maximum likelihood estimate for $\sigma^{2}$ given
$\theta$. However, this closed form solution involves forming and inverting
the full training covariance matrix $K_{\theta}(X,X)$, which is not possible
for large $n$. Our approximate estimate is the mean scale parameter within
each batched $k$ neighborhood or precisely

(13) $\widehat{\sigma^{2}}=\frac{1}{kb}\sum_{i\in
B}Y(X_{N_{i}})^{T}\Omega_{\theta}(X_{N_{i}},X_{N_{i}})^{-1}Y(X_{N_{i}}).$

Finally, predictions at testing (prediction) locations are computed using a
similarly local design. Define $X_{N_{i}^{\star}}$ as the set of training
observations nearest to the $i$th testing location $\mathbf{x}_{i}^{\star}$,
and let $Y(X_{N_{i}^{\star}})$ be their corresponding responses. Then,
approximate predictions and prediction errors are obtained via MuyGPs as

(14) $\displaystyle\widehat{Y}(X^{*}\mid X)$ $\displaystyle\approx
K_{\hat{\theta}}(X^{*},X_{N_{i}^{\star}})K_{\hat{\theta}}(X_{N_{i}^{\star}},X_{N_{i}^{\star}})^{-1}Y(X_{N_{i}^{\star}}),\text{
and}$ (15) $\displaystyle\text{Var}(\widehat{Y}(X^{*}\mid X))$
$\displaystyle\approx
K_{\hat{\theta}}(X^{*},X^{*})-K_{\hat{\theta}}(X^{*},X_{N_{i}^{\star}})K_{\hat{\theta}}(X_{N_{i}^{\star}},X_{N_{i}^{\star}})^{-1}K_{\hat{\theta}}(X_{N_{i}^{\star}},X^{*}),$

where $\hat{\theta}$ are the covariance parameters trained as aforementioned.

##  4 Numerical Studies

We demonstrate the effectiveness of our novel GP hyperparameter estimation on
a benchmark temperature dataset. The dataset considered is sourced from [17],
is pictured in Figure 1, and is available for download at
https://github.com/finnlindgren/heatoncomparison. The data were collected
using the Terra instrument onboard the MODIS satellite. It is composed of
daytime land surface temperatures for August 4, 2016 in longitudes -95.91153
to -91.28381 and latitudes 34.29519 to 37.06811 observed on a $500\times 300$
grid. This region is attractive due to its completeness, as 98.9% of locations
were observed in the area of interest on August 4, 2016. The cloud coverage
pattern from August 6, 2016 was applied to the data in order to create a
realistic missingness pattern for the training/testing data. This yielded
105,569 training observations, 42,740 testing observations, and 1,691 missing
(non-testing) observations, which are excluded from analysis. We normalize the
gridded observations to be in $[0,1]^{2}$ while maintaining the relationship
among latitude and longitude by

$\mathbf{x}_{i}=\frac{(latitude_{i},longitutde_{i})+218}{464},$

for $i=1,2,...105,569$ for all observations given in the training dataset.
This transformation creates a range of (0.2102833, 0.8078952) in latitude and
(0.001297902, 0.998651208) in longitude in the training observations. We apply
this same normalizing transformation to the testing locations.

We compute the same statistics to evaluate our methods used in the competition
in [17]. We consider mean absolute error (MAE), root mean squared error
(RMSE), continuous rank probability score (CRPS) [11], interval score (INT)
[11], and empirical coverage of 95% confidence intervals (COV). We computed
CRPS and INT source code drawn directly from [17] for consistency. The timings
of our method were obtained with a 2016 MacBook Pro with 2.9 CHz Quad-Core
Intel Core i7 with 16 GB of RAM. In [17], a machine with 256 GB of RAM and 2
Intel Xeon E5-2680 v4 2.40GHz CPUs with 14 coreseach and 2 threads per core -
totaling 56 possible threads for use in parallel computing was utilized for
comparison.

We implemented our newly-developed MuyGPs methodology as described in Section
3. Prior to fitting we must make modeling decisions including the selection of
a mean function and covariance function, the choice of batch size, and the
number of nearest neighbors to be employed. We realized models using a variety
of different choices in order to demonstrate the robustness of our results.
First, we assumed that the response function $Y$ has mean zero throughout the
entirety of this manuscript. Although this is a common assumption when fitting
stationary GPs, it is conventional and necessary in real-world data to first
remove mean trends prior to GP fitting. Therefore, define $Y^{obs}$ to be the
vector of non-zero mean temperature responses, and
$Y_{k}(x)=Y^{obs}(X)-\mu_{k}(X)$, where $Y_{k}(X)$ is the mean zero response
with mean function $\mu_{k}(X)$ removed (previously notated $Y(X)$). We
consider three simple mean functions so that $k=1,2,3$.

First, we consider a constant mean function. Let

$\mu_{1}(X)=c.$

This simple mean is estimated by the sample mean of the training data
$\hat{c}=\frac{1}{n}\sum_{i=1}^{n}Y^{obs}(\mathbf{x}_{i})$. Next, we consider
a linear mean function with an interaction term. Define $Z$ to be the $n\times
4$ matrix that contains a $n\times 1$ vector of ones, $X$, and the
multiplication of the rows elements of $X$ so that
$Z=[1_{n},X^{T},X^{T}[,1]*X^{T}[,2]]$. Then define

$\mu_{2}(X)=Z\beta,$

where $\beta$ is a $4\times 1$ vector of mean parameters. These parameters are
estimated using the ordinary least squares solution of the training data.

(a) $\mu_{1}$

(b) $Y_{1}$

(c) $\mu_{2}$

(d) $Y_{2}$

(e) $\mu_{3}$

(f) $Y_{3}$

Figure 2: Mean functions considered, and the residual observations that are
fit via MuyGPs using the testing data.

Finally, we consider a mean that is based on kernel smoothing. In particular,
we utilize the Nadaraya-Watson kernel smoother [25, 33]. Define

$\mu_{3}(X)=\frac{\sum_{i=1}^{n}G(X,\mathbf{x}_{i})Y(\mathbf{x}_{i})}{\sum_{i=1}^{n}G(X,\mathbf{x}_{i})},,$

where $G(X,\mathbf{x}_{i})$ is the Gaussian kernel observed for locations $X$
and $\mathbf{x}_{i}$
($\exp{-||X-\mathbf{x}_{i}||/\rho}$). We use the $smooth.2d$ function from the
package $fields$ in order to estimate implement this smoothing with length
scale $\rho=25$ on the $500\times 300$ grid [3]. This implementation utilizes
a 2-D fast Fourier transform (FFT) to efficiently estimate the smoother using
the partially gridded training data [3]. Employing this mean is similar to the
multi-resolution principle in [26, 19], where very smooth low frequency trends
are fit by this filtering mean, and the high-frequency trend are fit on
resulting residuals. We plot the three mean functions ($\mu_{k}$ ), and the
corresponding data fit in each mean ($Y_{k}$) is seen in Figure 2. Note that
because $\mu_{2}(X)$ and $\mu_{3}(X)$ depend on $X$, they lead to first order
non-stationary predictions.

MuyGPs can theoretically estimate the parameters of most kernel functions. For
its flexibility, we employ the Matérn covariance as stated in Equation 6 in
our experiments. Since the impact on prediction is low, we fix
$\tau^{2}=0.001$, and if otherwise not stated, we fix the length scale
$\rho=1.0$. In our experiments, $\nu$ and $\rho$ are difficult to
simultaneously estimate, and we found better performance when $\nu$ is allowed
to vary. Further, we explore other fixed values of $\rho$ in the results in
Table 1.

Figure 3: Empirical 90% confidence intervals of RMSE over 100 simulation
iterations for 50 exact nearest neighbors with $\mu_{3}$. Figure 4:
Performance in RMSE and computing time for 500 batched values using
approximate nearest neighbors and $\mu_{3}$.

Next, we must select a batch size for the MuyGPs estimation procedure. In all
of our experiments we sampled batch elements from our training data randomly
and without replacement. Using exact nearest neighbor sets of size 50 and
$\mu_{3}$, we demonstrate the variability in prediction performance using RMSE
for various batch sizes in Figure 3 for 100 independent simulation iterations
at each batch size. We then compute empirical 95% empirical confidence
intervals at each batch size to determine the variability of estimates from
that batch size. In even extremely small batch sizes of 25 members, the
variability in the RMSE is very small at a difference of about 0.014 length in
a 95% empirical confidence interval of the simulations. When considering batch
sizes of 500, there is less than 0.001 in variability, and at batch sizes of
2,000, there is nearly no variability in prediction RMSE. Note that this
dataset has 105,569 training observations, so a batch size of 2,000 offers
significant improvement where we can use less that 1.8% of the data in the
leave-one-out cross-validation optimization without impacting the results.

Figure 5: Comparison in RMSE and computing time in comparison of MuyGPs to all
methods in [17] and [4]. The multiple MuyGPs observations are the results from
Figure 4 with varying nearest neighbor set sizes between 25 and 200 in 25
incriments.

Following the tradeoffs associated with batch size selection, we demonstrate
selecting the number of nearest neighbors. As more neighbors are incorporated
into the model, the predictions become more accurate, but the computing time
is increased. We can achieve practical speed improvements by employing
approximate nearest neighbor algorithms such as the hierarchical navigable
small world (HNSW) algorithm [23]. HNSW offers further scaling benefits when
the number of variables in $\mathbf{x}_{i}$ is large. We implement this via
bindings for the ’hnswlib’ library available at
https://github.com/nmslib/hnswlib. We demonstrate this performance tradeoff in
including different numbers of approximate nearest neighbors in Figure 4 using
$\mu_{3}$ as an example. There is a large accuracy performance gap between 25
and 50 nearest neighbors, but all increases in number of nearest neighbors
improves the prediction accuracy, but increases total computing time due to
the $O(k^{3})$ term associated with solving the nearest neighbor systems of
linear equations.

Table 1: Numerical scoring for various versions of MuyGPs with competing method original results from [17] and [4] for comparison. Method | $\ell$ | MAE | RMSE | CRPS | INT | COV | Time (min)
---|---|---|---|---|---|---|---
| 0.1 | 1.14 | 1.66 | 0.86 | 9.25 | 0.95 | 0.70
| 0.25 | 1.15 | 1.64 | 0.84 | 8.40 | 0.95 | 0.68
MuyGPs, $\mu_{1}$ | 0.5 | 1.19 | 1.67 | 0.85 | 8.02 | 0.93 | 0.75
| 0.75 | 1.21 | 1.69 | 0.86 | 7.90 | 0.93 | 0.75
| 1.0 | 1.22 | 1.68 | 0.86 | 7.85 | 0.92 | 0.63
| 0.1 | 1.12 | 1.64 | 0.86 | 9.38 | 0.95 | 0.64
| 0.25 | 1.13 | 1.62 | 0.83 | 8.31 | 0.94 | 0.64
MuyGPs, $\mu_{2}$ | 0.5 | 1.15 | 1.62 | 0.83 | 8.00 | 0.94 | 0.64
| 0.75 | 1.19 | 1.65 | 0.85 | 7.85 | 0.93 | 0.63
| 1.0 | 1.19 | 1.64 | 0.84 | 7.80 | 0.93 | 0.63
| 0.1 | 1.07 | 1.54 | 0.84 | 9.35 | 0.95 | 0.76
| 0.25 | 1.08 | 1.53 | 0.80 | 8.24 | 0.94 | 0.66
MuyGPs, $\mu_{3}$ | 0.5 | 1.12 | 1.55 | 0.81 | 7.95 | 0.94 | 0.62
| 0.75 | 1.14 | 1.56 | 0.81 | 7.78 | 0.93 | 0.63
| 1.0 | 1.15 | 1.57 | 0.82 | 7.71 | 0.93 | 0.65
FRK |  | 1.96 | 2.44 | 1.44 | 14.08 | 0.79 | 2.32
Gapfill |  | 1.33 | 1.86 | 1.17 | 34.78 | 0.36 | 1.39
LatticeKrig |  | 1.22 | 1.68 | 0.87 | 7.55 | 0.96 | 27.92
LAGP |  | 1.65 | 2.08 | 1.17 | 10.81 | 0.83 | 2.27
Metakriging |  | 2.08 | 2.50 | 1.44 | 10.77 | 0.89 | 2888.52
MRA |  | 1.33 | 1.85 | 0.94 | 8.00 | 0.92 | 15.61
NNGP |  | 1.21 | 1.64 | 0.85 | 7.57 | 0.95 | 2.06
NNGP2 |  | 1.24 | 1.68 | 0.87 | 7.50 | 0.94 | 42.85
Partition |  | 1.41 | 1.80 | 1.02 | 10.49 | 0.86 | 79.98
Pred. Proc. |  | 2.15 | 2.64 | 1.55 | 15.51 | 0.83 | 160.24
SPDE |  | 1.10 | 1.53 | 0.83 | 8.85 | 0.97 | 120.33
Tapering |  | 1.87 | 2.45 | 1.32 | 10.31 | 0.93 | 133.26
Per Embed |  | 1.29 | 1.79 | 0.91 | 7.44 | 0.93 | 9.81
PALM |  | 1.59 | 1.93 | 1.15 | 11.78 | 0.78 | 4.64
Global + PALM2 |  | 1.44 | 1.76 | 1.03 | 9.28 | 0.84 | 4.64

Finally, we compare the performance of our method to all previous methods
computed on this data in [17] and [4]. We plot the tradeoffs between RMSE and
time visually in Figure 5, where we consider multiple numbers of nearest
neighbors in 4, with $\mu_{3}$ , 500 batched values, and $\rho=1$. Notably,
MuyGPs with a small number of nearest neighbors is the fastest method
presented and remains among the top five methods in terms of RMSE.
Additionally, MuyGPs with a large number of nearest neighbors is the most
accurate method considered in all manuscripts. Further, MuyGPs with number of
nearest neighbors greater than or equal to 50 ranks 2nd or better in all the
methods in accuracy. In [17], SPDE [21] is the most accurate method, but it
took approximately 120 minutes to compute, which is roughly more than 12 times
slower than the slowest presented MuyGPs setting. A similarly accurate
solution is obtained by MuyGPs in approximately only 5 minutes.

We consider the other statistics based on uncertainty in the results in Table
1. In all previous studies, the parameter $\ell$ from the Matérn covariance
has been fixed at 1.0. This is because the $\ell$ and $\nu$ parameters are
unidentifiable to estimate simultaneously. Given this were a blind competition
as the participants of [17], this value would likely have been fixed at a
reasonable value. To show the robustness of our method, we demonstrate the
performance of MuyGPs for several selections of fixed $\ell$ values. In fact,
our previous selection of $\ell=1.0$ is the least favorable reasonable value
for $\ell$ in the main results of Table 1. With all three of the mean
functions, even with only 50 approximate nearest neighbors, the accuracy of
MuyGPs is competitive. Further, all methods compared were under a minute in
computing time, which is faster than any method previously tested on this
data. In all cases, the coverage is near to that of the the expected 0.95, and
the other uncertainty interval statistics perform favorable in the field of
methods. In summary, the MuyGPs method offers more accurate and faster
computation time than any existing method, and no matter the modeling choices,
is a competitive with all known scalable GP computation methods.

##  5 Discussion

We have presented in this manuscript MuyGPs, a new method for stationary
Gaussian process hyperparameter estimation and prediction that is both
computationally efficient and highly performant when compared to state-of-the-
art approximate GP estimation methods. MuyGPs minimizes the leave-one-out
cross-validated mean squared error using only data from the $k$ nearest
neighbors for selected locations across the domain. Similarly, predictions
depend only upon their $k$ nearest neighbor training observations. We
demonstrate that although this method is simple, it is powerful and performs
well against state-of-the-art approximation methods.

Using this method requires the selection of a mean function, covariance
function, batch size, and number of nearest neighbors. We demonstrate that our
method performs well with very few batched points. In fact, when exact nearest
neighbors are employed, the accuracy has nearly no variance with a batch size
of about 2,000, which is approximately 1.8% the training data. Increasing the
number of nearest neighbors improves performance in terms of RMSE at the cost
of increasing computing time. Finally, we show that a non-stationary mean
function improves the accuracy of our model, but the magnitude of this effect
is small compared to changes in GP estimation methods. Although these
selections ultimately define a tradeoff between performance in accuracy and
computation time, our method performs favorably compared to existing methods
regardless of these modeling choices. Settings of our method offers both the
fastest, and the most accurate method among all known methods, and performs
among the top few methods in all other considered statistics.

To date we have optimized our method only with respect to MSE, but other
objective functions are possible. Using optimal parameters for this testing
dataset yields an RMSE of 1.42 with 0.95 coverage via our model with only 50
approximate nearest neighbors (although these are not the parameters selected
by our training method). This implies there may be room for improvement of
parameter estimation using a more complex objective function. One possible
idea could be to incorporate other statistics from Table 1 into the objective
function.

Although we achieve a first-order non-stationary predictions through
non-stationary mean functions, our hyperparameter estimation method is
amenable to estimation of second-order non-stationary models that have
previously been considered too expensive for large data. For example, one
could assume hierarchically that, for example, $\nu$ parameter follows a GP
over the domain. Although we found small batch sizes relatively effective with
a random sample under a stationarity assumption, more care may be needed in
order to select the batched observations and their size under this more
complex non-stationary model case.

## Acknowledgments

This work was performed under the auspices of the U.S. Department of Energy by
Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344 with
IM release number LLNL-JRNL-822013. Funding for this work was provided by LLNL
Laboratory Directed Research and Development grant 19-SI-004.

This document was prepared as an account of work sponsored by an agency of the
United States government. Neither the United States government nor Lawrence
Livermore National Security, LLC, nor any of their employees makes any
warranty, expressed or implied, or assumes any legal liability or
responsibility for the accuracy, completeness, or usefulness of any
information, apparatus, product, or process disclosed, or represents that its
use would not infringe privately owned rights. Reference herein to any
specific commercial product, process, or service by trade name, trademark,
manufacturer, or otherwise does not necessarily constitute or imply its
endorsement, recommendation, or favoring by the United States government or
Lawrence Livermore National Security, LLC. The views and opinions of authors
expressed herein do not necessarily state or reflect those of the United
States government or Lawrence Livermore National Security, LLC, and shall not
be used for advertising or product endorsement purposes.

## References

  * [1] S. Banerjee, A. E. Gelfand, A. O. Finley, and H. Sang, Gaussian predictive process models for large spatial data sets, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 70 (2008), pp. 825–848.
  * [2] A. Datta, S. Banerjee, A. O. Finley, and A. E. Gelfand, Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets, Journal of the American Statistical Association, 111 (2016), pp. 800–812.
  * [3] Douglas Nychka, Reinhard Furrer, John Paige, and Stephan Sain, fields: Tools for spatial data, 2017, https://doi.org/10.5065/D6W957CT, https://github.com/NCAR/Fields.  R package version 11.6.
  * [4] A. M. Edwards and R. B. Gramacy, Precision aggregated local models, 2020, https://arxiv.org/abs/2005.13375.
  * [5] M. Fuentes, Spectral methods for nonstationary spatial processes, Biometrika, 89 (2002), pp. 197–210.
  * [6] R. Furrer, M. G. Genton, and D. Nychka, Covariance tapering for interpolation of large spatial datasets, Journal of Computational and Graphical Statistics, 15 (2006), pp. 502–523.
  * [7] R. Furrer and S. R. Sain, spam: A sparse matrix r package with emphasis on mcmc methods for gaussian markov random fields, Journal of Statistical Software, 36 (2010), pp. 1–25.
  * [8] A. E. Gelfand, P. Diggle, P. Guttorp, and M. Fuentes, Handbook of spatial statistics, CRC press, 2010.
  * [9] F. Gerber, R. de Jong, M. E. Schaepman, G. Schaepman-Strub, and R. Furrer, Predicting missing values in spatio-temporal remote sensing data, IEEE Transactions on Geoscience and Remote Sensing, 56 (2018), pp. 2841–2853.
  * [10] M. N. Gibbs and D. J. MacKay, Variational gaussian process classifiers, IEEE Transactions on Neural Networks, 11 (2000), pp. 1458–1464.
  * [11] T. Gneiting and A. E. Raftery, Strictly proper scoring rules, prediction, and estimation, Journal of the American statistical Association, 102 (2007), pp. 359–378.
  * [12] R. B. Gramacy, lagp: Large-scale spatial modeling via local approximate gaussian processes in r, Journal of Statistical Software, 72 (2016), pp. 1–46.
  * [13] R. B. Gramacy and D. W. Apley, Local gaussian process approximation for large computer experiments, Journal of Computational and Graphical Statistics, 24 (2015), pp. 561–578.
  * [14] R. B. Gramacy et al., tgp: an r package for bayesian nonstationary, semiparametric nonlinear regression and design by treed gaussian process models, Journal of Statistical Software, 19 (2007), p. 6.
  * [15] R. Guhaniyogi and S. Banerjee, Meta-kriging: Scalable bayesian modeling and inference for massive spatial datasets, Technometrics, 60 (2018), pp. 430–444.
  * [16] J. Guinness, Spectral density estimation for random fields via periodic embeddings, Biometrika, 106 (2019), pp. 267–286.
  * [17] M. J. Heaton, A. Datta, A. O. Finley, R. Furrer, J. Guinness, R. Guhaniyogi, F. Gerber, R. B. Gramacy, D. Hammerling, M. Katzfuss, et al., A case study competition among methods for analyzing large spatial data, Journal of Agricultural, Biological and Environmental Statistics, 24 (2019), pp. 398–425.
  * [18] M. Jurek and M. Katzfuss, Multi-resolution filters for massive spatio-temporal data, 2019, https://arxiv.org/abs/1810.04200.
  * [19] M. Katzfuss, A multi-resolution approximation for massive spatial datasets, Journal of the American Statistical Association, 112 (2017), pp. 201–214.
  * [20] M. Lázaro-Gredilla and M. K. Titsias, Variational heteroscedastic gaussian process regression, in ICML, 2011.
  * [21] F. Lindgren, H. Rue, and J. Lindström, An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73 (2011), pp. 423–498.
  * [22] H. Liu, Y.-S. Ong, X. Shen, and J. Cai, When gaussian process meets big data: A review of scalable gps, IEEE transactions on neural networks and learning systems, 31 (2020), pp. 4405–4423.
  * [23] Y. A. Malkov and D. A. Yashunin, Efficient and robust approximate nearest neighbor search using hierarchical navigable small world graphs, IEEE transactions on pattern analysis and machine intelligence, 42 (2018), pp. 824–836.
  * [24] A. Muyskens, J. Guinness, and M. Fuentes, Non-stationary covariance estimation using the stochastic score approximation for large spatial data, arXiv preprint arXiv:1809.06904, (2018).
  * [25] E. A. Nadaraya, On estimating regression, Theory of Probability & Its Applications, 9 (1964), pp. 141–142.
  * [26] D. Nychka, S. Bandyopadhyay, D. Hammerling, F. Lindgren, and S. Sain, A multiresolution gaussian process model for the analysis of large spatial datasets, Journal of Computational and Graphical Statistics, 24 (2015), pp. 579–599.
  * [27] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning), The MIT Press, 2005\.
  * [28] H. Sang, M. Jun, and J. Z. Huang, Covariance approximation for large multivariate spatial data sets with an application to multiple climate model errors, The Annals of Applied Statistics, (2011), pp. 2519–2548.
  * [29] M. L. Stein, Interpolation of spatial data: some theory for kriging, Springer Science & Business Media, 2012.
  * [30] M. L. Stein, J. Chen, and M. Anitescu, Stochastic approximation of score functions for gaussian processes, The Annals of Applied Statistics, (2013), pp. 1162–1191.
  * [31] D. Tran, R. Ranganath, and D. M. Blei, The variational gaussian process, arXiv preprint arXiv:1511.06499, (2015).
  * [32] A. V. Vecchia, Estimation and model identification for continuous spatial processes, Journal of the Royal Statistical Society: Series B (Methodological), 50 (1988), pp. 297–312.
  * [33] G. S. Watson, Smooth regression analysis, Sankhyā: The Indian Journal of Statistics, Series A, (1964), pp. 359–372.
  * [34] P. Whittle, On stationary processes in the plane, Biometrika, (1954), pp. 434–449.
  * [35] A. Zammit-Mangion and N. Cressie, Frk: An r package for spatial and spatio-temporal prediction with large datasets, arXiv preprint arXiv:1705.08105, (2017).
  * [36] A. Zammit-Mangion and N. Cressie, Frk: An r package for spatial and spatio-temporal prediction with large datasets, 2018, https://arxiv.org/abs/1705.08105.
  * [37] C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization, ACM Transactions on Mathematical Software (TOMS), 23 (1997), pp. 550–560.