Sparse Phase Retrieval via Truncated Amplitude Flow

Introduction

In many fields of engineering and applied physics, one is often tasked with reconstructing a signal from the (squared) modulus of its Fourier (or any linear) transform, which is also known as phase retrieval (PR). Such a task arises naturally in applications such as X-ray crystallography, microscopy and ptychography, astronomy, optics, as well as array and coherent diffraction imaging. In these settings, optical sensors and detectors such as charge-coupled device cameras, photosensitive films, and human eyes record only the intensity (squared magnitude) of a light wave, but not the phase. In particular, solution to PR has led to significant accomplishments, including the discovery in $1953$ of DNA double helical structure from diffraction patterns, and the characterization of aberrations in the Hubble Space Telescope from measured point spread functions . Due to the absence of Fourier phase information, the one-dimensional ($1$D) Fourier PR problem is generally ill-posed. It can be shown that there are in fact exponentially many non-equivalent solutions beyond trivial ambiguities in the $1$D PR case . A common approach to overcome this ill-posedness is exploiting additional information on the unknown signal such as non-negativity, sparsity, or bounded magnitude . Other viable solutions consist of introducing redundancy into the measurement transforming system to obtain over-sampled and short-time Fourier transform (STFT) measurements , random Gaussian measurements , and coded diffraction patterns using structured illumination and random masks , just to name a few; see for contemporary reviews on the theory and practice of PR.

Past PR approaches can be mainly categorized as convex and nonconvex ones. A popular class of nonconvex approaches is based on alternating projections including the seminal works by Gerchberg-Saxton and Fienup , , , alternating minimization with re-sampling (AltMinPhase) , (stochastic) truncated amplitude flow (TAF) and the Wirtinger flow (WF) variants , trust-region , (stochastic) proximal linear algorithms . See also related discussion in . Specifically, the WF variants and the trust-region methods minimize the intensity (modulus squared) based empirical risk, while AltMinPhase and TAF cope with the amplitude-based empirical risk. The convex alternatives either rely on the so-called Shor’s relaxation to obtain semidefinite programming (SDP) based solvers abbreviated as PhaseLift and PhaseCut , or solve a basis pursuit problem in the dual domain as in PhaseMax .

Building on TWF and TAF, we propose here a novel sparse PR algorithm, which we call SPARse Truncated Amplitude flow (SPARTA). Adopting an amplitude-based nonconvex formulation of the sparse PR, SPARTA emerges as a two-stage iterative solver: In stage one, the support of the underlying signal is estimated first using a well-justified rule, and subsequently power iterations are employed to obtain an initialization restricted on the recovered support; while the second stage successively refines the initialization with a series of hard thresholding based truncated gradient iterations. Both stages are conceptually simple, scalable, and fast. Moreover, we demonstrate that SPARTA recovers any $k$-sparse $n$-dimensional real-/complex-valued signal $\bm{x}$ ($k\ll n$) with minimum nonzero entries (in modulus) on the order of $(1/\sqrt{k})\|\bm{x}\|_{2}$ from $\mathcal{O}(k^{2}\log n)$ measurements. Further, to reach any given solution accuracy $\epsilon>0$, SPARTA incurs total computational cost of $\mathcal{O}(k^{2}n\log n\log(1/\epsilon))$, which improves upon the state-of-the-art by at least a factor of $k$. This computational advantage is paramount in large-scale imaging applications, where the basis factor $n\log n$ is large, typically on the order of millions. In addition, SPARTA can be shown robust to additive noise of bounded support. Extensive simulated tests demonstrate markedly improved exact recovery performance (in the absence of noise), robustness to noise, and runtime speedups relative to the state-of-the-art algorithms.

The remainder of this paper is organized as follows. Section 2 reviews the sparse PR problem, and also presents known necessary and sufficient conditions for uniqueness. Section 3 details the two stages of the proposed algorithm, whose analytic performance analysis is the subject of Section 4. Finally, numerical tests are reported in Section 5, proof details are given in Section 6, and conclusions are drawn in Section 7. Supporting lemmas are presented in the Appendix.

Sparse Phase Retrieval

where $\{\psi_{i}\}_{i=1}^{m}$ are the observed modulus data, and $\{\bm{a}_{i}\}_{i=1}^{m}$ are known sensing (feature) vectors. The sparsity level $k\ll n$ is assumed known a priori for theoretical analysis purposes, while numerical implementations with unknown $k$ values will be tested as well. Alternatively, the data can be given in modulus squared (i.e., intensity) form as $\{y_{i}=|\langle\bm{a}_{i},\bm{x}\rangle|^{2}\}_{i=1}^{m}$. It has been established that $m=2k$ generic It is not within the scope of this paper to explain the meaning of generic vectors. Interested readers are referred to . (e.g., random Gaussian) measurements as in (1) are necessary and sufficient for uniquely determining a $k$-sparse solution in the real case, and $m\geq 4k-2$ are sufficient in the complex case . In the noisy scenario, stable compressive PR requires at least as many measurements as the corresponding compressive sensing problem since one is tasked with even less (no phase) information. Hence, stable sparse PR requires at least $\mathcal{O}(k\log(n/k))$ measurements as in compressive sensing . Indeed, it has been recently demonstrated that $\mathcal{O}(k\log(n/k))$ generic measurements also suffice for stable PR of a real-valued sparse signal .

Adopting the least-squares criterion (which coincides with the maximum likelihood one when assuming additive white Gaussian noise in (1)), the problem of recovering a $k$-sparse solution from phaseless quadratic equations naturally boils down to that of minimizing the ensuing amplitude-based empirical loss function

Broadening the TAF approach and the sparse PR solver TWF, the present paper puts forth a novel iterative solver for (2) that proceeds in two stages: S1) a sparse orthogonality-promoting initialization is obtained by solving a PCA-type problem with a few simple power iterations on an estimated support of the underlying sparse signal; and, S2) successive refinements of the initialization are effected by means of a series of truncated gradient iterations along with a hard thresholding per iteration to set all entries to zero, except for the $k$ ones of largest magnitudes. The two stages are presented in order next.

Algorithm: Sparse Truncated Amplitude Flow

Hereafter, assume $\bm{x}$ to be the fixed solution to problem (1) with $\phi(\bm{z})=0$; otherwise, one can replace $\bm{z}$ by $\bm{z}\text{e}^{i\phi}$, but the constant phase shift shall be dropped for notational brevity. Assume also without loss of generality that $\|\bm{x}\|_{2}=1$, which will be justified and generalized shortly.

When no sparsity is exploited, the orthogonality-promoting initialization proposed in starts with a popular folklore in stochastic geometry: High-dimensional random vectors are almost always nearly orthogonal to each other . The key idea is approximating the unknown $\bm{x}$ by another vector that is most orthogonal to a carefully chosen subset of sensing vectors $\{\bm{a}_{i}\}_{i\in\mathcal{I}^{0}}$, where $\mathcal{I}^{0}\subseteq[m]:=\{1,\,2,\,\ldots,\,m\}$ is some index set to be designed next. It is well known that the orthogonality between two vectors can be interpreted by their squared normalized inner-product $(\bm{a}_{i}^{\mathcal{T}}\bm{x})^{2}/(\|\bm{a}_{i}\|_{2}^{2}\|\bm{x}\|_{2}^{2})$. Intuitively, the smaller the squared normalized inner-product between two vectors $\bm{a}_{i}$ and $\bm{x}$ is, the more orthogonal they are to each other. Upon evaluating the inner-product between each $\bm{a}_{i}$ and $\bm{x}$ for all pairs $\{(\bm{a}_{i},\bm{x})\}_{i=1}^{m}$, one can construct $\mathcal{I}^{0}$ to include the indices of $\{\bm{a}_{i}\}$’s corresponding to the $|\mathcal{I}^{0}|$-smallest squared normalized inner-products with $\bm{x}$. Therefore, it is natural to approximate $\bm{x}$ by computing a vector $\bm{z}^{0}$ most orthogonal to the set $\mathcal{I}^{0}$ of sensing vectors . Mathematically, this is equivalent to solving a smallest eigenvector (defined to be the eigenvector associated with the smallest eigenvalue of a symmetric positive definite matrix) problem

The smallest eigenvalue (eigenvector) problem can be solved by fully eigen-decomposing the matrix $\frac{1}{|\mathcal{I}^{0}|}\sum_{i\in\mathcal{I}^{0}}\frac{\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}}{\|\bm{a}_{i}\|_{2}^{2}}$ at computational complexity $\mathcal{O}(n^{3})$ (assuming $|{\mathcal{I}}^{0}|$ to be on the order of $n$). Upon defining $\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}^{0}$ to be the complement of the set $\mathcal{I}^{0}$ in $[m]$, one can rewrite $\sum_{i\in\mathcal{I}^{0}}\frac{\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}}{\|\bm{a}_{i}\|_{2}^{2}}=\sum_{i\in[m]}\frac{\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}}{\|\bm{a}_{i}\|_{2}^{2}}-\sum_{i\in\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}^{0}}\frac{\bm{a}_{i}\bm{a}_{i}^{\mathcal{T}}}{\|\bm{a}_{i}\|_{2}^{2}}.$ Recall that for i.i.d. standard Gaussian sensing vectors $\{\bm{a}_{i}\sim\mathcal{N}(\bm{0},\bm{I}_{n})\}_{i=1}^{m}$, the following concentration result holds

The problem at hand is NP-hard in general due to the combinatorial constraint. Additionally, it can not be readily converted to a (sparse) PCA problem since the number of data samples available is much smaller than the signal dimension $n$, thus hardly validating the non-asymptotic result in (5). Although at much higher computational complexity than power iterations, semidefinite relaxation could be applied . Instead of coping with (7) directly, we shall take another route and develop our sparse orthogonality-promoting initialization approach to obtain a meaningful sparse initialization from the given limited number of measurements.

which will be shown to recover $\mathcal{S}$ exactly with high probability provided that $\mathcal{O}(k^{2}\log n)$ measurements are taken and the minimum nonzero entry $x_{\min}:=\min_{j\in S}|x_{j}|$ is on the order of $(1/\sqrt{k})\|\bm{x}\|_{2}$. The latter is postulated to guarantee such a separation between quantities having their indices belonging or not belonging to the support set. It is worth stressing that $k^{2}\log n\ll n$ when $k\ll n$, hence largely reducing the sampling size and also the computational complexity.

1.2 Orthogonality-promoting intialization

2 Thresholded Truncated Gradient Stage

Upon obtaining a sparse orthogonality-promoting initialization $\bm{z}^{0}$, our approach to solving (2) boils down to iteratively refining $\bm{z}^{0}$ by means of a series of $k$-sparse hard thresholding based truncated gradient iterations, namely,

for some $\gamma>0$ to be determined shortly, where $\{|\bm{a}_{i}^{\mathcal{T}}\bm{x}|=\psi_{i}\}$ are the given modulus data.

Main Results

The proposed sparse phase retrieval solver is summarized in Algorithm 1 along with default parameter values. Given data samples $\{(\bm{a}_{i};\psi_{i})\}_{i=1}^{m}$ generated from i.i.d. $\{\bm{a}_{i}\}_{i=1}^{m}\sim\mathcal{N}(\bm{0},\bm{I}_{n})$ sensing vectors, the following result establishes the statistical convergence rate for the proposed SPARTA algorithm in the case of $\gamma=+\infty$.

which holds with probability exceeding $1-c_{1}m{\rm e}^{-c_{0}k}-7/m$ provided that $m\geq C_{2}|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}^{0}|\geq C_{0}k^{2}\log(mn)$. Here, $c_{0},\,c_{1},\,C_{0},\,C_{2}$, and $0<\nu<1$ are some numerical constants.

Proof of Theorem 1 is deferred to Section 6 with supporting lemmas presented in the Appendix. We typically take parameters $|\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mathcal{I}}^{0}|=\lceil\frac{1}{6}m\rceil$, and $\mu=1$, which will also be validated by our analytical results on the feasible region of the step size. The constant $C_{0}$ depends on $C_{1}$, $\nu$ on $\mu$ and $C_{1}$, and $\underline{\mu}$ and $\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mu}$ rely on both $C_{1}$ and $C_{0}$. In the case of PR of unstructured signals, existing algorithms such as TAF ensures exact recovery when the number of measurements $m$ is about the number of unknowns $n$, i.e., $m\gtrsim n$. Hence, it would be more meaningful to study the sample complexity bound for PR of sparse signals when $m\lesssim n$. To this end, the sample complexity bound $m\geq C_{0}k^{2}\log(mn)$ in Theorem 1 can often be rewritten as $m\geq C_{0}^{\prime}k^{2}\log n$ for some constant $C_{0}^{\prime}>C_{0}$ and large enough $n$. Regarding Theorem 1, three observations are in order.

SPARTA recovers exactly any $k$-sparse signal $\bm{x}$ of minimum nonzero entries on the order of $(1/\sqrt{k})\|\bm{x}\|_{2}$ when there are about $k^{2}\log n$ magnitude-only measurements, which coincides with the number of measurements required by the state-of-the-art algorithms such as CPRL , sparse AltMinPhase , and TWF .

SPARTA converges at a linear rate to the globally optimal solution $\bm{x}$ with convergence rate independent of the signal dimension $n$. In other words, for any given solution accuracy $\epsilon>0$, after running at most $T=\log(1/\epsilon)$ SPARTA iterations (11), the returned estimate $\bm{z}^{T}$ is at most $\epsilon\|\bm{x}\|_{2}$ away from the global solution $\bm{x}$.

SPARTA enjoys a low computational complexity of $\mathcal{O}(k^{2}n\log n)$, and incurs a total runtime of $\mathcal{O}(k^{2}n\log n\log(1/\epsilon))$ to produce an $\epsilon$-accurate solution. The runtime is proportional to the time $\mathcal{O}(k^{2}n\log n)$ taken to read the data $\{(\bm{a}_{i},\psi_{i})\}_{i=1}^{m}$. To see this, recall that the support recovery incurs computational complexity $\mathcal{O}(k^{2}n\log n+n\log n)$, power iterations incur complexity $\mathcal{O}(k^{2}n\log n)$, and thresholded truncated gradient iterations have complexity $\mathcal{O}(k^{2}n\log n)$; hence, leading to a total complexity on the order of $k^{2}n\log n$. Given the linear convergence rate, SPARTA takes a total runtime of $\mathcal{O}(k^{2}n\log n\log(1/\epsilon))$ to achieve any fixed solution accuracy $\epsilon>0$.

Besides exact recovery guarantees in the case of noiseless measurements, it is worth mentioning that SPARTA exhibits robustness to additive noise, especially when the noise has bounded values. Numerical results using SPARTA for noisy sparse PR will be presented in the ensuing section.

Numerical Experiments

The first experiment evaluates the exact recovery performance of various approaches in terms of the empirical success rate over $100$ independent Monte Carlo trials, where the true signals are real-valued. A success is declared for a trial provided that the returned estimate incurs a relative mean-square error defined as

less than $10^{-5}$. We fixed the signal dimension to $n=1,000$, and the sparsity level at $k=10$, while the number of measurements $m/n$ increases from $0.1$ to $3$ by $0.1$. Curves in Fig. 1 clearly demonstrate markedly improved performance of SPARTA over state-of-the-art alternatives. Even when the exact number of nonzero elements in $\bm{x}$, namely, $k$ is unknown, setting $k$ in Algorithm 1 as an upper limit on the theoretically affordable sparsity level (e.g., $\lceil\sqrt{n}\;\rceil$ when $m$ is about $n$ according to Theorem 1) works well too (see the magenta curve, denoted SPARTA0). Comparison between TAF and SPARTA shows the advantage of exploiting sparsity in sparse PR settings.

The second experiment examines how SPARTA recovers real-valued signals of various sparsity levels given a fixed number of measurements. Figure 2 depicts the empirical success rate versus the sparsity level $k$, where $k$ equals the exact number of nonzero entries in $\bm{x}$. The results suggest that with a total of $m=n$ phaseless quadratic equations, TAF representing the state-of-the-art for PR of unstructured signals fails, as shown by the blue curve. Although TWF works in some cases, SPARTA significantly outperforms TWF, and it ensures exact recovery of sparse signals with up to about $25<\sqrt{n}\approx 32$ nonzero entries (due to existence of polylog factors in the sample complexity), hence justifying our analytical results.

The next experiment validates the robustness of SPARTA against additive noise present in the data. Postulating the noisy Gaussian data model $\psi_{i}=|\bm{a}_{i}^{\mathcal{T}}\bm{x}|+\eta_{i}$ , we generated i.i.d. Gaussian noise according to $\eta_{i}\sim\mathcal{N}(0,\,0.1^{2})$, $i=1,\,\ldots,\,m$. From Fig. 1, it is clear that to achieve exact recovery, SPARTA requires about $m=6k^{2}=600$ measurements, TAF about $3n=3,000$ measurements, and TWF much more than $3,000$. In this case, parameters were taken as $n=1,000$, $m=3,000$, and $k=10$, with the number of measurements large enough to guarantee that TWF and TAF also work. It is worth mentioning that SPARTA can work with a far smaller number of measurements than $m=3,000$. As seen from the plots, SPARTA performs only a few gradient iterations to achieve the most accurate solution among the three approaches, while its competing TAF and TWF require nearly an order more number of iterations to converge to less accurate estimates.

Regarding computation times, SPARTA converges much faster (both in time and in the number of iterations required to achieve certain solution accuracy) than TWF and TAF in all reported experiments. All numerical experiments were implemented with MATLAB R$2016$a on an Intel CPU @ $3.4$ GHz ($32$ GB RAM) computer.

Proof of Theorem 1

The proof of Theorem 1 will be provided in this section. To that end, we will first evaluate the performance of our sparse orthogonality-promoting initialization. The following result demonstrates that if the number of measurements is sufficiently large (on the order of $k^{2}$ within polylog factors), Step 3 of the SPARTA algorithm 1 reconstructs the support of $\bm{x}$ exactly with high probability.

Appealing to Lemma 4, one establishes for all $j\in\mathcal{S}$ that

Taking $\epsilon=x_{\min}^{2}:=\min_{j\in\mathcal{S}}x_{j}^{2}\leq x_{j}^{2}$ leads to

Recalling our assumption that $x_{\min}^{2}$ is on the order of $(1/k)\|\bm{x}\|_{2}^{2}$, i.e., $x_{\min}^{2}=(C_{1}/k)\|\bm{x}\|_{2}^{2}$ for certain constant $C_{1}>0$, the following holds with probability at least $1-1/m$ for all $j\in\mathcal{S}$

provided that $m\geq C_{0}k^{2}\log(mn)$ for some absolute constant $C_{0}>0$.

Now let us turn to the case of $j\notin\mathcal{S}$, in which $\sum_{i=1}^{m}Z_{i,j}=\sum_{i=1}^{m}\psi_{i}^{2}a_{i,j}^{2}$ is a weighted sum of $\chi_{1}^{2}$ random variables. According to Lemma 5, it holds that

In addition, for any constants $\epsilon^{\prime},\,\epsilon^{\prime\prime}>0$, Chebyshev’s inequality together with the union bound confirms that

Take $\epsilon:=\log(mn)$ in (20), $\epsilon^{\prime}:=\sqrt{m}$ and $\epsilon^{\prime\prime}:=4\log(mn)$ in (21). Then, with probability at least $1-4/m$, the next holds for all $j\notin\mathcal{S}$ and $m>C^{\prime}$

for some absolute constant $C^{\prime}>0$ depending on $n$.

On the other hand, the rotational invariance property of Gaussian distributions asserts that $\psi_{i}^{2}=|\bm{a}_{i}^{\mathcal{T}}\bm{x}|^{2}=|\bm{a}_{i,\mathcal{S}}^{\mathcal{T}}\bm{x}_{\mathcal{S}}|^{2}\overset{d}{=}a_{i,j}^{2}\|\bm{x}\|_{2}^{2}$ , in which the symbol $\overset{d}{=}$ means that terms involved on both sides of the equality enjoy the same distribution. Since the $\chi^{2}$ variables $a_{i,j}^{2}$ are sub-exponential, an application of Bernstein’s inequality produces the tail bound

for any $\epsilon\in(0,1)$, which can also be easily verified with a direct tail probability calculation from the tail probability of standard Gaussian distribution. Choosing $\epsilon:=\sqrt{16\log(m)/m}$ with $m>C^{\prime}$ gives rise to

which holds true with probability at least $1-{1}/{m}$ for all $j\in[m]$. Putting results in (6) and (24) together leads to

which holds with probability exceeding $1-5/m$ for large enough $m$.

The last inequality taken collectively with (19) suggests that there exists an event $E_{0}$ on which with probability at least $1-6/m$, the following holds

provided that $m\geq C_{0}k^{2}\log(mn)$ such that $C_{0}\geq 144/C_{1}^{2}$ with $x_{\min}^{2}=(C_{1}/k)\|\bm{x}\|_{2}^{2}$. ∎

provided that $m\geq C_{0}^{\prime}k$ for some absolute constant $C_{0}^{\prime}>0$.

The proof can be directly adapted from [9, Proposition 1], and hence it is omitted.

Take a constant learning parameter $\mu\in(\underline{\mu},\,\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mu})$. There exists an event of probability at least $1-c_{1}m^{-c_{0}k}$, such that on this event, starting from an initial estimate $\bm{z}^{0}$ satisfying ${\rm dist}(\bm{z}^{0},\bm{x})\leq(1/10)\|\bm{x}\|_{2}$, successive estimates by Step 5 with $\gamma=+\infty$ in Algorithm 1 obey

if $m\geq C_{0}^{\prime\prime}(3k)\log(n/(3k))$. Here, $\underline{\mu},\,\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\mu}_{0},\,c_{0},\,c_{1},\,C_{0}^{\prime\prime}>0$ are certain universal constants.

It is worth noting that Step 5 of Algorithm 1 guarantees linear convergence to the globally optimal solution $\bm{x}$ as long as the initial guess $\bm{z}^{0}$ lands within a small neighborhood of $\bm{x}$, regardless of whether $\bm{z}^{0}$ estimates exactly the support of $\bm{x}$ or not.

To start, let us establish a bit of notation, which will be used only in this section. Define for all $t\geq 0$

The proof of Lemma 3 will be mainly based on results in , and , . The former helps establishing the so-termed local regularity condition that will be key to proving linear convergence of iterative optimization algorithms to the globally optimal solutions of nonconvex optimization problems , while the latter two offer a standard approach to dealing with the nonlinear hard thresholding operator involved in our proposed SPARTA algorithm. Specifically, based on the triangle inequality of the vector $2$-norm, one arrives at

where in the last inequality the first term denotes the distance of $\bm{x}_{\Theta^{t+1}}$ to the estimate $\bm{d}^{t+1}_{\Theta^{t+1}}$ before hard thresholding, and the second denotes the distance between $\bm{d}^{t+1}_{\Theta^{t+1}}$ and its best $k$-approximation $\bm{z}^{t+1}_{\Theta^{t+1}}$ because $\bm{z}_{\Theta^{t+1}}^{t+1}$ has cardinality equal to $k$. The optimality of $\bm{z}^{t+1}_{\Theta^{t+1}}$ implies $\|\bm{z}^{t+1}_{\Theta^{t+1}}-\bm{d}_{\Theta^{t+1}}^{t+1}\|_{2}\leq\|\bm{x}_{\Theta^{t+1}}-\bm{d}_{\Theta^{t+1}}^{t+1}\|_{2}$. Plugging the latter inequality back into (29) yields

Define the estimation error $\bm{h}^{t}:=\bm{x}-\bm{z}^{t}$. Rewriting and substituting

$\|\bm{A}_{\mathcal{B}}^{\mathcal{T}}\bm{A}_{\mathcal{D}}\|_{2}\leq\delta_{3k}$, where $\mathcal{B}$ and $\mathcal{D}$ are disjoint sets of combined cardinality not exceeding $3k$;

$\|\bm{A}_{\mathcal{B}\cup\mathcal{D}}^{\mathcal{T}}\bm{A}_{\mathcal{B}\cup\mathcal{D}}-\bm{I}\|_{2}\leq\delta_{3k}$.

Having elaborated on the properties of RIP matrices, we are ready to derive bounds for the three terms on the right hand side of (31). Regarding the first term, it is easy to check that

where $\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\lambda},\,\underline{\lambda}>0$ are the largest and smallest eigenvalue of $(1/m)\sum_{i=1}^{m}\bm{a}_{i,\Theta^{t+1}}\bm{a}_{i,\Theta^{t+1}}^{\mathcal{T}}$, respectively. Specifically, the two inequalities in (33) are obtained based on the definition of the induced $2$-norm (i.e., the spectral norm) of matrices.

Next, we estimate the eigenvalues $\accentset{{\cc@style\underline{\mskip 9.5mu}}}{\lambda}$ and $\underline{\lambda}$. Using P2, it clearly holds that

due to $|\Theta^{t+1}|\leq 2k$. For the same reason, it further holds that

Taking the results in (34) and (35) into (33) yields

For the second term in (31), since $|\Theta^{t+1}\cup\Theta^{t}|\leq 3k$, the next holds with high probability

in which the first inequality arises again from the definition of the matrix $2$-norm. The last inequality can be obtained by appealing to P4.

Consider now the last term in (31). For convenience, define $\bm{A}_{\Theta^{t+1}}^{\mathcal{T}}:=[\bm{a}_{1,\Theta^{t+1}}~{}\cdots~{}\bm{a}_{m,\Theta^{t+1}}]$ with $|\Theta^{t+1}|\leq 2k$, and also $\bm{v}^{t}:=[v_{1}^{t}~{}\cdots~{}v_{m}^{t}]^{\mathcal{T}}$ with $v_{i}^{t}:=(\frac{\bm{a}_{i}^{\mathcal{T}}\bm{z}^{t}}{|\bm{a}_{i}^{\mathcal{T}}\bm{z}^{t}|}-\frac{\bm{a}_{i}^{\mathcal{T}}\bm{x}}{|\bm{a}_{i}^{\mathcal{T}}\bm{x}|})|\bm{a}_{i}^{\mathcal{T}}\bm{x}|$ for $i=1,\,\ldots,\,m$. Upon rearranging terms, the induced matrix $2$-norm definition implies that

Plugging the inequality in (6) into (39) leads to

Substituting the three bounds in (6), (6), and (6) into (31), we obtain

over disjoint sets $\Theta^{t+1}$ and $\Theta^{t}\setminus\Theta^{t+1}$. To ensure linear convergence, it suffices to choose a constant step size $\mu>0$ such that

For sufficiently small $\delta_{3k}>0$ and $\epsilon_{1}>0$, one has $\nu:=1-\rho\in(0,1)$, which justifies the linear convergence result in (14). ∎

Theorem 1 can be directly implied by combining Lemmas 1, 2, and 3. In fact, Lemma 1 ensures exact support recovery so that the orthogonality-promoting initialization can be effectively performed on the equivalent dimension-reduced data samples. Lemma 2 guarantees that the sparse initialization attained based on the dimensional-reduced data lands within a small neighborhood of the globally optimal solution (this region is also termed basin of attraction; see e.g., , , for more details) with high probability. Starting from any point within the basin of attraction, Lemma 3 confirms that successive iterates of SPARTA will be dragged toward the globally optimal solution at a linear rate provided that the step size and the truncation threshold are appropriately selected.

Concluding Remarks

This paper contributed a sparse truncated amplitude flow (SPARTA) algorithm for solving PR of sparse signals. SPARTA initially recovers the support of the underlying sparse signal, which is used to obtain a sparse orthogonality-promoting initialization using power iterations restricted on the estimated support; subsequently, SPARTA refines the initialization by means of hard thresholding based truncated gradient iterations to ensure overall simplicity and scalability. SPARTA enjoys provably exact recovery as soon as the number of noiseless Gaussian measurements exceeds a certain bound. In contrast to state-of-the-art algorithms, such as AltMinPhase and TWF, SPARTA requires the same sample size but can afford lower computational complexity. Simulated tests corroborate markedly improved recovery performance and computational efficiency of SPARTA relative to existing alternatives.

Appendix: Supporting Lemmas

for $\sigma^{2}:=m\max(b^{2},\,v^{2})$, and the cumulative distribution function of the standard normal distribution $\Phi(\cdot)$, where one can take $c_{0}=25$.

Let $X_{1},\,X_{2},\,\ldots,\,X_{m}$ be i.i.d. Gaussian random variables with zero mean and variance $1$, and $b_{1},\,b_{2},\,\ldots,b_{m}$ be nonnegative. The following inequality holds for any $\epsilon>0$

The proof of Lemma 6 can be found in [56, Page 30], which generalizes the result of [19, Lemma 3].

Acknowledgment

The authors would like to thank the anonymous reviewers for their thorough review and all constructive comments and suggestions, which helped to improve the quality of the manuscript. The authors also thank Prof. Xiaodong Li for sharing the codes of the thresholded Wirtinger flow algorithm.

References