On Fienup Methods for Regularized Phase Retrieval
Edouard Pauwels Amir Beck Yonina C. Eldar Shoham Sabach














 
 
 
 


todo
 
 TODO:







theoremTheorem[section]
lemma[theorem]Lemma
proposition[theorem]Proposition
corrolary[theorem]Corollary
definitionDefinition
exampleExample
assumptionAssumption
remarkRemark
corollary[theorem]Corollary
proof[1][]
  * Proof #1: 

IEEEtran































 
 
 

 
 
 




 
 




 


 














On Fienup Methods for Regularized Phase Retrieval
    Edouard Pauwels, Amir Beck, Yonina C. Eldar, Fellow, IEEE, Shoham Sabach
The work of Edouard Pauwels was partially supported by the Air Force Office of Scientific Research grant number FA9550-15-1-0500. The research of Amir Beck was partially supported by the Israel Science Foundation Grant 1821/16.
Edouard Pauwels is with the Informatics department (IRIT), Université Toulouse 3 Paul Sabatier, Toulouse 31062, France (e-mail epauwels@irit.fr).
A. Beck is with the department of Industrial Engineering, Technion–Israel Institute of Technology, Haifa, Israel 32000 (e-mail: becka@ie.technion.ac.il).
Y. C. Eldar is with the department of Electrical Engineering, Technion–Israel Institute of Technology, Haifa, Israel 32000 (e-mail: yonina@ee.technion.ac.il).
S. Sabach is with the department of Industrial Engineering, Technion–Israel Institute of Technology, Haifa, Israel 32000 (e-mail: ssabach@ie.technion.ac.il).

    Draft of December 30, 2023
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================


Alternating minimization, or Fienup methods, have a long history in phase retrieval. We provide new insights related to the empirical and theoretical analysis of these algorithms when used with Fourier measurements and combined with convex priors.  In particular, we show that Fienup methods can be viewed as performing alternating minimization on a regularized nonconvex least-squares problem with respect to amplitude measurements. We then prove that under mild additional structural assumptions on the prior (semi-algebraicity), the sequence of signal estimates has a smooth convergent behaviour towards a critical point of the nonconvex regularized least-squares objective. Finally, we propose an extension to Fienup techniques, based on a projected gradient descent interpretation and acceleration using inertial terms. We demonstrate experimentally that this modification combined with an ℓ_1 prior constitutes a competitive approach for sparse phase retrieval.



§ INTRODUCTION


Phase retrieval is an old and fundamental problem in a variety of areas within engineering and physics <cit.>. Many applications of the phase retrieval problem involve estimation of a signal from the modulus of its Fourier measurements. This problem is ill posed in general, so that uniqueness and recovery typically require prior knowledge on the input, particularly in one-dimensional problems. Here we focus on the estimation of real sparse signals from their Fourier magnitude, a problem which has been treated in several recent works <cit.>.

A longstanding line of algorithms to tackle the phase retrieval problem involve application of the alternating minimization method which alternate between the constraints in time and the Fourier magnitude constraints
<cit.>. These methods were pioneered by the work of Gerchberg
and Saxton and later extended by Fienup; see <cit.> for an optimization point of view on these techniques and a rich historical perspective.
Alternating minimization approaches have also been recently applied to phase retrieval from random measurements <cit.>. The main advantage of this class of algorithms is their simplicity and scalability.

A more recent approach to phase retrieval is to formulate the recovery as a smooth nonconvex least-squares estimation problem and use dedicated techniques to estimate the signal using continuous optimization algorithms that guarantee convergence to stationary points. The GESPAR algorithm <cit.> is an example of this approach which is based on the Gauss-Newton method coupled with sparsity priors.
For phase retrieval with random measurements, gradient descent methods have been proposed and analyzed such as Wirtinger flow <cit.> and truncated amplitude flow <cit.>. Both treat least-squares objectives where Wirtinger flow measures the loss with respect to the squared-magnitude of the measurements while the amplitude flow approach performs a truncated gradient descent on an amplitude objective.
Another line of work suggests the use of matrix lifting and semidefinite programming based relaxations <cit.>. These techniques are limited by the size of problems that can be tackled using available numerical solvers.

Our main contribution is to propose a new look at alternating minimization algorithms for phase retrieval in the context of Fourier measurements and convex priors. We refer collectively to these techniques as Fienup methods.
The use of Fourier measurements is less flexible than general measurements and is less suited for statistical analysis. On the other hand, the Fourier transform has very strong structure which allows for richer algorithmic constructions and analysis.

 As a first step we provide two new interpretations of Fienup algorithms. First we show that these techniques are naturally linked to a nonsmooth nonconvex least-squares problem with respect to an amplitude objective.
 Fienup approaches can then be understood as majorization-minimization methods for solving this problem. Second, we demonstrate that Fienup algorithms can be viewed as a projected gradient descent scheme to minimize a smooth convex objective function over a nonconvex constraint set. This observation  allows to characterize the behaviour of the algorithm and develop extensions based on known ideas for accelerating gradient methods using inertial terms <cit.>. We then specialize these results to the case of ℓ_0 and ℓ_1 priors, leading to a new inertial gradient scheme, which we refer to as FISTAPH: FISTA for PHase retrieval.

On the theoretical side, we show that if the convex prior is well structured (semi-algebraic or more generally representable), then the sequence of signal estimates produced by Fienup has a smooth convergence behaviour. Recall that, broadly speaking, an object is said to be semi-algebraic if it can be represented by systems of polynomial inequalities. The notion of smooth convergence is a very desirable property, even more in nonconvex settings where it is usually not possible to obtain global convergence estimates. The convergence analysis follows well established techniques from tame optimization <cit.>. These techniques build upon the Kurdyka-Łojasiewicz (KL) property which holds for many classes of functions <cit.>. We then provide numerical experiments based on synthetic problems to compare Fienup with ℓ_0 and ℓ_1 priors, GESPAR <cit.>, Wirtinger flow (or gradient) methods <cit.> with ℓ_0 and ℓ_1 priors and FISTAPH. Numerical results suggest that the latter combined with an ℓ_1 prior constitutes a very competitive alternative for sparse phase retrieval.

The rest of the paper is organized as follows. Section <ref> introduces our notation and states the problem of interest more formally. We also introduce several mathematical definitions that are required for the rest of the paper and review the numerical algorithms that are used in subsequent sections. Section <ref> describes our  characterization of Fienup methods in the context of phase retrieval from Fourier measurements with convex priors. We detail the relation of Fienup with a nonsmooth nonconvex least-squares problem as well as its interpretation as projected gradient descent. Our main convergence result and our new FISTAPH algorithm are presented in Section <ref>. Simulation results are provided in Section <ref>.



§ PROBLEM FORMULATION AND MATHEMATICAL BACKGROUND
 


 §.§ Notation

Throughout the paper vectors are denoted by boldface letters. For a vector ∈^n, [i] is the i-th entry of , i=1,2,…, n and () is the support of , namely, the set { i = 1, 2, … n; [i] ≠ 0 }. Furthermore, _0 counts the number of nonzero entries of the vector : _0 = |()| and _p denotes the ℓ_p norm of  for p ∈_+. The notations |·|, (·), (·) and ·̅ describe the modulus, real part, imaginary part and complex conjugate, respectively, defined over the field of complex numbers. If their argument is a vector, then they should be understood component-wise. Similarly, basic operations, e.g. powers, are taken component-wise when applied to vectors. For ∈^n and N ∈, (, N) ∈^N is the vector composed of the N first coefficients of the discrete Fourier transform of  (obtained by  zero padding if n <N). For simplicity, we use the shorthand notation () = (,n) to denote the standard discrete Fourier transform of ∈^n and ^-1 to denote its inverse. For a set S, δ_S  S →∪{+∞} is the indicator function of S (0 if its argument is in S, +∞ otherwise) and P_S denotes the Euclidean orthogonal projection onto the set S.


 §.§ Phase Retrieval

Given _0 ∈^n, we consider the data acquisition process

    = |(_0)| + ,

where ∈^n is an unknown vector of errors. In the rest of the paper, we actually assume that  has positive entries (it is always possible to set the potential negative entries of  to zero). The phase retrieval problem consists of producing an estimate ∈^n of _0 based solely on the knowledge of  given by (<ref>).

As mentioned in the introduction, phase retrieval of one-dimensional vectors from Fourier measurements requires the use of prior knowledge. We focus on support and sparsity inducing priors. For J ⊆{ 1,2,…,n }, we define the set X_J = {∈^n; () ⊆ J }. The prior function that we use will be denoted by g ^n →∪{+∞}. We focus on the following priors (for a given J):

	
  * g ↦_0 + δ_X_J(), or ℓ_0-based nonconvex prior.
	
  * g ↦_1 + δ_X_J(), or ℓ_1-based convex prior.

In the experimental section, we compare between these two classes of priors. The algorithmic derivations in this paper will be made under the assumption that g is proper and lower semicontinuous, and the main convergence result (c.f. Theorem <ref>) will require in addition convexity of g.
                In order to efficiently implement the proposed algorithm, we need to focus on priors for which the proximity operator <cit.> is easy to compute. We provide several examples of such priors in Section <ref>.

In the rest of the paper, ∈^n_+ denotes modulus measurements which are assumed to be given, fixed and obtained through (<ref>). Given ∈^n_+, we define Z_ = {∈^n;  |()| = } as the set of values  that could have produced  (ignoring the noise). To estimate _0, we consider the regularized least-squares problem

    min_∈^n, ∈ Z_1/2 - _2^2 + g(),

where g encodes our prior knowledge. Our algorithmic approach consists of employing an alternating minimization method, or one of its variants, to solve the above formulation.



 §.§ Prior Algorithms for Phase Retrieval


We briefly review several existing algorithms for phase retrieval that will be used in our experiments in Section <ref>.

One approach to sparse phase retrieval is the GESPAR algorithm which is based on the damped Gauss-Newton method in conjunction with an ℓ_0 prior <cit.>. Damped Gauss-Newton allows to solve smooth, nonlinear least-squares problems. The work of <cit.> is based on the notion of Wirtinger derivatives to treat the same smooth least-squares problem as GESPAR. The notion of Wirtinger derivative is needed since the objective is not differentiable (holomorphic) as a function of complex variables (see <cit.> for details). In the case of real valued functions of real variables, the Wirtinger derivative reduces to a standard gradient (up to a constant multiplicative factor).
An obvious extension of these types of methods is the use of proximal decomposition, or forward-backward methods which consist in alternating a gradient step on the smooth part of the objective with a proximal step on the nonsmooth part <cit.>. This is the approach that we use in the numerical experiments to treat phase retrieval with priors.

Finally, we consider alternating minimization methods that are the main focus of this work. This approach consists of solving (<ref>) by applying the alternating minimization algorithm. The special structure of the problem allows to perform each partial minimization efficiently. In particular, the projection onto Z_ is easy, as described below in (<ref>). These types of methods are also referred to as Fienup algorithms. A deeper interpretation of this approach is given in Section <ref>.




 §.§ Tools from Convex and Nonsmooth Analysis


Throughout the paper, our results will be based on tools from convex and nonsmooth analysis which we review here.

The gradient of a differentiable function f is denoted by ∇ f. This concept admits extensions to nonsmooth analysis; the subgradient of a nonsmooth function g is denoted by ∂ g. For convex functions, subgradients correspond to tangent affine lower bounds. This definition no longer holds for nonconvex functions. In this case, the proper understanding of subgradients involves much more machinery which will not be discussed here. We only consider the notion of a Fréchet critical point which generalizes classical first order criticality for differentiable functions (see <cit.>).

	
	Let S ⊆^n be a closed set and f^n → be a lower semicontinuous function. We say that ∈ S is a Fréchet critical point of the problem
	
    min_∈ S f()

	if

    c
    				x →x̅
    
    				x ≠x̅
    
    				x ∈ S
    			liminf     f(x) - f(x̅)/x - x̅≥ 0.





	In other words, the negative variations of f in S around x̅ are negligible at the first order.

We will also heavily use the notion of the proximity operator of a function.

	For a nonsmooth function g: ^n →∪{+∞}, the (potentially multivalued) proximity operator is denoted by _g and defined by
	
    _g() ≡_∈^n{1/2 - _2^2 + g()}.


Note that when g is proper lower semicontinuous and convex, this operator is single valued.

We next provide a few examples of such functions with their proximity operators; many more can be found, for example, in <cit.>.
[Proximity operators]

	
		
  * Support prior: If C ⊆^n is a closed convex set, then _δ_C is the Euclidean projection onto C. This can be used for example to encode knowledge about the support of the signal _0 by choosing C = X_J for some J ⊆{ 1, 2, …, n }. In this case, the projection simply consists in setting the coefficients [i] to 0 for i ∉J.
		
  * Sparsity prior: If g is the ℓ_1 norm, then the proximal operator is the soft thresholding operator.[ The soft thresholding operator is given by 𝒯_α()_i =  sgn(x_i)max{ |x_i|-α,0}. If g() = λ_1 for some λ>0, then _g() = 𝒯_λ().]  This can be combined with support information prior by first setting the coefficients outside of the support to 0 and then applying the soft thresholding operator.
		
  * Change of basis: Suppose that D is an n × n' real matrix such that its columns form an orthonormal family, that is D^T D is the identity in ^n'. Suppose that g̃^n'→ is a lower semicontinuous convex function and let g()=g̃(D^T). In this case, we have _g( )=  + D( _g̃(D^T) - D^T) (see <cit.>). This allows to express priors in different orthonormal bases, such as wavelets for example.
	

It is also worth mentioning that the proximity operator is efficiently computable for some nonconvex priors. For example, if g = δ_C where C = {∈^n; _0 ≤ k }, then the proximity operator is obtained by setting the n-k lowest coefficients (in absolute value) to 0. This can also be combined with support information.




§ FIENUP, MAJORIZATION-MINIMIZATION AND PROJECTED GRADIENT



In this section we expand on the alternating minimization approach to (<ref>) leading to the Fienup family of algorithms. For this section, the prior term g in (<ref>) is taken to be a general proper lower semicontinuous   function. We begin by describing the algorithm and then provide two interpretations of it.


 §.§ Alternating Minimization Algorithm

The alternating minimization algorithm applied to problem (<ref>) is explicitly written below.

7cmAlternating Minimization (Fienup) 

	Initialization. ^0∈^n

  General Step. For k ∈,
	
    ^k+1   ∈_∈ Z_1/2^k - _2^2,
    ^k+1   ∈_∈^n1/2 - ^k+1_2^2 + g().

	



The main interest in this scheme is that both partial minimization steps in (<ref>) can be carried out efficiently whenever g is “proximable", meaning that its prox  (or a member in its prox) is easily computed. First consider, in (<ref>), the partial minimization in  with ∈^n being arbitrary but fixed. This minimization amounts to computing P_Z_(^k), the orthogonal  projection of ^k onto Z_. For a given ∈^n, all the members in P_Z_() are of the form = ^-1(), where for j=1,2,…, n, we have (i=√(-1) in the equation below)

    [j] =
    	[j]()[j]/|()[j]|,		
    		   if  |()[j]| ≠ 0, 
    [j]e^iθ_j,    for an arbitrary θ_j  otherwise.

Next, we treat the subproblem in (<ref>) of minimizing with respect to  where ∈^n is arbitrary but fixed. The partial minimization in  is given by the expression

    argmin_∈^n{1/2 - _2^2 + g()} = _g( Re()),

where Re is the real part taken component-wise. We have used the definition of the proximity operator of g given in (<ref>). When this operator is easy to compute, each step of the algorithm can be carried out efficiently.

The iterations of the alternating minimization method are summarized as follows:

    ^k+1   ∈ P_Z_(_g( Re(^k))), 
    ^k+1   ∈_g( Re(P_Z_(^k))).

We now consider several special cases of (<ref>):
 
  * If g=0, then _g is the identity and we recover the original algorithm from Fienup <cit.>, or alternating projection <cit.>.

6cmFienup 

	Initialization. ^0∈^n.

  General Step. For k ∈,

    ^k+1 = Re(P_Z_(^k)).




The convergence result given in Theorem <ref> also holds in this case since constant functions are convex and continuous.

  * If g() = λ_1 for some λ>0, then _g = 𝒯_λ, where 𝒯_λ is the soft thresholding operator (see footnote on page footnote_soft). We refer to the resulting algorithm as “AM L1".


6cmAM L1 

	Initialization. ^0∈^n, λ>0

  General Step. For k ∈,

    ^k+1 = 𝒯_λ( Re(P_Z_(^k))).



  * If g=δ_C_K, where C_K is the set of all K-sparse vectors, C_K = {∈^n: _0 ≤ K}, then _g = P_C_K is the so-called hard thresholding operator. This operator outputs a vector which is all zeros except for the largest K components (in absolute values) of its input vector which are kept the same. The hard thresholding operator is multivalued and the resulting algorithm, which we term “AM L0" picks an arbitrary point in its range.


6cmAM L0 

	Initialization. ^0∈^n, K ∈

  General Step. For k ∈,

    ^k+1∈ P_C_K( Re(P_Z_(^k))).





 §.§ Majorization-Minimization Interpretation

In this section, we focus on partial minimization in . We show that the value of this partial minimization leads to a least-squares objective. This allows us to interpret the Fienup algorithm as a majorization-minimization process on this least-squares function. For the rest of this section, for any ∈^n, we denote by () an arbitrary but fixed member of P_Z_().


  §.§.§ Partial Minimization in 

The following lemma provides a connection between partial minimization in  and the evaluation of a nonsmooth least-squares objective.

	
	For any ∈^n, we have
	
    min_∈ Z_1/2 - _2^2 = 1/2n|()| - _2^2.



An optimal solution of the minimization problem is given by = ^-1() where  has the form (<ref>). Now,

    min_∈ Z_1/2 - _2^2    = 1/2 - ^-1()_2^2
       =1/2^-1(() -))_2^2
       =1/2n() - _2^2.

Using the expression of  in (<ref>), we have for all j=1, 2, …, n,

    | ()[j] - [j]| =
    	
    		|| ()[j]| - [j]|,		
    		   if  |()[j]| ≠ 0, 
    [j],     otherwise.

Putting everything together,

    min_∈ Z_1/2 - _2^2  =1/2n() - _2^2   = 1/2n|()| - _2^2,

which completes the proof.

A direct consequence of Lemma <ref> is the following corollary that connects between problem (<ref>) and a regularized nonlinear least-squares problem.
 The pair (,) is an optimal solution of problem (<ref>) if and only if  is an optimal solution of

    min{ F() ≡1/2n|()| - _2^2+g() },

and = ^-1(), where  is of the form given in (<ref>).

Note that in (<ref>), the least-squares objective is defined with respect to the amplitude |()| and not the magnitude-squared |()|^2.
For random measurements, it has been shown in <cit.> that the amplitude objective leads to superior performance over the standard magnitude-squared approach.


  §.§.§ Fienup as Majorization-Minimization


In order to understand further the connection with the Fienup algorithm, we define the following auxiliary function:

    h (,) ≡1/2 - ()_2^2 + g().

Now, for any ∈^n, using Lemma <ref>, we have the following properties (recalling the definition of F in (<ref>)):

    h(, )    = 1/2 - ()_2^2 + g() 
       ≥1/2 - ()_2^2 + g() = F(),  ∀∈^n,
    
    	h(, )    = 1/2 - ()_2^2 + g() = F().

In other words, using the convexity of g, we have that h(, ·) is a 1-strongly convex global upper bound on the objective F. Computing this upper bound amounts to performing partial minimization over  in (<ref>). Minimizing the upper bound h(, ) in  corresponds to partial minimization over  in (<ref>). The upper bound is tight in the sense that we recover the value of the objective at the current point, h(,) = F(). Therefore the alternating minimization algorithm is actually a majorization-minimization method for the nonsmooth least-squares problem

    min_∈^n1/2n|()| - _2^2 + g().

The steps presented in (<ref>) can then be summarized as follows:

    ^k+1   = _ h(^k, ) 
       = _g( Re((^k))) 
       =_g( Re(P_Z_(^k))),

which is exactly the mapping given in (<ref>).



 §.§ Projected Gradient Descent Interpretation


We now provide an additional interpretation of the alternating minimization algorithm as a projected gradient method for an optimization problem related to (<ref>) which consists of a smooth convex objective and a nonconvex constraint set. This interpretation is valid whenever g is assumed to be proper lower semicontinuous and convex.

For any ∈^n and ∈^n, we can write (<ref>) as - _2^2 =  - ()_2^2 + ()_2^2. To move from complex numbers to real numbers, we set _1 = () and _2=(). Defining a new constraint set Z̃_ = {(_1, _2) ∈^n ×^n: _1 + i _2 ∈ Z_}, problem (<ref>) can be  equivalently rewritten in the form

    min_∈^n, (_1, _2)∈Z̃_{1/2 - _1_2^2 + 1/2_2_2^2 + g()}.


Minimizing first w.r.t. , (<ref>) reduces to the following minimization problem in _1,_2:

    min_(_1, _2) ∈Z̃_{ H(_1,_2) ≡ G(_1)+1/2_2^2 },

where

    G(_1) ≡min_∈^n{1/2_1 - _2^2 + g()}.









The following result allows us to relate the gradient of H to the optimization primitives used in the alternating minimization method.

	
Assume that g is proper, lower semicontinuous and convex. Then the function H is continuously differentiable, its gradient is 1-Lipschitz and can be expressed as	

    ∇ 	H(_1,_2) =  ( _1 - _g(_1),  _2 ).











	From Moreau <cit.>, we know that G is differentiable and ∇ G() =  - _g () = _g^*(), where g^* is the conjugate function of g, which is convex. The computation of the gradient of H is then immediate. We can use the fact that proximity operators of convex functions are nonexpansive <cit.> to verify that ∇ H is 1-Lipschitz. Indeed, for any (_1, _2) and (_1, _2), we have
	
    ∇ 	H(_1,_2) - ∇ 	H(_1,_2)_2^2 
    
    		=    _g^*(_1) - _g^*(_1)_2^2 + _2 - _2_2^2 
    ≤    _1 - _1_2^2 + _2 - _2_2^2 
    
    		=    (_1,_2) - (_1,_2)_2^2,

completing the proof.





Consider applying projected gradient descent to solve (<ref>). From Lemma <ref>, we can use a step size of magnitude 1. In this case, taking into account the form of the gradient given in (<ref>), we obtain that the general update step takes the form






    (_1^k+1,_2^k+1)
       =    P_Z̃_((_1^k,_2^k) - ∇ H(_1^k,_2^k))
       =   P_Z̃_((_1^k,_2^k) - (_1^k-_g(_1^k),_2^k))
       =     P_Z̃_(_g(_1^k),0).

We now go back to the complex domain by setting = _1 + i _2. Note that projecting (_1, _2) onto Z̃_ is equivalent to projecting  onto Z_. With this notation, the iterations of projected gradient descent can be summarized by the following iteration mapping (on complex numbers):

    ^k+1   = P_Z_(_g((^k))),

which is exactly the same as (<ref>). Therefore, the Fienup algorithm is equivalent to projected gradient descent with unit stepsize applied to the formulation (<ref>). Note that from the point of view of nonsmooth analysis, problem (<ref>) is much better behaved than (<ref>).


§ CONSEQUENCES AND EXTENSIONS


The interpretations of Section <ref> can be used to analyze the convergence of alternating minimization applied to problem (<ref>) and to offer extensions of the method.


 §.§ Convergence Analysis

Our main convergence result is given in the following theorem. Recall that a function is semi-algebraic if its graph can be defined by combining systems of polynomial equalities and inequalities (for example, the ℓ_1 norm is semi-algebraic).

	
	Assume that g is proper, lower semicontinuous, convex and semi-algebraic. Then the sequence {^k,^k}_k ∈ generated by the alternating minimization algorithm satisfies the following:
	
		
  (i) It holds that ∑_k≥ 0^k+1 - ^k_2 < + ∞ and the sequence {^k}_k ∈ converges to a point ^* ∈^n.
		
  (ii) For any accumulation point ^* of {^k}_k ∈, (^*, ^*) is a Fréchet critical point of problem (<ref>) and (_1^*,_2^*)=((^*),(^*)) is a Fréchet critical point of problem (<ref>).
	

The proof is quite technical and is given in Appendix <ref>. The semi-algebraic assumption on g can be relaxed to representability in o-minimal structures over the real field, see <cit.>. Therefore, the proposed result actually applies to much more general regularizers. For example, using boundedness of the feasible set in (<ref>), the same result holds if g is analytic (see the dicussion in <cit.>). The arguments build upon a nonsmooth variant of the celebrated Kurdyka-Łojasiewicz (KL) property <cit.>. Note that direct application of the results of <cit.> to projected gradient descent or the results of <cit.> to majorization minimization is not possible here.

The most important implication of Theorem <ref> is that the sequence of estimated signals converges smoothly to a point which satisfies certain optimality conditions related to problems (<ref>), and (<ref>). This is a departure from standard convergence results that are only able to guarantee that accumulation points of the generated sequence of iterates satisfy certain optimality conditions. It is important to underline that the result is global: it holds for any initialization of the algorithm and does not require any regularity assumption beyond semi-algebraicity and convexity of g. This is in contrast with local convergence results which are typical for alternating projection methods <cit.> that are applicable when the prior term g is an indicator function.


 §.§ Acceleration and Momentum Term

A benefit of the interpretation of alternating minimization as a projected gradient method is that it allows to propose new variants inspired by known extensions for projected gradient algorithms. In this section we focus on the incorporation of an inertial term that results in an alternating minimization scheme that includes a momentum term. This line of research has a long history in optimization, starting with the development of the heavy-ball method <cit.> which inspired an optimal first order scheme for convex optimization developed by Nesterov <cit.>, and its extension to convex composite problems with the FISTA method <cit.>. Although this last technique was proposed and analyzed only in the context of convex optimization, we consider its application in our nonconvex constrained problem since it empirically provides interesting results. The resulting algorithm is referred to as FISTAPH, and is described as follows.

7cmFISTAPH: FISTA for Phase retrieval 

	Initialization. ^0∈ Z_ and α^k∈[0 , 1) for all
		k ∈. Set ^0 = ^0 and ^-1=^0.

  General Step. For k ∈,
	
		
  * ^k+1∈ P_Z_(_g((^k))).
		
  * ^k+1 = ^k + α^k( ^k - ^k-1).
	
	



If ^m is the last produced iteration, then the output of the algorithm is =_g((^m)). A typical choice for the weight sequence is α^k = k-1/k+2.
 The question of the convergence of the iterates produced by this method in nonconvex settings is an interesting topic to explore in future research. We may also further consider monotone variants of similar types of methods, see e.g. <cit.>.

In the numerical experiments we employ FISTAPH in the setting where g() = λ_1 for some λ>0. In this case, _g = 𝒯_λ with 𝒯_λ being the soft thresholding operator with parameter λ (see footnote on page footnote_soft).



§ EXPERIMENTS AND NUMERICAL RESULTS


In this section, we describe experiments and numerical results comparing the different algorithms introduced in Section <ref> on the task of phase retrieval.


 §.§ Experimental Setup

Given measurements  as in (<ref>), our problem consists of finding the corresponding _0. We focus on the setting in which _0 is known to be sparse. We vary  the signal size n (with J = {1, 2, …, n/2}), the sparsity level K and the signal to noise ratio (SNR). In the following discussion, we will refer to a recovery method ℳ which can be seen as a black box which takes as input a vector of measurements ∈_+^n, support information J, sparsity level K, an initial estimate  and outputs an estimate ∈^n with () ⊆ J and _0 ≤ K. One recovery experiment consists of the following:

	
  * Fix a recovery method ℳ, a signal length n, a support information set J = {1, 2, …, n/2}, a sparsity level K and an SNR.
	
  * Generate _0 ∈^n by the following procedure:
		
			
  * Choose K coordinates among J uniformly at random.
			
  * Set these coordinate values at random in [-4, -3] ∪ [3, 4].
			
  * Set all other coordinates to be 0.
		
	
  * Generate the measurements ^2 = |(_0)|^2 +, where  is white Gaussian noise according to the chosen SNR. Set negative entries of ^2 to be 0 in order to take square root.
	
  * Call method ℳ 100 times with data (, J, K) and randomly generate initial estimates to get 100 candidate solutions {_it}_it = 1, 2, …, 100.
	
  * Compute the best estimate _best with best = _it = 1,2…,100{|(_it)| - _2^2}.
	
  * Compare _best and _0 (modulo Fourier invariances) with the following metric (sign is understood coordinatewise with (0) = 0):
		
    recovery(_best, _0)    =
    				
    					1,   (_best) = (_0)
    
    					0,    otherwise.


This procedure was repeated 100 times. That is, for each method, signal length, sparsity level and SNR, we have 100 signal recovery experiments, each one associated with a support recovery status. We aggregate these results by considering the recovery probability (average of recovery(_best, _0)) and the median CPU usage for a single simulation (100 calls to the method with different initialization estimates). We use the same initialization for all methods by careful initialization of random seeds. All the experiments were performed on a desk station with two 3.2 GHz Quad Core Intel Xeon processors and 64GB of RAM.


 §.§ Implementation Details

In our numerical implementation, we used the following stopping criterion.

	
  * For alternating minimization and Wirtinger methods: the difference in successive objective value less than 10^-8.
	
  * For GESPAR: no swap improvement.
	
  * For FISTAPH: the norm of the gradient mapping less than 10^-8.

The tuning of these criteria allows to balance accuracy and computational time to some extent.

The ℓ_1 penalized problem includes a prior sparsity inducing term of the form g(·) = λ·_1. It is necessary to tune the λ parameter in order to obtain meaningful results. We considered the following strategies for different methods.

	
  * For alternating minimization, we set λ = 0.2 in all experiments.
	
  * For Wirtinger based method, λ is tuned a posteriori as a function of n and k. The experiment was conducted for λ = 1, 2.15, 4.64, 10, 21.5, 46.4, 100, 215, 464, and we report only the best experiment for each setting.

An interesting feature of alternating minimization based methods is that in our experiments, recovery performance was very consistent for different values of λ in different settings. As a result, we chose a single value of λ for all experiments. The tuning of λ for Wirtinger based algorithms is practically much more difficult. In particular, we found that the best λ was a highly dependant function of the sparsity level K.

Finally, we note that ℓ_0 based priors have the sparsity level of the estimate, K, as a parameter. On the other hand, ℓ_1 based priors will not necessarily produce K-sparse estimates. We therefore use truncation and keep the K largest entries in absolute value of the last iteration.



 §.§ Numerical Results

The performance in terms of support recovery are presented in Figure <ref> with the corresponding algorithm run time in Figure <ref>. Each point in these plots is an average over 100 simulations of the recovery process, each simulation consisting of 100 random initializations of the method considered. AM corresponds to Fienup methods with different priors, FISTA is the accelerated variant, and WIRT stands for Wirtinger.

We make the following observations from the numerical results:

	
  * For alternating minimization, there is a consistent increase in recovery performance by switching from ℓ_0 to ℓ_1 based regularization priors.
	
  * The ℓ_1 prior degrades the performance of Wirtinger based methods compared to the ℓ_0 prior.
	
  * FISTAPH consistently provides the best performance and is significantly faster than its competitors.
	
  * Fienup with ℓ_0 prior leads to lower performance compared to GESPAR, which was already reported in <cit.>.

As described in the experimental section, we added noise on the squared measurements rather than on the measurements themselves. This noise model is closer to the optimization model considered for GESPAR and Wirtinger flow than model (<ref>) which is related to problem (<ref>). We tried to change the noise model on a subset of experiments (additive noise on the measurements rather than squared measurements), however,  the performance of the different methods was very similar. Therefore, we only report results related to squared-measurement noise model.






§ CONCLUSION

The main theoretical contribution of this work is to provide a strong theoretical basis to the fact that Fienup-type methods, when used with Fourier transforms and convex priors, lead to smoothly converging sequences of estimates. This result holds under minimal assumptions and in particular, it holds globally, independently of the initialization point. Furthermore, we characterize the properties of the limiting point as Fréchet critical points of different optimization problems. These results shed light on important properties of one of the most well known algorithms used in the context of phase retrieval. Furthermore, based on an interpretation as a projected gradient method, we proposed a new variant of Fienup with the incorporation of a momentum term which we call FISTAPH.

On the practical side, we demonstrated based on numerical simulations that FISTAPH with ℓ_1 regularization constitutes a very competitive alternative to other methods in the context of sparse phase retrieval.


§ PROOF OF THEOREM <REF>


The proof involves many notions of nonsmooth analysis which can be found in <cit.>. Throughout the proof, we only consider subgradients of subdifferentially regular functions. Each subgradient can be interpreted as a Fréchet subgradient and the subgradient set valued mapping is closed.
	We adopt the notation of Section <ref>, letting = _1 + i _2 for two real vectors _1 and _2 and consider the constraint set Z̃_ = {(_1, _2) ∈^n ×^n; _1 + i _2 ∈ Z_}. We let K(, _1, _2) = 1/2 - _1_2^2 + 1/2_2_2^2 + g() be the objective function of problem (<ref>) which with this notation becomes
	
    min_∈^n, (_1, _2) ∈Z̃_ K(, _1, _2).

	We will denote by δ_Z̃_, the indicator function of the set Z̃_ (0 on the set and +∞ outside). We set K̃(, _1, _2) = K(, _1, _2) + δ_Z̃_(_1, _2) so that problem (<ref>) is equivalent to the (unconstrained) minimization of K̃.
	


	Proof of (i):
	Using <cit.>, the subgradient of this nonsmooth function is of the form
	
    ∂K̃(, _1, _2)
    		   = (
    		[         ∂_K̃(, _1, _2); ∂_(_1, _2)K̃(, _1, _2) ])
       = (
    		[                           - _1 + ∂ g(); (
    			[ _1 -;   _2 ]) + ∂δ_Z̃_ (_1, _2) ]).

	Partial minimization over iterations yields the following
	
    0    ∈^k+1 - ^k_1 + ∂ g(^k+1)
    
    		0    ∈(
    			[ ^k_1 - ^k;      ^k_2 ]) + ∂δ_Z̃_ (^k_1, ^k_2).

	Combining these, we have
	
    (
    		[                                                        0; (
    			[ ^k_1 - ^k+1;        ^k_2 ]) + ∂δ_Z̃_ (^k_1, ^k_2) ])
    ⊂    ∂K̃(^k+1, ^k_1, ^k_2).

	Using (<ref>),
	
    (
    		[         0; ^k - ^k+1;         0 ])
    		∈∂K̃(^k+1, ^k_1, ^k_2).

	Finally, from strong convexity of K̃ with respect to its first argument, we have
	
    K̃(^k+1, ^k_1, ^k_2) + 1/2^k+1 - ^k^2_2 
    ≤    K̃(^k, ^k_1, ^k_2) 
    ≤    K̃(^k, ^k-1_1, ^k-1_2).

	Since g is semi-algebraic, K̃ is also semi-algebraic. Any semi-algebraic function satisfies the nonsmooth Kurdyka-Łojasievicz property <cit.>. We can now use the now well established recipe <cit.> <cit.> with the two conditions (<ref>) and (<ref>) to obtain that the sequence {^k+1 - ^k_2 }_k ∈ is summable. This proves statement (i) (convergence holds by Cauchy criterion).
	


	Proof of (ii):
	Using the fact that K̃ has compact sublevel sets, the sequence { (^k+1, _1^k, _2^k) }_k∈ is bounded and hence has a converging subsequence. We fix an accumulation point (^*, _1^*, _2^*) of the sequence (note that ^* is given by (i)). We remark that, thanks to (<ref>) and the fact that ^k+1 - ^k→ 0, any accumulation point of the sequence is a critical point of K̃. Furthermore, since ^k →^*, we have using (<ref>) that
	
    -(
    		[ ^*_1 - prox_g(_1^*);                ^*_2 ]) ∈∂δ_Z̃_ (^*_1, ^*_2).

	This is actually the criticality condition for problem (<ref>) which proves statement (ii).







