Reversible and non-reversible Markov chain Monte Carlo algorithms for reservoir simulation problems

We compare numerically the performance of reversible and non-reversible Markov Chain Monte Carlo algorithms for high-dimensional oil reservoir problems; because of the nature of the problem at hand, the target measures from which we sample are supported on bounded domains. We compare two strategies to deal with bounded domains, namely reflecting proposals off the boundary and rejecting them when they fall outside of the domain. We observe that for complex high-dimensional problems, reflection mechanisms outperform rejection approaches and that the advantage of introducing non-reversibility in the Markov Chain employed for sampling is more and more visible as the dimension of the parameter space increases.


Introduction
Markov Chain Monte Carlo (MCMC) methods are popular algorithms which allow one to sample from a given target measure π on R N .In combination with the Bayesian inference approach, MCMC methods have been very successfully implemented in a vast range of problems in the applied sciences, and the literature about MCMC is extensive.The purpose of MCMC algorithms is to build a Markov chain {x k } k∈N which has the measure π as invariant measure.Traditionally this is obtained by ensuring that the chain satisfies the detailed balance condition with respect to the measure π, so that the resulting chains are reversible with respect to π.In recent years, non-reversible MCMC algorithms have attracted a lot of attention, because of their favourable convergence and mixing properties; the literature on the matter has rapidly become large, here we refer the reader to e.g.[14,15,2,3,7] and references therein; however, to the best of our knowledge most of the papers on non-reversible MCMC so far have tested this new class of algorithms only on relatively simple target measures.Furthermore, the performance of non-reversible algorithms has been discussed almost exclusively in the case in which the measure is supported on the whole of R N .However in many applications it is very important to be able to sample from measures supported on bounded domains.This is the case, for example, in applications to reservoir modelling and petroleum engineering, which we treat in this paper.The purpose of this paper is twofold: on the one hand we want to test the performance of non-reversible algorithms for complex, high-dimensional problems, which are completely out of reach for a full analytical treatment; on the other hand, we want to employ them for situations in which the target measure is supported in a bounded domain.The non-reversible algorithms that we consider in this paper are the so-called Horowitz algorithm, see [12], and the Second Order Langevin-Date: March 19, 2019.
Hamiltonian Monte Carlo (SOL-HMC) algorithm, introduced in [8].Both of them are nonreversible modifications of the well known Hamiltonian Monte Carlo (HMC) [13], which is reversible.More precisely, the Horowitz algorithm is a non-reversible version of HMC and the SOL-HMC algorithm is a modification of the Horowitz algorithm, well-posed in infinite dimensions and therefore well-adapted to sample from the high-dimensional target measures that we will treat here.
All the algorithms we discuss in this paper need in principle no modification in order to sample from measures supported on bounded domains.However, if they are not suitably modified, they will employ proposal moves which fall outside of the support of the target measure.For the problem we consider, this seems to give two major drawbacks, namely i) proposal moves that fall outside of the box are immediately rejected, so the algorithm wastes time rejecting moves which one knows a priori should not be made;1 ii) the likelihood function is calculated through the use of a simulator, which, further to being time-consuming to run, it will mostly fail to calculate correctly values that fall outside the support of the target.For this reason, we will consider two modifications of each one of the mentioned algorithms in which proposal moves that fall outside of the support of the target measures are either rejected or bounced off (or better, reflected off) the boundary of the support (see Section 2), so that the proposal is not automatically rejected as it will fall within the support.With these observations in mind, let us come to summarize the main conclusions of the paper: • We compare rejection and reflection strategies and test them on both low and high dimensional targets and conclude that, for the problems at hand, the two strategies perform similarly when implemented in low dimensions; however in high dimensions (and for more complex problems where a proxy is employed for the likelyhood function), reflections seem more advantageous • we compare the performance of HMC, Horowitz and SOL-HMC and conclude that, in high dimensions, the SOL-HMC algorithm is subtantially outperforming the other two.
Performance of all the algorithms is compared by using the normalized Effective Sample size (nESS) as a criterion for efficiency, see Section 4. We emphasize that one of the main purposes of this paper is to demonstrate how the SOL-HMC algorithm, while being a simple modification of the HMC algorithm, which requires truly minimal code adjustment with respect to HMC, can bring noticeable improvements with respect to the latter method; furthermore, such improvements are more noticeable when tackling high-dimensional complex target measures.The paper is organised as follows: in Section 2 we recall the HMC, SOL-HMC and Horowitz algorithms, present the numerical integrators that we use in order to implement such algorithms and introduce the versions of such methods which are adapted to sampling from measures with bounded support.In Section 3 we give details about the types of target measures used to compare the efficiency of these various algorithm and how they arise from reservoir simulation problems.This section explains mostly the mathematical structure of such target measures.Further details regarding the simulator and some basic background material about the reservoir model are deferred to Appendix B. In Section 4 we present numerical experiments.For completeness, we include Appendix A, containing some simple theoretical results regarding the modified algorithms.

Description of the algorithms
In this section we present the three main algorithms that we would like to compare, namely the Hamiltonian Monte Carlo (HMC) algorithm, the SOL-HMC algorithm and the Horowitz algorithm.With abuse of notation, throughout we will denote a probability measure and its density with the same letter, i.e. π(dx) = π(x)dx.
Suppose we wish to sample from a probability measure π defined on R N which has a density of the form π(x) ∝ e −V (x) e − x,C −1 x , i.e. the target density π is absolutely continuous with respect to a Gaussian measure with covariance matrix C (as customary, we assume that such a matrix is symmetric and positive definite).All three of our algorithms make use of the common trick of introducing an auxiliary variable p ∈ R N and sampling from the density π defined on the extended state space R N × R N as follows π(x, p) ∝ e −V (x) e − x,C −1 x e − 1 2 p 2 . (2.1) The original target measure π is therefore the marginal of π with respect to the variable x.
More precisely, the algorithms we will consider generate a chains which sample from the measure (2.1); because (2.1) is a product measure of our desired target times a standard Gaussian, if we consider just the first component of the chain {(x k , p k )} k , i.e. the chain {x k } k , then, for k large enough, such a chain will be sampling from the measure π.We now focus on explaining how the chain {(x k , p k )} ⊂ R N × R N is generated by the three algorithms we consider.
Let us introduce the Hamiltonian function then the associated Hamiltonian flow can be written as Let χ t denote a numerical integrator for the system (2.3) up to time t (we will comment below on our choices of integrator).The HMC algorithm then proceeds as follows: suppose that at time k the first component of the chain is x k and (1) pick p k ∼ N (0, I); (2) compute (x k+1 , pk+1 ) = χ δ (x k , p k ) and propose xk+1 as the next move; (3) calculate the acceptance probability α k , according to α k = min(1, e −(H(x k+1 ,p k+1 )−H(x k ,p k )) ; (4) Set q k+1 = qk+1 with probability α k , otherwise set q k+1 = q k .In principle any numerical integrator can be used in an HMC algorithm (see [17,13] for more detailed comments on this).In this paper we will consider two numerical integrators χ, which are two popular splitting schemes.The first is given by splitting the "momentum" and "position" equations, see e.g.[17] and references therein.That is, let M t denote the solution map at time t of the system (2.4) ẋ = 0, ṗ = −x − C∇V (x) and P t denote the solution map at time t of the system (2.5) ẋ = p ṗ = 0.
For HMC we shall always use the numerical integrator . Note that we can always write the maps M t and P t explicitly; indeed, The other splitting scheme that we will consider splits the Hamiltonian system (2.3) into its linear and nonlinear part.More precisely, let R t and Θ t be the flows associated with the following ODEs: (2.9) The resulting integrator is given by This is the integrator that we will use in the SOL-HMC algorithm (see step (1) of the SOL-HMC algorithm below); the use of this splitting scheme for high dimensional problems has been studied in [10].SOL-HMC is motivated as a time-discretisation of the SDE where {W t } t≥0 is a standard N -dimensional Brownian motion.Such an equation can be seen as a Hamiltonian dynamics perturbed by an Ornstein-Uhlenbeck process in the momentum variable.As is well known, the SDE (2.11) admits the measure (2.1) as unique invariant measure, see e.g.[18].With these observations in mind, define O ε to be the map which gives the solution at time ε of the system Note that we may solve this system explicitly, indeed where ξ is a standard normal random variable.In Section 4 we will set where i is a parameter we can tune; in which case we have The SOL-HMC algorithm is as follows: (1) Given (x k , p k ), let and propose (x k+1 , pk+1 ) = χ δ S (x k , pk ), where we recall that χ δ S is the integrator introduced in (2.10); (2) calculate the acceptance probability α k according to (3) set Finally, the algorithm that we will refer to as the Horowitz algorithm, is just the SOL-HMC algorithm when in step one, instead of using the integrator χ S , we use the integrator χ H (defined in (2.6)).
Remark 2.1.We do not give many details about HMC, SOL-HMC and the Horowitz algorithm here, and refer to the already cited literature.However we would like to stress the two following facts: • The chain {x k } k produced by the HMC algorithm is reversible with respect to the measure π, in the sense that it satisfies detailed balance with respect to π [13] -more precisely, the chain {(x k , p k )} k generated by HMC satisfies a generalised detailed balance condition with respect to π, see e.g.[17, Lemma 1] or [5].In contrast the chains generated by the Horowitz algorithm and the SOL-HMC do not satisfy any form of detailed balance with respect to π and they are therefore non-reversible, see [7,8].In Appendix A we will show that adding reflections to the algorithm does not alter this property.• The difference between the Horowitz algorithm and HMC may seem small, but in reality this small difference is crucial.Indeed, thanks to this choice of integrator, SOL-HMC is well-posed in infinite dimensions, while the Horowitz algorithm is not.For a discussion around this matter see [7,10].
As mentioned in the introduction, in this paper we will be interested in sampling from measures which are not necessarily supported on the whole space R N , but just on some box B = [−a, a] N .If this is the case then one may still use any one of the above algorithms and reject proposal moves that fall outside the box.We will briefly numerically analyse this possibility (see Section 4).Alternatively, one may want to simply make sure that all the proposed moves belong to the box B, so that the algorithm doesn't waste too much time rejecting the moves that fall outside the box.We therefore consider modified versions of the introduced algorithms by introducing appropriate reflections to ensure that all of the proposals belong to the box B. Because the proposal is defined through numerical integration of the Hamiltonian dynamics, we will need to modify the integrators χ H and χ S .
First consider the map P δ defined in (2.7); then we define map P δ bounce recursively as follows: (1 In which case P αδ (x, p) lies on the boundary of the box, so there exists some2 j ∈ {1, . . ., N } such that the jth component of P αδ (x, p) is either a or −a.Then we define (2.17) bounce (S j (P αδ (x, p)).Here S j is the reflection map S j (x, p) = (x, p 1 , . . ., p j−1 , −p j , p j+1 , . . ., p N ).Similarly we define R δ bounce by (1 In which case R αδ (x, p) lies on the boundary of the box, so there exists some j ∈ {1, . . ., N } such that the jth component of R αδ (x, p) is either a or −a.Then we define ).Note that it may occur that R δ (x, p) ∈ B however there is some point α ∈ [0, 1] such that R αδ (x, p) / ∈ B, in this case we still set R δ bounce (x, p) = R δ (x, p).Therefore the algorithm HMC-bounce (Horowitz-bounce, respectively) is defined like HMC (Horowitz, respectively), but the numerical χ δ H,Bounce = M δ/2 P δ bounce M δ/2 is employed in place of the integrator χ H ; analogously, the algorithm SOL-HMC-bounce is defined as the algorithm SOL-HMC with numerical integrator χ δ S,Bounce = Θ δ/2 R δ bounce Θ δ/2 in place of χ S .

Target measures
In this section we describe the three target measures that will be the object of our simulations.The first measure we consider, π Ros , is a change of measure from the popular 5D Rosenbrock, see (3.2).This is the most artificial example we consider.The other two target measures are posterior measures for parameter values in reservoir simulation models.Roughly speaking, the second target measure we describe, π f ull , is a posterior measure on the full set of parameters of interest in a reservoir model; for our reservoir model, which is quite realistic, we will be sampling from 338 parameters, hence this measure will be a measure supported on R 338 .The third measure, π light , is a measure on R 21 , which derives from considering a considerably simplified reservoir model.We will refer to the former as full reservoir model and to the latter as lightweight parametrization.In this section we explain the mathematical structure of π f ull and π light , without giving many details regarding the inverse problem related to the reservoir model.More details about the reservoir model and the simulator used to produce the likelihood function have been included in Appendix B for completeness.In the following we let I N denote the N × N identity matrix.Let us now come to describe our targets.
• Change of measure from 5D Rosenbrock (i.e.π Ros ).The first target measure we consider is a measure on R 5 and it is a change of measure from the 5D Rosenbrock measure; namely, the density of 5D Rosenbrock is given by The target we consider is given by where C is the prior covariance matrix.In all the numerical experiments (see Section 4) regarding π Ros we take C = 0.3 • I 5 .
• Full reservoir simulation.We study a Bayesian inverse problem for reservoir simulation.We consider a synthetic reservoir with 7 layers, and 40 producing/injecting wells.A diagrammatic representation of the reservoir is shown in Figure 1.Each layer is divided in blocks and, while the well goes through all the layers, it will not necessarily communicate through perforations with all the blocks it goes through (in Figure 1 we highlight in yellow the boxes containing a perforation of a well).We also assumed that, in each layer, the well goes through at most one block.In total our subsurface reservoir is made of 124 blocks: 38 blocks on the boundaries to represent the active aquifers, one block per layer per well, plus some additional blocks (which are neither acquifer blocks nor crossed by the wells).
The reservoir properties (i.e. the parameters that we will be interested in sampling from) are described by the pore volumes V of the blocks, ∈ {1, . . ., 124}, the transmissibilities T j between the interconnected blocks and j, and the perforation productivity coefficients J w for the well-block connections.We do not explain here the practical significance of such parameters and for more details on reservoir simulation we refer the reader to [1].Altogether the parameter space for this example is 338-dimensional.For the sake of clarity all (nonzero) T j are re-indexed with a single index as T p , p ∈ {1, . . ., 139}; and similarly J w are reindexed as J k , k ∈ {1 . . .75} and we denote by x ∈ R 338 the full vector of parameters, i.e. x = (V 1 , . . ., V 124 , T 1 , . . ., T 139 , J 1 , . . .J 75 ) T .There are 86 non-aquifer blocks in total, and we always assume an ordering of the parameters V such that the first 86 of them correspond to the non-aquifer blocks.In our Bayesian inverse problem for the parameters x the likelihood function is defined from the reservoir simulation output, and the prior is a Gaussian with covariance matrix C. The observed block pressure and the well bottom hole pressure (BHP) Yellow "v" stands for the block containing a perforation of a well.That is, the well goes through all the layers, but there is a hole for well-block communication only in correspondence of the yellow boxes.This figure does not show all the blocks -but only those perforated by the wells.In particular, it does not show the aquifer blocks located on the boundary.data are known for certain wells and time steps; we arrange such data into the vector d 0 .The likelihood L(d 0 |x), see equation (3.3) below, is defined using the simulator-modelled data d(x), the observed data d 0 , and the covariance matrix for data errors C d .The function d(x) is found by numerical solution of a system of ordinary differential equations, which we report in Appendix B, see (5.3) -(5.6); such a system describes the relation between the vector of reservoir properties x and the simulated pressures.The important thing for the time being is that such a system is high dimensional and the resulting posterior is analytically intractable. 3Finally, we seek to sample from the measure where the likelihood function is of the form In our numerical experiments we will always take the matrix C d to be diagonal, with the entries equal to either σ 2 BHP = 20 2 or σ 2 b = 9.We will give more details about this choice in Appendix B. The full parameterisation is further divided into three subcases denoted here as full-a, full-b, full-c, which have different min/max bounds for the parameters of interest or prior covariances.For the full-a case we define the minimum L i and maximum U i bounds of each parameter x i ∈ {V , T p , J k } as follows: let xi be the maximum likelihood value of the parameter x i , found approximately, by running an optimization algorithm on all the parameters; 4 we then take L i = 0.1x i , U i = 10x i , i = 1, ..., 338.Since the values of physical parameters x i may differ by several orders of magnitude, it makes sense to apply a transform to get a similar magnitude for all the parameters.Such a transform was done by function log 10 and a constant shift, mapping the parameters x i from the original range [L i , U i ] to [−1, 1].The prior covariance is taken as C full-a = 0.25 • I 338 .So, all the parameters in the transformed representation vary within the box [−1, 1] and have standard deviation 0.5.For the full-b case wider parameter bounds are taken: L i = 0.001x i , U i = 1000x i , i = 1, ..., 338.The parameters are transformed by log 10 function, and then mapped to the interval [−3, 3].The prior covariance is the same as in the full-a case, so all the parameters have standard deviation 0.5 in the transformed representation.Case full-c uses the same parameter bounds and the same transform as case full-b, but a wider prior covariance C full-c = 9 • C full-a , which means the prior standard deviation is 1.5 in the transformed representation.
• Lightweight parameterisation Here we consider a reduced, 21-dimensional, parameter space.Here we just fix the values of V 1 , . . ., V 86 (non-aquifer blocks), and we find the remaining V 87 , . . ., V 124 (aquifer blocks), T 1 , . . ., T 139 (all blocks), J 1 , . . .J 75 (all perforations), which are required by the simulator, using 21 new parameters.Such parameters essentially act as multipliers; namely, for each one of the seven layers n ∈ {A, . . .G} we introduce one pore volume multiplier for the aquifer blocks Ṽn , one transmissibility multiplier Tn , and one perforation productivity multiplier Jn .These parameters, collectively denoted by y ∈ R 21 , are those that we are interested in sampling from, by using the posterior measure where ρ(y) is a zero mean Gaussian with covariance matrix denoted by C 21 , described below.Because we are using the same simulator as for the full reservoir simulation, the likelihood function L is still the one defined in (3.3), hence necessarily we must have X(y) ∈ R 338 .To define the function X : R 21 → R 338 , we need to introduce some notation first.Denote by A n the number of aquifer blocks in layer n, P n the number of transmissibility coefficients in layer n, and K n the number of well perforations in this layer.Let V be the maximum likelihood value of the parameter V (similarly for Tp and Jk ), again found by running a maximum likelihood algorithm, and let the corresponding full vector denoted by x.The first 86 components of X(y) (corresponding to non-aquifer V ) are taken equal to V , = 1, . . .and so forth until the 7th column.The columns from 8 to 14 are built similarly, such that column n + 7 corresponds to layer n and has P n non-zero values equal to Tp (for appropriate indices p).The last seven columns are built in the same way, this time using the values Jk .
For the lightweight parameterisation the following minimum and maximum bounds were employed: [0.15,15] for all Ṽn , [0.07, 7] for all Tn ,and [0.11,11] for all Jn .As before, the physical parameters (multipliers) y i are mapped to the interval [−1, 1].The prior covariance C 21 , which acts in the transformed variables, was taken as essentially a diagonal matrix with the main diagonal equal to 0.25, however additional non-zero covariances equal to 0.1 were also specified between the transmissibility multiplier Tn and perforation productivity multiplier Jn for each layer n.A brief summary of the four discussed cases of the model parameterisation is presented in Table Table 1.Summary of the subcases for the reservoir simulation model.In physical representation the lower bounds are L i = l i b i , the upper bounds are U i = u i b i , where l i , u i are reported in the table, and b i are some base case parameter values (e.g. for all full parameterisations b i = xi ).

Numerics: sampling from measures supported on bounded domains
To compare efficiency of the algorithms we compute a normalised effective sample size (nESS), where the normalisation is by the number of samples N .Following [6], we define the Effective Sample Size ESS = N/τ int where N is the number of steps of the chain (after appropriate burn-in) and τ int is the integrated autocorrelation time, τ int := 1 + k γ(k), where γ(k) is the lag − k autocorrelation.Consistently, the normalised ESS, nESS, is just nESS := ESS/N .Notice that nESS can be bigger than one (when the samples are negatively correlated), and this is something that will appear in our simulations.As an estimator for τ int we will take the Bartlett window estimator (see for example [6,Section 6], and references therein) rather than the initial monotone sequence estimator (see again [6,Section 6]), as the former is more suited to include non-reversible chains.Since the nESS is itself a random quantity, we performed 10 runs of each case using different seeds, and our plots below show P10, P50, P90 percentiles of the nESS from these runs.4.1.Bounces vs Rejection.First we consider the performance of the two proposed methods for sampling from the box B. We illustrate these by comparing SOL-HMC-bounce and SOL-HMC-rej.
In Figure 2 we compare performance of SOL-HMC-bounce and SOL-HMC-rej for sampling from the 5D Rosenbrock target π Ros ; each one of the five parameters is taken to vary in the interval [−a, a], and Figure 2 shows how the performance varies as the size a of the box varies, a = 0.1, 0.2, . . ., 1.4.The target acceptance rate for both samplers was set to 0.9, and parameter i = 0.6 (defined in (2.14)).
As a "sanity test" the plots indicate that for the larger boxes (a ≥ 0.8) the two implementations SOL-HMC-bounce and SOL-HMC-rej are almost identical (in terms of nESS), which is natural as for large box sizes these two algorithms coincide.For small box sizes, the performance of the two samplers depends really on which coordinate is being sampled, so the performance of the two algorithms is substantially indistinguishable for this low dimensional problem.It is important to note the following practical drawback of SOL-HMC-rej (or indeed any other sampler which handles boundaries by the rejection mechanism) with respect to SOL-HMC-bounce: during the proposal step a trajectory may leave the box, and then return back inside the box.By construction of the algorithm such a trajectory is not rejected just because it escaped from the domain for a short while.The accept/reject decision is made only for the final point of the proposal trajectory, and thus every trajectory needs to be calculated till the end.However, if the trajectory is allowed to leave the box for the intermediate calculations, it may go to the extreme regions of the parameter space, where the simulator may suffer from severe numerical errors and abnormal behaviour.We illustrate this phenomena by comparing SOL-HMC-rej against HMC-bounce in Figure 3 for full-a in (A) and full-b in (B).(Here we think of HMC-bounce as sort of gold standard and for this reason we compare SOL-HMC-rej with HMC-bounce).We examine the ratio of nESS of SOL-HMC-rej and HMC-bounce and plot a histogram for the parameters.When the nESS ratio is bigger than one then SOL-HMC-rej is performing better than HMC-bounce.This is the case in (B) for full-b.However in (A) for full-a the boundary of B is encountered far more frequently, just because the size of the box for this target measure is smaller, see Table 1.Moreover a comparison of the histograms in Figure 3 (A) with the one in Figure 9 (A) shows better performance of SOL-HMC bounce with respect to SOL-HMC Rejections.From now on we consider SOL-HMC-bounce only.

4.2.
Comparison for 5D Rosenbrock.We consider the 5D Rosenbrock target π Ros where the minimum-to-maximum range for each one of the five parameters was taken as [−a, a], where a = 0.1, 0.2, . . ., 1.4.The plots in Figure 4 compare the performance of the HMCbounce, SOL-HMC-bounce, and Horowitz-bounce algorithms.The target acceptance rate is 0.9, and the parameter i = 0.6 for SOL-HMC-bounce and Horowitz-bounce.For this small dimensional problem we observe that SOL-HMC-bounce and Horowitz-bounce have similar nESS across the range of sizes for the box B. For smaller boxes B (e.g. a ≤ 0.5) all three algorithms have similar nESS.For larger box sizes we see for parameter x 1 an advantage in using SOL-HMC-bounce/Horowitz-bounce over the HMC-bounce however for x 2 there is a slight advantage to HMC-bounce.This corroborates the idea that in low dimension the advantage of introducing irreversibility in the sampler is hardly noticeable.

4.3.
Increasing parameter space and effectiveness of non-reversibility.We now consider our more realistic targets and increase the parameter space to 21 and then to 338.We now clearly see the advantage of the non-reversible algorithms SOL-HMC-bounce and Horowitz-bounce over HMC-bounce.
Figure 5 reports the nESS for the lightweight parameterisation of the reservoir simulation problem for the following four cases: HMC-bounce with an acceptance rate of 0.82 (target 0.8), SOL-HMC-bounce with acceptance rate of 0.77 (target 0.9), SOL-HMC-bounce with an acceptance rate of 0.69 (target 0.8) and Horowitz-bounce with an acceptance rate of 0.68 (target 0.8).All SOL-HMC-bounce and Horowitz-bounce algorithms took the parameter i = 0.5 and here we give results from a single MCMC run in each case.The plot clearly shows that the non-reversible algorithms outperform HMC for the majority of the parameters.We also observe the variability due to acceptance rate: for SOL-HMC-bounce a better nESS is achieved for the higher acceptance rate.As we further increase the dimension and complexity the advantage of the non-reversible algorithm becomes further apparent.In Figure 6 we compare for full-a SOL-HMC-bounce and HMC-bounce and observe a clear improved nESS for SOL-HMC-bounce across the whole parameter space.Finally we compare SOL-HMC-bounce and Horowitz against the benchmark of HMCbounce by examining the ratio of nESS.Recall that when the ratio is bigger than one then SOL-HMC-bounce (or Horowitz) has a larger nESS than HMC.We consider the targets fulla, full-b and full-c.In Figure 7 we compare for full-a SOL-HMC-bounce against Horowitzbounce.First note that in both bases the nESS ratio is > 1 for most parameters showing a clear improvement in the non-reversible algorithms over HMC.To aid comparison between SOL-HMC-bounce against Horowitz-bounce we plot on (A) and (B) a fit of the histogram from (A), this is the black dotted line.We see that the nESS for SOL-HMC-bounce over the parameters is larger than that for Horowitz-bounce and that there is an improvement using SOL-HMC-bounce.Here we took i = 0.5.Figure 8   Finally, in Figure 9, we examine SOL-HMC-bounce for full-a (A), full-b (B) and full-c (C) for the same value of i = 0.4.We see a clear improvement of the non-reversible SOL-HMCbounce over HMC in each case.We compare here to the SOL-HMC-bounce for full-b for the same value of i = 0.4 in (B).We observe a similar improvement for SOL-HMC-bounce over HMC in both cases.

Conclusion
We have investigated two different ways to deal with sampling measures on a bounded box B: rejection and bounces.This is crucial in many practical applications, for example to respect physical laws (such as porosity for reservoir modelling or pixel values in image reconstruction).We have explained and demonstrated why, for complex problems involving the use of a proxy, reflection algorithms should be preferred to rejection strategies.We have furthermore shown that when sampling from complex realistic target measures, such as those that arise in reservoir simulation, non-reversible algorithms such as SOL-HMC and Horowitz outperform standard reversible algorithms such as HMC.In addition, we see that as the problem size grows SOL-HMC is superior to Horowitz having larger nESS.
Q j = T j (P − P j − ρ g h j ), (5.3)where ρ is the known liquid density, h j is the known depth difference between the block centers, and g is the acceleration due to gravity.
The inflow q w into the perforation of well w in block is proportional to the difference of the bottom-hole pressure (BHP) and the block pressure: q w = J w (P − P BHP w − ρ g h w ), (5.4) where h w is the depth difference between the block center and the BHP reference depth.The total inflow into well w is obtained by summing up contributions related to this well; that is, q w = q w .(5.5) Finally, the volumetric inflows and outflows for block are balanced, with the excessive/deficient fluid volume leading to the block pressure change via the following compressibility equation: where t denotes time, and the first (second, respectively) summation on the right hand side is taken over all the blocks j connected to the block (all the wells w perforated in block , respectively).The compressibility c of the block is supposed to be known.The simulated reservoir time spans 12 years.

Figure 1 .
Figure 1.Perforations of the 40 wells (columns) in the seven layers (rows).Yellow "v" stands for the block containing a perforation of a well.That is, the well goes through all the layers, but there is a hole for well-block communication only in correspondence of the yellow boxes.This figure does not show all the blocks -but only those perforated by the wells.In particular, it does not show the aquifer blocks located on the boundary.

Figure 2 .
Figure 2. Normalised ESS for SOL-HMC-bounce (blue) and SOL-HMC-rej (black ), for different sizes of the box bounding the parameter space (X-axis).The two plots correspond to coordinates x 1 , x 3 only.The other three coordinates have nESS plots similar to x 3 .

Figure 3 .
Figure 3. Ratio of nESS (SOL-HMC-rej divided by HMC-bounce).Parameterisation full-a (A) is shown in green and full-b (B) in blue.

Figure 4 .
Figure 4. Normalised ESS for SOL-HMC-bounce (blue), HMC (black ), and Horowitz (orange), for different sizes of the box bounding the parameter space (X-axis).The two plots correspond to coordinates x 1 , x 2 only.The other three coordinates show a similar picture.

Figure 5 .
Figure 5. Normalised ESS (Y axis) for the reservoir simulation MCMC, lightweight parameterisation.X axis shows the 21 parameters.In the legend, the real acceptance rates are indicated.

Figure 6 .
Figure 6.Normalised ESS (Y axis) for the reservoir simulation MCMC, fulla parameterisation.X axis shows the 338 parameters.The samplers are HMC and SOL-HMC with i = 0.4.

Figure 7 .
Figure 7. Ratio of nESS for Horowitz-bounce by nESS for SOL-HMCbounce.The target measure here is the 338-dimensional full-a.Parameter i = 0.5.

Figure 9 .
Figure 9. Ratio of nESS for SOL-HMC-bounce for targets full-a (A), full-b (B) and full-c (C) and in each case i = 0.4. 1.