# Amplitude-Based Smoothed Sparse Phase Retrieval Algorithm

Authors: Zhuolei Xiao; Zhengyan Liu
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

Sparse phase retrieval mainly solves nonconvex and nonsmooth problems. Aiming at the nonsmooth problem in sparse phase retrieval, we propose a smoothing algorithm which is called sparse smoothed amplitude flow (SPSAF). The proposed SPSAF algorithm is an amplitude-based nonconvex sparse smoothing phase retrieval algorithm. First, the original phase retrieval loss function is smoothed without modifying the gradient in the gradient refinement stage, thereby reducing the computational complexity of the overall algorithm. Secondly, the support of the original signal is estimated by differential analysis of the gaps, and the initialization can be obtained through a carefully designed method based on this support. Finally, we get sparse estimates by gradient descent based on hard thresholding. Numerical experiments show that the proposed SPSAF algorithm has significant improvements in recovery performance, convergence speed, and sampling complexity. Further, the SPSAF algorithm is stable in noisy environments.


---

## Body

## 1. Introduction

In real-world applications, especially in the field of channel estimation and image processing, the primary signal is naturally sparse or accepts sparse representation after some known deterministic linear transformations [1–5]. When the measurements are undersampled, phase retrieval will introduce sparse assumption or direct sparse transformation for the signal [6–11]. For instance, the sparsely distributed stars observed in astronomy and the sparsely distributed atoms or distributions observed in crystallography can all be regarded as applications of sparse phase retrieval. This recovery method for solving sparse signals is called sparse phase retrieval or compressed phase retrieval, which has fundamental significance in signal recovery and reconstruction [12, 13].Sparse prior is essential for solving sparse phase retrieval problems, and many methods for sparse phase retrieval have been developed. Based on the PhaseLift algorithm, the CPRL algorithm [14] is designed by utilizing the ℓ1-norm regularization term, which improves the one-dimensional (1D) signal to N×N matrix, and finally obtains the sparse solution. Aiming at the application problem of the PhaseMax algorithm under sparse conditions, the SparsePhaseMax algorithm was proposed in [15, 16] by utilizing ℓ1-norm regularization. Based on the alternating minimization algorithm, [17] proposed the Sparse AltMinPhase algorithm, which solves the sparse solution of the equation by alternating between the missing phase information and the candidate solution. [18] developed a probabilistic phase retrieval algorithm based on generalized approximate message passing. [19] proposed a greedy sparse phase retrieval algorithm (GESPAR) by combining the damped Gauss-Newton method with the two-element optimization method. It dynamically updated the support of the signal by using the two-element optimization method and solved the optimal value of the current support by using the damped Gauss-Newton method. [20] proposed a soft threshold TWF algorithm for sparse phase retrieval. Qiu and Palomar developed C-PRIME and SC-PRIME algorithms based on the majorization-minimization method [21]. They replaced the objective function approximately with the convex function. [22] proposed the sparse truncated amplitude flow algorithm (SPARTA), a sparse phase retrieval algorithm based on the TAF algorithm. The SPARTA algorithm improves the performance of sparse phase retrieval by introducing truncation procedures. Zhang et al. proposed compressive reweighted amplitude flow (CRAF) based on the RAF algorithm in [23], making the gradient descent direction of sparse phase retrieval more effective and accurate by introducing a weighted program.When phase retrieval adopts a noise-free Gaussian random measurement model, the CPRL algorithm can accurately recover anyk-sparse signal from Ok2logn measurements, and the computational complexity is On3 [24]. Sparse AltMinPhase and soft threshold TWF algorithm need at least Ok2logn measurements, and the computational complexity is Ok2nlogn. The SPARTA algorithm accurately recovers the signal from about Ok2logn random Gaussian measurements with a computational complexity of about Ok2nlogn, and the total running time is a positive correlation to the time required to read the data. The CRAF algorithm can accurately recover the signal from Ok2logmn random Gaussian measurements.The performance of the amplitude-based phase retrieval method is better than that of the intensity-based phase retrieval algorithm in both numerical and experimental verification, which is the same in a sparse environment [25, 26]. The soft threshold TWF algorithm in [20] solves the smoothing problem in the sparse phase retrieval model based on the intensity yi=aiTz2i=1m. It adopts the adaptive hard threshold iterative algorithm based on compressed sensing. In the gradient optimization stage, each iteration only retains some maximum indexes.On the other hand, the Sparse AltMinPhase algorithm in [17] estimates the support of the original signal and solves the k-sparse phase retrieval problem. The Sparse AltMinPhase algorithm solves the sparse phase retrieval problem by alternating minimization and resampling conditions and performs matrix inversion in each iteration. Numerical experiments show that resampling is a necessary condition for the Sparse AltMinPhase algorithm; that is, many measurements are needed to estimate the support accurately, and the importance of the recovery of the support for signal recovery and even the whole phase retrieval is self-evident.Among the above methods, the illustrious amplitude-based sparse phase retrieval algorithms include the SPARTA algorithm [22] and CRAF algorithm [23]. The SPARTA algorithm applies the TAF algorithm in a sparse environment, divided into two stages: initialization and gradient refinement. Firstly, a reasonable rule is used to solve the support, and some simple power iterations are used to solve the initialization problem to obtain the sparse initialization. Through a series of truncated gradient iterations and the hard threshold of each iteration, all indexes except k maximum values are set to zero to realize the continuous update of the initial value, where k refers to the signal sparsity. The CRAF algorithm further considers the structured sparsity pattern based on the SPARTA and RAF algorithms. It proposes the amplitude-based (block) sparse phase retrieval problem. The CRAF algorithm developed a new sparse spectrum initialization method in the initialization phase, wisely assigning negative or positive weights to each sample. By this method, the mean value of the initialization matrix obtained by the CRAF algorithm has improved performance. Then, the CRAF algorithm uses the reweighted gradient in the gradient refinement stage to gradually improve the initialization of the hard threshold iteration.Based on the SAFPR algorithm and CRAF algorithm, we propose a new sparse phase retrieval algorithm, called sparse smoothed amplitude flow (SPSAF). The SPSAF algorithm is an amplitude-based nonconvex sparse smoothing phase retrieval algorithm with two stages: initialization and gradient refinement. In the initialization stage, we estimate the support by a reasonable rule. Then, the initial estimation is obtained by a carefully designed initialization method based on the support. The gradient descent method updates the initialization estimate based on the hard threshold in the gradient refinement stage. Numerical experiments show that the SPSAF algorithm is robust to additional noise in the finite support. Compared with the existing typical algorithms, the recovery performance and speed of the SPSAF algorithm are significantly improved.The rest of this paper is organized as follows. In Section2, the sparse smoothing phase retrieval problem is formulated. Section 3 introduces the SPSAF algorithm in initialization and gradient refinement two stages, respectively. Numerical experiments are provided in Section 4. Finally, Section 5 summarizes this paper.


## 2. Sparse Phase Retrieval Problem

### 2.1. Amplitude-Based Sparse Phase Retrieval

The amplitude-based sparse phase retrieval problem is a kind of problem to reconstruct the sparse signal from the phaseless measurements [19, 20, 22]. Mathematically, we can describe the sparse phase retrieval problem as a set of phaseless quadratic equations, namely,
(1)ψi=aiTx,1≤i≤msubjecttox0≤k,where k represents the level of sparsity. ·0 express a zero-norm operator, that is, the number of nonzero elements. The signal x is k-sparse, consisting of up to k nonzero elements and n−k zero elements. The goal of amplitude-based sparse phase retrieval is to reconstruct the sparse signal x∈ℝn or ℂn based on the given measurements ψi.For theoretical analysis, suppose sparsityk is prior known. In the real case, the number of measurements required to the k-sparse signal x recovery is at least m≥min2k,2n−1. [27] pointed out that k-sparse signals can be reconstructed by selecting m≥min2k,2n−1 measurements on general positions in real number space ℝn. Similarly, in the complex cases, the number of measurements required for k-sparse phase retrieval is at least m≥4k−2 [28]. Due to the lack of more phase information in the noise environment, stable sparse or compressed phase retrieval needs as many measurements as the relevant compressed sensing problem. Therefore, as with compressed sensing, stable sparse phase retrieval requires at least Oklogn/k measurements [29]. [30] proves that sparse phase retrieval requires Oklogn/k measurements to recover sparse signals stably in the real case.We adopt a real-valued Gaussian model to analyze sparse phase retrieval in this paper. The model assumes sparse signalx∈ℝn and i.i.d. standard Gaussian sensing vector ai~N0,In,i=1,⋯,m. Nevertheless, the proposed algorithm is also applicable to the complex-valued Gaussian phase retrieval model, namely, signal x∈ℂn and i.i.d. standard Gaussian sensing vector ai~CN0,In=N0,In/2+jN0,In/2,i=1,⋯,m. Given the data ai,ψii=1m and assuming that the equation has a uniquek-sparse solution, our goal is to develop a simple and efficient algorithm that can recover any k-sparse n-dimensional signal x from as few amplitude measurements as possible in (1).Using the least square criterion, the problem of reconstructingk-sparse solution from phaseless quadratic equation (1) can be naturally transformed into the problem of minimizing the amplitude-based empirical loss function:
(2)minzℓz∶=12m∑i=1mψi−aiTz2s.t.z0≤kSz⊆n≔1,2,⋯,nz∈ℝn,where Sz represents the support of z and k is the signal sparsity. The objective function ℓz in (2) is nonconvex and nonsmooth, which is subject to the combined constraint of zero norm z0≤k. Therefore, the sparse phase retrieval optimization problem is NP-hard, which is difficult to solve in the calculation [31–33]. The methods to solve these problems include the following:
(1)
In the initialization stage, the support of the signal needs to be accurately estimated, and the accurate support information can correctly restore the original sparse signal. The selection of initialization methods is also critical. At present, there are many initialization methods, including spectral initialization, orthogonal promotion initialization, and reweighted maximum correlation initialization, which can be applied to nonconvex sparse phase retrieval models(2)
In the gradient refinement stage, inspired by the iterative hard threshold algorithm of compressed sensing [34, 35], in the descent process, the adaptive hard threshold process that only retains some maximum indexes in each iteration is proved to be effective


### 2.2. Smoothing Scheme of Sparse Phase Retrieval

Due to the existence of modulus in the sparse phase retrieval objective function (2), there will be dramatic changes and discontinuous changes in the gradient function at point z when aiTz=0, which makes the normal gradient descent method ineffective. Therefore, the sparse optimization problem in (2) is nonconvex and nonsmooth. Usually, the way to solve this kind of nonsmooth function is to bypass this nonsmooth function, such as the SPARTA or CRAF algorithm. [22] points out that if aiTz/ψi is less than a certain threshold, the SPARTA algorithm will remove these “bad” gradient components by the truncation program. The CRAF algorithm [23] is based on the same principle, but it is not a simple removal of these “bad” gradient components, but according to the size of aiTz/ψi to reweight gradient components, the effect is better than the SPARTA algorithm. However, the above two methods will introduce additional calculation for the gradient, resulting in changes in the search direction and thus affecting the update of the entire gradient direction.In recent years, there have been many research on the nonconvex phase retrieval algorithm, among which the optimization of the nonsmooth term in the empirical loss function has made great progress. [36] proposed a mixed optimization method to modify the empirical loss function and solve the nonsmooth problem through the classical proximal method. Pinilla et al. [37] proposed to smooth the phase retrieval problem. They introduced a special smoothing function to replace the nonsmooth term in the original loss empirical function and solved the related problems using the projection conjugate gradient method. [38] adopted the same smoothing strategy. By introducing the smoothing approximation function, it replaced the nonsmooth term and used the same operation to put forward a new smoothing loss empirical function. [39] introduced two smoothing functions for the nonconvex phase retrieval loss function and proposed two smoothing algorithms. [40] further improves the smoothing operation for the loss function and proposes a faster stochastic smoothing phase retrieval algorithm.By introducing the smoothing function, we approximate and replace the absolute value functionaiTx in (2) to obtain an approximate sparse smoothing phase retrieval loss empirical function. Therefore, no additional operation on the gradient is needed in the iterative update process, which reduces the computational complexity. There are many substitution methods for the absolute value function. We mainly use the following method to smooth the phase retrieval loss function (2).We first define the concept of smoothing function to be used.Definition 1.
Define the absolute value functionfu=u. If gεu is smooth for any real constant ε≥0inℂn and for any fixed u∈ℂn(3)limε⟶0gεu=fu,then gεu:ℂn×ℝ⟶ℝ is the smoothing function of fu.By Definition1, the smooth substitution function of the absolute value function fu=u is defined as
(4)gε,pu∶=up+εp1/p,where the relaxation factor ε≥0 is a real constant and the smoothing parameter p>0.As a smooth approximation function of the absolute value functionfu, gε,pu must have the following properties:
(1)
gε,pu is a Lipschitz continuous function(2)
gε,pu is even, namely, gε,pu=gε,p−u(3)
gε,pu uniformly converges to fu in ℝProof.
(1)
By Definition1, gε,pu is smooth for any real constant ε≥0 in ℂn, so we can have(5)∇gε,pu=up−1up+εp1−1/p=uup+εp1/pp−1.
Since the real constantsε≥0, p>0, therefore,
(6)up+εp1/p≥up1/p=u.
According to (5) and (6), it can be obtained that
(7)∇gε,pu=gε,pu1−gε,pu2u1−u2≤1.
Therefore,gε,puis a Lipschitz continuous function with bounded first-order function, and its Lipschitz constant is 1.
(2)
By definition of the absolute value, it is obvious thatgε,pu=gε,p−u=up+εp1/p(3)
According to Definition1 and (4), we can have(8)gε,pu−fu=up+εp1/p−u.
According to the Minkowski inequality, it can be obtained that(9)gε,pu−fu=up+εp1/p−u≤up1/p+εp1/p−u=ε.
Therefore,gε,pu can converge uniformly to fu and only depends on the value of ε.Letu=aiTz, according to (2) and (4), the smoothing sparse phase retrieval loss function can be expressed as
(10)minz12m∑i=1mψi−aiTzε,p2s.t.z0≤kSz⊆n≔1,2,⋯,nz∈ℝn.To ensure that the introduction ofgε,paiTz does not affect the global minimum of the original nonconvex loss function (2), the measurement ψi adopts the same smoothing strategy, i.e.,
(11)ψiε,p∶=ψip+εip1/p.Due to the fact that smoothing vectorεi is a positive correlation to the amplitude ψi of signal x, we approximate it as
(12)εi∶=τψi,where the coefficient τ is a constant.Combining (10) and (11), we can obtain the amplitude-based sparse smoothing phase retrieval least-square optimization problem:
(13)minzℓε,pz∶=12m∑i=1maiTzp+εip1/p−ψip+εip1/p2s.t.z0≤kSz⊆n≔1,2,⋯,nz∈ℝn.Functionℓε,pz is an approximate smoothing function of the original loss function ℓz. Both of ℓε,pz and ℓz have the same global minimum. When the smoothing relaxation factor εi⟶0 or the smoothing parameter p⟶1, ℓε,pz will be reduced to the original function ℓz.


## 2.1. Amplitude-Based Sparse Phase Retrieval

The amplitude-based sparse phase retrieval problem is a kind of problem to reconstruct the sparse signal from the phaseless measurements [19, 20, 22]. Mathematically, we can describe the sparse phase retrieval problem as a set of phaseless quadratic equations, namely,
(1)ψi=aiTx,1≤i≤msubjecttox0≤k,where k represents the level of sparsity. ·0 express a zero-norm operator, that is, the number of nonzero elements. The signal x is k-sparse, consisting of up to k nonzero elements and n−k zero elements. The goal of amplitude-based sparse phase retrieval is to reconstruct the sparse signal x∈ℝn or ℂn based on the given measurements ψi.For theoretical analysis, suppose sparsityk is prior known. In the real case, the number of measurements required to the k-sparse signal x recovery is at least m≥min2k,2n−1. [27] pointed out that k-sparse signals can be reconstructed by selecting m≥min2k,2n−1 measurements on general positions in real number space ℝn. Similarly, in the complex cases, the number of measurements required for k-sparse phase retrieval is at least m≥4k−2 [28]. Due to the lack of more phase information in the noise environment, stable sparse or compressed phase retrieval needs as many measurements as the relevant compressed sensing problem. Therefore, as with compressed sensing, stable sparse phase retrieval requires at least Oklogn/k measurements [29]. [30] proves that sparse phase retrieval requires Oklogn/k measurements to recover sparse signals stably in the real case.We adopt a real-valued Gaussian model to analyze sparse phase retrieval in this paper. The model assumes sparse signalx∈ℝn and i.i.d. standard Gaussian sensing vector ai~N0,In,i=1,⋯,m. Nevertheless, the proposed algorithm is also applicable to the complex-valued Gaussian phase retrieval model, namely, signal x∈ℂn and i.i.d. standard Gaussian sensing vector ai~CN0,In=N0,In/2+jN0,In/2,i=1,⋯,m. Given the data ai,ψii=1m and assuming that the equation has a uniquek-sparse solution, our goal is to develop a simple and efficient algorithm that can recover any k-sparse n-dimensional signal x from as few amplitude measurements as possible in (1).Using the least square criterion, the problem of reconstructingk-sparse solution from phaseless quadratic equation (1) can be naturally transformed into the problem of minimizing the amplitude-based empirical loss function:
(2)minzℓz∶=12m∑i=1mψi−aiTz2s.t.z0≤kSz⊆n≔1,2,⋯,nz∈ℝn,where Sz represents the support of z and k is the signal sparsity. The objective function ℓz in (2) is nonconvex and nonsmooth, which is subject to the combined constraint of zero norm z0≤k. Therefore, the sparse phase retrieval optimization problem is NP-hard, which is difficult to solve in the calculation [31–33]. The methods to solve these problems include the following:
(1)
In the initialization stage, the support of the signal needs to be accurately estimated, and the accurate support information can correctly restore the original sparse signal. The selection of initialization methods is also critical. At present, there are many initialization methods, including spectral initialization, orthogonal promotion initialization, and reweighted maximum correlation initialization, which can be applied to nonconvex sparse phase retrieval models(2)
In the gradient refinement stage, inspired by the iterative hard threshold algorithm of compressed sensing [34, 35], in the descent process, the adaptive hard threshold process that only retains some maximum indexes in each iteration is proved to be effective


## 2.2. Smoothing Scheme of Sparse Phase Retrieval

Due to the existence of modulus in the sparse phase retrieval objective function (2), there will be dramatic changes and discontinuous changes in the gradient function at point z when aiTz=0, which makes the normal gradient descent method ineffective. Therefore, the sparse optimization problem in (2) is nonconvex and nonsmooth. Usually, the way to solve this kind of nonsmooth function is to bypass this nonsmooth function, such as the SPARTA or CRAF algorithm. [22] points out that if aiTz/ψi is less than a certain threshold, the SPARTA algorithm will remove these “bad” gradient components by the truncation program. The CRAF algorithm [23] is based on the same principle, but it is not a simple removal of these “bad” gradient components, but according to the size of aiTz/ψi to reweight gradient components, the effect is better than the SPARTA algorithm. However, the above two methods will introduce additional calculation for the gradient, resulting in changes in the search direction and thus affecting the update of the entire gradient direction.In recent years, there have been many research on the nonconvex phase retrieval algorithm, among which the optimization of the nonsmooth term in the empirical loss function has made great progress. [36] proposed a mixed optimization method to modify the empirical loss function and solve the nonsmooth problem through the classical proximal method. Pinilla et al. [37] proposed to smooth the phase retrieval problem. They introduced a special smoothing function to replace the nonsmooth term in the original loss empirical function and solved the related problems using the projection conjugate gradient method. [38] adopted the same smoothing strategy. By introducing the smoothing approximation function, it replaced the nonsmooth term and used the same operation to put forward a new smoothing loss empirical function. [39] introduced two smoothing functions for the nonconvex phase retrieval loss function and proposed two smoothing algorithms. [40] further improves the smoothing operation for the loss function and proposes a faster stochastic smoothing phase retrieval algorithm.By introducing the smoothing function, we approximate and replace the absolute value functionaiTx in (2) to obtain an approximate sparse smoothing phase retrieval loss empirical function. Therefore, no additional operation on the gradient is needed in the iterative update process, which reduces the computational complexity. There are many substitution methods for the absolute value function. We mainly use the following method to smooth the phase retrieval loss function (2).We first define the concept of smoothing function to be used.Definition 1.
Define the absolute value functionfu=u. If gεu is smooth for any real constant ε≥0inℂn and for any fixed u∈ℂn(3)limε⟶0gεu=fu,then gεu:ℂn×ℝ⟶ℝ is the smoothing function of fu.By Definition1, the smooth substitution function of the absolute value function fu=u is defined as
(4)gε,pu∶=up+εp1/p,where the relaxation factor ε≥0 is a real constant and the smoothing parameter p>0.As a smooth approximation function of the absolute value functionfu, gε,pu must have the following properties:
(1)
gε,pu is a Lipschitz continuous function(2)
gε,pu is even, namely, gε,pu=gε,p−u(3)
gε,pu uniformly converges to fu in ℝProof.
(1)
By Definition1, gε,pu is smooth for any real constant ε≥0 in ℂn, so we can have(5)∇gε,pu=up−1up+εp1−1/p=uup+εp1/pp−1.
Since the real constantsε≥0, p>0, therefore,
(6)up+εp1/p≥up1/p=u.
According to (5) and (6), it can be obtained that
(7)∇gε,pu=gε,pu1−gε,pu2u1−u2≤1.
Therefore,gε,puis a Lipschitz continuous function with bounded first-order function, and its Lipschitz constant is 1.
(2)
By definition of the absolute value, it is obvious thatgε,pu=gε,p−u=up+εp1/p(3)
According to Definition1 and (4), we can have(8)gε,pu−fu=up+εp1/p−u.
According to the Minkowski inequality, it can be obtained that(9)gε,pu−fu=up+εp1/p−u≤up1/p+εp1/p−u=ε.
Therefore,gε,pu can converge uniformly to fu and only depends on the value of ε.Letu=aiTz, according to (2) and (4), the smoothing sparse phase retrieval loss function can be expressed as
(10)minz12m∑i=1mψi−aiTzε,p2s.t.z0≤kSz⊆n≔1,2,⋯,nz∈ℝn.To ensure that the introduction ofgε,paiTz does not affect the global minimum of the original nonconvex loss function (2), the measurement ψi adopts the same smoothing strategy, i.e.,
(11)ψiε,p∶=ψip+εip1/p.Due to the fact that smoothing vectorεi is a positive correlation to the amplitude ψi of signal x, we approximate it as
(12)εi∶=τψi,where the coefficient τ is a constant.Combining (10) and (11), we can obtain the amplitude-based sparse smoothing phase retrieval least-square optimization problem:
(13)minzℓε,pz∶=12m∑i=1maiTzp+εip1/p−ψip+εip1/p2s.t.z0≤kSz⊆n≔1,2,⋯,nz∈ℝn.Functionℓε,pz is an approximate smoothing function of the original loss function ℓz. Both of ℓε,pz and ℓz have the same global minimum. When the smoothing relaxation factor εi⟶0 or the smoothing parameter p⟶1, ℓε,pz will be reduced to the original function ℓz.


## 3. The Proposed SPSAF Algorithm

To solve the smoothing sparse phase retrieval least-square optimization problem (13), we propose a new smoothing sparse phase retrieval algorithm, named SPSAF. It can be classified into two periods: initialization and gradient refinement, which are described in detail.


### 3.1. Sparse Weighted Spectrum Initialization

Based on the weighted initialization method proposed in [41], the SPSAF algorithm proposes an improved sparse weighted initialization estimation method. Compared with the orthogonal promotion initialization method used in the SPARTA algorithm, this method abandons the truncation procedure. It adopts the reweighted strategy so that the information of all samples is effectively applied. We divide the initialization algorithm of SPSAF into three parts, including the general initialization method, support recovery method, and sparse initialization method.


#### 3.1.1. General Initialization Method

First of all, assumex=1 without loss of generality. It can be seen from [25, 41] that there is a certain correlation between ai and x. The larger the amplitude ψi=aiTz of the inner product of ai and x is, the higher the correlation between ai and the unknown solution x is, so it carries more profitable direction information of signal x. Similar to the weighted maximum correlation method in [41], the corresponding sequence of ai can be obtained by sorting the amplitude ψi, and then, the unknown solution x can be estimated.Let setJ≔1,2,⋯,m denote the set of indexes that participate in the calculation initialization of the picking sensing vector ai. The vector ai corresponding to the largest J measurements ψi1≤i≤J have the largest correlation with the real solution x, where J is an integer of m. Therefore, we can estimate the real solution x by indexing ai with J, that is,
(14)z~0≔argmaxz=11J∑i∈JaiTz2=argmaxz=1zT1J∑i∈JaiaiTz.Secondly, to further strengthen the ability ofai to point to the real solution x, the initialization method uses the reweighted strategy to add weights to different aii∈J. Sorting index set J according to the size of ψi1≤i≤J. The larger J corresponds to the larger weight, and the smaller index corresponds to the smaller weight or even negative weight, namely,
(15)z~0≔argmaxz=1zT1J∑i∈JriaiaiTz,where ri∈ℝ is the weight corresponding to different indexes. According to the size of ψi, we define the weight
(16)ri=ψi1/2.Whenx≠1, utilizing the strong law of large numbers and rotation invariance of Gaussian sampling vector ai, we can get the initiation
(17)z0=1m∑i=1mψi2z~0.From Theorem 1 in [23], we can see that given a constant δ0∈0,1, there exist constants c0>0 and C0, and when m≥C0n, the error distz0,x between the initial estimate z0 and the real solution x meets
(18)distz0,x≤δ0x2,with probability exceeding 1−10exp−c0m.Finally, it is worth emphasizing that since the signalx is k-sparse and k≪n, the initialization method in (15) becomes
(19)z~0≔argmaxz=1zT1J∑i∈JriaiaiTzs.t.z0≤k,by utilizing ℓ0 regularization to represent sparse prior information.


#### 3.1.2. Accurate Support Recovery

When the measurements are undersampled, the number of measurements is less than the dimension of the signal; that is,m<n. At this time, the sparse assumption of the signal needs to be introduced.Without loss of generality, we set the supportS⊆n≔1,⋯,n of unknown k-sparse solution x, and S=k≪n. According to [22], we define random variables:
(20)Zi,j∶=ψi2ai,j2,j=1,⋯,n.For the normalized Gaussian sensing vectorai∼N0,In, let p be a positive integer; we can get
(21)Eai,j2p=2p−1!!,where !! denotes two-order multiplication. According to the rotation invariance of Gaussian distribution, it can be proved that for all 1≤j≤n,
(22)EZi,j=EaiTx2ai,j2=Eai,j4xj2+ai,/jTx/j2ai,j2=3xj2+x/j22=2xj2+x22,where ai,/jT,x/j∈ℝn−1 is aiT,x∈ℝn after removing the jth element.If indexj∈S which means the corresponding xj is a nonzero term, that is, xj≠0, then (22) can be converted to
(23)EZi,j=2xj2+x22,j∈S.Conversely, ifj∉S which means xj=0, we can get
(24)Zi,j=x22,j∉S.According to formulas (23) and (24), we can find that for j∉S and j∉S, the expected value of Zi,j has at least 2xj2 intervals, which we call the gap. As long as the gap 2xj2 is large enough, the support S can be restored accurately in this way [22]. That is, the index set corresponding to k maximum EZi,j can restore the support S of the original signal x. However, since EZi,j is actually unavailable, it is replaced by their independent implementation. At the same time, following the law of strong numbers, the sample mean should approach a whole Z~i,j.To estimate the supportS of the original signal x, we first calculate the empirical estimate of the so-called sample mean
(25)Z~i,j∶=1m∑i=1mZi,j=1m∑i=1mψi2ai,j2as the expected EZi,j. The larger Z~i,j, the greater probability of nonzero the corresponding element xj is. Therefore, we need to collect indexes corresponding to k maximums in Z~i,j, which form an estimated support S~ that can be expressed as
(26)S~∶=j∈nindexescorrespondingtokmaximumsinZ~i,j.[22] proves that to improve the probability of accurately recovering the support S, Ok2logn measurements are needed. At the same time, in order to ensure the separation of the index of the support, the smallest nonzero item
(27)xmin∶=minj∈Sxjof the signal is approximately 1/kx2; that is, for some constants C0>0, there are
(28)xmin=C0kx2.


#### 3.1.3. Sparse Initialization Method

After obtaining the support domain estimateS~ of the original signal, we can estimate the initialization according to S~. Specifically, for all 1≤i≤m, we rewrite the measurements ψi as
(29)ψi=aiTx=ai,S~TxS~,where ai,S~ includes the element of ai,j whose index j∈S~, and xS~ is the element of x whose index j∈S~.Apply the general weighted initialization method in (19) to the reduced data ai,S~,ψii=1m, that is,
(30)z~S~0≔argmaxz=1zT1J∑i∈Jriai,S~ai,S~Tzs.t.z0≤k.By zero-filling the element corresponding to the indexj∉S~ of z~S~0, we can construct k-sparse n-dimensional initialization z~0. When x≠1, the final k-sparse initialization z0 can be obtained through the norm estimation of x in (17).


### 3.2. Thresholded Gradient Stage

After obtaining an accuratek-sparse initialization z0, we put it into the gradient iteration for continuous refinement. Our SPSAF algorithm introduces a smoothing strategy for the amplitude-based phase retrieval loss function, avoiding the nondifferentiable objective function. It significantly simplifies the convergence analysis of SPSAF, eliminating the need to perform other operations on the gradient like other algorithms.The method we adopt is to iteratively refine through a series of gradient iterations based on thek-sparse hard threshold, that is,
(31)zt+1∶=Hkzt−μ∇ℓε,pzt,where t is the number of iterations and constant μ is the step size. Hku denotes the k-sparse hard threshold operator, which can transform all elements except k largest elements in u into zero, thus transforming n-dimensional vector u into k-sparse n-dimensional vector u¯, i.e.,
(32)u¯j=uj,ifj∈S,0,ifj∉S.Combining (13), the Wirtinger gradient ∇ℓε,pz of the SPSAF algorithm can be expressed as
(33)∇ℓε,pz=1m∑i=1maiTzp+εip1/p−ψip+εip1/paiTzp+εip1/p−1aiTzp−1aisgnaiTz.It is worth emphasizing that for the smoothed sparse phase retrieval model, the selection of differentp will affect the recovery performance and calculated efficiency of the algorithm. Specifically, the larger the p value, the better the restore performance of the algorithm, but the convergence rate will slow down; the smaller the p value, the faster the algorithm converges, but the recovery performance decreases. We set the p value selection according to the ratio of measurements to sparsity m/k. When the ratio m/k is less than 1.9 in the real case or 2.8 in the complex case, p=2 or 3; otherwise, p=4 or 5. The specific SPSAF algorithm is described as follows.Algorithm 1: The proposed SPSAF algorithm.
1.Input: Data ai,ψii=1m and sparsity k; learning rate μ; the maximum number of iterations T; subset cardinality J≔3m/13; adaptive smoothing factor p; smoothing vector εi≔τψi, and τ is a constant;2.Construct the support index set S~, which includes the index corresponding to the largest k inZ~i,ji=1m≔1/m∑i=1mψi2ai,j2i=1m;3.Calculate the sparse initial estimatez~S~0≔argmaxz=1zT1/J∑i∈Jriai,S~ai,S~Tz,where the weighting factorri=ψi1/2;4.Initializez0=∑i=1mψi2/mz~0wherez~0 is obtained by zero-filling z~S~0;5.fort=0toT−1dozt+1≔Hkzt−μ∇ℓε,pztwhereHk· is the hard threshold operator, and the gradient∇ℓε,pz=1/m∑i=1maiTzp+εip1/p−ψip+εip1/paiTzp+εip1/p−1aiTzp−1aisgnaiTz6.end for7.Output: zT.


## 3.1. Sparse Weighted Spectrum Initialization

Based on the weighted initialization method proposed in [41], the SPSAF algorithm proposes an improved sparse weighted initialization estimation method. Compared with the orthogonal promotion initialization method used in the SPARTA algorithm, this method abandons the truncation procedure. It adopts the reweighted strategy so that the information of all samples is effectively applied. We divide the initialization algorithm of SPSAF into three parts, including the general initialization method, support recovery method, and sparse initialization method.


### 3.1.1. General Initialization Method

First of all, assumex=1 without loss of generality. It can be seen from [25, 41] that there is a certain correlation between ai and x. The larger the amplitude ψi=aiTz of the inner product of ai and x is, the higher the correlation between ai and the unknown solution x is, so it carries more profitable direction information of signal x. Similar to the weighted maximum correlation method in [41], the corresponding sequence of ai can be obtained by sorting the amplitude ψi, and then, the unknown solution x can be estimated.Let setJ≔1,2,⋯,m denote the set of indexes that participate in the calculation initialization of the picking sensing vector ai. The vector ai corresponding to the largest J measurements ψi1≤i≤J have the largest correlation with the real solution x, where J is an integer of m. Therefore, we can estimate the real solution x by indexing ai with J, that is,
(14)z~0≔argmaxz=11J∑i∈JaiTz2=argmaxz=1zT1J∑i∈JaiaiTz.Secondly, to further strengthen the ability ofai to point to the real solution x, the initialization method uses the reweighted strategy to add weights to different aii∈J. Sorting index set J according to the size of ψi1≤i≤J. The larger J corresponds to the larger weight, and the smaller index corresponds to the smaller weight or even negative weight, namely,
(15)z~0≔argmaxz=1zT1J∑i∈JriaiaiTz,where ri∈ℝ is the weight corresponding to different indexes. According to the size of ψi, we define the weight
(16)ri=ψi1/2.Whenx≠1, utilizing the strong law of large numbers and rotation invariance of Gaussian sampling vector ai, we can get the initiation
(17)z0=1m∑i=1mψi2z~0.From Theorem 1 in [23], we can see that given a constant δ0∈0,1, there exist constants c0>0 and C0, and when m≥C0n, the error distz0,x between the initial estimate z0 and the real solution x meets
(18)distz0,x≤δ0x2,with probability exceeding 1−10exp−c0m.Finally, it is worth emphasizing that since the signalx is k-sparse and k≪n, the initialization method in (15) becomes
(19)z~0≔argmaxz=1zT1J∑i∈JriaiaiTzs.t.z0≤k,by utilizing ℓ0 regularization to represent sparse prior information.


### 3.1.2. Accurate Support Recovery

When the measurements are undersampled, the number of measurements is less than the dimension of the signal; that is,m<n. At this time, the sparse assumption of the signal needs to be introduced.Without loss of generality, we set the supportS⊆n≔1,⋯,n of unknown k-sparse solution x, and S=k≪n. According to [22], we define random variables:
(20)Zi,j∶=ψi2ai,j2,j=1,⋯,n.For the normalized Gaussian sensing vectorai∼N0,In, let p be a positive integer; we can get
(21)Eai,j2p=2p−1!!,where !! denotes two-order multiplication. According to the rotation invariance of Gaussian distribution, it can be proved that for all 1≤j≤n,
(22)EZi,j=EaiTx2ai,j2=Eai,j4xj2+ai,/jTx/j2ai,j2=3xj2+x/j22=2xj2+x22,where ai,/jT,x/j∈ℝn−1 is aiT,x∈ℝn after removing the jth element.If indexj∈S which means the corresponding xj is a nonzero term, that is, xj≠0, then (22) can be converted to
(23)EZi,j=2xj2+x22,j∈S.Conversely, ifj∉S which means xj=0, we can get
(24)Zi,j=x22,j∉S.According to formulas (23) and (24), we can find that for j∉S and j∉S, the expected value of Zi,j has at least 2xj2 intervals, which we call the gap. As long as the gap 2xj2 is large enough, the support S can be restored accurately in this way [22]. That is, the index set corresponding to k maximum EZi,j can restore the support S of the original signal x. However, since EZi,j is actually unavailable, it is replaced by their independent implementation. At the same time, following the law of strong numbers, the sample mean should approach a whole Z~i,j.To estimate the supportS of the original signal x, we first calculate the empirical estimate of the so-called sample mean
(25)Z~i,j∶=1m∑i=1mZi,j=1m∑i=1mψi2ai,j2as the expected EZi,j. The larger Z~i,j, the greater probability of nonzero the corresponding element xj is. Therefore, we need to collect indexes corresponding to k maximums in Z~i,j, which form an estimated support S~ that can be expressed as
(26)S~∶=j∈nindexescorrespondingtokmaximumsinZ~i,j.[22] proves that to improve the probability of accurately recovering the support S, Ok2logn measurements are needed. At the same time, in order to ensure the separation of the index of the support, the smallest nonzero item
(27)xmin∶=minj∈Sxjof the signal is approximately 1/kx2; that is, for some constants C0>0, there are
(28)xmin=C0kx2.


### 3.1.3. Sparse Initialization Method

After obtaining the support domain estimateS~ of the original signal, we can estimate the initialization according to S~. Specifically, for all 1≤i≤m, we rewrite the measurements ψi as
(29)ψi=aiTx=ai,S~TxS~,where ai,S~ includes the element of ai,j whose index j∈S~, and xS~ is the element of x whose index j∈S~.Apply the general weighted initialization method in (19) to the reduced data ai,S~,ψii=1m, that is,
(30)z~S~0≔argmaxz=1zT1J∑i∈Jriai,S~ai,S~Tzs.t.z0≤k.By zero-filling the element corresponding to the indexj∉S~ of z~S~0, we can construct k-sparse n-dimensional initialization z~0. When x≠1, the final k-sparse initialization z0 can be obtained through the norm estimation of x in (17).


## 3.1.1. General Initialization Method

First of all, assumex=1 without loss of generality. It can be seen from [25, 41] that there is a certain correlation between ai and x. The larger the amplitude ψi=aiTz of the inner product of ai and x is, the higher the correlation between ai and the unknown solution x is, so it carries more profitable direction information of signal x. Similar to the weighted maximum correlation method in [41], the corresponding sequence of ai can be obtained by sorting the amplitude ψi, and then, the unknown solution x can be estimated.Let setJ≔1,2,⋯,m denote the set of indexes that participate in the calculation initialization of the picking sensing vector ai. The vector ai corresponding to the largest J measurements ψi1≤i≤J have the largest correlation with the real solution x, where J is an integer of m. Therefore, we can estimate the real solution x by indexing ai with J, that is,
(14)z~0≔argmaxz=11J∑i∈JaiTz2=argmaxz=1zT1J∑i∈JaiaiTz.Secondly, to further strengthen the ability ofai to point to the real solution x, the initialization method uses the reweighted strategy to add weights to different aii∈J. Sorting index set J according to the size of ψi1≤i≤J. The larger J corresponds to the larger weight, and the smaller index corresponds to the smaller weight or even negative weight, namely,
(15)z~0≔argmaxz=1zT1J∑i∈JriaiaiTz,where ri∈ℝ is the weight corresponding to different indexes. According to the size of ψi, we define the weight
(16)ri=ψi1/2.Whenx≠1, utilizing the strong law of large numbers and rotation invariance of Gaussian sampling vector ai, we can get the initiation
(17)z0=1m∑i=1mψi2z~0.From Theorem 1 in [23], we can see that given a constant δ0∈0,1, there exist constants c0>0 and C0, and when m≥C0n, the error distz0,x between the initial estimate z0 and the real solution x meets
(18)distz0,x≤δ0x2,with probability exceeding 1−10exp−c0m.Finally, it is worth emphasizing that since the signalx is k-sparse and k≪n, the initialization method in (15) becomes
(19)z~0≔argmaxz=1zT1J∑i∈JriaiaiTzs.t.z0≤k,by utilizing ℓ0 regularization to represent sparse prior information.


## 3.1.2. Accurate Support Recovery

When the measurements are undersampled, the number of measurements is less than the dimension of the signal; that is,m<n. At this time, the sparse assumption of the signal needs to be introduced.Without loss of generality, we set the supportS⊆n≔1,⋯,n of unknown k-sparse solution x, and S=k≪n. According to [22], we define random variables:
(20)Zi,j∶=ψi2ai,j2,j=1,⋯,n.For the normalized Gaussian sensing vectorai∼N0,In, let p be a positive integer; we can get
(21)Eai,j2p=2p−1!!,where !! denotes two-order multiplication. According to the rotation invariance of Gaussian distribution, it can be proved that for all 1≤j≤n,
(22)EZi,j=EaiTx2ai,j2=Eai,j4xj2+ai,/jTx/j2ai,j2=3xj2+x/j22=2xj2+x22,where ai,/jT,x/j∈ℝn−1 is aiT,x∈ℝn after removing the jth element.If indexj∈S which means the corresponding xj is a nonzero term, that is, xj≠0, then (22) can be converted to
(23)EZi,j=2xj2+x22,j∈S.Conversely, ifj∉S which means xj=0, we can get
(24)Zi,j=x22,j∉S.According to formulas (23) and (24), we can find that for j∉S and j∉S, the expected value of Zi,j has at least 2xj2 intervals, which we call the gap. As long as the gap 2xj2 is large enough, the support S can be restored accurately in this way [22]. That is, the index set corresponding to k maximum EZi,j can restore the support S of the original signal x. However, since EZi,j is actually unavailable, it is replaced by their independent implementation. At the same time, following the law of strong numbers, the sample mean should approach a whole Z~i,j.To estimate the supportS of the original signal x, we first calculate the empirical estimate of the so-called sample mean
(25)Z~i,j∶=1m∑i=1mZi,j=1m∑i=1mψi2ai,j2as the expected EZi,j. The larger Z~i,j, the greater probability of nonzero the corresponding element xj is. Therefore, we need to collect indexes corresponding to k maximums in Z~i,j, which form an estimated support S~ that can be expressed as
(26)S~∶=j∈nindexescorrespondingtokmaximumsinZ~i,j.[22] proves that to improve the probability of accurately recovering the support S, Ok2logn measurements are needed. At the same time, in order to ensure the separation of the index of the support, the smallest nonzero item
(27)xmin∶=minj∈Sxjof the signal is approximately 1/kx2; that is, for some constants C0>0, there are
(28)xmin=C0kx2.


## 3.1.3. Sparse Initialization Method

After obtaining the support domain estimateS~ of the original signal, we can estimate the initialization according to S~. Specifically, for all 1≤i≤m, we rewrite the measurements ψi as
(29)ψi=aiTx=ai,S~TxS~,where ai,S~ includes the element of ai,j whose index j∈S~, and xS~ is the element of x whose index j∈S~.Apply the general weighted initialization method in (19) to the reduced data ai,S~,ψii=1m, that is,
(30)z~S~0≔argmaxz=1zT1J∑i∈Jriai,S~ai,S~Tzs.t.z0≤k.By zero-filling the element corresponding to the indexj∉S~ of z~S~0, we can construct k-sparse n-dimensional initialization z~0. When x≠1, the final k-sparse initialization z0 can be obtained through the norm estimation of x in (17).


## 3.2. Thresholded Gradient Stage

After obtaining an accuratek-sparse initialization z0, we put it into the gradient iteration for continuous refinement. Our SPSAF algorithm introduces a smoothing strategy for the amplitude-based phase retrieval loss function, avoiding the nondifferentiable objective function. It significantly simplifies the convergence analysis of SPSAF, eliminating the need to perform other operations on the gradient like other algorithms.The method we adopt is to iteratively refine through a series of gradient iterations based on thek-sparse hard threshold, that is,
(31)zt+1∶=Hkzt−μ∇ℓε,pzt,where t is the number of iterations and constant μ is the step size. Hku denotes the k-sparse hard threshold operator, which can transform all elements except k largest elements in u into zero, thus transforming n-dimensional vector u into k-sparse n-dimensional vector u¯, i.e.,
(32)u¯j=uj,ifj∈S,0,ifj∉S.Combining (13), the Wirtinger gradient ∇ℓε,pz of the SPSAF algorithm can be expressed as
(33)∇ℓε,pz=1m∑i=1maiTzp+εip1/p−ψip+εip1/paiTzp+εip1/p−1aiTzp−1aisgnaiTz.It is worth emphasizing that for the smoothed sparse phase retrieval model, the selection of differentp will affect the recovery performance and calculated efficiency of the algorithm. Specifically, the larger the p value, the better the restore performance of the algorithm, but the convergence rate will slow down; the smaller the p value, the faster the algorithm converges, but the recovery performance decreases. We set the p value selection according to the ratio of measurements to sparsity m/k. When the ratio m/k is less than 1.9 in the real case or 2.8 in the complex case, p=2 or 3; otherwise, p=4 or 5. The specific SPSAF algorithm is described as follows.Algorithm 1: The proposed SPSAF algorithm.
1.Input: Data ai,ψii=1m and sparsity k; learning rate μ; the maximum number of iterations T; subset cardinality J≔3m/13; adaptive smoothing factor p; smoothing vector εi≔τψi, and τ is a constant;2.Construct the support index set S~, which includes the index corresponding to the largest k inZ~i,ji=1m≔1/m∑i=1mψi2ai,j2i=1m;3.Calculate the sparse initial estimatez~S~0≔argmaxz=1zT1/J∑i∈Jriai,S~ai,S~Tz,where the weighting factorri=ψi1/2;4.Initializez0=∑i=1mψi2/mz~0wherez~0 is obtained by zero-filling z~S~0;5.fort=0toT−1dozt+1≔Hkzt−μ∇ℓε,pztwhereHk· is the hard threshold operator, and the gradient∇ℓε,pz=1/m∑i=1maiTzp+εip1/p−ψip+εip1/paiTzp+εip1/p−1aiTzp−1aisgnaiTz6.end for7.Output: zT.


## 4. Experimental Results

This section will introduce the relevant numerical results of the SPSAF algorithm. To reflect the superiority of the SPSAF algorithm, we compare it with the SPARTA algorithm [16] and CRAF algorithm [20], which are the latest methods of sparse phase retrieval. The parameters of all algorithms will use their recommended values. All simulation experiments were conducted with 100 independent Monte Carlo experiments. In each experiment, all algorithms’ initialization power iteration numbers and gradient refinement iteration numbers are set to 100.In all experiments, the signalx∈ℝ3000 or ℂ3000 to be recovered is k sparse signal in the real or complex case. When the sparse phase retrieval adopts the real-valued Gaussian model, the sparse signal is x~N0,I3000 and the sensing vector is ai∼N0,I3000,i=1,⋯,m. When sparse phase retrieval adopts the complex Gaussian model, the sparse signal is x~CN0,I3000=N0,1/2I3000+jN0,1/2I3000 and the sensing vector is ai=N0,1/2I3000+jN0,1/2I3000,i=1,⋯,m. In addition, other parameters of the SPSAF algorithm are selected according to experience: learning rate μ=0.8 in the real case or μ=1 in the complex case, smoothing parameter p=2.We use relative error as our performance index, that is,(34)Relativeerror∶=distz,xx2,where distz,x is the Euclidean distance from the estimated value z to the real solution x. When the relative error is less than 10−5, it can be considered that the original sparse signal has been successfully recovered. Namely, the current experiment is successful.


### 4.1. Recovery Success Rate with Unknown Sparsityk

In some practical applications, the sparsityk of the signal may be unknown, so it is necessary to test the running state of the algorithm in the case of unknown sparsity, namely, the algorithm’s stability. We define k~ as the estimation of sparsity k and set it to the upper bound of sparsity level in theory. Theorem 1 in [16] proves that when m is approximately equal to n, the upper limit of the sparsity level in this paper is about n=55.When the phase retrieval model is in the real case, the recovery success rates of SPSAF, SPARTA, and CRAF algorithms are compared as shown in Figure1, where the number of measurements m increases from 400 to 2400. Notably, these curves show that SPSAF has higher precision recovery performance than other comparison algorithms. When the sparsity k is unknown, the SPSAF algorithm shows higher stability, and the recovery rate reaches more than 90% when m=1000 and achieves accurate recovery when m=1400.Figure 1
The success rate of signal recovery with sparsityk=30 and estimated sparsity k~=55 in the real case.When the phase retrieval model is in the complex case, the recovery success rates of SPSAF, SPARTA, and CRAF algorithms are shown in Figure2, where the number of measurements m ranges from 800 to 2800. It can be found from Figure 2 that the SPSAF algorithm has a weak advantage in recovery ability compared with other algorithms. When the actual sparsity k=30, the estimated sparsity k~=55; the SPSAF algorithm is stable. It is worth emphasizing that when the number of measurements m=500, the SPSAF algorithm first reaches more than 80%, showing excellent recovery ability.Figure 2
The success rate of signal recovery with sparsityk=30 and estimated sparsity k~=55 in the complex case.


### 4.2. Recovery Success Rate with Known Sparsityk

In this section, we compare the signal recovery success rates of each algorithm when the sparsityk is known. The comparison of the algorithms in the real case is shown in Figure 3, where the number of measurements m increases from 400 to 2400. It can be seen from Figure 3 that the recovery performance of the SPSAF algorithm is significantly higher than the other two algorithms, showing excellent recovery ability. In addition, it is worth noting that compared with Figures 1 and 3, it can be found that the recovery performance of the SPSAF algorithm significantly improved when the sparse prior is known.Figure 3
The success rate of signal recovery with prior sparsityk=30 in the real case.Figure4 depicts the success rate of signal recovery in the complex case. We set the distribution of the number of measurements m from 800 to 2800, where the prior sparsity k=30 is known. It can be seen from Figure 4 that the SPSAF algorithm is significantly better than the other two algorithms. Moreover, by comparing Figures 2 and 4, we find that the recovery degree of the algorithm has been slightly improved when the sparsity k is known.Figure 4
The success rate of signal recovery with prior sparsityk=30 in the complex case.Figure5 depicts the relationship between the success rate of signal recovery and the sparsity k in the real case, where the step size parameters of each algorithm take the optimal value. It can be seen that with the increasing sparsity k, the SPSAF algorithm can still accurately restore the signal at the sparsity k=50 and can still maintain a success rate of more than 95% at k=100. In comparison, the success rates of the CRAF algorithm and SPARTA algorithm with the same sparsity are reduced to about 80%.Figure 5
Comparison of signal recovery success rate under different sparsityk selection.


### 4.3. Comparison of Convergence Consumption

This section compares the SPSAF algorithm with other algorithms in terms of convergence speed and time cost. We adopt the real Gaussian model when the number of measurements ism=1600 or the complex Gaussian model when m=2800, respectively. The prior sparsity of all models is k=30, and the relative error of the algorithm is less than 10−15 which is regarded as an accurate recovery. The convergence curves in the real or complex case are shown in Figures 6 and 7. It can be seen from these figures that the number of iterations required for the convergence of the SPSAF algorithm is significantly less than that of the other two algorithms.Figure 6
Comparison of convergence rate in the real case.Figure 7
Comparison of convergence rate in the complex case.On the other hand, Table1 compares convergence speed and time cost under the noise-free Gaussian model, where the coarse font is the current optimal value. It can be seen from the results that the convergence speed and time cost of the SPSAF algorithm are better than those of other algorithms in both real and complex cases, which show excellent performance.Table 1
Comparison of convergence speed and time cost.
AlgorithmsReal caseComplex caseIterationsTime (s)IterationsTime (s)SPSAF260.1035751.3525CRAF640.2013941.6509SPARTA270.1352961.8528


### 4.4. Noise Robustness

In order to prove the stability of the SPSAF algorithm in the presence of additional noise, we plot the relative mean square error (MSE) as a function of the signal-to-noise ratio (SNR) value in dB. Amplitude measurements with noise can be expressed as(35)ψi=aiTx+ηi,where the size of ηi~N0,σ2I is given by
(36)SNR∶=10log10∑i=1maiTx2mσ2.In this section, we adoptm=1600 real-valued Gaussian model with sparse prior k=30. Figure 8 depicts the relationship between the average relative error of the three algorithms when σ2=3. It can be seen that the SPSAF algorithm provides the most accurate estimation under the noise addition.Figure 8
Relative error of additional noise in the real case.We also describe the relative error as a function of SNR with differentm/n to verify the stability of the SPSAF algorithm. Figure 9 depicts the experimental results. We use a real-valued Gaussian model with sparse prior k=30, where SNR is between 5 dB and 50 dB. It can be seen that the SPSAF algorithm decreases approximately linearly with the increase of SNR in both real and complex cases, which proves that the proposed SPSAF algorithm is stable in noisy environments.Figure 9
Relationship between relative error and SNR.


## 4.1. Recovery Success Rate with Unknown Sparsityk

In some practical applications, the sparsityk of the signal may be unknown, so it is necessary to test the running state of the algorithm in the case of unknown sparsity, namely, the algorithm’s stability. We define k~ as the estimation of sparsity k and set it to the upper bound of sparsity level in theory. Theorem 1 in [16] proves that when m is approximately equal to n, the upper limit of the sparsity level in this paper is about n=55.When the phase retrieval model is in the real case, the recovery success rates of SPSAF, SPARTA, and CRAF algorithms are compared as shown in Figure1, where the number of measurements m increases from 400 to 2400. Notably, these curves show that SPSAF has higher precision recovery performance than other comparison algorithms. When the sparsity k is unknown, the SPSAF algorithm shows higher stability, and the recovery rate reaches more than 90% when m=1000 and achieves accurate recovery when m=1400.Figure 1
The success rate of signal recovery with sparsityk=30 and estimated sparsity k~=55 in the real case.When the phase retrieval model is in the complex case, the recovery success rates of SPSAF, SPARTA, and CRAF algorithms are shown in Figure2, where the number of measurements m ranges from 800 to 2800. It can be found from Figure 2 that the SPSAF algorithm has a weak advantage in recovery ability compared with other algorithms. When the actual sparsity k=30, the estimated sparsity k~=55; the SPSAF algorithm is stable. It is worth emphasizing that when the number of measurements m=500, the SPSAF algorithm first reaches more than 80%, showing excellent recovery ability.Figure 2
The success rate of signal recovery with sparsityk=30 and estimated sparsity k~=55 in the complex case.


## 4.2. Recovery Success Rate with Known Sparsityk

In this section, we compare the signal recovery success rates of each algorithm when the sparsityk is known. The comparison of the algorithms in the real case is shown in Figure 3, where the number of measurements m increases from 400 to 2400. It can be seen from Figure 3 that the recovery performance of the SPSAF algorithm is significantly higher than the other two algorithms, showing excellent recovery ability. In addition, it is worth noting that compared with Figures 1 and 3, it can be found that the recovery performance of the SPSAF algorithm significantly improved when the sparse prior is known.Figure 3
The success rate of signal recovery with prior sparsityk=30 in the real case.Figure4 depicts the success rate of signal recovery in the complex case. We set the distribution of the number of measurements m from 800 to 2800, where the prior sparsity k=30 is known. It can be seen from Figure 4 that the SPSAF algorithm is significantly better than the other two algorithms. Moreover, by comparing Figures 2 and 4, we find that the recovery degree of the algorithm has been slightly improved when the sparsity k is known.Figure 4
The success rate of signal recovery with prior sparsityk=30 in the complex case.Figure5 depicts the relationship between the success rate of signal recovery and the sparsity k in the real case, where the step size parameters of each algorithm take the optimal value. It can be seen that with the increasing sparsity k, the SPSAF algorithm can still accurately restore the signal at the sparsity k=50 and can still maintain a success rate of more than 95% at k=100. In comparison, the success rates of the CRAF algorithm and SPARTA algorithm with the same sparsity are reduced to about 80%.Figure 5
Comparison of signal recovery success rate under different sparsityk selection.


## 4.3. Comparison of Convergence Consumption

This section compares the SPSAF algorithm with other algorithms in terms of convergence speed and time cost. We adopt the real Gaussian model when the number of measurements ism=1600 or the complex Gaussian model when m=2800, respectively. The prior sparsity of all models is k=30, and the relative error of the algorithm is less than 10−15 which is regarded as an accurate recovery. The convergence curves in the real or complex case are shown in Figures 6 and 7. It can be seen from these figures that the number of iterations required for the convergence of the SPSAF algorithm is significantly less than that of the other two algorithms.Figure 6
Comparison of convergence rate in the real case.Figure 7
Comparison of convergence rate in the complex case.On the other hand, Table1 compares convergence speed and time cost under the noise-free Gaussian model, where the coarse font is the current optimal value. It can be seen from the results that the convergence speed and time cost of the SPSAF algorithm are better than those of other algorithms in both real and complex cases, which show excellent performance.Table 1
Comparison of convergence speed and time cost.
AlgorithmsReal caseComplex caseIterationsTime (s)IterationsTime (s)SPSAF260.1035751.3525CRAF640.2013941.6509SPARTA270.1352961.8528


## 4.4. Noise Robustness

In order to prove the stability of the SPSAF algorithm in the presence of additional noise, we plot the relative mean square error (MSE) as a function of the signal-to-noise ratio (SNR) value in dB. Amplitude measurements with noise can be expressed as(35)ψi=aiTx+ηi,where the size of ηi~N0,σ2I is given by
(36)SNR∶=10log10∑i=1maiTx2mσ2.In this section, we adoptm=1600 real-valued Gaussian model with sparse prior k=30. Figure 8 depicts the relationship between the average relative error of the three algorithms when σ2=3. It can be seen that the SPSAF algorithm provides the most accurate estimation under the noise addition.Figure 8
Relative error of additional noise in the real case.We also describe the relative error as a function of SNR with differentm/n to verify the stability of the SPSAF algorithm. Figure 9 depicts the experimental results. We use a real-valued Gaussian model with sparse prior k=30, where SNR is between 5 dB and 50 dB. It can be seen that the SPSAF algorithm decreases approximately linearly with the increase of SNR in both real and complex cases, which proves that the proposed SPSAF algorithm is stable in noisy environments.Figure 9
Relationship between relative error and SNR.


## 5. Conclusion

In this paper, we proposed the SPSAF algorithm to solve the problems of sparse phase retrieval. The proposed SPSAF algorithm is an amplitude-based nonconvex sparse smoothing phase retrieval algorithm divided into two stages: initialization and gradient refinement. The complexity of our algorithm is reduced by smoothing the phase retrieval loss function to avoid the truncation or weighting of the gradient. The SPSAF algorithm first estimates the support of the original signal by a reasonable rule, obtains the initial estimation by a carefully designed initialization method based on the support, and finally obtains the sparse estimation by a gradient descent method based on the hard threshold. Numerical experiments show that the SPSAF algorithm has significantly improved recovery performance and speed compared with the existing typical algorithms and has good robustness.



---

*Source: 7012261-2022-06-09.xml*