# Adaptive RKHS Fourier Features for Compositional Gaussian Process Models

Xinxing Shi
Department of Computer Science
The University of Manchester
<EMAIL_ADDRESS>
&Thomas Baldwin-McDonald
Department of Computer Science
The University of Manchester
<EMAIL_ADDRESS>&Mauricio A. Álvarez
Department of Computer Science
The University of Manchester
<EMAIL_ADDRESS>

###### Abstract

Deep Gaussian Processes (DGPs) leverage a compositional structure to model
non-stationary processes. DGPs typically rely on local inducing point
approximations across intermediate GP layers. Recent advances in DGP inference
have shown that incorporating global Fourier features from Reproducing Kernel
Hilbert Space (RKHS) can enhance the DGPs’ capability to capture complex non-
stationary patterns. This paper extends the use of these features to
compositional GPs involving linear transformations. In particular, we
introduce Ordinary Differential Equation (ODE) -based RKHS Fourier features
that allow for adaptive amplitude and phase modulation through convolution
operations. This convolutional formulation relates our work to recently
proposed deep latent force models, a multi-layer structure designed for
modelling nonlinear dynamical systems. By embedding these adjustable RKHS
Fourier features within a doubly stochastic variational inference framework,
our model exhibits improved predictive performance across various regression
tasks.

##  1 Introduction

Gaussian Processes (GPs) provide a principled Bayesian framework for function
approximation, making them particularly useful in many applications requiring
uncertainty calibration (Rasmussen and Williams, 2006), such as Bayesian
optimisation (Snoek et al., 2012) and time-series analysis (Roberts et al.,
2013). Despite offering reasonable uncertainty estimation, shallow GPs often
struggle to model complex, non-stationary processes present in practical
applications. To overcome this limitation, Deep Gaussian Processes (DGPs)
employ a compositional architecture by stacking multiple GP layers, thereby
enhancing representational power while preserving the model’s intrinsic
capability to quantify uncertainty (Damianou and Lawrence, 2013). However, the
conventional variational formulation of DGPs heavily depends on local inducing
point approximations across intermediate GP layers (Titsias, 2009; Salimbeni
and Deisenroth, 2017), which hinder the model from capturing the global
structures commonly found in real-world scenarios.

Incorporating _Fourier features_ into GP models has shown promise in
addressing this challenge in GP inference due to the periodic nature of these
features. A line of research uses Random Fourier Features (RFFs, (Rahimi and
Recht, 2007)) of stationary kernels to convert the original (deep) GPs into
Bayesian networks in weight space (Lázaro-Gredilla et al., 2010; Gal and
Turner, 2015; Cutajar et al., 2017). Building on this concept within a sparse
variational GP framework, recent advancements in inter-domain GPs (Lázaro-
Gredilla and Figueiras-Vidal, 2009a; Van der Wilk et al., 2020) directly
approximate the posterior of the original GPs by introducing fixed
_Variational Fourier Features_ (VFFs) through process projection onto a
Reproducing Kernel Hilbert Space (RKHS)(Hensman et al., 2018; Rudner et al.,
2020).

VFFs are derived by projecting GPs onto a different domain. The original GP
posterior that these VFFs attempt to approximate remains within the same
functional space as the original GP. In this setting, the VFFs produce a set
of static basis functions determined by a fixed set of frequency values. To
enhance these features and introduce greater flexibility, we propose a
generalised approach that incorporates Fourier features into inter-domain GPs
through linear transformations, such as convolution operations.

In this paper, we focus on a type of GP characterised as the output of a
convolution operation between a smoothing kernel and a latent GP. An example
of this construction is the Latent Force Model (LFM, (Alvarez et al., 2009)),
in which the smoothing kernel corresponds to the Green’s function associated
with an Ordinary Differential Equation (ODE). By incorporating RKHS Fourier
features into this framework, we derive adaptive global features inspired by
the ODE, allowing for the optimisation of amplitudes and phases. We name the
obtained features _Variational Fourier Response Features_ (VFRFs) since they
are derived from the output of a linear system. To enhance the capability of
our model, we further use these adaptive features in a compositional GP model
that stacks multiple LFMs, also known as Deep LFM (DLFM, (McDonald and
Álvarez, 2021)). This hierarchical structure facilitates more precise and
robust modelling of complex, non-stationary data. Our experimental results on
both synthetic and real-world data demonstrate that incorporating these ODE-
inspired RKHS Fourier features improves upon the standard practice of using
VFFs.

##  2 Background

This section reviews concepts and preliminaries relevant to this work,
including GPs, LFMs (Section 2.1), and the sparse variational framework with
Fourier features for GP inference (Section 2.2). We also establish the
notation used throughout the subsequent discussions.

###  2.1 Latent Force Models

##### Gaussian Processes

A Gaussian Process (GP)
$f(\cdot)\sim\mathcal{GP}(m(\cdot),k(\cdot,\cdot^{\prime}))$ places
probability measures on a function space
$\\{f:\mathbb{R}^{D}\rightarrow\mathbb{R}\\}$ (Rasmussen and Williams, 2006).
Its behaviour is characterised by the mean function
$m:\mathbb{R}^{D}\rightarrow\mathbb{R}$ and the covariance function
$k:\mathbb{R}^{D}\times\mathbb{R}^{D}\rightarrow\mathbb{R}$. The evaluation of
the function $f(\cdot)$ at an input of interest $\mathbf{x}$ is a random
variable denoted as $f(\mathbf{x})\in\mathbb{R}$. Given a dataset of inputs
$\mathbf{X}=[\mathbf{x}_{n}]_{n=1}^{N}\in\mathbb{R}^{N\times D}$ and
corresponding measurements $\mathbf{y}=[y_{n}]_{n=1}^{N}\in\mathbb{R}^{N}$, we
assume $\mathbf{y}$ is observed from a noise-corrupted GP:
$y_{n}=f(\mathbf{x}_{n})+\epsilon,\epsilon\sim\mathcal{N}(\epsilon\mid
0,\varepsilon^{2})$, where $\varepsilon^{2}$ is the noise variance. Then the
posterior distribution $p(f\mid\mathbf{y})$ at any input locations
$\mathbf{x},\mathbf{x}^{\prime}$ remains a GP
$f\mid\mathbf{y}\sim\mathcal{GP}(\mu(\mathbf{x}),\Sigma(\mathbf{x},\mathbf{x}^{\prime}))$,
with updated mean and covariance given by

$\displaystyle\mu(\mathbf{x})$
$\displaystyle=m(\mathbf{x})+k(\mathbf{x},\mathbf{X})(\mathbf{K}_{\mathbf{XX}}+\varepsilon^{2}\mathbb{I}_{N})^{-1}(\mathbf{y}-m(\mathbf{X})),$
(1) $\displaystyle\Sigma(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime})-k(\mathbf{x},\mathbf{X})(\mathbf{K}_{\mathbf{XX}}+\varepsilon^{2}\mathbb{I}_{N})^{-1}k(\mathbf{X},\mathbf{x}^{\prime}),$

where $[m(\mathbf{X})]_{n}=m(\mathbf{x}_{n})$,
$[\mathbf{K_{XX}}]_{i,j}=k(\mathbf{x}_{i},\mathbf{x}_{j})$. The kernel
parameters and the noise variance are typically optimised by maximising the
log-likelihood $\log p(\mathbf{y}\mid\mathbf{X})$.

##### Latent Force Model

A Latent Force Model (LFM, (Alvarez et al., 2009)) integrates differential
equations with GPs to probabilistically model dynamic physical systems. Based
on prior physical knowledge, a single-output LFM assumes the system’s output
$f(t)$ is influenced by $Q$ _latent forces_ $\\{u_{q}(t)\\}_{q=1}^{Q}$ through
differential equations. Commonly, a first-order LFM uses the following form of
ODE (Guarnizo and Álvarez, 2018)

$\frac{\mathop{}\\!\mathrm{d}f(t)}{\mathop{}\\!\mathrm{d}t}+\gamma
f(t)=\sum_{q=1}^{Q}S_{q}u_{q}(t),$ (2)

where $\gamma>0$ is a decay parameter, and $S_{q}\in\mathbb{R}$ is a
sensitivity parameter. The solution for the output process $f(t)$ takes the
form of weighted convolution integrals
$f(t)=\sum_{q=1}^{Q}S_{q}\int_{0}^{t}G(t-\tau)u_{q}(\tau)\mathop{}\\!\mathrm{d}\tau$,
where $G(\cdot)$ denotes the Green’s function associated with the ODE.

Latent forces are presumed to follow GP priors,
$u_{q}(t)\sim\mathcal{GP}(0,k_{q}(t,t^{\prime}))$, leading to a covariance
function for the outputs
$k_{f}(t,t^{\prime})=\sum_{q=1}^{Q}S^{2}_{q}\int_{0}^{t}\int_{0}^{t^{\prime}}G(t-\tau)k_{q}(\tau,\tau^{\prime})G(t^{\prime}-\tau^{\prime})\mathop{}\\!\mathrm{d}\tau\mathop{}\\!\mathrm{d}\tau^{\prime}.$
For some types of covariance functions $k_{q}(t,t^{\prime})$, e.g., the radial
basis function (RBF), the covariance $k_{f}(t,t^{\prime})$ can either be
computed explicitly (Lawrence et al., 2006) or approximated by the inner
product using convolved Random Fourier Features (RFFs) (Rahimi and Recht,
2007)

$\boldsymbol{\phi}_{q}(t)=\frac{S_{q}}{\sqrt{M}}\left[\int_{0}^{t}G(t-\tau)e^{j\omega_{m}\tau}\mathop{}\\!\mathrm{d}\tau\right]_{m=1}^{M},\quad\omega_{m}\overset{\mathrm{i.i.d.}}{\sim}p_{q}(\omega),$
(3)

where $j^{2}=-1$ and $\omega_{m}$ is sampled i.i.d. from the spectral density
$p_{q}(\omega)$ associated with $k_{q}$ (Guarnizo and Álvarez, 2018). By
plugging the physics-informed kernels into the GP posterior (1), LFMs embed
domain-specific knowledge into the learning process and can utilise the sparse
approximation techniques in GP inference.

###  2.2 Variational Fourier Features

##### Sparse Variational GPs

Exact GPs, as described in (1), suffer from $\mathcal{O}(N^{3})$ time
complexity and are limited to Gaussian likelihoods. Sparse Variational
Gaussian Processes (SVGPs, (Titsias, 2009; Hensman et al., 2013, 2015))
provide a scalable inference framework by introducing a small set of $M(\ll
N)$ inducing points $\mathbf{Z}=[\mathbf{z}_{m}]_{m}^{M}\in\mathbb{R}^{M\times
D}$ and corresponding inducing variables
$\mathbf{u}=[u(\mathbf{z}_{m})]_{m=1}^{M}\in\mathbb{R}^{M}$ from the GP prior,
i.e., $p(\mathbf{u})=\mathcal{N}(\mathbf{u}\mid
m(\mathbf{Z}),\mathbf{K_{ZZ}})$. A variational distribution
$q(\mathbf{u})=\mathcal{N}(\mathbf{u}\mid\mathbf{m},\mathbf{S})$ is employed
to approximate the posterior process $q(f(\mathbf{x}))=\int
p(f(\mathbf{x})\mid\mathbf{u})q(\mathbf{u})\mathop{}\\!\mathrm{d}\mathbf{u}=\mathcal{N}(f\mid\tilde{m},\tilde{\Sigma})$,
where

$\displaystyle\tilde{\mu}(\mathbf{x})$
$\displaystyle=m(\mathbf{x})+k(\mathbf{x},\mathbf{Z})\mathbf{K}^{-1}_{\mathbf{ZZ}}(\mathbf{m}-m(\mathbf{Z})),$
(4) $\displaystyle\tilde{\Sigma}(\mathbf{x},\mathbf{x}^{\prime})$
$\displaystyle=k(\mathbf{x},\mathbf{x}^{\prime})+k(\mathbf{x},\mathbf{Z})\mathbf{K}^{-1}_{\mathbf{ZZ}}(\mathbf{S}-\mathbf{K_{ZZ}})\mathbf{K}^{-1}_{\mathbf{ZZ}}k(\mathbf{Z},\mathbf{x}^{\prime}).$

SVGPs learn the optimal placement of the inducing points and the variational
distribution by maximising an Evidence Lower BOund (ELBO) of $\log
p(\mathbf{y}\mid\mathbf{X})$.

##### Variational Fourier Features

Inter-domain GPs (Lázaro-Gredilla and Figueiras-Vidal, 2009b; Álvarez et al.,
2010; Van der Wilk et al., 2020) extend the domain of inducing variables by
integrating the GP $f$ with a deterministic inducing function $g$:

$u(\mathbf{z})=\int_{\mathbb{R}^{D}}g(\mathbf{x},\mathbf{z})f(\mathbf{x})\mathop{}\\!\mathrm{d}\mathbf{x},\quad\mathbf{z}\in\mathbb{R}^{D^{\prime}}.$
(5)

This formulation allows for a redefinition of inducing variables
$\mathbf{u}=u(\mathbf{Z})$, which still share the GP prior, albeit with
alternative expressions of $k(\mathbf{x},\mathbf{Z})$ and $\mathbf{K_{ZZ}}$
used in (4). By choosing various functions for $g$, inter-domain GPs
facilitate the construction of vector-valued basis functions
$k(\cdot,\mathbf{Z})$ for more informative feature extraction while
maintaining the standard SVGP framework.

VFFs (Hensman et al., 2018) define each inter-domain inducing variable $u_{m}$
of $\mathbf{u}$ by projecting the original GP $f$ onto a Fourier basis:
$u_{m}=\langle\phi_{m},f\rangle_{\mathcal{H}}$, where
$\langle\cdot,\cdot\rangle_{\mathcal{H}}$ denotes the Matérn RKHS inner
product on an interval $[a,b]$, and $\phi_{m}$ is the $m$-th entry of a
truncated Fourier basis

$\boldsymbol{\phi}(x)=[1,\cos(z_{1}(x-a)),\cdots,\cos(z_{M}(x-a)),\sin(z_{1}(x-a)),\cdots,\sin(z_{M}(x-a))],$
(6)

with $x$ a scalar input. In this setting, $\mathbf{z}=[z_{m}]_{m=1}^{M}$ are
$M$ _inducing frequencies_ , analogous to inducing points in the SVGP context
111In Hensman et al. (2018), inducing frequencies are denoted as $\omega_{m}$.
In this work, we use $z_{m}$ instead to maintain a notation consistent with
the inducing points in SVGPs and inter-domain GPs.. The projection of $f$ onto
this basis results in sinusoidal terms in the cross-covariance, i.e.,
$\mathrm{Cov}[f(\cdot),u_{m}]=\phi_{m}(\cdot)$ due to the reproducing property
of the RKHS.

##  3 Methodology

This section describes the method of integrating RKHS Fourier features into
compositional inter-domain GPs, with a specific focus on LFMs within SVGP
framework. We start by adapting the conventional ODE used in LFMs to
incorporate the inter-domain GPs with VFFs as a special instance of our model
(Section 3.1). Details on VFRFs are provided in Section 3.2. We then extend
our model from a single-layer to a hierarchical structure in Section 3.3.

###  3.1 LFMs with Modified ODEs

In this work, we focus on a dynamical system modelled by a potential first-
order ODE without loss of generality

$\beta\frac{\mathop{}\\!\mathrm{d}f(t)}{\mathop{}\\!\mathrm{d}t}+\alpha
f(t)=u(t),$ (7)

where $\alpha,\beta$ are positive coefficients and
$u(t)\sim\mathcal{GP}(0,k(t,t^{\prime}))$ represents an unknown latent force
with a Matérn kernel with half-integer order. The Green’s function of (7) is
$G(t)=\frac{1}{\beta}\exp(-\frac{\alpha}{\beta}t)=\frac{1}{\beta}\exp(-\gamma
t)$. We introduce $\gamma=\frac{\alpha}{\beta}$ to remain consistent with the
decay parameter in (2). Unlike the conventional formulation (2), which
involves a weighted sum of multiple latent forces, the modified ODE (7)
simplifies it to a single process $u(t)$. It can be trivially decomposed into
distinct latent forces if necessary. Moreover, we will further show that, by
introducing coefficients $\alpha$ and $\beta$ in our model, the output process
$f$ can revert to a GP with VFFs as $\beta\rightarrow 0^{+}$. This formulation
enables practitioners to apply the proposed approach in scenarios where prior
knowledge of the system is limited and there is no prior knowledge indicating
if the dynamics encoded in the kernel accurately reflect the observed data.

A solution $f(t)$ can be expressed as a convolution integral

$f(t)=\int_{-\infty}^{t}G(t-\tau)u(\tau)\mathop{}\\!\mathrm{d}\tau=G\circ u,$
(8)

where the integral’s lower limit is extended to negative infinity to maintain
variance near the origin, though it can be adjusted based on data range or
prior knowledge in practice. The covariance function of the output process
$f(\cdot)$ is derived by applying the convolution operator to the kernel $k$’s
arguments, respectively:

$\text{Cov}[f(t),f(t^{\prime})]=\int_{-\infty}^{t^{\prime}}\int_{-\infty}^{t}G(t-\tau)k(\tau,\tau^{\prime})G(t^{\prime}-\tau^{\prime})\mathop{}\\!\mathrm{d}\tau\mathop{}\\!\mathrm{d}\tau^{\prime}.$
(9)

The covariance function (also LFM kernel) can be calculated analytically if
$k$ is a Matérn kernel with half-integer orders. We give the closed-form
covariances in Table 1 of Appx. 9.

##### Model interpretation

The Green’s function $G(\cdot)$, determined by the system’s dynamics, serves
as a signal filter. It effectively acts as a low-pass filter described by the
ODE (7), with $\gamma$ representing the “cutoff frequency”. Mathematically,
the convolution operator $G(\cdot)$ of the modified ODE will behave like the
Dirac delta function in (8) as $\alpha=1,\beta\rightarrow 0^{+}$
(i.e.,$\gamma\rightarrow+\infty$), causing $f(t)$ to closely replicate $u(t)$.
Fig. 1, left, illustrates this behaviour. We use the Matérn-$\tfrac{1}{2}$
kernel for the covariance $k(t,t^{\prime})$ of the latent force. The figure
shows this covariance function (dashed) and two LFMs covariance functions
(solid) $k_{f}(t,t^{\prime})$ with different $\gamma$ values. The LFM kernel
reverts to the latent force kernel as $\gamma$ increases. However, the
conventional LFM kernel without the ODE modification, i.e., (2) (Guarnizo and
Álvarez, 2018) will get flattened since the corresponding Green’s function
$\exp{(-\gamma t)}$ does not effectively mimic a valid Dirac delta function.

Figure 1: Covariance functions of LFMs (left) and Variational Fourier Response
Features (VFRFs) (right). The latent force $u(t)$ uses a Matérn-$\tfrac{1}{2}$
kernel with length-scale $l=0.2$ (left dashed). Left: The centred kernel of
the input latent force (dashed) and the output process $f(t)$ of LFMs with
different ODE parameters $\gamma$ (solid). Unlike the LFM kernel induced by
(9) (green), the modified LFM kernel from (7) can revert to the input Matérn-
$\tfrac{1}{2}$ kernel if increasing $\gamma$ (red to brown). Right: VFRFs
($G\circ\phi$, red solid) and VFFs ($\phi$, blue dashed) with different
inducing frequencies: $z_{m}=\tfrac{8\pi}{b-a}$ (upper) and
$\tfrac{28\pi}{b-a}$ (lower). The upper panel depicts the cosine basis with a
phase delay $\theta\approx\tfrac{\pi}{4}$ to the VFF, while the lower panel
displays the sine basis with a phase delay $\theta\approx\frac{5\pi}{12}$.

###  3.2 Variational Fourier Response Features for LFMs

Building upon the modified ODE described by (7), we introduce a spectral
approximation for the LFMs within the inter-domain GP framework. The latent
force $u$ is initially projected onto the Fourier basis entries $\phi_{m}$ as
defined in (6), yielding its spectral representations

$v_{m}=\langle\phi_{m},u\rangle_{\mathcal{H}},\quad m=0,1,\ldots,2M.$ (10)

The projected inducing variables $v_{m}$ are collected as
$\mathbf{v}=[v_{m}]_{m=0}^{2M}\in\mathbb{R}^{2M+1}$.

By the closure of GPs under linear operations, the output $f$ and the
projection $\mathbf{v}$ share a joint augmented GP prior. The covariance
matrix of inducing variables $\mathrm{Cov}[\mathbf{v},\mathbf{v}]$ has a low-
rank-plus-diagonal structure if inducing frequencies
$\mathbf{z}=[z_{m}]_{m=1}^{M}$ are harmonic on $[a,b]$, i.e.,
$z_{m}=\frac{2\pi m}{b-a}$, facilitating faster posterior computation (Hensman
et al., 2018). For a given input $t$, the cross-covariance of the output
process $f$ and the inducing variable $v_{m}$ is computed as

$\mathrm{Cov}[f(t),v_{m}]=\int_{-\infty}^{t}G(t-\tau)\cdot\langle
k(\tau,\cdot),\phi_{m}(\cdot)\rangle_{\mathcal{H}}\mathop{}\\!\mathrm{d}\tau,$
(11)

where we take advantage of the linearity and calculate the expectation over
$u(\cdot)$. The reproducing property of Matérn RKHS ensures that the inner
product $\langle k(\tau,\cdot),\phi_{m}(\cdot)\rangle_{\mathcal{H}}$ results
in well-defined sinusoidal functions within the interval $\tau\in[a,b]$.
Therefore, we can derive the RKHS Fourier features for LFMs on $t\in[a,b]$ as
follows:

$\mathrm{Cov}[f(t),v_{m}]=\int_{-\infty}^{t}G(t-\tau)\phi_{m}(\tau)\mathop{}\\!\mathrm{d}\tau=\begin{cases}\frac{\cos(z_{i}(t-a)+\theta)}{\beta\sqrt{z_{i}^{2}+\gamma^{2}}}+\xi_{i}&i=0,1,\ldots,M,\\\
\frac{\sin(z_{i}(t-a)+\theta)}{\beta\sqrt{z_{i}^{2}+\gamma^{2}}}+\xi_{i}&i=M+1,\ldots,2M,\end{cases}$
(12)

where $z_{0}=0$, cutoff frequency $\gamma=\frac{\alpha}{\beta}$, phase shift
$\theta=-\arctan(\frac{z_{i}}{\gamma})$, and $\xi_{i}$ represents an
exponential decay term. Since the integration variable $\tau$ ranges from
negative infinity and the inner product $\langle
k(\tau,\cdot),\phi(\cdot,z)\rangle_{\mathcal{H}}$ has different expressions
beyond $\tau\in[a,b]$, the covariance of $f$ and $v_{m}$ for $t\in\mathbb{R}$
emerges as a continuous piece-wise function (see Fig. 1 right). The detailed
derivation and complete expressions of $\text{Cov}[f(t),v(z)]$ for
Matérn-$\tfrac{1}{2}/\frac{3}{2}/\frac{5}{2}$ are given in Appx. 10 and
illustrated in Fig. 7 and 8.

The derived inter-domain features (12) reflect the filtering effect of the
system, i.e., how the ODE adaptively transforms the frequency components of
latent forces to the output through amplitude attenuation and phase shift. By
analogy with the frequency response of linear systems, we name the obtained
Fourier features from RKHS as _Variational Fourier Response Features_ (VFRFs).
Fig. 1, right, depicts the VFRFs of the LFM from the left subplot ($\gamma=4$,
solid red) and the VFFs of the corresponding latent force (dash blue). They
show that the VFRFs are learnable basis functions that adjust both the
amplitude and the phase according to the input frequencies and the ODE
parameters. Particularly, the system will allow nearly all frequency
components of the input process to pass through as $\gamma\rightarrow+\infty$.
Under this condition, VFRFs converge to the VFFs of inter-domain GPs. We would
like to emphasise here that the features derived from (11) can apply to more
general inter-domain GPs with other linear transformations $G(\cdot)$, not
just limited to LFMs. Moreover, for stable dynamical systems governed by
higher-order ODEs within the LFM framework, the derivation of VFRFs described
above can be readily extended by designing properly scaled Green’s functions
222The “stable” system here means its ODE has a general solution involving
exponentially decaying terms, in other words, the roots of the characteristic
equation lie in the left half of the complex plane. In this case, (11) may
still produce sinusoidal features..

###  3.3 Deep LFMs with VFRFs

DLFMs extend the concept of shallow LFMs by stacking them in a non-parametric
cascade, similar to DGPs. This hierarchical structure allows DLFMs to model
the non-stationarities present in nonlinear dynamical systems. In this
section, we detail the construction of a hierarchical composition of $L$ LFMs
within the framework of variational DGPs, each governed by the modified ODE
and enhanced with VFRFs for variational approximation of the posterior. We
leverage the layer-wise Monte Carlo technique in doubly stochastic variational
inference (Salimbeni and Deisenroth, 2017) to allow functional samples to
propagate through the compositional architecture efficiently.

The first layer of a DLFM processes a $D^{0}$-dimensional input
$\mathbf{t}=[t_{d}]_{d=1}^{D^{0}}$ and outputs a $D^{1}$-dimensional
independent process
$\boldsymbol{g}^{1}(\mathbf{t})=[g^{1}_{r}(\mathbf{t})]_{r=1}^{D^{1}}$ (the
superscripts indicate the layer index). To extend the application of VFRFs to
multidimensional inputs, we follow Hensman et al. (2018) to employ additive
LFM kernels for each output dimension, i.e., each output dimension
$g^{1}_{r}(\mathbf{t})$ is modelled as
$g^{1}_{r}(\mathbf{t})=\sum_{d=1}^{D^{0}}f^{1}_{d}(t_{d})$, where
$\\{f^{1}_{d}\\}_{d=1}^{D^{0}}$ are LFMs with ODE-induced covariance
functions. In this work, we assume the LFMs $f^{1}_{d}$ are independent, but
this assumption can be relaxed by allowing them to share the same latent
forces, which can lead to more complex kernels for the outputs.

Following the construction of a single-layer LFM, each $g^{1}_{r}(\mathbf{t})$
is equipped with a set of $M$ inducing frequencies
$\mathbf{Z}^{0}\in\mathbb{R}^{M\times D^{0}}$ and corresponding inducing
variables $\mathbf{V}^{1}=[v^{1}_{m,d}]^{2M+1,D^{0}}_{m=1,d=1}$. These
variables are created by the RKHS projection $v^{1}_{m,d}=\langle
u^{1}_{d},\phi_{m}\rangle_{\mathcal{H}}$. Therefore, the covariance functions
necessary for sparse variational inference are given by

$\mathrm{Cov}[g^{1}_{r}(\mathbf{t}),v^{1}_{m,d}]=\int_{-\infty}^{t_{d}}G^{1}_{d}(t_{d}-\tau)\phi_{m}(\tau)\mathop{}\\!\mathrm{d}\tau,\quad\mathrm{Cov}[v^{1}_{m,d},v^{1}_{m^{\prime},d^{\prime}}]=0\
\ (d\neq d^{\prime}).$ (13)

Assuming a variational distribution $q(\mathbf{V}^{1})$, the approximate
posterior $q(\boldsymbol{g}^{1}\mid\mathbf{t})$ of the first layer is derived
by substituting $k(\mathbf{x},\mathbf{Z})$ and $\mathbf{K_{ZZ}}$ in (4) with
the expressions of VFRFs. Samples from this approximate posterior are drawn
using the re-parameterisation trick (Kingma et al., 2015).

Given a training dataset with inputs $[\mathbf{t}_{i}]_{i=1}^{N}$ and targets
$[\mathbf{y}_{i}]_{i=1}^{N}$, the training of DLFMs involves maximising the
average ELBO over a mini-batch $\mathcal{B}$:

$\mathrm{ELBO}=\frac{1}{|\mathcal{B}|}\sum_{i\in\mathcal{B}}\mathbb{E}_{q(\mathbf{g}^{L}_{i}\mid\mathbf{t}_{i})}\log\left[p(\mathbf{y}_{i}\mid\mathbf{g}^{L}_{i})\right]-\frac{1}{N}\sum_{l=1}^{L}\mathrm{KL}\left[q(\mathbf{V}^{l})\mid\mid
p(\mathbf{V}^{l})\right],$ (14)

where $\mathbf{V}^{l}$ are the inducing variables of the $l$-th layer, and the
collection $\\{\mathbf{g}^{l}\\}^{L}_{l=1}$ denotes the output random
variables at hidden layers. The output of each layer serves as the input for
the subsequent layer, creating a chain of dependencies where the posterior of
each layer is computed based on the propagated samples
$\hat{\mathbf{g}}_{i}^{l}\sim
q(\mathbf{g}_{i}^{l}\mid\hat{\mathbf{g}}_{i}^{l-1})$. The predictive
distribution $q(\mathbf{y}_{*})=\int
p(\mathbf{y}_{*}\mid\mathbf{g}^{L}_{*})q(\mathbf{g}^{L}_{*})\mathop{}\\!\mathrm{d}\mathbf{g}^{L}_{*}$
at test location $\mathbf{t}_{*}$ follows a similar layer-wise procedure,
where $q(\mathbf{g}^{L}_{*})$ is a Gaussian mixture of $S$ hidden-layer
samples:
$q(\mathbf{g}_{*}^{L})\approx\frac{1}{S}\sum_{s=1}^{S}q(\mathbf{g}_{*}^{L}\mid\hat{\mathbf{g}}_{*}^{(s)^{L-1}})$.

##  4 Related Work

LFMs presented a physically-inspired approach to combining data-driven
modelling with differential equations (Alvarez et al., 2009). Álvarez et al.
(2010) further proposed variational inducing functions to handle non-smooth
latent processes within convolved GPs (Alvarez and Lawrence, 2011). Our model
builds upon LFMs and DGPs (Salimbeni and Deisenroth, 2017). Recently, various
approximate inference methods have been explored for DGP-based models, which
are generally categorised into variational inference techniques (Salimbeni and
Deisenroth, 2017; Salimbeni et al., 2019; Lindinger et al., 2020) and Monte
Carlo approaches(Havasi et al., 2018).

As outlined in Section 1, RFFs (Rahimi and Recht, 2007) and VFFs (Hensman et
al., 2018) have recently been incorporated into GP models. RFFs were used in
shallow LFMs models to approximate covariance matrices (Guarnizo and Álvarez,
2018) and expanded to a deeper architecture (McDonald and Álvarez, 2021,
2023). While related to these studies, our approach primarily uses features
similar to VFFs within the scope of inter-domain GPs (Van der Wilk et al.,
2020). Unlike RFF-based DGP models, which often alter the original covariance
functions by introducing a fully parametric variational distribution over
random frequencies, our model maintains the integrity of the original kernel
forms and approximates the DGP posterior directly. Another closely related
work is inter-domain DGPs (IDDGPs) (Rudner et al., 2020), which employ fixed
VFFs. In contrast, our model extends the compositional inter-domain GPs by
integrating ODEs to provide trainable, physics-informed RKHS Fourier features.
Appx. 7 illustrates a high-level comparison of our work with related works.

##  5 Experiments

This section presents experiments designed to evaluate our model using VFRFs.
We begin by examining the approximation quality of shallow LFMs with VFRFs and
RFFs on synthetic data. Additionally, we evaluate our model on a highly non-
stationary speech signal dataset, as well as on benchmark regression tasks,
comparing to a variety of baselines in both cases. We will make the
implementation publicly available after the review period.

##### Posterior Approximation of Shallow LFMs

VFRFs and RFFs both leverage Fourier features to facilitate approximate
inference in LFMs. In Fig. 2, we compare the approximation quality of VFRFs
and RFFs in a regression task using models with a Matérn-$\frac{1}{2}$ kernel.
The kernel’s parameters and the noise variance are initially optimised by
maximising the marginal likelihood of an exact LFM, then fixed across all
models. The frequencies of RFFs are sampled from the corresponding Cauchy
distribution (detailed in Appx. 9.2). Fig. 2 shows that the model with 20
VFRFs fits the data points reasonably well. Increasing VFRFs to 80 fills in
the details of the region with more observations, and the approximate
predictive posterior is quite close to the exact one. The same number of RFFs
yields a poorly fitted predictive mean and tends to underestimate the variance
in different regions of the input space, a phenomenon known as _variance
starvation_ (Wang et al., 2018). Achieving a comparable approximation requires
more RFFs due to the heavy-tailed spectral density of Matérn-$\frac{1}{2}$
kernels.

Figure 2: Illustrative example of Matérn-$\tfrac{1}{2}$ LFM posterior with
VFRFs / RFFs. The model’s feature is indicated at the lower right. Top row:
predictive posteriors of 20, 80, and 500 RFFs. Bottom row: predictive
posteriors of 20 and 80 inducing frequencies and an exact LFM. Noisy
observations are marked with red dots, posterior predictive means with blue
lines, and uncertainty (one or two standard deviations) with varying shades of
blue. In this example, VFRFs show a better approximation to the true
posterior, whereas RFRFs indicate variance underestimation with fewer
features. Figure 3: Posterior predictive distribution comparison of different
models on data points from a noisy multi-step function. The models and the
features used are noted at the bottom right of each subplot. The dashed lines
are samples from the predictive distributions. The experiment uses two layers
for deep models and Matérn-$\tfrac{3}{2}$ kernels except for the DGP (upper
left) and DLFM-RFF (lower left) that use RBF kernels. All models are trained
with 20 inducing points/Fourier features per layer. The DLFM models with VFRFs
perform best among the models.

##### Synthetic Data for Deep Structures

To further evaluate our DLFM-VFRF’s performance against other models, we
conduct tests on a synthetic multi-step function (Fig. 3). This task is
challenging for shallow GP models due to the need to capture global structures
in highly non-stationary data (Rudner et al., 2020). Although equipped with
VFRFs, our single-layer LFM (upper left) struggles to fit the non-stationarity
with a stationary kernel. In contrast, the models with compositional layers
exhibit better performance. The DLFM-RFF (lower left) (McDonald and Álvarez,
2021) generates high-frequency, wiggly posterior predictive samples, resulting
in an easily over-fitting model struggling to seize the slow-changing trend in
the data. DGP-RBF (upper left), IDDGP (lower middle), and our DLFM-VFRF (lower
right) offer smoother samples from the posterior distributions. Due to the
VFRFs’ flexibility to capture the global data structure, our DLFM-VFRF
outperforms both the DGP-RBF with local inducing points and the IDDGP,
especially inside $[0.8,1]$. Our model provides a more accurate predictive
mean throughout the steps and at abrupt step transitions, and demonstrates
narrower confidence intervals, indicating a better uncertainty calibration. In
the plot, GP models based on function-space inference tend to revert to prior
distributions outside the observed data range, displaying wide uncertainty
bands. In contrast, the DLFM-RFF yields relatively more confident non-zero
predictions in these areas.

##### TIMIT Speech Signals

We apply our model to a regression task on the TIMIT dataset, a speech
recognition resource previously used by Rudner et al. (2020), to explore the
capability of GP-based models in handling complex, non-stationary data. The
dataset features rapidly changing audio waves, posing significant challenges
for shallow GP models reliant on local approximation. Initially, we apply a
moving average filter to smooth the audio waves and select the first 10,000
data points for our analysis, reserving 30 % as test data. Our method uses the
Matérn-$\frac{3}{2}$ kernels.

(a) Learning progression of DLFMs and IDDGPs

(b) Standardised RMSE and NMLL

Figure 4: (a) Learning progression of DLFMs and IDDGPs with $M$ inducing
frequencies on the TIMIT dataset, presented in negative ELBO, test average
RMSE and NMLL. The DLFM in yellow maintains fixed $\beta=10^{-6}$ throughout
the first 2000 training epochs, after which $\alpha/\beta$ are allowed to
vary. The DLFM in red employs trainable ODE parameters from the start. The
DLFMs with VFRFs consistently outperform the IDDGPs; (b) Mean standardised
RMSE and NMLL with the standard deviations (over 10 random seeds) for models
employing varying numbers of inducing frequencies. The numbers following the
hyphen in the y-axis labels indicate the number of inducing
frequencies/points. A lower value (to the left) indicates better performance.

One of our goals is to evaluate how the performance of IDDGPs and DLFMs varies
with the number of global Fourier features and the effect of the ODE
parameters $\alpha$ and $\beta$ on the learning process. Fig. 4(a) illustrates
the progressions of performance metrics, e.g., test Root Mean Square Error
(RMSE) and Negative Marginal Log-Likelihood (NMLL) using 100 and 400 inducing
frequencies. The yellow lines of DLFMs align closely with the IDDGPs’ blue
line during the first 2000 iterations, where DLFMs maintain a fixed small
$\beta$ value ($\beta=10^{-6},\alpha=1$). We use this setting to illustrate
how our DLFM with VFRFs can replicate the original IDDGPs as expected. After
$2000$ iterations, we allow optimisation of $\alpha$ and $\beta$, leading to
subsequent improvements in the testing metrics, and suggesting continuous
learning with ODE-based Fourier features. Additionally, optimising all
parameters from the beginning (red lines) yields the best results across
various setups. Fig. 4(b) compares the performance of different models with
all parameters optimised from the beginning. It is unsurprising that
increasing the number of inducing frequencies typically results in better
performance. The results reveal that while RFF-based DLFMs exhibit the lowest
NMLL, they show the highest RMSE, reflecting a lack of precision. DLFMs
equipped with VFRFs consistently surpass both the DGP with local inducing
points and the IDDGPs in both RMSE and NMLL, highlighting our model’s enhanced
ability to accurately capture the global structure and non-stationarity of the
data.

##### UCI Regression Benchmarks

To demonstrate the versatility and effectiveness of our model on domain-
agnostic real-world data, we conduct evaluations on five diverse regression
datasets from the UCI Machine Learning Repository (Dua and Graff, 2019). These
datasets vary in size and feature dimensionality, allowing us to test the
model’s adaptability across different scenarios (see Fig. 5). Consistent with
standard practice (Salimbeni and Deisenroth, 2017), our regression tasks
involve multivariate inputs and a univariate target. We reserve 10% of each
dataset for testing, normalise the inputs to the range $[0,3]$, and
standardise the outputs based on the mean and standard deviation of the
training set (these transformations are then reversed for evaluation). All
models run with Matérn-$\frac{3}{2}$ kernels, employing 20 inducing points or
frequencies and three output dimensions per layer. We maintained the same
hyperparameters and initialisation across all tests. The figure illustrates
that our models achieve comparable performance to the baselines. Notably, our
models with two layers outperform IDDGP counterparts on the Energy, Power, and
Naval datasets. We also observe that increasing the number of layers generally
enhances the model’s representational capacity, resulting in improved
performance.

Figure 5: Regression test RMSE and NMLL results on UCI datasets, averaged over
10 random seeds. Lower values (to the left) indicate better performance. Model
names include the number of layers.

##  6 Conclusion

In this work, we proposed a variational approach to incorporating RKHS Fourier
features in compositional GP models. Particularly, we designed a set of
adaptive inter-domain Fourier features for hierarchical LFMs with modified
ODEs. Our empirical evaluations across various datasets demonstrate that our
model extends inter-domain GPs with RKHS Fourier features and has an enhanced
ability compared to conventional GPs with local approximation in fitting
complex data with non-stationary and global patterns. Future work will focus
on extending our model’s use to other challenging machine-learning tasks
requiring the integration of specific domain knowledge of higher-order ODEs.

## References

  * Alvarez et al. [2009] Mauricio Alvarez, David Luengo, and Neil D Lawrence.  Latent force models.  In _Artificial Intelligence and Statistics_ , pages 9–16. PMLR, 2009.
  * Álvarez et al. [2010] Mauricio Álvarez, David Luengo, Michalis Titsias, and Neil D Lawrence.  Efficient multioutput gaussian processes through variational inducing kernels.  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ , pages 25–32. JMLR Workshop and Conference Proceedings, 2010.
  * Alvarez and Lawrence [2011] Mauricio A Alvarez and Neil D Lawrence.  Computationally efficient convolved multiple output gaussian processes.  _The Journal of Machine Learning Research_ , 12:1459–1500, 2011.
  * Cutajar et al. [2017] Kurt Cutajar, Edwin V Bonilla, Pietro Michiardi, and Maurizio Filippone.  Random feature expansions for deep gaussian processes.  In _International Conference on Machine Learning_ , pages 884–893. PMLR, 2017.
  * Damianou and Lawrence [2013] Andreas Damianou and Neil D Lawrence.  Deep gaussian processes.  In _Artificial Intelligence and Statistics_ , pages 207–215. PMLR, 2013.
  * Dua and Graff [2019] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2019.  URL http://archive.ics.uci.edu/ml.  Irvine, CA: University of California, School of Information and Computer Sciences.
  * Durrande et al. [2016] Nicolas Durrande, James Hensman, Magnus Rattray, and Neil D Lawrence.  Detecting periodicities with gaussian processes.  _PeerJ Computer Science_ , 2:e50, 2016.
  * Duvenaud et al. [2014] David Duvenaud, Oren Rippel, Ryan Adams, and Zoubin Ghahramani.  Avoiding pathologies in very deep networks.  In Samuel Kaski and Jukka Corander, editors, _Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics_ , volume 33 of _Proceedings of Machine Learning Research_ , pages 202–210, Reykjavik, Iceland, 22–25 Apr 2014. PMLR.  URL https://proceedings.mlr.press/v33/duvenaud14.html.
  * Gal and Turner [2015] Yarin Gal and Richard Turner.  Improving the gaussian process sparse spectrum approximation by representing uncertainty in frequency inputs.  In _International Conference on Machine Learning_ , pages 655–664. PMLR, 2015.
  * Gardner et al. [2018] Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson.  Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Guarnizo and Álvarez [2018] Cristian Guarnizo and Mauricio A Álvarez.  Fast kernel approximations for latent force models and convolved multiple-output gaussian processes.  In _34th Conference on Uncertainty in Artificial Intelligence 2018, UAI 2018_ , pages 835–844. Association For Uncertainty in Artificial Intelligence (AUAI), 2018.
  * Havasi et al. [2018] Marton Havasi, José Miguel Hernández-Lobato, and Juan José Murillo-Fuentes.  Inference in deep gaussian processes using stochastic gradient hamiltonian monte carlo.  _Advances in Neural Information Processing Systems_ , 31, 2018.
  * Hensman et al. [2013] James Hensman, Nicolò Fusi, and Neil D. Lawrence.  Gaussian processes for big data.  In _Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence_ , UAI’13, page 282–290, Arlington, Virginia, USA, 2013. AUAI Press.
  * Hensman et al. [2015] James Hensman, Alexander Matthews, and Zoubin Ghahramani.  Scalable variational gaussian process classification.  In _Artificial Intelligence and Statistics_ , pages 351–360. PMLR, 2015.
  * Hensman et al. [2018] James Hensman, Nicolas Durrande, and Arno Solin.  Variational fourier features for gaussian processes.  _Journal of Machine Learning Research_ , 18(151):1–52, 2018.
  * Kingma et al. [2015] Durk P Kingma, Tim Salimans, and Max Welling.  Variational dropout and the local reparameterization trick.  _Advances in Neural Information Processing Systems_ , 28, 2015.
  * Lawrence et al. [2006] Neil Lawrence, Guido Sanguinetti, and Magnus Rattray.  Modelling transcriptional regulation using gaussian processes.  _Advances in Neural Information Processing Systems_ , 19, 2006.
  * Lázaro-Gredilla and Figueiras-Vidal [2009a] Miguel Lázaro-Gredilla and Anibal Figueiras-Vidal.  Inter-domain gaussian processes for sparse inference using inducing features.  _Advances in Neural Information Processing Systems_ , 22, 2009a.
  * Lázaro-Gredilla and Figueiras-Vidal [2009b] Miguel Lázaro-Gredilla and Anibal Figueiras-Vidal.  Inter-domain gaussian processes for sparse inference using inducing features.  _Advances in Neural Information Processing Systems_ , 22, 2009b.
  * Lázaro-Gredilla et al. [2010] Miguel Lázaro-Gredilla, Joaquin Quinonero-Candela, Carl Edward Rasmussen, and Aníbal R Figueiras-Vidal.  Sparse spectrum gaussian process regression.  _The Journal of Machine Learning Research_ , 11:1865–1881, 2010.
  * Lindinger et al. [2020] Jakob Lindinger, David Reeb, Christoph Lippert, and Barbara Rakitsch.  Beyond the mean-field: Structured deep gaussian processes improve the predictive uncertainties.  _Advances in Neural Information Processing Systems_ , 33:8498–8509, 2020.
  * McDonald and Álvarez [2021] Thomas McDonald and Mauricio Álvarez.  Compositional modeling of nonlinear dynamical systems with ode-based random features.  _Advances in Neural Information Processing Systems_ , 34:13809–13819, 2021.
  * McDonald and Álvarez [2023] Thomas McDonald and Mauricio Álvarez.  Deep latent force models: Ode-based process convolutions for bayesian deep learning.  _arXiv preprint arXiv:2311.14828_ , 2023.
  * Rahimi and Recht [2007] Ali Rahimi and Benjamin Recht.  Random features for large-scale kernel machines.  _Advances in Neural Information Processing Systems_ , 20, 2007.
  * Rasmussen and Williams [2006] Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  The MIT Press, 2006.  ISBN 026218253X.  URL http://www.gaussianprocess.org/gpml.
  * Roberts et al. [2013] Stephen Roberts, Michael Osborne, Mark Ebden, Steven Reece, Neale Gibson, and Suzanne Aigrain.  Gaussian processes for time-series modelling.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 371(1984):20110550, 2013.
  * Rudin [2017] Walter Rudin.  _Fourier analysis on groups_.  Courier Dover Publications, 2017.
  * Rudner et al. [2020] Tim GJ Rudner, Dino Sejdinovic, and Yarin Gal.  Inter-domain deep gaussian processes.  In _International Conference on Machine Learning_ , pages 8286–8294. PMLR, 2020.
  * Salimbeni and Deisenroth [2017] Hugh Salimbeni and Marc Deisenroth.  Doubly stochastic variational inference for deep gaussian processes.  _Advances in Neural Information Processing Systems_ , 30, 2017.
  * Salimbeni et al. [2019] Hugh Salimbeni, Vincent Dutordoir, James Hensman, and Marc Deisenroth.  Deep gaussian processes with importance-weighted variational inference.  In _International Conference on Machine Learning_ , pages 5589–5598. PMLR, 2019.
  * Snoek et al. [2012] Jasper Snoek, Hugo Larochelle, and Ryan P Adams.  Practical bayesian optimization of machine learning algorithms.  _Advances in neural information processing systems_ , 25, 2012.
  * Titsias [2009] Michalis Titsias.  Variational learning of inducing variables in sparse gaussian processes.  In David van Dyk and Max Welling, editors, _Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics_ , volume 5 of _Proceedings of Machine Learning Research_ , pages 567–574, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16–18 Apr 2009. PMLR.
  * Van der Wilk et al. [2020] Mark Van der Wilk, Vincent Dutordoir, ST John, Artem Artemev, Vincent Adam, and James Hensman.  A framework for interdomain and multioutput gaussian processes.  _arXiv preprint arXiv:2003.01115_ , 2020.
  * Wang et al. [2018] Zi Wang, Clement Gehring, Pushmeet Kohli, and Stefanie Jegelka.  Batched large-scale bayesian optimization in high-dimensional spaces.  In _International Conference on Artificial Intelligence and Statistics_ , pages 745–754. PMLR, 2018.

## Appendix

##  7 Model Comparison

We give a high-level comparison of our approach with other related models in
Fig. 6. Each point represents a corresponding model. The reference is attached
to the blue tag. The comparison dimensions include whether the model has a
multi-layer structure, whether it incorporates physics-informed modelling
involving ODEs and convolutions, and which kind of feature it uses. The red
points represent our work.

$\mathbf{Feature}$$\mathbf{Physics}$-$\mathbf{informed}$$\mathbf{Structure}$$Inducing\
Points$$RFF$$VFF$$no$$yes$$shallow$$deep$SVGP [Titsias, 2009]LFM [Alvarez et
al., 2009]DGP [Damianou and Lawrence, 2013][Salimbeni and Deisenroth,
2017]SSGP[Lázaro-Gredilla et al., 2010][Gal and Turner, 2015]DGP-RFF[Cutajar
et al., 2017]DLFM-RFF[McDonald and Álvarez, 2021]LFM-RFF[Guarnizo and Álvarez,
2018]IDDGP[Rudner et al., 2020]GP-VFF[Hensman et al., 2018]This work Figure 6:
An illustration on a comparison of our model with recent related work. The
comparison dimensions are the used feature, the structure depth, and whether
incorporating physical dynamics.

##  8 Variational Fourier Feature

Variational Fourier Features (VFFs, [15]) are designed on a Matérn Reproducing
Kernel Hilbert Space (RKHS). Specifically,
Matérn-$\frac{1}{2}/\frac{3}{2}/\frac{5}{2}$ kernels with inputs
$x,x^{\prime}\in\mathbb{R}$ are of particular interest:

$\displaystyle k_{1/2}(r)=\sigma^{2}e^{-\lambda r},\ \lambda=\frac{1}{l},$
(15) $\displaystyle k_{3/2}(r)=\sigma^{2}(1+\lambda r)e^{-\lambda r},\
\lambda=\frac{\sqrt{3}}{l},$ (16) $\displaystyle
k_{5/2}(r)=\sigma^{2}(1+\lambda r+\frac{1}{3}\lambda^{2}r^{2})e^{-\lambda r},\
\lambda=\frac{\sqrt{5}}{l},$ (17)

where $r=|x-x^{\prime}|$, $\sigma^{2}$ is the kernel’s output-scale (or
variance), and $l$ is the length-scale. We reiterate the closed-form RKHS
inner products for Matérn-$\frac{1}{2}$ on $[a,b]$ here ( for other Matérn
kernels and more details see Durrande et al. [7]) :

$\langle
g,h\rangle_{\mathcal{H}_{\frac{1}{2}}}=\frac{1}{2\lambda\sigma^{2}}\int_{a}^{b}(\lambda
g(x)+g^{\prime}(x))(\lambda
h(x)+h^{\prime}(x))\mathop{}\\!\mathrm{d}x+\frac{1}{\sigma^{2}}g(a)h(a).$ (18)

The explicit expressions of the Matérn RKHS not only allow us to verify the
reproducing property when $t\in[a,b]$

$\langle k(t,\cdot),h(\cdot)\rangle_{\mathcal{H}}=h(t),\ \forall
h\in\mathcal{H},t\in[a,b],$ (19)

but also make it feasible to complete the cross-covariance if $t$ is outside
$[a,b]$. In Appx. 10, we further utilise the conclusion of VFFs in [15] for
$t\in\mathbb{R}$ to calculate our _Variational Fourier Response Features_
(VFRFs).

Given a Matérn GP $f(t)\sim\mathcal{GP}(0,k(t,t^{\prime}))$, the explicit RKHS
inner product provides an alternative linear operator to construct an inter-
domain GP by $u_{m}=\langle f,\phi_{m}\rangle_{\mathcal{H}}$, where
$\phi_{m},m=0,\ldots,2M$ is from a set of truncated Fourier basis with
harmonic _inducing frequencies_

$\phi_{0}(\cdot)=1,\ \phi_{m}(\cdot)=\cos(z_{m}(\cdot-a)),\
\phi_{M+m}(\cdot)=\sin(z_{m}(\cdot-a)),\ z_{m}=\frac{2\pi m}{b-a}.$ (20)

$u(\cdot)$ is an _inter-domain_ GP sharing a joint Gaussian prior with
$f(\cdot)$. For $t\in[a,b]$, the covariances are

$\mathrm{Cov}[f(t),u_{m}]=\langle
k(t,\cdot),\phi_{m}(\cdot)\rangle_{\mathcal{H}}=\phi_{m}(t),\quad\mathrm{Cov}[u_{i},u_{j}]=\langle\phi_{i}(\cdot),\phi_{j}(\cdot)\rangle_{\mathcal{H}}.$
(21)

The VFFs approximate the posterior by replacing the covariance matrix in
Sparse Variational GPs (SVGPs) appropriately.

##  9 LFMs for First-Order Dynamical System

We recall in this work a dynamical system modelled by a first-order ODE

$\beta\frac{\mathop{}\\!\mathrm{d}f(t)}{\mathop{}\\!\mathrm{d}t}+\alpha
f(t)=u(t),\ u(t)\sim\mathcal{GP}(0,k(t,t^{\prime})),\ \alpha,\beta>0,$ (22)

where $u$ is an unobserved latent force with a Matérn kernel. The Green’s
function is $G(t)=\frac{1}{\beta}\exp\left(-\gamma
t\right),\gamma=\frac{\alpha}{\beta}$. We take the solution

$f(t)=\int_{-\infty}^{t}G(t-\tau)u(\tau)\mathop{}\\!\mathrm{d}\tau=G\circ u,$
(23)

where a convolutional operator ${G}$ acting on $u$ is represented as $f=G\circ
u$. Conventional LFMs establish a GP over $f\sim\mathcal{GP}(0,G\circ k\circ
G)$, where $k\circ G$ signifies $G$ operating the second argument of the
kernel. The LFM kernels are computed analytically in [17, 1], but their
expressions are based on RBF kernels. In the subsequent part, we will present
the closed-form covariance expressions $G\circ k\circ G$ for the
Matérn-$\frac{1}{2}/\frac{3}{2}/\frac{5}{2}$ kernels, respectively. All
analytical LFM covariance functions discussed in this work are summarised in
Table 1 and illustrated in the left panel of Fig. 1, 7 and 8. Furthermore, we
introduce the approximation of our LFM kernels using random Fourier features
in Appx. 9.2.

###  9.1 Analytical LFM Matérn Kernels

The LFM kernel of a Matérn-$\frac{1}{2}$ latent force when $t>t^{\prime}$ is
given by

$\displaystyle G\circ k\circ G$
$\displaystyle=\int_{-\infty}^{t}\int_{-\infty}^{t^{\prime}}\frac{1}{\beta}e^{-\gamma(t-\tau)}\cdot\sigma^{2}e^{-\lambda|\tau-\tau^{\prime}|}\cdot\frac{1}{\beta}e^{-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau\mathop{}\\!\mathrm{d}\tau^{\prime}$
$\displaystyle=\frac{\sigma^{2}}{\beta^{2}}\int_{-\infty}^{t^{\prime}}\int_{\tau^{\prime}}^{t}e^{-\gamma(t-\tau)-\lambda(\tau-\tau^{\prime})-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau\mathop{}\\!\mathrm{d}\tau^{\prime}$
$\displaystyle\quad+\frac{\sigma^{2}}{\beta^{2}}\int_{-\infty}^{t^{\prime}}\int_{-\infty}^{\tau^{\prime}}e^{-\gamma(t-\tau)-\lambda(\tau^{\prime}-\tau)-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau\mathop{}\\!\mathrm{d}\tau^{\prime}$
$\displaystyle=\frac{\sigma^{2}}{\beta^{2}\gamma(\gamma^{2}-\lambda^{2})}\left[\gamma
e^{-\lambda(t-t^{\prime})}-\lambda e^{-\gamma(t-t^{\prime})}\right].$ (24)

The derivation for $t<t^{\prime}$ is similar. As a result, we obtain a
stationary LFM kernel for $\forall t,t^{\prime}\in\mathbb{R}$,

$\displaystyle G\circ k\circ
G=\begin{cases}\frac{\sigma^{2}}{\beta^{2}\gamma(\gamma^{2}-\lambda^{2})}\left[\gamma
e^{-\lambda|t-t^{\prime}|}-\lambda e^{-\gamma|t-t^{\prime}|}\right]&\text{if}\
\gamma\neq\lambda,\\\
\frac{\sigma^{2}(1+\lambda|t-t^{\prime}|)}{2\beta^{2}\lambda^{2}}e^{-\lambda|t-t^{\prime}|}&\text{if}\
\gamma=\lambda.\end{cases}$ (25)

Likewise, the expressions of the other LFM kernels are present in Table 1. The
expressions exhibit continuity but non-differentiability at the point where
$\gamma=\lambda$.

Table 1: LFM kernels of Matérn-$\tfrac{1}{2}/\tfrac{3}{2}/\tfrac{5}{2}$ latent
forces

Latent Force $k(r)$ | LFM $G\circ k\circ G$ 1
---|---
Matérn-$\frac{1}{2}$ | $\frac{\sigma^{2}}{\beta^{2}\gamma(\gamma^{2}-\lambda^{2})}\left[\gamma e^{-\lambda r}-\lambda e^{-\gamma r}\right]$
Matérn-$\frac{3}{2}$ | $\frac{\sigma^{2}}{\beta^{2}}\left[\frac{\lambda r+1}{\gamma^{2}-\lambda^{2}}-\frac{2\lambda^{2}}{(\gamma^{2}-\lambda^{2})^{2}}\right]e^{-\lambda r}+\frac{2\lambda^{3}\sigma^{2}}{\beta^{2}\gamma(\gamma^{2}-\lambda^{2})^{2}}e^{-\gamma r}$
Matérn-$\frac{5}{2}$ | $\frac{\sigma^{2}}{3\beta^{2}}\left[\frac{\lambda^{2}r^{2}+3}{\gamma^{2}-\lambda^{2}}+\frac{\lambda(3\gamma^{2}-7\lambda^{2})r}{(\gamma^{2}-\lambda^{2})^{2}}+\frac{4\lambda^{2}(3\lambda^{2}-\gamma^{2})}{(\gamma^{2}-\lambda^{2})^{3}}\right]e^{-\lambda r}-\frac{8\lambda^{5}\sigma^{2}}{3\beta^{2}\gamma(\gamma^{2}-\lambda^{2})^{3}}e^{-\gamma r}$

  * 1

$r$ is the input distance $r=|t-t^{\prime}|$. $\\{\beta,\gamma\\}$ and
$\\{\sigma^{2},\lambda\\}$ are parameters from the ODE and the Matérn kernel,
respectively.

###  9.2 Random Fourier Approximation of Kernels

Stationary kernels can be approximated using random Fourier features [24]
using Bochner’s theorem [27]. In the illustrative experiment, we approximate
the Matérn LFM covariance proposed in our work using random Fourier features
(also termed _Random Fourier Response Features_ (RFRFs) by [11]). The random
Fourier features of the LFM with a Matérn-$\nu\
(\nu=\tfrac{1}{2}/\tfrac{3}{2}/\tfrac{5}{2})$ kernel is given by

$\varphi(t;\omega)=\int_{-\infty}^{t}e^{j\omega\tau}\cdot\frac{1}{\beta}e^{-\gamma(t-\tau)}\mathop{}\\!\mathrm{d}\tau=\frac{e^{j\omega
t}}{\beta(\gamma+j\omega)},\ \omega=\frac{\omega^{\prime}}{l},\
\omega^{\prime}\sim t_{2\nu}(\omega^{\prime}),$ (26)

where $j^{2}=-1$, $t_{2\nu}$ is a zero-mean Student’s $t$-distribution with
$2\nu$ degrees of freedom, and $l$ is the length-scale of the Matérn kernel.
Therefore,

$\displaystyle G\circ k$
$\displaystyle\approx\frac{\sigma^{2}}{M}\sum_{m=1}^{M}\varphi(t;\omega_{m})\cdot
e^{-j\omega_{m}t^{\prime}},$ (27) $\displaystyle G\circ k\circ G$
$\displaystyle\approx\frac{\sigma^{2}}{M}\sum_{m=1}^{M}\varphi(t;\omega_{m})\cdot\bar{\varphi}(t^{\prime};\omega_{m}),$
(28)

where $\sigma^{2}$ is Matérn kernel’s variance, $\\{\omega_{m}\\}_{m=1}^{M}$
are $M$ random Fourier frequencies sampled from the corresponding Student’s
$t$-distribution. $\bar{\varphi}$ denotes the complex conjugate of $\varphi$.

##  10 Variational Fourier Response Features for LFMs

We represent the projection of the latent force $u$ onto the truncated Fourier
basis $\phi$ as

$v(z)=P\circ u=\langle\phi(\cdot;z),u(\cdot)\rangle_{\mathcal{H}},$ (29)

where $v(z)$ is the projection process for the latent force, and $z$ is the
inducing frequency. For simplicity, the projection operator is denoted as $P$.
Consequently, the output process $f$ and the projection process $v$ share a
joint GP prior:

$\left[\begin{array}[]{c}f\\\
v\end{array}\right]\sim\mathcal{N}\left(\mathbf{0},\left[\begin{array}[]{cc}G\circ
k\circ G&G\circ k\circ P\\\ P\circ k\circ G&P\circ k\circ
P\end{array}\right]\right),$ (30)

where the covariance terms are given by

$\displaystyle\mathrm{Cov}[f(t),v(z)]$
$\displaystyle=\mathbb{E}\left[f(t)\langle
u(\cdot),\phi(\cdot;z)\rangle_{\mathcal{H}}\right]=\langle\mathbb{E}[f(t)u(\cdot)],\phi(\cdot;z)\rangle_{\mathcal{H}}$
$\displaystyle=\langle G\circ
k(t,\cdot),\phi(\cdot;z)\rangle_{\mathcal{H}}=G\circ\langle
k(t,\cdot),\phi(\cdot;z)\rangle_{\mathcal{H}}=G\circ k\circ P.$ (31)
$\displaystyle\mathrm{Cov}[v(z_{i}),v(z_{j})]$
$\displaystyle=\langle\phi(\cdot;z_{i}),\phi(\cdot;z_{j})\rangle_{\mathcal{H}}=P\circ
k\circ P.$ (32)

The rest of this section will specify the closed-form VFRF expressions of
$P\circ k\circ G$. With the integration lower limit going to negative
infinity, the input values outside the interval $[a,b]$ should be considered.

We detail the derivation of Matérn-$\frac{1}{2}$ (listed in Table 3 and 4) and
directly give the results of Matérn-$\frac{3}{2}$ (see Table 6 and 7) and
Matérn-$\frac{5}{2}$ (see Table 9 and 10). These expressions of VFRFs will
revert to VFFs under certain conditions. The LFM kernel and the VFRFs for
Matérn-$\frac{3}{2}$ and Matérn-$\frac{5}{2}$ with the same hyperparameters in
the main text are depicted in Fig. 7 and Fig. 8, respectively.

Figure 7: Covariance functions and VFRFs for LFMs with Matérn-$\frac{3}{2}$
kernel. Figure 8: Covariance functions and VFRFs for LFMs with
Matérn-$\frac{5}{2}$ kernel.

###  10.1 Matérn-$\frac{1}{2}$ Cross-covariance

We write the cross-covariance as

$P\circ k\circ
G=\langle\phi(\cdot;z),k(\cdot,\tau^{\prime})\rangle_{\mathcal{H}}\circ
G=\int_{-\infty}^{t^{\prime}}h(z,\tau^{\prime})\cdot
G(t^{\prime}-\tau^{\prime})\mathop{}\\!\mathrm{d}\tau^{\prime}.$ (33)

The location of $\tau$ determines the expression of
$h(z,\tau^{\prime})=\langle\phi(\cdot;z),k(\cdot,\tau^{\prime})\rangle_{\mathcal{H}}$
(i.e., $P\circ k$) [15], as summarized in the subsequent table (For
Matérn-$\frac{3}{2}/\frac{5}{2}$, see Table 5/8). The columns of the tables
indicate the VFFs for input $\tau^{\prime}$ located inside/outside $[a,b]$.
The cross-covariance can be derived by substituting $h(z,\tau^{\prime})$ in
(33). Table 3 and Table 4 collect the VFRFs for Matérn-$\frac{1}{2}$ LFMs with
cosine/sine projection basis functions at different locations of $t^{\prime}$
in (33).

Table 2: VFFs $h(z,\tau^{\prime})$ with Matérn-$\frac{1}{2}$ kernel $\phi(\cdot;z)$ | $\tau^{\prime}<a$ | $a\leq\tau^{\prime}\leq b$ | $\tau^{\prime}>b$
---|---|---|---
$\cos(z(\cdot-a))$ | $e^{-\lambda(a-\tau^{\prime})}$ | $\cos(z(\tau^{\prime}-a))$ | $e^{-\lambda(\tau^{\prime}-b)}$
$\sin(z(\cdot-a))$ | 0 | $\sin(z(\tau^{\prime}-a))$ | 0

####  10.1.1 Cosine Features ($\phi(t;z)=\cos(z(t-a))$)

##### Case 1:

$t^{\prime}<a$,

$P\circ k\circ
G=\int_{-\infty}^{t^{\prime}}e^{-\lambda(a-\tau^{\prime})}\cdot\frac{1}{\beta}e^{-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau^{\prime}=\frac{1}{\beta(\gamma+\lambda)}e^{-\lambda(a-t^{\prime})}.$
(34)

The covariance will converge to the VFF in Table 2 with a scaling coefficient
$\alpha$:

$\lim\limits_{\beta\to 0^{+}}P\circ k\circ
G=\frac{1}{\alpha}e^{-\lambda(a-t^{\prime})}.$ (35)

##### Case 2:

$a\leq t^{\prime}\leq b$,

$\displaystyle P\circ k\circ G$
$\displaystyle=\int_{-\infty}^{a}e^{-\lambda(a-\tau^{\prime})}\cdot\frac{1}{\beta}e^{-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau^{\prime}+\int_{a}^{t^{\prime}}\cos(z(\tau^{\prime}-a))\cdot\frac{1}{\beta}e^{-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau^{\prime}$
$\displaystyle=\frac{\gamma\cos(z(t^{\prime}-a))+z\sin(z(t^{\prime}-a))}{\beta(z^{2}+\gamma^{2})}+\frac{(z^{2}-\gamma\lambda)}{\beta(\gamma+\lambda)(z^{2}+\gamma^{2})}e^{-\gamma(t^{\prime}-a)}$
$\displaystyle=\frac{\cos(z(t^{\prime}-a)+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}+\xi_{\cos},$
(36)

where $\theta=-\arctan(\frac{z}{\gamma})$, and $\xi$ is a decay term.
Particularly, the cross-covariance extends the VFFs of the latent force since

$\lim\limits_{\beta\to 0_{+}}P\circ k\circ
G=\frac{1}{\alpha}\cos(z(t^{\prime}-a)),$ (37)

##### Case 3:

$t^{\prime}>b$,

$P\circ k\circ
G=\frac{e^{-\lambda(t^{\prime}-b)}}{\beta(\gamma-\lambda)}+\frac{(z^{2}-\gamma\lambda)e^{-\gamma(t^{\prime}-a)}}{\beta(\gamma+\lambda)(z^{2}+\gamma^{2})}-\frac{(z^{2}+\gamma\lambda)e^{-\gamma(t^{\prime}-b)}}{\beta(\gamma-\lambda)(z^{2}+\gamma^{2})},$
(38)

which utilizes harmonic $z=\frac{2\pi m}{b-a},m\in\mathbb{Z}_{+}$. The
covariance will also return to a scaled term in Table 2 when $\beta\rightarrow
0_{+}$.

####  10.1.2 Sine Features ($\phi(t;z)=\sin(z(t-a))$)

##### Case 1:

$t^{\prime}<a,\ P\circ k\circ G=0.$

##### Case 2:

$a\leq t^{\prime}\leq b$,

$\displaystyle P\circ k\circ G$
$\displaystyle=\int_{a}^{t^{\prime}}\sin(z(\tau^{\prime}-a))\cdot\frac{1}{\beta}e^{-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau^{\prime}$
$\displaystyle=\frac{-z\cos(z(t^{\prime}-a))+\gamma\sin(z(t^{\prime}-a))}{\beta(z^{2}+\gamma^{2})}+\frac{z}{\beta(z^{2}+\gamma^{2})}e^{-\gamma(t^{\prime}-a)}$
$\displaystyle=\frac{\sin(z(t^{\prime}-a)+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}+\xi_{\sin},\quad\theta=-\arctan(\frac{z}{\gamma})$
(39)

##### Case 3:

$t^{\prime}>b$,

$P\circ k\circ
G=\int_{a}^{b}\sin(z(\tau^{\prime}-a))\cdot\frac{1}{\beta}e^{-\gamma(t^{\prime}-\tau^{\prime})}\mathop{}\\!\mathrm{d}\tau^{\prime}=\frac{ze^{-\gamma(t^{\prime}-a)}-ze^{-\gamma(t^{\prime}-b)}}{\beta(\gamma^{2}+z^{2})}.$
(40)

The VFRFs of the LFM with a Matérn-$\frac{1}{2}$ latent force are summarised
in Table 3 and 4, where the absolute distances to the interval ends are
denoted as $r_{a}=|t^{\prime}-a|$ and $r_{b}=|t^{\prime}-b|$ and the phase
shift is $\theta=-\arctan(\frac{z}{\gamma})$. The features are continuous at
$\gamma=\lambda$ when $t^{\prime}>b$.

Table 3: Matérn-$\tfrac{1}{2}$ VFRFs on Fourier basis $\phi(x;z)=\cos(z(x-a))$.  $t^{\prime}\in\mathbb{R}$ | LFM Fourier Feature $P\circ k\circ G$ (cosine part)
---|---
$t^{\prime}<a$ | $\frac{1}{\beta(\gamma+\lambda)}e^{-\lambda r_{a}}$
$a\leq t^{\prime}\leq b$ | $\frac{\cos(zr_{a}+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}-\left[\frac{\gamma}{\beta(z^{2}+\gamma^{2})}-\frac{1}{\beta(\gamma+\lambda)}\right]e^{-\gamma r_{a}}$
$t^{\prime}>b\ (\gamma\neq\lambda)$ | $-\left[\frac{\gamma}{\beta(z^{2}+\gamma^{2})}-\frac{1}{\beta(\gamma+\lambda)}\right]e^{-\gamma r_{a}}+\left[\frac{\gamma}{\beta(z^{2}+\gamma^{2})}-\frac{1}{\beta(\gamma-\lambda)}\right]e^{-\gamma r_{b}}+\frac{1}{\beta(\gamma-\lambda)}e^{-\lambda r_{b}}$
$t^{\prime}>b\ (\gamma=\lambda)$ | $-\left[\frac{\lambda}{\beta(z^{2}+\lambda^{2})}-\frac{1}{2\beta\lambda}\right]e^{-\lambda r_{a}}+\left[\frac{\lambda}{\beta(z^{2}+\lambda^{2})}+\frac{r_{b}}{\beta}\right]e^{-\lambda r_{b}}$
Table 4: Matérn-$\tfrac{1}{2}$ VFRFs on Fourier basis $\phi(x;z)=\sin(z(x-a))$.  $t^{\prime}\in\mathbb{R}$ | LFM Fourier Feature $P\circ k\circ G$ (sine part)
---|---
$t^{\prime}<a$ |  $0$
$a\leq t^{\prime}\leq b$ | $\frac{\sin(zr_{a}+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}+\frac{z}{\beta(z^{2}+\gamma^{2})}e^{-\gamma r_{a}}$
$t^{\prime}>b$ | $\frac{z}{\beta(z^{2}+\gamma^{2})}e^{-\gamma r_{a}}-\frac{z}{\beta(z^{2}+\gamma^{2})}e^{-\gamma r_{b}}$

###  10.2 Matérn-$\frac{3}{2}/\frac{5}{2}$ Cross-covariance

Based on Table 5 and 8, we give the VFRFs for Matérn-$\frac{3}{2}$ and
Matérn-$\frac{5}{2}$ LFMs with $\theta=-\arctan(\frac{z}{\gamma})$ in Table
6,7,9 and 10. Also, the derived cross-covariance expressions can return to the
scaled VFFs of the latent force.

Table 5: VFFs $h(z,\tau^{\prime})$ with Matérn-$\frac{3}{2}$ kernel $\phi(\cdot;z)$ | $\tau^{\prime}<a$ | $a\leq\tau^{\prime}\leq b$ | $\tau^{\prime}>b$
---|---|---|---
$\cos(z(\cdot-a))$ | $(1+\lambda(a-\tau^{\prime}))e^{-\lambda(a-\tau^{\prime})}$ | $\cos(z(\tau^{\prime}-a))$ | $(1+\lambda(\tau^{\prime}-b))e^{-\lambda(\tau^{\prime}-b)}$
$\sin(z(\cdot-a))$ | $z(\tau^{\prime}-a)e^{-\lambda(a-\tau^{\prime})}$ | $\sin(z(\tau^{\prime}-a))$ | $z(\tau^{\prime}-b)e^{-\lambda(\tau^{\prime}-b)}$
Table 6: Matérn-$\tfrac{3}{2}$ VFRFs on Fourier basis $\phi(x;z)=\cos(z(x-a))$.  $t^{\prime}\in\mathbb{R}$ | LFM Fourier Feature $P\circ k\circ G$ (cosine part)
---|---
$t^{\prime}<a$ | $\left[\frac{\lambda r_{a}+1}{\beta(\gamma+\lambda)}+\frac{\lambda}{\beta(\gamma+\lambda)^{2}}\right]e^{-\lambda r_{a}}$
$a\leq t^{\prime}\leq b$ | $\frac{\cos(zr_{a}+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}-\left[\frac{\gamma}{\beta(z^{2}+\gamma^{2})}-\frac{\gamma+2\lambda}{\beta(\gamma+\lambda)^{2}}\right]e^{-\gamma r_{a}}$
$t^{\prime}>b\ (\gamma\neq\lambda)$ |  $-\left[\frac{\gamma}{\beta(z^{2}+\gamma^{2})}-\frac{\gamma+2\lambda}{\beta(\gamma+\lambda)^{2}}\right]e^{-\gamma r_{a}}+\left[\frac{\gamma}{\beta(z^{2}+\gamma^{2})}-\frac{\gamma-2\lambda}{\beta(\gamma-\lambda)^{2}}\right]e^{-\gamma r_{b}}$ $\quad+\left[\frac{\lambda r_{b}+1}{\beta(\gamma-\lambda)}-\frac{\lambda}{\beta(\gamma-\lambda)^{2}}\right]e^{-\lambda r_{b}}$
$t^{\prime}>b\ (\gamma=\lambda)$ | $-\left[\frac{\lambda}{\beta(z^{2}+\lambda^{2})}-\frac{3}{4\beta\lambda}\right]e^{-\lambda r_{a}}+\left[\frac{\lambda}{\beta(z^{2}+\lambda^{2})}+\frac{(\lambda r_{b}+2)r_{b}}{2\beta}\right]e^{-\lambda r_{b}}$
Table 7: Matérn-$\tfrac{3}{2}$ VFRFs on Fourier basis $\phi(x;z)=\sin(z(x-a))$.  $t^{\prime}\in\mathbb{R}$ | LFM Fourier Feature $P\circ k\circ G$ (sine part)
---|---
$t^{\prime}<a$ | $-\frac{z}{\beta}\left[\frac{r_{a}}{\gamma+\lambda}+\frac{1}{(\gamma+\lambda)^{2}}\right]e^{-\lambda r_{a}}$
$a\leq t^{\prime}\leq b$ | $\frac{\sin(zr_{a}+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}+\frac{z}{\beta}\left[\frac{1}{(z^{2}+\gamma^{2})}-\frac{1}{(\gamma+\lambda)^{2}}\right]e^{-\gamma r_{a}}$
$t^{\prime}>b\ (\gamma\neq\lambda)$ |  $\frac{z}{\beta}\left[\frac{1}{(z^{2}+\gamma^{2})}-\frac{1}{(\gamma+\lambda)^{2}}\right]e^{-\gamma r_{a}}-\frac{z}{\beta}\left[\frac{1}{(z^{2}+\gamma^{2})}-\frac{1}{(\gamma-\lambda)^{2}}\right]e^{-\gamma r_{b}}$ $\quad+\frac{z}{\beta}\left[\frac{r_{b}}{\gamma-\lambda}-\frac{1}{(\gamma-\lambda)^{2}}\right]e^{-\lambda r_{b}}$
$t^{\prime}>b\ (\gamma=\lambda)$ | $\frac{z}{\beta}\left[\frac{1}{(z^{2}+\lambda^{2})}-\frac{1}{4\lambda^{2}}\right]e^{-\lambda r_{a}}+\frac{z}{\beta}\left[-\frac{1}{(z^{2}+\lambda^{2})}+\frac{r_{b}^{2}}{2}\right]e^{-\lambda r_{b}}$
Table 8: VFFs $h(z,\tau^{\prime})$ with Matérn-$\frac{5}{2}$ kernel 1

$\phi(\cdot,z)$ | $\tau^{\prime}<a$ | $a\leq\tau^{\prime}\leq b$ | $\tau^{\prime}>b$
---|---|---|---
$\cos(z(\cdot-a))$ | $\left[1+\lambda r-\frac{(z^{2}-\lambda^{2})r^{2}}{2}\right]e^{-\lambda r}$ | $\cos(z(\tau^{\prime}-a))$ | $\left[1+\lambda r-\frac{(z^{2}-\lambda^{2})r^{2}}{2}\right]e^{-\lambda r}$
$\sin(z(\cdot-a))$ | $z(\tau^{\prime}-a)(1+\lambda r)e^{-\lambda r}$ | $\sin(z(\tau^{\prime}-a))$ | $z(\tau^{\prime}-b)(1+\lambda r)e^{-\lambda r}$

  * 1

$r=\min\\{|\tau^{\prime}-a|,|\tau^{\prime}-b|\\}.$

Table 9: Matérn-$\tfrac{5}{2}$ VFRFs on Fourier basis $\phi(x;z)=\cos(z(x-a))$.  $t^{\prime}\in\mathbb{R}$ | LFM Fourier Feature $P\circ k\circ G$ (cosine part)
---|---
$t^{\prime}<a$ | $-\left[\frac{(z^{2}-\lambda^{2})r_{a}^{2}}{2\beta(\gamma+\lambda)}+\frac{(z^{2}-\gamma\lambda-2\lambda^{2})r_{a}}{\beta(\gamma+\lambda)^{2}}+\frac{z^{2}-\gamma^{2}-3\gamma\lambda-3\lambda^{2}}{\beta(\gamma+\lambda)^{3}}\right]e^{-\lambda r_{a}}$
$a\leq t^{\prime}\leq b$ | $\frac{\cos(zr_{a}+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}-\left[\frac{z^{2}-\gamma^{2}-3\gamma\lambda-3\lambda^{2}}{\beta(\gamma+\lambda)^{3}}+\frac{\gamma}{\beta(z^{2}+\gamma^{2})}\right]e^{-\gamma r_{a}}$
$t^{\prime}>b\ (\gamma\neq\lambda)$ |  $-\left[\frac{z^{2}-\gamma^{2}-3\gamma\lambda-3\lambda^{2}}{\beta(\gamma+\lambda)^{3}}+\frac{\gamma}{\beta(z^{2}+\gamma^{2})}\right]e^{-\gamma r_{a}}+\left[\frac{z^{2}-\gamma^{2}+3\gamma\lambda-3\lambda^{2}}{\beta(\gamma-\lambda)^{3}}+\frac{\gamma}{\beta(z^{2}+\gamma^{2})}\right]e^{-\gamma r_{b}}$ $\quad\ -\left[\frac{(z^{2}-\lambda^{2})r_{b}^{2}}{2\beta(\gamma-\lambda)}-\frac{(z^{2}+\gamma\lambda-2\lambda^{2})r_{b}}{\beta(\gamma-\lambda)^{2}}+\frac{z^{2}-\gamma^{2}+3\gamma\lambda-3\lambda^{2}}{\beta(\gamma-\lambda)^{3}}\right]e^{-\lambda r_{b}}$
$t^{\prime}>b\ (\gamma=\lambda)$ | $-\left[\frac{z^{2}-7\lambda^{2}}{8\beta\lambda^{3}}+\frac{\lambda}{\beta(z^{2}+\lambda^{2})}\right]e^{-\lambda r_{a}}+\left[\frac{\lambda}{\beta(z^{2}+\lambda^{2})}-\frac{[(z^{2}-\lambda^{2})r_{b}^{2}-3\lambda r_{b}-6]r_{b}}{6\beta}\right]e^{-\lambda r_{b}}$
Table 10: Matérn-$\tfrac{5}{2}$ VFRFs on Fourier basis $\phi(x;z)=\sin(z(x-a))$.  $t^{\prime}\in\mathbb{R}$ | LFM Fourier Feature $P\circ k\circ G$ (sine part)
---|---
$t^{\prime}<a$ | $-\frac{z}{\beta}\left[\frac{\lambda r_{a}^{2}}{(\gamma+\lambda)}+\frac{(\gamma+3\lambda)r_{a}}{(\gamma+\lambda)^{2}}+\frac{\gamma+3\lambda}{(\gamma+\lambda)^{3}}\right]e^{-\lambda r_{a}}$
$a\leq t^{\prime}\leq b$ | $\frac{\sin(zr_{a}+\theta)}{\beta\sqrt{z^{2}+\gamma^{2}}}+\frac{z}{\beta}\left[\frac{1}{(z^{2}+\gamma^{2})}-\frac{(\gamma+3\lambda)}{(\gamma+\lambda)^{3}}\right]e^{-\gamma r_{a}}$
$t^{\prime}>b\ (\gamma\neq\lambda)$ |  $\frac{z}{\beta}\left[\frac{1}{(z^{2}+\gamma^{2})}-\frac{(\gamma+3\lambda)}{(\gamma+\lambda)^{3}}\right]e^{-\gamma r_{a}}-\frac{z}{\beta}\left[\frac{1}{(z^{2}+\gamma^{2})}-\frac{(\gamma-3\lambda)}{(\gamma-\lambda)^{3}}\right]e^{-\gamma r_{b}}$ $\quad\ +\frac{z}{\beta}\left[\frac{\lambda r_{b}^{2}}{(\gamma-\lambda)}+\frac{(\gamma-3\lambda)r_{b}}{(\gamma-\lambda)^{2}}-\frac{\gamma-3\lambda}{(\gamma-\lambda)^{3}}\right]e^{-\lambda r_{b}}$
$t^{\prime}>b\ (\gamma=\lambda)$ | $\frac{z}{\beta}\left[\frac{1}{(z^{2}+\lambda^{2})}-\frac{1}{2\lambda^{2}}\right]e^{-\lambda r_{a}}-\frac{z}{\beta}\left[\frac{1}{(z^{2}+\lambda^{2})}-\frac{(2\lambda r_{b}+3)r_{b}^{2}}{6}\right]e^{-\lambda r_{b}}$

##  11 Experimental Details

All models in the experimental section are implemented using GPyTorch [10],
trained by Adam Optimizer on an NVIDIA A100-SXM4 GPU (for TIMIT and UCI
datasets) or Apple Macbook CPUs (for illustrative examples), with a learning
rate of 0.01 and a batch size of 10,000. The models using doubly stochastic
variational inference, e.g., IDDGPs, DLFM-VFRF, employ five samples for layer-
wise sampling during training. We follow [29] to set up a linear mean function
for all the inner layers and a zero-mean function for the outer layer to avoid
pathological behaviour [8]. The weights of the linear mean function are fixed
and determined by SVD if the input and output dimensions are not equal. The
variational distributions over inducing variables are initialised to normal
distributions with zero mean and variances identity for the outer layers and
$10^{-5}$ for the inner layers. The inducing points are initialised with
K-means. All models, including RFF-based models, used 100 Monte Carlo samples
for evaluations on test data.

Unless specifically stated, the RKHS interval is set to $[a,b]=[-1,4]$, and
all input data are normalised to $[0,3]$. We initialise our model with length-
scale $l=0.1$ for the TIMIT dataset and $l=1$ for the UCI datasets, ODE
parameters $\alpha=1,\beta=0.01$, kernel variance $\sigma^{2}=0.1$ and noise
variance $\varepsilon^{2}=0.01$.

##### Intermediate Outputs for Synthetic Data

We present here the posterior distributions of the DGP, the IDDGP and our
model DLFM-VFRF on the synthetic data in Section 5.

Figure 9: Comparison of posterior distributions of different compositional GP
models on synthetic data. Top row: The output distributions of the
intermediate layers. Bottom row: The posterior predictive distributions.
Training points are marked with red dots, posterior means with blue lines, and
uncertainty with varying shades of blue. Each panel depicts 20 samples from
the posterior distribution. Although both the IDDGP and our model show better
fitting to the multi-step function, they have very different intermediate
posterior distributions.