PaperlayerView exact revision on arXiv

Source-linked AI summary

GESPAR: Efficient Phase Retrieval of Sparse Signals

Yoav Shechtman, Amir Beck, Yonina C. Eldar

arXiv:1301.1018v2cs.IT

TL;DR

Phase retrieval is ill-posed because Fourier magnitudes lose phase, so the paper targets recovery of sparse signals using prior sparsity information. It proposes GESPAR, a fast local-search method that avoids matrix lifting, and simulations report efficient, accurate recovery across varied settings.

  • Problem

    Fourier-magnitude phase retrieval is ill-posed because phase information is lost, motivating recovery methods that exploit signal sparsity.

  • Method

    GESPAR applies iterative 2-opt local search to a sparsity-constrained nonlinear optimization problem, updating the support and solving each support-restricted problem with damped Gauss–Newton iterations.

  • Results

    Simulations show that GESPAR is efficient and more accurate than current techniques, with robustness to noise and applicability to Fourier and more general quadratic measurements.

  • Takeaways & Limitations

    Because it avoids matrix lifting, GESPAR is potentially suitable for large-scale problems such as two-dimensional images.

Abstract

from arXiv · show

We consider the problem of phase retrieval, namely, recovery of a signal from the magnitude of its Fourier transform, or of any other linear transform. Due to the loss of the Fourier phase information, this problem is ill-posed. Therefore, prior information on the signal is needed in order to enable its recovery. In this work we consider the case in which the signal is known to be sparse, i.e., it consists of a small number of nonzero elements in an appropriate basis. We propose a fast local search method for recovering a sparse signal from measurements of its Fourier transform (or other linear transform) magnitude which we refer to as GESPAR: GrEedy Sparse PhAse Retrieval. Our algorithm does not require matrix lifting, unlike previous approaches, and therefore is potentially suitable for large scale problems such as images. Simulation results indicate that GESPAR is fast and more accurate than existing techniques in a variety of settings.

I. INTRODUCTION

Phase retrieval is ill-posed because Fourier magnitudes discard phase, motivating sparse-signal methods. GESPAR applies fast local search to sparse nonlinear optimization and is reported to outperform existing approaches in simulations.

  • I. INTRODUCTION: Phase retrieval recovers signals from Fourier-transform magnitudes, but the loss of phase makes one-dimensional recovery generally ill-posed.Prior information such as support, non-negativity, or signal magnitude is commonly used to address this issue.
  • I. INTRODUCTION: Existing sparse phase-retrieval methods include semidefinite-programming approaches and Fienup-type alternating-projection algorithms.The cited approaches involve matrix lifting or can suffer convergence problems, especially for one-dimensional signals.
  • I. INTRODUCTION: GESPAR applies a fast 2-opt local search method to a sparsity-constrained nonlinear optimization formulation.It iteratively updates the support and uses damped Gauss–Newton iterations to find a local minimum for the current support.
  • I. INTRODUCTION: GESPAR is reported to be efficient and more accurate than current techniques in numerical simulations.The simulations also examine robustness to noise and scalability, with reliable Fourier-magnitude recovery appearing to require measurements scaling like s^3.
  • I. INTRODUCTION: The formulation recovers a known s-sparse vector from magnitude-squared DFT measurements under a sparsity constraint, up to shift, global-phase, and mirroring ambiguities.The autocorrelation can provide support information when N ≥ 2n−1, but that information is not used for support estimation in noisy measurements.

B. Sparse Phase Retrieval: General Measurements

The paper extends its sparse phase-retrieval formulation beyond Fourier measurements. By changing the measurement matrices, the same approach accommodates arbitrary quadratic measurements and sparsity bases.

  • B. Sparse Phase Retrieval: General Measurements: The formulation applies to arbitrary quadratic measurements, not only Fourier measurements.The paper states that sparsity in a basis other than the identity basis is also supported.
  • B. Sparse Phase Retrieval: General Measurements: For arbitrary linear measurements, the quadratic formulation uses matrices Ai = φiφT_i derived from the measurement vectors.Here φi denotes the ith measurement vector.
  • B. Sparse Phase Retrieval: General Measurements: When x̄ = Dz is sparse in a basis D, the formulation uses Ai = DTφiφT_iD and therefore accommodates arbitrary sparsity bases.The same construction simultaneously represents general quadratic measurements.
  • B. Sparse Phase Retrieval: General Measurements: GESPAR solves the resulting formulation using an iterative local-search algorithm that assumes only symmetric measurement matrices.Fourier instances can be implemented more efficiently by exploiting Fourier structure.

III. GREEDY SPARSE PHASE RETRIEVAL (GESPAR)

The damped Gauss–Newton method solves the sparse-support nonlinear least-squares subproblem, while iterative updates and backtracking seek a stationary point under suitable conditions.

  • Damped Gauss–Newton Method: Damped Gauss–Newton minimizes the nonlinear least-squares objective over a fixed support.The method linearizes the terms in the objective, solves a linear least-squares problem, and updates with a backtracking-selected stepsize.
  • Damped Gauss–Newton Method: Under suitable conditions, the gradient norm converges to zero and every limit point is a stationary point of the objective.The convergence proof requires the Jacobian to have full column rank with uniformly bounded-below minimum eigenvalues.
  • Damped Gauss–Newton Method: The Gauss–Newton direction is computed from the Jacobian and least-squares solution, then applied with a backtracking stepsize.The implementation uses the update z_k = z_{k−1} − t_k d_k and stops when successive iterates are sufficiently close or the iteration limit is reached.
  • Support-Based Search: GESPAR initializes a random feasible support, invokes damped Gauss–Newton, and iteratively tests support swaps that improve the objective.Supports satisfy J1 ⊆ S0 ⊆ J2 and contain s indices; each accepted swap produces a new damped Gauss–Newton solution.
  • Support-Based Search: The algorithm accepts objective-decreasing updates and terminates when no tested swap improves the objective.Its implementation records the total number of swaps and returns the current solution when the stopping condition is met.

B. The 2-opt Local Search Method

The 2-opt local search repeatedly exchanges one support index with one off-support index, using damped Gauss–Newton to evaluate whether the objective improves.

  • 2-opt Local Search: The 2-opt method starts from a random support satisfying J1 ⊆ S ⊆ J2 and invokes damped Gauss–Newton.The support has size s, and the resulting solution initializes the local-search iterate.
  • 2-opt Local Search: Each iteration swaps the smallest-magnitude supported component with the largest-magnitude off-support gradient component.Only two support elements change, which gives the procedure its 2-opt designation.
  • 2-opt Local Search: A proposed swap is accepted when its damped Gauss–Newton solution decreases the objective function.The search advances through improving supports and stops when none of the tested swaps yields a better objective value.

C. The GESPAR Algorithm

GESPAR restarts the 2-opt search from different random supports to reduce the effect of local optima, stopping upon success or a swap budget limit.

  • GESPAR Algorithm: GESPAR repeatedly invokes 2-opt with different initial random support sets.The restart strategy addresses the tendency of 2-opt to become stuck at local optima.
  • GESPAR Algorithm: The algorithm stops successfully when the objective falls below a threshold and otherwise stops after the maximum total number of swaps.Random weights are regenerated for each damped Gauss–Newton invocation, and the authors observed fewer non-optimal local traps with this modification.

IV. FOURIER IMPLEMENTATION DETAILS

GESPAR can target sparse quadratic systems generally and becomes computationally efficient for Fourier measurements by using FFT-based calculations instead of explicitly forming measurement matrices.

  • General Sparse Quadratic Measurements: GESPAR seeks sparse solutions to systems of quadratic equations subject to an ℓ0 sparsity constraint.The formulation allows any symmetric measurement matrices representing the quadratic measurements.
  • Fourier Implementation: For efficiently implemented transforms, GESPAR avoids explicitly creating and storing the matrices representing the measurements.This implementation advantage is especially relevant for large vector sizes in two-dimensional Fourier phase retrieval.
  • Fourier Implementation: For Fourier measurements, the FFT computes the weighted objective, Jacobian-related quantities, and gradient used by the algorithm.The Jacobian computation involves only s columns of the Fourier matrix, and FFT2 can be used for two-dimensional problems.
  • Fourier Implementation: A 195 × 195 sparse image with s = 15 was recovered from 38,025 2D-Fourier magnitude measurements in 80 seconds.The example used 225 dictionary elements corresponding to non-overlapping circles on a 15 × 15 grid.

V. NUMERICAL SIMULATIONS

Numerical simulations evaluate GESPAR against SDP and Sparse-Fienup methods for recovery accuracy, runtime, and sensitivity to unknown sparsity. Across these tests, GESPAR achieves high recovery probability, faster runtime than SDP, and little sensitivity to using only an upper bound on sparsity.

  • A. Signal-recovery Accuracy: Over 90% successful recovery is achieved by GESPAR up to s = 15, compared with s = 8 and s = 7 for the other methods.Success probability is the ratio of correctly recovered signals across 100 simulations.
  • A. Signal-recovery Accuracy: GESPAR is faster than the SDP-based algorithm, while Fienup iterations are slightly slower and have a much lower success rate.Runtime is averaged over successful recoveries for n = 64 and N = 128.
  • B. Sensitivity to exact sparsity knowledge: Using only an upper limit of 25 for sparsity does not significantly reduce recovery probability compared with supplying the exact sparsity level.Performance is somewhat improved when more nonzero elements are allowed during the iterations.

C. Effect of the number of allowed swaps

GESPAR’s recovery probability improves as more support-index swaps are allowed, but computation time rises approximately linearly with swaps and eventually yields no additional recovery benefit. Simulations also examine oversampling, support information, noise robustness, and scalability with signal dimension.

  • Effect of the number of allowed swaps: Increasing the allowed number of index swaps increases recovery probability, but values above 6400 did not improve unsuccessful recoveries in the tested setting.The tested range was ITER = 100 to 25600, with N = 128 and s < 25 for the cited boundary.
  • Effect of oversampling and support information: Oversampling and autocorrelation-derived support information both improve GESPAR’s recovery performance for n = 64.The comparison used noiseless Fourier magnitude measurements with N = 64, 128, and 256, and separately used support information with N = 128.
  • Robustness to noise: GESPAR’s reconstruction performance improves as SNR increases and clearly outperforms sparse-Fienup in noise robustness.The comparison used normalized mean squared reconstruction error averaged over 100 random realizations.
  • Computation time and scalability: The maximal sparsity permitting successful recovery increases with vector length and appears to scale like n^1/3 for n ∈ [64, 2048].With n = 512, s = 35, N = 1024, and ITER = 6400, the mean reconstruction time was 33.5 seconds.
  • Computation time and scalability: Runtime scales approximately linearly with the number of swaps because each swap corresponds to one damped Gauss–Newton iteration.The matrix inversion process in the damped Gauss–Newton segment is the most time-consuming part of GESPAR.
  • Computation time and scalability: Beyond the successful recovery region, GESPAR uses the maximal 6400 swaps without obtaining a correct solution.The mean number of swaps depends on sparsity s and vector length n.

H. Two-Dimensional Fourier Phase Retrieval

The paper evaluates GESPAR for two-dimensional Fourier phase retrieval and large image-sized signals, comparing recovery probability and runtime with sparse-Fienup. GESPAR supports higher sparsity at large dimensions, while sparse-Fienup is faster for successful recoveries.

  • Two-Dimensional Fourier Phase Retrieval: GESPAR is applied to large-scale 2D Fourier phase retrieval problems with sparsities s ∈ [2, 82] and vector sizes n ∈ [256, 6400].The simulations use noiseless 2D DFT magnitudes without oversampling or autocorrelation-derived support information, with ITER = 6400.
  • Two-Dimensional Fourier Phase Retrieval: For n = 6400, GESPAR recovers signals with very high probability up to sparsity s = 57, compared with s = 42 for sparse-Fienup.The comparison is based on the recovery-probability curves for the two methods.
  • Two-Dimensional Fourier Phase Retrieval: Sparse-Fienup is faster than GESPAR for successful 2D recoveries as vector size n varies.The runtime comparison reports average computation time for a successful recovery as a function of n.
  • Conclusion: The paper concludes that GESPAR does not require matrix lifting and is potentially suitable for large-scale problems such as 2D images.Simulations also demonstrate successful treatment of phase-retrieval problems in one and two dimensions.

APPENDIX

The appendix develops the DGN method's least-squares step and establishes boundedness and convergence properties under smoothness and boundedness assumptions.

  • The vector-valued function h and its component functions h_i(z) define the least-squares formulation used by the method.
  • The least-squares solution yields a Gauss–Newton-style update involving (J(z_{k−1})^T J(z_{k−1}))^−1 J(z_{k−1})^T h(z_{k−1}).
  • The generated sequence is shown to be bounded because −d_k is a descent direction and the objective satisfies the stated descent inequalities.
  • Smoothness of g and continuity of J provide uniform bounds on the Hessian, J(z)^T J(z), and the gradient over the relevant ball.
  • The objective converges while ∇g(z_k) tends to zero, and all accumulation points of the sequence are stationary.
Loading 1301.1018v2…