  We present a highly efficient proximal Markov chain Monte Carlo methodology
to perform Bayesian computation in imaging problems. Similarly to previous
proximal Monte Carlo approaches, the proposed method is derived from an
approximation of the Langevin diffusion. However, instead of the conventional
Euler-Maruyama approximation that underpins existing proximal Monte Carlo
methods, here we use a state-of-the-art orthogonal Runge-Kutta-Chebyshev
stochastic approximation that combines several gradient evaluations to
significantly accelerate its convergence speed, similarly to accelerated
gradient optimisation methods. The proposed methodology is demonstrated via a
range of numerical experiments, including non-blind image deconvolution,
hyperspectral unmixing, and tomographic reconstruction, with total-variation
and $\ell_1$-type priors. Comparisons with Euler-type proximal Monte Carlo
methods confirm that the Markov chains generated with our method exhibit
significantly faster convergence speeds, achieve larger effective sample sizes,
and produce lower mean square estimation errors at equal computational budget.
