PaperlayerView exact revision on arXiv

Source-linked AI summary

Phase Retrieval via Wirtinger Flow: Theory and Algorithms

Emmanuel Candes, Xiaodong Li, Mahdi Soltanolkotabi

arXiv:1407.1065v3cs.ITmath.FAmath.NAmath.OCmath.ST

TL;DR

The paper addresses recovery of a complex signal from phaseless quadratic measurements, a nonconvex phase retrieval problem. It introduces Wirtinger Flow, combining spectral initialization with iterative gradient-like updates, and proves exact recovery for random measurements with geometric convergence. The analysis also supports a coded-diffraction variant and yields insights into nonconvex optimization.

  • Problem

    Phase retrieval asks how to reconstruct x ∈ C^n from magnitude-only measurements y_r = |⟨a_r,x⟩|^2, without directly observing measurement phases.

  • Method

    Wirtinger Flow combines a spectral initialization with iterative updates that refine the estimate in a gradient-descent-like scheme.

  • Results

    The algorithm's iterates converge geometrically to a solution under the paper's random sampling conditions, with initialization obtained from a leading-eigenvector spectral method.

  • Takeaways & Limitations

    The proposed approach enables exact phase retrieval from a nearly minimal number of random measurements and has a near-linear-time coded-diffraction variation.

  • Takeaways & Limitations

    The main theorem requires m ≥ c0·n log n measurements, leaving whether the result extends to the linear-sampling regime open.

Abstract

from arXiv · show

We study the problem of recovering the phase from magnitude measurements; specifically, we wish to reconstruct a complex-valued signal x of C^n about which we have phaseless samples of the form y_r = |< a_r,x >|^2, r = 1,2,...,m (knowledge of the phase of these samples would yield a linear system). This paper develops a non-convex formulation of the phase retrieval problem as well as a concrete solution algorithm. In a nutshell, this algorithm starts with a careful initialization obtained by means of a spectral method, and then refines this initial estimate by iteratively applying novel update rules, which have low computational complexity, much like in a gradient descent scheme. The main contribution is that this algorithm is shown to rigorously allow the exact retrieval of phase information from a nearly minimal number of random measurements. Indeed, the sequence of successive iterates provably converges to the solution at a geometric rate so that the proposed scheme is efficient both in terms of computational and data resources. In theory, a variation on this scheme leads to a near-linear time algorithm for a physically realizable model based on coded diffraction patterns. We illustrate the effectiveness of our methods with various experiments on image data. Underlying our analysis are insights for the analysis of non-convex optimization schemes that may have implications for computational problems beyond phase retrieval.

1 Introduction

The paper frames phase retrieval as recovering a complex signal from magnitude-only quadratic measurements, a problem arising broadly in science and engineering. Existing methods can work in practice but often rely on signal priors and lack fully resolved convergence theory.

  • Problem setting: The measurement model is a generalized nonconvex quadratic program with broad relevance to science and engineering.Related applications include X-ray crystallography, microscopy, astronomy, diffraction imaging, and optics.
  • Problem setting: Phase retrieval seeks to recover x ∈ C^n from measurements y_r = |⟨a_r,x⟩|^2, which record magnitudes rather than phases.In coherent diffraction imaging, the model represents an object's Fourier-magnitude data.
  • Existing methods: Error reduction methods iteratively alternate projections in data space and signal space to match observed Fourier magnitudes and impose prior structure.Typical priors include real-valuedness, nonnegativity, and spatial limitation.
  • Existing methods: These algorithms often work well empirically but appear to depend heavily on a priori signal information.The paper identifies this dependence as a limitation of the existing approach.
  • Existing methods: Because the Fourier-domain constraint is nonconvex, fundamental convergence questions for alternating-projection methods remain largely unresolved.The unresolved theory concerns algorithms cast as projections onto nonconvex sets.
  • Related optimization context: Several related combinatorial optimization problems are also quadratic programs, but their general forms are notoriously difficult and often addressed with heuristics or semidefinite relaxations.For certain random models, PhaseLift-type relaxations can provide exact solutions up to global phase.

2 Algorithm: Wirtinger Flow

Wirtinger Flow combines a spectral initialization with iterative Wirtinger-gradient updates to minimize a non-convex phase-retrieval objective. Under a random model and sufficiently accurate initialization, its iterates converge toward a solution, while the stochastic-gradient interpretation motivates adaptive step sizes.

  • Algorithm overview: The algorithm first computes an initial estimate using a spectral method, then refines it through novel iterative updates resembling gradient descent.The two stages together are called the Wirtinger Flow algorithm.
  • Non-convex objective: The formulation minimizes a loss-based non-convex objective over z ∈ C^n, using quadratic loss ℓ(x,y) = (x − y)^2.The loss measures scalar misfit, but the resulting objective is not convex and can have many stationary points.
  • Iterative updates: The WF update is a steepest-descent step based on Wirtinger derivatives, with µτ+1 interpreted as a step size.For complex-valued variables, the objective is not holomorphic, so ordinary complex differentiation does not apply.
  • Initialization: The spectral initializer is the leading eigenvector of a positive semidefinite Hermitian matrix constructed from sampling vectors and observations.With A formed from the sampling vectors, the initializer can be computed by repeatedly applying A, entrywise multiplication by y, and A* via the power method.
  • Convergence and step sizes: The main result guarantees convergence toward a solution when the initialization is sufficiently accurate, while WF can be interpreted as a stochastic-gradient scheme.Gradient-estimate uncertainty is larger farther from the solution, motivating smaller early step sizes and larger later ones.
  • Initialization: In the Gaussian model, the population initialization matrix is I + 2xx*, whose leading eigenvector recovers x up to a global sign or phase with infinitely many samples.Finite data yields an initializer that is not perfect but is hoped to remain sufficiently correlated with x.

3 Main Results

The paper establishes recovery guarantees for Wirtinger flow under Gaussian measurements and compares its sampling, convergence, and computational properties with alternative non-convex schemes.

  • Exact phase retrieval via Wirtinger flow: Under the Gaussian model, Wirtinger flow is proved correct, and a diffraction-model extension provides exact recovery results as well.The theorem requires Gaussian sampling; the paper later develops exact recovery results for a physically inspired diffraction model.
  • Exact phase retrieval via Wirtinger flow: m ≥ c0·n log n measurements suffice for the theorem, so WF loses at most a logarithmic factor relative to the injectivity threshold 4n.The paper notes that 2n measurements are necessary in principle and that injectivity holds for m ≥ 4n in the sampling model.
  • Exact phase retrieval via Wirtinger flow: With µ = c1/n, WF reaches relative accuracy ϵ in O(n log 1/ϵ) iterations and has computational complexity O(mn^2 log 1/ϵ).A theoretical modification retains the same sampling complexity while reducing computational complexity to O(mn log 1/ϵ).
  • Comparison with other non-convex schemes: Gerchberg-Saxton lacks a clear convergence rate and initialization theory, whereas WF is shown to converge geometrically in the Gaussian model without solving a least-squares problem at each iteration.The comparison also identifies reduced per-iteration computational complexity for WF.
  • Comparison with other non-convex schemes: AltMinPhase uses disjoint sample blocks and least-squares updates, requiring c·n log n·(log^2 n + log 1/ϵ) samples for ϵ accuracy under its stated guarantees.The paper argues that partitioning samples and using each measurement only once makes this design of little practical value.

4 Numerical Experiments

The experiments evaluate Wirtinger flow on Gaussian and coded diffraction measurements, random signals, and natural images. The algorithm achieves high-probability recovery with few measurements and high precision while remaining computationally practical.

  • Experimental setup: Experiments primarily assess Wirtinger flow using coded diffraction patterns, with additional Gaussian measurement-model tests.The coded diffraction model is physically realizable in optical applications, particularly microscopy.
  • Experimental setup: The coded diffraction model uses random modulation patterns, including octanary patterns with eight possible values.For each modulation pattern, the signal is modulated before collecting a diffraction-pattern magnitude, producing multiple coded views.
  • Experimental setup: The tests use random low-pass and random Gaussian complex signals, with n = 128 and one signal of each type reused across experiments.The low-pass model is a bandlimited version of the Gaussian model with matched expected signal power.
  • Random-signal recovery: Around 4.5n Gaussian phaseless measurements suffice for exact recovery with high probability, while about six octanary patterns are sufficient.Success is declared when the relative reconstruction error falls below 10^-5, based on 100 random trials.
  • Natural-image recovery: On color photographs, 20 random octanary patterns yield 12 or 13 digits of precision within minutes.The image experiments use 50 power-method initialization iterations and 300 update iterations for each RGB image.
  • Computational cost: For about 13 digits of accuracy, the image experiments require 21,000–42,000 FFT units, within a factor of 1.5–3 of the estimated optimum.The workload is roughly 350 applications of the sensing matrix and its adjoint.
  • Computational cost: Semidefinite-programming methods such as PhaseLift and PhaseCut are impractical on a laptop because lifting an n-pixel image produces an approximately n^2/2-sized array.The smallest cited example requires approximately 85 GB of storage, while the Milky Way image would require about 17 TB.

5 Theory for the Coded Diffraction Model

For coded diffraction patterns, a resampled Wirtinger Flow variation achieves exact recovery with high probability under logarithmic pattern and sample requirements, while attaining near-linear computational complexity.

  • Algorithm: The coded diffraction variation uses resampled Wirtinger Flow iterations, partitioning observations into independent groups for initialization and refinement.Each group corresponds to a distinct random code realization, producing stochastic independence across groups.
  • Recovery guarantees: L ≥ c0(log n)^4 admissible coded diffraction patterns suffice for the theorem’s initialization and convergence guarantees.The theorem assumes a sufficiently large numerical constant c0 and gives high-probability events for initialization and subsequent convergence.
  • Recovery guarantees: The method achieves perfect recovery from on the order of n(log n)^4 coded diffraction samples.The paper compares this rate with contemporaneous PhaseLift guarantees and identifies further sampling reduction as future work.
  • Computational complexity: Setting µ = c1 yields ϵ accuracy in O(log 1/ϵ) iterations and overall computational work at most O(nL log n log 1/ϵ) in the CDP model.The per-iteration work is dominated by matrix-vector products involving A and A∗.

6 Wirtinger Derivatives

The paper develops Wirtinger calculus for real-valued functions of complex variables, using conjugate coordinates to express gradients, Hessians, and descent updates in a form analogous to real-variable calculus.

  • Complex coordinates: Complex variables can be represented through conjugate coordinates z and z̄, which formally substitute for real coordinates x and y.This representation applies to differentiable real-valued component functions even when the complex function is not holomorphic.
  • Derivatives: The complex gradient and Hessian are defined in the conjugate-coordinate system, with gradients represented as column vectors and derivatives as row vectors.The Hessian is organized into blocks involving z and z̄ derivatives.
  • Gradient updates: The resulting update can be interpreted as gradient descent in conjugate coordinates for real-valued functions over complex variables.The paper uses this interpretation to connect the Wirtinger update to steepest descent.
  • Gradient updates: For real-valued objectives, the second set of gradient updates is the conjugate of the first, so only one update needs to be tracked.This reduces the bookkeeping required for gradient-based optimization over complex variables.
  • Complex coordinates: Wirtinger derivatives simplify derivative expressions for real-valued functions of complex variables and make them resemble derivatives in the real case.The paper contrasts this with directly representing z = x + iy as a vector in R^2n.

7 Proofs

The proofs establish that local curvature and smoothness control the Wirtinger Flow dynamics, while concentration results verify these properties for Gaussian and coded diffraction measurements.

  • Proof normalization: The proof analysis uses a rescaled update equivalent to the WF update, allowing convergence results to transfer after adjusting the learning-parameter bound.The rescaled update is invariant to the Euclidean norm of x, so the proofs assume ∥x∥ = 1 without loss of generality.
  • Concentration results: Hessian concentration holds with high probability when m ≥ c(δ)n log n in the Gaussian model or L ≥ c(δ)log^3 n in the CDP model.The corresponding probability bounds are stated separately for the Gaussian and CDP settings.
  • Concentration results: Gradient concentration is established for fixed vectors near the solution, with high-probability bounds in both Gaussian and admissible CDP models.The lemma applies when dist(z,x) ≤ 1 and z is independent of the sampling vectors.
  • Convergence framework: Under the regularity condition, Wirtinger Flow initialized near the global minimizer set converges to a global optimizer at a geometric rate.The regularity condition is imposed on E(ϵ), the set of points close to the path of global minimizers.
  • Local conditions: The local curvature condition requires the objective to curve sufficiently upward near the global minimizer path.This condition is combined with local smoothness to establish the regularity condition.
  • Local conditions: The local smoothness condition requires the gradient to remain well behaved near the curve of global optimizers.Together with local curvature, it supplies the two ingredients used to derive the regularity condition.

A.1 Proof of Lemma 7.1

The appendix proof for Lemma 7.1 supplies Gaussian-model identities used in the main concentration and convergence analysis.

  • Gaussian identities: For the Gaussian model, the proof invokes two auxiliary lemmas for fixed vectors x.These lemmas provide model-specific identities whose proofs are omitted in the cited passage.
  • Gaussian identities: One resulting Gaussian identity gives E[ar ar^T] = 2xx^T.The identity is stated as part of the appendix calculations for the Gaussian sampling model.

A.2 Proof of Lemma 7.2

The proof invokes separate prior lemmas for the coded diffraction pattern and Gaussian measurement models.

  • The CDP model is handled by applying Lemma 3.1 from reference [16].
  • The Gaussian model is handled by applying Lemma A.1 stated above.
  • These lemma applications provide the stated conclusion.

A.3 Proof of Lemma 7.3

The proof analyzes Gaussian measurements using unitary invariance, moment identities, concentration inequalities, and net arguments to establish high-probability bounds.

  • Unitary invariance reduces the Gaussian analysis to vectors aligned with the first coordinate.The proof takes v=e1 and represents u using its first two coordinates.
  • Standard normal moment calculations establish the identity used in the proof.
  • Chebyshev’s, Hoeffding’s, and Bernstein-type inequalities provide the intermediate high-probability bounds.The resulting bounds include probabilities at least 1−3n−2, 1−3e−2γn, and 1−2e−2γn.
  • The proof controls expressions uniformly over unit vectors by decomposing vectors into first-coordinate and remaining components.
  • For sufficiently large C, m≥C⋅nlog n yields (A.2) with probability at least 1−5e−γn−4n−2.The argument uses an ϵ-net, a union bound, and suitable choices of δ0, ϵ, and γ.
  • The Hessian bound follows by combining concentration around the mean with an upper semidefinite inequality.The proof uses ∥∇2f(x)−E[∇2f(x)]∥≤δ to obtain ∇2f(x)⪯E[∇2f(x)]+δI.

A.6 Proof of Corollary 7.6

The proof of the corollary combines previously established concentration and derivative bounds, with the remaining argument obtained by symmetry.

  • The proof begins from a high-probability statement established in the proof of Lemma 7.4.
  • Gradient concentration is established by bounding directional deviations ⟨u,∇f(z)−E∇f(z)⟩.
  • Lemma 7.2 supplies an additional bound used in the argument.
  • Expressions (A.7) and (A.8) are combined using the triangle inequality and Lemma 7.4.

A.8 Proof of Lemma 7.8

The result is established separately for the coded diffraction pattern and Gaussian models using model-specific or standard covariance-deviation results.

  • For the CDP model, the result follows from Lemma 3.3 in reference [16].
  • For the Gaussian model, the result follows from standard sample-covariance deviation results.The cited example is Theorem 5.39 in reference.
  • The Gaussian argument concerns deviation of the sample covariance matrix from its mean.

B The Power Method

The power method initializes phase retrieval by computing the leading eigenvector of Y, repeatedly applying Y and normalizing. Its accuracy improves geometrically, and a constant eigenvalue gap implies few iterations suffice.

  • The method uses a randomly initialized power iteration to compute the first eigenvector of Y = Adiag{y}A∗.The input is matrix Y, and the initial vector is sampled uniformly from the unit sphere of C^n.
  • In the Gaussian model, each power-method iteration costs 2mn multiplications, while the CDP model costs 2L n-point FFTs.
  • Each iteration multiplies the current vector by Y and normalizes the result, returning the final iterate as the initialization ˜z0.
  • After k iterations, eigenvector accuracy is O(tanθ0(λ2/λ1)^k), where λ1 and λ2 are the two largest eigenvalues and θ0 is the initial angle.
  • A constant lower bound on the eigenvalue gap means that roughly log(n/ϵ)/log(λ1/λ2) iterations yield ϵ-accurate initialization.
Loading 1407.1065v3…