{'fraction_non_alphanumeric': 0.04838812233738157, 'fraction_numerical': 0.01611877662618427, 'mean_word_length': 4.456201214223764, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The split Gibbs sampler revisited: improvements to its algorithmic structure and augmented target distribution

This paper proposes a new accelerated proximal Markov chain Monte Carlo (MCMC) methodology to perform Bayesian computation efficiently in imaging inverse problems. The proposed methodology is derived from the Langevin diffusion process and stems from tightly integrating two state-of-the-art proximal Langevin MCMC samplers, SK-ROCK and split Gibbs sampling (SGS), which employ distinctively different strategies to improve convergence speed. More precisely, we show how to integrate, at the level of the Langevin diffusion process, the proximal SK-ROCK sampler which is based on a stochastic Runge-Kutta-Chebyshev approximation of the diffusion, with the model augmentation and relaxation strategy that SGS exploits to speed up Bayesian computation at the expense of asymptotic bias. This leads to a new and faster proximal SK-ROCK sampler that combines the accelerated quality of the original SK-ROCK sampler with the computational benefits of augmentation and relaxation. Moreover, rather than viewing the augmented and relaxed model as an approximation of the target model, positioning relaxation in a bias-variance trade-off, we propose to regard the augmented and relaxed model as a generalisation of the target model. This then allows us to carefully calibrate the amount of relaxation in order to simultaneously improve the accuracy of the model (as measured by the model evidence) and the sampler's convergence speed. To achieve this, we derive an empirical Bayesian method to automatically estimate the optimal amount of relaxation by maximum marginal likelihood estimation. The proposed methodology is demonstrated with a range of numerical experiments related to image deblurring and inpainting, as well as with comparisons with alternative approaches from the state of the art.


Introduction
The problem of estimating an unknown image from noisy and/or incomplete data is central to imaging sciences [14,7]. Canonical examples include, for example, noise removal [29], image inpainting [50], image deblurring [18], medical imaging [27,3,32], astronomical imaging [12,13]. Estimation by direct inversion of the foward model relating the unknown image to the data is not usually possible, inasmuch the inverse problem is often severely ill-conditioned or ill-posed. The literature describes a range of mathematical frameworks to incorporate regularisation and formulate well-posed solutions (see, e.g., [26,43,7]).
We consider Bayesian statistical solutions to imaging inverse problems [26]. In this case, the solution is a probability distribution characterising our knowledge about the value of the unknown image of interest given the observed data. This distribution can then be used to derive image estimators, calibrate unknown model parameters, perform uncertainty quantification analyses, or Bayesian model selection (see, e.g., [22,45,37,11]).
There are three main strategies to perform Bayesian computation in imaging inverse problems. One strategy relies on optimisation methods, which typically scale efficiently to large models and offer detailed convergence guarantees, but can only support maximum-a-posteriori (MAP) estimation and a very limited range of other inferences [14,13,39]. Another strategy is to iteratively The remainder of the paper is organised as follows: In Section 2 we introduce the models considered in this work, and recall the state-of-the-art MCMC methods to sample from them. In Section 3 we revisit the augmented model and show empirically that there is a subset of hyperparemeter values that enhances this model, while we also present a method to automatically compute an optimal choice of these hyperparameters. In Section 4 we establish a formal connection between MYULA and SGS which allows us to propose two novel and more efficient sampling algorithms for the augmented model. Section 5 illustrates the proposed methodologies with two experiments related to image deblurring and image inpainting, where we report detailed comparisons with state-of-the-art algorithms. Conclusions and perspectives for future work are reported in Section 6. In this work, we pay special attention to problems where the estimation of x given y is ill-posed or ill-conditioned 1 . We address this difficulty by considering the Bayesian framework, where we regularize the estimation problem by specifying a prior distribution, given for any x and θ ∈ (0, +∞) p by p(x|θ) ∝ e −θ T g(x) , for some vector of statistics g : R d → R p , where θ parametrises this prior distribution and controls the level of imposed regularity. Using Bayes' theorem, we derive the posterior distribution p(x|y, θ) = p(y|x)p(x|θ) R d p(y|x)p(x|θ)dx ∝ exp −fy(x) − θ T g(x) , (2.1) which models the knowledge we have about x given the observed data y. We focus on Bayesian computational methodology for log-concave models of the form (2.1), where fy and g satisfy the following conditions: 1. fy : R d → R is convex and Lipschitz continuously differentiable with constant L f 2 .
Models with these characteristics are widely adopted in the imaging community due to the variety of currently available Bayesian optimization tools that exploit the convexity properties of p(x|y, θ) [39,33], such as the computation of the MAP estimate, which can be formulated as a convex optimization problem [35] given bŷ xMAP = argmax x p(x|y, θ) = argmin x fy(x) + θ T g(x), (2.2) and can be solved by using state-of-the-art convex optimization algorithms [14]. However, there are other more complicated Bayesian analyses beyond point estimation, such as model calibration, Bayesian model selection and hypothesis testing [41], which cannot be addressed by using optimization algorithms and typically require the calculation of probabilities and expectations w.r.t p(x|y, θ). This is challenging in imaging problems since it requires calculating intractable integrals on R d . In this case, the application of MCMC methods [40,6] and, in particular, proximal MCMC methods [34,22,37,46,48], specialised for non-smooth log-concave distributions, is the preferred approach.

Sampling via the Langevin diffusion
We consider proximal MCMC methods derived from the discretization of the overdamped Langevin diffusion process, which we discuss below. Assume that we are interested in calculating probabilities w.r.t. a smooth distribution with density π(x), and consider the stochastic differential equation where Wt is a d-dimensional Brownian motion. Under mild assumptions on π(x), this SDE has a unique strong solution and admits π(x) as its unique invariant density [42, Theorem 2.1]. However, in imaging applications it is usually not possible to solve (2.3) exactly, so a numerical approximation needs to be employed. Most works consider the Euler-Maruyama (EM) scheme given by where δ > 0 is a given step-size and (Zn+1) n≥0 is an i.i.d. sequence of d-dimensional standard Gaussian random vectors. This recursion is known as the unadjusted Langevin algorithm (ULA) and has been shown to be a highly efficient method for highdimensional Bayesian inference when π(x) is log-concave and smooth with ∇ log π(x) L-Lipschitz continuous and δ < 1/L [21,20]. However, in many imaging models, π(x) is not smooth and hence appropriate adjustments need to be made to ULA.
Applying the ULA scheme to the smooth posterior approximation p λ (x|y, θ) leads to the recursion which is known as the Moreau-Yosida unadjusted Langevin algorithm (MYULA) [22]. The main benefit of the MYULA is that now since p λ (x|y, θ) is smooth and preserves log-concavity, the results from [21,20] apply hence providing an efficient method for imaging applications. One of the main computational bottlenecks of ULA and MYULA is the fact that, in order to converge, one needs to choose δ ≤ 1/L where L is the Lipschitz constant of ∇ log p λ (x|y, θ), given by L = L f + p/λ (see [22,Theorem 3.2]). This step-size restriction is not problematic for moderate values of L, for example in denoising problems. However, in problems of the form y = Ax + ξ with ξ ∼ N (0, σ 2 Im) and σ > 0 where the forward operator A is very poorly conditioned, the MYULA sampler will have poor mixing properties, particularly on the subspace of R d where x has high uncertainty (this is related to the slow components of the Markov chain). A similar situation occurs when the problem requires a high level of accuracy in the approximated prior term g λ i (x) (i.e., λ is very small), for example in problems where one has to enforce domain constraints on the solution space.

SK-ROCK
A natural way of overcoming the step-size limitation of the MYULA is to adopt a discretization scheme for the Langevin Diffusion (2.3) with better numerical stability properties. In particular, an explicit stochastic Runge-Kutta-Chebyshev discretization of the Langevin SDE was proposed in [2] called SK-ROCK, and a proximal variant suitable for computational imaging problems was recently proposed in [37]. This highly advanced Runge-Kutta stochastic integration scheme extends the deterministic Chebyshev method [1] to sampling processes and has been shown to require fewer gradient evaluations than ULA or MYULA to reach a given level of accuracy. In particular, it has been shown in [37, Proposition 3.1] that, in order to be ε close to a multivariate Gaussian target distribution (in the 2-Wasserstein distance), SK-ROCK requires O( √ κ) gradient evaluations (where κ is the condition number of ∇ log π(x)) instead of O(κ) required by the ULA, similar to accelerated optimization methods. This is possible due to the fact that SK-ROCK uses s gradient evaluations per iteration at carefully chosen extrapolated points, which allows using a much larger step-size than the MYULA and thus having a faster decorrelation of the Markov chain in a stable manner. This is one of the great advantages of SK-ROCK in very ill-posed and ill-conditioned models, compared to MYULA, accelerating its convergence while maintaining its stability, thus improving the quality of the generated samples. The SK-ROCK scheme is shown in Algorithm 1.

Sampling via augmentation: Split Gibbs sampler (SGS)
A separate line of research seeks to address the limitation of MYULA by introducing an auxiliary variable z ∈ R d and operating on the augmented state-space (x, z). This allows to relax the original model (2.1) and instead uses the following augmented posterior
where ρ 2 controls the correlation between the variable of interest x and the auxiliary variable z, and fy, g are the same as in Section 2.1. If we now consider the marginal posterior distribution it is possible to show that it converges in total variation to the original posterior p(x|y, θ) as ρ 2 → 0. This approach was first introduced several decades ago as a way of calculating maximum likelihood estimates from incomplete data [17], and as an efficient method for sampling from posterior distributions [44] (see [23] for a review of these techniques). In the current literature, this model was revisited by [38] in the context of consensus Monte Carlo in distributed settings and applied to imaging inverse problems in [46], where its similarities to the algorithmic structure of the Alternating Direction Method of Multipliers (ADMM) optimization algorithm [10] were also discussed.
From a computational point of view, as in the case of MYULA, because g is not differentiable one needs to approximate p(x, z|y, θ, ρ 2 ) by To sample (2.9), [46,48] proposed a Gibbs-like splitting strategy scheme, applied on the following conditional distributions This method is known as the split Gibbs sampler (SGS). See Algorithm 2.
In the case where the likelihood is Gaussian one can exactly sample from (2.10) [25] (for a review and comparison of existing Gaussian sampling approaches, see [47]). A main benefit of this splitting approach is that the step-size one needs to set for the proximal MCMC method used for sampling (2.11) will be independent of the Lipschitz constant associated with the likelihood distribution, and will only depend on the parameters λ and ρ 2 (i.e., the step-size now depends on the Lipschitz constant indirectly via λ and ρ 2 , there is still some dependence w.r.t. the Lipschitz constant through the bias incurred). This can lead to faster sampling algorithms compared to MYULA for suitably chosen values of the parameter ρ 2 [46], albeit for a biased posterior distribution.

Enhancing Bayesian imaging models by smoothing
As discussed previously, the augmented model (2.7) was originally proposed as a relaxation of (2.1) that allows for a faster exploration of the target distribution, at the expense of some additional bias when compared to the original model. One then might think that ρ 2 = 0 represents the best model for inference (at the expense of higher computing cost). However, we have found empirically that this is not the case.
As an illustration, Figure 1(a) shows the estimation mean-squared error (MSE) for a Bayesian image debluring problem (the details of this experiment will be explained in Section 5.1). The error is computed w.r.t. the posterior mean, as estimated by an adaptation of the SK-ROCK method to target (2.9) (see Section 4.3 for details), using a value of θ = 4.4 × 10 −2 estimated by [45,Algorithm 1], and by using different values for ρ 2 . Recalling that increasing ρ 2 improves convergence speed, one can clearly identify a regime of small values of ρ 2 for which convergence speed improves without a deterioration in estimation accuracy (in fact, there is a mild improvement). Beyond this range, the estimation MSE deteriorates dramatically. This suggests the need for a method to automatically set the value of ρ 2 .
We propose an empirical Bayesian method to estimate optimal values for θ and ρ 2 directly from y by maximum marginal likelihood estimation (MMLE) where Θ ⊂ (0, +∞) p and Ω ⊂ (0, +∞) are compact convex sets. To solve (3.1) we modify the stochastic approximation proximal gradient (SAPG) algorithm of [45]. Notice that by maximising the model evidence, (3.1) seeks to select the best model to perform inference within the class of posterior distributions parametrised by θ ∈ Θ, ρ 2 ∈ Ω [49].

Computing the optimal values for θ and ρ 2
We adopt the approach of [45] to solve (3.1) and estimate optimal values for θ and ρ 2 in (2.7). The method [45] was proposed for models of the form (2.1), so we will now adapt it to the augmented model (2.7). We are interested in estimating the parameters θ ∈ Θ, ρ 2 ∈ Ω by MMLE (3.1), where p(y|θ, ρ 2 ) is defined for (2.7) as and If we had access to the gradients ∇ ρ 2 log p(y|θ, ρ 2 ) and ∇ θ log p(y|θ, ρ 2 ), then we could construct an interative algorithm that converges to the solution of (3.1) by using the projected gradient algorithm [30] where ΠΩ and ΠΘ are the projection onto Ω and Θ respectively, and (γn) n∈N is a sequence of non-increasing step-sizes. However, due to the complexity of the model, ∇ ρ 2 log p(y|θ, ρ 2 ) and ∇ θ log p(y|θ, ρ 2 ) are intractable. As shown in [45], one can construct carefully designed stochastic estimates of these gradients that satisfy the conditions for the solution to converge to (3.1). To build these stochastic estimators, we are going to express the gradients as expectations by applying the Fisher's identity [19,Proposition D.4] which we can then approximate using MCMC. We will see that in fact one MCMC sample will suffice to obtain an estimate of the gradient accurate enough to converge asymptotically to (3.1). More precisely, we have that As p(x, z, y|θ, ρ 2 ) = p(y|x)p(x, z|θ, ρ 2 ) = p(y|x)p(x|z, ρ 2 )p(z|θ), we have and similarly, replacing (3.2) in p(z|θ) gives Because of the complexity of the model, A θ,ρ 2 (y) and B θ,ρ 2 (y) are not available analytically and need to be approximated by MCMC computation (e.g., by using the methods we develop in Section 4). With respect to C θ,ρ 2 (y), and more precisely, the integral between brackets, we can follow a similar procedure as in [45, Section 3.2.1]. In particular, if we consider the case where each gi(z) is αi positively homogeneous 4 , which is the case for many regularizers such as ℓ1, ℓ2 or TV, we have that (See [45] for more details and [45, Section 3.2] for the case of inhomogeneous regularizers).
Following on from this, and by using Monte Carlo approximations of A θ,ρ 2 (y) and B θ,ρ 2 (y), we construct an SAPG algorithm [24,45] to solve (3.1) and produce optimal estimates of θ and ρ 2 . This method is presented in Algorithm 3. We refer the reader to [16,15] for details about the convergence properties of this kind of SAPG algorithm. To illustrate Algorithm 3 in action, Figure 1 shows the value of ρ 2 estimated by the algorithm for the image deblurring problem. Observe that the MMLE estimate is close to the value that produces the best estimation MSE in this case. This is in agreement with the results reported in [45] for other problems.

Noisy MYULA
We now proceed to show that the SGS algorithm 2 can be viewed as a noisy version of MYULA. This link will be crucial in allowing us to write it as a noisy discretisation of an SDE, which will help us to propose more efficient MCMC methods for sampling (2.7).
First, note that the marginal of z computed from (2.9) can be written as follows is α positively homogeneous if, for any x ∈ R d and t > 0, g(tx) = t α g(x).
As p(x, y|z, ρ 2 ) = p(y|x)p(x|z, ρ 2 ), we have ∇z log p(x, y|z, ρ 2 ) = E x|y,z,ρ 2 ∇z log p(x|z, ρ 2 ) Using this expression in (4.1) we obtain We are now ready to explicitly establish the connection to SGS. SGS stems from dealing with the presence of the expectation in this algorithm by replace it by a Monte Carlo empirical average, i.e., More precisely, to recover SGS we take N = 1 and substitute in (4.2) to obtain Since X (1) is an exact sample from p(x|y, Zn, ρ 2 ), (4.4) corresponds to the fourth line of Algorithm 2. This establishes that SGS is equivalent to a noisy version of MYULA that relies on one sample from p(x|y, z, ρ 2 ) to compute a stochastic estimate of the gradient ∇z log p(y|z, ρ 2 ) via (4.3). Using multiple samples from p(x|y, Zn; ρ) would improve the estimation of the expectation (4.3) and hence the behaviour of the algorithm. Alternatively, in the experiments considered in this paper p(x|y, Zn; ρ) is Gaussian, and hence this expectation can be calculated exactly. This is exploited in the MCMC methods proposed below.

Latent space MYULA
We established that SGS is equivalent to MYULA targeting the marginal of z with an inexact (i.e., stochastic) estimate of the gradient. Replacing this stochastic estimate with its exact value in Algorithm 2 produces the following recursion where ζi+1 ∼ N (0, I d ).
We now discuss how to use samples {Zi} m i≥1 to compute expectations w.r.t. the marginal of interest x|y, θ, ρ 2 . More precisely, consider the computation of an expectation E x|y,θ,ρ 2 [h(x)] for some function h w.r.t. the posterior distribution p λ (x|y, θ, ρ 2 ) defined in (2.8) by using (4.5). Formally, Using the fact that p(x, z|y, θ, ρ 2 ) = p(x|y, z, ρ 2 )p(z|θ) we have that In cases where E x|y,θ,ρ 2 [h(x)] is available analytically, we suggest using a Rao-Blackwellised estimator of the form [40] The computation of E x|y,Z i ,ρ 2 [h(x)] can be done as a postprocessing step, or alternatively within the iterations of the sampler. If E x|y,θ,ρ 2 [h(x)] is not available analytically, we would draw samples from the conditional x|y, Zi, θ, ρ 2 and apply a standard Monte Carlo estimator.
We are now ready to present our first new MCMC method, summarised in Algorithm 4 below. We henceforth refer to this method as latent space MYULA (ls-MYULA), since it corresponds to MYULA applied to the marginal of the latent variable z.

Latent space SK-ROCK
In the same way that an exact MYULA discretization is more beneficial than the stochastic MYULA discretization used in SGS, we can further improve results by using an exact SK-ROCK discretization which, as we described in Section 2.2.2, has many important advantages compared to MYULA. In particular, we present this method in Algorithm 5, and we will refer to it as latent space SK-ROCK (ls-SK-ROCK). The main difference between this algorithm and Algorithm 1 is that the conditional expectation E x|y,Z j ,ρ 2 [x] is computed on each internal stage s.

Implementation guidelines
Setting λ As the priors of the experiments performed in this work are non-differentiable, we will use the Moreau-Yosida envelope defined in We chose λ = L −1 f in our numerical experiments, however, we have found numerically that values of λ = 5L −1 f or λ = 10L −1 f lead to faster convergence at the cost of additional bias.
Setting γ i , γ ′ i and n With respect to Algorithm 3, it is suggested in [45] to set γi = C0i −p and γ ′ i = C ′ 0 i −p where p ∈ [0.6, 0.9] (in the experiments performed in this paper, we have set p = 0.8), C0 and C ′ 0 starting with (θ0d) −1 and (ρ 2 0 d) −1 respectively, and then readjusting it if necessary. With respect to n, we have followed the recommendation in [45] and used a single sample (i.e., n = 1) on each iteration (we did not observe significant difference for larger values of n).

Setting an stopping criteria
It is recommended to supervise the evolution of |θm+1 − θm|/θm and |ρ 2 m+1 − ρ 2 m |/ρ 2 m in the execution of Algorithm 3 until they reach a tolerance level β to stop the algorithm execution. In our imaging experiments, we set β = 10 −4 but we have observed that β = 10 −3 is often enough to reach an acceptable estimate of the hyper-parameters in small computational times.

Other implementation considerations
In the implementation of the SAPG method, it is important to update the step-size of the MCMC method to sample Xi and Zi within each iteration of the SAPG scheme, as the maximum step-size depends on the value of ρ 2 i . Regarding the Lipschitz constant L one needs to compute the step-size for MYULA and SK-ROCK algorithms, the model (2.5) has L = λ −1 + L f . With respect to the augmented model (2.8), the Lipschitz constant is La = λ −1 + (ρ 2 + L −1 f ) −1 which we use to implement SGS, ls-MYULA and ls-SK-ROCK. Therefore, we set the step-size of the MCMC methods to 1/L for MYULA, to 1/La for SGS and ls-MYULA,, and to δ max

Numerical experiments
We now illustrate the improvement that can be obtained by sampling the augmented model (2.7) using Algorithms 4 and 5 together with an optimal estimate of ρ 2 using Algorithm 3. To evaluate the performance of the methods in a variety of situations, we perform two imaging experiments related to image deblurring (whose model is strongly log-concave) and image inpainting (whose model is weakly log-concave). We implement these algorithms as described in the implementation guidelines (see Section 4.4). For a fair comparison the results we show have been plotted as a function of the number of gradient evaluations, i.e., the number of times ∇ log p λ (x|, y, θ) and ∇z log p λ (z|y, θ, ρ 2 ) are computed in our algorithms 6 . The plots we show include the evolution of the MCMC samples in the burn-in stage using the scalar statistic log p(Xn|y, θ) for MYULA and SK-ROCK, and log p(X grad n |y, θ) for SGS, ls-MYULA and ls-SKROCK. We have also plotted the progression of the mean-squared error (MSE) between the posterior mean and the true image, when all the algorithms have reached stationarity, including the MAP estimate defined in (2.2) and computed using a highly efficient optimization algorithm called SALSA [5,4] for the image deblurring and image inpainting experiments.
We also provide pixel-wise standard deviation plots as a way of quantifying the uncertainty in the delivered solution. We have also computed standard deviation plots performing a downsampling by averaging the samples by a factor of 2 × j where j = {1, 2, 4}, which allows us to observe the uncertainty in image structures at different scales. Finally, we also show autocorrelation plots of the slowest component of the samples produced by each of the methods 7 , applying a 1-in-s thinning to the MYULA, SGS and ls-MYULA chains to equal the number of gradient evaluations between the mentioned methods (one gradient evaluations per iteration) and SK-ROCK/ls-SKROCK methods (s gradient evaluations per iteration). With this quantity, we have also compute effective sample sizes (ESS) 8 of the five algorithms discussed in this paper, where the sum is truncated at lag k when the lag-k autocorrelation reaches a value less than 0.05.
For completeness, in Table 7 we have also provided computing times of all the experiments. These results have been obtained on an Intel core i5-8350U@1.70GHz workstation running MATLAB R2018a.

Image deblurring
To examine the performance of the MCMC methods in different scenarios, we consider a deblurring problem with three test images (cameraman, boat and man) of size d = 256 × 256 pixels. An uniform blur operator H of size 5 × 5 is applied to the true image x ∈ R d and then additive Gaussian noise η is added with a sigma-to-noise level of 40dB, to produce an observation y ∈ R d related to the true image by y = Hx + η. As H is nearly singular, the problem becomes ill-conditioned. So, to promote regularity, we have used the TV norm as a prior which leads to the following posterior distributions where fy(x) = y − Hx 2 /2σ 2 and g(x) = TV(x).  We begin estimating optimal values for θ and ρ 2 for the given models implementing Algorithm 3 setting γi = γ ′ i = 10 × i −0.8 /d, θ0 = 0.04, ρ 2 0 = L −1 f = σ 2 and X0 = Z0 = H ⊺ y. The corresponding results for the parameters estimation are given in Table 1, together with the Lipschitz constants L and La required to sample (5.1) and (5.2) respectively. We then generate 5 × 10 6 samples using MYULA and 5 × 10 6 /s samples using SK-ROCK (with s = 15) from (5.1), and 5 × 10 6 samples using SGS and ls-MYULA and 5 × 10 6 /s samples using ls-SK-ROCK (with s = 15) from (5.2). The results of these three experiments are plotted in Figure 3. In particular, we note from the evolution of the MSE (when the chains have reached the typical set of the target distributions) that ls-MYULA and ls-SK-ROCK outperform SGS, as we are using an exact MYULA implementation rather than a noisy one, as shown in Section 4. We have also reported the step-size used by each method in Table 2.
We also plot the autocorrelation function of the slowest component from the chains of the MCMC algorithms, this is shown in Figure 3(c),(f),(i) and, as can be seen, ls-SK-ROCK presents the fastest decay. In addition, Table 3 shows the effective sample sizes (ESS) associated with these autocorrelation plots, and one can notice that ls-SK-ROCK reaches the largest ESS. In addition, we also illustrated in Figure 4 the minimum mean-square estimator (MMSE) of all the MCMC methods for all three deblurring experiments.
Finally, Figure 5 shows the marginal posterior standard deviation of the cameraman deblurring experiment at different scales. One can notice that edges show higher uncertainty, which is expected due to the nature of the forward operator. As can be seen, ls-MYULA and, in particular, ls-SK-ROCK outperform SGS in terms of delivering accurate (i.e., less noisy) estimates, showing the benefit of using these algorithms to sample in a more efficient way the augmented posterior distribution.

Image inpainting
To conclude, we perform an image inpainting experiment. It consists of randomly selecting 60% of the image pixels x ∈ R d to form the observation vector y ∈ R m (m < d) and then adding Gaussian noise with a SNR level of 40dB (the observation images y are illustrated in Figure 6). To test the MCMC methods in different regimes, we will use the same three test images given in Section 5.1 (cameraman, boat and man) and illustrated in Figure 2(a)-(c). For this experiment, we consider the following models where fy(x) = y − Ax 2 /2σ 2 , A ∈ R m×d is a rectangular matrix obtained by taking a random subset of rows from the identity matrix in dimension d, and g(x) = TV(x).
With the generated samples, we proceed to plot the results of these experiments in Figure 7. We first notice the acceleration one can get from ls-SK-ROCK from the convergence to equilibrium of the MCMC samples in the burn-in stage represented by the evolution of the scalar estimate log p(Xn|y, θ). Then, we illustrate the evolution of the mean-squared error (MSE) between the  mean of the samples and the true image x in stationarity and, as can be seen, ls-SK-ROCK is computationally efficient in being the fastest method to reach the MMSE in all three experiments, even outperforming the MAP estimate in terms of accuracy in all three experiments; moreover, the improvement of ls-MYULA over SGS , in terms of accuracy is evident similar to our previous results. We also plot the autocorrelation function of the pixel values for the slowest component in Figure 7(c),(f),(i) and, as can be seen, the ACF of the ls-SK-ROCK samples decays faster than all the other MCMC methods. In addition, we also illustrate in Figure 8 the MMSE of all the MCMC methods for all three inpainting experiments and, as in previous numerical results, we can see in Figure  7(b),(e),(h) that ls-SK-ROCK is the fastest method in compute this estimate.
Finally, Figure 9 presents uncertainty quantification plots by showing pixel-wise standard deviation estimates for the cameraman inpainting experiment. In this case, the uncertainty is concentrated on the unobserved pixels, which is expected given the nature of the inpainting problem. One can notice that ls-MYULA and ls-SK-ROCK is slightly less noisy than SGS, showing the good performance of these algorithms in sampling the augmented posterior distribution.

Conclusions
We presented a strategy to combine MYULA or proximal SK-ROCK with augmentation and relaxation in the manner of SGS. This was achieved by first establishing that SGS is equivalent to a noisy ULA scheme applied to the marginal distribution of the     latent variable z in an augmented Bayesian model x, z|y, θ, ρ 2 . This then naturally led to two new samplers that apply MYULA and SK-ROCK to the latent marginal distribution z|y, θ, ρ 2 . Probabilities and expectations w.r.t. the primal marginal x|y, θ, ρ 2 are then straightforwardly computed by using a Rao-Blackwellised Monte Carlo estimator. Moreover, we also observed empirically that there is a range of values for ρ 2 for which convergence speed and model quality improve (in the sense of the model evidence). Increasing ρ 2 beyond this range leads to improvements in convergence speed at the expense of significant estimation bias. We therefore proposed to adopt an empirical Bayesian approach and set ρ 2 , together with the regularisation parameter θ, by maximum marginal likelihood estimation from y. This was achieved by using an SAPG scheme that convergences in very few iterations. We illustrated the benefits from adopting the proposed methodology with two experiments, image deblurring and image inpainting. The results showed that the new proximal SK-ROCK algorithm that benefits from augmentation and relaxation outperforms the other methods from the state of the art in terms of computational efficiency. Future work will focus on extending the proposed approach to plug-and-play priors encoded by neural network denoisers [28], and on establishing non-asymptotic convergence results for the methods.