Manon Lamy

Kamal Hammani
email: kamal.hammani@u-bourgogne.fr

Juan Arocas

Christophe Finot

Jean-Claude Weeber

Broadband etching-free metal grating couplers embedded in titanium dioxide waveguides

Embedded metal grating couplers into titanium dioxide layer are proposed. A coupling efficiency better than 20 % is experimentally demonstrated with a 3-dB bandwidth of 86 nm in agreement with simulation results. This allowed us to perform error-free transmissions of 10 Gbit/s wavelength multiplexed signals in the C-band. ©2017 Optical Society of America

The efficient coupling of light from an optical fiber to photonic waveguides remains an essential need for the emergence of integrated optical devices. The excitation of an integrated waveguide by an optical fiber is basically achieved by two ways known as end-fire coupling (or butt-coupling) or grating coupling. In the case of end-fire coupling, the optical fiber must be carefully aligned with respect to the (clean) waveguide input and sophisticated approaches such as cantilever adiabatic couplers [1] are needed to account for the large mismatch between the optical fiber and the waveguide mode size. On the other hand, grating couplers (GC) lead to very high coupling efficiencies [2] without the requirement of a highly accurate positioning of the optical fiber but at the cost of tedious fabrication processes. In particular, the etching depth of the grating couplers turns to be a critical parameter that can fully jeopardized their efficiency. To circumvent those limitations, couplers relying on metallic gratings have been recently proposed [3,4]. In those previous works, the grating couplers were deposited either on top of or beneath the photonic layer limiting the number of structural parameters available for the optimization of the couplers.

In this study, we investigate the situation of a metallic grating coupler embedded within the photonic guiding layer comprised of titanium dioxide deposited on a glass substrate. In particular, we optimize the location of the grating into the photonic waveguide and we show experimentally a factor two improvement of the coupling efficiency for the embedded grating as compared to the situation of the coupler deposited on top of the waveguide. Beyond their simple coupling efficiency characterizations, we also show that the embedded couplers are fully convenient for error-free transmission of 10 Gbit/s WDM signals. A diagrammatic view of the embedded GCs we consider is depicted in Fig. 1(a). The GC is comprised of gold ridges (width w1, height hAu) periodically distributed (period Λ) and embedded into a titanium dioxide (TiO2) layer lying on a glass substrate. In our configuration, gold has been selected for its chemical stability to prevent alteration of the ridges during the encapsulation process. However, it is well-know that gold ridges can support localized surface plasmon resonances (LSPRs). The interaction of photonic modes with periodic arrangement of metallic nanostructures sustaining LSPRs has attracted much interest recently [5][6][7]. The spectral shaping of LSPRs and the generation of surface lattice resonances originating from the hybridization of the LSPR and the photonic modes have been demonstrated in this context. In the present case, the cross-sections of the gold ridges are such that they do not support LSPR in the telecom frequency range we consider. Hence, although of great interest, the contribution LSPR in metallic grating couplers is out of scope of the present study. Given that our GCs do not exploit the plasmonic properties of the metal, one can anticipate that cost effective CMOS compatible metals such as titanium nitride or aluminum could be conveniently used in place of gold. Figures 1(b) and 1(c) show an atomic force microscope image and a SEM image of the fabricated GC respectively. The GCs are sandwiched in between two layers of TiO2 with a total fixed height of 255nm throughout the study. TiO2 has indeed been recognized [8][9][10][11] as a promising material for linear and non-linear integrated optics applications due to its high refractive index (n = 2.45 at 1.55 µm), broad transparency window and negligible two photon absorption in the telecom C-band. TiO2 can be readily deposited at low temperature making this material very convenient for the encapsulation of metallic GCs.

For a fixed thickness of the photonic waveguide, the period, the filling factor, the height and the distance hbottom [see Fig. 1(a)] have to be considered for an optimum design. As discussed in more details in the following, the encapsulation of the GC is obtained by a conformal deposition of a TiO2 layer onto the metal grating in such a way that a replica of the metal grating is created at the top surface of the upper TiO2 layer. We carried out an optimization of the structural parameters of two-dimensional GCs by using a commercial finite-element software. Assuming a Gaussian beam (free-space wavelength 1550nm) falling at an angle of incidence of Θ = 30° (corresponding to our experimental configuration), the best set of parameters (giving the maximum coupling efficiencies) for a fixed total thickness of 255 nm have been obtained thanks to Monte-Carlo routine and are summarized in Table1.

Table 1. Best set of parameters and coupling efficiency (CE) depending on the configuration and on the initial polarized

Gaussian beam

CE ~20 % (~-7dB) ~ 15 % (~-8dB)
The embedded configuration allows to reach higher coupling efficiencies for both TE and TM polarization compared to the top configuration. Moreover, while the difference between TE and TM polarization in top configuration is not significant, the TM polarization is clearly more efficient than the TE one with the embedded configuration. Therefore, in the following, we focus our attention on the TM polarization.

Figure 2 provides an overview of the impact of each parameter on the coupling efficiency per facet. When the period increases, as expected the maximum efficiency is shifted toward longer wavelength [Fig. 2(a)] and the highest coupling is obtained for filling factors close to 50 % [Fig. 2(b)]. Regarding the impact of the position of the embedded grating, it is clear from Fig. 2(c) that the coupling efficiency is higher when the grating is not directly on the substrate [3]. Actually, the depth of embedded grating is related to the mode field distribution of the unperturbed waveguide. More precisely, the grating has to be positioned where the power flow (i.e. the Poynting vector) is the highest into the layer. Another variable parameter is the height of the grating which is not critical for the coupling efficiency beyond a minimal thickness. However, when the height of the grating increases, the optimal wavelength is shifted towards higher wavelength [Fig. 2(d)] due to an increase of the effective index. Moreover, beyond a certain value (here 56nm), the height of the grating becomes deleterious and leads to lower efficiency. Based on those results, embedded GCs have been obtained using the following process flow. First, a TiO2 layer (thickness 69nm) is deposited on a glass substrate by DC magnetron sputtering of a 99.9% pure titanium target. Next, the gold gratings are fabricated on top of the TiO2 layer by standard electron-beam lithography, thermal metal evaporation (3nm of chromium used as an adhesion layer+56nm of gold) followed by a lift-off process. Then the metal gratings are encapsulated by depositing a second TiO2 layer (thickness 193nm). An overlay electron lithography using alignment marks and eventually a dry-etching of the TiO2 layers are conducted to fabricate the photonic structures (waveguides and tapers). We have fabricated several devices with the period varying from 1300 nm to 1500 nm by step of 25 nm. Finally, the width of the grating is about 610 nm (targeted 666 nm). Note that two coupling tapers, where the gratings are embedded, are used to achieve beam spot conversion into a 1.5 µm-wide rib waveguide. The total length of the device including these 85-µm long tapers is 575 µm.

The experimental setup is sketched in Fig. 3(a). Optical fibers equipped with focusers are used to couple light in and out the waveguide. The focusers generate a spot with a diameter around 10 µm at a working distance of 8 mm and are tilted at an angle of 30° with respect to the normal to the substrate. An IR and a visible camera mounted on a vertically oriented microscope are operated simultaneously to align accurately the focuser onto the GCs (see Figs 3(b) and 3(c) respectively). The light source we use is generated from a high peak power picosecond laser, centered at 1550 nm injected into a highly nonlinear fiber with anomalous dispersion. This apparatus generates a supercontinuum [12] with a spectrum ranging from 1300 nm to 1900 nm. The polarization of the light injected into the input focuser is adjusted by polarization controllers and is oriented to be TM. At the distal end of the waveguide, the transmitted light decoupled by the grating is collected by the output focuser and analyzed with an optical spectral analyzer (OSA) or an accurate powermeter (PM). For the quantitative measurement of the coupling efficiencies, a thick gold layer on a glass substrate is used as a highly reflecting surface. The light collected by the output focuser after specular reflection of the incident light onto the gold mirror is optimized and a reference spectrum is recorded. By normalizing the transmission spectra along the waveguide to the reference spectrum, the spectral dispersion of the coupling efficiency for each grating coupler is readily obtained. Taking into account the symmetry of the coupling in and the coupling out configuration in our setup and neglecting the confinement and material losses experienced upon propagation, the coupling efficiency per facet is defined as half the transmission loss along the entire waveguide.

Figure 4(a) displays the coupling efficiency per facet for GCs with period around 1400 nm (optimal period in simulations). As expected from Fig. 2(a), we observe a red shift of the optimal wavelength with the increase of the GC period. Figure 4(b) compares the numerical and experimental results obtained for the optimal experimental period Λ = 1375 nm. The dashed curve shows the optimized coupling efficiency that can be obtained after a full optimization of the GC parameters. The dotted grey curve is obtained by using for the numerical modeling the actual parameters (hbottom = 69 nm, htop=193 nm, w = 610 nm, hAu= 56 nm) measured after fabrication of the GC whereas the solid line displays the experimental coupling efficiencies obtained for this GC. The experimental coupling efficiency per facet features a maximum at 1550 nm around -7dB (20 %) and a 3 dB-bandwidth of 86 nm (1 dB bandwidth= 50 nm) in fairly good agreement with the numerical results obtained with the actual parameters of the GC. We note however a decrease of about 2dB between the experimental and the numerical results (black dashed curve) attributed to the propagation losses along the waveguide and mainly related to fabrication defects [see scattering spots in Fig. 3(b)]. In addition, by referring to the optimum performances of the GC, we note that an improvement by about 2 dB of the coupling efficiency per facet is expected by a further optimization of the GC parameters. This waveguide is polarization-sensitive, with a measured extinction ratio of 12.9 dB. To go further, we have fabricated GCs deposited on top of the photonic waveguides [Fig. 5(a)] at the same time (with the same lithography design). Thus the period and the filling factor of the GC are identical as in the case of the embedded gratings. Therefore, the TM mode is not optimally coupled due to a wrong filling factor. More precisely, while the optimal coupling (see Table 1) is expected to be -8 dB, with our set of experimental parameters, the coupling efficiency should drop down to -13 dB per facet. In our experiment, we actually measured a total loss close to -30 dB which means a coupling efficiency being about -15 dB per facet, validating qualitatively the predictions of the numerical simulations. Regarding the TE mode, the same grating with a filling factor close to 50 % can be optimal but of course not at 1550 nm. Indeed, the effective index is significantly higher and leads to a shift towards higher wavelength. Thus by operating the supercontinuum source, the spectra displayed in Fig. 5(b) were obtained. The best result was about -10dB (where the simulation gives -7 dB) which is still lower than the efficiency with the embedded grating coupler (expected to be -5dB). This demonstrates the potential of such embedded structures to couple light in photonic waveguides and confirms that the embedded GCs are a particularly powerful way to couple the TM-polarized light.

The superior performances of embedded GCs as compared to deposited GCs can be understood qualitatively from the fact that the grating is located at the maximum of the power density of the unperturbed waveguide but also from the interaction of the incident light with the grating located at the upper surface of the photonic waveguide. The role of this upper grating can be evaluated on the basis of numerical modeling. Indeed, for the optimal GC we have checked that the coupling efficiency of the embedded coupler is reduced by a factor 2 if the top grating is removed. The interaction of the top grating and embedded grating is thus a key ingredient for the optimum performance of the embedded GCs. The intrinsic properties of the embedded GCs have been discussed so far but it is also worthy to demonstrate those GCs in the context of practical applications such as a high-bit rate traffic. To demonstrate the efficiency of the embedded coupler device, we have transmitted a 10 Gbit/s signal at standard telecommunications wavelengths. More precisely, we used enhanced small form-factor pluggable (SFP+) to transmit and receive light. Those transceivers are commonly used as part of wavelength-division multiplexing optical networks. The 10 Gbit/s signal is coded through a 2 31 pseudo random bit sequence (PRBS) generator. At the output, a bit error rate tester (BERT) and a high speed sampling oscilloscope are used to evaluate the bit error rate (BER) and monitor the eye diagram. For SFP+ typical of coarse wavelength multiplexing emitting at 1530, 1550 and 1570 nm, BER lower than 10 -12 are achieved, indicating error-free transmission. For SFP+ centered at 1510 nm, the BER decreases to 1. 6 10 -11 due to 2 dB extra-loss that prevents the integrated amplifier and reshaper to perfectly work. More precisely, at 1550 nm, the transmission is still error-free whereas at 1570 nm the BER is about 7.5 .10 -7 . The degradation between single and multi channel configurations is explained by the insertion of two optical wavelength multiplexers leading to 2 dB of extra-loss.

To conclude, we have shown, for the first time to our knowledge, the fabrication of metal grating couplers in gold embedded into a titanium dioxide waveguide. We demonstrate coupling efficiency greater than 20 % at 1550 nm for a TM polarization. The bandwidth of this device is about 86 nm at 3 dB. The comparison with a metal grating coupler on top of the waveguide confirms the higher potential of the embedded device. Our simulations are in good agreement with our numerical results even though small discrepancies can be explained by fabrication defects. However, simulations demonstrated that efficiencies close to 45 % could be reached. The device has been tested in real configuration using integrated transceivers (SFP+) to transmit coded signal at 10 Gbit/s without any error. To increase the performances of such device, several ways should be explored with in particular a metallic pad on the substrate that acts as a mirror or the optimization of coupler shape by using parabolic grating [4]. Such coupling efficiency will allow to consider the titanium dioxide platform as a serious candidate for nonlinear applications in near infrared [9].

Fig. 1 .

 1 Fig. 1. (a) Embedded metal grating in titanium Dioxide layer (side view) (b) Atomic Force Microscope picture of the top surface. (c) SEM image of a fabricated device (top view)

Fig. 2 .

 2 Fig. 2. Coupling efficiency as a function of (a) the period, (b) the filling factor, (c) the position and (d) the thickness of the metal embedded grating. For non-varying parameters, the values are Λ = 1400 nm, w1 = 666 nm, hAu= 56 nm, hbottom= 65 nm and htop= 255 nm -hbottom.

Fig. 3 .

 3 Fig. 3. (a) Experimental Setup. PC: polarization controller, OSA: optical spectrum analyzer, PM: powermeter, CWL: continuous wave laser, HNLF: highly nonlinear fiber. Top view of under test 1.5 µm-wide waveguide with (b) IR or (c) visible camera.

Fig. 4 .

 4 Fig. 4. Coupling efficiency per facet of embedded grating coupler as a function of the injected wavelength (a) for several period. (b).The best experimental results (=1375nm, solid black curve) are compared with numerical simulations taking into account real parameters (dotted grey curve) and optimum parameters (dashed black curve).

Fig. 5 .

 5 Fig. 5. (a) Top metal grating on titanium dioxide (side view). (b) Coupling efficiency per facet of embedded grating coupler as a function of the injected wavelength for several periods.

Fig. 6 (

 6 Fig. 6 (a) Example of experimental setup for 10 Gbit/s transmission. Each SFP+ is composed of a transmitter Tx and a receiver Rx. WDM = wavelength division multiplexer. The square called 'set-up' corresponds to the setup shown in Fig. 3(a). (b) Optical spectrum in copropagation configuration (c) Eye diagrams at 1550 nm (c1) and 1570 nm (c2) in counter-propagation configuration. To go beyond, two WDM experiments have been performed. First, a co-propagating configuration [Fig. 6(a)] has been considered with a signal centered at 1550 nm and another

Funding.

  PARI PHOTCOM Région Bourgogne, Labex Action (ANR-11-LABX-01-01) and FEDER-FSE Bourgogne 2014/2020.

Acknowledgment. We thank J. Fatome and M. G. Nielsen for advices and fruitful discussion on experimental and numerical parts.


Renliang Gu
email: renliang@iastate.edu

Senior Member, IEEE Aleksandar Dogandžić

Projected Nesterov's Proximal-Gradient Algorithm for Sparse Signal Recovery

Keywords: Convex optimization, Nesterov's momentum acceleration, sparse signal reconstruction, Poisson compressed sensing, proximal-gradient methods

We develop a projected Nesterov's proximal-gradient (PNPG) approach for sparse signal reconstruction that combines adaptive step size with Nesterov's momentum acceleration. The objective function that we wish to minimize is the sum of a convex differentiable data-fidelity (negative log-likelihood (NLL)) term and a convex regularization term. We apply sparse signal regularization where the signal belongs to a closed convex set within the closure of the domain of the NLL; the convex-set constraint facilitates flexible NLL domains and accurate signal recovery. Signal sparsity is imposed using the `1-norm penalty on the signal's linear transform coefficients. The PNPG approach employs a projected Nesterov's acceleration step with restart and a duality-based inner iteration to compute the proximal mapping. We propose an adaptive step-size selection scheme to obtain a good local majorizing function of the NLL and reduce the time spent backtracking. Thanks to step-size adaptation, PNPG converges faster than the methods that do not adjust to the local curvature of the NLL. We present an integrated derivation of the momentum acceleration and proofs of O.k 2 / objective function convergence rate and convergence of the iterates, which account for adaptive step size, inexactness of the iterative proximal mapping, and the convex-set constraint. The tuning of PNPG is largely application independent. Tomographic and compressed-sensing reconstruction experiments with Poisson generalized linear and Gaussian linear measurement models demonstrate the performance of the proposed approach.

I. Introduction

Most natural signals are well described by only a few significant coefficients in an appropriate transform domain, with the number of significant coefficients much smaller than the signal size. Therefore, for a vector x 2 R p that represents the signal and an appropriate sparsifying dictionary matrix ‰, ‰ H x is a signal transform-coefficient vector with most elements having negligible magnitudes. Real-valued ‰ 2 R p p 0 can accommodate discrete wavelet transform (DWT) or gradient-map sparsity with anisotropic total-variation (TV) sparsifying transform (with ‰ D OE‰ v ‰ h ); a complex-valued ‰ D ‰ v C j‰ h 2 C p p 0 can accommodate gradient-map sparsity and the 2D isotropic TV sparsifying transform; here ‰ v ; ‰ h 2 R p p 0 are the ver- tical and horizontal difference matrices similar to those in [START_REF] Boyd | Vectors, Matrices, and Least Squares[END_REF]Sec. 15.3.3]. The idea behind compressed sensing [START_REF] Candes | Near-optimal signal recovery from random projections: Universal encoding strategies?[END_REF] is to sense the significant components of ‰ H x using a small number of measurements; here, " H " denotes the conjugate transpose.

We use the negative log-likelihood (NLL) (data-fidelity) function L.x/ to describe the noisy measurement process. Consider signals x that belong to a closed convex set C and assume

C Â cl.dom L/ (1)
which ensures that L.x/ is computable for all x 2 int C . If C n dom L is not empty, then L.x/ is not computable in it, which needs special attention; see Section III. The nonnegative signal scenario with

C D R p C ( 2
)
is of significant practical interest and applicable to X-ray computed tomography (CT), single photon emission computed tomography (SPECT), positron emission tomography (PET), and magnetic resonance imaging (MRI), where the pixel values correspond to inherently nonnegative density or concentration maps [START_REF] Prince | Medical Imaging Signals and Systems[END_REF]. Harmany et al. consider such a nonnegative sparse signal model and develop in [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms-theory and practice[END_REF] and [START_REF] Harmany | Gradient projection for linearly constrained convex optimization in sparse signal recovery[END_REF] a convex-relaxation sparse Poisson-intensity reconstruction algorithm (SPIRAL) and a linearly constrained gradient projection method for Poisson and Gaussian linear measurements, respectively. In addition to signal nonnegativity, other convex-set constraints have been considered in the literature: prescribed value in the Fourier domain; box, geometric, and total-energy constraints; intersections of these sets [START_REF] Youla | Image restoration by the method of convex projections: Part 1-theory[END_REF]; and unit simplex [START_REF] Hansen | Jackknife model averaging[END_REF].

We adopt the analysis regularization framework and minimize f .x/ D L.x/ C ur.x/ (3a)

with respect to the signal x, where L.x/ is a differentiable convex NLL and

r.x/ D I C .x/ C .x/ (3b)
is a convex regularization term that imposes convex-set constraint on x, x 2 C , and sparsity of an appropriate transformed x through the convex penalty .x/ [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms-theory and practice[END_REF][START_REF] Gu | Projected Nesterov's proximalgradient signal recovery from compressive Poisson measurements[END_REF][START_REF] Gu | Projected Nesterov's proximal-gradient algorithm for sparse signal reconstruction with a convex constraint[END_REF][START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF][START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF]. Here, u > 0 is a scalar tuning constant that quantifies the weight of the regularization term, and I C .x/ , ( 0;

x 2 C C1; otherwise is the indicator function. The penalty .x/ is often selected as the `1-norm of the signal transform-coefficient vector [START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF]:

.x/ D k‰ H xk 1 :

Define the proximal operator for a function r.x/ scaled by > 0 at argument a 2 R p : prox r a D arg min x2R p 1 2 kx ak 2 2 C r.x/:

(5) 0000-0000/00$00.00 © 2017 IEEE

In this paper (see also [START_REF] Gu | Projected Nesterov's proximalgradient signal recovery from compressive Poisson measurements[END_REF][START_REF] Gu | Projected Nesterov's proximal-gradient algorithm for sparse signal reconstruction with a convex constraint[END_REF]), we develop a projected Nesterov's proximal-gradient (PNPG) method whose momentum acceleration accommodates adaptive step-size selection and convex-set constraint on the signal x. Computing the proximal operator with respect to r.x/ in (3b) needs iteration and is therefore inexact [START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF][START_REF] Schmidt | Convergence rates of inexact proximal-gradient methods for convex optimization[END_REF][START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF]. We establish conditions for the O.k 2 / convergence rate of the objective function as well as the convergence of PNPG iterates. These results are the first for an accelerated proximal-gradient (PG) method with step-size adaptation (and, therefore, adjustment to the local curvature of the NLL) that establish convergence of the iterates (Theorem 2) and incorporate inexact proximal operators into objective function convergence rate and convergence of the iterates analyses (Theorems 1 and 2). We modify the original Nesterov's acceleration [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O.1=k 2[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] so that we can establish these results when the step size is adaptive and adjusts to the local curvature of the NLL. (Local-curvature adjustments of the NLL by step-size adaptation have also been used in other algorithms under different contexts in [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF][START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF][START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF]; see also the following and discussion in Section IV-A1.) Our integration of the adaptive step size and convex-set constraint extends the application of the Nesterov-type acceleration to more general measurement models than those used previously, such as the Poisson compressed-sensing scenario described in Section II-A. Furthermore, a convex-set constraint can bring significant improvement to signal reconstructions compared with imposing signal sparsity only, as illustrated in Section V-B. See Section IV-A for further discussion of O.k 2 / acceleration approaches [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF][START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF][START_REF] Auslender | Interior gradient and proximal methods for convex and conic optimization[END_REF][START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF].

Optimization problems (3a) with composite penalty-term structure in (3b) have been considered in [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms-theory and practice[END_REF][START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF][START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF][START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF], which use PG (forward-backward)-type methods with nested inner iterations. The general optimization approach in these references is close to ours. Unlike PNPG, these methods approximate the NLLs whose gradients are not Lipschitz continuous and [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms-theory and practice[END_REF][START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF][START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] do not have fast O.k 2 / convergence-rate guarantees; [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF] observes the benefits of larger step sizes and stepsize adaptation. The nested forward-backward splitting iteration in [START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF] applies fast iterative shrinkage-thresholding algorithm (FISTA) [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] in both the outer and inner loops using duality to formulate the inner iteration; however, it does not employ step-size adaptation or analyze effects of inexact proximalmapping computations. References [START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF][START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF][START_REF] Raguet | A generalized forward-backward splitting[END_REF][START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF] describe splitting schemes to minimize (3a), where [START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF][START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF] are inspired by the parallel proximal algorithm (PPXA) [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Some splitting schemes, e.g., [START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF][START_REF] Pustelnik | Parallel proximal algorithm for image restoration using hybrid regularization[END_REF], apply proximal operations on individual summands L.x/, u .x/, and I C .x/, which is useful if all individual proximal operators are easy to compute. Both [START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF] and generalized forward-backward (GFB) splitting [START_REF] Raguet | A generalized forward-backward splitting[END_REF] require inner iterations to solve prox a for .x/ in (4) in the general case where the sparsifying matrix ‰ is not orthogonal. Reference [START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF] applies the primal-dual approach by Chambolle and Pock [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], which allows solving its Poisson reconstruction problems without approximating the NLL: (3a) is split into L.x/ and r.x/ and also into L.x/ C .x/ and I C .x/, where the second approach (termed CP) does not require nested iterations. The primal-dual splitting (PDS) method in [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF] does not require inner iterations for general L.x/ and sparsifying matrix. GFB and PDS need Lipschitz-continuous gradient of L and the value of the Lipschitz constant is important for tuning their parameters. The convergence rate of both GFB and PDS methods can be upper-bounded by C =k where k is the number of iterations and the constant C is determined by values of the tuning proximal and relaxation constants [START_REF] Liang | Convergence rates with inexact non-expansive operators[END_REF][START_REF] Davis | Convergence rate analysis of primal-dual splitting schemes[END_REF]. In Section V, we show the performances of CP, GFB, and PDS.

Variable-metric methods with problem-specific diagonal scaling matrices have been considered in [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF][START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF][START_REF] Salzo | The variable metric forwardbackward splitting algorithm under mild differentiability assumptions[END_REF]; [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF] applies Barzilai-Borwein (BB) step size and an Armijo line search for the overrelaxation parameter. It accounts for inexact proximal operator and establishes convergence of iterates but does not employ acceleration or provide fast convergence-rate guarantees. [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF] does not require the Lipschitz continuity of the gradient of the NLL in general, except for proving the convergence rate of the objective function. Salzo [START_REF] Salzo | The variable metric forwardbackward splitting algorithm under mild differentiability assumptions[END_REF] analyzes variable-metric algorithms without acceleration (of the type [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF]) and relies on the uniform continuity of rL for the convergence analysis of both objective function and iterates; however, [START_REF] Salzo | The variable metric forwardbackward splitting algorithm under mild differentiability assumptions[END_REF] does not account for inexact proximal operators. In practice, special care is needed in selecting a good scaling matrix, and no clear guidelines are given in [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF][START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF][START_REF] Salzo | The variable metric forwardbackward splitting algorithm under mild differentiability assumptions[END_REF] for this selection. Setting the overrelaxation parameter in [START_REF] Salzo | The variable metric forwardbackward splitting algorithm under mild differentiability assumptions[END_REF] to unity leads to a variablemetric/scaling scheme with an adaptive step size; further, setting the scaling matrix to identity leads to a PG iteration with adaptive step size.

Similar to templates for first-order conic solvers (TFOCS) [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF], PNPG code is easy to maintain: for example, the proximalmapping computation can be easily replaced as a module by the latest state-of-the-art solver. Furthermore, PNPG requires minimal application-independent tuning; indeed, we use the same set of tuning parameters in two different application examples. This is in contrast with the existing splitting methods, which require problem-dependent (NLL-and u-dependent) tuning, with convergence speed sensitive to the choice of tuning constants.

We review the notation: 0, 1, I , denoting the vectors of zeros and ones and identity matrix, respectively; " " is the elementwise version of " "; " T " and " H " are transpose and conjugate transpose, respectively. For a vector a D .a i / N i D1 2 R N , define the projection and soft-thresholding operators:

P C .a/ D arg min x2C kx ak 2 2 (6a) OET .a/ i D sgn.a i / max ja i j ; 0 (6b)
and the elementwise exponential function OEexp ı a i D exp a i . The projection onto R N C and the proximal operator (5) for the `1-norm kxk 1 can be computed in closed form:

P R N C .a/ i D max.a i ; 0/; prox k k 1 a D T .a/: (6c)

A. Preliminary Results

Define the "-subgradient [START_REF] Bertsekas | Convex Optimization Algorithms[END_REF]Sec. 3.3] (" > 0):

@ " r.x/ , ˚g 2 R p j r.z/ r.x/ C .z x/ T g "; 8z 2 R p « (7)
and an inexact proximal operator [START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF]: (11) To simplify the notation, we omit the dependence of G.p/; y x.p/; S.p/ and y p on a and . We will add the subscripts " a; " to these quantities when we wish to emphasize their dependence on a and .

The following proposition extends the result in [14, Sec. We introduce representative NLL functions (Section II), describe the proposed PNPG signal reconstruction algorithm (Section III), establish its convergence properties (Section IV), present numerical examples (Section V), and make concluding remarks (Section VI).

II. Probabilistic Measurement Models

For numerical stability, we normalize the likelihood function so that the corresponding NLL L.x/ is lower-bounded by zero.

A. Poisson Generalized Linear Model

Generalized linear models (GLMs) with Poisson observations are often adopted in astronomic, optical, hyperspectral, and tomographic imaging [START_REF] Prince | Medical Imaging Signals and Systems[END_REF][START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms-theory and practice[END_REF][START_REF] Ollinger | Positron-emission tomography[END_REF] [38, Sec. I-A2], see also [START_REF] Gu | Blind X-ray CT image reconstruction from polychromatic Poisson measurements[END_REF]) and identity:

g. / D b; .x/ D ˆx C b (17)
used for modeling the photon count in optical imaging and radiation activity in emission tomography [3, Ch. 9.2], as well as for astronomical image deconvolution. Here, ˆ2 R N p C and b 2 R N1 C are the known sensing matrix and intercept term, respectively; the intercept b models background radiation and scattering estimated, e.g., by calibration before the measurements y have been collected. The nonnegative set C in (2) satisfies [START_REF] Boyd | Vectors, Matrices, and Least Squares[END_REF], where we have used the fact that the elements of ˆare nonnegative. If b has zero components, C n dom L is not empty and the NLL does not have a Lipschitz-continuous gradient.

B. Linear Model with Gaussian Noise

The linear measurement model with zero-mean additive white Gaussian noise (AWGN) leads to the following scaled NLL:

L.x/ D 1 2 ky ˆxk 2 2 ( 18
)
where y 2 R N is the measurement vector, and constant terms (not functions of x) have been ignored. This NLL belongs to the Gaussian GLM with identity link without intercept: g. / D .

Algorithm 1: PNPG iteration

Input: x .0/ , u, , b, n, m, , Á, and threshold Output: arg min x f .x/ Initialization: x . 1/ 0, i 0, Ä 0, ˇ. [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF] and penalty (3b) with .x/ in ( 4) is an analysis basis pursuit denoising (BPDN) problem with a convex signal constraint.

ˇ.i/ > ˇ.i 1/ then // increase n n n C m ˇ.i/ ˇ.i/ and Ä 0 if f .x .i / / > f .x .i 1/ / then if (21) holds then if x x .i/ ¤ x .i 1/ and f .x .i/ / Ä f .x x .i / / then // function restart  .i 1/ 1, i i 1, and continue if Á > Á MIN and f .x .i/ / > f .x x .i/ /

III. Reconstruction Algorithm

We propose a PNPG approach for minimizing (3a) that combines convex-set projection with Nesterov acceleration [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O.1=k 2[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] and applies adaptive step size to adapt to the local curvature of the NLL and restart to ensure monotonicity of the resulting iteration. The pseudo code in Algorithm 1 summarizes our PNPG method.

Define the quadratic approximation of the NLL L.x/ as

Q ˇ.x j x x/ D L.x x/ C .x x x/ T rL.x x/ C 1 2ˇk x x xk 2 2 (19)
with step-size tuning constant ˇ> 0. Iteration i of the PNPG method proceeds as follows:

B .i/ D ˇ.i 1/ =ˇ. i / (20a)  .i/ D ( 1; i Ä 1 1 C p b C B .i/ . .i 1/ / 2 ; i > 1 (20b) ' .i/ D . .i 1/ 1/= .i / (20c) x x .i/ D P C x .i 1/ C ' .i / .x .i 1/ x .i 2/ / (20d) x .i/ D prox ˇ.i/ ur x x .i/ ˇ.i/ rL.x x .i / / (20e)
where ˇ.i/ > 0 is an adaptive step size chosen to satisfy the majorization condition

L.x .i/ / Ä Q ˇ.i/ .x .i / j x x .i / / (21)
using a simple adaptation scheme that aims at keeping ˇ.i/ as large as possible; see also Section III-B and Algorithm 1. Here,

2; b 2 OE0; 1=4 (22)
in (20b) are momentum tuning constants. We will denote  .i/ as  .i/ ;b when we wish to emphasize its dependence on and b. We declare convergence when

p ı .i/ Ä kx .i / k 2 (23a)
where > 0 is the convergence threshold and ı .i/ is the local variation of signal iterates:

ı .i / , kx .i / x .i 1/ k 2 2 : (23b)
We need B .i / in (20a) to derive the theoretical guarantee for the convergence speed of the PNPG iteration and its sequence convergence. A similar idea for handling the increasing step size in its TFOCS framework is seen in [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF]. However, [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF] does not address this modification in detail or establish convergence properties of the corresponding method.

The acceleration step (20d) extrapolates the two latest iteration points in the direction of their difference x .i 1/

x .i 2/ , followed by the projection onto the convex set C , which has also been proposed in our preliminary work [START_REF] Gu | Projected Nesterov's proximalgradient signal recovery from compressive Poisson measurements[END_REF] and in the variable-metric/scaling method [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF]. For nonnegative C in (2), this projection has closed form; see (6c). If C is an intersection of convex sets with a simple individual projection operator for each, we can apply projections onto convex sets (POCS) [START_REF] Youla | Image restoration by the method of convex projections: Part 1-theory[END_REF].

For .x/ in (4), we compute the proximal mapping (20e) using the dual formulation in [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF] and a simpler version of PNPG, Nesterov's projected-gradient algorithm, because the proximal step in this case reduces to projection onto H in (10c); for the TV penalty, this method is similar to the TV denoising scheme in [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF]. Because of its iterative nature, (20e) is inexact; this inexactness can be modeled as

x .i/ Ñ " .i/ prox ˇ.i/ ur x x .i / ˇ.i/ rL.x x .i / / (24)
where " .i/ quantifies the precision of the PG step in Iteration i .

If we remove the convex-set constraint by setting C D R p , iteration (20a)-(20e) reduces to the Nesterov's proximal-gradient iteration with adaptive step size that imposes signal sparsity only in the analysis form (termed NPG S ); see Section V-B for an illustrative comparison between NPG S and PNPG.

We now extend [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF]Lemma 2.3] to the inexact proximal operation:

Lemma 1: Assume convex and differentiable NLL L.x/ and convex .x/, and consider an inexact PG step (24) with step size ˇ.i/ that satisfies the majorization condition [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF]. Then,

f .x/ f .x .i / / 1 2ˇ. i/ kx .i/ xk 2 2 kx x .i / xk 2 2 ." .i/ / 2 (25)
for all i 1 and any x 2 R p .

Proof: See Appendix A. Lemma 1 is general and algorithm independent, because x

x .i / can be any value in dom L and we have used only the fact that step size ˇ.i/ satisfies the majorization condition [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF], rather than depending on specific details of the step-size selection. We use this result to establish the monotonicity property in Proposition 3 and to derive and analyze our accelerated PG scheme.

A. Restart and Monotonicity

If f .x x .i/ / > f .x .i/ / > f .x .i 1/ / or x x .i/ 2 C n dom L, set  .i 1/ D 1 (restart), (26)
re-evaluate (20b)-(20e), and refer to this action as function restart [START_REF] O'donoghue | Adaptive restart for accelerated gradient schemes[END_REF] or domain restart, respectively; see Algorithm 1. Function and domain restarts ensure that the PNPG iteration is monotonic and x x .i/ remains within dom f as long as the projected initial value is within dom f : f .P C .x .0/ // < C1. In this paper, we employ PNPG iteration with restart, unless specified otherwise (e.g., in Theorems 1 and 2 in Section IV).

Proposition 3 (Monotonicity): The inexact PG step ( 24) is monotonic:

f .x .i/ / Ä f .x x .i / / (27a)
if it is sufficiently accurate such that

" .i / Ä kx x .i/ x .i/ k 2 : (27b)
Hence, the PNPG iteration with restart and inexact PG steps ( 24) is non-increasing:

f .x .i/ / Ä f .x .i 1/ / (28)
if (27b) holds for all i. Proof: (27a) is straightforward by plugging x D x x .i / and (27b) into [START_REF] Raguet | A generalized forward-backward splitting[END_REF].

If there is no restart in Iteration i , the objective function has not increased. If there is a restart, Â .i 1/ D 1, (20d) simplifies to x

x .i / D P C .x .i 1/ / D x .i 1/ , and monotonicity follows due to x

x .i/ D x .i 1/ . To establish the monotonicity in Proposition 3, the step size ˇ.i/ need satisfy only the majorization condition [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF].

B. Adaptive Step Size

Define the step-size adaptation parameter

2 .0; 1/: (29)
We propose the following adaptive scheme for selecting ˇ.i/ : i) if there have been no step-size backtracking events or increase attempts for n consecutive iterations (i n to i 1), start with a larger step size:

x ˇ.i/ D ˇ.i 1/ = (increase attempt); (30a) otherwise start with x ˇ.i/ D ˇ.i 1/ I (30b)
ii) (backtracking search) select

ˇ.i/ D t i x ˇ.i/ (30c)
where 0 Ä t i Ä t MAX is the smallest integer such that (30c) satisfies [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF]; backtracking event corresponds to

t i > 0. iii) if max.ˇ. i/ ; ˇ.i 1/ / < x ˇ.i/ , increase n by a nonnegative integer m: n n C m: (30d)
We select the initial step size x ˇ.1/ using the BB method [START_REF] Barzilai | Two-point step size gradient methods[END_REF]. If there has been an attempt to change the step size in any of the previous n consecutive iterations, we start the backtracking search ii) with the step size from the latest completed iteration. Consequently, ˇ.i/ will be approximately piecewise constant as a function of the iteration index i ; see Fig. 1, which shows the evolutions of ˇ.i/ for measurements following the Poisson generalized linear and Gaussian linear models corresponding to Figs. 4a and 6b in Sections V-A and V-B. To reduce sensitivity to the choice of the tuning constant n, we increase its value by m if there is a failed attempt to increase the step size in Iteration i ; i.e., x ˇ.i/ > ˇ.i 1/ and ˇ.i/ < x ˇ.i/ . The adaptive step-size strategy keeps ˇ.i/ as large as possible subject to [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF], which is important not only because the signal iterate may reach regions of L.x/ with different local Lipschitz constants, but also because of the varying curvature of L.x/ in different updating directions. For example, a (backtrackingonly) PG-type algorithm with non-adaptive step size would fail or converge very slowly if the local Lipschitz constant of rL.x/ decreases as the algorithm iterates, because the step size will not adjust and track this decrease; see also Section V, which demonstrates the benefits of step-size adaptation.

Setting n D C1 corresponds to step-size backtracking only. A step-size adaptation scheme with n D m D 0 initializes the step-size search aggressively, with an increase attempt (30a) in each iteration.

C. Inner-Iteration Warm Start and Convergence Criteria

For .x/ in (4), the inner iteration solves (20e) using the dual problem (10b). Denote by p .i;j / the iterates of the dual variable p in the j th inner iteration step within Iteration i ; this inner iteration solves (20e) using (10b). The initial p .i;0/ is the latest p from Iteration i 1, which is referred to in [START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF] as the warm restart. (The variable metric inexact line-search algorithm (VMILA) [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF] also uses warm restart.)

We consider two convergence criteria. The first tracks local variation of the signal iterates (23b):

kx .i;j / x .i;j 1/ k Ä Á p ı .i 1/ (31a)
where Á is a tuning constant. The second duality-gap-based criterion relies on the result in Proposition 2 to guarantee that . .k/ " .k/ / 2 decreases at a rate of O.k q / within each iteration segment without restart; this guarantee allows us to control the decrease of the convergencerate upper bound in Section IV. Denote by à i the iteration index of the latest restart prior to (and excluding) Iteration i (i 1); set its initial value à 1 D 0. We select the duality-gap based inner-iteration convergence criterion as (see also (10e) and ( 11))

G .i;j / ˇ.i/ u y x .i;j / Ä Á .i à i / q . .i / / 2 (31b)
where Á is a tuning constant and q is the accuracy rate [START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF]. Here, G .i;j / and y x .i;j / are the duality gap G a; .p .i;j / / and y x a; .p .i;j / / in [START_REF] Dupé | Deconvolution under Poisson noise using exact data fidelity and synthesis or analysis sparsity priors[END_REF] and (10e) (respectively) for a D x x .i/ ˇ.i/ rL.x x .i/ / and D ˇ.i/ u. Without restart (i.e., Ã i Á 0) and step-size adaptation, (31b) reduces to the inneriteration convergence criterion in [START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF]Sec. 6.1].

Adjusting Á. We use Á in (31a)-(31b) to trade off accuracy and speed of the inner iterations. If f .x .i/ / > max f .x .i 1/ /; f .x x .i/ / indicating that the monotonicity condition (27b) does not hold, we decrease Á by an order of magnitude (10 times) and re-evaluate (20a)-(20e).

IV. Convergence Analysis

We now bound the convergence rate of the PNPG method without restart.

Theorem 1 (Convergence of the Objective Function): Assume that the NLL L.x/ is convex and differentiable, .x/ is convex, the closed convex set C satisfies

C Â dom L ( 32
)
(implying no need for domain restart), and the momentum tuning constants are within the range [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. Consider the PNPG iteration without restart with the inexact PG step [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] in place of (20e).

The convergence rate of the PNPG iteration is bounded as follows: for k 1,

.k/ Ä kx .0/ x ? k 2 2 C E .k/ 2ˇ. k/ . .k/ / 2 (33a) Ä 2 kx .0/ x ? k 2 2 C E .k/ 2 p ˇ.1/ C P k iD1 p ˇ.i/ 2 (33b)
where x ? is a minimum point of f .x/ and

.k/ , f .x .k/ / f .x ? / (34a)

E .k/ , k X i D1 . .i/ " .i/ / 2 (34b)
are the centered objective function and the cumulative error term, which accounts for the inexact PG steps, respectively.

Proof: See Appendix A for the proof of (33a); then, to obtain (33b), use

 .k/ q ˇ.k/ 1 q ˇ.k/ C  .k 1/ q ˇ.k 1/ (35a) 1 k X i D2 q ˇ.i/ C  .1/ q ˇ.1/ (35b)
for all k > 1, where (35a) follows from the definitions of B .k/ and  .k/ in (20a) and (20b), and (35b) follows by repeated application of the inequality (35a) with k replaced by k 1; k 2; : : : ; 2.

Theorem 1 shows that better initialization, smaller proximalmapping approximation error, and larger step sizes .ˇ. i / / k iD1 help lower the convergence-rate upper bounds in [START_REF] Salzo | The variable metric forwardbackward splitting algorithm under mild differentiability assumptions[END_REF]. This result motivates our step-size adaptation with the goal of maintaining large .ˇ. i / / k i D1 ; see Section III-B. To derive this theorem, we have used only the fact that the step size ˇ.i/ satisfies the majorization condition [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF], rather than taking advantage of specific details of the step-size selection.

To minimize the upper bound in (33a), we can select  .i/ to satisfy (A16b) with equality, which corresponds to  .i/ 2;1=4 in (20b), on the boundary of the feasible region in [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. By (35a), p ˇ.k/  .k/ and the denominator of the bound in (33a) are strictly increasing sequences. The upper bound in (33b) is not a function of b and is minimized with respect to for D 2, given the fixed step sizes .ˇ. i/ / C1 i D0 . Corollary 1: Under the assumptions of Theorem 1, the convergence of PNPG iteration x .k/ without restart is bounded as follows:

.k/ Ä 2 kx .0/ x ? k 2 2 C E .k/ 2.k C 1/ 2 ˇmin (36a) for k 1, provided that ˇmin , C1 min kD1 ˇ.k/ > 0: (36b)
Proof: Use (33b) and the fact that

p ˇ.1/ C P k iD1 p ˇ.i/ .k C 1/ p ˇmin .
The assumption (36b) is implied by, and weaker than, the Lipschitz continuity of rL.x/; indeed, ˇmin > =L if rL.x/ has a Lipschitz constant L; see also [START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF].

According to Corollary 1, the PNPG iteration attains the O.k2 / convergence rate as long as the step size ˇ.i/ is bounded away from zero (see (36b)) and the cumulative error term (34b) converges:

E .C1/ , lim k!C1 E .k/ < C1 (37)
which requires that . .k/ " .k/ / 2 decreases at a rate of O.k q / with q > 1. This condition, also key for establishing convergence of iterates in Theorem 2, motivates us to use decreasing convergence criteria (31a)-(31b) for the inner proximal-mapping iterations, where (31b) guarantees [START_REF] Mccullagh | Generalized Linear Models[END_REF] upon choosing an appropriate q.

We now contrast our result in Theorem 1 with existing work on accommodating inexact proximal mappings in PG schemes. By recursively generating a function sequence that approximates the objective function, [START_REF] Villa | Accelerated and inexact forward-backward algorithms[END_REF] gives an asymptotic analysis of the effect of " .i/ on the convergence rate of accelerated PG methods with inexact proximal mapping. However, no explicit upper bound is provided for .k/ . Schmidt et al. [START_REF] Schmidt | Convergence rates of inexact proximal-gradient methods for convex optimization[END_REF] provide convergence-rate analysis and an upper bound on .k/ , but their analysis does not apply here because it relies on fixed step-size assumption, uses different form of acceleration [13, Prop. 2], and has no convex-set constraint. Bonettini et al. [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF] analyze the inexactness of proximal mapping but for proximal variablemetric/scaling methods with O.k 1 / convergence rate for the objective function.

We now establish convergence of the PNPG iterates. Theorem 2 (Convergence of Iterates): Assume that 1) the conditions of Theorem 1 hold;

2) E .C1/ exists: (37) holds;

3) the momentum tuning constants . ; b/ satisfy

> 2; b 2 OE0; 1= 2 I (38)
4) the step-size sequence .ˇ. i/ / C1 i D1 is bounded within the range OEˇm in ; ˇmax , for ˇmin > 0.

Consider the PNPG iteration without restart with the inexact PG step [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF] in place of (20e). Then, the sequence of PNPG iterates x .i/ converges to a minimizer of f .x/.

Proof: See Appendix B. Observe that Assumption 3 requires a narrower range of . ; b/ than (22): indeed (38) is a strict subset of [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. The intuition is to leave a sufficient gap between the two sides of (A16a) so that their difference becomes a quantity that is roughly proportional to the growth of  .i / , which is important for proving the convergence of signal iterates [START_REF] Chambolle | On the convergence of the iterates of the 'fast iterative shrinkage/thresholding algorithm[END_REF]. Although the momentum term (20b) with D 2 is optimal in terms of minimizing the upper bound on the convergence rate (see Theorem 1), it appears difficult or impossible to prove convergence of the signal iterates

x .i / for this choice of because, in this case, the gap between the two sides of (A16a) is upper-bounded by a constant.

Iterate convergence results in [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF][START_REF] Chambolle | On the convergence of the iterates of the 'fast iterative shrinkage/thresholding algorithm[END_REF][START_REF] Aujol | Stability of over-relaxations for the forward-backward algorithm, application to FISTA[END_REF] apply to momentumaccelerated methods that require non-increasing step-size sequences and do not adjust to the local curvature of the NLL. Aujol and Dossal [START_REF] Aujol | Stability of over-relaxations for the forward-backward algorithm, application to FISTA[END_REF] analyze both the convergence of the objective function and the iterates with inexact proximal operator for B .1/ D 1 and n D 1, i.e., with decreasing step size only, and for a different (less aggressive) Â .i/ than ours in (20b). Bonettini et al. use the ideas from [START_REF] Chambolle | On the convergence of the iterates of the 'fast iterative shrinkage/thresholding algorithm[END_REF] to establish convergence of iterates for their variable-metric/scaling approach in [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF], but this analysis does not account for inexact proximal steps.

A. O.k 2 / Convergence Acceleration Approaches

A few variants that accelerate the PG method achieve the O.k 2 / convergence rate [18, Sec. 5.2]. One competitor proposed by Auslender and Teboulle in [START_REF] Auslender | Interior gradient and proximal methods for convex and conic optimization[END_REF]Sec. 5] and restated in [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF] where it was referred to as AT, replaces (20d)-(20e) with

x x .i/ D 1 1  .i/ Á x .i 1/ C 1  .i / z x .i 1/ (39a) z x .i/ D prox  .i / ˇ.i/ ur z x .i 1/  .i/ ˇ.i/ rL.x x .i / / (39b) x .i/ D 1 1  .i/ Á x .i 1/ C 1  .i / z x .i/ (39c)
where  .i/ D  .i / 2;1=4 in (20b). Here, ˇ.i/ in the TFOCS implementation [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF] is selected using the aggressive search with n D m D 0.

All intermediate signals in (39a)-(39c) belong to C and do not require projections onto C . However, as  .i/ increases with i , step (39b) becomes unstable, especially when an iterative solver is needed for its proximal operation. To stabilize its convergence, AT relies on periodic restart by resetting  .i/ using (26) [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF]. However, the period of restart is a tuning parameter that is not easy to select. For a linear Gaussian model, this period varies according to the condition number of the sensing matrix ˆ [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF], which is generally unavailable and not easy to compute for largescale problems. For other models, there are no guidelines how to select the restart period.

In Section V, we show that AT converges slowly compared with PNPG, which justifies the use of projection onto C in (20d) and (20d)-(20e) instead of (39a)-(39c). PNPG usually runs uninterrupted (without restart) over long stretches and benefits from Nesterov's acceleration within these stretches, which may explain its better convergence properties compared with AT. PNPG may also be less sensitive than AT to proximal-step inaccuracies; we have established convergence-rate bounds for PNPG that account for inexact proximal steps (see [START_REF] Salzo | The variable metric forwardbackward splitting algorithm under mild differentiability assumptions[END_REF] and (36a)), whereas AT does not yet have such bounds, to the best of our knowledge.

1) Relationship with FISTA:

The PNPG method is a generalized FISTA [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] that accommodates convex constraints, more general NLLs, 2 and (increasing) adaptive step size. In contrast with PNPG, FISTA has a non-increasing step size ˇ.i/ , which allows for setting B .i / D 1 in (20b) for all i (see Appendix A-II); upon setting . ; b/ D .2; 1=4/, this choice yields the standard FISTA (and Nesterov's [START_REF] Nesterov | A method of solving a convex programming problem with convergence rate O.1=k 2[END_REF]) update. Convergence of signal iterates has not been established for FISTA with . ; b/ D .2; 1=4/ The method in [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF] is a variable-metric/scaling version of FISTA with projection of the extrapolation step to account for the convex constraints. [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF] analyzes a version of FISTA where the (decreasing) step size is adjusted using a condition in [START_REF] Tseng | A modified forward-backward splitting method for maximal monotone mappings[END_REF] different from the majorization condition ( 21), and establishes objective-function convergence under the assumption that the step size is lower bounded. As FISTA, [START_REF] Bonettini | A variable metric forward-backward method with extrapolation[END_REF] and [START_REF] Bello Cruz | On the convergence of the forward-backward splitting method with linesearches[END_REF] do not adapt the step size and hence do not adjust to the local curvature of the NLL.

V. Numerical Examples

We now evaluate our proposed algorithm by means of numerical simulations. We adopt the nonnegative C in (2) and the `1-norm sparsifying penalty in [START_REF] Harmany | This is SPIRAL-TAP: Sparse Poisson intensity reconstruction algorithms-theory and practice[END_REF]. The PNPG iterations with the local-variation and duality-gap inner convergence criteria (31a) and (31b) are labeled PNPG and PNPG d , respectively.

All iterative methods that we compare use the convergence criterion (23a) with

D 10 9 (40)
and have the maximum number of iterations I max D 10 4 . In the presented examples, PNPG uses momentum tuning constants . ; b/ D .2; 1=4/ and adaptive step-size parameters n D 4 (unless specified otherwise), m D n, D 0:8, inner-iteration convergence constants Á D 10 2 and .Á; q/ D .1; 1:0001/ for PNPG and PNPG d (respectively), and maximum number of inner iterations J max D 1000. Here, PNPG d uses q D 1:0001 with goal to guarantee [START_REF] Mccullagh | Generalized Linear Models[END_REF]. We apply the AT method (39) implemented in the TFOCS package [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF] with adaptive step size and a periodic restart every 200 iterations (tuned for its best performance) and our proximal mapping. Our inner convergence criteria (31b) cannot be implemented in the TFOCS package (i.e., it require editing its code). Hence, we select the proximal mapping that has a relative-error inner convergence criterion kp .i;j / p .i;j 1/ k 2 Ä 0 kp .i;j / k 2 ;

(41a) where p .i;j / is the dual variable employed by the inner iterations. This relative-error inner convergence criterion is easy to incorporate into the TFOCS software package [START_REF] Becker | Templates for convex cone problems with applications to sparse signal recovery[END_REF] and is already used by the SPIRAL package; see [START_REF] Harmany | The sparse Poisson intensity reconstruction algorithms (SPIRAL) toolbox[END_REF]. Here, we select

0 D 10 6 (41b)
for both AT and SPIRAL and set their maximum number of inner iterations to 100.

We apply the CP approach based on [23, Sec. 7.5]: [START_REF] Chambolle | On the convergence of the iterates of the 'fast iterative shrinkage/thresholding algorithm[END_REF]. In this paper, we adopt the version of CP in [START_REF] Chambolle | On the convergence of the iterates of the 'fast iterative shrinkage/thresholding algorithm[END_REF].

z prox 1 F z C 1 ˆx x (42a) p P uH .p C 2 ‰ H x x/ (42b) x x 2P C x ˆT z C Re.‰p/ Á x (42c) x .x x C x/=2 (
All the numerical examples were performed on a Linux workstation with an Intel Xeon CPU E31245 (3.30 GHz) and 8 GB memory. The operating system is Ubuntu 14.04 LTS (64bit). The Matlab implementation of the proposed algorithms and numerical examples is available [START_REF] Gu | Projected Nesterov's proximal-gradient algorithm source code[END_REF].

A. PET Image Reconstruction from Poisson Measurements

In this example, we adopt the Poisson GLM (16a) with identity link in [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF]. Consider PET reconstruction of the 128 128 concentration map x in Fig. 2a, which represents simulated radiotracer activity in the human chest. Assume that the corresponding 128 128 attenuation map Ä is known, which is needed to model the attenuation of the gamma rays [START_REF] Ollinger | Positron-emission tomography[END_REF] and compute the sensing matrix ˆin this application. We collect the photons from 90 equally spaced directions over 180°, with 128 radial samples in each direction. Here, we adopt the parallel strip-integral matrix [START_REF] Fessler | Image Reconstruction[END_REF]Ch. 25.2] and use its implementation in the Image Reconstruction Toolbox (IRT) [START_REF] Fessler | Image reconstruction toolbox[END_REF] with sensing matrix ˆD w diag exp ı . Ä C c/ (43) where c is a known vector generated using a zero-mean independent, identically distributed (i.i.d.) Gaussian sequence with variance 0:3 to model the detector efficiency variation; w > 0 is a scaling constant controlling the expected total number of detected photons due to electron-positron annihilation; and 1 T E.y b/ D 1 T ˆx, which is a signal-to-noise ratio (SNR) measure. Assume that the background radiation, scattering effect, and accidental coincidence combined lead to a known (generally nonzero) intercept term b in the Poisson GLM [START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF]. The elements of the intercept term have been set to a constant equal to 10 % of the sample mean of ˆxtrue : b D .1 T ˆxtrue /=.10N /1.

The above model, choices of parameters in the PET system setup, and concentration map have been adopted from IRT [49, emission/em_test_setup.m].

Here, we consider the DWT and isotropic TV sparsifying transforms. We use the 2-D Haar DWT with 6 decomposition levels and a full circular mask [START_REF] Dogandžić | Mask iterative hard thresholding algorithms for sparse image reconstruction of objects with known contour[END_REF] to construct a sparsifying dictionary matrix ‰ 2 R 12 449 14 056 with orthonormal rows, i.e., ‰‰ T D I , which allows efficient inner iteration. We compare the filtered backprojection (FBP) [START_REF] Ollinger | Positron-emission tomography[END_REF] and PG methods that aim at minimizing (3) with nonnegative C in (2) and DWT and TV sparsifying transforms.

We implemented SPIRAL with TV penalty using the centered NLL term (16a), which improves the numerical stability compared with the original code in [START_REF] Harmany | The sparse Poisson intensity reconstruction algorithms (SPIRAL) toolbox[END_REF]. We do not compare with SPIRAL that uses DWT penalty because its inner iteration for the proximal step requires a square orthogonal ‰ (see [4, Sec. II-B]), which is not the case here. We also compare with VMILA [START_REF] Bonettini | Variable metric inexact line-search-based methods for nonsmooth optimization[END_REF][START_REF] Bonettini | Variable metric inexact line-search algorithm (VMILA)[END_REF] with both DWT and TV penalties and its default tuning constants, which yield good performance; hence VMILA is insensitive to tuning.

In this example, we adopt the following form of the regularization constant u:

u D 10 a ; (44)
vary a in the range OE 6; 3 with a grid size of 0.5, and search for the reconstructions with the best average relative square error (RSE) performance; here, RSE D ky x x true k 2 2 =kx true k 2 2 , where x true and y

x are the true and reconstructed signals, respectively. All iterative methods were initialized by FBP reconstructions implemented by IRT [START_REF] Fessler | Image reconstruction toolbox[END_REF].

Figs. 2b-2d show reconstructions and corresponding RSEs for one random realization of the noise and detector variation c, with the expected total annihilation photon count (SNR) equal to 10 8 ; the optimal a is 0.5. All sparse reconstruction methods (PNPG, AT, CP, SPIRAL, and VMILA) perform similarly as long as they employ the same penalty: the TV sparsity penalty is superior to its DWT counterpart; see [9, Fig. 6] which shows average RSEs of different methods as functions of 1 T ˆxtrue .

Figs. 3 and4 show the normalized centered objectives .i/ =f .x ? / as functions of the number of iterations and CPU time for the DWT and TV signal sparsity regularizations and two random realizations of the noise and detector variation with different total expected photon counts. The legends in Figs. 3b and 4b apply to Figs. 3a and 4a as well. Fig. 3 examines the convergence of PNPG as a function of the momentum tuning constants . ; b/ in (22), using 2 f2; 5; 15g and b 2 f0; 1=4g. For small Ä 5, there is no significant difference between different selections and no choice is uniformly the best, consistent with [START_REF] Chambolle | On the convergence of the iterates of the 'fast iterative shrinkage/thresholding algorithm[END_REF] which considers only b D 0 and nonadaptive step size. As we increase further ( D 15), we observe slower convergence. In the remainder of this section, we use . ; b/ D .2; 1=4/.

To illustrate the benefits of step-size adaptation, we present in Fig. 4 the performance of PNPG (n D 1), which does not adapt to the local curvature of the NLL and has monotonically non-increasing step sizes ˇ.i/ , similar to FISTA. PNPG (n D 4) outperforms PNPG (n D 1) because it uses step-size adaptation; see also Fig. 1a, which corresponds to Fig. 4a and shows that the step size of PNPG (n D 4) is consistently larger than that of PNPG (n D 1). Initializing PNPG iterations by a vector close to 0 (rather than with FBP) will lead to an even larger difference in convergence speed between PNPG (n D 1) and PNPG (n D 4). The advantage of PNPG (n D 4) over the aggressive PNPG (n D 0) scheme is due to the patient nature of its step-size adaptation, which leads to a better local majorization function of the NLL and reduces time spent backtracking. Indeed, if we do not account for the time spent on each iteration and only compare the objectives as functions of the iteration index, then PNPG (n D 4) and PNPG (n D 0) perform similarly; see [START_REF] Gu | Projected Nesterov's proximalgradient signal recovery from compressive Poisson measurements[END_REF]Fig. 4 further outperformed by PNPG (n D 4). Fig. 4b shows that SPIRAL, which does not employ PG step acceleration, takes at least three times longer than PNPG (n D 4) to reach the same objective function.

In Fig. 4b, AT and SPIRAL reach the performance floor due to their fixed inner convergence criterion in (41a); we observe a performance floor for AT in Fig. 4a as well. Reducing 0 in (41b) will lower this floor, at the cost of slowing down the two algorithms. This result justifies our convex-set projection in (20d) for the Nesterov's acceleration step, shows the superiority of (20d) over AT's acceleration in (39a) and (39c), and is consistent with the results in Section V-B. PNPG d (n D 4) employs the duality-gap-based inner convergence criterion (31b) with q D 1:0001. Since the goal of using (31b) is to guarantee [START_REF] Mccullagh | Generalized Linear Models[END_REF], this inner criterion is more stringent and leads to slower overall performance of PNPG d (n D 4) compared to PNPG (n D 4). Indeed, if we do not account for the time spent on each iteration and only compare the objectives as functions of the iteration index, then PNPG d (n D 4) and PNPG (n D 4) perform similarly with the former slightly better.

The CP method uses the following tuning constants carefully selected for this particular problem: . 1 ; 2 ; / D .10 6 ; 1; 103 / and . 1 ; 2 ; / D .104 ; 1; 10 2 / for the DWT and TV penalties, respectively. CP is sensitive to tuning and a different selection of . 1 ; 2 ; / can significantly slow down its convergence. Initially, CP converges quickly and then slows down as it approaches the optimum.

Considering its O.k 1 / theoretical convergence rate, VMILA performs quite well, thanks to its use of the variablemetric/scaling approach.

B. Skyline Signal Reconstruction from Linear Measurements

We adopt the DWT sparsifying transform and linear measurement model with Gaussian noise in Section II-B where each column of the sensing matrix ˆare i.i.d. and drawn from the uniform distribution on unit sphere. Due to the widespread use of this measurement model, we can compare a wider range of methods than in the Poisson PET example in Section V-A.

We have designed a "skyline" signal of length p D 1024 by overlapping magnified and shifted triangle, rectangle, sinusoid, and parabola functions; see Fig. 5a. We generate the noiseless measurements using y D ˆxtrue . The DWT matrix ‰ is constructed using the Daubechies-4 wavelet with three decomposition levels whose approximation by the 5 % largestmagnitude wavelet coefficients achieves RSE D 98 %. We compare:

AT, PNPG, and PNPG d ; CP with the parameters 2 D 1 [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] and D 1 with 1 tuned separately for best performance in each experiment; 3 linearly constrained gradient projection method [START_REF] Harmany | Gradient projection for linearly constrained convex optimization in sparse signal recovery[END_REF], part of the SPIRAL toolbox [START_REF] Harmany | The sparse Poisson intensity reconstruction algorithms (SPIRAL) toolbox[END_REF] and labeled SPIRAL herein; the GFB method [START_REF] Raguet | A generalized forward-backward splitting[END_REF] (see ( 4)):

z 1 z 1 C prox .ru=w/ .2x z 1 rrL.x// x (45a) z 2 z 2 C OEP C .2x z 2 rrL.x// x (45b) x wz 1 C .1 w/z 2 (45c)
with r D 1:9=kˆk 2 2 , D 1, and w D 0:5 tuned for best performance; and the PDS method [START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF]:

N z P OE u;u p .z C ‰ T x/ (46a) N x P C x rL.x/ ‰.2N z z/ (46b) z z C r.N z z/ (46c) x x C r. N x x/ (46d)
where we choose D 1=. C kˆk 2 2 =2/ and r D 2 0:5kˆk 2 2 . 1 / 1 with tuned for best performance, 4 all of which aim to solve the generalized analysis BPDN problem with a convex signal constraint. Here, p 0 D p, ‰ 2 R p p is an orthogonal matrix (‰‰ T D ‰ T ‰ D I ), and prox a D ‰ T .‰ T a/ has a closed-form solution (see (6c)), which simplifies the implementation of the GFB method ((45a), in particular); see the discussion in Section I. The other tuning options for SPIRAL, and AT are kept to their default values, unless specified otherwise.

where a is an integer selected from the interval OE 9; 1 and U is an upper bound on u of interest. Indeed, the minimum point

x ? reduces to 0 if u U [38, Sec. II-D].
As before, PNPG (n D 4) and PNPG (n D 0) converge at similar rates as functions of the number of iterations. However, due to the excessive attempts to increase the step size at every iteration, PNPG (n D 0) spends more time backtracking and converges at a slower rate as a function of CPU time compared with PNPG (n D 4); see also Fig. 6 presents the normalized centered objectives .i/ =f .x ? / as functions of CPU time for a random realization of the sensing matrix ˆwith normalized numbers of measurements N=p D 0:34 and several different regularization constants a. (For the GFB method, we compute the normalized centered objectives using P C .x .i / / instead of x .i/ in (34a) because its x .i/ may be outside C .) The legend in Fig. 6c applies also to Figs. 6a and6b. To achieve good performance, CP and PDS need to be manually tuned for each a. CP and PDS have optimal 1 D 2 equal to 0:01; 0:1; 1 and 0:0026; 0:026; 2:6 for a D 5; 4; 3, respectively.

PNPG and PNPG d have the steepest descent rate, followed by PNPG d . AT and SPIRAL reach the performance floor around the relative precision of 10 6 due to their fixed inner convergence criterion in (41a). The GFB and primal-dual methods, PDS and CP, are sensitive to the selection of the tuning constants. After a careful selection of the tuning constants, CP performs exceptionally well in Figs. 6a and6b. The performance of GFB is affected significantly by the value of the regularization parameter a.

VI. Conclusion

We developed a fast algorithm for reconstructing sparse signals that belong to a closed convex set by employing a projected proximal-gradient scheme with Nesterov's acceleration, restart, and adaptive step size. We applied the PNPG method to construct one of the first Nesterov-accelerated proximal-gradient reconstruction algorithm for Poisson compressed sensing. We presented integrated derivation of the proposed algorithm and convergence-rate upper bound that accounts for inexactness of the proximal operator and also proved convergence of iterates.

Our PNPG approach is computationally efficient compared with other state-of-the-art methods. and b .i/ D w 2 .a .i/ / 2 a .i/ ; see also (A5a). Upon defining

 .i/ , w 2 a .i/ (A11a)
we have

w 2 c .i/ D . .i / / 2 I w 2 b .i/ D . .i / / 2  .i / : (A11b)
Plug (A11) into (A10) and reorganize to obtain the following form of momentum acceleration:

y x .i/ D x .i 1/ C ' .i/ .x .i 1/ x .i 2/ /: (A12) Although x x .i/ D y
x .i / satisfies (A6a), it is not guaranteed to be within dom L; consequently, the proximal-mapping step for this selection may not be computable.

b) Selecting x x .i/ 2 C that satisfies (A6a): We now seek x x .i/ within C that satisfies the inequality (A6a). Since x .i 1/ and x ? are in C , z .i / 2 C by the convexity of C ; see (A5c). According to Lemma 2, projecting (A12) onto C preserves or reduces the distance between points. Therefore, x x .i/ D P C .y x .i/ / (A13) belongs to C and satisfies the condition (A6a):

c .i 1/ t .i 1/ D c .i/ ky x .i/ z .i / k 2

Finally, let k ! C1 and use [START_REF] Mccullagh | Generalized Linear Models[END_REF] and (B2) to conclude (B1). where we obtain the inequality (B7a) by combining the terms on the right-hand size and using (35a) and (B7b) holds because p ˇ.k/ Â .k/ is an increasing sequence (see Section IV). Now,

… j Ä C1 X kDj k Y `Dj
ˇ.` 2/ . .` 2/ / 2 ˇ.` 1/  .` 1/  .`/ D C1 X kDj ˇ.j 2/ . .j 2/ / 2  .j 1/ ˇ.k 1/ . .k 1/ / 2  .k/ (B8a) Ä ˇ.j 2/ . .j 2/ / 2  .j 1/ p ˇ.j 2/  .j 2/ p ˇ.j 1/  .j 1/ D p B .j 1/  .j 2/ (B8b)

where (B8a) follows by using (20c), (A16a) with i D ` 1, and fraction-term cancellation; (B8b) is obtained by substituting (B7b) into (B8a) and canceling summation terms. (B8b) implies (B6) by using (35a) with k D j 1. Define

.i/ , kx .i / x ? k 2 2 ; ƒ .i/ , .i/ .i 1/ : (B9) Since f .x .i/ / converges to f .x ? / D min x f .x/ as the iteration index i grows and x ? is a minimizer, it is sufficient to prove the convergence of .i/ ; see [START_REF] Chambolle | On the convergence of the iterates of the 'fast iterative shrinkage/thresholding algorithm[END_REF]Th. 4 where (B13b) is due to 1 < 1= (see [START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF]) and the following: where we have used (20c), the fact that p ˇ.i/ Â .i/ is an increasing sequence, ˇ.i/ =ˇ. i 1/ 1= (see Section III-B), and [START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF].

' .i / < Â .i
According to (35b) and the fact that the sequence .ˇ. i / / is bounded (by Assumption 4), there exists an integer J such that  .j 1/ 2; ' .j / 1  .j / > 0 (B15) for all j J , where the second inequality follows from the first and the definition of ' .j / ; see (20c). Then .i/ , max.0; ƒ .i / / Ä ' .i / Ä

. ' .`/ (B16b) for i J , where the inequality in (B16a) follows by combining the inequalities (B13b) and .i 1/ ƒ .i 1/ , and (B16b) follows by recursively applying inequality (B16a) with i replaced by i 1; i 2; : : : ; J . Now, sum the inequalities (B16b) over i D J; J C 1; : : : ; C1 and exchange the order of summation over i and j on the right-hand side (see also (B6)):

C1 X i DJ .i/ Ä C1 X j DJ … j Ä 2ı .j 1/ C
." .j / / 2 ' .j / C … J .J 1/ : (B17) For j J 3, .2 .j 1/ 1/ … j . .j 1/ 1/ C 1 > 0 (B18a) 2 . .j 1/ 1/ … j . .j 1/ 2/ C 1 > 0 (B18b)

where the first and second inequalities in (B18) follow by applying Lemma 4 and (B15), respectively; consequently, C1 X j DJ … j ı .j 1/ Ä C1 X j DJ .2Â .j / 1/ı .j / < C1 (B19a) C1 X j DJ … j ." .j / / 2 ' .j / Ä 2 C1 X j DJ

." .j / / 2 Â .j 1/ 1 ' .j / (B19b)

D 2 C1 X j DJ
." .j / / 2 Â .j / (B19c)

Ä 2 C1 X j DJ
. .j / " .j / / 2 (B19d)

where (B19a) follows from (B18a) and Lemma 3 (for the second inequality) and (B19b) follows by using (B18b); (B19c) and (B19d) are due to (20c) and (B15), respectively. Combine (B19a) and (B19d) with (B17) to conclude that C1 X i D1

.i/ < C1: (B20)

The remainder of the proof uses the technique employed by Chambolle and Dossal to conclude the proof of [42, Th. 4.1, p. 978], which we repeat for completeness. Define X .i/ , .i/ P i j D1

.j / , which is lower bounded because .i/ and P i j D1

.j / are lower and upper bounded, respectively; see (B9) and (B20). Furthermore, .X .i/ / is an non-increasing sequence:

X .i C1/ D .i C1/ .i C1/ i X j D1
.j / Ä X .i/ ; (B21)

where we used the fact that .i C1/ ƒ .i C1/ D .i C1/ .i/ . Hence, .X .i/ / converges as i ! C1. Since P i j D1

.j / converges, . .i/ / also converges.

Fig. 1 :

 1 Fig.1:Step sizes ˇ.i/ as functions of the number of iterations for Poisson and Gaussian linear models.

Fig. 2 :

 2 Fig. 2: (a) True emission image and (b)-(d) the reconstructions of the emission concentration map.

Fig. 3 :

 3 Fig. 3: Normalized centered objectives as functions of the number of iterations for (a) DWT and (b) TV regularizations.

  ]. Although PNPG (n D 0) and AT have the same step-size selection strategy and O.k 2 / convergence-rate guarantees, PNPG (n D 0) converges faster; both schemes are

Fig. 4 :

 4 Fig. 4: Normalized centered objectives as functions of the CPU time for (a) DWT and (b) TV regularizations.

Fig. 5 : 3 Fig. 6 :

 536 Fig. 5: Nonnegative skyline signal and its PNPG and NPG S reconstructions for N=p D 0:34.

  Fig. 1b which corresponds to Fig. 6b and shows the step sizes as functions of the number of iterations for a D 4 and N=p D 0:34. Hence, we present only the performances of PNPG with n D 4 in this section. Fig. 5 shows the advantage brought by the convex-set nonnegativity signal constraints (2). Figs. 5b and 5c present the PNPG .a D 5/ and NPG S .a D 4/ reconstructions from one realization of the linear measurements with N=p D 0:34 and a tuned for the best RSE performance. Recall that NPG S imposes signal sparsity only. Here, imposing signal nonnegativity significantly improves the overall reconstruction and does not simply rectify the signal values close to zero.

Lemma 4 :j 3 , 1 /

 431 For  .k 1/  .k/ Ä p ˇ.k 1/  .k 1/ p ˇ.k/  .k/

  Proposition 1: x Ñ " prox ur a implies kx prox ur ak 2 Ä ".

	2u if .a x/=u 2 @ " 2	r.x/.
	Proof: By Definition 1, the following holds for any z:
	ur.z/ ur.x/ C .z x/ T .a x/ 0:5" 2	(9a)
	which is equivalent to
	0:5kz ak 2 2 C ur.z/ 0:5kx ak 2 2 C ur.x/ C 0:5kz xk 2 2 0:5" 2 : (9b)
	Since z D prox ur a minimizes the left-hand side of (9b), substituting it into (9b) completes the proof.
	Now, we adapt the results in [24, Sec. IV-A] and [23,
	Sec. 5.2.4] to complex ‰ using the fact that, for complex y
	and p, kyk 1 D max kpk 1 Ä1 Re.p H y/. The proximal operator (5) with .x/ in (4) can be rewritten as
	x.y p/ prox r a D y p 2 C p 0 solves the dual problem [23, 24]: where y y p D arg min p2H 1 2 kS.p/k 2 2 1 2 kS.p/ y x.p/k 2 2	(10a) (10b)
	and
		H , ˚w 2 C p 0 S.p/ , a Re.‰p/ j kwk 1 Ä 1 «	(10c) (10d)
		y x.p/ , P C S.p/ 2 R p :			(10e)
	When p 2 H , the objective function in (10b) is differentiable with respect to the real and imaginary parts of p. When ‰
	p must be real-valued and hence is real-valued, the optimal y (10b) reduces to optimization with respect to p over the unit
	hypercube.
	The duality gap for the optimization problem (10b) is
	G.p/ D	˚	y x.p/	y x T .p/ Re.‰p/	«	C I H .p/:
							with
							"-precision, denoted
							x Ñ " prox ur a	(8)

Definition 1:

We say that x approximates prox ur .a/

  W R p 7 ! R C is a convex function of the signal x.Here, the relationship between the linear predictor ˆx and the expected value .x/ of the measurements y is summarized by the link function g. / W R N 7 ! R N [37]:

		and are used to model event
	counts, e.g., numbers of particles hitting a detector. Assume
	that the measurements y D .y n / N nD1 2 N N 0 are independent Poisson-distributed 1 with means OE .x/ n .
	Upon normalization, we obtain the generalized Kullback-
	Leibler divergence form of the NLL [36]
	L.x/ D 1 T OE .x/ y C	X n;y n ¤0	y n ln	y n OE .x/ n	: (16a)
	(16b) Two typical link functions in the Poisson GLM are log The NLL L.x/ E.y/ D .x/ D g 1 .ˆx/: Note that cl.dom L/ D fx 2 R p j .x/ 0g.
	(described in

  Another version of CP can be obtained by associating the regularization parameter u with ‰ instead of the `1-norm function, which leads to replacing uH with H and ‰; ‰ H with u‰; u‰ H , respectively, in

42d) obtained by splitting the objective function (3a) into the sum of F .ˆx/ C uk‰ H xk 1 and I C .x/, where the first summand is a convex lower semicontinuous function of OEˆH ‰ H x and F .ˆx/ D L.x/.

In our examples in Sections V-A and V-B, F .y/ and its convex conjugate F .z/ have analytical proximal operators

[START_REF] Anthoine | Some proximal methods for Poisson intensity CBCT and PET[END_REF] Sec. 7.4]

; hence, CP does not require an inner iteration. In the original CP algorithm in

[START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]

, Chambolle and Pock select 1 D 2 D . However, when the difference between kˆk 2 and k‰k 2 becomes large, it is hard to find tuning constants of the form . 1 ; 2 ; / D . ; ; / that ensure fast convergence of the CP algorithm, which is why we do not impose 1 D 2 D here.

  .1]. Proof of Theorem 2: Use (A3a) and .i / 0 to obtain C .' .i / / 2 ı .i 1/ C 2' .i/ .x .i 1/ x ? / T .x .i 1/ x .i 2/ / (B11a)Ä .i 1/ C .' .i/ / 2 ı .i 1/ C ' .i / .ƒ .i 1/ C ı .i 1/ / (B11b)where (B11a) and (B11b) follow by using the nonexpansiveness of the projection operator (see also (A12)) and the identity 2.a b/ T .a c/ D ka bk 2 2 C ka ck 2

			0	.i /	kx x .i/ x ?	k 2 2	." .i / / 2 :	(B10)
	Now,
	kx x .i/ x ?	k 2 2 Ä ky x .i/ x ?	k 2 2 D .i 1/
							2	kb ck 2 2 (B12)
	respectively. Combine the inequalities (B11b) and (B10) to get
	ƒ .i/	Ä ' .i/ ƒ .i 1/ Ä ' .i/ .ƒ .i 1/	C .' .i/ C 2ı .i 1/ = / C ." .i / / 2 C 1/ı .i 1/ C ." .i/ / 2 (B13a) (B13b)

Here, we use the extended Poisson probability mass function (pmf) Poisson.y j / D . y =yŠ/ e for all 0 by defining 0 0 D 1 to accommodate the identity-link model.

FISTA has been developed for the linear Gaussian model in Section II-B.

We select 1 D 2 as in[START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] because kˆk 2 and k‰k 2 have approximately the same scale in this example.

These choices of and r are at the boundary of the convergence region in[START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF] Th. 3.1]. We have searched for and r inside this convergence region as well, but found that the boundary choices that we select are the best, or close to the best.

Renliang Gu was born in Nantong, China. He re-

Acknowledgment

The authors are grateful to Prof. M. Hong, Iowa State University, for valuable discussions and to the anonymous reviewers for their helpful comments.

This work was supported by the U.S. National Science Foundation under Grant CCF-1421480.

Appendix A Derivation of Acceleration (20a)-(20d) and

Proofs of Lemma 1 and Theorem 1

We first prove Lemma 1 and then derive the acceleration (20a)-(20d) and prove Theorem 1.

Proof of Lemma 1: According to Definition 1 and (24), ur.x/ ur.x .i/ / C .x x .i/ / T h x x .i / x .i/ ˇ.i/ rL.x x .i / / i ." .i/ / 2 2ˇ. i / (A1a)

for any x 2 R p . Moreover, due to the convexity of L.x/, L.x/ L.x x .i/ / C .x x x .i / / T rL.x x .i/ /: (A1b) Summing (A1a), (A1b), and ( 21) completes the proof.

The following result from [START_REF] Bertsekas | Convex Optimization Algorithms[END_REF]Prop. 3.2.1 in Sec. 3.2] states that the distance between x and y can be reduced by projecting them onto a closed convex set C .

Lemma 2 (Projection theorem): The projection mapping onto a nonempty closed convex set C Â R p is nonexpansive

for all x; y 2 R p . We now derive the projected Nesterov's acceleration step (20b)-(20d) with the goal of selecting the x

x .i/ in the proximal step (20e) that achieves the convergence rate of O.k 2 /. This derivation and convergence-rate proof are inspired by-but are more general than- [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF]. We start from [START_REF] Raguet | A generalized forward-backward splitting[END_REF] with x replaced by x D x ? and x D x .i 1/ ,

and design two coefficient sequences a .i/ > 0 and b .i/ > 0 that multiply (A3a) and (A3b), respectively, which ultimately leads to (20a)-(20d) and the convergence-rate guarantee in (33a). Consider sequences a .i/ > 0 and b .i/ > 0. Multiply them by (A3a) and (A3b), respectively, add the resulting expressions, and multiply by ˇ.i/ to obtain

where

We arranged (A4) using completion of squares so that the first two summands are similar (but with opposite signs), with the goal of facilitating cancellations as we sum over i . Since we have control over the sequences a .i / and b .i/ , we impose the following boundary conditions for i 1:

where

Now, apply the inequality (A6a) to the right-hand side of (A4):

and sum (A8a) over i D 1; 2; : : : ; k, which leads to summand cancellations and

where (A8b) follows by discarding a nonnegative term c .k/ t .k/ . Now, due to .i/ .i/ 0 (see (34a) and (A6b)), the inequality (A8b) leads to

with simple upper bound on the right-hand side, thanks to summand cancellations facilitated by the assumptions (A6).

As long as ˇ.k/ c .k/ grows at a rate of k 2 and the inexactness of the proximal mappings leads to bounded P k i D1 c .i/ ." .i/ / 2 , the centered objective function .k/ can achieve the desired bound decrease rate of 1=k 2 .

In the following section, we show how to satisfy (A6a), which will lead to the projected momentum acceleration step (20d). We approach the constraints (A6a) by first aiming to meet them with equality, which is possible in the absence of the convexset constraint (C D R p ). We then use the nonexpansiveness of the convex-set projection to construct a .i / and b .i/ that satisfy (A6a) with inequality in the general case where the convex-set constraint is present. Finally, we show how to satisfy (A6b), which will allow us to construct the recursive update of  .i/ in (20b) and verify the allowed range of momentum tuning constants in [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF].

I Satisfying Conditions (A6)

a) Imposing equality in (A6a): (A6a) holds with equality for all i and any x ? when we choose x

x .i / D y x .i / that satisfies p c .i 1/ .x .i 1/ z .i 1/ / D p c .i/ .y x .i/ z .i/ /: (A10)

Now, (A10) requires equal coefficients multiplying x ? on both sides; thus a .i / = p c .i/ D 1=w for all i, where w > 0 is a constant (not a function of i), which implies c .i/ D w 2 .a .i / / 2 (A14a)

where (A14a) and (A14b) follow from (A10) and by using Lemma 2, respectively; see also (A5b). Without loss of generality, set w D 1 and rewrite and modify (A7), (A5b), and (A8b) using (A11) to obtain .i/ D ˇ.i/ . .i / / 2 ˇ.iC1/  .i C1/ . .iC1/ 1/; i 1 (A15a)

. .i/ / 2 t .i/ D  .i/ x .i/ . .i / 1/x .i 1/ x ? 2 2 (A15b)

where (A15c) is obtained by discarding the negative term 2ˇ. k/ . .k/ / 2 .k/ and the zero term ˇ.1/  .1/ . .1/ 1/ .0/ (because  .1/ D 1) on the left-hand side of (A8b). Now, (33a) follows from (A9) by using  .0/ D  .1/ D 1 (see (20b)), (A11), and (A15b) with i D 0. c) Satisfying (A6b): By substituting (A15a) into (A6b), we obtain the conditions

and interpret . .i/ / C1 iD1 as the sequence of gaps between the two sides of (A16a); (A16a) implies

Comparing (20b) with (A16b) justifies the constraints in [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF]. [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF] If the step-size sequence .ˇ. i/ / is non-increasing (e.g., in the backtracking-only scenario with n D C1), (20b) with B .i/ D 1 also satisfies the inequality (A16b). In this case, (33a) still holds but (33b) does not because [START_REF] Ollinger | Positron-emission tomography[END_REF] no longer holds. However, because

II Connection to Convergence-Rate Analysis of FISTA in

which generalizes [START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF]Th. 4.4] to include the inexactness of the proximal operator and the convex-set projection.

Appendix B Convergence of Iterates

To prove convergence of iterates, we need to show that the centered objective function .k/ decreases faster than the righthand side of (33b). We introduce Lemmas 3 and 4 and then use them to prove Theorem 2. Throughout this Appendix, we assume that Assumption 1 of Theorem 2 holds, which justifies (A3) and (A16) as well as results from Appendix A that we use in the proofs.

Lemma 3: Under Assumptions 1-3 of Theorem 2,

Proof: By letting k ! C1 in (A15c) and using (37), we obtain

For i 1, rewrite (A15a) using  .i/ expressed in terms of  .i C1/ (based on (20b)):

where the inequality in (B3) is due to b 2 1 < 0; see Assumption 3. Apply nonexpansiveness of the projection operator to (A3b) and use (A12) to obtain 2ˇ. i/ . .i 1/ .i/ / ı .i/ .' .i/ / 2 ı .i 1/ ." .i / / 2 I (B4) then multiply both sides of (B4) by . .i/ / 2 , sum over i D 1; 2; : : : ; k and reorganize:

where (see (A15a))

and we drop the zero term . .0/ 1/ 2 ı .0/ and the negative term . .k/ / 2 ı .k/ 2ˇ. k/ . .k/ / 2 .k/ from (B5a) and use the fact that % .i / Ä OE2=.

2/ .i / implied by (B3) to obtain (B5b).