A Fast and Provable Algorithm for Sparse Phase Retrieval
Authors: Jian-Feng CAI
Yu Long
Ruixue WEN
Jiaxi Ying
Venue: ICLR 2024 poster (2024)
Keywords: Sparse Phase Retrieval
Quadratic Convergence
Newton-type Method
Hard Thresholding
Nonconvex Optimization
Primary area: optimization

Abstract: We study the sparse phase retrieval problem, which seeks to recover a sparse signal from a limited set of magnitude-only measurements. In contrast to prevalent sparse phase retrieval algorithms that primarily use first-order methods, we propose an innovative second-order algorithm that employs a Newton-type method with hard thresholding. This algorithm overcomes the linear convergence limitations of first-order methods while preserving their hallmark per-iteration computational efficiency. We provide theoretical guarantees that our algorithm converges to the $s$-sparse ground truth signal $\boldsymbol{x}^{\natural} \in \mathbb{R}^n$ (up to a global sign) at a quadratic convergence rate after at most $O(\log (\Vert\boldsymbol{x}^{\natural} \Vert /x_{\min}^{\natural}))$ iterations, using $\Omega(s^2\log n)$ Gaussian random samples. Numerical experiments show that our algorithm achieves a significantly faster convergence rate than state-of-the-art methods.

Published as a conference paper at ICLR 2024
A FAST AND PROVABLE ALGORITHM FOR SPARSE PHASE RETRIEVAL
Jian-Feng Cai Yu Long" Ruixue Wen', Jiaxi Ying
1 Hong Kong University of Science and Technology
2 HKUSS Shenzhen-Hore Kong Collaborative Innovation Research Institute
3 Guangxi University
ABSTRACT
We study the sparse phase retrieval problem, which seeks to recover a sparse signal from limited set magnitude-only measurements. In contrast prevalent sparse phase retrieval algorithms that primarily use first-order methods, we propose an innovative second-order algorithm that employs a Newton-type method with hard thresholding. This algorithm overcomes the linear convergence limitations of first-order methods while preserving their hallmark per-iteration computational efficiency. We provide theoretical guarantees that our algorithm converges to the s-sparse ground truth signal æ² E R." (up to global sign) at a quadratic convergence rate after at most ||= iterations, using s2ss log n) Gaussian random samples. Numerical experiments show that our algorithm achieves significantly faster convergence rate than ttate-of-the-ee the-art methods.
INTRODUCTION
We study the phase retrieval problem, which involves reconstructing an n-dimensional signal ø using its magnitude-only measurements;
y/- Ka,,, x 1 i 1,2.......
(1)
where each y; represents a measurement, a; denotes a sensing vector, æt is the unknown signal to be recovered, and m the total number measurements. The phase retrieval problem arises in various applications, including diffraction imaging (Maiden & Rodenburg, 2009), X-ray crystallography (Miao et al., 2008), and optics (Shechtman et al., 2015).
Although the phase retrieval problem ill-posed and even NP-hard (Fickus et al., 2014), various algorithms can recover target signals. These are broadly categorized into convex and nonconvex approaches. Convex methods, such PhaseLift (Candes et al., 2015; 2013), PhaseCut (Waldspurger et al., 2015), and PhaseMax (Goldstein & Studer, 2018; Hand & Voroninski, 2018), offer optimal sample complexity but are computationally challenging in high-dimensional cases. To improve computational efficiency, nonconvex approaches explored such as alternating minimization (Netrapalli et al., 2013; Cai et al., 2022b), Wirtinger flow (Candes et al., 2015), truncated amplitude flow (Wang et al., 2017a), Riemannian optimization (Cai & Wei, 2023), Gauss-Newton (Gao & Xu, 2017; Ma et al., 2018), Kaczmarz (Wei, 2015; Chi & Lu, 2016), and unregularized gradient descent (Ma et al., 2020). Despite the nonconvex nature of its objective function, the global geometric landscape lacks spurious local minima (Sun et al., 2018; Li et al., 2019; Cai et al., 2023b), allowing algorithms with random initialization work effectively (Chen et al., 2019; Waldspurger, 2018).
The nonconvex approaches previously mentioned can guarantee successful recovery of the ground truth (up to global phase) using m 22 log" n) measurements, where a 2 0. This complexity is nearly optimal, as the phase retrieval problem requires m 2 2n for real signals and m 2 4n - 4 for complex signals (Conca al., 2015). However, in practical situations, especially in high-dimensional cases, the number of available measurements is often less than the signal dimension (i.e., m < n), leading to a need for further reduction in sample complexity.
"Corresponding authors.

---

Published as a conference paper at ICLR 2024
In this paper, we focus on the sparse phase retrieval problem, which aims to recover sparse signal from a limited number of phaseless measurements, It has been established that the minimal sample complexity required to ensure s-sparse phase retrievability in the real case is only 2s for generic sensing vectors (Wang & Xu, 2014). Several algorithms have been proposed to address the sparse phase retrieval problem (Ohlsson et al., 2012; Cai et al., 2016; Wang et al., 2017b; Jagatap & Hegde, 2019; Cai et al., 2022c). These approaches have been demonstrated to effectively reconstruct the ground truth using (2 log n) Gaussian measurements. While this complexity is not optimal, is significantly smaller than that in general phase retrieval.
1.1 CONTRIBUTIONS
Existing algorithms for sparse phase retrieval primarily employ first-order methods with linear convergence. Recent work (Cai et al., 2022c) presented second-order method, while fails to obtain quadratic convergence. The main contributions of this paper can be summarized in three points:
1. We propose second-order algorithm for sparse phase retrieval that maintains the same per- iteration computational complexity popular first-order methods. Our algorithm enhances convergence by integrating second-orde derivative information from intensity- based empirical loss into the search direction. To ensure computational efficiency, the Newton step is applied only a subset of variables identified through the amplitude-bast bbased empirical loss.
2.
We establish a non-asymptotic quadratic convergence rate for our proposed algorithm and provide the iteration complexity. Specifically, we prove that the algorithm converges to the ground truth (up to global sign) at a quadratic rate after at most ([ =0 /// (l/rr())))) )) iterations, with m = O(s² log n) measurements. To the best of our knowledge, this is the first algorithm to establish a quadratic convergence rate for sparse phase retrieval.
3. Numerical experiments demonstrate that the proposed algorithm achieves significantly faster convergence rate in comparison to state-of-the-art methods. The experiments also reveal that our algorithm attains higher success rates in the exact recovery of signals from noise-free measurements and offers improved signal reconstruction in the presence of noise.'
Notation: ||æ||| denotes the number of nonzero entries of æ. and ||>|| denotes the f--nomm. For
matrix A E Rmxx || is the spectral norm of A. For any Q1 2 and 92 2 1, | denotes the induced operator norm from the Banach space (RR" |1-Iloz) ) to (I ||l|-1 lg.)) ). Amin (A) and Amax (A) denote the smallest and largest eigenvalues of A. |S] denotes the number of elements in S. a @ b denotes the entrywise product of a and b. We write f(n = Ogg(n) to indicate (n) I cig(n) for some constant ei > 0, and (n) (2g(n)) to denote ((n) 2 c2g(n) for some constant C2 > 0. For ar a² RR", the distance between æ: and 2 is defined as dist((x, .²) = min {||æ: I || a= []]
Lmin denotes the smallest nonzero entry in magnitude of æ².
2
PROBLEM FORMULATION AND RELATED WORK
We first present the problem formulation for sparse phase retrieval, and then review related work and provide a comparative overview of state-of-the-art algorithms and our proposed algorithm.
2.1
PROBLEM FORMULATION
The standard sparse phase retrieval problem can concisely expressed as finding æ that satisfies
( x 2 = y Vi li=.....,m and x llo S s,
where {a,)m, and {y: represent known and fixed sensing vectors and phaseless measurements, respectively. Each y/1 (a, 12 with æ² the ground truth (|=²1lo S s). While sparsity level is assumed known priori for theoretical analysis, To address Problem (2), various problem reformulations have been explored. Convex formulations, such as the l reguullard PhaseLift method (Ohlsson et al., 2012), often use the lifting technique and solve the problem in the n X matrix space, resulting in high computational costs.
experiments will explore cases with unknown
Our codes are available http://giitaa.ccmm/yyyin//uune
2

---

Published as a conference paper at ICLR 2024
Table Overview of per-iteration computational cost, iteration complexity, and loss function types for ThWF (Cai et al., 2016), SPARTA (Wang et al., 2017b), CoPRAM (Jagatap & Hegde, 2019), HTP (Cai et al., 2022c), and our proposed algorithm. Here, x' represents the ground truth with dimension and sparsity and zinin denotes the smallest nonzero entry in magnitude of x².
Methods
Per-iteration computation
Iteration complexity
Loss function
ThWF
O(n logn)
O(log(1/e))
fx
SPARTA
( log m)
O[log(1/e))
fa(w)
CoPRAM
( logn) m)
O[log(L/e))
fa(x)
HTP
O((n + s²)s² logn)
O(log(log(n")
fa(x)
Proposed
O((n + 5
O(log(log(1/e) log[ (/)mm))))
fx()), fa(x)
To enhance computational efficiicc nonconvex approaches (Cai et al., 2016; Wang et al., 2017b; Cai et al., 2022c; Soltanolkotabi, 2019) are explored, which can be formulated as: minimize f(z), subject to ||2]||0 S s.
Both the loss function f(x)) and the lo-norm constraint in Problem (3) are nonconvex, making it challenging to solve. Two prevalent loss functions are investigated: intensity- empirical loss
m
fr(x) = 4m
(Iaa; x)1² yi)
and amplitude-basee based empirical loss
(a
1
Σ (I az)| 21)²
2m
where 2; VUi. 1,..., m. The intensity- based loss fr( is smooth, while the amplitude-base -base
loss fa(x) is non-smooth because of the modulus.
2.2 RELATED WORK
Existing nonconvex sparse phase retrieval algorithms can be broadly classified into two categories: gradient projection methods and alternating minimization methods. Gradient projection methods, such as ThWF (Cai et al., 2016) and SPARTA (Wang et al., 2017b), employ thresholded gradient descent and iterative hard thresholding, respectively. On the other hand, alternating minimization methods, including CoPRAM (Jagatap & Hegde, 2019) and HTP (Cai et al., 2022c), alternate between updating the signal and phase. When updating the signal, formulated squares problem, CoPRAM leverages the cosamp method (Needell & Tropp, 2009), while HTP applies the hard thresholding pursuit algorithm (Foucart, 2011).
sparsity-constraing least
Contrary to previously discussed algorithms, our algorithm is rooted in a Newton-type method with hard thresholding and incorporates second-order information to accelerate convergence. Unlike alternating minimization methods, our algorithm eliminates the need to update the signal and phase separately. A sample complexity of 12s² log under Gaussian measurements is required for successful recovery across all discussed algorithms. ThWF employs intensity- loss as the objective function, while SPARTA, CoPRAM, and HTP utilize an amplitude-bare loss. Our algorithm, distinctively, adopts both loss types: It uses intensity based loss as the objective function and amplitude-based loss to determine the support for the Newton update.
Discussions: Table provides comparative overview of various algorithms. ThWF, SPARTA, and CoPRAM, as first-order methods, exhibit linear convergence. In contrast, our algorithm achieves quadratic convergence rate. Although HTP is a second-order method that converges in a finite number of iterations, it fails to establish quadratic convergence and has higher iteration complexity than our algorithm. Empirical evidence further shows that our algorithm converges faster than HTP. This can attributed to our algorithm's 's more effective exploitation of the second-order information of the objective function when constructing the search direction.
3

---

Published as a conference paper at ICLR 2024
3 MAIN RESULTS
In this section, we present our proposed algorithm for sparse phase retrieval. Generally, nonconvex methods comprise two stages: initialization and refinement. The first stage generates an initial guess close the target signal, while the second stage refines the initial guess. Our proposed algorithm adheres to this two-stage strategy. In the first stage, we employ existing effective method to generate an initial point. Our primary focus is on the second stage, wherein we propose an efficient second-order algorithm refine the initial guess.
3.1 PROPOSED ALGORITHM
We introduce our proposed second-order algorithm for sparse phase retrieval, which adopts a dual loss strategy: It employs the intensitybbased loss defined in (4) for the objective function and uses the amplitude-based loss from (5) to identify the support for the Newton update. In alignment with established Newton-type methods for sparse optimization (Cai et al., 2023a; Zhou et al., 2021; Meng Zhao, 2020), our algorithm involves two primary steps: (1) identifying the sets of free and fixed variables, and (2) computing the search direction.
3.1.1 IDENTIFYING FREE AND FIXED VARIABLES
In Newton-type methods, computing the Newton direction each iteration typically requires solving a linear system, process that often incurs significant computational cost of O(n"). Here, n denotes the dimension of the problem space, and w represents the matrix multiplication constant, which could be less than if fast matrix multiplication employed. Unfortunately, this complexity can still render the algorithm impractical for high-dimensional scenarios where n is large.
To address this challenge, we categorize variables into two groups at each iteration: free and fixed, and update them separately. The free variables, which consist of at most s variables, are updated according to the (approximate) Newton direction. while the fixed variables are set to zero. This strategy substantially cuts the computational cost from O(n²) to O(s*) by only requiring the solution of a linear system size s x 8, and ensures s-sparsity at each iteration.
We identify free variables using one-step iterative hard thresholding (IHT) of the loss fa(x) in (5):
Sk+1 supp (H, fax"))
where vj the stepsize, and H. denotes the s-sparse hard-thresholding operator, defined as follows:
H.(w) == arg min | - w| --wl| subject to ||0||| Ss,
(6)
This operator picks the largest magnitude entries of the input vector and sets the rest to zero. Therefore, there at most free variables. Since fa is non-smooth, we adopt the generalized gradient (Zhang et al., 2017) as fA. The set of fixed variables is the complement of S We only update free variables along the approximate Newton direction and set others to zero.
While the objective function of algorithm is the intensity- loss fr(w), it's noteworthy that we use the amplitude-base loss fA (æ:) in place of fx( when identifying free variables. Empirical evidence shows that this approach yields faster convergence compared to using fz
3.1.2 COMPUTING SEARCH DIRECTION
We update the free variables in Sk+1 by solving the following support-constrainet problem:
minimize (x)),, subject to supp(æ:) C Sk+11
where (x) is an approximation of the intensity-based loss ft at , To ensure fast convergence and efficient computation, we choose (x) in (7) as the second-order Taylor expansion of fr at k
the (z)) fr(('* (V/r (xx x w²) ( x², fr -=="=") x²)).
For notational simplicity, define g [V]r(z>)]] )s.." which denotes the sub-vector of V V (x+ indexed by Sk+1. H Seta Sets (x**)) S Sepr which represents the principle sub matrix
4

---

Published as a conference paper at ICLR 2024
of the Hessian indexed by Sk+1, and H H
fi("))
14 denoting the sub-matrix
whose rows and columns are indexed by Sk-1 and S respectively. Let x denote the minimizer
of Problem (7). Following from the first-order optimality condition of Problem (7), we obtain that S and T5a41 satisfies
H Set1 (IS - 2 x'S.1)) = H SE E'S
(8)
As a result, we obtain the next iterate wht1 +1 by
a ++1 S441
x SS4t
P and
x 1
k+1 O,
(9)
where psa.. represents the approximate Newton direction over Sk+1, which can calculated by
H Se Ps I x J211 + g
(10)
where Je41 ;= S \ Sk+1 with |J++1| S 8. In contrast to Eq. (8), we replace sse with æ7 in
(10), as Jk+1 captures all nonzero elements in st as follows:
(Z
(
x s nsp
nS
1-
x Su+1
1-
Je+1
(11)
0
where operator 9 arranges all nonzero elements of vector to appear first, followed by zero elements. The first equality in (11) follows from the fact that supp(x") C Sy.
By calculating H H rather than H St as in (10), the computational cost substantially reduced from (smmn) to Os)²m). The costs for computing S....S., and solving the linear system in (10) are (s)²m) and O(s"), respectively. Thus, the cost for Step 2 is O(s²m). The cost for Step amounts to ((m)). which involves calculating fa(" Therefore, the overall cost per iteration is (n + s² )m, with on the order of 3² log n, which is generally necessary for theoretical recovery guarantee. Assuming for sparse phase retrieval would reduce to that of general methods), our algorithm's per-iteration complexity matches that of popular first-order methods O(ns logn).
O(vπ)) (otherwise, the complexity Q(s log
Algorithm 1 Proposed algorithm
Input: Data {a,, | sparsity s, initial estimate x., and stepsize 17.
1: for k == 0, 2, do
2: Identify the set of free variables Sk+1 supp(H,, - TA (x)))));;
3: Compute the search direction pls... over Sk+1 by solving (10);
4: Update k+1,
k+1 Sut1
5
and
s $++1
k+1
0.
5: end for
Output: k+1
3.1.3
LEVERAGING TWO TYPES OF LOSSES
We provide an intuitive explanation of our algorithm, with a particular focus on the utilization of two types of loss functions: the intensity-based loss and the amplitude-basee based loss.
Our algorithm employs the intensity-based loss fr as the objective function; its smoothness facilitates the computation of the Hessian and the construction of the Newton direction. While one might intuitively use the same loss function to determine variables for the Newton update, numerical experiments indicate that the amplitude-based loss fA is more effective for this purpose.
One potential reason is the superior curvature exhibited by fA around the true solution, which often behaves more similarly that of quadratic least-squares loss (Chi al., 2019; Zhang al., 2017; Wang et al., 2017a). This indicates that the gradient of / could offer more effective search direction than that of fr. Furthermore, studies (Zhang et al., 2018; Cai et al., 2022a) have demonstrated that algorithms founded on fa consistently outperform those based on fi in numerical results.
5

---

Published as a conference paper at ICLR 2024
3.2 INITIALIZATION
The nonconvex nature of phase retrieval problems often requires good initial guess to find the target signal. Spectral initialization (Candes et al., 2015) is a common approach. We adopt a sparse variant of the spectral initialization method obtain a favorable initial guess. would be intriguing to consider other well-established initialization methods, such sparse orthogonality- tattogonaityy-romoiing promoting initialization (Wang et al., (b)), diagonal thresholding initialization (Cai et al., 2016), and modified spectral initialization method (Cai et al., 2023c).
Assuming {a: m independently drawn from Gaussian distribution (O,, I,), the expectation of the matrix 1 Ee y:a.a? is M == || I I, + 2x (x+)² The leading eigenvector of M is precisely = Hence, the leading eigenvector of J E== y;a,a can be close to == (Candes et al., 2015). However, this method requires the sample complexity of at least 2(n), which is excessively high for sparse phase retrieval. Leveraging the sparsity of æ4 crucial to lower this complexity.
We adopt the sparse spectral initialization method proposed in (Jagatap & Hegde, 2019). Specifically, first collect the indices of the largest 8 values from { [ Σx 3h [a:1)) and obtain the set S. serving as estimate of the support of the true signal x Next, construct the initial guess x as follows: z is the leading eigenvector of J. Ei= Bi [alall [aal.... and zz-- O. Finally, we scale x such that ||= m E= Mc. ensuring the power of x closely aligns with the power of x
The study (Jagatap & Hegde, 2019) demonstrates that, given a sample complexity m = S2 2 log n)) the aforementioned sparse spectral initialization method can produce initial estimate x0 that falls within small-sized neighborhood around the ground truth. Specifically, it holds dist(x x²) < >||>0|| for any Y E (0, 1), with a probability of at least Sm 11.-1. Our theoretical analysis demonstrates that, once the estimate enters this neighborhood,
iterative updates in the refinement stage will consistently stay within this neighborhood and be attracted towards the ground truth, exhibiting at least linear convergence and achieving quadratic convergence after iterations.
3.3 THEORETICAL RESULTS
Given the nonconvex nature of the objective function and constraint set in the sparse phase retrieval problem, a theoretical analysis crucial to guarantee the convergence of our algorithm. In this subsection, we provide comprehensive convergence analysis for both noise-free and noisy cases.
3.3.1 NOISE-FREE CASE
We begin by the noise-free free case, in which each measurement y, = ((a,, )² Starting with an initial guess obtained via the sparse spectral initialization method, the following theorem shows that our algorithm exhibits quadratic convergence rate after at most Olog( [|=* (im))) iterations.
Theorem 3.1. Let {a, ) be i.i.d. random vectors distributed as .N(0, In), and x² E R" be any signal with ||=e" llo < s. Let {=*] be the sequence generated by Algoritium I with I,...,m, and the initial guess x² generated by the sparse spectral initialization method mentioned earlier. There exist positive constants p, m., 12, C1, C2, Ca, C4, C₁, such that if the stepsize 1) E [1,, and m 2 Cis² log then with probability at least (C2K + Cz)m= the sequence {x* J&21 converges to the ground truth x² at quadratic rate after at most O (log([]=>\//==_.)
the input measurements y/s (a, {2,
iterations, i.e.
dist(x**1 x²4) S dist x²)), Vkz K,
where K S C4 log (1 Cs. and wi... is the smallest nonzero entry in magnitude of wh
The proof of Theorem 3.1 available in Appendix B.3.1. Theorem 3.1 establishes the non-asymptotic quadratic convergence rate of our algorithm it converges to the ground truth, leading to an iteration complexity of (loog(lo(()/e)) + log(I] (1==+ for achieving an e-accurate solution. While superlinear convergence for Newton-type methods has been widely recognized in literature, typically holds only asymptotically, Consequently, the actual iteration complexity remains undefined, highlighting the importance of establishing a non-asymptotic superlinear convergence rate. For the first K iterations, our algorithm exhibits at least linear convergence.

---

Published as a conference paper at ICLR 2024
3.3.2 NOISY CASE
We demonstrate the robustness of our algorithm under noisy measurement conditions. Building upon (Cai et al., 2016; Chen & Candès, 2017); we assume that the noisy measurements are given by:
y: = x + G.. for i= 1....,..
where represents vector of stochastic noise that is independent of {a. Throughout this paper, we ussume, without loss of generality, that the expected value of e is O.
Theorem 3.2. Ler {a, [a.);=1111 be i.i.d. random vectors distributed as (N,,, and x² E R be any signal with ||==1lo $ s. Let {ø"})>> JR21 be the sequence generated by Algorithm I with noisy input y4 Ka, 2²)12 + Ee,
L...., m. There exists positive constants n, 72, Co Cr, Cg. and E (0, 1/8), such that if the stepsize E M2], m 2 Ces² log n. and the initial guess x" obeys dist(x", a² ) S >|1=01| with ||00||o S 8, then with probability at least (OyK' + Ca)m m1
dist(x**+ x²)) S p dist(æ", (²") + vllell,
AO k S K'
where o' (0, V E (0, 1). and K' is a positive integer.
The proof of Theorem 3.2 is provided in Appendix B.3.2. Theorem 3.2 validates the robustness of our algorithm, demonstrating its ability to effectively recover the signal from noisy measurements, It reveals a linear convergence rate for our algorithm in the presence of noise. However, our algorithm incorporates the second-order information when determining the search direction, which always leads to faster convergence than first-order algorithms. In addition, Theorem 3.2 does not suggest that the established inequality fails to hold when k 2 K"
4 EXPERIMENTAL RESULTS
In this section, we present series of numerical experiments designed to validate the efficiency and accuracy of our proposed algorithm. All experiments were conducted on a 2 GHz Intel Core i5 processor with 16 GB of RAM. and all compared methods were implemented using MATLAB.
Unless explicitly specified, the sensing vectors {a: were generated by the standard Gaussian distribution. The true signal æ² has nonzero entries, where the support is selected uniformly from all subsets of [x] with cardinality s, and their values are independently generated from the standard Gaussian distribution (0, 1). In the of noisy measurements, have:
y/ = (a,, >² des, for i 1=....,mm. ,,,
(12)
where {s:)m follow i.i.d standard Gaussian distribution, and > 0 determines the noise level.
We compare the convergence speed of our algorithm with state-of-the-art methods, including ThWF (Cai et al., 2016), SPARTA (Wang et al., 2017b), CoPRAM (Jagatap & Hegde, 2019), and HTP 0.7 for ThWF; 7 0.5, μ 1 and |Z|= = [m/6] for SPARTA; !) 0.95 for both HTP and our algorithm. The maximum number of iterations for each algorithm is 1000. The Relative Error (RE) between the estimated signal & and the ground truth d defined as RE == x
(Cai et al., 2022c). We fine-tune the parameters and set:
di A recovery is deemed successful if RE 10 10-3
Figure compares the number of iterations required for convergence across various algorithms. We 5000, and the
set the number of measurements m 3000, the dimension of the signal sparsity levels to 80 and 100. We consider both noise-free and noisy measurements with a noise level of a = 0.03. We observe that all algorithms perform well under both noise-free and noisy conditions: however, our algorithm converges with significantly fewer iterations.
Table 2 presents a comparison of the convergence running times for various algorithms, corresponding to the experiments depicted in Figure 1. For noise-free measurements, all algorithms are set to terminate when the iterate satisfies the following condition: dist( _x²') 10 10-3, which indicates |
a successful recovery. In the case of noisy measurements, the termination criterion is set as d 10-3 As evidenced by the results in Table 2, our algorithm consistently outperforms state-of-the-arr methods terms of running time, for both noise-free and noisy cases, highlighting its superior efficiency for sparse phase retrieval applications.
|
7

---

Published as a conference paper at ICLR 2024
16
Saura
104
Bekinve
5.
Baratons
Barutoms
(a) Noise free, sparsity 80
(b) Noise free, sparsity
100
16
CAPRAN
16
10
107
10
J00
10
25
100
Iasalams
(c) Noise level σ = 0.03, sparsity 8 80
(d) Noise level 0 0.03, sparsity s 100
Figure Relative error versus iterations for various algorithms, with fixed signal dimension n 5000 and sample size 3000. The results represent the average of 100 independent trial runs.
Table 2: Comparison of running times (in seconds) for various algorithms in the recovery of signals with sparsity levels of 80 and 100 for both noise-free and noisy scenarios.
Methods
ThWF
SPARTA
CoPRAM
HTP
Proposed
Sparsity 80
0.3630
1.0059
0.9762
0.0813
0.0530
Noise (a
Sparsity 100
0.6262
1.2966
3.3326
0.2212
0.1024
Sparsity 80
0.2820
1.1082
1.3426
0.1134
0.0803
Noisy (a 0.03)
Sparsity 100
0.4039
1.6368
4.1006
0.2213
0.1187
Probabiliy success
0.6
0.6
Probabilly
0.4
COPRAM
0.4
COPRAM
HTP
HTP
ThWF
ThWF
PATTA
SparTA
Proposed
Procosed
500
1000
1500
2000
2500
3000
500
1000
1500
2000
2500 3000
Measurements
(a) 8
25
(b) 8 = 50
Figure 2: Phase transition performance of various algorithms for signals of dimension 3000 with
sparsity levels 25 and 50. The results represent the average of 200 independent trial runs. The parameter settings for ThWF, SPARTA, CoPRAM, and HTP in experiments consistent with those in the studies by Jagatap & Hegde (2019); Cai et al. (2022c).
Figure depicts the phase transitions different algorithms, with the true signal dimension fixed 3000 and sparsity levels set to 25 and 50. The phase transition graph generated by evaluating the successful recovery rate of each algorithm over 200 independent trial runs. Figure shows that the probability of successful recovery for each algorithm transitions from zero to one as sample size m increases. Furthermore, our algorithm consistently outperforms ofate-of-the-art methods, achieving higher successful recovery rate across various measurement counts.
8

---

Published as a conference paper at ICLR 2024
In practical upplications, natural signals may not be inherently sparse; however, their wavelet coefficients often exhibit sparsity. Figure 3 illustrates the reconstruction performance of signal from noisy phaseless measurements, where the true signal, with dimension of 30,000, exhibits sparsity and contains 208 nonzero entries under the wavelet transform, using 20,000 samples. The sampling matrix A E 220,000030,00 is constructed from random Gaussian matrix and an inverse wavelet transform generated using four levels of Daubechies wavelet. The noise level is set to a = 0.03.
To evaluate recovery accuracy, we use the Peak Signal-to-Noise Ratio (PSNR), defined as PSNR V2 10 log MSE. where V represents the maximum fluctuation in the ground truth signal, and MSE denotes the mean squared error of the reconstruction. higher PSNR value generally indicates better reconstruction quality. As depicted in Figure 3, our proposed algorithm outperforms cate-of-the-an art methods in terms of both reconstruction time and PSNR. It achieves higher PSNR while requiring considerably less time for reconstruction. In the experiments, the sparsity level assumed to be unknown, and the hard thresholding sparsity level is set to 300 for various algorithms.
(a) Original signal
ThWF: PSNR = 71, Time(s) 11.586
SPARTA: PSNR = 66, Time(s) 38.379
(c)
CoPRAM: PSNR 61, (d) Time(s) = 117.685
HTP: PSNR 62, (c) Time(s) = 6.689
Proposed: PSNR 78,
(f)
Time(s) 2.629
Figure 3: Reconstruction of the signal with dimension of 30,000 from noisy phaseless measurements by various algorithms. Time(s) is the running time in seconds.
5 CONCLUSIONS AND DISCUSSIONS
In this paper, have introduced an efficient second-order algorithm for sparse phase retrieval. Our algorithm attains non-asymptotic quadratic convergence rate while maintaining the same per- iteration computational complexity as popular first-order methods, which exhibit linear convergence limitations. Empirical results have demonstrated a significant improvement in the convergence rate of our algorithm. Furthermore, experiments have revealed that our algorithm excels in attaining higher success rate for exact signal recovery.
Finally, we discuss the limitations of our paper, which also serve as potential avenues for future research. Both algorithm and state-of-the-aar methods share the same sample complexity of [(s2logn) n) for successful recovery; however, our algorithm requires this complexity in both the initialization and refinement stages, while state-of-the-arr methods require ( log for initialization and O(s log n//) for refinement. Investigating ways achieve tighter complexity in our algorithm's 's refinement stage is a worthwhile pursuit for future studies.
Currently, the initialization methods for spanse phase retrieval exhibit a sub-optimal sample complexity of 9(s² logn). n). A key challenge involves reducing its quadratic dependence on s. Recent study (Jagatap & Hegde, 2019) attained complexity of ((s log nt), closer to the information-theowd limit, but relied on the strong assumption of the power law decay for sparse signals. Developing an initialization method that offers optimal sample complexity for a broader range of sparse signals an engaging direction for future research.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
This work was supported by the Hong Kong Research Grant Council GRFs 16310620, 16306821, and 16307023, and Hong Kong Innovation and Technology Fund MHP//09/20, and the Project of Hetao Shenzhen-Howe Kong Science and Technology Innovation Cooperation Zone under Grant We would also like to thank the anonymous reviewers for their valuable feedback on the manuscript.
REFERENCES
Thomas Blumensath and Mike E. Davies. Iterative hard thresholding for compressed sensing. Applied and Computational Harmonic Analysis, 2763):265-274, 2009.
Jian-Feng Cai and Ke Wei. Solving systems of phaseless equations via Riemannian optimization with optimal sampling complexity. Journal of Computational Mathematics, 2023.
Fang Cai, Meng Huang. Dong Li, and Yang Wang. Solving phase retrieval with random initial guess nearly as good as by spectral initialization. Applied and Computational Harmonic Analysis, 58:60-84, 84, 2022a.
Jian-Feng Cai, Yuling Jiao, Xiliang Lu, and Juntao You. Sample-effciient sparse phase retrieval via stochastic alternating minimization. IEEE Transactions Signal Processing, 70:4951-4966. 2022b.
Hlan-Feng Cai, Jingzhi Li, Xiliang Lu, and Juntao You. Sparse signal recovery from phaseless measurements via hard thresholding pursuit. Applied and Computational Harmonic Analysis, 56: 367-390, 2022c.
Jian-Feng Cai, José Vinícius de Miranda Cardoso, Daniel P. Palomar, and Jiaxi Ying. Fast projected Newton-like method for precision matrix estimation under total positivity. In Advances Neural Information Processing Systems, volume 36, pp. 73348-73370, 2023a.
Jian-Feng Cai, Meng Huang, Dong Li, and Yang Wang. Nearly optimal bounds for the global geometric landscape of phase retrieval. Inverse Problems, 39(7)077001,, 2023b.
Jian-Feng Cai, Jingyang Li, and Juntao You. Provable sample-efficient sparse phase retrieval initialized by truncated power method. Inverse Problems, 39(7):075008, jun 2023c.
T. Tony Cai, Xiaodong Li, and Zongming Ma. Optimal rates of convergence for noisy sparse phase retrieval via thresholded Wirtinger flow. The Annals of Statistics, 44(5):2221-225 2016.
Emmanuel Candes and Terence Tao. Decoding by linear programming. IEEE Transactions on Information Theory, (1(12):2203-42) 2005.
Emmanuel Candes, Thomas Strohmer, and Vladislav Voroninski. Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming. Communications Pure and Applied Mathematics, 66(8): 12744124,, 2013.
Emmanuel J Candes, Xiaodong Li, and Mahdi Soltanolkotabi. Phase retrieval via Wirtinger flow: Theory and algorithms. IEEE Transactions on Information Theory, 6185) 2015.
Yuxin Chen and Emmanuel J Candès. Solving random quadratic systems of equations is nearly as easy as solving linear systems. Communications on Pure and Applied Mathematics, 70(5): 822-883, 2017.
Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma. Gradient descent with random initialization: Fast global convergence for nonconvex phase retrieval. Mathematical Programming, 176(1):5-37, 2019.
Yuejie Chi and Yue Lu. Kaczmarz method for solving quadratic equations. IEEE Signal Processing Letters, 23(9):1 188-1187, 2016.
10

---

Published as a conference paper at ICLR 2024
Yuejie Chi, Yue M Lu, and Yuxin Chen. Nonconvex optimization meets low-rank matrix factorization: An overview. IEEE Transactions on Signal Processing, 67(20))5239-522) 2019.
Aldo Conca, Dan Edidin, Milena Hering, and Cynthia Vinzant. An algebraic characterization of injectivity in phase retrieval. Applied and Computational Harmonic Analysis, 38(2)334-336, 2015.
Matthew Fickus, Dustin Mixon, Aaron Nelson, and Yang Wang. Phase retrieval from very few measurements. Linear Algebra and its Applications, 449:475-499, 2014.
Simon Foucart. Hard thresholding pursuit: An algorithm for compressive sensing. SIAM Journal on Numerical Analysis, 49(6):2543-2566 2011.
Simon Foucart and Holger Rauhut. An invitation compressive sensing. In A Mathernatical Introduction to Compressive Sensing, pp. 1-39. 2013.
Bing Gao and Zhiqiang Xu. Phaseless recovery using the Gauss-Newton method. IEEE Transactions on Signal Processing, 65(22):5885-5886, 2017.
Tom Goldstein and Christoph Studer. Phasemax: Convex phase retrieval via basis pursuit. IEEE Transactions on Information Theory, 64(4))2675-2689, 2018.
Paul Hand and Vladislav Voroninski. An elementary proof of convex phase retrieval in the natural parameter space via the linear program phasemax. Communications in Mathematical Sciences, 16 (7):2047-2051, 2018.
Gauri Jagatap and Chinmay Hegde. Sample-efficient algorithms for recovering structured signals from magnitude-only measurements. IEEE Transactions Information Theory, 15(7):443444555 2019.
Zhenzhen Li, Jian-Feng Cai, and Ke Wei. Toward the optimal construction a loss function without spurious local minima for solving quadratic equations. IEEE Transactions Information Theory, 66(5):3242-3266, 2019.
Chao Ma, Xin Liu, and Zaiwen Wen. Globally convergent Gevvenberg-Marquuatee method for phase retrieval. IEEE Transactions on Information Theory, 65(4))2343-234) 2018.
Cong Ma, Kaizheng Wang. Yuejie Chi, and Yuxin Chen. Implicit regularization in nonconvex statistical estimation: Gradient descent converges linearly for phase retrieval, matrix completion, and blind deconvolution. Foundations of Computational Mathematics, 20:451-632, 2020.
Andrew M Maiden and John M Rodenburg. An improved ptychographical phase retrieval algorithm for diffractive imaging. Ultraiccroscopy, 109(10): 1262, 2009.
Nan Meng and Yun-Bin Zhao. Newton-step-based hard thresholding algorithms for sparse signal recovery. IEEE Transactions Signal Processing, 68:6594-6606. 2020.
Jianwei Miao, Tetsuya Ishikawa, Qun Shen, and Thomas Eareest. Extending X-ray crystallography to allow the imaging of noncrystalline materials, cells, and single protein complexes. Annual Review of Physical Chemistry, 59 59(1)3877410, 2008
Deanna Needell and Joel A Tropp. Cosamp: Iterative signal recovery from incomplete and inaccurate samples. Applied and Computational Harmonic Analysis, 26(3):301-321, 2009,
Praneeth Netrapalli, Prateek Jain, and Sujay Sanghavi. Phase retrieval using alternating minimization. Advances in Neural Information Processing Systems, 26, 2013.
Henrik Ohlsson, Allen Yang, Roy Dong, and Shankar Sastry. CPRL-An extension of compressive sensing to the phase retrieval problem. Advances Neural Information Processing Systems, 25, 2012.
11

---

Published as a conference paper at ICLR 2024
Yoav Shechtman, Yonina C Eldar, Oren Cohen, Henry Nicholas Chapman, Jianwei Miao, and Mordechai Segev. Phase retrieval with application to optical imaging: A contemporary overview. IEFE Signal Processing Magazine, 32(3):87-109, 2015.
Mahdi Soltanolkotabi. Structured signal recovery from quadratic measurements: Breaking sample complexity barriers via nonconvex optimization. IEEE Transactions on Information Theory, 65(4): 2374-2400, 2019.
Ju Sun, Qing Qu, and John Wright. A geometric analysis of phase retrieval. Foundations of Computational Mathematics, 18(5):1 2018.
Irene Waldspurger. Phase retrieval with random Gaussian sensing vectors by alternating projections.
IEEE Transactions on Information Theory, 64(5):330113312, 2018.
Irene Waldspurger, Alexandre d' Aspremont, and Stéphane Mallat. Phase recovery, MaxCut and
complex semidefinite programming. Mathematical Programming, 149(1):47-81, 2015.
Gang Wang.
Georgios B Giannakis, and Yonina C Eldar. Solving systems of random quadratic equations via truncated amplitude flow. IEEE Transactions on Information Theory, 64(2)773--99, 2017a.
Gang Wang, Liang Zhang, Georgios Giannakis, Mehmet Akçakaya, and Jie Chen. Sparse phase retrieval via truncated amplitude flow. IEEE Transactions on Signal Processing, 66(2):479-491 2017b.
Yang Wang and Zhiqiang Xu. Phase retrieval for sparse signals. Applied and Computational Harmonic Analysis, 37(3)5531-544 2014.
Ke Wei. Solving systems of phaseless equations via Kaczmarz methods: A proof of concept study. Inverse Problems, 31 31(12):22508, 2015.
Huishuai Zhang, Yi Zhou, Yingbin Liang, and Yuejie Chi. A nonconvex approach for phase retrieval: Reshaped wirtinger flow and incremental algorithms. Journal of Machine Learning Research, 18, 2017.
Liang Zhang, Gang Wang, Georgios B Giannakis, and Jie Chen. Compressive phase retrieval via reweighted amplitude flow. IEEE Transactions on Signal Processing, 66(19):5029-5044 2018.
Shenglong Zhou, Naihua Xiu, and Hou-Duo Qi. Global and quadratic convergence of Newton hard-thresholding pursuit. The Journal of Machine Learning Research, 22(1):599-66 643, 2021.
12

---

Published as a conference paper at ICLR 2024
Additional experimental results are provided in Appendix A, while the proofs for Theorems 3.1 and 3.2 can be found in Appendix B.
A ADDITIONAL EXPERIMENTS
In this section, we conduct a series of experiments to evaluate the scalability of the proposed algorithm, its phase transition characteristics during signal recovery, its resilience when confronted with various levels of input sparsity, and its robustmess in the face of noisy measurements.
A.1 SCALABILITY ACROSS VARYING DIMENSIONS
We first investigate the scalability of our algorithm in terms of varying dimensions. In particular, we extend the range of signal dimensions from 10000 to 50000 and adjust the sample size ratio (m/n) from 0.3 to 0.7. We define successful recovery as the termination condition for the algorithm, specifically when the iterate fulfills: dist(m=" _x² < 10-3, Table offers an in-depth view of the aas
efficiency and scalability of our proposed algorithm.
Table Running time comparison (in seconds) for the proposed algorithm while recovering signals with dimensions ranging from 10000 to 50000 and sample size ratios (m/n) from 0.3 to 0.7. The underlying signal has sparsity level of 100. The reported results represent the average of 200 independent trial runs.
Dimension (10")
1.5
2
2.5
3
3.5
4
4.5
5
0.3
0.212
0.420
0.458
0.675
0.987
1,214
1.472
1.854
2.147
Sample size ratio m/n
0.5
0.242
0.497
0.585
0.912
1.382
1.788
2.237
2.709
3.022
0.7
0.278
0.545
0.833
1.179
1.783
2.224
2.789
3.581
4.252
A.2 PHASE TRANSITIONS FOR BL.OCK-SPARSE, SIGNALS
We present the phase transitions of various algorithms in Figure 4, focusing on block-sparse signals with fixed dimension of n 3000 and sparsity levels set to s = 20 and 30. The signal generation process aligns with the experiments depicted in Figure of the study by Jagatap & Hegde (2019). We generated the phase transition graph by assessing the success rate of each algorithm's 's recovery across 200 independent trial runs. As shown in Figure 4, our algorithm always achieves a higher successful recovery rate than the state-of-the-an methods across various measurement counts.
Probablity of success
0.8
Probabiity of success
0.8
0.6
0.6
0.4
CoP RAM
0.4
COPRAM RAM
0.2
0.2
SparTA
SparTA
Proposed
Proposed
500
1000
1500
2000
2500
3000
500
1000
1500
2000
2500
3000
Measurements [m]
Measurements (m)
(a) 8 20
(b) 8 = 30
Figure 4: Phase transition comparison of various algorithms applied to block-sparse signals with dimension of n = 3000 and sparsity levels s = 20 and 30. The signal generation process aligns with the experiments depicted in Figure 2 of the study by Jagatap Hegde (2019). The parameter settings for ThWF, SPARTA, CoPRAM, and HTP are consistent with those used in the studies by Jagatap & Hegde (2019); 9); Cai et al. (2022c). The results reflect the average of 200 independent trial runs. A recovery considered successful if the relative error was less than 10-3
13

---

Published as a conference paper at ICLR 2024
A.3 PHASE TRANSITIONS ACROSS VARYING SPARSITY LEVELS
Figure 5 presents the success rate of different algorithms as function of varying sparsity levels and the number measurements m. With a fixed signal dimension of n = 3000, we vary the signal sparsity s from 120 and the number of measurements from 150 3000. A signal recovery is considered successful if the relative error dist((b.a=)) success rate: black corresponds and gray to rate between 0% and 100%. As Figure 5 demonstrates, our algorithm outperforms the state-of-ther tart methods at higher values s. For lower values 8, our algorithm achieves slightly better results compared ThWF, similar the performances of SPARTA, CoPRAM, and HTP.
10 The gray level of block represents the 0% successful recovery. white 100% successful recovery.
|
Spmittyo
90
00
Spanlyy 10
20
20
303
1000
1500
2002 2000 2000
500
2000
1500
2002
2300 3000
500
2000
1500 2000 2300 2000
Neasuremenss OT))
Neasurements
Meas.rement (T))
(a) ThWF
(b) SPARTA
(c) CoPRAM
320
120
00
40
1596 2000 Mussuremens [m]
1000 2500
Messurements (m)
(d) HTP
Proposed
Figure 5: Comparing phase transitions among various algorithms for a signal dimension of n = 3000, across different sparsity levels and numbers of measurements. The successful recovery rates are indicated by varying grey levels in the corresponding block. Black signifies 0% successful recovery rate, white indicates 100%, and grey represents values between 0% and 100%. A signal recovery is considered successful relative error is less than 10-3
A.4 PERFORMANCE COMPARISON WITH UNKNOWN SPARSITY
In Table 4, we consider scenarios with unknown sparsity. We input various sparsity levels into each algorithm and compare the success rates of various algorithms in recovering the signal. In these experiments, the underlying signal has sparsity of 30, signal dimension of 3000, and the number of measurements is 2000. We excluded ThWF from the comparison because it does not require input sparsity. As demonstrated in Table 4, our proposed algorithm exhibits significant robustness to changes in input sparsity levels.
Table 4: Comparison of success rates of various algorithms with unknown signal sparsity. The underlying signal has sparsity of 30, a dimension of 3000, and a total of 2000 samples.
Input sparsity
10
20
30
50
70
100
150
200
250
300
CoPRAM
0
0
1
1
1
0.75
0.09
0
0
HTP
0
0
1
1
0.71
0.22
0.02
0.01
SPARTA
0
0
1
0.09
0
0
Proposed
1
0.93
0.85
0.76
0.66
14

---

Published as a conference paper at ICLR 2024
A.5 NOISE ROBUSTNESS
We investigate the impact of the noise level of measurements on the recovery enror of our proposed algorithm. Noise levels are represented by signal-to to-se ratios (SNR), i.e., | where x is the ground truth signal and a is parameter determining the standard deviation of Gaussian noise, as defined in (12). We set the dimension of the ground truth signal to 2000, the sparsity level to ===00 20, and the number of measurements to m 1500.
Figure depicts the relative error of the proposed algorithm as a function of signal-to-noise ratios (SNR) in dB. We observe a nearly linear decrease in the relative error as the SNR increases, implying that its recovery error can be controlled by a multiple of the measurement noise level. This result demonstrates the robustness of our algorithm against noise in measurements. Additionally, the small error bars shown in Figure 6 emphasize the low variability of recovery errors of our algorithm.
12
122
5MR (6h1
Figure 6: Robustness of the proposed algorithm against additive Gaussian noise. The y-axis represents the relative error of the proposed algorithm, while the r-axis corresponds to the signal-to-noise ratios (SNR) of the measurements. The results are averaged over 100 independent trial runs, with error bars indicating the standard deviation of the recovery error. We set = 2000, m 1500, and = 20.
B PROOFS
We provide proofs for Theorem 3.1, the recovery guarantee for foise-free measurements, and Theorem 3.2, the recovery guarantee for noisy measurements. Technical lemmas are presented in Appendix B. which serve as the foundation for proving Theorems 3.1 and 3.2, and their proof's are available in Appendix B.2. Subsequently, the proofs of Theorems 3.1 and 3.2 can be found in Appendix B.3.
For more concise representation, arrange the sampling vectors and observations as follows:
:= [a 02 ... aml y y [3/ y/2 ... yml?. z := [2] 82 ... 2m T
(13)
where VT.. I...., m. We denote a,,s as the vector that represents the i-th row of A. retaining only the entries indexed by S.
B. TECHNICAL LEMMAS
The following lemma serves as an extension of Lemma 7.4 found in (Candes et al., 2015).
Lemma B.1. For any s-saase vector " F" with support S supp(æ*),, let {a be identically and independently distributed as (N In) and define the matrix A as in (13). For any subset S C [n] such that supp(æ*) C S and |S| < r for some integer r S n. With probability at least m C, 2m Cb exp (-c,82m log m), we have
(1 (2nnsss 2a3((a²)²) ||
a aaaiaasa'sa.
61|200112
(14)
provided m 2 C(t)r log(n//)///))) Here C(8) is a constant depending on d, ca. Cb and e. are positive absolute constants, and a, represents the i-th row of A, retaining only the entries indexed by S.
15

---

Published as a conference paper at ICLR 2024
The following lemma is derived through the application of Hölder's inequality and Lemma B.9:
Lemma B.2. Given two vectors æ, 2 E R each with sparsity no larger than and two subsets
S. S [n]. Define the subset R == SUT SLTTUsupp supp(()) Usupp(z). Under the event (20) with support
R, there holds
- -Vo
3 - - m
(33m)144 + ((R1) + 201ogm | + z =||| z||. (15)
The following inequalities in Lemma B.3 can be derived using Lemmas B.1, B.2, B.11, and B.12.
Lemma B.3. Let x² E IR be any signal with sparsity || I s and support S Let {a, I be random Gaussian vectors identically and independently distributed as (N, Define A, y and fi (x) as in (13) and (4) respectively. Given two subsets S,T = [n] satisfyyng [S] S s and [T| S s. Then if m 2 30s for any s-sparse vector x Rº with supp(æ) S T and obeying
dist(a x²) S >|| with 0 < S 0.1, under events (20) and (14), the following inequalities hold:
(i)
4 ||| - Vsssfi(z))) S 12/1u/l. Vu TuERo 1
(16)
where 1 = (2 (2-28- (2 + 7)) 202 and L2 == (6+ 28 + 10y(2 + 7)) | (ii)
Ssssusssiii (z) s = (2 26 + 10 + 7)) |001|
(17)
The next lemma is adapted from Lemma 3 in (Cai et al., 2022c).
Lemma B.4. Let fa= h221 be the sequence generated by the Algorithm 1. Let x² := 2 O
Assume 1| xd S > x=1) that
Then under the event (18) with s, 2s and the event (19), holds
r Set
(**
Aui) s / C1 + 5,) 10'11
4 where Cy= = (1-92) (co + vV ) with co 10-3, and S supp(H,.(** (/
The following lemma provides an upper bound on the estimation error for the vector obtained after one iteration of IHT, as described in (Blumensath & Davies, 2009). To make this paper self-contained, we include the details of the proof for the reader's 's convenience.
Lemma B.5. Given an s-sparse estimate x² satisfying || - -x'll S 7 ||= Define the vector obtained by one iteration of IHT with stepsize 1) to be
t == H. (xx
V/V (x²)).
Under the RIP event (18), it holds that
||u² x:ll S <s $111r 2r'1|
where Ç < 2 22 (22 max 77 (1 )2.)3 + TV Cr (1 + 622)
The following lemma asserts that, given a sufficiently large m. the normalized random Gaussian matrix A obeys the restricted isometry property (RIP) with overwhelming probability. This conclusion well-established in compressive sensing literature (Candes & Tao, 2005; Foucart & Rauhut, 2013).
Lemma B.6. (Foucart & Rauhut, 2013, Theorem 9.27) Let A be defined as in (13) with each vector a; distributed as the random Gaussian vector a ~ (0, Ia) independently for i = 1,2,.... m. There exists some universal positive constants C² C² such that for any positive integer r $ and any 8, E (0, 1), m > 8, rog (n/r), then vm A satisfies r-IIP with constunt o,, namely
1
(1 (1-66) 61]11]1
Ax
(1 (1+6,)]<*]1
Vx E R 1ithth/ S r,
(18)
mm
with probability at least 1 -0mm
16

---

Published as a conference paper at ICLR 2024
The results presented below have been previously established in (Cai et al., 2016).
Lemma B.7. (Cai et al., 2016, Lemma A.6) On event with probability at least m we have
m
a s
x
aa.ssaa (1=~~1) + 22 (mig)))
- 6611=111
m
provided m 2 C())s log Here C(8) is constant depending on 6.
The subsequent lemma, direct outcome from (Soltanolkotabi, 2019), plays a crucial role bounding the term ||
|| Te+1 (z Ax
Lemma B.8. (Soltanolkotabi, 2019. Lemma 25) Let {a, 1= be i.i.d. Gaussian random vectors with mean and variance matrix Let y be any constant in (0, l.. Fixing any Eo > 0, then for any 8-sparse vector æ satisfying dist(æx,"* x²4) S >oo|=>"|, with probability at least e C. there holds that
m a 2²12 1 ((aFa)(afat)50) S (1 2o)
21
Σ
6D + Ao
| - 201
(19)
m
20
provided m log(ns).. Here C²2 and C/2 are some universal positive constants.
The next lemma, derived from (Cai et al., 2016), asserts that the induced norm of the submatrix of random Gaussian matrix A bounded by its size.
Lemma B.9. (Cai et al., 2016, Lemma A.S) Let {a,);- be identically and independently distributed as (0, In) and define the matrix A as in (13). Given support S S [m] with cardinality with probability at least 4m there holds
Ass12-s4 S (3m) 8 2 loggogo m
(20)
where As represents the sub-matrix of A, retaining only the columns indexed by S.
J.A. A can be deduced by The subsequent result concerning the spectral norm of submatrices of vm
employing the restricted isometry property of matrix J. A.
Lemma B.10. (Needell Tropp, 2009, Proposition 3. Suppose the matrix satisfies the inequality
(18) for values of specified and s'. Then for any disjoint subsets S and T of 2,, m) satisfying (S] < and [T] I s', the following inequalities hold:
1
A
I v +
(21)
771
A AT
S 6xta'
(22)
1 - 8, 5
= A As
3 11+60 11+
(23)
The following lemma provides the expectation of the sub-Hessian of the intensity-based loss fi at at. As this result can be easily derived through basic calculations, we will not delve into the details here.
Lemma B.11. For any subset S C [m]] such that supp(x*) S, the expectation of V2 (
E (Vsss(((a') (1==*)2((nssss 20r (ab)
and it has one eigenvalue of 6| and all other eigenvalues are 2 2|
The next lemma is the so-called Weyl Theorem, which a classical linear algebra result.
17

---

Published as a conference paper at ICLR 2024
Lemma B.12. Suppose M, N E Rux* are symmetric matrices. The eigenvalues of M are denoted as X1 2 X2 2 ... An and the eigenvalues of N are denoted as μ1 2 P2 2 ... Pa. Then we have
k4 A | - N Vi 1,2,
Lemma B.13. (Hoeffding--yy) prppedinng--eeeeeeee inequality) Let X1, .N be independent centered sub- Gaussian random variables, and let K- maX4e[N] || wi. where the nub-Gaussian norm
|x == supp 11/2 (E (1
P21
Then for every b [bi;... ;b,] E INN and every 2 O, we have
T
Σ Xx
20 S Seexp
cf2
Here e universal constant.
Lemma B.14. (Bernstein-type inequality) Let X,.... N be independent centered exb-exnnnetiaal random variables, and let K max, [| where the nub-exxonential norm
|| supp (III
P21
Then for every b = [b13--. ;bN] E E and every 2 0, we have
2 t
(-emin
P
by Xx
exp
K2 bl12 Kbbll
Here is universal constant.
Lemma B.15. (Bernstein's inequality for bounded distributions) Let XI, NN be independent mean zero random variables, such that [x|] < K for all Then, for every 2 O, have
N
(- 121
X
zt
S 2exp
7-2 K1/3
Here a² = ENN EX? the variance of the sum.
B.2 PROOFS OF TECHNICAL LEMMAS
In this subsection, we provide proof's for the technical lemmas: Lemmas B.1, B.2, B.3, B.4, and B.5.
B.2.1 PROOF OF LEMMA B.1
Proof of Lemma B.I. If [S] = | then S = S and the result is established by Lemma B.7. We
now focus on the case where := S o We first define two matrices as follows:
a; a Ss T
ai
(In)) 2x (x
H
a
ag
a 1T
x (I
Then we rephrase the term on the left-hand side of (14) follows:
T
H
<
1
Σ
a
an
12 Paassaaaa (1)=2)1(((n)))) + 2x (2r (2ng.))
1=1
a
x
a
a iT
+
a
a ot
-1
m
i=1
18

---

Published as a conference paper at ICLR 2024
The first term can be bounded with overwhelming probability via direct application of Lemma B.7 as below whenever m 2 c log(ns)))
a i x sa - (11z5|)(In)s=+ + 2x (2)))) = 6/4.
m
For the other two terms, it is enough to consider æ² because the unitary invariance of the Gaussian measure and rescaling. Then for prefixed subset S (thus T is also fixed), there exists unit vector u E In with sup S such that the second term is equal to
Σ m
1
m
(a.(1))
m m Σ (a)(1)) (af u).
We emphasize that af 2 a,,- uy is a random Gaussian variable distributed as (0,,) for all E IR", || = 1, supp(u) C T and independent of (1; for all i E [m]. So for one realization of (a:(1)), Lemma B.13 (Hoeffding-typp inequality) implies
-
Σ (a)(1))
af > S eeep
com2t 2t2
T2 (a:(1)16)
for any / >> 0. Define the set W to be a collection of all the index set in [n]] with cardinality no larger than r, i.e. WT = (S S [n] : |S] sr. Note that the cardinality of WP can be bounded by Σド= ()) < em) Taking t 5/8, together with union bound the set Wh obtain for any
subset = St,St.
1 Σ (a)) P
of
> 8/4] S e eexp (- (-64
64 E [1:(1] + 6r log(n/r))
Com 82
Now with application of Chebyshev's inequality, we have Di=1 |a;(1)] < 20m with probability at least - Cym Substituting this into the above, we conclude that for any subset S S Wp. whenever m 2-r log(n/r) for some sufficiently large c4.
m
(1)( of.
< 6/4,
with probability at least Cym e expl - m))
Assuming a² e₁ and fixing subset S C [n] (thus T is also fixed), we can obtain:
(1)11) (a.,7a ITT) 1 (1)
a.
<
for which there exists unit vector u E R with supp(u) S T such that
(a)(1))² (a.,af, I[T]]
Σ (a.(1))) uat -1/-
Note that (a)(1)} is independent of {a an application of Bernstein's 's inequality implies
m
(1))) (aT - 2 2 exp C6 min (d1, dz)):
P
m
where di
o Em ( and d2 maxieim] 2 Taking t 8/8, together with a union bound on the subset WP we obtain for any subset S E wm.
P
m Σ
(1)P (a.Tar TT)
> 5/4
< 2 2xx (-c6 min (dg, da) + Gr (og(n/r))
19

---

Published as a conference paper at ICLR 2024
where d3
6 - (1)l and d4 the union bound, we obtain:
m 6
Bor maxxe(m) (1) Applying Chebyshev's inequality and
md
m
Σ
la:(1)1* S 10m and max [1, < 10 m KEim
hold with probability at least 1 cgm -4 To conclude, for any subset S W., when m 2 cr00rrrroppp (ogg(n/r) for some sufficiently large constant Co.
(1)12 (ai,TaT 1rri) 8/4
2 2exp (-c1002m/lol m)). Finally, an application of with probability at least 1 - cam Chebyshev's 's inequality implies
(a;(1)) S 5/4
with probability at least 0110-2mmm -1, The proof is finished by combining the above bounds and probabilities.
B.2.2 PROOF OF LEMMA B.2
Proof of Lemma B.2. A simple calculation yields
3 A D( Amr |Az 2 A) V s.) V(sTf((()=
m
There exist unit vectors u, v E R" with support supp(u) S S and supp(v) CT such that
| V.T//(x)) - -V3 (r))\(
3
| Ax |==1))AT]
m
3
(laTx - -1a x)) (a)
m
3
(la - 10 2र ()) (a un ( (
m
T (x + 2R)) (att (x: =R)) (a)
m
1=1
3
| 1 + =||||= = 2||
m
1A
3
(3m)+ + 1/2 + 2 log m
|| + =||||= 2|1-
where the last inequality is implied by Lemma B.9.
B.2.3 PROOF OF LEMMA B.3
Proof of Lemma B.3. We begin by proving result (i) in Lemma B.3. Define R TUSUS Then [R| < 3s. Applying Lemma B.2 yields that
3 m
(3m) 1/4 + (3s) + 2 V1og || + æ x
S 10 x + x
||| - 11
20

---

Published as a conference paper at ICLR 2024
provided m 30s Furthermore, following from Lemmas B.1, B.12, and B.11, and the interlacing inequality, we obtain that:
Amin (.))(()))) min (V (
Amin (R. (x*)) R (x) - - (at)
min
(E [ RI (z*))) V R (x
RR. (29)] 10 |0|| + x || x²)
2 | | - 26 | - 10 x x
|
x²1|
(2 28 (0
||
where the last inequality is implied by dist(æ x²) S >||=7||- Similarly, the upper bound for the largest eigenvalue of (x) can be bounded as follows:
Amax
V (æ)) I Amax R (a))
max RR (z)) + V x((((()) Ri (x*)
max (E) RR (z*) + R (x*) - E RR (z*)]] + 10 10|l= + x²1 ||| x²)
< | 2 + 26 ||=
2 + 10 | | -
- (6 + 28 + 10y((2 + 7)) ||
Now we turn to proving result (ii) in Lemma B.3. Define R = TUSUS St. Thus | S 3s. For the disjoint subsets and S*,S, we consider sf (x), which a submatrix V (x)
1||==11/2. We note that the spectral norm of submatrix never exceeds the norm of the entire matrix. By employing the result from part (i), can deduce that
V <f,((())) - |||=>1|0]||
< max (6 + 25 + 10 + - 4,4 (2- 25 10y(2 + +))} ||
28 + 10 7)) ||
completing the proof.
B.2.4 PROOF OF LEMMA B.4
Proof of Lemma B.4. Define the sets (G follows:
GK = ii sgn x²) sgn(a xt), S m).
With ** == sgn(Ax") Ax*), where z is defined (13), we deduce that
1 A
1
sgn (a xt) - sgn (af' a x²
a x
1 (aT (aT
m
vEgg
(24)
4
21 TV
| x
6o +
7)2
20
C
where the first inequality follows from [sgn(a? x*) sgn(a, x sgn(a,"a
2 and sgn(a, x*)
[Taa 2²) = 0 on Gk. and the second inequality follows from Lemma B.8. Together with (21) in Lemma B.10, (24) leads to
AF Afa (**
Auti| s CC,(1 (c) x x²11
completing the proof.
21

---

Published as a conference paper at ICLR 2024
B.2.5 PROOF OF LEMMA B.5
Proof of Lemma B.5. Define S := supp (at)), Th+1 == Sk+1 U.S,, and ob = x - -1/V SA (2x
Since u is the best s-term approximation of o², we have
|| 20 S || - --01l
which implies
/ute ot Te111
|
S 1
because the relation supp (x*) - Tk+1 and supp (x+) S Te+11 Then, follows from the triangle inequality and the inequality above that
|| 1 z011 = [ 1/ x Tats
Thes + Vete 20 -
U Tx+1
U
Tut1
Tu+1
+
Tets
Tett
(25)
ut Thi1
2
Ti/e+1
Define z* = 3 sgn (Ax*)) with in (13). Using the definition of wht direct calculation yields
2 ||
-x T1+1
"
A Th-1 A (x - x²) + ATim (z" - -==a)
T2+1
<
AT (2# - -Axe))
+ (
AT Tets ATez)) ) (atti 2 atte))
11
+
AT
Te4s ATATeH
(oth x0 TeiT
We will now proceed to estimate T,, I2, and I3 sequentially. For l An application of Lemma B.4 yields that
A (z* - --+
s | C, (1 820) | x²
(26)
For I2: Let 1) E (0, ) follows from (23) in Lemma B.10 as well as Weyl's inequality that
(1 - (1 + )2.)) |ell s |
"7 A A S (1 17(1 82-)) ||/a|],
m
for any Te+1l, which deducts
(1
1 AT... Tet A m
(atie - x = )
S (1 - 17(1 82x)) / - a ata
For Ig: (22) in Lemma B.10 implies that
7 AT Tat1
AT Tu+1
[x - x² TAYT = S 17830
|x0 - -xa =
Combining all terms together, we obtain
k++ Tu-1
2 2
< J + T2 + I3 S 2(1 + 12) T3
< 2 max {1d83s, 1 - 17 (1 )2.)) 1| - x'1| + C, (1 + 820) 1|0 =
22 max (nōss, 1 1/ 82.)) + 'n ///C
C, (1 + )2.) 1100 2x
We complete the proof by using (25).
22

---

Published as a conference paper at ICLR 2024
B.3 PROOF OF THEOREMS
We present the proofs of Theorems 3.1 and 3.2 in this subsection.
B.3.1 PROOF OF THEOREM 3.1
Proof of Theorem 3.1. In this proof, consider the case where || 2ل Consequently, the distance between the initial estimate and the true signal is given by dist(x") "est((x"2"" 2²)) || =²I. Note that the case where ||x0 + æ$Il ||x8 - x= can be addressed in a similar manner. For the purpose of this proof, assume, without loss of generality, that the true signal has unit norm, i... ||> =1.
||x00 + xa||
Let x represent the k-th iterate generated by Algorithm 1. Given an 8-sparse estimate æ with support Sk, which is close the target signal, i.e. dist a²) & 7 /||=>]]l For any 0 S ( S 1, denote a=(() == 2² 1(x* 2.). It is evident that supp(ar()) S supp(æ*) Uuupp(#*)) and the size of the support of at) is at most 2s, Le., (supp(x(()) S 2s. Furthermore, we obtain
| xt))|| - (1 t) 1 x
s | -=*-1"
and
x + stt]| [1 t) (1 t)z² (1++) ||~*|| + (1 ||~i||
$ (1 + (1 + 7) 1 + (1 - t) ||~i||| = 2(7 + 1) x
where the last inequality holds because 1 1.
Assume events (20) and (14) hold, then by (15) in Lemma B.2 with x², = an(),
Sk+1, T Sy U S there holds
S (x*) V S S.us:
3
(3m) + (3s) 2 logm)
22
1 x(t)| ||a x(t)|
(27)
m
<10
x + a+|||| 1 x
<In 1| - 2il
where La :- 20(r + 1)||==*||. Also, Eq. (16) in Lemma B.3 with = x SS+1 indicates that
( (V fo (x*)) | (Amin ( (2*1 S
1
(28)
T
Moreover, by the mean value theorem, one has
fi (a*) f
fr (2+)
42 (z(1)) (s+- dt.
(29)
We also have the following chain of equalities
| 1
2i1
1-t Sa12
112
|
112)
112] 1/2
Se11
SE
sta.
|
1|2
p 112 1/2
(30)
where p is the search direction that can be obtained by (10).
We first bound the term I in (30). Note that z's.., is subvector of œt nt, Based on this observation and applying Lemma B.5, we deduce that
I = | ||ar* - u²
S < <||a0* - -x011
(31)
where (= (22 max max
- n1(1 - -d24]} + 11 C1 62s) with < I48
23

---

Published as a conference paper at ICLR 2024
We now proceed to bound the term Ig in (30). By plugging the expression of P we obtain that
T2- SeLr St(x"))
(sens + S x Se+1
(32)
We further can deduce that:
1
v Se
S (r (ass) - Vsendr(("))) Vit. Se+i fr(#*)) (x - 2bse
1
1
V fi - V (i (x* (a*==")) dt
+ L x
4
11
1
V (( V Suorfi | dt
S/a
||00 x11
11
l
$ 1 1|, x2/1? ( (1-
Lh
(i-t) + {la || a11
h
= Pi || x
where the first inequality follows from (28), the second inequality is based on Lemma B.3 and (29)
together with the fact that V (r((=)) = O, and the third inequality is derived from (31), and the last inequality is obtained from (27). The equality includes Pr, defined as follows:
Thp 21'1
C23
20(1 + 7)n
C(2 + 25 10> 7))
x 24
P1
4
2(2 - 25 - (2
2 26 - 10 + 7)
2C(1 + 10 + + ((2 ?))
- (2 + 7)
By substituting the upper bounds of terms 1n and I2 into (30), we obtain:
- s V pi + <²
|
æ''|
(33)
Let / pi ç², then م 1 can be ensured by properly choosing parameters. For example, when da. 0.05 and 8 = 0.001, and set 1) = 0.95, then S 0.6 < 1 provided that y S 0.01. Therefore, || - x²1 S +|| a²1 S m>| for some p E (0, 1).
Let R2 S U S 1 U S Assume that the event (20) with S = Rr holds for k iterations. As stated in Theorem IV.1 of (Jagatap & Hegde, 2019), the initial guess x" is guaranteed to satisfy dist(x", x² < > || By applying mathematical induction, we can show that for any integer O k k1, there exists a constant p E (0, 1) such that dist( w4) S p dist(x" 2*).
Let K be the minimum integer such that it holds
|001 pK < min.
(34)
We then assert that S C S₁. for all 2 K. This because if it were not the case, there would exist an index i S 0, such that ||ac* 20 2 |a:|| - However, this contradicts our assumption ||z I=1 S >1|=011 S <>110")1pt < xmin- Consequently, based on (34), we derive:
K 1000011===1////
+1 Ca log
+ C..
log(p-)
Note that S₁- = S for all ts 2 K implies R2 R++1 for all ls > K. As result, the probability of event (20) occurring for all 1: 2 can be bounded by 1 4Km To conclude, with probability at least 1 (4K C. C.)m-L Cge- C.m/ log there holds
h+1 x²) S p dist(x* )²).
with some (0, 1) provided m 2 Cys² log(n/s). Moreover, for ls 2 K, utilizing the result the obtain:
result Sk+1 St and consequently x s
2/2/2 =
|| | +1 - 2i11 =
x
1-1 Set1
x
S.
||
+
k+1 S441
x
a -I - 23301 + pPS...
St+1
24

---

Published as a conference paper at ICLR 2024
We can further obtain
|
2
- 11
V
S Supr ( (x*)x* fr(")f(("") V Sap (x)))()s S
V S (x (x) - xd) V Vi -
1
/ (i - V ( (xt - - 2r)jde
$ 11
1 /
1
V
Sets Scust V uss 1 - - dt
La 1/=t
-02113 / (1- t)dt
La 1|at*
2x:1?
241
where the first inequality follows from (28) and (32), and the last inequality follows from (27). Consequently, the sequence {**} converges quadratically.
L
B.3.2 PROOF OF THEOREM 3.2
Proof of Theorem 3.2. Owing to the independence between the sensing matrix and noise, along with
the assumption E(c) 0, Lemma B.1 remains valid in the noisy setting. Therefore, by applying Lemma B.I, m = (2(s log // s). then with overwhelming probability, it holds that
V (x*)
E
[V (=*)] s (5 + c) ||
Under the assumptions in Lemma B.3, it follows that
[u]l < V (x))ul s E111ell Vu RR
(35)
and
| f7..fi( ( (z) S 1/z
(36)
where = (2 25 2c 10 + 7)) x |l. 1/2 = (6 + 25 + 2e 10>(2 + 7)) ||, and
1/2 = (2 26 + 2c + 10y(2 + 7)) ||
In the noisy case, ( 122 is given by
z² = (1 + c) @ sgg((z*)
Then using the same argument as the proof of the inequality in Lemma B.4, have
7 AT Awi)
(=*
"
AT Tk-1
(|Ax*) @ sgn( (Ax Axt))
m Te+.
m
C'n
+ m
AT Te11 (co ((o(((Aa))))
(37)
C'n / VI+522lle 82. |c]l,
< <ny
C1
82a 1|2
x'1l= +
Vm
where the last inequality follows from Lemma B.4 and (21) in Lemma B.10. Then, we modify Lemma B.5 to estimate ||wu - æ0I| in the noisy case. All arguments Lemma B.5 go except that the estimation of in (26) should be replaced by (37). Thus obtain
C'n Vm
lu k
a:
S
1|00 - x0ll +
/ 12.lelll
(38)
where tt 7) the as those in Lemma B.5.
25

---

Published as a conference paper at ICLR 2024
We now proceed prove the convergence property for the noisy case, employing similar argument as in the proof of Theorem 3.1. Notably, equality (30) remains valid in the presence of noise. As for the first term in (30), since " is a subvector of x² u², it follows from (38) that
|| || =21l c110 x +
C' / 82el|el1,
Vm
where Ç 2 /2 max((nss n11 82s)} } + ny C1 d22))
with rj < 1/(1 + 626).
Furthermore, the second term in (30) can now be estimated as
| - I's.. + +
V fi Si Ja S St fi (a )
1
ViLtifilz' f(a)) + m
A (EO (Azi))
V fi - x0 Fi - =z7)d()
1
V 1. fi a -1 imm
1
A5. (c @ (Ax*)) ||
+ 1/
vI + 6g Axx 1l
V2 F V ((((t) (x* - =xa)dtt
St
+
|||||
7 I/ Seenssuusrf V
s(ssss | ctt
5122 4
||+ 2+||+ C'GgnV++ + 82 + vI 0= A 11 |100|1)) |||||
+
l'mm
La P 110 - -0113 J/' (1 t)dt + 51/2 ||0t
O'GynV1 + 82 + /1 + 8, Ax
-x--l x'1|
||||
livm
=p2l|=* ++||+ C'GnnV + 82 + v1 8g A
|||||
(39)
l//m
where the first inequality follows from (35), together with the equality ffr(x +A" (e (Ax*)) =
O, and the third inequality is based on (36) and the inequality a @ bll < ||a||>>||0] for any two vectors a, b. The last equality includes P2. defined as follows:
Ln
5h m <
20(1 + 7)m
((2 + 10y(2 7))
1 - x²1|+
P2 =
2m's
2(2 2 2E 10 +7))
2 2 28 2e 10>(2
2c11 + (1+ <(2 + +))
2 - 2c - 10 +
Following from (a² + b2y2122 S a for ab 2 0 and putting the two terms together yields
| 2 S p' | 2 + e||e|].
where p¹ = =p2 + < and v= G'U +//m//111110 tvI+A A Noticing that
l/ym
1
log(ms)
Ax²
max aa 2 < [E(m]
max |روو| - ||0+||| = 3 m ie JESt
|
Vm
v
Jest
m
with probability least 1 - (ms) 2, Consequently, J ||Ax loo can be quite small. As result, properly setting parameters can lead to u E (0, 1) and p' E (0, 1).
26


Reviews:
[{"review_id": "MJ2ipDCVKO", "reviewer": "Reviewer_s5s7", "summary": "In this paper, the authors address the sparse phase retrieval problem, which focuses on the recovery of a sparse signal from a limited number of phaseless measurements. The authors propose an efficient second-order algorithm based on Newton's projection, which maintains the same per-iteration computational complexity as popular first-order methods. The paper's primary focus is on the sparse phase retrieval problem, where the goal is to recover a sparse signal from phaseless measurements. The proposed algorithm is theoretically guaranteed to converge to the ground truth at a quadratic convergence rate after some upper-bounded iterations. Numerical experiments demonstrate that their algorithm outperforms state-of-the-art methods in terms of achieving a significantly faster convergence rate.", "questions": "See weakness.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The authors have introduced an efficient second-order algorithm rooted in Newton projection that maintains an equivalent per-iteration computational complexity to widely used first-order methods. The algorithm is accompanied by a theoretical guarantee, ensuring a quadratic convergence rate when converging to the ground truth. Their algorithm's superior performance is substantiated through numerical experiments, demonstrating a notably faster convergence rate compared to existing state-of-the-art methods.", "weaknesses": "The experimental section could provide a more comprehensive perspective on the algorithm's efficiency by considering different dimensions. For example, it would be beneficial to explore why the experiments were consistently conducted with fixed dimensions, such as signal dimension (n) at 5000 and sample size (m) at 3000. This choice of fixed dimensions raises questions about the scalability of efficiency when the dimension varies.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "9ETQQvY8vM", "reviewer": "Reviewer_M25C", "summary": "This paper proposes a second-order Newton-type algorithm for solving the sparse phase retrieval problems. It is shown that the algorithm enjoys a quadratic convergence rate with similar sample complexies relative to existing algorithms. Numerical tests verify the theoretical findings.", "questions": "In different papers of sparse PR, SPARTA performs better and faster than e.g., ThWF; see e.g., the paper CoPRAM. It is strange why SPARTA is the worst performing method in almost all tests provided here. Please explain and replicate an experimental setting in SPARAT or CoPRAM paper using all methods.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "+ The design of the algorithm is interesting and delicate which is different than the existing sparse PR algorithms based on hard/soft thresholding. \n+ The paper is easy to follow.\n+ Extensive convincing numerical results.", "weaknesses": "- The first stage of finding the initialization is not reinvented. \n- Known sparsity level s a priori.\n- When noise is considered, the convergence rate recovers that of the first-order algorithms. In practice, all measurements are noisy. Please explain what is the advantage of the proposed algorithm for practical settings then.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wzVRjBajXY", "reviewer": "Reviewer_qMEw", "summary": "This paper proposes a new (second order) algorithm for the sparse phase retrieval problem, which is based on Newton projection. The algorithm is proven to converge to the ground truth at a quadratic convergence rate. This paper also provides various experiments to demonstrate its performance against SOTA algorithm (e.g., with/without noise and phase transition plots).", "questions": "N/A", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper is easy to follow for audience with phase retrieval background.\n2. Theoretical results for both with/without noise are provided.\n3. Experiments are sufficient, and demonstrate its performance against SOTA algorithms. I really like the phase transition plot, which will be asked to do so if it is not provided.", "weaknesses": "Although it is easy to follow for audience with phase retrieval background, I suggest authors to provide more details such as intuitions and figures of geometry to improve the readability for a broader audience. A good example will be ``Yuxin Chen and Emmanuel J Candès. Solving random quadratic systems of equations is nearly as easy as solving linear systems. Communications on pure and applied mathematics, 70(5):822–883, 2017.'' I will raise my rating if it could be realized.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "Qnqg8nvOKL", "reviewer": "Reviewer_5Ggq", "summary": "This paper investigates the sparse phase retrieval problem, which seeks to reconstruct an unknown $s$-sparse vector $x$ from $m$ measurements $y_i=\\langle a_i, x\\rangle^2$. This problem has been extensively studied, and several first- and second-order methods have been developed. The main contribution of this paper is to propose a second-order algorithm for sparse phase retrieval. The cost per iteration is $\\tilde{O}(s^4+s^2n^2)$, and the number of iterations is $\\log \\log (1/\\epsilon)+\\log (\\|x\\|/x_{\\min})$. The algorithm also considers noisy measurements and proves a linear convergence rate. Additionally, the paper conducts several experiments in different parameter regimes and compares the new algorithm's performance with previous algorithms. The experimental results show that the new algorithm has significant advantages in efficiency, accuracy, and robustness.", "questions": "1. Page 2, bullet 2: what is |x|?\n2. Page 2, bullet 2: “It is important to emphasize…” Note that the non-asymptotic convergence of Newton's methods has been studied in many optimization problems in theoretical computer science and optimization. It is indeed interesting to obtain an iteration bound for the sparse phase retrieval problem. But this emphasis seems to be unnecessary.\n3. Page 3, Sec. 2.1: the formulation of the problem is unclear. For example, are the sensing vectors $a_i$ fixed or chosen by the algorithm?\n4. Table 1: it would be better to add references to these methods.\n5. Is there any lower bound on the computation complexity of this problem? \n6. Page 4, Sec. 3.1: $O(n^3)$ can be replaced with $n^{\\omega}$ using fast matrix multiplication since you are referring to the complexity. \n7. Page 4: “We identify the set of free variables using one-step IHT”, IHT is not defined.\n8. Eq. (4): is there a closed-form expression for the s-sparse hard thresholding? Does it just pick the s largest magnitude entries of the input vector?\n9. Bottom of Page 4: is there any intuition as to why using f_A for free variable identification can improve the algorithm? \n10. Thm. 3.1: “ There exists positive constant” -> “ exist”\n11. Thm. 3.1: what is the convergence rate of the first K iterations?\n12. Thm. 3.1: what’s the scale of $\\rho$? Is it smaller than 1 or bigger than 1?\n13. Thm 3.2: in this theorem, the convergence rate is characterized for the first K iteration. Why is this different from the noiseless case? Does it mean the algorithm is not robust in a long time?\n14. Lem. B.2: Does $A_S$ mean the principal sub-matrix indexed by $S$? Or does it only include the rows in $S$? How tight is this lemma?\n15. Lem. B.4: why $r=s$ and $s’$? Is it a typo here? Also, what is $\\delta_{s+s’}$?\n16. Lem. B.7: $\\mathcal{T}$ is not defined here.\n17. Eq. (29): it would be better to recall what is $p$ here.\n18. Proof of Thm 3.1: It actually shows a linear convergence rate for the first K iterations. After that, the support of the ground-truth vector is fixed, and the quadratic convergence appears. It seems that in the fast-convergence phase, the algorithm uses Newton’s method to solve a general phase retrieval problem (without sparsity) of dimension $s$. Is the quadratic convergence rate known before for this problem?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The problem studied in this paper is of great significance in both theory and practice. The new fast algorithm proposed here can be utilized for a variety of real-world applications. Theoretically, the time complexity for each iteration of the algorithm matches the previous state-of-the-art first-order method, while the number of iterations improves on all previous first- and second-order methods. This is an essential theoretical improvement. Empirically, the proposed algorithm can be efficiently implemented and runs very fast and robustly. Furthermore, this paper is well-written, and most statements appear to be mathematically sound to me.", "weaknesses": "The super-linear convergence rate of the iterative method requires a warm start, meaning that it only holds after $K$ iterations. The main text of this paper does not provide enough intuition for why this initialization is necessary. Second, the initial point of the iterative method is generated via the sparse spectral initialization method by Jagatap and Hegde (2019). This makes the whole algorithm quite complicated. Third, it is unclear what are the technical differences between this paper and the papers by Cai et al. (2016) and Cai et al. (2022). It seems that the proofs in this paper rely on several technical lemmas from previous works.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)