Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
41 changes: 23 additions & 18 deletions cuqi/implicitprior/_restorator.py
Original file line number Diff line number Diff line change
Expand Up @@ -98,50 +98,55 @@ def get_conditioning_variables(self):


class MoreauYoshidaPrior(Distribution):
"""
r"""
This class defines (implicit) smoothed priors for which we can apply
gradient-based algorithms. The smoothing is performed using
the Moreau-Yoshida envelope of the target prior potential.

In the following we give a detailed explanation of the
Moreau-Yoshida smoothing.

We consider a density such that - \log\pi(x) = -g(x) with g convex, lsc,
We consider a density such that :math:`- \log\pi(x) = -g(x)` with g convex, lsc,
proper but not differentiable. Consequently, we cannot apply any
algorithm requiring the gradient of g.
Idea:
We consider the Moreau envelope of g defined as

g_{smoothing_strength} (x) = inf_z 0.5*\| x-z \|_2^2/smoothing_strength + g(z).
.. math::

g_{\mathrm{smoothing\_strength}} (x) = \inf_z 0.5 \| x-z \|_2^2/\mathrm{smoothing\_strength} + g(z).

g_{smoothing_strength} has some nice properties
- g_{smoothing_strength}(x)-->g(x) as smoothing_strength-->0 for all x
- \nabla g_{smoothing_strength} is 1/smoothing_strength-Lipschitz
- \nabla g_{smoothing_strength}(x) = (x - prox_g^{smoothing_strength}(x))/smoothing_strength for all x with
:math:`g_{\mathrm{smoothing\_strength}}` has some nice properties

- :math:`g_{\mathrm{smoothing\_strength}}(x) \to g(x)` as :math:`\mathrm{smoothing\_strength} \to 0` for all x
- :math:`\nabla g_{\mathrm{smoothing\_strength}}` is :math:`1/\mathrm{smoothing\_strength}`-Lipschitz
- :math:`\nabla g_{\mathrm{smoothing\_strength}}(x) = (x - \mathrm{prox}_g^{\mathrm{smoothing\_strength}}(x))/\mathrm{smoothing\_strength}` for all x with

prox_g^{smoothing_strength}(x) = argmin_z 0.5*\| x-z \|_2^2/smoothing_strength + g(z) .
.. math::

\mathrm{prox}_g^{\mathrm{smoothing\_strength}}(x) = \operatorname{argmin}_z 0.5 \| x-z \|_2^2/\mathrm{smoothing\_strength} + g(z) .

Consequently, we can apply any gradient-based algorithm with
g_{smoothing_strength} in lieu of g. These algorithms do not require the
full knowledge of g_{smoothing_strength} but only its gradient. The gradient
of g_{smoothing_strength} is fully determined by prox_g^{smoothing_strength}
and smoothing_strength.
:math:`g_{\mathrm{smoothing\_strength}}` in lieu of g. These algorithms do not require the
full knowledge of :math:`g_{\mathrm{smoothing\_strength}}` but only its gradient. The gradient
of :math:`g_{\mathrm{smoothing\_strength}}` is fully determined by :math:`\mathrm{prox}_g^{\mathrm{smoothing\_strength}}`
and :math:`\mathrm{smoothing\_strength}`.
It is important as, although there exists an explicit formula for
g_{smoothing_strength}, it is rarely used in practice, as it would require
:math:`g_{\mathrm{smoothing\_strength}}`, it is rarely used in practice, as it would require
us to solve an optimization problem each time we want to
estimate g_{smoothing_strength}. Furthermore, there exist cases where we dont't
the regularization g with which the mapping prox_g^{smoothing_strength} is
estimate :math:`g_{\mathrm{smoothing\_strength}}`. Furthermore, there exist cases where we dont't
the regularization g with which the mapping :math:`\mathrm{prox}_g^{\mathrm{smoothing\_strength}}` is
associated.

Remark (Proximal operators are denoisers):
We consider the denoising inverse problem x = u + n, with
n \sim \mathcal{N}(0, smoothing_strength I).
:math:`n \sim \mathcal{N}(0, \mathrm{smoothing\_strength} \, I)`.
A mapping solving a denoising inverse problem is called denoiser. It takes
the noisy observation x as an input and returns a less noisy version of x
which is an estimate of u.
We assume a prior density \pi(u) \propto exp(- g(u)).
We assume a prior density :math:`\pi(u) \propto \exp(- g(u))`.
Then the MAP estimate is given by
x_MAP = \argmin_z 0.5 \| x - z \|_2^2/smoothing_strength + g(z) = prox_g^smoothing_strength(x)
:math:`x_{\mathrm{MAP}} = \operatorname{argmin}_z 0.5 \| x - z \|_2^2/\mathrm{smoothing\_strength} + g(z) = \mathrm{prox}_g^{\mathrm{smoothing\_strength}}(x)`
Then proximal operators are denoisers.

Remark (Denoisers are not necessarily proximal operators): Data-driven
Expand Down
44 changes: 17 additions & 27 deletions demos/howtos/defining_posterior.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,8 @@
# Then consider the following Bayesian model
#
# .. math::
# \begin{align*}
# \mathbf{x} &\sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
# \mathbf{y} &\sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\,\mathbf{I})
# \end{align*}
# \mathbf{x} \sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
# \mathbf{y} \sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\,\mathbf{I})
#
# which can be written in CUQIpy as

Expand Down Expand Up @@ -64,10 +62,8 @@
# Suppose we had two forward models :math:`\mathbf{A}` and :math:`\mathbf{B}`:
#
# .. math::
# \begin{align*}
# \mathbf{y} &= \mathbf{A}\mathbf{x}\\
# \mathbf{d} &= \mathbf{B}\mathbf{x}\\
# \end{align*}
# \mathbf{y} = \mathbf{A}\mathbf{x}\\
# \mathbf{d} = \mathbf{B}\mathbf{x}

# Both observations come from the same unknown x
A, y_obs, _ = cuqi.testproblem.Deconvolution1D().get_components()
Expand All @@ -77,11 +73,9 @@
# Then consider the following Bayesian model
#
# .. math::
# \begin{align*}
# \mathbf{x} &\sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
# \mathbf{y} &\sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\mathbf{I})\\
# \mathbf{d} &\sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01^2\mathbf{I})
# \end{align*}
# \mathbf{x} \sim \mathcal{N}(\mathbf{0}, 0.1\,\mathbf{I})\\
# \mathbf{y} \sim \mathcal{N}(\mathbf{A}\mathbf{x}, 0.05^2\mathbf{I})\\
# \mathbf{d} \sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01^2\mathbf{I})

x = cuqi.distribution.Gaussian(np.zeros(A.domain_dim), 0.1)
y = cuqi.distribution.Gaussian(A(x), 0.05**2)
Expand Down Expand Up @@ -115,11 +109,9 @@
# we have the following 3 forward models
#
# .. math::
# \begin{align*}
# \mathbf{y} &= \mathbf{A}\mathbf{x}\\
# \mathbf{d} &= \mathbf{B}\mathbf{x}\\
# \mathbf{b} &= C(\mathbf{x})
# \end{align*}
# \mathbf{y} = \mathbf{A}\mathbf{x}\\
# \mathbf{d} = \mathbf{B}\mathbf{x}\\
# \mathbf{b} = C(\mathbf{x})
#
# where :math:`C` is a nonlinear function.

Expand All @@ -133,15 +125,13 @@
# Then consider the following Bayesian model
#
# .. math::
# \begin{align*}
# q &\sim \mathcal{U}(0.1, 10)\\
# l &\sim \mathrm{Gamma}(1, 1)\\
# s &\sim \mathrm{Gamma}(1, 10^{-2})\\
# \mathbf{x} &\sim \mathcal{N}(\mathbf{0}, l^{-1}\mathbf{I})\\
# \mathbf{y} &\sim \mathcal{N}(\mathbf{A}\mathbf{x}, s^{-1}\mathbf{I})\\
# \mathbf{d} &\sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01\mathbf{I})\\
# \mathbf{b} &\sim \mathcal{L}(\mathbf{C}(\mathbf{x}), q)
# \end{align*}
# q \sim \mathcal{U}(0.1, 10)\\
# l \sim \mathrm{Gamma}(1, 1)\\
# s \sim \mathrm{Gamma}(1, 10^{-2})\\
# \mathbf{x} \sim \mathcal{N}(\mathbf{0}, l^{-1}\mathbf{I})\\
# \mathbf{y} \sim \mathcal{N}(\mathbf{A}\mathbf{x}, s^{-1}\mathbf{I})\\
# \mathbf{d} \sim \mathcal{N}(\mathbf{B}\mathbf{x}, 0.01\mathbf{I})\\
# \mathbf{b} \sim \mathcal{L}(\mathbf{C}(\mathbf{x}), q)

q = cuqi.distribution.Uniform(0.1, 10)
l = cuqi.distribution.Gamma(1, 1)
Expand Down
54 changes: 26 additions & 28 deletions demos/howtos/gibbs.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,28 +2,28 @@
Gibbs sampling
==============

This tutorial shows how to use CUQIpy to perform Gibbs sampling.
Gibbs sampling is a Markov chain Monte Carlo (MCMC) method for
sampling a joint probability distribution.

Opposed to jointly sampling the distribution simultaneously, Gibbs
sampling samples the variables of the distribution sequentially,
one variable at a time. When a variable represents a random vector, the
whole vector is sampled simultaneously.
The sampling of each variable is done by sampling from the conditional
distribution of that variable given (fixed, previously sampled) values
of the other variables.

This is often a very efficient way of sampling from a joint
distribution if the conditional distributions are easy to sample
from. This is one way to exploit the structure of the joint
distribution. On the other hand, if the conditional distributions
are highly correlated and/or are difficult to sample from, then
Gibbs sampling can be very inefficient.

For these reasons, Gibbs sampling is often a double-edged sword,
that needs to be used in the right context.
This tutorial shows how to use CUQIpy to perform Gibbs sampling.
Gibbs sampling is a Markov chain Monte Carlo (MCMC) method for
sampling a joint probability distribution.

Opposed to jointly sampling the distribution simultaneously, Gibbs
sampling samples the variables of the distribution sequentially,
one variable at a time. When a variable represents a random vector, the
whole vector is sampled simultaneously.

The sampling of each variable is done by sampling from the conditional
distribution of that variable given (fixed, previously sampled) values
of the other variables.

This is often a very efficient way of sampling from a joint
distribution if the conditional distributions are easy to sample
from. This is one way to exploit the structure of the joint
distribution. On the other hand, if the conditional distributions
are highly correlated and/or are difficult to sample from, then
Gibbs sampling can be very inefficient.

For these reasons, Gibbs sampling is often a double-edged sword,
that needs to be used in the right context.

"""
# %%
Expand Down Expand Up @@ -77,12 +77,10 @@
#
# .. math::
#
# \begin{align}
# d &\sim \mathrm{Gamma}(1, 10^{-4}) \\
# l &\sim \mathrm{Gamma}(1, 10^{-4}) \\
# \mathbf{x} &\sim \mathrm{GMRF}(\mathbf{0}, d) \\
# \mathbf{y} &\sim \mathcal{N}(\mathbf{A} \mathbf{x}, l^{-1} \mathbf{I}_m)
# \end{align}
# d \sim \mathrm{Gamma}(1, 10^{-4}) \\
# l \sim \mathrm{Gamma}(1, 10^{-4}) \\
# \mathbf{x} \sim \mathrm{GMRF}(\mathbf{0}, d) \\
# \mathbf{y} \sim \mathcal{N}(\mathbf{A} \mathbf{x}, l^{-1} \mathbf{I}_m)
#
# where :math:`\mathbf{y}` is the observed data, and :math:`\mathbf{x}`
# is the unknown signal. The hyperparameters :math:`d` and :math:`l` are
Expand Down
2 changes: 1 addition & 1 deletion docs/api/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ The main modules are:
- :doc:`cuqi.distribution <_autosummary/cuqi.distribution>` for defining probability distributions.
- :doc:`cuqi.model <_autosummary/cuqi.model>` for defining deterministic models.
- :doc:`cuqi.sampler <_autosummary/cuqi.sampler>` for sampling from probability distributions.
- :doc:`cuqi.inference <_autosummary/cuqi.implicitprior>` for defining implicit priors.
- :doc:`cuqi.implicitprior <_autosummary/cuqi.implicitprior>` for defining implicit priors.

The following modules provide higher-level interfaces:
- :doc:`cuqi.testproblem <_autosummary/cuqi.testproblem>` for defining specific test problems.
Expand Down
Loading