Franc ¸ois-Xavier Dupé

Sandrine Anthoine

GENERALIZED SUBSPACE PURSUIT AND AN APPLICATION TO SPARSE POISSON DENOISING

Keywords: Sparse regularization, Subspace Pursuit, Poisson denoising, Greedy algorithm

We present a generalization of Subspace Pursuit, which seeks the k-sparse vector that minimizes a generic cost function. We introduce the Restricted Diagonal Property, which much like RIP in the classical setting, enables to control the convergence of Generalized Subspace Pursuit (GSP). To tackle the problem of Poisson denoising, we propose to use GSP together with the Moreau-Yosida approximation of the Poisson likelihood. Experiments were conducted on synthetic, exact sparse and natural images corrupted by Poisson noise. We study the influence of the different parameters and show that our approach performs better than Subspace Pursuit or ℓ 1 -relaxed methods and compares favorably to state-of-art methods.

INTRODUCTION

Classical greedy algorithms such as MP [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] or CoSaMP [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF] that minimize a quadratic cost under a l 0 sparsity constraint have recently been generalized to other costs. [START_REF] Zhang | Sparse recovery with orthogonal matching pursuit under RIP[END_REF][START_REF] Bahmani | Greedy sparsity-constrained optimization[END_REF] for example show that a k-sparse minimizer of a function can be found if its gradient satisfies specific properties. Furthermore these analyses lead to a better understanding of the behavior of these greedy methods. Such generalization open the doors to the use of greedy algorithms for general sparse optimization problems for which the theory is developed so far almost only when the ℓ 0 -semi-norm is relaxed to a ℓ 1 -norm.

In the context of Poisson denoising, we recently suggested a way to fill the theoretical gap by providing a CoSaMP-like greedy algorithm and its analysis [START_REF] Dupé | A greedy approach to sparse poisson denoising[END_REF]. Such effort is relevant because while some Poisson denoising algorithms have a solid theoretical ground, e.g. variance stabilization combined with Gaussian methods [START_REF] Zhang | Wavelets, ridgelets, and curvelets for Poisson noise removal[END_REF][START_REF] Makitalo | Optimal inversion of the anscombe transformation in low-count poisson image denoising[END_REF] or convex minimization with ℓ 1 -norm for sparsity [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF], others are transposition from the Gaussian world and still need a thorough analysis, e.g. non-local mean [START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF] or greedy (hybrid) algorithms [START_REF] Giryes | Sparsity based Poisson denoising[END_REF]. In this paper we pursue these efforts by providing a new greedy algorithm for sparse Poisson denoising and exhibiting sufficient conditions for its convergence. Our contributions are twofold. [START_REF] Zhang | Sparse recovery with orthogonal matching pursuit under RIP[END_REF] We propose the first (to our knowledge) generalization of the Subspace Pursuit (called GSP) to generic cost with a convergence theorem. 2) We devise a new algorithm for Poisson denoising using GSP and the same regularization procedure as in [START_REF] Dupé | A greedy approach to sparse poisson denoising[END_REF] and demonstrate its efficiency.

The structure of this paper is the following. In Section 2, we present Generalized Subspace Pursuit and introduce the Restricted Diagonal Property which ensures its convergence for a wide class of cost function (including non-convex ones). In Section 3, we use GSP with a Moreau-Yosida regularization to devise a new Poisson denoising algorithm. Experiments are conducted to see the behavior of this new method depending on the sparsity, the noise level and the regularization parameter. Finally, we compare it to state-of-art methods like SAFIR [START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF], MSVST [START_REF] Zhang | Wavelets, ridgelets, and curvelets for Poisson noise removal[END_REF], the classical ℓ 1 -norm relaxation method and Subspace Pursuit (SP) [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF].

GENERALIZED SUBSPACE PURSUIT

Classical greedy algorithms such as MP, CoSaMP or SP attempt to find the "best" k-sparse vector x that explains the observed measured vector y under a linear model: y = Φx, where "best" is understood in the ℓ 2 sense which means that these algorithms solve the following problem:

Find x ∈ argmin x y -Φx 2 2 s.t. x 0 ≤ k . (P1) (||x|| 0 = Card{i : x i = 0}
is the number of non-zero entries of x). This can also be viewed as solving the linear inverse problem y = Φx under the assumption of Gaussian noise and with a hard sparsity constraint as regularization.

In many inverse problems, the quadratic cost y -Φx 2 is not appropriate and one rather wishes to solve:

Find x ∈ argmin f (x) s.t. x 0 k , ( P2
)
for a more complex f . A vast literature gives solutions for a relaxed sparsity regularization (i.e. using the ℓ 1 -norm instead of the ℓ 0 -semi-norm) in the case of convex f for example. In contrast, only few attempts have successfully extended greedy algorithms to Problem (P2). Indeed the convergence proofs of MP, SP etc. heavily rely on linear algebra useful because the gradient of y -Φx 2 2 is affine (see [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF][START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF][START_REF] Jain | Orthogonal matching pursuit with replacement[END_REF]), which is not the case for a generic f .

Algorithm

We propose an extension of Subspace Pursuit [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] designed to solve Problem (P2) for a cost function f : R n → R differentiable. This algorithm named Generalized Subspace Pursuit (GSP) is detailed in Algorithm 1.

Here, we denote by supp(x) the support of x: supp(x) = {i = 1, . . . , n | x i = 0} and x k the restriction of x to its k largest coefficients in magnitude.

As SP, GSP seeks the support of the solution of size k (denoted by T ). It updates it at each iteration by adding a set of k candidate directions (step 3-5), obtaining an extended support set S of size 2k. The problem is solved on this extended support (step 6). T is updated by the best support of size k extracted from this solution (step 7). Finally, one seeks the best solution of support T (step 8).

GSP exactly mimics SP, except for steps 3, 6 and 8 where it accounts for f instead of y -Φx 2 2 . Steps 6 and 8 consist in finding an optimal solution to the original problem on a predefined support set, which can usually be done efficiently (e.g. with a projected gradient descent).

Algorithm 1 Generalized Subspace Pursuit 1: initialization: â0 ← 0 . 2: for t = 0 . . . N -1 do 3: g ← ∇f (â t ) , ⊲ Compute local directions. 4:
G ← g k , ⊲ Select new directions.

5:

S ← G ∪ supp(â t ) , ⊲ Set extended support. 9: end for 10: outputs: âN .

Theoretical guarantees

In the classical case (f = ||Φx -y|| 2 ), a sufficient condition for MP, SP or CoSaMP to converge is that the matrix Φ satisfies the RIP condition with a small constant δ k . This essentially ensures that Φ acts is almost isometry on the set of k-sparse vectors so that Φ * Φ which controls the increments ∇f (x)-∇f (x ′ ) is close to the identity. It turns out that in the generalized case, ensuring the increments of ∇f are close to the identity on the set of k-sparse vectors also guarantees the convergence of GSP. This condition can even be weakened by replacing the identity by any diagonal operator bounded away from zero and having this property hold only locally.

To state the convergence of GSP, let us introduce D 1 the set of diagonal operators bounded away by 1:

D 1 = {D : R n → R n , D diagonal, ∀x Dx x }. (1)
Definition 1 (Restricted Diagonal Property). T is said to have the Restricted Diagonal Property (RDP) of order k if there exists

δ k > 0 such that if a, b ∈ R n then, Card{supp(a) ∪ supp(b)} k ⇒ ∃ D ab ∈ D 1 , T(a) -T(b) -D ab (a -b) δ k a -b . ( 2
)
Remark 1. δ 2k < 1 implies that T is injective on the set of k-sparse vectors since T(a) -T(b) (1 -δ 2k ) x -y .
Remark 2. If D ab is the Identity for all a and b, there is a connection with RIP since then (1

-δ 2k ) x -y T(x) -T(y) (1 + δ 2k ) x -y .
However the Restricted Diagonal Property is a more local property. Let us now state the convergence theorem.

Theorem 2. Assume f : R n → R, differentiable and (P2) has at least one solution. Denote by a ⋆ any solution of (P2) and δ 0 the unique real root of g(x) = x 3 + x 2 + 7x -1.

If there exists β > 0 such that ∇f β has the Restricted Diagonal Property of order 3k with δ 3k δ 0 . Then GSP verifies

âN -a ⋆ 1 2 N â0 -a ⋆ + 15 β ∇f (a ⋆ ) 2k . ( 3
)
The distance between the iterates of GSP and a solution of (P2) are thus controlled by two terms. The first one decreases exponentially fast to zero while the second expresses the error made, and is controlled by the gradient. Theorem 2 thus shows that if f admits global minimizer that is k-sparse, GSP recovers it exactly (since in that case ∇f (a ⋆ ) = 0): this is the analog to the exact recovery case for SP, MP and the like. Otherwise, the quality of the iterates is controlled by the error term ∇f (a ⋆ ) 2k , which is analog the the uncompressible error term exhibited in the classical case (see [START_REF] Needell | CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[END_REF]).

The proof of this theorem (not given here due to lack of space) is very similar to that of [START_REF] Dupé | A greedy approach to sparse poisson denoising[END_REF], where a generalization of CoSaMP is proposed and analyzed using a stronger property.

Let us pinpoint the pros and cons of the proposed analysis. First, this is the first generalization of SP to non-quadratic costs, and we provide a convergence theorem. An advantage of extending SP versus CoSaMP (as in GRASP [START_REF] Bahmani | Greedy sparsity-constrained optimization[END_REF][START_REF] Dupé | A greedy approach to sparse poisson denoising[END_REF]) lies in GSP's last step (step 8), which consists in finding the best k-sparse vector with the updated support T of size k. This is more costly than the simple pruning used in CoSaMP or GRASP, but has been empirically proved to speed up the convergence in the classical setting (SP vs CoSaMP). More importantly, this enables to consider cost functions f incorporating hard constraints since those will be naturally respected at all steps (not shown here).

Second, although the RDP is not easy to check in practice, it is clear that it does not require f to be convex (for example f (x) = x T (D + γI)x with D = diag(±1, . . . , ±1), so that D ab = D and δ k = γ). So the proposed analysis extends the scope of greedy algorithms to non-convex costs, as opposed to previous analyses like [START_REF] Bahmani | Greedy sparsity-constrained optimization[END_REF].

In the next Section we shall apply GSP to the problem of Poisson denoising, which is a typical example where the quadratic cost function is not adapted.

APPLICATION TO POISSON DENOISING

Poisson denoising under a hard sparsity constraint

Let us assume that we observe y ∈ R n , a Poisson corrupted version of the true image x ∈ R n , both containing n pixels. In our model, we also assume that x has a k-sparse representation on the dictionary Φ = (ϕ 1 , . . . , ϕ m ) ∈ R n×m :

x = Φα = α j ϕ j with ||α|| 0 = k ≪ n ,
where the atoms are normalized ( ϕ j = 1).

Our goal is to reconstruct x given the data y, the sparsity k and the dictionary Φ, which may be done by solving:

x = Φα, where α = argmin

α∈R m s.t. ||α||0 k F y (Φα) , (P3)
where F y (x) is a data fidelity term that quantifies how well an estimated image x fits the observed data y. A natural fidelity term is the negative-log-likelihood F y (x) =log P(y|x) which reads in the case of Poisson noise

F y (x) = -log P(y|x) = n i=1 f (x i , y i ), with f (ξ, η) =      -η log(ξ) + ξ if η > 0 and ξ > 0, ξ if η = 0 and ξ 0, +∞ otherwise. (4)
Notice that F y (x) is finite only when x complies with the data, which implies x ∈ R n + and x i > 0 if y i > 0. Moreover, due to the logarithm, it gradient is not defined on its all domain. Hence GSP may not be applied on F y • Φ. These technical difficulties reflect a deeper conceptual problem in solving Problem (P3): in fact the hard constraints imposed by the Poisson negative-log-likelihood on Φα ((Φα) i 0 for all i and (Φα) i > 0 if y i > 0) may not be compatible with hard ℓ 0 constraint α 0 k.

To alleviate this problem, we substitute in (P3) F y • Φ for its Moreau-Yosida regularization [START_REF] Lemaréchal | Practical Aspects of the Moreau-Yosida Regularization: Theoretical Preliminaries[END_REF] of parameter λ:

M λ,Fy•Φ (α) = inf s∈R m 1 2λ ||s -α|| 2 + F y • Φ(s) . (5)
The Moreau-Yosida envelope of a convex l.s.c. function has full domain, its gradient is Lipschitz. Moreover, its minimizer is also minimizer of the original function. It does provides a good candidate for relaxing F y • Φ and using it with GSP. To do so, we need the gradient of this function, which is computable when the dictionary Φ is a tight frame.

Proposition 3 (Gradient of the Moreau-Yosida regularization of the Poisson neg-log-likelihood). If Φ is a tight frame of constant ν > 0, then the gradient of M λ,Fy•Φ is [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF]:

∇M λ,Fy•Φ = 1 νλ Φ * • (I -prox νλFy ) • Φ with prox νλFy (x) i = x i -νλ + |x i -νλ| 2 + 4νλy i 2 . (6)
In the following, we thus denoise y using GSP to solve

x = Φα, where α = argmin α s.t. ||α||0 k M λ,Fy•Φ (α) . (P4)
It is not clear whether M λ,Fy•Φ has the RDP, however we empirically show in the next section that the proposed algorithm converges.

Experiments

The strengths and the weaknesses of our approach are assessed on three kinds of experiments. For each one, the quantitative values are averaged over 10 realizations of the noise. The dictionary Φ is fixed to cycle spinning wavelets with a redundancy of 9. In the first experiment, we study the evolution of the support recovery with the noise and the sparsity levels. Next, we study the influence of the parameter of the Moreau-Yosida regularization on the recovered support. Finally, we compare the performance of our method to that of state-of-the-art ones.

Support recovery

In this experiment, we generate random k-sparse images of size 64 × 64 with different sparsities and peak intensities. They are subsequently corrupted by Poisson noise. Each image is denoised using our method with the regularization parameter

λ = 1. Fig 1(a)
shows the phase diagram of support recovery (number of atoms of the true support recovered divided by k) as a function of sparsity level (k/n) and peak intensity (higher intensity means lower noise level). As expected when the image is sparser or less degraded (lower left), the support is easier to find.

Influence of the regularization

Let us asses the influence of the Moreau-Yosida parameter (λ) on support recovery. Fig. 1(b) pictures the evolution of support recovery with λ for a sparse version of the Cameraman (k/n = 0.15) with a peak intensity at 11 obtained by Subspace Pursuit. This illustrates the following facts. First, the lower the λ is, the better we reconstruct the support. This is expected as when λ tends to 0, the Moreau-Yosida regularized function tends to the true function. Second, the quality degrades for high λ, eventually becoming worse than that of Subspace Pursuit (0.79), which shows that too much regularization degrades the estimation. Finally, the plateau reached for small λ indicates that λ does not need to be specifically tuned but only small enough. For the next experiments, we choose λ to obtain the best support recovery.

Performances

We compare our approach (GSP) with other states-of-art methods: Subspace Pursuit (SP) [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF], SAFIR [START_REF] Boulanger | Patch-based nonlocal functional for denoising fluorescence microscopy image sequences[END_REF] (with the parameters from [START_REF] Makitalo | Optimal inversion of the anscombe transformation in low-count poisson image denoising[END_REF]), MSVST [START_REF] Zhang | Wavelets, ridgelets, and curvelets for Poisson noise removal[END_REF] and a convex ℓ 1 -relaxation of (P3) i.e. a procedure minimizing the Poisson likelihood on a ℓ 1 -ball (using a projection onto the ℓ 1 -ball [START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators[END_REF][START_REF] Chierchia | Epigraphical projection and proximal tools for solving constrained convex optimization problems: Part i[END_REF]). We apply these methods on two images, a sparse version of the Cameraman (see 3.2.2) and the NGC 2997 galaxy (peak intensity at 5, see Fig. 2). We use the exact parameters (sparsity for GSP and SP and ℓ 1 -norm for the ℓ 1 -relaxation) for the sparse cameraman and tune them for the Galaxy. For each method we compute the PSNR, the mean absolute deviation (MAE) and the SSIM [START_REF] Wang | Image quality assessment: From error visibility to structural similarity[END_REF] and display them in Table 1.

The MSVST gives the best results in terms of PSNR and MAE and SAFIR does so for the SSIM. GSP is competitive with both these state-of-the-arts methods in Poisson denoising. Moreover GSP is in both cases better than SP and the ℓ 1 -relaxation. While the performance over SP shows the benefits of using non-quadratic cost in the case of high Poisson noise level, the performance over the ℓ 1 -relaxation emphasizes that using the exact sparse ℓ 0 a-priori enables to recover a better dynamic that its convex relaxation counterpart. These remarks are visually confirmed for the galaxy (Fig. 2).

CONCLUSION

In this paper we proposed a greedy algorithm, Generalized Subspace Pursuit, seeking a k-sparse minimizer for a generic cost function f . We exhibited a property on the gradient of f , the Restricted Diagonal Property that controls its convergence. Although hard to verify in practice, this property opens the way for using greedy algorithms on non-convex functions. A method of Poisson denoising is proposed using GSP on the Moreau-Yosida regularization of the Poisson likelihood. The experiments show that it is competitive to state-of-thearts methods and improves over classical sparse approaches.

Future work on the theoretical side include improving the constants in the theorem, which seem quite conservative in view of the conducted experiments; studying possible links to notions used in the linear case, such as coherence; and adapting the analysis to the co-sparse [START_REF] Nam | The cosparse analysis model and algorithms[END_REF] case. One may also want to reduce the complexity of GSP by removing one of the minimization steps, preserving the last one (step 8), as is done in OMPR [START_REF] Jain | Orthogonal matching pursuit with replacement[END_REF].

A further analysis of the regularization of the Poisson likelihood needs to be conducted to ensure the convergence of the proposed denoising method. Furthermore, it would be interesting to consider learning a dictionary to improve the quality of the results. However, such dictionaries are not likely to be tight frames or even frames and would make the previous analysis even more complex.

6 : 8 :

 68 b ∈ argmin {a:supp(a)⊆S} f (a) , ⊲ Solve on extended support. 7: T ← supp(b k ) , ⊲ Set support. ât+1 ∈ argmin {a:supp(a)⊆T } f (a) , ⊲ Solve on the support.

Fig. 1 .

 1 Fig. 1. Support recovery as a function of sparsity and noise level (a) and of the Moreau-Yosida regularization (b).

Fig. 2 .

 2 Fig. 2. NGC 2997 Galaxy image (a), a noisy version (b) and several denoising results (c-g).

Table 1 .

 1 Comparison of denoising methods on a sparse version of Cameraman (k/n = 0.15) and the NGC 2997 Galaxy.

		Sparse Cameraman		Galaxy
		PSNR	MAE SSIM	PSNR MAE	SSIM
	Noisy	19.3	1.57	0.32	22.9	0.63	0.19
	GSP	28.6	0.32	0.87	26.0	0.17	0.71
	SP [11]	26.1	0.55	0.63	27.6	0.28	0.55
	SAFIR [9]	27.6	0.36	0.86	29.9	0.15	0.84
	MSVST [6]	29.5	0.31	0.84	30.5	0.12	0.83
	ℓ 1 -relaxation	22.7	0.64	0.73	21.5	0.32	0.50

This work is partially supported by the French GIP ANR under contract ANR GRETA 12-BS02-004-01 Greediness: theory and algorithms.