Reshaped Wirtinger Flow and Incremental Algorithm for Solving Quadratic System of Equations

Introduction

Many problems in machine learning and signal processing can be reduced to solve a quadratic system of equations. For instance, in phase retrieval applications, i.e., X-ray crystallography and coherent diffraction imaging , the structure of an object is to be recovered from the measured far field diffracted intensity when an object is illuminated by a source light. Mathematically, such a problem amounts to recover the signal from only measurements of magnitudes. Specifically, the problem is formulated below.

Various algorithms have been proposed to solve this problem since 1970s. The error-reduction methods proposed in work well empirically but lack theoretical guarantees. More recently, convex relaxation of the problem has been formulated, for example, via PhaseLift and PhaseCut , and the correspondingly developed algorithms typically come with performance guarantee. The reader can refer to the review paper to learn more about applications and algorithms of the phase retrieval problem.

While with good theoretical guarantee, these convex methods often suffer from computational complexity particularly when the signal dimension is large. On the other hand, more efficient nonconvex approaches have been proposed and shown to recover the true signal as long as initialization is good enough. proposed AltMinPhase algorithm, which alternatively updates the phase and the signal with each signal update solving a least-squares problem, and showed that AltMinPhase converges linearly and recovers the true signal with $\mathcal{O}(n\log^{3}n)$ Gaussian measurements. More recently, introduces Wirtinger flow (WF) algorithm, which guarantees signal recovery via a simple gradient algorithm with only $\mathcal{O}(n\log n)$ Gaussian measurements and attains $\epsilon-$accuracy within $\mathcal{O}(mn^{2}\log 1/\epsilon)$ flops. More specifically, WF obtains good initialization by the spectral method, and then minimizes the following nonconvex loss function

WF is further improved by truncated Wirtinger flow (TWF) algorithm proposed in , which adopts a Poisson loss function of $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|^{2}$, and keeps only well-behaved measurements based on carefully designed truncation thresholds for calculating the initial seed and every step of gradient. Such truncation assists to yield linear convergence with certain fixed step size and reduces both the sample complexity to $\mathcal{O}(n)$ and the convergence time to $\mathcal{O}(mn\log 1/\epsilon)$.

It can be observed that WF uses the quadratic loss of $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|^{2}$ so that the optimization objective is a smooth function of $\boldsymbol{a}_{i}^{T}\boldsymbol{z}$ and the gradient step becomes simple. But this comes with a cost of increasing the order of $\boldsymbol{a}_{i}^{T}\boldsymbol{z}$ to be four in the loss function. In this paper, we adopt the quadratic loss of $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|$. Although the loss function is not smooth everywhere, it reduces the order of $\boldsymbol{a}_{i}^{T}\boldsymbol{z}$ to be two, and the general curvature can be more amenable to convergence of the gradient method. The goal of this paper is to explore potential advantages of such a nonsmooth lower-order loss function.

Furthermore, incremental/stochastic methods have been proposed to solve Problem 1. Specifically, Kaczmarz method for phase retrieval (Kaczmarz-PR) is shown to have excellent empirical performance, but no global convergence guarantee was established. Incremental truncated Wirtinger flow (ITWF) is a stochastic algorithm developed based on TWF and exhibits linear convergence to the true signal once initialized properly. In this paper, we consider the incremental/stochastic version of RWF (IRWF) and compare its performance with Kaczmarz-PR and ITWF, in order to further demonstrate the advantage of the lower-order loss function.

This paper adopts the following loss functionThe loss function (3) was also used in to derive a gradient-like update for the phase retrieval problem with Fourier magnitude measurements. However, the focus of this paper is to characterize global convergence guarantee for such an algorithm with appropriate initialization, which was not studied in .

Compared to the loss function (2) in WF that adopts $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|^{2}$, the above loss function adopts the absolute value/magnitude $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|$ and hence has lower-order variables. For such a nonconvex and nonsmooth loss function, we develop a gradient descent-like algorithm, which sets zero for the “gradient” component corresponding to nonsmooth samples. We refer to such an algorithm together with an initialization using a new spectral method (different from that employed in TWF or WF) as reshaped Wirtinger flow (RWF). We show that the lower-order loss function has great advantage in both statistical and computational efficiency, although scarifying smoothness. In fact, the curvature of such a loss function behaves similarly to that of a least-squares loss function in the neighborhood of global optimums (see Section 2.2), and hence RWF converges fast. The nonsmoothness does not significantly affect the convergence of the algorithm because only with negligible probability the algorithm encounters nonsmooth points for some samples, which furthermore are set not to contribute to the gradient direction by the algorithm.

We further exploit the loss function (3) to design the incremental/stochastic reshaped Wirtinger flow (IRWF), and we show that IRWF also enjoys the advantageous local curvature of such a loss function and achieves excellent statistical and computation performance. In particular, IRWF performs better than other competitive incremental methods (ITWF and Kaczmarz-PR) as well as batch algorithms (RWF, TWF and WF) numerically.

We summarize our main results as follows.

Statistically, we show that RWF recovers the true signal with $\mathcal{O}(n)$ samples, when the design vectors consist of independently and identically distributed (i.i.d.) Gaussian entries, which is optimal in the order sense. Thus, even without truncation in gradient steps (truncation only in initialization stage), RWF improves the sample complexity $\mathcal{O}(n\log n)$ of WF, and achieves the same sample complexity as TWF with truncation in gradient step. It is thus more robust to random measurements.

Computationally, RWF converges geometrically to the true signal, requiring $\mathcal{O}(mn\log(1/\epsilon))$ flops to reach $\epsilon-$accuracy. Again, without truncation in gradient steps, RWF improves computational cost $\mathcal{O}(mn^{2}\log(1/\epsilon)$ of WF and achieves the same computational cost as TWF.

Numerically, RWF adopts fixed step size and does not require truncation, which avoids the trouble to set truncation thresholds in practice. It is generally two times faster than TWF and four to six times faster than WF in terms of the number of iterations and time cost.

We also show that RWF is robust to bounded additive noise. The estimation error is shown to diminish geometrically to the power of bounded noise up to a certain coefficient. Experiments on Poisson noise further corroborate the stability guarantee.

We further propose incremental/stochastic algorithm based on RWF (referred to as IRWF) and show that IRWF converges to the true signal geometrically under an appropriate initialization. More interestingly, we show that randomized Kaczmarz-PR (i.e., Kaczmarz method adapted for phase retrieval ) can be viewed as IRWF under a specific way of choosing step size, via which we further established geometric convergence guarantee for randomized Kaczmarz-PR.

Compared to WF and TWF, the new form of the gradient step due to nonsmoothness of the loss function, in terms of technical analysis, requires new developments of bounding techniques. On the other hand, our technical proof of performance guarantee is much simpler, because the lower-order loss function allows to bypass higher-order moments of variables and truncation in gradient steps. We also anticipate that such analysis is more easily extendable.

2 Connection to Related Work

Along the line of developing nonconvex algorithms with global performance guarantee for the phase retrieval problem, developed alternating minimization algorithm, developed/studied first-order gradient-like algorithms, and a recent study characterized geometric structure of the nonconvex objective and designed a second-order trust-region algorithm. This paper is most closely related to , but develops a new gradient-like algorithm based on a lower-order nonsmooth (as well as nonconvex) loss function that yields advantageous statistical/computational efficiency.

Stochastic algorithms are also developed for the phase retrieval problem. studied the incremental truncated Wirtinger flow (ITWF) and showed that ITWF needs much fewer passes of data than TWF to reach the same accuracy. adapted the Kaczmarz method to solve the phase retrieval problem and demonstrated its fast empirical convergence. We propose a stochastic algorithm based on RWF (IRWF) and show that it has close connection with Kaczmarz-PR. We also show that IRWF runs faster than ITWF due to the benefit of low-order loss function.

After our work was posted on arXiv, an independent work was subsequently posted, which also adopts the same loss function but develops a slightly different algorithm TAF (i.e., truncated amplitude flow). One major difference of our algorithm RWF from TAF is that RWF does not require truncation in gradient loops while TAF employs truncation. Hence, RWF has fewer parameters to tune, and is easier to implement than TAF in practice. Furthermore, RWF demonstrates the performance advantage of adopting a lower-order loss function even without truncation, which cannot be observed from TAF. Moreover, we analyze stochastic algorithm based on new loss function while does not.

More generally, various problems have been studied by minimizing nonconvex loss functions. For example, a partial list of these studies include matrix completion , low-rank matrix recovery , robust PCA , robust tensor decomposition , dictionary learning , community detection , phase synchronization, blind deconvolution , etc.

For minimizing a general nonconvex nonsmooth objective, various algorithms have been proposed, such as gradient sampling algorithm and majorization-minimization method . These algorithms were often shown to convergence to critical points which may be local minimizers or saddle points, without explicit characterization of convergence rate. In contrast, our algorithm is specifically designed for the phase retrieval problem, and can be shown to converge linearly to global optimum under appropriate initialization.

The advantage of nonsmooth loss function exhibiting in our study is analogous in spirit to that of the rectifier activation function (of the form $\max\{0,\cdot\}$) in neural networks. It has been shown that rectified linear unit (ReLU) enjoys superb advantage in reducing the training time and promoting sparsity over its counterparts of sigmoid and hyperbolic tangent functions, in spite of non-linearity and non-differentiability at zero. Our result in fact also demonstrates that a nonsmooth but simpler loss function yields improved performance.

3 Paper Organization and Notations

The rest of this paper is organized as follows. Section 2 describes RWF algorithm in detail and establishes its performance guarantee. Section 3 introduces the IRWF algorithm and establishes its performance guarantee and compares it with existing stochastic algorithms. Section 4 compares RWF and IRWF with other competitive algorithms numerically. Finally, Section 5 concludes the paper with comments on future directions.

Throughout the paper, boldface lowercase letters such as $\boldsymbol{a}_{i},\boldsymbol{x},\boldsymbol{z}$ denote vectors, and boldface capital letters such as $\boldsymbol{A},\boldsymbol{Y}$ denote matrices. For two matrices, $\boldsymbol{A}\preceq\boldsymbol{B}$ means that $\boldsymbol{B}-\boldsymbol{A}$ is positive definite. For a complex matrix or vector, $\boldsymbol{A}^{*}$ and $\boldsymbol{z}^{*}$ denote conjugate transposes of $\boldsymbol{A}$ and $\boldsymbol{z}$ respectively. For a real matrix or vector, $\boldsymbol{A}^{T}$ and $\boldsymbol{z}^{T}$ denote transposes of $\boldsymbol{A}$ and $\boldsymbol{z}$ respectively. The indicator function $\boldsymbol{1}_{A}=1$ if the event $A$ is true, and $\boldsymbol{1}_{A}=0$ otherwise.

Reshaped Wirtinger Flow

where it is simply $\min\|\boldsymbol{z}\pm\boldsymbol{x}\|$ for real case. We focus on the real-valued case in analysis, but the algorithm designed below is applicable to the complex-valued case and the case with coded diffraction pattern (CDP) as we demonstrate via numerical results in Section 4.

We design RWF (see Algorithm 1) for solving the above problem, which contains two stages: spectral initialization and gradient loop. Suggested values for parameters are $\alpha_{l}=1,\alpha_{u}=5$ and $\mu=0.8$For complex Gaussian case, we suggest $\mu=1.2$.. The scaling parameter in $\lambda_{0}$ and the conjugate transpose $\boldsymbol{a}_{i}^{*}$ allow the algorithm readily applicable to complex and CDP cases. We next describe the two stages of the algorithm in detail in Sections 2.1 and 2.2, respectively, and establish the convergence of the algorithm in Section 2.3. Finally, we provide the stability result of RWF in Section 2.4.

Differently from the spectral initialization methods for WF in and for TWF in , both of which are based on $|\boldsymbol{a}_{i}^{*}\boldsymbol{x}|^{2}$, we propose an alternative initialization in Algorithm 1 that uses magnitude $|\boldsymbol{a}_{i}^{*}\boldsymbol{x}|$ instead, and truncates samples with both lower and upper thresholds as in (6). We show that such initialization achieves smaller sample complexity than WF and the same order-level sample complexity as TWF, and furthermore, performs better than both WF and TWF numerically.

Our initialization consists of estimation of both the norm and direction of $\boldsymbol{x}$. The norm estimation of $\boldsymbol{x}$ is given by $\lambda_{0}$ in Algorithm 1 with mathematical justification in Appendix A. Intuitively, with real Gaussian measurements, the scaling coefficient $\frac{mn}{\sum_{i=1}^{m}\|\boldsymbol{a}_{i}\|_{1}}\approx\sqrt{\frac{\pi}{2}}$. Moreover, $y_{i}=|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|$ are independent sub-Gaussian random variables for $i=1,\ldots,m$ with mean $\sqrt{\frac{2}{\pi}}\|\boldsymbol{x}\|$, and thus $\frac{1}{m}\sum_{i=1}^{m}y_{i}\approx\sqrt{\frac{2}{\pi}}\|\boldsymbol{x}\|$. Combining these two facts yields the desired argument.

We next provide the formal statement of the performance guarantee for the initialization step that we propose.

Fix $\delta>0$. The initialization step in Algorithm 1 yields $\boldsymbol{z}^{(0)}$ satisfying $\text{dist}(\boldsymbol{z}^{(0)},\boldsymbol{x})\leq\delta\|\boldsymbol{x}\|$ with probability at least $1-\exp(-cm\epsilon^{2})$, if $m>C(\delta,\epsilon)n$, where $c$ is some positive constant and $C$ is a positive number only affected by $\delta$ and $\epsilon$.

Finally, Figure 1 demonstrates that RWF achieves better initialization accuracy in terms of the relative error $\frac{\text{dist}(\boldsymbol{z}^{(0)},\boldsymbol{x})}{\|\boldsymbol{x}\|}$ than WF and TWF. Furthermore, we also include the orthonormal promoting initialization method proposed for truncated amplitude flow (TAF) in the independent work , in the comparison. It can be seen that our initialization is slightly better.

2 Gradient Loop and Why RWF is Fast

The gradient loop of Algorithm 1 is based on the loss function (3), which is rewritten below

We let the update direction be as follows:

We next provide the intuition about why reshaped WF is fast. Suppose that the spectral method sets an initial point in the neighborhood of the ground truth $\boldsymbol{x}$. We compare RWF with the following problem of solving $\boldsymbol{x}$ from linear equations $y_{i}=\langle\boldsymbol{a}_{i},\mathbf{x}\rangle$ with $y_{i}$ and $\boldsymbol{a}_{i}$ for $i=1,\ldots,m$ given. In particular, we note that this problem has both magnitude and sign of the measurements observed. In this case, it is natural to use the least-squares loss. Assume that the gradient descent is applied to solve this problem. Then the gradient is given by

We now argue intuitively that the gradient (9) of RWF behaves similarly to the least-squares gradient (10). For each $i$, $y_{i}\frac{\boldsymbol{a}_{i}^{T}\boldsymbol{z}}{|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|}=|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|\text{sgn}(\boldsymbol{a}_{i}^{T}\boldsymbol{z})$, and hence the two gradient components are close if $|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|\cdot\text{sgn}(\boldsymbol{a}_{i}^{T}\boldsymbol{z})$ is viewed as an estimate of $\boldsymbol{a}_{i}^{T}\boldsymbol{x}$. The following lemma shows that if $\text{dist}(\boldsymbol{z},\boldsymbol{x})$ is small (guaranteed by initialization), then $\boldsymbol{a}_{i}^{T}\boldsymbol{z}$ has the same sign as $\boldsymbol{a}_{i}^{T}\boldsymbol{x}$ for large $|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|$.

Let $\boldsymbol{a}_{i}\sim\mathcal{N}(0,\boldsymbol{I}_{n\times n})$. For any given $\boldsymbol{x}$ and $\boldsymbol{z}$, independent from $\{\boldsymbol{a}_{i}\}_{i=1}^{m}$, satisfying $\|\boldsymbol{x}-\boldsymbol{z}\|<\frac{\sqrt{2}-1}{\sqrt{2}}\|\boldsymbol{x}\|$, we have

where $\text{erfc}(u):=\frac{2}{\sqrt{\pi}}\int_{u}^{\infty}\exp(-\tau^{2})d\tau$.

It is easy to observe in (11) that large $\boldsymbol{a}_{i}^{T}\boldsymbol{x}$ is likely to have the same sign as $\boldsymbol{a}_{i}^{T}\boldsymbol{z}$ so that the corresponding gradient components in (9) and (10) are likely equal, whereas small $\boldsymbol{a}_{i}^{T}\boldsymbol{x}$ may have different sign as $\boldsymbol{a}_{i}^{T}\boldsymbol{z}$ but contributes less to the gradient. Hence, overall the two gradients (9) and (10) should be close to each other with a large probability.

This fact can be further verified numerically. Figure 2(a) illustrates that RWF takes almost the same number of iterations for recovering a signal (with only magnitude information) as the least-squares gradient descent method for recovering a signal (with both magnitude and sign information).

𝒛𝒙𝑇𝒛𝒙f(z)=\min\{(\boldsymbol{z}-\boldsymbol{x})^{T}(\boldsymbol{z}-\boldsymbol{x}),(\boldsymbol{z}+\boldsymbol{x})^{T}(\boldsymbol{z}+\boldsymbol{x})\} with $\boldsymbol{x}=^{T}$. (b) Expected loss function of RWF for $\boldsymbol{x}=^{T}$. (c) Expected loss function of WF for $\boldsymbol{x}=^{T}$. Furthermore, Figure 3 illustrates a quadratic function $f(\cdot)$ (Figure 3(a)), the expected loss surface of RWF (Figure 3(b), and see Appendix B for expression), and the expected loss surface for WF (Figure 3(c), and see Appendix B for expression). It can be seen that the loss of RWF rather than the loss of WF has a similar curvature to the quadratic function around the global optimums, which justifies its better performance than WF.

3 Geometric Convergence of RWF

We characterize the convergence of RWF in the following theorem.

We outline the proof here with details relegated to Appendix C. Compared to WF and TWF, our proof requires new development of bounding techniques to deal with nonsmoothness, but is much simpler due to the lower-order loss function that RWF relies on.

We first introduce a global phase notation for real case as follows:

For the sake of simplicity, we let $\boldsymbol{z}$ be $e^{-j\Phi(\boldsymbol{z})}\boldsymbol{z}$, which indicates that $\boldsymbol{z}$ is always in the neighborhood of $\boldsymbol{x}$.

Here, the central idea is to show that within the neighborhood of global optimums, RWF satisfies the Regularity Condition $\mathsf{RC}(\mu,\lambda,c)$ , i.e.,

for all $\boldsymbol{z}$ obeying $\|\boldsymbol{z}-\boldsymbol{x}\|\leq c\|\boldsymbol{x}\|$, where $0<c<1$ is some constant. Then, as shown in , once the initialization lands into this neighborhood, geometric convergence can be guaranteed, i.e.,

for any $\boldsymbol{z}$ with $\|\boldsymbol{z}-\boldsymbol{x}\|\leq c\|\boldsymbol{x}\|$.

And Lemma 4 in Appendix C further yields that

Therefore, the above two bounds imply that Regularity Condition (14) holds for $\mu$ and $\lambda$ satisfying

We note that (17) implies an upper bound $\mu\leq\frac{0.74}{2}=0.37$, by taking $\epsilon$ and $\delta$ to be sufficiently small. This suggests a range to set the step size in Algorithm 1. However, in practice, $\mu$ can be set much larger than such a bound, say $0.8$, while still keeping the algorithm convergent. This is because the coefficients in the proof are set for convenience of proof rather than being tightly chosen.

Theorem 1 indicates that RWF recovers the true signal with $\mathcal{O}(n)$ samples, which is order-level optimal. Such an algorithm improves the sample complexity $\mathcal{O}(n\log n)$ of WF. Furthermore, RWF does not require truncation of weak samples in the gradient step to achieve the same sample complexity as TWF. This is mainly because RWF benefits from the lower-order loss function given in (8), the curvature of which behaves similarly to the least-squares loss function locally as we explain in Section 2.2.

Theorem 1 also suggests that RWF converges geometrically at a constant step size. To reach $\epsilon-$accuracy, it requires computational cost of $\mathcal{O}(mn\log 1/\epsilon)$ flops, which is better than WF ($\mathcal{O}(mn^{2}\log(1/\epsilon)$). Furthermore, it does not require truncation in gradient steps to reach the same computational cost as TWF. Numerically, as we demonstrate in Section 4, RWF is two times faster than TWF and four to six times faster than WF in terms of both iteration count and time cost in various examples.

4 Stability to Bounded Noise

We have established that RWF guarantees exact recovery with geometric convergence for noise-free case. We now study RWF in the presence of noise. Suppose the measurements are corrupted by bounded noise, and are given by

where $\{w_{i}\}_{i=1}^{m}$ denote the additive noise. Then the following theorem shows that RWF is robust under such noise corruption.

Consider the model (18). Suppose that the measurement vectors are independently Gaussian, i.e., $\boldsymbol{a}_{i}\sim\mathcal{N}(0,\boldsymbol{I})$ for $1\leq i\leq m$, and the noise is bounded, i.e., $\|\boldsymbol{w}\|/\sqrt{m}\leq c\|\boldsymbol{x}\|$. Then there exist some universal constants $\mu_{0}>0$ ($\mu_{0}$ can be set as $0.8$ in practice), $0<\rho<1$ and $c_{0},c_{1},c_{2}>0$ such that if $m\geq c_{0}n$ and $\mu<\mu_{0}$, then with probability at least $1-c_{1}\exp(-c_{2}m)$, Algorithm 1 yields

The numerical result under the Poisson noise model in Section 4 further corroborates the stability of RWF.

Incremental Reshaped Wirtinger Flow

In large sample size and online scenarios, stochastic algorithm is preferred due to its potential advantage of fast convergence and low computational complexity. Thus, in this section, we develop the stochastic algorithm based on RWF, named incremental reshaped Wirtinger flow (IRWF). We show that IRWF guarantees exact recovery with linear convergence rate. We further draw the connection between IRWF and the Kaczmarz-PR algorithm recently developed for phase retrieval .

We describe IRWF in Algorithm 2. More specifically, IRWF applies the same initialization step as in RWF, but calculates each gradient update using only one sample selected randomly.

We characterize the convergence of IRWF in the following theorem.

We suggest the step size $\rho_{0}=1$ in practice.

The proof is relegated to Appendix E.1, which uses technical lemmas established for proving Theorem 1. ∎

Theorem 3 establishes that IRWF achieves linear convergence to the global optimum. For general objectives, it is not anticipated that incremental/stochastic first-order method achieves linear convergence due to the variance of stochastic gradient. However, for our specific problem, the variance of stochastic gradient reduces as the estimate approaches the true solution, and hence a fixed step size can be employed and linear convergence can be established. Such a result is also established in for the stochastic algorithm based on TWF (ITWF). We comment further on comparison between our algorithm and ITWF in Section 3.4. Another explanation may be due to the fact that stochastic gradient method yields linear convergence to the minimizer $\boldsymbol{x}_{\star}$ when the objective $F(\boldsymbol{x})=\sum_{i}f_{i}(\boldsymbol{x})$ is a smooth and strongly convex function and $\boldsymbol{x}_{\star}$ minimizes all components $f_{i}(\boldsymbol{x})$. This may also hold for our objective (3) whose summands share a same minimizer, although it is neither convex nor smooth.

2 Minibatch IRWF: Algorithm and Convergence

In oder to fully exploit the processing throughput of CPU/GPU, we develop a minibatch version of IRWF, described in Algorithm 3. The minibatch IRWF applies the same initialization step as in RWF, but uses a minibatch of data for each gradient update in contrast to IRWF that uses only a single sample.

We characterize the convergence of minibatch IRWF in the following theorem.

We suggest that $\rho_{0}=1$ in practice.

3 Connection to Kaczmarz Method for Phase Retrieval

Kaczmarz method was originally developed for solving the linear equation systems $\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b}$ . In , it was adapted to solve phase retrieval problem, which we refer to as Kaczmarz-PR. It has been demonstrated in that Kaczmarz-PR exhibits better empirical performance than error reduction (ER) and Wirtinger flow (WF) . However, theoretical guarantee of Kaczmarz-PR has not been well established yet although Kaczmarz method for least-squares problem achieves linear convergence guarantee . For instance, obtained a bound on the estimation error which can be as large as the signal energy no matter how many iterations are taken. requires infinite number of samples to establish the asymptotic convergence.

In this section, we draw the connection between IRWF and Kaczmarz-PR, and thus the theoretical guarantee of Kaczmarz-PR can be established by adapting that of IRWF. This is analogous to the connection established in between Kaczmarz method and stochastic gradient method when solving the least-squares problem. Here, the connection is interesting because RWF rather than WF and TWF connects to Kaczmarz-PR due to the lower-order loss function that RWF adopts.

To be more specific, the Kaczmarz-PR (Algorithm 3 in ) employs the following update rule

where $i_{t}$ is selected either in a deterministic manner or randomly. We focus on the randomized Kaczmarz-PR where $i_{t}$ is selected uniformly at random.

Comparing (24) and (20), the update rule of randomized Kaczmarz is a special case of IRWF with step size $\mu$ replaced by $\frac{1}{\|\boldsymbol{a}_{i_{t}}\|^{2}}$. Moreover, these two update rules are close if $\mu$ is set as suggested, i.e., $\mu=\frac{1}{n}$, because $\|\boldsymbol{a}_{i_{t}}\|^{2}$ concentrates around $n$ by law of large numbers. As we demonstrate in empirical results (see Table 1), these two methods have similar performance as anticipated. Thus, following the convergence result Theorem 3 for IRWF, we have the convergence guarantee for the randomized Kaczmarz-PR as follows.

Assume the measurement vectors are independent and each $\boldsymbol{a}_{i}\sim\mathcal{N}(0,\boldsymbol{I})$. There exist some universal constants $0<\rho<1$ and $c_{0},c_{1},c_{2}>0$ such that if $m\geq c_{0}n$, then with probability at least $1-c_{1}m\exp(-c_{2}n)$, the randomized Kaczmarz update rule (24) yields

holds for all $\boldsymbol{z}^{(t)}$ satisfying $\frac{\text{dist}(\boldsymbol{z}^{(t)},\boldsymbol{x})}{\|\boldsymbol{z}\|}\leq\frac{1}{10}$.

The above theorem implies that once the estimate $\boldsymbol{z}^{(t)}$ enters the neighborhood of true solutions (often referred as to basin of attraction), the error diminishes geometrically by each update in expectation.

Furthermore, also provided a block Kaczmarz-PR (similar to the minibatch version), whose update rule is given by

where $\odot$ is element-wise product, $\text{Ph}(\boldsymbol{a})$ denotes a phase vector of $\boldsymbol{a}$ (stacking phase of each element of $\boldsymbol{a}$ together), $\boldsymbol{A}$ is a matrix with each row being measurement vector $\boldsymbol{a}_{i}^{*}$, and $\Gamma_{t}$ is a selected block at iterate $t$ containing row indices. Moreover, $\dagger$ represents Moore-Penrose pseudoinverse, which can be computed as follows:

Comparing (26) and the minibatch IRWF update in (22), these two update rules are similar to each other if $\boldsymbol{A}_{\Gamma_{t}}\boldsymbol{A}_{\Gamma_{t}}^{*}$ approaches $\frac{n}{\rho_{0}}\boldsymbol{I}_{|\Gamma_{t}|}$. For the case with Gaussian measurements, $\boldsymbol{A}_{\Gamma_{t}}$ has linearly independent rows with high probability if $|\Gamma_{t}|\leq n$ and hence $\boldsymbol{A}_{\Gamma_{t}}\boldsymbol{A}_{\Gamma_{t}}^{*}$ is not far from $n\boldsymbol{I}_{|\Gamma_{t}|}$. Our empirical results (see Table 1) further suggest similar convergence rate for these two algorithms with the same block size.

Next, we argue that for the CDP setting, block Kaczmarz-PR is the same as the minibatch IRWF with $\mu=\frac{1}{n}$. The CDP measurements are collected in the following form

where $\boldsymbol{F}$ represents the discrete Fourier transform (DFT) matrix, and $\boldsymbol{D}^{(l)}$ is a diagonal matrix (mask). We choose the block size $|\Gamma_{t}|$ to be $n$, the dimension of the signal, for the convenience of Fourier transform. Then $\boldsymbol{A}_{\Gamma_{t}}$ becomes a Fourier transform composed with $\boldsymbol{D}^{(l)}$ (mask effect) and $\boldsymbol{A}_{\Gamma_{t}}^{*}$ becomes $\boldsymbol{D}^{(l)*}$ multiplied by inverse Fourier transform. Therefore, $(\boldsymbol{A}_{\Gamma_{t}}\boldsymbol{A}_{\Gamma_{t}}^{*})=\boldsymbol{I}$ if the diagonal elements of $\boldsymbol{D}^{(l)}$ have unit magnitude. Taking the step size $\mu=1$, the two algorithms are identical.

On the other hand, since the block Kaczmarz-PR needs to calculate the matrix inverse or to solve an inverse problem, the block size cannot be too large. However, minibatch IRWF works well for a wide range of batch sizes which can even vary with signal dimension $n$ as long as a batch of data is loadable into memory.

4 Comparison with Incremental Truncated Wirtinger Flow (ITWF)

Recently, designed and analyzed an incremental algorithm based on TWF, which is referred to as ITWF. More specifically, ITWF employs the same initialization procedure as TWF and randomly chooses one sample for gradient update as follows.

where $\boldsymbol{1}_{\mathcal{E}_{1,t}^{i_{t}}\cap\mathcal{E}_{3}^{i_{t}}}$ is the truncation rule determined by two events $\mathcal{E}_{1,t}^{i_{t}}$ and $\mathcal{E}_{3}^{i_{t}}$. Compared to ITWF developed based on TWF, our IRWF uses lower-order variable $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|$ rather than $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|^{2}$ used in ITWF. Moreover, IRWF does not employ any truncation in gradient loops and hence has fewer parameters to tune, which is easier to implement in practice.

proved that ITWF converges linearly to the true signal as long as $m/n$ (sample size/ signal dimension) is large enough. Compared to ITWF, IRWF also achieves the same linear convergence, but runs faster than ITWF numerically as demonstrated in Section 4. The proof of IRWF requires different bounding techniques, but is simpler than ITWF due to the lower-order loss function that RWF adopts and the avoidance of truncation in the gradient update.

Numerical Results

In this section, we demonstrate the numerical efficiency of RWF and IRWF by comparing their performances with other competitive algorithms. Our experiments are conducted not only for real Gaussian case but also for complex Gaussian and the CDP cases. All the experiments are implemented in Matlab 2015b and carried out on a computer equipped with Intel Core i7 3.4GHz CPU and 12GB RAM.

We first compare the sample complexity of RWF and IRWF with those of TWF and WF via the empirical successful recovery rate versus the number of measurements. For RWF, we follow Algorithm 1 with suggested parameters. For IRWF, we adopt a block size 64 for efficiency and set the step size $\rho_{0}=1$. For WF, TWF, we use the codes provided in the original papers with the suggested parameters. For ITWF, we also adopt a block size 64 and set the step size $\rho_{0}=0.6$ (optimal step size). We conduct the experiment for real Gaussian, complex Gaussian and CDP cases respectively. For real and complex cases, we set the signal dimension $n$ to be 1000, and set the ratio $m/n$ to take values from $2$ to $6$ by a step size $0.1$. For each $m$, we run $100$ trials and count the number of successful trials. For each trial, we run a fixed number of iterations/passes $T=1000$ for all algorithms. A trial is declared to be successful if $\boldsymbol{z}^{(T)}$, the output of the algorithm, satisfies $\text{dist}(\boldsymbol{z}^{(T)},\boldsymbol{x})/\|\boldsymbol{x}\|\leq 10^{-5}$. For the real Gaussian case, we generate signal $\boldsymbol{x}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{n\times n})$, and the measurement vectors $\boldsymbol{a}_{i}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{I}_{n\times n})$ i.i.d. for $i=1,\ldots,m$. For the complex Gaussian case, we generate signal $\boldsymbol{x}\sim\mathcal{N}(0,\boldsymbol{I}_{n\times n})+j\mathcal{N}(0,\boldsymbol{I}_{n\times n})$ and measurements $\boldsymbol{a}_{i}\sim\frac{1}{2}\mathcal{N}(0,\boldsymbol{I}_{n\times n})+j\frac{1}{2}\mathcal{N}(0,\boldsymbol{I}_{n\times n})$ i.i.d. for $i=1,\ldots,m$. For the CDP case (28), we set $n=1024$ for convenience of FFT and $m/n=L=1,2,\ldots,8$. All other settings are the same as those for the real case.

Figure 4 plots the fraction of successful trials out of 100 trials for all algorithms, with respect to $m$. It can be further seen that IRWF exhibits the best sample complexity for all three cases, which is close to the theoretical limits . It can be seen that the two incremental methods (IRWF and ITWF) outperform batch methods (RWF, TWF and WF). This can be due to the inherent noise in incremental methods helps to escape bad local minimums, which is extremely helpful in the regime of small number of samples where local minimums do exist near the global ones. Comparing among the three batch methods (RWF, TWF and WF), it can be seen that although RWF outperforms only WF (not TWF) for the real Gaussian case, it outperforms both WF and TWF for complex Gaussian and CDP cases. An intuitive explanation for the real case is that a substantial number of samples with small $|\boldsymbol{a}_{i}^{T}\boldsymbol{z}|$ can deviate gradient so that truncation indeed helps to stabilize the algorithm if the number of measurements is not large. Furthermore, RWF exhibits sharper transition than TWF and WF.

We next compare the convergence rate of RWF, IRWF with those of TWF, ITWF, WF and AltMinPhase. We run all of the algorithms with suggested parameter settings in the original codes. We generate signal and measurements in the same way as those in the first experiment with $n=1024,m=8n$. All algorithms are seeded with RWF initialization. In Table 1, we list the number of passes and time cost for those algorithms to achieve the relative error of $10^{-14}$ averaged over 10 trials. For incremental methods, one update passes $k$ samples and one pass amounts to $m/k$ updates. Clearly, IRWF with minibatch size 64 runs fastest for both real and complex cases. Moreover, among batch (deterministic) algorithms, RWF takes much less number of passes as well as running much faster than TWF and WF. Although RWF takes more iterations than AltMinPhase, it runs much faster than AltMinPhase due to the fact that each iteration of AltMinPhase needs to solve a least-squares problem that takes much longer time than a simple gradient update in RWF.

We also compare the performance of the above algorithms on the recovery of a real image from the Fourier intensity measurements (two dimensional CDP case). The image (see Figure 5) is the Milky Way Galaxy with resolution $1920\times 1080$. Table 2 lists the number of passes and the time cost of the above six algorithms to achieve the relative error of $10^{-15}$ for one R/G/B channel. All algorithms are seeded with RWF initialization. To explore the advantage of FFT, we run the incremental/stochastic methods with minibatch size of the one R/G/B channel. We note that with such a minibatch size IRWF is equivalent to block Kaczmarz-PR from the discussion in Section 3.3. It can be seen that in general, the incremental/stochastic methods (IRWF/Kaczmarz-PR and ITWF) run faster than the batch methods (RWF, TWF, WF and AltMinPhase). Moreover, among batch methods, RWF outperforms other three algorithms in both number of passes and the computational time cost. In particular, RWF runs two times faster than TWF and six times faster than WF in terms of both the number of iterations and computational time cost.

We next demonstrate the robustness of RWF to noise corruption and compare it with TWF. We consider the phase retrieval problem in imaging applications, where random Poisson noises are often used to model the sensor and electronic noise . Specifically, the noisy measurements of intensity can be expressed as y${}_{i}=\sqrt{\alpha\cdot\text{Poisson}\left(|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|^{2}/\alpha\right)},\quad\text{for }i=1,2,...m$ where $\alpha$ denotes the level of input noise, and $\text{Poisson}(\lambda)$ denotes a random sample generated by the Poisson distribution with mean $\lambda$. It can be observed from Figure 6 that RWF performs better than TWF in terms of recovery accuracy under different noise levels.

Conclusion

In this paper, we proposed RWF and its stochastic version IRWF to recover a signal from a quadratic systems of equations, based on a nonconvex and nonsmooth quadratic loss function of absolute values of measurements. This loss function sacrifices the smoothness but enjoys advantages in statistical and computational efficiency. It has potential to be extended in various scenarios. One interesting direction is to extend such an algorithm to exploit signal structures (e.g., non-negativity, sparsity, etc) to assist the recovery. The lower-order loss function may offer great simplicity to prove performance guarantee in such cases.

Another interesting direction is to study the convergence of algorithms from random initialization. In the regime of large sample size ($m\gg n$), the empirical loss surface approaches the asymptotic loss (Figure 3(b)) and hence has no spurious local minimums. Due to , it is conceivable that gradient descent converges from random starting point. Similar phenomenons have been observed in . However, under moderate number of measurements ($m<10n$), authentic local minimums do exist and often locate not far from the global ones. In this regime, the batch gradient method often fails with random initialization. As always believed, stochastic algorithms are efficient in escaping bad local minimums or saddle points in nonconvex optimization because of the inherent noise . We observe numerically that IRWF and block IRWF from random starting point still converge to global minimum even with very small sample size which is close to the theoretical limits . It is of interest to analyze theoretically that stochastic methods escape these local minimums (not just saddle points) efficiently.

We let $\|\cdot\|_{1}$ and $\|\cdot\|$ denote the $l_{1}$ norm and $l_{2}$ norm of a vector, respectively. Moreover, let $\|\cdot\|_{F}$ and $\|\cdot\|$ denote the Frobenius norm and the spectral norm of a matrix, respectively. We note that the constants $c,C,c_{0},c_{1},c_{2}$ may be different from line to line, for the sake of notational simplicity.

Appendix A Proof of Proposition 1: Performance Guarantee for Initialization

Compared to the proof for TWF , this proof has new technical developments to address the magnitude measurements and truncation from both sides.

We first estimate the norm of $\boldsymbol{x}$ as

Since $\boldsymbol{a}_{i}\sim\mathcal{N}(0,\boldsymbol{I}_{n\times n})$, by Hoeffding-type inequality, it can be shown that

holds with probability at least $1-2\exp(-c_{1}mn\epsilon^{2})$ for some constant $c_{1}>0$.

Moreover, given $\boldsymbol{x}$, $y_{i}$’s are independent sub-Gaussian random variables. Thus, by Hoeffding-type inequality, it can be shown that

holds with probability at least $1-2\exp(-c_{1}m\epsilon^{2})$ for some constant $c_{1}>0$.

On the event $E_{1}=\{\text{both }\eqref{eq:scalingest}\text{ and }\eqref{eq:normestbound}\text{ hold}\}$, it can be argued that

Without loss of generality, we let $\|\boldsymbol{x}\|=1$. Then on the event $E_{1}$, the truncation function satisfies the following bounds

we have $\boldsymbol{Y}_{1}\prec\boldsymbol{Y}\prec\boldsymbol{Y}_{2}$. We further compute the expectations of $\boldsymbol{Y}_{1}$ and $\boldsymbol{Y}_{2}$ and obtain

where $\xi\sim\mathcal{N}(0,1)$. For given $\alpha_{l}$ and $\alpha_{u}$, small value of $\epsilon$ yields arbitrarily close $\beta_{1}$ and $\beta_{3}$, as well as arbitrarily close $\beta_{2}$ and $\beta_{4}$. For example, taking $\alpha_{l}=1,\alpha_{u}=5$ and $\epsilon=0.01$, we have $\beta_{1}=0.9678,\beta_{2}=0.4791,\beta_{3}=0.9678,\beta_{4}=0.4888$.

Now applying standard results on random matrices with non-isotropic sub-Gaussian rows [61, equation (5.26)] and noticing that $\boldsymbol{a}_{i}\boldsymbol{a}_{i}^{T}|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|\boldsymbol{1}_{\{\alpha_{l}(1+\epsilon)<|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|<\alpha_{u}(1-\epsilon)\}}$ can be rewritten as $\boldsymbol{b}_{i}\boldsymbol{b}_{i}^{T}$ for sub-Gaussian vector $\boldsymbol{b}_{i}:=\boldsymbol{a}_{i}\sqrt{|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|}\boldsymbol{1}_{\{\alpha_{l}(1+\epsilon)<|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|<\alpha_{u}(1-\epsilon)\}}$, one can derive

Appendix B Supporting Arguments for Section 2.2

The expectation of the loss function (2) of WF is given by as

We next show that the expectation of the loss function (3) of RWF has the following form:

where $\rho=\frac{\boldsymbol{z}^{T}\boldsymbol{x}}{\|\boldsymbol{x}\|\|\boldsymbol{z}\|}$ and $K_{0}(\cdot)$ is the modified Bessel function of the second kind.

Following , the density function of $u\cdot v$ is given by

for $|\rho|<1$. Therefore, if $|\rho|<1$, then

where the last step follows by changing variables.

If $|\rho|=1$, then $|uv|$ becomes a $\chi_{1}^{2}$ random variable, with the density

B.2 Proof of Lemma 1

We evaluate the conditional probability as follows.

where $\text{erfc}(z):=\frac{2}{\sqrt{\pi}}\int_{z}^{\infty}\exp(-t^{2})dt$, and the last inequality holds if $\|\boldsymbol{h}\|<(1-1/\sqrt{2})\|\boldsymbol{x}\|$.

Appendix C Proof of Theorem 1: Geometric Convergence of RWF

The general structure of the proof follows that for WF in and TWF in . However, the proof requires development of new bounds due to the nonsmoothness of the loss function and absolute value based measurements. On the other hand the proof is much simpler due to the lower-order loss function adopted in RWF.

The idea of the proof is to show that within the neighborhood of global optimums, RWF satisfies the Regularity Condition $\mathsf{RC}(\mu,\lambda,c)$ , i.e.,

for all $\boldsymbol{z}$ and $\boldsymbol{h}=\boldsymbol{z}-\boldsymbol{x}$ obeying $\|\boldsymbol{h}\|\leq c\|\boldsymbol{x}\|$, where $0<c<1$ is some constant. Then, as shown in , once the initialization lands into this neighborhood, geometric convergence can be guaranteed, i.e.,

for any $\boldsymbol{z}$ with $\|\boldsymbol{z}-\boldsymbol{x}\|\leq c\|\boldsymbol{x}\|$.

To show the regularity condition, we first define a set $\mathcal{S}:=\{i:1\leq i\leq m,(\boldsymbol{a}_{i}^{T}\boldsymbol{z})(\boldsymbol{a}_{i}^{T}\boldsymbol{x})<0\}$, and then derive the following bound:

The first term in (44) can be bounded using Lemma 3.1 in , which we state below.

For any $0<\epsilon<1$, if $m>c_{0}n\epsilon^{-2}$, then with probability at least $1-2\exp(-c_{1}\epsilon^{2}m)$,

The above equation can be further upper bounded by the following lemma.

For any $\epsilon>0$, if $m>c_{0}n\epsilon^{-2}\log\epsilon^{-1}$, then with probability at least $1-C\exp(-c_{1}\epsilon^{2}m)$,

Therefore, combining Lemmas 2 and 3 with (44) yields

Fix $\delta>0$, and assume $y_{i}=|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|$. Suppose that $m\geq c_{0}n$ for a certain constant $c_{0}>0$. There exist some universal constants $c,C>0$ such that with probability at least $1-C\exp(-cm)$,

Thus, applying Lemma 4 to (48), we conclude that Regularity Condition (42) holds for $\mu$ and $\lambda$ satisfying

which concludes the proof. The proofs of two major lemmas are provided in the following two subsections.

We first prove bounds for any fixed $\boldsymbol{h}\leq\frac{1}{10}\|\boldsymbol{x}\|$, and then develop a uniform bound later on. We introduce a series of auxiliary random Lipschitz functions to approximate the indicator functions. For $i=1,\ldots,m$, define

and then $\chi_{i}(t)$’s are random Lipschitz functions with Lipschitz constant $\frac{1}{\delta}$. We further have

For convenience, we denote $\gamma_{i}:=\frac{|\boldsymbol{a}_{i}^{T}\boldsymbol{h}|^{2}}{\|\boldsymbol{h}\|^{2}}\boldsymbol{1}_{\{(1-\delta)|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|^{2}<|\boldsymbol{a}_{i}^{T}\boldsymbol{h}|^{2}\}}$ and $\theta:=\|\boldsymbol{h}\|/\|\boldsymbol{x}\|$. We next estimate the expectation of $\gamma_{i}$, by conditional expectation,

where $f(\tau_{1},\tau_{2})$ is the density of two joint Gaussian random variables with correlation $\rho=\frac{\boldsymbol{h}^{T}\boldsymbol{x}}{\|\boldsymbol{h}\|\|\boldsymbol{x}\|}\neq\pm 1$. We then continue to derive

Furthermore, $\chi_{i}(|\boldsymbol{a}_{i}^{T}\boldsymbol{h}|^{2})$’s are sub-exponential with sub-exponential norm $\mathcal{O}(\|\boldsymbol{h}\|^{2})$. By the sub-exponential tail bound (Bernstein type) , we have

for some universal constant $c$, as long as $\|\boldsymbol{h}\|\leq\frac{1}{10}\|\boldsymbol{x}\|$.

holds with probability at least $1-(1+2/\epsilon)^{n}\exp(-cm\epsilon^{2})$.

Since $\chi_{i}(t)$’s are Lipschitz functions with constant $1/\delta$, we have the following bound

holds with probability at least $1-C\exp(-c_{1}m)$ as long as $m>c_{0}n$ for some constants $C,c_{0},c_{1},c_{2}>0$. Consequently, on the event that (60) holds, we have

where the last inequality is due to the Lemma 2 in .

On the intersection of events that (58) and (60) hold, we have

for all $\boldsymbol{h}$ with $\|\boldsymbol{h}\|=\frac{1}{10}\|\boldsymbol{x}\|$.

For the case when $\|\boldsymbol{h}^{\prime}\|<\frac{1}{10}\|\boldsymbol{x}\|$, $\boldsymbol{h}^{\prime}=\omega\boldsymbol{h}$ for some $\boldsymbol{h}$ satisfying $\|\boldsymbol{h}\|=\frac{1}{10}\|\boldsymbol{x}\|$ and $0<\omega<1$. By the definition of $\chi_{i}(\cdot)$, it can be verified that

Applying (61), on the same event that (58) and (60) hold, we have

for all $\|\boldsymbol{h}^{\prime}\|<\frac{1}{10}\|\boldsymbol{x}\|$. Since $\epsilon$ can be arbitrarily small, the proof is completed.

C.2 Proof of Lemma 4

Denote $v_{i}:=\boldsymbol{a}_{i}^{T}\boldsymbol{z}-|\boldsymbol{a}_{i}^{T}\boldsymbol{x}|\text{sgn}(\boldsymbol{a}_{i}^{T}\boldsymbol{z})$. Then

where $\boldsymbol{A}$ is a matrix with each row being $\boldsymbol{a}_{i}^{T}$ and $\boldsymbol{v}$ is a $m-$dimensional vector with each entry being $v_{i}$. Thus,

as long as $m\geq c_{1}n$ for some sufficiently large $c_{1}>0$, where the spectral norm bound $\|\boldsymbol{A}\|\leq\sqrt{m}(1+\delta)$ follows from [61, Theorem 5.32].

We next bound $\|\boldsymbol{v}\|$. Let $\boldsymbol{v}=\boldsymbol{v}^{(1)}+\boldsymbol{v}^{(2)}$, where $v_{i}^{(1)}=\boldsymbol{a}_{i}^{T}\boldsymbol{h}$ and $v_{i}^{(2)}=2\boldsymbol{a}_{i}^{T}\boldsymbol{x}\boldsymbol{1}_{\{(\boldsymbol{a}_{i}^{T}\boldsymbol{z})(\boldsymbol{a}_{i}^{T}\boldsymbol{x})<0\}}$. By triangle inequality, we have $\|\boldsymbol{v}\|\leq\|\boldsymbol{v}^{(1)}\|+\|\boldsymbol{v}^{(2)}\|$. Furthermore, given $m>c_{0}n$, by [9, Lemma 3.1] with probability $1-\exp(-cm)$, we have

By Lemma 3, we have with probability $1-C\exp(-c_{1}m)$

Appendix D Proof of Theorem 2: Stability to Bounded Noise

We consider the model (18) with bounded noise, i.e., $y_{i}=\left|\langle\boldsymbol{a}_{i},\boldsymbol{x}\rangle\right|+w_{i}$ for $i=1,\cdots,m$. The initialization analysis is similar to Appendix A. To analyze the gradient loop, we consider two regimes.

$\bullet$ Regime 1: $c_{4}\|\boldsymbol{z}\|\geq\|\boldsymbol{h}\|\geq c_{3}\frac{\|\boldsymbol{w}\|}{\sqrt{m}}$. In this regime, error contraction by each gradient step is given by

where the second inequality is because the spectral norm bound $\|\boldsymbol{A}\|\leq\sqrt{m}(1+\delta)$ following from [61, Theorem 5.32]. Given the regime condition $\|\boldsymbol{h}\|\geq c_{3}\frac{\|\boldsymbol{w}\|}{\sqrt{m}}$, we further have

The $\mathsf{RC}$ is guaranteed if $\mu,\lambda,\epsilon$ are chosen properly, $c_{3}$ is sufficiently large, and $s$ is sufficiently small.

for some constant $c_{5}$. As long as $\|\boldsymbol{w}\|/\sqrt{m}$ is sufficiently small, it is guaranteed that $c_{5}\frac{\|\boldsymbol{w}\|}{\sqrt{m}}\leq c_{4}\|\boldsymbol{x}\|$. If the iterate jumps out of Regime 2, it falls into Regime 1.

Appendix E Convergence of IRWF, Minibatch IRWF and Kaczmarz-PR

Since the initialization is the same as that in Algorithm 1, it suffices to show the convergence of gradient loops given that the initial point lands into the neighborhood of global minimums. To prove Theorem 3, the major step is to prove the following Proposition 2 which characterizes how the error of an estimate decays upon one iteration of Algorithm 2. Once Proposition 2 is established, we take expectation on both sides of (80) with respect to $i_{t-1}$, and apply Proposition 2 one more time to obtain

Continuing this process until the initialization point $\boldsymbol{z}^{(0)}$ yields Theorem 3. We next focus on proving Proposition 2 stated bellow

Assume the measurement vectors are independent and each $\boldsymbol{a}_{i}\sim\mathcal{N}(0,\boldsymbol{I})$. There exist some universal constants $0<\rho,\rho_{0}<1$ and $c_{0},c_{1},c_{2}>0$ such that if $m\geq c_{0}n$ and $\mu=\frac{\rho_{0}}{n}$, then with probability at least $1-c_{1}\exp(-c_{2}m)$, we have

to hold for all $\boldsymbol{z}^{(t)}$ satisfying $\frac{\text{dist}(\boldsymbol{z}^{(t)},\boldsymbol{x})}{\|\boldsymbol{z}\|}\leq\frac{1}{10}$.

Without loss of generality, we assume $\boldsymbol{z}^{(t)}$ is in the neighborhood of $\boldsymbol{x}$ (otherwise it is in the neighborhood of $-\boldsymbol{x}$). Let $\boldsymbol{h}=\boldsymbol{z}^{(t)}-\boldsymbol{x}$.

We follow the notations in Appendix C and let $\mathcal{S}=\{i:(\boldsymbol{a}_{i}^{T}\boldsymbol{x})(\boldsymbol{a}_{i}^{T}\boldsymbol{z}^{(t)})<0\}$. Then we have

where $(a)$ is due to the fact $i_{t}$ is sampled uniformly at random from $\{1,2,\cdots,m\}$. Applying Lemma 2, we have that if $m\geq c_{0}\epsilon^{-2}n$, then with probability $1-2\exp(-c_{1}m\epsilon^{2})$

holds for all vectors $\boldsymbol{h}$. Furthermore, by Lemma 3, we have that with probability $1-C\exp(-c_{1}m\epsilon^{2})$,

holds for all $\boldsymbol{h}$ satisfying $\|\boldsymbol{h}\|/\|\boldsymbol{x}\|\leq\frac{1}{10}$.

By choosing the step size $\mu\leq\frac{0.24}{n}$, the proposition is proved. ∎

E.2 Proof of Theorem 4

As argued in Appendix E.1, it suffices to show that one iteration of Algorithm 2 satisfies the following property.

Assume the measurement vectors are independent and each $\boldsymbol{a}_{i}\sim\mathcal{N}(0,\boldsymbol{I})$. There exist some universal constants $0<\rho,\rho_{0}<1$ and $c_{0},c_{1},c_{2}>0$ such that if $m\geq c_{0}n$ and $\mu=\rho_{0}/n$ for the update rule (22), then with probability at least $1-c_{1}\exp(-c_{2}m)$, we have

to hold for all $\boldsymbol{z}^{(t)}$ satisfying $\frac{\text{dist}(\boldsymbol{z}^{(t)},\boldsymbol{x})}{\|\boldsymbol{z}\|}\leq\frac{1}{10}$.

Without loss of generality, we assume $\boldsymbol{z}^{(t)}$ is in the neighborhood of $\boldsymbol{x}$ (otherwise it is in the neighborhood of $-\boldsymbol{x}$). Let $\boldsymbol{h}=\boldsymbol{z}^{(t)}-\boldsymbol{x}$.

We follow the notations in Appendix C and let $\mathcal{S}=\{i:(\boldsymbol{a}_{i}^{T}\boldsymbol{x})(\boldsymbol{a}_{i}^{T}\boldsymbol{z}^{(t)})<0\}$. Then we have

where $(a)$ is due to the fact that $\Gamma_{t}$ is uniformly chosen from all subsets of $\{1,2,\ldots,m\}$ with cardinality $k$.

Following the arguments for obtaining (82), we have

By choosing the step size $\mu\leq\frac{0.24}{n}$, the proposition is proved. ∎

E.3 Proof of Theorem 5

Without loss of generality, we assume $\boldsymbol{z}^{(t)}$ is in the neighborhood of $\boldsymbol{x}$ (otherwise it is in the neighborhood of $-\boldsymbol{x}$). Let $\boldsymbol{h}=\boldsymbol{z}^{(t)}-\boldsymbol{x}$.

We follow the notations in Appendix C and let $\mathcal{S}=\{i:(\boldsymbol{a}_{i}^{T}\boldsymbol{x})(\boldsymbol{a}_{i}^{T}\boldsymbol{z}^{(t)})<0\}$. Then we have

Moreover, $\{\sqrt{n}\frac{\boldsymbol{a}_{i}}{\|\boldsymbol{a}_{i}\|}\}_{i=1}^{m}$ are sub-Gaussian and the sub-Gaussian norm is bounded by an absolute constant. Thus, we have that if $m\geq c_{0}\epsilon^{-2}n$, then with probability $1-2\exp(-c_{1}m\epsilon^{2})$,

holds for all vectors $\boldsymbol{h}$. By Lemma 3, we have that with probability $1-C\exp(-c_{1}m\epsilon^{2})$

holds for all $\boldsymbol{h}$ satisfying $\|\boldsymbol{h}\|/\|\boldsymbol{x}\|\leq\frac{1}{10}$.

References