PROBABILISTIC NUMERIC CONVOLUTIONAL NEURAL NETWORKS
Year: 2021

Abstract:
Continuous input signals like images and time series that are irregularly sampled or have missing values are challenging for existing deep learning methods. Coherently defined feature representations must depend on the values in unobserved regions of the input. Drawing from the work in probabilistic numerics, we propose Probabilistic Numeric Convolutional Neural Networks which represent features as Gaussian processes (GPs), providing a probabilistic description of discretization error. We then define a convolutional layer as the evolution of a PDE defined on this GP, followed by a nonlinearity. This approach also naturally admits steerable equivariant convolutions under e.g. the rotation group. In experiments we show that our approach yields a 3× reduction of error from the previous state of the art on the SuperPixel-MNIST dataset and competitive performance on the medical time series dataset PhysioNet2012.

heading: 1 INTRODUCTION
text: Standard convolutional neural networks are defined on a regular input grid. For continuous signals like time series and images, these elements correspond to regular samples of an underlying function f defined on a continuous domain. In this case, the standard convolutional layer of a neural network is a numerical approximation of a continuous convolution operator A. Coherently defined networks on continuous functions should only depend on the input function f , and not on spurious shortcut features (Geirhos et al., 2020) such as the sampling locations or sampling density, which enable overfitting and reduce robustness to changes in the sampling procedure. Each application of A in a standard neural network incurs some discretization error which is determined by the sampling resolution. In some sense, this error is unavoidable because the features f (`) at the layers ` depend on the values of the input function f at regions that have not been observed. For input signals which are sampled at a low resolution, or even sampled irregularly such as with the sporadic measurements of patient vitals data in ICUs or dispersed sensors for measuring ocean currents, this discretization error cannot be neglected. Simply filling in the missing data with zeros or imputing the values is not sufficient since many different imputations are possible, each of which can affect the outcomes of the network.
Probabilistic numerics is an emergent field that studies discretization errors in numerical algorithms using probability theory (Cockayne et al., 2019). Here we build upon these ideas to quantify the dependence of the network on the regions in the input which are unknown, and integrate this uncertainty into the computation of the network. To do so, we replace the discretely evaluated feature maps {f (`)(xi)}Ni=1 with Gaussian processes: distributions over the continuous function f (`) that track the most likely values as well as the uncertainty. On this Gaussian process feature representation, we need not resort to discretizing the convolution operator A as in a standard convnet, but instead we can apply the continuous convolution operator directly. If a given feature is a Gaussian process, then applying linear operators yields a new Gaussian process with transformed mean and covariance functions. The dependence of Af on regions of f which are not known translates into ∗Work done during internship at Qualcomm AI Research †Qualcomm AI Research is an initiative of Qualcomm Technologies, Inc.
the uncertainty represented in the transformed covariance function, the analogue of the discretization error in a CNN, which is now tracked explicitly. We call the resulting model Probabilistic Numeric Convolutional Neural Network (PNCNN).

heading: 2 RELATED WORK
text: Over the years there have been many successful convolutional approaches for ungridded data such as GCN (Kipf and Welling, 2016), PointNet (Qi et al., 2017), Transformer (Vaswani et al., 2017), Deep Sets (Zaheer et al., 2017), SplineCNN (Fey et al., 2018), PCNN (Atzmon et al., 2018), PointConv (Wu et al., 2019), KPConv (Thomas et al., 2019) and many others (de Haan et al., 2020; Finzi et al., 2020; Schütt et al., 2017; Wang et al., 2018). However, the target domains of sets, graphs, and point clouds are intrinsically discrete and for continuous data each of these methods fail to take full advantage of the assumption that the underlying signal is continuous. Furthermore, none of these approaches reason about the underlying signal probabilistically.
In a separate line of work there are several approaches tackling irregularly spaced time series with RNNs (Che et al., 2018), Neural ODEs (Rubanova et al., 2019), imputation to a regular grid (Li and Marlin, 2016; Futoma et al., 2017; Shukla and Marlin, 2019; Fortuin et al., 2020), set functions (Horn et al., 2019) and attention (Shukla and Marlin, 2020). Additionally there are several works exploring reconstruction of images from incomplete observations for downstream classification (Huijben et al., 2019; Li and Marlin, 2020).
Most similar to our method are the end-to-end Gaussian process adapter (Li and Marlin, 2016) and the multi-task Gaussian process RNN classifier (Futoma et al., 2017). In these two works, a Gaussian process is fit to an irregularly spaced time series and sampled imputations from this process are fed into a separate RNN classifier. Unlike our approach where the classifier operates directly on a continuous and probabilistic signal, in these works the classifier operates on a deterministic signal on a regular grid and cannot reason probabilistically about discretization errors.
Finally, while superficially similar to Deep GPs (Damianou and Lawrence, 2013) or Deep Differential Gaussian Process Flows (Hegde et al., 2018), our PNCNNs tackle fundamentally different kinds of problems like image classification1, and our GPs represent epistemic uncertainty over the values of the feature maps rather than the parameters of the network.

heading: 3 BACKGROUND
text: Probabilistic Numerics: We draw inspiration for our approach from the community of probabilistic numerics where the error in numerical algorithms are modeled probabilistically, and typically with a Gaussian process. In this framework, only a finite number of input function calls can be made, and therefore the numerical algorithm can be viewed as an autonomous agent which has epistemic uncertainty over the values of the input. A well known example is Bayesian Monte Carlo where a Gaussian process is used to model the error in the numerical estimation of an integral and optimally select a rule for its computation (Minka, 2000; Rasmussen and Ghahramani, 2003). Probabilistic numerics has been applied widely to numerical problems such as the inversion of a matrix (Hennig, 2015), the solution of an ODE (Schober et al., 2019), a meshless solution to boundary value PDEs (Cockayne et al., 2016), and other numerical problems (Cockayne et al., 2019). To our knowledge, we are the first to construct a probabilistic numeric method for convolutional neural networks.
Gaussian Processes: We are interested in operating on the continuous function f(x) underlying the input, but in practice we have access only to a collection of the values of that function sampled on a finite number of points {xi}Ni=1. Classical interpolation theory reconstructs f deterministically by assuming a certain structure of the signal in the frequency domain. Gaussian processes give a way of modeling our beliefs about values that have not been observed (Rasmussen et al., 2006), as reviewed in appendix A. These beliefs are encoded into a prior covariance k of the GP f ∼ GP(0, k) and updated upon seeing data with Bayesian inference. Explicitly, given a set of sampling locations x = {xi}Ni=1 and noisy observations y = {yi}Ni=1 sampled yi ∼ N (f(xi), σ2i ), using Bayes rule
1While GPs can be applied directly to image classification, they are not well suited to this task even with convolutional structure baked in, as shown in Kumar et al. (2018).
Input Data
one can compute the posterior distribution f |y,x ∼ GP(µp, kp), which captures our epistemic uncertainty about the values between observations. The posterior mean and covariance are given by
µp(x) = k(x) >[K + S]−1y , kp(x, x ′) = k(x, x′)− k(x)>[K + S]−1k(x′) , (1)
where Kij = k(xi, xj), k(x)i = k(x, xi) and S = diag(σ2i ). Below we shall choose the RBF kernel2 as prior covariance, due to its convenient analytical properties: kRBF(x, x′) = aN (x;x′, l2I) = a ( 2πl2
)− d2 exp(− 12l2 ||x− x′||2). In typical applications of GPs to machine learning tasks such as regression, the function f that we want to predict is already the regression model. In contrast, here we use GPs as a way of representing our beliefs and epistemic uncertainty about the values of both the input function and the intermediate feature maps of a neural network.

heading: 4 PROBABILISTIC NUMERIC CONVOLUTIONAL NEURAL NETWORKS

heading: 4.1 OVERVIEW
text: Given an input signal f : X → Rc, we define a network with layers that act directly on this continuous input signal. We define our neural network recursively from the input f (0) = f , as a series of L continuous convolutions A(`) with pointwise ReLU nonlinearities and weight matrices which mix only channels (known as 1× 1 convolutions) M ∈ Rc×c:
f (`+1) = M (`)ReLU[A(`)f (`)], (2)
and a final global average pooling layer P which acts channel-wise as natural generalization of the discrete case: P(f (L))α = ∫ f (L) α (x)dx for each α = 1, 2, . . . , c. Denoting the space of functions on X with c channels byHc, the convolution operatorsA(`) are linear operators fromHc` toHc`+1 . Like in ordinary convolutional neural networks, the layers build up increasingly more expressive spatial features and depend on the parameters in A(`) and M (`). Unlike ordinary convolutional networks, these layers are well defined operations on the underlying continuous signal.
While it is clear that such a network can be defined abstractly, the exact values of the function f (L) cannot be computed as the operators depend on unknown values of the input. However, by adapting a probabilistic description we can formulate our ignorance of f (0) with a Gaussian process and see how the uncertainties propagate through the layers of the network, yielding a probabilistic output. Before delving into details, we outline the key components of equation 2 that make this possible.
Continuous Convolutional Layers: Crucially, we consider continuous convolution operators A that can be applied to input Gaussian process f ∼ GP(µp, kp) in closed form. The output is another Gaussian process with a transformed mean and covariance Af ∼ GP(Aµp,AkpA′) where A′ acts
2For convenience, we include the additional scale factor (2πl2)d/2 relative to the usual definition.
to the left on the primed argument of kp(x, x′).3 In section 4.2 we show how to parametrize these continuous convolutions in terms of the flow of a PDE and show how they can be applied to the RBF kernel exactly in closed form.
Probabilistic ReLUs: Applying the ReLU nonlinearity to the GP yields a new non Gaussian stochastic process h(`) = ReLU[A(`)f (`)], and we show in section 4.5 that the mean and covariance of this process has a closed form solution which can be computed.
Channel Mixing and Central Limit Theorem: The activations h(`) are not Gaussian; however, for a large number of weakly dependent channels we argue that f (`+1) = M (`)h(`) is approximately distributed as a Gaussian Process in section 4.5.
Measurement and Projection to RBF Gaussian Process: While f (`+1) is approximately a Gaussian process, the mean and covariance functions have a complicated form. Instead of using these functions directly, we take measurements of the mean and variance of this process and feed them in as noisy observations to a fresh RBF kernel GP, allowing us to repeat the process and build up multiple layers without increasing complexity.
The Gaussian process feature maps in the final layer f (L) are aggregated spatially by the integral pooling P which can also be applied in closed form (see appendix D), to yield a Gaussian output. Assembling these components, we implement the end to end trainable Probabilistic Numeric Convolutional Neural Network which integrates a probabilistic description of missing data and discretization error inherent to continuous signals. The layers of the network are shown in figure 1.

heading: 4.2 CONTINUOUS CONVOLUTIONAL LAYERS
text: On a discrete domain such as the lattice X = Zd, all translation equivariant linear operators A are convolutions, a fact which we review in appendix B. In general, these convolutions can be written in terms of a linear combination of powers of the generators of the translation group: the shift operators τi, i = 1, . . . , d shift all elements by one unit along the i-th axis of the grid. For a one dimensional grid, one can always write A = ∑ kWkτ
k where the weight matrices Wk ∈ Rc×c act only on the channels and the shift operator τ acts on functions on the lattice. In d dimensions, A = ∑ k1,...,kd Wk1,...,kdτ k1 1 · · · τ kd d for some set of integer coefficients k1, . . . , kd. For example when d = 2, we can take k1, k2 ∈ {−1, 0, 1} to fill out a 3× 3 neighborhood. On the continuous domain X = Rd we similarly parametrize convolutions with A = ∑ kWke
Dk , where Dk is given by powers of the partial derivatives ∂i, i = 1, . . . , d which generate infinitesimal translations along the i-th axes. Setting d = 1 for simplicity, we can indeed verify that the operator exponential τa = ea∂x applied to a function g(x) is a translation: ea∂xg(x) = g(x) + ag′(x) + 1 2a
2g′′(x)+ · · · = g(x+a), which is the Taylor series expansion of g(x+a) around x. Exponentials of operators can be defined similarly in terms of the formal Taylor series eD = ∑∞ k=0Dk/k! or more broadly as the solution to the PDE:
∂tg(t, x) = (Dg)(t, x) , g(0, x) = g(x) (3)
at time t = 1: eDg(x) = g(t = 1, x).
Following the discussion in the discrete case, translation invariance of Dk imposes that it is expressed in terms of powers of the generators. Collecting the derivatives into the gradient ∇, we can write the general form of Dk as αk + β>k ∇ + 12∇
>Σk∇ + ... for any constants αk, vectors βk, matrices Σk etc. For simplicity, we truncate the series at second order to get
Dk = β>k ∇+ 12∇ >Σk∇ , (4)
where we omit the constants αk that can be absorbed into the definition of Wk. For this choice of D, the PDE in equation 3 is nothing but the diffusion equation with drift βk and diffusion Σk. When discussing rotational equivariance in section 4.4, we also consider a more general form of D. The diffusion layer can also be viewed in another way as the infinitesimal generator of an Ito diffusion (a stochastic process). Given an Ito process with constant drift and diffusion dXt =
3More generally neural networks have affine layers including both convolutions and biases. An affine transformation Af + b of a Gaussian process is also a Gaussian process f ∼ GP(Aµp + b,AkpA′), and we include biases in our network but omit them from the derivations for simplicity.
βdt + Σ1/2dBt where Bt is a d dimensional Brownian motion, the time evolution operator can be written via the Feynman-Kac formula as etDf(x) = E[f(Xt)] where X0 = x. In other words, the operator layerA = etD is the expectation under a parametrized Neural Stochastic Differential equation (Li et al., 2020; Tzen and Raginsky, 2019) that is homogeneous and therefore shift invariant. The flow of this SDE depends on the drift and diffusion parameters β and Σ.
To recap, we define our convolution operator through the general form A = ∑ kWke
Dk where the weight matrices Wk ∈ Rc×c mix only channels and eDk is the forward evolution by one unit of time of the diffusion equation with drift βk and diffusion Σk containing learnable parameters {(Wk, βk,Σk)}Kk=1. The translation equivariance of A follows directly from the fact that the generators commute ∀k, i : [Dk,∇i] = 0 and therefore [A, τi] = 0 (the bracket [a, b] = ab − ba is the commutator of the two operators). In appendix B we show that our definition of A reduces to the usual one in the discrete case and is thus a principled generalization to the continuous domain.

heading: 4.3 EXACT APPLICATION ON RBF GPS
text: Although the application of the linear operator A = ∑ kWke
Dk involves the time evolution of a PDE, owing to properties of the RBF kernel we fortuitously can apply the operator to an input GP in closed form! Gaussian processes are closed under linear transformations: given f ∼ GP(µp, kp), we need only compute the action of A on the mean and covariance: Af ∼ GP(Aµp,AkpA′), where A′ is the adjoint w.r.t. the L2(X ) inner product. The application of time evolution eDk is a convolution with a Green’s function Gk, so Af = ∑ kWke Dkf = ∑ kWkGk ∗ f . As we derive in appendix C, the Green’s function for Dk = β>k ∇ + (1/2)∇>Σk∇, is nothing but the multivariate Gaussian density Gk(x) = N (x;−βk,Σk):
Af = ∑ k Wke Dkf = ∑ k WkGk ∗ f = ∑ k WkN (−βk,Σk) ∗ f . (5)
In order to apply etD to the posterior GP, we need only to be able to apply the operator to the posterior mean and covariance. This posterior mean and covariance in equation 1 are expressed in terms of kRBF = aN (x;x′, `2I) and the computation boils down to a convolution of two Gaussians:
etDkRBF(x, x ′) = N (x;−tβ, tΣ) ∗ aN (x;x′, `2I) = aN (x;x′ − tβ, `2I + tΣ) (6)
etD1kRBF(x, x ′)etD ′ 2 = aN (x;x′ − t(β1 − β2), `2I + tΣ1 + tΣ2) . (7)
The application of the channel mixing matrices Wk and summation is also straightforward through matrix multiplication for the mean and covariance. To summarize, because of the closed form action on the RBF kernel, the layer can be implemented efficiently and exactly with no discretization or approximations.
We note that with the Green’s function above, the action ofA encompasses the ordinary convolution operator on the 2d lattice as a special case. Given drift βk ∈ {−1, 0, 1}×2, k = 1, . . . , 9 filling out the 9 elements of a 3× 3 grid and as the diffusion Σk → 0, the Green’s function is a Dirac delta, so that: Af(x) = ∑ kWkδ(x− βk) ∗ f(x) = ∑ i,j=−1,0,1Wijf(x1 − i, x2 − j) = W ∗Z2 f(x).

heading: 4.4 GENERAL EQUIVARIANCE
text: The convolutional layers discussed so far are translation equivariant. We discuss how to extend the continuous linear operator layers to more general symmetries such as rotations. Feature fields in this more general case are described by tensor fields, where the symmetry group acts not only on the input space X but also on the vector space attached to each point x ∈ X . A linear layer A is equivariant if its action commutes with that of the symmetry. In appendix E we derive constraints for general linear operators and symmetries, which generalize those appearing in the steerableCNN literature (Weiler and Cesa, 2019; Cohen et al., 2019). Then we show how equivariance under continuous roto-translations in 2d constrains the form of a convolutional layer by solving the equivariance constraint. Non-trivial solutions require that the operator D in the PDE of equation 3 has a non-trivial matrix structure.

heading: 4.5 PROBABILISTIC NONLINEARITIES AND RECTIFIED GAUSSIAN PROCESSES
text: Gast and Roth (2018) derive the mean and variance for a univariate rectified Gaussian distribution for use in a neural network. We generalize these results to the full covariance function (and higher moments) of a rectified Gaussian process in appendix J and present the results here. For the input GP A(`)f (`)(x) ∼ GP(µ(x), k(x, x′)), we denote σ(x) = √ k(x, x), Σ the matrix with components Σij = k(xi, xj) for i, j = 1, 2 and µ = [µ(x1), µ(x2)]. We use notation Φ(z) for the univariate standard normal CDF, and Φ(z; Σ) for (two dimensional) multivariate CDF of N (0,Σ) at z. Σ1 and Σ2 are the column vectors of Σ. The first and second moments of h = ReLU[Af ] are:
E[h(x)] = µ(x)Φ(µ(x)/σ(x)) + σ(x)Φ′(µ(x)/σ(x)) , (8) E[h(x1)h(x2))] = (k(x1, x2) + µ(x1)µ(x2))Φ(µ; Σ) (9)
+ (µ(x1)Σ > 2 + µ(x2)Σ > 1 )∇Φ(µ; Σ) + Σ>1 ∇∇>Φ(µ; Σ)Σ2,
where∇∇>Φ denotes the Hessian of Φ with respect to the first argument. The first and higher order derivatives of the Normal CDF are just the PDF and products of the PDF with Hermite polynomials. Note that the mean and covariance interact through the nonlinearity.

heading: 4.6 CHANNEL MIXING AND CENTRAL LIMIT THEOREM
text: After the non-linearity the process is no longer Gaussian. To overcome this issue we introduce a channel mixing matrix M (`) ∈ Rc`+1×c` and define the feature map in the following layer by f (`+1) = M (`)h(`), where h(`) = ReLU[A(`)f (`)]. So long as the channels of h(`) are only weakly dependent, we can apply the central limit theorem (CLT) to each function f (`+1)α = ∑c` β=1M (`) αβh (`) β so that in the limit of large c`, the statistics of the f (`+1) α ’s converge to a GP with first and second moments given by:
E[f (`+1)(x)] = ME[h(`)(x)], E[f (`+1)(x)f (`+1)(x′)>] = ME[h(`)(x)h(`)(x′)>]M> . (10)
The convergence to a Gaussian process here is reminiscent of the well known infinite width limit of Bayesian neural networks Neal (1996); de G. Matthews et al. (2018); Yang (2019). However the setting here is fundamentally different. Unlike the Bayesian case where the distribution of M is given by a prior or posterior, in PNCNN M is a deterministic quantity and instead the uncertainty is about the input. PNCNN is not a Bayesian method in the sense of representing uncertainty about the parameters of the model, but instead it is Bayesian in representing and propagating the uncertainty in the value of the inputs. We go through the central limit theorem argument in appendix I, and go over some potential caveats such if the weights were to converge to 0 or the independence assumptions are violated. In the following section, we evaluate the Gaussianity empirically.

heading: 4.7 MEASUREMENT AND PROJECTION TO RBF GAUSSIAN PROCESS
text: As a last step we simplify the mean and covariance functions of the approximate GP f (`+1). While we can readily compute the values of these functions, unlike in the RBF kernel case, we cannot apply the convolution operator etD in closed form. In order to circumvent this challenge, we model the (approximately) Gaussian process f (`+1) with an RBF Gaussian process as follows: we evaluate the mean yi = E[f (`+1)(xi)] and variance σ2i = Var[f (`+1)(xi)] of the approximate Gaussian process f (`+1) at a collection of points {xi}Ni=1 using equations 8, 9 and 10. These values yi are treated as measurements of the underlying signal with a heteroscedastic noise σ2i that varies from point to point. We can then compute the RBF-based posterior GP of this signal f̂ (`+1)|{(xi, yi, σi)}Ni=1 ∼ GP(µp, kp) with posterior mean and covariance given by equation 1 for the heteroschedastic noise model. The uncertainty in the input f (`) is propagated through to the RBF posterior f̂ (`+1)|{(xi, yi, σi)}Ni=1 via the measurement noise σi. Crucially, this Gaussian process mean and covariance functions are written in terms of the RBF kernel and we can therefore continue applying convolutions in closed form in future layers.
As we describe in the following section, the RBF kernel in each layer is trained to maximize the marginal likelihood of the data that it sees, and thereby minimize the discrepancy with the underlying generating distribution f (`+1). While this measurement/projection approach is effective in
many scenarios, in networks with many layers or a very large number of observations uncertainty information can get attenuated as it passes through the layers, a phenomenon which we investigate in appendix H.
With a network that is trained on a version of MNIST that is randomly subsampled to 75 pixels, in figure 2 (left) we evaluate the mean and uncertainty of the internal feature maps as we vary the number of pixels of the inputs at test time. As expected, the mean functions for the feature maps slowly converge and the predicted uncertainties decrease in magnitude as the input resolution is increased. In figure 2 (middle) we show that in early layers the uncertainties decrease at a similar rate to the O(1/ √ N) of discretization error that we would expect from a standard convolutional layer which is discretized to a square grid.4 Despite the fact that these resolutions differ substantially from those seen at training time and the fact that there are no explicit uncertainty targets for these internal layers, the predictions are reasonably well calibrated as demonstrated in figure 2 (right). While the prediction residuals have fatter tails than a standard Gaussian (as evidenced by the deviations from the line in the QQ plot), the mean and standard deviation are close to the theoretically optimal 0 and 1 values across a range of resolutions.

heading: 4.8 TRAINING PROCEDURE
text: Our neural network has two sets of parameters: the channel mixing and diffusion parameters, {(M (`),W (`),β(`),Σ(`))}L`=1, as well as kernel hyperparameters of the Gaussian Processes {(l(`), a(`))}L`=1. We train all parameters jointly on the loss Ltask + λLGP , where Ltask is the cross entropy with logits given by the mean µP of the pooled features P(f (L)) ∼ N (µP ,ΣP ) and LGP is the marginal log likelihoods of the GP feature maps:
LGP(f) = 1
2 L∑ `=1 c∑̀ α=1 [( fTα [KXX + Sα] −1fα ) + log det [KXX + Sα] +N log 2π ](`) , (11)
where for each layer `, fα = [ fα(x1), ..., fα(xN ) ] ∈ RN are the observed values for channel α at locations X = [x1, . . . , xN ], KXX is the covariance of the RBF kernel and Sα = diag(σ2α) the measurement noise for each channel α and spatial location. Notably the GP marginal likelihood is independent of the class labels.

heading: 5 EXPERIMENTAL RESULTS
text: We evaluate the Probabilistic Numeric CNN on two different problems which have incomplete and irregular observations.
4A 2D discrete convolution layer usingN = m2 points can be interpreted as a Riemann sum approximation of the continuous integral and will therefore have an error of O(1/m) = O(1/ √ N) which is the same rate as would be achieved through Monte Carlo sampling.
Superpixel MNIST is an adaptation of the MNIST dataset where the 784 pixels of the original images are replaced by 75 salient superpixels that are non uniformly spread throughout the domain and are different for each image (Monti et al., 2017). Despite the simplicity of the underlying images, the lack of grid structure and high fraction of missing values make this a challenging task. Example inputs are visualized at the left of figure 1. We compare to Monet (Monti et al., 2017), SplineCNN (Fey et al., 2018), Graph Convolutional Gaussian Processes (GCGP) (Walker and Glocker, 2019), and Graph Attention Networks (GAT) (Avelar et al., 2020). As shown in table 5, the probabilistic numeric CNN greatly outperforms the competing methods reducing the classification error rate by more than 3× over the previous state of the art.
We conduct an ablation study where uncertainty propagation is removed: the probabilistic ReLU is replaced with the deterministic one applied to the mean, and the uncertainties in each layer are set to 0. This form of the network still makes use of the fact that the input function is continuous by use of the GP interpolation of the means, but crucially it does not integrate the uncertainty in the computation resulting from the missing data. While this variant (PNCNN w/o σ) with a 3.03% error rate outperforms existing methods from the literature, it is substantially worse than the PNCNN that integrates the uncertainty at 1.24% error. This validates both that the underlying architecture (using the continuous convolution operators) has good inductive biases and that reasoning about discretization errors probabilistically can improve performance directly.
In figure 3 we evaluate the performance on zero shot generalization to a different test resolution for a variety of training resolutions. In order to compare to an ordinary CNN we sample MNIST on a regular square grid: PNCNN with uncertainty is the most robust to this train test sampling distribution shift, followed by PNCNN w/o uncertainty, and finally the ordinary CNN which is quite sensitive to these changes.
Irregularly Spaced Time Series For the second task, we evaluate our model on the irregularly spaced time series dataset PhysioNet2012 (Silva et al., 2012) for predicting mortality from ICU vitals signs. This dataset is particularly challenging because different vital sign channels are observed at different times, even within a single patient record. This means that we cannot compute the GP inference formula of equation 1 efficiently for all channels simultaneously because the observation
points {xi} and hence the matrices K in that formula differ between the channels, increasing computational complexity. To circumvent this difficulty, we employ a stochastic diagonal estimator to compute the variances as described in appendix G. We compare against IP-Nets (Shukla and Marlin, 2019), SEFT-ATTN (Horn et al., 2019), and GRU-D (Che et al., 2018) as reported in Horn et al. (2019). PNCNN performs competitively, although not a breakout performance as in the image dataset which we attribute to the use of the stochastic variance estimates over an exact calculation.

heading: 6 CONCLUSION
text: Based on the ideas of probabilistic numerics, we have introduced a new class of neural networks which model missing values and discretization errors probabilistically within the layers of a CNN. The layers of our network are a series of operators defined on continuous functions, removing dependence on shortcut features like the sampling locations and distribution. On irregularly sampled and incomplete spatial data we show improved generalization and robustness. We believe that generative modeling of the internal feature maps of a neural network for representing uncertainty from incomplete observations is an underexplored research direction and has significant promise for other modalities of data, which we leave to future work.

heading: A REVIEW OF GAUSSIAN PROCESSES
text: We briefly review here the main ideas of Gaussian Processes for machine learning, see Stein (2012); Rasmussen et al. (2006) for more details. We start to explain how to use stochastic processes for Bayesian inference. We see the stochastic process as prior over functions p(f) and as we are given samples x = (x1, . . . , xN ),y = (y1, . . . , yN ), yi ≡ f(xi), we update our beliefs about the function by constructing the posterior via Bayes rule: p(f |y,x) = p(y|f,x)p(f)/p(y|x). Here we need to specify the likelihood of the data with our model p(y|f,x) = ∏N i=1 p(yi|f, xi) and the denominator, called the evidence or marginal likelihood, follows: p(y|x) = Ef∼p(f)[p(y|f,x)]. The power of this approach is that the value of the signal y at an unseen point x has an uncertainty which depends on our knowledge of its neighbourhood: p(y|x,y,x) = Ef∼p(f |y,x)[p(y|f, x)] and allows us to reason probabilistically about the underlying signal.
A particular convenient class of random function is Gaussian processes (GPs) for which inference can be done exactly. A stochastic process can be presented in terms of the finite distributions of the random variables {f(xi)}Mi=1 at points {xi}Mi=1. For a GP these distributions are Gaussian and can be defined uniquely by specifying means and covariances, and so a GP is specified entirely by its mean function µ(x) and covariance kernel k(x, x′). We shall write f ∼ GP(µ, k). Let us assume a Gaussian likelihood model as well, i.e. p(yi|f, xi) = N (f(xi), σ2i ), where σn represents aleatoric uncertainty on the measurement. (For simplicity we take here the function to be scalar valued but the reasoning can be easily generalized.) Then properties of the Gaussian distribution (see (Rasmussen et al., 2006, Chap. 2) for a detailed derivation of the formulas) lead to the following posterior distributions after seeing data y,x: p(f |y,x) = GP(µp, kp), with
µp(x) = k(x) T [K + S]−1y , kp(x, x ′) = k(x, x′)− k(x)T [K + S]−1k(x′) . (12)
where Kij = k(xi, xj), k(x)i = k(x, xi) and S = diag(σ2i ).

heading: B FROM DISCRETE TO CONTINUOUS CONVOLUTIONAL LAYERS
text: We here show that the general formula A = ∑ k Wke Dk (13)
with Dk a function of spatial derivatives, reduces in the case of discrete input space X to the usual convolution we encounter in deep learning.
For simplicity we shall assume a 1d grid as input space X = {1, . . . , N}. Let us start by recalling the form of the classical discrete convolution when C` = C`+1 = 1. We define a convolutional layer as a linear map that commutes with the translation operator. To make the symmetry exact, we need assume periodic boundaries. Then in the standard basis of RN , {ei}Ni=1 of vectors localized at site i, the translation operator τ acts as τei = ei+1 mod N . An N ×N matrix B is translation invariant iff τB = Bτ . Since τ is diagonal in Fourier space, the most general solution is B = F−1diag(b̂)F , where Fjk = e− 2πi N jk is the discrete Fourier transform. Such matrices are called circulant and can
be written alternatively as B = ∑N−1 i=0 bN−iτ i, b = F−1b̂. Explicitly:
B =  b0 b1 . . . bN−2 bN−1 bN−1 b0 b1 bN−2 ... bN−1 b0 . . . ...
b2 . . . . . . b1 b1 b2 . . . bN−1 b0
 . (14)
This shows that the most general convolutional layer is a circulant matrix. E.g. if bi = 0 unless i = 0, 1, N −1, B coincides with the matrix representing a periodic convolution of filter size 3. The matrix B is invertible as long as b̂k 6= 0 for all k. In a convolutional network the parameters bi are random variables and the measure of the set where B is not invertible is zero. Thus the role of eD is replaced in the discrete case by the B.
The discrete analog of A is then:
A = ∑ i Wi ⊗Bi . (15)
Introducing the unit matrices Eα,β which have 1 at the row α and column β and 0 otherwise, we can rewrite it as:
A = ∑ j,α,β Eα,β ⊗ τ jWα,βj , W α,β j = ∑ i Wα,βi bi,N−j . (16)
Since Eα,β⊗ τ j is a linear basis of the space of convolutional layers, we see that equation 13 indeed reduces to the usual one when discretizing the input domain and is a principled generalization to the continuous domain.

heading: C GREENS FUNCTION
text: Given the operator D = β>∇ + 12∇ >Σ∇ we can compute the action of etD in terms of convolutions. Using the d dimensional Fourier transforms F [h](k) = (2π)− d2 ∫ h(x)e−ik
>xdx and F−1 = F†, we can rewrite the derivative operator D in terms of elementwise multiplication in the Fourier domain, which diagonalizes D. Since ∇ = F−1(ik)F ,
D = F−1(iβ>k − 12k >Σk)F . (17)
Using the series definition etD = ∑∞ n=0(tD)n/n!, we have:
etD = F−1et(iβ >k− 12k >Σk)F . (18)
Applying this operator to a test function h(x) yields
etDh = F−1[et(iβ >k− 12k >Σk)F [h](k)] = F−1[F [Gt] · F [h]] = Gt ∗ h, (19)
where the final step follows from the Fourier convolution theorem, and we define the function Gt = F−1[et(iβ >k− 12k >Σk)]. Directly applying the Fourier integral yields a Gaussian integral
Gt(x) = (2π) − d2 ∫ eik >(x+tβ)− 12k >tΣkdk = e− 1 2 (x+tβ) >(tΣ)−1(x+tβ)det(2πtΣ)−1/2. (20)
This function Gt(x) = N (x;−tβ, tΣ) is nothing but a multivariate heat kernel, the Greens function (also known as the fundamental solution or time propagator) for the diffusion equation ∂tGt(x − x′) = DGt(x− x′), and indeed limt→0Gt(x− x′) = δ(x− x′).
D INTEGRAL POOLING The integral pooling operator P[f ] = ∫ Rd f(x)dx can be applied to the Gaussian process just like any other linear operator. Given f (L) ∼ GP(µ, k), we have that
Pf (L) ∼ GP(Pµ,PkP ′) = N (Pµ,PkP ′). (21)
Again, computing the mean µP = Pµ and covariance matrix ΣP = PkP ′ we need just to be able to apply P to the RBF kernel.
PkRBF(x′) = ∫ Rd kRBF(x, x ′)dx = a (22)
PkRBFP ′ = ∫ Rd×Rd kRBF(x, x ′)dxdx′ =∞ (23)
For many applications such as image classification using the mean logit value, we require only the predictive mean, so an unbounded covariance matrix ΣP is acceptable. We use this form for all of our experiments.
However for some applications an output uncertainty can be useful, so we also provide a variant that integrates over a finite region [0, 1]d, Pf = ∫ [0,1]d f(x)dx.
PkRBF(x′) = ∫
[0,1]d kRBF(x, x
′)dx = a d∏ i=1 [ Φ( x′i ` )− Φ( x′i−1 ` ) ]
(24)
PkRBFP ′ = ∫
[0,1]d×[0,1]d kRBF(x, x
′)dxdx′ = a [ ` √
2 π (e −1/2`2 − 1) + 2Φ( 1` )− 1
]d (25)
where Φ is again the univariate standard normal CDF.

heading: E EQUIVARIANCE
text: E.1 RELATED WORK
We note that there has been considerable research effort in the development of equivariant CNNs which we build on top of. The group equivariant CNN was introduced by Cohen and Welling (2016a) for discrete groups on lattices. This work has been extended for continuous groups (Worrall et al., 2017; Zhou et al., 2017) and with steerable equivariance (Cohen and Welling, 2016b; Weiler and Cesa, 2019) where other group representations are used. There have also been group equivariant networks designed for point clouds and other irregularly spaced data (Thomas et al., 2018; Finzi et al., 2020; Fuchs et al., 2020; de Haan et al., 2020). In Shen et al. (2020), layers using finite difference estimation of derivative operators are used for defining equivariant layers in an equivariant CNN, effectively a change of basis.
Most closely related to our PDE operator approach to equivariance is work by Smets et al. (2020). In this work, the authors define layers of their convolutional network through the time evolution of a PDE which is a nonlinear generalization of the diffusion equation, and which additionally includes pooling like behaviour. Smets et al. (2020) parametrize their equivariant PDEs in terms of left invariant vector fields. These PDEs are applied to scalar functions on the group, unlike our work where feature fields (of various representations) live on the base space Rn. These two approaches reflect the two major approaches to equivariance in the literature: using regular representations and using irreducible/tensor representations. Both have advantages and disadvantages and we focus on the second approach in this work. It may be possible that there is a deeper relationship between the solutions we find and the left invariant vector fields in Smets et al. (2020), but this is beyond the scope of our work.
E.2 TRANSLATION EQUIVARIANCE
A key factor in the generalization of convolutional neural networks is their translation equivariance. Patterns in different parts of an input signal can be seen in the same way because convolution is translation equivariant. Our learnable linear operators A are equivariant to continuous transformations. Two linear operators eC and eB commute [eC , eB] = 0 if and only if their generators commute: [C,B] = 0. Since the generator of diffusions Di is a sum of derivative operators, and the generators of translations are just ∇ as mentioned in section 4.2, the two commute: [Dk,∇] = 0. Therefore [A, τa] = [ ∑ kWke Dk , τa] = ∑ kWk[e Dk , ea >∇] = 0 and A is translation equivariant.
E.3 STEERABLE EQUIVARIANCE FOR LINEAR OPERATORS
For some tasks like medical segmentation, aerial imaging, and chemical property prediction there are additional symmetries in the data it makes sense to exploit other than mere translation equivariance. Below we show how to enforce equivariance of the Linear operator A to other symmetry groups G such as the group of continuous rotations SO(d) in Rd. Applying equivariance constraints separately on each of the components of eDi on top of translation equivariance yields very restricted set of operators. For example, enforcing equivariance to continuous rotations G = SO(d), the operator must be an isotropic heat kernel: Dk = ck∇>∇. The reason for this apparent restriction is a result of considering the different channels independently, as scalar fields.
The alternative is to use features fields which transform under more general representations of the symmetry group, introduced in steerable-CNNs (Cohen and Welling, 2016b) and used in (Worrall
et al., 2017; Thomas et al., 2018; Weiler et al., 2018; Weiler and Cesa, 2019) and others. In this way, the symmetry transformation acts not only on the spatial domain X , but also transforms the channels. The way that the group acts on Rc (i.e. the channels) is formalized by a representation matrix ρ(g) ∈ Rc×c for each element g ∈ G in the transformation group that satisfies ∀g, h ∈ G : ρ(gh) = ρ(g)ρ(h). Choosing the type of each intermediate feature map is equivalent to choosing their representations, and we describe a simple way of doing this with tensor representations in the later section.
Operator Equivariance Constraint: Returning to linear operators, we derive the equivariance constraint and show how to use constructs from the previous sections to implement steerable rotation equivariance. Equivariance of a linear operatorA : (Rd → Rcin)→ (Rd → Rcout) requires that, for any input function, transforming the input function first (both argument and channels) and applying A is equivalent to first applying A and then transforming the output: Aρin(g)Lgf = ρout(g)LgAf where Lgf(x) = f(g−1x). Rearranging the terms, one sees that the equivariance constraint on the linear operator A is: ρout(g)LgALg−1ρin(g−1) = A, (26) where the operators Lg and L−1g are understood not to act on the representation matrices ρ (although implicitly a function of g). As shown in Appendix E.4, eq. 26 is a direct generalization of the equivariance constraint for convolutions ∀x : ρout(g)K(g−1x)ρin(g−1) = K(x) described in the literature (Weiler and Cesa, 2019; Cohen et al., 2019).
As shown in Appendix E.6, the equivariance constraint for continuous rotations applied to the diffusion operators A = ∑ kWke
Dk has only the trivial solutions of isotropic diffusion without any drift. For this reason we instead consider a more general form of diffusion operator where the PDE itself couples the different channels. For the coupled PDE:
∂f ∂t = ∑ k WkDkf (27)
the time evolution contains the matrices Wk in the exponential A = e ∑ kWkDk . Like with the
example of translation above, this operator is equivariant if and only if the infinitesmal generator∑ kWkDk is equivariant.
Because equation 26 applies generally to linear operators and not just convolutions, we can compute the equivariance constraint for these derivative operators. We can simplify the summation∑ kWkDk = ∑ kWk(β T k ∇ + (1/2)∇TΣk∇) by writing it in terms of the collection of matrices
Bi = ∑ kWkβki and Sij = (1/2) ∑ kWkΣkij to express Aderiv = ∑ iBi∂i + ∑ i,j Sij∂i∂j where the indices i, j = 1, 2..., d enumerate the spatial dimensions of each vector βk and each matrix Σk. As we derive in appendix E.5, the necessary and sufficient conditions for the equivariance of ∑ kWkDk and therefore A is that ∀g ∈ G : [ρout ⊗ ρ∗in ⊗ ρ(1,0)](g)vec(B) = vec(B) and ∀g ∈ G : [ρout⊗ρ∗in⊗ρ(2,0)](g)vec(S) = vec(S) where vec(·) denotes flattening the elements into a single vector and ρ(r,s) is the tensor representation with r covariant and s contravariant indices.
E.4 GENERALIZATION OF EQUIVARIANCE CONSTRAINT FOR CONVOLUTIONS
This equivariance constraint is a direct generalization of the equivariance constraint for convolution kernels as described in Weiler and Cesa (2019); Cohen et al. (2019). In fact, whenA is a convolution operator,Af = K ∗f , the action of Lg by conjugationA is equivalent to transforming the argument of the kernel K:
Lg(K∗)Lg−1f(x) = ∫ K(g−1x− x′)f(gx′)dµ(x′)
= ∫ K(g−1(x− x′′))f(x′′)dµ(x′′) = (Lg[K]) ∗ f.
Letting both sides of eq 26 act on the product of a constant unit vector ei and a delta function, f = eiδ the expression ∀ei : ρout(g)Lg[K]ρin(g−1) ∗ eiδ = K ∗ eiδ can be rewritten as ∀x : ρout(g)K(g −1x)ρin(g −1) = K(x) which is precisely the constraint for steerable equivariance for convolution described in the literature. 5
5This assumes as is typically done that measure µ over which the convolution is performed is left invariant. For the more general case, see the discussion in Bekkers (2019).
E.5 EQUIVARIANT DIFFUSIONS WITH MATRIX EXPONENTIAL
Below we solve for the necessary and sufficient conditions for the equivariance of the operator Aderiv. We will use tensor representations for their convenience, but the approach is general to allow other kinds of representations. A rank (p, q) tensor t is an element of the vector space T(p,q) := V ⊗p ⊗ (V ∗)⊗q where V is some underlying vector space, V ∗ is its dual and (·)⊗p is the tensor product iterated p times. In common language T(0,0) are scalars, T(1,0) are vectors, and T(1,1) are matrices. Given the action of a group G on the vector space V , the representation on T(p,q) is ρ(p,q)(g) = g
⊗p ⊗ (g−>)⊗q where −> is inverse transpose and ⊗ on the matrices is the tensor product (Kronecker product) of matrices. Composite representations can be formed by stacking different tensor ranks together, such as a representation of 50 scalars, 25 vectors, 10 matrices and 5 higher order tensors: T 50(0,0) ⊕ T 25 (1,0) ⊕ T 10 (1,1) ⊕ T 5 (1,2), where ⊕ in this context is the same as
the Cartesian product. For a composite representation U = ⊕
i T(pi,qi) the group representation is similarly ρU (g) = ⊕ i ρ(pi,qi)(g) where ⊕ concatenates matrices as blocks on the diagonal.
Noting that the operator Lg that acts only on the argument and the matrix ρin(g) acts only on the components, the two commute and we can rewrite the constraint for Aderiv as∑
i
ρout(g)Biρin(g −1)Lg∂iLg−1 + ∑ ij ρout(g)Sijρin(g −1)Lg∂i∂jLg−1 = Aderiv (28)
We can simplify the expression Lg∂iLg−1 by seeing how it acts on a function. For any differentiable function ∂iLg−1f(x) = ∂∂xi [f(gx)] = ∑ j gji[∂jf ](gx) = Lg−1 ∑ j gji∂jf(x) where gij are the components of the matrix g. Since this holds for any f , we find that Lg∇Lg−1 = gT∇ and therefore Lg∇∇TLg−1 = Lg∇Lg−1Lg∇TLg−1 = gT∇∇T g. Since equation 28 holds as an operator equation, it must be true separately for each component ∂i and ∂i∂j . This means that the constraint separates into a constraint for B and a constraint for S:
1. ∀g, i : ∑ j gijρout(g)Bjρin(g −1) = Bi
2. ∀g, i, j : ∑ kl gi`gikρout(g)S`kρin(g −1) = Sij .
These relationships can be expressed more succinctly by flattening the elements of B and S into vectors: [ρout(g) ⊗ ρin(g−T ) ⊗ ρ(1,0)(g)]vec(B) = vec(B) and [ρout(g) ⊗ ρin(g−T ) ⊗ ρ(2,0)(g)]vec(S) = vec(S).
E.6 ROTATION EQUIVARIANCE CONSTRAINT FOR SCALAR DIFFUSIONS HAS ONLY TRIVIAL SOLUTIONS
The diffusion operator A = ∑ kWke
Dk leads to only trivial βk = 0 and Σk ∝ I if it satisfies the continuous rotation equivariance constraint.
Proof:
The application of eDk is just a convolution with the Greens function∑ k Wke Dkf = ∑ k Wk[e − 12 (x+βk) >Σ−1k (x+βk)det(2πΣk) −1/2] ∗ f = ∑ k WkGk ∗ f (29)
where the Greens function is the multivariate Gaussian density: Gk(x) = N (x;−βk,Σk). As shown in appendix E.4, for convolutions the operator constraint is equivalent to the kernel equivariance constraint ρout(g)K(g−1x)ρin(g−1) = K(x) from (Weiler and Cesa, 2019). With K(x) = ∑ kWkGk(x) this reads:
∀x ∈ Rd, g ∈ G : ∑ k ρout(g)WkN (g−1x;−βk,Σk)ρin(g−1) = ∑ k WkN (x;−βk,Σk),
For rotations g ∈ SO(2) where we can parametrize gθ = eθJ in terms of the antisymmetric matrix J = [[0, 1], [−1, 0]] ∈ R2×2 and the translation operator can be written Lg = e−θx T JT∇, we can
take derivatives with respect to θ to get (now with double sums implicit): ∀x ∈ Rd : ∑ k [ dρoutWkN (x;−βk,Σk)−WkN (x;−βk,Σk)dρin−Wk(xTJT∇)N (x;−βk,Σk) ] = 0.
Here the Lie Algebra representation of J is dρ := ∂∂θρ(gθ)|θ=0. Factoring out the normal density: ∀x ∈ Rd : ∑ k [ dρoutWk −Wkdρin −Wk(xTJTΣ−1k (x+ βk)) ] N (x;−βk,Σk) = 0.
Without loss of generality we may assume that each of the Gaussians βk,Σk pairs are distinct since if they were not then we could replace the collection with a single element. Since the (finite) sum of distinct Gaussian densities is never a Gaussian density, and monomials of order > 0 multiplied by a Gaussian density cannot be formed with sums of Gaussian densities or sums multiplied by monomials of a different order and Gaussian densities are never 0, this constraint separates out into several independent constraints.
1. ∀i : dρoutWk = Wkdρin 2. ∀i, x : Wk(xTJTΣ−1k βk) = 0 3. ∀i, x : Wk(xTJTΣ−1k x) = 0
We may assume w.l.o.g. that Wk is not 0 for all components of the matrix (otherwise we could have deleted this element of k and continue). Therefore there is some component which is nonzero, and the expressions in parentheses in equations 2 and 3 must be 0. Given that this holds for all x, eq 3 implies: JTΣ−1k = 0 or equivalently Σ −1 k J = 0 because Σk is symmetric, and since J = −JT this can be expressed concisely as [Σ−1k , J ] = 0 for which the only symmetric solution is proportional to the identity Σk = ckI . Since both Σk and J are invertible, equation 2 yields βk = 0. Therefore there are no nontrivial solutions for β,Σ in A = ∑ kWke Dk for continuous rotation equivariance.

heading: F DATASET AND TRAINING DETAILS
text: In this section we elaborate on some of the details regarding hyperparameters, network architecture, and the datasets.
As described in the main text, the PNCNN is composed of a chain of convolutional blocks containing a convolution layer, a probabilistic ReLUs, and linear channel mixing layer (analogue of the colloquial 1 × 1 convolution). In each of these convolutional blocks, the input is a collection of points and feature mean value at those points along with the feature elementwise standard deviation at those points: {(xi, µ(xi), σ(xi)}Ni=1. These observations seed the GP layer, and the block is evaluated at the same collection of points for the output (although it can be evaluated elsewhere since it is a continuous process, and we make use this fact to visualize the features in figures 1 and 2).
Hyperparameters: For the PNCNN on the Superpixel MNIST dataset, we use 4 PNCNN convolution blocks with c = 128 channels and with K = 9 basis elements for the different drift and diffusion parameters in ∑K k=1Wke
Dk . We train for 20 epochs using the Adam optimizer (Kingma and Ba, 2014) with lr = 310−3 with batch size 50.
For the PNCNN on the PhysioNet2012 dataset, we use the variant of the PNCNN convolution layer that uses the stochastic diagonal estimator described in appendix G with P = 20 probes. In the convolution blocks we use c = 96 channels, K = 5 basis elements and we train for 10 epochs using the same optimizer settings above. For both datasets we tuned hyperparameters on a validation set of size 10% before folding the validation set back into the training set for the final runs. Both models take about 2 hours to train.
SuperPixel-MNIST We source the SuperPixel MNIST dataset (Monti et al., 2017) from Fey and Lenssen (2019) consisting of 60k training examples and 10k test represented as collections of positions and grayscale values {(xi, f(xi))}75i=1 at the N = 75 super pixel centroids. PhysioNet2012 We follow the data preprocessing from Horn et al. (2019) and the 10k-2k train test split. The individual data points consist of 42 irregularly spaced vital sign time series signals as
well as 5 static variables: Gender, ICU Type, Age, Height, Weight. We use one hot embeddings for the first two categoric variables, and we treat each of these static signals as fully observed constant time series signals. As the binary classification task exhibits a strong label imbalance, 14% positive signals, we apply an inverse frequency weighting of 1/.14 to the binary cross entropy loss.

heading: G STOCHASTIC DIAGONAL ESTIMATION FOR PHYSIONET2012
text: In order to compute the mean and variance of the rectified Gaussian process, the activations of the probabilistic ReLU, we need compute the diagonal of AkpA′(xn, xn) for the relevant points {xn}Nn=1. In the usual case where each of the channels α = 1, 2, ..., c are observed at the same locations this can be done efficiently. First one computes the application of eDi on the left and eD ′ j on the right onto the posterior kp:
Nij = (e Dikpe D′j )(xn, xn) = (e DikeD ′ j )(xn, xn)− (eDik>)(xn)[K + S]−1(keD ′ j )(xn)
where k is the RBF kernel and we have reused the notation from appendix A. Notably, this quantity is the same for each of the channels, and the elementwise variance is just:
vα(xn) = (AkpA′)αα(xn, xn) = ∑ i,j,β Wαβi NijW αβ j (30)
where the α, β index the channels of each of the matrices Wi. Because N is the same for all channels, we can compute this quantity efficiently with a reasonable memory cost and compute.
For the PhysioNet2012 dataset where the observation points differ between the channels we must consider a different observation set {xβn}Nn=1 for each channel β. This means that evaluated kernel depends on the channel and we have the objects: kβ , Kβ and Sβ . As a result, we have an additional index for Nβij and the desired computation is
vα(x α n) = (AkpA′)αα(xαn, xαn) = ∑ i,j,β Wαβi N β ijW αβ j . (31)
While each of the terms in the computation can be computed without much difficulty, performing the summation explicitly requires an unreasonably large memory cost and also compute.
However, by the same approach we can consider the full covariance matrix B(αn)(βm) = (AkpA′)αβ(xαn, xβm), and while it would not be feasible to compute this matrix directly we can define matrix vector multiplies onto vectors of size RcN implicitly using the sequence of operations that define it. Crucially, this sequence of operations has much more modest memory consumption (and compute cost) over the direct expression in equation 31. These implicit matrix vector multiplies can then be used to compute a stochastic diagonal estimator (Bekas et al., 2007) given by:
v̂α(x α n) = 1 P P∑ p=1 zp Bzp (32)
with Gaussian probe vectors zp ∼ N (0, I), and where is elementwise multiplication (see Bekas et al. (2007) for more details on this stochastic diagonal estimator). We use this estimator with P = 20 probes for computing the variances for PhysioNet. We note that with P = 20 the variance estimates are still quite noisy, however without the estimator cannot readily apply the PNCNN to PhysioNet. We leave a better approach for handling this kind of data to future work.

heading: H PATHOLOGIES IN PROJECTION TO RBF GAUSSIAN PROCESS
text: In section 4.7 describe an approach by which a Gaussian process with a complex mean and covariance function is projected down to the posterior of a (simpler) RBF kernel GP from a set of observations. We know given the representation capacity of the RBF kernel that with the right set of observations, a complex function can be well approximated in principle. However, the relationship for uncertainty is less straightforward.
The properties of the input Gaussian process must be conveyed to the output Gaussian process by only the (uncorrelated) noisy observations {(xi, µ(xi), σ(xi))}Ni=1. As the uncertainty in original GP increases, so do the measurement uncertainties in the transmission, and therefore the output GP also has a higher uncertainty. However, the uncertainty in the input GP is in the form of a full covariance kernel k(x, x′) and it seems that individual observations will not easily be able to communicate the covariance of the values of the GP function at different spatial locations despite the heterogeneous noise model.
Fundamentally, the problem is that the observation values are treated as independent, an incorrect assumption which has other knock-on effects when the number of observations is large. With some fixed measurement error no matter how high but a large enough set of independent observations, the mean value can be pinned down precisely. If in contrast the observations are not independent, then there may be a situation where the mean value cannot be known more precisely than some limiting uncertainty. This effect leads the output GP to have less uncertainty and be more confident in the values that it should be given the input GP.
If the observations are sparse, then the effective sample size of the estimator for the mean of the GP at any given location is small, and then the amount by which uncertainty is underestimated is small. However, if there are very many observations then this kind of observation transmission of information with the independence assumption will attenuate the uncertainty. We would also expect that over the course of many layers, this attenuation can accumulate. We believe that this is what causes the poorer uncertainty calibration in layers 3 and 4 of the PNCNN shown in figure 2. We hope that this problem can be resolved perhaps by removing the independence assumption or providing an alternative projection method in future work.

heading: I CENTRAL LIMIT THEOREM FOR STOCHASTIC PROCESSES
text: In this section we show that each of the outputs features fα = ∑ βMαβhβ are marginally distributed as Gaussian processes in the limit as cin →∞, provided that hβ are independent and have bounded moments, and that the fraction of zeros in each row of M does approach 1. The result follows from an application of the multivariate Lyapunov CLT to the joint distribution of the functions evaluated at a finite collection of points.
Since we are only making the case about the marginal distribution of each of the output channels, we can separate them out individually and apply a CLT argument on each one. For any given output channel g, the output is a sum g = ∑ βmβhβ for values mβ given by the particular row of M . For
each coefficient mβ that is nonzero (with index β = i), we can define an element h̃i = mihi. The output is the sum of such elements g = ∑ i h̃i.
Now choose any finite collection of points x1, x2, . . . , xN . The random vector [g(x1), g(x2), ..., g(xN )] = ∑ i[h̃i(x1), h̃i(x2), ..., h̃i(xN )]. In the limit as the number of channels cin goes to infinity, we can apply the Lyapunov central limit theorem to this sum of random vectors, so long as the number of elements in the sum also goes to infinity. This means that the vector [g(x1), g(x2), ..., g(xN )] converges multivariate Gaussian with mean µ having entries µn = ∑ i E[h̃i(xn)] and covariance Σ with entries Σnm = ∑ i E[h̃i(xn)h̃i(xm)]. Since we can repeat the same argument for any collection of points x1, x2, . . . , xN , we must deduce that g is a Gaussian process according to the definition. The mean and covariance of ḡ(x)→ GP(µ, k) follow from the mean and covariance of these finite collections, µ(x) = ∑ i E[h̃i(x)] = ∑ imiE[hi(x)]
and k(x, x′) = ∑ i E[h̃i(x)h̃i(x′)] = ∑ im 2 iE[hi(x)hi(x′)].
Conceivably the CLT could be violated during training if a large enough fraction of Mαβ converged to 0 during training, or also if the separate channels became more strongly correlated. We investigate the non-Gaussianity of the features in section 4.7 and find that the features are mostly Gaussian but have somewhat longer tails indicating deviation from the theoretical behavior.

heading: J MOMENTS OF RECTIFIED GAUSSIAN RANDOM VARIABLES
text: Let f ∼ N (µ,Σ) be a d dimensional random Gaussian vector. We compute here
E(ReLU(f1) · · ·ReLU(fd)) = 1
NΣ ∫ f>0 ddf f1 · · · fd exp− 12 (f − µ) TΣ−1(f − µ) (33)
NΣ = (2π) d/2 det(Σ)1/2 . (34)
We use the generating function technique. Define
Z(b) = 1
NΣ ∫ f>0 ddf exp[− 12 (f − µ) TΣ−1(f − µ) + bTf ] (35)
= 1 NΣ eb Tµ ∫ f>−µ ddf exp[− 12f TΣ−1f + bTf ] . (36)
Then
E(ReLU(f1) · · ·ReLU(fd)) = ∂ ∂b1 · · · ∂ ∂bd Z(b) ∣∣∣∣ b=0 . (37)
To compute Z(b) we proceed as in Gaussian case. We change variables to f = Σb+ g , (38)
and define z = µ+ Σb to get:
Z(b) = eb Tµ+ 1 2b TΣb 1
NΣ ∫ g<+z ddg exp[− 12g TΣ−1g] (39)
= eS(b) Φ(d)(z; 0,Σ) , S(b) = bTµ+ 12b TΣb . (40)
Φ being the multivariate standard Normal CDF:
Φ(d)(z;µ,Σ) = ∫ g<+z ddgψ(d)(g;µ,Σ) , (41)
ψ(d)(g;µ,Σ) = 1
NΣ exp[− 12 (g − µ) TΣ−1(g − µ)] . (42)
Now we compute the first two derivatives. Note that in d = 1, denoting σ2 = Σ: ∂
∂z Φ(1)(z, 0, σ2) = ψ(1)(z, 0, σ2) . (43)
In d = 2, we can use the conditional probability decomposition to get the required derivatives:
ψ(2)(g; 0,Σ) = ψ(1)(g1;α1g2, β1) · ψ(1)(g2; 0,Σ22) (44) α1 = Σ12Σ −1 22 , β1 = Σ11 − Σ12Σ −1 22 Σ21 (45)
∂z2Φ (2)(z, 0,Σ) = ψ(1)(z2; 0,Σ22) ∫ z1 −∞ dg1ψ (1)(g1;α1z2, β1) , (46)
= ψ(1)(z2; 0,Σ22)Φ (1)(z1, α1z2, β1) (47)
∂2z2Φ (2)(z, 0,Σ) = − z2
Σ22 ψ(1)(z2; 0,Σ22)Φ (1)(z1, α1z2, β1) (48)
+ ψ(1)(z2; 0,Σ22)∂z2Φ (1)(z1, α1z2, β1) (49)
∂z1∂z2Φ (2)(z, 0,Σ) = ψ(2)(z; 0,Σ) . (50)
So denoting ∂i = ∂∂bi , we get: ∂iZ(b) = (µi + ∑ j Σijbj)Z(b) + e S(b) ∑ `
∂z`Φ (d)(z, 0,Σ)Σ`,i︸ ︷︷ ︸
mi(b)
(51)
∂k∂iZ(b) = ΣikZ(b) + (µi + ∑ j Σijbj)∂kZ(b) (52)
+ (µk + ∑ j Σkjbj)mi(b) + e S(b) ∑ `,q ∂z`∂zqΦ (d)(z, 0,Σ)Σ`,iΣq,k . (53)
In particular, for the first moment d = 1 we have:
E(ReLU(f)) = µΦ(1)(µ; 0, σ2) + ψ(1)(µ, 0, σ2)σ2 , (54)
which coincides with equation 8. Note that ψ(1)(µ, 0, σ2) = 1σψ(µ/σ; 0, 1) because of the normalization factor. For the second moments d = 2 we have:
E(ReLU(f1)ReLU(f2)) = Σ12Φ(2)(µ; 0,Σ) + µ1µ2Φ(2)(µ; 0,Σ) + µ1m2(0) + µ2m1(0) (55) + ∑
`,q=1,2
Σ`,1Σq,2∂z`∂zqΦ (2)(z, 0,Σ)|b=0 . (56)
which can be rewritten in the form of equation 9.

Review summaries:
This work presents an uncertainty aware continuous convolutional layers for learning from continuous signals like time series/images. This work is most useful in the setting of irregularly sampled data. Gaussian processes (GP) are used to represent the irregularly sampled input. The proposed continuous convolutional layers can be directly applied to input Gaussian Process in a closed form, which subsequently outputs another GP with transformed mean and variance. Finally, the continuous convolutional layers are parameterized in terms of the flow of a PDE. The use of GP for feature representation and the ability of continuous convolutional layers to take GP as input provides the model with the ability to propagate uncertainty.

The paper considers the problem of using CNNs on an irregularly sampled grid. It proposes to incorporate the numerical uncertainty related to the disretisation of the domain using an approach inspired by probabilistic numerics. This involves defining continuous convolutional layers as affine transformations of the input GP and use rectified GPs to include the nonlinearity in the mapping. The resulting non-Gaussian distribution is then approximated with another GP in each layer. These are then approximated with a GP with an RBF kernel using Monte Carlo to produce observations with input dependent noise. The model is trained using MAP estimation.