OSA's Digital Library

Journal of the Optical Society of America A

Journal of the Optical Society of America A


  • Vol. 19, Iss. 7 — Jul. 1, 2002
  • pp: 1334–1345

Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization

Heinz H. Bauschke, Patrick L. Combettes, and D. Russell Luke  »View Author Affiliations

JOSA A, Vol. 19, Issue 7, pp. 1334-1345 (2002)

View Full Text Article

Acrobat PDF (235 KB)

Browse Journals / Lookup Meetings

Browse by Journal and Year


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools



The phase retrieval problem is of paramount importance in various areas of applied physics and engineering. The state of the art for solving this problem in two dimensions relies heavily on the pioneering work of Gerchberg, Saxton, and Fienup. Despite the widespread use of the algorithms proposed by these three researchers, current mathematical theory cannot explain their remarkable success. Nevertheless, great insight can be gained into the behavior, the shortcomings, and the performance of these algorithms from their possible counterparts in convex optimization theory. An important step in this direction was made two decades ago when the error reduction algorithm was identified as a nonconvex alternating projection algorithm. Our purpose is to formulate the phase retrieval problem with mathematical care and to establish new connections between well-established numerical phase retrieval schemes and classical convex optimization methods. Specifically, it is shown that Fienup’s basic input–output algorithm corresponds to Dykstra’s algorithm and that Fienup’s hybrid input–output algorithm can be viewed as an instance of the Douglas–Rachford algorithm. We provide a theoretical framework to better understand and, potentially, to improve existing phase recovery algorithms.

© 2002 Optical Society of America

OCIS Codes
(100.3010) Image processing : Image reconstruction techniques
(100.3020) Image processing : Image reconstruction-restoration
(100.5070) Image processing : Phase retrieval

Heinz H. Bauschke, Patrick L. Combettes, and D. Russell Luke, "Phase retrieval, error reduction algorithm, and Fienup variants: a view from convex optimization," J. Opt. Soc. Am. A 19, 1334-1345 (2002)

Sort:  Author  |  Year  |  Journal  |  Reset


  1. A. Barty, D. Paganin, and K. Nugent, “Phase retrieval in Lorentz microscopy,” Exp. Methods Phys. Sci. 36, 137–166 (2001).
  2. J.-F. Brun, D. de Sousa Meneses, B. Rousseau, and P. Echegut, “Dispersion relations and phase retrieval in infrared reflection spectra analysis,” Appl. Spectrosc. 55, 774–780 (2001).
  3. J. C. Dainty and J. R. Fienup, “Phase retrieval and image reconstruction for astronomy,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 231–275.
  4. L. D. Marks, W. Sinkler, and E. Landree, “A feasible set approach to the crystallographic phase problem,” Acta Crystallogr. Sect. A 55, 601–612 (1999).
  5. R. P. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Am. A 7, 394–411 (1990).
  6. P. Jaming, “Phase retrieval techniques for radar ambiguity problems,” J. Fourier Anal. Appl. 5, 309–329 (1999).
  7. W. O. Saxton, Computer Techniques for Image Processing in Electron Microscopy (Academic, New York, 1978).
  8. L. S. Taylor, “The phase retrieval problem,” IEEE Trans. Antennas Propag. AP-29, 386–391 (1981).
  9. P. Roman and A. S. Marathay, “Analyticity and phase retrieval,” Nuovo Cimento 30, 1452–1464 (1963).
  10. A. Walther, “The question of phase retrieval in optics,” Opt. Acta 10, 41–49 (1963).
  11. E. Wolf, “Is a complete determination of the energy spectrum of light possible from measurements of degree of coherence?” Proc. Phys. Soc. London 80, 1269–1272 (1962).
  12. Lord Rayleigh (J. W. Strutt), “On the interference bands of approximately homogeneous light; in a letter to Prof. A. Michelson,” Philos. Mag. 34, 407–411 (1892).
  13. R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik (Stuttgart) 35, 237–246 (1972).
  14. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982).
  15. D. C. Youla and H. Webb, “Image restoration by the method of convex projections: Part I—theory,” IEEE Trans. Med. Imaging MI-1, 81–94 (1982).
  16. A. Lent and H. Tuy, “An iterative method for the extrapolation of band-limited functions,” J. Math. Anal. Appl. 83, 554–565 (1981).
  17. A. Levi and H. Stark, “Signal reconstruction from phase by projection onto convex sets,” J. Opt. Soc. Am. 73, 810–822 (1983).
  18. H. J. Trussell and M. R. Civanlar, “The feasible solution in signal restoration,” IEEE Trans. Acoust. Speech Signal Process. ASP-32, 201–212 (1984).
  19. L. M. Brègman, “The method of successive projection for finding a common point of convex sets,” Sov. Math. Dokl. 6, 688–692 (1965).
  20. P. L. Combettes, “The foundations of set theoretic estimation,” Proc. IEEE 81, 182–208 (1993).
  21. P. L. Combettes, “The convex feasibility problem in image recovery,” in Advances in Imaging and Electron Physics, P. W. Hawkes, ed. (Academic, Orlando, Fla., 1996), Vol. 95, pp. 155–270.
  22. H. Stark, ed., Image Recovery: Theory and Application (Academic, Orlando, Fla., 1987).
  23. H. Stark and Y. Yang, Vector Space Projections: A Numerical Approach to Signal and Image Processing, Neural Nets, and Optics (Wiley, New York, 1998).
  24. P. L. Combettes and H. J. Trussell, “Stability of the linear prediction filter: a set theoretic approach,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (Institute of Electrical and Electronics Engineers, New York, 1988), pp. 2288–2291.
  25. M. Hedley, H. Yan, and D. Rosenfeld, “Motion artifact correction in MRI using generalized projections,” IEEE Trans. Med. Imaging 10, 40–46 (1991).
  26. R. G. Hoptroff, P. W. McOwan, T. J. Hall, W. J. Hossak, and R. E. Burge, “Two optimization approaches to cohoe design,” Opt. Commun. 73, 188–194 (1989).
  27. B. E. A. Saleh and K. M. Nashold, “Image construction: optimum amplitude and phase masks in photolithography,” Appl. Opt. 24, 1432–1437 (1985).
  28. A. Levi and H. Stark, “Image restoration by the method of generalized projections with application to restoration from magnitude,” J. Opt. Soc. Am. A 1, 932–943 (1984).
  29. A. Levi and H. Stark, “Restoration from phase and magnitude by generalized projections,” in Image Recovery: Theory and Applications, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp 277–320.
  30. J. A. Cadzow, “Signal enhancement—a composite property mapping algorithm,” IEEE Trans. Acoust. Speech Signal Process. 36, 49–62 (1988).
  31. P. L. Combettes and H. J. Trussell, “Method of successive projections for finding a common point of sets in metric spaces,” J. Optim. Theory Appl. 67, 487–507 (1990).
  32. N. E. Hurt, “Signal enhancement and the method of successive projections,” Acta Appl. Math. 23, 145–162 (1991).
  33. S. Chrétien and P. Bondon, “Cyclic projection methods on a class of nonconvex sets,” Numer. Funct. Anal. Optim. 17, 37–56 (1996).
  34. R. Barakat and G. Newsam, “Algorithms for reconstruction of partially known, band-limited Fourier-transform pairs from noisy data,” J. Opt. Soc. Am. A 2, 2027–2039 (1985).
  35. R. Barakat and G. Newsam, “Algorithms for reconstruction of partially known, band-limited Fourier-transform pairs from noisy data. II. The nonlinear problem of phase retrieval,” J. Integral Eq. 9, 77–125 (1985).
  36. The issue of nonuniqueness of the projection is not to be confused with the uniqueness of solutions to the phase problem. The results surveyed, for instance, in Ref. 37, are not affected by the multivaluedness of the projection operators.
  37. M. H. Hayes, “The unique reconstruction of multidimensional sequences from Fourier transform magnitude or phase,” in Image Recovery: Theory and Application, H. Stark, ed. (Academic, Orlando, Fla., 1987), pp. 195–230.
  38. D. R. Luke, “Analysis of wavefront reconstruction and deconvolution in adaptive optics,” Ph.D. thesis (University of Washington, Seattle, Wash., June 2001), ftp://amath.washington.edu/pub/russell/Dissertation.ps.gz.
  39. D. R. Luke, J. V. Burke, and R. Lyon, “Optical wavefront reconstruction: theory and numerical methods,” SIAM Rev. 44 (2) (2002), ftp://amath.washington.edu/pub/russell/Luke_Burke_Lyon_01.ps.gz.
  40. J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. 3, 27–29 (1978).
  41. J. R. Fienup, “Space object imaging through the turbulent atmosphere,” Opt. Eng. 18, 529–534 (1979).
  42. J. R. Fienup, “Iterative method applied to image reconstruction and to computer-generated holograms,” Opt. Eng. 19, 297–305 (1980).
  43. J. W. Goodman, Statistical Optics (Wiley Interscience, New York, 1985).
  44. “a.e.” stands for “almost everywhere” in the sense of measure theory, since, strictly speaking, the elements of L are classes of equivalence of signals that may differ on a set of zero measure. For technical details on L, see, for instance, Ref. 45.
  45. C. D. Aliprantis and K. C. Border, Infinite-Dimensional Analysis, 2nd ed. (Springer-Verlag, Berlin, 1999).
  46. J. N. Cederquist, J. R. Fienup, C. C. Wackerman, S. R. Robinson, and D. Kryskowski, “Wave-front phase estimation from Fourier intensity measurements,” J. Opt. Soc. Am. A 6, 1020–1026 (1989).
  47. P. L. Combettes, “Inconsistent signal feasibility problems: least-squares solutions in a product space,” IEEE Trans. Signal Process. 42, 2955–2966 (1994).
  48. G. T. Herman, Image Reconstruction from Projections—The Fundamentals of Computerized Tomography (Academic, New York, 1980).
  49. J. H. Seldin and J. R. Fienup, “Numerical investigation of the uniqueness of phase retrieval,” J. Opt. Soc. Am. A 7, 412–427 (1990).
  50. H. H. Bauschke, “The composition of finitely many projections onto closed convex sets in Hilbert space is asymptotically regular,” Proc. Am. Math. Soc. (to be published).
  51. H. H. Bauschke and J. M. Borwein, “On the convergence of von Neumann’s alternating projection algorithm for two sets,” Set-Valued Anal. 1, 185–212 (1993).
  52. H. H. Bauschke, J. M. Borwein, and A. S. Lewis, “The method of cyclic projections for closed convex sets in Hilbert space,” in Recent Developments in Optimization Theory and Nonlinear Analysis (American Mathematical Society, Providence, R.I., 1997), pp. 1–38.
  53. P. L. Combettes and P. Bondon, “Hard-constrained inconsistent signal feasibility problems,” IEEE Trans. Signal Process. 47, 2460–2468 (1999).
  54. L. G. Gubin, B. T. Polyak, and E. V. Raik, “The method of projections for finding the common point of convex sets,” USSR Comput. Math. Math. Phys. 7, 1–24 (1967).
  55. D. C. Youla and V. Velasco, “Extensions of a result on the synthesis of signals in the presence of inconsistent constraints,” IEEE Trans. Circuits Syst. CAS-33, 465–468 (1986).
  56. Iterative signal recovery projection algorithms have also been implemented optically without sampling the continuous waveforms (e.g., Ref. 57). In such instances, the underlying signal space is L itself.
  57. R. J. Marks II, “Coherent optical extrapolation of 2-D band-limited signals: processor theory,” Appl. Opt. 19, 1670–1672 (1980).
  58. Let R be a set of real numbers. If R≠⌀, then inf (R) stands for the infimum of R, i.e., the largest number in [−∞, +∞] that is smaller than or equal to all elements of R. By convention, inf(⌀)=+∞.
  59. J.-P. Aubin and H. Frankowska, Set-Valued Analysis (Birkhäuser, Boston, Mass., 1990).
  60. For theoretical reasons, the sets (and functions) that we deal with must be “measurable”—this is not the same as being “physically measurable” or “observable”! For our purposes, measurable sets and functions constitute a sufficiently large class to work with; thus all closed and open subsets (and all continuous functions) are measurable, as well as various combinations of those.
  61. Mathematically, this set is assumed to have nonzero measure.
  62. The complex Hilbert space L from the phase retrieval problem is also a real Hilbert space, provided that we use the real part of the inner product as the new inner product.
  63. C. Tisseron, Notions de Topologie—Introduction aux Espaces Fonctionnels (Hermann, Paris, 1985).
  64. Recall the notation from Remark 3.4.
  65. F. Deutsch, Best Approximation in Inner Product Spaces (Springer-Verlag, New York, 2001).
  66. In our setting, this means that the set {t∈R N :A(t)⋂Z≠ ⌀} is measurable for every closed (or, equivalently, open) set Z in C ; see Section 8.1 in Ref. 59 and Section 14.A in Ref. 67.
  67. R. T. Rockafellar and R. J.-B. Wets, Variational Analysis (Springer-Verlag, Berlin, 1998).
  68. If x∈L, then Re (x) denotes the function defined by Re (x): t ↦ Re [x(t)], and [Re (x)]+ subsequently denotes the positive part: [Re (x)]+: t ↦ max {0, Re [x(t)]}.
  69. R. W. Gerchberg, “Super-resolution through error energy reduction,” Opt. Acta 21, 709–720 (1974).
  70. D. C. Youla, “Generalized image restoration by the method of alternating orthogonal projections,” IEEE Trans. Circuits Syst. CAS-25, 694–702 (1978).
  71. A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst. CAS-22, 735–742 (1975).
  72. A. E. Çetin and R. Ansari, “Convolution-based framework for signal recovery and applications,” J. Opt. Soc. Am. A 5, 1193–1200 (1988).
  73. K. Goebel and W. A. Kirk, Topics in Metric Fixed Point Theory (Cambridge U. Press, Cambridge, UK, 1990).
  74. Z. Opial, “Weak convergence of the sequence of successive approximations for nonexpansive mappings,” Bull. Am. Math. Soc. 73, 591–597 (1967).
  75. A. Genel and J. Lindenstrauss, “An example concerning fixed points,” Isr. J. Math. 22, 81–86 (1975).
  76. Alternative descriptions of these algorithms have been proposed; see, for instance, Ref. 77.
  77. G. T. Herman and D.-W. Ro, “Image recovery using iterative data refinement with relaxation,” Opt. Eng. 29, 513–523 (1990).
  78. Another formulation of the algorithm also includes a nonnegativity constraint.14
  79. The reason for this difference is that Fienup, on page 2763 of Ref. 14, defines γ as the set of all points where (in our notation) PM (xn) violates the object domain constraints. Hence γ={t∈∁D: (PM (xn)) (t)≠0}, or t∈γ if and only if t∈D and PM (xn)(t)≠0. It follows that t belongs to the complement of γ if and only if t∈D or PM (xn)(t)=0. The latter condition then leads to this different interpretation of the HIO algorithm. Sticking with this interpretation for another moment, we could set D(n)=D∪{t∈ ∁D: PM (xn)(t)=0} and S(n)={z∈L: z ⋅ 1∁D(n) = 0} and obtain analogously xn+1 =(PS(n) (2PM −I)+(I−PM)) (xn). In practical experiments for problem (5), however, this ambiguity has hardly an impact, as the sets γ and ∁D almost always coincide.
  80. H. Takajo, T. Shizuma, T. Takahashi, and S. Takahata, “Reconstruction of an object from its noisy Fourier modulus: ideal estimate of the object to be constructed and a method that attempts to find that object,” Appl. Opt. 38, 5568–5576 (1999).
  81. H. Takajo, T. Takahashi, and T. Shizuma, “Further study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A 16, 2163–2168 (1999).
  82. H. Takajo, T. Takahashi, R. Ueda, and M. Taninaka, “Study on the convergence property of the hybrid input–output algorithm used for phase retrieval,” J. Opt. Soc. Am. A 15, 2849–2861 (1998).
  83. The corresponding mask is certainly much easier to implement.
  84. J. R. Fienup and C. C. Wackerman, “Phase retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3, 1897–1907 (1986).
  85. The algorithms discussed here for solving problem (25) can be viewed in the broader context of finding a zero of the sum of two maximal monotone operators. Good starting points are Refs. 86 87 88.
  86. P. L. Combettes, “Fejér-monotonicity in convex optimization,” in Encyclopedia of Optimization, C. A. Floudas and P. M. Pardalos, eds. (Kluwer, New York, 2001), Vol. 2, pp. 106–114.
  87. J. Eckstein, “Splitting methods for monotone operators with applications to parallel optimization,” Ph.D. thesis (Department of Civil Engineering, Massachusetts Institute of Technology, Cambridge, Mass., 1989), available as Rep. LIDS-TH-1877 (Laboratory for Information and Decision Sciences, MIT).
  88. J. Eckstein and D. P. Bertsekas, “On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators,” Math. Program. (Ser. A) 55, 293–318 (1992).
  89. Unfortunately, PA PB is generally not firmly nonexpansive. However, it is strongly nonexpansive (see Ref. 90 for a precise definition and further information), and, for this class of mappings, a result corresponding to Fact 3.20 does exist.90
  90. R. E. Bruck and S. Reich, “Nonexpansive projections and resolvents of accretive operators in Banach spaces,” Houston J. Math. 3, 459–470 (1977).
  91. H. H. Bauschke and J. M. Borwein, “On projection algorithms for solving convex feasibility problems,” SIAM Rev. 38, 367–426 (1996).
  92. H. H. Bauschke, “Projection algorithms: results and open problems,” in Inherently Parallel Algorithms in Feasibility and Optimization and Their Applications, Studies in Computational Mathematics, D. Butnariu, Y. Censor, and S. Reich, eds. (Elsevier, Amsterdam, 2001), Vol. 8, pp. 11–22.
  93. P. L. Combettes, “Quasi-Fejérian analysis of some optimization algorithms,” in Inherently Parallel Algorithms in Feasibility and Optimization and Their Applications, Studies in Computational Mathematics, D. Butnariu, Y. Censor, and S. Reich, eds. (Elsevier, Amsterdam, 2001), Vol. 8, pp. 115–152.
  94. R. L. Dykstra, “An algorithm for restricted least squares regression,” J. Am. Stat. Assoc. 78, 837–842 (1983).
  95. J. P. Boyle and R. L. Dykstra, “A method for finding projections onto the intersection of convex sets in Hilbert spaces,” in Advances in Order Restricted Statistical Inference (Springer, Berlin, 1986), pp. 28–47.
  96. In the aforementioned context of maximal monotone operators, 85 Dykstra’s algorithm can be interpreted as a tight version of the Peaceman–Rachford algorithm. See page 77 in Ref. 87 for further information. Let us also note that in the standard linear case, the Peaceman–Rachford and Douglas–Rachford algorithms can be viewed from a unifying standpoint (see Section 7.4 in Ref. 97).
  97. R. S. Varga, Matrix Iterative Analysis, 2nd ed. (Springer-Verlag, New York, 2000).
  98. H. H. Bauschke and J. M. Borwein, “Dykstra’s alternating projection algorithm for two sets,” J. Approx. Theory 79, 418–443 (1994).
  99. H. H. Bauschke and A. S. Lewis, “Dykstra’s algorithm with Bregman projections: a convergence proof,” Optimization 48, 409–427 (2000).
  100. N. Gaffke and R. Mathar, “A cyclic projection algorithm via duality,” Metrika 36, 29–54 (1989).
  101. S.-P. Han, “A successive projection method,” Math. Program. (Ser. A) 40, 1–14 (1988).
  102. H. Hundal and F. Deutsch, “Two generalizations of Dykstra’s cyclic projections algorithm,” Math. Program. (Ser. A) 77, 335–355 (1997).
  103. P. L. Combettes, “Signal recovery by best feasible approximation,” IEEE Trans. Image Process. 2, 269–271 (1993).
  104. The Douglas–Rachford algorithm was originally developed as a linear implicit iterative method to solve partial differential equations in Ref. 105 (see also Chap. 7 in Ref. 97). It was extended to an operator splitting method for finding a zero of the sum of two maximal monotone operators by Lions and Mercier in Ref. 106. When it is applied to the normal cone maps of the constraint sets, one obtains a method for solving problem (25). See Refs. 86 87 88 and 106 for further information.
  105. J. Douglas and H. H. Rachford, “On the numerical solution of heat conduction problems in two or three space variables,” Trans. Am. Math. Soc. 82, 421–439 (1956).
  106. P.-L. Lions and B. Mercier, “Splitting algorithms for the sum of two nonlinear operators,” SIAM (Soc. Ind. Appl. Math.) J. Numer. Anal. 16, 964–979 (1979).
  107. u is a weak cluster point of a sequence (un) if there exists a subsequence (ukn) such that ukn ⇀u.
  108. If we had used the literal update rule for the HIO algorithm, the present observation would change only in one respect: the set A would be replaced with S(n) (see Remark 4.1 and Ref. 79) and hence vary with n.
  109. E. H. Zarantonello, “Projections on convex sets in Hilbert space and spectral theory,” in Contributions to Nonlinear Functional Analysis, E. H. Zarantonello, ed. (Academic, New York, 1971), pp. 237–424.
  110. L. Debnath and P. Mikusiński, Introduction to Hilbert Spaces with Applications, 2nd ed. (Academic, San Diego, Calif., 1999).
  111. However, as shown in Property 4.1 in Ref. 39, the set M is not weakly closed; i.e., if a sequence (xn) of points in M converges weakly to a point x, then x may not be in M.
  112. While PB is nonexpansive and therefore Lipschitz continuous, this property is not sufficient to draw the conclusion advertised in Corollary 6.1 in Ref. 88, namely (in our context), that (PB xn) converges weakly to a point in A⋂B. Such a conclusion requires additional assumptions, e.g., that PB be weakly continuous (if so, then PB xn ⇀PB x), as is the case when dim H <+∞ (or when B is a closed affine subspace). Note, however, that the projector onto a closed convex set may fail to be weakly continuous. An example is on page 245 in Ref. 109.

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited