Diffusion models (DMs) have shown great potential in modeling complex image distributions, making them promising priors for solving Bayesian inverse problems. However, existing methods often rely on approximations to make DMs generic across different inverse problems, leading to inaccurate posterior distributions. To address this, the authors propose a Markov chain Monte Carlo (MCMC) algorithm called Plug-and-Play Diffusion Models (PnP-DM), which reduces posterior sampling to sampling the posterior of a Gaussian denoising problem. The key innovation is leveraging a general DM formulation to rigorously solve the denoising problem with various state-of-the-art DMs. The method is demonstrated on six inverse problems, including a real-world black hole imaging problem, showing superior performance in terms of accuracy and posterior estimation compared to existing DM-based methods. The PnP-DM algorithm is based on the Split Gibbs Sampler, which alternates between likelihood and prior steps, with an annealing schedule to improve mixing and avoid local minima. Theoretical analysis provides insights into the non-asymptotic behavior of PnP-DM, showing that it converges to the target posterior at a rate of \(O(1/K)\). Experimental results on synthetic and real-world data validate the effectiveness of PnP-DM, highlighting its ability to handle both linear and nonlinear inverse problems, including highly ill-posed cases.Diffusion models (DMs) have shown great potential in modeling complex image distributions, making them promising priors for solving Bayesian inverse problems. However, existing methods often rely on approximations to make DMs generic across different inverse problems, leading to inaccurate posterior distributions. To address this, the authors propose a Markov chain Monte Carlo (MCMC) algorithm called Plug-and-Play Diffusion Models (PnP-DM), which reduces posterior sampling to sampling the posterior of a Gaussian denoising problem. The key innovation is leveraging a general DM formulation to rigorously solve the denoising problem with various state-of-the-art DMs. The method is demonstrated on six inverse problems, including a real-world black hole imaging problem, showing superior performance in terms of accuracy and posterior estimation compared to existing DM-based methods. The PnP-DM algorithm is based on the Split Gibbs Sampler, which alternates between likelihood and prior steps, with an annealing schedule to improve mixing and avoid local minima. Theoretical analysis provides insights into the non-asymptotic behavior of PnP-DM, showing that it converges to the target posterior at a rate of \(O(1/K)\). Experimental results on synthetic and real-world data validate the effectiveness of PnP-DM, highlighting its ability to handle both linear and nonlinear inverse problems, including highly ill-posed cases.