# Spectral Diffusion Processes

Angus Phillips , Thomas Seror11footnotemark: 1 , Michael
Hutchinson22footnotemark: 2 ,
Valentin De Bortoli33footnotemark: 3 , Arnaud Doucet22footnotemark: 2 , Émile
Mathieu Equal contribution. Correspondence to:
<EMAIL_ADDRESS>of Statistics, University of Oxford,
Oxford, UK.Dept. of Computer Science ENS, CNRS, PSL University Paris,
France.Dept. of Engineering, University of Cambridge, Cambridge, UK.

###### Abstract

Score-based generative modelling (SGM) has proven to be a very effective
method for modelling densities on finite-dimensional spaces. In this work we
propose to extend this methodology to learn generative models over functional
spaces. To do so, we represent functional data in spectral space to dissociate
the stochastic part of the processes from their space-time part. Using
dimensionality reduction techniques we then sample from their stochastic
component using finite dimensional SGM. We demonstrate our method’s
effectiveness for modelling various multimodal datasets.

##  1 Introduction

Score-based Generative Models (SGMs) (Song and Ermon,, 2019; Song et al.,,
2021; Ho et al.,, 2020; Dhariwal and Nichol,, 2021) are a powerful class of
generative models. Noise is progressively added to data following a Stochastic
Differential Equation (SDE)—the forward _noising_ process—until it is
approximately Gaussian. The generative model is given by an approximation of
the associated time-reversed process called the backward _denoising_ process.
This process is an SDE whose drift depends on the gradient of the logarithm of
the densities of the forward process. The drift is estimated by leveraging
techniques from deep learning and score matching (Hyvärinen,, 2005; Vincent,,
2011).

Probabilistic modelling over functional spaces is of high importance in the
physical sciences where objects of interest are often _fields_ , mapping
(space-time) inputs to tensors like scalars or vectors. Aleatoric and
epistemic uncertainties motivate for a _probabilistic_ treatment of these
phenomena, like weather forecasting (Ravuri et al.,, 2021) or electric
potential prediction (Yang et al.,, 2020).

A well-established framework to model distributions over functions are
Gaussian processes (GPs) (Rasmussen,, 2003). To alleviate their scalability
and flexibility limitations, neural processes (NPs) have been introduced and
successfully applied on a wide variety of probabilistic prediction tasks. Yet,
their Gaussian likelihood assumption limits their flexibility.

Figure 1: Illustration of our methodology. SGM is performed in a spectral
space.

In this work, we circumvent this problem by casting the modelling problem in a
spectral space. Doing so, we decouple the _stochastic_ part of the process
from its _space-time_ part. This approach turns an _uncountable_ dimensional
problem into a _countable_ dimensional one, which is turned into a _finite_
dimensional one by truncation. Finally, we deal with this last problem using
standard SGM techniques. Our methodology is illustrated in Fig. 1 and
presented in more detail in Sec. 3.

##  2 Background

#### Score-based generative modelling.

We briefly recall the concepts behind SGMs on the Euclidean space
$\mathbb{R}^{d}$. Let $(\mathbf{Y}_{t})_{t\geq 0}$ be a _noising_ process
defined by the following SDE

$\mathrm{d}\mathbf{Y}_{t}=-\mathbf{Y}_{t}\mathrm{d}t+\sqrt{2}\mathrm{d}\mathbf{B}_{t},\quad\mathbf{Y}_{0}\sim
p_{0},$ (1)

with $(\mathbf{B}_{t})_{t\geq 0}$ a $d$-dimensional Brownian motion and
$p_{0}$ the data distribution on $\mathbb{R}^{d}$. This process converges
towards the unit Gaussian distribution $\mathrm{N}(0,\operatorname{Id})$.
Under conditions on $p_{0}$, the time-reversal
$(\bar{\mathbf{Y}}_{t})_{t\in\left[0,T\right]}=(\mathbf{Y}_{T-t})_{t\in\left[0,T\right]}$
is given by (Cattiaux et al.,, 2021; Haussmann and Pardoux,, 1986)

$\mathrm{d}\bar{\mathbf{Y}}_{t}=\\{\bar{\mathbf{Y}}_{t}+2\nabla\log
p_{T-t}(\bar{\mathbf{Y}}_{t})\\}\mathrm{d}t+\sqrt{2}\mathrm{d}\mathbf{B}_{t},\quad\bar{\mathbf{Y}}_{0}\sim
p_{T},$ (2)

with $p_{t}$ the density of $\bm{Y}_{t}$. Denoising diffusion models are
defined to sample approximately from Eq. 2. First, we replace $p_{T}$ by
$\mathrm{N}(0,\operatorname{Id})$. Then $\nabla\log p_{t}$ is approximated by
a parametric function $\bm{s}_{\theta}(t,\cdot)$ so as to minimise
$\ell_{t}(\mathbf{s})=\mathbb{E}[\|\mathbf{s}_{\theta}(t,\mathbf{Y}_{t})-\nabla_{y_{t}}\log
p_{t|0}(\mathbf{Y}_{t}|\mathbf{Y}_{0})\|^{2}]$ (Hyvärinen,, 2005). Finally, an
Euler–Maruyama discretisation of Eq. 2 is considered using a discretisation
step $\gamma$ and $N=T/\gamma$

$\bar{Y}_{n+1}=\bar{Y}_{n}+\gamma\\{\bar{Y}_{n}+2\mathbf{s}_{\theta}(T-n\gamma,\bar{Y}_{n})\\}+\sqrt{2\gamma}G_{n+1},\quad\bar{Y}_{0}\sim\mathrm{N}(0,\operatorname{Id}),\quad
G_{n+1}\overset{\textup{i.i.d.}}{\sim}\mathrm{N}(0,\operatorname{Id}).$ (3)

#### Stochastic processes.

Here we are interested in modelling distributions over functions. Given a
compact input space $\mathcal{X}$, a $\mathbb{R}^{d}$-valued stochastic
process $(\mathbf{Y}_{x})_{x\in\mathcal{X}}$ is a collection of random
variables $\mathbf{Y}_{x}$ 111In this work, we focus on _continuous_
stochastic processes, i.e. any trajectory of
$(\mathbf{Y}_{x})_{x\in\mathcal{X}}$ belongs to
$\mathrm{C}(\mathcal{X},\mathbb{R}^{d})$. In that case
$\omega\mapsto(\mathbf{Y}_{x}(\omega))_{x\in\mathcal{X}}$ is measurable and we
denote $\pi\in\mathscr{P}(\mathrm{C}(\mathcal{X},\mathbb{R}^{d}))$ its
distribution.. Most existing methods which model stochastic processes rely on
the parametrisation of finite dimensional marginals
$\\{\mathbf{Y}_{x_{i}}\,:i\in\\{1,\dots,n\\},\ x_{i}\in\mathcal{X}\\}$ for
every $n\in\mathbb{N}$. However, it is not clear _a priori_ that a collection
of probability distributions
$\mathrm{S}=\\{\pi_{x_{1},\dots,x_{n}}\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}$ defines a distribution
$\pi\in\mathscr{P}((\mathbb{R}^{d})^{\mathcal{X}})$222$\mathscr{P}((\mathbb{R}^{d})^{\mathcal{X}})$
is the space of distributions over the product space
$(\mathbb{R}^{d})^{\mathcal{X}}$.. This is the focus of the Kolmogorov
extension theorem (Charalambos and Aliprantis,, 2013, Thm 15.26) which ensures
the existence of such a distribution under a set of conditions. The conditions
required of the system are exchangeability and consistency. A system
$\mathrm{S}$ is _exchangeable_ , if for any $n\in\mathbb{N}$, $n$-permutation
$\sigma\in\mathfrak{S}_{n}$, $x_{1},\dots,x_{n}\in\mathcal{X}$ and continuous
and bounded function $f\in\mathrm{C}_{b}((\mathbb{R}^{d})^{n})$, we have $\int
f(y_{\sigma^{-1}(1)},\dots,y_{\sigma^{-1}(n)})\mathrm{d}\pi_{x_{\sigma(1)},\dots,x_{\sigma(n)}}(y_{1},\dots,y_{n})=\int
f(y_{1},\dots,y_{n})\mathrm{d}\pi_{x_{1},\dots,x_{n}}(y_{1},\dots,y_{n})$. A
system $\mathrm{S}$ is _consistent_ , if for any $n_{1}\leq n_{2}$,
$x_{1},\dots,x_{n_{2}}\in\mathcal{X}$ and
$f\in\mathrm{C}_{b}((\mathbb{R}^{d})^{n_{1}})$ we have $\int
f(y_{1},\dots,y_{n_{1}})\mathrm{d}\pi_{x_{1},\dots,x_{n_{1}}}(y_{1},\dots,y_{n_{1}})=\int
f(y_{1},\dots,y_{n_{1}})\mathrm{d}\pi_{x_{1},\dots,x_{n_{2}}}(y_{1},\dots,y_{n_{2}})$.
Existing generative modelling approaches often satisfy the exchangeability
criterion but fall short on the consistency. In this work we aim at bridging
this gap by developing a _spectral_ method which immediately yields the
consistency of the obtained system.

##  3 Spectral Score-Based Processes

#### Karhunen-Loève theorem.

In order to derive a consistent method, our key insight is to dissociate the
_stochastic_ part of the process from its _space-time_ component using a
spectral decomposition. This approach stems from the Karhunen-Loève theorem
which we recall below.

######  Theorem 1.

Let $\mathcal{X}$ be a compact space and $(\mathbf{Y}_{x})_{x\in\mathcal{X}}$
a continuous stochastic process such that for any $x\in\mathcal{X}$,
$\mathbb{E}[\mathbf{Y}_{x}]=0$ and
$\mathbb{E}[\|\mathbf{Y}_{x}\|^{2}]<+\infty$. For any
$x_{1},x_{2}\in\mathcal{X}$, we denote the covariance kernel as
$\mathrm{K}_{\mathbf{Y}}(x_{1},x_{2})=\mathbb{E}[\langle\mathbf{Y}_{x_{1}},\mathbf{Y}_{x_{2}}\rangle]$
and define a $\mathrm{L}^{2}(\mathcal{X})$ operator
$T_{\mathrm{K}_{\mathbf{Y}}}$ given for any $f\in\mathrm{L}^{2}(\mathcal{X})$
by
$T_{\mathrm{K}_{\mathbf{Y}}}f(x_{2})=\int_{\mathcal{X}}\mathrm{K}_{\mathbf{Y}}(x_{1},x_{2})f(x_{1})\mathrm{d}x_{1}$.
Denoting by
$(e_{m})_{m\in\mathbb{N}}\in(\mathrm{L}^{2}(\mathcal{X}))^{\mathbb{N}}$,
$(\lambda_{m})_{m\in\mathbb{N}}\in\mathbb{R}^{\mathbb{N}}$ the eigenfunctions
and eigenvalues of $T_{\mathrm{K}_{\mathbf{Y}}}$, we have

$\textstyle{\lim_{M\to+\infty}\mathbb{E}[\sup_{x\in\mathcal{X}}\|\mathbf{Y}_{x}-\sum_{m=0}^{M}{\lambda_{m}}^{1/2}Z_{m}e_{m}(x)\|^{2}]=0,}$
(4)

where for any $m\in\mathbb{N}$,
$\textstyle{Z_{m}=\lambda_{m}^{-1/2}\int_{\mathcal{X}}\langle\mathbf{Y}_{x},e_{m}(x)\rangle\mathrm{d}x}$.

Therefore, a natural approximation of the process
$(\mathbf{Y}_{x})_{x\in\mathcal{X}}$ 333We assume that
$\mathbb{E}[\|\mathbf{Y}_{x}\|^{2}]<+\infty$ and w.l.o.g.
$\mathbb{E}[\mathbf{Y}_{x}]=0$ as the mean can be subtracted otherwise. is
given by $F_{M}\triangleq\sum_{m=0}^{M}{\lambda_{m}}^{1/2}Z_{m}e_{m}$. In the
previous theorem, we dissociated the _stochastic_ part of the process—embedded
in $\\{Z_{m}\\}_{m=0}^{M}$—from its _space-time_ part—embedded in
$\\{e_{m}\\}_{m=0}^{M}$. In fact, even though our approach is originally based
on the Karhunen-Loève theorem, in practice we consider _arbitrary_ kernels
$\mathrm{K}$ (not necessarily based on the process $\mathbf{Y}$) using
Mercer’s theorem (Ferreira and Menegatto,, 2009, Thm 1.1).

######  Proposition 1.1.

We denote by $\pi^{M}$ the distribution of
$(\sum_{m=0}^{M}{\lambda_{m}}^{1/2}Z_{m}e_{m}(x))_{x\in\mathcal{X}}$ and
$\mathrm{S}=\\{\pi^{M}_{x_{1},\dots,x_{n}}\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}$. Then we have that $\mathrm{S}$ is _exchangeable_ and
_consistent_.

In order to sample from
$(\sum_{m=0}^{M}{\lambda_{m}}^{1/2}Z_{m}e_{m}(x))_{x\in\mathcal{X}}$ one only
needs to sample $\\{Z_{m}\\}_{m=0}^{M}$, effectively approximating the
infinite dimensional modelling problem by a finite-dimensional one. We are now
ready to present our method: Spectral Process Score-Based Generative Modelling
(SP-SGM).

#### Spectral Process Score-Based Generative modelling.

SP-SGM works by first choosing a kernel $\mathrm{K}$ and computing the
associated system of orthonormal eigenfunctions
$\\{e_{m}\\}_{m\in\mathbb{N}}$—see Sec. B.3. We assume we have access to a
dataset $\mathcal{D}=\\{\mathbf{Y}^{i}\\}_{i=1}^{L}$ of _functions_ 444In
practice we have access to these functions via a finite number of
observations. from the target distribution $\pi$. We set $M\in\mathbb{N}$ and
compute a new $M+1$-dimensional dataset
$\mathcal{D}^{M}=\\{\\{Z_{m}^{i}\\}_{m=0}^{M}\,:i\in\\{1,\dots,L\\}\\}$ where
$Z_{m}^{i}=\lambda_{m}^{-1/2}\int_{\mathcal{X}}\langle\mathbf{Y}_{x}^{i},e_{m}(x)\rangle\mathrm{d}x$.
Then, we train a $M+1$-dimensional standard SGM using $\mathcal{D}^{M}$. The
SGM model is initialised with
$\\{\bar{Z}_{m,0}\\}_{m=0}^{M}\sim\mathrm{N}(0,\operatorname{Id})$ and outputs
$\\{\bar{Z}_{m}\\}_{m=0}^{M}$ which is approximately distributed as
$\\{Z_{m}\\}_{m=0}^{M}$. Finally, SP-SGM outputs the _functional_ model
$\bar{F}_{M}=\sum_{m=0}^{M}\lambda_{m}^{1/2}\bar{Z}_{m}e_{m}$. We remark that
standard SGMs would only be able to deal with _finite_ input space. Note that
$\bar{F}_{M,0}=\sum_{m=0}^{M}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ is a
Gaussian process (GP) and in particular—when
$\mathrm{K}=\mathrm{K}_{\mathbf{Y}}$— it is the closest GP to the target
distribution $\pi$ in the following sense.

######  Proposition 1.2.

Let $\pi^{0}$ be the distribution of
$\sum_{m=0}^{+\infty}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ and $\pi$ the target
distribution. Denote $\mathrm{GP}(\mathcal{X})$ the space of Gaussian
processes on $\mathcal{X}$ and assume that $\mathrm{K}$ is the covariance
kernel. Then,
$\pi^{0}\in\operatornamewithlimits{argmin}_{\pi_{\mathrm{GP}}\in\mathrm{GP}(\mathcal{X})}\operatorname{KL}(\pi|\pi_{\mathrm{GP}})$.
In addition, $\sum_{m=0}^{M}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ is the
projection of $\sum_{m=0}^{+\infty}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ on the
subspace of $\mathrm{L}^{2}(\mathcal{X})$ spanned by $\\{e_{m}\\}_{m=0}^{M}$.

##  4 Related work

#### Gaussian processes and the Neural processes family.

One standard and powerful framework to construct distributions over functional
spaces are Gaussian processes (Rasmussen,, 2003). However, they scale badly
with the number of datapoints and are restricted in their modelling capacity.
These problems can be partially alleviated by using neural processes (Kim et
al.,, 2019; Garnelo et al., 2018b, ; Garnelo et al., 2018a, ; Jha et al.,,
2022; Louizos et al.,, 2019; Singh et al.,, 2019), although they also assume a
Gaussian likelihood and only satisfy predictive consistency, meaning one must
condition on a context set to achieve consistency. Recently, denoising
diffusion models have been combined with attention network architectures to
build consistent models of stochastic processes (Dutordoir et al.,, 2022).
Finally, Dupont et al., (2022) model weights of implicit neural representation
using SGM.

#### Spectral approaches.

Another possibility to circumvent the consistency problem is to consider
spectral decompositions—such as the Karhunen-Loéve one— which shows that
stochastic processes can be represented over a functional basis with random
coefficients. Some prior works have explored this avenue using variational
auto-encoders to model these coefficients (Mishra et al.,, 2020) whereas Lim
et al., (2022) considered energy-based models. Another related line of work is
‘inter-domain Gaussian processes’ (Lázaro-Gredilla and Figueiras-Vidal,,
2009), which rely on inducing variables that are obtained via inducing
functions such as elements of a Fourier basis (Hensman et al.,, 2018).

#### Spatial structure in score-based generative models

A variety of approaches have also been proposed to incorporate spatial
correlation in the noising process of finite-dimensional diffusion models
leveraging the multiscale structure of data (Jing et al.,, 2022; Guth et al.,,
2022; Ho et al.,, 2022; Saharia et al.,, 2021; Hoogeboom and Salimans,, 2022;
Rissanen et al.,, 2022). Our methodology can also be seen as a principled way
to modify the forward dynamics in classical denoising diffusion models. Indeed
by applying the diffusion in the spectral space the destruction process
_blurs_ the samples instead of _noising_ them, see Fig. 3 for an illustration.
Therefore, our contribution can be understood in the light of recent advances
in generative modelling on cold and soft denoising diffusion models (Daras et
al.,, 2022; Bansal et al.,, 2022; Hoogeboom and Salimans,, 2022).

##  5 Experimental results

First, we assess our method’s capacity to model bi-modality by qualitative
comparison of samples on the synthetic Quadratic dataset. Figure 3 shows that
our method is capable of expressing bi-modality while NPs are not. While SGMs
alone exhibit multimodality, they only apply to finite input spaces, hence
expressing multimodality in stochastic processes is indeed a feature of our
approach. To provide a quantitative assessment we compute the power of a two-
sample hypothesis test for functional data (Wynne and Duncan,, 2022) which
tries to discriminate between samples from the model and the dataset. Table 3
shows that our method outperforms GPs and NPs.

Then, we illustrate the behaviour of our proposed parametrisation of
stochastic processes. We plot the forward noising process of both SP-SGM
(reconstructed in the original image space) and of a standard SGM in Figure 3.
We see that by performing the diffusion in the spectral space, we capture and
incorporate spatial correlations governed by the spatial component of the
Karhunen-Loève theorem.

Finally we conduct ablation studies on the MNIST dataset. We consider the
effect of the truncation order $M$ in Figure 5, where we find best performance
for intermediate values of $M$. We decompose this effect into the
reconstruction between spectral and spatial domain, which improves with $M$,
and the quality of the stochastic model in the spectral domain, which
decreases with $M$. These conflicting objectives are also visualised in Figure
5. We also assess different kernels in Figure 5. Qualitatively, eigenfunctions
of the covariance kernel capture spatial correlations effectively and thus
generates the best samples, while the RBF kernel imposes a smoothness which
tends to produce blurry samples.

|
---|---
Figure 2:  Samples from the Quadratic dataset (orange), and from a trained NP
[Left, blue] and a trained SP-SGM [Right, blue].

| SP-SGM | NP | GP
---|---|---|---
Quadratic | $\bf{{5.4}_{\pm 0.7}}$ | ${8.6}_{\pm 1.5}$ | ${100.0}_{\pm 0.0}$
Melbourne | $\bf{{5.3}_{\pm 0.7}}$ | ${10.1}_{\pm 1.9}$ | ${20.1}_{\pm 4.0}$
Gridwatch | $\bf{{4.7}_{\pm 0.5}}$ | ${51.8}_{\pm 15.1}$ | ${29.2}_{\pm 5.5}$

tablePower (percent) of a kernel two-sample hypothesis test on 1D datasets.
Lower is better. Statistically significant best result is in bold.

Figure 3: Forward process in SP-SGM [Top] vs standard SGM [Bottom] on MNIST
digits. Pixel-wise mean and standard deviation of reference measure in
rightmost columns respectively.

---

figure Eigenfunctions [row 1, 3] and samples from SP-SGM [row 2, 4] for
covariance and RBF kernels.

---

$M=20$ | $M=100$ | $M=200$
---|---|---
|  |
|  |

figureOn the left, functional MMD vs $M$. Lower is better. On the right,
Karhunen-Loève recompositions of MNIST samples [top row] and the distribution
of the first spectral component from the dataset [orange] and from the SGM in
spectral space [blue].

## Acknowledgements

We thank the anonymous reviewers for their insightful and constructive
comments. AP and MH are funded through the StatML CDT through grant
EP/S023151/1. AD acknowledges support of the UK Defence Science and Technology
Laboratory (Dstl) and Engineering and Physical Research Council (EPSRC) under
grant EP/R013616/1. This is part of the collaboration between US DOD, UK MOD
and UK EPSRC under the Multidisciplinary University Research Initiative. AD is
also partially supported by the EPSRC grant EP/R034710/1 CoSines. EM is
supported by an EPSRC Prosperity Partnership EP/T005386/1 between Microsoft
Research and the University of Cambridge.

## References

  * Ambrosio et al., (2008) Ambrosio, L., Gigli, N., and Savaré, G. (2008).  Gradient Flows in Metric Spaces and in the Space of Probability Measures.  Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, second edition.
  * Baker, (1979) Baker, C. T. H. (1979).  Numerical Integration in the Treatment of Integral Equations.  In Hämmerlin, G., editor, Numerische Integration: Tagung Im Mathematischen Forschungsinstitut Oberwolfach Vom 1. Bis 7. Oktober 1978, pages 44–53. Birkhäuser Basel, Basel.
  * Bansal et al., (2022) Bansal, A., Borgnia, E., Chu, H.-M., Li, J. S., Kazemi, H., Huang, F., Goldblum, M., Geiping, J., and Goldstein, T. (2022).  Cold diffusion: Inverting arbitrary image transforms without noise.  arXiv preprint arXiv:2208.09392.
  * Cattiaux et al., (2021) Cattiaux, P., Conforti, G., Gentil, I., and Léonard, C. (2021).  Time reversal of diffusion processes under a finite entropy condition.  arXiv preprint arXiv:2104.07708.
  * Charalambos and Aliprantis, (2013) Charalambos, D. and Aliprantis, B. (2013).  Infinite Dimensional Analysis: A Hitchhiker’s Guide.  Springer-Verlag Berlin and Heidelberg GmbH & Company KG.
  * Daras et al., (2022) Daras, G., Delbracio, M., Talebi, H., Dimakis, A. G., and Milanfar, P. (2022).  Soft diffusion: Score matching for general corruptions.  arXiv preprint arXiv:2209.05442.
  * Dhariwal and Nichol, (2021) Dhariwal, P. and Nichol, A. (2021).  Diffusion models beat GAN on image synthesis.  arXiv preprint arXiv:2105.05233.
  * Dupont et al., (2022) Dupont, E., Kim, H., Eslami, S., Rezende, D., and Rosenbaum, D. (2022).  From data to functa: Your data point is a function and you should treat it like one.  arXiv preprint arXiv:2201.12204.
  * Dupuis and Ellis, (2011) Dupuis, P. and Ellis, R. S. (2011).  A weak convergence approach to the theory of large deviations.  John Wiley & Sons.
  * Dutordoir et al., (2022) Dutordoir, V., Saul, A., Ghahramani, Z., and Simpson, F. (2022).  Neural diffusion processes.  arXiv preprint arXiv:2206.03992.
  * Fernique, (1975) Fernique, X. (1975).  Regularité des trajectoires des fonctions aléatoires gaussiennes.  In Ecole d’Eté de Probabilités de Saint-Flour IV—1974, pages 1–96. Springer.
  * Ferreira and Menegatto, (2009) Ferreira, J. and Menegatto, V. (2009).  Eigenvalues of integral operators defined by smooth positive definite kernels.  Integral Equations and Operator Theory, 64(1):61–81.
  * (13) Garnelo, M., Rosenbaum, D., Maddison, C., Ramalho, T., Saxton, D., Shanahan, M., Teh, Y. W., Rezende, D., and Eslami, S. M. A. (2018a).  Conditional neural processes.  In Dy, J. and Krause, A., editors, Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pages 1704–1713. PMLR.
  * (14) Garnelo, M., Schwarz, J., Rosenbaum, D., Viola, F., Rezende, D. J., Eslami, S., and Teh, Y. W. (2018b).  Neural processes.  arXiv preprint arXiv:1807.01622.
  * Guth et al., (2022) Guth, F., Coste, S., De Bortoli, V., and Mallat, S. (2022).  Wavelet score-based generative modeling.  arXiv preprint arXiv:2208.05003.
  * Haussmann and Pardoux, (1986) Haussmann, U. G. and Pardoux, E. (1986).  Time reversal of diffusions.  The Annals of Probability, 14(4):1188–1205.
  * Hensman et al., (2018) Hensman, J., Durrande, N., and Solin, A. (2018).  Variational fourier features for gaussian processes.  Journal of Machine Learning Research, Volume 18, issue 1:52.
  * Ho et al., (2020) Ho, J., Jain, A., and Abbeel, P. (2020).  Denoising diffusion probabilistic models.  Advances in Neural Information Processing Systems.
  * Ho et al., (2022) Ho, J., Saharia, C., Chan, W., Fleet, D. J., Norouzi, M., and Salimans, T. (2022).  Cascaded diffusion models for high fidelity image generation.  J. Mach. Learn. Res., 23:47–1.
  * Hoogeboom and Salimans, (2022) Hoogeboom, E. and Salimans, T. (2022).  Blurring diffusion models.  arXiv preprint arXiv:2209.05557.
  * Hyvärinen, (2005) Hyvärinen, A. (2005).  Estimation of non-normalized statistical models by score matching.  Journal of Machine Learning Research, 6(4).
  * Jha et al., (2022) Jha, S., Gong, D., Wang, X., Turner, R. E., and Yao, L. (2022).  The neural process family: Survey, applications and perspectives.  arXiv preprint arXiv:2209.00517.
  * Jing et al., (2022) Jing, B., Corso, G., Berlinghieri, R., and Jaakkola, T. (2022).  Subspace diffusion generative models.  arXiv preprint arXiv:2205.01490.
  * Kim et al., (2019) Kim, H., Mnih, A., Schwarz, J., Garnelo, M., Eslami, A., Rosenbaum, D., Vinyals, O., and Teh, Y. W. (2019).  Attentive neural processes.  In International Conference on Learning Representations.
  * Lázaro-Gredilla and Figueiras-Vidal, (2009) Lázaro-Gredilla, M. and Figueiras-Vidal, A. (2009).  Inter-domain gaussian processes for sparse inference using inducing features.  In Bengio, Y., Schuurmans, D., Lafferty, J., Williams, C., and Culotta, A., editors, Advances in Neural Information Processing Systems, volume 22. Curran Associates, Inc.
  * Lim et al., (2022) Lim, J. N., Vollmer, S., Wolf, L., and Duncan, A. (2022).  $\mathcal{F}$-ebm: Energy based learning of functional data.  arXiv preprint arXiv:2202.01929.
  * Louizos et al., (2019) Louizos, C., Shi, X., Schutte, K., and Welling, M. (2019).  The functional neural process.  Advances in Neural Information Processing Systems, 32.
  * Mishra et al., (2020) Mishra, S., Flaxman, S., Berah, T., Pakkanen, M., Zhu, H., and Bhatt, S. (2020).  $\pi$vae: Encoding stochastic process priors with variational autoencoders.  arXiv preprint arXiv:2002.06873.
  * Rasmussen, (2003) Rasmussen, C. E. (2003).  Gaussian processes in machine learning.  In Summer school on machine learning, pages 63–71. Springer.
  * Ravuri et al., (2021) Ravuri, S., Lenc, K., Willson, M., Kangin, D., Lam, R., Mirowski, P., Fitzsimons, M., Athanassiadou, M., Kashem, S., Madge, S., Prudden, R., Mandhane, A., Clark, A., Brock, A., Simonyan, K., Hadsell, R., Robinson, N., Clancy, E., Arribas, A., and Mohamed, S. (2021).  Skilful precipitation nowcasting using deep generative models of radar.  Nature, 597(7878):672–677.
  * Rissanen et al., (2022) Rissanen, S., Heinonen, M., and Solin, A. (2022).  Generative Modelling With Inverse Heat Dissipation.
  * Saharia et al., (2021) Saharia, C., Ho, J., Chan, W., Salimans, T., Fleet, D. J., and Norouzi, M. (2021).  Image super-resolution via iterative refinement.  arXiv preprint arXiv:2104.07636.
  * Singh et al., (2019) Singh, G., Yoon, J., Son, Y., and Ahn, S. (2019).  Sequential neural processes.  Advances in Neural Information Processing Systems, 32.
  * Song and Ermon, (2019) Song, Y. and Ermon, S. (2019).  Generative modeling by estimating gradients of the data distribution.  In Advances in Neural Information Processing Systems.
  * Song et al., (2021) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B. (2021).  Score-based generative modeling through stochastic differential equations.  In International Conference on Learning Representations.
  * Sun et al., (2019) Sun, S., Zhang, G., Shi, J., and Grosse, R. (2019).  Functional variational bayesian neural networks.  arXiv preprint arXiv:1903.05779.
  * Vincent, (2011) Vincent, P. (2011).  A connection between score matching and denoising autoencoders.  Neural Computation, 23(7):1661–1674.
  * Wynne and Duncan, (2022) Wynne, G. and Duncan, A. B. (2022).  A kernel two-sample test for functional data.  Journal of Machine Learning Research, 23(73):1–51.
  * Yang et al., (2020) Yang, L., Zhang, D., and Karniadakis, G. E. (2020).  Physics-informed generative adversarial networks for stochastic differential equations.  SIAM Journal on Scientific Computing, 42(1):A292–A317.

##  Appendix A Introduction to the supplementary material

We provide details on the implementation of SP-SGM in Appendix B. The proof of
Proposition 1.1 is presented in Appendix D. The proof of Proposition 1.2 is
presented in Appendix E. The likelihood computation is given in Appendix C.
Experimental details are provided in Appendix F.

##  Appendix B SP-SGM in practice

###  B.1 SP-SGM algorithm

Algorithm 1 Spectral Process Score-Based Generative Model (SP-SGM)

1:$T,\mathcal{D},\theta_{0},N_{\mathrm{iter}},\varepsilon,K,\\{(\lambda_{m},e_{m})\\}_{m=0}^{M}$

2:/// TRAINING ///

3:Get $\mathcal{D}^{M}$ from $\mathcal{D}$ $\triangleright$ Dataset projection

4:for $n\in\\{0,\dots,N_{\textrm{iter}}-1\\}$ do

5: Get $\\{\mathbf{Y}_{m,0}\\}_{m=0}^{M}$ mini-batch from $\mathcal{D}^{M}$

6: $t\sim U(\left[\varepsilon,T\right])$ $\triangleright$ Uniform sampling
between $\varepsilon$ and $T$

7:
$\mathbf{Y}_{m,t}=\mathrm{e}^{-t}\mathbf{Y}_{m,0}+(1-\mathrm{e}^{-2t})^{1/2}G,\quad
G\sim\mathrm{N}(0,\operatorname{Id})$ $\triangleright$ Diffuse

8: Get DSM loss $\ell(\theta_{n})$ $\triangleright$ Compute score matching
loss

9: $\theta_{n+1}=\verb|optimiser_update|(\theta_{n},\ell(\theta_{n}))$
$\triangleright$ ADAM optimiser step

10:end for

11:$\theta^{\star}=\theta_{N_{\textrm{iter}}}$

12:/// SAMPLING ///

13:$\\{\bar{Y}_{m,0}\\}_{m=0}^{M}\sim\mathrm{N}(0,\operatorname{Id})$
$\triangleright$ Sample from Gaussian distribution

14:$b_{\theta}^{\star}(t,y)=\mathbf{s}_{\theta^{\star}}(T-t,\bar{y})$ for any
$t\in\left[0,T\right]$, $\bar{y}\in\mathbb{R}^{M+1}$ $\triangleright$ Reverse
process drift

15:$\\{\bar{Y}_{m,n}\\}_{m=0,n=0}^{M,N}$ Euler-Maruyama with drift
$b_{\theta}^{\star}$ $\triangleright$ Approximate reverse diffusion

16:return
$\theta^{\star},x\mapsto\sum_{m=0}^{M}\lambda_{m}^{1/2}\bar{Y}_{m,N}e_{m}(x)$

###  B.2 Karhunen-Loève coefficients

In this section we show that the random coefficients from the KL basis
expansion, given for any $m\in\mathbb{N}$ by

$\textstyle{Z_{m}=\langle\mathbf{Y},\lambda_{m}^{-1/2}e_{m}\rangle_{\mathrm{L}^{2}(\mathcal{X})}=\int_{\mathcal{X}}\langle(\mathbf{Y}(x)-\mu(x)),\lambda_{m}^{-1/2}e_{m}(x)\rangle\mathrm{d}x,}$
(5)

satisfy $\mathbb{E}[Z_{m}]=0$ and furthermore if
$\mathrm{K}=\mathrm{K}_{\mathbf{Y}}$ then
$\mathbb{E}[Z_{m}Z_{m^{\prime}}]=\updelta_{m,m^{\prime}}$. First, we have that
for any $m\in\mathbb{N}$

$\textstyle{\mathbb{E}[Z_{m}]=\int_{\mathcal{X}}\langle\mathbb{E}[\mathbf{Y}(x)-\mu(x)],\lambda_{m}^{-1/2}e_{m}(x)\rangle\mathrm{d}x=0.}$
(6)

Then, we have for any $m,m^{\prime}\in\mathbb{N}$

$\displaystyle\mathbb{E}[Z_{m}Z_{m^{\prime}}]$
$\displaystyle=\textstyle{\mathbb{E}[\int_{\mathcal{X}}\langle(\mathbf{Y}(x)-\mu(x)),\lambda_{m}^{-1/2}e_{m}(x)\rangle\mathrm{d}x\int_{\mathcal{X}}\langle(\mathbf{Y}(x^{\prime})-\mu(x^{\prime})),\lambda_{m^{\prime}}^{-1/2}e_{m^{\prime}}(x^{\prime})\rangle\mathrm{d}x^{\prime}]}$
(7a)
$\displaystyle=\lambda_{m}^{-1/2}\lambda_{m^{\prime}}^{-1/2}\textstyle{\int_{\mathcal{X}}\int_{\mathcal{X}}\operatorname{Tr}(\mathbb{E}[(\mathbf{Y}(x)-\mu(x))(\mathbf{Y}(x^{\prime})-\mu(x^{\prime}))^{\top}]e_{m}(x)e_{m^{\prime}}^{\top}(x^{\prime}))\mathrm{d}x\mathrm{d}x^{\prime}}$
(7b)
$\displaystyle=\lambda_{m}^{-1/2}\lambda_{m^{\prime}}^{-1/2}\textstyle{\int_{\mathcal{X}}e_{m}(x)^{\top}\mathrm{K}_{\mathbf{Y}}(x,x^{\prime})e_{m^{\prime}}(x^{\prime})\mathrm{d}x\mathrm{d}x^{\prime}}$
(7c)
$\displaystyle=\lambda_{m}^{-1/2}\lambda_{m^{\prime}}^{-1/2}\textstyle{\int_{\mathcal{X}}\lambda_{m}\langle
e_{m}(x^{\prime}),e_{m^{\prime}}(x^{\prime})\rangle\mathrm{d}x^{\prime}}=\lambda_{m}^{-1/2}\lambda_{m^{\prime}}^{-1/2}\lambda_{m}\updelta_{m,m^{\prime}}=\updelta_{m,m^{\prime}}.$
(7d)

###  B.3 Computation of eigenfunctions

We focus on the case where $\mathbf{Y}$ takes values in $\mathbb{R}$ but the
computations can be extended to the case where $\mathbf{Y}$ takes values in
$\mathbb{R}^{d}$ with $d\in\mathbb{N}$. Recall from Theorem 1 the
$\mathrm{L}^{2}(\mathcal{X})$ operator $\mathrm{T}_{\mathrm{K}}$ associated
with kernel $\mathrm{K}$, given by
$\mathrm{T}_{\mathrm{K}}f(x_{2})=\int_{\mathcal{X}}\mathrm{K}(x_{1},x_{2})f(x_{1})\mathrm{d}x_{1}$.
In the following we describe how we approximate the orthonormal system of
eigenfunctions $\\{e_{m}\\}_{m\in\mathbb{N}}$ of $\mathrm{T}_{\mathrm{K}}$
from a dataset $\mathcal{D}=\\{\mathbf{Y}^{i}\\}_{i=1}^{L}$.

Firstly, we note that in practice we have access to the functions in our
dataset at a finite set of query points, that is
$\mathcal{D}=\\{\\{\mathbf{Y}^{i}_{x_{i,n}}\\}_{n=1}^{N_{i}}\\}_{i=1}^{L}$.
For ease of notation, in the following we denote
$\\{x_{s}\\}_{s=1}^{S}=\\{x_{i,n}\,:i\in\\{1,\dots,L\\},\ n\in\\{1,\dots
N_{i}\\}\\}$, with $s=\sum_{i=1}^{L}N_{i}$.

The eigensystem $\\{e_{m}\\}_{m\in\mathbb{N}}$ satisfies
$\mathrm{T}_{\mathrm{K}}e_{m}(x)=\lambda_{m}e_{m}(x)$ for all
$x\in\mathcal{X}$, which we approximate by

$\textstyle{\lambda_{m}e_{m}(x)=\mathrm{T}_{\mathrm{K}}e_{m}(x)\approx\frac{1}{S}\sum_{s=1}^{S}\mathrm{K}(x,x_{s})e_{m}(x_{s}).}$
(8)

Evaluating the above at the entire set of evaluation points
$x=\\{x_{s}\\}_{s=1}^{S}$ results in the matrix eigenproblem

$\textstyle{\frac{1}{S}Ku_{m}=\lambda_{m}^{\text{mat}}u_{m},}$

where $K$ is an $S\times S$ matrix and has entries
$K_{i,j}=\mathrm{K}(x_{i},x_{j})$, $\lambda_{m}^{\text{mat}}\in\mathbb{R}$ is
the matrix eigenvalue and $u_{m}\in\mathbb{R}^{S}$ is the corresponding matrix
eigenvector, normalised such that $u_{m}^{\top}u_{m}=1$. By Baker’s theorem
(Baker,, 1979, Theorem 3.4), $\lambda_{m}^{\text{mat}}\to\lambda_{m}$ as
$S\to\infty$. Note that if we take $e_{m}(x_{s})=(u_{m})_{s}$ then we have

$\textstyle{\int_{\mathcal{X}}e_{m}(x)e_{m}(x)p(x)\mathrm{d}x\approx\frac{1}{S}\sum_{s=1}^{S}e_{m}(x_{s})e_{m}(x_{s})=\frac{1}{S}u_{m}^{\top}u_{m}=\frac{1}{S}}$

where $p(x)$ is the density of observation points on $\mathcal{X}$. So, in
order to have an orthonormal eigensystem
$\int_{\mathcal{X}}e_{m}(x)e_{m^{\prime}}(x)p(x)\mathrm{d}x=\updelta_{mm^{\prime}}$
we must rescale according to $e_{m}(x_{s})=\sqrt{S}(u_{m})_{s}$. Finally,
rearranging Eq. 8 we get

$\textstyle{\hat{e}_{j}(x)\approx(S\lambda_{j}^{\text{mat}})^{-1}\sum_{s=1}^{S}k(x,x_{s})e_{j}(x_{s})\approx(\sqrt{S}\lambda_{j}^{\text{mat}})^{-1}\sum_{s=1}^{S}k(x,x_{s})u_{j}(x_{s}).}$
(9)

Having obtained an approximate eigensystem as above, it is simple to obtain
the spectral dataset
$\mathcal{D}^{M}=\\{\\{Z_{m}^{i}\\}_{m=0}^{M}\,:i\in\\{1,\dots,L\\}\\}$ using

$\textstyle{Z_{m}^{i}\approx\tfrac{1}{N_{i}}\sum_{n=1}^{N_{i}}(\mathbf{Y}_{x_{n}}-\mu(x_{n}))(\lambda_{m}^{\text{mat}})^{-1/2}\hat{e}_{m}(x_{n})}.$

where $\mu(x)=\mathbb{E}[\mathbf{Y}_{x}]$ is the process mean.

###  B.4 Truncation order

For a given approximate eigensystem
$\\{(\hat{\lambda}_{m},\hat{e}_{m})\\}^{L}_{m=1}$, we define the truncation
order $M$ as

$\textstyle{M={\arg\min}_{m\in\mathbb{N}}\sum_{j=1}^{m}\hat{\lambda}_{j}/\sum_{j=1}^{L}\hat{\lambda}_{j}\geq\eta,}$

where $\eta\in[0,1]$ is a threshold parameter, fixed to $\eta=0.99$ for our
experimental results on the 1D datasets. Alternatively, as we do on the MNIST
dataset, one can treat the truncation order $M$ as an hyperparameter.

###  B.5 Eigenvalue weighted objective

We propose the following weighted DSM loss to train the score network
$\bm{s}_{\theta}:[0,T]\times\mathbb{R}^{M+1}\rightarrow\mathbb{R}^{M+1}$:

$\textstyle{\mathbb{E}_{t}\left\\{\mathbb{E}_{\mathbf{Z}_{0},\mathbf{Z}_{t}}\left[\left\|\tilde{\lambda}^{\alpha}\odot\left(\bm{s}_{\theta}(\mathbf{Z}_{t},t)-\nabla_{\mathbf{Z}_{t}}\log
p(\mathbf{Z}_{t}|\mathbf{Z}_{0})\right)\right\|^{2}_{2}\right]\right\\},}$
(10)

with
$\tilde{\lambda}=[\lambda_{0}/\Lambda,\dots,\lambda_{M}/\Lambda]^{\top}$555The
exponentiation in Eq. 10 is to be understood in a pointwise fashion.,
$\Lambda=\sum_{m=0}^{M}\lambda_{m}$, the vector of normalised eigenvalues,
$\alpha>0$ a tuneable hyperparameter, $t\sim U([0,T])$,
$\mathbf{Z}_{0}=\\{\lambda_{m}^{-1/2}\int_{\mathcal{X}}\langle\mathbf{Y}_{x},e_{m}(x)\rangle\mathrm{d}x\\}_{m=0}^{M}$
and

$\mathrm{d}\mathbf{Z}_{t}=-(1/2)\mathbf{Z}_{t}\mathrm{d}t+\beta_{t}\mathrm{d}\mathbf{B}_{t},$
(11)

where $t\mapsto\beta_{t}$ is an hyperparameter. The motivation behind Eq. 10
is to put more importance on fitting well the lower frequency (i.e. higher
eigenvalue) components of the decomposition
$\sum_{m=0}^{M}{\lambda_{m}}^{1/2}Z_{m}e_{m}(x)$ since these matter the most
in terms of quality of reconstruction as per Theorem 1. With $\alpha=0$, we
have that $\lambda^{\alpha}=\bm{1}$, thus that all dimensions of the spectral
space are weighted equally and we recover the standard DSM loss.

##  Appendix C Likelihood evaluation

Let’s denote the Karhunen-Loeve expansion with
$\phi:\mathbb{R}^{M+1}\rightarrow\mathbb{R}^{d}$ such that
$(Z_{1},\dots,Z_{M})\mapsto\mu(\cdot)+\sum_{m=0}^{M}Z_{m}\sqrt{\lambda_{m}}e_{m}(\cdot)$.
We have that $\frac{\partial\phi}{\partial Z_{m}}=\sqrt{\lambda_{m}}e_{m}$,
hence with $J_{z}\phi=\left(\frac{\partial\phi}{\partial
z_{1}},\dots,\frac{\partial\phi}{\partial z_{N}}\right)$ we get
$J_{z}\phi^{\top}J_{z}\phi=\left(\sqrt{\lambda_{m}}\sqrt{\lambda_{n}}\langle
e_{m},e_{n}\rangle\right)_{mn}=\left(\sqrt{\lambda_{m}}\sqrt{\lambda_{n}}\delta_{mn}\right)_{mn}=\text{diag}(\lambda_{m})$.
Hence, denoting a realisation $F=\phi(\mathbf{Z})$, applying the change of
variable formula we have

$\textstyle{\log p(F)=\sum_{m}\log
p(Z_{m})-\log|J_{z}\phi^{\top}J_{z}\phi|^{1/2}=\sum_{m}\left(\log
p(Z_{m})-1/2\log\lambda_{m}\right).}$ (12)

##  Appendix D Proof of Proposition 1.1

We recall Proposition 1.1. See 1.1

###### Proof.

Note that this result is the trivial direction of the Kolmogorov extension
theorem (Charalambos and Aliprantis,, 2013, Theorem 15.26) since $F_{M}$ is by
definition a stochastic process. For completeness, we provide a proof of this
result in our setting. We denote $Z=\\{Z_{m}\\}_{m=0}^{M}$ and $\pi$ the model
distribution on $Z$. For any $n\in\mathbb{N}$ and
$f\in\mathrm{C}_{b}((\mathbb{R}^{d})^{n},\mathbb{R})$ we have

$\displaystyle\textstyle{\int
f(y_{1},\dots,y_{n})\mathrm{d}\pi^{M}_{x_{1},\dots,x_{n}}(y_{1},\dots,y_{n})}$
$\displaystyle=\textstyle{\int\int_{\mathbb{R}^{M+1}}f(y_{1},\dots,y_{n})\mathrm{d}\pi^{M}_{x_{1},\dots,x_{n}}(y_{1},\dots,y_{n},Z=z)\mathrm{d}z}$
(13)
$\displaystyle=\textstyle{\int\int_{\mathbb{R}^{M+1}}f(y_{1},\dots,y_{n})\mathrm{d}\pi^{M}_{x_{1},\dots,x_{n}}(y_{1},\dots,y_{n},|Z=z)\mathrm{d}\pi(z)}$
(14)
$\displaystyle=\textstyle{\int\int_{\mathbb{R}^{M+1}}f(y_{1},\dots,y_{n})\prod_{i=1}^{n}\mathrm{d}\pi^{M}_{x_{i}}(y_{i}|Z=z)\mathrm{d}\pi(z).}$
(15)

Hence, for any $n\in\mathbb{N}$, $n$-permutation $\sigma\in\mathfrak{S}_{n}$,
$x_{1},\dots,x_{n}\in\mathcal{X}$, and continuous and bounded function
$f\in\mathrm{C}_{b}((\mathbb{R}^{d})^{n})$ we have

$\textstyle{\int
f(y_{\sigma^{-1}(1)},\dots,y_{\sigma^{-1}(n)})~{}\mathrm{d}\pi^{M}_{x_{\sigma(1)},\dots,x_{\sigma(n)}}(y_{1},\dots,y_{n})=\int
f(y_{1},\dots,y_{n})~{}\mathrm{d}\pi^{M}_{x_{1},\dots,x_{n}}(y_{1},\dots,y_{n}).}$

Similarly, for any $n_{1}\leq n_{2}$, $x_{1},\dots,x_{n_{2}}\in\mathcal{X}$
and $f\in\mathrm{C}_{b}((\mathbb{R}^{d})^{n_{1}})$ we have

$\displaystyle\textstyle{\int
f(y_{1},\dots,y_{n_{1}})~{}\mathrm{d}\pi^{M}_{x_{1},\dots,x_{n_{2}}}(y_{1},\dots,y_{n_{2}})}$
(16) $\displaystyle\qquad=\textstyle{\int
f(y_{1},\dots,y_{n_{1}})\int_{\mathbb{R}^{M+1}}~{}\mathrm{d}\pi^{M}_{x_{1},\dots,x_{n_{2}}}(y_{1},\dots,y_{n_{2}}|Z=z)\mathrm{d}\pi(z)}$
(17) $\displaystyle\qquad=\textstyle{\int
f(y_{1},\dots,y_{n_{1}})\int_{\mathbb{R}^{M+1}}~{}\prod_{i=1}^{n_{2}}\mathrm{d}\pi^{M}_{x_{i}}(y_{i}|Z=z)\mathrm{d}\pi(z)}$
(18) $\displaystyle\qquad=\textstyle{\int
f(y_{1},\dots,y_{n_{1}})\int_{\mathbb{R}^{M+1}}~{}\prod_{i=1}^{n_{1}}\mathrm{d}\pi^{M}_{x_{i}}(y_{i}|Z=z)\left(\prod_{i=n_{1}+1}^{n_{2}}\mathrm{d}\pi^{M}_{x_{i}}(y_{i}|Z=z)\right)\mathrm{d}\pi(z)}$
(19) $\displaystyle\qquad=\textstyle{\int
f(y_{1},\dots,y_{n_{1}})\int_{\mathbb{R}^{M+1}}~{}\prod_{i=1}^{n_{1}}\mathrm{d}\pi^{M}_{x_{i}}(y_{i}|Z=z)\mathrm{d}\pi(z)}$
(20) $\displaystyle\qquad=\textstyle{\int
f(y_{1},\dots,y_{n_{1}})~{}\mathrm{d}\pi^{M}_{x_{1},\dots,x_{n_{1}}}(y_{1},\dots,y_{n_{1}}).}$
(21)

In this setting, the random variable $Z$ is finite dimensional as it is
supported in $\mathbb{R}^{M}$, but more generally this result is still true
with an (infinite dimensional) process $F$. ∎

##  Appendix E Proof of Proposition 1.2

In this section, we provide two proofs of Proposition 1.2.

###  E.1 Partition definition

We recall Proposition 1.2.

######  Proposition 1.3.

Let $\pi^{0}$ be the distribution of
$\sum_{m=0}^{+\infty}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ and $\pi$ the target
distribution. Denote $\mathrm{GP}(\mathcal{X})$ the space of Gaussian
processes on $\mathcal{X}$ and assume that $\mathrm{K}$ is the covariance
kernel. Then,
$\pi^{0}\in\operatornamewithlimits{argmin}_{\pi_{\mathrm{GP}}\in\mathrm{GP}(\mathcal{X})}\operatorname{KL}(\pi|\pi_{\mathrm{GP}})$.
In addition, $\sum_{m=0}^{M}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ is the
projection of $\sum_{m=0}^{+\infty}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ on the
subspace of $\mathrm{L}^{2}(\mathcal{X})$ spanned by $\\{e_{m}\\}_{m=0}^{M}$.

Note that in this first approach we do not assume any regularity on the
samples of the process. This proof is based on (Sun et al.,, 2019, Theorem 1)
which considered a modified definition of the Kullback-Leibler divergence.
Namely, for any $\mu,\nu$ probability measures over a probability space
$(\Omega,\mathcal{F})$ we define

$\textstyle{\operatorname{KL}(\mu|\nu)=\sup\\{\operatorname{KL}(\mu_{P}|\nu_{P})\,:P\
\text{ measurable finite partition of $\Omega$}\\}},$ (22)

where for any measurable finite partition
$P=\\{\Omega_{1},\dots,\Omega_{N}\\}$ we define
$\mu_{P}=\\{\mu(\Omega_{1}),\dots,\mu(\Omega_{N})\\}$, similarly for
$\nu_{P}$. With that definition (Sun et al.,, 2019, Theorem 1) holds. Namely,
for any $\mu,\nu\in\mathscr{P}((\mathbb{R}^{d})^{\mathcal{X}})$666The product
space $(\mathbb{R}^{d})^{\mathcal{X}}$ is a measurable space with the
cylindrical sigma algebra, see Sun et al., (2019) for more details in that
context., we have

$\textstyle{\operatorname{KL}(\mu|\nu)=\sup\\{\operatorname{KL}(\mu_{x_{1},\dots,x_{n}}|\nu_{x_{1},\dots,x_{n}})\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}.}$ (23)

Let $\pi$ be the target distribution of the stochastic process. We define
$\mathrm{GP}(\mathcal{X})$ the space of Gaussian processes over the input
space $\mathcal{X}$. We have that

$\displaystyle\inf\\{\sup\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}\,:\pi^{0}\in\mathrm{GP}(\mathcal{X})\\}$ (24)
$\displaystyle\qquad\geq\sup\\{\inf\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:\pi^{0}\in\mathrm{GP}(\mathcal{X})\\}\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}.$ (25)

In addition, we have

$\displaystyle\sup\\{\inf\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:\pi^{0}\in\mathrm{GP}(\mathcal{X})\\}\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}$ (26)
$\displaystyle\qquad=\sup\\{\inf\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:\pi^{0}_{x_{1},\dots,x_{n}}\in\mathrm{GP}(\\{x_{1},\dots,x_{n}\\})\\}\,:$
(27) $\displaystyle\qquad\qquad(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\},$

where we emphasise that the set $\mathrm{GP}(\\{x_{1},\dots,x_{n}\\})$ is
simply the set of $nd$-dimensional Gaussian probability measures. Hence, for
any $\\{x_{1},\dots,x_{n}\\}\in\mathcal{X}^{n}$ and $n\in\mathbb{N}$, we have
that

$\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0,\star}_{x_{1},\dots,x_{n}})=\inf\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:\pi^{0}_{x_{1},\dots,x_{n}}\in\mathrm{GP}(\\{x_{1},\dots,x_{n}\\})\\},$
(28)

where $\pi^{0,\star}_{x_{1},\dots,x_{n}}$ is the Gaussian measure with same
mean and covariance matrix as $\pi_{x_{1},\dots,x_{n}}$. Using the Kolmogorov
extension theorem (Charalambos and Aliprantis,, 2013, Theorem 15.26), there
exists $\pi^{\star}\in\mathscr{P}((\mathbb{R}^{d})^{\mathcal{X}})$ such that
for any $\\{x_{1},\dots,x_{n}\\}\in\mathcal{X}^{n}$ and $n\in\mathbb{N}$,
$\pi^{\star}_{x_{1},\dots,x_{n}}=\pi^{0,\star}_{x_{1},\dots,x_{n}}$.
Therefore, we have that

$\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{\star}_{x_{1},\dots,x_{n}})=\inf\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:\pi_{x_{1},\dots,x_{n}}\in\mathrm{GP}(\mathcal{X})\\}.$
(29)

Hence, using Eq. 23 we get that

$\displaystyle\sup\\{\inf\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:\pi^{0}\in\mathrm{GP}(\mathcal{X})\\}\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}$ (30)
$\displaystyle\qquad=\sup\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{\star}_{x_{1},\dots,x_{n}})\,:\\{x_{1},\dots,x_{n}\\}\in\mathcal{X}^{n}\\}=\operatorname{KL}(\pi|\pi^{\star}).$
(31)

Therefore, we have that

$\inf\\{\sup\\{\operatorname{KL}(\pi_{x_{1},\dots,x_{n}}|\pi^{0}_{x_{1},\dots,x_{n}})\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}\,:\pi^{0}\in\mathrm{GP}(\mathcal{X})\\}\geq\operatorname{KL}(\pi|\pi^{\star}),$
(32)

which implies using Eq. 23 that

$\operatorname{KL}(\pi|\pi^{\star})\leq\inf\\{\operatorname{KL}(\pi|\pi^{0})\,:\pi^{0}\in\mathrm{GP}(\mathcal{X})\\}.$
(33)

The equality holds since $\pi^{\star}\in\mathrm{GP}(\mathcal{X})$. Finally,
since $\pi^{\star}$ and $\pi$ have the same covariance kernels, they have the
same Karhunen-Loève eigensystems. Therefore, there exists
$\\{Z_{m}\\}_{m\in\mathbb{N}}$ such that
$\mathbf{Y}^{\star}=\sum_{m\in\mathbb{N}}\lambda_{m}^{1/2}Z_{m}e_{m}$ has
distribution $\pi^{\star}$, which concludes the proof.

###  E.2 Sample continuous

In our second approach, we restrict ourselves to the case of sample continuous
processes. Namely, we now longer consider
$\mu\in\mathscr{P}((\mathbb{R}^{d})^{\mathcal{X}})$ but
$\mu\in\mathscr{P}(\mathrm{C}(\mathcal{X},\mathbb{R}^{d}))$.

######  Proposition 1.4.

Assume that $\pi\in\mathrm{C}(\mathcal{X},\mathbb{R}^{d})$ and that there
exists $\phi:\ \left[0,+\infty\right)$ such that for any
$x_{1},x_{2}\in\mathcal{X}$,

$\textstyle{\mathbb{E}[\|\mathbf{Y}_{x_{1}}-\mathbf{Y}_{x_{2}}\|^{2}]\leq\phi(\|x_{1}-x_{2}\|),}$
(34)

such that $\int_{0}^{+\infty}\phi(\exp[-t^{2}])\mathrm{d}t<+\infty$. Let
$\pi^{0}$ be the distribution of
$\sum_{m=0}^{+\infty}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$. Denote
$\mathrm{GP}(\mathcal{X})$ the space of Gaussian processes on $\mathcal{X}$
and assume that $\mathrm{K}$ is the covariance kernel. Then,
$\pi^{0}\in\operatornamewithlimits{argmin}_{\pi_{\mathrm{GP}}\in\mathrm{GP}(\mathcal{X})}\operatorname{KL}(\pi|\pi_{\mathrm{GP}})$.
In addition, $\sum_{m=0}^{M}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ is the
projection of $\sum_{m=0}^{+\infty}\lambda_{m}^{1/2}\bar{Z}_{m,0}e_{m}$ on the
subspace of $\mathrm{L}^{2}(\mathcal{X})$ spanned by $\\{e_{m}\\}_{m=0}^{M}$.

In that case, we show that

$\textstyle{\operatorname{KL}(\mu|\nu)=\sup\\{\operatorname{KL}(\mu_{x_{1},\dots,x_{n}}|\nu_{x_{1},\dots,x_{n}})\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}.}$ (35)

More precisely, we have the following proposition.

######  Proposition 1.5.

Let $\mu,\nu\in\mathscr{P}(\mathrm{C}(\mathcal{X},\mathbb{R}^{d}))$

$\textstyle{\operatorname{KL}(\mu|\nu)=\sup\\{\operatorname{KL}(\mu_{x_{1},\dots,x_{n}}|\nu_{x_{1},\dots,x_{n}})\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}.}$ (36)

###### Proof.

First, note that for any $\\{x_{1},\dots,x_{n}\\}\in\mathcal{X}^{n}$,
$n\in\mathbb{N}$ we have using the data processing theorem (Ambrosio et al.,,
2008, Lemma 9.4.5)

$\operatorname{KL}(\mu|\nu)\geq\operatorname{KL}(\mu_{x_{1},\dots,x_{n}}|\nu_{x_{1},\dots,x_{n}}).$
(37)

Since $\mathcal{X}$ is compact, for any $n\in\mathbb{N}$, there exists
$\\{x_{1},\dots,x_{n}\\}$ such that
$\mathcal{X}\subset\cup_{k=1}^{n}\operatorname{B}(x_{k},1/(n+1))$. For any
$n\in\mathbb{N}$, let $\\{\varphi_{k}\\}_{k=1}^{n}$ be a smooth partition of
unity associated with $\\{\operatorname{B}(x_{k},1/(n+1))\\}_{k=1}^{n}$. For
any $\mathrm{C}(\mathcal{X},\mathbb{R}^{d})$-valued random variable
$\mathbf{Y}$ and $n\in\mathbb{N}$, denote $\mathbf{Y}^{n}$ such that for any
$x\in\mathcal{X}$

$\textstyle{\mathbf{Y}^{n}_{x}=\sum_{k=1}^{n}\varphi_{k}(x)\mathbf{Y}_{x_{k}}.}$
(38)

Since $\mathbf{Y}$ is continuous and $\mathcal{X}$ is compact we have that
$\mathbf{Y}$ is uniformly continuous. Hence, for any $\varepsilon>0$, there
exists $n\in\mathbb{N}$ such that for any $x_{1},x_{2}$, $\|x_{1}-x_{2}\|\leq
1/(n+1)$, $\|\mathbf{Y}_{x_{1}}-\mathbf{Y}_{x_{2}}\|\leq\varepsilon$.
Therefore, we have that

$\textstyle{\|\mathbf{Y}^{n}_{x}-\mathbf{Y}_{x}\|\leq\sum_{k=1}^{n}\varphi_{k}(x)\|\mathbf{Y}_{x_{k}}-\mathbf{Y}_{x}\|\leq\varepsilon.}$
(39)

Therefore, we have that
$\lim_{n\to+\infty}\sup_{x\in\mathcal{X}}\|\mathbf{Y}^{n}_{x}-\mathbf{Y}_{x}\|=0$.
Therefore, for any $n\in\mathbb{N}$, denoting $\mu_{n}$ the distribution of
$\mathbf{Y}^{n}$, we get that $(\mu_{n})_{n\in\mathbb{N}}$ converges to $\mu$
in $\mathscr{P}(\mathrm{C}(\mathcal{X},\mathbb{R}^{d}))$. For any
$n\in\mathbb{N}$, let
$\mathbf{a}=\\{a_{k}\\}_{k=1}^{n}\in(\mathbb{R}^{d})^{n}$ and
$f_{n}^{\mathbf{a}}:\ \mathcal{X}\to\mathbb{R}^{d}$ such that for any
$x\in\mathcal{X}$, $f_{n}^{\mathbf{a}}(x)=\sum_{k=1}^{n}a_{k}\varphi_{k}(x)$.
Denote
$\mathrm{C}_{n}(\mathcal{X},\mathbb{R}^{d})=\\{f_{n}^{\mathbf{a}}\,:\mathbf{a}\in(\mathbb{R}^{d})^{n}\\}$.
Define $\varphi_{n}:\
(\mathbb{R}^{d})^{n}\to\mathrm{C}_{n}(\mathcal{X},\mathbb{R}^{d})$ such that
for any $\mathbf{a}\in(\mathbb{R}^{d})^{n}$,
$\varphi_{n}(\mathbf{a})=f_{n}^{\mathbf{a}}$. We have that $\varphi_{n}$ is a
bijection. In addition, we have that
$(\varphi_{n})_{\\#}\mu_{n}=\mu_{x_{1},\dots,x_{n}}$. Therefore, using the
data processing inequality we have that for any $n\in\mathbb{N}$

$\operatorname{KL}(\mu_{n}|\nu_{n})=\operatorname{KL}(\mu_{x_{1},\dots,x_{n}}|\nu_{x_{1},\dots,x_{n}}).$
(40)

In addition, since $(\mu_{n},\nu_{n})\to(\mu,\nu)$ we using that the Kullback-
Leibler divergence is lower semi-continuous (Dupuis and Ellis,, 2011, Lemma
1.4.3)

$\textstyle{\lim_{n\to+\infty}\operatorname{KL}(\mu_{n}|\nu_{n})=\lim_{n\to+\infty}\operatorname{KL}(\mu_{x_{1},\dots,x_{n}}|\nu_{x_{1},\dots,x_{n}})\geq\operatorname{KL}(\mu|\nu).}$
(41)

Therefore, we get that

$\textstyle{\operatorname{KL}(\mu|\nu)\leq\sup\\{\operatorname{KL}(\mu_{x_{1},\dots,x_{n}}|\nu_{x_{1},\dots,x_{n}})\,:(x_{1},\dots,x_{n})\in\mathcal{X}^{n},\
n\in\mathbb{N}\\}.}$ (42)

Combining this result and Eq. 37, we conclude the proof. ∎

The rest of the proof is similar to Sec. E.1, expect that one needs to check
that the obtained Gaussian process $\pi^{\star}$ is sample continuous (up to a
modification). This is done using that there exists $\phi:\
\left[0,+\infty\right)$ such that for any $x_{1},x_{2}\in\mathcal{X}$,

$\textstyle{\mathbb{E}[\|\mathbf{Y}_{x_{1}}-\mathbf{Y}_{x_{2}}\|^{2}]\leq\phi(\|x_{1}-x_{2}\|),}$
(43)

such that $\int_{0}^{+\infty}\phi(\exp[-t^{2}])\mathrm{d}t<+\infty$, see
Fernique, (1975).

##  Appendix F Experimental details

In the following sections we provide details on our experimental procedures.

###  F.1 1D Datasets

#### Quadratic

A synthetic dataset constructed to exhibit clear bi-modality. Samples consist
of function evaluations of $f(x):=ax^{2}+\epsilon$ with
$a\sim\text{Unif}(\\{-1,1\\})$ and $\epsilon\sim\text{N}(0,1)$ on a uniformly
sampled grid of 100 evaluation points in $[-10,10]$. We allow unlimited
samples during training, validation and testing phases.

#### Melbourne

A real dataset recording the number of pedestrians on streets in Melbourne,
each sample being a period of 24 hours with readings taken hourly (giving a
grid of 24 evaluation points per sample). The dataset is sourced from
http://www.timeseriesclassification.com/description.php?Dataset=MelbournePedestrian.
We pre-processed the dataset by removing rows with unobserved values and
rescaling the dataset to have unit variance. The pre-processed dataset
contains 2319 samples and we split the dataset in the ratio of $[0.8,0.1,0.1]$
for training, validation and testing respectively.

#### Gridwatch

A real dataset recording the energy demand of the UK grid, each sample being a
period of 24 hours with readings taken every five minutes (giving a grid of
288 evaluation points per sample). The dataset is sourced from
https://www.gridwatch.templar.co.uk/download.php by selecting the ‘demand’
field. We performed a number of pre-processing steps:

  * •

Select days where the readings were taken exactly on the fifth minute,

  * •

Remove days where the the difference between two consecutive readings was
beyond the 99.5th percentile of differences,

  * •

Remove days where consecutive readings did not vary for half an hour or more,

  * •

Centered the dataset so each sample has zero mean,

  * •

Scaled the dataset so each sample has unit variance.

This resulted in 1013 samples which we split in the ratio of $[0.8,0.1,0.1]$
for training, validation and testing respectively.

###  F.2 Baselines

#### Gaussian Processes

We used a Matern-1/2 kernel with lengthscale and noise variance learnt by
maximising the marginal likelihood.

#### Neural Processes

We modified a PyTorch implementation of Neural Proccesses available online at
https://github.com/EmilienDupont/neural-processes. Our modifications adapted
the implementation to unconditional training/sampling by implementing Equation
7 of Garnelo et al., 2018b . We set the dimension of the representation of
context points and the dimension of the latent variable to be 512 (1024 for
Gridwatch dataset). The encoder and decoder were three-layered FCNs with 512
(1024 for Gridwatch dataset) units and ReLU activations, trained to 1000
epochs.

###  F.3 Functional MMD and two-sample test power

Our quantitative results are based off the functional MMD and kernel two-
sample hypothesis test of Wynne and Duncan, (2022).

The hypothesis test considered uses the null hypothesis that samples from the
model and dataset are from the same distribution, while the alternative states
that they are not. We know that the null is false as the distribution of
samples from the model will never truly match the data distribution.
Therefore, we expect a perfectly powerful test to always reject the null
hypothesis. By choosing a specific test with finite power (for example by
restricting the number of samples used when computing the test statistic), we
can compare the quality of samples from different models by computing the
power of the specific test on each of the sets of samples. If the test
exhibits a lower power on a set of samples, it indicates those samples are
harder to distinguish from the data distribution and thus can be considered of
higher quality.

The specific test we choose is the kernel two-sample test with the ID kernel
computed on 10 samples. We perform 1000 tests to estimate the power.
Confidence intervals are obtained by repeating the procedure on models trained
from different seeds.

###  F.4 Uncurated samples

We visualise samples from all models on the 1D datasets and from SP-SGM on
MNIST.

Dataset | SP-SGM | NP | GP
---|---|---|---
|  |  |
|  |  |
|  |  |
Table 1: Samples from trained models on 1D datasets; Quadratic [top],
Melbourne [middle] and Gridwatch [bottom]. Figure 4: Uncurated MNIST samples,
covariance kernel, $M=100$

###  F.5 Further ablation studies

$M=20$ | $M=100$ | $M=200$
---|---|---
|  |
|  |

figureEffect of increasing truncation order $M$. Increasing reconstruction
quality on MNIST digits [top], decreasing quality of fit in spectral space
[bottom].

|
---|---
|
|
|
|

figureEigenfunctions [left] and samples from SP-SGM [right] for RBF kernel
with lengthscale in [1,2,3,4,5] [top to bottom respectively]

Figure 5: Forward noising process of SP-SGM on MNIST digits using RBF kernel
with lengthscales in [1,2,3,4,5] [top to bottom respectively]. Rightmost
columns show the pixel-wise mean and standard deviation of the reference
measure respectively.

|
---|---
|
|
|
|

figureEigenfunctions [left] and samples from SP-SGM [right] for Matern-3/2
kernel with lengthscale in [1,2,3,4,5] [top to bottom respectively]

Figure 6: Forward noising process of SP-SGM on MNIST digits using Matern-3/2
kernel with lengthscales in [1,2,3,4,5] [top to bottom respectively].
Rightmost columns show the pixel-wise mean and standard deviation of the
reference measure respectively.


11institutetext: APC, Univ Paris Diderot, CNRS/IN2P3, CEA/Irfu, Obs de Paris,
Sorbonne Paris CitÃ, France 22institutetext: IFIC - Instituto de FÃsica
Corpuscular (CSIC - Universitat de ValÃncia) c/ CatedrÃ tico JosÃ BeltrÃ n, 2
E-46980 Paterna, Valencia, Spain 33institutetext: Institut Universitaire de
France, 75005 Paris, France 44institutetext: DÃpartment de physique nuclÃaire
et corpusculaire, UniversitÃ de GenÃve, GenÃve, Switzerland. 55institutetext:
RWTH Aachen University, III. Physikalisches Institut B, Aachen, Germany.
66institutetext: Dipartimento di Fisica "Aldo Pontremoli", UniversitÃ degli
Studi di Milano e INFN, Milan, Italy 77institutetext: Institut de Physique
Nucléaire d’Orsay (IPNO), CNRS-IN2P3, Orsay, France 88institutetext: Centro
Brasileiro de Pesquisas FÃsicas, Rio de Janeiro, RJ, Brazil 99institutetext:
Centro AtÃmico Bariloche, S.C de Bariloche, Argentina 1010institutetext:
Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba, Japan
1111institutetext: Service de Physique ThÃorique, UniversitÃ Libre de
Bruxelles, Ixelles, Belgium

# Search for a correlation between the UHECRs measured by the Pierre Auger
Observatory and the Telescope Array and the neutrino candidate events from
IceCube and ANTARES

J. Aublin, A. Coleiro, A. Kouchner for the ANTARES Collaboration
I. Al Samarai, A. Barbano, T. Montaruli L. Schumacher C. Wiebusch for the
IceCube Collaboration
L. Caccianiga, P.L. Ghia, U. Giaccari, G. Golup for the Pierre Auger
Collaboration
H. Sagawa, P. Tinyakov for the Telescope Array Collaboration
1111221133444444555566 lorenzo.caccianiga(at)unimi.it 77889910101111

###### Abstract

High-energy neutrinos are expected to be produced by the interaction of
accelerated particles near the acceleration sites. For this reason, it is
interesting to search for correlation in the arrival directions of ultra-high
energy cosmic rays (UHECRs) and HE neutrinos. We present here the results of a
search for correlations between UHECR events measured by the Pierre Auger
Observatory and Telescope Array and high-energy neutrino candidate events from
IceCube and ANTARES. We perform a cross-correlation analysis, where the
angular separation between the arrival directions of UHECRs and neutrinos is
scanned. When comparing the results with the expectations from a null
hypothesis contemplating an isotropic distribution of neutrinos or of UHECR we
obtain post-trial p-values of the order of $\sim 10^{-2}$.

##  1 Introduction

The sources of ultra-high energy cosmic rays (UHECRs) are still unknown,
mostly because of the difficulties in tracking their origin after they are
deflected by magnetic fields. On the other hand, neutrinos, being neutral, do
not suffer deflections and carry directional information on their origin.
High-energy neutrinos are expected to be produced at the same sites of UHECRs,
when the latter interact with the source itself or in its vicinity. It would
be then natural to search for a correlation in the arrival directions of
UHECRs and neutrinos. Previous searches performed by the Pierre Auger
Observatory, Telescope Array and IceCube collaborations showed an interesting
hint of correlation at a $3\sigma$ level jcap  which was later reduced with
additional data ICRC2017 . In these proceedings we show the results of an
update of this analysis, with new data from the original collaborations and
with the addition of data from the ANTARES neutrino telescope.

##  2 The UHECR observatories and datasets

Ultra-high energy cosmic rays are observed via the extensive air showers (EAS)
they induce in the atmosphere. The EAS can be observed in their longitudinal
development through telescopes that detect the fluorescence emission of the
atmosphere particles or sampled at the ground with arrays of particle
detectors. The first measurement gives a direct calorimetric estimate of the
energy of the primary particle but can be performed only during clear moonless
nights (low duty cycle). The hybrid concept that both the Pierre Auger
Observatory and Telescope Array exploit relies on a cross-calibration of the
surface detector with the events observed through fluorescence. In this way,
an estimate of the energy of the primary particle can be performed with
surface data only ($\sim 100\%$ duty cycle) without relying on simulations. A
summary of UHECR data used in this and previous works can be found in Table 1.

###  2.1 The Pierre Auger Observatory

The Pierre Auger Observatory augernim ; NIM2015 , located in Argentina at a
latitude of $\sim-35^{\circ}$, is a $3000\;\text{km}^{2}$ cosmic-ray detector.
It is a hybrid detector, using both a surface detector array (SD) composed of
1660 water-Cherenkov detectors and an atmospheric fluorescence detector (FD)
made of 27 fluorescence telescopes to observe extensive air showers induced in
the atmosphere by primary particles with energies
$E>0.1\;\text{EeV}$111$1\;\text{EeV}=10^{18}\;\text{eV}$. For this study, 324
events with $E>52\;$EeV recorded by SD from January 2004 to April 2017 were
used, resulting in 90 events more (corresponding to three years) with respect
to previous analyses ICRC2017 . Events are observed up to a zenith angle of
$80^{\circ}$, which translates into a field of view ranging from $-90^{\circ}$
to $+45^{\circ}$ in declination. The systematic uncertainty on the energy
scale is 14%, and the statistical uncertainty in the energy is smaller than
12% for these events VerziICRC2013 , while the angular uncertainty is less
than $0.9^{\circ}$ at these energies.

###  2.2 Telescope Array

Telescope Array (TA) TASDNIM ; TAFDNIM  is a UHECR detector located in Utah,
USA, at a latitude of $+39^{\circ}$. It is composed of 507 plastic
scintillators distributed on a square grid covering a $700\;\text{km}^{2}$
surface, overlooked by three fluorescence detector stations housing 38
telescopes. The events used for this analysis have an energy $E\geq 57\>$EeV
and zenith angles up to $55^{\circ}$. A total of 143 events recorded from May
2008 to May 2017 is used in this work, that corresponds to 34 more (two years)
than the ones used in the latest update ICRC2017 . These events have about
$1.5^{\circ}$ angular resolution, $\sim 20\%$ energy resolution is and a $\sim
22\%$ systematic uncertainty on the energy scaleTAspec .

##  3 The neutrino observatories and datasets

High-energy astrophysical neutrinos are observed by instrumenting large
volumes ($\sim 0.01-1\;\text{km}^{3}$) of water or ice with photomultipliers
(PMT) in order to use them as Cherenkov detectors. The shape of a neutrino
event depends on the neutrino flavor and on the type of interaction it
undergoes. Neutrinos of all flavors via neutral current (NC) and $\nu_{e}$ via
charged current (CC) produce electromagnetic showers that are usually well-
contained within the detector, offering a good energy estimate but a poor
angular resolution (shower-like events, or cascades). Vice-versa, muons
produced in $\nu_{\mu}$ CC interactions tend to cross through the whole
detector resulting in a good angular resolution and a poor energy estimate
(track-like events, or tracks). A summary of neutrino candidate events used in
this and previous works can be found in Table 1.

###  3.1 The ANTARES neutrino telescope

The ANTARES neutrino telescope AntNIM  is located at $2475\;$m of depth in the
Mediterranean sea, off the coast of Toulon, France (latitude: $+43^{\circ}$.
It is composed of 12 $450\>$m-long strings, each equipped with 75 optical
modules, at a depth of about $2500\>$m. The events used for this work are
selected from the diffuse nine-year sample AntDiffuse  by requiring a
âsignalnessâ $>40\%$. The signalness is defined as the ratio of the number of
expected astrophysical events over the sum of the expected atmospheric and
astrophysical events for a given energy proxy. The selection resulted in two
track-like events and no cascade events. The same signalness cut applied to
the point-source sample AntPS  resulted in an additional track event selected,
bringing the ANTARES dataset used in this analysis to a total of three tracks
and no cascades. For these events, the median angular resolution is $\sim
0.4^{\circ}$.

###  3.2 The IceCube South Pole Neutrino Observatory

IceCube icecubeA ; icecube  is a high-energy neutrino detector of $\sim
1\;\text{km}^{3}$ composed by 86 strings instrumented by 5160 PMTs installed
deep ($1.5-2.5\;$km) in the ice at the geographic South Pole (latitude:
$-90^{\circ}$). The dataset used for this analysis is composed of events from
two different samples. The main one is the six-year High-Energy Starting
Events (HESE) dataset IC1ICRC2017 , which includes events generated by
neutrinos of all flavors interacting inside the detection volume (starting
events) with deposited energies ranging from 60 TeV up to 2 PeV. It includes
58 cascades and 15 tracks. The track sample is integrated with data from the
dataset of through-going muons induced by charged current interactions of
candidate $\nu_{\mu}$ from the Northern sky (diffuse through up-going tracks )
IceCubeMu . This dataset adds 35 tracks corresponding to seven years of data
from the eight-year data sample published in IC2ICRC2017 . Only events with
âsignalnessâ $>50\%$ are used (as described in Section 3.1), corresponding to
a threshold of $\sim 200\>$TeV on the muon energy proxy.

Cascade events in the HESE sample have an angular resolution of $\sim
15^{\circ}$ above 100 TeV. Tracks are reconstructed with an angular resolution
of $\lesssim 1^{\circ}$. The resolution of the deposited energy for tracks and
cascades is around 15% IceCubeEnergy , but cascades have a better resolution
for the reconstructed neutrino energy since they deposit most of their energy
inside the detector, while tracks do not.

##  4 Analysis methods

The search for correlations between UHECR and neutrino arrival directions is
performed using different methodsjcap : a cross-correlation analysis, a
likelihood analysis using neutrinos as stacked sources and a likelihood
analysis using UHECRs as stacked sources to search for correlation with point-
source neutrino samples (not described in paragraph 3). A recent update of the
third analysis with an improved method is described in Schumacher2018 . In
these proceedings an update of the cross-correlation analysis is presented.

The cross-correlation analysis is a simple yet completely model-independent
analysis. It does not require assumptions on the deflection of UHECR in
magnetic fields. This means that a part from the global energy cut chosen for
the dataset by each collaboration, no energy information is needed on the
single UHECR event. For this reason, no re-scaling of the UHECR energy scales
to match the TA and Auger flux in the common field of view will be applied in
the following, as it was done in previous works. Information on the flux
matching of the two UHECR datasets can, however, be found in other joint works
presented at this conference, in particular Ivanov:2018  and Biteau:2018 .

The principle of the cross-correlation analysis is to count the number
$n_{OBS}$ of UHECR-neutrino pairs separated by less than a certain angular
distance $\psi$. This number is then compared to the expected number of pairs
$n_{EXP}$ separated by less than $\psi$ in the null-hypothesis scenario. For
this analysis, two different null-hypotheses were investigated. The first one
is an isotropic distribution of UHECR, obtained by generating isotropic
datasets according to the UHECR exposure of the two experiments. The second
null-hypothesis is the complementary one: an isotropic distribution of
neutrinos, obtained by scrambling the real neutrino datasets, generating
randomly right ascension values but keeping the original declination, thus
preserving the declination-dependent acceptance of the neutrino observatories.
The two null-hypotheses are investigated separately and thus when neutrino
positions are randomized, UHECR ones are kept fixed and vice-versa. A large
number of samples satisfying the null hypothesis are generated, and the
fraction of generated samples where $n_{EXP}\geq n_{OBS}$ is observed is
considered the local p-value. The analysis was performed at various $\psi$
values, ranging from $1^{\circ}$ to $30^{\circ}$ in $1^{\circ}$ steps in order
to find the maximum departure from isotropy and the minimum local p-value.
This p-value is then penalized taking into account the scan on $\psi$ and
resulting in the final global p-value that will be reported in the following.
The analysis was performed independently on the cascade and track neutrino
datasets, thus obtaining four different p-values, one for each neutrino sample
and each null hypothesis.

##  5 Results

###  5.1 Previous results

In jcap , the first results of the analysis described in section 4 were
reported. No significant result was found using neutrino tracks, while with
cascades the penalized p-values were $5\times 10^{-4}$ with respect to an
isotropic flux of UHECR and $8.5\times 10^{-3}$ with respect to an isotropic
flux of neutrinos. In ICRC2017 , with additional neutrino data from IceCube
and one year more of UHECR events from TA, the correlation with tracks
continued to be not significant, while with cascades the p-values went up
roughly one order of magnitude, respectively to $5.4\times 10^{-3}$ and
$1\times 10^{-2}$.

###  5.2 Updated results

We performed the analysis on the updated datasets as described in paragraph 2
and 3. The results of the scan in angle are shown in figure 1, where the
relative excess of pairs with respect to the expected value from an isotropic
distribution of UHECRs is shown for tracks (left) and cascades (right). In
these conditions, the maximum departure from isotropy is found at $1^{\circ}$
for tracks, where four pairs are observed and at $16^{\circ}$ for tracks, with
observed 623 pairs. The post-trial p-values are 0.45 for tracks and $2.7\times
10^{-2}$ for cascades. The most significant angles are found to be the same
for the null hypothesis of an isotropic distribution of neutrinos. The
corresponding penalized p-values are 0.49 and $2.6\times 10^{-2}$ for tracks
and cascades respectively. A summary of the updated results compared to the
previous ones can be found in table 2.

##  6 Discussion

In jcap , a possible hint of a correlation between the arrival directions of
UHECR and HE neutrinos was reported, using a simple and model-independent
method. A first update published in ICRC2017  lowered the significance of the
finding. This work uses new UHECR data from the Pierre Auger Observatory and
Telescope Array experiment and new neutrino tracks from the ANTARES neutrino
telescope. With these updated datasets, the post-trial p-values are now higher
than 0.01 both considering as a null hypothesis an isotropic distribution of
UHECRs or neutrinos.

This result can be explained by a number of considerations. First of all, the
composition of ultra-high energy cosmic rays is not well known at these
energies. This, together with the poor knowledge of the Galactic and extra-
galactic magnetic field, implies that we cannot distinguish between UHECRs
that suffer small or large deflections. Secondly, as already pointed out in
jcap , the propagation of UHECRs is limited by the interactions with the
cosmic backgrounds (GZK effects and similar): depending on the composition, it
is expected that cosmic rays at these energies cannot arrive from further than
10-100 Mpc. On the other hand, neutrinos can reach us from cosmological
distances so that in principle only a few percent of neutrinos would be
expected from the sources within this horizon and thus in correlation with the
arrival directions of the detected UHECRs. Finally, neutrinos originated at
cosmic rays acceleration sites are expected to carry few percents of the
energy of the original cosmic-ray. Thus, the neutrinos observed by IceCube and
ANTARES would have been produced by cosmic rays of much lower energy than the
ones in the datasets by Auger and TA. For this reason, it is possible that
only a fraction of the neutrinos in our data could come from a site that is
actually capable of accelerating cosmic rays up to $\sim 10^{20}\;$eV. We
expect to carry on this work by soon updating the other two analyses with the
latest data.

Figure 1: Results of the cross-correlation analysis applied to the neutrino tracks (left) and cascades (right) samples. The relative excess of pairs is shown as a function of the maximum separation angle for data (dots) and expectation from an isotropic distribution of UHECR (colored area, with one, two and three-$\sigma$ dispersion).  Table 1: Summary of the data used in this work compared to previous ones. | JCAP 2016 jcap  | ICRC 2017 ICRC2017  |  | UHECR 2018
---
(This work)
Auger UHECR |  | 231 events
---
1/1/2004 - 31/3/2014
| 324 events
---
1/1/2004 - 30/4/2017
TA UHECR |  | 87 events
---
11/5/2008 - 1/5/2014
| 109 events
---
11/5/2008 - 1/5/2015
| 143 events
---
11/5/2008 - 1/5/2017
IceCube HE $\nu$ |  | 39 cascades
---
16 tracks
(4 yrs HESE + diffuse up-going)
| 58 cascades
---
49 tracks
(6 yrs HESE + diffuse up-going)
ANTARES HE $\nu$ | - |  | 3 tracks
---
(9-years diffuse + point source )
Table 2: Summary of the post-trial p-values obtained in this work compared to previous ones. | JCAP 2016 jcap  | ICRC 2017 ICRC2017  |  | UHECR 2018
---
(This work)
tracks wrt an isotropic flux of UHECR | 0.28 | 0.48 | 0.45
tracks wrt an isotropic flux of neutrinos |  | 0.52 | 0.49
cascades wrt an isotropic flux of UHECR | $5\times 10^{-4}$ | $5.4\times 10^{-3}$ | $\boldsymbol{2.7\times 10^{-2}}$
cascades wrt an isotropic flux of neutrinos $\nu$ | $8.5\times 10^{-3}$ | $1.0\times 10^{-2}$ | $\boldsymbol{2.6\times 10^{-2}}$

## References

  * (1) M.G. Aartsen et al. (IceCube, Pierre Auger, Telescope Array), JCAP 01(2016)037 (2016)
  * (2) I. Al Samarai for the IceCube, Telescope Array and Pierre Auger Collaborations, PoS(ICRC2017)961 (2017)
  * (3) The Pierre Auger Collaboration, Nucl. Instrum. Meth. A 523, 50 (2004)
  * (4) The Pierre Auger Collaboration, Nucl. Instrum. Meth. A 798, 172 (2015)
  * (5) V. Verzi, the Pierre Auger Collaboration, Proceedings of the 33rd ICRC, Rio de Janeiro, Brasil (2013)
  * (6) T. Abu-Zayyad et al. (Telescope Array), Nucl. Instrum. Meth. A689, 87 (2013), 1201.4964
  * (7) H. Tokuno et al., Nucl. Instrum. Meth. A676, 54 (2012), 1201.0002
  * (8) T. Abu-Zayyad, R. Aida, M. Allen, R. Anderson, R. Azuma, E. Barcikowski, J.W. Belz, D.R. Bergman, S.A. Blake, R. Cady et al., The Astrophysical Journal Letters 768, L1 (2013)
  * (9) M. Ageron et al. (ANTARES), Nucl. Instrum. Meth. A656, 11 (2011), 1104.1607
  * (10) A. Albert et al. (ANTARES), Astrophys. J. 853, L7 (2018), 1711.07212
  * (11) A. Albert et al. (ANTARES), Phys. Rev. D96, 082001 (2017), 1706.01857
  * (12) A. Achterberg et al. (IceCube), Astropart. Phys. 26, 155 (2006), astro-ph/0604450
  * (13) T. Deyoung, Advances in Space Research 51, 238 (2013)
  * (14) C. Kopper for the IceCube Collaboration, PoS(ICRC2017)981 (2017)
  * (15) M.G. Aartsen et al. (IceCube), Astrophys. J. 833, 3 (2016), 1607.08006
  * (16) C. Haack, C. Wiebusch on behalf of the IceCube Collaboration, PoS(ICRC2017)1005 (2017)
  * (17) M.G. Aartsen et al. (IceCube), JINST 9, P03009 (2014), 1311.4767
  * (18) L.Schumacher for the ANTARES, IceCube, Pierre Auger and Telescope Array Collaborations, Proceedings of VLVnT2018, EPJ Web of Conferences (in press) (2019)
  * (19) A. Ivanov for the Pierre Auger and Telescope Array Collaborations, Proceedings of UHECR2018, to appear in EPJ Web of Conferences (2019)
  * (20) J. Biteau for the Pierre Auger and Telescope Array Collaborations, Proceedings of UHECR2018, to appear in EPJ Web of Conferences (2019)