Gaussian Process Predictions with Uncertain Inputs Enabled by Uncertainty-Tracking Microprocessors

Venue: NeurIPS_cc_2024_Workshop_MLNCP
Published: 2024-09-12 12:59:05

Abstract: Gaussian Processes (GPs) are a theoretically grounded mechanism for modelling with both aleatoric and epistemic uncertainty, but, the well-known solutions of the predictive GP posterior distribution apply only for deterministic inputs. If however, the input is uncertain, closed-form solutions aren't generally available and approximation schemes, such as moment-matching and Monte Carlo simulation, must be used. Moment-matching is only available under restricted conditions on the input distribution and the GP prior, and will miss the nuances of the posterior distribution; Monte Carlo simulation can be computationally expensive. In this article, we present a general method that uses a novel uncertainty-tracking microprocessor, Laplace, to implicitly calculate the posterior predictive distribution with uncertain inputs. We show that our method captures the nuances of the posterior distribution while being orders of magnitude faster than Monte Carlo simulation.

Main paper: 
 introduction: 
 In the domain of Machine Learning (ML), Gaussian Processes (GPs) are a class of non-parametric models that can intrinsically capture model uncertainty. They are also able to make predictions that capture the impact of uncertainty in inputs on the predicted output. The prediction of a GP, for a unobserved input, is referred to as the predictive posterior distribution since it is a distribution over the predictive output conditioned on the observed data.For deterministic inputs, there is a closed-form solution for the predictive posterior distribution of a Gaussian Process; this posterior is itself also a Gaussian distribution. However, when there is uncertainty in the input, the predictive posterior distribution of a Gaussian process is in general arbitrary, and no general closed-form solution exists. A typical solution to handle uncertain inputs is to carry out moment-matching, where the first n moments of the true predictive posterior distribution are calculated. Often, only the first two moments (mean and variance) are calculated [1, 2] . Closedform solutions for the first two moments are only available in special cases, and even then, such a Gaussian approximation loses key statistical properties of the true posterior distribution, such as any additional modes. To obtain a better approximation, one would carry out a (potentially) computationally expensive Monte Carlo simulation. We note that uncertain inputs are not an obscure case! They occur whenever a GP is used to model a dynamical system, where the output of one time step (always uncertain for a GP) is used as the input to the next time step.In this article, we present a simple algorithm for implicitly computing the predictive posterior distribution of a GP with an uncertain input using recent advances in uncertainty-tracking microprocessors [3, 4] . This method does not rely on Monte Carlo simulation but rather exploits the internal distributional representation of the uncertainty-tracking microprocessor. Our method is agnostic to the distribution of the input and the Gaussian Process prior.We evaluate the effectiveness of our method against the Monte Carlo method by comparing the accuracies of the solutions to the ground-truth (obtained by a Monte Carlo simulation with a large number of samples). We also compare the run time of each method. We show in Section 6 that using Laplace can achieve orders of magnitude faster performance while achieving the same accuracy. 
gaussian processes: 
 A Gaussian Process is a collection of random variables (Definition A.2 in Appendix A) where any finite subcollection of those random variables is jointly distributed as a multivariate Gaussian. It is the distribution over the infinitely-many random variables that satisfy this condition on any finite subcollection of them. We can define a Gaussian Process implicitly by specifying a mechanism for deriving the mean vector and covariance matrix of the multivariate Gaussian distribution that distributes any given subset of random variables. We can do so using a mean function m and covariance kernel k to derive the mean vector and covariance matrix of the multivariate Gaussian distribution over a finite subset of random variables using the indexes of the random variables. Definition 2.1 (Implicit definition of a Gaussian Process). Let X be an index set. Then, for a finite subset of X, {x 1 , ..., x n } ∈ X, let the random variables Y x1 , ..., Y xn be jointly distributed as a multivariate Gaussian, N (µ, Σ), where the elements µ i of the mean vector µ and Σ ij of the covariance matrix Σ are implicitly derived using:µ i = m(x i ),andΣ ij = k(x i , x j ).The function m : X → R is called the mean function and the function k : X × X → R is the covariance kernel. We denote an implicitly defined Gaussian Process by GP(m, k).We denoted the index set with X and the random variable with Y xi , for x i ∈ X to highlight that a Gaussian Process is a distribution over functions where X is the input space and Y = Y xi∈X is the output space. Thus the sample space of each Y xi is Y . Let X = R and Y = R be an input and output space respectively and F denote a space of functions, where (f ∈ F) : X → Y denotes a possible function from X to Y . We can place a Gaussian Process prior over F by writing F ∼ GP(m, k). In this view, a sample from a Gaussian Process is a sample from all the random variables Y x = f (x), ∀x ∈ X.The function space F over which a Gaussian Process is defined is implicitly described by the chosen mean function and the covariance kernel; we have placed a GP prior over F. Popular kernels include the squared exponential (this is a basic 'all-rounder' type kernel) and the Matern class of kernels [5] .The squared exponential kernel is given by:EQUATIONwhere c ∈ R + and λ ∈ R are hyperparameters of the kernel. c is typically called the signal variance, and the λ is called the characteristic lengthscale 1 . The squared exponential kernel describes a function space of smooth functions that have a variability that is described by the signal variance and the lengthscale; the signal variable loosely specifies how much the function can vary from the mean function, and the lengthscale specifies how quickly the function can vary from the mean function.Let X = R d denote an input space, where d denotes the dimension of X, and let Y = R denote an output space. Let f : X → Y be an arbitrary function from X to Y from a function space F. We define a Gaussian process prior over F as F ∼ GP(m, k). The mean function m and covariance kernel k implicitly define the function space F. Further, let the likelihood function be ỹ = f (x) + ϵ, where ϵ ∈ N (0, σ 2 ϵ ) and ỹ ∈ ( Ỹ = Y ) denotes a noisy observation. Finally, let X ∈ R k×d be a matrix of k input points from X, and let Y ∈ R k×1 be a matrix of k noisy output points from Ỹ .We are interested in the posterior Gaussian Process over F given the data (X, Y). However, rather than finding an explicit distribution over F|X, Y, we express the posterior distribution of the posterior Gaussian Process evaluated at particular point x * ∈ X, which according to the definition of a Gaussian Process is a Gaussian random variable. Then, for an arbitrary point x * ∈ X, the predictive posterior density of the random variable Y x * = f (x * ) is given by:EQUATIONwhere,EQUATIONσ 2 * : X → R x * → σ 2 * (x * ) = k(x * , x * ) -k(x * , X)(K + σ 2 ϵ I) -1 k(x * , X) T .For a deterministic input, Gaussian Process prediction can be thought of as a function from the input to a mean and a variance, as Equation 3shows, that defines the predictive posterior distribution. 
gaussian process prediction with uncertain inputs: 
 A useful enhancement is to consider a distributional input. Let the input to the Gaussian Process now be a random variable X * ; the sample space of X * = X. Thus, X * ∼ p X * ( • ). We are still interested in the predictive distribution of Y x * ∼ p f (X * ) (y|X, Y), where x * ∈ X * . However, we must now marginalize over X * because it is a random variable:EQUATIONIn general, this integral is intractable. Therefore, we must usually resort to approximation methods.We discuss two approximation methods used in the literature below. In Section 4, we introduce a novel method for implicitly computing the posterior distribution with uncertain input using a state-of-the-art uncertainty-tracking processor, Laplace [3, 4] .When the input random variable is a Gaussian random variable we can approximate the predictive posterior distribution as a Gaussian distribution. This is done by defining a Gaussian distribution using the first and second moments of the predictive posterior distribution [1, 2] . This mean and variance can be computed analytically when using the squared exponential kernel and the input is a Gaussian random variable [1, 2] . However, for more complicated input distributions and GP priors, it is difficult to find closed-form solutions. Furthermore, a Gaussian approximation to a non-Gaussian random variable would lack important statistical information.The Monte Carlo method can be used to obtain a more informative approximation of p f (X * ) (y * |X, Y).Monte Carlo simulation [6, 7, 8] applies the Monte Carlo method to find an approximation to the probability density function of transformed random variable. That is, if we want the probability density function of Y = f (X), then we first take a n samples of x i ∼ X, transform each sample by appylying f to obtain n samples y i = f (x i ) of Y , and create a density estimation (such as a histogram) using the transformed samples. The key variable that is changed to obtain better approximations is the number of samples n.Girard et al [2] suggests estimating the integral in Equation 4 using Monte Carlo integration, where,EQUATIONFor this method, we would need to evaluate the sum in Equation 5for each value that f (x * ) can take. Instead, we will obtain samples of f (x * ) under p f (x * ) (y * |X, Y) by taking samples of x * i ∼ X * , calculating the mean µ * (x * i ) and variance σ 2 * (x * i ) according to Equation 3 for p f (X * ) (y * |X, Y, x * ),and taking samplesy * i ∼ N (µ * (x * i ), σ 2 * (x * i )).The Monte Carlo method can be used to obtain an arbitrarily good approximation of p f (X * ) (y * |X, Y).However, doing so can be prohibitively expensive due to the slow convergence rate of the method [8] . 
uncertain gaussian process prediction using laplace: 
 To overcome the accuracy limitations of moment-matching and the computational costs of Monte Carlo simulations, we present a novel method that is enabled by recent advances in uncertaintytracking microprocessors. Laplace [3, 4] is a computer microarchitecture that can represent distributional information in its microarchitectural state and track how these distributions evolve under arithmetic operations, transparently to the applications running on it. Laplace does this by providing an in-processor representation of the probability distribution, as well as a suite of arithmetic operations that can be used to manipulate these representations. Laplace carries out deterministic computations on probability distributions and does not repl on any sampling methods. Therefore, unlike Monte Carlo methods that require convergence to a good solution by increasing the number of Monte Carlo iterations, Laplace is convergence-oblivious up-to the fidelity of its representation [9] .Laplace has an in-processor distribution representation and this representation has an associated representation size, r, which describes the precision of the probability distribution representation.Larger representation sizes result in more accurate representations. A useful analogy is the IEEE-754 standard for representing floating-point numbers [10, 11] for approximately representing the infinite set of real numbers in a computer that can only represent a finite number of values.A practical way of thinking of what Laplace does is to think of it as a way of carrying out uncertainty propagation. That is, if we have a random variable X and a function f : X → Y , then Laplace directly computes the distribution of Y = f (X) without resorting to Monte Carlo sampling. If X is represented in the Laplace's representation for distributions, then we can compute the representation of the distribution of Y by directly computing it through the function f .Since Laplace can represent any distribution for X, it would be useful if we were able to use Laplace to compute the predictive posterior distribution of a GP with an uncertain input. A Gaussian Process is practically two functions that map some input x * to the mean µ * (x * ) and variance σ 2 * (x * ) of the output distribution at x * (Equation 3), and not to the output random variable itself. Simply plugging a Laplace-represented distribution to Equation 3 would result in a distribution of means and variances, and not a distribution of the output random variable that we desire.We can express a transformation of random variables that results in an output random variable that has the desired posterior distribution p f (X * ) (y * |X, Y) by using the "reparametrization trick" [12] , EQUATIONp f (X * ) (y * |X, Y) = p f (x * ) (y * |X, Y, x * )p X (x * )dx * = N µ * (x * ), σ 2 * (x * ) p X (x * )dx * , where p f (X * ) (y * |X, Y, x *) is the conditional probability density function of f (x * ) given X, Y, and a particular value of x * ∈ X * .Proof. See Appendix B. ■Using Algorithm 1, we can easily find the uncertain Gaussian Process predictive posterior when using Laplace. Here, X * , ϵ, and y are distributional variables that are represented in Laplace.Require:X * ⇐= DistributionalVariable, X, Y m ← µ * (X * ) ▷ From Equation 3 v ← σ 2 * (X * ) ▷ From Equation 3 ϵ ← GaussianDistributionalVariable(0, 1) ▷ A distributional variable in Laplace y ← m + √ vϵ return y 5 MethodsWe compare the performance of Monte Carlo simulation and Algorithm 1 for computing the predictive posterior distribution. Our methodology was closely inspired by Petangoda et al [9] . We measure the timing performance in terms of the wall-clock run time of each method and we measure accuracyWasserstein distance [13] of the result of each method compared to a ground-truth predictive posterior distribution obtained by a large (n = 1, 000, 000) iteration Monte Carlo simulation.We compare the run times and the Wasserstein distances between the Monte Carlo simulations and the Laplace-based method across different number of Monte Carlo iterations n and representation sizes r respectively to compare the trade-offs between converging to the output distribution and the extra time required to do so. For each configuration of n or r, we carry out 30 repetitions to account for variation in the process of sampling 2 .We set the prior mean function to m(x) = 0. We set X = R and Y = R, and generate a dataset (X, Y) of k = 10 data points by setting x 0 = -3π, x i = x i-1 + 0.6π 3 , and y i as,EQUATIONwhere each ϵ i ∼ N (0, 0.1) 4 . For each case, we use the squared exponential kernel k sq c=1,λ= √ 1.7 and we set the uncertain input distribution to X * ∼ N (0, 4).For the Monte Carlo method, we evaluate on n ∈ {4, 256, 1152, 2048, 4096, 8192, 16000, 32000, 128000, 512000}. For Laplace, we evaluate on r ∈ {16, 32, 64, 128, 256, 512, 2048}. To run our application on Laplace, we use a multi-dimensional variation of the Laplace [3, 4] microprocessor with automatic correlation tracking between arbitrary random variables, as provided by the commercial implementation of Laplace [14] . We refer to this variant of the microprocessor as Laplace (CT).We implemented the experiments described above using the C programming language and a custom tensor library implemented in C 5 . Please see Appendix C for more detail on our method. 
results: 
 Figure 1 summarizes our results; Table 2 in the Appendix shows numerical results for the other configurations. Figure 1a shows that the results from Laplace are almost always at the Pareto frontier.Therefore, for any required accuracy except for the highest, in the ranges that we had tested, Laplace produces an output orders of magnitude faster. Table 1 shows that for r = 128, the Monte Carlo 2 We calculate the Wasserstein distance for Laplace's representation of the output distributions by generating 1,000,000 samples from the representation. Therefore, we also repeat the Laplace experiments 30 times for each representation size even though Laplace's uncertainty-tracking methods are deterministic and convergenceoblivious. 3 Essentially, we have equally divided the range [-3π, 3π] to n = 10 points. 4 We set the random seed to 42. 5 We will release all code necessary to reproduce our experiments on Github method with n = 128, 000 shows a nearly identical Wasserstein distance (albeit with twice the standard deviation) while taking approximately 151.15× longer to compute.Furthermore, Laplace does not suffer from large variances in the accuracy since it carries out deterministic computation on probability distributions. We can see this from the constant variance shown by the Laplace result. This constant uncertainty is due to the fact that we calculated the Wasserstein distance by taking finite samples from the output distribution. It is possible to calculate the Wasserstein distance directly from Laplace's representation, but we did not do so at this time.After r = 128, the accuracy of the Laplace result does not increase, and the trade-off with run time increases. Therefore, if larger accuracy is required, the Monte Carlo method must be used.Figures 1b and 1c show that the distribution of the output is more complex than a Gaussian. A moment-matched Gaussian solution would, for example, miss the fact that this posterior distribution contains several modes. 
conclusions: 
 Existing methods for approximating the often-intractable Gaussian Process predictive posterior distribution with uncertain inputs can be inaccurate and constrained in their use (moment-matching), or computationally expensive (Monte Carlo simulation). We present a novel and simple method for computing the desired posterior distribution on an uncertainty-tracking microprocessor by making use of the reparametrization trick. Our method can be used on any input distribution and GP prior unlike moment-matching, and we show experimentally that our method can be accurate and orders of magnitude faster compared to Monte Carlo simulation. 
 
 Appendix: In this section, we define some key mathematical details. We follow conventions and techniques as presented by Petangoda et al [9] .A.1 Random Variables Definition A.1 (Probability density function). Let X be a set called the sample space and p X be a map from X to R + ,p X : X → R + , that satisfies: X p X (x) dx = 1.We define p X to be the probability density function on X.Definition A.2 (Random Variable). Let X be a set and p X a probability density function on X. We define the tuple (X, p X ) as a random variable 6 .In the rest of this article, we suppress the notation (X, p X ) and simply write X to denote a random variable. We use the simple definition of a random variable above rather than the typical measuretheoretic definition to capture what we think of as a random variable in practice; we often think of a random variable in terms of its probability density function (for example, when we think of a Gaussian random variable, we often think of the bell curve drawn by the Gaussian probability density function).Often, the set X = R d , where R d is the d-dimensional real space. In the integration above (and below), we have suppressed the fact that it contains multiple integrals over each dimension of X.The probability distribution of a random variable is a set function P X : Ω → R + that tells us about the probability of the random variable taking on the values inside the given set, where Ω is a set of subsets of X. Ω must technically be a σ-algebra, which for our purposes, can be thought of as a set of sets each of which can be expressed as a countable union of disjoint sets that also belong to Ω7 .Since we are only considering random variables that contain probability density functions, we define the distribution of a random variable as:Definition A.3 (Distribution of a Random Variable). Given a random variable X with probability density function p X , the distribution of X is given by the set functionEQUATIONwhere U i is a collection of disjoint sets whose union equals the input set ω.Definition A.4 (Cumulative Distribution Function of a Random Variable). The cumulative distribution function F X of the random variable X with probability density function p X is given by,F X : X → [0, 1]x → F X (x) = P({t|t ∈ X and t < x}) = {t|t∈X and t<x} p X (t) dt, (9) Note that in the integrations above, we have suppressed the fact that the integrations need to occur over dimentions of X. By applying the Fundamental Theorem of Calculus to each such integral, we obtain that p X (x) = ∂ d F X ∂x1...∂x d (x) This also gives us a way of going from a distribution of a random variable to its probability density function. Therefore, in this article, we often refer to the distribution and the probability density function interchangeably.An example of a high dimensional random variable is the multivariate Gaussian random variable:Definition A.5 (Multivariate Gaussian random variable). Let µ ∈ R d and Σ ∈ R d×d be a mean vector and a covariance matrix respectively. The multivariate Gaussian random variable (X N , p N X ) is defined as:X N = R d , p N X : R d → R + , x → 1 (2π) d det Σ exp - 1 2 (x -µ) T Σ -1 (x -µ) .p N X is called the multivariate Gaussian distribution. A Gaussian distribution is fully specifed by its mean vector µ and its covariance matrix Σ. Therefore, we can also write that a random variable X N is distributed as a multivariate Gaussian by writingX N ∼ N (µ, Σ).The mean vector µ is the expectation of the random variable X N , which is defined for any random variable X as, Definition A.6 (Expectation of a Random Variable). Let (X, p X ) be a random variable. The expectation of X is defined as:E[X] = X xp X (x) dx.An alternative view of a multivariate Gaussian random variable, or any other random variable over R d , is that it is a collection of univariate (single dimensional) random variables that are jointly distributed as the multivariate distribution. The covariance matrix Σ is a matrix the elements Σij = cov[X i , X j ],where X i and X j are the i-th and j-th univariate random variables of X respectively. Generally, the covariance of two random variables is defined as:Definition A.7 (Covariance of Random Variables). Let (X, p X ) and (Y, p Y ) be two random variables.The covariance of X and Y is defined as:cov[X, Y ] = E[(X -E[X])(Y -E[Y ])].We note that the covariance between the same random variable is its variance:cov[X, X] = E[(X -E[X]) 2 ] = var[X].In this article, we write cov[X] to denote the variance of X.Give a collection of random variables {X i } i=1,...,n , we can construct a covariance matrix Σ:Definition A.8 (Covariance Matrix). Let X be an index set and {Y i }x = 1 n be a collection of random variables indexed by x ∈ X. Then, the covariance matrix Σ ∈ R n×n of the collection of random variables that contains the pairwise covariances of the random variables:Σ ij = cov[Y i , Y j ].Let f : X → Y denote a transformation from a random variable X to a random variable Y (i.e., Y = f (X)) 8 . We can also apply f to an instance value x of X to obtain an instance value y= f (x) of Y .If f is invertible and once-differentiable, then Theorem A. Two matrix inputs: k(X, X ′ ) is the kernel applied to two sets of matrices obtained from sets {x i } i=1,...,n and {x ′ j } i=1,...,m collected into matrices X ∈ R n×d and X ′ ∈ R m×d respectively. The output is a matrix k(X, X ′ ) ∈ R n×m , where each element k(X, X ′ ) ij = k(x i , x ′ j ). Our methods follow closely from Petangoda et al [9] . We measure the run time as the sum of the time taken to generate samples incurred during the sampling step of the Monte Carlo method or the initializing step of Laplace, and the time taken for the evaluation step. For both methods, we measure time using the gettimeofday function from the Standard C library [17] . We measure the time from the start of the main entry point until the end of key computations. The reported times omit any time spent by the programs on saving and reporting the results.In order to explicitly quantify the post-processing step of the Monte Carlo method, we compute the mean and the variance of the samples obtained from the Monte-Carlo-based experiments. Such a step is not necessary with Laplace, because Laplace already provides a usable representation of the output distribution. We note that we are being generous to the Monte Carlo method, since the mean and the variance alone does not fully capture the shape of a non-Gaussian distribution. In contrast, Laplace captures the full distribution in its representation.It is possible to write the squared exponential kernel in terms of a specific lengthscale for each dimension of the input space. For simplicity, will be using a single, common lengthscale for all dimensions in this article, as in Equation1.There can be measure theoretic random variables, such as the random variable that has the Cantor distribution, that do not admit probability density functions[15]. In this article, we do not consider such exotic random variables.The sets in a σ-algebra must also satisfy that countable unions and intersections also belong to Ω. 



Reviews: review: title: title: Review of Paper Submission 50
main_review: main_review: This paper describes an application of uncertainty-tracking microprocessor on gaussian process predication.
It shows its advantages of high prediction accuracy and faster runtime against Monte Carlo simulation.
The writing is clear and easy to follow. Here is a list of pros and cons:

**Pros:**
* This paper identifies the limitation and challenges of exisiting methods of uncertain gaussian process predication.
* It theoretically demonstrates how this problem can be addressed by uncertainty-tracking microprocessors (Laplace) with proof.
* The results shows strong improvements in terms of runtime and prediction quality against Monte Carlo simulation, which implies its promising applications in practice.

**Cons:**
* The application of Laplace to solve uncertain gaussian process prediction needs more details to explain how it works. The current Section 4 only highlights the unique properties of Laplace to convey the high-level idea. As a reader, I have to go through the paper of reference [3] to have a much better understanding. 

* The implmentation setup and runtime measurement details of the experiment results needs more clarification. In Appendix, this paper mentions that it uses Apple M1 Pro as the processor to run Mont Carlo simulation without parallel execution. 
    * Which type of M1-Pro's core is used? Efficient core or Performance core?
    * Is M1-Pro run with battery mode or power plug-in mode?
    * What software language is used to implement Mont-Carlo simulation? C/C++, Python, or Matlab? (A reference to [9] is not helpful at all)

* The details of Laplace hardware is not clear. The referenced website [14] is not telling what kind of implementation of this microprocessor is used. Is it an emulated processor in RISC-V CPUs, x86 CPUs or FPGAs? 

* It will be more insightful to shows the comparison results on some real-world uncertain Gaussian process predication problems. The current problem is too small to stress any processor cache/memory sub-system bound. In particular, the paper will be stronger to show results on larger problems that need to store long chain of ancestor addresses.
title: title: Review for Gaussian Process Predictions with Uncertain Inputs Enabled by Uncertainty-Tracking Microprocessors
main_review: main_review: The paper discusses how uncertainty-tracking microprocessors can be used to efficiently compute posterior distributions with a Gaussian process prior applied to uncertain inputs. Typical previous approaches using Gaussian process priors required the input to be deterministic. For uncertain inputs, some approaches only approximated the posterior with moment-matching for the first few moments, but this often cannot represent many statistical properties of posterior distributions including additional modes. Alternatively, computationally-expensive Monte Carlo simulation can provide more accurate estimates of the posterior. In this work, the authors proposed to use a microarchitecture inherently suited for modelling distributions to make accurate estimates of the posterior distribution for uncertain inputs computationally tractable.

Before the methods section, the paper provides a detailed and very helpful derivation of the previous approaches for computing posterior distributions with Gaussian processes for deterministic and stochastic outputs. This section is well-written and effectively defines the vocabulary and related work in the field. It also clearly highlights the challenges with the existing moment-matching and Monte Carlo estimation methods.

Comparatively, the discussion on the Laplace microprocessor is lacking in details. While the cited papers do discuss how the microprocessor implements the distribution-to-distribution operations, it would be good to have at least a couple examples of arithmetic operations and distribution representations in this paper to provide a bit more specificity to the discussion on how the Laplace microprocessor works. As is, the current discussion on the Laplace processor seems a bit too vague.

It would also be useful to provide more detail on why the hardware configuration for the Monte Carlo method (Apple M1 Pro) is a suitable comparison for the Laplace microprocessor. Would a desktop GPU be a better estimate of a hardware that is well suited for Monte-Carlo simulation in a way that is similar to how the Laplace microprocessor is well-suited for the proposed approach? Providing timing for the Monte Carlo approach across a variety of hardware configurations would be another way to mitigate this concern. 

Overall, this paper provides a meaningful contribution by showcasing how uncertainty-tracking microprocessors could be very useful for communities interested in Gaussian Process-based uncertainty estimation.