Solving Random Quadratic Systems of Equations Is Nearly as Easy as Solving Linear Systems

Yuxin Chen; Emmanuel J. Candes

Source-linked AI summary

Solving Random Quadratic Systems of Equations Is Nearly as Easy as Solving Linear Systems

Yuxin Chen, Emmanuel J. Candes

arXiv:1505.05114v2May 19, 2015cs.ITcs.LGmath.NAmath.STstat.ML

TL;DR

The paper addresses solving quadratic intensity-measurement systems arising in phase retrieval, where sensors record intensities rather than diffraction phases. It develops truncated Wirtinger Flow with spectral initialization and adaptive updates, proving exact recovery and linear-time computational cost under Gaussian measurements, with stability in noise.

Problem
Quadratic intensity measurements create phase-retrieval systems because optical sensors record y = |Ax|^2 rather than the diffraction pattern Ax.
Method
The method combines spectral initialization with a nonconvex Poisson log-likelihood optimization whose gradient updates truncate excessively influential components.
Results
TWF achieves exact recovery from O(n) quadratic equations with total cost O(mn log(1/ϵ)) and provides stability guarantees in the noisy case.
Takeaways & Limitations
Under independent Gaussian designs, solving random quadratic systems can have linear computational cost and nearly minimal mean-square loss in noise.
Takeaways & Limitations
The theory leaves the CDP measurement model for future research because it offers less randomness for the mathematical analysis.

Abstract

from arXiv · show

We consider the fundamental problem of solving quadratic systems of equations in $n$ variables, where $y_i = |\langle \boldsymbol{a}_i, \boldsymbol{x} \rangle|^2$, $i = 1, \ldots, m$ and $\boldsymbol{x} \in \mathbb{R}^n$ is unknown. We propose a novel method, which starting with an initial guess computed by means of a spectral method, proceeds by minimizing a nonconvex functional as in the Wirtinger flow approach. There are several key distinguishing features, most notably, a distinct objective functional and novel update rules, which operate in an adaptive fashion and drop terms bearing too much influence on the search direction. These careful selection rules provide a tighter initial guess, better descent directions, and thus enhanced practical performance. On the theoretical side, we prove that for certain unstructured models of quadratic systems, our algorithms return the correct solution in linear time, i.e. in time proportional to reading the data $\{\boldsymbol{a}_i\}$ and $\{y_i\}$ as soon as the ratio $m/n$ between the number of equations and unknowns exceeds a fixed numerical constant. We extend the theory to deal with noisy systems in which we only have $y_i \approx |\langle \boldsymbol{a}_i, \boldsymbol{x} \rangle|^2$ and prove that our algorithms achieve a statistical accuracy, which is nearly un-improvable. We complement our theoretical study with numerical examples showing that solving random quadratic systems is both computationally and statistically not much harder than solving linear systems of the same size---hence the title of this paper. For instance, we demonstrate empirically that the computational cost of our algorithm is about four times that of solving a least-squares problem of the same size.

1 Introduction

The paper introduces Truncated Wirtinger Flow (TWF), a vector-domain nonconvex method for recovering signals from magnitude-only quadratic measurements. Under independent Gaussian designs, it combines exact noiseless recovery, linear-time complexity, and near-optimal noisy-data accuracy with favorable empirical performance.

Problem and approach: TWF solves magnitude-only quadratic systems by preserving vector iterates and directly optimizing a nonconvex objective instead of lifting to matrix variables.The approach follows Wirtinger Flow but avoids the higher-dimensional matrix formulation.
Problem and approach: Adaptive truncation trims overly influential measurements during initialization and gradient updates, producing tighter initial guesses and more stable search directions.The iteration-dependent subsets Tt exclude components that are excessively large.
Numerical performance: 4 times slower than conjugate gradient is the reported empirical cost of solving a random quadratic system relative to a least-squares problem of the same size.The comparison uses m = 8n, Gaussian data, and TWF initialized with 10 power iterations.
Theoretical guarantees: O(n) quadratic measurements suffice for exact recovery, and TWF reaches ε-accuracy in O(log(1/ε)) iterations with total cost O(mn log(1/ε)).These guarantees improve on WF's O(n log n) sample complexity and O(mn^2 log(1/ε)) computational cost.
Theoretical guarantees: TWF's noisy-data error shrinks geometrically to a noise-dependent floor, while its Poisson-model accuracy matches the minimax lower bound.The lower bound states that no estimator can achieve vanishing error in the proportional-measurement regime under the stated signal-energy condition.

2 Algorithm: Truncated Wirtinger Flow

Truncated Wirtinger Flow (TWF) combines spectral initialization with adaptive trimming of gradient and initialization components. The procedure is designed to stabilize nonconvex updates and improve initialization in random quadratic systems.

Regularized gradient stage: Large gradient components can have excessive leverage and cause WF updates to overshoot, especially when |a_i^T z| is unusually small.The resulting weights can become as large as m||x||^2/||z|| in the real-valued setting.
Regularized gradient stage: TWF discards high-leverage gradient components, accepting a small bias to reduce variance and stabilize the descent direction.The trimming rules regularize both numerator and denominator terms in the component weights.
Truncated spectral initialization: The original spectral method can require m ≳ n log n for a good initialization because heavy-tailed observations can dominate its empirical sum.A truncated spectral method instead discards observations several times larger than the mean before computing the leading eigenvector.
Truncated spectral initialization: Under the reported Gaussian and coded-diffraction experiments, truncated spectral initialization maintains nearly identical accuracy as dimension grows, unlike the original spectral method.Figure 6 averages relative error over 50 Monte Carlo runs with m = 6n for Gaussian measurements and m = 12n for coded diffraction patterns.
Implementation choices: A fixed step size works under the stated conditions for any constant 0 < μ < 0.28, while backtracking line search is another supported option.The convergence rate depends on the fixed step size.

3 Why TWF works?

TWF works by controlling gradient components so that the update direction remains aligned with the error near the planted solution. Its regularity condition yields local contraction without requiring strong convexity.

Gradient intuition: The untrimmed gradient contains a useful component proportional to −4(a_i^T h)a_i and a potentially unbounded residual term.The residual becomes problematic when a_i^T z is arbitrarily small.
Gradient intuition: Adaptive trimming removes abnormal gradient components, leaving most of the regularized gradient aligned with the deviation from the current iterate to the solution.The trimming introduces a small bias when only a small fraction of observations is discarded.
Regularity condition: The regularity condition states a local relationship between the gradient and the error within an ε-ball around the planted solution.This neighborhood forms a basin of attraction for the iterative procedure.
Regularity condition: When the regularity condition holds, the TWF update is locally contractive around the planted solution.The contraction applies for iterates satisfying the stated neighborhood condition.
Regularity condition: Conditions (37) and (38) imply regularity with μ and λ of order one when ||h||/||z|| is sufficiently small, enabling a constant contraction rate 1 − μλ.This supports use of a constant step size.
Regularity versus convexity: Regularity does not necessarily imply strong convexity because it constrains pairs relative to the fixed planted solution rather than all local point pairs.The paper illustrates this distinction with a non-convex function satisfying the regularity condition.

4 Numerical experiments

Numerical experiments evaluate TWF on noiseless and noisy random quadratic systems across Gaussian, complex Gaussian, and coded-diffraction measurements. They report strong recovery behavior, improved performance over WF, and inverse-SNR error scaling.

Noiseless recovery: For real Gaussian measurements, TWF’s empirical success rate is evaluated at n = 1000 as m varies from 2n to 6n, using a relative-error threshold of 10^-5.Each value is averaged over 100 Monte Carlo trials.
Noiseless recovery: m > 4.5n and m ≥ 6n are often sufficient for complex Gaussian and coded-diffraction models, respectively.The reported success rates average 100 Monte Carlo trials.
Comparison with WF: TWF’s empirical success rates outperform WF after 1000 iterations under the reported Gaussian settings.The comparison suggests faster convergence or improved phase-transition behavior, without distinguishing between them.
Noisy recovery: Relative MSE scales inversely with SNR across the tested numbers of equations, matching the stability guarantee’s approximately −1 slope on the dB scale.The noisy experiments use real Gaussian designs and Poisson noise, averaging over 100 independent runs.

5 Exact recovery from noiseless data

The noiseless analysis establishes that truncated Wirtinger flow has a locally contractive update and achieves exact recovery under random measurements once the equation-to-variable ratio is sufficiently large.

A truncated spectral initializer is accurate with high probability when m/n exceeds a numerical constant, enabling local convergence analysis.
The local error contraction proposition guarantees that iterates remain near x and approach it at a geometric rate.
The regularity-condition proof controls the descent direction by separating dominant, higher-order, and rare high-influence terms.

6 Stability

The stability analysis shows geometric error reduction until a noise-limited regime, where iterates remain within an error scale determined by the noise magnitude.

The analysis establishes stability guarantees for the iterative gradient stage while deferring initialization analysis to the appendix.
Regime 1: In Regime 1, the estimation error is contractive and decreases geometrically under suitable regularity conditions.
Regime 2: In Regime 2, the error may cease to contract, but each update has bounded size and cannot increase the error excessively.
Regime 2: If an iterate leaves Regime 2 while the noise is sufficiently small, it returns to Regime 1 rather than escaping the controlled neighborhood.
Overall, the error shrinks geometrically before entering Regime 2 and thereafter remains on the order of ∥η∥√m∥x∥.

7 Minimax lower bound

The minimax analysis constructs many separated hypotheses near the true signal, bounds their information divergence, and derives a lower bound on estimation accuracy.

The proof targets a minimax lower bound using conditional KL-divergence inequalities and a Fano-type argument.
A collection of exp(n/30) separated hypotheses is constructed near x under Gaussian design vectors.
The hypotheses are rescaled so their pairwise distances are about δ while remaining clustered around x.
The KL-divergence bound yields a lower-bound scale involving n/m after choosing δ appropriately.
Because the constructed hypotheses lie much closer to x than its norm, Euclidean estimation error agrees with the phase-invariant distance.

8 Discussion

The discussion identifies extensions beyond the analyzed Poisson and Gaussian-measurement settings and points to low-rank recovery applications for modified TWF schemes.

More general objective functions: The analysis framework is stated to extend from the Poisson log-likelihood to a broad class of nonconvex objectives, including the Gaussian log-likelihood.
Sub-Gaussian measurements: Sub-Gaussian measurement vectors are within the theory's scope, although the truncation threshold may need adjustment for their sub-Gaussian norm.
Sub-Gaussian measurements: The coded diffraction pattern case is identified as more challenging because it offers less randomness for mathematical analysis and is left for future research.
Low-rank extensions: A modified TWF scheme is proposed as a possible linear-time approach for recovering low-rank positive semidefinite matrices from rank-one measurements.
Low-rank extensions: Potential applications include matrix completion, robust PCA, latent-variable models, and alignment or matching problems.

A Proofs for Section 5

The proof solves the relevant quadratic inequality separately across sign regimes and derives matching inner and outer bounds. Homogeneity reduces the argument to unit vectors before the final case analysis.

The argument first simplifies the quadratic inequality defining D_i before deriving equivalent interval descriptions.
Under assumption (112), the proof obtains inner and outer bounds from the preceding observations.
The same inner and outer bounds are derived when a_i^⊤z < 0, completing the sign-dependent analysis.

A.2 Proof of Lemma 4

The proof establishes uniform concentration through Lipschitz truncations, tail inequalities, epsilon-net discretization, and union bounds. These arguments yield high-probability control over unit vectors and parameter ranges when the sample size is sufficiently large relative to dimension.

Hoeffding- and Bernstein-type inequalities control sub-Gaussian and sub-exponential terms arising from the truncated functions.
An epsilon-net discretizes the unit sphere, after which union bounds extend fixed-vector concentration estimates uniformly.
The proof combines Lemmas 1 and 2 to obtain high-probability bounds for all unit vectors.
The proof restricts attention to unit vectors and replaces discontinuous indicators with bounded Lipschitz auxiliary functions.
For the relevant parameter range, the resulting uniform estimates hold once m exceeds a sufficiently large constant multiple of n.
The remaining terms are bounded using the inclusion (63), a probability estimate of 1 − exp(−Ω(m)), and the relation γ4 = 3αh.

B.1 Proof of Lemma 8

The proof constructs exponentially many candidate vectors near x using random packing, then shows that sufficiently many satisfy the required inequalities. Independence, concentration, union bounds, and conditional counting establish the collection required by Lemma 8.

A random packing constructs an exponentially large set M1 whose pairwise distances concentrate in a common range.
The vectors are decomposed into parallel and orthogonal components, and indicator variables divide the constraints into two control groups.
The first indicator group imposes stringent per-equation conditions, yet retains exponentially many candidates because M1 is exponentially large.
Conditioning on the first group preserves independence from the remaining sensing vectors, enabling a union-bound argument for the second group.
With probability 1 − o(1), at least a specified collection of vectors in M2 satisfies all required indicator conditions and forms a valid set M.
The auxiliary Gaussian matrix argument uses k = m/(4 log m) and concentration for sufficiently large m.

C Initialization via truncated spectral Method

The truncated spectral method is shown to work when m is proportional to n. Concentration of the truncated data matrix and Davis–Kahan perturbation control yield an accurate initial spectral estimate.

The section demonstrates that the truncated spectral method works when m ≍ n.
For the stated quadratic model, the truncated spectral initializer returns z(0) close to the target direction with probability exceeding 1 − exp(−Ω(m)).
The guarantee requires m > c0n for a constant c0 > 0.
Concentration for isotropic sub-Gaussian rows controls the empirical matrices and their expectation when m/n exceeds a sufficiently large constant.
Davis–Kahan perturbation bounds convert matrix approximation into an eigenvector accuracy guarantee for the leading eigenvector of Y.
The Poisson observation model satisfies the required concentration condition with high probability through Chernoff bounds and a union bound.

D Local error contraction with backtracking line search

The section sketches a proof that backtracking line search yields local error contraction, with the noiseless argument extending to noisy systems. Under stated parameter conditions, the admissible step size is bounded above by 0.384, while the proof does not optimize constants.

Proof strategy: The proof establishes local error contraction for a backtracking line search, presenting the noiseless case while stating that the argument extends to noisy systems.The section explicitly treats the proof as a concise sketch and does not seek an optimized constant.
Parameter conditions: αh ≥ 6 is sufficient for Proposition 1 to remain valid, according to Proposition 4.The proposition passage also includes an additional condition involving αub, but that expression is truncated in the supplied text.
Step-size control: 0.384 bounds the backtracking step size from above when δ and ˜ϵ are sufficiently small.The bound follows from τ ≤ τub ≤ 0.384.
Parameter conditions: The step-size analysis assumes z ≥ 5, αh ≥ 6, and αp ≥ 5, together with a small-δ regime.The supplied proof also introduces an assumption and sets δ through auxiliary quantities h, b, and related identities.
Step-size control: A complementary minorization argument, omitted for conciseness, indicates that the step size is bounded away from zero.The statement combines the minorization with the stopping criterion but does not provide the resulting lower bound.