Source-linked AI summary
This is SPIRAL-TAP: Sparse Poisson Intensity Reconstruction ALgorithms - Theory and Practice
Zachary T. Harmany, Roummel F. Marcia, Rebecca M. Willett
TL;DR
The paper addresses sparse Poisson inverse problems with nonnegative intensities and potentially undersampled observations. It proposes SPIRAL, a flexible penalized-likelihood framework using separable quadratic approximations and several regularizers. Experiments report strong reconstruction performance, with SPIRAL-TV outperforming state-of-the-art emission-tomography approaches and producing fewer spurious artifacts than wavelet-regularized methods.
Problem
The paper asks how to reconstruct sparse or compressible nonnegative intensities from Poisson data when unknowns may outnumber observations.
Method
SPIRAL minimizes a penalized negative Poisson log-likelihood under nonnegativity constraints using separable quadratic approximations and ℓ1, total-variation, or partition-based penalties.
Results
SPIRAL-TV outperforms state-of-the-art emission-tomography approaches and produces fewer spurious artifacts than wavelet-regularized methods.
Takeaways & Limitations
The framework supports flexible regularization choices, including convex total-variation reconstruction, for Poisson inverse problems.
Abstract
from arXiv · showhide
The observations in many applications consist of counts of discrete events, such as photons hitting a detector, which cannot be effectively modeled using an additive bounded or Gaussian noise model, and instead require a Poisson noise model. As a result, accurate reconstruction of a spatially or temporally distributed phenomenon (f*) from Poisson data (y) cannot be effectively accomplished by minimizing a conventional penalized least-squares objective function. The problem addressed in this paper is the estimation of f* from y in an inverse problem setting, where (a) the number of unknowns may potentially be larger than the number of observations and (b) f* admits a sparse approximation. The optimization formulation considered in this paper uses a penalized negative Poisson log-likelihood objective function with nonnegativity constraints (since Poisson intensities are naturally nonnegative). In particular, the proposed approach incorporates key ideas of using separable quadratic approximations to the objective function at each iteration and penalization terms related to l1 norms of coefficient vectors, total variation seminorms, and partition-based multiscale estimation methods.
I. INTRODUCTION
The paper studies sparse or compressible Poisson inverse problems with potentially more unknowns than observations, using a nonnegative penalized likelihood formulation and flexible sparsity penalties.
- Applications include nuclear medicine, night vision, astronomy, traffic analysis, and microscopy, where data are collected by counting discrete events.
- The goal is to estimate a sparse or compressible intensity f⋆ from Poisson observations when the number of unknowns may exceed the observations.
- Poisson compressed sensing differs from conventional compressed sensing because the sensing matrix and estimated intensity must be nonnegative.
- The proposed optimization minimizes a penalized negative Poisson log-likelihood subject to nonnegativity constraints.
- SPIRAL supports penalties based on canonical sparsity, arbitrary-basis coefficients, total variation, and recursive dyadic partition complexity.
B. Related work
Prior Poisson reconstruction methods use variance stabilization, Bayesian or multiscale penalties, and total variation, while SPIRAL provides a broadly applicable constrained optimization framework with convergence guarantees.
- Existing approaches include variance-stabilizing transforms, multiscale Bayesian and complexity-penalized likelihood methods, and total-variation regularization.
- Poisson inverse problems are difficult because the negative log-likelihood requires sophisticated optimization and nonnegativity creates constrained feasible sets.
- The proposed framework requires no special structure in the sensing matrix A and can exploit fast matrix-vector multiplications when available.
- The analysis establishes global convergence under mild practical assumptions and conditions for uniqueness, supported by limited-angle emission-tomography simulations.
- SPIRAL uses separable quadratic approximations with scaled-identity Hessians and constrained subproblems resembling nonnegative ℓ2 denoising.
- A modified Barzilai-Borwein step captures current Poisson curvature, while acceptance checks permit controlled nonmonotonic objective behavior.
A. Canonical basis with sparsity penalty
The paper solves sparsity-penalized subproblems in canonical and non-canonical bases, using direct thresholding in the canonical case and a feasible dual method otherwise.
- A. Canonical basis with sparsity penalty: With a canonical-basis ℓ1 penalty, the nonnegative minimization subproblem has an analytic solution, making the constrained update straightforward.
- B. Non-canonical basis: For an arbitrary orthonormal basis W, nonnegativity of f̂ = Wθ̂ means the subproblem generally lacks an analytic solution.
- B. Non-canonical basis: The non-canonical problem is reformulated with θ = u − v, doubling parameters while making the objective differentiable.
- B. Non-canonical basis: Its Lagrangian dual is solved by alternating updates for λ and γ under componentwise bounds.
- B. Non-canonical basis: Early termination of the dual iterations preserves feasibility but may produce a less accurate estimate; this variant is called SPIRAL-ℓ1.
C. Total variation penalty
The total-variation variant uses TV regularization to favor images with few abrupt intensity changes and homogeneous regions, reducing the problem to a nonnegative denoising subproblem.
- C. Total variation penalty: The anisotropic TV seminorm measures the magnitude of vertical and horizontal first-order pixel differences.
- C. Total variation penalty: Images with small TV are sparse in an overcomplete first-difference representation and tend to contain few abrupt changes with many homogeneous regions.
- C. Total variation penalty: Choosing the identity blur operator and nonnegative constraint converts the TV problem into the constrained ℓ2 denoising subproblem required by SPIRAL.
- C. Total variation penalty: The resulting total-variation algorithm is called SPIRAL-TV.
D. Partition-based methods
SPIRAL-RDP uses recursive dyadic partitions to construct nonnegative, multiresolution intensity estimates, with efficient dynamic-programming optimization and a translation-invariant variant for smoother results.
- D. Partition-based methods: Recursive dyadic partitions model intensities by fitting a nonnegative model within each partition cell.The paper illustrates piecewise constant fitting on an emission image.
- D. Partition-based methods: Large homogeneous regions are represented by large cells, while edges are approximated through finer partitioning.
- D. Partition-based methods: The partition search is computed efficiently using dynamic programming, with constant-cell fits obtained non-iteratively from clipped local averages.
- D. Partition-based methods: SPIRAL-RDP-TI applies cycle-spun translation invariance to produce a smoother estimator while retaining high computational efficiency.
- D. Partition-based methods: RDP penalties impose a hereditary wavelet-like constraint but are nonconvex, so the method guarantees convergence only to a local minimizer.With sufficiently accurate initialization, the RDP method is reported to perform competitively with convex penalties in speed and accuracy.
III. ALGORITHMIC DETAILS
SPIRAL uses acceptance-controlled iterative subproblem updates, with practical initialization and termination criteria designed for constrained Poisson reconstruction. The approach balances convergence behavior against computational cost and subproblem accuracy.
- III. ALGORITHMIC DETAILS: Checking acceptance criteria can dominate computation because each candidate requires recomputing Af^k+1, especially for dense unstructured matrices.
- III. ALGORITHMIC DETAILS: Initialization matters because finite runtimes make approximate solutions dependent on the starting point, while nonconvex RDP penalties can introduce difficult local optima.
- III. ALGORITHMIC DETAILS: Scaled ATy, Fourier inversion, and filtered back-projection are described as low-cost initialization strategies in different application settings.
- III. ALGORITHMIC DETAILS: The SPIRAL-ℓ1 subproblem uses a duality-gap termination criterion because its transformed objective is convex with affine constraints.
- III. ALGORITHMIC DETAILS: Early subproblem solves need not be highly accurate when the current iterate remains far from the optimum.
- III. ALGORITHMIC DETAILS: General termination checks monitor small changes in consecutive iterates or objective values, but may stop before reaching the true solution.The authors report not observing this premature termination in practice.
D. Convergence proof
The convergence analysis establishes global convergence of SPIRAL under mild assumptions for the constrained penalized Poisson objective. Objective values converge sublinearly in general and R-linearly when F is strongly convex.
- D. Convergence proof: The analysis assumes proper convexity and Lipschitz continuous differentiability of F, coercivity of Φ, and continuity of the penalty over the nonnegative domain.
- D. Convergence proof: The negative Poisson log-likelihood is Lipschitz continuously differentiable, with a computable gradient-Lipschitz bound based on A and the observations.The bound avoids requiring direct computation of the spectral norm of A.
- D. Convergence proof: Accepted iterates remain in the initial compact sublevel set, ensuring at least one convergent subsequence.
- D. Convergence proof: The step-size acceptance rule and supporting lemmas imply that iterate differences vanish and objective values converge to a finite limit.
- D. Convergence proof: All accumulation points are critical points, so SPIRAL converges to a minimizer of the convex objective under the stated assumptions.
- D. Convergence proof: Objective values converge sublinearly in general and R-linearly when F is strongly convex.
E. Uniqueness of the minimizer
The minimizer need not be unique in underdetermined problems with common nonsmooth penalties. Uniqueness can be established using properties of the forward operator and penalty, including an RIP condition for the ℓ1 formulation.
- E. Uniqueness of the minimizer: Under the general assumptions, the minimizers form a compact convex set that may contain more than one point.
- E. Uniqueness of the minimizer: In underdetermined settings, the likelihood is constant along directions in the kernel of A, while ℓ1 and TV penalties generally are not strictly convex.
- E. Uniqueness of the minimizer: The paper analyzes uniqueness for ℓ1 and TV penalties but excludes RDP because global-optimum convergence cannot be explicitly guaranteed.
- E. Uniqueness of the minimizer: If two solutions exist, they have identical projections through A and equal penalty values, making them equally faithful to the data and equally sparse or smooth under the penalty.
- E. Uniqueness of the minimizer: For ℓ1 regularization, an RIP condition on AW for sparsity s = 2|S| implies uniqueness of the coefficient solution.
A. Simulation Setup
The simulation reconstructs a known axial emission map from limited-angle Poisson tomographic data under known attenuation. Ten noisy data realizations support average performance evaluation, while figures show the image components and initialization.
- Simulation Setup: The experiment reconstructs a true axial emission map from simulated limited-angle emission tomography data.The emission and attenuation images are standard Image Reconstruction Toolbox test images.
- Simulation Setup: Figure 2 presents the true emission image, attenuation map, noisy projection data, and initialization estimate.The panels organize the simulated inputs and starting point for reconstruction.
- Simulation Setup: The known attenuation map modifies the projection model through A ≡ diag[exp(−Rµ)]R.The projection geometry uses 128 radial and 128 angular samples over 135 degrees.
- Simulation Setup: Ten Poisson-noise realizations are simulated to examine ten-trial average reconstruction performance.Only images from one particular realization are shown.
B. Algorithm Setup
The evaluation compares SPIRAL variants and competing Poisson or least-squares methods under a common initialization and iteration protocol. Algorithm parameters are tuned for best-case RMSE, with cross-validation identified as a practical alternative.
- Algorithm Setup: The evaluation includes SPIRAL-ℓ1, TV, RDP, and RDP-TI, two competing Poisson methods, and SpaRSA.SPIRAL-ℓ1 uses the DB-6 wavelet basis in the reported results.
- Algorithm Setup: SPS-OS uses a separable paraboloidal surrogate with ordered subsets, whereas EPL-INC-3 uses an incremental penalized Poisson-likelihood EM algorithm.Both are treated as competing Poisson reconstruction methods.
- Algorithm Setup: SpaRSA solves an ℓ2-ℓ1 regularized least-squares problem, and its output is thresholded to obtain a nonnegative feasible solution.This comparison evaluates the formulation using the Poisson likelihood.
- Algorithm Setup: All methods share the same initialization and run for at least 50 iterations, with global convergence declared below tolP = 5 × 10−4 relative iterate change.The initialization comes from two iterations of EPL-INC-3 starting from filtered back-projection.
- Algorithm Setup: Algorithm parameters are selected to minimize reconstruction RMSE for best-case comparisons, although practical settings can use cross-validation.The reported RMSE is 100 · ∥bf − f⋆∥2/∥f⋆∥2.
C. Results analysis
The results show that SPIRAL-TV provides the strongest reconstruction quality, while relaxed convergence settings can substantially improve computational speed with limited RMSE loss. Convergence behavior and regularization choice create a clear accuracy–efficiency trade-off across methods.
- Reconstruction quality: SPIRAL-TV achieves the lowest RMSE and highest visual quality among the proposed approaches, followed by SPIRAL-RDP-TI and SPIRAL-ℓ1.The comparison also indicates that image structure beyond basis sparsity can improve reconstruction quality.
- Reconstruction quality: SPIRAL-TV and EPL-INC-3 best preserve strong edges while smoothing homogeneous regions, with SPIRAL-TV more effectively suppressing homogeneous-region variation.SPIRAL-RDP-TI and ℓ1-based methods tend to over-smooth, while EPL-INC-3 exhibits intensity variations in homogeneous regions.
- Convergence and computation: SPIRAL-ℓ1 and RDP-TI terminate in under eight seconds after sharp initial RMSE decreases, whereas EPL-INC-3 reaches termination after 20 seconds with slower improvement.SPIRAL-TV also decreases sharply at first but then enters a prolonged period of little RMSE change.
- Convergence and computation: Relaxing near-monotonicity accelerates SPIRAL-TV convergence but can cause large transient RMSE increases during execution.The nonmonotonic variants typically show early sharp RMSE decreases, whereas the monotonic variant achieves a more accurate solution at termination in ten-trial averages.
- Convergence and computation: Combining nonmonotonic updates with loose subproblem convergence makes SPIRAL-TV approximately four times faster while nearly matching the best EPL-INC-3 method in RMSE.The relaxed criteria support faster convergence, with little loss in final RMSE according to Table I.
- Overall comparison: The total-variation regularized method outperforms current state-of-the-art emission-tomography approaches, producing fewer spurious artifacts than wavelet-regularized methods.Unlike partition-regularized methods, it is obtained through a convex optimization procedure.