Inverse Problems with Poisson noise: Primal and Primal-Dual Splitting1footnote 11footnote 1Submitted to ICIP 2011 on the 01/21/11.

# Inverse Problems with Poisson noise: Primal and Primal-Dual Splitting111Submitted to ICIP 2011 on the 01/21/11.

###### Abstract

In this paper, we propose two algorithms for solving linear inverse problems when the observations are corrupted by Poisson noise. A proper data fidelity term (log-likelihood) is introduced to reflect the Poisson statistics of the noise. On the other hand, as a prior, the images to restore are assumed to be positive and sparsely represented in a dictionary of waveforms. Piecing together the data fidelity and the prior terms, the solution to the inverse problem is cast as the minimization of a non-smooth convex functional. We establish the well-posedness of the optimization problem, characterize the corresponding minimizers, and solve it by means of primal and primal-dual proximal splitting algorithms originating from the field of non-smooth convex optimization theory. Experimental results on deconvolution and comparison to prior methods are also reported.

Inverse Problems with Poisson noise: Primal and Primal-Dual SplittingSubmitted to ICIP 2011 on the 01/21/11.

F.-X. Dupé, M.J. Fadili and J.-L. Starck
 a AIM UMR CNRS - CEA 91191 Gif-sur-Yvette France b GREYC CNRS-ENSICAEN-Université de Caen 14050 Caen France

Index Terms—  Inverse Problems, Poisson noise, Duality, Proximity operator, Sparsity.

## 1 Introduction

Linear inverse problems in presence of Poisson noise have attracted less interest in the literature than their Gaussian counterpart, presumably because the noise properties are more complicated to handle. Such inverse problems have however important applications in imaging such as restoration (e.g. deconvolution in medical and astronomical imaging), or reconstruction (e.g. computerized tomography). For instance, the well-known Richardson-Lucy has been proposed for deconvolution. The RL algorithm, however, amplifies noise after a few iterations, which can be avoided by introducing regularization. In [1], the authors presented a Total Variation (TV)-regularized RL algorithm, and [2] advocated a wavelet-regularized RL algorithm.

In the context of Poisson linear inverse problems using sparsity-promoting regularization, a few recent algorithms have been proposed. For example, [3] stabilize the noise and proposed a family of nested schemes relying upon proximal splitting algorithms (Forward-Backward and Douglas-Rachford) to solve the corresponding optimization problem. The work of [4] is in the same vein. However, nested algorithms are time-consuming since they necessitate to sub-iterate. Using the augmented Lagrangian method with the alternating method of multipliers algorithm (ADMM), which is nothing but the Douglas-Rachford splitting applied to the Fenchel-Rockafellar dual problem, [5] presented a deconvolution algorithm with TV and sparsity regularization. This scheme however necessitates to solve a least-square problem which can be done explicitly only in some cases.

In this paper, we propose a framework for solving linear inverse problems when the observations are corrupted by Poisson noise. In order to form the data fidelity term, we take the exact Poisson likelihood. As a prior, the images to restore are assumed to be positive and sparsely represented in a dictionary of atoms. The solution to the inverse problem is cast as the minimization of a non-smooth convex functional, for which we prove well-posedness of the optimization problem, characterize the corresponding minimizers, and solve them by means of primal and primal-dual proximal splitting algorithms originating from the realm of non-smooth convex optimization theory. Convergence of the algorithms is also shown. Experimental results and comparison to other algorithms on deconvolution are finally conducted.

### Notation and terminology

Let a real Hilbert space, here a finite dimensional vector subspace of . We denote by the norm associated with the inner product in , and is the identity operator on . is the norm. and are respectively reordered vectors of image samples and transform coefficients. We denote by the relative interior of a convex set . A real-valued function is coercive, if , and is proper if its domain is non-empty . is the class of all proper lower semicontinuous (lsc) convex functions from to . We denote by the spectral norm of the linear operator , and its kernel.

Let be an image. can be written as the superposition of elementary atoms parameterized by such that . We denote by the dictionary (typically a frame of ), whose columns are the atoms all normalized to a unit -norm

## 2 Problem statement

Consider the image formation model where an input image of pixels is indirectly observed through the action of a bounded linear operator , and contaminated by Poisson noise. The observed image is then a discrete collection of counts which are bounded, i.e. . Each count is a realization of an independent Poisson random variable with a mean . Formally, this writes in a vector form as

 y∼P(Hx)\leavevmode\nobreak . (1)

The linear inverse problem at hand is to reconstruct from the observed count image .

A natural way to attack this problem would be to adopt a maximum a posteriori (MAP) bayesian framework with an appropriate likelihood function — the distribution of the observed data given an original — reflecting the Poisson statistics of the noise. As a prior, the image is supposed to be economically (sparsely) represented in a pre-chosen dictionary as measured by a sparsity-promoting penalty supposed throughout to be convex but non-smooth, e.g. the norm.

From the probability density function of a Poisson random variable, the likelihood writes:

 p(y|x)=∏i((Hx)[i])y[i]exp(−(Hx)[i])y[i]!\leavevmode\nobreak . (2)

Taking the negative log-likelihood, we arrive at the following data fidelity term:

 f1 :η∈Rn↦n∑i=1fpoisson(η[i]), (3) if y[i]>0,fpoisson(η[i]) ={−y[i]log(η[i])+η[i]if η[i]>0,+∞otherwise, if y[i]=0,fpoisson(η[i]) ={η[i]if η[i]∈[0,+∞),+∞otherwise.

Our aim is then to solve the following optimization problems, under a synthesis-type sparsity prior111Our framework and algorithms extend to an analysis-type prior just as well, though we omit this for obvious space limitation reasons.,

 argminα∈H′J(α),J : α↦f1∘H∘Φ(α)+γΨ(α)+ıC∘Φ(α)\leavevmode\nobreak . (Pγ,ψ)

The penalty function is positive, additive, and chosen to enforce sparsity, is a regularization parameter and is the indicator function of the convex set . In our case, is the positive orthant since we are fitting Poisson intensities, which are positive by nature.

From the objective in (), we get the following,

###### Proposition 1.

1. is a convex function and so are and .

2. is strictly convex if . remains strictly convex if is an orthobasis and .

3. Suppose that . Then .

### 2.1 Well-posedness of (Pγ,ψ)

Let be the set of minimizers of problem (). Suppose that is coercive. Thus is coercive. Therefore, the following holds:

###### Proposition 2.

1. Existence: () has at least one solution, i.e. .

2. Uniqueness: () has a unique solution if is strictly convex, or under (ii) of Proposition 1.

## 3 Iterative Minimization Algorithms

### 3.1 Proximal calculus

We are now ready to describe the proximal splitting algorithms to solve (). At the heart of the splitting framework is the notion of proximity operator.

###### Definition 3 ([6]).

Let . Then, for every , the function achieves its infimum at a unique point denoted by . The operator thus defined is the proximity operator of .

Then, the proximity operator of the indicator function of a convex set is merely its orthogonal projector. One important property of this operator is the separability property:

###### Lemma 4 ([7]).

Let and let . Then .

The following result can be proved easily by solving the proximal optimization problem in Definition 3 with as defined in (3), see also [8].

###### Lemma 5.

Let be the count map (i.e. the observations), the proximity operator associated to (i.e. the Poisson anti log-likelihood) is,

 proxβf1x=(x[i]−β+√(x[i]−β)2+4βy[i]2)1⩽i⩽n\leavevmode\nobreak . (4)

We now turn to which is given by Lemma 4 and the following result:

###### Theorem 6 ([9]).

Suppose that : (i) is convex even-symmetric, non-negative and non-decreasing on , and ; (ii) is twice differentiable on ; (iii) is continuous on , and admits a positive right derivative at zero . Then, the proximity operator has exactly one continuous solution decoupled in each coordinate :

 ^α[i]={0if |β[i]|⩽δψ′i+(0)βi−δψ′i(^α[i])if |β[i]|>δψ′i+(0) (5)

Among the most popular penalty functions satisfying the above requirements, we have , in which case the associated proximity operator is soft-thresholding, denoted in the sequel.

### 3.2 Splitting on the primal problem

#### 3.2.1 Splitting for sums of convex functions

Suppose that the objective to be minimized can be expressed as the sum of functions in , verifying domain qualification conditions:

 argminx∈H\leavevmode\nobreak (F(x)=K∑k=1Fk(x))\leavevmode\nobreak . (6)

Proximal splitting methods for solving (6) are iterative algorithms which may evaluate the individual proximity operators , supposed to have an explicit convenient structure, but never proximity operators of sums of the .

Splitting algorithms have an extensive literature since the 1970’s, where the case predominates. Usually, splitting algorithms handling have either explicitly or implicitly relied on reduction of (10) to the case in the product space . For instance, applying the Douglas-Rachford splitting to the reduced form produces Spingarn’s method, which performs independent proximal steps on each , and then computes the next iterate by essentially averaging the individual proximity operators. The scheme described in [10] is very similar in spirit to Spingarn’s method, with some refinements.

#### 3.2.2 Application to Poisson noise inverse problems

Problem () is amenable to the form (6), by wisely introducing auxiliary variables. As () involves two linear operators ( and ), we need two of them, that we define as and . The idea is to get rid of the composition of and . Let the two linear operators and . Then, the optimization problem () can be equivalently written:

 argmin(x1,x2,α)∈H×K×H′f1(x2)+ıC(x1)+γΨ(α)G(x1,x2,α)+ (7) ıkerL1(x1,x2,α)+ıkerL2(x1,x2,α)\leavevmode\nobreak . (8)

Notice that in our case by virtue of separability of the proximity operator of in , and ; see Lemma 4.

The proximity operators of and are easily accessible through Lemma 5 and 6. The projector onto the positive orthant is also trivial. It remains now to compute the projector on , , which by well-known linear algebra arguments, is obtained from the projector onto the image of .

###### Lemma 7.

The proximity operator associated to is

 (9)

The inverse in the expression of is can be computed efficiently when is a tight frame. Similarly, for , the inverse writes , and its computation can be done in the domain where is diagonal; e.g. Fourier for convolution.

Finally, the main steps of our primal scheme are summarized in Algorithm 1. Its convergence is a corollary of [10][Theorem 3.4].

###### Proposition 8.

Let be a sequence generated by Algorithm 1. Suppose that Proposition 1-(iii) is verified, and . Then converges to a (non-strict) global minimizer of ().

### 3.3 Splitting on the dual: Primal-dual algorithm

Our problem () can also be rewritten in the form,

 argminα∈H′F∘K(α)+γΨ(α) (10)

where now and . Again, one may notice that the proximity operator of can be directly computed using the separability in and .

Recently, a primal-dual scheme, which turns to be a pre-conditioned version of ADMM, to minimize objectives of the form (10) was proposed in [11]. Transposed to our setting, this scheme gives the steps summarized in Algorithm 2.

Adapting the arguments of [11], convergence of the sequence generated by Algorithm 2 is ensured.

###### Proposition 9.

Suppose that Proposition 1-(iii) holds. Let , choose and such that , and let as defined by Algorithm 2. Then, converges to a (non-strict) global minimizer () at the rate on the restricted duality gap.

### 3.4 Discussion

Algorithm 1 and 2 share some similarities, but exhibit also important differences. For instance, the primal-dual algorithm enjoys a convergence rate that is not known for the primal algorithm. Furthermore, the latter necessitates two operator inversions that can only be done efficiently for some and , while the former involves only application of these linear operators and their adjoints. Consequently, Algorithm 2 can virtually handle any inverse problem with a bounded linear . In case where the inverses can be done efficiently, e.g. deconvolution with a tight frame, both algorithms have comparable computational burden. In general, if other regularizations/constraints are imposed on the solution, in the form of additional proper lsc convex terms that would appear in (), both algorithms still apply by introducing wisely chosen auxiliary variables.

## 4 Experimental results

Our algorithms were applied to deconvolution. In all experiments, was the -norm. Table 1 summarizes the mean absolute error (MAE) and the execution times for an astronomical image, where the dictionary consisted of the wavelet transform and the PSF was that of the Hubble telescope. Our algorithms were compared to state-of-the-art alternatives in the literature. In summary, flexibility of our framework and the fact that Poisson noise was handled properly, demonstrate the capabilities of our approach, and allow our algorithms to compare very favorably with other competitors. The computational burden of our approaches is also among the lowest, typically faster than the PIDAL algorithm. Fig. 1 displays the objective as a function of the iteration number and time (in s). We can clearly see that Algorithm 2 converges faster than Algorithm 1.

## 5 Conclusion

In this paper, we proposed two provably convergent algorithms for solving the Poisson inverse problems with a sparsity prior. The primal-dual proximal splitting algorithm seems to perform better in terms of convergence speed than the primal one. Moreover, its computational burden is lower than most comparable of state-of-art methods.

## References

• [1] N. Dey et al., “A deconvolution method for confocal microscopy with total variation regularization,” in ISBI 2004. 2004, pp. 1223–1226, IEEE.
• [2] J.-L. Starck and F. Murtagh, Astronomical Image and Data Analysis, Springer, 2006.
• [3] F.-X. Dupé, M.J. Fadili, and J.-L.Starck, “A proximal iteration for deconvolving poisson noisy images using sparse representations,” IEEE Trans. on Im. Pro., vol. 18, no. 2, pp. 310–321, 2009.
• [4] C. Chaux, J.-C. Pesquet, and N. Pustelnik, “Nested iterative algorithms for convex constrained image recovery problems,” SIAM Journal on Imaging Sciences, vol. 2, no. 2, pp. 730–762, 2009.
• [5] M. Figueiredo and J. Bioucas-Dias, “Restoration of poissonian images using alternating direction optimization,” IEEE Transactions on Image Processing, 2010, , Submitted.
• [6] J.-J. Moreau, “Fonctions convexes duales et points proximaux dans un espace hilbertien,” CRAS Sér. A Math., vol. 255, pp. 2897–2899, 1962.
• [7] P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,” SIAM Multiscale Model. Simul., vol. 4, no. 4, pp. 1168–1200, 2005.
• [8] P. L. Combettes and J.-. Pesquet, “A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery,” IEEE J. Selec. Top. Sig. Pro., vol. 1, no. 4, pp. 564–574, 2007.
• [9] M.J. Fadili, J.-L. Starck, and F. Murtagh, “Inpainting and zooming using sparse representations,” The Computer Journal, 2006.
• [10] P. L. Combettes and J.-C. Pesquet, “A proximal decomposition method for solving convex variational inverse problems,” Inv. Prob., vol. 24, no. 6, 2008.
• [11] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” Tech. Rep., CMAP, Ecole Polytechnique, 2010.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters