# Unrolled Variational Bayesian Algorithm for Image Blind Deconvolution

Yunshi Huang, Emilie Chouzenoux, and Jean-Christophe Pesquet The authors
acknowledge support from the _Agence Nationale de la Recherche_ of France
under MAJIC (ANR-17-CE40-0004-01) project, and from the European Research
Council Starting Grant MAJORIS ERC-2019-STG-850925.

###### Abstract

In this paper, we introduce a variational Bayesian algorithm (VBA) for image
blind deconvolution. Our generic framework incorporates smoothness priors on
the unknown blur/image and possible affine constraints (e.g., sum to one) on
the blur kernel. One of our main contributions is the integration of VBA
within a neural network paradigm, following an unrolling methodology. The
proposed architecture is trained in a supervised fashion, which allows us to
optimally set two key hyperparameters of the VBA model and lead to further
improvements in terms of resulting visual quality. Various experiments
involving grayscale/color images and diverse kernel shapes, are performed. The
numerical examples illustrate the high performance of our approach when
compared to state-of-the-art techniques based on optimization, Bayesian
estimation, or deep learning.

###### Index Terms:

Variational Bayesian approach, Kullback-Leibler divergence, Majorization-
Minimization, blind deconvolution, image restoration, neural network,
unrolling, deep learning.

##  I Introduction

Image blind deconvolution problem appears in many fields of image processing
such as astronomy [1], biology [2] and medical imaging [3]. Given a degraded,
blurred and noisy image, the aim is to restore a clean image along with an
estimate of the blur kernel. Blind deconvolution is a severely ill-posed
problem as there exists an infinite number of pairs (image/blur) that lead to
the same observed image. Blind deconvolution methods adopt either a sequential
identification process [4], or a joint estimation approach [5]. In the former,
the blur kernel is identified first, possibly through a calibration step [6,
7, 8]. Then the unknown image is inferred using a non-blind image restoration
method. In the latter, the blur kernel and unknown image are simultaneously
estimated. Since the problem is highly ill-posed, it is mandatory to
incorporate prior knowledge on the sought unknowns. The retained prior
strongly influences the choice for the solver. Let us distinguish three main
classes of joint blind deconvolution approaches. A first option consists of
formulating the problem as the minimization of a cost function gathering a
data fidelity term (e.g., least-squares discrepancy) and penalties/constraints
acting on the image and kernel variables. In such a way, it is quite standard
to impose normalization and sparsity enhancing constraints on the kernel
coefficients to avoid scale ambiguity inherent to the blind deconvolution
model [9, 10, 11]. One can also easily impose the smoothness of the image, by
adopting total-variation based regularization [12]. Several other efficient
choices have been proposed in the literature, along with suitable iterative
optimization methods to solve the resulting problems [13, 14, 15, 16, 17]. The
main advantage of this family of methods is probably its flexibility. But this
comes at the price of heavy parameter tuning. The second option is to resort
to a Bayesian formulation to express the model and a priori knowledge on the
variables. The estimates are then defined from the estimation of the moments
(typically, the mean) of a posterior distribution given the observed data and
prior. As this typically involves the evaluation of intractable integrals,
sampling [18, 19] or approximation [20] strategies are used. Markov chain
Monte Carlo (MCMC) methods have been widely used for blind deconvolution
involving 1D sparse signals [21, 22, 23], but it is up to our knowledge
scarcely employed in large-scale problems [24], probably for computational
time reasons. Another family of approach consists in adopting the so-called
variational Bayesian approximation paradigm [25, 26]. Then, a simpler (usually
separable) approximation to the posterior is built through the minimization of
a suitable divergence. This approach leads to fast Bayesian-based algorithms,
whose great performance has been assessed in the context of non-blind [27, 28]
and blind [29] image restoration. Bayesian-based techniques usually require
less parameters than optimization-based ones. Moreover, they can provide
higher-order moments estimates, such as covariance matrices, which are of high
interest for assessing probabilistically the uncertainty of the results.
However, dealing with complex noise models and priors in such methods may be
tricky, and the algorithms may be quite computationally heavy. A recent trend
is to insert optimization-based steps in Bayesian sampling/approximation
methods for a more versatile modality and faster computations. See, for
example, [27, 30, 31] for applications of such ideas in the context of large-
scale image processing. The third category of methods is more recent, as the
references only go back to the last decade [32, 33, 34]. These methods rely on
the deep learning methodology. More precisely, a supervised learning strategy
is adopted to learn (implicitly) some prior information on the image/kernel
from a so-called training set. A highly non-linear and multi-layers
architecture is built, and its parameters (i.e., neuron weights) are estimated
by back-propagation to minimize a given loss function associated with the task
at hand (e.g., image visual quality). Several recent works propose neural
network architectures dedicated to the problem of image blind deconvolution.
Let us mention DeblurGAN [32], based on conditional generative adversarial
networks and a multi-component loss function, SRCNN [33] and its extended
version, DBSRCNN [34], relying on a CNN architecture, SelfDeblur [35]
combining an optimization-based method with two GAN networks. These methods
can reach very good performance, as long as the training set is large and
representative enough. Moreover, they are well suited to GPU-based
implementation. However, they have traditionally suffered from lack of
interpretability and robustness [36]. An emerging set of methods, in the field
of inverse problems in signal/image processing, performs algorithm unrolling
[37], where an iterative method (e.g., an optimization algorithm) is unrolled
as layers of a neural network. The reduced set of parameters of this network
are learnt by supervised training. Promising results have been obtained in the
context of image deconvolution in [38, 39, 40]. Theoretical results assessing
the stability and robustness of unrolling techniques can be found in [40, 41,
42]. These methods are also closely related to plug-and-play techniques where
a trained neural network is employed as the denoiser [40, 43].

In this paper, we propose a novel approach for blind image deconvolution, that
aims at gathering the best of the three aforementioned methods. We first
introduce a variational Bayesian algorithm (VBA) enhanced by optimization-
based ideas from [27], with the advantages to cope with a large set of priors
on the kernel and the image, and to present a reduced computational cost.
Then, we apply the unrolling paradigm to create a deep neural network
architecture, where VBA iterations are integrated as layers. This allows us to
(i) learn the hyperparameters (in particular, the noise level) of VBA in an
automatic supervised fashion, (ii) improve further the quality of the results
by choosing a dedicated loss in the training phase, (iii) implement the method
by taking full advantage of possible GPU resources, thus considerably reducing
the processing time during the test phase. In contrast to standard deep
learning methods for blind deconvolution, all these benefits come along with a
preservation of the interpretability of the method, thanks to the unrolling
technique. Let us emphasize that variational Bayesian methods often appear in
deep learning context. Indeed, they are backbones of variational autoencoders
[44] and also constitute methods of choice for training Bayesian neural
networks [45]. However, up to our knowledge, our work is the first to
investigate the unrolling of a variational Bayesian technique.

The rest of our paper is organized as follows. In Section II, we introduce the
image degradation model and introduce our Bayesian modeling, and provide the
background for deriving our algorithm. Section III explicitly describes the
iterative updates of the proposed VBA. The unrolling of VBA is presented in
Section IV. Numerical results, including comparisons with various methods, are
presented in Section V. Section VI concludes this paper.

##  II Problem statement

###  II-A Observation model

We focus on the restoration of an original image
$\widetilde{\textbf{x}}\in\mathbb{R}^{N}$, from a degraded version of it
$\textbf{y}\in\mathbb{R}^{N}$, related to $\widetilde{\textbf{x}}$ according
to the following model:

${\textbf{y}=\widetilde{\textbf{H}}\widetilde{\textbf{x}}+\textbf{n}.}$ (1)

Hereabove, $\textbf{n}\in\mathbb{R}^{N}$ models some additive random
perturbation on the observation. Moreover,
$\widetilde{\textbf{H}}\in\mathbb{R}^{N\times N}$ is a linear operator
modeling the effect of a blur kernel
$\widetilde{\textbf{h}}\in\mathbb{R}^{M}$. In this work, we focus on the
generalized blind deconvolution problem where the matrix associated with a
given kernel $\textbf{h}=[h_{1},\ldots,h_{M}]^{\top}$ reads

$\textbf{H}=\sum^{M}_{m=1}h_{m}\textbf{S}_{m},$ (2)

with $\\{\textbf{S}_{1},\ldots,\textbf{S}_{M}\\}$ is a set of $M$ known sparse
$N\times N$ real-valued matrices. This model allows to retrieve the standard
image deblurring model, as a special case when H identifies with a 2D discrete
convolution matrix with suitable padding. The considered problem amounts to
retrieving an estimate ($\widehat{\textbf{x}},\widehat{\textbf{h}}$) of the
pair of variables $(\widetilde{\textbf{x}},\widetilde{\textbf{h}})$ given y.
Due to the ill-posedness of this inverse problem, assumptions are required on
the sought image / kernel and on the noise statistics to reach satisfying
results. In the sequel, we will assume that n is a realization of an additive
Gaussian noise with zero mean and standard deviation $\sigma$. In the
remainder of the paper, it will be convenient to set $\beta=\sigma^{-2}$.
Furthermore, we introduce a linear equality constraint on the blur kernel
estimate h. A general expression of such a constraint is as follows:

$\textbf{h}=\textbf{T}\textbf{z}+\textbf{t},$ (3)

where $\textbf{T}=(T_{m,p})_{1\leq m\leq M,1\leq p\leq
P}\in\mathbb{R}^{M\times P}$ is a matrix of rank $P\in\\{1,\ldots,M\\}$ and
$\textbf{t}=[t_{1},\dots,t_{M}]^{\top}\in\mathbb{R}^{M}$ is some vector of
predefined constants. Vector
$\textbf{z}=[z_{1},\ldots,z_{P}]^{\top}\in\mathbb{R}^{P}$ becomes the new
unknown of the problem, along with the image x. A typical linear equality
constraint in such context is the sum-to-one constraint, i.e.
$\sum_{m=1}^{M}h_{m}=1$. Other examples will be provided in the experimental
section. We can thus rewrite (2) as

$\textbf{H}=\sum^{P}_{p=1}z_{p}\textbf{K}_{p}+\textbf{K}_{0}=\mathcal{H}(\textbf{z}),$
(4)

with

$(\forall
p\in\\{1,\ldots,P\\})\quad\textbf{K}_{p}=\sum_{m=1}^{M}T_{m,p}\textbf{S}_{m}\in\mathbb{R}^{N\times
N},$ (5)

and

$\textbf{K}_{0}=\sum_{m=1}^{M}t_{m}\textbf{S}_{m}\in\mathbb{R}^{N\times N}.$
(6)

###  II-B Hierarchical Bayesian Modeling

Let us now introduce the hierarchical Bayesian model on which our VBA method
will be grounded.

####  II-B1 Likelihood

First, we express the likelihood $p(\textbf{y}|\textbf{x},\textbf{z})$ of the
observed data, given the unknowns $(\textbf{x},\textbf{z})$. Since the noise
is assumed to be Gaussian distributed, the likelihood can be expressed as
follows:

$\displaystyle p(\textbf{y}|\textbf{x},\textbf{z})$
$\displaystyle=\beta^{\frac{N}{2}}\textrm{exp}\left(-\frac{\beta}{2}||\textbf{y}-\mathcal{H}(\textbf{z})\textbf{x}||^{2}\right),$
(7)

where we recall that $\beta$ denotes the inverse of the noise variance.

####  II-B2 Prior

As already mentioned, it is necessary to incorporate suitable prior knowledge
on the sought quantities to limit the problem ill-posedness. We here consider
a wide range of sparsity enhancing prior for the image x, by adopting the
generic model,

$p(\textbf{x}|\gamma)\propto\gamma^{\frac{N}{2\kappa}}\exp\Big{(}-\gamma\sum_{j=1}^{J}||\textbf{D}_{j}\textbf{x}||^{2\kappa}\Big{)},$
(8)

with $\kappa\in(0,1]$ a scale parameter and $(\textbf{D}_{j})_{1\leq j\leq
J}\in(\mathbb{R}^{S\times N})^{J}$ both assumed to be known. For instance, an
isotropic total variation prior is obtained by setting $\kappa=1/2$, $S=2$,
$J=N$ and for every $j\in\\{1,\ldots,N\\}$,
$\textbf{D}_{j}\textbf{x}=[[\nabla^{h}\textbf{x}]_{j},[\nabla^{v}{\textbf{x}}]_{j}]\in\mathbb{R}^{2}$
gathers the horizontal and vertical gradients of x at pixel $j$. Other
relevant choices are discussed in [27]. Hereabove, $\gamma>0$ is a
regularization hyperparameter that we incorporate in our hierarchical model.
We assume a Gamma distribution on $\gamma$,

$p(\gamma)\propto\gamma^{\alpha-1}\exp(-\eta\gamma),$ (9)

where $\alpha\geq 0$ and $\eta\geq 0$ are the (known) shape and inverse scale
parameters of the Gamma distribution. Such choice for the hyperparameter is
rather standard in the context of Bayesian image restoration.
Regarding the blur h, we adopt the so-called SAR model, successfully used for
Bayesian-based blind deconvolution in [29]. The model relies on the following
Gaussian model,

$p(\textbf{h}|\xi)\propto\xi^{\frac{M}{2}}\exp\Big{(}-\frac{\xi}{2}||\textbf{A}(\textbf{h}-\textbf{m})||^{2}\Big{)},$
(10)

where $\textbf{A}\in\mathbb{R}^{Q\times M}$ with
$Q\in\mathbb{N}\setminus\\{0\\}$ denotes a matrix of rank $M$.
$\textbf{m}\in\mathbb{R}^{M}$ is the mean of the underlying Gaussian
distribution, and $\xi>0$ is such that $\xi\textbf{A}^{\top}\textbf{A}$ is its
inverse covariance matrix. If h follows this distribution, the projection of h
onto the affine space defined by (3) is also Gaussian as well as the vector z
associated with each projected vector. More precisely, z follows a Gaussian
distibution with mean $\bm{\mu}=\textbf{T}^{-1}(\textbf{m}-\textbf{t})$ and
covariance matrix
$\xi^{-1}\textbf{T}^{-1}(\textbf{A}^{\top}\textbf{A})^{-1}(\textbf{T}^{-1})^{\top}$
where $\textbf{T}^{-1}$ is the left inverse of T, i.e.
$\textbf{T}^{-1}=(\textbf{T}^{\top}\textbf{T})^{-1}\textbf{T}^{\top}$. This
yields the following prior for the variable of interest z:

$p(\textbf{z}|\xi)\propto\xi^{\frac{P}{2}}\exp\Big{(}-\frac{\xi}{2}(\textbf{z}-\bm{\mu})^{\top}\textbf{L}(\textbf{z}-\bm{\mu})\Big{)},$
(11)

where
$\textbf{L}=\textbf{T}^{\top}\textbf{T}(\textbf{T}^{\top}\big{(}\textbf{A}^{\top}\textbf{A})^{-1}\textbf{T}\big{)}^{-1}\textbf{T}^{\top}\textbf{T}$.
We will consider $(\textbf{L},\bm{\mu})$ to be predefined by the user, so as
to be adapted to the sought properties of the blur kernel to estimate. The
hyperparameter $\xi$ will be learned during a training phase, as we will
explain in Section IV.

####  II-B3 Hierarchical model

Let us assume that $(\textbf{x},\gamma)$ and z are mutually independent.
According to Bayes formula, the posterior distribution of the unknowns
$\Theta=(\textbf{x},\textbf{z},\gamma)$ given the observed data y is defined
as

$p(\Theta|\textbf{y})\propto
p(\textbf{y}|\textbf{x},\textbf{z})p(\textbf{x}|\gamma)p(\textbf{z}|\xi)p(\gamma),$
(12)

where the four factors on the right side have been defined above.

###  II-C Variational Bayesian Inference

The Bayesian inference paradigm seeks for solving the blind restoration
problem through the exploration of the posterior $p(\Theta|\textbf{y})$.
Typically, one would be interested in the posterior mean, its covariance, or
its modes (i.e., maxima). Let us make (12) explicit:

$\displaystyle
p(\Theta|\textbf{y})\propto\textrm{exp}\left(-\gamma\sum_{j=1}^{J}||\textbf{D}_{j}\textbf{x}||^{2\kappa}-\frac{\beta}{2}||\textbf{y}-\mathcal{H}(\textbf{z})\textbf{x}||^{2}\right)$
$\displaystyle\times\gamma^{\frac{N}{2\kappa}+\alpha-1}\exp(-\eta\gamma)\xi^{\frac{P}{2}}\exp\Big{(}-\frac{\xi}{2}(\textbf{z}-\bm{\mu})^{\top}\textbf{L}(\textbf{z}-\bm{\mu})\Big{)}.$
(13)

Unfortunately, neither $p(\Theta|\textbf{y})$, nor its moments (e.g., mean,
covariance), nor its mode positions have a closed form. In particular
$p(\textbf{y})$, which acts as a normalization constant, cannot be calculated
analytically. We thus resort to the variational Bayesian framework to
approximate this distribution by a more tractable one, denoted by $q(\Theta)$,
for which the estimators are easier to compute. The approximation is computed
with the aim to minimize the Kullback-Leibler (KL) divergence between the
target posterior and its approximation, which amounts to determining

$\displaystyle q^{\textrm{opt}}(\Theta)$ $\displaystyle=\textrm{argmin}_{q}\
\mathcal{KL}(q(\Theta)||p(\Theta|\textbf{y})),$
$\displaystyle=\textrm{argmin}_{q}\ \int
q(\Theta)\textrm{ln}\left(\frac{q(\Theta)}{p(\Theta|\textbf{y})}\right)\textrm{d}\Theta,$
(14)

where the equality holds only when $q(\Theta)=p(\Theta|\textbf{y})$. In order
to make the solution of the above minimization problem tractable, a typical
strategy is to make use of a variational Bayesian algorithm (VBA) based on a
so-called mean field approximation of the posterior, combined with an
alternating minimization procedure.

The mean field approximation reads as a factorized structure
$q(\Theta)=\prod_{r=1}^{R}q_{r}(\Theta_{r})$, which is assumed for the
distribution $q$. Each of the $R$ factors are then obtained by minimizing the
KL divergence by iterative update of a given factor $q_{r}$ while holding the
others unchanged. This procedure takes advantage of the property that the
minimizer of the KL divergence with respect to each factor can be expressed as

$\displaystyle(\forall r\in\\{1,...,R\\})\quad
q_{r}^{\textrm{opt}}(\Theta_{r})\propto\exp\left(<\textrm{ln}p(\textbf{y},\Theta)>_{\prod_{i\neq
r}q_{i}^{\textrm{opt}}(\Theta_{i})}\right)$ (15)

where $<\cdot>_{\prod_{i\neq r}q_{i}(\Theta_{i})}=\int\,\prod_{i\neq
r}q_{i}(\Theta_{i})\textrm{d}\Theta_{i}$. Here, we will consider the following
factorization:

$q(\Theta)=q_{\textbf{X}}(\textbf{x})q_{\textbf{Z}}(\textbf{z})q_{\Gamma}(\gamma).$
(16)

In Section III, we describe the steps of VBA for this particular choice. Due
to the intricate form of the chosen prior on the image, we introduce an extra
approximation step, relying on a majoration-minimization (MM) strategy,
reminescent from [27]. In addition, we propose a strategy to reduce the time
complexity of VBA, so as to deal with medium to large size images. As we will
emphasize, the method requires the setting of two cumbersome hyperparameters,
namely the regularization weight $\xi$ and the noise level $\beta$. Then, in
Section IV, we show how to unroll the VBA method as a neural network
structure, so as to learn the parameters $(\xi,\beta)$ in a supervised
fashion.

##  III VBA for blind image deconvolution

We now describe our proposed implementation of the VBA when applied to the
approximation to the posterior in (13). We first present an MM-based procedure
to handle the complicated form of the prior term on variable x. Then, we give
the explicit expressions of the updates performed in the alternating
minimization method.

###  III-A MM-based approximation

Let us focus on the prior term in (8). This distribution is difficult to deal
with as soon as $\kappa$ is different from $1$ (in which case a Gaussian
distribution is retrieved). We thus propose to construct a surrogate for the
prior on x. We use the tangent inequality for concave functions, which yields
the following majorant function for the $\ell_{\kappa}$-function with
$\kappa\in(0,1]$:

$(\forall u>0)(\forall v\geq 0)\quad v^{\kappa}\leq(1-\kappa)u^{\kappa}+\kappa
u^{\kappa-1}v.$ (17)

Let us introduce the vector of auxiliary positive variables
$\bm{\lambda}=(\lambda_{j})_{1\leq j\leq J}$. From the previous inequality, we
then deduce the following majorant function for the negative logarithm of the
prior distribution:

$(\forall\textbf{x}\in\mathbb{R}^{N})\quad\gamma\sum_{j=1}^{J}||\textbf{D}_{j}\textbf{x}||^{2\kappa}\leq\sum_{j=1}^{J}F_{j}(\textbf{D}_{j}\textbf{x},\lambda_{j};\gamma),$
(18)

where, for every $j\in\\{1,\ldots,J\\}$,

$F_{j}(\textbf{D}_{j}\textbf{x},\lambda_{j};\gamma)=\gamma\frac{\kappa||\textbf{D}_{j}\textbf{x}||^{2}+(1-\kappa)\lambda_{j}}{\lambda_{j}^{1-\kappa}}.$
(19)

This majorant function can be understood as a Gaussian lower bound on the
prior distribution on x, which will appear more tractable in the VBA
implementation. We will also show that the update of the auxiliary variables
remains rather simple, thus not impacting the complexity of the whole
procedure.

In a nutshell, using (13), and (18), we obtain the following inequality:

$p(\Theta|\textbf{y})\geq\mathcal{F}(\Theta|\textbf{y};\bm{\lambda})$ (20)

where the lower bound on the posterior distribution is

$\displaystyle\mathcal{F}(\Theta|\textbf{y};\bm{\lambda})=$ $\displaystyle
C\gamma^{\frac{N}{2\kappa}}\exp\left(-\frac{\beta}{2}||\textbf{y}-\mathcal{H}(\textbf{z})\textbf{x}||^{2}-F(\textbf{x},\bm{\lambda};\gamma)\right)p(\gamma)p(\textbf{z}|\xi).$
(21)

Hereabove we have introduced the shorter notation

$F(\textbf{x},\bm{\lambda};\gamma)=\sum^{J}_{j=1}F_{j}(\textbf{D}_{j}\textbf{x},\lambda_{j};\gamma)$
(22)

and $C$ is a multiplicative constant independent from $\Theta$. Inequality
(20) leads to the following majorization of the $\mathcal{KL}$ divergence
involved in (14):

$\mathcal{KL}(q(\Theta)||p(\Theta|\textbf{y}))\leq\mathcal{KL}(q(\Theta)||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda})).$
(23)

By minimizing the upper bound in (23) with respect to $\bm{\lambda}$, we can
keep it as tight as possible, so as to guarantee the good performance of the
VBA. To summarize, we propose to solve Problem (14) through the following four
iterative steps:

  1. 1.

Minimizing
$\mathcal{KL}(q(\Theta)||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}))$ w.r.t.
$q_{\textbf{X}}(\textbf{x})$.

  2. 2.

Minimizing the upper bound
$\mathcal{KL}(q(\Theta)||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}))$ in (23)
w.r.t. $q_{\textbf{Z}}(\textbf{z})$.

  3. 3.

Update the auxiliary variables $(\lambda_{j})_{1\leq j\leq J}$ to minimize
$\mathcal{KL}(q(\Theta)||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}))$.

  4. 4.

Minimizing
$\mathcal{KL}(q(\Theta)||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}))$ w.r.t.
$q_{\Gamma}(\gamma)$.

Subsequently, at a given iteration $k$ of the proposed algorithm, the
corresponding estimated variables will be indexed by $k$.

###  III-B VBA updates

Let us now describe the four steps of the proposed VBA, starting from a given
iteration $k$ associated with the current approximated distributions
$q^{k}_{\textbf{X}}(\textbf{x}),q^{k}_{\textbf{Z}}(\textbf{z})$, and
$q^{k}_{\Gamma}(\gamma)$, and the auxiliary parameter estimate
$\bm{\lambda}^{k}$. We also denote by
$(\textbf{x}^{k},\textbf{z}^{k},\gamma^{k})$ the estimates of the means of
$q^{k}_{\textbf{X}}$, $q^{k}_{\textbf{Z}}$, and $q^{k}_{\Gamma}$, and
$(\textbf{C}_{\textbf{x}}^{k},\textbf{C}_{\textbf{z}}^{k})$ the covariance
estimates for $q^{k}_{\textbf{X}}$ and $q^{k}_{\textbf{Z}}$.

####  III-B1 Update of $q_{\textbf{X}}(\textbf{x})$

By definition,

$q_{\textbf{X}}^{k+1}(\textbf{x})\\\
=\textrm{argmin}_{q_{\textbf{X}}}\mathcal{KL}(q_{\textbf{X}}(\textbf{x})q^{k}_{\Gamma}(\gamma)q^{k}_{\textbf{Z}}(\textbf{z})||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}^{k})).$
(24)

The standard solution provided by (15) remains valid, by replacing the joint
distribution by a lower bound chosen proportional to
$\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}^{k})$:

$\displaystyle q_{\textbf{x}}^{k+1}(\textbf{x})$
$\displaystyle\propto\exp\left(<\textrm{ln}\,\mathcal{F}(\textbf{x},\textbf{z},\gamma\mid\textbf{y};\bm{\lambda}^{k})>_{q^{k}_{\Gamma}(\gamma),q^{k}_{\textbf{Z}}(\textbf{z})}\right)$
$\displaystyle\propto\exp\left(\int\int\textrm{ln}\,\mathcal{F}(\textbf{x},\textbf{z},\gamma\mid\textbf{y};\bm{\lambda}^{k})q^{k}_{\Gamma}(\gamma)q^{k}_{\textbf{Z}}(\textbf{z})\textrm{d}\gamma\textrm{d}\textbf{z}\right).$
(25)

By decomposing the different terms and using (4),

$\displaystyle q_{\textbf{x}}^{k+1}(\textbf{x})\propto$
$\displaystyle\exp\Biggl{\\{}-\frac{1}{2}\textbf{x}^{\top}\Biggr{(}\beta\Big{(}\mathbb{E}_{q_{\textbf{Z}}^{k}(\textbf{z})}(\textbf{H})^{\top}\mathbb{E}_{q_{\textbf{Z}}^{k}(\textbf{z})}(\textbf{H})$
$\displaystyle+\sum^{P}_{p=1}\sum^{P}_{q=1}\textbf{e}^{\top}_{p}\textrm{cov}_{q^{k}_{\textbf{Z}}(\textbf{z})}(\textbf{z})\textbf{e}_{q}\textbf{K}^{\top}_{p}\textbf{K}_{q}\Big{)}$
$\displaystyle+2\mathbb{E}_{q_{\Gamma}^{k}(\gamma)}(\gamma)\textbf{D}^{\top}\Lambda^{k}\textbf{D}\Biggr{)}\textbf{x}+\beta\textbf{x}^{\top}\mathbb{E}_{q_{\textbf{Z}}^{k}(\textbf{z})}(\textbf{H})^{\top}\textbf{y}\Biggl{\\}}$
(26)

where

$\displaystyle\mathbb{E}_{q_{\textbf{Z}}^{k}(\textbf{z})}(\textbf{H})$
$\displaystyle=\sum^{P}_{p=1}\textbf{e}^{\top}_{p}\mathbb{E}_{q^{k}_{\textbf{Z}}(\textbf{z})}(\textbf{z})\textbf{K}_{p}+\textbf{K}_{0},$
(27) D
$\displaystyle=[\textbf{D}_{1}^{\top},\ldots,\textbf{D}_{J}^{\top}]^{\top},$
(28)

$\Lambda^{k}$ is the block diagonal matrix whose diagonal elements are
$(\kappa(\lambda_{j}^{k})^{\kappa-1}\textbf{I}_{S})_{1\leq j\leq J}$, and
$(\textbf{e}_{1},\ldots,\textbf{e}_{P})$ is the canonical basis of
$\mathbb{R}^{P}$. We thus obtain a Gaussian distribution:

$q_{\textbf{X}}^{k+1}(\textbf{x})=\mathcal{N}(\textbf{x};\check{\textbf{x}}^{k+1},\check{\textbf{C}_{\textbf{x}}}^{k+1}),$
(29)

parametrized by

$\displaystyle{(\check{\textbf{C}_{\textbf{x}}}^{k+1})^{-1}}=\;$
$\displaystyle\beta\Biggr{(}(\textbf{H}^{k})^{\top}\textbf{H}^{k}+\sum^{P}_{p=1}\sum^{P}_{q=1}\textbf{e}^{\top}_{p}{\textbf{C}_{\textbf{z}}^{k}}\textbf{e}_{q}\textbf{K}^{\top}_{p}\textbf{K}_{q}\Biggr{)}$
$\displaystyle+2\bm{\gamma}^{k}\textbf{D}^{\textrm{T}}\bm{\Lambda}^{k}\textbf{D},$
(30) $\displaystyle\check{\textbf{x}}^{k+1}=\;$
$\displaystyle\beta{\check{\textbf{C}_{\textbf{x}}}^{k+1}}(\textbf{H}^{k})^{\top}\textbf{y},$
(31)

with $\textbf{H}^{k}=\mathcal{H}(\textbf{z}^{k})$.

In image restoration applications, dimension $N$ can be rather large
(typically greater than $10^{6}$ variables), so that the storage of the full
covariance matrix $\check{\textbf{C}_{\textbf{x}}}^{k+1}$ is neither desirable
nor usually possible. We thus propose to resort to a diagonal approximation to
this matrix when required, so that the update finally reads:

$q_{\textbf{X}}^{k+1}(\textbf{x})=\mathcal{N}(\textbf{x};{\textbf{x}^{k+1}},\textbf{C}_{\textbf{x}}^{k+1}),$
(32)

with

$\displaystyle\textbf{C}_{\textbf{x}}^{k+1}$
$\displaystyle=\textrm{Diag}\left(\delta_{\textbf{x}}^{k+1}\right)$ (33)
$\displaystyle{\textbf{x}^{k+1}}$
$\displaystyle=\textrm{CG}\Big{(}{(\check{\textbf{C}_{\textbf{x}}}^{k+1})^{-1}},\beta(\textbf{H}^{k})^{\top}\textbf{y}\Big{)},$
(34)

where $\delta_{\textbf{x}}^{k+1}\in\mathbb{R}^{N}$ is the vector of the
inverses of the diagonal elements of
$(\check{\textbf{C}_{\textbf{x}}}^{k+1})^{-1}$, and
$\textrm{CG}(\texttt{A,b})$ denotes the application of a linear conjugate
gradient solver to the linear system Ax=b.

####  III-B2 Update of $q_{\textbf{Z}}(\textbf{z})$

According to the VBA principle,

$q_{\textbf{Z}}^{k+1}(\textbf{z})\\\
=\textrm{argmin}_{q_{\textbf{Z}}}\mathcal{KL}(q^{k+1}_{\textbf{X}}(\textbf{x})q^{k}_{\Gamma}(\gamma)q_{\textbf{Z}}(\textbf{z})||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}^{k})).$
(35)

Using (15) and the previously introduced bound
$\mathcal{F}(\Theta\mid\textbf{y};\bm{\lambda}^{k})$, we have

$\displaystyle q_{\textbf{Z}}^{k+1}(\textbf{z})$
$\displaystyle\propto\exp\left(\int\int\textrm{ln}\,\mathcal{F}(\textbf{x},\textbf{z},\gamma\mid\textbf{y};\bm{\lambda}^{k})q^{k}_{\Gamma}(\gamma)q^{k+1}_{\textbf{X}}(\textbf{x})\textrm{d}\gamma\textrm{d}\textbf{x}\right).$
(36)

Replacing the involved quantities by their expression yields

$q_{\textbf{Z}}^{k+1}(\textbf{z})\propto\exp\Biggl{\\{}-\frac{1}{2}\textbf{z}^{\top}\left(\beta\textbf{B}^{k+1}+\xi\textbf{L}\right)\textbf{z}\\\
+\textbf{z}^{\top}\left(\beta\textbf{a}^{k+1}+\xi\textbf{L}\bm{\mu}\right)\Biggl{\\}},$
(37)

where $\textbf{a}^{k+1}=(a^{k+1}_{p})_{1\leq p\leq P}\in\mathbb{R}^{P}$ and
$\textbf{B}^{k+1}=(B^{k+1}_{p,q})_{1\leq p,q\leq P}\in\mathbb{R}^{P\times P}$
are such that, for every $(p,q)\in\\{1,\ldots,P\\}^{2}$,

$\displaystyle a^{k+1}_{p}=$
$\displaystyle\mathbb{E}_{q_{\textbf{X}}^{k+1}(\textbf{x})}(\textbf{x})^{\top}\textbf{K}^{\top}_{p}\textbf{y}-\mathbb{E}_{q_{\textbf{X}}^{k+1}}(\textbf{x}^{\top}\textbf{K}^{\top}_{p}\textbf{K}_{0}\textbf{x})$
$\displaystyle=(\textbf{x}^{k+1})^{\top}\textbf{K}^{\top}_{p}\textbf{y}-\textbf{B}^{k+1}_{p,0},$
(38) $\displaystyle B^{k+1}_{p,q}$
$\displaystyle=\mathbb{E}_{q_{\textbf{X}}^{k+1}}(\textbf{x}^{\top}\textbf{K}^{\top}_{p}\textbf{K}_{q}\textbf{x})$
$\displaystyle=\textrm{trace}\left(\textbf{K}_{p}\textbf{C}_{\textbf{x}}^{k+1}\textbf{K}^{\top}_{q}\right)+({\textbf{x}^{k+1}})^{\top}\textbf{K}^{\top}_{p}\textbf{K}_{q}\textbf{x}^{k+1}$
(39)

with

$\displaystyle B^{k+1}_{p,0}$
$\displaystyle=\mathbb{E}_{q_{\textbf{X}}^{k+1}}(\textbf{x}^{\top}\textbf{K}^{\top}_{p}\textbf{K}_{0}\textbf{x})$
$\displaystyle=\textrm{trace}\left(\textbf{K}_{p}\textbf{C}_{\textbf{x}}^{k+1}\textbf{K}^{\top}_{0}\right)+({\textbf{x}^{k+1}})^{\top}\textbf{K}^{\top}_{p}\textbf{K}_{0}\textbf{x}^{k+1}.$
(40)

Thus, the update for the distribution $q_{\textbf{Z}}$ reads

$q^{k+1}_{\textbf{Z}}(\textbf{z})=\mathcal{N}(\textbf{z};{\textbf{z}^{k+1}},\textbf{C}_{\textbf{z}}^{k+1}),$
(41)

with

$\displaystyle{(\textbf{C}_{\textbf{z}}^{k+1})^{-1}}$
$\displaystyle=\beta\textbf{B}^{k+1}+\xi\textbf{L},$ (42)
$\displaystyle{\textbf{z}^{k+1}}$
$\displaystyle={\textbf{C}_{\textbf{z}}^{k+1}}\left(\beta\textbf{a}^{k+1}+\xi\textbf{L}\bm{\mu}\right).$
(43)

####  III-B3 Update of $\bm{\lambda}$

Let us now express the update of the auxiliary variable. We aim at finding

$\bm{\lambda}^{k+1}\\\
=\textrm{argmin}_{\bm{\lambda}}\mathcal{KL}(q^{k+1}_{\textbf{X}}(\textbf{x})q^{k}_{\Gamma}(\gamma)q^{k+1}_{\textbf{Z}}(\textbf{z})||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda})).$
(44)

This amounts to finding, for every $j\in\\{1,\ldots,J\\}$,

$\displaystyle\lambda^{k+1}_{j}$
$\displaystyle=\textrm{argmin}_{\lambda_{j}\in[0,+\infty)}\int
q^{k+1}_{\textbf{X}}(\textbf{x})q^{k}_{\Gamma}(\gamma)q^{k+1}_{\textbf{Z}}(\textbf{z})$
$\displaystyle\qquad\quad\times\textrm{log}\frac{q^{k+1}_{\textbf{X}}(\textbf{x})q^{k}_{\Gamma}(\gamma)q^{k+1}_{\textbf{Z}}(\textbf{z})}{\mathcal{F}(\Theta|\textbf{y},\bm{\lambda})}\textrm{d}\Theta,$
$\displaystyle=\textrm{argmin}_{\lambda_{j}\in[0,+\infty)}\sum^{J}_{j=1}\int\int
q^{k+1}_{\textbf{X}}(\textbf{x})q^{k}_{\Gamma}(\gamma)$
$\displaystyle\qquad\quad\times
F_{j}(\textbf{D}_{j}\textbf{x},\lambda_{j};\gamma)\textrm{d}\textbf{x}\textrm{d}\gamma,$
$\displaystyle=\textrm{argmin}_{\lambda_{j}\in[0,+\infty)}\frac{\kappa\mathbb{E}_{q^{k+1}_{\textbf{x}}(\textbf{x})}\left[||\textbf{D}_{j}\textbf{x}||^{2}\right]+(1-\kappa)\lambda_{j}}{\lambda_{j}^{1-\kappa}}.$
(45)

The explicit solution to the above minimization problem yields the following
update:

$\displaystyle\lambda^{k+1}_{j}$
$\displaystyle=\mathbb{E}_{q^{k+1}_{\textbf{x}}(\textbf{x})}\left[||\textbf{D}_{j}\textbf{x}||^{2}\right]$
$\displaystyle=||\textbf{D}_{j}\textbf{x}^{k+1}||^{2}+\textrm{trace}\left(\textbf{D}^{\textrm{T}}_{j}\textbf{D}_{j}\textbf{C}_{\textbf{x}}^{k+1}\right).$
(46)

####  III-B4 Update of $q_{\Gamma}(\gamma)$

Finally, the update related to the hyperparameter $\gamma$ is expressed as

$q_{\Gamma}^{k+1}(\gamma)\\\
=\textrm{argmin}_{q_{\Gamma}}\mathcal{KL}(q^{k+1}_{\textbf{X}}(\textbf{x})q_{\Gamma}(\gamma)q^{k+1}_{\textbf{Z}}(\textbf{z})||\mathcal{F}(\Theta|\textbf{y};\bm{\lambda}^{k+1})).$
(47)

Using (15), we have

$q_{\Gamma}^{k+1}(\gamma)\propto\exp\Biggr{(}\int\int\textrm{ln}\,\mathcal{F}(\textbf{x},\textbf{z},\gamma\mid\textbf{y};\bm{\lambda}^{k+1})\\\
\times
q^{k+1}_{\textbf{X}}(\textbf{x})q^{k+1}_{\textbf{Z}}(\textbf{z})\,\textrm{d}\textbf{x}\textrm{d}\textbf{z}\Biggr{)}.$
(48)

The above integral has the following closed form expression:

$\displaystyle q_{\Gamma}^{k+1}(\gamma)$
$\displaystyle\propto\gamma^{\frac{N}{2\kappa}+\alpha-1}\exp(-\eta\gamma)$
$\displaystyle\times\exp\left(-\gamma\sum^{J}_{j=1}\frac{\kappa\,\mathbb{E}_{q^{k+1}_{\textbf{x}}(\textbf{x})}\left[||\textbf{D}_{j}\textbf{x}||^{2}\right]+(1-\kappa)\lambda^{k+1}_{j}}{(\lambda^{k+1}_{j})^{1-\kappa}}\right).$
(49)

It thus follows from (46) that the update of $q_{\Gamma}$ is

$\displaystyle q_{\Gamma}^{k+1}(\gamma)=\Gamma(d,b^{k+1}),$ (50)

that is the Gamma distribution with parameters

$\displaystyle d=\frac{N}{2\kappa}+\alpha,\ \
b^{k+1}=\sum^{J}_{j=1}(\lambda^{k+1}_{j})^{\kappa}+\eta.$ (51)

The mean of $q_{\Gamma}^{k+1}$ is finally given by

$\gamma^{k+1}=\frac{d}{b^{k+1}}.$ (52)

Note that parameter $d$ is not iteration dependent and can thus be precomputed
from the beginning of the VBA.

###  III-C Overview of VBA

Algorithm 1 provides a summary of the resulting VBA for solving the blind
deconvolution problem introduced in Section II. We also specify our
initialization strategy. More practical details about the latter will be
discussed in the experimental section. As a result, the optimal posterior
distributions for both variables x and z will be approximated as Gaussian
distributions, while the one for hyperparameter $\gamma$ is approximated by a
Gamma distribution. In particular, after $K$ iterations, it is direct to
extract from VBA outputs an estimate for the posterior mean of the image and
the kernel, through variable $\textbf{x}^{K}$ and
$\textbf{T}\textbf{z}^{K}+\textbf{t}$. The associated covariance matrices are
given by $\textbf{C}_{\textbf{x}}^{K}$ and
$\textbf{T}\textbf{C}_{\textbf{z}}^{K}\textbf{T}^{\top}$. These matrices can
be useful to perform uncertainty quantification of the results. The VBA also
allows us to estimate easily the hyperparameter $\gamma$ involved in the image
prior. Nonetheless, it appears difficult to find an efficient manner to
estimate the hyperparameter $\xi$ using a variational Bayesian approach, as
this value highly fluctuates from one image/kernel pair to the other so that a
simple prior modeling of does not appear obvious. Moreover, the VBA requires
the knowledge of the noise level, through the parameter $\beta$. This is
limitating, and one might prefer to have this quantity estimated in an
automatic manner. Thus, we propose in the next section, to resort to a
supervised learning strategy to learn both $\xi$ and $\beta$ along the
iterates of VBA, in the spirit of recent works [39] on the unrolling (also
called unfolding) of iterative algorithms.

Algorithm 1 VBA approach for image blind deconvolution

1:Set hyperparameters $(\xi,\beta,\alpha,\eta)$. Define initial values for
$(\textbf{x}^{0},\textbf{C}_{\textbf{x}}^{0},\textbf{z}^{0},\textbf{C}_{\textbf{z}}^{0})$.
Compute $\bm{\lambda}^{0}$ and $\gamma^{0}$ using (46) and (52), respectively.

2:For $k=0,1,\ldots,K$:

3:Update the mean $\textbf{x}^{k+1}$ and the covariance matrix
$\textbf{C}_{\textbf{x}}^{k+1}$ of $q_{\textbf{X}}^{k+1}(\textbf{x})$ using
(33)-(34).

4:Update the mean $\textbf{z}^{k+1}$ and the covariance matrix
$\textbf{C}_{\textbf{z}}^{k+1}$ of $q_{\textbf{Z}}^{k+1}(\textbf{z})$ using
(42)-(43).

5:Update $\lambda^{k+1}_{j}$ using (46), for every $j\in\\{1,\ldots,J\\}$.

6:Update the mean $\gamma^{k+1}$ of $q_{\Gamma}^{k+1}(\gamma)$ using
(51)-(52).

##  IV Supervised learning of VBA hyperparameters

###  IV-A Overview

We introduce a supervised learning strategy to estimate the hyperparameter
$\xi$ and the inverse of the noise variance $\beta$, that are required to run
VBA. We adopt the so-called _unrolling_ (or unfolding) methodology [37]. The
idea is to view each iteration of an iterative algorithm as one layer of a
neural network structure. Each layer can be parametrized by some quantities
that are learned from a training database so as to minimize a task-oriented
loss function. The advantage of the unrolling approach is threefold: (i) each
layer mimics one iteration of the algorithm and thus it is highly
interpretable, (ii) the choice of the loss is directly related to the task at
the end, which is beneficial for the quality of the results, (iii) once
trained, the network can be applied easily and rapidly on a large set of test
data without any further tuning. In particular, its implementation can make
use of GPU-accelerated frameworks. Several recent examples in the field of
image processing have shown the benefits of unrolling [46, 47, 48, 49] when
compared to standard black-box deep learning techniques or more classical
restoration methods based on Bayesian or optimization tools. Let us in
particular mention the works [38, 50] for the application of unrolling in the
context of blind image restoration.

Let us now specify the unrolling procedure in the context of VBA. Let $K>0$ be
the number of iterations of the VBA described in Algorithm 1, thus
corresponding to $K$ layers of a neural network architecture. Iteration
$k\in\\{0,\ldots,K-1\\}$ of our unrolled VBA can be conceptually expressed as

$\displaystyle(\textbf{x}^{k+1},\textbf{C}_{\textbf{x}}^{k+1},\textbf{z}^{k+1},\textbf{C}_{\textbf{z}}^{k+1},\bm{\lambda}^{k+1},\bm{\gamma}^{k+1})$
$\displaystyle=\mathcal{A}(\textbf{x}^{k},\textbf{C}_{\textbf{x}}^{k},\textbf{z}^{k},\textbf{C}_{\textbf{z}}^{k},\bm{\lambda}^{k},\bm{\gamma}^{k},\xi^{k},\beta^{k}).$
(53)

The initialization procedure for
$(\textbf{x}^{0},\textbf{C}_{\textbf{x}}^{0},\textbf{z}^{0},\textbf{C}_{\textbf{z}}^{0},\bm{\lambda}^{0},\bm{\gamma}^{0})$
is detailed in Algorithm 1. For $k\in\\{0,\ldots,K-1\\}$, the expressions of
$(\textbf{x}^{k+1},\textbf{C}_{\textbf{x}}^{k+1},\textbf{z}^{k+1},\textbf{C}_{\textbf{z}}^{k+1},\bm{\gamma}^{k+1},\bm{\lambda}^{k+1})$
as a function of the input arguments of $\mathcal{A}(\cdot)$ are given
respectively by (33)-(34), (42)-(43), (46), and (51)-(52). Furthermore,
$(\xi^{k},\beta^{k})_{0\leq k\leq K-1}$ are now learned, instead of being
constant and preset by the user. This leads to the _unfoldedVBA_ architecture
depicted in Fig. 1, which can be summarized into the composition of $K$ layers
$\mathcal{L}_{K-1}\circ\cdots\circ\mathcal{L}_{0}$. Each layer
$\mathcal{L}_{k}$ with $k\in\\{0,\ldots,K-1\\}$ is made of three main blocks,
that are two neural networks, namely $\textrm{NN}_{\sigma}^{k}$ and
$\textrm{NN}_{\xi}^{k}$, and the core VBA block $\mathcal{A}(\cdot)$. There
remains to specify our strategy for building the two inner networks, with the
aim to learn $(\xi^{k},\beta^{k})_{0\leq k\leq K-1}$.

Figure 1: Architecture of _unfoldedVBA_ network.

###  IV-B Learning hyperparameter $\xi$

For every $k\in\\{0,...,K-1\\}$, neural network $\textrm{NN}_{\xi}^{k}$ takes
as input the current kernel estimate
$\textbf{h}^{k}=\textbf{T}\textbf{z}^{k}+\textbf{t}$ and delivers $\xi^{k}$ as
an output. The architecture of the neural network is shown in Fig. 2. Note
that the Softplus function, defined as

$(\forall x\in\mathbb{R})\quad\textrm{Softplus}(x)=\textrm{ln}(1+\exp{(x)}),$
(54)

is used as a last layer, in order to enforce the strict positivity of the
output hyperparameter $\xi^{k}$.

Figure 2: Neural network architecture $\textrm{NN}_{\xi}^{k}$ for estimating
$\xi^{k}$, for $k\in\\{0,...,K-1\\}$.

###  IV-C Learning noise parameter $\beta$

When the noise parameter $\beta$ is unknown, it might be useful to include a
procedure to learn it automatically, again in a supervised fashion. In this
case, we propose to introduce simple nonlinear mappings such that, for every
$k\in\\{0,\dots,K-1\\}$,

$\displaystyle{{\sigma}^{k}}$
$\displaystyle=\textrm{NN}_{\sigma}^{k}(\textbf{y}),$
$\displaystyle=\textrm{Softplus}(\rho_{k})\widehat{\sigma}(\textbf{y})+\textrm{Softplus}(\tau_{k}),$
(55)

and $\beta^{k}=({{\sigma}^{k}})^{-2}$. Hereabove, y is the observed degraded
image, from which we deduce the wavelet-based variance estimator (also used in
[39]),

$\displaystyle\widehat{\sigma}(\textbf{y})=\frac{\textrm{median}(|\mathbf{W}_{H}\textbf{y}|)}{0.6745},$
(56)

where $|\mathbf{W}_{H}\textbf{y}|$ gathers the absolute value of the diagonal
coefficients of the first level Haar wavelet decomposition of the degraded
image y. Moreover, $(\rho_{k},\tau_{k})_{0\leq k\leq K-1}$ are two scalar
parameters to be learned during the training phase.

###  IV-D Complete architecture

We now present our complete blind deconvolution architecture for grayscale
images and color images in Fig. 3. First, let us notice that VBA and its
unrolled variant is designed for grayscale images. We thus generalized the
architecture from Fig. 3(top), to process color images. To this end, we first
transform the input RGB image to its YUV representation, which takes human
perception into consideration. The network $\textrm{NN}_{\sigma}^{k}$ is first
applied to the luminance part $\textbf{y}_{Y}$ of the image. After applying
the _unfoldedVBA_ network (see Fig. 1), we obtain $\textbf{z}_{\text{VBA}}$
and $\textbf{x}_{\text{VBA}}$ as outputs. The latter is a restored version of
the luminance channel. The remaining (U,V) color channels are simply obtained
by median filtering of $(\textbf{y}_{U},\textbf{y}_{V})$. Both architectures
in Fig. 3 additionally involve post-processing layers. More precisely, we
first include a linear layer so as to encode the linear transformation (3),
and then deduce the estimated blur kernel $\widehat{\textbf{h}}$. Second, we
also allow a post-processing layer $\mathcal{L}_{\text{pp}}$ acting on the
image, so as to reduce possible residual artifacts, finally yielding
$\widehat{\textbf{x}}$. In the case of color images, the post-processing is
applied on the RGB representation to avoid chromatic artifacts.

###  IV-E Training procedure

The training of both proposed architectures from Fig. 3 requires to define a
loss function, measuring the discrepancy between the output
$(\widehat{\textbf{x}},\widehat{\textbf{h}})$ and the ground truth
$(\overline{\textbf{x}},\overline{\textbf{h}})$, that we denote hereafter by
$\ell(\widehat{\textbf{x}},\widehat{\textbf{h}},\overline{\textbf{x}},\overline{\textbf{h}})$.
In the blind deconvolution application, one can for instance consider a loss
function related to the error reconstruction on the kernel, or to the image
quality, or a combination of both. Two training procedures will be
distinguished and discussed in our experimental section, namely:
Greedy training The parameters of the unfolded VBA are learned in a greedy
fashion so as to minimize the kernel reconstruction error at each layer. Then,
the post-processing network is learned in a second step, so as to maximize an
image quality metric such as the SSIM [51].
End-to-end training The parameters of the complete architecture are learned
end-to-end so as to maximize the image quality metric.

Whatever the chosen training procedure, it is necessary to make use of a back-
propagation step, that is to differentiate the loss function with respect to
all the parameters of the network. Most operations involved in Fig. 3 can be
differentiated efficiently using standard auto-differentiation tools. However,
we observed in our experiments that it is beneficial (and sometimes even
necessary) for a stable training phase to avoid using such tools for
differentiating the VBA layer $\mathcal{A}(\cdot)$ involved in Fig. 1. In
practice, we used the explicit expressions for the partial derivatives of it.
Note that we followed the approach from [52] to obtain the expression of the
derivatives for the CG solver.

---
Figure 3: Proposed blind image restoration pipeline for grayscale (top) and
color (bottom) images.

##  V Numerical experiments

###  V-A Problem formulation and settings

####  V-A1 Problem overview

We focus on the resolution of the blind image deconvolution problem, where
$\widetilde{\textbf{x}}\in\mathbb{R}^{N}$ is an original image, either
grayscale or colored one. We come back to the model presented in Section II-A,
where the linear operator $\widetilde{\textbf{H}}\in\mathbb{R}^{N\times N}$
models the application of a blur kernel
$\widetilde{\textbf{h}}\in\mathbb{R}^{M}$ to the image. The noise n is assumed
to be an additive white Gaussian noise with zero mean and standard deviation
$\sigma$. In the case of color images, we assume that the same kernel, and the
same noise level, is applied to each of the three RGB channels.

####  V-A2 VBA settings

Let us now specify our practical choices, for the implementation of the VBA
step. In all our experiments, we seek for kernels whose entries satisfy two
equality constraints, namely a sum to one constraint, and an axial symmetry
along the main diagonal axis. This can be easily translated into the affine
constraint (3). In such case, the degree of freedom of the kernel model is
equal to $P=\frac{(\sqrt{M}+1)\sqrt{M}}{2}-1$. Regarding the choice for the
prior, we set $\textbf{A}\in\mathbb{R}^{(2M+1)\times M}$ as the matrix that
computes the horizontal and vertical differences between pixels, augmented
with an additional first row corresponding to an averaging operation, which
takes the form $[1,\ldots,1]/M$. This choice allows to promote smooth
variations in the kernels, while satisfying the required full column rank
assumption on A. A constant vector with entries equal to $\frac{1}{M}$ is set
for the prior mean m. Matrices $(\textbf{D}_{j})_{1\leq j\leq J}$ and
parameter $\kappa$, involved in (8), are set in such a way that the chosen
prior on the image yields an isotropic total-variation regularization (see our
comment in Sec. II-B2). We must now specify the initialization for VBA
iterates/layers. Our initial guess $\textbf{x}^{0}$ for the image is the
degraded image. The associated covariance matrix $\textbf{C}_{\textbf{x}}^{0}$
is initialized using the identity matrix. The blur is initialized with a
uniform kernel with size $5\times 5$, from which we deduce the corresponding
$\textbf{z}^{0}$, and the covariance matrix $\textbf{C}_{\textbf{z}}^{0}$ is
set to a multiple of identity matrix. The hyperparameters $(\alpha,\eta)$
involved in the prior law on parameter $\gamma$are set to zero in practice
which is equivalent to impose a non-informative Jeffrey improper prior.
Finally, the conjugate gradient solver used for the update of the image is run
over 10 iterations which appears sufficient to reach practical stability. The
solver is initialized with the degraded image.

####  V-A3 Datasets

Let us now introduce the two datasets we employ to train and test our network,
and compare it to state-of-the-art techniques. In both cases, the training set
is made of 100 images from the COCO training set. The validation set contains
40 images from the BSD500 validation set. The test set consists of 30 images
from the Flickr30 test set. Each image is center-cropped using a window of
size $N=256\times 256$. Each original image $\widetilde{\textbf{x}}$ is
associated to a degraded version of it, y, built from Model (1). Various blur
kernels and noise levels are used, as detailed hereafter.
Dataset 1: All involved images are converted in grayscale. Each image of the
database is blurred with 10 randomly generated Gaussian blurs, and then
corrupted by additive noise. Thus in total, we have 1000 ($=100\times 10$)
training images, 400 ($=40\times 10$) validation images, and 300 ($=30\times
10$) test images for Dataset 1. The Gaussian blurs are of size $9\times 9$.
Two of them are isotropic with standard deviation randomly generated following
a uniform law within $[0.2,0.4]$. Eight of them are anisotropic with
orientation either $\pi/4$ or $3\pi/4$ (with equal probability) and
vertical/horizontal widths (i.e., standard deviations of the 2D Gaussian
shape) uniformly drawn within $[0.15,0.4]$. On this dataset, the noise
standard deviation is set to $\sigma=0.01$, and assumed to be known (so that
blocks $(\text{NN}^{k}_{\sigma})_{1\leq k\leq K}$ of our architecture are
overlooked).
Dataset 2: All the images are then colored ones. We degraded each of them with
15 different blurs, namely 10 Gaussian blurs (simulated using the same
procedure as above), two uniform blurs with width $5\times 5$ and $7\times 7$
pixels, and 3 out-of-focus blurs. For the latter, the vertical and horizontal
widths are set randomly within $[0.2,0.5]$, and the orientation is either
$\pi/4$ or $3\pi/4$ (with equal probability). Furthermore, for each blurred
image, zero-mean Gaussian noise is added, with standard deviation $\sigma$
randomly chosen, with uniform distribution over $[0.005,0.05]$. The noise
level is not assumed to be known and is estimated using the proposed
$\text{NN}_{\sigma}$ architecture. In total, we have 1500 ($=100\times 15$)
training images, 600 ($=40\times 15$) validation images, and 450 ($=30\times
15$) test images, on this dataset.

Examples of blurs involved in _Dataset 2_ are depicted in Fig. 4.

---
Figure 4: Examples of blur kernels used to construct _Dataset 2_.

####  V-A4 Training specifications

We present results obtained by adopting the two training strategies described
in Section IV-E. For the _greedy training_ , we make use of the mean squared
error on the estimated kernel, as a loss function for the _unfoldedVBA_
layers, defined as
$\ell(\widehat{\textbf{x}},\widehat{\textbf{h}},\widetilde{\textbf{x}},\widetilde{\textbf{h}})=\|\widetilde{\textbf{h}}-\widehat{\textbf{h}}\|^{2}$.
The SSIM loss ([51]), between the output image $\widehat{\textbf{x}}$ and the
ground truth $\widetilde{\textbf{x}}$ is used to train the post-processing
layer $\mathcal{L}_{\rm pp}$. For the _end-to-end training_ , we use again
SSIM between $\widehat{\textbf{x}}$ and $\widetilde{\textbf{x}}$. We use warm
initialization for _end-to-end training_ , that is we initialize with the
weights learnt during the _greedy training_ phase, associated with a weight
decay procedure. ADAM optimizer, with mini-batch size equal to 10, is used for
all the training procedures. Its parameters such as learning rate (lr), weight
decay (wd) and epochs number are finetuned, so as to obtain stable performance
on each validation set. The number of layers $K$ (i.e., number of VBA
iterations) is set during the _greedy training_ , and kept the same for the
_end-to-end training_. In practice, we increase $K$ as long as a significant
decrease in the averaged MSE over the training set was observed. We specify in
Table I all the retained settings. The train/validation/test phase are
conducted with a code implemented in Pytorch (version 1.7.0) under Python
(version 3.6.10) environment, and run on an Nvidia DGX wokstation using one
Tesla V100 SXM2 GPU (1290 MHz frequency, 32GB of RAM). Our code is made
available at https://github.com/yunshihuang/unfoldedVBA.

|  Dataset 1 |  Dataset 2
---|---|---
Greedy training | UnfoldedVBA
$K=6$, epoch $=10$ |  $K=21$, epoch $=10$
lr $=5\times 10^{-3}$ |  lr $=5\times 10^{-3}$ (for $\mathcal{L}_{0}$), lr $=10^{-3}$ (for other layers)
Post-processing $\mathcal{L}_{pp}$
U-net [53] |  Residual network [39, Fig.4]
epoch $=30$, lr $=10^{-3}$ |  epoch $=200$, lr $=10^{-3}$
End-to-end training |  $K=6$ |  $K=21$
epoch $=6$ |  epoch $=6$
|  lr $=5\times 10^{-5}$ |  lr $=5\times 10^{-5}$
|  wd $=10^{-4}$ |  wd $=10^{-4}$
TABLE I: Settings for the training phases in our experiments

####  V-A5 Comparison to other methods

The proposed method is compared to several blind deconvolution approaches
available in the literature:
Optimization-based methods: We first evaluate the VBA described in Section
III, in the favorable situation where the noise level $\sigma$ is assumed to
be known, and parameter $\xi$ is finetuned empirically (see more details
hereafter). VBA is run until reaching practical convergence, i.e. when the
relative squared distance between two consecutive image iterates gets lower
than $10^{-5}$. We also test two optimization-based approaches for blind
deconvolution. The first one is called _deconv2D_. It makes use of the
proximal alternating algorithm from [17], to minimize a least-squares data
fidelity term combined with various priors, namely total variation and
positivity constraint on the image, sum-to-one and quadratic constraint on the
kernel. This method is implemented in Matlab, and inherits some of the
software accelerations discussed in [16] for blind video deconvolution. The
second competitor in this category is the _blinddeconv_ approach 111Matlab
code: https://dilipkay.wordpress.com/blind-deconvolution/ from [11]. For the
sake of fair comparisons, for both datasets, we finetune the hyperparameters
of these three methods on 40% of the training set and apply an average of the
found values on the test set. Moreover, following the use of these three
methods, we perform a non-blind deconvolution step BM3D-DEB 222Matlab code:
https://webpages.tuni.fi/foi/GCF-BM3D/index.html#ref_software [54], which uses
their respective estimated blur kernel to restore the image.
Deep learning methods: We perform comparisons with three recent deep learning
architectures for blind deconvolution. SelfDeblur 333Python/Pytorch code:
https://github.com/csdwren/SelfDeblur [35] is an unsupervised approach able to
jointly perform the image restoration and kernel estimation tasks. DBSRCNN
444Python/Pytorch code: https://github.com/Fatma-ALbluwi/DBSRCNN [34] and
DeblurGAN 555Python/Pytorch (training) and Matlab C-mex (testing) codes:
https://github.com/KupynOrest/DeblurGAN [32] are two supervised deep learning
techniques. In contrast with SelfDeblur, they both only provide the estimated
image, but do not estimate the kernel. Both these methods have been retrained
on each of our datasets, using the same settings as in their initial
implementation. Moreover, we adapted DBSRCNN to color images using the same
pipeline as for our method, that is applying DBSRCNN on the luminance channel
while simply nonlinearly filtering the chrominance ones.

####  V-A6 Evaluation metrics

All the methods are evaluated in terms of their performance on the blur kernel
estimation (when available) and on the image restoration. Different metrics
are used. For the blur kernels, we evaluate (i) the MSE, (ii) the so-called
$\mathcal{H}_{\infty}$ error defined as the $\ell_{\infty}$ norm of the
difference between the 2D discrete Fourier coefficients (with suitable
padding) of the estimated and of the true kernel, and (iii) the mean absolute
error (MAE) defined as the $\ell_{1}$ norm of the difference between
$\widetilde{\textbf{h}}$ and $\widehat{\textbf{h}}$. For evaluating the image
quality, we compute (i) the SSIM, (ii) the PSNR (Peak-Signal-to-Noise Ratio),
and (iii) the PieAPP value [55], between the estimated image
$\widehat{\textbf{x}}$ and the ground truth $\widetilde{\textbf{x}}$.

###  V-B Experimental results

####  V-B1 Dataset 1

In Tables II and III, we report the results of kernel estimation and image
restoration, computed on the test set, using the different methods. As could
be expected, the greedy approach tends to give more weight to the kernel
quality than the end-to-end training. Our two training approaches yield great
performance, when compared to all the other tested approaches. One can notice
that the VBA with finetuned value for $\xi$ performs quite well, showing the
validity of our Bayesian formulation. The proposed unrolled VBA technique
allows us to avoid a manual tuning of this parameter, and further increases
the resulting quality. This is a direct outcome of the supervised training
procedure aiming at maximizing quality scores, and also to the introduction of
a post-processing step on the images. DBSRCNN has a good performance in terms
of image restoration in this dataset. However, it is not capable of estimating
the blur kernel, which might be useful for various applications. We display
two examples of results in Fig 6, extracted from our test set. One can notice,
by visual inspection of these results, the high quality of the restored
images. No artifacts can be observed, which is confirmed by a low average
value of the PieAPP index on the test set. Moreover, the kernels are generally
estimated quite accurately, as shown by the low MSE score and the good
retrieval of their general structure. In the few cases when the unfolded VBA
algorithm fails to give a perfect recovery of the blur kernel as in Fig.
6(bottom), the estimation is still accurate enough to yield a good recovery of
the image whatever _greedy training_ or _end-to-end training_ is used. One can
also notice that our method tends to provide better contrasted images,
compared to its closest competitor in the image restoration task that is
DBSRCNN. We display in Fig. 5(left) the evolution of the SSIM loss during the
end-to-end training of the proposed architecture, showing the increase of the
loss, then its stabilization, for both training and validation set, thus
confirming an appropriate setting of ADAM optimizer parameters. Finally, Table
IV(left) displays the average test time for each methods, that is the
computational time required to restore one example of the dataset, once the
method is finetuned/trained. We displayed CPU time for a fair comparison
between methods, for codes ran on a Dell workstation equipped with an Xeon(R)
W-2135 processor (3.7 GHz clock frequency and 12 GB of RAM). GPU time is also
indicated when available. The fastest method is DBSRCNN, though we must
emphasize that, in contrast with all the other methods based on Matlab/Python
softwares, DBSRCNN makes use of an optimized C implementation, for its test
phase on CPU. DeblurGAN is also very fast, but the resulting quality was quite
poor in our experiments. Our method reaches a reasonable computational time on
CPU. It becomes quite competitive when making use of GPU implementation, as
the unrolled architecture is well suited for that purpose. This allows to drop
the test time per image to few seconds, making it advantageous, with the
addition benefit of better quality results in average, and of an available
kernel estimate.

---
Figure 5: Evolution of SSIM loss along epochs of _end-to-end training_ phase, averaged either on training or on validation sets of _Dataset 1_ (left) and _Dataset 2_ (right).  Method | MSE |  $\mathcal{H}_{\infty}$ error | MAE
---|---|---|---
VBA | 0.0017 (0.0022) | 0.1674 (0.1003) | 0.0472 (0.0317)
deconv2D | 0.0025 (0.0035) | 0.1483 (0.1037) | 0.0489 (0.0395)
blinddeconv | 0.0013 (0.0011) | 0.1553 (0.0660) | 0.0417 (0.0203)
SelfDeblur | 8.9165 (4.0668) | 15.0213 (6.8490) | 3.5314 (0.8998)
Proposed (greedy) |  0.0008 (0.0012) |  0.1165 (0.0677) |  0.0281 (0.0168)
Proposed (end-to-end) | 0.0009 (0.0013) | 0.1188 (0.0672) | 0.0289 (0.0170)
TABLE II: Quantitative assessment of the restored kernels. Mean (standard deviation) values computed over the test sets of _Dataset 1_. Method | SSIM | PSNR | PieAPP
---|---|---|---
Blurred | 0.6542 (0.1072) | 22.2254 (2.3779) | 4.1794 (0.9005)
VBA | 0.7603 (0.0814) | 23.7332 (2.5672) | 1.5109 (0.6184)
deconv2D | 0.7668 (0.0912) | 24.5970 (2.8656) | 1.9289 (0.4959)
blinddeconv | 0.7528 (0.0963) | 23.9347 (2.4299) | 1.9170 (0.6630)
SelfDeblur | 0.6948 (0.1006) | 22.2704 (2.1255) | 3.3178 (0.7291)
DBSRCNN | 0.7780 (0.0895) |  24.9561 (2.9800) | 1.5959 (0.6463)
DeblurGAN | 0.6613 (0.0731) | 22.4388 (2.4074) | 1.8937 (0.7630)
Proposed (greedy) | 0.7945 (0.0890) |  24.7093 (2.9351) | 1.4047 (0.6437)
Proposed (end-to-end) |  0.7989 (0.0886) |  24.6638 (3.0711) |  1.1976 (0.5433)
TABLE III: Quantitative assessment of the restored images. Mean (standard deviation) values computed over the test sets of _Dataset 1_. |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/1_blur9_x0.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/1_blur9_h0.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/1_blur9_VBA_golden_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/1_blur9_unfoldedVBA_golden.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/1_blur9_deconv2D_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/1_blur9_deconv2D.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/1_blur9_norm_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/1_blur9_norm_sparsity_average.eps}} \end{overpic}
---|---|---|---|---
Degraded | Original | VBA | deconv2D | blinddeconv
|  | MSE = 0.0010 | MSE = 0.0016 |  MSE = $5.3788\times 10^{-4}$
| PieAPP = 4.7468 | PieAPP = 1.6914 | PieAPP = 2.7365 | PieAPP = 1.7290
\begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/1_blur9_selfdeblur_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/1_blur9_selfdeblur.eps}} \end{overpic} |  |  |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/1_blur9_unfoldedVBA_lpp_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/1_blur9_unfoldedVBA_block7.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/1_blur9_unfoldedVBA_lpp_ssim_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/1_blur9_unfoldedVBA_lpp_ssim.eps}} \end{overpic}
SelfDeblur | DBSRCNN | DeblurGAN | proposed (greedy) | proposed (end-to-end)
MSE = 6.2415 |  |  |  MSE = 1.2555$\mathbf{\times 10^{-4}}$ |  MSE =1.1817$\mathbf{\times 10^{-4}}$
PieAPP = 4.1130 | PieAPP = 1.8130 | PieAPP = 1.9762 | PieAPP = 1.2950 | PieAPP = 1.2088
|  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/29_blur2_x0.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/29_blur2_h0.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/29_blur2_VBA_golden_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/29_blur2_unfoldedVBA_golden.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/29_blur2_deconv2D_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/29_blur2_deconv2D.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/29_blur2_norm_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/29_blur2_norm_sparsity_average.eps}} \end{overpic}
Degraded | Original | VBA | deconv2D | blinddeconv
|  | MSE = 0.0054 | MSE = 0.0058 | MSE = 0.0040
PieAPP = 3.4024 |  | PieAPP = 1.6356 | PieAPP = 1.9079 | PieAPP = 1.8397
\begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/29_blur2_selfdeblur_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/29_blur2_selfdeblur.eps}} \end{overpic} |  |  |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/29_blur2_unfoldedVBA_lpp_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/29_blur2_unfoldedVBA_block7.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/29_blur2_unfoldedVBA_lpp_ssim_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/29_blur2_unfoldedVBA_lpp_ssim.eps}} \end{overpic}
SelfDeblur | DBSRCNN | DeblurGAN | proposed (greedy) | proposed (end-to-end)
MSE = 7.0703 |  |  | MSE = 0.0035 |  MSE = 0.0034
PieAPP = 2.5108 | PieAPP = 1.4002 | PieAPP = 1.5206 | PieAPP = 1.3922 |  PieAPP =1.2468
Figure 6: Ground-truth image/blur, degraded image, restored images (with PieAPP index) and estimated blurs (with MSE score) when available, for various methods, on two examples in the test set of _Dataset 1_. Method | Dataset 1 | Dataset 2
---|---|---
VBA | 153s (15s) | 156s (18s)
deconv2D | 16s | 19s
blinddeconv | 19s | 22s
SelfDeblur | 452s (51s) | 455s (55s)
DBSRCNN | 1s | 2s
DeblurGAN | 2s (1s) | 3s (2s)
Proposed | 36s (4s) | 113s (12s)
TABLE IV: Average test time per image, using CPU (resp. GPU).

####  V-B2 Dataset 2

The results of kernel estimation and image restoration on Dataset 2 using the
various methods are shown in Tables V and VI, respectively. This dataset is
more challenging, as it includes color images, various blur shapes, and
various noise levels. The latter are not assumed to be known anymore.
Hereagain, we can observe that the _greedy training_ yields the best
performance in terms of kernel estimation for the three considered metrics. In
contrast, _end-to-end training_ tends to favor the restored image quality
while still providing a good kernel quality compared to other methods. In this
more complicated context, standard VBA does not perform very well, as setting
$\xi$ becomes tedious for such an heterogeneous dataset. Let us note that the
noise level is assumed to be known for this particular method, putting it in a
quite favorable situation, compared to the other competitors, including our
proposed approach. DBSRCNN provides again a good image recovery, but our
proposed approach still outperforms it for both SSIM and PieAPP metrics. We
display two examples of restoration in Fig 7, when the sought blur is uniform,
and out-of-focus, respectively. Such blur shapes are challenging and the MSE
on the estimated blur might appear not excellent. Nevertheless, our method
remains the best among the compared ones. The visual quality of the image
generated by the proposed method is also very satisfying. We display in Fig.
5(right) the evolution of the SSIM loss during the _end-to-end training_ ,
witnessing the absence of any overfitting issue. Moreover, we present in Fig 8
the evolution of the MSE loss on the kernel estimate, along the $K=21$ layers
of the architecture trained in an _end-to-end_ manner. The MSE was averaged on
test set examples associated to either Gaussian or out-of-focus blurs,
respectively. These plots show that, for our choice of $K$ (finetuned on the
validation set), the MSE values are close to minimal. Larger $K$ implied an
increase of memory and training time, while not necessarily improving the
results quality. One can also notice more fluctuations in the case of out-of-
focus blur, which turns out to be more challenging to restore. A similar curve
was obtained for uniform blurs, not shown by lack of space. Finally, Table
IV(right) presents the average test time of the different methods. Again, our
method appears competitive in terms of running time.

Method | MSE |  $\mathcal{H}_{\infty}$ error | MAE
---|---|---|---
VBA | 0.0148 (0.0139) | 0.4492 (0.1638) | 0.1339 (0.0627)
deconv2D | 0.0099 (0.0160) | 0.2796 (0.1692) | 0.0869 (0.0576)
blinddeconv | 0.0245 (0.0264) | 0.3113 (0.1409) | 0.1596 (0.1106)
SelfDeblur | 1.7533 (1.4455) | 2.5647 (2.7609) | 1.3752 (0.5132)
Proposed (greedy) |  0.0037 (0.0079) |  0.1888 (0.1061) |  0.0570 (0.0414)
Proposed (end-to-end) | 0.0039 (0.0079) | 0.1960 (0.1056) | 0.0588 (0.0411)
TABLE V: Quantitative assessment of the restored kernels. Mean (standard deviation) values computed over the test sets of _Dataset 2_. Method | SSIM | PSNR | PieAPP
---|---|---|---
Blurred | 0.5427 (0.1150) | 21.7994 (2.1679) | 4.2378 (0.8539)
VBA | 0.4024 (0.1571) | 16.0371 (4.1798) | 2.4218 (0.5545)
deconv2D | 0.6880 (0.1065) | 23.1940 (2.8986) | 2.2245 (0.6721)
blinddeconv | 0.6961 (0.1034) | 23.2663 (2.7229) | 2.3259 (0.8080)
SelfDeblur | 0.5107 (0.1305) | 19.9943 (2.1467) | 5.9269 (1.4066)
DBSRCNN | 0.6948 (0.1688) |  23.6041 (4.2073) | 1.9474 (0.7171)
DeblurGAN | 0.3370 (0.0740) | 17.2781 (1.2909) | 3.6581 (1.0040)
Proposed (greedy) | 0.7454 (0.1015) | 23.2169 (2.4442) |  1.7250 (0.5324)
Proposed (end-to-end) |  0.7518 (0.1025) |  23.5631 (2.5959) |  1.7681 (0.5502)
TABLE VI: Quantitative assessment of the restored images. Mean (standard deviation) values computed over the test sets of _Dataset 2_. |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/7_blur10_x0.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/7_blur10_h0.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/7_blur10_VBA_golden_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/7_blur10_VBA_golden.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/7_blur10_deconv2D_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/7_blur10_deconv2D.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/7_blur10_norm_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/7_blur10_norm_sparsity_average.eps}} \end{overpic}
---|---|---|---|---
Degraded | Original | VBA | deconv2D | blinddeconv
|  | MSE = 0.0115 | MSE = 0.0190 | MSE = 0.0075
PieAPP = 3.8507 |  |  PieAPP = 0.8278 | PieAPP = 1.6273 | PieAPP = 1.0216
\begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/7_blur10_selfdeblur_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/7_blur10_selfdeblur.eps}} \end{overpic} |  |  |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/7_blur10_unfoldedVBA_lpp_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/7_blur10_unfoldedVBA_lpp.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/7_blur10_unfoldedVBA_lpp_ssim_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/7_blur10_unfoldedVBA_lpp_ssim.eps}} \end{overpic}
SelfDeblur | DBSRCNN | DeblurGAN | proposed (greedy) | proposed (end-to-end)
MSE = 3.6152 |  |  |  MSE = 0.0015 | MSE = 0.0020
PieAPP = 8.1479 | PieAPP = 1.2058 | PieAPP = 2.6971 | PieAPP = 1.2793 | PieAPP = 1.1663
|  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/30_blur12_x0.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/30_blur12_h0.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/30_blur12_VBA_golden_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/30_blur12_VBA_golden.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/30_blur12_deconv2D_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/30_blur12_deconv2D.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/30_blur12_norm_BM3D_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/30_blur12_norm_sparsity_average.eps}} \end{overpic}
Degraded | Original | VBA | deconv2D | blinddeconv
|  | MSE= 0.0293 | MSE = 0.0184 | MSE = 0.0388
PieAPP = 2.8997 |  | PieAPP = 1.3625 | PieAPP = 1.0608 | PieAPP = 1.0318
\begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/30_blur12_selfdeblur_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/30_blur12_selfdeblur.eps}} \end{overpic} |  |  |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/30_blur12_unfoldedVBA_lpp_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/30_blur12_unfoldedVBA_lpp.eps}} \end{overpic} |  \begin{overpic}[height=85.35826pt,trim=28.45274pt 28.45274pt 28.45274pt 28.45274pt,clip]{figures/30_blur12_unfoldedVBA_lpp_ssim_mk.eps} \put(22.0,6.0){\includegraphics[height=28.45274pt,trim=113.81102pt 42.67912pt 96.73918pt 28.45274pt,clip]{figures/30_blur12_unfoldedVBA_lpp_ssim.eps}} \end{overpic}
SelfDeblur | DBSRCNN | DeblurGAN | proposed (greedy) | proposed (end-to-end)
MSE = 0.9536 |  |  | MSE = 0.0049 |  MSE = 0.0046
PieAPP = 2.0619 | PieAPP = 0.8841 | PieAPP = 2.7761 | PieAPP = 0.7790 |  PieAPP = 0.6987
Figure 7: Ground-truth image/blur, degraded image, restored images (with
PieAPP index) and estimated blurs (with MSE score) when available, for various
methods, on two examples in the test set of _Dataset 2_.
---
Figure 8: MSE loss along the layers for proposed method using _end-to-end
training_. Average over _Dataset 2_ test examples involving either Gaussian
(left), or out-of-focus (right) blur shapes.

##  VI Conclusion

This paper proposes a novel method for blind image deconvolution that combines
a variational Bayesian algorithm with a neural network architecture. Our
experiments illustrate the excellent performance of this new method on two
datasets, comprising grayscale and color images, and degraded with various
kernel types. Compared to state-of-the-art variational and deep learning
approaches, our method delivers a more accurate estimation of both the image
and the blur kernels. It also includes an automatic noise estimation step, so
that it requires little hyperparameter tuning. The proposed method is very
competitive in terms of computational time during the test phase, while
showing similar train time to its deep learning competitors. The main core of
the proposed architecture is highly interpretable, as it implements unrolled
iterates of a well sounded Bayesian-based blind deconvolution method. As a
byproduct, it also outputs estimates for the covariance matrices of both
sought quantities (image/kernel). This information could be of interest for
uncertainty quantification and model selection tasks (see for instance [7,
56]). More generally, our work demonstrates that unrolling VBA algorithms
constitutes a promising research direction for solving challenging problems
arising in Data Science.

## References

  * [1] F. Murtagh, E. Pantin, and J.-L. Starck, Deconvolution and blind deconvolution in astronomy, pp. 277–316.  Taylor and Francis, 2007.
  * [2] T. J. Holmes, D. Biggs, and A. Abu-Tarif, Blind Deconvolution, pp. 468–487.  Boston, MA: Springer US, 2006.
  * [3] O. Michailovich and A. Tannenbaum, “Blind deconvolution of medical ultrasound images: A parametric inverse filtering approach,” IEEE Transactions on Image Processing, vol. 16, no. 12, pp. 3005–3019, 2007.
  * [4] A. S. Carasso, “Direct blind deconvolution,” SIAM Journal on Applied Mathematics, vol. 61, no. 6, pp. 1980–2007, 2001.
  * [5] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman, “Understanding and evaluating blind deconvolution algorithms,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2009), (Miami, FL, USA), pp. 1964–1971, 20-25 Jun. 2009.
  * [6] T. Bell, J. Xu, and S. Zhang, “Method for out-of-focus camera calibration,” Applied Optics, vol. 55, no. 9, pp. 2346–2352, 2016.
  * [7] Y. Huang, E. Chouzenoux, and V. Elvira, “Probabilistic modeling and inference for sequential space-varying blur identification,” IEEE Transactions on Computational Imaging, vol. 7, pp. 531–546, 2021.
  * [8] E. Chouzenoux, T.-K. Lau, C. Lefort, and J.-C. Pesquet, “Optimal multivariate Gaussian fitting with applications to psf modeling in two-photon microscopy imaging,” Journal of Mathematical Imaging and Vision, vol. 61, no. 7, pp. 1037–1050, 2019.
  * [9] N. Komodakis and N. Paragios, “MRF-based blind image deconvolution,” in Proceedings of the Asian Conference on Computer Vision (ACCV 2012), (Daejeon, Korea), pp. 361–374, 5-9 Nov. 2012.
  * [10] J. Kotera, F. Šroubek, and P. Milanfar, “Blind deconvolution using alternating maximum a posteriori estimation with heavy-tailed priors,” Computer Analysis of Images and Patterns, vol. 8048, pp. 59–66, 2013.
  * [11] D. Krishnan, T. Tay, and R. Fergus, “Blind deconvolution using a normalized sparsity measure,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2011), (Colorado Springs, CO, USA), pp. 233 – 240, 20-25 Jun. 2011.
  * [12] T. Chan and C.-K. Wong, “Total variation blind deconvolution,” IEEE transactions on image processing, vol. 7, no. 3, pp. 370–375, 1998.
  * [13] A. Levin, “Blind motion deblurring using image statistics,” in Proceedings of the Annual Conference on Neural Information Processing Systems (NIPS 2006), vol. 19, (Vancouver, British Columbia, Canada), pp. 841–848, 4-7 Dec. 2006.
  * [14] Q. Shan, J. Jia, and A. Agarwala, “High-quality motion deblurring from a single image,” ACM Transactions on Graphics, vol. 27, pp. 1–10, Aug. 2008\.
  * [15] N. Joshi, R. Szeliski, and D. J. Kriegman, “PSF estimation using sharp edge prediction,” in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR 2008), (Anchorage, AK, USA), pp. 1–8, 23-28 Jun. 2008.
  * [16] F. Abboud, Émilie Chouzenoux, J.-C. Pesquet, J.-H. Chenot, and L. Laborelli, “An alternating proximal approach for blind video deconvolution,” Signal Processing: Image Communication, vol. 70, pp. 21–36, 2019.
  * [17] J. Bolte, P. L. Combettes, and J.-C. Pesquet, “Alternating proximal algorithm for blind image recovery,” in Proceedings of the IEEE International Conference on Image Processing (ICIP 2010), (Honk Kong, China), pp. 1673–1676, 26-29 Sep. 2010.
  * [18] C. P. Robert and G. Casella, Monte Carlo Statistical Methods.  Springer, 2004.
  * [19] D. Luengo, L. Martino, M. Bugallo, V. Elvira, and S. Särkkä, “A survey of monte carlo methods for parameter estimation,” EURASIP Journal on Advances in Signal Processing, vol. 2020, pp. 1–62, 2020.
  * [20] M. Evans and T. Swartz, “Methods for approximating integrals in statistics with special emphasis on bayesian integration problems,” Statistical Science, vol. 10, no. 3, pp. 254–272, 1995.
  * [21] O. Rosec, J.-M. Boucher, B. Nsiri, and T. Chonavel, “Blind marine seismic deconvolution using statistical MCMC methods,” IEEE Journal of Oceanic Engineering, vol. 28, no. 3, pp. 502–512, 2003.
  * [22] D. Ge, J. Idier, and É. Le Carpentier, “Enhanced sampling schemes for MCMC based blind Bernoulli-Gaussian deconvolution,” Signal Processing, vol. 91, pp. 759–772, Apr. 2011.
  * [23] G. Kail, J.-Y. Tourneret, F. Hlawatsch, and N. Dobigeon, “Blind deconvolution of sparse pulse sequences under a minimum distance constraint: A partially collapsed gibbs sampler method,” IEEE Transactions on Signal Processing, vol. 60, no. 6, pp. 2727–2743, 2012.
  * [24] T. E. Bishop, R. Molina, and J. R. Hopgood, “Blind restoration of blurred photographs via ar modelling and mcmc,” in Proceedings of the 15th IEEE International Conference on Image Processing (ICIP 2008), pp. 669–672, 2008\.
  * [25] C. Fox and S. Roberts, “A tutorial on variational bayesian inference,” Artificial Intelligence Review, vol. 38, p. 85–95, 2012.
  * [26] D. Blei, A. Kucukelbir, and J. McAuliffe, “Variational inference: A review for statisticians,” Journal of the American Statistical Association, vol. 112, Jan. 2016.
  * [27] Y. Marnissi, Y. Zheng, E. Chouzenoux, and J.-C. Pesquet, “A variational Bayesian approach for image restoration—application to image deblurring with Poisson–Gaussian noise,” IEEE Transactions on Computational Imaging, vol. 3, no. 4, pp. 722–737, 2017.
  * [28] Y. Zheng, A. Fraysse, and T. Rodet, “Efficient variational Bayesian approximation method based on subspace optimization,” IEEE Transactions on Image Processing, vol. 24, no. 2, pp. 681–693, 2015.
  * [29] S. Babacan, R. Molina, and A. Katsaggelos, “Variational Bayesian blind deconvolution using a total variation prior,” IEEE Transactions On Image Processing, vol. 18, pp. 12–26, Feb. 2009.
  * [30] M. Pereyra, P. Schniter, E. Chouzenoux, J.-C. Pesquet, J.-Y. Tourneret, A. O. Hero, and S. McLaughlin, “A survey of stochastic simulation and optimization methods in signal processing,” IEEE Journal of Selected Topics in Signal Processing, vol. 10, p. 224–241, Mar. 2016.
  * [31] Y. Marnissi, E. Chouzenoux, A. Benazza-Benyahia, and J.-C. Pesquet, “Majorize–minimize adapted Metropolis–Hastings algorithm,” IEEE Transactions on Signal Processing, vol. 68, pp. 2356–2369, 2020.
  * [32] O. Kupyn, V. Budzan, M. Mykhailych, D. Mishkin, and J. Matas, “DeblurGAN: Blind motion deblurring using conditional adversarial networks,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2018), (Salt Lake City, UT, USA), pp. 8183–8192, 18-22 Jun. 2018.
  * [33] C. Dong, C. C. Loy, K. He, and X. Tang, “Image super-resolution using deep convolutional networks,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 38, no. 2, pp. 295–307, 2016.
  * [34] F. Albluwi, V. A. Krylov, and R. Dahyot, “Image deblurring and super-resolution using deep convolutional neural networks,” in Proceedings of the IEEE 28th International Workshop on Machine Learning for Signal Processing (MLSP 2018), (Aalborg, Denmark), pp. 1–6, 17-20 Sep. 2018\.
  * [35] D. Ren, K. Zhang, Q. Wang, Q. Hu, and W. Zuo, “Neural blind deconvolution using deep priors,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2020), (Seattle, WA, USA), pp. 3338–3347, 16-18 Jun. 2020.
  * [36] A. Bietti, G. Mialon, and J. Mairal, “On regularization and robustness of deep neural networks.” arXiv:1810.00363, 2018.
  * [37] V. Monga, Y. Li, and Y. C. Eldar, “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,” IEEE Signal Processing Magazine, vol. 38, no. 2, pp. 18–44, 2021.
  * [38] Y. Li, M. Tofighi, J. Geng, V. Monga, and Y. C. Eldar, “Efficient and interpretable deep blind image deblurring via algorithm unrolling,” IEEE Transactions on Computational Imaging, vol. 6, pp. 666–681, 2020.
  * [39] C. Bertocchi, E. Chouzenoux, M.-C. Corbineau, J.-C. Pesquet, and M. Prato, “Deep unfolding of a proximal interior point method for image restoration,” Inverse Problems, vol. 36, Sep. 2019.
  * [40] J.-C. Pesquet, A. Repetti, M. Terris, and Y. Wiaux, “Learning maximally monotone operators for image recovery,” SIAM Journal on Imaging Sciences, vol. 14, no. 3, pp. 1206–1237, 2021.
  * [41] P. L. Combettes and J.-C. Pesquet, “Deep neural network structures solving variational inequalities,” Set-Valued and Variational Analysis, vol. 28, pp. 491–518, Sep. 2020.
  * [42] E. Chouzenoux, C. D. Valle, and J.-C. Pesquet, “Inversion of integral models: a neural network approach.” arXiv:2105.15044, 2021.
  * [43] Y. Sun, Z. Wu, X. Xu, B. Wohlberg, and U. S. Kamilov, “Scalable plug-and-play ADMM with convergence guarantees,” IEEE Transactions on Computational Imaging, vol. 7, pp. 849–863, 2021.
  * [44] D. Kingma and M. Welling, “Auto-encoding variational Bayes,” in Proceedings of the International Conference on Learning Representations (ICLR 2014), (Banff, AB, Canada), 14-16 Apr. 2014.
  * [45] L. V. Jospin, W. L. Buntine, F. Boussaid, H. Laga, and M. Bennamoun, “Hands-on Bayesian neural networks - a tutorial for deep learning users,” ArXiv, vol. abs/2007.06823, 2020.
  * [46] D. Gilton, G. Ongie, and R. Willett, “Model adaptation for inverse problems in imaging,” IEEE Transactions on Computational Imaging, vol. 7, pp. 661–674, 2021.
  * [47] Y. Li, O. Bar-Shira, V. Monga, and Y. C. Eldar, “Deep algorithm unrolling for biomedical imaging,” Computing Research Repository, vol. abs/2108.06637, 2021.
  * [48] B. Tolooshams, S. Mulleti, D. Ba, and Y. C. Eldar, “Unfolding neural networks for compressive multichannel blind deconvolution.” arXiv:2010.1139 [Online], 2021\.
  * [49] Y. Nan, Y. Quan, and H. Ji, “Variational-EM-based deep learning for noise-blind image deblurring,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2020), (Seattle, WA, USA), 16-18 Jun. 2020.
  * [50] Y. Li, M. Tofighi, V. Monga, and Y. C. Eldar, “An algorithm unrolling approach to deep image deblurring,” in Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2019), (Brighton, UK), pp. 7675–7679, 12-17 May 2019.
  * [51] Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Transactions on Image Processing, vol. 13, no. 4, pp. 600–612, 2004.
  * [52] B. Charlier, J. Feydy, J. Glaunès, F.-D. Collin, and G. Durif, “Kernel operations on the GPU, with autodiff, without memory overflows,” Journal of Machine Learning Research, vol. 22, no. 74, pp. 1–6, 2021.
  * [53] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI 2015), vol. 9351, (Munich, Germany), pp. 234–241, 5-9 Oct. 2015.
  * [54] M. Lebrun, “An analysis and implementation of the BM3D image denoising method,” Image Processing On Line, vol. 2, pp. 175–213, Aug. 2012.
  * [55] E. Prashnani, H. Cai, Y. Mostofi, and P. Sen, “PieAPP: Perceptual image-error assessment through pairwise preference,” in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR 2018), (Salt Lake City, UT, USA), pp. 1808–1817, 18-22 Jun. 2018.
  * [56] A. Repetti, M. Pereyra, and Y. Wiaux, “Scalable Bayesian uncertainty quantification in imaging inverse problems via convex optimization,” SIAM Journal on Imaging Sciences, vol. 12, no. 1, pp. 87–118, 2019.