Absolute Uniqueness of Phase Retrieval with Random Illumination
Albert Fannjiang
Introduction
Phase retrieval is a fundamental problem in many areas of physical sciences such as X-ray crystallography, astronomy, electron microscopy, coherent light microscopy, quantum state tomography and remote sensing. Because of loss of the phase information a central question of phase retrieval is the uniqueness of solution which is the focus of the present work.
Researchers in phase retrieval, however, have long settled with the notion of relative uniqueness (i.e. irreducibility) for generic (i.e. random) objects, without a practical means for deciding the reducibility of a given (i.e. deterministic) object, and searched for various ad hoc strategies to circumvent problems with stagnation and error in reconstruction. The common problem of stagnation may be due to the possibility of the iterative process to approach the object and its twin or shifted image, the support not tight enough or the boundary not sharp enough . Besides the uniqueness issue, phase retrieval is also inherently nonconvex and many researchers have believed the lack of convexity in the Fourier magnitude constraint to be a main, if not the dominant, source of numerical problems with the standard phasing algorithms . While there have been dazzling advances in applications of phase retrieval in the past decades , we still do not know just how much of the error and stagnation problems is attributable to to the lack of uniqueness or convexity.
We propose here to refocus on the issue of uniqueness as uniqueness is undoubtedly the first foundational issue of any inverse problem, including phase retrieval. Specifically we will first establish uniqueness in the absolute sense with random illumination under general, physically reasonable object constraints (Figure 1) and secondly demonstrate that random illumination practically alleviates most numerical problems and drastically improves the quality of reconstruction.
The Fourier transform can be obtained from the -transform as
by some abuse of notation. The discrete phase retrieval problem is to determine from the knowledge of the Fourier magnitude .
The question of uniqueness was partially answered in which says that in dimension two or higher and with the exception of a measure zero set of finite sequences phase retrieval has a unique solution up to the equivalence class of “trivial associates” (i.e. relative uniqueness). These trivial, but omnipresent, ambiguities include constant global phase,
Conjugate inversion produces the so-called twin image.
This landmark uniqueness result, however, does not address the following issues. First, a given object array, there is no way of deciding a priori the irreducibility of the corresponding -transform and the relative uniqueness of the phasing problem. Secondly, although visually no different from the true image the trivial associates (particularly spatial shift and conjugate inversion) nevertheless “confuse” the standard numerical iterative processes and cause serious stagnation .
In this paper, we study the notation of absolute uniqueness: if two finite objects and give rise to the same Fourier magnitude data, then unequivocally. More importantly, we present the approach of random (phase or amplitude) illumination to the absolute uniqueness of phase retrieval. The idea of random illumination is related to coded-aperture imaging whose utility in other imaging contexts than phase retrieval has been established experimentally as well as mathematically .
Our basic tool is an improved version (Theorem 2) of the irreducibility result of with, however, a completely different perspective and important practical implications. The main difference is that while the classical result works with generic (thus random) objects from a certain ensemble Theorem 2 can deal with a given, deterministic object whose support has rank . This improvement is achieved by endowing the probability measure on the ensemble of illuminations, which we can manipulate, instead of the space of objects, which we can not control, as in the classical setting.
On the basis of almost sure irreducibility, the mere assumption that the phases or magnitudes of the object at two arbitrary points lie in a countable set enforces uniqueness, up to a global phase, in phase retrieval with a single random illumination (Theorem 3). The absolute uniqueness can be enforced then by imposing the positivity constraint (Corollary 1). For objects satisfying a tight sector condition, absolute uniqueness is valid with high probability depending on the object sparsity for either phase or amplitude illumination (Theorem 4). For complex-valued objects under a magnitude constraint, uniqueness up to a global phase is valid with high probability (Theorem 5). For general complex-valued objects, almost sure uniqueness, up to global phase, is proved for phasing with two independent illuminations (Theorem 6).
The paper is organized as follows. In Section 2 we discuss various sources of ambiguity. In Section 3 we prove the almost sure irreducibility (Theorem 2 and Appendix). In Section 4 we derive the uniqueness results (Theorem 3, 4, 5, 6 and Corollary 1). We demonstrate phasing with random illumination in Section 5. We conclude in Section 6.
Sources of ambiguity
As commented before the phase retrieval problem does not have a unique solution. Nevertheless, the possible solutions are constrained as stated in the following theorem .
Let the -transform of a finite complex-valued sequence be given by
where are nontrivial irreducible polynomials. Let be the -transform of another finite sequence . Suppose . Then must have the form
where is a subset of .
is the autocorrelation function of . Note the symmetry .
The theorem then follows straightforwardly from the equality between the autocorrelation functions of and , because , and the unique factorization of polynomials (see for more details).
If the finite array is known a priori to vanish outside the lattice , then by Shannon’s sampling theorem for band-limited functions the sampling domain for can be limited to the finite regular grid
since is band-limited to the set .
There are three sources of ambiguity. First, the linear phase term in (1) remain undetermined because the autocorrelation operation destroys information about spatial shift. The unspecified constant phase is another source of ambiguity.
To understand the physical meaning of the operation
which is the -transform of the conjugate space-inversed array . The same is true in multi-dimensions.
The subtlest form of ambiguity is caused by partial conjugate inversion on some, but not all, factors of a factorable object, with a reducible -transform, without which the conjugate inversion, like spatial shift and global phase, is global in nature and considered “trivial” in the literature (even though the twin image may have an opposite orientation).
In this paper, we consider both types, trivial and nontrivial, of ambiguity, as they both can degrade the performance of phasing schemes. Our main purpose is to show by rigorous analysis that with random illumination it is possible to eliminate all ambiguities at once.
Irreducibility
Random illumination amounts to replacing the original object by
where , representing the incident field, is a known array of samples of random variables (r.v.s). The idea is to first modify the object by the encoding array so that phase retrieval has unique solution and then use the prior knowledge of to recover .
Nearly independent random illumination can be produced by a diffuser placed near the object, cf. Figure 1. The illumination field can be randomly modulated in phase only with the use computer generated holograms , random phase plates and liquid crystal phase-only panels . One of the best known amplitude masks is uniformly redundant array and its variants . The advantage of phase mask, compared to amplitude mask, is the lossless energy transmission of an incident wavefront through the mask. By placing either phase or amplitude mask at a distance from the object, one can create an illumination field modulated in both amplitude and phase in a way dependent on the distance .
If the object support does not touch all the coordinate hyperplanes, then the the irreducibility holds true, up to some monomial of . In view of Theorem 1 this is sufficient for our purpose.
The theorem does not hold if the rank-2 condition fails. For example, let be any monomial and consider
which is factorable by, again, the fundamental theorem of algebra.
The proof of Theorem 2 is given in the Appendix.
Theorem 2 improves in several aspects on the classical result that the set of the reducible polynomials has zero measure in the space of multivariate polynomials with real-valued coefficients . The main improvement is that while the classical result works with generic (thus random) objects Theorem 2 deals with any deterministic object with minimum (and necessary) conditions on its support set. By definition, deterministic objects belong to the measure zero set excluded in the classical setting of . It is both theoretically and practically important that Theorem 2 places the probability measure on the ensemble of illuminations, which we can manipulate, instead of the space of objects, which we can not control.
In the next section, we go further to show that with additional, but for all practical purposes sufficiently general, constraints on the values of the object, we can essentially remove all ambiguities with the only possible exception of global phase factor. This decisive step distinguishes our method from the standard approach.
Uniqueness
Without additional a priori knowledge on the object Theorem 2, however, does not preclude the trivial ambiguities such as global phase, spatial shift and conjugate inversion. For example, we can produce another finite array that yields the same measurement data by setting
Suppose the object support has rank . Suppose either of the following cases holds:
Then is determined uniquely, up to a global phase, by the Fourier magnitude measurement on the lattice with probability one.
For the two-point constraint in case (i) to be convex, it is necessary for the constraint set to be a singleton, namely the phases of the object at two nonzero points must take on a single known value. On the other hand, the amplitude constraint in case (ii) can never be convex unless the set is a singleton and the object phases are the same at the two points.
By Theorem 2 the -transform of is irreducible with probability one. We prove the theorem case by case.
Case (i): Suppose the phases of and belong to the coutable set . Let us show the probability that the phase of as given by (9) with takes on a value in at two distinct points is zero.
Since and the phases of are independent, continuous r.v.s on , the phase of is continuously distributed on for all .
Now suppose the phase of for some lies in the set . This implies that must belong to the countable set which is shifted by the negative phase of . The phase of at a different location , however, almost surely does not take on any value in the set for any fixed unless . Since a countable union of measure-zero sets has zero measure, the probability that the phases of at two points lie in is zero if .
Likewise, has a random phase that is continuously distributed on and by the same argument the probability that the phases of as given by (10) at two points lie in is zero.
The global phase , however, can not be determined uniquely in either case.
The global phase factor can be determined uniquely by additional constraint on the values of the object. For example, the following result follows immediately from Theorem 3 (i).
With a real, positive object, the countable set for phase is the singleton and the global phase is uniquely fixed. ∎
2. Sector constraint
More generally, we consider the sector constraint that the phases of belong to . For example, the class of complex-valued objects relevant to -ray diffraction typically have nonnegative real and imaginary parts where the real part is the effective number of electrons coherently diffracting photons, and the imaginary part represents the attenuation . For such objects, .
Generalizing the argument for Theorem 3 we can prove the following.
Suppose the object support has rank . Let the finite object satisfy the sector constraint that the phases of belong to . Let be the sparsity (the number of nonzero elements) of the object.
In both cases, the global phase is uniquely determined if the sector is tight in the sense that no proper interval of contains all the phases of the object.
Case (i): Consider first the expression (9) with any and the independently distributed r.v.s of corresponding to nonoverlapping pairs of points . The probability for every such the phase of to lie in the sector is for any and hence the probability for all with to lie in the sector is at most . The union over of these events has probability at most .
Likewise the probability for all given by (10) to lie in the first quadrant for any is at most .
Case (ii): For (9) with any the independently distributed random variables corresponding to nonoverlapping pairs of points , satisfy the sector constraint with probability at most if . Hence the probability that all with satisfy the sector constraint is at most
For (10) with and any , at and hence lies in the first quadrant with probability one. For , satisfies the sector constraint with probability if . Now the independently distributed r.v.s corresponding to nonoverlapping pairs of points satisfy the sector constraint with probability at most if . Hence the probability that all given by (10) with arbitrary satisfy the sector constraint is at most . ∎
3. Magnitude constraint
Likewise if the object satisfies a magnitude constraint then we can use random amplitude illumination to enforce uniqueness (up to a global phase).
The proof is similar to that for Theorem 4(ii).
For (9) with any the independently distributed random variables corresponding to nonoverlapping pairs of points satisfy with probability less than for any . Hence the probability that with satisfy the magnitude constraint at or more points is at most .
For (10) with any , at and hence satisfies the magnitude constraint with probability one. For , there is at most probability for to satisfy the magnitude constraint. By independence, the independently distributed r.v.s corresponding to nonoverlapping pairs of points satisfy the magnitude constraint with probability at most . Hence the probability that given by (10) with arbitrary satisfy the magnitude constraint at or more points is at most .
The global phase factor is clearly undetermined. ∎
As in Theorem 3 case (ii) the magnitude constraint here, however, is not convex.
4. Complex objects without constraint
For general complex-valued objects without any constraint, we consider two sets of Fourier magnitude data produced with two independent random illuminations and obtain almost sure uniqueness modulo global phase.
Let be a finite complex-valued array whose support has rank . Let and be two independent arrays of r.v.s satisfying the assumptions in Theorem 2.
Then with probability one is uniquely determined, up to a global phase, by the Fourier magnitude measurements on with two illuminations and .
If the second illumination is deterministic and results in an irreducible -transform while is random as above, then the same conclusion holds.
Let be another array that vanishes outside and produces the same data. By Theorem 1, 2 and Remark 1
Four scenarios of ambiguity exist but because of the independence of none can arise.
This almost surely can not occur unless in which case equals up to a global phase factor.
The other possibilities can be similarly ruled out:
for any .
The same argument above applies to the case of deterministic if the resulting -transform is irreducible.
Numerical examples
We test the case of random phase illumination on a real, positive image consisting of the original Cameraman in the middle, surrounded by a black margin (zero padding) of 13 pixels in width (Figure 2(a)). We synthesize and sample the Fourier magnitudes at the Nyquist rate (Remark 1) and implement the standard Error Reduction (ER) and Hybrid-Input-Output (HIO) algorithms in the framework of the oversampling method . By Corollary 1, absolute uniqueness holds with a random phase illumination.
Let and be the Fourier transform and the diagonal matrix representing the illumination. For the uniform illumination . The ER and HIO algorithms are described below.
In the original version of HIO , the hard thresholding is replaced by
where the feedback parameter is used in the simulations. ER has the desirable property that the residual is reduced after each iteration under either uniform or random phase illumination . When absolute uniqueness holds, a vanishing residual then implies a vanishing reconstruction error.
Figure 2 shows the results of ER reconstruction with random phase illumination. The ER iteration converges to the true image after 40 iterations. For HIO reconstruciton we apply 50 ER iterations after 100 HIO iterations as suggested in . HIO has essentially the same performance as ER (Figure 3 (a), (b)). The relative residual curve, Figure 3(d), and the relative error curve, Figure 4, however, indicate a small improvement by HIO. The close proximity between the vanishing residual curve and the vanishing error curve for ER and HIO reflects the absolute uniqueness under random illumination.
With uniform illumination, ER produces a poor result (Figure 5 (a)), resulting a error (Figure 5 (b)) after more than 1000 iterations. The relative change curve, Figure 5(c), indicates stagnation or convergence to a fixed point after 100 iterations and the relative residual plot, Figure 5(d), shows non-convergence to the true image. For HIO reconstruction, we augment it with 50 ER iterations at the end of 1000 HIO iterations. While HIO improves the performance of ER but still leads to a shifted, inverted image which is also severely distorted (Figure 6 (a), (b)). The ripples and stripes in Figure 6 (a) are a well known artifact of HIO reconstruction . As expected, HIO reduces the residual and does not stagnate as much as ER (Figure 6 (c), (d)) but its error is greater than that of ER due to the interferences from shifted and twin images present under the uniform illumination (Figure 7).
To summarize, under a random phase illumination, the problems of stagnation and error disappear and phasing with ER/HIO achieves accurate, high-quality recovery. These experiments confirm our belief that a central barrier to stable and accurate phasing by the standard methods is the lack of absolute uniqueness.
Conclusions
In conclusion, we have proposed random illumination to address the uniqueness problem of phase retrieval. For general random illumination we have proved almost sure irreducibility for any complex-valued object whose support has rank (Theorem 2). We have proved the almost sure uniqueness, up to a global phase, under the two-point assumption (Theorem 3). The absolute uniqueness is then enforced by the positivity constraint (Corollary 1). Under the tight sector constraint, we have proved the absolute uniqueness with probability exponentially close to unity as the object sparsity increases (Theorem 4). Under the magnitude constraint, we have proved uniqueness up to a global phase with probability exponentially close to unity (Theorem 5). For general complex-valued objects without any constraint, we have established almost sure uniqueness modulo global phase with two independent illuminations (Theorem 6).
Numerical experiments reveal that phasing with random illumination drastically reduces the reconstruction error, the number of Fourier magnitude data and removes the stagnation problem commonly associated with the ER and HIO algorithms. Enforcement of absolute uniqueness therefore appears to have a profound effect on the performance of the standard phasing algorithms.
Systematic and detailed study of phasing in the presence of (additive or multiplicative) noise with low-resolution random illuminations and sub-Nyquist sampling rates will be presented in the forthcoming paper .
Appendix A Proof of Theorem 2
Our argument is based on and can be extended to the case of more than two independent variables. For simplicity of notation, we present the proof for the case of two independent variables.
First we state an elementary result from algebraic geometry (see, e.g., , page 65).
If a homogeneous polynomial of (total) degree is irreducible, then is also irreducible with degree .
For a polynomial of degree , the expression
defines a homogeneous polynomial of degree with the property . The process from to is called homogenization while the reverse process is called dehomogenization. Homogenization, in conjunction with Proposition 1, is a useful tool for studying the question of irreducibility.
A subset of a projective space is closed in the Zariski topology if and only if it is an algebraic variety, i.e. the common zero set
of a finite number of homogeneous polynomials of the homogeneous coordinates of the projective space. The Zariski topology is much cruder than the metric topology. Indeed, a Zariski closed set is either the whole space or a measure-zero, nowhere-dense closed set in the metric topology as stated in the following (see, e.g. , page 115).
Any Zariski closed proper subset of a (real or complex) projective variety has measure zero with respect to the standard measure on the projective variety.
Now the multiplication of two polynomials of degree and determines a regular (i.e. polynomial) mapping
in the following way. Let and be homogeneous polynomials of degrees and , respectively. Let and be the coefficients of and , respectively. Then the coefficients of the image point are given by
In other words is bilinear in and and thus is regular. Clearly we have
Let and let the ensemble of polynomials corresponding to be identified with
Following the suggestion in , we now prove
The argument is based on two observations. First the polynomial
Secondly, for any satisfying the assumptions of Theorem 2, there exists a set of three points which can be transformed into , the support of (16), under a rational map.
We separate the analysis of the second observation into two cases.
Case 1: . Then there are at least two other points, say , belonging to . Without loss of generality, we assume . Because has rank 2, .
to . This amounts to a linear transformation from the set of independent vectors
to the set . This transformation can be accomplished by the following matrix
where the divisor is nonzero. To ensure integer entries in (18) we set
Case 2: for some positive integers . Then there is at least another point such that are not collinear, which means .
Suppose . Consider the polynomial
By the same analysis above the form (16) can be achieved by the transformation matrix
which has integer entries if is a multiple of . With the choice
Suppose . Consider the polynomial
The form (16) can be achieved by the transformation matrix
Suppose . Consider the polynomial
The form (16) can be achieved by the transformation matrix
To conclude the proof of Proposition 4, in any above case, if the polynomial is reducible (i.e. has a non-monomial factor), then we can write and
where are non-monomial factors. Let be the lowest (possibly negative) power in of . If , then the factorization (19) implies that has a non-monomial factor. If , then the factorization (19) implies that has a non-monomial factor. Either case contradicts the fact that is irreducible. So is irreducible. The proof of Proposition 4 is complete.
Acknowledgements. I am grateful to my colleagues Greg Kuperberg and Brian Osserman for inspiring discussions on the proof of Theorem 2, an improvement of the earlier version which assumes convexity of the support. I thank my student Wenjing Liao for performing simulations and producing the figures.