Probabilistic Numeric Convolutional Neural Networks
Authors: Marc Anton Finzi
Roberto Bondesan
Max Welling
Venue: ICLR 2021 Poster (2021)
Keywords: probabilistic numerics
gaussian processes
discretization error
pde
superpixel
irregularly spaced time series
misssing data
spatial uncertainty

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, 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\times$ 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.

PROBABILISTIC NUMERIC
CONVOLUTIONAL NEURAL NETWORKS
Mare Finzi*
Roberto Bondesan & Max Welling
Qualcomm AI Research
Qualcomm AI Research
New York University ma www.0ennyueeeee
(rbondesa, mwelling) @qt i. qua comm. com
ABSTRACT
Continuous input signals like images and time series that are irregularly sampled or have missing values are challenging for existing deep learning methods. Co- herently defined feature representations must depend the values in unobserved regions the input. Drawing from the work in probabilistic numerics, we propose Probabilistic Numerie Convolutional Neural Networks which represent features Gaussian processes (GPs), providing a probabilistic description of discretization error. We then define convolutional layer as the evolution of 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 that our approach yields a 3 reduction of error from the previous state of the art on the SuperPixel-MNNI dataset and competitive performance on the medical time series dataset PhysioNet2012.
show
INTRODUCTION
Standard convolutional neural networks are defined regular input grid. For continuous signals like time series and images, these elements correspond to regular samples of an underlying function defined on a continuous domain. In this case, the standard convolutional layer of neural network is a numerical approximation of a continuous convolution operator A.
Coherently defined networks continuous functions should only depend on the input function f. and not on spurious shortcut features (Geirhos et al. (020) such as the sampling locations or sam- pling 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 deter- mined by the sampling resolution. In some sense, this error is unavoidable because the features f(() the layers l 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 dispersed sensors for measuring currents, this discretization етгоr cannot be neglected. Simply filling in the missing data with zeros or imputing the values is not sufficient since many different imputations 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 un- certainty into the computation of the network. To do so, we replace the discretely evaluated feature maps (x;;), with Gaussian processes: distributions over the continuous function (() that track the most likely values as well as the uncertainty. On this Gaussian process feature represen- tation, 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 given feature is a Gaussian process, then applying linear operators yields new Gaussian process with transformed mean and covariance functions. The dependence of Af on regions of f which are not known translates into
"Wor done during internship at Qualcomm Al Research 1 "Qualcomm AI Research is 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).
2 RELATED WORK
Over the years there have been many successful convolutional approaches for ungridded data such as
GCN (Kipi 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 SShhitt al. 2017 Wang al.) 2018). However, the target domains of sets, graphs, and point clouds are intrinsically ly 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 regular grid (Li and Marlin Futoma et al., [017 Shukla and Marlin (20 19 Fortuin al., 22020). set functions (Horn
et al. 2019) and attention (Shukla and Marlin (020)). Additionally there are several works exploring
reconstruction of images from incomplete observations for downstream classification (Huijhen 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, Gaussian process is fit to an irregularly spaced time series and sampled imputations from this process are fed into separate RNN classifier. Unlike our approach where the classifier operates directly on continuous and probabilistic signal, in these works the classifier operates on a deterministic signal on regular grid and cannot reason prohabilistically about discretization errors.
Finally, while superficially similar to Deep GPs (Damianou and Lawrence, 2013) or Deep Differen- tial Gaussian Process Flows (Hegde al. 2018), our PNCNNs tackle fundamentally different kinds of problems like image classificatior and our GPs represent epistemic uncertainty over the values of the feature maps rather than the parameters of the network.
3 BACKGROUND
Probabilistic Numerics: We draw inspiration for our approach from the community of proba- bilistic numerics where the error in numerical algorithms are modeled probabilistically, and typically with 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 Gaussian process is used to model the error in the numerical estimation of an integral and optimally select rule for its computation (Minkal 2000 Rasmussen and Ghahramani 2003). Probabilistic numerics has been applied widely to numerical problems such as the inversion of matrix (Hennig. 2015), the solution of an ODE (Schober al. 2019). 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 fxx 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 {at J.V. 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 reviewed in appendix A These beliefs are encoded into a prior covariance i of the GP f ~ CP k) and updated upon seeing data with Bayesian inference. Explicitly, given set of sampling locations x = (xr N. and noisy observations y [yr sampled ~ (x.())) of), using Bayes rule
beliefs about values that have not been observed (Rasmussen al.. 2006), as
'hile GPs can be applied directly to image classification, they not well suited to this task even with convolutional structure baked in, shown in kumar al. (2018).
2

---

Input Data
GP
GP
GP
GP
o
MD
M
M
M
M
J f((x)xx
O
N 2p)
ы(н) o(H)
p(()
ы{x) o(x)
Probabilistic Rel.U
Probabilistic
Probabilistic ReLU
ReLU
Figure 1: The PNCNN operating on SuperPixel-MNIST images shown on the left. The mean and elementwise uncertainty of the Gaussian process feature maps are shown as they are transformed through the network by the convolution layers. Observation points shown as green dots (xx)
one can compute the posterior distribution fly. x ~ GP(pp. kp), which captures our epistemic uncertainty about the values between observations. The posterior mean and covariance given by
(x(xx) = k(x)) [K + 51 y; kip = k( ((x)) T + 5)='k((a'))
(1)
where K, k(,,, 22), k(z); k(x 24) and S diag(o??. Below we shall choose the RBF ker- neflas prior covariance, due its convenient analytical properties: KRBF (Ir, x²) = aN ( ((1)
(2x(²) exp(- / x 2²) In typical applications of GPs 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 way of representing our beliefs and epistemic uncertainty about the values of both the input function and the intermediate feature maps of neural network.
PROBABILISTIC NUMERIC CONVOLUTIONAL NEURAL NETWORKS
4.1 OVERVIEW
Given an input signal f X' Re, 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) = ʃ. as series of continuous convolutions ( with pointwise ReLU nonlinearities and weight matrices which mix only channels (known as x convolutions) M E Reexe
p(++ 2
(2)
and final global average pooling layer P which acts channel-wise as natural generalization of the discrete case: P(f(2)) = f f((" (z)dr for each a 1, .... Denoting the space of functions on X' with e channels by He, the convolution operators ( are linear operators from He, to Heet1 Like in ordinary convolutional neural networks, the layers build up increasingly more expressive spatial features and depend on the parameters in ( and M Unlike ordinary convolutional networks, these layers are well defined operations on the underlying continuous signal.
While it is clear that such network can be defined abstractly, the exact values of the function (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 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 ky) in closed form. The output is another Gaussian process with a transformed mean and covariance Af GP (App, Ak) where A' acts
For convenience, include the additional scale factor (2m12)d//2 relative to the usual definition.
3

---

to the left the primed argument of ky x²) In section 4.2 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.
Applying the ReLU nonlinearity to the GP yields a new non Gaussian stochastic process h (e) ReLU ((0)], and we show in section 4.5 that the mean and covariance of this process has closed form solution which can be computed.
Probabilistic ReLUs:
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(l+1) M(0 is approximately distributed as a Gaussian Process in section 4.5
Measurement and Projection RBF Gaussian Process: While f((+1) is approximately a Gaus- sian process, the mean and covariance functions have 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 Gaussian out- put. Assembling these components, we implement the end to end trainable Probabilistic Numeric Convolutional Neural Network which integrates a probabilistic description of missing data and dis- cretization error inherent to continuous signals. The layers of the network are shown in figure
4.2 CONTINUOUS CONVOLUTIONAL LAYERS
On a discrete domain such as the lattice X' Zdd 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 linear combination of powers of the generators of the translation group: the shift operators Ti i = .....,dd 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 Σn WAT where the weight matrices W4 e gexe act only on the channels and the shift operator T acts on functions on the lattice. In d dimensions, for some set of integer coefficients ky...., kg. For example when d 2, we can take ,,, h2 E (, 1} to fill out 3 x 3 neighborhood.
A = Ekn....ha WA Wk......e A
T1 ... -d ki
On the continuous domain X' = Rd we similarly parametrize convolutions with A= Σ. Wee where Dr is given by powers of the partial derivatives 0, i l...., d which generate intinitesimal translations along the i-th axes. Setting d = 1 for simplicity, we can indeed verify that the operator exponential 7" ead applied to function g(z) is translation: end= g(x) g(x) + ag' ag'(z) +
}a' + 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 e. =>x D /k! or more broadly as the solution to the PDE:
Oig(1, x) = (I)) g(0,x) = g(x)
at time t 1: e g(t 1,
Following the discussion in the discrete case, translation invariance of Dz imposes that it is ex- pressed in terms of powers of the generators. Collecting the derivatives into the gradient V, we can write the general form of Dr as a + 3 V + 2 ! Σ, V +.. for any constants a*₁, vectors Sk. matrices Σg etc. For simplicity, we truncate the series at second order get
D3 = 1 A E
(4)
where we omit the constants ay that can be absorbed into the definition of Wt. For this choice of D. the PDE in equation iss nothing but the diffusion equation with drift By and diffusion Σg. When discussing rotational equivariance in section 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 d.,,
MMore generally neural networks have affine layers including both convolutions and biases. An affine transformation Af + b of a Gaussian process is also Gaussian process f ~ GP(AAt + b. Akp ) and we include biases in our network hut omit them from the derivations for simplicity.
4

---

3dt EB where B, is a d dimensional Brownian motion, the time evolution operator can be written via the Feynman- Kac formula as f'D((²) = [X)(X,) where Xu = 2. In other words, the operator layer A= e¹ is the expectation under parametrized Neural Stochastic Differential equa-
tion (Li et al. 2020 Tzen and Raginsky 2019) that is homogeneous and therefore shift invariant.
The flow of this SDE depends the drift and diffusion parameters B and Σ.
To recap, we define our convolution operator through the general form = En Wg where the weight matrices Ws E IRex mix only channels and e is the forward evolution by one unit of time of the diffusion equation with drift 3x and diffusion E.. containing learnable parameters {(W,, Se, Ex) H=) The translation equivariance of A follows directly from the fact that the gen- erators commute Ak, : (Px: V.] and therefore [A, 7.] (the bracket [a, ab ba the commutator of the two operators). In appendix B we show that our definition of AA reduces to the usual one in the discrete case and is thus a principled generalization the continuous domain.
4.3 EXACT APPLICATION ON RBF GPs
Although the application of the linear operator A = = Wxe involves the time evolution of 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 we need only compute the action of AA on the mean and covariance: Af ~ (Atpp A') where A' is the adjoint w.r.t. the (2()) inner product. The application of time evolution e is convolution with Green's function Gk, so Af Σn Wa = Σn W&G2 * f. As we derive in appendix C the Green's 's function for Dr P V + (1/2) En V, nothing but the multivariate Gaussian density G& = NNx;; - Bk, Ex):
GP kp),
Af = Σ Whe Whe Puf= f= Σ W
Σ E.
(5)
In order to apply eD 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 [ are expressed in terms of KRBF = =N x', (1 ((())) and the computation boils down to convolution of two Gaussians:
knaF (x, x²) = N(;; -2- aN (n) = - -1, c KRBr KREr (f aW t(1
(6)
t22)
The application of the channel mixing matrices W; 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 's funetion above, the action of encompasses the ordinary convolution
operator the 2d lattice as special case. Given drift 8g E {-1, 1Jx2 k l,..., 9 filling out the 9 elements of a x 3 grid and as the diffusion Σς - 0, the Green's function a Dirac delta, so that: Az Σx Wh5(1 3x) f(z) Σxj=-1.0.1 Wiy f -i, i, j) W *z2 f(r).
4.4 GENERAL EQUIVARIANCE
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 z E X'. A linear layer A is cquivariant if its action commutes with that the symmetry. In appendix E we derive constraints for general linear operators and symmetries, which generalize those appearing in the steerable- CNN literature (Weiler and Cesa 2019 Cohen et al.) 2019). Then we show how equivariance under continuous roto-translations in 2d constrains the form of convolutional layer by solving the equivariance constraint. Non-trivial solutions require that the operator D in the PDE of equation 3 has non-trivial matrix structure.
5

---

4.5 PROBABILISTIC NONLINEARITIES AND RECTIFIED GAUSSIAN PROCESSES
Gast and Roth (0018) derive the mean and variance for a univariate rectified Gaussian distribution for use in neural network. We generalize these results to the full covariance function (and higher moments) of a rectified Gaussian process in appendix D and present the results here. For the input GP (())())) ~ (P(u(z), k(x, x')). we denote (xx
k(x,,) x), Σ the matrix with components Eij k(zi, Zg) for i.j 1.2 and μ = [u(a)), (x22) We use notation ((2) for the univariate standard normal CDF, and (2; Σ) for (two dimensional) multivariate CDF of N Σ) at z. Σ and Σg are the column vectors of E. The first and second moments of RcLU[Af] are:
B[h(=)] = (( (( (x())
(8)
B(h(m)) = + (M(2)) + ((2)) E (w Σ) + E VV (1;; Σ)Σ2,
(11; + ((()))())))))))))) 4 (11; Σ)
(9)
where VV o 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.
4.6 CHANNEL MIXING AND CENTRAL LIMIT THEOREM
After the non-linearity the process is no longer Gaussian. To overcome this issue we introduce a channel mixing matrix M( E TR xee and define the feature map in the following layer by f(e+1) = (6) where (0 = ReLUU f(O], So long as the channels of (0) are only weakly dependent, we can apply the central limit theorem (CLT) to each function f!((t1) Ei= M Magh so that in the limit of large ct, the statistics of the fa f(f+1) converge to GP with first and second moments given by:
(2)
[ (E+))) (x)) MET (c)); (x)
(x) (TM (10)
The convergence to Gaussian process here is reminiscent of the well known infinite width limit of Bayesian neural networks Neal (196);; 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 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]] [[ 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.
4.7 MEASUREMENT AND PROJECTION TO RBF GAUSSIAN PROCESS
As last step we simplify the mean and covariance functions of the approximate GP f((+1), While we can readily compute the values of these functiinn, unlike in the RBF kernel case, we can- not apply the convolution operator e!D in closed form. In order to circumvent this challenge, we model the (approximately) Gaussian process ((+1) with an RBF Gaussian process as follows: we evaluate the mean y: (ff()))(xz)) and variance o² = (arlf(+1)) (x:)) of the approximate Gaussian process f((+1) at a collection of points {a,).M using equations g 9 and 10 These values Yi are treated as measurements of the underlying signal with a heteroscedastic noise of that varies from point to point. We can then compute the RBF-baser posterior GP of this signal f(++)) yr, ) GP(Mp, kp) with posterior mean and covariance given by equation | for the heteroschedasttcc noise model. The uncertainty in the input f(()) is propagated through to the RBF posterior f((+1) [[r, yi. x via the measurement noise a, Crucially, this Gaussian pro- cess 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 kemel in each layer is trained to maximize the marginal likelihood of the data that sees, and thereby minimize the discrepancy with the under- lying generating distribution (E+1). While this measurement/projjetioon approach is effective in

---

Uncertainty Prediction
101
layer
(uDx) -gt]aaxx) Quantiles
5.0
layer
256
Resolution (Pixels)
+
2.5
100
+
0.0
(x) u (P))))
OCI/VR)
-2.5
(x)())
10
10a
10a
Train Res N-75
-5.0
Test Pixels 14 29
58
110
167
239
Test Resolution (pixels)
N(O,1) Quantiles (offset)
Figure 2: Left: Qualitative convergence of the mean and uncertainty of the first 3 channels of the feature maps is shown in RGB color as the input test resolution is increased. Middle: Median predicted uncertainties over spatial locations as function of the test resolution. Right: Using the predictions of the highest resolution model ground truth, the distribution of prediction residuals is shown in Q-Q plot for each layer (shifted horizontally for clarity) with the black lines showing the theoretical relationship, and the overall distribution histogram shown on the right.
many scenarios, in networks with many layers 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 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 (middle) we show that in early layers the uncertainties decrease at similar rate to the (()/NN) of discretization error that we would expect from a standard convolutional layer which is discretized to square grid 4l Despite the fact that these resolutions differ substantially from those seen at training time and the fact that there are 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 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 range of resolutions.
4.8 TRAINING PROCEDURE
Our neural network has two sets of parameters: the channel mixing and diffusion parame- ters, {(M()) W( ,(0),, H as well as kernel hyperparameters of the Gaussian Processes {((() ( We train all parameters jointly on the loss Lrask + ALGp, where Lrask is the cross entropy with logits given by the mean pP of the pooled features P(f(t)) ~ (NP,,, Ep) and Lgp the marginal log likelihoods of the GP feature maps:
L
Lop(f)
1
ΣΣ t=1
1Ixx Sal fa) + log det [xxx + So] + N log 2n
(1I)
where for each layer E, fa (f. [aN)]] E RN are the observed values for channel at locations X = ..,,..., Kxx is the covariance of the RBF kernel and S. = diag(oā) the measurement noise for each channel a and spatial location. Notably the GP marginal likelihood is independent of the class labels.
5 EXPERIMENTAL RESULTS
We evaluate the Probabilistic Numeric CNN on two different problems which have incomplete and irregular observations.
2D discrete convolution layer using N m²² points can be interpreted as Riemann sum approximation of the continuous integral and will therefore have error of O(1/m) O(1//N) which is the same rate as would be achieved through Monte Carlo sampling.
7

---

CNN
PNCNN without
PNCNN with
Train Resolution
16
Test Error
36
0.5
64
100
144
0.0
196
256
0
100
200
0
100
200
0
100
200
Test Resolution (pixels)
Figure 3: Zero shot generalization to other resolutions: Having trained on MNIST at given training resolution shown by the color, we evaluate the performance an PNCNN, PNCNN without uncer- lainty, and an ordinary CNN on varying test resolutions. Notably, PNCNN with uncertainty is the most robust.
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 the underlying images, the lack of grid structure and high fraction of missing values make this challenging task. Example inputs are visualized at the left figure ! We compare to Monet (Monti et al., 2017), SplineCNN (Fey et al, 20188) Graph Convolutional Gaussian Processes (GCGP) (Walker and Glocker 2019), and Graph Attention Networks (GAT) (Avelar et al., (020). As shown table the probabilistic numeric CNN greatly outperforms the competing methods reducing the classification eπrоr rate by more than 3x over the previous state of the art.
MNIST Superpixel 75
Error Rate (L)
PhysIONet2012
AP (t)
AUROC (1)
Monet
8.89
SplineCNN
4.78
IP-Nets
51.86
86.24
GCGP
4.2
SEFT-ATTN
53.67
85.14
GAT
3.81
GRU-D
54.97
86.99
PNCNN
0.24++0.12
PNCNN
53.92+t.17
86.13.+
PNCNN w/o a
3.03
Table Left: Classification Error on the 75-SuperPixel MNIST problem. Right: Average Precision and Area Under ROC curve metrics for PhysioNET2012. Mean and standard deviation are computed over trials.
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 O. This form of the network still makes use of the fact that the input function is continuous by 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 σr) with 3.03% error rate outperforms existing methods from the literature. it is substantially worse than the PNCNN that integrates the uncertainty at 1.24% eтror. 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 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 regular square grid: PNCNN with uncertainty 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 model on the irregularly spaced time series dataset PhysioNet2012 (Silva et al., 2012) for predicting mortality from ICU vi- tals 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 Jefficiently for all channels simultaneously because the observation
8

---

points } and hence the matrices K in that formula differ between the channels, increasing com- putational complexity. To circumvent this difficulty, we employ stochastic diagonal estimator to compute the variances described in appendix G We compare against IP-Nets (Shukla and Mar- lin 2019), SEFT-ATTN (Horn et al. 2019), and GRU-D (Che et al.) 2018 as reported in Horn al. (2019). PNCNN performs competitively, although not breakout performance in the image dataset which we attribute to the use of the stochastic variance estimates over an exact calculation.
6 CONCLUSION
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 CNN. The layers of our network are a series of operators defined on continuous functions, removing de- pendence 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 gen- erative modeling of the internal feature maps of 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.

---

REFERENCES
Matan Atzmon, Haggai Maron, and Yaron Lipman. Point convolutional neural networks by exten- sion operators. arXiv preprint arXiv: 10091, 2018.
Pedro HC Avelar, Anderson R Tavares, Thiago LT da Silveira, Cláudio R Jung, and Luís C Lamb. Superpixel image classification with graph attention networks. arXiv preprint 0rX55:2002.0554,, 2020.
Costas Bekas, Effrosyni Kokiopoulou, and Yousef Saad. An estimator for the diagonal of matrix. Applied numerical mathematics, 5701112):1214-2225 2007.
Erik Bekkers. B-spline cnns on lie groups. arXiv preprint arXiv:1909, 12057, 2019.
Zhengping Che, Sanjay Purushotham, Kyunghyun Cho, David Sontag, and Yan Liu. Recurrent neural networks for multivariate time series with missing values. Scientific reports, 8(1):1-12, 2018.
Jon Cockayne, Chris Oates, Tim Sullivan, and Mark Girolami. Probabilistic numerical methods for partial differential equations and bayesian inverse problems. arXiv preprint arXiv: 1605.07811, 2016.
Jon Cockayne, Chris J Oates, Timothy John Sullivan, and Mark Girolami. Bayesian probabilistic numerical methods. SIAM Review, 61(4):756-789 2019.
Taco Cohen and Max Welling. Group equivariant convolutional networks. In International confer- on machine learning, pages 2990-2999, 2016a.
Taco S Cohen and Max Welling. Steerable enns. arXiv preprint arXiv:6612.08495 2016b.
Taco S Cohen, Mario Geiger, and Maurice Weiler. A general theory of equivariant cnns on ho- mogeneous spaces. In Advances in Neural Information Processing Systems, pages 9145-9156, 2019.
Andreas Damianou and Neil Lawrence. Deep gaussian processes. In Artificial Intelligence and Statistics, pages 207-215, 2013.
Alexander G. de Matthews, Mark Rowland, Jiri Hron, Richard E. Turner, and Zoubin Ghahra- mani. Gaussian process behaviour in wide deep neural networks, 2018.
Pim de Haan, Maurice Weiler, Taco Cohen, and Max Welling. Gauge equivariant mesh cnns: Anisotropic convolutions on geometric graphs. arXiv preprint arXiv:2003.05425, 2020.
Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019.
Matthias Fey, Jan Eric Lenssen, Frank Weichert, and Heinrich Müller. Splinecnn: Fast geomet- ric deep learning with continuous b-spline kernels. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 869-877. 2018.
Mare Finzi, Samuel Stanton, Pavel Izmailov, and Andrew Gordon Wilson. Generalizing convolu- tional neural networks for equivariance to lie groups arbitrary continuous data. arXiv preprint arXiv:2002.12889. 2020.
Vincent Fortuin, Dmitry Baranchuk, Gunnar Rätsch, and Stephan Mandt. Gp-vae: tic time series imputation. In International Conference on Artificial Intelligence and Statistics, pages 1651-1661, PMLR, 2020.
Deep probabilis-
Fabian B Fuchs, Daniel E Worrall, Volker Fischer, and Max Welling, Se (3)ntffsforrm: 3d roto- translation equivariant attention networks. arXiv preprint arXiv:2006. 10503, 2020.
Joseph Futoma, Sanjay Hariharan, and Katherine Heller. Learning to detect sepsis with a multitask gaussian process mnn classifier. arXiv preprint arXiv: 1706.04152, 2017.
10

---

Jochen Gast and Stefan Roth. Lightweight probabilistic deep networks. In Proceedings the IEEE Conference on Computer Vision and Pattern Recognition, pages 3369-3378, 2018.
Robert Geirhos, Jörn- Henrikk Jacobsen, Claudio Michaelis, Richard Zemel, Wieland Brendel, Matthias Bethge, and Felix A Wichmann. Shortcut learning in deep neural networks. arXiv preprinl arXiv:2004.07788 2020.
Pashupati Hegde, Markus Heinonen, Harri Lähdesmäki, and Samuel Kaski. Deep learning with differential gaussian process flows. arXiv preprint ftth::881004066. 2018.
Philipp Hennig. Probabilistic interpretation of linear solvers. SIAM Journal on Optimization. 25(1): 234-260, 2015.
Max Horn, Michael Moor, Christian Bock, Bastian Rieck, and Karsten Borgwardt. Set functions for time series. arXiv preprint arXiv: 1909,1004,, 2019.
Iris AM Huijben, Bastiaan S Veeling, and Ruud JG van Sloun. Deep probabilistic subsampling for task aask-adaptive compressed sensing. In International Conference on Learning Representations, 2019.
Diederik Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv: (412.6422..) 2014.
Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional net- works. arXiv preprint arXiv: 2016.
Vinayak Kumar, Vaibhav Singh, PK Srijith, and Andreas Damianou. Deep gaussian processes with 2018.
convolutional kernels. arXiv preprint /
Steven Cheng-Xian Xian Li and Benjamin M Marlin. A scalable end-to-end gaussian process adapter for irregularly sampled time series classification. In Advances in neural information processing systems, pages
2016.
Steven Cheng-Xiaa Xian Li and Benjamin M Marlin. Learning from irregularly-sampled time series: A missing data perspective. arXiv preprint arXiv:2008.075: 07599, 2020.
Xuechen Li. Ting-Kan Kam Leonard Wong, Ricky TQ Chen, and David Duvenaud. Scalable gradients for stochastic differential equations. arXiv preprint arXiv:2001.014 3288 2020.
Thomas P Minka. Deriving quadrature rules from gaussian processes. Technical report, Technical report, Statistics Department, Carnegie Mellon University, 2000.
Federico Monti, Davide Boscaini, Jonathan Masci, Emanuele Rodola, Jan Svoboda, and Michael M Bronstein. Geometric deep learning on graphs and manifolds using mixture model enns. In Proceedings of the IEEE Conference Computer Vision and Pattern Recognition, pages 5115- 5124, 2017.
Radford M Neal. Priors for infinite networks. In Bayesian Learning for Neural Networks, pages 29-53. Springer, 1996.
Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas. Pointnet: Deep learning point sets for 3d classification and segmentation. In Proceedings the IEEE conference on computer vision and pattern recognition, pages 652-660, 2017.
Carl Edward Rasmussen and Zoubin Ghahramani. Bayesian monte carlo. Advances neural information processing systems, pages 505-512,
C.E. Rasmussen, C.K.I. Williams, M.I.T. Press, F. Bach, and ProQuest (Firm). Gaussian Processes for Machine Learning. Adaptive computation and machine learning. MIT Press, 2006. ISBN 9780262182539, URL hatL:ttts://ooo
Yulia Rubanova, Ricky TQ Chen, and David Duvenaud. Latent odes for ireeularly-samplled time series. arXiv preprint arXiv: arXiv: (907:039077077 2019.

---

Michael Schober, Simo Särkkä, and Philipp Hennig. A probabilistic model for the numerical solu- tion of initial value problems. Statistics and Computing, 29(1):99-122, 2019.
Kristof Schütt, Pieter-Jan Kindermans, Huziel Enoc Sauceda Felix, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller. Schnet: continuous-filter convolutional neural network for modeling quantum interactions. In Advances in neural information processing systems, pages 991-1001, 2017.
Zhengyang Shen, Lingshen He, Zhouchen Lin, and Jinwen Ma. Pdo-econvs: Partial differential operator based equivariant convolutions. arXiv preprint arXiv:2007. 2020.
Satya Narayan Shukla and Benjamin Marlin. Multi-time attention networks for irregularly sampled time series. 2020.
Satya Narayan Shukla and Benjamin M Marlin. Incerpolation-precictioon networks for irregularly sampled time series. arXiv preprint arXiv: 190.00790 2019.
Ikaro Silva, George Moody, Daniel Scott, Leo Celi, and Roger G Mark. Predicting in-hospital mortality of icu patients: The physionet/compuuing in cardiology challenge 2012. In 2012 Com- puting in Cardiology, pages 245 248. IEEE, 2012.
Bart Smets, Jim Portegies, Erik Bekkers, and Remco Duits. Pde-based group equivariant convolu- tional neural networks. arXiv preprint arXiv:2001.09046 2020.
Michael L Stein. Interpolation of spatial data: some theory for kriging. Springer Science & Business Media, 2012.
Hugues Thomas, Charles R Qi, Jean- -eaannmma Deschaud, Beatriz Marcotegui, François Goulette, and Leonidas J Guibas. Kpconv: Flexible and deformable convolution for point clouds. In the IEEE International Conference on Computer Vision, pages 6411 2019.
Proceedings the
Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-and translation-qquivarian neural networks for 3d point clouds. arXiv preprint vrXiv:1802.082/99 2018.
Belinda Tzen and Maxim Raginsky. Neural stochastic differential equations: Deep latent gaussian models in the diffusion limit. arXiv preprint 1905. 008833, 2019.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Fukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in neural information processing systems, pages 5998-6008, 2017.
Ian Walker and Ben Glocker. arXiv:1905.05739, 2019.
Graph convolutional
gaussian processes. arXiv preprint
Shenlong Wang, Simon Suo, Wei-Chiu Ma, Andrei Pokrovsky, and Raquel Urtasun. Deep para- metric continuous convolutional neural networks. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 2589-2599 2597, 2018.
Maurice Weiler and Gabriele Cesa. General e (2)-equivrr steerable cnns. In Advances Neural Information Processing Systems, pages 14345, 2019.
Maurice Weiler, Fred Hamprecht, and Martin Storath. Learning steerable filters for rotation equiv- ariant cnns. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 849-858, 2018
Daniel E Worrall, Stephan J Garbin, Daniyar Turmukhambetov, and Gabriel J Brostow. Harmonic networks: Deep translation and rotation equivariance. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5028-5037, 2017.
Wenxuan Wu, Zhongang Qi, and Li Fuxin. Pointcomv: Deep convolutional networks on 3d point clouds. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 9621-9630 2019.
12

---

Greg Yang. Wide feedforward recurrent neural networks of any architecture gaussian pro- cesses. H. Wallach, H. Larochelle, A. Beygelzimer, F. dAlché-Buc, E. Fox, and R. Gar- nett, editors, Advances in Neural Information Processing Systems. volume 32. Curran Asso- ciates, Inc., 2019. URL https::tpp..........t.. (76aaaaaaaaa22200116777766666130
neurips httpa//scc./aaeee/////////// pdf
Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola. Deep sets. In Advances in neural information processing systems, pages 3391-3401,
Yanzhao Zhou, Qixiang Ye, Qiang Qiu, and Jianbin Jiao. Oriented response networks. In Pro- ceedings of the IEEE Conference Computer Vision and Pattern Recognition, pages 519-528, 2017.
13

---

A REVIEW OF GAUSSIAN PROCESSES
We briefly review here the main ideas of Gaussian Processes for machine learning, see Stein (012)); Rasmussen et al. (2006) for more details. We start 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 æ = ..I,...,,N) y = (y)..... 1/N); y1 = (xx), we update our belief's about the function by constructing the posterior via Bayes rule: pf y. x) = p(yl,, a))([))/((U
Here we need to specify the likelihood of the data with our model p(y) f. = II;_ p(silf, and the denominator, called the evidence or marginal likelihood, follows: p(y (z) = s-((()) (pty) (p(yl...f...... x))) The power of this approach 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 2r y, x) = (-pi (ppl:, 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. stochastic process can be presented in terms of the finite distributions of the random variables {f [+ points [ 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 ju(x;) and covariance kernel k(x, x²). We shall write ~ GP k). Let assume Gaussian likelihood model as well, i.e. f( xx) = (x;) o??) where on represents aleatoric uncertainty on the measurement. (For simplicity take here the function to be scalar valued but the reasoning can be easily generalized.) Then properties of the Gaussian distribution (see (Rasmussen al. 2006 Chap. 2) for detailed derivation of the formulas) lead the following posterior distributions after seeing data x: p(flp x) = (P kp), with
(x(xx) = k<t) + 5) 'y
y: hp( x') k(x, x²) k(x)"(K + -)-1)(a').
(12)
where Kig = k( k(z); = k(æ,, æ.) and S = diag(o?).
B FROM DISCRETE TO CONTINUOUS CONVOLUTIONAL LAYERS
We here show that the general formula
A Dr
(13)
with D& a function of spatial derivatives, reduces in the case of discrete input space X' the usual convolution we encounter in deep learning.
For simplicity we shall assume 1d grid as input space X' {1, .I,...,,A J. Let us start by recalling the form of the classical discrete convolution when C₁ Cc+1 = 1. We define convolutional layer as linear map that commutes with the translation operator. To make the symmetry exact, we need assume periodic boundaries. Then in the standard basis of R.V. (e:).z of vectors localized at site i, the translation operator т acts as TE, = 6s+1 mod N. An N N matrix B is translation invariant iff rB = Br. Since T is diagonal Fourier space, the most general solution is B = diag where Fjx- = Ag is the discrete Fourier transform. Such matrices are called circulant and can be written alternatively as B
DN'1 bN ; b F-16. Explicitly:
60
0i
bN- bN-
bz-
bo
b
bN
B
bx
bo
(14)
b2
b1
bi
b2
by
bo
This shows that the most general convolutional layer circulant matrix. F.g. if b. = unless i == 0, N, - 1, B coincides with the matrix representing periodic convolution of filter size 3. The matrix B is invertible long as b # for all k. In a convolutional network the parameters by are random variables and the measure of the set where is not invertible is zero. Thus the role of e is replaced in the discrete case by the B.
14

---

The discrete analog of A is then:
A= Σ W.BB,
(15)
Introducing the unit matrices E which have at the row a and column 3 and otherwise, we rewrite as:
A
Σ Ea,s & r' W Σ W.a.bbNN) j,008
(16)
Since E & rd is a linear basis of the space of convolutional layers, we see that equation 13 indeed reduces the usual one when discretizing the input domain and is principled generalization the continuous domain.
C GREENS FUNCTION
Given the operator D = 3 V + EV we can compute the action of e in terms of con- volutions. Using the d dimensional Fourier transforms [[h](k) (2x) h(z)e=ik dd and F-1 Ft, we can rewrite the derivative operator 2D in terms of elementwise multiplication in the Fourier domain, which diagonalizes D. Since V F (ik)F.
D F (iB k 1 EE)F.
(17)
Using the series definition eD
Em=o(t)) (tD) /n!, we have:
1 = |(3 1 Tx)
Ib)
(18)
Applying this operator to test function h(x)) yields
eD =h= [|(()
F[h]] = G2 ht
(19)
where the final step follows from the Fourier convolution theorem, and we define the function G. F (G 2- T Directly applying the Fourier integral yields Gaussian integral
G,(z) = (2x) eth (0+18)- N (Ek dk (*+++)) (1 (r++ det(2 (20)
This function (x(x) Nx: -13:-1,,,13 15) is nothing but multivariate heat kernel, the Greens function (also known as the fundamental solution or time propagator) for the diffusion equation 0,G₂ x') DG,(I - x'), and indeed lim,- Gz(z z²) = ō(x -')..
D INTEGRAL POOLING
The integral pooling operator Pif] = Sgo fx))dx can be applied to the Gaussian process just like
any other linear operator. Given ((L) (), k), we have that
P GP(P,, PkP' = N(P,, PkP').
(21)
Again, computing the mean pp Py and covariance matrix E, PkP' we need just to able to apply P to the RBF kernel.
J KRBE (z, x' x')xx =a
PKRBF(I') =
(22)
-f knaF (x,, x' )dada'
PkuEr
(23)
For many applications such image classification using the mean logit value, require only the predictive mean, so an unbounded covariance matrix Ep is acceptable. We use this form for all of experiments.
15

---

However for some applications output uncertainty can be useful, also provide variant that integrates over a finite region [0, 1| Pf = [oo11- fz))dz.
PkRBr((')
kRB x' dd II + (4) - $
(24)
PKRBFP's [0.11"x(0,1'"
knnp (x, a' dd alev 2 (e-1/22 - -1)+20(1)) -1]"
(25)
where 4 is again the univariate standard normal CDF.
E EQUIVARIANCE
E. 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 (20 (2016a) 6aa for discrete groups on lattices. This work has been extended for continuous groups (Worrall 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 al., [2020] [Fuchs al., [2200] de Haan et al.) (2020). In [Shen 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 SSett et al. 22200) In this work, the authors define layers of their convolutional network through the time evolution of a PDE which is 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 applied to scalar functions the group, unlike work where feature fields (of various representations) live on the base space R". These two approaches reflect the major approaches 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 seen in the same way because convolution translation equivariant. Our learnable linear operators A are equivariant to continuous transforma- tions. Two linear operators el and cB commute [eC, c = if and only if their generators commute: [C, = O. Since the generator of diffusions D; is sum of derivative operators, and the generators of translations are just V as mentioned section 4.2 the two commute: [Ds, VJ = 0. Therefore
[A, Ta]
E WeePs, Ta] = Στ We(ePi ea
0 and A is translation equivariant.
E.3 STEERABLE EQUIVARIANCE FOR LINBAR 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 such as the group of continuous rotations SO(d) ) in Rd Applying equivariance constraints separately on each of the components of ePi on top of translation equivariance yields very restricted set of operators. For example, enforeing equivariance to continuous rotations G SO(d), the operator must be an isotropic heat kernel: Dj = c V V. The reason for this apparent restriction is a result of considering the different channels independently, scalar fields.
the Linear operator to other symmetry groups G
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
16

---

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 IS (i.c. the channels) is formalized by representation matrix p(g) E Rexe for each element g E G in the transformation group that satisfies Vg, h E G p(gh) (()). Choosing the type of each intermediate feature map is equivalent choosing their representations, and we describe a simple way of doing this with tensor representations in the later section.
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 operator A (IR" + R + (IR - R ) 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: Apin (g) Pout (g) Ay where Lg (xr) f(g Rearranging the terms, one sees that the equivariance constraint on the linear operator is:
Operator Equivariance Constraint:
Pout(g)
= A,
(26)
where the operators Lg and L,2 are understood not to act on the representation matrices p (although
implicitly a function of g). As shown in Appendix E.4
eq.
26 is direct generalization of the Pin ( 1 = K (xx) described in the
equivariance constraint for convolutions Vor Poot g
literature (Weiler and Cesa
2019; Cohen et al. 2019).
As shown Appendix E.6 the equivariance constraint for continuous rotations applied to the dif- fusion operators A En WyeP* 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:
0f = --EWn WhDaff
(27)
Øl
the time evolution contains the matrices Wt in the exponential A = eE. WeDs, Like with the example of translation above, this operator is equivariant and only if the infinitesmal generator >n WaDx is equivariant.
Because equation 26 applies generally to linear operators and not just convolutions, we can com- pute the equivariance constraint for these derivative operators. We can simplify the summation Wx = En (
V V) by writing it in terms of the collection of matrices B. = E. W WRB and Sij (1/2) Ex WaEkij to express Aderiv where the indices j = 1,2.., d enumerate the spatial dimensions of each vector 3r and each ma- trix Σx. As we derive in appendix EE.5, the necessary and sufficient conditions for the equivariance of Σx W&Dx and therefore A is that Yg E G [Pout & Pin (g), vec(S') where vec(-) denotes flattening the elements into a single vector and P(r.s) is the tensor representation with r covariant and s contravariant indices.
B.
E, B,ü; + Σis S,y0,0,
vec( () and
Vgg E G [Pout & phn & P(2.,) (g)ve(())
E.4 GENERALIZATION OF EQUIVARIANCE CONSTRAINT FOR CONVOLUTIONS
This equivariance constraint is direct generalization of the equivariance constraint for convolution kernels as described in Weiler and Cesa (2019); Cohen et al. (019)) In fact, when A is a convolution operator, Af K f, the action of Lg by conjugation A is equivalent to transforming the argument of the kernel K:
/ 1I -2))x==))()))'d
Lo)
= (x- )d((=") (Lg
Letting both sides of eq 26 act the product of a constant unit vector e; and delta function, = =e, the expression Ve, (out(g) [ (in((ggg)))))) * e, = K * e,ō can be rewritten as Vz :
= (xx) which is precisely the constraint for steerable equivariance for convolution described in the literature.
Porat (g) K K(g-tr)pmn(
TThis assumes as is typically done that measure /a over which the convolution is performed is left invariant. For the more general case, see the discussion in Bekkers (209).
17

---

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 é an element of the vector space T(p..) V@p (V*)&9 where V is some underlying vector space, V* its dual and (-)or is the tensor product iterated times. In common language T(0,0) are scalars, T(1.0) are vectors, and T(,1) are matrices. Given the action of group G on the vector space V, the representation on (p-9) is P(p.g) = gop & Tg= where -T 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% @ (1.0) (1..)) @ T(, where @ in this context is the same as the Cartesian product. For a composite representation = Dj Tip.g.)) the group representation is similarly (v D P(p..a)) (g) where ① concatenates matrices blocks on the diagonal.
Noting that the operator Tng that acts only the argument and the matrix Pin (g) acts only the components, the two commute and we can rewrite the constraint for Aderiv as
Σ Pout Bipin ( Logot gg + Σ Port (ow((g) SipPin )) L.000, Aderiv
(28)
We can simplify the expression L, by seeing how it acts function. For any differentiable function f()) ăz, [f(g))] = E, flll L E, (x) where 9ij the components of the matrix g. Since this holds for any f. we find that Lg V = g" and therefore VVV...
Lg VVVVL = Lg L =
Since equation 28 holds as operator equation, it must be true separately for each component a and 0,0,. This means that the constraint separates into a constraint for B and a constraint for St
Vg, i: Pis Bg ()
Bi
2. Vg,ú,j Σn Puri (9) SekPinng Sij-
These relationships can be expressed more succinctly by flattening the elements of and S into vectors: (g) @ pin(g-1) (e) vec(B) and (ooa( & Pin (g "") @
(9))ee))) vec(S).
E.6 ROTATION EQUIVARIANCE CONSTRAINT FOR SCALAR DIFFUSIONS HAS ONLY TRIVIAL SOLUTIONS
The diffusion operator A Σn We leads to only trivial 8x 0 and Eg 7 if it satisfies the continuous rotation equivariance constraint.
Proof:
The application of eD is jsst a convolution with the Greens function
Σ Whee Σ Wale ((=-A)TT
(x+(++ "x, det( ( 1/2 1/21 f= Σ W +
(29)
where the Greens function is the multivariate Gaussian density: Gx(x) = - En).
As shown in appendix E.4 for convolutions the operator constraint is equivalent the kernel equivariance constraint (o) ((gg()))eoo
K(x) from (Weiler and Cesa 2019). With
(x)) E. (x) this reads:
Vx E Rd EEG: G: Σ Pout Powelo)W,N(1 x; -Bx E (mn = Σ WaN((: -k, Ex),
k
For rotations E SO(2) where we can parametrize go = 0o in terms of the antisymmetric matrix
J = [00 1), 1, 0] x2x2 and the translation operator be written I.g e 0xx -0=="_"
y²
18

---

take derivatives with respect to & to get (now with double sums implicit):
Vo E IRd Σ [doma Ww: -x;- Ww - -xx, Ee )dpi - Wu( J (1; Ex) - 0.
Here the Lie Algebra representation of is dp == d0 p I0=0 Factoring out the normal density:
Va E ER
E [dpons W Wadpin - Wa(a" (x (o) Be, Ex) =
Without loss of generality we may assume that each of the Gaussians 3k, Ex pairs are distinct since if they were not then we could replace the collection with single element. Since the (finite) sum of distinct Gaussian densities never a Gaussian density, and monomials of order > multiplied by Gaussian density cannot be formed with sums of Gaussian densities sums multiplied by monomials of a different order and Gaussian densities are never 0, this constraint separates out into several independent constraints.
1. Vi : dpout WA= = Wadpin
2. vi W& (1= 3x) =
3. vi Welx J' x)==
We may assume w.l.o.g. that W₂ is not O 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 O. Given that this holds for all eq 3 implies: J' 0 or equivalently 2xp = because MB is symmetric, and since J = -J" this can be expressed concisely as [E J = 0 for which the only symmetric solution is proportional the identity Σg C& J. Since both Eg and J are invertible, equation 2 yields 3g 0. Therefore there are no nontrivial solutions for 3, Σ A = Σn Wge for continuous rotation equivariance.
F DATASET AND TRAINING DETAILS
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 chain of convolutional blocks contain- ing convolution layer, probabilistic ReLUs, and linear channel mixing layer (analogue of the colloquial 1 x 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: {xx; juz;), (x;)) ) These observations seed the GP layer, and the block is eval- uated the same collection of points for the output (although can be evaluated elsewhere since
a continuous process, and we make use this fact to visualize the features in figures I [
Hyperparameters: For the PNCNN on the Superpixel MNIST dataset, we use 4 PNCNN convo- lution blocks with 128 channels and with K 9 basis elements for the different drift and diffusion parameters in E WxePs, We train for 20 epochs using the Adam optimizer KKinmma and Ba 2014) with Ir 310-3 with batch size 50.
For the PNCNN on the PhysioNet2012 dataset, we use the variant of the PNCNN convolution layer 20 probes. In the
that uses the stochastic diagonal estimator described in appendix G with convolution blocks we use 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 validation set of size 10% before folding the validation set back into the training set for the final runs Both models take about hours to train.
SuperPixel-MNIST We source the SuperPixel MNIST dataset (Monti et al. 2017 from Fey and Lenssen (2019) consisting of 60% training examples and 10k test represented as collections of posi- tions and grayscale values {x,, f (a.))) ?o the N 75 super pixel centroids.
PhysioNet2012 We follow the data preprocessing from Horm 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
19

---

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 these static signals as fully observed constant time series signals. As the binary classification task exhibits strong label imbalance, 14% positive signals, we apply an inverse frequency weighting of 1/.14 to the binary cross entropy loss.
G STOCHASTIC DIAGONAL ESTIMATION FOR PHYSIONET2012
In order to compute the mean and variance of the rectified Gaussian process, the activations of the probabilistic ReLU, need compute the diagonal of (In Ix) for the relevant points {
In the usual case where each of the channels a 1,2,..., e are observed at the same locations this can be done efficiently. First one computes the application c:: on the left and eP's on the right onto the posterior kip
N₂.- = (e kpp (i (Zn, = (e (e ke (2n, In) (ePxkT)(ww)IA + 5T- (In)
where k is the RBF kernel and we have reused the notation from appendix A Notably, this quantity the same for each of the channels, and the elementwise variance is just:
Va (IT) (AAlipA') (2n, Fn.)
Σ Wi Wop 6.3.8
(30)
where the at 33 index the channels of each of the matrices W. 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 { for each channel 3. This means that evaluated kernel depends on the channel and we have the objects: k2, KP and S As result, we have an additional index for N and the desired computation is
Va (2 (AAiiaA)) Jan (xn mn)
w Ew. N2
(31)
$J..3
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 (cn)( (asp ), and while it would not be feasible to compute this matrix directly we can define matrix vector multiplies onto vectors of size ReN 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 These implicit matrix vector multiplies can then be used to compute stochastic diagonal estimator (Bekas et al. 2007) given by:
(2
1
Σ 2p B
(32)
with Gaussian probe vectors Zp NO, 7), and where @ is elementwise multiplication (see Bekas ct 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 better approach for handling this kind of data to future work.
H PATHOLOGIES IN PROJECTION TO RBF GAUSSIAN PROCESS
In section 4.7
describe an approach by which a Gaussian process with complex mean and co- variance 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.
20

---

The properties of the input Gaussian process must be conveyed to the output Gaussian process by only the (uncorrelated) noisy observations {(z,, (lzz), a (x 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 a full covariance kernel k(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 eтror no matter how high but large enough set of independent observations, the mean value can he pinned down precisely. 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 We hope that this problem can be resolved perhaps by removing the independence assumption or providing an alternative projection method in future work.
I CENTRAL LIMIT THEOREM FOR STOCHASTIC PROCESSES
In this section we show that each of the outputs features fa = E Mashg are marginally distributed as Gaussian processes in the limit as can. + oo, 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 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 a sum g E. mahg for values ma given by the particular row of M. For each coefficient my that is nonzero (with index 3 i), we can define element h, m,hy. The output the sum of such elements g = E,
Now choose any finite collection of points T1, Ø2. EN. The random vector [ø(x1), (1,2),..,
(z)h( hı(22), (=N)) In the limit as the number of channels Gin 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 (((2))) g(.,))(())))) converges multivariate Gaussian with mean μ having entries Since we can
Nn = =; E,,, repeat the same argument for any collection of points I1, T22 .....--. ,N,, we must deduce that g is 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, {u(x) = Σ, [[) = Σ, m. [T)())) and k(x,x') Σ, (z)), = Σ,
and covariance E with entries Σnm = E. (x,, hy
Conceivably the CLT could be violated during training if a large enough fraction of Ma converged to 0 during training, or also if the separate channels became more strongly correlated. We investigate the non- Gon-siasiyyy of the features in section an7 and find that the features are mostly Gaussian but have somewhat longer tails indicating deviation from the theoretical behavior.
21

---

J. MOMENTS OF RECTIFIED GAUSSIAN RANDOM VARIABLES
Let f ~~ (N, Σ) be d dimensional random Gaussian vector. We compute here
1
[ I' fi ... -f
B(ReLLU(A))-. ... (eel(((a))))) Nv f>o
14)
(33)
Ny = (2 det det(E) 1//
2/2
(34)
We use the generating function technique. Define
1 [ d' ((f-1)² 11)EE'1(f - p1) + fl
Z(b) Nx $20
(35)
1 cb"
d' d'fe 11--"""''' + fb
(36)
No
1-2-1
Then
B(ReLLU(I)) RelLU(a)) 0
0
Db, abd %(b)
(37)
b 0-0
To compute Z(b) we proceed as in Gaussian case. We change variables to f= + += g,
(38)
and define μ + Eb to get:
Z(b)= .b #+ 6 Eb 1
e'geepll33 1/- E gg
(39)
NE gK+z
e (()) (()) (z; 0, Σ), S(b) = bb p + 2b Σb.
(40)
i being the multivariate standard Normal CDF:
(d) (zi p4, E)- = [
(41)
gg+z
(4) (g) p4, Σ)
expl- -00 - 11) " (g - [+]]-
(42)
NE
Now we compute the first two derivatives. Note that in d 1, denoting 0² E: (2, 0²²))
(1)) (2,0,0²) = (0))
(43)
dz
In d 2, we can use the conditional probability decomposition to get the required derivatives:
(e) 0, Σ) و6 (g1: a192, 81) ws (g2: O, E22)
(44)
a1 Σng ,31 = En - En2
(45)
0_ $ (2) (z, 0, Σ) 1 (22; 0, M222)
(46)
=201 O. Exe) $ (21, 01 11)
(47)
8. 2.000 (2) (2,0,)
22 (( (22) O, Eee 1 (21, 31)
(48)
+ +b(()) (22; 0, Eee 0s, + (2), 2122
(49)
02 0,,0 (2) (2) (z, E) (1;) (2;
(50)
2. So denoting a Tb. we get:
0, = (wi + Σ 2i)bb))(()) (60(6))++55) +c Edn dz (1)) 0, E)>ie
(51)
mb((b)
0g0,Z(b) = Exx + (H4 + Σ (bbb))0aZ(b)
(52)
+ (H4 + Σ (ms + -e s(b)
Σ
Dr 2, (2)00 E,),,
(53)
lgg
22

---

In partieular, for the first moment d = 1 we have:
E(ReLU(f)) -= (w: (,:0, + + (,,00000 20)
(54)
which coincides with equation 8 Note that æb()) (m,, 0,02) = 1 1/4//0 0; 1) because of the normal-
ization factor. For the second moments = 2 we have:
((ReL))((())eelU((3)) E12 (2) (1; Σ) + (1 (1:: 0, Σ) (1m2(0) + (2mm1 (0) (55)
+ Σ Σe, 2, 0= (2 (2, 0, E)bb=
(56)
1,0112
which can be rewritten in the form of equation
23


Reviews:
[{"review_id": "suIt9BqTLEh", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary: 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.\n\nPositives:\n1. The problem addressed is very important in many domains. \n2. The proposed approach is novel for learning from irregularly sampled data.\n3. Empirical results show that proposed model performs well.  \n\nConcerns:\n1. Given the time complexity of Gaussian processes, could the authors comment on the run-time of the proposed method as compared to standard CNN?\n2. This work is missing comparisons with one of the recent work focusing on the same problem (Rubanova et al, 2019). It would be interesting to also compare the run-time of both these methods. \n3. In Section 2 para 3, the authors mention that Gaussian process adapter (Li and Marlin, 2016) operates on a deterministic signal\non a regular grid and cannot reason probabilistically about discretization errors. I don't think this is true. Infact, GP adapter is an uncertainty aware classification framework which has the ability to pass the uncertainty through to classification framework similar to the proposed approach. This makes these two approach more similar in what they are trying to achieve and I would expect the authors to compare the performance of these methods.\n5. The authors repeatedly discuss that their approach does not require any discretization unlike some the previous approaches based on GP(mentioned in Section 2). But in Section 4.7, they mention that their approach requires to evaluate Gaussian Process at a set of points which is similar to what was done in Li and Marlin(2016). Could the authors clarify this?\n\nAdditional Comments:\n1. Could the authors comment on the possibility of using the proposed approach for interpolation in irregularly sampled time series?\n\nI would be happy to increase the score once my concerns are addressed.", "extra_scores": {}, "extra_text": {}}, {"review_id": "qHYyjuKlM3O", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper is way out of my comfort zone so I can only provide very general high level comments.\n\nFirst, the motivation for the paper seems very clear. For deep learning to work, we often have to discretize, and have to do so at lower resolution, which induces errors we typically ignore.\n\nSecond, the architecture is well described and illustrated. \n\nSome references (e.g. p1 Cockayne 2019) are in the wrong format. \n\nI didn't check the math in all details, but the three experiments with superpixel resolution, generalization to new resultions and and irregular time series appear convincing.", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ui7KPwMf08U", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "SUMMARY:\n\nThe 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.\n\n\n##################################################################\n\nREASON FOR SCORE:\n\nI like the motivation for this work and the paper was generally a pleasure to read. Naturally, many technical challenges arise when trying to approach this problem from a probabilistic perspective; they are discussed in a fair bit of detail and generally seem reasonable. However, I still have a number of questions that are listed below that would help me better understand the main ideas and contributions.\n\n\n##################################################################\n\nPROS:\n+ The motivation for the work seems sound.\n+ The main technical contribution seems to be the definition of a continuous convolution on GPs with an RBF kernel. This is presented clearly and in a fair bit of detail.\n+ The writing quality is very good. The appendix was very informative as well.\n\n\n##################################################################\n\nCONS:\n- The paper makes quite a few approximations (which mostly make sense) but doesn't give a thorough discussion of the effect of these choices.\n- The paper doesn't seem to make full use of the setup (see questions below).\n- I found Fig.2, especially 1 and 3, hard to interpret and the figure caption wasn't very informative. It appears from the left figure that the mean is very smooth in all layers. Is that desirable? In the right figure, what is gt on the y-axis?\n\n\n##################################################################\n\nQUESTIONS and COMMENTS:\n- My understanding is that PN approaches are typically used for tasks like, (a) to impose specific priors, (b) to improve the design of the problem or (c) to learn something about the implicit priors in established deterministic methods.\nIt would be interesting to think about these in the context of your model. For example, (a) is it reasonable to use a GP with a smooth kernel as an interpolant of the inputs? If the kernel choice was not limited by the analytic tractability, would you still choose an RBF kernel? For (b), could you use the model to decide which points in the input domain are the most informative for the end task? I suppose that you would need the full posterior distribution rather than just a point estimate, which would be very tricky given the multi-layer setup. Also, given that your approach seems to outperform existing (deterministic?) models or some tasks, how could you change the existing models to gain a similar advantage? \n\n- It would be interesting to know what is the effect of the various approximations, for example, the Gaussian approximation in 4.6 and the MC in 4.7, on the final task.\n\n- The CLT argument in Appendix I is not very clear to me. The Lyapunov CLT holds in cases of RVs being independent, but not necessarily identically distributed. Is it correct that in your case the RVs are weakly-dependent and identically distributed? What is c in the inequality |c - c'| >> 1? More generally, is the argument in this section some well-known result? It shares some ideas with the argument of Neal [1996] of an infinite limit of NNs being equivalent to GPs but it would be good to either give a more thorough discussion of this argument or provide links to related results in the literature.\n\n- It would be informative to discuss the computational and memory complexity of the proposed approach and other related questions, e.g. the number of samples in the MC steps.\n\n- I'm a little surprised that removing the noise gives worse results in your MNIST experiment. I would guess that this may be more related to some optimisation phenomenon than to actual model. For example, if instead of using the data-dependent noise parameter as described in Sec. 4.7 you used a spherical Gaussian noise with fixed sigma for all inputs, would that give better results than discarding the noise entirely? And how would it compare with your data-dependent noise? Using spherical Gaussian noise relates with some of the practices in DGP models where the input to each layer is assumed to be a noisy version of the outputs of the previous layer. \n\n- Related literature: It might be interesting to compare to some DGP models that also use convolutions, e.g. Deep Gaussian Processes with Convolutional Kernels [2018] by Kumar et al and Deep convolutional Gaussian processes [2019] by Blomqvist et al. \n\n##################################################################\n\nMINOR:\n\n- p. 1 - Probalistic Numeric\n- p. 6 - the number the number\n- p. 15 - infinitesmal\n- p. 18 - many many observations\n- p. 19 - a variant a variant", "extra_scores": {}, "extra_text": {}}, {"review_id": "Kl6Ds3MPRZA", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "# Summary\n\nTo work with data that is not sampled on a grid, this work represents the activations of a neural network using a Gaussian process with RBF covariance. Convolutions can be applied in closed form, but nonlinearities get projected to another RBF-posterior GP by\n\n1) calculating the first two moments of the post-nonlinearity activation, and\n2) projecting the moment-matched GP to a finite set of noisy (independent) observations,\n\nwhich can then be fitted with an RBF GP again, and so on until the end. The proposed method works very well for irregularly-sampled MNIST, and fairly well for the PhysioNet2012 data set for ICU mortality.\n\n# Comments\n\nThe paper references related work, makes a meaningful contribution, and I think the empirical methodology is sound. Thus, I think it should be accepted. As negatives, there are a few small issues with the paper, and it is hard to follow.\n\n## Central limit theorem\n\nThe CLT argument referenced in section 4.6, only works if all the elements of M are equal. (As a counterexample: suppose only one column of M is 1, and all the other ones are 0. Then, if the corresponding input process is non-Gaussian, all the outputs will be non-Gaussian).\n\nIn fact, the argument in the appendix forgoes coefficients for every $g_c$ altogether. It is correct, but it is not the same object in section 4.6.\n\nI don't buy your CLT argument, but I think the rest of the paper is still great without it. Simply, you're doing a moment-matching approximation, like Expectation Propagation (Minka, 2001b).\n\nI would consider dropping the CLT claims. Or, you could argue that in practice the elements of M are similar.\n\nYou could also argue that the elements of M are close to independent with mean zero. If M is also fairly large, you can argue that $f$ is close to a GP, following the theoretical arguments from the infinitely wide NN literature (Matthews et al., 2018; or generalizations thereof such as Yang, 2019).\n\nMinka, T. P. (2001). Expectation propagation for approximate Bayesian inference. UAI (pp. 362-369)\nMatthews, A G de G, et al. (2018). Gaussian process behaviour in wide deep networks. https://arxiv.org/abs/1804.11271\nYang, G., (2019). http://papers.neurips.cc/paper/9186-wide-feedforward-or-recurrent-neural-networks-of-any-architecture-are-gaussian-processes\n\n## Some experiment and figure clarification\n\nIn section 4.7 you talk about a collection of points ${x_i}_{i=1}^N$ to which you project the activation. Are these points the same for every layer? Are they always the points at which the network is currently being evaluated (I understand these are points for one input point, i.e. sampling locations for an MNIST image)?\n\nFor figure 2, how many test images do you use? Is it one, with varying number of sampling locations per image?\n\nFor Figure 3, how do you change the test-time resolution of a normal CNN, one that has discretely specified pixels? Do you linearly interpolate their values? This should be explained.\n\n## Other clarification requests and suggestions\n\nI found the paper hard to follow, and some things are unclear. This is understandable, because it uses mathematical tools unusual for an ML paper, but you can still make it easier for the reader. The rest of my review is some requests for clarity and writing suggestions.\n\n*Taylor expansion above eq. 3:* it is an expansion around $x$ of $g(x + a)$, no? I believe stating this simple fact after \"Verify by Taylor expansion [around $x$ of $g(x + a)$] would make it much less confusing.\n\nWhat is the \"[a, b]\" operator at the top of page 5?\n\n*Eq 9*:  I think $\\nabla\\nabla^T {\\boldsymbol \\Phi}({\\boldsymbol \\mu}, {\\boldsymbol \\Sigma})$ is supposed to represent the Hessian of ${\\boldsymbol \\Phi}({\\boldsymbol \\mu}, {\\boldsymbol \\Sigma})$, but I'm not sure. It's not quite standard notation and should be mentioned in the paragraph above.\n\nIn section 4.7, when projecting the current activations to a set of points, you lose the covariance. It's good that you comment on this in Appendix H, and it would be better if you pointed the reader to it.\n\nI think specifying that you use Adam in section 4.8 is helpful to the reader and doesn't take much space.\n\nAppendix F: $3 \\cdot 10^{-3}$, not $310^{-3}$.\n\nAre you going to release the code?", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This is a fairly technical paper bridging deep learning with uncertainty propagation in computations (i.e. probabilistic numerics). It is well structured, but it could benefit from further improvements in readability given that there are only very few researchers that are experts in all sub-domains associated with this work. Given the above, as well as low overall confidence by the reviewers, I attempted a more thorough reading of the paper (even if not an expert myself), and I was also happy to see that the discussion clarified important points. Overall, the idea is novel, convincing and seems well executed, with good results. The technical advancements needed to make the idea work are fairly complicated and are appreciated as contributions, because they are expected to be useful in other applications too (beyond irregular sampled data) where uncertainty propagation matters.
Decision: Accept (Poster)