{'fraction_non_alphanumeric': 0.032978069276491086, 'fraction_numerical': 0.005943840951014552, 'mean_word_length': 4.912627241880756, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Phase Retrieval: From Computational Imaging to Machine Learning

Phase retrieval consists in the recovery of a complex-valued signal from intensity-only measurements. As it pervades a broad variety of applications, many researchers have striven to develop phase-retrieval algorithms. Classical approaches involve techniques as varied as generic gradient-descent routines or specialized spectral methods, to name a few. Yet, the phase-recovery problem remains a challenge to this day. Recently, however, advances in machine learning have revitalized the study of phase retrieval in two ways: significant theoretical advances have emerged from the analogy between phase retrieval and single-layer neural networks; practical breakthroughs have been obtained thanks to deep-learning regularization. In this tutorial, we review phase retrieval under a unifying framework that encompasses classical and machine-learning methods. We focus on three key elements: applications, overview of recent reconstruction algorithms, and the latest theoretical results.


I. INTRODUCTION
Phase retrieval is a longstanding computational problem that is simple to define yet difficult to solve. It consists in the d-dimensional search for x * ∈ C d such that There, one assumes that the measurements y ∈ R n and the matrix A ∈ C n×d are known, while |·| is the elementwise modulus operator (Figure 1(a)). Compared to the linear equation y = Ax * , the missing phase in (1) makes the computational problem more challenging. Phase retrieval is a problem that is encountered in many fields. In the context of imaging with electromagnetic waves, it arises whenever one wants to image an object from intensity-only measurements (Figure 1(b)). It has been described in the physics literature since the 1950s, first in crystallography [1], followed by intense studies in astronomy since the '80s [2]. More recently, it has been intimately linked with the next generation of image-reconstruction algorithms which solve this computational problem to push the limits of conventional imaging [3,4,5,6]. Today, phase retrieval encompasses vastly different imaging settings, from X-ray to optical or THz imaging, to observe objects from molecular to astronomical scales. This equation is also encountered outside of computational imaging, for example in computergenerated holography [7] or optical computing [8], to name a few. Despite the wide range of applications, the underlying physics-based models lead to the same universal phase-retrieval problem.
The phase-retrieval problem has stimulated the development of a variety of novel imaging modalities. Likewise, a plethora of reconstruction algorithms have been designed [9], from the first alternating-projection algorithms to convex relaxations [10] and spectral methods [11,12]. Given the rich palette of new reconstruction strategies, it may be unclear for practitioners which one to choose for a particular application. A description of their performance and range of applicability would help in this regard.
Phase retrieval raises theoretical questions such as the uniqueness of the solution [13] and the performance guarantees of algorithms [12,14]. For instance, strong results have been recently obtained for random models where the elements of A are drawn as i.i.d. random variables. In this case, a clear picture has emerged on regions of solvability [15,16] and the best algorithm to use. While recent theoretical progress has pushed our understanding to more structured matrices A [16], tackling other models related to real-world settings remains a challenge.
Phase retrieval is also deeply linked with machine learning. The underlying forward model can be interpreted as a single-layer neural network with quadratic activation function (Figure 1(c)). Thus, the solution to the phase-retrieval problem amounts to the learning of the network weights, which is a typical machine-learning optimization problem. While phase retrieval certainly does not expose one to the full complexity of deep learning, it is already rich enough to provide an entry point to the difficult questions debated in the machine-learning community today. In particular, it exhibits local minima, strongly relies on the quality of the initialization, and behaves differently depending on the amount of training data. As a result, it has sparked a strong interest in the theoretical machine-learning community. In particular, phase retrieval could help us understand better the nonconvex optimization procedures typical of real-world machine learning. Ultimately, phase retrieval represents a transdisciplinary topic that impacts fields ranging from imaging applications to fundamental mathematics.
In this tutorial, we present an overview of phase retrieval today, from imaging applications to algorithms and theory. We introduce a unifying framework for the various problems in which a phase-retrieval equation arises: diverse and scattered applications can actually be grouped into four main categories depending on the matrix A. We also present the latest algorithmic and theoretical advances, in particular, the recent results developed by the theoretical machine-learning community during the past few years. The intuition and principles behind these concepts are described without diving into the technical points. These are described in previous works or dedicated reviews we reference for the interested reader. We hope that our unified overview will foster interaction between research fields and lead to a better understanding of the phase-retrieval equation.

A. Mathematical model
The phase-retrieval equation relates three mathematical objects. In the context of imaging, one wants to reconstruct an object of interest x * from the knowledge of the captured intensities y and the imaging system modeled by the linear operator A. Beyond the field of imaging, this equation often arises over the estimation of a physical quantity x * from energy measurements on a sensor, as detectors are often sensitive to the squared modulus of a field. It will be useful to rewrite this equation as where A is the vertical concatenation of the rows a H i (H denoting Hermitian conjugation), and x, y = x H y is the complex inner product. Thus, each measurement is obtained from the squared modulus of a scalar product with a sampling vector a i ∈ C d . The phase-retrieval problem can then be restated as follows: How many phaseless scalar products are required to reconstruct the unknown vector x * ?
In order to make algorithms robust to experimental measurements, the model can be enriched by incorporating a general noise channel p out such that y i ∼ p out · |a H i x * | 2 , which typically is additive Gaussian or Poissonian. The noise statistics can indeed be leveraged in maximum-likelihood and Bayesian algorithms.
An important quantity to determine the difficulty of the problem is the oversampling ratio α = n/d. It compares the number of measurements to the number of parameters we want to estimate. The greater the oversampling ratio α, the easier the computational problem becomes. This sampling ratio also plays a role in the simpler problem of linear regression y = Ax * . There, the solution is not unique-the set of linear equations is underdetermined-for α < 1 and can be solved exactly for α ≥ 1, in the absence of noise and assuming A to be full-rank. Similar considerations apply to phase retrieval, although with a more complicated behavior. In this tutorial, we define and present the latest results regarding injectivity, weak recovery, full recovery, and asymptotic algorithmic performance. These discussions are highlighted in standalone theory boxes.
Let us emphasize that phase retrieval is, at first sight, a non-convex problem. Indeed, if some x = 0 is in the solution set S of the phase-retrieval equation (1), then the whole set {e iθ x, θ ∈ [0, 2π)} is included in S, while the middle point 1 2 (x + (−x)) = 0 is not in S. This global phase ambiguity can be removed for instance by a quotient of the solution space, or by considering instead the estimation of invariant quantities such as x * (x * ) H . Notwithstanding, the solution set might still be non-convex. Thus, the general phase-retrieval problem cannot be analyzed directly through the prism of convex optimization.

B. Neural-network representation
As (1) describes a linear operation followed by an element-wise nonlinearity, phase retrieval is equivalent to the training of a single-layer neural network with a quadratic activation function. More precisely, for each i = 1, . . . , n, the input data a i is presented to a network with weights x, with a target output y, as depicted in Figure 1(c). Training the network for this regression task amounts to the recovery of x * in the phase-retrieval equation.
This corresponds to the simplest neural network architecture with a nontrivial polynomial nonlinearity. Still, it captures the difficulties of training a shallow neural network with a non-convex loss functional. In particular, convergence is not assured with gradient-descent methods as they may get trapped into local minima or saddle points. More generally, single-layer neural networks belong to the class of generalized linear models. This simple yet rich architecture has attracted the attention of theoreticians, which explains the recent profusion of studies around phase retrieval. This parallel enriches the discussion between physics and machine learning, as we shall further illustrate throughout this tutorial.

A. Unifying Framework for Phase-Retrieval Applications
The phase-retrieval problem is central to a large variety of applications in science and technology. Although many of these approaches were developed independently over the years, they share the same phase-retrieval equation (1). Diversity comes notably from the set of computations performed by the matrix A that depends on the physical model of the imaging system. This domain knowledge sometimes makes it difficult for non-experts to navigate the landscape of phase-retrieval applications.
To offer a unifying perspective, we are proposing to classify phase retrieval into four major groups: Fourier, coded-illumination, coded-detection, and random models (Figure 1(d)). Targeted toward imaging, this classification has been designed to facilitate interaction between fields. It may be used by practitioners to draw links between applications and benefit from the expertise of other fields. For theoreticians, it may give directions to follow when defining new theoretical models to study.

B. Fourier Phase Retrieval
Fourier phase retrieval has been historically the first occurrence of this nonlinear equation. Wave propagation over a long distance, referred to as far-field propagation, is modeled by a Fourier transform [17], as in A = F, where F is a 2D Fourier-transform matrix. The first phase-retrieval problems appeared in modalities where optical components are difficult to manufacture, most notably for X-rays or electron imaging [3]. In these cases, the imaging system is often minimalist and only involves far-field propagation. From the determination of molecular structures using X-rays to aberration correction in astronomy [2], Fourier phase retrieval has been associated with major scientific discoveries in the past century.
Since the Fourier transform is a unitary operation (i.e. F H = F −1 ), we necessarily have that n = d, leading to the oversampling ratio α = 1. This problem is fundamentally ill-posed as we only know amplitudes with no phase information about the Fourier transform of x * . Additional information on x * is required to lift this degeneracy. The famous work of Gerchberg and Saxton [18] uses another intensity measurement of x * directly so that, in (1), A becomes a vertical concatenation of the identity and Fourier operators. Another common assumption exploits a support and non-negativity constraint on x * [2], which is often valid in astronomy. Other more involved assumptions may be invoked, and they will be described in Section V dedicated on regularization. New applications of Fourier phase retrieval have recently appeared, such as speckle correlation imaging for nonline-of-sight imaging or complex media imaging [19,20]. Ultimately, while Fourier phase retrieval has been the first physical model of phase retrieval, it remains a challenging computational problem.

C. Coded-illumination Phase Retrieval
More advanced imaging techniques were developed to overcome the overly restrictive assumptions of Fourier phase retrieval. A solution to lift the phase ambiguity consists in shining different illuminations onto the object of interest in order to capture multiple images after far-field propagation. Each image corresponds to a slightly different view of the same object. This provides the algorithm with complementary and partially redundant information and improves its convergence, at the expense of a more involved experimental setup and longer measurement time. This is the general principle underlying coded-illumination phase retrieval. Here, the linear operator A is a vertical concatenation of operators A l : A l = FD l for l = 1, . . . , L, where D l is a diagonal matrix that corresponds to the coded illumination, n = Ld, and the oversampling ratio α is equal to the number L of images. A variety of imaging applications based on coded illuminations are listed below.
• Coded diffraction imaging relies on a modulation near the object plane [21]. This modality allows ample flexibility on the choice and implementation of the modulation scheme, which explains its application in many different imaging settings. Gratings, masks, as well as tilted wavefronts all represent solutions for structured illumination. • Ptychography is arguably the most successful implementation of coded-illumination phase retrieval [22]. It records the diffraction pattern of the object while it is scanned by an illumination probe with finite support. Oversampling is linked with the overlap between two consecutive illuminated areas-the higher the overlap, the easier it is to solve the phaseretrieval problem. It has been successfully used with a broad range of electromagnetic waves, from Xrays to the THz regime, as well as for electron microscopy. • Fourier ptychography is a variant of ptychography in Fourier space where a sequence of tilted plane waves are sent on the sample [4,5]. In this case, the object x * is defined as the Fourier transform of the object transmission function, and the diagonal matrix D i corresponds to a finite-support probe that is translated in Fourier space. This technique is used in optical microscopy to achieve high-resolution on a large field-of-view.
Imaging scientists have developed specialized algorithms for each of these applications. For example, the Ptychographic Iterative Engine [22] is a popular projection algorithm, while gradient-based methods are also routinely used. However, the structure present in these sensing matrices prevents a thorough establishment of theoretical results for now.

D. Coded-detection Phase Retrieval
Acting on the illumination in the proximity of the object may not always be feasible, for instance in astronomical observations. Nevertheless, multiple images can be obtained upon modulation on the detection side, leading to coded-detection phase retrieval. This can be done by taking several shots at different defocused positions-a technique derived again from astronomy and called phase diversity [23], but applied to other imaging fields as well, such as lensless imaging [6].
There are many ways to describe defocus via freespace wave propagation, and all rely on various degrees of approximation of Maxwell's equations. Using the angular-spectrum method, the operator A is now the vertical concatenation of where u denotes normalized coordinates in Fourier space and the unitless curvature z l is proportional to the defocus distance. Interestingly, this class of phase retrieval also appears outside of imaging for computer-generated holography [7]. There, one wants to produce a 3D hologram-a volumetric intensity pattern. This target defines a set of equations at multiple planes, such that a phase mask is optimized to satisfy these constraints and applied with a wavefront-shaping device.
More general coded-detection strategies may be designed to simplify the phase recovery by using the diversity of the diagonal matrices D l in Fourier space. For example, one can add higher-order optical aberrations beyond defocus or a physical aperture. Once again, a thorough theoretical understanding of this phase-retrieval class is still lacking. It would especially benefit from more results regarding solvability and algorithm choice.

E. Random Phase Retrieval
In random phase retrieval, A is a random matrix. A common setting is the Gaussian model, in which each element of A is drawn from an i.i.d. complex Gaussian distribution. Such an operator A may be seen as a generic linear operation: each sampling vector is independent from the others and adds a little bit of information on the unknown x * . Contrary to expectations, random phase retrieval arises in practical applications too, such as imaging in complex media [24]. There, the random linear operator comes from multiple light scattering or random input patterns synthesized by a programmable device.
Distinct from the previous models, randomness greatly facilitates theoretical studies and allows one to derive asymptotic results on recovery guarantees and to characterize the typical behavior of algorithms. These advances are recent and will be described in more details below. The bottomline is that phase retrieval in the random setting is quite well-understood with active research still ongoing.
Some imaging applications exhibit pseudo-random models, such as coded diffraction imaging with random masks or ptychography with a random probe [22]. They emulate random models, but residual correlations are present between rows of A. This variable amount of randomness makes it challenging to derive precise theoretical guarantees on performance but, to some extent, motivates the use of the latest algorithms designed for the fully random setting.
These observations motivate further work to bridge the gap between experiments and theory. Indeed, to illustrate current efforts, a few theoretical results have recently been generalized beyond the i.i.d. random setting [16]. We can distinguish between two classes of random models, for which strong theoretical results have been reported ∼ P , with P a complex (centered) probability distribution with finite moments of all order (for simplicity). R.2 (right unitarily-invariant) A is a random matrix with probability distribution P (A), and P (A) = P (AU) for every unitary matrix U. Equivalently, the right eigenvectors of A are completely delocalized. As a final remark, the observed robustness of random phase retrieval may inspire new ideas to introduce randomness in experimental settings.

IV. RECONSTRUCTION ALGORITHMS
A large number of reconstruction algorithms have been proposed to tackle the diversity of phase-retrieval models. We present them here together with some intuitive explanation and indications on when to use them in practice. The purpose of this overview is to illustrate the inspiring creativity that researchers deployed while designing phase retrieval algorithms. Dedicated surveys cover more technical descriptions of the recent developments [9,27].

A. Projection Algorithms
Projection algorithms were the first strategy proposed to solve the phase-retrieval problem. They provide a direct and intuitive way to exploit prior information or several intensity measurements. Thanks to their ease of implementation and flexibility, projection algorithms are still widely used today.
More precisely, each intensity image or prior knowledge corresponds to a constraint labelled from 1 to N . The earliest implementation is the Gerchberg-Saxton algorithm [18], which uses the intensity at two measurement planes as constraints in the object and Fourier planes. The error-reduction algorithm resorts to a modulus constraint in the Fourier plane and a support constraint in the object plane [2]. Other techniques such as ptychography [28,22], Fourier ptychography [4,5], or phase diversity perform projections onto more than two sets to exploit even more information for to recover the missing phase.
Injectivity: When is the solution unique? In order to build a mathematically sound theory of phase retrieval, it is important to understand if the solution to the (noiseless) phase-retrieval problem is unique. Stated differently, for a given measurement matrix A ∈ C n×d , is the map A A : x ∈ C d → |Ax| 2 ∈ R n injective a ? Since 2014, the following picture has gradually emerged on a complete characterization of injectivity. Conjecture 1: For any d ≥ 2, and any full-rank A ∈ C n×d the following holds First described as the "(4d − 4) conjecture" [13], point (ii) was later proven in [25]. Meanwhile, point (i) has been proven for small dimensions d ≤ 3 in [13], disproved for d = 4 and refined for large dimensions in a probabilistic sense c [26].
In particular, for large d 1, the threshold n/d 4 should indeed distinguish between "typically injective" and "typically non-injective" behaviors. In practice, in order to efficiently recover in the not-too-oversampled regime n/d 4, this suggests that one should leverage a priori information on the object (e.g. regularization methods).
a While we do not mention it explicitly, the reader should remember that injectivity is always meant up to a global phase. b We refer to the description of the conjecture in [13] for a more detailed definition of the notion of "generic" A. c Let n = 4d − 5. Let im(A) be drawn at random from the set of d-dimensional subspaces of C n , and denote p d the probability that A A is injective. Then: (i) for all d ≥ 2, one has p d < 1, and (ii) lim d→∞ p d = 0.
Each constraint is satisfied by a set of vectors x ∈ C d , and one is aiming at finding the intersection of such sets. One thus projects on all the sets sequentially, to hopefully converge to a suitable solution. The initial estimate isx (0) and the iterations proceed according tô where Π l denotes the projection operators on the solution set of constraint l = 1, . . . , N . For example, we can write the projector of a modulus constraint for an intensity image y l = |A l x * | 2 that corresponds to an invertible linear operator A l as This operator effectively replaces the square modulus of A l u by the one obtained from the measurements y l . Similarly, the projector of a support constraint γ sets all the values outside the support to zero: Π(u) = Diag(1 γ )u, where 1 γ denotes the indicator function of γ. While this class of algorithms is widely used in practice, it unfortunately does not necessarily converge to the set of solutions. More advanced algorithms have been suggested for better convergence, such as the hybrid input-output algorithm [2] that uses relaxed projections for Fourier measurements with support constraint. Another strategy replaces the sequential projection with an average to avoid being stuck in loops. When applied to ptychographic algorithms, these concepts gave rise to the difference map and relaxed averaged alternating reflections algorithms [22].

B. Gradient-based Optimization
Gradient-based iterative methods are the methods of choice to optimize nonlinear objective functions in machine learning. This class of algorithms is very flexible, with many variations and acceleration strategies. They have been used to solve the phase-retrieval problem [29,11], with applications ranging from ptychography and coded-diffraction imaging to defocus-based implementations.
One wants to find an estimatex that minimizes a loss function L(x, y), for instance the square loss L 2 (x, y) = y − |Ax| 2 2 . At each iteration k, a step of size η k toward the negative gradient of L is taken to refine the solution, leading to the procesŝ The choice of the loss function is of crucial importance since it determines the performance of the algorithm. An educated guess of a suited loss function can be made if the noise statistics are known. Typical noise models are Poissonian (for shot-noise limited measurements) or additive Gaussian [5]. In this case, the loss function corresponds to the negative log-likelihood function − log(p out y |Ax * | 2 ). Typically, the loss function L 2 corresponds to the negative log-likelihood of theadditive-Gaussian noise case. As the optimization is non-convex, gradient descent may get stuck in local minima or saddle points. Many acceleration strategies are at our disposal to circumvent these difficulties, from conjugate gradient and Nesterov acceleration to second-order methods. Another classical technique in machine learning to escape local minima is stochastic gradient descent, in which the gradient of the loss function at each iteration is computed using a partial dataset only. Intensity images are often used one at a time, akin to sequential projection algorithms, although this introduces noise in the gradient-descent updates.
Gradient descent represents a preferred tool to solve the phase-retrieval equation. Nevertheless, the theoretical bases of gradient-based algorithms still represent an active research topic. For example, the choice of the initial estimate greatly impacts the performance: in Section IV-D we shall describe algorithms that yield informed initializations.

C. Convex Relaxations
In contrast to gradient descent, which is a universal optimization technique, convex relaxations exploit the quadratic form of the phase-retrieval equation.
Convex relaxations recast the phase-retrieval problem using auxiliary variables. We seek a slightly different or relaxed formulation of the problem, which makes it convex. Convex optimization then provides a large variety of algorithms, with stronger convergence guarantees. For instance, PhaseLift [10] replaces the unknown vector x * ∈ C d by the matrix X * = x * (x * ) H ∈ C d×d . The quadratic equation (2) then becomes linear in the components of X * since one now has that y i = a H i X * a i for i = 1, . . . , n.
In its original form, the phase-retrieval equation would require one to seek a rank-one solution. As this rank constraint is itself non-convex, PhaseLift looks for a minimal-trace solution instead, leading to the convex program minimize Tr(X), subject to X ≥ 0 and y i = a H i Xa i for i = 1, . . . , n. Written in this form, the problem falls into the category of semidefinite programming, where one can readily use an off-the-shelf solver. The output is then given by the leading eigenvector of the final estimate X. Other convex relaxations have been developed [30,31], a detailed list being reviewed in [9].
The issue associated to semidefinite programming lies in its computational complexity, as one needs to build a matrix in C d×d . To illustrate this issue, consider a 100 × 100 image, which would correspond to d = 10 4 . Standard convex-relaxation techniques would then require one to perform optimization on a d × d = 10 8 -dimensional space. For this reason, semidefiniteprogramming strategies have only been used as a proof of concept in real imaging settings, although promising sketching methods to keep the computation tractable have been proposed recently [32].

D. Spectral Methods
To enhance the performance of the possibly expensive methods described in the previous sections, practitioners often require an "informed" initializationx (0) , preferably one that is close to optimum. Recently proposed, the privileged class of algorithms to obtain such initializations at a low computational cost are spectral methods [11]. These methods build an estimatex as the principal eigenvector of a weighted covariance matrix, defined from the intensities {y i } and the sampling vectors {a i } as Intuitively, since y i = |a H i x * | 2 , this scalar measurement is minimal when a i and x * are orthogonal, while it increases when the two vectors are correlated. For any increasing preprocessing function T , the matrix in (10) thus gives "more weight" to those sampling vectors a i that are more aligned with x * .
The performance of a given spectral method heavily depends on the choice of the preprocessing function T . Various choices have been considered in the literature, for example the trimming scheme T (y) = y1{|y| ≤ t} (with possibly t = ∞, meaning in this case that T is the identity). Interestingly, an optimal preprocessing method has been derived for the random models R.1 and R.2 and any type of noise [12,33]. In particular, in a noiseless setting it reads 1 T (y) = 1 − y −1 .
In practice, the simplicity of spectral methods ensures that their computational cost remains very low. The leading eigenvector can basically be recovered with power iterations or its accelerated variants 2 . Importantly, the convergence toward the leading eigenvector is ensured, a guarantee that is missing in many other nonlinear optimization strategies. Unlike convex relaxation methods, spectral methods do not involve a lifting process either, thus simplifying their adoption in imaging settings such as ptychography [34]. The solution estimated via spectral methods is then typically refined with an iterative optimization algorithm such as gradient descent.
In Figure 2 we illustrate the performance of some implementations of spectral methods, in the noiseless phase-retrieval problem with A being a complex Gaussian matrix. The optimal spectral method M(T ) yields almost perfect correlation with the true solution x * , even without subsequent refinement by gradient descent. Results are also compared with bare gradient-descent reconstructions. They show how all these theoretical findings transfer very well to the recovery of a real image.

E. Bayesian Algorithms
Bayesian estimators are a widely used and powerful ensemble of algorithms that leverage the properties of the posterior distribution, which is given by Bayes' law: Here, y are the observations, p out is the likelihood, and p 0 a prior distribution, typically Gaussian if one does not want to enforce any particular structure on the reconstructed image. Bayesian algorithms are not based on the minimization of a loss function, but rather on the computation of estimates defined from the posterior distribution. Thus, they can potentially exploit more efficiently the available information in the measurements. In this regard, an important class of iterative Bayesian algorithms is approximate message passing (AMP) [35], which is based on high-dimensional asymptotics of the belief-propagation algorithm [36]. The power of this class of algorithms arises from several insights. First, AMP is known to be optimal among a large class of polynomial-time algorithms for many random models [37] and is conjectured to be so for many more. Moreover, it is well understood, in the sense that one can analytically track its performance in the high-dimensional limit. However, AMP is computationally heavy compared to simpler algorithms, such spectral methods or gradient-based optimization procedures. More generally, the main limitation of Bayesian algorithms is that they assume the knowledge of the generative model (i.e. of the model of the noise present in the observations, and of the distribution of the random sensing matrix). While algorithmic procedures such as AMP can be used off-theshelf beyond these assumptions (e.g. for a non-random A), their performance and convergence properties are no longer guaranteed 3 .
We show in Figure 3 the performance of AMP in noiseless phase retrieval with a complex Gaussian matrix. We obtain a very good agreement with the theory (box on asymptotic guarantees) and, as for the spectral methods [33]), a very consistent performance for the recovery of real images.

A. Conventional Regularization
When faced with a difficult phase-retrieval problem, an appropriate regularization can be an efficient solution to enhance recovery of the desired images. It introduces prior information on the class of images to recover, which will help the reconstruction algorithm. This concept has long been associated with phase retrieval, since the first projection algorithms were based on support and positivity constraints.
Most regularization approaches can be formulated as the optimization problem arg min The first term enforces consistency with the intensity measurements y (data-fidelity term). It corresponds to the metric in the gradient-based reconstruction algorithms discussed in Section IV-B. The second term penalizes unrealistic estimates to favor suitable properties on x. In the Bayesian setting of (12) with maximum-likelihood estimation, the two terms correspond to p out and p 0 , respectively. In Section IV-E, the prior was set to a generic Gaussian distribution, but it can be modified to introduce other prior information on the recovered image.
Indicator functions may be used to enforce nonnegativity or support constraints. Another frequently studied regularization term is x 1 , the L 1 norm of x. It promotes sparsity in the final reconstruction and bears deep links with the field of compressed sensing. Sparsity in a transform domain can often be interesting for images, such as sparsity of Lx with L the discrete gradient operator for total-variation regularization. To tackle these potentially non-differentiable compositeoptimization problems, proximal algorithms alternatively minimize the loss L and the regularization R using gradients or projection operators. A large variety of algorithms have been designed for this task. We refer the Finally, a comparison with randomly initialized gradient descent is shown, which does not achieve a positive correlation at any α ≤ 5. In the last two columns, the optimal spectral method was chosen.
Weak Recovery Threshold: When do reconstructions start to be possible?
Complementary to injectivity studies, many recent works have focused on the question of weak recovery, in other words when does it become possible to nontrivially estimate the signal x * ? In mathematical terms, this means that there exists a procedure that is polynomial-time in d and that returns an estimatorx satisfying Note that anx uniformly sampled on the sphere (i.e. completely random) satisfies with high probability | x, A series of precise results on the weak recovery of phase retrieval under a random design assumption were obtained in the recent years. In particular, under R.1 or R.2, [15,16] showed that if n, d → ∞ with α = n/d, there exists a sharp threshold α WR distinguishing between impossible (α < α WR ) and possible (α > α WR ) weak recovery. For example, in noiseless complex phase retrieval with a Gaussian measurement matrix, this yields α WR = 1, as seen in Figures 2 and 3. a Note that the "optimal" spectral method we described in Section IV-D has been shown to reach the asymptotic weak-recovery threshold, i.e. it satisfies Eq. (11) for any α > α WR [15,33]. a [16] provides an explicit equation solved by αWR, solely as a function of the singular-value spectrum of the measurement matrix A. For example, in complex noiseless phase retrieval, αWR is the only solution to α = 2(E(λ)) 2 /E(λ) 2 , in which E(·) is taken with respect to the eigenvalues {λ} of AA H /d. Note that the right-hand side depends on α in general, making this an implicit equation.
reader to the tutorial [27] on conventional regularization and sparsity for more details.
Machine learning algorithms and, in particular, convolutional neural networks (CNNs), have recently emerged as powerful new tools to analyze natural images. Remarkably, they achieve state-of-the-art performance on many image benchmarks such as classification, object detection, denoising, or image generation. Not surprisingly, CNNs have also been deployed for image reconstruction as they can provide a very strong prior for natu-ral images. This prior consistently outperforms classical regularization on natural images. There are several ways to introduce it in our phase-retrieval algorithms.

B. Plugging-in a Denoiser
CNNs are remarkably effective for denoising, in other words, for the retrieval of an image from noisy observations. They are typically trained on a large dataset made of ground-truth images along with their corrupted Correlation | x, x * | achieved by AMP in noiseless phase retrieval with a complex Gaussian matrix. We compare the theoretical prediction, known as "state evolution" (blue curve), with actual runs of the algorithm. The results are very consistent whether the signal x * is synthetic (i.e. random) or an actual image, achieving perfect recovery for α 2.03. The inset of the left column compares the performance of AMP with the information-theoretic (IT) performance, which is the optimal error that any algorithm can reach, no matter its running time. We uncover a computational gap (or "hard phase") for 2 < α 2.03, for which AMP is not able to reach the IT performance [16]. We use an efficient modular implementation of AMP in a wide class of inference procedures (including phase retrieval), publicly available in [38].
Perfect recovery: When can we reconstruct "almost exactly?" Crucially, injectivity studies inherently consider "worst-case" scenarios on the choice of the signal x * to recover and the measurement matrix A, and do not leverage the high dimensionality of the problem. In order to go beyond this framework, the perfect recovery (or full recovery) phenomenon also attracted interest in the past years. This question is fundamentally high-dimensional and can be formulated as follows: for d 1, when is α = n/d large enough so that, for any typical x * and A, there exists an estimatorx of x * (with x = 1) that satisfies | x, x * | 1? Leveraging prior knowledge on the structure of the signal x * (i.e. giving a sense to the typicality), the question above has been answered in the very large class of random models R.1-R.2 [16]. In particular, perfect recovery can be achieved much before the typical injectivity threshold α = 4, e.g. at α 2 for Gaussian matrices A and a uniformly-sampled x * , cf. Fig. 3. Perfect-recovery studies turn out to be more relevant to practical phase-retrieval setups than injectivity analyses. In these setups, the measurement matrix and signal are not worst-case but rather constrained by the properties of the physical problem (e.g. a natural image). Then, one is often willing to accept vanishingly small errors.
versions, which corresponds to a self-supervised learning scheme.
Denoising networks can then be deployed at each iteration of an optimization framework. This applies in particular to Plug-and-Play (PnP) methods, which are a class of proximal algorithms commonly used to solve (13). These solvers alternate between two distinct steps: (a) decrease a data-fidelity loss and (b) correct the current guess by enforcing regularization. This class of algorithms is very flexible. For example, step (a) can represent a proximal operator of the data-fidelity term or a gradient-descent iteration, leading to the alternatingdirection method of multipliers (ADMM) and fast iter-ative shrinkage-thresholding algorithm (FISTA) respectively. As discussed in V-A, step (b) has traditionally been a proximal operator of a well-defined regularization term R. PnP methods replace this step by the denoising CNN.
Even though convergence is in general not assured and can only be proven in restricted settings, PnP with a denoising network is extremely flexible and has been used with great success in many applications. For instance, such an approach enables computationally efficient large-scale phase retrieval with state-of-the-art performance [41].
Likewise, regularization by denoising (RED) [42] also Asymptotic guarantees: Can we characterize the performance of reconstruction algorithms? A considerable effort has been devoted to characterize the performance of the diverse algorithms presented in Section IV, with a particular emphasis on the theoretical analysis of the error for large n, d, often under randomdesign assumptions (cf. Section III-E). For instance, the Bayesian algorithms of Section IV-E and, in particular, AMP (which conjecturally reaches the optimal polynomial-time performance) can be analytically studied in the high-dimensional limit with a set of tools known as State Evolution [16,35]. On one hand, these allow us to establish the theoretical solid lines in Figure 3, which agree very well with finite-size simulations of AMP. On the other hand, the spectral methods described in Section IV-D can be tackled more directly using tools of random-matrix theory and the authors of [12] derive analytical formulas for the asymptotic performance achieved by a wide class of spectral methods. Finally, gradientbased algorithms (Section IV-B) can also be analyzed in high dimension. For instance, [39] strikingly explained how over-parameterization greatly enhances their performance in phase retrieval. Another avenue to understand the trajectories of (stochastic) gradient descent in general settings is Dynamical Mean Field Theory [14,40]. We summarize in Table I these performance analyses in the specific case of noiseless complex-valued phase retrieval with the Gaussian/i.i.d. model R.1. Note that many of these results have been extended to more general noise and measurement models. Despite missing pieces and several simplifying assumptions, theoreticians have reached a very precise understanding of the performance of phase-retrieval algorithms. Nevertheless, almost all these studies rely on high-dimensional random models which, even though they can be quite refined, fail to describe many practical setups such as coded diffractive imaging with structured Fourier matrices or low-dimensional configurations. In particular, the theory of Bayesian methods and of spectral methods has been empirically shown to fail when the eigenspaces of A are very structured (e.g. with vanilla Fourier matrices). In this sense, theoretical asymptotic guarantees for structured matrices are still an open problem. exploits inner-loop denoisers to improve the quality of the reconstructions. However, the regularization term is now where f : C d → C d is the denoiser. This formulation makes the regularization term differentiable, which offers an interesting alternative to PnP. See [43] for a more indepth discussion.

C. Generative models
Generative Adversarial Networks (GANs) oppose two networks, a generator G and a discriminator D. The goal of the generator is to generate realistic images that follow the distribution of a dataset of natural images. It receives as input a random vector z also called latent variable and outputs x = G(z). Then, the task of the discriminator is to choose whether a given image is from the real dataset or some output of the generator. In this adversarial scheme, the two networks are trained together. Throughout training, the generator learns a manifold of realistic images while the discriminator gets better and better at discerning real and fake images.
To use such a generator for phase retrieval, one constrains the image x to be an output of the generator x = G(z) [44]. Gradient descent is performed directly on the latent variable z. This re-parameterization constrains the search space to the learned manifold of generator outputs. With this point of view, one can not only produce a point estimate of the unknown signal, but also perform a Bayesian sampling of the posterior distribution. This enables the computation of statistical estimates such as the spatial distribution of the variance of the final image [45].

D. Deep Image Prior
The machine-learning approaches of Sections V-B-V-C rely on a dataset of images for training. This raises questions about overfitting: by design the final image will resemble a realistic image from the training set while crucial distinctive information may be discarded. Deep Image Prior (DIP) has been proposed to remove any dependence on a training set. Here, one assumes that the estimates are generated from a CNN x = G w (z), with w the weights of the network and z a random input. Optimization is then directly performed on the weights w of the CNN. DIP has been implemented in several phase-retrieval applications [46,47].
While both the GAN and DIP approaches rely on a generative network, their re-parameterization differ considerably. The former optimizes a low-dimensional latent vector z, while the second tunes the CNN weights w. The dimension of w is typically larger than the dimension of x. As a result, DIP can overfit any input x,  even random noise. Nevertheless, it has been observed that it naturally converges at first toward smooth natural images. Regularization comes from the architecture of the CNN itself but requires early stopping to avoid fitting noise. In the end, the exact principle behind DIP regularization remains an open question.

E. Discussion
Historically, the first attempts at the application of deep learning targeted a direct inversion of the phase retrieval problem. From the intensity measurements y, they trained a network to output the desired image x and applied it to ptychography [48] or holography [49]. This strategy corresponds to an usual supervisedlearning scheme. However, the algorithms following this strategy are specific to one particular imaging application with a specific forward model, noise level, or camera sampling. This greatly limits their range of applicability. The methods described in Sections V-B-V-D belong to a second generation of algorithms with increased robustness and flexibility.
All these various strategies demonstrate how deeplearning methods offer a powerful alternative to classical regularization. They open new perspectives for solving the difficult nonlinear equation of phase retrieval. While our understanding is progressing fast, many of these methods remain black boxes that still need to be deciphered. In practice, one also needs to be aware of the notorious instabilities of neural networks [50], which can sometimes create artefacts (hallucinations) or remove certain features of interest. These issues are particularly troublesome and difficult to detect. The final result may look deceptively realistic (in reason of the regularization inherent with deep learning) even though the algorithm has actually failed.
From their inception, phase-retrieval algorithms have relied on prior information about the object to reconstruct. Regularization has gone a long way and evolved substantially, from nonnegativity and support constraints to deep learning, but the full potential and limits of the most recent tools at our disposal are not yet clearly established. The majority of the network-based proposed works require a training dataset, which may impede certain applications, typically, in medical setups. DIP and other unsupervised schemes offer an interesting alternative and may pave the way to powerful and controlled regularization strategies based on learning.

VI. CONCLUSION
We have described the many facets of phase retrieval, a problem which is central to a broad range of applications in computational optics. Various experimental configurations and reconstruction algorithms have been proposed to solve the missing-phase problem. To classify such a diversity of applications, we have outlined a unifying framework of many phase-retrieval models. We have also provided a concise overview of the latest phase-retrieval algorithms, to help practitioners choose one that fits their own application.
Viewed through the prism of machine learning, the links between phase retrieval and neural networks are twofold. First, the parallel between the defining equation (1) and single-layer neural networks has enabled the application of many theoretical machine-learning studies, pushing forward our understanding of this nonlinear equation. Second, the unrivaled performance of neural networks when dealing with images brings strong promises for future applications, but also comes with challenges on how to counteract unstability issues.
Regarding future developments, the introduction of more randomness in the illumination and/or detection could facilitate phase recovery in experiments, by benefiting from the strong theoretical results obtained for the random setting. However, many realistic phase-retrieval models still remain unexplored from a theoretical viewpoint, calling for more studies. With this document, we hope to promote a fruitful exchange between different disciplines, from applied physics to theoretical computer science.