#### NOTE
New to DeepInverse? Get started with the basics with the
[5 minute quickstart tutorial](https://deepinv.org/auto_examples/basics/demo_quickstart.html.md#sphx-glr-auto-examples-basics-demo-quickstart-py)..
# Uncertainty quantification with PnP-ULA.
This code shows you how to use sampling algorithms to quantify uncertainty of a reconstruction
from incomplete and noisy measurements.
ULA obtains samples by running the following iteration:
$$
x_{k+1} = x_k + \alpha \eta \nabla \log p_{\sigma}(x_k) + \eta \nabla \log p(y|x_k) + \sqrt{2 \eta} z_k
$$
where $z_k \sim \mathcal{N}(0, I)$ is a Gaussian random variable, $\eta$ is the step size and
$\alpha$ is a parameter controlling the regularization.
The PnP-ULA method is described in the paper Laumont *et al.*[1](#footcite-laumont2022bayesian).
```Python
import deepinv as dinv
from deepinv.utils.plotting import plot
import torch
from deepinv.utils import load_example
```
## Load image from the internet
This example uses an image of Messi.
```Python
device = dinv.utils.get_device()
x = load_example("messi.jpg", img_size=32).to(device)
```
```none
Selected GPU 0 with 24085.25 MiB free memory
```
## Define forward operator and noise model
This example uses inpainting as the forward operator and Gaussian noise as the noise model.
```Python
sigma = 0.1 # noise level
physics = dinv.physics.Inpainting(mask=0.5, img_size=x.shape[1:], device=device)
physics.noise_model = dinv.physics.GaussianNoise(sigma=sigma)
# Set the global random seed from pytorch to ensure reproducibility of the example.
torch.manual_seed(0)
```
```none
```
## Define the likelihood
Since the noise model is Gaussian, the negative log-likelihood is the L2 loss.
$$
-\log p(y|x) \propto \frac{1}{2\sigma^2} \|y-Ax\|^2
$$
```Python
# load Gaussian Likelihood
likelihood = dinv.optim.data_fidelity.L2(sigma=sigma)
```
## Define the prior
The score a distribution can be approximated using Tweedie’s formula via the
[`deepinv.optim.ScorePrior`](https://deepinv.org/api/stubs/deepinv.optim.ScorePrior.html.md#deepinv.optim.ScorePrior) class.
$$
\nabla \log p_{\sigma}(x) \approx \frac{1}{\sigma^2} \left(D(x,\sigma)-x\right)
$$
This example uses a pretrained DnCNN model.
From a Bayesian point of view, the score plays the role of the gradient of the
negative log prior
The hyperparameter `sigma_denoiser` ($sigma$) controls the strength of the prior.
In this example, we use a pretrained DnCNN model using the [`deepinv.loss.FNEJacobianSpectralNorm`](https://deepinv.org/api/stubs/deepinv.loss.FNEJacobianSpectralNorm.html.md#deepinv.loss.FNEJacobianSpectralNorm) loss,
which makes sure that the denoiser is firmly non-expansive (see Terris *et al.*[2](#footcite-terris2020building)), and helps to
stabilize the sampling algorithm.
```Python
sigma_denoiser = 2 / 255
prior = dinv.optim.ScorePrior(
denoiser=dinv.models.DnCNN(pretrained="download_lipschitz")
).to(device)
```
## Create the MCMC sampler
Here we use the Unadjusted Langevin Algorithm (ULA) to sample from the posterior defined in
[`deepinv.sampling.ULAIterator`](https://deepinv.org/api/stubs/deepinv.sampling.ULAIterator.html.md#deepinv.sampling.ULAIterator).
The hyperparameter `step_size` controls the step size of the MCMC sampler,
`regularization` controls the strength of the prior and
`iterations` controls the number of iterations of the sampler.
```Python
regularization = 0.9
step_size = 0.01 * (sigma**2)
iterations = int(5e3) if torch.cuda.is_available() else 10
params = {
"step_size": step_size,
"alpha": regularization,
"sigma": sigma_denoiser,
}
f = dinv.sampling.sampling_builder(
"ULA",
prior=prior,
data_fidelity=likelihood,
max_iter=iterations,
params_algo=params,
thinning=1,
verbose=True,
)
```
## Generate the measurement
We apply the forward model to generate the noisy measurement.
```Python
y = physics(x)
```
## Run sampling algorithm and plot results
The sampling algorithm returns the posterior mean and variance.
We compare the posterior mean with a simple linear reconstruction.
```Python
mean, var = f.sample(y, physics)
# compute linear inverse
x_lin = physics.A_adjoint(y)
# compute PSNR
print(f"Linear reconstruction PSNR: {dinv.metric.PSNR()(x, x_lin).item():.2f} dB")
print(f"Posterior mean PSNR: {dinv.metric.PSNR()(x, mean).item():.2f} dB")
# plot results
error = (mean - x).abs().sum(dim=1).unsqueeze(1) # per pixel average abs. error
std = var.sum(dim=1).unsqueeze(1).sqrt() # per pixel average standard dev.
imgs = [x_lin, x, mean, std / std.flatten().max(), error / error.flatten().max()]
plot(
imgs,
titles=["measurement", "ground truth", "post. mean", "post. std", "abs. error"],
)
```
```none
0%| | 0/5000 [00:00
* **References:**
* Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra. Bayesian imaging using plug & play priors: when langevin meets tweedie. *SIAM Journal on Imaging Sciences*, 15(2):701–737, 2022.
* Matthieu Terris, Audrey Repetti, Jean-Christophe Pesquet, and Yves Wiaux. Building firmly nonexpansive convolutional neural networks. In *ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*, 8658–8662. IEEE, 2020.
**Total running time of the script:** (0 minutes 3.571 seconds)
[`Download Jupyter notebook: demo_sampling.ipynb`](https://deepinv.org/auto_examples/sampling/_downloads/5a89e1000770b922d57eed794f82857e/demo_sampling.ipynb)
[`Download Python source code: demo_sampling.py`](https://deepinv.org/auto_examples/sampling/_downloads/5596c4a992481456ea481739b5152618/demo_sampling.py)
[`Download zipped: demo_sampling.zip`](https://deepinv.org/auto_examples/sampling/_downloads/199de052057d8e7a8b39ef1dbe23d4c9/demo_sampling.zip)
[Gallery generated by Sphinx-Gallery](https://sphinx-gallery.github.io)