Robust phase retrieval with the swept approximate message passing (prSAMP) algorithm
In phase retrieval, the goal is to recover a complex signal from the magnitude of its linear measurements. While many well-known algorithms guarantee deterministic recovery of the unknown signal using i.i.d. random measurement matrices, they suffer serious convergence issues some ill-conditioned matrices. As an example, this happens in optical imagers using binary intensity-only spatial light modulators to shape the input wavefront. The problem of ill-conditioned measurement matrices has also been a topic of interest for compressed sensing researchers during the past decade.
In this paper, using recent advances in generic compressed sensing, we propose a new phase retrieval algorithm that well-adopts for both Gaussian i.i.d. and binary matrices using both sparse and dense input signals. This algorithm is also robust to the strong noise levels found in some imaging applications.
This paper considers the fundamental problem of recovering a complex signal, , from magnitude of its linear projections. This problem is called phase retrieval (PR). Indeed, in many imaging setups, detectors (for instance, CCD cameras) are fundamentally intensity-only. Getting access to phase measurements may not be possible, or may involve a significantly more complex physical setup, e.g. with interferometry. Some of these applications include X-ray crystallography harrison93 (), X-ray diffraction imaging bunk07 (), optical imagers walther63 (); liutkus14 () and astronomical imaging fienup87 (). PR problems in the presence of additive noise be may formulated dremeau15 () as:
where known (measured) output, is the known complex projection matrix, is the unknown input, and is the ”noise” term - upon which some statistical assmuptions are made. Many methods have been reported in phase recovery while is the Fourier transform or a random matrix with iid Gaussian coefficients. These methods include, but are not limited to, convex relaxation algorithms such as phaseLift candes13 () and phaseCut waldspurger15 (), error reduction algorithms such as Gerchberg and Saxton gerchberg72 () and Fienup fienup78 () and several variants of them marchesini07 (); netrapalli15 () and spectral recovery method alexeev14 ().
Here, we are interested in the more challenging problem of recovering using a binary projection matrix, . This is the situation we face in real imaging applications using binary intensity spatial light modulators (SLM) such as digital micromirror devices (DMD) dremeau15 (). Using these ill-conditioned matrices, one is often faced with convergence issues with most of the afore-mentioned algorithms.
Signal recovery using ill-conditioned is also a challenging problem in other signal processing fields. Recently, in compressed sensing there have been attempts to reconstruct a sparse signal using generic matrices vila14 (); cakmak14 (); manoel14 (). Based on a Bayesian approach to a well-known compressed sensing approach which is called approximate message passing algorithm (AMP) maleki10 (), in manoel14 () the authors develop an idea which demonstrates good convergence properties over ill-conditioned noisy matrices. The extension is called swept AMP (SwAMP).
Following the Bayesian method, there are also a few phase retrieval methods such as phase retrieval generalized AMP (prGAMP) schniter15 (). By utilizing a magnitude only output prior, prGAMP reaches near optimal results to the classic PR problem with a smaller number of measurements. Phase retrieval variational Bayes expectation maximization (prVBEM) dremeau15 () is a mean-field variational Bayes phase retrieval technique that was developed for the task of calibrating the transmission matrix of a strongly scattering material, using binary measurements dremeau15 (). Although prVBEM has both small complexity and robustness to strong noise, its application has only been demonstrated in the context of light focusing rajaei15 ().
In this paper, we mix the idea of SwAMP with phase retrieval strategies, in order to solve (1) over binary projection matrices. The algorithm is called prSAMP and is already partially presented in rajaei15 ()111C source code to produce the same results as this paper is accessible at corresponding IPOL web page http://www.ipol.im.. In the context of compressive imaging though scattering material liutkus14 (), we here show that prSAMP can effectively deal with the phase-less recovery problem using intensity-only SLM. This yields to these two different problems in calibration and recovery steps: 1) complex input and binary measurement matrix and 2) binary input and complex measurement matrix. Obviously, the special case of complex input and matrix as addressed by most PR algorithms is also solvable by the algorithm.
In this section a brief summary of the notations that is used throughout the paper, is provided. As usual, scalars, vectors and matrices are written in small regular-face, small bold-face and capital bold-face letters, respectively. The th entry of a vector is denoted by and the th column of matrix by . The and operators stand for vector and element-wise multiplication. We also use and for element-wise square power and division. In algorithms, a function is represented as that is defined either in text or in another algorithm.
Iii prSAMP Algorithm
The proposed phase retrieval swept approximate message passing (prSAMP) algorithm is a mixture of two ideas in compressed sensing (CS) and phase retrieval, in addition to some modifications to work in 2D image recovery. The first part is swept approximate message passing algorithm (SwAMP) manoel14 (), which is one of many variants of approximate message passing (AMP) method maleki10 () in compressed sensing. In the context of CS, AMP is an iterative algorithm that reconstructs a sparse signal from a set of under-determined linear noisy measurements, , where . Figure 1 shows the statistical approach in the AMP method krzakala12 (), where the algorithm starts from initial posterior estimates of signal average and variance and . It then follows three main steps iteratively; 1) calculate output mean and variance variables, and , 2) calculate input maximum likelihood terms, and , which are also called AMP Gaussian fields and 3) use AMP denoisers to update input signal mean and variance, and . The AMP denoisers carry prior knowledge of the input unknown signal. Later in this paper, we define two denoisers for binary and Gaussian signals. AMP has been shown to converge to optimal solution while working with i.i.d. gaussian projection matrices bayati11 (). However, it does not necessarily converge for generic matrices caltagirone14 (). This is where the idea of SwAMP brings up.
SwAMP is a simple change in step 2 of AMP. Instead of standard parallel calculation, a sequential, or swept, random update of AMP maximum likelihood variables is suggested which shows significant stabilization of the AMP loop while working with different non i.i.d. and/or ill-conditioned projection matrices manoel14 ().
To extend AMP framework to our PR problem of eqn. (1), we first use the generalized AMP (GAMP) rangan11 () which is an extension of AMP for arbitrary output channels, i.e. . This adds an output function, , which is dependent on the stochastic description of . Normal CS problems follow an additive white Gaussian noise (AWGN) channel as output prior. For the PR case, we follow what is proposed in a GAMP-based phase retrieval algorithm, which is called prGAMP schniter15 () and formulates for and .
Algorithm III describes phase retrieval version of SwAMP, denoted as prSAMP, which combines the swept update ordering and the phase retrieval output channel in the AMP iteration. Beside the intensity measurements, , and the projection matrix, , the algorithm has a few other input parameters. These include the two stopping parameters, maximum number of iterations, , and the precision threshold, . The algorithm stops if it reaches iterations or if the difference between two successive estimations is less than , . The other parameter is . During prSAMP iterations, variance terms may become negative or very small. This prevents the algorithm to improve its current estimation which happens often during the first iterations. In these cases the bad variance values replace by . There are also two damping parameters, and . Damping is necessary in case of ill-conditioned matrices. We use the first damp factor for and variables in step 2 and the second for 2D signals. If the input signal is actually vectorized version of a 2D image, after step 3 we add one step to take into account the 2D relation between elements. Here, we employ a simple damping with respect to which is local median in iteration but more sophisticated 2D priors may establish for better smoothness in recovered signal. Finally, we have input and output prior functions and their associated parameters which are defined in separate algorithms.
In the main loop, the algorithms starts by estimating output average and variance terms, and . Then we have output prior applies over these variables to calculate and mean and variance terms. In the case of phase retrieval experiment, these are defined distinctly in Algorithm III but in normal compressed sensing with additive white Gaussian noise (AWGN) the calculation is straightforward. AWGN output prior indicates and for and variables, respectively.
In the second step, we have the sequential swept iteration for maximum likelihood terms, and . It has been claimed in the SwAMP original paper that random computation of involved variables result in better convergence therefore we also follow the same method. After each index is calculated from input signal, the updates should apply over output channel variables. Finally, the estimate of unknown input signal is returned as variable.
Considering a circular Gaussian additive white noise in measurements, follows a Rician probability density function which is the basis for a PR output channel derivation in prGAMP paper schniter15 (). We also follow the same formulation. Algorithm III explains the PR output prior. Here, and functions are respectively and -order modified Bessel functions of first kind.
As we mentioned earlier, in this paper we are interested in solving PR problem in two cases: 1) calibration step with and and 2) recovery step with and . Therefore, a Gaussian input prior for the calibration phase is a reasonable choice as it is described in Algorithm III. Furthermore, a possible binary prior is explained in Algorithm III for the reconstruction step.
To reconstruct the complex signal (up to a global phase) using its intensity-only projections, the size of the measurement vector should be at least - it has been established recently that, in a generic case, measurements are required bodmann15 () to recover a unique . This means that prSAMP follows a computational complexity of which is a bottleneck for real-time imaging. Due to sequential nature of swept loop we can not solve this scaling issue directly but there are two possibilities to alleviate it: 1) in calibration phase since different rows of projection matrix are inherently independent, the algorithm is fully parallel. In the supplementary files, two extensions of prSAMP using OMP and MPI parallel tools, are provided. 2) the other enhancement option is an idea we call block-based phase retrieval rajaei16 (). This block-based PR method starts by splitting the input problem into , sub-problems, where , and . The sub-problems are then solved in parallel. Finally, all the partial results are merged with a few extra global measurements, by applying a low-dimension global phase tuning step. In this way the order of prSAMP algorithm breaks down into . This comes at a price of being able to design the measurement matrix in a general block-diagonal manner which is the case in any physical systems where one can probe the whole object by parts. Block-based prSAMP is extended in the supplementary files using Matlab.
V Parameter study
There are two groups of parameters: first, the main prSAMP parameters and second, priors parameters. Depending on prior knowledge of input and output signals, and , we may require to provide parameters like, noise variance in measurements, , an estimation of input sparsity level, , or input mean and variance, and , in case of Gaussian input prior. Beside these obvious prior-dependent parameters, there are a few main parameters that they play important role in algorithm convergence. The initial estimation of unknown input signal, , and the damping factor, , are the two most important ones.
As it is well-known the compressive phase retrieval problem generally suffers from convergence to local minima vila14 (). Empirical studies show that the situation is worse while working with ill-conditioned non-gaussian i.i.d. random matrices manoel14 (). In case of Gaussian input signals, like what we have in the calibration phase, using a pseudo random generator to initialize seems a reasonable choice. Afterwards if the algorithm diverges, a complete restart with a new random initial vector is necessary. Multiple restarts was first suggested in prGAMP paper schniter15 (). The solution that yields to lowest normalized residual, , is selected as the algorithm output. In case of other types of input signals, a good initial point would guarantee the algorithm convergence. For example, in recovery phase of our optical imager, we employ a low resolution (LR) version of input image. This LR signal may be gathered from negative output of DMD array or numerically estimated based on specific image database. For signals of length , LR version of is a good start point. Depending on initial point confidence, variance vector is selected from interval. In calibration and recovery steps we empirically set values at 0.5 and 0.1, respectively.
The other important parameter is the damping factor, . Damping slows down the convergence rate of the algorithm, and hence prevents being stuck into a possibly wrong local minima, while still keeping information from previous iterations. Here, is a scalar from interval where 0 indicates no damping situation. In case of ill-conditioned projection matrices we need more damping. In our experiments we use 0.9 and 0.2 values for calibration and recovery steps, respectively.
In case of 2D input signals we have another damping parameter, . In this paper as a 2D prior we used a simple damping step which mixes the current solution,, with a representative of its neighborhood, . The representative is median over a block centered at element . This may improve by taking into account learned priors like the RBM prior as it is proposed in tramel15 (). The more sophisticated priors usually come with the price of an offline learning step. Hence, since our simple damping prior provides satisfactory results, as it is shown in the next section, we left further improvements to the interested reader.
The other parameters include: maximum number of iterations, , precision factor, , and negative variance factor, . Number of iterations is usually a factor of number of nonzero elements in input unknown signal, . In calibration step, with a full rank input vector, we set at empirically. But for small , it is necessary to let algorithm pass the initial oscillations. In small , we may use .
Precision factor, , is another measure of convergence which ensures a minimum difference between two successive solutions. A difference less than indicates the algorithm is iterating around a local minima and, Hence, there is no progress.
Finally, negative variance factor, , is employed in case of resulting a negative variance term. There are various variance variables in prSAMP algorithm like: and variables in the main algorithm and , and variables in priors. These terms have to be positive and not extremely small. Therefore, in case of negative or very small variance terms, is used as a replacement value. This parameter should be sufficiently large and in the range of because negative variance indicates a bad situation in prSAMP iteration and we should set the variance at a large value to let algorithm converge to another mean point.
Vi Experimental results
In this section we investigate the application of using prSAMP algorithm to solve phase retrieval problem (1) in two different situations; first, and and second, and . Using binary transition matrix to recover complex input signal, figure 2 shows the phase transition plot for and snr equal to 30 dB (). The error is measured in terms of normalized mean square error (NMSE) between original and the recovered signal after compensating the global phase shift. Each point in the plot is the lowest NMSE obtained by prSAMP in 50 distinct trials. As a result of ill-conditioned binary measurement matrix, the damping factor, is set to 0.9. A phase transition curve is generated by applying a NMSE threshold of 0.2. The plot confirms the recently established rate of to reconstruct dense signal in phase retrieval regime. The effect of increasing the number of measurements on recovered signal is shown in figure 3 using the same settings as the previous experiment and a random input. Here, we have .
For , reconstruction performance of prSAMP is studied at four different noise levels (snr equal to 30, 20, 10 and 5 dB) and in figure 4. According to this experiment in case of strong noisy measurements, after adding more samples does not improve the results significantly.
In the second experiment, prSAMP is applied to the problem of reconstructing binary random input signals. Figure 5 shows the corresponding phase transition plot. Except which is set to 0.2 the other parameters are similar to the first experiment. The recovery error is measured in terms of best correlation out of 50 runs. Here, the number of necessary measurements for complete recovery is decreased significantly at different sparsity levels probably due to binary input prior. This fact also has been shown in figure 6 which plots a random binary signal with and its two reconstructions at and . Comparing to figure 3 complete recovery is happened using significantly less measurements ().
Vii Computational complexity
Finally, it would be interesting to have a brief discussion on computational complexity of prSAMP algorithm. Even though, as our experiments show, the algorithm performs well for ill-conditioned matrices and strong noise situations, it does not scale well as the size of input signal increases. In Algorithm III the number of iterations, , and measurements, , grow linearly with input size . Therefore, prSAMP follows a cubic computational complexity. In addition to this, the amount of data that the algorithm has to handle at least scales with . This is challenging at large inputs. In rajaei16 (), a block-based version of prSAMP has been proposed that can reduce computational cost and also memory requirements of original prSAMP by several orders of magnitude.
In this study, a new phase retrieval algorithm has been proposed, called phase retrieval swept AMP (prSAMP). prSWAMP is here numerically evaluated in two situations inspired by real imaging setups. In particular, prSWAMP solves the challenging problem of estimating a complex input signal using binary patterns. In reverse, we also show that prSAMP accurately estimates a binary unknown signal using a complex transmission matrix.
This research has received funding from PSL Research University under contract CSI:PSL, from the European Research Council under the EU’s 7th Framework Programme (FP/2007- 2013/ERC Grant Agreement 307087-SPARCS and 278025- COMEDIA) ; and from LABEX WIFI under references ANR-10-LABX-24 and ANR-10-IDEX-0001-02-PSL.
This work was also granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche.
- (1) B. Alexeev, A.S. Bandeira, M. Fickus, and D.G. Mixon, Phase retrieval with polarization, SIAM Journal on Imaging Sciences, 7 (2014), pp. 35–66.
- (2) M. Bayati and A. Montanari, The dynamics of message passing on dense graphs, with applications to compressed sensing, IEEE Trans. on Info. Theory, 57 (2011), pp. 764–785.
- (3) B.G. Bodmann and N. Hammen, Stable phase retrieval with low-redundancy frames, Advances in computational mathematics, 41 (2015), pp. 317–331.
- (4) O. Bunk, A. Diaz, F. Pfeiffer, C. David, B. Schmitt, D. K. Satapathy, and J. F. Veen, Diffractive imaging for periodic samples: Retrieving one-dimensional concentration profiles across microfluidic channels, Acta Crystallograph. Sect. A: Found. Crystallogr., 63 (2007), pp. 306–314.
- (5) B. Çakmak, O. Winther, and B. Fleury, S-AMP: Approximate message passing for general matrix ensembles, in IEEE Info. Theory Work., 2014, pp. 192–196.
- (6) F. Caltagirone, L. Zdeborová, and F. Krzakala, On convergence of approximate message passing, in IEEE Int. Symp. on Info. Theory, 2014, pp. 1812–1816.
- (7) E.J. Candes, T. Strohmer, and V. Voroninski, Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming, Communications on Pure and Applied Mathematics, 66 (2013), pp. 1241–1274.
- (8) A. Drémeau, A. Liutkus, D. Martina, O. Katz, C. Schülke, F. Krzakala, S. Gigan, and L. Daudet, Reference-less measurement of the transmission matrix of a highly scattering material using a DMD and phase retrieval techniques, Opt. Express, 23 (2015), pp. 11898–11911.
- (9) C. Fienup and J. Dainty, Phase retrieval and image reconstruction for astronomy, Image Recovery: Theory and Application, (1987), pp. 231–275.
- (10) J.R. Fienup, Reconstruction of an object from the modulus of its fourier transform, Optics letters, 3 (1978), pp. 27–29.
- (11) R.W. Gerchberg, A practical algorithm for the determination of phase from image and diffraction plane pictures, Optik, 35 (1972), p. 237.
- (12) R.W. Harrison, Phase problem in crystallography, J. Opt. Soc. Amer. A, 10 (1993), pp. 1046–1055.
- (13) F. Krzakala, M. Mézard, F. Sausset, Y. Sun, and L. Zdeborová, Probabilistic reconstruction in compressed sensing: Algorithms, phase diagrams, and threshold achieving matrices, J. Stat. Mech.: Theory & Exp., 2012 (2012), p. P08009.
- (14) A. Liutkus, D. Martina, S. Popoff, G. Chardon, O. Katz, G. Lerosey, S. Gigan, L. Daudet, and I. Carron, Imaging with nature: Compressive imaging using a multiply scattering medium, Sci. Rep., 4 (2014).
- (15) A. Maleki, Approximate Message Passing Algorithms for Compressed Sensing, PhD thesis, Stanford Uni., 2010.
- (16) A. Manoel, F. Krzakala, E.W. Tramel, and L. Zdeborová, Swept approximate message passing for sparse estimation, in Int. Conf. on Machine Learning, 2014.
- (17) S. Marchesini, Invited article: A unified evaluation of iterative projection algorithms for phase retrieval, Review of scientific instruments, 78 (2007), p. 011301.
- (18) P. Netrapalli, P. Jain, and S. Sanghavi, Phase retrieval using alternating minimization, Signal Processing, IEEE Transactions on, 63 (2015), pp. 4814–4826.
- (19) B. Rajaei, S. Gigan, F. Krzakala, and L. Daudet, Fast phase retrieval for high dimensions: A block-based approach, arXiv preprint arXiv:1602.02944, (2016).
- (20) B. Rajaei, E.W. Tramel, S. Gigan, F. Krzakala, and L. Daudet, Intensity-only optical compressive imaging using a multiply scattering material and a double phase retrieval approach, arXiv preprint arXiv:1510.01098, (2015).
- (21) S. Rangan, Generalized approximate message passing for estimation with random linear mixing, in IEEE Int. Symp. Inf. Theory, IEEE, 2011, pp. 2168–2172.
- (22) P. Schniter and S. Rangan, Compressive phase retrieval via generalized approximate message passing, IEEE Trans. on Signal Proc., 63 (2015), pp. 1043–1055.
- (23) E.W. Tramel, A. Drémeau, and F. Krzakala, Approximate message passing with restricted Boltzmann machine priors, J. Stat. Mech.: Theory & Exp., (2016). to appear.
- (24) J. Vila, P. Schniter, S. Rangan, F. Krzakala, and L. Zdeborová, Adaptive damping and mean removal for the generalized approximate message passing algorithm, in IEEE Int. Conf. on Acoustics, Speech, & Signal Proc., 2015, pp. 2021–2025.
- (25) I. Waldspurger, A. dâAspremont, and S. Mallat, Phase recovery, maxcut and complex semidefinite programming, Mathematical Programming, 149 (2015), pp. 47–81.
- (26) A. Walther, The question of phase retrieval in optics, Opt. Acta., 10 (1963), pp. 41–49.