An SVD Approach for Estimating the Dimension of Phaseless Data on Multiple Arcs in Fresnel Zone

: In this article, we tackle the question of evaluating the dimension of the data space in the phase retrieval problem. With the aim to achieve this task, we ﬁrst exploit the lifting technique to recast the quadratic model as a linear one. After that, we evaluate analytically the singular values of the lifting operator, and we quantify the dimension of the data space by counting the number of “signiﬁcant” singular values. In the last part of the article, we show some numerical results in order to corroborate our analytical prediction on the singular values’ behavior of the lifting operator and on the dimension of the data space. The analysis is performed for a 2D scalar geometry consisting of an electric current strip whose square magnitude of the radiated ﬁeld is observed on multiple arcs of circumference in Fresnel zone.


Introduction
Inverse problems cover a wide range of applications in the electromagnetic and optical literature [1][2][3][4][5][6]. However, especially when the frequency increases, accurate phase measurements may become difficult; hence, the need to perform reconstructions from only intensity data may arise. Such a task falls into the realm of the phase retrieval problem.
From a mathematical perspective, phase retrieval is an inverse problem that consists of retrieving the unknown from the square magnitude (or from the magnitude) of its transformation. In its quadratic formulation, it can be formalized as recovering the unknown function f from the equation: where T is an operator such that f ∈ X −→ g ∈ Y with X and Y denoting two functional spaces. As can be seen from (1), the absence of the phase knowledge makes the problem nonlinear also when the operator T is linear. The nonlinearity of the problem implies that the corresponding least squares problem involves a quartic cost functional [25] that generally is non-convex. Indeed, the latter may contain trap points like local minima and saddle points, which make the task of finding the actual solution harder.
In this framework, deterministic optimization techniques do not ensure the recovery of the actual solution, since they converge to the closest stationary point. To overcome this drawback, globally convergent techniques like genetic algorithms, particle swarm optimization, and simulated annealing can be exploited. Despite this, such approaches exhibit a very slow converge speed that makes them unsuitable for large-scale problems.
To overcome the traps issue, the lifting technique was introduced [26]. The latter exploits a redefinition of the unknown space, which allows recasting the phase retrieval problem as a linear one. However, since in all lifting-based methods, the new number of unknowns is equal to the square of the original number of unknowns, their range of applicability is restricted to small-size problems.
In light of the previous considerations, the reduced computational burden and the low memory requirements of the approach based on the minimization of the quartic functional have led to a revival of such a method and have prompted researchers to investigate the conditions under which it is possible to ensure the convergence. In particular, in these studies, it has been shown that the convergence to the actual solution can be ensured by exploiting two different strategies. The first one is that of using a special initialization, which provides a starting point in the attraction basin of the cost functional, if some mathematical conditions involving the ratio between the dimension of the data space (M) and the dimension of the unknowns space (N) are satisfied [27][28][29]. The second strategy is based on the observation that a sufficiently high value of the ratio M/N makes the quartic functional free from traps [30][31][32]. In such a case, no special initialization is required, and the starting point can be chosen at random.
From the previous discussion, it is evident that the dimension of the data space (M) plays a key role for convergence guarantees; hence, after recalling its definition, it is worth investigating how to estimate it from an analytical point of view.
The dimension of the data space can be defined as the number of independent functions that allow representing the data with a given degree of accuracy. It can be estimated in two steps. The first one is that of exploiting the lifting technique to obtain a linear representation of the data. Once such a task has been achieved, the dimension of the data space can be computed by employing different tools like Singular-Value Decomposition (SVD) [1], the Gram-Schmidt orthogonalization (GSO) [25], or Principal Component Analysis (PCA) [33,34]. The main difference between such methods is that the SVD and GSO estimate the dimension of the data space studying the lifting operator, instead, PCA reaches the same goal studying the data of the problem.
In this article, we choose to evaluate the dimension of the data space by exploiting the approach introduced in [35]; hence, we estimate such a parameter by counting the number of "significant" singular values of the lifting operator. For such a reason, with reference to a 2D scalar geometry concerning a strip current observed on multiple arcs in Fresnel zone, we will provide an approximate expression of the singular values of the pertinent lifting operator.
Let us remark that analytical results regarding the singular values of the linear operator that links the source current with the radiated field in amplitude and phase are already available in the literature [36][37][38]. Conversely, to the best of our knowledge, this is the first time that analytical results concerning the singular values' behavior of the lifting operator are provided.
The paper is structured as follows. In Section 2, we sketch the geometry of the problem and recall some preliminary results on the corresponding problem with the data in amplitude and phase. In Sections 3 and 4, we study respectively the singular values' behavior of the lifting operator in the case of one and more observation arcs. In Section 5, we validate the theoretical results by means of some numerical experiments. The conclusions follow.

Geometry of the Problem and Preliminary Results
In this article, we consider the 2D scalar geometry depicted in Figure 1 where the y-axis represents the axis of invariance. An electric current J(x) = J(x)î y supported on the set [−a, a] of the x-axis radiates within a homogeneous medium with wavenumber β.
The electric field E radiated by such a strip source has one component directed along the y-axis; hence, E(x, z) = E(x, z)î y . The square amplitude of the electric field is observed in the Fresnel zone on P arcs whose radii are denoted by r 1 , r 2 , ... , r P . Such arcs subtend the same angle; indeed, they extend along the polar coordinate θ on the same set [−θ max , θ max ].
For the geometry at hand, the radiated electric field on the i-th arc of circumference can by expressed in the variable u = sin(θ) by the equation: where ∀i = {1, 2, .., P}, T i is the linear integral operator that realizes the mapping: with L 2 (SD) denoting the space of square-integrable functions on the set SD = [−a, a] and L 4 (OD i ) indicating the space of functions whose amplitude to the fourth power is integrable on the set OD i = [−u max , u max ] = [−sin θ max , sin θ max ]. Under the paraxial Fresnel approximation, the operator T i can be explicitly written in the form: where: With the aim to simplify the discussion of the next sections, let us recall some results on the linear integral operator T i T † i where T † i stands for the adjoint of T i . Such an operator can be expressed as below: where: Accordingly, the operator T i T † i is convolution, and its kernel is a band-limited function of the sinc type.

Singular Values of the Lifting Operator in the Case of One Observation Arc
Before addressing the case of P observation arcs, let us analyze the configuration where the observation domain is made up by one arc of radius r 1 .
The square amplitude of the electric field over such an arc, |E 1 (u)| 2 , is given by: Starting from the quadratic model in (8), at first, we provide a linear representation of |E 1 | 2 . After, by studying the singular values of such linear operator, we find the dimension of the data space.
In order to obtain a linear representation of |E 1 | 2 , let us rewrite (8) as: where (·) * denotes the conjugate of the function (·). From Equation (9), it is evident that if we redefine the unknown space and consider as unknown the function: then the operator A 1 that links the unknown function F(x, x) with the data function |E 1 (u)| 2 is linear. The latter is known in the literature as the lifting operator, and it is defined as: with: Accordingly, the square amplitude distribution |E 1 | 2 can be written through the following linear model: By Equation (12), it follows that the adjoint operator A † 1 is given by: where (·) denotes the function of the variable u on which the adjoint operator acts. In order to evaluate the singular values σ m of the lifting operator A 1 , we find the eigenvalues λ m of the operator A 1 A † 1 , which are related to each other by the equation: The operator A 1 A † 1 can be expressed as: where: Hence, the kernel of A 1 A † 1 is the square of the kernel of T 1 T † 1 . Taking into account (7), H 11 (u, u) can be written as: Accordingly, the operator A 1 A † 1 can be expressed in the following form: The eigenvalues of a convolution operator with a a sinc-squared kernel are known in the closed-form [39]. Hence, the eigenvalues of A 1 A † 1 can be analytically evaluated, and they are given by: where: with [ 4 π βa u max ] denoting the integer nearest to 4 π βa u max . Consequently, the singular values of A 1 are given by: The value of M 1 provides the number of relevant singular values of the lifting operator; hence, it can be taken as an evaluation of the dimension of the data space.

Singular Values of the Lifting Operator in the Case of P Observation Arcs
In this section, we consider a configuration where the observation domain is an ensemble of P observation arcs in the Fresnel zone. In such a case, the square amplitude distributions |E 1 | 2 , |E 2 | 2 , ..., |E P | 2 are related to the current distribution J by the quadratic model: With reference to such a configuration, at first, we introduce a linear operator that represents the phaseless data |E 1 | 2 , ..., |E P | 2 . After, we evaluate analytically the singular values of such an operator and the dimension of the data space.
As already done in Section 3, here, we consider as unknown the function F(x, x) = J(x) J * (x). This redefinition of the unknown function allows expressing the square amplitude distributions of the radiated field through the following linear model: where ∀i ∈ {1, ..., P} Accordingly, the linear operator A, which represents the phaseless data in the case of P observation arcs, is defined as: Such an operator is still called the lifting operator, and it can be expressed as: The adjoint operator A † is given by: where ∀i = {1, ..., P} To estimate the singular values σ m of the lifting operator A, let us find the eigenvalues λ m of the operator AA † . The latter is given by: where ∀(i, j) ∈ {1, ... , P}: with: With the aim to evaluate H ij (u, u), let us distinguish the case i = j from the case i = j. For i = j, the integral in (32) can be easily evaluated, and it results that: Differently, for i = j, the kernel H ij cannot be expressed in terms of elementary functions. In such a case, a closed-form expression of the integral (32) can be provided by resorting to the imaginary error function, which will is denoted by erfi (for further details, see Appendix A).
In particular, as shown in [40], it results that: Accordingly, for i = j, H ij (u, u) can be expressed as: where: Let us point out that: • the operators along the main diagonal of (30) and those off-diagonal are both convolutions. • ∀(i, j) ∈ {1, ..., P} the kernel of A i A † j becomes more and more similar to A i A † i when r j approaches r i .
As said before, our goal is that of evaluating the singular values of the lifting operator A by computing the eigenvalues of AA † . Naturally, the behavior of such eigenvalues depends on the distance between the arcs on which the square amplitude of the electric field is observed.
In particular, if the radii of the arcs are very similar (r 1 ≈ r 2 ≈ ... ≈ r P ), then all the operators in (30) have approximately the same kernel. Accordingly, the number of relevant eigenvalues of AA † is approximately equal to the case of one observation arc.
The situation changes completely when the radii of the observation arcs are sufficiently different from each other. Indeed, in this case, each scanning increases the number of relevant eigenvalues and, consequently, the dimension of the data space.
With the aim to evaluate the number and the value of the relevant eigenvalues of AA † , let us remark that if the norm of the off-diagonal terms in (30) is negligible with respect to the norm of the terms along the main diagonal, then the eigenvalues of AA † are approximately given by the union of the eigenvalues of A 1 A † 1 , A 2 A † 2 , ..., A P A † P . From the previous discussion, it is evident that ∀(i, j) ∈ {1, ..., P}, the norm of A i A † j plays a key role; hence, it is worth evaluating such a quantity. Since ∀(i, j) ∈ {1, ..., P}, the operator A i A † j is self-adjoint, and its norm can be computed according to the equation: Hence, to compute the norm of A i A † j , it is necessary to estimate its eigenvalues. For each i = j, the eigenvalues of A i A † j can be analytically evaluated, since the eigenvalues of a convolution operator with a sinc squared kernel are provided in [39].
Conversely, for each i = j, the eigenvalues of A i A † j are not known in the closed-form. Despite this, the following considerations can be made to understand the eigenvalue behavior of such an operator.
The kernel of A i A † j is provided by (35). The latter exhibits an oscillating behavior with a main lobe and a number of side lobes. The parameters p and q that appear in the arguments of the erfi functions are related respectively to the width of the main lobe and to the level of the side lobes of H ij (u, u). In particular, it happens that if the distance between r i and r j increases, then: • the parameter p decreases, • the parameter q increases.
Consequently, the main lobe of H ij (u, u) enlarges, while the level of the side lobes rises. Such changes make the original main lobe indistinguishable by the first side lobes; hence, the total effect of increasing the distance |r i − r j | is to raise the width of the main lobe of H ij (u, u). Since the eigenvalue behavior of a convolution operator is asymptotically given by the Fourier transform of its kernel, when the main lobe of the kernel becomes wider, the support of its Fourier transform is restricted, and the number of relevant eigenvalues decreases. From the previous discussion, it follows that if the distance between the radii r i and r j increases, then the eigenvalues of A i A † j decay more quickly, and consequently, the norm of A i A † j is reduced.
Hence, for each i = j, a minimum value of the distance |r i − r j | exists such that ||A i A † j || is negligible with respect to ||A i A † i || and ||A j A † j ||. In such a condition, the operator AA † in (30) can be approximated as: Now, it is relevant to quantify the minimum distance between the radii such that the operator AA † can be approximated as in (38). To this end, it is worth noting that if we introduce a weighted adjoint A † and we repeat for the considered geometry the same reasoning made in [41], we obtain that the (i, j) entry of the operator AA † can be approximated as below: From Equation (39), it is evident that the off-diagonal terms of AA † are negligible with respect to the terms on the main diagonal when the condition: is fulfilled. Let us remark that the use of the weighted adjoint changes only the dynamics of the eigenvalues, but not the critical index at which they become negligible. For such a reason, the last equation also works for the operator A i A † j . Accordingly, we can state that if Condition (40) is satisfied, then ||A i A † j || is negligible with respect to ||A i A † i || and ||A j A † j ||. From Equation (40), it is clear that the fulfillment of such an equation is more critical when j = i + 1. Therefore, we can conclude that if the quantity ∆r i = r i+1 − r i satisfies the inequality: ∀i = 1, .., P − 1, then the operator AA † can be approximated as in (38). In such a case, the eigenvalues of AA † are given by: where ∀i ∈ {1, 2, ..., P}: with M i = 4 π βa u max . Consequently, the singular values of the lifting operator A can be computed by exploiting the equation: where ∀i ∈ {1, 2, ..., P}: From the previous discussion, it follows that if the spacing ∆r i between the observation arcs satisfies the inequality (40) ∀i = {1, 2, .., P}, then the number of relevant singular values of the lifting operator A or, in other words, the dimension of the data space is approximately given by:

Numerical Verification of the Theoretical Results
In this section, we check the results on the singular values' behavior provided in Sections 3 and 4 by means of some numerical simulations.

Numerical Verification in the Case of One Observation Arc
In this section, we check the analytical results provided in Section 3 for the case of one observation arc of radius r 1 = 25λ. To this end, in Figure 2, we compare the numerical evaluation of the singular values of the lifting operator A 1 with the theoretical prediction provided by (22). As can be seen from Figure 2, the two diagrams exactly overlap until the index M 1 = 40, which provides an accurate estimation of the dimension of the data space.

Numerical Verification in the Case of Two Observation Arcs
In this section, we verify the analytical results of Section 4 in the case of two observation arcs. In such a case, the operator AA † is a two by two matrix given by: In Figures 3 and 4, with reference to the case where r 1 = 25λ, we show respectively: • a cut of the kernel of A 1 A † 2 , • the eigenvalues of A 1 A † 2 , for three different values of the radius r 2 . In particular, in the two figures, the black diagram refers to the case r 2 = 27λ, while the red and the blue diagrams refer respectively to the configurations r 2 = 30λ and r 2 = 35λ.  In perfect agreement with the theoretical discussion, from Figures 3 and 4, it is evident that when the distance ∆r 1 = r 2 − r 1 increases, then: 1.
the width of the main lobe of H 12 is larger, 2.
the eigenvalues of A 2 A † 1 decay more quickly, and consequently, ||A 2 A † 1 || is lower. For the considered configuration (a = 10λ, u max = 0.5, r 1 = 25λ), the inequality (41) is satisfied if ∆r 1 ≥ 8.33λ. In Figure 5, we sketch the singular values of A numerically computed and the theoretical prediction of such singular values provided by (44) in the case where r 2 = 27λ (∆r 1 = 2λ). Instead, in Figure 6, we show the same quantity with reference to the case r 2 = 35λ (∆r 1 = 10λ).  As can be seen from Figure 5, when ∆r 1 = 2λ, the inequality (41) is not satisfied, and the theoretical prediction of (44) does not work. Conversely, when ∆r = 10λ (see Figure 6), the inequality (41) is satisfied, and the theoretical estimation of the singular values (44) works very well. In this last case, also the dimension of the data space is correctly provided by (46); indeed, as can be seen from Figure 6, the singular values of A are relevant until the index M = [ 4P π βa u max ] = 80.

Numerical Verification in the Case of Three Observation Arcs
In this last section, we verify the analytical results of Section 4 in the case of three observation arcs.
As regards the value of r 1 and r 2 , we choose r 1 = 25λ, r 2 = 35λ. In such a case, it results that ∆r 1 = r 2 − r 1 = 10λ, and the condition (41) is fulfilled.
As concerns r 3 , we choose it in such a way that also ∆r 2 = r 3 − r 2 satisfies the condition on the spacing of the arcs. In particular, the minimum value of ∆r 2 = r 3 − r 2 such that condition (41) is satisfied is 18.85λ, and this implies that r 3 ≥ 53.85λ. In our numerical experiment, we choose r 3 = 55λ.
In Figure 7, we compare the singular values of A numerically computed with their theoretical prediction provided by (44) with reference to the case a = 10λ, r 1 = 25λ, r 2 = 27λ, r 3 = 55λ, u max = 0.5. Figure 7. Numerical evaluation and theoretical prediction of the singular values of A for the configuration a = 10λ, r 1 = 25λ, r 2 = 35λ, r 3 = 55λ. The two diagrams are in dB, and they are normalized with respect to the maximum of the numerical evaluation.
As illustrated in Figure 7, the singular value behavior of the lifting operator and the critical index at which the singular values become negligible are both predicted with sufficient accuracy by (44) and (46), respectively.
The main difference between the two diagrams arises in correspondence with the knee of the curve of singular values. Indeed, as can be seen from the figure, the index at which the knee occurs in the black diagram is lower than the index at which the knee occurs in the red one. This means that when the number of arcs increases, the theoretical estimation of the data space dimension provided by (46) is not exactly equal to the actual value of the dimension of the data space, but it represents a good upper-bound.

Conclusions
In this paper, we introduce a strategy that allows finding the dimension of the data by means of the singular-value decomposition of the lifting operator.
The study was done with reference to a 2D scalar configuration made up by a strip source whose square magnitude of the radiated field is observed on multiple arcs of circumference in the Fresnel zone.
In particular, for the considered geometry, we first introduced the lifting operator. After that, we found its singular values, and finally, we established the dimension of the data space by finding the index at which the singular values of the lifting operator become negligible.
Before concluding, it is worth point out how the noise in the data affects the dimension of the data space. To clarify this aspect, let us recall that the lifting operator is compact; accordingly, its singular values approach zero as the index increases. Furthermore, as shown in Sections 3 and 4, the singular values of the lifting operator are not step-like, but they exhibit a dynamics before becoming negligible in correspondence with the critical index [ 4P π βa u max ]. This implies that the dimension of the data space is not independent of the noise level, but weakly dependent on it. In other words, if the SNR (Signal-to-Noise Ratio) is high, then the data space dimension M is essentially equal to [ 4P π βa u max ]. Conversely, if the SNR is low, then the dimension of the data space M is smaller than [ 4P π βa u max ]. Naturally, the noise in data affects not only the dimension of the data space, but also the inversion process, making the use of a regularized inversion scheme mandatory. However, this last issue is out of the scope of this paper.
As a final remark, let us point out that our SVD-based approach for the estimation of the data space dimension can also be extended to near zone configurations and more realistic scenarios involving 3D geometries.

Acknowledgments:
The authors would to thank Raffaele Solimene for the useful advice during the revision process and Ehsan Akbari Sekehravani for proofreading the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this Appendix, we provide some details on the erfifunction. The imaginary error function er f i(z) is an entire function defined by the equation: with er f (z) denoting the erf function (also called the error function). The derivative and the integral of the imaginary error function are given by: d dz er f i(z) = 2 √ π e z 2 (A2) er f i(z)dz = z er f i(z) − 1 √ π e z 2 (A3)