11institutetext: Université Côte d’Azur, CNRS, LJAD, France
11email<EMAIL_ADDRESS>
# Python client for Isabelle server††thanks: This work has been supported by
the French government, through the 3IA Côte d’Azur Investments in the Future
project managed by the National Research Agency (ANR) with the reference
number ANR-19-P3IA-0002.

Boris Shminke 11 0000-0002-1291-9896

###### Abstract

We contribute a Python client for the Isabelle server, which gives researchers
and students using Python as their primary programming language an opportunity
to communicate with the Isabelle server through TCP directly from a Python
script. Such an approach helps avoid the complexities of integrating the
existing Python script with languages used for Isabelle development (ML and
Scala). We also describe new features that appeared since the announcement of
the first version of the client a year ago. Finally, we give examples of the
client’s applications in research and education and discuss known limitations
and possible directions for future development.

###### Keywords:

Isabelle Python client-server architecture

##  1 Introduction

Isabelle [8] interactive theorem prover (ITP) has included the Isabelle server
as part of its standard distribution since 2018 [15]. The Isabelle server
enables users to run multiple sessions and manage concurrent tasks to process
Isabelle theory files through TCP. It makes, in principle, possible to
communicate with the Isabelle server using any popular programming language
[12], including Python. Python clients already exist for other major ITPs, for
example, one [13] for Lean [7] or another one [9] for Coq [11]. To our best
knowledge, the previous version (0.2.0) of the isabelle-client announced last
year [6] was the first Python client for Isabelle. This paper describes
version 0.3.5 of the client (the digital artefact is available on Zenodo
[14]).

##  2 How to use

A typical scenario of isabelle-client application is represented by Fig. 2.
First, one should start the Isabelle server, for example, using a utility
function start_isabelle_server from the package. Second, one creates an
instance of IsabelleClient object, e.g. using the factory function
get_isabelle_client. Finally, one can issue any command supported by the
Isabelle server (see section 4.4 in [16] for a full list) using IsabelleClient
object, which implements all these commands as methods. Usually, these
commands rely on existing Isabelle theory files, for example, generated by
third-party Python scripts (not by isabelle-client). See also listing 1 for a
basic code snippet.

⬇

""" An example of the “isabelle-client“ usage """

from isabelle_client import get_isabelle_client, start_isabelle_server

\par# first, we start Isabelle server

server_info, _ = start_isabelle_server(

name="test", port=9999, log_file="server.log"

)

# then we create an “IsabelleClient“ instance

isabelle = get_isabelle_client(server_info)

# now we can send theory files to the server and get a response

isabelle.use_theories(theories=["Example"], master_dir=".")

# or we can build a session document using ROOT and root.tex files

isabelle.session_build(dirs=["."], session="examples")

isabelle.shutdown()

Figure 1: How to use isabelle-client.

Interested potential users of the client can also follow the package homepage
111https://pypi.org/project/isabelle-client/ for more detailed information.

Figure 2: In a project where Python scripts are generating Isabelle theory
files, isabelle-client can be used to start the Isabelle server, sending these
files for parallel processing and getting back the results.

##  3 New features since previous version

The first public version of this client [6], which appeared in March 2021,
worked only for Python 3.7 on GNU/Linux and was supposed to be installed only
with the pip package manager. The current version is available for any Python
3.6+ on GNU/Linux and Windows. Every new build is tested in a continuous
integration workflow against each supported Python version. The package is
hosted now not only on the Python Package Index (PyPI), but also on Conda
Forge [1], which enables its installation with both pip and conda package
managers. In addition, one can run the client inside a Docker container, for
example, in a cloud using Binder [5]. This option provides the client coupled
with the Isabelle server and is particularly useful for students specialising
in logic but not necessarily having much experience in information
technologies.

The current version of the client is tested to work with the latest Isabelle
2021-1 (released in December 2021). The last client version returns all server
replies as a Python list, not only the last one as it was in the previous
version giving more flexibility to the end-user. Last but not least, the
current version arrives with detailed documentation pages and is nearly 100%
covered with unit tests using fixtures for emulating a working Isabelle server
behaviour.

##  4 Known applications in research and education

The original system entry for isabelle-client expressed hope that it would be
helpful to other projects. After the first year, we saw several applications
confirming these aspirations. The most impressive one was a discovery of a
proof [2] for an algebraic problem which stood open for two years despite the
efforts of specialists in the field. For this project, the isabelle-client
played a role of a bridge connecting Python scripts which generated relevant
Isabelle theory files with the Isabelle server processing them.

Another use case comes from higher education. The isabelle-client running in a
Docker container on Binder was used during the practical sessions of the
Advanced Logic course taught at the Université Côte d’Azur in the autumn of
the 2021-2022 academic year. The client helped students not trained in
functional programming languages used for Isabelle development to concentrate
on understanding the Isabelle language syntax and consequently generating
theory files in it with Python scripts without installing and running the
Isabelle GUI on their laptops.

Also, Fabian Huch used the isabelle-client for debugging the ’Proving for Fun’
backend [4]. The ’Proving for Fun’ project already had in its codebase parts
resembling the isabelle-client, but using the socket Python package that
provides access to the BSD socket interface. The isabelle-client uses asyncio,
providing a more high-level interface for TCP communications, and exists as a
stand-alone piece of software striving to be easily reusable for different
projects.

##  5 Known limitations, workarounds, and future work

The first thing one would want from a client-server application is running the
server and the client on different machines (probably, with various operating
systems) or, at least, in separate Docker containers. The attempts to do that
with isabelle-client demonstrated that the current implementation of the
Isabelle server doesn’t support binding its IP address to anything but
localhost for security reasons. A desirable solution is implementing an
authentication mechanism for the remote Isabelle server (for example, with
access tokens, as is done for the Jupyter [3] server). The current workaround
is to use SSH tunnels or package the isabelle-client together with the
Isabelle server in the same Docker container.

Although the isabelle-client works perfectly on Windows for communication with
the Isabelle server running on the same machine, it’s impossible to start the
Isabelle server on Windows with the utility function from the isabelle-client
package. The reason is related to the fact that the Isabelle server on Windows
runs under Cygwin and not as a proper Windows application. A possible solution
could be to rewrite the Windows CMD file, which starts the Isabelle server, to
make it pass the started server info from the Cygwin subsystem to the parent
Windows process. The current workaround is to launch the Isabelle server
manually.

At some point, it became known that it’s impossible to start a clean build of
a session from the isabelle-client. The reason is that the Isabelle server
(unlike isabelle build tool) lacks -c parameter in its API. Thus, achieving
relative feature parity of the Isabelle server and other Isabelle tools seems
to be an important direction in making the isabelle-client more applicable and
interesting to its end users.

##  6 Conclusion

A year after its appearance isabelle-client became a more mature piece of
software and found its applications in both scientific research and higher
education. As shown in this paper, there are still many ways to improve it,
and we hope for many possible occasions for its further application, not
limited to but including machine learning (as a domain dominated by Python
[10]) for theorem proving.

## References

  * [1] conda-forge community: The conda-forge Project: Community-based Software Distribution Built on the conda Package Format and Ecosystem (Jul 2015). https://doi.org/10.5281/zenodo.4774216, https://doi.org/10.5281/zenodo.4774216
  * [2] Fussner, W., Shminke, B.: Mining counterexamples for wide-signature algebras with an isabelle server. CoRR abs/2109.05264 (2021), https://arxiv.org/abs/2109.05264
  * [3] Granger, B.E., Pérez, F.: Jupyter: Thinking and storytelling with code and data. Comput. Sci. Eng. 23(2), 7–14 (2021). https://doi.org/10.1109/MCSE.2021.3059263, https://doi.org/10.1109/MCSE.2021.3059263
  * [4] Haslbeck, M.P., Wimmer, S.: Competitive proving for fun. Kalpa Publications in Computing 10, 9–14 (2019)
  * [5] Project Jupyter, Matthias Bussonnier, Jessica Forde, Jeremy Freeman, Brian Granger, Tim Head, Chris Holdgraf, Kyle Kelley, Gladys Nalvarte, Andrew Osheroff, Pacer, M., Yuvi Panda, Fernando Perez, Benjamin Ragan Kelley, Carol Willing: Binder 2.0 - Reproducible, interactive, sharable environments for science at scale. In: Fatih Akici, David Lippa, Dillon Niederhut, Pacer, M. (eds.) Proceedings of the 17th Python in Science Conference. pp. 113 – 120 (2018). https://doi.org/10.25080/Majora-4af1f417-011
  * [6] Líska, M., Lupták, D., Novotný, V., Ruzicka, M., Shminke, B., Sojka, P., Stefánik, M., Wenzel, M.: Cicm’21 systems entries. In: Kamareddine, F., Coen, C.S. (eds.) Intelligent Computer Mathematics - 14th International Conference, CICM 2021, Timisoara, Romania, July 26-31, 2021, Proceedings. Lecture Notes in Computer Science, vol. 12833, pp. 245–248. Springer (2021). https://doi.org/10.1007/978-3-030-81097-9_20, https://doi.org/10.1007/978-3-030-81097-9_20
  * [7] Moura, L.d., Ullrich, S.: The lean 4 theorem prover and programming language. In: Platzer, A., Sutcliffe, G. (eds.) Automated Deduction – CADE 28. pp. 625–635. Springer International Publishing, Cham (2021)
  * [8] Nipkow, T., Paulson, L.C., Wenzel, M.: Isabelle/HOL - A Proof Assistant for Higher-Order Logic, Lecture Notes in Computer Science, vol. 2283. Springer (2002). https://doi.org/10.1007/3-540-45949-9, https://doi.org/10.1007/3-540-45949-9
  * [9] Gallego Arias, E.J., Martinez, T.: Pycoq: Access coq from python! (Jan 2022), https://github.com/ejgallego/pycoq
  * [10] Kaggle Inc.: State of data science and machine learning 2021 (Oct 2021), https://www.kaggle.com/kaggle-survey-2021
  * [11] The Coq Development Team: The coq proof assistant (Jan 2022). https://doi.org/10.5281/zenodo.5846982, https://doi.org/10.5281/zenodo.5846982
  * [12] TIOBE Software BV: The tiobe programming community index (Apr 2022), https://www.tiobe.com/tiobe-index/
  * [13] Rute, J., Massot, P., Berman, J., Roux, F.L.: Lean client for python (Aug 2021), https://github.com/leanprover-community/lean-client-python
  * [14] Shminke, B.: Python client for isabelle server (Apr 2022). https://doi.org/10.5281/zenodo.6490275, https://doi.org/10.5281/zenodo.6490275
  * [15] Wenzel, M.: Isabelle/pide after 10 years of development. In: 13th International Workshop on User Interfaces for Theorem Provers (UITP 2018). Federated Logic Conference 2018 (2018), https://sketis.net/wp-content/uploads/2018/08/isabelle-pide-uitp2018.pdf
  * [16] Wenzel, M.: The isabelle system manual (Dec 2021), https://isabelle.in.tum.de/dist/Isabelle2021-1/doc/system.pdf


# Denoising Score-Matching for Uncertainty Quantification in Inverse Problems

Zaccharie Ramzi
CEA (Neurospin and Cosmostat), Inria (Parietal)
Gif-sur-Yvette, France
<EMAIL_ADDRESS>
&Benjamin Remy, François Lanusse, Jean-Luc Starck
AIM, CEA, CNRS
Université Paris-Saclay
Université Paris Diderot
Sorbonne Paris Cité
{benjamin.remy, francois.lanusse<EMAIL_ADDRESS>&Philippe Ciuciu
CEA (Neurospin), Inria (Parietal)
Gif-sur-Yvette, France
<EMAIL_ADDRESS>https://www.cosmostat.org/people/zaccharie-ramzi

###### Abstract

Deep neural networks have proven extremely efficient at solving a wide range
of inverse problems, but most often the uncertainty on the solution they
provide is hard to quantify. In this work, we propose a generic Bayesian
framework for solving inverse problems, in which we limit the use of deep
neural networks to learning a prior distribution on the signals to recover. We
adopt recent denoising score matching techniques to learn this prior from
data, and subsequently use it as part of an annealed Hamiltonian Monte-Carlo
scheme to sample the full posterior of image inverse problems. We apply this
framework to Magnetic Resonance Image (MRI) reconstruction and illustrate how
this approach not only yields high quality reconstructions but can also be
used to assess the uncertainty on particular features of a reconstructed
image.

##  1 Introduction

Generative modeling has enjoyed remarkable success in recent years with models
such as Generative Adversarial Networks (GANs) (Goodfellow et al., 2014)
reaching extremely high quality results on complex high resolution images
(Karras et al., 2020). Yet, GANs are still prone to issues including unstable
training and mode collapse, i.e. a lack of diversity in generated images. In
addition, GANs only provide a convenient way to sample from the learned
distribution, they do not give access to the density function itself. Another
class of commonly used models, Variational Auto-Encoders (VAEs) (Kingma and
Welling, 2014) have also been able to reach high quality sample resolutions
(Vahdat and Kautz, 2020), but can only provide a bound on the density function
of the model.

In some practical applications of these generative models however, the most
desirable aspect is not necessarily to be able to sample from the model, or
directly have access to its density functions, but have access to gradients of
the log density function, i.e. its score. In particular, in this work we focus
on the generic problem of performing Bayesian inference to solve an inverse
problem (i.e. deconvolution, denoising, inpainting, etc), using a learned
generative model as a Bayesian prior. Such problems can typically be solved by
gradient-based inference methods such as Variational Inference (VI), and
Langevin or Hamiltonian Monte-Carlo (HMC), but all these methods have in
common that they only rely on the score of the generative model.

Recent work (Lim et al., 2020; Song and Ermon, 2019) has investigated various
approaches to directly train a deep neural network (coined, a score network)
to estimate the score of a density function instead of the density itself. One
of the most promising and scalable approach relies on denoising score-
matching, where the score network is trained using a denoising loss on an
augmented dataset where various amounts of Gaussian noise is added to the
data. Not only has this approach already been demonstrated to reach convincing
results in score estimation, but Song and Ermon (2020) for instance
demonstrated that these score networks can be used to sample high resolution
natural images through Langevin dynamics.

In this work, we develop the first application of deep denoising score
matching to Magnetic Resonance Image (MRI) reconstruction and its use for
Uncertainty Quantification (UQ) in the context of imaging inverse problems.

##  2 Related Works

#### Score Matching

Score Matching was introduced by Hyvärinen (2005). It was later revisited in
Vincent (2011) and Alain and Bengio (2013) in the context of Deep Learning and
Denoising Auto-Encoders (DAE) in particular. Recently, Song and Ermon (2019)
applied these findings to images showing the potential of this approach
combined with sampling techniques. They refined their work in Song and Ermon
(2020) to take it to high dimension and provide some techniques for hyper-
parameter setting.

#### Plug-and-Play Priors

The idea of using a denoising model as a prior for inverse problems solving
was introduced by Venkatakrishnan et al. (2013) who used non-learned denoisers
to replace the proximity operator in the ADMM algorithm (Boyd, 2011). More
recently, Meinhardt et al. (2017) have made use of learned denoiser networks
in Plug-and-Play inverse problems solving. Annealed HMC sampling can be seen
as a fuzzy-version of these approaches with theoretical grounds.

#### Deep Uncertainty Quantification in Inverse Problems

The idea of using generative models for inverse problems such as data
imputation and denoising was for instance suggested in an early VAE work
(Rezende et al., 2014). In this approach, a generative network is first
trained on high quality data, and inference is then performed in the latent
space of the model typically using VI. The main advantage comes from the
reduced dimensionality of the parameter space. Recent examples of this
approach include (Wu et al., 2018; Böhm et al., 2019), and an application to
MRI can for instance be found in (Edupuganti et al., 2019). An important
aspect of this approach is that it relies on the known likelihood of
observations.

A different approach was proposed in Adler and Öktem (2018) based on a
conditional Wasserstein Generative Adversarial Network (cWGAN). In this
formulation, a generative model is trained to sample high quality images
conditioned on degraded observations. At test time, independent samples from
the posterior distribution induced by the GAN are obtaind by simply sampling
different latent space variables. Similar frameworks based on other
conditional generative models have alse been proposed, such as (Denker et al.,
2020) using conditional Normalizing Flows, or (Tonolini et al., 2019) using
conditional VAEs. While this conditional generative model approach allows for
fast inference, thanks to ancestral sampling, it is worth noting that none of
these models include an explicit data consistency step, nor do they make use
of a test-time likelihood. As a direct consequence, for instance, these models
cannot be used on data with different noise levels as seen during training.

##  3 Deep Denoising Score Matching for Posterior Inference

Our main objective in this work is to perform probabilistic inference over a
Bayesian model typically represented as

$p({\mathb{x}}|{\mathb{y}})\propto p({\mathb{y}}|{\mathb{x}})p({\mathb{x}})$
(1)

where ${\mathb{y}}$ are some measurements and the posterior
$p({\mathb{x}}|{\mathb{y}})$ is the distribution of possible solutions
${\mathb{x}}$ compatible with observations and $p({\mathb{x}})$ a prior
knowledge. The problem-specific likelihood $p({\mathb{y}}|{\mathb{x}})$
encodes the forward process of the model and accounts for observational noise,
while the prior $p({\mathb{x}})$ encapsulates any a priori information we have
on the solution of the problem. Our goal for UQ is to sample solutions
${\mathb{x}}$ belonging to that posterior. Multiple inference techniques can
be leveraged for sampling from this posterior $p({\mathb{x}}|{\mathb{y}})$,
but for high-dimensional problems modern techniques rely on gradient-based
methods, including Variational Inference (Hoffman et al., 2013) and Langevin
Diffusion, or HMC (Neal, 2011). All of these techniques have in common that
they only require having access to the score $\nabla_{x}\log p(x)$ of the
target distribution. We illustrate in particular in Appendix B how Metropolis
Hastings-calibrated HMC can be conducted with only knowledge of the score.
Two terms will enter in the score of the posterior distribution in Equation 1,
the score of the likelihood, and the score of the prior. In many problems,
such as in the MRI problem presented later in this work, the likelihood score
can be derived analytically, only the score of the prior remains unknown but
can be learned from data by score matching.

#### Deep Denoising Score Matching

As originally identified by Vincent (2011) and Alain and Bengio (2013), the
score of a given target distribution $P$ can be modeled using a DAE, i.e. by
introducing an auto-encoding function
$\boldsymbol{r}:\mathbb{R}^{n}\times\mathbb{R}\mapsto\mathbb{R}^{n}$ trained
to reconstruct under an $\ell_{2}$ loss a true ${\mathb{x}}\sim P$ given a
noisy version ${\mathb{x}}^{\prime}={\mathb{x}}+\boldsymbol{n}$ with
$\boldsymbol{n}\sim\mathcal{N}(0,\sigma^{2}\boldsymbol{I})$. An optimal
denoiser $\boldsymbol{r}^{\star}$ would then be achieved for:

$\boldsymbol{r}^{\star}({\mathb{x}}^{\prime},\sigma)={\mathb{x}}^{\prime}+\sigma^{2}\nabla_{{\mathb{x}}}\log
p_{\sigma^{2}}({\mathb{x}}^{\prime})$ (2)

where $p_{\sigma^{2}}=p\ast\mathcal{N}(0,\sigma^{2})$. In other words, the
optimal denoiser is closely related to the score we wish to learn and when the
noise variance $\sigma^{2}$ tends to zero, should exactly match the score of
the target density. In practice to learn this score efficiently, we adopt the
residual noise-conditional denoising score matching technique proposed by Lim
et al. (2020), and train a model to minimize the following Denoising Score
Matching loss:

$\mathcal{L}_{DSM}=\underset{{\mathb{x}}\sim
P}{\mathbb{E}}\underset{\begin{subarray}{c}\boldsymbol{u}\sim\mathcal{N}(0,I)\\\
\sigma_{s}\sim\mathcal{N}(0,s^{2})\end{subarray}}{\mathbb{E}}\parallel\boldsymbol{u}+\sigma_{s}\boldsymbol{r}_{\theta}({\mathb{x}}+\sigma_{s}\boldsymbol{u},\sigma_{s})\parallel_{2}^{2}$
(3)

In this formulation, the network $\boldsymbol{r}_{\theta}({\mathb{x}},\sigma)$
is now directly modeling the score $\nabla_{{\mathb{x}}}\log
p_{\sigma^{2}}({\mathb{x}})$ of the Gaussian-convolved target distribution.

#### Annealed Hamiltonian-Monte Carlo Sampling

Given the noise-conditional neural scores learned with the procedure described
above, it is now possible to use a variety of inference methods to access the
Bayesian posterior. In this work, we adopt an annealed HMC procedure which
provides an efficient way to obtain parallel batches of independent samples
from the target posterior despite the high dimensionality of the problem. This
is a sampling procedure closely related to the Annealed Langevin Diffusion
proposed in (Song and Ermon, 2019, 2020), but benefits from the faster
Hamiltonian dynamics and Metropolis-Hastings calibration (see Appendix B).

To build our procedure, we consider a Gaussian-convolved version of our target
density:

$\log p_{\sigma^{2}}({\mathb{x}}|{\mathb{y}})=\log
p_{\sigma^{2}}({\mathb{y}}|{\mathb{x}})+\log p_{\sigma^{2}}({\mathb{x}})+cst$
(4)

where $\sigma^{2}$ plays the role of the inverse temperature found in
classical annealing. The likelihood can be obtained analytically, and in the
case of a Gaussian likelihood takes the following form: $\log
p_{\sigma^{2}}({\mathb{y}}|{\mathb{x}})=-\frac{\parallel{\mathb{x}}-f({\mathb{x}})\parallel_{2}^{2}}{2(\sigma_{n}^{2}+\sigma^{2})}+cst$,
where $\sigma_{n}^{2}$ is the noise variance in the measurements. As for the
prior term, the noise-conditional score network introduced above already
models the score of $\log p_{\sigma^{2}}$. This distribution is gradually
annealed to low temperatures and the chain progressively moves towards a point
in the target distribution.

##  4 Application to Bayesian Inverse Problems

#### The MRI Problem

Magnetic Resonance Imaging (MRI) is a non-invasive modality used to probe soft
tissues. Compressed sensing, introduced for MRI by Lustig et al. (2007) is
used to reduce its significant acquisition time, and recently, deep learning
approaches(see for example (Schlemper et al., 2018; Hammernik et al., 2018;
Pezzotti et al., 2020)) have been shown to perform extremely well on the
reconstruction problem. The idealized reconstruction problem is usually
formalized in the following way, for the single coil setting, with uniform
acquisition:

${\mathb{y}}=M_{\Omega}{\mathb{F}}{\mathb{x}}+\boldsymbol{n}$ (5)

where ${\mathb{y}}$ is the acquired Fourier coefficients, also called the
k-space data, $M_{\Omega}$ is a mask, ${\mathb{F}}$ is the classical 2D
Fourier transform, ${\mathb{x}}$ is the real anatomical image, and
$\boldsymbol{n}\sim\mathcal{N}(0,\sigma_{n}^{2})$ is measurement noise (we set
$\sigma_{n}=0.1$ in our experiments). The data we consider is the single coil
data with the Proton-Density (PD) contrast from the fastMRI dataset introduced
in (Zbontar et al., 2018). For comparison with deep learning approaches, we
use the state-of-the-art Primal-Dual net enhanced with a U-net for image
correction (UPDNet) from (Ramzi et al., 2020). The under-sampling was done
retrospectively using an acceleration factor of $4$ and a random mask as
described for knee images in (Zbontar et al., 2018).

#### Results

We use a simple residual U-net inspired by Ronneberger et al. (2015) with
ResNet building blocks from (He et al., 2016) as our score network. In order
to promote the regularity of the learned scores, we regularize the spectral
norm of each convolutional layer, and we find that setting the spectral norm
to $\simeq 2$ yields the best results. More details about the network
architecture and training can be found in appendix A. We then sample from the
posterior following our annealed HMC procedure down to relatively low
temperatures, and apply one final denoising step on the last sample from the
chain, following the Expected Denoised Sample (EDS) scheme detailed in
(Jolicoeur-Martineau et al., 2020).

Figure 1:  Bayesian posterior sampling for MRI reconstruction. The top
leftmost image is the ground truth image. The top second to the left image is
the zero-filled retrospectively under-sampled image
${\mathb{F}}^{T}{\mathb{y}}$. The top third to the left image is the
reconstruction of the under-sampled image by the UPDNet. All the other images
are denoised samples from the estimated posterior distribution obtained by a
tempered HMC. The zero-filling achieves a PSNR of 25.55 dB, each sample 27.63
dB on average, the mean of the samples 30.04 dB and the neural network 32.15
dB. A zoom of the region in the red square is provided in appendix A.4. An
animation of posterior samples is available online .

Figure 1 compares samples from the MRI posterior to the ground truth (top
left), zero-filled image (top second left), and UPDNet reconstruction (top
center). We see that although individual samples carry slightly less details
than the neural network reconstruction, confidence in any particular part of
the image can be gauged by looking for variability or stability across
multiple independent posterior samples. We highlight in particular the red
region, where posterior samples show significant variability, indicating that
this part of the image is poorly constrained by data. In contrast, a direct
neural network reconstruction (top center) does not match the ground truth in
that region, and does not provide an estimate of uncertainty which may lead a
physician to misinterpret the image.

##  5 Conclusions and Discussions

We have presented in this work the first (to the best of our knowledge)
instance of a framework for Bayesian inverse problems based on Deep Denoising
Score Matching and applied to MRI reconstruction and Uncertainty
Quantification. We illustrated the merits of this approach on an MRI example
where the ability to sample from the full posterior highlights what features
present in a reconstruction are not actually constrained by data.
We also showed that the Deep Denoising Score Matching approach can
effortlessly be used for dimensions as high as $320\times 320\times 2$. It is
therefore much more scalable than GAN- or VAE-based models. We note however
that finding an optimal HMC annealing schedule and temperature-adaptive step-
size proved difficult. In general we did not manage to lower the temperature
below a certain level, prompting us to resort to a denoising step on the last
samples from our HMC chains. This is a direction for further research.

## Acknowledgments and Disclosure of Funding

We acknowledge the financial support of the Cross-Disciplinary Program on
Numerical Simulation of CEA for the project entitled SILICOSMIC. We also
acknowledge the French Institute of development and ressources in scientific
computing (IDRIS) for their AI program allowing us to use the Jean Zay
supercomputer’s GPU partitions.

## References

  * Goodfellow et al. (2014) Ian J. Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Mirza, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative Adversarial Nets.  _Advances in neural information processing systems_ , 2014.  URL https://arxiv.org/pdf/1406.2661.pdf.
  * Karras et al. (2020) Tero Karras, Samuli Laine, Miika Aittala, Janne Hellsten, Jaakko Lehtinen, and Timo Aila.  Analyzing and Improving the Image Quality of StyleGAN.  pages 8107–8116, 2020.  doi: 10.1109/cvpr42600.2020.00813.
  * Kingma and Welling (2014) Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  _2nd International Conference on Learning Representations, ICLR 2014 - Conference Track Proceedings_ , (Ml):1–14, 2014.
  * Vahdat and Kautz (2020) Arash Vahdat and Jan Kautz.  NVAE: A Deep Hierarchical Variational Autoencoder.  pages 1–20, 7 2020.  URL http://arxiv.org/abs/2007.03898.
  * Lim et al. (2020) Jae Hyun Lim, Aaron Courville, Christopher Pal, and Chin-Wei Huang.  AR-DAE: Towards Unbiased Neural Entropy Gradient Estimation.  2020\.  URL http://arxiv.org/abs/2006.05164.
  * Song and Ermon (2019) Yang Song and Stefano Ermon.  Generative Modeling by Estimating Gradients of the Data Distribution.  (NeurIPS), 2019.  URL http://arxiv.org/abs/1907.05600.
  * Song and Ermon (2020) Yang Song and Stefano Ermon.  Improved Techniques for Training Score-Based Generative Models.  2020\.  URL http://arxiv.org/abs/2006.09011.
  * Hyvärinen (2005) Aapo Hyvärinen.  Estimation of Non-Normalized Statistical Models by Score Matching.  _Journal of Machine Learning Research_ , 6:695–708, 2005\.  URL http://jmlr.org/papers/volume6/hyvarinen05a/old.pdf.
  * Vincent (2011) Pascal Vincent.  A connection between scorematching and denoising autoencoders.  _Neural Computation_ , 23(7):1661–1674, 2011\.  ISSN 08997667.  doi: 10.1162/NECO{\\_}a{\\_}00142.
  * Alain and Bengio (2013) Guillaume Alain and Yoshua Bengio.  What regularized auto-encoders learn from the data generating distribution.  _1st International Conference on Learning Representations, ICLR 2013 - Conference Track Proceedings_ , 15:3743–3773, 2013.
  * Venkatakrishnan et al. (2013) Singanallur V. Venkatakrishnan, Charles A. Bouman, and Brendt Wohlberg.  Plug-and-Play priors for model based reconstruction.  _2013 IEEE Global Conference on Signal and Information Processing, GlobalSIP 2013 - Proceedings_ , pages 945–948, 2013.  doi: 10.1109/GlobalSIP.2013.6737048.
  * Boyd (2011) Stephen Boyd.  Alternating Direction Method of Multipliers.  _IEEE Conference on Decision and Control_ , 3:1–44, 2011\.  ISSN 1935-8237.  doi: 10.1145/2783258.2783400.
  * Meinhardt et al. (2017) Tim Meinhardt, Michael Moeller, and Daniel Cremers.  Learning Proximal Operators : Using Denoising Networks for Regularizing Inverse Imaging Problems.  (12):1781–1790, 2017.
  * Rezende et al. (2014) Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.  Stochastic Backpropagation and Approximate Inference in Deep Generative Models.  _31st International Conference on Machine Learning, ICML 2014_ , 4:3057–3070, 1 2014.  URL http://arxiv.org/abs/1401.4082.
  * Wu et al. (2018) Ga Wu, Justin Domke, and Scott Sanner.  Conditional Inference in Pre-trained Variational Autoencoders via Cross-coding.  5 2018.  URL http://arxiv.org/abs/1805.07785.
  * Böhm et al. (2019) Vanessa Böhm, François Lanusse, and Uroš Seljak.  Uncertainty Quantification with Generative Models.  (NeurIPS):1–9, 10 2019.  URL http://arxiv.org/abs/1910.10046.
  * Edupuganti et al. (2019) Vineet Edupuganti, Morteza Mardani, Shreyas Vasanawala, and John Pauly.  Uncertainty Quantification in Deep MRI Reconstruction.  pages 1–10, 2019.  URL http://arxiv.org/abs/1901.11228.
  * Adler and Öktem (2018) Jonas Adler and Ozan Öktem.  Deep Bayesian Inversion.  11 2018.  URL http://arxiv.org/abs/1811.05910.
  * Denker et al. (2020) Alexander Denker, Maximilian Schmidt, Johannes Leuschner, Peter Maass, and Jens Behrmann.  Conditional Normalizing Flows for Low-Dose Computed Tomography Image Reconstruction.  6 2020.  URL http://arxiv.org/abs/2006.06270.
  * Tonolini et al. (2019) Francesco Tonolini, Jack Radford, Alex Turpin, Daniele Faccio, and Roderick Murray-Smith.  Variational Inference for Computational Imaging Inverse Problems.  21:1–46, 4 2019.  URL http://arxiv.org/abs/1904.06264.
  * Hoffman et al. (2013) Matthew D. Hoffman, David M. Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 14:1303–1347, 2013\.  ISSN 15324435.
  * Neal (2011) Radford M. Neal.  MCMC using hamiltonian dynamics.  _Handbook of Markov Chain Monte Carlo_ , pages 113–162, 2011.  doi: 10.1201/b10905-6.
  * Lustig et al. (2007) Michael Lustig, David Donoho, and John M. Pauly.  Sparse MRI: The Application of Compressed Sensing for Rapid MR Imaging Michael.  _Magnetic Resonance in Medicine_ , 2007.  ISSN 10014861.  doi: 10.1002/mrm.21391.
  * Schlemper et al. (2018) Jo Schlemper, Jose Caballero, Joseph V Hajnal, Anthony Price, and Daniel Rueckert.  A Deep Cascade of Convolutional Neural Networks for MR Image Reconstruction.  _IEEE Transactions on Medical Imaging_ , 37(2):491–503, 2018.  doi: 10.1109/TMI.2017.2760978.
  * Hammernik et al. (2018) Kerstin Hammernik, Teresa Klatzer, Erich Kobler, Michael P Recht, Daniel K Sodickson, Thomas Pock, and Florian Knoll.  Learning a Variational Network for Reconstruction of Accelerated MRI Data.  _Magnetic Resonance in Medicine_ , 3071(September 2017):3055–3071, 2018.  doi: 10.1002/mrm.26977.
  * Pezzotti et al. (2020) Nicola Pezzotti, Sahar Yousefi, Mohamed S Elmahdy, Jeroen Van Gemert, Christophe Schülke, Mariya Doneva, Tim Nielsen, Sergey Kastryulin, Boudewijn P F Lelieveldt, Matthias J P Van Osch, Elwin De Weerdt, and Marius Staring.  An Adaptive Intelligence Algorithm for Undersampled Knee MRI Reconstruction: Application to the 2019 fastMRI Challenge.  Technical report, 2020.
  * Zbontar et al. (2018) Jure Zbontar, Florian Knoll, Anuroop Sriram, Matthew J Muckley, Mary Bruno, Aaron Defazio, Marc Parente, Krzysztof J Geras, Joe Katsnelson, Hersh Chandarana, Zizhao Zhang, Michal Drozdzal, Adriana Romero, Michael Rabbat, Pascal Vincent, James Pinkerton, Duo Wang, Nafissa Yakubova, Erich Owens, C Lawrence Zitnick, Michael P Recht, Daniel K Sodickson, and Yvonne W Lui.  fastMRI: An Open Dataset and Benchmarks for Accelerated MRI.  Technical report, 2018.  URL https://arxiv.org/pdf/1811.08839.pdf.
  * Ramzi et al. (2020) Zaccharie Ramzi, Philippe Ciuciu, and Jean Luc Starck.  Benchmarking MRI reconstruction neural networks on large public datasets.  _Applied Sciences (Switzerland)_ , 10(5), 2020.  ISSN 20763417.  doi: 10.3390/app10051816.
  * Ronneberger et al. (2015) Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-Net: Convolutional Networks for Biomedical Image Segmentation.  Technical report, 2015.  URL http://lmb.informatik.uni-freiburg.de/.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  _Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , 2016-Decem:770–778, 2016.  ISSN 10636919.  doi: 10.1109/CVPR.2016.90.
  * Jolicoeur-Martineau et al. (2020) Alexia Jolicoeur-Martineau, Rémi Piché-Taillefer, Rémi Tachet des Combes, and Ioannis Mitliagkas.  Adversarial score matching and improved sampling for image generation.  pages 1–29, 2020.  URL http://arxiv.org/abs/2009.05475.
  * Yoshida and Miyato (2017) Yuichi Yoshida and Takeru Miyato.  Spectral Norm Regularization for Improving the Generalizability of Deep Learning.  2017\.  URL http://arxiv.org/abs/1705.10941.
  * Ryu et al. (2019) Ernest K. Ryu, Jialin Liu, Sicheng Wang, Xiaohan Chen, Zhangyang Wang, and Wotao Yin.  Plug-and-play methods provably converge with properly trained denoisers.  _36th International Conference on Machine Learning, ICML 2019_ , 2019-June:9764–9775, 2019.
  * Metropolis et al. (1953) Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller.  Equation of state calculations by fast computing machines.  _The Journal of Chemical Physics_ , 21(6):1087–1092, 1953.  ISSN 00219606.  doi: 10.1063/1.1699114.
  * Hastings (1970) W. K. Hastings.  Monte carlo sampling methods using Markov chains and their applications.  _Biometrika_ , 57(1):97–109, 1970.  ISSN 00063444.  doi: 10.1093/biomet/57.1.97.

##  Appendix A MRI experiments

###  A.1 Network architecture

The network is a 3-scale U-net with residual blocks composed of 3 convolutions
followed by a batch normalisation. Each batch normalisation is followed by a
ReLU non-linearity except the last one. A projection is used for the input of
residual blocks whose number of input channels is not the same as that of
output channels. Each scale uses 2 residual blocks for the downsampling path
and 2 residual blocks for the upsampling path. Downsampling is performed via
average pooling and upsampling is performed via up-convolution (as designed in
[Ronneberger et al., 2015]), in order to avoid checkerboard artefacts. We use
the following sequence of number of channels: [32, 64, 128]. In order to deal
with MRI images, that are complex in nature, we concatenate the real and
imaginary part of the image at the input of the network, forming effectively a
2-channel $320\times 320$ image. For input noise level conditioning, we
concatenate a noise standard deviation map to the input and in the lower scale
of the network. Following the recommendation of Song and Ermon [2020], we also
divide the output of the network by the absolute value of the noise power.
Finally, we use Spectral Norm regularisation from [Yoshida and Miyato, 2017]
in order to make sure the score does not take inconsistent values in low-
density regions. The Spectral Norm regularisation indeed forces the spectral
norm (maximum eigenvalue) of each layer to a certain value which in turn
lowers the Lipschitz constant of the network, preventing it to blow in unseen
regions. The influence of the Spectral Norm has been studied in the case of
Plug-and-Play approaches by Ryu et al. [2019].

###  A.2 Network training

We use the Adam optimizer for network training, with a learning rate of
$10^{-4}$. We add a white Gaussian noise of variance $\sigma_{s}^{2}$ to the
images scaled by a factor of $10^{6}$. $\sigma_{s}$ is drawn on-the-fly from a
Gaussian distribution of variance $s=50$. This means that at training time,
the standard deviation of the noise can be negative, following the
recommendation of Lim et al. [2020] to go from extrapolation to interpolation.

###  A.3 HMC procedure

The HMC procedure starts from a zero-filled reconstruction of the image with
added white Gaussian noise of variance $\sigma_{init}=100$. For the reduction
of noise standard deviation at reduction step, we use a factor $\gamma=0.995$,
following the recommendation of Song and Ermon [2020]. We take a step size
$\alpha$ dependent of the sampling temperature at step $i$,
$\alpha=\epsilon\left(\frac{\sigma_{i}}{\sigma_{0}}\right)^{1.5}$ with
$\epsilon=10$.

###  A.4 Zoom on Reconstruction

Figure 2:  Zoom of Bayesian posterior sampling for MRI reconstruction. The
ordering is the same as in Figure 1.

The zoom on this part of the reconstruction shown in Figure 2, illustrates a
failure of the Neural network reconstruction. It indeed produces a sharp
artifact that might hinder the examiner’s conclusions. The bayesian sampling
allows us to see that this region is not necessarily very sharp.

##  Appendix B Score-Based Hamiltonian Monte Carlo

We used HMC for sampling from the posterior distribution. The sampling
procedure follows the leapfrog integrator presented bellow:

$\begin{array}[]{lcl}\boldsymbol{m}_{t+\frac{\alpha}{2}}&=&\boldsymbol{m}_{t}+\dfrac{\alpha}{2}\nabla_{{\mathb{x}}}\log
p({\mathb{x}}_{t}|{\mathb{y}})\\\
{\mathb{x}}_{t+\alpha}&=&{\mathb{x}}_{t}+\alpha\boldsymbol{\mathrm{M}}^{-1}\boldsymbol{m}_{t+\frac{\alpha}{2}}\\\
\boldsymbol{m}_{t+\alpha}&=&\boldsymbol{m}_{t+\frac{\alpha}{2}}+\dfrac{\alpha}{2}\nabla_{{\mathb{x}}}\log
p({\mathb{x}}_{t+\alpha}|{\mathb{y}})\end{array}$ (6)

where $\alpha$ is the step size, $\boldsymbol{m}$ the momentum and
$\boldsymbol{\mathrm{M}}$ is a preconditioning matrix taking into account the
space metric, but in our case the identity matrix. The gradient of the log
posterior is computed as $\nabla_{{\mathb{x}}}\log
p({\mathb{x}}_{t}|{\mathb{y}})=\nabla_{{\mathb{x}}}\log
p({\mathb{y}}|{\mathb{x}}_{t})+\nabla_{{\mathb{x}}}\log
p({\mathb{y}}|{\mathb{x}}_{t})$, using the inverse Bayes theorem in equation
1.

It is worth noting that the discretized HMC with final step size $\epsilon$ is
not guaranteed to converge to the desired target distribution. This is the
motivation for the Metropolis-Hastings algorithm [Metropolis et al., 1953,
Hastings, 1970], which corrects for the discretization errors with a
correction step. Concretely, a new step is proposed using Equation 6, but it
is accepted with probability of the form
$\min\\{1,p({\mathb{x}}^{*})q({\mathb{x}}_{n}|{\mathb{x}}^{*})/p({\mathb{x}}_{n})q({\mathb{x}}^{*}|{\mathb{x}}_{n})\\}$
where ${\mathb{x}}_{n}$ is the last accepted point in the chain,
${\mathb{x}}^{*}$ is the proposed update, $p$ is the target density, and $q$
is the proposal corresponding to taking a step of Equation 6. In the context
of a score-only approach, we do not have access to $p$ and therefore cannot
compute directly this acceptance probability. We can nonetheless estimate from
the scores themselves by using a path integral:

$\log p({\mathb{x}}_{*})-\log
p({\mathb{x}}_{n})=\int_{0}^{1}\nabla_{{\mathb{x}}}\log
p(t*({\mathb{x}}_{*}-{\mathb{x}}_{n})+{\mathb{x}}_{n})\cdot({\mathb{x}}_{*}-{\mathb{x}}_{n})dt$
(7)

This integral can be evaluated numerically to arbitrary precision (at the cost
of additional score evaluations), but we find excellent agreement with the
theoretical value on toy model densities like the two-moons distribution with
a simple 4 points Simpson integration rule.

This illustrates that for gradient-based MCMC samplers, which also includes
Metropolis Adjusted Langevin Algorithm, these algorithms can fully be
implemented with only knowledge of the score.