Phase retrieval is an inverse problem of recovering the phase of a complex signal from its measured amplitude. It appears in various modifications in many scientific and engineering fields, including astronomy imaging DaiFie87 ; HarTho00 , X-ray crystallography Mil90 ; Har93 , microscopy Arr99 ; KimZhoGodPop16 and adaptive optics Mugnier2006 ; AntVer15 ; VisVer13 ; VisBruVer16 . An important application of phase retrieval in optics is to quantify the properties of an imaging system via its generalized pupil function Jan02 ; BraDirJan02 ; DoeThaVer18 ; PisGupSolVer18 . The fundamental advantage of this approach compared to those using intensity point spread functions (PSFs) or intensity optical transfer functions is that it is modifiable and automatically takes the specific characteristics of the imaging system under investigation into account. In adaptive optics, one needs to know the phase of the optical field in the system aperture to be able to compensate for an optical aberration, and the phase retrieval is a basis for a wide class of focal-plane based wavefront sensors.
Since the fundamental work Say52 of Sayre in 1952, which reveals that the phase of a scattered wave can be recovered from the recorded images at and between Bragg peaks of a diffracted wavefront, a wide variety of solution methods for phase retrieval has been proposed and developed. For an overview of phase retrieval algorithms, we refer the reader to the papers Fie13 ; SheEldCohChaMiaSeg15 ; Luk17 ; LukSabTeb19 . Direct methods usually require insights about the crystallographic structure to recover the missing phase Hau86 . Such a structural information is not only costly in terms of computational complexity but also sensitive to noise and approximation, for example, due to physical limitation or model deviation. As a consequence, this approach lacks practicability and becomes less popular in practice. The second class of solution algorithms relies on the fact that phase retrieval problems can be reformulated as linear equations with rank and positive semidefinite constraints in higher dimensional spaces. Well known examples of this algorithm class are MaxCut GoeWil95 , PhaseCut WalAspMal15 and PhaseLift CanEldStrVor13 ; CanStrVor13 . This convex relaxation approach requires the matrix lifting step which is computationally demanding and hence not suitable for large-scale problems. The most popular class of phase retrieval methods is based on projections and pioneered by the work of Gerchberg and Saxton GerSax72 , which deals with phase retrieval given a single PSF image and the amplitude of the complex signal, which in the sequel will be referred to as the amplitude constraint in order to clearly differentiate it from the intensity constraints determined by data images. The need to deal with more and more phase retrieval models, for example, incorporating various types of a priori constraint Fie82 , being given multiple images and involving regularization schemes, has given rise to a wide range of solution methods in this class. It was recently observed by Luke et al. LukSabTeb19 that this class of methods actually outperforms the other classes of phase retrieval algorithms.
In light of LevSta84 ; BauComLuk02 ; LukBurLyo02 , phase retrieval can be interpreted as mathematical feasibility problems and, as a consequence, all algorithmic schemes for set feasibility can be adapted for phase retrieval. The current research is devoted to that topic. The alternating projection (AP) and the Douglas–Rachford (DR) algorithms are perhaps the most widely known solution methods for set feasibility and have served as a basis for a wide range of modifications and regularizations, see, for example, BauMou17 ; KruLukNgu18 . It has been observed that AP is stable, always convergent and to some extent able to suppress noise but it may get stuck at undesired local minima and the convergence speed can be very slow Fie82 . In contrast, DR can be faster in convergence and better in escaping from bad local minima but less robust against noise and model deviation Luk08 . As a result, this algorithm can not be naively applied to practical problems which intrinsically involve noise and model approximation. This fact has motivated a number of its efficient relaxation schemes such as the usage of the Krasnoselski–Mann relaxation, the Fienup’s hybrid input-output (HIO) algorithm Fie82 , the relaxed averaged alternating reflections (RAAR) algorithm Luk05 ; Luk08 and the DRAP algorithm Tha18 .
In this paper, we analyze the DRAP algorithm for solving the phase retrieval problem for the first time after having observed that it appears to be the most efficient algorithm for the problem setting under consideration, see Section 5. Interestingly, DRAP mathematically coincides with the convex combination of the AP and DR operators in the phase retrieval setting. As a result, DRAP admits two mathematically equivalent descriptions (see (19) and (20) in Section 3). The first one ensures that its computational complexity is only approximate to that of each of the constituent operators and thus it is used for numerical implementation. The second description as a convex combination of the AP and the DR operators exhibits a concrete connection to the fundamental projection algorithms and hence it is intuitively better situated on the map of projection methods (see Remark 6).
The main contribution of this paper is the convergence analysis of the DRAP algorithm for solving the phase retrieval problem. First, using the analysis approach initiated by Chen and Fannjiang CheFan18 , we establish a convergence criterion for DRAP (Theorem 4.1), which extends the convergence result of the DR algorithm formulated in that paper. It is worth mentioning here that extending a convergence criterion for DR to a corresponding one for its relaxations such as HIO, RAAR and DRAP algorithms is not trivial111For example, similar criterion for RAAR was proved in LiZho17 while the one for HIO remains unknown.. Proposition 2 extends the applicable scope of this type of convergence results222Including the criteria for DR and RAAR algorithms formulated in CheFan18 and LiZho17 , respectively. to cover also phase retrieval problems with amplitude constraint. Second, applying the analysis scheme developed by Luke et al. LukNguTam18 , we establish another convergence criterion for the DRAP algorithm (Theorem 4.2) by integrating the physical properties of the phase retrieval problem Luk08 into the earlier known results for DRAP Tha18 . Recall that the analysis of the latter article involves only abstract mathematical notions in the general setting of set feasibility. As a comparison, we make an attempt on connecting the two convergence criteria by linking their key mathematical assumptions to a single physical condition on the phase diversities which are the almost only adjustable figures of the phase retrieval problem (see Remark 15).
The paper is organized as follows. In the last part of this introductory section, we introduce the mathematical notation used in the paper. Section 2 is devoted to formulating the phase retrieval problem and addressing in details the key steps towards its solutions using projection algorithms. A discussion on projection methods for phase retrieval is presented in Section 3. In Section 4, convergence results of the DRAP algorithm are established using two different analysis approaches: 1) spectral analysis in Section 4.1 and 2) variational analysis in Section 4.2. Numerical simulation is presented in Section 5.
Mathematical notation. The underlying space in this paper is a finite dimensional Hilbert space denoted by . The element-wise multiplication is denoted by . The element-wise division , absolute value , square and square root operations are also frequently used but without need for extra notation. and denote the real and the imaginary parts of a complex object in th brackets, respectively. The imaginary unit is . denotes the identity mapping while denotes the identity matrix of size . The distance to a set is defined by
and the set-valued mapping
is the projector on . A selection is called a projection of on . The reflection operator associated with is accordingly defined by . Given a subset , the Fréchet and limiting normal cones to at a point are defined, respectively, as follows:
where means that and . The set of fixed points of an operator is defined by . Our other basic notation is standard; cf. Mor06.1 ; VA . stands for the open ball with radius and center . For a linear subspace of ,
is the orthogonal complement subspace of .
2 Problem formulation
2.1 Phase retrieval
Phase diversities and the Fourier transform are key ingredients of the phase retrieval problem studied in this paper. Recall that adding a phase diversity to the phase of a complex signal is a unitary transform and the (discrete) Fourier transform is also a unitary operator. Since unitary transforms are one-to-one represented as unitary matrices, the phase retrieval problem can be formulated in the form of matrix-vector-multiplication as follows. For an unknown complex object, let be the propagation matrix which is normalized to be isometric, and be the measured data of . The phase retrieval problem is to find an (approximate) solution to the equation:
where represents unknown noise333Dimension corresponds to the pixel totality of one image..
To formulate the phase retrieval problem in the matrix-vector-multiplication form (2) or any feasibility model in Section 2.2, we need to vectorize all array objects in a consistent manner and rewrite all linear mappings as matrix multiplication operations in higher dimensional spaces, see, for example, (DoeThaVer18, , section 2A). This one-to-one conversion allows us to do the theoretical analysis in the simple matrix-vector-multiplication formulation without loss of generality.
In this paper, we study the phase retrieval setting with several phase diversities, and the propagation matrix takes the following form:
where is the number of data images,
is the unitary matrix representing the discrete Fourier transform, andare unitary matrices representing the phase diversities which will be denoted by in the sequel . Note that .
Remark 2 (phase modulators versus out-of-focus measurements)
There are two widely used techniques of acquiring the PSF images for phase-diversity phase retrieval. First, a phase modulator is used for introducing phase diversities in the pupil plane corresponding to which the images are measured in the focal plane. Second, the images are registered in out-of-focus planes along the optical axis (i.e, parallel to the focal plane at some known distances) without the use of phase modulator. It is well known that the two techniques are mathematically equivalent Goodman05 .
When a priori knowledge of the solutions is available, that is for some known subset , one can expect more accurate phase retrieval. The formulation (2) is naturally modified as follows:
2.2 Feasibility models
Several feasibility models of phase retrieval have been formulated in either the physical domain444The unknown variable is the signal in the pupil plane. LukBurLyo02 ; LevSta84 or the Fourier domain555The unknown variable relates to the signal in the pupil plane via the Fourier transform. CheFan18 . Viewing the Fourier transform and phase-diversity addition as unitary transforms, we clarify the relationship between various feasibility models of the phase retrieval problem.
Then, the problem (4) can be approached via the following feasibility problem involving multiple sets:
where captures a priori knowledge of the solutions.
Remark 3 (nonconvexity feasibility)
All the problem models appearing in this paper are nonconvex due to the nonconvexity of the intensity constraints defined in (5).
When addressing the phase retrieval problem with noise and model deviation, an appropriate averaging process is essential for suppressing noise. For this, we consider the following feasibility model in the product space:
The equivalence between (6) and (7) in the general setting of set feasibility finds its root in Pie84 . Without a priori constraint, i.e., , the set is the (-dimensional subspace) diagonal of the product space . The counterpart of (7) in the Fourier domain is as follows:
Proposition 1 (equivalences of feasibility models)
Remark 4 (inconsistent feasibility)
In practical circumstances, for example, due to the presence of noise and model deviation, the intersection in (6), (7) and (9) is likely to be empty. There are natural interpretations of inconsistent feasibility in terms of minimization involving indicator and distance functions. For example, let us interpret the AP method for solving the (possibly inconsistent) feasibility (9) in terms of classical algorithms for minimization. The worrisome issue regarding the emptiness of the intersection would be eased when one associates (9) with the following minimization problem:
In view of Proposition 2 (which is proved later in Section 4.1), the set defined in (10) can be assumed to be convex, and hence the objective function in (11) is differentiable with the gradient given by for every point PolRocThi00 . Then, alternating projection for solving (9) is precisely the projected gradient method for solving (11).
The decisive step of solving the feasibility problem (9) by projection algorithms is to calculate the two projectors on the sets and defined in (10). Since is geometrically the product of a number of circles of the complex number plane, an explicit form of the projector , which is in general a set-valued mapping, is available BauComLuk02 ; LukBurLyo02 :
with the convention that whenever , where denotes the complex unit circle666The subscript indicates the th entry of the object.. In numerical computation, the (single-valued) selection of corresponding whenever is sufficient.
Remark 5 (projector on regularized sets)
In view of Remark 4, the set can have no common point with the set . For ways of handling such a feasibility gap, one can think of regularizing or approximating the set . For example, Luke Luk12 proposed to enlarge the set to
where can be viewed as the radius of enlargement, is the Bregman distance, associated with a strictly convex function which is differentiable on the interior of its domain, given by
The function should be chosen in accordance with the statistical model of the noise in (4). More specifically, let us consider the Gaussian and Poisson models of noise, which are perhaps the most relevant to phase retrieval. The Bregman distance associated with the half energy kernel operator is simply the Euclidean norm, and it is appropriate for Gaussian noise. Let us define the function by
The Bregman distance associated with the function given by (13) is the Kullback-Leibler divergence, and it is appropriate for Poisson noise. The projector on the regularized set can be viewed as an approximation of the projector on , and hence it can be used in the framework of projection methods. The cyclic projection algorithm using approximate projectors of this type has been analyzed by Luke Luk12 , and in fact his idea can also be extended to other projection methods. However, since the projector on a regularized set is often much more complicated to be computed than the one on the original set, we can instead treat the latter one as an approximation of the former one (LukSabTeb19, , page 22). This insight about approximate projectors for inconsistent feasibility allows us to simply use the formula (12) for both analytical and numerical purposes without any worrisome issue.
The projector on the set can also be explicitly described777Note that convexity of is not required in Lemma 1.. We make use of the following notation:
For the propagation matrix given in (3), it holds that
Let us first define the unitary matrix based on the matrix as follows:
This block diagonal matrix is unitary since all of its constituent blocks are so. By the structure of and , we have that
Since is unitary, it holds that
Since is a subspace containing , by the properties of the metric projection, we have that
We next calculate . Note that is also a block diagonal matrix whose blocks are the conjugate transpose of the corresponding blocks of . Let us denote the column vector whose entries taken from correspond to the block of , . We have that
Since is the -dimensional diagonal of the product space , we obtain by solving the minimizing problem (1) that
The proof is complete. ∎
The formula (14) shows that the complexity of heavily depends on that of .
3 Projection algorithms
Projection algorithms for phase retrieval can be considered as descendants of the well known Gerchberg–Saxton (GS) algorithm GerSax72 which deals with phase retrieval given the amplitude constraint and a single PSF image. Their introduction has been motivated by the rapidly growing application of phase retrieval originated from a wide variety of physical settings. For example, the famous input-output, output-output and hybrid-input-output algorithms Fie82 arose up when dealing with the support and the real and nonnegative constraints instead of the amplitude constraint as the GS method. Extensions for solving problems given multiple images and for obtaining better restoration have been among the main objectives of this class of phase retrieval algorithms. In light of the groundwork BauComLuk02 , in Section 2.2 we have interpreted the phase retrieval problem (4) as a feasibility problem in one of the equivalent forms (6), (7) and (9). Having calculated the projectors and in Section 2.3, we are now ready to discuss algorithmic schemes for the solutions. From now on, we analyze the feasibility model (9).
AP and DR are perhaps the most widely known solution methods for feasibility and have been the basis for a wide variety of modification and regularization schemes. We refer the reader to, for example, KruLukNgu18 ; BauMou17 for an overview of these basic methods in the setting of set feasibility. For an early discussion in the context of phase retrieval, we refer the reader to the surveys LukBurLyo02 ; BauComLuk02 . It has been observed that AP is stable, always convergent and to some extent able to suppress noise, but the convergence speed can be very slow Fie82 . In contrast, DR can be fast in convergence, but sensitive to noise and model deviation Luk08 . Indeed, only relaxations of DR can be used for problems in the presence of noise and model mismatch.
The use of the Krasnoselski–Mann relaxation is perhaps the most widely known. Mathematically, it is the convex combination of the DR operator and the identity mapping:
where is the relaxation parameter. The Fienup’s hybrid-input-output (HIO) method Fie82 can be viewed as a relaxation of DR:
where is the relaxation parameter. Another relaxation of DR known as the relaxed averaged alternating reflections (RAAR) algorithm was proposed and analyzed in Luk05 ; Luk08 for phase retrieval. It is the convex combination of the DR operator and one of the projectors:
where is the relaxation parameter888Relaxation parameter zero is not allowed for KMDR, HIO and RAAR..
Interestingly, in the phase retrieval setting (9), coincides with the convex combination of the AP and DR operators provided that is an affine set. The latter condition implies that the set given by (10) is affine. Hence, the projector is linear and we obtain the the following expression:
where is the AP operator.
The two expressions (19) and (20) play their own role in explaining interesting features of DRAP999They do differ in general settings.. On the one hand, only two projections are required for computing an iteration of (19) ( once and once) compared to three projections for (20) ( once and twice). This means that the computational complexity of DRAP is at the same level as that of the other projection methods if (19) is used in numerical implementation. On the other hand, the expression (20) as a convex combination of and explains better the idea leading to the introduction of DRAP as a relaxation of DR compared to the less intuitive form (19).
where and stand for the two consecutive iterations of DRAP. In the case , (21) further reduces to
In the remainder of this paper, we analyze the DRAP algorithm in the phase retrieval setting (9) and demonstrate its advantages over the other algorithms.
4 Convergence analysis
In this section, we study convergence properties of DRAP using two different analysis schemes. Since the problem (9) is nonconvex, we can only obtain local convergence criteria though it is observed from numerical results that the quality of phase retrieval is not affected by the starting point for the algorithm.
4.1 A result from spectral analysis
The analysis in this section is based on the observation that the projector given by (14) is linear, and the projector given by (12) also has a good first order approximation around any solution of (9). We follow the analysis approach initiated by Chen and Fannjiang CheFan18 where they established a local linear convergence result for the DR algorithm. The mentioned result of CheFan18 was later extended for the RAAR algorithm in LiZho17 . We will show that DRAP also enjoys that kind of convergence result101010Similar results for the HIO algorithm are unknown.. In this section, we assume that the lowest intensity of the images is strictly positive:
where is a solution to (9) and denotes the diagonal matrix with elements on its diagonal taken from the vector in the brackets. Since vanishes nowhere111111Recall that the square amplitude is element-wise. by (23), for all sufficiently close to , also vanishes nowhere. In particular, for a fixed vector , the vector vanishes nowhere provided that is sufficiently small. The next lemma establishes the first order approximation of as a complex vector valued function around in a given direction.
Lemma 2 (first order approximation of )
For a vector and a sufficiently small number , we have
where and .
The next step is to analyze the spectrum of the real decomposition of the complex matrix as follows:
Note that is isometric since is so. Define also the real decomposition of a complex vector by
be the singular values ofwith the corresponding right singular vectors and the left singular vectors
. We have by the definition of the singular value decomposition (SVD) that
The next technical result regarding the spectrum of is crucial.
(CheFan18, , Proposition 5.6) There holds that , , , and .
Thanks to Lemma 3 and the definition of the SVD, one has the following expression of the second largest singular value of :
The following theorem establishes linear convergence of the DRAP algorithm for solving (9). Since phase retrieval is ambiguous (at least) up to a global phase shift121212That is the first element of the orthogonal basis of Zernike polynomials., the following distance between two complex vectors is of interest:
Theorem 4.1 (linear convergence of DRAP)
Let us denote . Thanks to Lemma 2, we have that
Multiplying both sides of the above equality by and taking the isometry property of into account, we obtain that
Due to (30) and the fact that we have
In other words, . By basic properties of the Hermitian inner product, one has . As a result, . Taking Lemma 3 into account, we have just shown that is orthogonal to which is the first left singular vector of . This together with the expression (27) of implies that
Since by assumption (29), there exists a number such that for all with sufficiently small, it holds that
The proof is complete. ∎
Remark 9 (region of convergence)
Since the algorithm operates in the underlying space , for the sake of brevity, let us speak of region around instead of . In view of Theorem 4.1, such a convergence region, if exists, is mutually dependent on the constant . More specifically, given a number , it is the region in which the first order approximation (24) of around is valid and condition (34) is satisfied for all . Note that the latter involves not only and but also the sequence itself. The intersection of the regions over all possible sequences complied with can be taken as the region of convergence. Obviously, such a statement is not informative and hence it has not ever been an objective of local convergence analysis.
Remark 10 (influence of on convergence)
In view of Theorem 4.1, the relaxation parameter obviously has influence on the region in which the first order approximation of (Lemma 2) is valid and condition (34) is satisfied, however, its influences on the convergence speed of DRAP is unclear131313Of course, the influence is clearly observed from numerical computation..
We have analyzed the DRAP algorithm in the phase retrieval setting (9) with . The latter condition limits the effectiveness of Theorem 4.1 to phase retrieval without a priori constraint. In the remainder of this section, we will show that the convergence criterion can also be applicable to phase retrieval problems with an amplitude constraint, which is a helpful prior information and often available in practice141414This is because the light distribution in the pupil plane is often known, for example, it can be uniform or truncated Gaussian..
The amplitude constraint is described by
where is the known amplitude of the complex signal. The next result shows that the problem (9) with an amplitude constraint can equivalently be reformulated as a problem without a priori constraint in a higher dimensional space.