Sparse Phase Retrieval of One-Dimensional Signals by Prony's Method

In this paper, we show that sparse signals f representable as a linear combination of a finite number N of spikes at arbitrary real locations or as a finite linear combination of B-splines of order m with arbitrary real knots can be almost surely recovered from O(N^2) Fourier intensity measurements up to trivial ambiguities. The constructive proof consists of two steps, where in the first step the Prony method is applied to recover all parameters of the autocorrelation function and in the second step the parameters of f are derived. Moreover, we present an algorithm to evaluate f from its Fourier intensities and illustrate it at different numerical examples.


. Introduction
Phase retrieval problems occur in many scienti c elds, particularly in optics and communications. They have a long history with rich literature regarding uniqueness of solutions and existence of reliable algorithms for signal reconstruction, see e.g. [SEC + ] and references therein. Usually, the challenge in solving one-dimensional phase retrieval problems is to overcome the strong ambiguousness by determining appropriate further information on the solution signal. Previous literature on characterization of ambiguities of the phase retrieval problem with given Fourier intensities is often concerned with the discrete problem, where a signal x in R N or C N has to be recovered. For an overview on the complete characterization of nontrivial ambiguities is this discrete case as well as on appropriate additional signal information we refer to our survey [BP a] and further recent results in [BP , Bei a, Bei b]. Contribution of this paper. In this paper, we consider the continuous one-dimensional sparse phase retrieval problem to reconstruct a complex-valued signal from the modulus of its Fourier transform. Applications of this problem occur in electron microscopy, wave front sensing, laser optics [SST , SSD + ] as well as in X-ray crystallography and speckle imaging [RCLV ]. For the posed problem, we will show that for sparse signals the given Fourier intensities are already su cient for an almost sure unique recovery, and we will give a construction algorithm to recover f .

B G P
We assume that the sparse signal is either of the form ) or, for m > , with c (m) j ∈ C, T j ∈ R for j = , . . . , N , where δ denotes the Delta distribution, and B j,m is the B-spline of order m being determined by the (real) knots T j < T j+ < . . . < T j+m . We want to recover these signals from the Fourier intensities | f (ω)| and will show that only O(N ) samples are needed to recover f , i.e. all coe cients c (m) j , j = , . . . , N and knots T j , j = , . . . , N + m, almost surely up to trivial ambiguities. The proposed procedure is constructive and consists in two steps. In a rst step, we employ Prony's method in oder to recover all parameters of the (squared) Fourier intensity function F[f ](ω) . In a second step, we recover the parameters T j and the complex coe cients c j that determine the desired signal. Related work on sparse phase retrieval. While the general phase retrieval problem has been extensively studied for a long tome, the special case of sparse phase retrieval grew to a strongly emerging eld of research only recently, particularly often connected with ideas from compressed sensing. Most of the papers consider a discrete setting, where the N -dimensional real or complex k-sparse vector x has to be reconstructed from measurements of the more general form | a j , x | with vectors a j forming the rows of a measurement matrix A ∈ C M ×N . The needed number M of measurements depends on the sparsity k.
If A presents rows of a Fourier matrix, this setting is close to the sparse phase retrieval problem considered in optics, see e.g. [JOH ]. Here the problem is rst rewritten as (non-convex) rank minimization problem, then a tight convex relaxation is applied and the optimization problem is solved by a re-weighted l -minimization method. The related approach in [ESM + ] employs the magnitudes of the short-time Fourier transform and applies the occurring redundancy for unique recovery of the desired signal. A corresponding reconstruction algorithm is then based on an adaptation of the GESPAR algorithm in [SBE ].
In [LV ], the measurement matrix A is taken with random rows and the PhaseLift approach [CSV ] leads to a convex optimization problem that recovers the sparse solution with high probability. Employing a thresholded gradient descent algorithm to a non-convex empirical risk minimization problem that is derived from the phase retrieval problem, Cai et al. [CLM ] have established the minimax optimal rates of convergence for noisy sparse phase retrieval under sub-exponential noise.
Other papers rely on the compressed sensing approach to construct special frame vectors a j to ensure uniqueness of the phase retrieval problem with high probability, where the number of needed vectors is O(k ), see e.g. [WX , OE , IVW ]. We would like to emphasize that all approaches employing general or random measurement matrices in phase retrieval are quite di erent in nature from our phase retrieval problem based on Fourier intensity measurements. In this paper, we want to stick on considering Fourier intensity measurements because of their particular relevance in practice.
Early attempts to exploit sparsity of a discrete signal for unique recovery using Fourier intensities go back to unpublished manuscripts by Yagle [Yaga, Yagb], where a variation of Prony's method is applied in a non-iterative algorithm to sparse signal and image reconstruction. Unfortunately, the algorithm proposed there not always determines the signal support correctly.
The continuous one-dimensional phase retrieval problem has been rarely discussed in the literature, see [Wal , Hof , RCLV , Bei b, BP b]. In the preprint [RCLV ], the authors also considered the recovery of sparse continuous signals of the form ( . ). However, in that paper the sparse phase retrieval problem is in turn transferred into a turnpike problem that is computationally expensive to solve. Moreover there exist cases, where a unique solution cannot be found, see [Blo ]. Our method circumvents this problem by proposing an iterative procedure to x the signal support (resp. the knots of the signal represented as a B-spline function) where the corresponding signal coe cients are evaluated simultaneously. Organization of this paper. In Section , we shortly recall the mathematical formulation of the considered sparse phase retrieval problem and the notion of trivial ambiguities of the phase retrieval problem that always occur.
Section is devoted to the special case of phase retrieval for signals of the form ( . ). Using Prony's method, we give a constructive proof for the unique recovery of the Nsparse signal f up to trivial ambiguities using / N (N − ) + Fourier intensity measurements. Here we have to assume that the knot di erences T j − T k are pairwise di erent.
In Section , the ansatz is generalized to the unique recovery of spline functions of the form ( . ) where we need to employ / (N + m)(N + m − ) + Fourier intensity measurements. In Section , we present an explicit algorithm for the considered sparse phase retrieval problem and illustrate it at di erent examples.
. Trivial ambiguities of the phase retrieval problem We wish to recover an unknown complex-valued signal f : R → C of the form ( . ) or ( . ) with compact support from its Fourier intensity | F[f ] | given by Unfortunately, the recovery of the signal f is complicated because of the well-known ambiguousness of the phase retrieval problem. Transferring [BP a, Proposition . ] to our setting, we can recover f only up to the following ambiguities.
Proposition . . Let f be of a signal of the form ( . ) or a non-uniform spline function of the form ( . ). Then (i) the rotated signal e iα f for α ∈ R, (iii) and the conjugated and re ected signal f (−·) have the same Fourier intensity | F[f ] |.
Proof. Applying the properties of the Fourier transform, we have Considering the absolute value of each equation yields the assertion.
Although the ambiguities in Proposition . always occur, they are of minor interest because of their close relation to the original signal. For this reason, we call ambiguities caused by rotation, time shift, conjugation and re ection, or by combinations of these trivial. In the following, we will show that for the considered sparse signals the remaining non-trivial ambiguities only occur in rare cases.
. Phase retrieval for distributions with discrete support Initially, we restrict ourselves to the recovery of signals f of the form ( . ) with complexvalued coe cients c ( ) j and spike locations T < · · · < T N .
and the known squared Fourier intensity | F[f ] | can be represented by Thus, in order to recover f being determined by the coe cients c ( ) j ∈ C and the knots T j ∈ R, j = , . . . , N , we will recover all parameters of the exponential sum in ( . ) in a rst step and then derive the desired parameters of f in a second step.
. . First step: Parameter recovery by Prony's method Assuming that the non-zero knot di erences T j − T k with j k are pairwise di erent, and denoting the distinct frequencies T j − T k in increasing order by In order to recover the parameters τ ℓ and the unknown coe cients γ ℓ from the exponential sum ( . ) we employ Prony's method [Hil , PT ]. Let h > be chosen such that hτ ℓ < π for all ℓ = , . . . , N (N − ) / .
Using the de nition of the Prony polynomial Λ(z) in ( . ), we observe that Consequently, the vector of remaining coe cients λ ≔ (λ , . . . , λ N (N − ) ) T of the Prony polynomial Λ(z) can be determined by solving the linear equation system . Since the Hankel matrix H can be written as This assumption has been ensured by choosing an h such that hτ ℓ ∈ (−π, π), since the τ ℓ had been supposed to be pairwise di erent.
Knowing the coe cients λ k of Λ(z), we can determine the unknown frequencies τ ℓ by evaluating the roots of the Prony polynomial ( . ). The coe cients γ ℓ can now be computed by solving the over-determined equation system with a Vandermonde-type system matrix.
The procedure summarized above is the usual Prony method, adapted to the nonnegative exponential sum P (ω) in ( . ). In the numerical experiments in Section , we will apply the approximate Prony method (APM) in [PT ]. APM is based on the above considerations but it is numerically more stable and exploits the special properties Let us now investigate the question, how many intensity values are at least necessary for the recovery of P (ω) in ( . ). Counting the number of unknowns of P (ω) in ( . ), we only need to recover the / N (N − ) + real values γ and Re γ ℓ , Im γ ℓ , τ ℓ , for ℓ = , . . . N (N − ) / . We will show now that using the special structure of the real polynomial P (ω) in ( . ) and of the Prony polynomial Λ(z) in ( . ), we indeed need only / N (N − ) + exact equidistant real measurements P (kh), k = , . . . , / N (N − ) to recover all parameters determining P (ω). This can be seen as follows.

. . Second step: Unique signal recovery
Having determined the parameters τ ℓ as well as the corresponding coe cients γ ℓ of ( . ), we want to reconstruct the parameters T j and c ( ) j , j = , . . . , N , of f in ( . ) in a second step.
Theorem . . Let f be a signal of the form ( . ), whose knot di erences T j − T k di er pairwise for j, k ∈ { , . . . , N } with j k, and whose coe cients satisfy |c ( ) | |c ( ) N |. Further, let h be a step size such that h(T j − T k ) ∈ (−π, π) for all j, k. Then f can be uniquely recovered from its Fourier intensities | F[f ](hℓ) | with ℓ = , . . . , / N (N − ) up to trivial ambiguities.
Proof. Applying Prony's method to the given data | F[f ](hℓ) |, we can compute the frequencies τ ℓ and the related coe cientsγ ℓ of the squared Fourier intensity ( . ). Again, we assume that the frequencies τ ℓ occur in increasing order and, further, denote the list of positive frequencies by T ≔ {τ ℓ } N (N − ) / ℓ= . Obviously, the maximal distance τN (N − ) / now corresponds to the length T N −T of the unknown f in ( . ). Due to the trivial shift ambiguity, we can assume without loss of generality that T = and T N = τN (N − ) / . Further, the second largest distance τ ( N (N − ) / )− corresponds either to T N − −T or to T N −T . Due to the trivial re ection and conjugation ambiguity, we can assume that T N − = τ ( N (N − ) / )− . By de nition, there exists a τ ℓ * > in our sequence of parameters T such that τ ℓ * + τ ( N (N − ) / )− = τN (N − ) / , and τ ℓ * hence corresponds to the knot di erence T N − T N − . Thus, we obtain These equations lead us to and thus to Since f can only be recovered up to a global rotation, we can assume that c ( ) is real and non-negative, which allows us to determine the coe cients c ( ) , c ( ) N , and c ( ) N − in a unique way.
To determine the remaining coe cients and support knots of f , we notice that the third largest distance τ ( N (N − ) / )− corresponds either to T N −T or to T N − −T . As before, we always nd a frequency τ ℓ * such that τ ( with the related coe cient γ ℓ * = c ( ) c ( ) . Moreover, we have γ ( However, only one of the two equalities in ( . ) and ( . ) can be true, since if both were true then γ ℓ * c ( ) contradicting the assumption that |c N | |c |. Consequently, either the equation in ( . ) or the equation in ( . ) holds true and we can either determine T with c ( ) Removing all parameters τ ℓ from the sequence of distances T that correspond to the di erences T j − T k of the recovered knots, we can repeat this approach to nd the remaining coe cients and knots of f inductively.
If we identify the space of complex-valued signals of the form ( . ) with the real space R N , the condition that two knot di erences T j − T k and T j − T k are equal for xed indices j , j , k , and k de nes a hyper plane with Lebesgue measure zero. An analogous observation follows for the condition |c ( ) | = |c Remark . . . Since the proof of Theorem . is constructive, it can be used to recover an unknown signal ( . ) analytically and numerically. If the number N of spikes is known beforehand then the assumption of Theorem . can be simply checked during the computation. If the assumption regarding pairwise di erent distances T j − T k is not satis ed, then the application of Prony's method in the rst step yields less than B G P N (N − ) + pairwise distinct parameters τ ℓ . The second assumption |c N | |c | can be veri ed in the second step, where c ( ) , c ( ) N − , and c ( ) N are evaluated. . A similar phase retrieval problem had been transferred to a turnpike problem in [RCLV ]. The turnpike problem deals with the recovery of the knots T j from an unlabeled set of distances. Although this problem is solvable under certain conditions, a backtracing algorithm can have exponential complexity in the worst case, see [LSS ].

. Retrieval of spline functions with arbitrary knots
In this section, we generalize our ndings to spline functions of order m ≥ . Let us recall that the B-splines B j,m in ( . ) being generated by the knot sequence T < · · · < T N +m are recursively de ned by see for instance [Boo , p. ]. Further, we notice that for ≤ k ≤ m − the kth derivative of the spline f in ( . ) is given by where the coe cients c (m−k ) j are recursively de ned by with the convention that c (m−k+ ) = c (m−k+ ) N +k = , see [Boo , p. ]. For k = m − , equation ( . ) coincides with a step function, i.e., with the right derivative of the linear spline f (m− ) . Further, in a distributional manner, the mth derivative of f is given by with the coe cients and the Dirac delta distribution δ.
Applying the Fourier transform to ( . ), we now obtain and thus Since the exponential sum on the right-hand side of ( . ) has exactly the same structure as the exponential sum in ( . ), we can immediately generalize Theorem . by considering Theorem . . Let f be a spline function of the form ( . ) of order m, whose knot distances T j − T k di er pairwise for j, k ∈ { , . . . , N + m} with j k, and whose coe cients satisfy |c ( ) | |c ( ) N +m |. Further, let h be a step size such that h(T j −T k ) ∈ (−π, π) for all j, k. Then f can be uniquely recovered from its Fourier intensities | with c ( ) ≔ and c (m−k+ ) ≔ , which nishes the proof.
Corollary . . Almost all spline functions f of order m in ( . ) can be uniquely recovered from their Fourier intensities F[f ] up to trivial ambiguities.

. Numerical experiments
Since the proofs of Theorem . and Theorem . are constructive, they can be straightforwardly transferred to numerical algorithms to recover a spline function from its Fourier intensity. However, the classical Prony method introduced in subsection . is numerically unstable with respect to inexact measurements and to frequencies lying close together. For this reason, there are numerous approaches to improve the classical method. In order to verify Theorem . and Theorem . numerically, we apply the so-called approximate Prony method (APM) proposed by Potts and Tasche in [PT , Algorithm . ] for recovery of parameters of an exponential sum of the form with τ ℓ = −τ −ℓ and γ ℓ = γ −ℓ . The algorithm can be summarized as follows, where the exact number M + of the occurring frequencies in ( . ) needs not be known beforehand.
. Repeat step with respect to the remaining frequencies τ ℓ .
A second adaption of the proof of Theorem . concerns the reconstruction of the coe cients c (m) j from the recovered coe cients c ( ) j . In order to describe the relation between the coe cients as linear equation system, we de ne the rectangular matrices Then, the recursion between the coe cients c (m−k+ ) j and c (m−k ) j can be stated as where we use the coe cient vectors c (m−k ) ≔ (c (m−k ) j ) N +k j= . Instead of computing the coe cients stepwise from left to right, we can determine the coe cients c With these modi cations, we recover a spline function of order m from its Fourier intensity by the following algorithm. . Compute the squared Fourier intensity of the mth derivative of the spline at the given points by . Apply the approximate Prony method (Algorithm . ) to determine the knot distances τ ℓ with ℓ = − (N +m)(N +m− ) / , . . . , (N +m)(N +m− ) / in increasing order and the corresponding coe cients γ ℓ .
. Update the reconstructed distances and coe cients by . Set T ≔ , T N +m ≔ τ(N+m)(N+m− ) / , and T N +m− ≔ τ ( (N +m)(N +m− ) / )− ; nd the index ℓ * with |τ ℓ * − T N +m + T M +m− | ≤ ε; and compute the corresponding coe cients by Initialize the lists of recovered knots and coe cients by and remove the used knot distances from the set T : . For the maximal remaining distance τ k * in T, determine the index ℓ * with |τ k * + τ ℓ * − T M +n | ≤ ε.
Remove all distances between the new knot and the already recovered knots from T and repeat step until the set T is empty.
-. i Table : Knots T j and coe icients c ( ) j of the spike function in Example . -.
-. i Table : Knots T j and coe icients c ( ) j of the spline function in Example .