Phase Retrieval for Partially Coherent Observations
Jonas Kornprobst, Alexander Paulus, Josef Knapp, Thomas F. Eibert
IIEEE TRANSACTIONS ON SIGNAL PROCESSING 1
Phase Retrieval for Partially Coherent Observations
Jonas Kornprobst,
Graduate Student Member, IEEE , Alexander Paulus,
Graduate Student Member, IEEE ,Josef Knapp,
Graduate Student Member, IEEE , and Thomas F. Eibert,
Senior Member, IEEE
Abstract —Phase retrieval is in general a non-convex and non-linear task and the corresponding algorithms struggle with theissue of local minima. We consider the case where the mea-surement samples within typically very small and disconnectedsubsets are coherently linked to each other — which is a rea-sonable assumption for our objective of antenna measurements.Two classes of measurement setups are discussed which canprovide this kind of extra information: multi-probe systemsand holographic measurements with multiple reference signals.We propose several formulations of the corresponding phaseretrieval problem. The simplest of these formulations poses alinear system of equations similar to an eigenvalue problemwhere a unique non-trivial null-space vector needs to be found.Accurate phase reconstruction for partially coherent observationsis, thus, possible by a reliable solution process and with judgmentof the solution quality. Under ideal, noise-free conditions, therequired sampling density is less than two times the numberof unknowns. Noise and other observation errors increase thisvalue slightly. Simulations for Gaussian random matrices and forantenna measurement scenarios demonstrate that reliable phasereconstruction is possible with the presented approach.
Index Terms —linear phase retrieval, non-convex non-linearcost function minimization, phaseless/magnitude-only near-fieldfar-field transformation, equivalent source reconstruction.
I. I
NTRODUCTION P HASE RETRIEVAL is a type of inverse problem arisingin many research fields, including optics [1], [2], X-raycrystallography [3], [4], high-frequency engineering [5]–[8],transmission electron microscopy [9], [10], coherent diffrac-tion imaging [11]–[13] and ptychography [14]–[16]. From amathematical point of view, the behavior of phase retrievalalgorithms is often studied for abstract scenarios, e.g., Gaus-sian random matrices, to develop new solution strategies, sincethis allows to (statistically) analyze the convergence behaviorof the algorithms [17]–[24].For some other than random Gaussian matrices — e.g.,Fourier-base magnitude-only measurements —, phase retrievalcan not be proven to work with absolute certainty, or absolutecertainty can only be attained for unrealistic measurementaccuracies and with quadratic oversampling rates [8]. Hence,a rather popular strategy is to integrate additional informationinto the problem formulation, going beyond magnitude-onlymeasurements. If feasible, one may change the observationkernel (masking [25], exploitation of multiple measurementdistances [5], [6], [26]) or enforce sparsity [27]–[31]. Anotherway is to relax the original assumption of magnitude-onlyobservations to some extent.
The authors are with the Chair of High-Frequency Engineering, Departmentof Electrical and Computer Engineering, Technical University of Munich,80290 Munich, Germany (e-mail: [email protected], [email protected]).This work was supported by the German Federal Ministry for EconomicAffairs and Energy under Grant 50RK1923.
In this work, we consider the case where (possibly small)subsets of the total observations are captured coherently. Inantenna measurements, this can be achieved with specialmulti-probe or multi-frequency measurements [7], [32]–[38].In optics, measurements with specialized masks containingstructured modulations can be seen as multi-probe measure-ments [39]. A similar problem arises in holography for thestitching of holographic images, where the relation betweenreference signal (from a reference antenna or beam) and mea-surement signal is changed. Examples of holographic antennameasurement techniques are found in [40]–[47] and of opticalholography in [48]–[50]. Angular synchronization [51]–[53]or phase unwrapping [54], [55] are similar but different.The phase differences obtained from such measurementsare typically disconnected and many phase differences neededfor straightforward phase reconstruction by concatenation aremissing. Thus, the determination of the absolute phases atall the measurement locations remains a non-trivial phaseretrieval problem. In order to solve the phase-retrieval problemby the concatenation of phase differences, the observationsneed to be collected in an appropriate manner [33]. Due topractical restrictions related to positioning accuracy and errorconcatenation, such an approach is only of limited utility.We do not restrict our investigations to the special caseof near-field (NF) antenna measurements and propose severalgeneral formulations of the magnitude-only inverse problemwith partially coherent observations together with relatedsolution strategies. With sufficiently oversampled observationsand suitable forward operator properties, we show that thephase retrieval problem even can be formulated as a linearsystem of equations which contains the observation data ina non-linear manner. The formulations are applicable to anykind of phase retrieval where partially coherent observationsare available. The great advantage is that the search for theglobal minimum — even under the influence of measurementerrors — becomes feasible. Phase retrieval results for syntheticdata, noisy and ideal, demonstrate that a solution as close aspossible to the true one can be reconstructed with certainty,once the corresponding sampling limit is reached — given thatmild conditions on the forward operator of the inverse problemare fulfilled.The paper is structured as follows. Section II introducesthe magnitude-only inverse problem. Several formulations forphase retrieval with partially coherent observations are pro-posed in Section III. Section IV demonstrates the applicabilityto random complex matrices and illustrates the limitationsof all variants. Furthermore, we investigate a possible wayof how to incorporate the method into NF far-field (FF)transformations (NFFFTs) with special multi-element probes. a r X i v : . [ ee ss . SP ] J a n IEEE TRANSACTIONS ON SIGNAL PROCESSING [ x ] [ x ] [ x ] [ x ] [ b ] [ b ] [ b ] [ b ] [ b ] [ b ] ⇒ [| b |] [| b |] [| b |] [| b |] [| b |] [| b |] [| b |] [| b |] [| b |] [| b |] e j Δ 𝜙 , [| b |] e j Δ 𝜙 , [| b |] e j Δ 𝜙 , ⇒ | · | | · | | · | A ∈ C 𝑀 × 𝑁 A ∈ C 𝑀 × 𝑁 A ∈ C 𝐶𝑀 × 𝑁 x ∈ C 𝑁 b ∈ C 𝐶𝑀 | b | ∈ R 𝐶𝑀 B ∈ C 𝐶𝑀 × 𝑀 | b |∈ R 𝑀 | b |∈ R 𝑀 diag B ∈ R 𝑀 diag B ∈ C 𝑀 B ∈ R 𝑀 × 𝑀 B ∈ C 𝑀 × 𝑀 Δ 𝜙 , Δ 𝜙 , Δ 𝜙 , ⇒ Figure 1. An illustration of the ordering of partially coherent observations in the observation vector and observation matrices for the case 𝑁 = , 𝑀 = and 𝐶 = . From left to right: The complex matrix-vector product A x = b , the magnitude-only observations | b | which are split into the sub-vectors | b | and | b | , the unnamed “partially-coherent” observation vector with entry-wise coherence between the first and second halves of the observations, which areemployed as the diagonal entries of the diagonal matrices B and B , and, eventually, the “partially-coherent” observation matrix B = [ B B ] T . II. P
HASE R ETRIEVAL — P
ROBLEM S TATEMENT
A. The Complex Inverse Problem
In order to introduce the phase retrieval problem statement,we start with a discrete inverse problem
A x = b written as aleast-squares optimization problem min x (cid:13)(cid:13) A x − b (cid:13)(cid:13) (1)for retrieving the unknown column vector x ∈ C 𝑁 from theobservation column vector b ∈ C 𝑀 . The relation betweenobservations and unknowns is established by the forward-operator matrix A ∈ C 𝑀 × 𝑁 with rk A ≤ min ( 𝑁, 𝑀 ) . We con-sider rk A as the number of degrees of freedom (DoFs) of theinverse problem. We restrict our theoretical investigations tothe case of a uniquely defined solution x , with 𝑀 ≥ rk A = 𝑁 . B. The Classical Magnitude-Only Inverse Problem
The phase-retrieval problem min x (cid:13)(cid:13) | A x | − | b | (cid:13)(cid:13) (2)enforces magnitude-only equality between the reconstruction Ax and the observations | b | ∈ R 𝑀 , where | · | is the element-wise absolute-value operator. An alternative formulation of thesame problem reads min 𝝓 , x (cid:13)(cid:13) A x − diag ( 𝝓 )| b | (cid:13)(cid:13) s . t . |[ 𝝓 ] 𝑚 | = 𝑚 ∈ N [ ,𝑀 ] , (3)where N [ ,𝑀 ] represents the natural numbers { , , . . . , 𝑀 } and diag (·) creates a diagonal matrix from a vector (or a vectorfrom the diagonal of a matrix). An additional unknown vector,the phase vector 𝝓 ∈ C 𝑀 with the 𝑚 th entry [ 𝝓 ] 𝑚 = e j 𝜙 𝑚 (4)has been introduced in (3). For reasonably small problems, a rank-revealing decomposition can beemployed to reduce the singular case (with
𝑁 , 𝑀 > rk A and 𝑁 ≠ 𝑀 ) tothe case discussed here. This may not be feasible for large 𝑁 or 𝑀 . III. T HE M AGNITUDE -O NLY P ROBLEM WITH P ARTIALLY C OHERENT O BSERVATIONS
A. Basic Assumptions
Let us assume that two specific observations – the 𝑚 th and 𝑘 th ones – are measured coherently. Then, the phase differencebetween these two observations is known. According to (4),we are able to introduce the additional constraint [ 𝝓 ] 𝑘 /[ 𝝓 ] 𝑚 = e j ( 𝜙 𝑘 − 𝜙 𝑚 ) = e j Δ 𝜙 𝑘,𝑚 (5)including the observed quantity Δ 𝜙 𝑘,𝑚 . The observed phasedifference Δ 𝜙 𝑘,𝑚 is only one out of many, where most mayremain unknown. Hence, the problem has to be reformulatedwith all these remaining unknown phase differences in mind.In order to study the effect of such partially coherentobservations, we constrain the way of how these observationsare taken. A special observation probe shall be able to capture 𝐶 independent observations coherently whenever it performsa measurement. For an illustrating example, see Fig. 1, where 𝑁 = unkowns and 𝑀 = sets of observations takenby a special probe with 𝐶 = elements are considered.The three phase differences e j Δ 𝜙 , , e j Δ 𝜙 , , and e j Δ 𝜙 , withinthree distinct pairs of measurement samples are observed inaddition to the magnitude-only measurement samples. Withthree of these observation pairs — i.e., 𝑀 = —, we observe 𝐶 𝑀 = magnitude samples but only ( 𝐶 − ) 𝑀 = phasedifferences. 𝑀 = phase unknowns remain to be found; thisis a compromise between retrieving just a single global phase(fully coherent measurements) and retrieving all 𝐶 𝑀 = phases (not coherent at all). Phase retrieval for partially coherent observations is asimpler task than the general phase retrieval problem sinceadditional information is available. Nevertheless, the algo-rithms found in literature, which tackle this particular problem, This restriction is not necessary in order to benefit from the proposedphase-retrieval method, but it helps to simplify the notation and to predict atwhich oversampling ratio reliable phase retrieval is possible. In some kindsof more general measurements, 𝐶 may change rather arbitrarily from oneobservation to another and the derived formulations do also hold in that case.However, the notation and, more importantly, the derivation founded on thisnotation, is simplified considerably by the regular structure of the observations. ORPROBST et al. : PHASE RETRIEVAL FOR PARTIALLY COHERENT OBSERVATIONS 3 unkownobject [ | b |] 𝑘 [ | b |] 𝑚 Δ 𝜙 𝑘,𝑚 (a) unkownobjectreference 2 𝐶 observationsreference 1 𝐶 observations (b)Figure 2. Electromagnetic measurement setups for partially coherent obser-vations. (a) A multi-probe approach, 𝐶 = . (b) A holographic approach withtwo coherent data sets, 𝑀 = and 𝐶 is large. are limited to solving non-convex non-linear minimizationproblems [7], [36] or they utilize restrictive — unrealistic oreven unfeasible — sampling strategies [33], [36], [56].In the general case, the forward matrix is A ∈ C 𝐶𝑀 × 𝑁 ,the magnitude-only observation vector is | b | ∈ R 𝐶𝑀 , and thephase unknowns vector is 𝝓 ∈ C 𝐶𝑀 . The inverse problem min 𝝓 , x (cid:13)(cid:13) A x − diag ( 𝝓 )| b | (cid:13)(cid:13) s . t . |[ 𝝓 ] 𝑚 | = 𝑚 ∈ N [ ,𝐶𝑀 ] [ 𝝓 ] 𝑚 + 𝑐𝑀 [ 𝝓 ] 𝑚 = e j Δ 𝜙 𝑚 + 𝑐𝑀,𝑚 for 𝑚 ∈ N [ ,𝑀 ] and for 𝑐 ∈ N [ ,𝐶 − ] (6)is now additionally constrained by the observed phase differ-ences Δ 𝜙 𝑚 + 𝑐𝑀,𝑚 . Alternatively, we can state that the phasedifferences Δ 𝜙 𝑘,𝑚 according to (5) are only observed if mod ( 𝑘 − 𝑚, 𝑀 ) = . These phase differences do not carrythe same information as the standard complex data, sincethey are fewer in number and concatenation is not possible —for instance, the phase differences from the first subset of 𝐶 observations 𝑚 = { , 𝑀 + , 𝑀 + , . . . , ( 𝐶 − ) 𝑀 + } to any other observation outside this subset are missing.Let us emphasize the key aspect that there are now 𝐶 𝑀 magnitude-only observations and ( 𝐶 − ) 𝑀 coherent andlinearly independent phase-difference observations. Simpleconcatenation of the phase differences is not feasible as 𝑀 − phase differences plus a global phase are still unknown. Aug-menting the phase retrieval problem (3) in (6), we expect theminimum number of observations for successful reconstructionin the range 𝑁 ≤ 𝐶 𝑀 ≤ 𝑁 , with the approximate empiricalupper bound 𝑁 for standard phase retrieval [19]. B. Possible Phase-Difference Measurement Techniques
Let us consider two basic types of measurement scenarios:i) multi-probe systems, see Fig. 2(a), and ii) holographicsystems with multiple reference signals, see Fig. 2(b), wherereciprocal permutations are possible (the reference may be onthe transmitting or receiving side). Fig. 2 focuses on antennameasurement systems; comparable systems in optics have beenmentioned in the introduction.The first case of multi-probe observations seems morechallenging since 𝐶 is typically rather small, e.g., 𝐶 = in [7], [32]–[35]. In such a scenario, very small (localized) and pos-sibly unconnected “isles” of coherent observations are takeneverywhere on the measurement surface of interest, which isfeasible in two ways. The particular phase differences of inter-est are observed either directly by multi-channel receivers withshared oscillator signals or via distinct magnitude observationsin the form of [| b |] 𝑘 , [| b |] 𝑚 , |[ b ] 𝑘 + [ b ] 𝑚 | , |[ b ] 𝑘 + j [ b ] 𝑚 | ,which allows to reconstruct the phase differences as Δ 𝜙 𝑘,𝑚 = atan |[ b ] 𝑘 + [ b ] 𝑚 | − [| b |] 𝑘 − [| b |] 𝑚 |[ b ] 𝑘 + j [ b ] 𝑚 | − [| b |] 𝑘 − [| b |] 𝑚 . (7)Other sets of at least four linear combinations of [ b ] 𝑘 and [ b ] 𝑚 may yield the same reconstructed phase differences.Essential for the remainder of the paper is only that we can as-sume the phase differences of the measurement samples withinthe small “isles” to be known. Note that, however, the influenceof noise differs for different hardware implementations.In the second case, the relevant task is to stitch holographicsub-images which are, in themselves, fully coherent — thenumber of coherent observations 𝐶 is here typically ratherlarge. One can think of myriads of variations of such mea-surement setups, where the constant behavior of the referencesignal source has to be known a priori or may even remainunknown but constant. One possible antenna measurementscenario, see Fig. 2(b), is to record two coherent sub-imagesfor two installation locations of a receiving reference antenna. C. Phase Retrieval Formulations with Coherence Constraints1) A non-linear minimization problem:
The structuring ofthe 𝐶 similar blocks is not yet visible in (6). Hence, wewill explain the various newly introduced variables as shownin Fig. 1. The observation vector | b 𝑐 | ∈ R 𝑀 represents the 𝑐 th block of the magnitude-only observations, composing thecomplete observation vector as, see Fig. 1, | b | = (cid:2) | b T1 | . . . | b T 𝑐 | . . . | b T 𝐶 | (cid:3) T (8)with the corresponding forward-operators A 𝑐 ∈ C 𝑀 × 𝑁 . Here,the transpose is denoted by (·) T . Furthermore, we employa reduced phase unknowns vector 𝝍 ∈ C 𝑀 for the phaseunknowns of | b | only. The phase differences to the 𝑐 th blockare implemented together with the observed magnitudes as adiagonal matrix with entries in the 𝑚 th row and column [ B ] 𝑚𝑚 = [| b |] 𝑚 , B ∈ R 𝑀 × 𝑀 , (9) [ B 𝑐 ] 𝑚𝑚 = [| b 𝑐 |] 𝑚 e j Δ 𝜙 𝑚 + 𝑐𝑀,𝑚 , B 𝑐 ∈ C 𝑀 × 𝑀 ,𝑐 ∈ N [ ,𝐶 − ] , (10)leading to the overall block-structured matrix, see Fig. 1, B = (cid:2) B B . . . B 𝑐 . . . B 𝐶 (cid:3) T , B ∈ C 𝐶𝑀 × 𝑀 (11)including all observed phase differences and magnitudes. The 𝑀 subgroups of observations are assembled in B , whereeach column contains 𝐶 entries including the observed phasedifferences. Furthermore, the 𝑚 th column has an unknownphase [ 𝝍 ] 𝑚 for its complex entries. This choice is arbitrary and, due to the flexibility of the phase vector, hasno influence on the solution or on the solution process at all.
IEEE TRANSACTIONS ON SIGNAL PROCESSING
These auxiliary quantities allow to rewrite (6) as min 𝝍 , x (cid:107) A x − B 𝝍 (cid:107) s . t . |[ 𝝍 ] 𝑚 | = 𝑚 ∈ N [ ,𝑀 ] . (12)This formulation is already much easier to implementthan (6), since only one non-convex side constraint withreduced dimension remains. One simple trick to get rid ofthis side constraint is to replace the phase unknowns by thereconstructed observations in the manner of min x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) | A x | − | b | A ... A 𝑐 ... A 𝐶 x − B ... B 𝑐 ... B 𝐶 B − A x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) . (13)Employing 𝝍 = B − A x instead of |[ 𝝍 ] 𝑚 | = as in (12) is analternative restriction. Its implications remain to be studied.
2) A linear formulation for the original unknowns:
Thenon-linear side constraints (dependent on the type of solver,including a yet-to-determine weighting) and the additional 𝝍 unknowns might be obstructive to deal with. Hence, we ana-lyze under which circumstances a unique 𝝍 exists — renderingthe magnitude-one side constraint obsolete.We proceed by looking at the 𝐶 sub-equations of (12) A 𝑐 x = B 𝑐 𝝍 . (14)The matrices B 𝑐 are diagonal with non-zero entries and, thus,invertible. Solving for 𝝍 yields 𝝍 = B − A x = . . . = B − 𝑐 A 𝑐 x = . . . = B − 𝐶 A 𝐶 x , (15)which resembles a concatenation of generalized (pseudo-)eigenvalue problems [57]–[62] (with a known eigenvalue ofvalue ) and enables us to eliminate the phase unknownsfrom the problem by simultaneously solving 𝐶 − generalizedeigenvalue problems (with 𝜆 = ) in the form of B − A x = 𝜆 B − 𝑐 A 𝑐 x for 𝑐 ∈ N [ ,𝐶 ] . (16)Multiplying by the diagonal matrices B B 𝑐 leads to B 𝑐 A x = B A 𝑐 x for 𝑐 ∈ N [ ,𝐶 ] , (17)and allows to recast the overall problem as Q x = B A − B A ... B 𝑐 A − B A 𝑐 ... B 𝐶 A − B A 𝐶 x = (18)with the matrix Q ∈ C 𝑀 ( 𝐶 − )× 𝑁 .Equation (18) is just the linear and homogeneous partof (13) multiplied by B in order to avoid the divisionby potentially small B entries and without the magnitudeconstraint | A x | = | b | . Having dropped the non-linear partof the constraints implies that not all information is employedin the source reconstruction. This may lead to a somewhatsuboptimal solution, but with the benefit of solving a linearsystem of equations — implying an improvement in reliability. ,
000 1 ,
200 1 , − − −
202 10 . . singular value index 𝑖 s i ngu l a r v a l u e s l og 𝜎 𝑖 SVD spectrum of (18)SVD spectrum of (21)
Figure 3. SVD spectra for a random matrix and random right-hand side, with 𝑁 = , 𝑀 = , 𝐶 = , and 𝑛 = − . The properties of (18) are influenced by the properties of A 𝑐 .We have assumed that rk A 𝑐 = min { 𝑀, 𝑁 } . Hence, with aknown x , 𝝍 follows immediately. However, it is still unclearunder which circumstances x is unique.Due to the block subtractions B 𝑐 A − B A 𝑐 in the matrix Q ∈ C 𝑀 ( 𝐶 − )× 𝑁 , the true x has to be in ker Q . We deduce acondition for a unique reconstruction: There has to be a non-trivial ker Q with dim ker Q = . To fulfill this, we recall that Q has 𝑁 columns. Hence, for dim ker Q = , it is requiredthat rk Q = 𝑁 − . In order to yield this rank, a necessary butinsufficient condition on the number of rows of Q reads 𝑀 ( 𝐶 − ) ≥ 𝑁 − . (19)A further requirement is that the eigenvalue 𝜆 = in the gen-eral eigenvalue problems (16) is unique, i.e., it is not degener-ate. For the subtraction of two random Gaussian matrices, theprobability of having a degenerate eigenvalue approaches zero.For other matrices, this becomes more difficult to achieve, seethe application example in Section V. If noisy observations arecontained in the matrices B 𝑐 , degenerate or near-degenerateeigenvectors may prevent a reliable reconstruction. This isfurther discussed when measurement errors are considered.The task is now to determine the unique non-trivial vectorin ker A 𝑐 . We consider the example of a Gaussian random A and b with 𝑁 = , 𝑀 = , 𝐶 = , and a noise-to-signalratio in the observation vector 𝑛 = − according to 𝑛 = (cid:107) b (cid:48) − b (cid:107) (cid:107) b (cid:107) , (20)where the primed vector b (cid:48) contains the noise-contaminatedobservations. The spectrum of the singular value decomposi-tion (SVD) of Q = USV H is shown in Fig. 3. Exactly onenoise-limited singular value with a magnitude of about − is observed. In a magnitude-ordered SVD, we refer to thisvale as 𝜎 𝑁 . The corresponding singular vector v 𝑁 solves thephase-retrieval problem to a comparable accuracy level.The task of finding the non-trivial vector in the null spacecan be tackled in different ways. The obvious one is to performan SVD and pick the vector for the smallest singular value, butthis is computationally rather expensive. From a complexitypoint of view, we can employ an iterative approach (preferablya Krylov-subspace method, e.g., Arnoldi iterations [63]) toestimate the required SVD vector and scale the retrieved x appropriately afterwards. Such an iteratively attained solutionis also unique once the discussed conditions on Q are fulfilled. ORPROBST et al. : PHASE RETRIEVAL FOR PARTIALLY COHERENT OBSERVATIONS 5
Furthermore, it does not suffer from local minima, i.e., itis independent from the initial guess. This is explained bythe linearity (and, thus, convexity) of the formulation. Sincethe nullity of Q is one, and (18) is a homogeneous linearsystem of equations, even standard solvers for linear systemsof equations may be employed if the trivial solution is avoided.The formulation offers two possibilities to judge whetherthe reconstruction was successful. Firstly, a drop in the SVDspectrum between the smallest and the second smallest sin-gular value should be observable. Secondly, the reconstructedphase vector 𝝍 according to (15) is required to have entrieswith constant magnitude. Unit-magnitude entries of the phasevector are to be created by a suitable scaling. If the magnitudesof the vector entries fluctuate a lot, the reconstruction failed.
3) A linear formulation for the phase unknowns:
In the stepfrom (12) to (13), the phase unknowns were replaced by x .However, it is also possible the other way round. If 𝝍 is uniqueaccording to (19), x may be replaced by the Moore-Penrosepseudo-inverse A + applied to the “complex” observations B 𝝍 yielding again a linear null-space equation R 𝝍 = (cid:2) AA + B − B (cid:3) 𝝍 = (21)with the matrix R ∈ C 𝐶𝑀 × 𝑀 .In Fig. 3, the spectrum of R = USV H is shown for arandom matrix with 𝑁 = and 𝑛 = − . Again, onesingular vector v 𝑀 of a non-trivial noise-limited null-space(with a singular value 𝜎 𝑀 ) appears. While the null-spaceseems to have a numerically decreased dynamical range, thereconstruction quality is the same in this particular example.Whether one of the two versions is superior to the other isstudied in the results sections. The main difference to (18) isfound in the fact that the forward operator appears only in theform of the projector AA + , removing most of the influence ofthe spectrum or null-space of A .After the reconstruction of 𝝍 by (21), the intuitive approachfor the solution of the phaseless problem is to solve A x = B 𝝍 (22)in a subsequent step. This is done in the following unless statedotherwise. Note, however, that the magnitude-one constraintin the reconstruction of 𝝍 was neglected in order to obtaina linear system of equations. Since not all information isconsidered in the reconstruction process, the retrieved solutionmay be suboptimal. Instead of B 𝝍 , the complex measurementvector may be reconstructed as A x = B diag (| 𝝍 |) − 𝝍 (23)enforcing the non-linear magnitude-one constraint. D. Discussion of the Linearized Reconstruction Algorithms1) Suboptimality of the Retrieved Solution:
An importantquestion is what happens if the conditions on Q (or, equiv-alently, R ) are not met. This might happen if the number ofobservations is not sufficiently large or if (near-) degenerate The initial guess is an important concept for non-convex minimizationproblems. Due to the possibility of local minima on the path to the globalsolution, it matters at which initially guessed vector the minimization starts. − − − − − − − −
202 log ( (cid:107) δ B (cid:107) /(cid:107) B (cid:107) ) l og ( (cid:107) · (cid:107) / (cid:107) 𝝍 (cid:107) ) solution (cid:107) 𝝍 (cid:48) (cid:107) /(cid:107) 𝝍 (cid:107) bound for (cid:107) δ 𝝍 (cid:107) /(cid:107) 𝝍 (cid:107) error (cid:107) δ 𝝍 (cid:107) /(cid:107) 𝝍 (cid:107) (a) − − − − ( (cid:107) δ B (cid:107) /(cid:107) B (cid:107) ) (b)Figure 4. Simulation results for the first-order bound (26), 𝐶 𝑀 = , 𝑁 = , and random realizations of A , b , and δ b . (a) 𝐶 = . (b) 𝐶 = . eigenvalues appear due to observation errors. Then, the nu-merically determined nullity dim ker Q is greater than one. Westill know that the true solution x ∈ ker Q , but (18) or (21)alone are not sufficient anymore. This offers two strategies.Either the minimization problem is constraint by choosingonly search vectors in ker Q , or the null-space equation isaugmented by additional constraints. For instance, the phasevector constraints |[ 𝝍 ] 𝑚 | = may be included again. Anotherway to get rid of false solutions in ker Q for the homogeneousequation (18) is to fix the phase and the magnitude of up to 𝐶 coherent (and hence complex) observations. As part of thefollowing error analysis, a single 𝑖 th entry in 𝝍 is fixed in orderto obtain the inhomogeneous and invertible linear system R ★ 𝝍 = u 𝐶𝑀 + with R ★ = (cid:2) R T u 𝑖 (cid:3) T (24)from (21), where u 𝑖 refers to the 𝑖 th unit vector.
2) Influence of Measurement Errors:
We assume that anoise-contaminated observation vector b (cid:48) = b + δ b , see (20),leads to perturbations δ B and δ R via a non-linear relation.From the perturbation theory of linear systems [64] appliedto (24), we state the first-order bound (cid:107) δ 𝝍 (cid:107)(cid:107) 𝝍 (cid:107) ≈ (cid:13)(cid:13) R + ★ δ R ★ 𝝍 (cid:13)(cid:13) (cid:107) 𝝍 (cid:107) ≤ (cid:13)(cid:13) R + ★ (cid:13)(cid:13) (cid:107) δ R ★ (cid:107) . (25)With (cid:107) δ R ★ (cid:107) ≤ (cid:107) δ B (cid:107) , we have (cid:107) δ 𝝍 (cid:107)(cid:107) 𝝍 (cid:107) (cid:46) 𝜅 (cid:107) δ B (cid:107)(cid:107) B (cid:107) , (26)which holds for any submultiplicative and unitarily invariantmatrix norm (cid:107)·(cid:107) . Since factoring out B from R + ★ is prohibitedby the pseudo-inversion, the dependence on the problem sen-sitivy 𝜅 = (cid:13)(cid:13) R + ★ (cid:13)(cid:13) (cid:107) B (cid:107) remains. Relating the bound to (cid:107) δ b (cid:107)/(cid:107) b (cid:107) is not easy to achieve, since (cid:107) δ b (cid:107) (cid:46) (cid:107) δ B (cid:107) .As observed in Fig. 4, noise-to-signal ratios above 𝑛 ≈ lead to a failure of the algorithm since the one-dimensionalnullspace of R disappears due to noise. It is not possibleanymore to find a suitable phase vector and, hence, a zerovector from the null space of R is retrieved. The other wayround this means that a zero singular value, which is unique inthe noiseless case, does not become degenerate in the presenceof small enough noise. The error bound proves (24) to bestable since limited measurement noise has limited impact on IEEE TRANSACTIONS ON SIGNAL PROCESSING the solution vector. Note that the error is generally lower forlarge values of 𝐶 than for smaller values of 𝐶 . Due to thesimilarities in the approaches, we expect a similar upper boundto hold true for (18).
3) Observations with Zero Magnitude:
At first sight, onemight be afraid that measurement samples with zero mag-nitude cause problems. Indeed, it is impossible to define aphase difference between two complex numbers with one ofthem being zero. However, a (magnitude-only) observationwith zero magnitude constitutes a complex measurement withzero real and imaginary parts. Since the entries of B as definedin (10) contain plain zeros for cases where the phase differencecan neither be observed nor defined, no problem arises forthe proposed algorithms. For a complete tuple of 𝐶 zeroobservations, we can even rewrite the condition in (19) as 𝐶 𝑀 − rk B ≥ 𝑁 − , (27)where rk B and the required number of measurements decreasefrom 𝑀 by one for each complete tuple of zero-observations.
4) On the Computational Complexity:
The computationalcomplexity of the proposed linear formulations is dominatedby the inversion of the linear forward operator A . The mostchallenging case is (21). If a direct solver is utilized, it isemployed twice: for the calculation of the pseudo-inverse A + and for the inversion of R . If an iterative solver is utilized —possibly even with an efficient evaluation of A —, it may alsobe employed in a nested manner for each evaluation of R . Inany case, only the prefactor of the computational complexityof a direct or an iterative solver is increased, where this mayvary dependent on the formulation. E. State of the Art: Algorithms for Comparison
For the standard phase-retrieval problem (2), we considerthe algorithms provided by PhasePack [2], [17], [19], [65]–[69]. The only generally applicable method for incorporatingpartial knowledge about phase differences, which may beemployed to compare it with the proposed algorithm, is foundin [7]. There, the structure of the forward operator is changed,leading to the non-linear non-convex minimization problem min x (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A A A + A A + j A x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − | b || b || b + b || b + j b | (cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) , (28)here given for the case 𝐶 = . The phase-difference constraintis included in the cost functional — and not written as a sideconstraint. Of course, (28) can be rewritten for larger 𝐶 .IV. R ESULTS FOR G AUSSIAN R ANDOM M ATRICES
The linear systems of equations are solved as describedin Section III. Non-linear minimization problems are solvedwith the cost function minimizers provided by Matlab [70],where the active-set method is employed for the minimizationwith equality side constraints. Custom implementations forvarious solvers capable of handling problems of larger sizehave also been realized, based on the memory limited 𝐿 -BFGS method [71], [72], by Broyden, Fletcher, Goldfarb, and − − local minima / outlierssuccess thresholdoversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b (1)(2)(12), 𝐶 = (a) − − success thresholdoversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b (1)(2)(13), 𝐶 = (b)Figure 5. RD of the inverse problem solved for Gaussian random matrices, 𝑁 = , 𝐶 = , 𝑛 = − . Each 𝐶 𝑀 / 𝑁 -ratio with simulations.(a) Magnitude cost-functional minimization with reduced phase unknowns.(b) Magnitude cost-functional minimization with eliminated phase unknowns. Shannon. The initial guess x for the non-linear solvers (6),(12), (13), and (28) is obtained by a spectral method as [19] x = v max | b T | A v max (cid:107)
A v max (cid:107) , (29)with the eigenvector v max corresponding to the largest eigen-value of A H B H BA . The non-linear solvers rely on the initialguess to avoid local minima. The solvers based on the null-space search do not require an initial guess. A. An Extensive Solver Comparison for 𝑁 = Considering 𝑁 = , the phase reconstruction is performedfor 200 random picks of A , each with a randomly picked truesolution 𝝃 and a corresponding right-hand side b = A 𝝃 . Afterthe solution process, the true reconstruction deviation (RD) 𝜖 b = (cid:107) Ax − A 𝝃 (cid:107) /(cid:107) A 𝝃 (cid:107) (30)is evaluated, where, however, the solution x is obtained for anoise-contaminated vector b (cid:48) with the noise-to-signal ratio 𝑛 .For 𝑀 ∈ N [ , ] and 𝑛 = − , the results are shown inFig. 5(a) for the standard solver with phase (1), the standardmagnitude-only solver (2), and the proposed solver with a 𝐶 = phase-differences side constraint according to (12). Thescatter plot provides the insight that the fully-coherent complexsolver always works and the magnitude-only versions requirea certain oversampling for a reliable reconstruction. We furtherobserve that the solver with 𝐶 = converges with fewer 𝐶 𝑀 .In Fig. 5(b), the fusion of magnitude-minimization and null-space condition (13) shows a better convergence than (12),with the main advantage of avoiding to get stuck in localminima for this scenario. In Fig. 6(a), the linear formu-
ORPROBST et al. : PHASE RETRIEVAL FOR PARTIALLY COHERENT OBSERVATIONS 7 − − success thresholdoversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b (1)(2)SVD (18), 𝐶 = (a) − − success thresholdoversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b (1)(2)SVD (21), 𝐶 = (b) − − success thresholdoversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b (1)(2)SVD (21), 𝐶 = (c)Figure 6. RD of the inverse problem solved for Gaussian random matrices, 𝑁 = , 𝐶 = , 𝑛 = − . Each 𝐶 𝑀 / 𝑁 -ratio with simulations. (a) Firstnull-space search with SVD. (b) Second null-space search with SVD andfully linear reconstruction (22). (c) Second null-space search with SVD andenforced magnitude-one constraint (23). lation, i.e., the null-space vector search, is included. Theconvergence behavior is slightly different. Almost exactly at 𝐶 𝑀 / 𝑁 = / ≈ . as expected, we observe a certainconvergence. In contrast to Fig. 5(a), there are no outliers(i.e., local minima) above 𝐶 𝑀 / 𝑁 > . Unfortunately, thelimit for successful reconstruction is a bit higher than in theother two cases. The second null-space equation (21) shows acomparable behavior in Fig. 6(b) with slightly lower RDs. InFig. 6(c), the magnitude-one constraint for the reconstructedphase unknowns is enforced according to (23). This seems toimprove the RD if the reconstruction is working, i.e., if thesuccess threshold 𝐶 𝑀 / 𝑁 ≥ / is fulfilled.We introduce a threshold of 𝑛 with the noise-to-signal ratio 𝑛 = − as defined in (20) and call a reconstruction with a RDbelow this limit successful and above failed . As seen in Figs. 5and 6, this is a rather demanding definition of a successfulreconstruction which excludes three kinds of solutions: globalfalse solutions due to insufficient sampling, wrong solutionsdue to local minima, and almost acceptable solutions, where,e.g., the solver convergence was too slow.This allows us to introduce a success rate for the recon- % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (2)WirtFlowAmplitudeFlowGerchbergSaxtonFienupPhaseMaxPhaseLiftPhaseLamp Figure 7. Success rate for Gaussian matrices, 𝑁 = , noise 𝑛 = − , andfor the comparison 𝐶 𝑀 / 𝑁 = 𝑀 / 𝑁 , PhasePack solvers [67]. struction. In Fig. 7, the minimization of (2) with the 𝐿 -BFGSmethod is compared to many methods provided by PhaseP-ack [67] and the simple cost function minimization is amongthe best-performing solvers. It is a reasonable choice that (28)(the comparison method for partially coherent phase retrieval)is minimized with the L-BFGS method in the following.In Fig. 8(a), the required oversampling ratio 𝐶 𝑀 / 𝑁 for ahigh chance of success for the phase-difference solvers (6)moves towards a value of 𝐶 𝑀 / 𝑁 = with increasing 𝐶 . Thesame is observed for (12) in Fig. 8(b). The standard phaselesssolver (2) performs worst since it has the smallest knowledgeabout the inverse problem. The solver (13), which solves onlyfor x , performs better than the two previous versions.So far, the magnitude-only solvers, including the versionswith partially coherent observations, have shown a rather goodconvergence rate — with a slight advantage of (12) over (6)and a great advantage for (13). However, for increasing 𝑀 ,local minima possibly prevent a reconstruction rate. InFig. 8(d), the SVD is employed to identify the vector for thesmallest singular value in (18). As expected, the transitionfrom failed to successful reconstruction occurs at around 𝐶 𝑀𝑁 ≥ 𝐶 ( 𝑁 − ) 𝑁 ( 𝐶 − ) , (31)which is a recasted version of (19). The SVD-based solver (21)performs marginally better in the transition to a certain recon-struction in Fig. 8(e).The minimizations according to (6) and (12) never reach acertain reconstruction since they strongly depend on the initialguess, and sometimes the spectral method fails in this respect.At an oversampling ratio where the minimizers come closeto a success rate of , the SVD null-space solutionsare able to reach a certain reconstruction. The comparisonmethod (28) is non-convex without guaranteed convergenceand may get stuck in local minima, but empirically it per-forms here quite well. The downsides of this method becomeapparent later. Also, the minimization of (13) achieves certainreconstruction in the hitherto presented results.Finally, we investigate the noise-free case, i.e., 𝑛 = , forthe second null-space solver (21) in Fig. 9. The achievableaccuracy is below − once the necessary oversampling cri-teria are met, e.g., at 𝑀 / 𝑁 = / for 𝐶 = . All proposedphase retrieval algorithms gain reconstruction accuracy oncethe ideal noise-free case is considered. The accuracy of the IEEE TRANSACTIONS ON SIGNAL PROCESSING % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (1) (2)(6), 𝐶 = (6), 𝐶 = (6), 𝐶 = (6), 𝐶 = (6), 𝐶 = (28), 𝐶 = (a) % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (1) (2)(12), 𝐶 = (12), 𝐶 = (12), 𝐶 = (12), 𝐶 = (12), 𝐶 = (28), 𝐶 = (b) % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (1) (2)(13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (28), 𝐶 = (c) % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (1) (2)SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = (28), 𝐶 = (d) % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (1) (2)SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = (28), 𝐶 = (e)Figure 8. Success rate for Gaussian matrices, 𝑁 = , noise level 𝑛 = − ,and for the comparison 𝐶 𝑀 / 𝑁 . (a) Different non-convex solvers (6) in-cluding phase unknowns. (b) Different non-convex solvers (12) including areduced number of phase unknowns. (c) Different non-convex solvers (13)with eliminated phase unknowns. (d) Linearized solvers (18) with SVD.(e) Linearized solvers (21) with SVD. − − − − − − − oversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b SVD (21), 𝐶 = SVD (21), 𝐶 = SVD (21), 𝐶 = SVD (21), 𝐶 = SVD (21), 𝐶 = Figure 9. RD of the inverse problem solved for Gaussian random matrices, 𝑁 = , 𝑛 = . Each 𝐶 𝑀 / 𝑁 ratio with simulations. % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (2)(13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (28), 𝐶 = (a) % % % oversampling ratio 𝐶 𝑀 / 𝑁 r ec on s t r u c ti on s u cce ss r a t e (2)SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = SVD, 𝐶 = (28), 𝐶 = (b)Figure 10. Success rate for Gaussian matrices, 𝑁 = and 𝑛 = − , over 𝐶 𝑀 / 𝑁 . (a) Non-linear cost function solver (13). (b) SVD-solver (18). cost-function minimizations of course depends on the stoppingcriteria for the iterative solver. B. A Larger Scenario with 𝑁 = Now, the noise level is set to 𝑛 = − and the number ofunknowns is increased to 𝑁 = . Among the non-linearminimization techniques, we investigate only the hitherto bestone, which is (13). The two SVD-based solvers were on parso far and we investigate just one of them. The success rates,again for a threshold of 𝑛 , are depicted for two solvers inFig. 10. Two differences to the 𝑁 = case are observed.On the one hand, the SVD-based solver requires a slightlylarger oversampling ratio for convergence than expected fromthe threshold according to (19). E.g., for 𝐶 = , the thresholdfor a successful reconstruction is at 𝐶 𝑀 / 𝑁 = . , butsuccess is only observed at 𝐶 𝑀 / 𝑁 > . The reason is theinteraction of a noise-induced transition period, compare Fig. 6to Fig. 9, and the demanding success threshold of 𝑛 .On the other hand, the cost-function minimization (13)shows a worse success rate than for the case 𝑁 = . Inparticular, it suffers from local minima for this example; thesame happens for the noise-free case and large 𝑁 . If our goal ORPROBST et al. : PHASE RETRIEVAL FOR PARTIALLY COHERENT OBSERVATIONS 9 − − − − signal-to-noise ratio − log 𝑛 R D l og 𝜖 b (2) (28), 𝐶 = (13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (13), 𝐶 = (1) (a) − − − − (23)(22) signal-to-noise ratio − log 𝑛 R D l og 𝜖 b (2) (28), 𝐶 = (21), 𝐶 = (21), 𝐶 = (21), 𝐶 = (21), 𝐶 = (21), 𝐶 = (1) (b)Figure 11. Signal-to-noise ratio analysis, 100 random simulations, 𝑁 = , 𝐶 𝑀 = . 𝑁 , 𝑛 ∈ { − , − , − } . (a) Non-linear solver (13). (b) Linearsolver (21) with phase reconstruction according to (22) and (23). is a certain reconstruction rate of , only the linear solverin Fig. 10 beats the magnitude-only solver and the state-of-the-art solver (28). It becomes clear that the non-linear non-convex minimizations fail to ensure a correct reconstruction.The SVD-based solvers are able to provide absolute certaintyat the theoretically determined thresholds dependent on 𝐶 ,within a noise-caused margin over the expected threshold.This influence of noise on the behavior of the non-linearsolver (13) and the linear solver (21) is analyzed in Fig. 11for 𝑁 = , 𝐶 𝑀 = . 𝑁 , random simulations, and 𝑛 ∈ { − , − , − } . Comparison methods comprise thestandard complex (best-performing) and phaseless (worst-performing) solvers and the multi-probe comparison algo-rithm for 𝐶 = (often suffering from local minima). Theperformance of the partially-coherent phaseless solvers for 𝐶 ∈ N [ , ] depends linearly on the chosen signal-to-noiseratio, except for local minima. The linear method (21) withenforced magnitude-one constraints (23) reliably provides agood solution.V. A N E XEMPLARY A PPLICATION : P
HASE R ETRIEVAL FOR S YNTHETIC A NTENNA N EAR -F IELD M EASUREMENT D ATA
We now consider a synthetic antenna near-field measure-ment setup as illustrated in Fig. 12. For more informationon the idea of NFFFTs, for instance refer to [73]–[80].As part of a phaseless NFFFT, the phases of the observedNFs are to be reconstructed. In the context of the setup inFig. 12, the measurement vector | b | corresponds to the signalsreceived by known probe antennas placed at sample locations(blue diamonds) on a closed hull surrounding the antennaunder test (AUT), which here is a horn antenna. In order to − −
20 0 20 40 − − 𝑥 in cm 𝑦 i n c m Figure 12. Synthetic antenna near-field measurement setup. An L-shapedprobe array (red arrows) around the measurement reference position (greendot) is used to acquire magnitude and local phase difference information (here: 𝐶 = ) at the sample locations (blue diamonds). Probe locations (red arrows)at different sampling locations do not coincide; phase concatenation is notpossible. The solution vector x corresponds to surface currents densities placedtangentially on the smallest sphere (orange sphere) enclosing the horn antenna. model the electromagnetic radiation of the AUT, the unknowncoefficients x of equivalent sources on an enclosing surface(orange sphere) are introduced. The received probe signalsand the coefficients of the equivalent sources are linked via theelectromagnetic radiation operator A , which typically does notfeature Gaussian distributed rows. Skipping the details — werefer the interested reader to [79]–[82] —, the spectrum of A isstrongly decaying and typically exhibits a non-trivial kernel,just to name two major differences to Gaussian matrices.Assuming a probe antenna array and coherence betweenthe probe elements, we are able to apply the presented phaseretrieval algorithms. To pick reasonable multi-antenna probes,we recapitulate that this measurement setup and its fielddistributions are three-dimensional, but a two-dimensionaldescription on the measurement surface is sufficient. Linkingphases in two dimensions is possible in the case of a three-element probe array with an L-shape: Two linearly indepen-dent phase differences are acquired at every measurementlocation, resulting in 𝐶 = . The L-probe (red arrows) placedat an exemplary measurement location (green dot) is illustratedin Fig. 12. For a comparison case with 𝐶 = , we pick thetwo diagonal probe elements only. Real-world measurementsetups for such a scenario are found in [7], [32]–[36]. Westick to synthetic data since i) only simulated results offer theknowledge of the true solution and ii) studies on a large set ofdifferent antennas are not really feasible with measurements.In the considered synthetic measurement setup, theequivalent-source sphere enclosing the AUT and the measure-ment sphere exhibit diameters of and wavelengths, respec-tively. As equivalent sources, 𝑁 = tangential Hertziandipoles are utilized and the horizontal as well as the verticalspacing between the probe-array elements is one wavelength.Each probe element is modeled as a single Hertzian dipole.The obtainable RDs for the two known formulations (2)and (28), as well as the two proposed ones (21) and (18), isdepicted in Fig. 13. The described cases of 𝐶 = and 𝐶 = . . . − − − 𝐶 = oversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b (2)(28)(18)(21) (a) . . . − − − 𝐶 = oversampling ratio 𝐶 𝑀 / 𝑁 R D l og 𝜖 b (2)(28)(18)(21) (b)Figure 13. RD of the inverse problem solved for a synthetic antenna near-fieldmeasurement setup. The average RDs are indicated by solid lines. (a) Utilizinga probe with two diagonal array elements ( 𝐶 = ). (b) Employing a three-element, L-shaped probe array ( 𝐶 = ). are given in Fig. 13(a) and (b), respectively. For every ratio of 𝐶 𝑀 / 𝑁 , random orientations of the AUT were simulated,resulting in different measurement vectors. All results wereobtained for a noise-to-signal ratio of 𝑛 = − .In both cases, the best results are obtained with formula-tion (21), which is observed to reliably yield accurate resultsabove a certain ratio of 𝐶 𝑀 / 𝑁 . The existing formulations (2)and (28) are observed to either fail completely or they arenot guaranteed to find a satisfactory solution. All formulationsexploiting the phase differences are observed to yield better re-sults for the case of the full L-shaped probe ( 𝐶 = ) comparedto the two-element diagonal probe ( 𝐶 = ). Especially formu-lation (28) is observed to benefit significantly. It yields similarresults as (18) in the case of 𝐶 = . The difference betweenthe two proposed null-space formulations, (18) and (21), canbe explained by looking at the spectrum of the singular valuesof Q and R . Considering a noise-free setup with 𝐶 𝑀 / 𝑁 = and the L-shaped probe, the spectra of Q and R are given inFig. 14. The null-space is more distinct for (21), i.e., the ratioof the second smallest singular value to the smallest one issignificantly greater for R than for Q .Whenever perturbations affect the observations, and thusthe spectrum of the singular values, it is more difficult toavoid false solutions and to maintain the desired null-spaceof Q . Since R features a more pronounced separation betweenits non-trivial null-space vector (smallest singular value) andfalse solutions (any other singular value), this formulation ismore robust with respect to noise for the considered exemplaryscenario and its forward operator.VI. C ONCLUSION
Various formulations of the phase retrieval problem withadditional knowledge of phase differences within subsets of ,
000 1 ,
500 2 ,
000 2 , − − −
20 10 . . . . singular value index 𝑖 s i ngu l a r v a l u e s l og 𝜎 𝑖 (18), 𝐶 = (18), 𝐶 = (21), 𝐶 = (21), 𝐶 = (a) ,
000 1 ,
500 2 ,
000 2 , − − −
20 10 . . . . singular value index 𝑖 s i ngu l a r v a l u e s l og 𝜎 𝑖 (18), 𝐶 = (18), 𝐶 = (21), 𝐶 = (21), 𝐶 = (b)Figure 14. Normalized SVD spectra for Q of (18) and R of (21) for thecases 𝐶 = and 𝐶 = of the synthetic measurement data. (a) Noise-freecase 𝑛 = . (b) Noise-contaminated observations with 𝑛 = − . observations have been presented. Two of the formulationsare based on homogeneous linear systems of equations andrequire only defined oversampling and mild conditions to workwith certainty. Since the non-linear magnitude constraints areonly implicitly fulfilled in these linear equations, they do notexploit the full available information and may be seen as suboptimal yet simply reachable solutions. In order to improvethe solution obtained in particular for cases with insufficientoversampling, a minimization problem for a non-linear costfunctional has also been presented, which provides betterresults than comparable algorithms found in literature.The presented results demonstrate that reliable phase re-trieval for partially coherent observations is possible if certainconditions on the phase retrieval equations are fulfilled. In thecase of a well-conditioned random Gaussian forward opera-tor, the performance of the two approaches is very similar.However, significant performance advantages are found forthe phase reconstruction with projection-based method (21),in particular in the presence of noise. Our explanation isthat the projection matrix (21) removes the influence of theconditioning of the forward operator (i.e., its decaying singularvalue spectrum). R EFERENCES[1] Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, andM. Segev, “Phase retrieval with application to optical imaging: Acontemporary overview,”
IEEE Signal Process. Mag. , vol. 32, no. 3,pp. 87–109, May 2015.[2] R. W. Gerchberg and W. O. Saxton, “A practical algorithm for thedetermination of phase from image and diffraction plane pictures,”
Optik ,vol. 35, no. 2, pp. 237–246, 1972.[3] J. Miao, R. L. Sandberg, and C. Song, “Coherent X-ray diffractionimaging,”
IEEE J. Sel. Topics Quantum Electron. , vol. 18, no. 1, pp.399–410, 2012.
ORPROBST et al. : PHASE RETRIEVAL FOR PARTIALLY COHERENT OBSERVATIONS 11 [4] F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval anddifferential phase-contrast imaging with low-brilliance X-ray sources,”
Nat. Phys. , vol. 2, no. 4, pp. 258–261, Apr. 2006.[5] T. Isernia, G. Leone, and R. Pierri, “Radiation pattern evaluation fromnear-field intensities on planes,”
IEEE Trans. Antennas Propag. , vol. 44,no. 5, pp. 701–710, May 1996.[6] R. G. Yaccarino and Y. Rahmat-Samii, “Phaseless bi-polar planar near-field measurements and diagnostics of array antennas,”
IEEE Trans.Antennas Propag. , vol. 47, no. 3, pp. 574–583, Mar. 1999.[7] A. Paulus, J. Knapp, and T. F. Eibert, “Phaseless near-field far-fieldtransformation utilizing combinations of probe signals,”
IEEE Trans.Antennas Propag. , vol. 65, no. 10, pp. 5492–5502, Oct. 2017.[8] J. Knapp, A. Paulus, and T. F. Eibert, “Reconstruction of squaredfield magnitudes and relative phases from magnitude-only near-fieldmeasurements,”
IEEE Trans. Antennas Propag. , vol. 67, no. 5, pp. 3397–3409, May 2019.[9] W. Coene, G. Janssen, M. O. de Beeck, and D. van Dyck, “Phaseretrieval through focus variation for ultra-resolution in field-emissiontransmission electron microscopy,”
Phys. Rev. Lett. , vol. 69, no. 26, pp.3743–3746, Dec. 1992.[10] H. M. L. Faulkner and J. M. Rodenburg, “Movable aperture lenslesstransmission microscopy: A novel phase retrieval algorithm,”
Phys. Rev.Lett. , vol. 93, no. 2, p. 023903, Jul. 2004.[11] M. Guizar-Sicairos and J. R. Fienup, “Phase retrieval with transversetranslation diversity: A nonlinear optimization approach,”
Opt. Express ,vol. 16, no. 10, pp. 7264–7278, May 2008.[12] E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval from codeddiffraction patterns,”
Appl. Comput. Harmonic Anal. , vol. 39, no. 2, pp.277–299, Sep. 2015.[13] J. Bacca, S. Pinilla, and H. Arguello, “Super-resolution phase retrievalfrom designed coded diffraction patterns,”
IEEE Trans. Image Process. ,vol. 29, pp. 2598–2609, 2020.[14] M. A. Iwen, A. Viswanathan, and Y. Wang, “Fast phase retrieval fromlocal correlation measurements,”
SIAM J. Imaging Sci. , vol. 9, no. 4,pp. 1655–1688, 2016.[15] T. Ramos, B. E. Grønager, M. S. Andersen, and J. W. Andreasen, “Directthree-dimensional tomographic reconstruction and phase retrieval of far-field coherent diffraction patterns,”
Phys. Rev. A , vol. 99, no. 2, p.023801, Feb. 2019.[16] N. Sissouno, F. Boßmann, F. Filbir, M. Iwen, M. Kahnt, R. Saab,C. Schroer, and W. zu Castell, “A direct solver for the phase retrievalproblem in ptychographic imaging,”
Math. Comput. Simulation , vol. 176,pp. 292–300, Oct. 2019.[17] E. J. Candes, T. Strohmer, and V. Voroninski, “Phaselift: Exact andstable signal recovery from magnitude measurements via convex pro-gramming,”
Commun. Pure Appl. Math. , vol. 66, no. 8, pp. 1241–1274,2013.[18] E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski, “Phaseretrieval via matrix completion,”
SIAM Rev. , vol. 57, no. 2, pp. 225–251, 2015.[19] E. J. Candes, X. Li, and M. Soltanolkotabi, “Phase retrieval via Wirtingerflow: Theory and algorithms,”
IEEE Trans. Inf. Theory , vol. 61, no. 4,pp. 1985–2007, Apr. 2015.[20] P. Netrapalli, P. Jain, and S. Sanghavi, “Phase retrieval using alternatingminimization,”
IEEE Trans. Signal Process. , vol. 63, no. 18, pp. 4814–4826, Sep. 2015.[21] P. Grohs and M. Rathmair, “Stable gabor phase retrieval for multivariatefunctions,” arXiv preprint arXiv:1903.01104 , 2019.[22] M. A. Iwen, S. Merhi, and M. Perlmutter, “Lower lipschitz bounds forphase retrieval from locally supported measurements,”
Appl. Comput.Harmonic Anal. , vol. 47, no. 2, pp. 526–538, 2019.[23] C. Cheng, I. Daubechies, N. Dym, and J. Lu, “Stable phase retrievalfrom locally stable and conditionally connected measurements,” arXivpreprint arXiv:2006.11709 , 2020.[24] P. Grohs, S. Koppensteiner, and M. Rathmair, “Phase retrieval: Unique-ness and stability,”
SIAM Rev. , vol. 62, no. 2, pp. 301–350, 2020.[25] V. Pohl, F. Yang, and H. Boche, “Phaseless signal recovery in infinitedimensional spaces using structured modulations,”
J. Fourier Anal.Appl. , vol. 20, no. 6, pp. 1212–1233, 2014.[26] C. H. Schmidt and Y. Rahmat-Samii, “Phaseless spherical near-fieldantenna measurements: Concept, algorithm and simulation,” in
Proc.IEEE Antennas Propag. Soc. Int. Symp. (APSURSI) , Charleston, SC,USA, Jun. 2009.[27] K. Jaganathan, S. Oymak, and B. Hassibi, “Sparse phase retrieval:Uniqueness guarantees and recovery algorithms,”
IEEE Trans. SignalProcess. , vol. 65, no. 9, pp. 2402–2410, May 2017. [28] E. J. R. Pauwels, A. Beck, Y. C. Eldar, and S. Sabach, “On Fienupmethods for sparse phase retrieval,”
IEEE Trans. Signal Process. , vol. 66,no. 4, pp. 982–991, Feb. 2018.[29] T. Qiu and D. P. Palomar, “Undersampled sparse phase retrievalvia majorization–minimization,”
IEEE Trans. Signal Process. , vol. 65,no. 22, pp. 5957–5969, Nov. 2017.[30] G. Wang, L. Zhang, G. B. Giannakis, M. Akcakaya, and J. Chen,“Sparse phase retrieval via truncated amplitude flow,”
IEEE Trans.Signal Process. , vol. 66, no. 2, pp. 479–491, Jan. 2018.[31] G. Baechler, M. Krekovic, J. Ranieri, A. Chebira, Y. M. Lu, andM. Vetterli, “Super resolution phase retrieval for sparse signals,”
IEEETrans. Signal Process. , vol. 67, no. 18, pp. 4839–4854, Sep. 2019.[32] S. Costanzo, G. Di Massa, and M. D. Migliore, “Integrated microstripprobe for phaseless near-field measurements on plane-polar geometry,”
Electron. Lett. , vol. 37, no. 16, pp. 1018–1020, Aug. 2001.[33] S. Costanzo and G. Di Massa, “An integrated probe for phaseless plane-polar near-field measurements,”
Microw. Opt. Technol. Lett. , vol. 30,no. 5, pp. 293–295, 2001.[34] S. Costanzo, G. Di Massa, and M. D. Migliore, “A novel hybridapproach for far-field characterization from near-field amplitude-onlymeasurements on arbitrary scanning surfaces,”
IEEE Trans. AntennasPropag. , vol. 53, no. 6, pp. 1866–1874, Jun. 2005.[35] S. Costanzo and G. Di Massa, “Wideband phase retrieval technique fromamplitude-only near-field data,”
Radioengineering , vol. 17, no. 4, pp. 8–12, 2008.[36] R. Tena-Sánchez, M. Sierra-Castañer, and L. J. Foged, “Relative phasereconstruction based on multiprobe solutions and post-processing tech-niques,” in
Proc. Eur. Conf. Antennas Propag. , Copenhagen, DK, Mar.2020.[37] A. Paulus, J. Knapp, J. Kornprobst, and T. F. Eibert, “Improved-reliability phase-retrieval with broadband antenna measurements,” in
Proc. Eur. Conf. Antennas Propag. , Copenhagen, DK, Mar. 2020.[38] J. Knapp, A. Paulus, J. Kornprobst, U. Siart, and T. F. Eibert, “Multi-frequency phase retrieval for antenna measurements,”
IEEE Trans.Antennas Propag. , 2020 (early access).[39] V. Pohl, F. Yang, and H. Boche, “Phase retrieval from low-rate samples,”
Sampling Theory Signal Image Process. , vol. 14, pp. 71–99, 2015.[40] G. Junkin, T. Huang, and J. C. Bennett, “Holographic testing of terahertzantennas,”
IEEE Trans. Antennas Propag. , vol. 48, no. 3, pp. 409–417,Mar. 2000.[41] G. Castaldi and I. M. Pinto, “Well-posed well-conditioned phase retrievaltechnique using a known reference source,” in
Proc. IEEE AntennasPropag. Soc. Int. Symp. USNC/URSI Nat. Radio Sci. Meeting , Salt LakeCity, UT, USA, Jul. 2000, pp. 1780–1782.[42] J. Laviada and F. Las-Heras, “Phaseless antenna measurement on non-redundant sample points via Leith-Upatnieks holography,”
IEEE Trans.Antennas Propag. , vol. 61, no. 8, pp. 4036–4044, Aug. 2013.[43] J. Laviada Martinez, A. Arboleya, Y. Alvarez-Lopez, C. Garcia-Gonzalez, and F. Las-Heras, “Phaseless antenna diagnostics based onoff-axis holography with synthetic reference wave,”
IEEE AntennasWireless Propag. Lett. , vol. 13, pp. 43–46, 2014.[44] A. Arboleya, J. Laviada, J. Ala-Laurinaho, Y. Alvarez, F. Las-Heras,and A. V. Räisänen, “Phaseless characterization of broadband antennas,”
IEEE Trans. Antennas Propag. , vol. 64, no. 2, pp. 484–495, Feb. 2015.[45] A. Arboleya, J. Laviada, and F. Las-Heras, “Scalar calibration forbroadband phaseless antenna measurements based on indirect off-axisholography,”
IEEE Trans. Antennas Propag. , no. 6, pp. 3241–3246, Jun.2018.[46] R. Tena-Sánchez and M. Sierra-Castafier, “Evaluation of software de-fined radio receiver for phaseless near-field measurements,” in
Proc. 40thAnn. Meeting Symp. Antenna Meas. Techn. Assoc. (AMTA) , Williams-burg, VA, USA, Nov. 2018.[47] P. Berlt, C. Bornkessel, and M. A. Hein, “Accurate 3D phase recovery ofautomotive antennas through LTE power measurements on a cylindricalsurface,” in
Proc. Eur. Conf. Antennas Propag. , Copenhagen, DK, Apr.2020.[48] D. Gabor, “Microscopy by reconstructed wave-fronts,” in
Proc. R. Soc.Lond. A , vol. 197, 1949, pp. 454–487.[49] E. N. Leith and J. Upatnieks, “Reconstructed wavefronts and commu-nication theory,”
J. Opt. Soc. America , vol. 52, no. 10, pp. 1123–1130,Oct. 1962.[50] D. A. Barmherzig, J. Sun, P.-N. Li, T. J. Lane, and E. J. Candès,“Holographic phase retrieval and reference design,”
Inverse Problems ,vol. 35, no. 9, p. 094001, 2019.[51] A. Viswanathan and M. Iwen, “Fast angular synchronization for phaseretrieval via incomplete information,” in
Wavelets Sparsity XVI , vol.9597. International Society for Optics and Photonics, 2015, p. 959718. [52] T. Gao and Z. Zhao, “Multi-frequency phase synchronization,” arXivpreprint arXiv:1901.08235 , 2019.[53] M. A. Iwen, B. Preskitt, R. Saab, and A. Viswanathan, “Phase retrievalfrom local measurements: Improved robustness via eigenvector-basedangular synchronization,”
Appl. Comput. Harmonic Anal. , vol. 48, no. 1,pp. 415–444, 2020.[54] R. M. Goldstein, H. A. Zebker, and C. L. Werner, “Satellite radar in-terferometry: Two-dimensional phase unwrapping,”
Radio Sci. , vol. 23,no. 4, pp. 713–720, 1988.[55] D. C. Ghiglia and M. D. Pritt, “Two-dimensional phase unwrapping:Theory, algorithms, and software,”
Hoboken, NJ, US: John Wiley &Sons , Apr. 1998.[56] V. Pohl, H. Boche, and C. Yapar, “Modulated wideband converter forcompressive phase retrieval of bandlimited signals,” in
Proc. Int. ITGConf. Syst. Commun. Coding , Hamburg, Germany, Feb. 2017.[57] K.-C. Toh and L. N. Trefethen, “Calculation of pseudospectra by thearnoldi iteration,”
SIAM J. Sci. Comput. , vol. 17, no. 1, pp. 1–15, 1996.[58] L. N. Trefethen, “Pseudospectra of linear operators,”
SIAM Rev. , vol. 39,no. 3, pp. 383–406, 1997.[59] ——, “Pseudospectra of matrices,”
Numer. Anal. , vol. 91, pp. 234–266,1991.[60] G. Boutry, M. Elad, G. H. Golub, and P. Milanfar, “The generalizedeigenvalue problem for nonsquare pencils using a minimal perturbationapproach,”
SIAM J. Matrix Anal. Appl. , vol. 27, no. 2, pp. 582–601,2005.[61] D. Chu and G. H. Golub, “On a generalized eigenvalue problem fornonsquare pencils,”
SIAM J. Matrix Anal. Appl. , vol. 28, no. 3, pp.770–787, 2006.[62] D. Kressner, E. Mengi, I. Naki´c, and N. Truhar, “Generalized eigenvalueproblems with specified eigenvalues,”
IMA J. Numer. Anal. , vol. 34,no. 2, pp. 480–501, Apr. 2014.[63] W. E. Arnoldi, “The principle of minimized iterations in the solutionof the matrix eigenvalue problem,”
Q. Appl. Math. , vol. 9, no. 1, pp.17–29, Apr. 1951.[64] J. W. Demmel,
Applied Numerical Linear Algebra . Society forIndustrial and Applied Mathematics, 1997.[65] J. R. Fienup, “Phase retrieval algorithms: A comparison,”
Appl. Opt. ,vol. 21, no. 15, pp. 2758–2769, 1982.[66] H. Zhang and Y. Liang, “Reshaped Wirtinger flow for solving quadraticsystem of equations,” in
Proc. Adv. Neural Inf. Process. Syst. 29 ,Barcelona, Spain, Dec. 2016, pp. 2622–2630.[67] R. Chandra, Z. Zhong, J. Hontz, V. McCulloch, C. Studer, and T. Gold-stein, “PhasePack: A phase retrieval library,” in
Proc. 51st AsilomarConf. Signals Syst. Comput. , Pacific Grove, CA, USA, Oct. 2017, pp.1617–1621.[68] O. Dhifallah, C. Thrampoulidis, and Y. M. Lu, “Phase retrieval via poly-tope optimization: Geometry, phase transitions, and new algorithms,” arXiv preprint arXiv:1805.09555 , 2018.[69] T. Goldstein and C. Studer, “Phasemax: Convex phase retrieval via basispursuit,”
IEEE Trans. Inf. Theory , vol. 64, no. 4, pp. 2675–2689, Apr.2018.[70] MATLAB, version 9.7.0 (R2019b) Update 2
Math. Program. , vol. 45, no. 1, pp. 503–528,Aug. 1989.[72] J. Nocedal and S. J. Wright,
Numerical Optimization , 2nd ed. NewYork, NY, USA: Springer Science+Business Media, 2006.[73] A. Ludwig, “Near-field far-field transformations using spherical-waveexpansions,”
IEEE Trans. Antennas Propag. , vol. 19, no. 2, pp. 214–220, Feb. 1971.[74] J. E. Hansen,
Spherical Near-Field Antenna Measurements . London,UK: Peter Pelegrinus Ltd., 1988.[75] T. K. Sarkar and A. Taaghol, “Near-field to near/far-field transformationfor arbitrary near-field geometry utilizing an equivalent electric currentand MoM,”
IEEE Trans. Antennas Propag. , vol. 47, no. 3, pp. 566–573,Mar. 1999.[76] Y. Álvarez, F. Las-Heras, and M. R. Pino, “Reconstruction of equivalentcurrents distribution over arbitrary three-dimensional surfaces based onintegral equation algorithms,”
IEEE Trans. Antennas Propag. , vol. 55,no. 12, pp. 3460–3468, Dec. 2007.[77] C. H. Schmidt, M. M. Leibfritz, and T. F. Eibert, “Fully probe-correctednear-field far-field transformation employing plane wave expansion anddiagonal translation operators,”
IEEE Trans. Antennas Propag. , vol. 56,no. 3, pp. 737–746, Mar. 2008. [78] M. A. Qureshi, C. H. Schmidt, and T. F. Eibert, “Efficient near-fieldfar-field transformation for nonredundant sampling representation onarbitrary surfaces in near-field antenna measurements,”
IEEE Trans.Antennas Propag. , vol. 61, no. 4, pp. 2025–2033, Apr. 2013.[79] J. L. Araque Quijano and G. Vecchi, “Field and source equivalence insource reconstruction on 3D surfaces,”
Prog. Electromagn. Res. , vol.103, pp. 67–100, 2010.[80] J. Kornprobst, R. A. M. Mauermayer, O. Neitz, J. Knapp, and T. F.Eibert, “On the solution of inverse equivalent surface-source problems,”
Prog. Electromagn. Res. , vol. 165, pp. 47–65, 2019.[81] N. Bleistein and J. K. Cohen, “Nonuniqueness in the inverse sourceproblem in acoustics and electromagnetics,”
J. Math. Physics , vol. 18,no. 2, pp. 194–201, 1977.[82] T. B. Hansen, A. Paulus, and T. F. Eibert, “On the condition number ofa normal matrix in near-field to far-field transformations,”
IEEE Trans.Antennas Propag. , vol. 67, no. 3, pp. 2028–2033, Mar. 2019.
Jonas Kornprobst (Graduate Student Member,IEEE) received the B.Eng. degree in electrical en-gineering and information technology from the Uni-versity of Applied Sciences Rosenheim, Rosenheim,Germany, in 2014, and the M.Sc. degree in electricalengineering and information technology from theTechnical University of Munich, Munich, Germany,in 2016. Since 2016, he has been a Research Assis-tant with the Chair of High-Frequency Engineering,Department of Electrical and Computer Engineering,Technical University of Munich. His current researchinterests include numerical electromagnetics, in particular integral equationmethods, antenna measurement techniques, antenna and antenna array design,as well as microwave circuits.
Alexander Paulus (Graduate Student Member,IEEE) received the M.Sc. degree in electrical engi-neering and information technology from the Tech-nical University of Munich, Munich, Germany, in2015. Since 2015, he has been a Research Assistantat the Chair of High-Frequency Engineering, Depart-ment of Electrical and Computer Engineering, Tech-nical University of Munich. His research interestsinclude inverse electromagnetic problems, compu-tational electromagnetics and antenna measurementtechniques.
Josef Knapp (Graduate Student Member, IEEE)received the M.Sc. degree in electrical engineer-ing and information technology from the TechnicalUniversity of Munich, Munich, Germany, in 2016.Since 2016, he has been a Research Assistant at theChair of High-Frequency Engineering, Departmentof Electrical and Computer Engineering, TechnicalUniversity of Munich. His research interests in-clude inverse electromagnetic problems, computa-tional electromagnetics, antenna measurement tech-niques in unusual environments, and field transfor-mation techniques.