Deep Random Features for Scalable Interpolation of Spatiotemporal Data
Authors: Weibin Chen
Azhir Mahmood
Michel Tsamados
So Takao
Venue: ICLR 2025 Poster (2025)
Keywords: Random Features
Deep Gaussian Processes
Bayesian Deep Learning
Remote Sensing
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We propose a scalable Bayesian deep learning framework to interpolate remote sensing data for increased accuracy and flexibility.

Abstract: The rapid growth of earth observation systems calls for a scalable approach to interpolate remote-sensing observations. These methods in principle, should acquire more information about the observed field as data grows. Gaussian processes (GPs) are candidate model choices for interpolation. However, due to their poor scalability, they usually rely on inducing points for inference, which restricts their expressivity. Moreover, commonly imposed assumptions such as stationarity prevents them from capturing complex patterns in the data. While deep GPs can overcome this issue, training and making inference with them are difficult, again requiring crude approximations via inducing points. In this work, we instead approach the problem through Bayesian deep learning, where spatiotemporal fields are represented by deep neural networks, whose layers share the inductive bias of stationary GPs on the plane/sphere via random feature expansions. This allows one to (1) capture high frequency patterns in the data, and (2) use mini-batched gradient descent for large scale training. We experiment on various remote sensing data at local/global scales, showing that our approach produce competitive or superior results to existing methods, with well-calibrated uncertainties.

Published as a conference paper at ICLR 2025
DEEP RANDOM FEATURES FOR SCALABLE INTERPOLATION OF SPATIOTEMPORAL DATA
Weibin Chen', Azhir Mahmood 1. Michel Tsamados' & So Takao
Universityy College London, UK 2PhysicsX, UK
CCliformia Institute of Technology, USA
ABSTRACT
The rapid growth of earth observation systems calls for a scalable approach interpolate remote-sensing observations. These methods in principle, should ac- quire more information about the observed field data grows. Gaussian processes (GPs) are candidate model choices for interpolation. However, due to their poor scalability, they usually rely on inducing points for inference, which restricts their expressivity. Moreover, commonly imposed assumptions such as stationarity pre- vents them from capturing complex patterns the data. While deep GPs can overcome this issue, training and making inference with them are difficult, again requiring crude approximations via inducing points. In this work, we instead ap- proach the problem through Bayesian deep learning, where spatiotemporal fields are represented by deep neural networks, whose layers share the inductive bias stationary GPs on the plane/sphere via random feature expansions. This allows one to (1) capture high frequency patterns in the data, and (2) use mini-batched gradient descent for large scale training. We experiment
various remote sens- ing data at local/global scales, showing that our approach produce competitive or superior results to existing methods, with well-calibrated uncertainties.
INTRODUCTION
The advent of earth observation systems have made it possible to monitor virtually all of earth's 's atmosphere and the ocean at unprecedented scales. This development has been pivotal to the un- derstanding of anthropogenic impact on the environment, including global warming and rise in sea level. Hence, it is crucial that we are able to process the voluminous data effectively and extract maximal information from it to make better informed decisions in our path to achieving sustainable development goals.
However, observations from satellite products inherently sparse in space-time, requiring methods to effectively fill in the gap at unobserved locations (Le Traon et al.) 1998). This typically relies on data assimilation techniques such as the ensemble Kalman filter (Evensen 2003). which requires one to have access to a physical model that describes the evolution of the field. While this can produce detailed and accurate reconstructions of the field, the physical models are typically expensive to run at high resolutions, often requiring access to high performance compute clusters. This can be challenging when one does have the expertise nor the resources gain access and/or run the models. On the other hand, statistical methods such as Gaussian process regression (GPR, Williams & Rasmussen] (2006)) can be deployed. However, GPR scales poorly to large data sets, necessitating approximate inference schemes such sparse Gaussian processes (Titsias 2009). which may result crude approximations if the underlying process does not have sufficiently large lengthscale or smoothness (Burt et al. 2019). Moreover, kernels used for GPR are often too simplistic, which can prevent learning of detailed fluctuations in the underlying non-stationary and multi-scale field. Deep Gaussian processes (DGPs) (Damianou & Lawrence 2013) have emerged as attractive solution to the latter problem. However, they still suffer from the difficulty of computing the posterior, again requiring variational inference learn only crude approximation to the true posterior.
In recent years, Bayesian deep learning (BDL) have emerged as an alternative paradigm for sta- tistical modelling, which combines the flexibility and scalability of deep learning methods with

---

Published as a conference paper at ICLR 2025
SVGP
DRF
0.2
0.1
-0.1
0.00
0.04
0.0
Figure 1: We propose deep random features (DRF) for accurate and flexible interpolation of satellite measurements of the earth's surface (Left). Compared to sparse variational GPs (SVGP, Centre), an ensemble of DRFs is able to achieve more detailed reconstructions of the field with sensible uncertainty estimates (Right).
Bayesian modelling principles (Papamarkou et al. 2024). In our current setting, we can approach spatiotemporal interpolation using BDL, by representing the ground truth underlying field fl by Bayesian neural network fo X' - V, (here, X denotes the spatiotemporal input space and y the output signals) and training input-output pairs D {(z.. N for In E X' and yn V, corresponding to earth observations. However, naïve design choices for fo can lead to poor recon- structions of ft; for example, vanilla deep ReLUU network is bound perform poorly it fails to learn high frequency features Tancik et al. (202)) On the other hand, deep neural networks with trigonometric activations (Sitzmann et al. 2020 Lu & Shafto (022) have emerged as an effective model for representing high-frequency spatiotemporal signals. However, they are not designed for interpolation of sparse data mind and are therefore prone to overfitting.
Taken altogether, we propose to design fa inspired by DGPs, such that it retains the learning capac- ity of DNNs, while having the interpretability and inductive biases of GPs. Our main contributions are follows: We propose the use of kernel-derived random features (Rahimi & Recht, 2007) as building blocks for BNNs to model spatiotemporal fields. We demonstrate through extensive exper- iments that they are capable of capturing fine-scale information in data, while being able to quantify uncertainty accurately by considering deep ensembles. Furthermore, motivated by recent develop- ments in geometric probabilistic modelling (Borovitskiy et al. 2020), also consider analogous random features on the sphere, leading to a novel DNN architecture with Gegenbauer polynomial activation functions that model global weather fields that adapted to the sphere. Our models are easily implementable in modern deep learning frameworks such as PyTorch and scale up to large datasets exceeding millions of data points through mini-batched gradient-based optimisation, pushing the boundary of what is currently possible with statistical interpolation.
1.1 RELATED WORKS
In Cutajar et al. 2017). trigonometric feature expansion of DGPs similar to ours have been consid- ered, with the intent of proposing a tractable variational inference (VT) scheme for DGPs. They show superior performance to mean-field VI (Damianou & Lawrence] 2013), however, have been largely overlooked due part to the adoption of doubly stochastic VI (Sllimeeni Deisenroth (017). as
the de facto standard method for DGP inference. [Jiang al. (2024) similarly leverages random
Fourier features to approximate kernel machines, emphasising composite kernels for incorporating prior knowledge into neural networks. DNNs with trigonometric activations have resurfaced as object of interest more recently, with the emergence of neural radiance fields (Mildenhall et al. 202TT
and subsequent work on implicit neural representations Tancik et al. 20200 Sitzmann et al. 2020
2

---

Published as a conference paper at ICLR 2025
Rigorous study of trigonometric networks and their connection to DGPs have been considered in Lu & Shafto (202)). and more general investigation of wide DNNs with bottlenecks in relation to DGPs have been considered in Agrawal et al. (2020); Pleiss & Cunningham (02I). Other closely related works include Meronen et al. (2020 2021). who study calibration of shallow networks with periodic activations, Garnelo al. [2018 proposes different approach combining aspects of GPs with
DNNs, and the works Sun et al. 2020);; Dutordoir et al. (2021) establish connections between neural
network layers and inducing points for GPs/DGPs.
2
BACKGROUND
2.1
GAUSSIAN PROCESSES AND DEEP GAUSSIAN PROCESSES
A Gaussian process (GP) is a random function f R' R such that for any N 0 and any set of points In E IR' for = I,..., NN, we have that (f(.))..... f(xx)) E R.V is Gaus- sian. GPs are characterised by a mean function m IR. + R and a kernel such that E [(x))
R x R - R.
m(x) and Cov (ax), ((x') = k(z, for all æ, x' E R' Willlams &
Rasmussen 2006 Extending these to have vector outputs f R' - R is made possible by
considering vector-valued means m TRd + RD and matrix-valued kemels k :R' x TR. T RxOxO satisfying E [f](a]]] (x) and Cov [f(( fi(x')] = kay(a x²), Vi = 1,...,0 O. We write GP m, k) to denote that f is GP with mean and kernel k. A deep GP (DGP) R extends GPs by considering compositions f(x) = f" a a.... fx)) where f f" R " R= for f 2.......
f R f : RR.
and f : The intermediate states RB are referred to as the bottleneck. We note that DGPs are more flex- ible class of models than GPs. However, due their compositional structure, they are longer Gaussian and therefore require approximate methods for inference, e.g. using variational Bayes.
R are vector-GPs.
2.2 RANDOM FOURIER FEATURES
Consider a zero-mean scalar GP f GP(O, k) for some kernel k. We say that ks is stationary if there exists a function k IR - IR such that k(rr x') = <(x - x²). In Rahimi & Recht 2007). it is shown that any stationary kernel on R. can be expressed as an expectation k(x, x²)
Ew.D [cos(w] x + b) cos(w> T x' +)]
(1)
20
cos( ba) cos(w, x' bn), Wh. p(w), by U([0 2x)
for some a > 0, where p(w) is the normalised Fourier transform of the function and U[0, 2m)) denotes the uniform distribution in the interval [o, 2n). From the weight-space viewpoint of GPs, equation 2 implies that we have
H
f(x)) Σ (n))(r), On ~ N(0, 1),
(3)
h=1
where (a(x)) 20 /H cos(w)) bu), 1= H,
(4)
with w/ ~ p(w) and by ~ U(0 2x. Extension to vector-valued GPs f R RO with inddpne- dent output components is straightforward, leading to random Fourier feature representation of the form f(x) @Φ(#) for O E goxH with O.j ~N 1), i.i.d. Vi, j. For details and examples of random Fourier features, we refer the readers to Appendix
DEEP RANDOM FEATURES FOR SPATIOTEMPORAL MODELLING
GPss are commonly used in spatiotemporal modelling due to their interpretability and smoothness inductive biases that are appealing to many geostatistical applications (Wikle et al. 2019). However, they are limited by their poor scalability and Gaussian assumptions. On the other hand, deep neural networks (DNN) offer a scalable, flexible modelling framework, however, do not have the desirable inductive bias of GPs. Motivated by this, we consider deep random features (Figure 2). which use random features corresponding to stationary GPs as building blocks for a larger neural network model tailored for spatiotemporal modelling, combining the benefits of both GPs and DNNs.
3

---

Published as a conference paper at ICLR 2025
Observations
Ensembles
Loss
Deep Random Features
Input coordinates
Spatial
Output
Temporal
Ski connections
(+) oos(uz 8), x nu p(s), 07
Figure 2: Illustration of spatiotemporal modelling with deep random features.
3.1 DEEP RANDOM FEATURES
we have seen that shallow GP can be approximated by a combination random features according to equation 3 In similar fashion, Cutajar et al. (2017) propose random feature expansion deep GPs, by replacing each GP layer by their corresponding random features, yielding DNN architecture that mimic the behaviour of the original DGP. In further details, let Φb R H be a random feature (equation which may be viewed equivalently as a single layer of a neural network with weights {wim H - blases {bm, 1 /-m= and cosine activation. The first hidden layer in deep GP, which itself is a vector-valued GP ft R + RB, can be approximated by linear model
In Section 2.2
h h'(x)) θ (1)'
(5)
where Q E BxH with 0 ~ NN 1), I,...,,
- B,j H. Similarly, given random
features o' RD + R for - 2.... L, Gaussian weights O E R axff for l ...... L and @ R we may consider DNN fe R' Eº of the form
fe = h- ..... h² o (x), where h(()) = O O')((x), 1= L,
(6)
which we refer to as the random feature expansion
DGP f R R Generally, we may consider building DNNs independently of a DGP by using layers of the form h(x) +Φ(x)) (x) as building blocks for a neural network. We refer to such models as deep random features. Note that while each layer h (æz) approximates stationary GP, its compositions
longer stationary.
3.1.1 TRAINING
In contrast to standard neural networks, when we train deep random features, we opt to alternate between trainable and fixed layers, where the parameters wf bf, in the layers Φe(-) are fixed upon initialisation, but the parameters
:= (e "...... @²} are trained. This is to mimic training of DGPs from the weight-space perspective, where inference should only be made with respect θ. Given a dataset D = {(X,, yn) ) ) and an arbitrary loss f. TRO x R > IR, we minimise
Lerain (0:
1 N
Σ
( n) vnd 0113
for some regularisation parameter 3 > From generalised Bayes' perspective, this is equivalent to maximum a priori estimation with the generalised posterior p( (D) x exp(-é(fe) (X), y)) ( (issiri et al. 2016). Using mean-squared eтror as the loss and considering a shallow network, minimising equation via gradient descent be neural tangent kernel limit (Lee et al.) 2019 He et al. 2020). Using other losses such as the Huber loss, this becomes akin to robust GP regression (Algikar & Mili 2023 Altamirano et al. 2024
sampling from the GP posterior in the
4

---

Published as a conference paper at ICLR 2025
3.1.2 SPHERICAL RANDOM FEATURES
When modelling signals over the sphere, which arises when need interpolate global satellite measurements (see Figure MT we require an analogous notion of random features defined over the sphere. In Borovitskky et al. (2020). commonly used kernels such the Matérn kernels are extended to be defined over general Riemannian manifolds, including the two-sphere S?. In general, such kernels can be approximated by the Mercer sum
(ssss)
(()))1((8)9)(6)) s, s².
C$ j=0
for some 4 : RR--- + R and constant Cq determined from the kernel, and {yy)j=o- {p)) the top eigenvalues and eigenfunctions respectively of the negative Laplace- Belttamm operator; on s2, the latter is precisely the spherical harmonics. Furthermore. on S", by making use of the addition theorem for spherical harmonics and the result (Azangulov et al., 2024, Proposition 7), we get an alternative expression for the kernel (see Appendix[AA. [.2 for the derivation)
k(ss s') zs Eui [< (dgz (s) (dgx (d, E.S'
(9) where (.) are the Gegenbauer polynomials of order n and weight parameter a, dga (:.)) denotes the geodesic distance on S². Cw is appropriate scaling constant, and the ex- pectation is taken over b ~ USS2), the uniform distribution over the sphere, and w ~ ( ((A))...., CT (()))) Then, by considering Monte Carlo approximation of the expectation in equation this implies random feature maps of the form
m m. M, where Wm ~ Multinomial ((A))....... " ((A)), bm ~ ~S
Pe (s) V M- (cmm 22/2 (s, (m,b,)), $ E
(10)
(11)
This gives us an analogous notion of random features (equation 4 on the sphere, which we can use as a component in our deep random feature model when our input is spherical.
Remark We may also consider the deterministic features (s)= C = (D (s), derived
=
from equation g which is analogous the regular Fourier features (Hensman et al. 2018 Solin &
Särkkä 2020 in the planar case. However in practice, we find that working with random features (equation 10p produce more stable results when using single precision arithmetic.
3.1.3 SPATIOTEMPORAL MODELLING WITH DEEP RANDOM FEATURES
So far, we have only discussed how to process spatial inputs. In order to deal with the temporal components in our data, we first consider deep random features in the spatial domain h!!"²) X + R (X' = Red or 8²), and temporal domain h, (L,) R - BB separately, before combining them as
f(x,1) @ (h (x)) h (()))
(12)
where θ ROx2B are learnable weights initialised with i.i.d. standard Gaussians. At short timescales, geospatial fields are approximately stationary, hence we
single layer network model the temporal component h((")) (i.e.. L, = 1). To introduce more complex spatiotempo- ral dependence, we can replace the linear output layer in equation 12 with deep random features.
However we find that in most applications this is unnecessary, only introducing extra cost.
3.1.4 SKIP CONNECTIONS
To prevent pathological behaviour from emerging as increase the network depth, add skip connections to the inputs (Duvenaud et al. 2014 Dunlop et al.) 2018). In the planar case, this takes
((+))()) O ++1
(concat (9 (x),
(13)
in the (e + 1)-th layer of the network, where + BB+ RM In the spherical case, this is not straightforward as the outputs of each layer will be Euclidean while the input is spherical. To this end, consider tf+1 to be additive random features corresponding to the sum kernel k((x, s), (x',s') = kgg (x, x²) + kga (s, for æ, x' E R s E S², where Kzg (., kga ().)) stationary kernels their respective spaces (see Appendix B.T for details).
5

---

Published as a conference paper at ICLR 2025
3.2 UNCERTAINTY QUANTIFICATION
Uncertainty quantification (UQ) with deep random features achieved using standard Bayesian deep learning techniques. In particular, consider the following methods in our experiments:
Variational inference. This considers Gaussian approximation to the posterior p() D) ~ ( m. C). Here, the moments of the variational distribution (()) =N(OOlm C) are learned by maximising the evidence lower bound (ELBO)
LELBO = Ey (fl - (C(g))o))
(14)
where Pe ()) =N I) is the prior 6 and (CC(-)):) is the Kullback-Leibler divergence.
Dropout at test time. A simple heuristic for obtaining uncertainty estimates is to apply dropout not only training time but also at test time. This yields an ensemble of random outputs, whose empirical distribution informs us of the model uncertainty. In fact, one can understand this as form of variational inference, as shown in Gal & Ghahramani (2016).
Deep ensembles. Deep ensembles [Lakshminaraymaan al. 2017) obtain uncertainty estimates by training an ensemble of models, initialised from different random seeds. The ensemble of outputs then used estimate uncertainty, similar to the dropout method for UQ. While being simple method, this has been shown to be surprisingly effective at obtaining uncertainty estimates. More- over, provided the model is small enough (which is often the case for deep random features), the ensembles can be trained in parallel on a single GPU.
3.3 HYPERPARAMETER SELECTION
Our deep random features model contain several hyperparameters X, including those of the kernel (e.g. lengthscales) that use to derive random features. If variational inference is used for UQ. we can take the ELBO (equation (4] for model comparison, as it may be viewed as a surro- gate for the log model evidence log p A[)), being its lower bound. When using ensemble based methods, rely performance hyperparameters. In particular, we consider the following validation loss to select X
held-out validation set D* { ";))"", select our
C (A)
(15)
N
where {f are the ensembles. Note that this can be viewed as approximating the negative log-predictive density (see Remark 3 Appendix I.II) We minimise this loss using Bayesian opti- misation. In practice, we find that learning À from equation alone may still lead to overfitting models. Therefore to prevent this, may opt to add an extra functional regularisation term
ト
Lealtrog (A) (1 a) (A) + V Vf(-:
(16)
for a E [0. 1), which helps to penalise those that give rise to functions f with sharp gradients (see Remark 2 in Appendix B.2 why we do not consider hyperpriors for regularisation). Here, 11; denotes the appropriate - norm depending the input space and V the gradient. For spherical inputs, we refer the readers to Appendix B.2 for more details.
EXPERIMENTS
We evaluate the spatiotemporal deep random features (DRF) model on various remote sensing datasets and compare against various baselines assess its ability to make predictions and quantify uncertainty. In our first experiment, we consider interpolation model's 's ability to recover the ground truth. In our second and third experiments, we consider inter- polation of real satellite data at local and global scales to test the robustness our method. Details can be found in Appendix C All experiments are performed using the NVIDIA L4 GPU.
synthetic data, and evaluate our

---

Published as a conference paper at ICLR 2025
4.1 BASELINE MODELS
Throughout this section, we consider several models as baselines to compare our model against. We consider both GP-based baselines and DNN-based baselines. In the former category, we consider the sparse variational GP (SVGP) model in Hensman et al.] (2013), deep GPs (DGP) using doubly stochastic variational inference (Salimbeni & Deisenroth (2117) and a mixture model of local GPs using the Sat library (Gregory et al.) 2024b). In the latter caregory, we consider deep ensembles of multilayer perceptrons (MLP) with ReLU activations, Fourier features network (FFN, [Tancik (FFN,fTancik al.) (020)), which use random Fourier features in the first layer only, MLP with sinusoidal activations (SIREN, Sitzmann et al. )0200) and conditional neural processes Garnelo et al. (2018).
4.2 EVALUATION ON SYNTHETIC DATA
The purpose of our first experiment is use synthetic observations from ground truth field to evaluate our model's ability to reconstruct the field. We use mean sea surface height (MSS) in the arctic as our ground truth, synthesised from 12 years of altimetry readings of Sentinel-3A, 3B (S3A, 3B) and CryoSat-2 (CS2) satellites. We then generate artificial measurements along S3A, 3B and CS2 tracks between the dates March 1st-10th 2020, taking the MSS values along the tracks and adding i.i.d. Gaussian noise to mimic measurement noise. Our final dataset comprise 1,158,500 505 datapoints; we select 80% of these randomly for training and the remaining 20% for validation. We train all models using the mean-squared error loss (for GP baselines, this corresponds to a Gaussian likelihood) with fixed weight decay parameter matching the observation noise variance. Visual comparison of predictions from all models can be found in Appendix C.5.I
4.2. EFFECT OF DEPTH
In Figure we show the effect of depth on our model's ' ability to reconstruct the true MSS field (in terms of RMSE) and corresponding computation time of the entire workflow, including the time to tune the kernel hyperparameters (see Appendix C.2.2). Generally, we find that deep networks outperform the shallow network on the RMSE, with the four layer model performing the best on this example. With layers, we start to some overfitting, which explains the higher RMSE for the 10 and 20 layer models. In Figure 4. we display mean results for models with two and four layers. We that the deeper network is able to capture higher frequency details, resulting in the improved RMSE. The time it takes to train and tune models with 1-4 layers are not significantly different.
RMSE and Time vs Number of Layers
2-layer DRF
4-layer DRF
0.18
RMSE
Time (minutes)
0.4
45
0.17
40
Time (minutes)
0.2
RMSE
0.16
0.0
0.15
30
-0.2
25
0.14
-0.4
0.13
1234
10
20
Figure 4: Comparison of predictions from DRF with two layers (left) and four layers (right). The four layer model is able to capture finer details compared to the two layer model.
Number of Layers
Figure 3: Comparison of RMSE and computation time vs. number of layers.
4.2.2 UQ COMPARISONS
Next, we compare various UQ methods applied to DRF model with four layers. In Table we dis- play comparisons with respect to the root mean squared error (RMSE), the negative lig-likelihood (NLL), and the continuous ranked probability score (CRPS) (see Appendix C. for details our evaluation method). The latter two evaluate the quality of uncertainties produced. Overall, we find that deep ensembles produce the best result in the CRPS and the RMSE, whereas variational infer- ence (VI) yielded the best NLL. The lower NLL using VI may be due to the fact that its predictions typically underconfident (see Figure 8 Appendix C.5.1 and NLL penalises them more lightly
7

---

Published as a conference paper at ICLR 2025
Model
CRPS
NLL
RMSE
Time (minutes)
DRF (Ensembles)
0.046 + 0.005
13.590 + 4.899
0.135 0.006
19.7 0.4
DRF (VI)
0.071 + 0.019
-0.407 + 0.756
0.166 + 0.021
6.40 + 0.06
DRF (Dropout)
0.174 + 0.001
425.987 + 208.969
0.238 + 0.001
48.6 + 2.5
SVGP
0.230 + 0.001
320.811 + 52.960
0.155 + 0.002
14.6 + 0.0
DGP
0.058 + 0.001
1614.069 #32..517
0.135 0.002
42.3
GPSat
0.045 0.007
74.738 + 15.622
0.126 = 0.001
6..6
ReLU MLP
0.062 0.000
7.877
0.146 0.000
10.07
FFN
0.072 + 0.008
126.869 + 111.178
0.153 + 0.005
31.3
SIREN
0.066 + 0.000
1 0.393
0.155 - 0.000
2.55 = 0.002
CNP
0.238 = 0.070
2.525 + 0.459
0.202 = 0.010
21.0 + 0.1
Table 1: Comparison of the CRPS, NLL and RMSE scores for a four-layer DRF (with different UQ methods) against various baselines the synthetic experiment. Best performing model in bold, second best performing in blue and third best performing in orange. We display the mean and standard deviation over five experiments.
than overconfident ones. However, the results in Figure suggest that the results from deep en- sembles are better calibrated to the observations, which explains the lower CRPS. Dropout does not perform well in neither the mean prediction nor uncertainty estimation.
4.2.3 BASELINE COMPARISONS
In T we also display results for the other baseline models described in Section 4.1 Regarding computation times, for DRF and FFN, we include the time tune the kernel hyperparameters using Bayesian optimisation (Section 3.3 to make fair comparison with the GP-based baselines, where recorded. However,
the total time for training and inference such as number of layers and hidden units to be fixed (L = 4, B = 128, H = 1000). For the other DNNebased baselines, we assume the architecture tuned ahead of time and fixed.
assume other hyperparameters.
Comparing with the GP baselines, we find that DGP and the GPSat model to be closest competitors to the DRF deep ensembles, with GPSat surpassing its performance on the RMSE. However, the time taken to train the DGP and CPSat model are two to three times longer than the time taken to train and tune the ensemble DRE For example, Ga trains 1225 local GP models on this example, which makes computation heavy. DGP has overall low predictive variances (see Figure 9 Appendix C. which results in high NLL values. In contrast, the uncertainty estimates of DRF and GP Sat well-calibrated to the satellite tracks. Qualitatively, all three models recover the ground truth well, with GPSat and DRF reconstructing it almost perfectly.
Comparing to other DNN basclines, SIREN's performance is noteworthy, being similar to DRF in that it uses trigonometric activations and differing only in the way the weights are initialised and whether it has bottleneck layers with fixed preactivations. The DRF ensemble is better able to capture the spatiotemporal patterns of the field, reinforcing the importance of the subtle architectural differences. The ReLU network, FFN and CNP all lead to oversmoothing results, similar SVGP.
4.3 FREEBOARD ESTIMATION FROM SATELLITE ALTIMETRY DATA
Here, we consider the interpolation of real altimeter readings of sea-ice freeboard taken along the Sentinel-3A, 3B and CryoSat-2 satellites (Gregory et al. 2024a). Real satellite altimetry measure- ments are typically noisy with heavy-tailed statistics (see Figure 7 Appendix C.3). hence they present more challenging setting than previous synthetic experiment. Our goal here is to test the robustness of our approach in comparison to other methods in this more realistic setting. Ex- perimental details can be found in Appendix and visual comparison of all results can be found in Figure 10, Appendix C.5.2
We compare two-layer DRF against the same baselines as before and display the root mean abso- lute error (RMAE), CRPS and negative log predictive density (NLPD) on a separately held out test data comprising 15% of the entire data in Table 2 We use the RMAE instead of the RMSE here
8

---

Published as a conference paper at ICLR 2025
as it is more robust to the heavy-tailed statistics of the measurement error, and therefore provides more reliable performance metric in this setting. The hyperparameter search for DRF was performed
with functional regularisation (see Section 3.3) we found that without it, optimising only the validation loss lead to overfitting models (see Figure Here, we used penalty weight of a = 0.9.
We find that out of our GP-based baselines, SVGP and the GPSat model give comparable perfor- mance to DRF ensembles, with CP Sat giving the best performance quantitatively. However, when examine the outputs from all three models (Figure(5), we find that the results from GPSat con- tain spurious patterns resulting from unstable hyperparameter optimisation at several local expert locations. This issue occurs since the GPSat local experts only see data in local regions, making them more sensitive to the heavy noise present in the data. On the other hand, DRF sees data glob- ally, which helps them to identify the larger structures the data, while simultaneously capturing the finer details owing to their deep architecture. We find that qualitatively, SVGP performs well on this example, due to the larger prominence of low frequency features in the underlying field that extend across the basin. We see that DRF provides middle ground between the two, being neither "too local" as we see in GP Sat, nor "too global", demonstrating its ability to adapt to the character- istics of the field. The other DNN baselines are found to perform poorly, with the ReLU MLP and CNP showing especially poor fit. The quantitative metrics for FFN and SIREN are actually decent, with FFN performing best RMAE. However, the qualitative results in Figure TO show signs of heavy overfitting with both models, which is especially prominent in SIREN.
Model
CRPS
NLPD
RMAE
DRF (Ensembles)
0.077 + 0.000
-0.944 + 0.020
0.322 + 0.000
SVGP
0.079 + 0.001
-1300 + 0.004
0.322 + 0.001
DGP
0.208 + 0.000
-0.159 + 0.000
0.339 + 0.001
GPSat
0.076 + 0.001
-1.167 0.025
0.318 + 0.000
ReLU MLP
0.076 1.578
0.751 + 0.554
FFN
0.080 0.001
2.145 = 1.811
0.316 + 0.004
SIREN
0.088 = 0.000
-1.009 0.010
0.345 0.001
CNP
0.101 + 0.001
-0.898 + 0.030
0.328 + 0.002
Table 2: Comparison of the CRPS, NLPD and RMAE scores for a two-layer DRF against various baselines sea-ice freeboard interpolation from S3A, 3B and CS2 satellite altimetry readings. Best performing model in bold, second best performing in blue and third best performing in orange.
DRF reg)
DRF(reg)
GPSat
SVGP
0.3
0.2
0.1
0.0
-0.1
Figure 5: Mean results of DRF, GPSat and SVGP freeboard interpolation. For DRF, plot results obtained both with and without functional regularisation during hyperparameter search.
4.4 LARGE SCALE INTERPOLATION OF GLOBAL SEA LEVEL ANOMALY
In this final experiment, we investigate the potential of DRF to interpolate global fields using spheri- cal random features (Section 3.1iiinn3.1... 3.2 in the spatial inputs. For this experiment, we use real data of sea level anomaly measurements collected from the Sentinel 3A, 3B satellites (Copernicus Data Space Ecosystem 2024). By considering four days of measurements, our final data consists of 8,094,569 datapoints. We use 80% for training, and 20% for validation. Similar to our previous data ob- tained from real satellite measurements, this data contains many outliers, making it a challenge to interpolate the data robustly, let alone whilst being consistent with the geometry of the sphere.

---

Published as a conference paper at ICLR 2025
Our goal is to fit spatiotemporal field f & x R R. To this end, we consider DRF model whose first layer in the spatial component is given by the spherical random feature Øga : S² + R (equation TO). The subsequent layers are given by Euclidean random feautures. To train model, we opted to use the Huher loss instead of MSE, which gave rise to slightly more robust results, likely due to the large number of outliers. We also used functional regularisation with penalty weight of = 0.95 when tuning hyperparameters.
In Figure we compare the mean predictions of the spherical DRF model with predictions from (1) SVGP using the spherical Matérn kernel of [Borovitskiy et al.] (2020), and (2) the Euclidean DRF model, taking the longitude and latitude coordinates of the satellite tracks as spatial puts in R²2 We use the spherical Matérn kernel implementation in the geometric-keznel package Mostowsky et al. (2024) to model the spatial component of our SVGP baseline. The temporal component is included by modelling the GP with a product kernel (((2r, f), (x',t') kga (2, (t,) Comparing the SVGP output with DRF, we see that they are both able to capture the larger patterns in the data. However, SVGP fails to capture some of the finer fluctuations (as also indicated by quantitative metries in Table g in Appendix C.4.)). for instance those around the Antarctic circumpolar current, known for its intense ocean activities.
For the Euclidean DRF, while it admits a deep structure, we find that it is not flexible enough to adapt to the spherical geometry of the imput space. For example, there spurious distortions around the poles cause by the stereographic projection, Figure [ in Appendix (C..3). Perhaps more interestingly, the Euclidean DRF is not able to learn the fine scale fluctuations that the spherical DRF is able to pick up, only being able to learn large scale trends in the data, similar SVGP, This highlights the importance of explicitly incorporating the spherical inductive bias into the model when modelling global fields.
addition to a discontinuity at longitude = 0 (see
Spherical SVGP
Euclidean DRF
Spherical DRF
0.3
0.2
0.1
-0.1
aoo
-0.2
-0.3
Figure 6: Mean results for global sea level anomaly interpolation. From left to right: SVGP using the spherical Matérn kernel, Euclidean DRF, and Spherical DRF. Spherical DRF able learn more intricate details compared to the other two baselines.
CONCLUSION
In this paper, we propose to model spatiotemporal fields using deep neural networks, whose layers derived from random feature expansions of stationary kernels. This neural representation can be trained on observations to effectively fill in the gaps between remote sensing observations of the earth's surface. Our experiments on various remote sensing data demonstrate that the deep ensemble model is able to flexibly adapt to the data, being able to learn both low and high-frequency structures that exist in the underlying field. current limitation of our approach the difficulty of tuning kernel hyperparameters: we use Bayesian optimisation (BO) on the validation loss to achieve this, which require knowledge of the ranges each hyperparameter may take. This is not clear due to the deep architecture, making the hyperparameters less interpretable than the shallow case. Additionally, we observe that it is sometimes necessary to add functional regularisation to reduce BO variance, necessitating hand tuning of the penalty weight a, hyper-hyperpe hyperrhhprrrrremtter Despite this, our promising results suggest the potential for deep learning methods to pave the way for more accurate and flexible reconstructions of spatiotemporal fields from remote sensing data.
10

---

Published as a conference paper at ICLR 2025
ACKNOWLEDGMENTS
WC and MT acknowledge support from ESA (Clev2er: CRISTAL LE procEssor prototype and R&D); MT acknowledges support from PESSA/AD/1-113//1//N//A #ESS//A0/-1006// 1/99//E SIN' XS: Sea Ice and Iceberg and Sea-ice Thickness Products cnter-comparison Exercise) and NERC (#NE/T000546/1 761 and #NE/X004643/1), AM is supported by the EPSRC-funded UCL Centre for Doctoral Training in Intelligent, Integrated Imaging in Healthcare (i4health) (EP/S021930/1) ST is supported by Department of Defense Vannevar Bush Faculty Fellowship held by Prof. An- drew Stuart, and by the SciAl Center, funded by the Office of Naval Research (ONR), under Grant Number N00014-23-1-2722 MSS data, William Gregory for support on the GPSat baseline experiments, Daniel Augusto de Souza and Mare Peter Deisenroth for useful initial discussions, Viacheslav Borovitskiy for tutoring ST about spherical feature maps, Michalis Michaelides for useful feedbacks on the first draft of our manuscript, and PhysicsX for providing additional GPU support.
We also thank Ronald MacEachern for generating and sharing the
CODE AVAILABILITY
The code for reproducing our experiments is available in the following github repository
httegs///gttta..ccoo//ccoor/r//r eaturee
REFERENCES
Devanshu Agrawal, Theodore Papamarkou, and Jacob Hinkle. Wide neural networks with bottle- necks are deep Gaussian processes. Journal of Machine Learning Research, 2020.
Pooja Algikar and Lamine Mili. Robust Gaussian process regression with Huber likelihood. arXiv preprint arXiv:2301.07858, 2023.
Matias Altamirano, Francois- Xavier Briol, and Jeremias Knoblauch. Robust and conjugate Gaussian process regression. In International Conference on Machine Learning. PMLR, 2024.
Iskander Azangulov, Andrei Smolensky, Alexander Terenin, and Viacheslav Borovitskiy. Stationary kernels and Gaussian processes on Lie groups and their homogeneous spaces I: the compact case. Journal of Machine Learning Research, 2024.
Maximilian Balandat, Brian Karrer, Daniel R. Jiang, Samuel Daulton, Benjamin Letham, An- drew Gordon Wilson, and Eytan Bakshy. BoTorch: A Framework for Efficient Monte-Carle Bayesian Optimization. Advances in Neural Information Processing Systems 33, 2020.
Pier Giovanni Bissiri, Chris Holmes, and Stephen G Walker. A general framework for updating belief distributions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2016.
Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, et al. Matérn Gaussian processes on Riemannian manifolds. In Advances in Neural Information Processing Systems, 2020.
David Burt, Carl Edward Rasmussen, and Mark Van Der Wilk. Rates of convergence for sparse vari- ational Gaussian process regression. In International Conference on Machine Learning. PMLR, 2019.
Copernicus Data Space Ecosystem. Sentinel-3 sea level anomaly data, 2024. URL https:// dataspace. copernicus
Kurt Cutajar, Edwin V Bonilla, Pietro Michiardi, and Maurizio Filippone. Random feature expan- sions for deep Gaussian processes. International Conference on Machine Learning. PMLR, 2017.
Andreas Damianou and Neil D Lawrence. Deep Gaussian processes. In Artificial Intelligence and Statistics. PMLR, 2013.
Matthew M Dunlop, Mark A Girolami, Andrew M Stuart, and Aretha L Teckentrup. How deep are deep Gaussian processes? Journal of Machine Learning Research, 2018.
11

---

Published as a conference paper at ICLR 2025
Vincent Dutordoir, James Hensman, Mark van der Wilk, Carl Henrik Ek, Zoubin Ghahramani, and Nicolas Durrande. Deep neural networks as point estimates for deep Gaussian processes. In Advances in Neural Information Processing Systems, 2021.
David Duvenaud, Oren Rippel, Ryan Adams, and Zoubin Ghahramani. Avoiding pathologies in very deep networks. In Artificial Intelligence and Statistics. PMLR, 2014.
David Duvenaud, Hannes Nickisch, and Carl Rasmussen. Additive Gaussian processes. Ad- vances in Neural Information Processing Systems, 2011.
Geir Evensen. The ensemble Kalman filter: Theoretical formulation and practical implementation. Ocean dynamics, 2003. doi: https://doi.org//00 1007/s 10236-003-000000
Yarin Gal and Zoubin Ghahramani. Dropout as a Bayesian approximation: Representing model uncertainty in deep learning. In International Conference on Machine Learning. PMLR, 2016.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPy- Torch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems, 2018.
Marta Garnelo, Dan Rosenbaum, Christopher Maddison, Tiago Ramalho, David Saxton, Murray Shanahan, Yee Whye Teh, Danilo Rezende, and SM Ali Eslami. Conditional neural processes. In International Conference on Machine Learning. PMLR, 2018.
Tilmann Gneiting and Adrian Raftery. Strictly proper scoring rules, prediction, and estima- tion. Journal of the American Statistical Association, 2007. doi: https://doi.or//10. 198/
016214506000004437
William Gregory, Isobel R Lawrence, and Michel Tsamados. A Bayesian approach towards daily pan-Arctic sea ice freeboard estimates from combined CryoSat-2 and Sentinel-3 satellite obser- vations. The Cryosphere, 2021. doi: https://doi.org/10.. 10./////001//-----/ 2022
William Gregory, Ronald MacEachern, So Takao, Isobel R Lawrence, Carmen Nab, Marc Peter Deisenroth, and Michel Tsamados. Datasets for "Scalable interpolation of satellite altimetry data with probabilistic machine learning", 2024a. URL htt.:/ttt:::/eeod/ee org/ doi/ 10. 52B1/ zenodo. 13218421444
William Gregory, Ronald MacEachern, So Takao, Isobel R Lawrence, Carmen Nab, Mare Pe- Deisenroth, and Michel Tsamados. Scalable interpolation of satellite altimetry data with probabilistic machine learning. Nature Communications, 2024b. doi: https://doi.org/10 1033/ 441467-024-51900-x
Bobby He, Balaji Lakshminarayanan, and Yee Whye Teh. Bayesian deep ensembles via the neural tangent kernel. Advances in Neural Information Processing Systems, 2020.
James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence, 2013.
James Hensman, Nicolas Durrande, and Arno Solin. Variational Fourier features for Gaussian pro- cesses. Journal of Machine Learning Research, 2018
Ziyang Jiang, Tongshu Zheng, Yiling Liu, and David Carlson. Incorporating prior knowledge into neural networks through an implicit composite kernel. Transactions on Machine Learning Re- search, 2024. URL https: httttp::///peneevieeet net wttp://ttoee/o..tll......... Saawwav
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell, Simple and scalable predictive uncertainty estimation using deep ensembles. In Advances in Neural Information Processing Systems, 2017.
PY Le Traon, F Nadal, and N Ducet. An improved mapping method of multisatellite altime- ter data. Journal of Ammosperri and Oceanic Technology, 1998. doi: https:/ddo..or/100.113
12

---

Published as a conference paper at ICLR 2025
Jachoon Lee, Lechao Xiao, Samuel Schoenholz, Yasaman Bahri, Roman Novak, Jascha Sohl- Dickstein, and Jeffrey Pennington. Wide neural networks of any depth evolve as linear models under gradient descent. In Advances in Neural Information Processing Systems, 2019.
Nils Lehmann. Lightning UQ Box, 2024.
URL batp:://Liggttiing--q-bbo
readthedocs io/en/latesta
Chi-Ken Lu and Patrick Shafto. On connecting deep trigonometric networks with deep
Gaussian processes: Covariance, Expressivity, and Neural Tangent Kernel. arXiv preprint
1.Xi::2003.07411, 2022. doi: http://d...rrr11005555550x2....2
Lassi Meronen, Christabella Irwanto, and Arno Solin. Stationary activations for uncertainty calibra- tion in deep learning. In Advances in Neural Information Processing Systems, 2020.
Lassi Meronen, Martin Trapp, and Amn Solin. Periodic activation functions induce stationarity. In Advances in Neural Information Processing Systems, 2021.
Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. NeRF: Representing scenes as neural radiance fields for view synthesis. Communications of the ACM, 2021. doi: https://oi.orr//1.1114 1.55//503500
Peter Mostowsky, Vincent Dutordoir, Iskander Azangulov, Noémie Jaquier, Michael John Hutchin- son, Aditya Ravuri, Leonel Rozo, Alexander Terenin, and Viacheslav Borovitskiy. The Geomet- ricKernels package: Heat and Matérn kernels for geometric learning on manifolds, meshes, and graphs. arXiv:2407.0808% 2024.
Theodore Papamarkou, Maria Skoularidou, Konstantina Palla, Laurence Aitchison, Julyan Arbel, David Dunson, Maurizio Filippone, Vincent Fortuin, Philipp Hennig, José Miguel Hernández- Lobato, al. Position: Bayesian deep learning is needed in the age of large-scale AI. In Interna- tional Conference on Machine Learning, 2024.
Geoff Pleiss and John P Cunningham. The limitations of large width in neural networks: A deep Gaussian process perspective. In Advances in Neural Information Processing Systems, 2021.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, 2007.
Hugh Salimbeni and Mare Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. In Advances in Neural Information Processing Systems, 2017.
Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell, and Gordon Wetzstein. Im- plicit neural representations with periodic activation functions. Advances in Neural Information Processing Systems, 2020.
Henriette Skourup, Sinéad Louise Farrell, Stefan Hendricks, Robert Ricker, Thomas WK Ar- mitage. Andy Ridout, Ole Baltazar Andersen, Christian Haas, and Steven Baker. An assess- ment of state-of-the-art mean sea surface and geoid models of the Arctic ocean: Implica- tions for sea ice freeboard retrieval. Journal of Geophysical Research: Oceans, 2017. doi: https://doi.org/1..oo 10tp//000000///1111113377
Arno Solin and Simo Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. Statistics and Computing, 2020. doi: https://doi.or/// 1007/s 1 1222-019-09886-w
Shengyang Sun, Jiaxin Shi, and Roger Baker Grosse. Neural networks as inter-domain inducing points. In Third Symposium on Advances in Approximate Bayesian Inference, 2020.
Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains. In Advances in Neural Information Pro- cessing Systems, 2020.
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics. PMLR, 2009.
13

---

Published as a conference paper at ICLR 2025
Christopher K Wikle, Andrew Sammit-Mangion and Noel Cressie. Spatio-temporal statistics with R. Chapman and Hall/CRC, 2019.
Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning- MIT press Cambridge, MA, 2006.
A RANDOM FEATURES FOR STATIONARY GAUSSIAN PROCESSES
A.1 RANDOM FOURIER FEATURES FOR PLANAR GAUSSIAN PROCESSES
Let k Rd x Rd + R be a stationary kernel. That is, there exists a function & R + such that k( x²) k(x== x¹). Then Bochner's theorem states that there exists spectral density 8 : Rd > R. such that
<(x x') / s(w)ee (m-w')
(17)
=wwmml) (Sw (co
(18)
where (x)) == سi 0² == -ge sdw) and p(w) = s(w)/o² so that p(w) probability density function. Expanding this further, have
1
o² "wwpp() [cos(w> T Φ)) cos(w x') + sin( x) سin x')l i Bumplw) /cos(w T x)) sin(w" x') + sin xt) cos( x')
(19)
(20)
=0 (since function is odd)
Jo 20 20 (wwn((() [Enmer(0,27] [cos(w To x + b) cos(w> T x' + bb]]
Ewwppi)
cos(w T x + b) cos(w x' + b)db
(21)
π
(22)
20 cos(w by) cos(w o bm),
(23)
N
where wn p(w), bn ~ U10, 2m.. The final expression gives us the random features
(24)
20 M cos(w x + bm))
V
1,..., M,
(25)
corresponding to the stationary kernel kt The only information that changes as we change the kernel is the (normalised) spectral density p(w), used to sample the weights Below, we give examples of such p(w) for the squared-exponential and Matérn kernels.
Example (Squared-exponential kernel) The squared-exponeee kernel is given by
| -
<(α x') = 0² o'eep exp
(26)
where E is the lengthscale hyperparameter and 0² is the kernel variance. The Fourier transform of (up to normalisation constant) can be checked to be given by
gd
||0|| 2
p(w)
(2m)a/2 exp
(27)
Note that this is precisely the probability density function of multivariate Gaussian N(0, 2 Thus, we obtain a random features upproximation to the squared--xponential GP by sampling wm NO, 2 in equation 24
I).
Example 2 (Matérn kernel) The Matérn kernel is given by
<(x x²) a 2'-"
(v20
x
xo
K, (v22
xx
(28)
T(v)
14

---

Published as a conference paper at ICLR 2025
where is the smoothness hyperparameter and K, is the modified Bessel function of the second kind. It is well-known that the Fourier transform of the Matérn kernel equation 28 reads
p(w) (V(V + a) (2V) [(u)
(20 4
(29)
and in fact, this can be identified as the probability density function of the multivariate t-distribution (0, I). Thus, we obtain a random features approximation to the Matérn GP by sampling Wm (0,
1) 21) in equation 24
A.2 MATÉRN GAUSSIAN PROCESSES THE SPHERE
Following Borovitskiy et al. 2020), we define a Matérn Gaussian process on the two-sphere & with lengthscale and smoothness parameter be a stochastic process f defined by the solution the stochastic partial differential equation
2v 224 c -g
Wo.
(30)
Here, V. denotes the space-time white noise process over (S2 with spectral density and Aga denotes the Laplace-Beltrami operator on the two-sphere (see (eeovviiiss et al. (2020) for the notion of solution to this equation). By (Borovitskiy et al. 2020) Theorem 5), the corresponding kernel has an explicit expression of the form (follows from Mercer's theorem)
0² 00 2j+1
2v
(s, s') )
Σ
+2,)
(sA()))sN(ss), Vs, E s?,
(31)
C₂
E2
j=U
where { are the eigenvalues of the positive definite operator - Aga, (pga) 1 100 the spherical harmonics (these are also eigenfunctions of -gga), and
2v
C2
C +
(32)
1=0
is a normalisng constant for the Matérn spectral density on the sphere p(i) x ( +  From equation 31 the weight- space view of Gaussian processes implies that f is equivalent general linear model with deterministic feature maps of the form
1.1 88 (s) (8)=1 V
a²
2v
+1) Ø
-
(s),s),, j= 1_1...,,00 = 2j +1.
(33)
That is,
00 2j+1
f(s) = Σ
0 fj
Jink
(34)
det
for fi ~ No, 1) i.idd. for all 1,.... DO and & = 1 2j + 1. This set of features can cause issues when using them in practice, since spherical harmonics have increasingly high fluctuations as j + oo, eventually not being realisable numerically due to aliasing issues. For example, using single precision floating point arithmetics, observe that such issues oceur around the level j = 20. We observe that random features, which we will consider next, are better suited if we wish keep floating arithmetics single precision, in order to save memory.
A.2.1 RANDOM FEATURES ON THE SPHERE
We now derive an alternative feature map representation of Matérn GPs over S² that is analogous to the random features equation 25 in the Euclidean setting. First, the addition theorem for spherical harmonics states that
2j+1
2j+1 + 0/ (dga (s,5)),
(s)(s)))))
(35)
4m
15

---

Published as a conference paper at ICLR 2025
where dga (..)) denotes the geodesic distance the sphere and 9% (:) are the Gegenbauer polyno- mials. Furthermore, (Azangulov et al. 2024 Proposition 7) gives us that
(de (8,6)) 2j + 4x
[ (dss (dex (s) (dez (6',u) (sa)) (du);
(36)
where uga is the invariant measure on the sphere. Putting all of this together gives
kgo (s, 80) = o² L (247
(
( (s,u)) (dga (s', u) usz (du) Vs, E e
4T
(37)
where pM (dw) is the discrete measure
1
2v
1-
n (dw)
+ +
(d()).
(38)
C2
(2
In practice, consider truncated series for equation discrete measure
giving us equation 37 with the finite
(dw)
= Cv.J j=1
Σ
J
2v
C
+2,
1--
dy(dw), C2. =
Σ
J
2v
12
-P-1
(39)
5=0
which is equivalent the multinomial distribution
(C,))
2v 82 + t
2v
)
-1-1
w ~
Multinomial
C.,.
tis
(40)
f2
Thus, by Monte Carlo approximation, we have
2wmm +
kg- 8') Re
G1l2 (dge (s, um) ) (d (dge (d,um)) Va,s' E S²,
(41)
4T
with b E US?) the uniform distribution on the sphere, and is sampled according equation 40 This implies the random spherical features
/ lcae G2/2 (8, bm)), 8 EES? m = 1,...,,,, M, Co = (o(2w + 1)
4x
(42)
We may also conisder the limiting kernel // + oo, giving us the so-called heat kernel (Borovitskiy al. 2020). which is analogous to the gguared-exponentn kernel in the Euclidean case. This the same as before, except now the measure to sample w reads
(dw) = C j=1
Σ e
C =
(43)
j=0
From a more general perspective, on the two-sphere S?, we have an analogous notion of a kernel 8 x + IR to be stationary if it satisfies
kxx, v) k(Oz, Oy).
(44)
for all I, 82 and O E 5O(3). Then general result in (Azangulov et al. 2024) states that for any stationary kernel on 8 there exist a random feature representation, extending the classic result in Euclidean space by Rahimi & Recht (2007). In fact, this construction can be further generalised to any kernels over homogeneous spaces, which g is merely a special case (Azangulov et al., 2024). However, for the purpose of this paper, do not need to extend beyond the spherical setting.
16

---

Published as a conference paper at ICLR 2025
B ADDITIONAL DETAILS FOR DEEP RANDOM FEATURES
B.1 ADDITIVE RANDOM FEATURES
In our deep neural network architecture, we add skip connections from the input layer to each bottleneck layer in order to prevent the emergence of pathological behaviours as we increase depth (Duvenaud al.
2014). When considering spherical inputs, this requires building layers R We achieve this by taking Φf to be the H random features of GP R with the additive kernel k((2,, 5), (a' s')) = kg (ær, x') - Kgg (s. 8'), where kgn :R x RB RBxB and Agz S x g - B are matrix-valuee stationary kernels on respective spaces. This corresponds to an additive GP (Duvenaud et al. (01TT) f fR² + Sgx, where fee R Rd and fs² S + I are independent GPs with kernels kgs and kga, respectively. In particular, this implics layers ø' of the form
R
&² $ Φ
(45)
where on &bz are the random features corresponding to the GPs fR", fsz, respectively.
Putting this together, our final deep neural network architecture fe : S² BO on the sphere with skip connections reads, for all E S
(s) @ ())) h' s) == O² (0%0 (x) + (s)), 2 x x R (46) fe(s) (s), where (s) (s), (s) h (hh t= 2 L, (47)
where of E BBxH for = I,..,, and O E R BxO are trainable weights.
B.2 FUNCTIONAL REGULARISATION
Here, we provide details on the functional regularisation that we use regularise training of the model hyperparameters, We assume the case O = 1 (one output dimension) for simplicity, but this can be extended easily to multiple output dimensions. In the planar model, the regularisation term reads
V / v (xx,t) fax l)dæ.
(48)
The gradient is computed using Py Porch''s automatic differentiation and the integral is computed using the trapezoidal rule (note: our input dimesion is typically small, e.g. two or three).
In the spherical setting, must adapt the computations account for the curvature of the sphere. First, we consider parameterisation of the sphere in the following spherical coordinates:
(Longitude) (Latitude)
to E (0,2
(49)
E [-T/2,
(50)
We consider 7.2-inner product with respect to the Haar measure on the sphere μ = cos d0 dyp and consider the Laplace- -Bettacm operator the sphere, which reads
1
Of) 00
A
cos
(51)
cos 0 00
cos2 0 aPp
By integration-by-parte we get the following expression for functional regularisation the sphere
r/2
a
as 00 )
(.(AAs))) =
((0,4)
cos 0
cos 0 d0 dep
cos 0 de
cos? A d.p
1/2
(52)
2x
x//
05 2
1
[5 2
cos 0 do dyp
(53)
J_n/2
a0
cos2
0
aΦp)
Again, the gradients are computed via automatic differentiation and the integral is computed using the trapezoidal rule.
17

---

Published as a conference paper at ICLR 2025
Remark 2 We consider functional regularisation regularise our hyperparameterss X instead of considering hyperpriors. since in deep GPs. the kernel lyperparaneterss in each layer are less in- terpretable compared to shallow GP hyperaparameters. Therefore, it becomes questionable prior on the hyperparameters, when we do not really have good knowledge of them. The functional regulrisation that we here do not require prior knowledge it assumes prior knowledge on the corresponding model fe(; X), which is more realistic; for in- stance, we might have some prior knowledge, such as smoothness, of the ground truth field f that we are trying to model. This knowledge can then be used to regularise the search for X, without directly imposing prior on X.
the hyperparameters. Instead,
C EXPERIMENT DETAILS
C.I EVALUATION METRICS
Denote by ((X,, yo)) the training set and D = {(X,,, yn)) Xn (In, (n) denotes spatiotemporal coordinates. We consider the following metrics to evaluate
1N2 be test set, where
our model f. trained on the set D.
Root Mean Squared Error (RMSE). We first consider the root mean squared error, given by
(RMSE(f:D) = B x lll
E Ef
(54)
yy E [X,
(55)
where we denote by [[ the conditional expectation with respect the event of observing the training set D. For ensemble based models, we approximate the conditional expectation by the empirical mean of the trained models.
Root Mean Absolute Error (RMAE). Similar to the RMSE score, we consider the root mean absolute error, computed
CRMAE(f:D) E(x.1) v) [ll Eff(() II]
(56)
N Σ -(f n=1
(57)
N
Compared to the RMSE, the RMAE is more robust to outlier values, making it a more reliable metric when the data generating model has heavier tails than Gaussian.
The RMSE and RMAE only evaluates the quality a single statistic the distribution f (here, used the predictive mean). Below, introduce metries that also evaluates the quality of the predictive uncertainties.
Negative Log- Like-iield (NLL). When have access the ground truth field fi, we the negative log-likelihood score to evaluate how well our model's predictive distribution describes the ground truth. This is computed for each input coordinate X
( (XX), D) log (f
(58)
where PDD denotes the probability density function of the random variable f(X) In partic- ular, assuming that this is Gaussian, P):x)p = (ylux Σxx), we we get
58
(f' -xx) Σxx (fl (X) /x) log Σxxl const.
(59)
We then use the following empirical risk to assess the overall quality of uncertainty of our model f
(NLL(f: D) = Ex [f (XX)
(60)
1
NLL(F(X.)): f (X.).D)
(61)
N
18

---

Published as a conference paper at ICLR 2025
For ensemble based models, use its Gaussian approximation constructed from the empirical means and variances to evaluate the NLL
Negative Log-Predictive Density (NLPD). On the other hand, when we do not have access to the ground truth field f), we can use the negative log-predictive density to evaluate the uncertainty estimates of our model. This is given by
NLPD(f:D.X.u))
log B pllyf.X)
(62)
In the case of Gaussian likelihoods, i.e., plylf. x) (x)) and Gaussian posteriors f|D (p, E). we have the following closed-form expression for the NLPD
02
log. N2y px. Σxx 1
(63)
y p
[xx + o²T)- (リー px) + log Σxx + o² + const.
(64)
2
2
Again, we use the corresponding risk to evaluate the predictive uncertaintics of our model
[N E() NNLPD(f: D, x))
N x,,,)
(65)
DL yn).
(66)
N ==1
As with the NLL, for ensemble based models, use the Gaussian approximation constructed from its empirical means and variances to evaluate equation 64
Remark In non-conjugate models, we may consider the following upper bound of the NLPD
66
Bip |- poy xx)]
(67)
J lo ~ p(pl))
(68)
where we used Jensen's 's inequality in the second line and Monte Carlo approximation in the last line. The resulting emprical risk reads
1
CNLPD
log Degluvile ogg(vilff, p(f(D.
(69)
Continuous Ranked Probability Score (CRPS). The continuous ranked probability score is other metric used to evaluate the quality of predicted uncertainties. Compared to the NLL or NLPD, CRPS is more robust to outliers and measures holistic calibration rather than pointwise calibration (Gneiting & Rattery 2007).
Given data y E R. the CRPS of random variable X R for modelling the data y is computed as
CRPS(X;
(Px < <) H(x < y)) dx.
(70)
where Px denotes the law X. Alternativly, one can write this as
CRPS(X: = E, || - -ll En. [| x']],
(71)
where E,,, Es,.r' denotes expectation with respect to the laws Ix,, Px @ Fx. This formulation can be used to evaluate the CRPS empirically using i.i.d. samples ar, x' ~ IPx. Furthermore, when X is Gaussian with mean and standard deviation a,, we have a closed form expression of the CRPS of the form
CRPS(X:v) =o
-0(y/ 24
y - )
-1)
+ 20
T -1/2 (
Y-
(72)
where @ and & are the probability density function and cumulative density function respectively of the standard Gaussian random variable.
19

---

Published as a conference paper at ICLR 2025
Assuming that the output dimension of our model is (i.e. O 1), get the following metric for assessing model calibration based on the CRPS
LcRrS (J;D) Exx,v a) [CRPS([(X))P:o))
(73)
N
(x, (D;
(74)
C.2 EXPERIMENT 1: EVALUATION ON A SYNTHETIC DATA
To generate our synthetic ground truth field ft for this experiment, we create a MSS using years of altimeter readings of sea surface height (SSH) from CS2 (2010-2022) covering the polar region densely. To de-slope the MSS we subtract the FGM2008 geoid (Skourup et al. 2017). We bin measurements 5x5km grid in the arctic and take the average in each bin to get estimate of a typical polar mean SSH field. The resulting field is suitable for our assessment since exhibits spatial non-stationarity and contains several high frequency features, which would be challenge to recover accurately. To generate the artificial measurements, we first take nine days of space-time time coordinates (between March 1st and 10th 2020) of time-evolving tracks from Sentinel-3A, 3B and CryoSat-2 satellites, then extract the value of our artificially generated mean SSH field at those spatial coordinates, and finally add i.i.d. Gaussian noise to mimic epistemic uncertainty. That is, for space-time coordinates (In, (n) of the satellite track, we take
yn = ( (n) + Ens ~ (O,
(75)
where we set Fg 0.01. Our final dataset comprise 1,158,505 datapoints; we select 80% of these randomly for training and the remaining 20% for validation.
C.2.1 TRAINING DETAILS
As usual, let = (X,,, yn yr) > denote our training data, where X, = -xx,, f..) for 7 = 1,..., N spatiotemporal coordinates of the satellite tracks. We train model using the MSE loss, giving us the regularised empirical risk
N Liaii (0:))
= N Σ (z),))) - yn 22 81|11
(76)
which we minimise find the optimal neural network weights O In this experiment, we assume knowledge of the data generating process in equation 75 giving us the negative log posterior
N
1 - 1012 const.
- log D)
fe( b.) yn
(77)
20
n=1
Comparing equation 76 with equation
77
we have the relation
3 = ==o o
(78)
Assuming knowledge of the observation noise oy we the weight decay parameter / accordingly.
We summarise the training details as follows:
Optimiser: Adam
Learning Rate: 0.001
Loss Function: Mean Squared Error (MSH) with 1.2 regularisation (equation 76
Number of Epochs:
Batch Size: 1024
C.2.2 DRF MODEL DETAILS
For the DRF model used in this experiment, we use the following model configurations:
Bottleneck size: B 128 for both spatial and temporal layers.
20

---

Published as a conference paper at ICLR 2025
Hidden unit size: H 1000 for both spatial and temporal layers.
Number of spatial layers: L2 1, 2, 4, 10, 20. For results in Table 1 fix Lu 4.
Number of temporal layers: Le 1.
We use random features corresponding to the Matérn-3/2 kernel for both spatial and temporal layers (see Example 2 Appendix A.I)). I). The lengthscale and amplitude hyperparameters for both the spatial and temporal layers tuned with Bayesian optimisation using the Python library BoTorch Balandat et al. 2020) with the following configurations:
Spatial Lengthscale (Aspaiiat): Optimised over the range [0.1, 10]
Temporal Lengthscale (Cremporat): Optimised over the range [0.1, 10]
Amplitude (o): Optimized over the range [0.1, 1]
Observation Noise (ay)): Not optimised; set to 1 x 10-2 in the loss function
Total Iterations: 10
Initial Samples: 10 (generated via Latin Hypercube Sampling)
Acquisition Function: Expected Improvement (EI)
Surrogate Model: Gaussian Process Regression (SingleTaskGF)
Number of Restarts for Optimisation: 100
Number of Raw Samples for Initialisation: 100
We used the same Bayesian optimisation configurations for DRF with different UQ methods (Deep Ensembles and MC Dropout).
C.2.3 UQ DETAILS
Deep Ensembles We consider 10 models to form an ensemble in this experiment. Each model in the ensemble uses the same architecture 4 Matérn random feature layers with skip connections, each trained for epoch per Bayesian optimisation iteration. The training and evaluation processes for the ensemble are parallelised using the joblib library, allowing efficient computation across multiple cores.
Monte-Carlo Dropout We consider the same architecture (4 Matérn RFF layers) for this UQ method, however we trained for 10 epochs to make sure the model is well-converged before dropout can be used effectively during inference. The implementation is based on the Python pack- age 1 liigghiinnu---ebae (Lehmann) 2024).
Variational Inference We consider the same architecture (4 Matérn RFF layers) for this UQ method, however we use ELBO maximisation (equation 14) for hyperparameter tuning. The im- plementation is based on Lightning-uq-box.
C.2.4 BASELINE DETAILS
GPSat. GPSat (Gregory al. 2024b) Python package for making predictions based on mixture of local GP models, developed specifically for satellite interpolation. Individual local expert models are trained on data subregions of the domain of interest and combined together in post-processing step. We use expert locations spaced evenly Scalable Earth (EASE) grid covering the arctic region. This results in 1225 expert locations. Each local expert model uses the Sparse Gaussian Process Regression (SGPR) model in Titsias (2009) for inference, each with 500 inducing points and trained with the L-BFGS-B optimiser.
200km resolution Equal-Area
Sparse Variational Gaussian Process (SVGP). We consider the sparse variational Gaussian pro- cess model in Hensman et al. 2013) with 3000 inducing points. The model is trained using the Adam optimiser with a learning rate 0.01, a minibatch size of 1024, and 5 epochs for training. Implemented using the Python package Torto (Gardner et al. 2018)
21

---

Published as a conference paper at ICLR 2025
Deep Gaussian Processes (DGP). We consider deep Gaussian process with two layers in this experiment. We use the doubly stochastic variational inference (Salimbeni Deisenroth 2017) for training with 512 inducing points per layer, Monte Carlo sample size of 10 and minibatch size of 1024. We found that increasing the number of layers did not lead to significant improvements. The model is trained using the Adam optimiser with learning of 0.01 and trained for 10 epochs. Implemented using the Python package GPyTorch.
ReLU MLP. We consider a vanilla RELU MLP with hidden layers, 1024 hidden units and trained with learning rate of 0.01 for 5 epochs. We used deep ensemble of 10 models to obtain the uncertainty estimates.
Fourier Features Network (FFN). Similar to SIREN, consider FFN with 4 layers (1 RBF RFF layer and fully connected linear layers). 512 hidden units for each layer and trained with learning rate of 1 x 10-5 for 10 epochs. We used deep ensemble of 10 models to obtain the uncertainty estimates. For the positional encoding, we used random Fourier features (equation 25 and employed Bayesian Optimisation tune the kernel hyperparameters & and σ.
We consider a SIREN network with 4 hidden layers, 512 hidden units and trained with a learning rate of 1 x 10-5 for 10 epochs. We used deep ensemble of 10 models to obtain the uncertainty estimates.
SIREN.
Conditional Neural Processes (CNP), The encoder consists fully connected layers with di- mensions [Fdim + Mdim: hamm...... -diml to extract the context information, and the decoder maps this context to the target values through a series of layers [zamm + Tdim: ham- .., 2 x 3dmm]
,
Encoder dimensions: [ddim + Udim = 4, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000]
Decoder dimensions:
[diim +Tdim 1003, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 2
Ydim
Remark 4 The architectures for the neural network baselines were selected by optimising (via grid- search) the mean NLPD score (equation 62) on the held-out validation
C.3 EXPERIMENT 2: FREEBOARD ESTIMATION FROM SATELLITE ALTIMETRY DATA
In this experiment, we use along-track data from the Sentinel-3A+B (S3A, S3B) and CryoSat-2 (CS2) satellites measuring sea-ice freeboard (in fact, the satellite measurement is the radar freeboard, which is distinct from sea-ice freeboard. This converted to sea-ice freeboard a postprocessing step). Freeboard measures the height of ice surface relative sea level, which can be used to compute sea ice thickness (Gregory et al. 12021); the latter is important for monitoring the effect of climate change the polar climate system.
Freeboard measurements
0.4
Sea-ice freeboard (cm)
0.3
0.2
0.1
0.0
0.1
80
82
84
86
88
Longitude (*E)
Figure 7: Freeboard measurements along a sample satellite track segment (highlighted in orange in the left figure). Real satellite measurements typically very noisy.
22

---

Published as a conference paper at ICLR 2025
We use the freeboard observations along S3A, 3B and CS2 satellite tracks downloaded from Gregory et al. (2024a). Our final dataset contains 1,174,848 data points. Since we do not have access to ground truth field, hold out separate test set that use 70% of the data (randomly sampled) for training, 15% for validation and the remaining 15% for testing. The satellite measurements challenge to interpolate robustly, picking out the signal from the noise (see for example Figure for freeboard measurements along an example track).
make our final evaluations. We
extremely noisy and contain many outliers, making it
C.3.1 DRF MODEL DETAILS
We use the same model architecture and Bayesian Optimisation settings as in our first experiment. For training, we also use the same procedure as our first experiment, except we now treat the observation noise by tunable hyperparameter, learned jointy with the kernel hyperparameters via Bayesian optimisation.
C.3.2 BASELINE DETAILS
GPSat.
We the same configuration as in first experiment.
Sparse Variational Gaussian Process (SVGP). We use the same configuration in our first experiment.
Deep Gaussian Processes (DGP). We use the same configuration as in first experiment.
ReLU MLP. We consider a vanilla RELU MLP with hidden layers, 1024 hidden units and trained with learning rate of 0.01 for 5 epochs. We used deep ensemble of 10 models to obtain the uncertainty estimates.
Fourier Features Network (FFN). We consider FFN with 4 layers (1 RBF RFF layer and fully connected linear layers), 256 hidden units for each layer and trained with learning rate of 0.01 for epochs. We used a deep ensemble of 10 models to obtain the uncertainty estimates. For the positional encoding, we used random Fourier features (equation 25 and employed Bayesian Optimisation to tune the kernel hyperparameters i and σ.
We consider SIREN network with 6 hidden layers, 512 hidden units and trained with learning rate of 0.01 for 5 epochs. We used a deep ensemble of models to obtain the uncertainty estimates.
SIREN.
Conditional Neural Processes (CNP). We use the same configuration as in our first experiment.
C.4 EXPERIMENT 3: LARGE SCALE INTERPOLATION OF GLOBAL SEA LEVEL ANOMALY
In our final experiment, we used global sea level anomaly measurements from the Sentinel 3A satellite for the period March 1st-4th, 2020, resulting 8,094,569 data points. We use 80% of these for training and 20% for validation.
C.4.1
DRF MODEL DETAILS
For the DRF model used in this experiment, we use the following model configurations:
Number of spatial layers: One spherical random feature layer + two Euclidean random
feature layers.
Number of temporal layers: L, = 1.
Bottleneck size: B = 128 for both spatial and temporal layers.
Hidden unit size: H 1000 for both spatial and temporal layers.
23

---

Published as a conference paper at ICLR 2025
We use random features corresponding to the Matérn-3/2 kernel for both spatial (spherical and Eu- clidean) and temporal layers. The lengthscale and amplitude hyperparameters for both the spatial and temporal layers tuned with Bayesian optimisation with the following configurations:
Spatial Lengthscale for Spherical Layer (Sspherica)): 10-$,0.1]
Optimised over the range x
Spatial Lengthscale for Euclidean Layer (Cenclidean): Optimised over the range [1 x
10-9,10]
Amplitude for Spherical Layer (ospheriall:: Optimised over the range [1 10
Temporal Lengthscale (Creporaa): Optimised over the range 10 10]
Amplitude for Euclidean Layer (uuclidam): Optimised over the range [1 x 10 10-5,
Total Iterations: 10
Initial Samples: 15 (generated via Latin Hypercube Sampling)
Acquisiiinn Function: Expected Improvement (EI)
Surrogate Model: Gaussian Process Regression (SingleTaskGP)
Number of Restarts for Optimisation: 100
Number of Raw Samples for Initialisation: 100
We train the model with a Huber loss
[l(( yl² ((f(() - y 6/2)
Luber (f: X, (,),
(79)
with 5 = 0.1. This loss is more robust than the mean squared error to outlier values.
C.4.2 BASELINE DETAILS
Spherical SVGP. We used the spherical Matérn-3/2 GP (Borovitskiy et al. 2020) with 3000 in- ducing points, trained using the Adam optimiser with learning rate of 0.01 and minibatch size of 1024. Implemented using the Python package GPyTorch with the spherical kernel implemented with the geometric senetrrrelekkkeeeeee package.
Euclidean DRE. We used a Euclidean DRF model with two layers, each with bottleneck dimen- sion B = 128 and hidden dimension H 1000. We also utilised Bayesian Optimisation to tune the kernel hyperparameters Expatial (temporrll and a.
C.4.3 QUANTITATIVE RESULTS
Model
NLPD
CRPS
Spherical DRF
0.0063 0.0002
0.0918 0.0021
Euclidean DRF SVGP
0.0387 = 0.0001
0.1158 + 0.0008
0.0085 = 0.0001
0.1210 + 0.0010
Table 3: NLPD and CRPS comparison for Spherical DRF, Fuclidean DRF and SVGP Models
In Table we present the guantitative results for experiment 3, comparing the performance of spherical DRF, Euclidean DRF and spherical SVGP models using the NLPD and CRPS metries. For the NLPD computation, we use its approximation (equation 69 with the likelihood defined by the Huber loss. That is, p(y f X) exp (- X, y, 5)).
Consistent with what we see in the qualitative performance (Figure 6). the spherical DRF signifi- cantly outperforms the other baselines, achieving the lowest NL.D and CRPS.
24

---

Published as a conference paper at ICLR 2025
C.5 GALLERY OF PREDICTIONS
C.5.1 GALLERY OF PREDICTIONS: EXPERIMENT
Mean: DRF (Deep ensemble)
Variance: DRF (Deep ensemble)
Ground Truth
Mean: DRF (Dropout)
Variance: DRF (Dropout)
Ground Truth
Mean: DRF (VI)
Variance: DRF (VI)
Ground Truth
-0.50 -0.25 0.00 0.25 0.50
0.000
0.005
0.010
-0.50
-0.25
0.00
0.25
0.50
Figure 8: Comparison of predictive means and variances of the DRF model using various UQ meth- ods the synthetic experiment. Left column: predictive means, Middle column: predictive vari- ances, Right column: ground truth MSS field.
25

---

Published as a conference paper at ICLR 2025
Mean: SVGP
Veriance: SVGP
Ground Truth
Mean: Deep
Variance: Deep GP
Ground Truth
Moon: GPSat
Variance: GPSot
Ground Truth
Varience: ReLU
Ground
Mean: FFN
Variance: FFN
Ground Truth
Mean: SIREN
Variance: SIREN
Ground Truth
Mean: CNP
Variance: CNP
Ground Truth
0.50 -0.25
0.00
0.25
0.50
0 0.000 O 0.00 0.006 O. 0.0000 0.000
-0.50
-0.25
0.00
0.25
0.50
Figure 9: Comparison of the baseline models on the synthetic experiment. Left column: predictive
means, Middle column: predictive variances, Right column: ground truth MSS field.
26

---

Published as a conference paper at ICLR 2025
C.5.2 GALLERY PREDICTIONS; EXPERIMENT
Mean
SVGP Mean
DeepGP Mean
GPSat Mean
0.3
0.2
0.1
0.0
DRF Variance
SVGP Variance
DeepGP Variance
GPSat Variance
0.010
0.008
0.006
0.004
0.002
0.000
RELU Mean
FFN Mean
SIREN Mean
CNP Mean
0.3
0.2
0.1
0.0
-0.1
Variance
Variance
SIREN Variance
CNP Variance
0.010
0.008
0.006
0.004
Figure 10: Comparison of DRF and other baseline models on freeboard interpolation. We display both predictive means and variances for each model.
27

---

Published as a conference paper at ICLR 2025
C.5.3 GALLERY PREDICTIONS: EXPERIMENT
Spherical DRF (mean)
Spherical DRF (varlance)
30
10
-10
-10
30
-30
50
550
-70
-70
90
90
-
180-140
-000
-60
-20
20
60
100
140
180
-180
140
-100
-60
-20
20
100
140
180
Euclidean DRF (mean)
Euclidean DRF (variance)
90
90
70
70
50
50
30
30
10
10
10
30
50
-70
-180-140- -100 -60
100
140
180
-180-140
-60
-20
180
Spherical SVGP (mean)
Spherical SVGP (variance)
90
100
-10
30
33
50
50
70
-70
90
90
180
140
100
20
20
60
100
140
180
180
-40
1100
-20
20
60
100
140
180
0.15
0.10
0.05
0.00
0.05
0.10
0.15
0.000
0.025 0.050 0.075
0.100
0.125 0.150 0.175 0.200
Figure 11: Result maps (mean and variance) for global sea level anomaly interpolation. From top to bottom: Spherical DRF, Euclidean DRF and Spherical SVGP.
28

---

Published as a conference paper at ICLR 2025
C.6 FURTHER ABLATION STUDY
C.6.1 DATA SCALABILITY
Here, we conducted an experiment where we trained on 10% and 1% subsampled data of our global sea level anomaly interpolation experiment. The dataset partition follows the following ratios: 70% training, 15% validation and 15% testing splits to compute the metrics. The results can be found in Table where we compare our DRF model against SVGP (3, 000 inducing points, trained for 2 epochs).
This show that, as expected, the DRF model achieves better performance with increasing data size, as evidenced by the lower NLPD and CRPS values. We also find that for larger data (10% and 100%), DRF performs better than SVGP due to its ability to capture finer and finer scale features as we increase data resolution (note that the performance of SVGP does not change much increasing the data from 10% to 100%). On the other hand, for smaller data proportions (1%), the simpler inductive bias in the SVGP helps generalise better, resulting in better performance compared to DRF. See figure 12 for comparisons.
(a) DRF trained with 100% dataset
(b) SVGP trained with 100% dataset
(c) DRF trained with 10% datasct
(d) SVGP trained with 10% dataset
(e) DRF trained with 1% dataset
(f) SVGP trained with 1% dataset
Figure 12: Results maps shows data scalability and corresponding performances for DRF (left col- umn) and SVGP (right column)
C.6.2 SKIP CONNECTIONS
Below, we show how the addition the skip connections affect our global sea level anomaly in- terpolation results. We find that without the skip connections, get unstable results, sometimes
29

---

Published as a conference paper at ICLR 2025
Dataset
NLPD
CRPS
Time (mins)
Full Dataset-DRF
0.0063
0.0927
291.15
Subsampled-DRF (1 in 10)
0.0079
0.1114
45.61
Subsampled-DRF (1 in 100)
0.1314
0.6490
17.82
Full Dataset-SVCG
0.0085
0.1205
345.04
Subsampled-SVGG (1 in 10)
0.0086
0.1200
47.46
Subsampled-SVGG (I in 100)
0.0161
0.1435
18.10
Table 4: Comparison NLPD, CRPS, and Time for Various Dataset Sizes for DRF and SVGP
leading to poor reconstructions of the field as illustrated in the left figure 13a responding metrics in Table 5 where we see that the model with the skip connection produce better scores result.
We display the cor-
For our other experiments, we did not see major differences between the model with or without skip connections.
(a) Example of poor result when using DRF model without skip connections
(b) Typical prediction of DRF model with skip con- nections
Figure 13: Performances of DRF models with and without skip connections
Model
NLPD
CRPS
Model with skip connections Model without skip connections
0.0063 = 0.0002
0.0918 + 0.0021
0.0078 # 0.0007
0.1092 # 0.0101
Table 5: Comparison of NLPD, CRPS, for model with or without skip connections
C.6.3 SCALABILITY WITH MODEL WIDTH
Here, we present the results of ablation experiment with respect to model width defined as the number of random feature per layer. In the first part, fixed the bottleneck size B per layer and varied the number of hidden units between the bottlenecks, while keeping the model depth fixed 3. The number of hidden units H directly corresponds the number of random features per layer. Our findings indicate that, in general, increasing the number of random features (i.c., enlarging H) leads improved performance. shown in T
DRF Model
NLPD
CRPS
1000, 128
0.0063 +: 0.0002
0.0918 + 0.0021
H 500, B = 128
0.0064 0.0001
0.0944 # 0.0009
= 250, B = 128
0.0065 = 0.0001
0.0951 + 0.0006
Table 6: NLPD and CRPS Comparison for DRF models with different hidden size
30

---

Published as a conference paper at ICLR 2025
In the second part. have conduced ablation experiment with respect the bottleneck size B while hidden size held fixed (H 1000).
DRF Model
NLPD
CRPS
H 1000, B 32
0.0066 = 0.0001
0.0969 + 0.0006
1000, 64
0.0067 0.0003
0.0969 0.0030
H 1000, B 128
0.0063 + 0.0002
0.0918 + 0.0021
1000,
1000
0.0065 0.0003
0.0932 0.0035
Table 7: NLPD and CRPS Comparison for DRF models with different bottleneck size
We also observe that increasing the bottleneck size generally improves performance, but this only holds true up to certain point. Beyond that, such as when we begin to deteriorate, and the hyperparameter tuning process becomes less stable. Figure [4 shows an example poor results we got when we set H 1000, B 1000.
1000, the results actually
Figure 14: Example of poor result of DRF model with H/ 1000, B 1000
31


Reviews:
[{"review_id": "u7boXzOpYa", "reviewer": "Reviewer_VtAL", "summary": "The paper introduces a model leveraging deep random features (DRF) to achieve scalable interpolation for large-scale, high-dimensional spatiotemporal data. By integrating random Fourier feature layers with DNNs, the model retains the inductive biases of GPs, while with potential for scaling efficiently with the dataset size. The novelty lies primarily in adapting random feature layers to produce an architecture capable of high-resolution spatiotemporal interpolation with uncertainty quantification, moving beyond the limitations of traditional GPs and recent DGPs in this domain.", "questions": "- To me, the proposed method closely resembles the Fourier features network and SIREN in various aspects. Specifically, it feels like a generalization of Tancik et al.’s repeated Gaussian mapping followed by linear transformations, applied across layers rather than in just the input. If random features were only used in the first layer, would the network’s performance still mirror that of the Fourier features network? Exploring this could clarify if deeper layers with random features add unique value.\n- The model also feels akin to a generalized SIREN, with deterministic parameter layers replaced by additional linear transformations and stochastic components. In this context, I'd appreciate a breakdown of the number of trainable parameters in this model versus SIREN and other baselines. This would help me interpret differences in the quantitative results presented in Tables 1 and 2.\n- Since the model’s roots are claimed to be in deep Gaussian processes, I think it would be helpful to see an ablation study comparing deep GPs approximated with random features against this proposed model. This doesn’t need to be overly complex—a simple synthetic dataset might be sufficient to highlight the supposed benefits of deep random features combined with learnable linear layers.\n- I find the skip connections interesting but somewhat confusing. From what I see, after the projection by $\\phi$, the features are concatenated back to a space seemingly compatible with the original coordinates. However, I’m unsure if this compatibility holds without explicit regularization. An ablation study on removing or adjusting these skip connections would help clarify their impact, making it easier to interpret how well they function as intended.\n- The paper highlights scalability as a major benefit, but it lacks empirical validation in terms of how the model performs with increasing data points or the number of random features. Scalability results as a function of data volume and feature count would substantiate claims about the model’s efficiency and approximation performance.\n- Lastly, given the reliance on random sampling for the random Fourier features, I want to see the variability in the model’s performance. Not on a UQ sense but evaluating the model across multiple runs with varied random seeds.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- The paper integrates multiple uncertainty quantification techniques (variational inference, dropout, deep ensembles), offering flexibility in obtaining uncertainty estimates for the model outputs.\n- The architecture performs well across varied datasets, including synthetic, local, and global satellite measurements, and the spherical adaptation for global data fields demonstrates its versatility.", "weaknesses": "I don’t see any major concerns or weaknesses in this paper, but I believe it would benefit from some ablation studies, which are currently missing, to better illustrate the contributions of the proposed model. Please refer to the questions for details.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "NhvmlNaZbh", "reviewer": "Reviewer_pDLc", "summary": "The paper presents a deep Gaussian process (GP) with deep random features for modeling spatiotemporal data where the inputs are spatial/temporal coordinates. Deep random features are utilized to address the high-frequency limitations of neural networks, and mini-batched gradient descent is used for large-scale training. Experiments are conducted on the synthetic data, satellite data, and sea level anomaly data.", "questions": "See above.", "limitations": "", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "- The paper addresses a clear and timely topic, focusing on coordinate-based neural representation for sparse observations. \n- Figures 1 and 2 effectively illustrate both the task and the model structure. \n- The theoretical background on random features is covered in sufficient depth.", "weaknesses": "- The model does not appear to outperform baselines in terms of accuracy or computation time. \n- Although the primary contribution is the application of random features, the paper does not cover advanced Fourier, Wavelet features, such as [1]. It would be helpful for the authors to consider incorporating these features or at least to discuss how their approach might be extended to include them. \n- The assumption of stationary kernel seems limit model performance, and further clarification on this assumption would be beneficial.\n- It would be beneficial for the authors to include tasks for both time interpolation and extrapolation, like [2].\n- The paper also overlooks many recent works that apply implicit neural representations for spatiotemporal data, which should be included for a more comprehensive comparison. For the existing works, temporal dimensions are often treated separately from spatial dimensions, typically using autoregressive time-stepping, which has shown efficiency in modeling physical spatiotemporal data. For example, [2], and [3] treat the spatial dimension continuously and employ neural ODEs to model time as a continuous variable. The author should justify why the proposed method has advantages over existing continuous spatial-temporal models. \n\n[1] Multiplicative Filter Networks, ICLR, 2021.\n\n[2] Continuous PDE Dynamics Forecasting with Implicit Neural Representations, ICLR, 2023.\n\n[3] Operator Learning with Neural Fields: Tackling PDEs on General Geometrics, Neurips, 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "aW8zkoCS44", "reviewer": "Reviewer_jnCy", "summary": "The paper proposes to blend the strengths of neural networks (NNs) and Gaussian processes (GPs) by introducing the deep random features (DRF) framework, where the weights of the Fourier expansion $\\boldsymbol{\\Theta}$ are viewed as trainable parameters (similar to weights in NN) and the random features themselves are viewed as neurons. To perform interpolation on sphere, the authors also construct random feature maps based on a kernel constructed by a Mercer sum. The hyperparameters of the model are tuned either by comparing the evidence lower bound (ELBO) or by using a held-out validation set, depending on the choice of uncertainty quantification method. The authors test DRF on both synthetic and remote sensing datasets to show that DRF is able to learn both coarse-scale and fine-scale data patterns with well-calibrated uncertainty.", "questions": "•\tSince $\\boldsymbol{\\omega} \\sim p(\\boldsymbol{\\omega})$ and $p(\\boldsymbol{\\omega})$ is the normalized Fourier transform of the kernel function, can you elaborate the procedure to efficiently sample from $p(\\boldsymbol{\\omega})$?\n\n•\tWhy does DRF with dropout have much higher computational time than an ensemble of DRFs?\n\n•\tDoes it really make sense to compare with neural processes in Tables 1 and 2 as NPs are mainly used for meta-learning?\n\n•\tAre ground truth values available for the sea level anomaly measurements? If so, it would be beneficial to present the interpolation and uncertainty calibration errors to validate the importance of learning fine-scale fluctuations.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "•\tThe authors present the methodology clearly, supplemented by informative visualizations, and provide a comprehensive discussion of the experimental results.\n\n•\tThe proposed framework integrates the flexibility of NNs with the inductive bias of GPs, offering a scalable solution for large datasets.\n\n•\tDetecting high-frequency patterns is important in environmental data analysis, as such patterns are often closely associated with extreme events.\n\n•\tWhile there exists some amount of literature on constructing random features from kernel machines, the application of random features for interpolating spatial and temporal data is relatively underexplored and can be considered novel.", "weaknesses": "•\tFor clarity, it would be better to differentiate $H$ from $B$ in Section 3.1, as $H$ plays a similar role as $M$ in Section 2.2 which represents the dimension of the random features, while the bottleneck dimension $B$ is analogous to the hidden layer dimension in NNs.\n\n•\tA more comprehensive discussion of related work would enhance the paper’s contribution. For instance, [1] investigated representation learning through the approximation of kernels using random Fourier features. Additionally, [2] explored the connection between deep ensembles trained with squared-error loss and GP posteriors.\n\n•\tAs noted in Section 2.2, random Fourier features are applicable only to stationary kernel functions. However, spatiotemporal data may exhibit non-stationary trends, which presents a potential limitation of the proposed method. A discussion on this point (e.g., extensions to handle non-stationarity) will be appreciated.\n\nReferences:\n\n[1]. Jiang, Z., Zheng, T., Liu, Y., & Carlson, D. (2022). Incorporating prior knowledge into neural networks through an implicit composite kernel. arXiv preprint arXiv:2205.07384.\n\n[2]. He, B., Lakshminarayanan, B., & Teh, Y. W. (2020). Bayesian deep ensembles via the neural tangent kernel. Advances in neural information processing systems, 33, 1010-1022.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)