Optimal Rates of Convergence for Noisy Sparse Phase Retrieval via Thresholded Wirtinger Flow

Introduction

In a range of fields in science and engineering, researchers face the problem of recovering a $p$-dimensional signal of interest $\bm{x}$ by probing the signal via a set of $p$-dimensional sensing vectors $\bm{a}_{j}$, $j=1,\dots,m$, and hence the observations are the $(\bm{a}_{j}^{\prime}\bm{x})$’s contaminated with noise. This gives rise to the linear regression model in statistical terminology where $\bm{x}$ is the regression coefficient vector and $\bm{A}=[\bm{a}_{1},\dots,\bm{a}_{m}]^{\prime}$ is the design matrix. There is an extensive literature on the theory and methods for the estimation/recovery of $\bm{x}$ under such a linear model. However, in many important applications, including X-ray crystallography, microscopy, astronomy, diffraction and array imaging, interferometry, and quantum information, it is sometimes impossible to observe $\bm{a}_{j}^{\prime}\bm{x}$ directly and the measurements that one is able to obtain are the magnitude/energy of $\bm{a}_{j}^{\prime}\bm{x}$ contaminated with noise. In other words, the observations are generated by the following phase retrieval model:

Efficient computational methods for phase retrieval have been proposed in the community of optics, and they are mostly based on the seminal work by Gerchberg, Saxton, and Fienup . The effectiveness of these methods relies on careful exploration of prior information of the signal in the spatial domain. Moreover, these methods were revealed later as non-convex successive projection algorithms . This provides insight for occasional observation of stagnation of iterates and failure of convergence.

Recently, inspired by multiple illumination, novel computational methods were proposed for phase retrieval without exploring and employing a priori information of the signal. These methods include semidefinite programming , polarization , alternating minimization , gradient methods , alternating projection , etc. More importantly, profound and remarkable theoretical guarantees for these methods have also been established. As for noiseless sparse phase retrieval, semidefinite programming has been proven to be effective with theoretical guarantees . Other empirical methods for sparse phase retrieval include belief propagation and greedy methods .

Regarding noisy phase retrieval, some stability results have been established in the literature; See . In particular, stability results have been established in for noisy sparse phase retrieval by semidefinite programming, though the authors did not study the optimal dependence of the convergence rates on the sparsity of the signal and the sample size. Nearly minimax convergence rates for sparse phase retrieval with Gaussian noise have been established in under sub-gaussian design matrices. However, the optimal rates are achieved by empirical risk minimization under sparsity constraints, in which both the objective function and the constraint are non-convex, implying that the procedure is not computationally feasible.

Methodology

The major component of the our method is a thresholded gradient descent algorithm to obtain a sparse solution to a given non-convex empirical risk minimization problem. Due to the non-convex nature of the problem, in order to avoid any local optimum that is far away from the truth, the initialization step is crucial. Thus, we also provide a candidate method which can be justified theoretically for yielding a good initializer. The methodology is proposed assuming that $\bm{A}$ has standard Gaussian entries, though it could potentially also be used when such an assumption does not necessarily hold.

Given the sensing vectors $\bm{a}_{j}$ and the noisy magnitude measurements $y_{j}$ as in (1.1) for $j=1,\dots,m$, one can consider estimating $\bm{x}$ by minimizing the following empirical risk function

Statistically speaking, in the low-dimensional setup with fixed $p$ and $m\rightarrow\infty$, if the additive noises are heavy-tailed, least-absolute-deviations (LAD) methods might be more robust than least-squares methods. However, recent progress in modern linear regression analysis shows that least-squares could be preferable to LAD when $p$ and $m$ are proportional, even the noises are sub-exponential . Due to this surprising phenomenon, we simply take the least-squares empirical risk in (2.1), although phase retrieval is a nonlinear regression problem, which could be very different from linear regression. More importantly, close-form gradient methods can be induced from the empirical risk function in (2.1), which is computationally convenient. To be specific, at any current value of $\bm{z}$, one updates the estimator by taking a step along the gradient direction

until a stationary point is reached. Indeed, Candès et al. showed that under appropriate conditions, initialized by an appropriate spectral method, a gradient method, referred to as Wirtinger flow, leads to accurate recovery of $\bm{x}$ up to a global phase in the complex domain and noiseless setting.

However, the direct application of gradient descent is not ideal for noisy sparse phase retrieval since it does not utilize the knowledge that the true signal $\bm{x}$ is sparse in order to mitigate the contamination of the noise. To incorporate this a priori knowledge, it makes sense to seek a “sparse minimizer” of (2.1). To this end, suppose we have a sparse initial guess $\bm{x}^{(0)}$ for $\bm{x}$. To update $\bm{x}^{(0)}$ to another sparse vector, we may take a step along $\nabla f(\bm{x}^{(0)})$, and then sparsify the result by thresholding.

Indeed, if we were given the oracle knowledge of the support $S$ of $\bm{x}$, then we can reduce the problem to recovering $\bm{x}_{S}$ based on the $\{y_{j},a_{jS}\}_{j=1}^{m}$. By avoiding estimating any coordinate of $\bm{x}$ in $S^{c}$, we could greatly reduce variance of the resulting estimator of $\bm{x}$. In reality, we do not have such oracle knowledge and the additional thresholding step added on top of gradient descent is intended to mimic the oracle behavior by hopefully restricting all the updated coordinates on $S$.

Let $\mathcal{T}_{\tau}$ be any thresholding function satisfying

For any vector $\bm{b}=(b_{1},\dots,b_{p})^{\prime}$, let $\mathcal{T}_{\tau}(\bm{b})=(\mathcal{T}_{\tau}(b_{1}),\dots,\mathcal{T}_{\tau}(b_{p}))^{\prime}$. With the foregoing definition, the proposed thresholded gradient descent method can be summarized as Algorithm 1. In view of the Wirtinger flow method for noiseless phase retrieval , we name our approach the “Thresholded Wirtinger Flow” method. The data-driven choice of the threshold level in (2.3) is motivated by the following intuition. Recall that we assume the sensing vectors $\{\bm{a}_{j}:j=1,\ldots,m\}$ are independent standard Gaussian vectors. For a fixed $\bm{z}$, if we act as if each $(|\bm{a}_{j}^{\prime}\bm{z}|^{2}-y_{j})(\bm{a}_{j}^{\prime}\bm{z})$ is a fixed constant, then the gradient in (2.2) is a linear combination of Gaussian vectors and hence has i.i.d. Gaussian entries with mean zero and variance $\frac{1}{m^{2}}\sum_{j=1}^{m}(|\bm{a}_{j}^{\prime}\bm{z}|^{2}-y_{j})^{2}(\bm{a}_{j}^{\prime}\bm{z})^{2}$. Therefore, the threshold $\tau(\bm{z})$ is simply $\sqrt{\beta\log(mp)}$ times the standard deviation of these Gaussian random variables, which is essentially the universal thresholding in the Gaussian sequence model literature . Although the above intuition is not exactly true, the resulting thresholds in (2.3) are indeed the right choices as justified later in Section 3, and illustrated in Section 4. Notice that there are two tuning parameters $\mu$ and $\beta$, which should be treated as absolute constants. We will validate some theoretical choices and also provide practical choices later.

2 Initialization

It is worth noting that the success of Algorithm 1 depends crucially on the initial estimator for two reasons. First, the empirical risk (2.1) is a non-convex function of $\bm{z}$ and hence it could have multiple local minimizers. Hence the success of a gradient descent based approach depends naturally on the starting point. Moreover, an accurate initializer can reduce the required number of iterations in the thresholded Wirtinger flow algorithm. In view of its crucial rule, we propose in Algorithm 2 an initialization method which can be proven to yield a decent starting point for Algorithm 1 under our modeling assumption.

The motivation of the algorithm is similar to that of diagonal thresholding for sparse PCA: we want to identify a small collection of coordinates with big marginal signals and then compute an estimator of $\bm{x}$ by focusing only on these coordinates. In particular, the quantity $I_{l}$ in (2.7) captures the marginal signal strength of the $l$-th coordinate and $\widehat{S}_{0}$ (2.8) selects all coordinates with big marginal signals. Last but not least, (2.9) and (2.10) computes the initial estimator by focusing only on the coordinates in $\widehat{S}_{0}$. There is a tuning parameter $\alpha$ needed as input of the algorithm, which can be treated as an absolute constant. We will provide some justified theoretical choice later.

Theory

Suppose $\beta=4$ in (2.3), and $\alpha=K\left(1+\frac{\sigma}{\|\bm{x}\|_{2}^{2}}\right)$ in (2.8) for some absolute constant $K$. Suppose $\mu\leq\mu_{0}$ in (2.4) and $m\geq C\left(1+\frac{\sigma^{2}}{\|\bm{x}\|_{2}^{4}}\right)k^{2}\log(mp)$. For all $t=1,2,3,\ldots$, there holds

where $\mu_{0}$, $C$, and $C_{0}$ are some absolute constants.

The proof is given in Section 6. Lemma 6.3 guarantees the efficacy of the initialization step Algorithm 2, and Lemmas 6.4 and 6.5 explain why the thresholded Wirtinger flow method leads to accurate estimation. Here $\beta=4$ and $\alpha=K\left(1+\frac{\sigma}{\|\bm{x}\|_{2}^{2}}\right)$ are chosen for analytical convenience. The discussion of empirical choices of $\beta$, $\alpha$, and $\mu$ are deferred to Section 4.

Let us interpret Theorem 3.1 by considering the following cases. In the noiseless case, with high probability, we obtain $\min\limits_{i=0,1}\|\widehat{\bm{x}}^{(t)}-(-1)^{i}\bm{x}\|_{2}\leq\frac{1}{6}\left(1-\frac{\mu}{16}\right)^{t}\|\bm{x}\|_{2}$. This implies that thresholded gradient descent method leads to linear convergence to the original signal up to a global sign.

In the noisy case, if $\mu>0$ is an absolute constant, by letting $t\asymp\log\left(1/\delta\right)$ where $\delta=\frac{\sigma}{\|\bm{x}\|_{2}^{2}}\sqrt{\frac{k\log p}{m}}$, we obtain $\min\limits_{i=0,1}\|\widehat{\bm{x}}^{(t)}-(-1)^{i}\bm{x}\|_{2}\precsim\frac{\sigma}{\|\bm{x}\|_{2}}\sqrt{\frac{k\log p}{m}}$ with high probability. If the knowledge of $\delta$ is not available, by choosing $t=O(\log p)$, we can obtain $\min\limits_{i=0,1}\|\widehat{\bm{x}}^{(t)}-(-1)^{i}\bm{x}\|_{2}\precsim\frac{\sigma}{\|\bm{x}\|_{2}}\sqrt{\frac{k\log p}{m}}+\frac{1}{p^{c}}$ for any predetermined $c>0$. The convergence rate $\frac{\sigma}{\|\bm{x}\|_{2}}\sqrt{\frac{k\log p}{m}}$ is better than the upper bound result established in , which is achieved by the intractable sparsity constrained empirical risk minimization. Our contribution is to show that this rate can be obtained tractably by a fast algorithm.

Ignoring any polylog factor, the above convenient properties of thresholded Wirtinger flow are guaranteed by the sample size condition $m\gtrsim k^{2}$. When $m\ll p$, this condition is crucial for the effectiveness of initialization Algorithm 2. An immediate question is whether such a minimum sample size condition is in some sense necessary for any computationally efficient algorithm, if the sensing matrix is random and structureless? A similar phenomenon has been previously observed in the related but different problem of sparse principal component analysis. Assuming the hardness of the planted clique problem , a series of papers have shown that a comparable minimum sample size condition is necessary for any estimator computable in polynomial time complexity to achieve consistency and optimal convergence rates uniformly over a parameter space of interest. In particular, it was shown in that this is the case even for the most restrictive parameter space in sparse principal component analysis – (discretized) Gaussian single spiked model with a sparse leading eigenvector. Establishing comparable computational lower bounds for sparse phase retrieval, especially under the Gaussian design, is an interesting project for future research.

In the case when $m\gtrsim p$ ignoring any log factor, it is well-known that a consistent initializer can be obtained by spectral methods , no matter whether $\bm{x}$ is sparse or not. In other words, the diagonal thresholding idea in Algorithm 2 is not as crucial as in the case $m\ll p$. It is interesting to investigate whether $m\gtrsim k^{2}$ can be relaxed such that the optimal converge rates can still be achieved by thresholded Wirtinger flow.

The convergence rate $\frac{\sigma}{\|\bm{x}\|_{2}}\sqrt{\frac{k\log p}{m}}$ is essentially optimal. The following lower bound result has been essentially proven in :

provided $m\geq C\left(\frac{\sigma^{2}}{\|\bm{x}\|_{2}^{4}}+1\right)k\log(ep/k)$, where both $C$ and $C_{0}$ are some absolute constants.

Notice that for a standard Gaussian variable with variance $\sigma^{2}$, its sub-exponential norm is a constant multiple of $\sigma$. For brevity, we do not scale the Gaussian noises such that their sub-exponential norms are strictly less than or equal to $\sigma$.

Numerical Simulation

In this section, we report numerical simulation results to demonstrate how the relative estimation error depends on the thresholding parameter $\beta$, the noise-to-signal ratio (NSR) $\sigma/\|\bm{x}\|_{2}^{2}$, the sample size $m$, and the sparsity $k$. To guarantee fair comparison, we always fix the length of the signal $p=1000$ and the initialization parameter $\alpha=0.1$ (except for the first case on thresholding effect). Moreover, in each numerical experiment, we conservatively choose gradient parameter $\mu=0.01$, and the number of iterations $T=1000$ for thresholded Wirtinger flow. The resulting estimator is denoted as $\widehat{\bm{x}}=\widehat{\bm{x}}^{(1000)}$. With each fixed $k$, the support of $\bm{x}$ is uniformly distributed at random. The nonzero entries of $\bm{x}$ are i.i.d. $\sim\mathcal{N}(0,1)$. The noise $\epsilon\sim\mathcal{N}(\bm{0},\sigma^{2}\bm{I}_{m})$, where $\sigma$ is determined by $\|\bm{x}\|_{2}$ and the choice of NSR $\sigma/\|\bm{x}\|_{2}^{2}$. As discussed before, the design matrix $\bm{A}$ consists of independent standard Gaussian random variables.

Thresholding effect: Fix $\alpha=0.1$, $m=7000$, $k=100$, and $\sigma/\|\bm{x}\|_{2}^{2}=1$. For each $\beta=0,0.25,0.5,\ldots,3$, we implement the algorithm for $10$ times with independently generated $\bm{A}$, $\bm{x}$, and $\bm{\epsilon}$. and then take the average of the $10$ independent relative errors $\min(\|\widehat{\bm{x}}-\bm{x}\|_{2},\|\widehat{\bm{x}}+\bm{x}\|_{2})/\|\bm{x}\|_{2}$. The relation between the average relative error and the choice of $\beta$ is plotted as the red curve in Figure 1. The result shows that the average relative error essentially decreases from $0.2365$ to $0.1151$ as the thresholding parameter increases from to $0.75$, and then increases slowly up to $0.1684$ as $\beta$ continues to increase to $3$.

We implement the above experiments again with the only difference $\alpha=0.5$. The relation curve between the relative estimation error and $\beta$ is plotted as the blue curve in Fig. 1. It is clear that the performance of the algorithm is very close to the case $\alpha=0.1$.

Noise effect: Fix $m=7000$, $k=100$, and $\beta=1$. In each choice of NSR $\sigma/\|\bm{x}\|_{2}^{2}=0,0.1,\ldots,1$, with $5$ instances of $(\bm{A},\bm{x},\bm{\epsilon})$ generated independently, we take the average of the relative error $\min(\|\widehat{\bm{x}}-\bm{x}\|_{2},\|\widehat{\bm{x}}+\bm{x}\|_{2})/\|\bm{x}\|_{2}$. In Figure 2, it shows how the average relative error depends on NSR. The average relative error strictly increases from $0.0000$ to $0.1219$, as the NSR increases from to $1$.

Sample size effect: Fix $k=100$, $\sigma/\|\bm{x}\|_{2}^{2}=1$, and $\beta=1$. In each choice of $m=2000,3000,\ldots,11000$, with $5$ instances of $(\bm{A},\bm{x},\bm{\epsilon})$ generated independently, we take the average of the relative error $\min(\|\widehat{\bm{x}}-\bm{x}\|_{2},\|\widehat{\bm{x}}+\bm{x}\|_{2})/\|\bm{x}\|_{2}$. In Figure 3, it shows how the average relative error depends on the sample size. When the sample sizes are $2000$ and $3000$, i.e., twice and three times as large as $p$, the average relative errors are $0.8444$ and $0.3651$ respectively. In these cases, the thresholded gradient descent method leads to poor recovery of the original signal. When the sample size increases from $4000$ to $11000$, the average relative error decreases steadily from $0.1692$ to $0.0956$.

Sparsity effect: Fix $m=7000$, $\sigma/\|\bm{x}\|_{2}^{2}=1$, and $\beta=1$. In each choice of sparsity $k=25,50,\ldots,200$, with $10$ instances of $(\bm{A},\bm{x},\bm{\epsilon})$ generated independently, we take the average of the relative error $\min(\|\widehat{\bm{x}}-\bm{x}\|_{2},\|\widehat{\bm{x}}+\bm{x}\|_{2})/\|\bm{x}\|_{2}$. Figure 4 demonstrates the relation between the average relative error and the sparsity. The average relative error essentially increases from $0.1059$ to $0.1666$, as the sparsity increases from $25$ to $200$.

Discussion

In this paper, we established the optimal rates of convergence for noisy sparse phase retrieval under the Gaussian design in the presence of sub-exponential noise, provided that the sample size is sufficiently large. Furthermore, a thresholded gradient descent method called “Thresholded Wirtinger Flow” was introduced and shown to achieve the optimal rates.

Iterative thresholding has been employed in a variety of problems in high-dimensional statistics, machine learning, and signal processing, under the assumption that the signal or parameter vector/matrix satisfies a sparse or low-rank constraint. Examples include compressed sensing/sparse approximation , sparse principal component analysis , high-dimensional regression , and low-rank recovery .

Regarding the application of iterative thresholding and projected gradient methods in high-dimensional $M$-estimation, their statistical optimality has been established when the empirical risk function satisfies certain properties, such as restrictive strong convexity and smoothness (RSC and RSS) . Although our thresholded gradient method aims to solve (2.1) for a sparse solution, the existing analytical framework for high-dimensional $M$-estimation does not apply to the sparse phase retrieval problem, since the empirical risk function in (2.1) does not satisfy RSC in general, no matter how large the sample size is. Instead, we have shown that thresholded gradient methods can achieve optimal statistical precision for signal recovery, even when the empirical risk function does not satisfy the common assumption of RSC.

Besides thresholded gradient methods, convexly and non-convexly regularized methods are also widely-used for high-dimensional $M$-estimation. In fact, some iterative thresholding methods are induced by regularizations; See, e.g., . Therefore, an alternative candidate method for solving the noisy sparse phase retrieval problem is to penalize the empirical risk function in (2.1) before taking the minimum, in order to promote a sparse solution. The major difficulty is apparently the non-convexity of the empirical risk function. An interesting result in guarantees the statistical precision of all local optima, as long as the non-convex penalty satisfies certain regularity conditions, and the empirical risk function, possibly non-convex, satisfies the restricted strong convexity. A similar result appeared in , in which the empirical risk function is required to satisfy a sparse eigenvalue (SE) condition. However, back to noisy sparse phase retrieval, the empirical risk function in (2.1) satisfies neither RSC nor SE in general, so there is no guarantee that all local optima are consistent. A natural question is whether some penalized version of (2.1) is strongly convex in a sufficiently large neighborhood of its global minimum, such that a tractable initializer lies in this neighborhood provided the sample size is sufficiently large. Another interesting question is whether the global minimizer of such penalized version of (2.1) is a rate-optimal estimator of the original sparse signal. We leave these questions for future research.

Proof of Theorem 3.1

For any two two random variables/vectors/matrices/sets $X$ and $Y$, we denote by X\rotatebox[origin={c}]{90.0}{\models}Y if $X$ and $Y$ are independent.

From the model (1.1), we have \bm{y}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}. Moreover, we have \{I_{1},\ldots,I_{k}\}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}} and \phi\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}, where $\phi$ and $\{I_{1},\ldots,I_{k}\}$ are defined in (2.6) and (2.7), respectively.

Proof The fact $\bm{y}=|\bm{A}\bm{x}|^{2}+\bm{\epsilon}=|\bm{A}_{S}\bm{x}_{S}|^{2}+\bm{\epsilon}$ implies straightforwardly that \bm{y}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}. By (2.7), we know for all $l=1,\ldots,k$, $I_{l}$ are defined by $\bm{y}$ and $\bm{A}_{S}$, which implies that I_{l}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}} for all $l=1,\ldots,k$. Finally, by (2.6), we know $\phi$ is determined uniquely by $\bm{y}$, which implies that \phi\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}.

On an event $\widetilde{E}_{0}$ with probability at least $1-\frac{3}{m}$,

for some numerical constant $C_{0}>0$. As a consequence, as long as $\frac{m}{\log m}\geq C(\delta)\left(1+\frac{\sigma^{2}}{\|\bm{x}\|_{2}^{4}}\right)$, there holds

Proof By the definition of $\phi^{2}$ and $y_{j},j=1,\ldots,m$, we have

As shown in Lemma A.7, with probability at least $1-\frac{1}{m}$,

for some numerical constant $C_{0}>0$. Moreover, since $\bm{x}$ is fixed, there holds

By Lemma 4.1 of , with probability at least $1-\frac{2}{m}$, we have

Let $\alpha=K\left(1+\frac{\sigma}{\|\bm{x}\|_{2}^{2}}\right)$ for some large enough absolute constant $K$, and $\widehat{\bm{x}}^{(0)}$ be defined in Algorithm 2. There exists a random vector $\bm{x}^{(0)}$ satisfying \bm{x}^{(0)}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}} and $\operatorname{supp}(\bm{x}^{(0)})\subset S$, such that on an event $E_{01}$ with probability at least $1-\frac{16}{m}-2e^{-k}$, we have

provided $m\geq C\left(1+\frac{\sigma^{2}}{\|\bm{x}\|_{2}^{4}}\right)k^{2}\log(mp)$. Here $C$ is an absolute constant.

Proof Recall that $S=\{1,\ldots,k\}$ and $I_{l}=\frac{1}{m}\sum_{j=1}^{m}y_{j}a_{jl}^{2}$ for $l=1,\ldots,p$. Define

with $2$-norm $\phi$. This easily implies $\operatorname{supp}(\bm{x}^{(0)})\subset S_{0}\subset S$. Since \{\bm{W}_{S_{0}S_{0}},\phi\}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}, we also have \bm{x}^{(0)}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}.

To simplify notation, let us write for any $j\in[m]$, $\widetilde{y}_{j}:=(\bm{a}_{j}^{\prime}\bm{x})^{2}=(\bm{a}_{jS}^{\prime}\bm{x})^{2}$, which implies $y_{j}=\widetilde{y}_{j}+\epsilon_{j}$. Notice that

in which we will first control the second term. For a given $l\in[p]$, we know $a_{1l}^{2}-1,\ldots,a_{ml}^{2}-1$ are i.i.d. centered sub-exponential random variables with sub-exponential norms being an absolute constant. Then, by Bernstein inequality (see, e.g., Proposition 16 in ), we have with probability at least $1-\frac{2}{mp}$,

for some absolute constant $C_{0}$. Then by Lemma A.7, with probability at least $1-4/m$, we have

provided $m\geq C(\log p)$ for some absolute constant $C$.

Next, we prove that with high probability $\bm{x}^{(0)}=\widehat{\bm{x}}^{(0)}$. It suffices to prove $\widehat{S}_{0}=S_{0}$, i.e., $\widehat{S}_{0}\subset S$. For any $l\in S^{c}$, $a_{jl}$ and $\widetilde{y}_{j}$ are independent, and so conditional on $\{\widetilde{y}_{j},j\in[m]\}$, $\sum_{j=1}^{m}\widetilde{y}_{j}a_{jl}^{2}$ is a weighted sum of $\chi^{2}_{1}$ variables. By Lemma 4.1 of ,

Moreover, Chebyshev’s inequality, the Gaussian tail bound and the union bound lead to

Thus, with probability at least $1-\frac{4}{m}$, for all $l\in S^{c}$,

Here the last inequality holds when $m\geq C$ for some absolute constant $C$.

Since $\alpha=K\left(1+\frac{\sigma}{\|\bm{x}\|_{2}^{2}}\right)$ with large enough $K$, by (6.3), (6.5), (6.4) and Lemma 6.2, we obtain that with probability at least $1-\frac{11}{m}$, for all $l\in S^{c}$,

which implies that $\widehat{S}_{0}\subset S$.

Next, Lemma 4.1 of leads to with probability at least $1-\frac{1}{m}$,

The last two inequalities, together with (6.4) and (6.3), imply that with probability at least $1-\frac{6}{m}$, for all $l\in S$,

Define $S_{-}=\left\{l\in S:x_{l}^{2}\geq\left(11+\frac{3}{5}\alpha\right)\|\bm{x}\|_{2}^{2}\sqrt{\log(mp)\over m}\right\}$. Then, for all $l\in S_{-}$ we have

Since $\alpha=K\left(1+\frac{\sigma}{\|\bm{x}\|_{2}^{2}}\right)$ with sufficiently large absolute constant $K$, by lemma 6.2, we have or all $l\in S_{-}$,

with probability at least $1-9/m$. This implies $S_{-}\subset S_{0}$.

By Lemma A.6, with probability at least $1-1/m$, we have

provided $m\geq C(\delta)k\log p$. Moreover, by Lemma A.7 and Lemma A.8, with probability at least $1-2/m-2e^{-k}$, we have ${\left\lVert{\sum_{j=1}^{m}\epsilon_{j}\bm{a}_{jS}\bm{a}_{jS}^{\prime}}\right\rVert}\leq C_{0}\sigma\sqrt{m(k+\log m)}$. By assuming $m\geq C(\delta)\frac{\sigma^{2}}{\|\bm{x}\|_{2}^{4}}k\log(mp)$, we have $\frac{1}{m}{\left\lVert{\sum_{j=1}^{m}\epsilon_{j}\bm{a}_{jS}\bm{a}_{jS}^{\prime}}\right\rVert}\leq\frac{\delta}{2}\|\bm{x}\|_{2}^{2}$. This implies that

By Lemma 6.2, we have $1+\delta\geq\phi/\|\bm{x}\|_{2}\geq 1-\delta$. Together with $1\geq\|\bm{x}_{S_{0}}\|_{2}/\|\bm{x}\|_{2}\geq 1-\delta$, we can easily obtain that $\min(\|\bm{x}^{(0)}-\bm{x}\|_{2},\|\bm{x}^{(0)}+\bm{x}\|_{2})\leq C_{0}\delta\|\bm{x}\|_{2}$ for some absolute constant $C_{0}$. By letting $\delta$ be small enough, we have $\min(\|\bm{x}^{(0)}-\bm{x}\|_{2},\|\bm{x}^{(0)}+\bm{x}\|_{2})\leq 1/6\|\bm{x}\|_{2}$.

Proof For $\bm{z}$ supported on $S$, define

Since $\operatorname{supp}(\bm{z})\subset S=\{1,\ldots,k\}$, we have

It suffices to bound $T_{1}$, $T_{2}$ and $\tau(\bm{z})$.

In what follows, we derive lower bound for $T_{11}$ and upper bound for $T_{12}$ separately.

First, by Lemma A.6 with probability at least $1-1/m$, we have

By Lemma A.5, with probability at least $1-2/m$, we have

provided $m\geq Ck^{2}$ for some sufficiently large numerical constant $C$. This implies

As to the upper bound for $T_{12}$, we can find $\|\bm{w}\|_{2}=1$, such that

By Holder’s inequality and Lemma A.5, we have

provided $m\geq Ck^{2}$, with sufficiently large constants $C_{0}$ and $C$. To summarize, with probability at least $1-3/m$,

By Lemma 6.2, letting $\delta$ small enough, we have with probability at least $1-6/m$,

provided $\mu\leq\mu_{0}$ with sufficiently small absolute constant $\mu_{0}>0$.

By Lemma A.7 and Lemma A.8, with probability at least $1-2/m-2e^{-k}$, we have

provided $m/\log m\geq k$. In summary, by Lemma 6.2, we have that with probability at least $1-5/m-2e^{-k}$,

Bound for τ​(𝒛)𝜏𝒛\tau(\bm{z})

By Holder’s inequality and Lemma A.5, with probability at least $1-2/m$, we have

for some numerical constant $C_{0}$. By Lemma A.7 and Lemma A.8, with probability at least $1-2/m-2e^{-k}$, we have,

for some numerical constant $C_{0}$, provided $\frac{m}{\log^{2}m}\geq k$. In summary,

provided $m\geq C\max(k\log p,k^{2}\sqrt{\log p})$.

Summary

We can guarantee that, with probability at least $1-\frac{15}{m}-4e^{-k}$,

for some absolute constant $C_{0}>0$, provided $m\geq Ck^{2}\log(mp)$ and $\mu\leq\mu_{0}$.

Suppose $E_{0}$ is the intersection of the events $E_{01}$ and $E_{02}$ described by Lemmas 6.3 and 6.4, respectively. Then we have

The following induction argument guarantees the effectiveness of thresholded Wirtinger flow:

provided $m\geq C\left(1+\frac{\sigma^{2}}{\|\bm{x}\|_{2}^{4}}\right)k^{2}\log(mp)$ for sufficiently large $C$.

Proof The improved estimation is defined as

where $\mathcal{T}_{\tau}$ is the soft-thresholding operator. We now define

By the definition of $\nabla f$, $\tau$ and $\phi$, as well as the assumption that \bm{x}^{(n)}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}\text{~{}and~{}}\operatorname{supp}(\bm{x}^{(n)})\subset S, we can prove $\operatorname{supp}(\bm{x}^{(n+1)})\subset S$ as well as \bm{x}^{(n+1)}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}. In fact, by the definition (2.3), we know if $\bm{x}^{(n)}$ is supported on $S$ and independent of $\bm{A}_{S^{c}}$, then $\tau(\bm{x}^{(n)})$ is independent of $\bm{A}_{S^{c}}$. Moreover, by the definition of the gradient (2.2), we know $\left(\nabla f(\bm{x}^{(n)})\right)_{S}$ is supported on $S$ and independent of $\bm{A}_{S^{c}}$. The assertion is established by the obvious fact \phi\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}} shown in Lemma 6.1.

In the following, we will construct $E_{n+1}\subset E_{n}$ such that $\widehat{\bm{x}}^{(n+1)}=\bm{x}^{(n+1)}$ on $E_{n+1}$. For any $i=k+1,k+2,\ldots,p$, with probability $1-\frac{1}{m^{2}p^{2}}$,

The first inequality is due to $\operatorname{supp}(\bm{x}^{(n)})\subset S$ and \bm{x}^{(n)}\rotatebox[origin={c}]{90.0}{\models}\bm{A}_{S^{c}}, and the second inequality is due to $\beta=4$. Then with probability at least $1-\frac{1}{m^{2}p}$,

Notice that on the event $E_{n}$, we have $\widehat{\bm{x}}^{(n)}=\bm{x}^{(n)}$, and hence

Since $E_{n}\subset E_{0}$ and $\bm{x}^{(n+1)}=\eta(\bm{x}^{(n)})$, by Lemma 6.4, we have

provided $m\geq C(\sigma^{2}/\|\bm{x}\|_{2}^{4})k\log p$ for a sufficiently large absolute constant $C$. Since $E_{n+1}\subset E_{n}$, and $\widehat{\bm{x}}^{(n+1)}=\bm{x}^{(n+1)}$ on $E_{n+1}$, we have

Theorem 3.1 can be directly implied by Lemma 6.5. In fact, by Lemma 6.3, we know the initial condition in 6.5 holds. For all $t=1,2,3,\ldots$, straight forward calculation yields

Appendix A Preliminaries and supporting lemmas

(Proposition 33 ) Consider two centered Gaussian processes $(X_{t})_{t\in T}$ and $(Y_{t})_{t\in T}$ whose increments satisfy the inequality

Proof The proof follows that of Theorem 32 in step by step. Define $X_{\bm{u},\bm{v}}=\langle\bm{A}_{S}\bm{u},\bm{v}\rangle$ on

Then $\|\bm{A}_{S}\|_{2\rightarrow 6}=\max_{(\bm{u},\bm{v})\in T}X_{\bm{u},\bm{v}}$. Define

For any $(\bm{u},\bm{v}),(\bm{u}^{\prime},\bm{v}^{\prime})\in T$, we have

due to $\|\bm{u}\|_{2}=\|\bm{u}^{\prime}\|_{2}=1$, $\|\bm{v}\|_{2}\leq\|\bm{v}\|_{6/5}=1$, and $\|\bm{v}^{\prime}\|_{2}\leq\|\bm{v}^{\prime}\|_{6/5}=1$. Then by Lemma A.3, we have

Since $\|\cdot\|_{2\rightarrow 6}$ is a $1$-Lipschitz function, by Lemma A.2, there holds with probability at least $1-2\exp(-t^{2}/2)$

Similarly, with probability at least $1-2\exp(-t^{2}/2)$

On an event with probability at least $1-1/m$, we have

provided $m\geq C(\delta)k\log k$, where $C(\delta)$ is constant only depending on $\delta$. Here $(\bm{I}_{p})_{S}$ by definition is a diagonal matrix with first $k$ diagonal entries equal to $1$, whereas other entries being . Furthermore, it implies that

The proof of this lemma is the same as that of Lemma 7.4 in .

Suppose $\epsilon_{1},\ldots,\epsilon_{m}$ are independent zero-mean sub-exponential random variables with

Then with probability at least $1-\frac{3}{m}$, we have

provided $m\geq m_{0}$ for some numerical constants $C_{0}$ and $m_{0}$.

This implies that with probability at least $1-\frac{2}{m^{10}}$, we have

provided $m\geq m_{0}$. This implies that

By the basic properties of sub-exponential random variables, for each $j=1,\ldots,m$, we have

which implies that $|\epsilon_{j}|\leq C_{0}\sigma\log m$ with probability at least $1-{e}/{m^{11}}$. This implies that

with probability at least $1-e/m^{10}$. Since

By letting $t=(4\sigma)^{2}$, we obtain that with probability at least $1-1/m$, we have $|X|\leq 20\sigma^{2}$. Similarly, with probability at least $1-1/m$, we have $\left|\frac{1}{m}\sum_{j=1}^{m}\epsilon_{j}^{4}\right|\leq C_{0}\sigma^{4}$ for some absolute constant $C_{0}$.

where $\mathcal{N}_{\frac{1}{4}}$ is the $1/4$-net of the unit sphere $\mathcal{T}^{k-1}$.

For fixed $\bm{x}\in\mathcal{N}_{\frac{1}{4}}$, let $y_{j}=|\bm{z}_{j}^{\prime}\bm{x}|^{2}-1$. Then

Notice that $y_{j}$, $j=1,\ldots,m$ are IID sub-exponential variables with $\|y_{j}\|_{\psi_{1}}\leq K$ where $K$ is an absolute constant. By Bernstein inequality (see, e.g., Proposition 16 in ), we have with probability at least $1-2\exp(-4k)$,

Since $|\mathcal{N}_{\frac{1}{4}}|\leq 9^{k}$, we know with probability at least $1-2e^{-k}$, we have

References