Plug-and-Play Quantum Adaptive Denoiser for Deconvolving Poisson Noisy Images

A new Plug-and-Play (PnP) alternating direction of multipliers (ADMM) scheme is proposed in this paper, by embedding a recently introduced adaptive denoiser using the Schroedinger equation's solutions of quantum physics. The potential of the proposed model is studied for Poisson image deconvolution, which is a common problem occurring in number of imaging applications, such as limited photon acquisition or X-ray computed tomography. Numerical results show the efficiency and good adaptability of the proposed scheme compared to recent state-of-the-art techniques, for both high and low signal-to-noise ratio scenarios. This performance gain regardless of the amount of noise affecting the observations is explained by the flexibility of the embedded quantum denoiser constructed without anticipating any prior statistics about the noise, which is one of the main advantages of this method. The main novelty of this work resided in the integration of a modified quantum denoiser into the PnP-ADMM framework and the numerical proof of convergence of the resulting algorithm.


I. INTRODUCTION
R ESTORATION of a distorted image is one of the most fundamental tasks in inverse problems related to imaging applications such as denoising, deblurring, superresolution, compression or compressed sensing. In number of applications such as limited photon acquisition, X-ray computed tomography, positron emission tomography, etc., the noise degrading the acquired data follows a Poisson distribution. These Poissonian models have been extensively studied in the fields of astronomical [1]- [3], photographic [4], [5] or biomedical [6]- [11] imaging. The inversion process is expressed as the estimation of a clean image x ∈ R n from observed degraded image y ∈ R m . The estimation of the underlying hidden image from this distorted observation is often formulated as the optimization of a cost function implementing the idea of the maximum a posteriori (MAP) estimator [12], i.e., the maximization of the posterior proba-bility, defined asx where P (x|y) is the posterior probability density function that defines x for a given measurement y andx represents the estimation of the unobserved image x. Taking −log(·) element wise and applying the Bayes' theorem, the maximization problem above becomeŝ x = arg min x − log (P (y|x)) − log (P (x)) + log (P (y)) . (2) f (x) = −log (P (y|x)) is the negative log-likelihood function whose expression depends on the observation (degradation) model, and g(x) = −log (P (x)) is the a priori logdistribution of x, that only depends on some prior knowledge on the image to estimate and is also called regularization VOLUME 00, 0000 1 arXiv:2107.00407v3 [eess.IV] 20 Oct 2021 function. Note that P (y) does not depend on x and is usually ignored in the estimation ofx. With these notations, the optimization problem to solve can be expressed aŝ Using a suitable choice of the regularization function, based for example on the a priori statistics of the image to estimate, proximal operator- [13] based iterative schemes have been extensively studied to solve (3) [14]- [24]. In particular, the alternating direction method of multipliers (ADMM) [18]- [24] has been largely used, by redefining the optimization problem (3) into a constrained optimization framework. During the last decade, a new approach was proposed in the literature, enabling the use of state-of-theart denoisers instead of the proximal operator, known as the plug-and-play (PnP) scheme [25]. PnP paves the way of using a wide range of state-of-the-art denoisers such as patch-based dictionary learning methods [26], block-matching 3D filtering (BM3D) [27], non-local means (NLM) [28], high-order variational models [29], etc. The interest of PnP schemes in image restoration have been shown by number of studies, e.g., [30]- [45]. Interestingly, these PnP-ADMM methods do not require any prior information about the hidden image x, as a consequence of the intrinsic association between the regularizer and the external denoiser.
More recently, alternative learning-based approaches were developed in the literature using deep learning (convolutional neural network (CNN)) techniques for tackling inverse problems [46]- [50]. During the past few years the implementation of these Deep-CNN networks has been introduced for image denoising [51], [52] and further extended to the PnP schemes [53]. These Deep-CNN networks give several advantages such as reconstruction accuracy and convergence speed [54]. However, more often they suffer from some drawbacks. First, such denoisers should be trained using the noise variance in each iteration. Hence, during the iterative process of the PnP framework, the noise variance is usually unknown since it varies at each iteration, and leads to a divergence of the algorithm for a pre-trained Deep-CNN architecture [55]. Second, the training procedure is very costly since Deep-CNN denoisers require expensive retraining whenever the noise level or noise type change. Also, each iteration involves a Deep-CNN denoising process, so using a large neural network and/or too many iterative operations leads to a time consuming task. Third, the theoretical aspects of Deep-CNN denoiser-based PnP models are still not clear.
This work focuses on PnP-ADMM algorithms applied to Poisson deconvolution problems, i.e., recover an image from a blurred observation contaminated by Poisson noise. Since the state-of-the-art denoisers (e.g., BM3D [27]) used within PnP schemes were primarily designed for additive Gaussian noise, they consequently exhibit inconsistency with a non-Gaussian model. Furthermore, decoupling the restoration and denoising steps within PnP frameworks alternatively converts the noise distribution affecting the observed distorted image into a possibly different noise model, and in particularly into a non-Gaussian noise. To mitigate this limitation, a variance stabilizing transformation (VST) [56]- [59], known as the Anscombe transformation, was embedded in several PnP-ADMM algorithms to adapt them to a data-dependent model. Indeed, VST was designed to remodel approximately a random data-dependent noise into an additive Gaussian noise, before processing through a Gaussian denoiser. Although these refined VST-based PnP schemes exhibit very good performance for low-intensity noise [30]- [32] and outperform existing state-of-the-art prior based models, they are less accurate while dealing with high-intensity noise (i.e., low SNR) [60]. Furthermore, the nonuniform nature of the convolution operator under a VST leads to fundamental flaws in the deconvolution algorithms [31], [32], [61].
In this paper, we address these shortcomings by embedding into a PnP-ADMM scheme a new adaptive denoiser [62], [63] designed by borrowing tools from quantum mechanics. The adaptive nature of this denoiser makes it highly efficient at selectively eliminating noise from higher intensity pixels, without relying on any statistical assumption about the noise [64]. Its efficiency regardless of the assumption of Gaussian noise represents the main motivation of its interest in Poisson deconvolution PnP-ADMM algorithms, discarding the necessity of a VST. To summarize, the main novelty of the paper is the use of quantum mechanical concepts in the field of image restoration. The primary contributions are the quantum denoiser, its integration into a PnP-ADMM scheme, and the experimental proof of convergence of the final algorithm.
The remainder of the paper is organized as follows. After a brief discussion on PnP-ADMM algorithms in Section II, the construction of the proposed method referred to as QAB-PnP is illustrated for Poisson inverse problems in Section III. Section IV regroups the numerical experiments and Section V draws the conclusions and the perspectives.

A. ALTERNATING DIRECTION METHOD OF MULTIPLIERS
ADMM is an iterative convex optimization algorithm, resulting from the fusion of the dual decomposition method with the method of multipliers [65]- [70]. Several developments have been proposed during the last few decades, resulting into a rapidly growing literature [16]- [20]. ADMM algorithm is able to solve constrained optimization problems of the form minimize x,z f (x) + g(z) where f and g are assumed to be closed convex functions of variables x ∈ R n and z ∈ R m , with A ∈ R p×n , B ∈ R p×m and c ∈ R p . The associated augmented Lagrangian function is defined as where v ∈ R p is the Lagrangian multiplier, and λ ∈ R + is the penalty parameter of the augmented Lagrangian. An equivalent expression of the augmented Lagrangian L λ (x, z, v) can be obtained by scaling the Lagrangian multiplier u = (1/λ)v, as follows: ADMM algorithm decouples the augmented Lagrangian into three iterative steps as follows: The convergence of this iterative scheme has been widely discussed in the literature of convex programming and within various statistical problems [71]- [73]. ADMM technique has a broad spectrum of applications in the context of signal and image restoration applications [74]- [78].

B. ADMM APPLICATION TO IMAGE RESTORATION
Let us consider the following general image restoration problem, characterized by the forward model where y is the observed image related to the underlying image x through the degradation operator O. ADMM can be used to estimate the MAP solution of such an image restoration task by reformulating it as (4) using the following parameterization: z = x, thus A = −B = I n×n , c = 0 n , where I n×n is the identity matrix of size n × n and 0 n is a zero vector of size n. The associated augmented Lagrangian is given by where f (x) = −log (P (y|x)) is the data fidelity term depending on O and g(z) the regularization function. To accelerate the convergence, the penalty parameter λ is usually increased at each iteration, by multiplication by a factor of γ > 1 [39], instead of using a fixed value. At each iteration, ADMM performs the following steps:

C. PLUG-AND-PLAY (PNP) FRAMEWORK
Since its initial development, the PnP scheme [25] is largely accepted for signal and image restoration problems due to its extremely promising performance [30]- [45]. The primary goal of PnP is to consider a state-of-the-art denoiser as the prior of a constrained optimization process. Interestingly, no prior knowledge is required about the image to estimate to derive the regularization function g, since g is intrinsically defined through the external denoiser used. The efficiency of ADMM algorithm mainly reposes on its ability of decoupling the optimization processes over each variable, as shown in the previous section. ADMM steps performed at each iteration, (12), (13) and (14), can be interpreted as follows. (12) is originally an inversion step to get the best possible primary image satisfying the data through the data fidelity function f (x), while the third step (14) updates the Lagrangian multiplier. The second step (13) can be rewritten as The expression on the right hand side of (16) fundamentally intends to find the solution that optimizes the compromise between the difference between z and (x k+1 + u k ) and the regularization function g(z). Thus, it can be associated to a denoising problem designed to denoise (x k+1 + u k ). Therefore it is possible to rewrite this step as where D(·) is a denoising operator. Hence it is feasible to implement a state-of-the-art denoiser to handle the denoising operation as proposed in [25]. The most interesting feature representing the key benefit of this approach is that this PnP model does not require the prior term g(z) explicitly, rather it is indirectly related to the choice of the denoiser D(·) (see, e.g., [27], [28], [79], [80]). Despite its interest shown in number of imaging applications, PnP-ADMM still presents important theoretical challenges while dealing with Poisson deconvolution. Indeed, most advanced denoisers available in the literature generally consider additive Gaussian models and cannot be implemented directly for other noise removal processes which do not follow Gaussian statistics. Furthermore, despite observing an image degraded by a specific noise model (e.g., Poisson in our case), the image (x k+1 + u k ) to be denoised VOLUME 00, 0000 at each iteration in (17) does not necessarily follow the same noise distribution. Therefore, handling an inverse problem using the PnP-ADMM algorithm requires to transform the unknown noise distribution of the noisy image (x k+1 + u k ) into an additive Gaussian distribution before implementing a Gaussian denoiser. In this context VST-like [56] transformations propose an efficient way of estimating approximately a Gaussian distribution from other types of data-dependent models. The convolution product is however not invariant under this VST and consequently leads to theoretical flaws. Therefore, a versatile denoiser adapted to different noise models, without a priori hypothesis about the noise statistics, is desirable to be efficient regardless of the prior noise distribution in this PnP framework.
In this work, our primary focus will be on the formulation of a PnP-ADMM model using an adaptive denoiser, constructed from the principles of quantum mechanics [62], [63], and its implementation into Poisson deconvolution processes. This quantum adaptive basis (QAB)-based denoiser does not require any explicit noise model. Therefore, while included in an PnP-ADMM scheme, it does not need the use of a VST before denoising and mitigates this theoretical limitation.

D. CONVERGENCE OF PNP-ADMM ALGORITHMS
One major challenge of PnP-ADMM algorithms is to prove their convergence, due to the implicit relation between the regularization function g(z) and the denoising operator D(·). Note that the convergence of conventional ADMM has been largely discussed in the literature, primarily in [71] and [16] and more recently in [18] based on the proximal operator [81] or in [82]. The proof of global convergence of PnP-ADMM algorithm [38] has been shown in the case of non-expansive denoisers belonging to the family of symmetric smoothing filters [83]- [86]. Yet these conditions are too restrictive for generalisation to all the denoisers. To overcome this issue, a series of works has been published during the last few years showing the fixed point convergence of PnP-ADMM algorithms for bounded denoisers not necessarily symmetric and non-expansive [39]- [45], but we stress that all these algorithms were constructed for Gaussian noise.

A. POISSONIAN DECONVOLUTION MODEL
Let us denote by x ∈ R n 2 the image to be recovered from the observation y ∈ R n 2 , a degraded version by a point spread function (PSF) and Poisson process denoted by P(·). Without loss of generality, we consider herein square images of size n × n, written as vectors in lexicographical order. The resulting image formation model is where H ∈ R n 2 ×n 2 is a block circulant with circulant blocks (BCCB) matrix acounting for 2D circulant convolution with the PSF. The pixels of the observed blurry and noisy image y are denoted by y[i], i = 1, 2, · · · , n 2 , and are contemplated as the independent realizations of a Poisson process with parameter (Hx)[i] ≥ 0 given by where (·)[i] represents the i-th component of a vectorized image. The restoration of x from the noisy-blurred observation y is the primary objective of Poisson deconvolution methods. One standard way to estimate x from the observation model (18) is to use the MAP estimator in (1). The Poisson noise probability density function is defined as Thus, the log-likelihood term, i.e., the data fidelity term f (x) used within the MAP estimator, is given by where 1 is a vector of length n 2 with all elements equal to 1. As explained previously, the function g(x) in (3), a prior of x, depends on some prior knowledge on the image to estimate. In a PnP framework, this prior is intrinsically defined through the external denoiser, removing the fact of defining the prior term g(x) explicitly. Hence, using the data fidelity term f (x) in (21), the PnP-ADMM steps depicted in (12), (14), (15) and (17) become: where D(·) is the denoising operator considered within the PnP-ADMM algorithm. In this work, following [87], a gradient descent algorithm is used to solve the minimization problem (22), that requires the use of the gradient of the augmented Lagrangian L λ given by where ∇ x represents the derivative with respect to x and y/(Hx) stands for element-wise division.
The following subsection describes the Poisson denoiser inspired from quantum mechanics used within the proposed PnP-ADMM algorithm for Poisson image deconvolution, to solve the step in (23).
The denoiser embedded in the proposed method is based on the construction of an adaptive basis inspired from quantum mechanics, as originally proposed in [62], [63]. An illustration of the adaptive basis construction is given in Fig. 2. It displays the relationship between a clean and a noisy image in the quantum mechanical framework. The basic idea is to use the image as a potential of a quantum system, where the height of the potential is determined by the pixel intensity. For illustration purpose, we considered the Boat image with half of it contaminated by Gaussian noise. Two patches, one clean and one noisy, are extracted from the image and plotted as 3D surfaces, which will ultimately act as the potentials of the system. In this system, the wave function governs the probability of presence of a quantum particle with energy E at some position on the surface. For a clean image, the wave function uses a broad range of frequencies to probe the surface. In presence of random noise, the wave function collapses and becomes localized at some particular position on the surface, as highlighted in Fig. 2  the adaptive basis is the fact that the pixel intensity is directly linked to the local frequency of the wave. The localization property in the presence of noise is actually a hindrance, cured by performing a pre-smoothing of the noisy potential in order to create an adaptive basis extended over the whole image. For more details on the construction of the basis, we refer the reader to [62]. For self-consistency, we recall hereafter the main steps of the QAB (quantum adaptive basis) technique.

1) Background on the adaptive QAB transform
In the non-relativistic quantum mechanics, the timeindependent Schroedinger equation yields an equation for the stationary wave solution ψ(a), given by where is the Planck constant and ψ(a) characterizes the energy state E of the particle with mass m in a potential V . The probability amplitude of the particle is given by |ψ(a)| 2 , normalized under |ψ(a)| 2 da = 1. The wave function ψ(a) is an element of the Hilbert space of square-integrable functions. It is possible to rewrite the equation (27) as where H QAB = − 2 2m ∇ 2 + V is the Hamiltonian operator. One can conclude from (28) that the solution ψ(a) of the equation (27) represents an eigen-state of the system described by the Hamiltonian operator. These eigen-states of (28) are oscillatory functions and primarily have two properties: i) the oscillation frequency increases with energy and ii) for the same eigen-function, the local frequency depends on the local value of the potential, and this dependence is regulated by the value of 2 /2m which acts as a hyperparameter herein.
In the perspective of designing an adaptive transformation for image processing, one may consider the image pixels' values as the potential V in the Schroedinger equation (27) for a discretized space. The stationary solutions of (27) can be obtained by computing the eigen-pairs of the discretized Hamiltonian operator defined as: where x ∈ R n 2 is an image (i.e., V = x), vectorized in lexicographical order and H QAB [i, j] represents the (i, j)th element of the operator H QAB ∈ R n 2 ×n 2 . Note that zero padding is used to handle the boundary conditions. As a consequence some violations of the rule (29) can be observed. More precisely, .., n, n 2 − n + 1, n 2 − n + 2, ..., n 2 } in order to respect the boundary conditions, and H QAB [i, i + 1] = H QAB [i + 1, i] = 0 for any i multiple of n apart from n 2 . More details about the construction of the Hamiltonian operator associated to an image can be found in [62].
The corresponding eigenbasis of the Hamiltonian operator (29) represents the adaptive transform. In the seminal works [62], [63], it was shown that this adaptive basis gives an efficient way of image denoising, especially in the presence of Gaussian, Poisson or speckle noise. In this work, this adaptive basis, referred to as quantum adaptive basis (QAB), is used to construct the denoiser D QAB (·) embedded in the proposed PnP-ADMM scheme.
These basis vectors belong to the family of oscillating functions along with the Fourier and wavelet bases, but with a local frequency depending on the local value of 2m(E − V )/ . Due to its dependence on the difference between the energy E and potential V , in the same basis vector the lower values of the potential are associated with oscillations of higher frequency. Thus, the property of these adaptive basis vectors able to describe different image pixels' values using different frequency levels, makes it fundamentally distinct from the Fourier and wavelet bases. From the above discussion it is understandable that the local frequency depends on the value of 2 /2m, which is a hyperparameter. Apart from that, the level of noise also has an impact on the basis vectors. Indeed, the presence of random noise in the system leads to a subtle quantum phenomenon [96] which makes these vectors localize exponentially at different positions of the potential in the system. To mitigate this phenomenon which degrades the denoising, it is important to low-pass the corrupted image using, for example, a Gaussian filter with suitable standard deviation σ QAB , before the computation of the QAB from the Hamiltonian operator (29). The reader may refer to [62] for an in-depth discussion about the QAB vector localization in the presence of noise.
The QAB explained above is used to denoise an image, as suggested in [62], as follows: project the noisy image onto the QAB to identify the valuable information and the noise, followed by a soft-thresholding of the projection coefficients, before taking the inverse projection of the modified coefficients to recover the noise-free image. The denoised image is retrieved as following: with where α i are the coefficients representing the image x in QAB, whose basis vectors are ψ i . s and ρ are two thresholding hyperparameters. The denoising process thus corresponds to expanding the signal in the adaptative basis and thresholding the coefficients according to an energy criterion (see [62] for a detailed discussion of this procedure).

C. QAB-PNP ALGORITHM
This section illustrates, in the context of Poisson image deconvolution, the proposed PnP-ADMM algorithm, denoted as QAB-PnP, incorporating the QAB denoiser introduced in the previous section. In this particular context, various stateof-the-art denoisers have been introduced in the literature, such as Gaussian denoisers (e.g., BM3D [27], etc) fused with VST-like transforms or not. Using QAB D QAB instead of a classical denoiser is the main contribution of this work. It consists in including a modified version of the QAB denoiser into the deconvolution PnP-ADMM method from Section III-A, more precisely to solve (23). The denoising process integrated in the proposed QAB-PnP algorithm requires the computation of the coefficients α i , obtained by projecting the noisy image onto the QAB. This is a time consuming task for a large image and affects the computational load of the deconvolution algorithm given that the denoising process is performed at each iteration. However, one may note that most of the α i are not used for reconstructing the denoised image given that they are discarded by the threshohlding operation. To increase the computational efficiency of the proposed algorithm, only the coefficients which contribute the most in the restoration process are computed. To this end, let us focus on T basis vectors α i from D QAB , corresponding to an energy level below E, assuming that higher energy levels naturally correspond to higher frequencies, where E is considered as a free hyperparameter. The corresponding T coefficients will be the most significant for the reconstruction of the clean image, and can be computed using the orthogonal matching pursuit (OMP) algorithm [97]- [100].
The OMP algorithm was fundamentally designed to obtain a sparse approximationα i with sparsity T of the corresponding coefficients α i while projecting the noisy image, say v ∈ R n 2 onto the denoising basis D QAB . Therefore VOLUME 00, 0000 Algorithm 1: Modified Orthogonal Matching Pursuit algorithm.
Compute the sparse coefficientsα i with sparsity T by using the measurement data z and the operator D QAB following the modified orthogonal matching pursuit method as illustrated in Algorithm 1.
Output:ẑ the primary goal of OMP is to recover coefficientsα i with T non-zero elements, such that v D QABαi . To get the best possible approximation, it is important to identify the columns ψ i ∈ D QAB which contribute in the reconstruction of v. The basic idea is to choose the column of D QAB which is mostly correlated with v, followed by subtracting its contribution and repeat the step on the residual. After T iterations one can have the desired set of basis vectors and projection coefficients. Within the adaptive basis D QAB , the basis eigenvectors are organized in ascending order, the first T basis vectors with energy less than E being the most correlated with v. Therefore, the OMP algorithm is modified herein so that it estimates only the projection coefficients onto the subspace formed by these T basis vectors. This modified OMP algorithm is detailed in Algorithm 1.
The sparse coefficientsα i estimated by Algorithm 1 are further used by the denoising method detailed in Algorithm 2, integrated in the proposed QAB-PnP deconvolution method in Fig. 1 and Algorithm 3 1 . 1 The Matlab code of the proposed Plug-and-Play-ADMM algorithm using the quantum-adaptive-basis denoiser is [Online]. Available: https://github.com/SayantanDutta95/QAB-PnP-ADMM-Deconvolution.git Algorithm 3: Poisson deconvolution using QAB-PnP algorithm.
, σ QAB , N Initialization: x 0 , z 0 , u 0 Compute a smooth version of y by low-pass Gaussian filter with standard deviation σ QAB Form the Hamiltonian matrix H QAB based on the smoothed version of y using (29) Calculate the eigen-pairs of H QAB Construct D QAB using the eigenvectors ψ i of H QAB Find the total number of eigenvalue T , less than the energy level E begin ADMM process: for k from 0 to N − 1 do Step 1: Step 2: z k+1 = D QAB (x k+1 + u k ), following QAB denoising Algorithm 2 Step 3: The computational complexity of the algorithm is dominated by the eigendecomposition of the high dimensional Hamiltonian matrix and the QAB image projection. For a n×n image, the Hamiltonian matrix is of size n 2 ×n 2 . Usual textbook diagonalization methods would require O(n 6 ) operations (time complexity) and O(n 4 ) storage space. However, the Hamiltonian matrix is extremely sparse, and is more efficiently diagonalized by iterative methods such as the Lanczos method (as we actually did). In this case the computational complexity would be O(n 4 ) if we compute all eigenvalues and eigenvectors (and still O(n 4 ) in storage space). If we compute only T of these eigenvalues and eigenvectors (with T ≤ n 2 ), the time complexity becomes O(T n 2 ) and the storage space (space complexity) also O(T n 2 ). The QAB (a) Performed on the image in Fig. 3(a).
(b) Performed on the image in Fig. 3(b).
(c) Performed on the image in Fig. 3(c). ≤ σ k M for any x k ∈ R n 2 , performed on the sample images in Fig. 3(a-c).
image projection is O(n 4 ) with the simplest algorithm, and becomes O(T n 2 ) in time and space with the OMP algorithm. We thus conclude that our algorithm requires O(T n 2 ) time and space resources, with T ≤ n 2 , for a n × n image.
To further decrease the complexity, a block-wise approach could be used as proposed in [62], where a large image is divided into smaller patches denoised independently by the QAB denoiser. In this the complexity is O(T P m 2 ) for P patches of size m (<< n). Moreover, such a patch-based architecture can be improved by considering the dependence between neighboring patches by borrowing tools from the quantum interaction theory as suggested in [101].

D. CONVERGENCE ANALYSIS OF QAB-PNP ALGORITHM
Despite their popularity during the last decade, the proof of convergence of PnP-ADMM algorithms may still be an issue. Some interesting developments have been proposed during the last few years on global [38] and fixed point [39]- [45] convergence of these algorithms, while imposing restrictions on the denoising operator. In this section, our goal is to analyse the fixed point convergence of the proposed QAB-PnP algorithm.
To enable the fixed point convergence and in particular to avoid the issue of unbounded gradient in (26) for pixel values equal to 0, i.e., to overcome the singularity problem at x = 0, we slightly modify the observation model (18) by introducing a small positive constant 1, as suggested in [102]: Therefore the negative Poisson log-likelihood (21) becomes and the corresponding gradient ∇f (x) = −H T (y/(Hx + 1)) + H T 1.
One should note that within practical experiments, is much smaller than any background value, so that its influence on the final output is negligible [102].
Proof: Since is the lower bound of (Hx + 1), therefore 1/ is the upper bound of 1/(Hx + 1). Since y and H are constants, they are bounded. Hence one can write: where δ 1 , δ 2 , L ∈ R + .

Remark 2.
Denoiser D QAB is a bounded denoising operator with a parameter σ k .
We cannot offer a general proof of this statement, also it intuitively appears highly likely. The denoising process denoted by D QAB certainly reduces the level of noise at each iteration and gets D QAB (x k ) closer and closer to x k . It is therefore fair to consider that D QAB (x k ) − x k 2 decreases with k. It is also bounded by x k 2 since D QAB is a projection operator.
The rate of decrease is not a priori easy to bound, but we offer numerical evidence that the decrease is fast. Indeed, in all three examples shown in Fig. 4 the decrease is very fast. In particular, it is much faster that the rate of decrease of σ k def = 1/λ k . We thus generalize this result and take as generic
* Second condition: The second condition should hold generically as discussed in Remark 2.
Given that the two conditions are satisfied within the proposed framework, let us move to the proof of the fixed point convergence in Remark 3. We start by proving the following statements: where C 1 , C 2 and C 3 are constants and λ k is the penalty parameter with λ k+1 = γλ k , where γ > 1. * First step: Proof of condition (36). From (12), we have The first order optimality implies Since the minimizer is obtained in x = x k+1 , replacing x by x k+1 and using the boundedness property of ∇f (x), we have Furthermore, since the denoiser D QAB is bounded and z k+1 = D QAB (x k+1 + u k ), one can write One also has Finally, using (41) and (42), we obtain ) * Second step: Proof of condition (38).
From (14), we get Using (45), we have (46) * Third step: Proof of condition (37). (14) can be written as Using (47), we have Hence all three conditions (36), (37) and (38) are true. Next, we aim at proving that {x k } ∞ k=1 is a Cauchy sequence. Therefore, one has to show that for all integer n > k, x n − x k 2 → 0 as n → ∞ and k → ∞. For any finite n and k, one can write using the condition (37) Therefore, as n → ∞ and k → ∞, x n − x k 2 → 0, since γ > 1, so {x k } ∞ k=1 is a Cauchy sequence. Hence, the sequence {x k } ∞ k=1 is convergent, thus there exits x * ∈ [0, 1] n 2 such that x k − x * 2 → 0 as k → ∞. Similarly, one can show that the sequence {z k } ∞ k=1 and {u k } ∞ k=1 are convergent, so there exit z * , u * ∈ [0, 1] n 2 such that z k − z * 2 → 0 and u k − u * 2 → 0 as k → ∞. Therefore we can conclude that the proposed QAB-PnP algorithm converges to a fixed point.
The proof we propose is not a convergence proof in the mathematical sense, since it reposes on Remark 2 for which we only have plausibility arguments and numerical evidence. Nevertheless, the discussion above and the numerical results in Fig. 4 for three very different images, indicate that with high confidence the algorithm should converge in practice for any image.

IV. SIMULATION RESULTS
This section illustrates the efficiency of the proposed QAB-PnP algorithm for Poisson image deconvolution. An analysis of the influence of the hyperparameters on the deconvolution accuracy is first provided in Subsection IV-A, before comparing its performance to several state-of-the-art methods in Subsection IV-B. In [62] we already performed a detailed analysis of the hyperparameters σ QAB , s and ρ for the efficiency of the denoiser. We recall that these hyperparameters control respectively the smoothing of the potential to avoid localization effects in the expansion basis, and the cutoff in energy which leads to denoising. We therefore chose these hyperparameters to be optimal according to the study in [62]. However, the computational method used in the present work (OMP algorithm) introduces a new hyperparameter E which controls the accuracy and efficiency of the OMP process. The accuracy of OMP increases for increasing E, but at the cost of higher computational time. A trade-off is thus necessary, and we will show that the optimal value of E is also influenced by the value of the hyperparameter 2 /2m, which fixes how the local frequencies of the basis vectors vary as a function of pixels' amplitudes.
The simulations are conducted on three images, shown in Fig. 3. Two of them represent cropped versions of the standard Lena and fruits images. The third one was synthetically constructed so that it contains high frequencies for low gray levels and, vice versa, low frequencies for high intensity pixels. Its purpose is to illustrate the ability of the proposed deconvolution method, and in particular of the embedded quantum-based denoiser, to handle such images. All the sample images are distorted with two Gaussian blurring kernel h 4×4 σ of size 4 × 4 and standard deviation σ = 3 and σ = 5 respectively. The study was conducted with three different Poisson noise levels corresponding to SNRs of 20, 15 and 10 dB. Note that the noise was image-dependent Poisson distributed and that the SNRs of the observations was computed a posteriori to emphasize the amount of noise. Finally, Subsection IV-C shows the abbility of the the proposed method to enhance experimental fluorescence microscopy images.

A. HYPERPARAMETER ANALYSIS
This subsection presents a detailed analysis on the influence of the hyperparameters on the proposed method. In particular, the role of the hyperparameter E will be evaluated, given its important impact on the compromise between accuracy and computational time, and its relationship with the hyperparameter 2 /2m will be assessed. It is important to mention that in general the hyperparameter 2 /2m and the number of significant wave vectors T vary in an opposite way, one of them increasing when the other one decreases. In addition, there is a linear relation between T and the processing time. Therefore, to achieve an optimal behaviour of the algorithm, a good balance between the hyperparameters 2 /2m and E needs to be achieved. We will also discuss the choice of the hyperparameter λ 0 which controls the iterations of the ADMM algorithm described in Section II.
From this perspective, we first show that considering the wave vectors up to the energy level E and evaluating only the corresponding coefficients α i following the modified OMP algorithm in Algo. 1 helps reducing the computation time with minimal accuracy loss. Quantitative results showing the influence of E on the simulations performed over the three sample images in Fig. 3, distorted by a Gaussian blurring kernel h 4×4 σ of size 4 × 4 and standard deviation σ = 3, and corrupted by Poisson noise corresponding to a SNR of 20 dB, 15 dB, and 10 dB, have been regrouped in Table 1, where the best results have been highlighted in bold. Similarly, the average peak signal to noise ratios (PSNR) values for different SNR, obtained with the proposed deconvolution method with and without the modified OMP algorithm, are shown in Fig. 5. The results in Fig. 5 and Table 1 prove that the accuracy loss, caused by the use of the parameter E within the modified OMP algorithm, is very limited. This accuracy loss is caused by the denoising process that reconstructs the denoised image only from the wave functions associated with an energy level lower than E. Indeed, although wave functions associated with higher energies are dominated by noise, they may still carry information about certain features of the clean image. The average computation time for different images obtained with a Matlab implementation on a desktop computer, with and without E, given in Table 2, confirms the computational efficiency gain enabled by the modified OMP algorithm embedded in QAB-PnP method.
In addition to E, as stated previously, 2 /2m is also an important hyperparameter of the proposed deconvolution  (a) Performed on the image in Fig. 3(a) with E = 3.9, 2 /2m = 4 (b) Performed on the image in Fig. 3(b) with E = 4.1, 2 /2m = 4 (c) Performed on the image in Fig. 3(c) with E = 4.5, 2 /2m = 4.3 technique. The hyperparameter 2 /2m dictates how the local frequencies of the basis vectors vary with the amplitude of the image pixel values. On the other hand, E is associated with the sparsity. Given their mutual dependence, Fig. 6(a) shows the accuracy of QAB-PnP algorithm for different couple values of these two hyperparameters over an acceptable range. This experiment consisted in recovering the image in Fig. 3(a) from a degraded version blurred by a 4×4 Gaussian kernel with standard deviation equal to 3 and Poisson noise corresponding to a SNR of 20 dB.
Similarly, Figs. 6(b) and (c) show the variation of the number of the significant wave vectors T and of the computation time. These results also justify the linear proportionality of T and processing time. Note that as explained previously, the other hyperparameters, σ QAB , s and ρ, were chosen as suggested in [62].
Finally, the choice of the hyperparameter λ 0 used within the iterations of the ADMM algorithm described in Section II is important to accelerate the convergence. The curves in Fig. 7 show, within a logarithmic scale, the evolution of the root mean square error (RMSE) over the iterations of the proposed deconvolution method, for different values of λ 0 . These simulations were performed for the three images in Fig. 3, distorted by a Gaussian blurring kernel h 4×4 σ of size 4×4 and standard deviation σ = 3, and corrupted by Poisson process corresponding to a SNR of 20 dB.
The studies performed in this subsection show that a certain range of optimal choice of the hyperparameters considered is possible. Without a priori knowledge, it should be possible to use values in this range for arbitrary images, taking care to choose E and 2 /2m in a correlated way. As a further note, keeping the hyperparameters constant to the same values for all the images considered hereafter leads to a very low PSNR degradation of about 0.1 dB. From the discussions above, one may note that the hyperparameters 2 /2m and E are primarily associated with the construction of the quantum adaptive basis and the sparsity of the clean image in this basis, both related to the denoising process. In contrast, λ 0 , the penalty parameter, regulates the restoration process by accelerating the convergence. Therefore, the optimal choice of 2 /2m and E discussed above is independent of the value of λ 0 .

B. POISSON DECONVOLUTION RESULTS
Poisson deconvolution is a well discussed domain in the literature where PnP algorithms implanting a Gaussian denoiser with or without a VST transformation have exhibited promising outcomes [31], [32]. The proposed method is intrinsically adaptive, which makes it well-adapted to different noise statistics for the problem addressed and does not require using any additional transformation in the denoising step.
This subsection regroups image deconvolution results obtained with the proposed method and five approaches from the literature. The experiments consisted in recovering the   Fig. 3 from degraded versions by Gaussian blurring kernels with different variances and Poisson noise at different SNRs. The first comparative method is a standard Poisson deconvolution method that consists in estimating the image that minimizes a cost function formed by the data fidelity term in (21) and the classical total variation regularization [11]. This method will be denoted by TV-ADMM hereafter. The second method denoted by ADMM+BM3D is an inte-gration of the BM3D denoiser in the PnP-ADMM algorithm. Similarly, a deep learning denoiser trained on natural images was integrated into the PnP-ADMM scheme and used as comparison method. In particular, the CNN-based flexible learning method, known as the trainable nonlinear reaction diffusion (TNRD) [103], was used given its efficiency within regularization by denoising approaches [104]. Finally, a PnP-ADMM algorithm coupled with an Anscombe transformation (VST) and a BM3D denoiser, denoted by P 4 IP in [32] was used for comparison. Note that TNRD has been also used with and without VST. The resulting algorithms are denoted by ADMM+TNRD and ADMM+VST+TNRD. It is important to mention that the methods used for comparisons such as TV-ADMM, P 4 IP and ADMM+VST+TNRD are particularly designed for handling data degraded by Poisson noise, and are therefore appropriate choices as comparative methods to the proposed Poisson deconvolution algorithm.
As explained previsouly, the proposed method does not require such a VST-like transformation due to the adaptive nature of the embedded denoiser [62]. Therefore, the proposed algorithm is expected to present better generic convergence properties compared to P 4 IP. In the example in Fig. 8, where P 4 IP had fast convergence, the rate of convergence of QAB-PnP is similar to P 4 IP and faster than TV-ADMM, ADMM+BM3D, ADMM+TNRD and ADMM+VST+TNRD. To evaluate the computational complexity of the proposed algorithm in comparison with other standard techniques, the average computational time and required number of iterations before convergence are given in Table 2 with respect to different images. The results confirm the faster convergence of the proposed method, albeit, at the cost of higher computational time per iteration.
The deconvolution results obtained with the six methods can be visually appreciated in Figs. 9, 10 and 11. The  PSNR and the structure similarity (SSIM) [105] were used to evaluate the deconvolution accuracy. The resulting numerical results, for two different blurring kernels and three different SNRs, are regroupped in Table 3. In particular, average and standard deviation values are reported for 200 noise realizations. For further investigation, the quantitative results obtained with the proposed method in presence of very highintensity noise, in particular, with SNRs close to 5 dB and 0 dB, are provided in Table 4. One may observe that the proposed scheme is capable to adapt both to low and high level of noise and outperforms the five other methods in almost all the simulations. It is important to note that QAB-PnP not only provides the best average values, but also the lowest standard deviations, in particular compared to P 4 IP. This observation is confirmed by the results in Fig. 12, that displays, for a given simulation, the best, the worst and an intermediate result over 200 noise realizations. While the difference between these three results is barely observable for the proposed method, this is not the case for P 4 IP. Finally, one may observe the big accuracy difference between the proposed method and the five others for the synthetic image.

C. APPLICATION TO FLUORESCENCE MICROSCOPY IMAGING
This section highlights the applicability of the proposed deconvolution method to real-life imaging applications, in particular to fluorescence microscopy imaging using, e.g., confocal [106] or two-photon [107] microscopes. Fluorescence microscopy images are intrinsically noisy, contaminated by Poisson-Gaussian noise. Poisson noise is the dominating source of noise [11], [108], [109], due to a limited number (∼ 10 2 per pixel) of quantized photons captured by a microscopic detector compared to normal photography (∼ 10 5 per pixel). Therefore, enhancing such contaminated fluorescence images is of interest for many modern biological studies.        Herein, we used three microscopy images from the online available data-set 2 to illustrate the potential of the proposed 2 http://tinyurl.com/y6mwqcjs method. Fig. 13 regroups the observed distorted images, their corresponding ground truth, and the deblurred images estimated by the six methods. PSNR and SSIM values comparing the observed and the deblurred images to the clean ones are given in Table 5. These results clearly show the efficiency of the proposed algorithm in real fluorescence microscopy image enhancement.

V. CONCLUSIONS
This paper proposed a new PnP-ADMM scheme to handle Poisson deconvolution problems. Although Gaussian denoiser-based PnP-ADMM algorithms have achieved enormous success in this domain of image restoration, they are still facing a theoretical limitation related to the Anscombe transformation used to approximately transform the Poisson noise into additive Gaussian noise. Under this transformation, the convolution operation is not invariant. To overcome this drawback, we proposed in this work the QAB denoiser derived from principles of quantum mechanics, whose architecture makes it well adapted to different noise statistics, explaining its good behavious as denoiser embedded in a PnP-ADMM algorithm.
The simulation results allowed to provide an in-depth analysis of the impact of the hyperparameters on the accuracy and computation efficiency of the proposed method. They also allowed to show its interest compared to five existing methods. An issue of our proposal is the computational burden. The use of the OMP algorithm already dramatically decreases this time compared to earlier implementations [62], but other improvements are certainly possible. As shown in our previous work [62] in the proposed quantum adaptive basis is equally efficient for Gaussian, Poisson and speckle noise removal problems without considering any prior information about the noise statistics. Therefore, the proposed deconvolution method could be suitable for other noise degradation than Poisson, and its evaluation in such conditions represents an interesting perspective. As another future perspective of this work one may think of implementing a more advanced inversion algorithm for a Poissonian model (e.g., SPIRAL-TAP [102]) instead of using a gradient descent method. Moreover, blind deconvolution is also an interesting perspective for future study, by coupling the proposed deconvolution algorithm with a PSF estimation method [110], [111]. Finally, such a PnP scheme can be further extended to other reconstruction problems, such as compressed sensing or super-resolution, using more efficient quantum mechanics based algorithms or by absorbing the patch-based procedure to the proposed framework, using for example the many-body quantum theory.

VI. ACKNOWLEDGMENTS
We would like to thank Dr Yoann Altmann, Assistant Professor at Heriot-Watt University, Edinburgh, Scotland, for his valuable inputs in this work. We also thank CNRS for funding through the 80 prime program.