Inverse Problems with Poisson noise: Primal and PrimalDual Splitting^{1}^{1}1Submitted 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 (loglikelihood) 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 nonsmooth convex functional. We establish the wellposedness of the optimization problem, characterize the corresponding minimizers, and solve it by means of primal and primaldual proximal splitting algorithms originating from the field of nonsmooth convex optimization theory. Experimental results on deconvolution and comparison to prior methods are also reported.
Inverse Problems with Poisson noise: Primal and PrimalDual Splitting^{†}^{†}Submitted to ICIP 2011 on the 01/21/11.
F.X. Dupé, M.J. Fadili and J.L. Starck  

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 wellknown RichardsonLucy 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 waveletregularized RL algorithm.
In the context of Poisson linear inverse problems using sparsitypromoting 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 (ForwardBackward and DouglasRachford) to solve the corresponding optimization problem. The work of [4] is in the same vein. However, nested algorithms are timeconsuming since they necessitate to subiterate. Using the augmented Lagrangian method with the alternating method of multipliers algorithm (ADMM), which is nothing but the DouglasRachford splitting applied to the FenchelRockafellar dual problem, [5] presented a deconvolution algorithm with TV and sparsity regularization. This scheme however necessitates to solve a leastsquare 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 nonsmooth convex functional, for which we prove wellposedness of the optimization problem, characterize the corresponding minimizers, and solve them by means of primal and primaldual proximal splitting algorithms originating from the realm of nonsmooth 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 realvalued function is coercive, if , and is proper if its domain is nonempty . 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
(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 prechosen dictionary as measured by a sparsitypromoting penalty supposed throughout to be convex but nonsmooth, e.g. the norm.
From the probability density function of a Poisson random variable, the likelihood writes:
(2) 
Taking the negative loglikelihood, we arrive at the following data fidelity term:
(3)  
Our aim is then to solve the following optimization problems, under a synthesistype sparsity prior^{1}^{1}1Our framework and algorithms extend to an analysistype prior just as well, though we omit this for obvious space limitation reasons.,
() 
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.
Proposition 1.

is a convex function and so are and .

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

Suppose that . Then .
2.1 Wellposedness of ()
Let be the set of minimizers of problem (). Suppose that is coercive. Thus is coercive. Therefore, the following holds:
Proposition 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 loglikelihood) is,
(4) 
We now turn to which is given by Lemma 4 and the following result:
Theorem 6 ([9]).
Suppose that : (i) is convex evensymmetric, nonnegative and nondecreasing 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 :
(5) 
Among the most popular penalty functions satisfying the above requirements, we have , in which case the associated proximity operator is softthresholding, 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:
(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 DouglasRachford 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:
(7)  
(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 wellknown 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.
3.3 Splitting on the dual: Primaldual algorithm
Our problem () can also be rewritten in the form,
(10) 
where now and . Again, one may notice that the proximity operator of can be directly computed using the separability in and .
3.4 Discussion
Algorithm 1 and 2 share some similarities, but exhibit also important differences. For instance, the primaldual 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 stateoftheart 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.
RLMRS [2]  RLTV [1]  StabG [3]  PIDALFS [5]  Alg. 1  Alg. 2  
MAE  63.5  52.8  43  43.6  46  43.6 
Times  230s  4.3s  311s  342s  183s  154s 
5 Conclusion
In this paper, we proposed two provably convergent algorithms for solving the Poisson inverse problems with a sparsity prior. The primaldual 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 stateofart 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. BioucasDias, “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 forwardbackward splitting,” SIAM Multiscale Model. Simul., vol. 4, no. 4, pp. 1168–1200, 2005.
 [8] P. L. Combettes and J.. Pesquet, “A DouglasRachford 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 firstorder primaldual algorithm for convex problems with applications to imaging,” Tech. Rep., CMAP, Ecole Polytechnique, 2010.