Phase retrieval with a multivariate Von Mises prior: from a Bayesian formulation to a lifting solution

TL;DR

This paper introduces a Bayesian approach for phase retrieval that incorporates prior phase information via a multivariate Von Mises distribution, leading to a lifting optimization solution.

Contribution

It presents a novel Bayesian formulation for phase retrieval using a multivariate Von Mises prior, connecting it to a lifting optimization method.

Findings

Effective incorporation of phase priors improves retrieval accuracy

Formulation using Mahalanobis distance enables convex optimization

Lifting method provides a practical solution for phase recovery

Abstract

In this paper, we investigate a new method for phase recovery when prior information on the missing phases is available. In particular, we propose to take into account this information in a generic fashion by means of a multivariate Von Mises dis- tribution. Building on a Bayesian formulation (a Maximum A Posteriori estimation), we show that the problem can be expressed using a Mahalanobis distance and be solved by a lifting optimization procedure.

Equations

$$\displaystyle\mbox{\boldmath$y$}=|\mbox{\bf A}\mbox{\boldmath$x$}|,$$

$$\mbox{\boldmath$y$}=\mbox{\boldmath$\operatorname{Diag}$}\{\mbox{\boldmath$\phi$}\}^{H}\mbox{\bf A}\mbox{\boldmath$x$}+\mbox{\boldmath$n$}$$

$$\displaystyle p(\theta_{m})=\frac{1}{2\pi I_{0}(\kappa_{m})}\exp\big{(}\kappa_{m}\cos(\theta_{m}-\mu_{m})\big{)},$$

$$\displaystyle p(\mbox{\boldmath$\theta$})=$$

$$\displaystyle-\mbox{\boldmath$s$}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\mu$})^{\top}\mbox{\boldmath$\Delta$}\mbox{\boldmath$s$}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\mu$})-\mbox{\boldmath$c$}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\mu$})^{\top}\mbox{\boldmath$\Delta$}\mbox{\boldmath$c$}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\mu$})\big{)},$$

$$\displaystyle c_{m}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\mu$})\negmedspace=\negmedspace\cos(\theta_{m}\negmedspace-\negmedspace\mu_{m}),\quad s_{m}(\mbox{\boldmath$\theta$},\mbox{\boldmath$\mu$})\negmedspace=\negmedspace\sin(\theta_{m}\negmedspace-\negmedspace\mu_{m}).$$

$$\hat{\mbox{\boldmath$\phi$}}\triangleq[e^{j\hat{\theta}_{1}},\dots,e^{j\hat{\theta}_{M}}]^{\top}=\operatorname*{argmin}_{\begin{subarray}{c}\mbox{\boldmath$\phi$}\\ |\phi_{i}|^{2}=1\;\forall i\end{subarray}}||\mbox{\boldmath$\phi$}-\mathbf{1}_{M}||_{\Gamma_{\phi}}^{2}$$

$$\displaystyle(\Gamma_{\phi}^{-1})_{ik}=\left\{\begin{array}[]{ll}\Delta_{ik}&\text{ if }k\neq i,\\ \frac{1}{2}\kappa_{i}-\sum_{l\neq i}\Delta_{il}&\text{ if }k=i.\end{array}\right.$$

$$\displaystyle\hat{\mbox{\boldmath$\phi$}}_{\text{\tiny{MAP}}}$$

$$\displaystyle=\operatorname*{argmin}_{\begin{subarray}{c}\mbox{\boldmath$\phi$}\\ |\phi_{i}|^{2}=1\;\forall i\end{subarray}}\negmedspace\frac{1}{\sigma_{n}^{2}}||\mbox{\boldmath$y$}\negmedspace-\negmedspace\mbox{\boldmath$\operatorname{Diag}$}\{\mbox{\boldmath$\phi$}\}^{H}\mbox{\bf A}\mbox{\boldmath$x$}||^{2}_{2}$$

$$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad+||\mbox{\boldmath$\phi$}-\mathbf{1}_{M}||_{\Gamma_{\phi}}^{2}.$$

$$\displaystyle\hat{\mbox{\boldmath$x$}}_{\text{\tiny{ML}}}$$

$$\displaystyle=\operatorname*{argmin}_{\mbox{\boldmath$x$}}||\mbox{\boldmath$y$}\negmedspace-\negmedspace\mbox{\boldmath$\operatorname{Diag}$}\{\mbox{\boldmath$\phi$}\}^{H}\mbox{\bf A}\mbox{\boldmath$x$}||^{2}_{2},$$

$$\displaystyle=\mbox{\bf A}^{+}\mbox{\boldmath$\operatorname{Diag}$}\{\mbox{\boldmath$\phi$}\}\mbox{\boldmath$y$},$$

$$\displaystyle\hat{\mbox{\boldmath$\phi$}}_{\text{\tiny{MAP}}}\negmedspace$$

$$\displaystyle\quad\quad\quad\quad\quad\quad+\left\|(\mbox{\bf I}_{M}\quad-\mathbf{1}_{M})\cdot\left(\begin{array}[]{c}\mbox{\boldmath$\phi$}\\ 1\end{array}\right)\right\|_{\Gamma_{\phi}}^{2},$$

$$\hat{\mbox{\boldmath$u$}}=\operatorname*{argmin}_{\scriptstyle{\begin{subarray}{c}\mbox{\boldmath$u$}\\ |u_{i}|^{2}=1\;\forall i\end{subarray}}}\mbox{\boldmath$u$}^{H}\mbox{\bf Q}\mbox{\boldmath$u$},\;\textrm{where}$$

$$\text{mbox}\mathbf{Q}=\left(\begin{array}{cc}\text{mbox}\mathbf{M}+{\sigma }_n^2{\Gamma }_{\varphi }^{-1}& -{\sigma }_n^2{\Gamma }_{\varphi }^{-1}{1}_M\\ -{\sigma }_n^2{1}_M^{\top }{\Gamma }_{\varphi }^{-1}& {\sigma }_n^2{\sum }_{i,k}({\Gamma }_{\varphi }^{-1}{)}_{ik}\end{array}\right)\in {\mathbb{C}}^{(M+1{)}^2},$$

$$\begin{array}[]{cc}\operatorname{argmin}&\mbox{trace}\>\{\mbox{\bf Q}\mbox{\bf U}\}\\ \mbox{\bf U}\succeq\mbox{\bf 0}&\\ \operatorname{diag}\{\mbox{\bf U}\}=\mbox{\bf 1}_{M}&\\ \end{array}$$

$$\frac{|{\hat{\mathbf{x}}}^H\mathbf{x}|}{\Vert \hat{\mathbf{x}}{\Vert }_2\Vert \mathbf{x}{\Vert }_2},$$

$$\displaystyle p(\mbox{\boldmath$\theta$})\simeq\frac{1}{C(\mbox{\boldmath$\kappa$},\mbox{\boldmath$\Delta$})}$$

$$\exp \left(-\sum _i\left(({\Gamma }_{\varphi }^{-1}{)}_{ii}{\theta }_i^2-2({\Gamma }_{\varphi }^{-1}{)}_{i(i-1)}{\theta }_i{\theta }_{i-1}+o({\theta }_i^2)\right)\right),$$

$$\displaystyle(\Gamma_{\phi}^{-1})_{ik}=\left\{\begin{array}[]{ll}-\frac{a}{2\sigma_{\theta}^{2}}&\text{ if }k=i+1\text{ or }k=i-1,\\ \frac{1+a^{2}}{2\sigma_{\theta}^{2}}&\text{ if }k=i\neq M,\\ \frac{1}{2\sigma_{\theta}^{2}}&\text{ if }k=i=M,\\ 0&\text{elsewhere.}\end{array}\right.$$

$$\displaystyle||\mbox{\boldmath$\phi$}-\mathbf{1}_{M}||_{\Gamma_{\phi}}^{2}$$

$$\displaystyle=(\mbox{\boldmath$\phi$}-\mathbf{1}_{M})^{H}\Gamma_{\phi}^{-1}(\mbox{\boldmath$\phi$}-\mathbf{1}_{M})$$

$$=\text{mbox}trace\{{\Gamma }_{\varphi }^{-1}\}+\sum _{i,k}({\Gamma }_{\varphi }^{-1}{)}_{ik}-2\sum _i({\Gamma }_{\varphi }^{-1}{)}_{ii}\cos ({\theta }_i)$$

$$-\sum _i\sum _{k\ne i}({\Gamma }_{\varphi }^{-1}{)}_{ik}{e}^{-j({\theta }_i-{\theta }_k)}-\sum _i\sum _{k\ne i}({\Gamma }_{\varphi }^{-1}{)}_{ik}{e}^{-j{\theta }_i}$$

$$-\sum _i\sum _{k\ne i}({\Gamma }_{\varphi }^{-1}{)}_{ki}{e}^{j{\theta }_i},$$

$$=\text{mbox}trace\{{\Gamma }_{\varphi }^{-1}\}+\sum _{i,k}({\Gamma }_{\varphi }^{-1}{)}_{ik}-2\sum _i\left(\sum _k({\Gamma }_{\varphi }^{-1}{)}_{ik}\right)\cos ({\theta }_i)$$

$$+\sum _i\sum _{k\ne i}({\Gamma }_{\varphi }^{-1}{)}_{ik}\cos ({\theta }_i-{\theta }_k),$$

$$||\mbox{\boldmath$\phi$}-\mathbf{1}_{M}||_{\Gamma_{\phi}}^{2}\propto-\log p(\mbox{\boldmath$\theta$})$$

Full text

**Index Terms— ** Phase retrieval, multivariate Von Mises distribution, Mahalanobis distance, lifting.

1 Introduction

Since more than twenty years, phase retrieval has been a constantly filled topic. This is because the problem interests numerous application domains, from crystallography [1] to optical imaging [2]. Formally, it can be written as follows: given \mbox{\boldmathy}\in\mathbb{R}^{M}, recover \mbox{\boldmathx}\in\mathbb{C}^{K} such as

[TABLE]

where A is a $M\times K$ known complex-measurement matrix. Several answers to this non-convex optimization problem have been proposed, that we can roughly divide into three families: i) alternating-projection algorithms, where we can find the works of Gerchberg & Saxton [3], Fienup [4] or Griffin & Lim [5], which alternate projections on the span of the measurement matrix and on the object domain, ii) algorithms based on convex relaxations, such as the recent PhaseLift [6] and PhaseCut [7], which replace the phase recovery problem by relaxed problems that can be efficiently solved by standard optimization procedures, and iii) Bayesian approaches, which express the phase recovery problem as the solution of a Bayesian inference problem and apply statistical tools to solve it, such as variational approximations [8, 9].

In the above procedures, the phases are completely missing from the observations: only intensities or amplitudes are acquired. In this paper, we are interested in phase retrieval problems where phases are observed but marred by noise. At the interface between the last two above families, we propose a Bayesian formulation of the problem and resort to a lifting optimization procedure to solve it. A priori knowledge over observed phases through various probabilistic laws have been exploited in previous works [10, 11]. Compared to them, our approach presents two appealing novelties: i) it is generic in the sense that it can handle multivariate phase priors and thus arbitrary dependencies; ii) the proposed Bayesian optimization problem is cast into a generalization of the recently proposed PhaseCut problem [7], for which a number of efficient estimation procedures readily exist, including convex relaxations. The last point is made possible by exploiting a previously unseen connection between a multivariate generalization of the Von Mises distribution and the Mahalanobis distance.

2 Bayesian formulation

In this section, we introduce the Bayesian modeling that we propose to exploit in the following and discuss its link to the Mahalanobis distance, particularly interesting for the optimization procedure.

2.1 Observation model

Let $M$ sensors record $K$ complex signals through linear instantaneous mixing, in the presence of both additive noise and multiplicative phase noise. The noisy observation \mbox{\boldmathy}\in\mathbb{C}^{M} is then expressed as

[TABLE]

where $\mbox{\bf A}\in\mathbb{C}^{M\times K}$ is the mixing matrix, \mbox{\boldmathx}\in\mathbb{C}^{K} is the source signal, \mbox{\boldmathn}\in\mathbb{C}^{M} is the noise vector, \mbox{\boldmath\phi}=[e^{j\theta_{1}},\dots,e^{j\theta_{M}}]^{\top} is the phase vector with \mbox{\boldmath\theta}\triangleq[\theta_{1},\dots,\theta_{M}]^{\top} $\in]-\pi,\pi]^{M}$, the operator \mbox{\boldmath\operatorname{Diag}}\{.\} transforms row- or column-vectors into diagonal matrices and $\cdot^{H}$ denotes the complex conjugate transpose. For simplicity, we assume that the additive noise is zero-mean i.i.d. circular complex Gaussian with variance $\sigma_{n}^{2}$. Note that generalizing subsequent derivations to an arbitrary noise covariance matrix $\Gamma_{n}$ instead is straightforward with appropriate changes of variable.

2.2 Von Mises prior

In the literature, model (2) has been already considered in phase retrieval problems with a uniform prior on the phases $\theta$ (see e.g. [8, 9]). Here, we look for a more informative model enforcing uncertain structures on and between phases.

Considering phases naturally leads to directional statistics. Among them, the most familiar one is probably the Von-Mises distribution, defined independently for each variable $\theta_{m}$, $m\in\{1,\ldots,M\}$ as

[TABLE]

where $\kappa_{m}\in\mathbb{R}$ and $\mu_{m}\in]-\pi,\pi]$ are parameters of the distribution, and $I_{0}(.)$ is the modified Bessel function of the first kind of order 0. This distribution has been considered in the literature, e.g. in [10]. In practice, it is well-adapted to situations where we want to take into account prior information (such as the mean through parameter $\mu_{m}$ or the variance through $\kappa_{m}$) on the phases independently of one another. Its extension to the multivariate case is not straightforward and can take different forms [12]. In this paper, we assume $\theta$ to be distributed according to

[TABLE]

where C(\mbox{\boldmath\kappa},\mbox{\boldmath\Delta}) is a normalizing constant and functions $c$ and $s$ are respectively defined by, $\forall m\in\{1,\ldots,M\}$,

[TABLE]

The matrix $\Delta$ is real-symmetric with zeros on its diagonal and captures dependencies between phases. Without loss of generality111Assuming \mbox{\boldmath\mu}\neq[0,\ldots,0]^{\top} amounts to considering the observation model \mbox{\boldmathy}=\mbox{\boldmath\operatorname{Diag}}\{\widetilde{\mbox{\boldmath\phi}}\}^{H}\widetilde{\mbox{\bf A}}\mbox{\boldmathx} with \widetilde{\mbox{\boldmath\phi}}=\mbox{\boldmath\operatorname{Diag}}\{\mbox{\boldmathu}\}^{H}\mbox{\boldmath\phi}, \widetilde{\mbox{\bf A}}=\mbox{\boldmath\operatorname{Diag}}\{\mbox{\boldmathu}\}\mbox{\bf A} and \mbox{\boldmathu}=[e^{j\mu_{1}},\ldots,e^{j\mu_{M}}]., we will assume in the sequel that \mbox{\boldmath\mu}=[0,\ldots,0]^{\top}\triangleq\mbox{\boldmath0}_{M}. This multivariate extension of the Von-Mises distribution was suggested at the end of [12], but does not seem to have been extensively studied or used. We prefer it here over other alternatives due to the following result (proof in appendix A):

Lemma 1

Let \mbox{\boldmath\mu}=\mbox{\boldmath0}_{M}, and let \hat{\mbox{\boldmath\theta}}=[\hat{\theta}_{1},\ldots,\hat{\theta}_{M}]^{\top} maximize the multivariate Von Mises distribution (4). We have:

[TABLE]

where $||\cdot||_{\Gamma_{\phi}}$ denotes the Mahalanobis distance with covariance $\Gamma_{\phi}$, $\mathbf{1}_{M}\triangleq[1,\ldots,1]^{\top}$, and $\forall(i,k)\in\{1,\ldots,M\}^{2}$

[TABLE]

In other words, maximizing the density (4) can be cast as a quadratically-constrained norm-minimization problem. This type of problem is central in the classical phase retrieval literature (see [7]), but does not seem to appear when using other multivariate generalizations of the Von Mises distribution as phase priors, e.g., the one studied in [12].

3 Phase and Signal Estimation

3.1 Maximum a posteriori

Using Lemma 1, it follows that the Maximum A Posteriori (MAP) estimate of $\phi$ within model (2) and (4) writes:

[TABLE]

Following a similar idea as in [7], we couple this MAP estimation of the phase vector $\phi$ with a Maximum Likelihood estimation of the source signal $x$:

[TABLE]

where $\mbox{\bf A}^{+}$ stands for the Moore-Penrose pseudo-inversion of matrix A. Reinjecting this estimate in the MAP problem (10) leads to

[TABLE]

where $\mbox{\bf I}_{M}$ stands for the identity matrix and the Mahalanobis distance term has been re-written to be homogeneous in \mbox{\boldmathu}=[\mbox{\boldmath\phi}\quad 1]^{\top}\in\mathbb{C}^{M+1}. Using this trick, it follows that solving (3.1) is equivalent to solving the $M+1$-dimensional problem

[TABLE]

[TABLE]

with \mbox{\bf M}=\mbox{\boldmath\operatorname{Diag}}\{\mbox{\boldmathy}^{H}\}(\mbox{\bf I}_{M}\negmedspace-\negmedspace\mbox{\bf A}\mbox{\bf A}^{+})\mbox{\boldmath\operatorname{Diag}}\{\mbox{\boldmathy}\}. It is easily verified that if \hat{\mbox{\boldmathu}} is solution of (17), then \hat{\mbox{\boldmath\phi}}_{\text{\tiny{MAP}}}=\hat{\mbox{\boldmathu}}_{1:M}/\hat{u}_{M+1} is solution of (3.1).

Interestingly, when \Gamma_{\phi}^{-1}=\mbox{\boldmath0} (uninformative prior on phases), (17) is equivalent to the program proposed by Waldspurger et al. [7] for classical phase retrieval. They refer to this complex quadratically-constrained quadratic program as PhaseCut, in reference to its real counterpart which is known to be equivalent to the classical graph-partition problem MaxCut [13]. These non-convex problems are NP-hard in general, difficult to solve in practice, and have been extensively studied, yielding a number of efficient optimization schemes for particular instances. The most straightforward approach consists in iteratively minimizing (17) with respect to each $u_{i}$ alternatively, which can be done in closed-form [7]. Since the problem is non-convex, this method is bound to converge to a local minimum which depends on the initialization.

3.2 Lifting solution

A particularly popular alternative to solve (17) is referred to as Lifting, and consists in solving the following convex semi-definite program (SDP) instead :

[TABLE]

where $\succeq\mbox{\bf 0}$ denotes positive semi-definiteness. Note that (19) is a relaxation of (17), in the sense that if \hat{\mbox{\bf U}}=\hat{\mbox{\boldmathu}}\hat{\mbox{\boldmathu}}^{H} is a rank-1 solution of (19), then \hat{\mbox{\boldmathu}} is a solution of (17). However, $\hat{\mbox{\bf U}}$ may not always be rank-1 in practice. In the classical prior-less phase retrieval case where \Gamma_{\phi}^{-1}\!=\!\mbox{\boldmath0}, the combined extensive research efforts in [6] and [7] lay theoretical grounds providing conditions on A for which solving (19) enables stable recovery of the phase vector $\phi$ and signal vector $x$ with high probability. Extending these theories to the proposed Bayesian generalization necessitates a deep research investigation, which cannot be tackled within this short paper. Rather, an experimental validation of the lifting approach in the multivariate Von-Mises phase retrieval setting is conducted in Section 4.

3.3 Algorithms

A large number of efficient generic SDP solvers are available, including interior-point methods [14] or augmented Lagrangian methods [15]. As mentioned in [7], the block-coordinate descent (BCD) method proposed in [16] is particularly simple and efficient for problems of the form (19), and is therefore used here. In practice, when the obtained solution $\hat{\mbox{\bf U}}$ is not rank-1, a natural approach consists in selecting the leading eigenvector of U.

4 Experiments

In this section, we propose two different experimental setups to assess the relevance of the above procedure. More precisely, we consider two particular cases of the multivariate Von-Mises prior (4) : the 1D Von-Mises distribution and the Markov chain.

For both setups, we confront it to two state-of-the-art phase retrieval algorithms, namely PhaseCut [7] and prVBEM [9]. The first one relies on the same optimization procedure as the one proposed here, but does not exploit any information on the phases to recover. The second one shares the same Bayesian formulation (2) as the algorithm proposed here but considers a non-informative, uniform distribution on the phases. In the sequel, we will refer to our approach as “informed PhaseCut”.

We consider the following general experimental setup. Observations are generated according to model (2) with $M=256$ and $K=64$. The elements of the dictionary A (resp. vector $x$) are i.i.d. realizations of a zero-mean circular Gaussian distribution with variance $M^{-1}$ (resp. $1$). We assess the performance in terms of the reconstruction of the signal $x$. In particular, we consider the correlation between the estimated signal and the one used to generate the data,

[TABLE]

as a function of the noise variance $\sigma_{n}^{2}$. This figure of merit is evaluated from $50$ trials for each simulation points.

4.1 1D Von Mises prior

As a first experimental setup, we consider the case where the phase noise is distributed on each sensor independently of one another according to the Von Mises law (3) with parameter $\mu_{i}=0$ and $\kappa_{i}=1$, $\forall i\in\{1,\ldots,M\}$.

Figure 1 presents the performance of the three algorithms with this particular prior distribution. As expected, informed PhaseCut outperforms the other algorithms, proving a good inclusion of the prior additional information. More particularly, the gap between them increases with the noise variance: for $\sigma_{n}^{2}=0.6$, informed PhaseCut achieves a correlation around $0.7$ against $0.3$ for PhaseCut and prVBEM.

4.2 Markov chain

In a second experimental setup, we consider the particular case where only the first two subdiagonals of $\Delta$ are non-zero. Considering a small variance of the phases $\theta$, straightforward calculus leads to

[TABLE]

where $\Gamma_{\phi}^{-1}$ is linked to the parameters \mbox{\boldmath\Delta},\mbox{\boldmath\kappa} through (8). This expression can be directly identified to a Markov chain such as $\forall i\in\{2,\ldots,M\}$, $\theta_{i}=a\;\theta_{i-1}+\omega_{i},$ where $\omega_{i}\sim\mathcal{N}(0,\sigma^{2}_{\theta})$, $\forall i$ and $\theta_{1}\sim\mathcal{N}(0,\sigma^{2}_{\theta})$, providing that

[TABLE]

We suppose here that $a=0.8$, $\sigma_{\theta}^{2}=0.1$.

Figure 2 confirms the good behavior of informed PhaseCut observed in the first experiment setup: taking into account the structure of the missing phases, it allows a better estimation (in the sense of the correlation) of the signal of interest $x$. The advantage brought by such prior inclusion increases with the noise variance: informed PhaseCut reveals here again more robustness.

5 Conclusion

In this paper, we have presented a novel algorithm able to solve the phase recovery problem with a multivariate Von Mises prior distribution. To that end, we have showed that this particular prior information can be efficiently integrated into a Maximum A Posteriori estimation by means of a Mahalanobis distance. The proposed solution relies on a lifting procedure and, to the extent of our experiments, reveals a coherent behavior with regard to non-informed state-of-the-art algorithms.

Appendix A Proof of Lemma 1

We have:

[TABLE]

where we have assumed that $(\Gamma_{\phi}^{-1})_{ki}=(\Gamma_{\phi}^{-1})_{ik}$, $\forall(i,k)\in\{1,\ldots,M\}^{2}$, or, in other words, $(\Gamma_{\phi}^{-1})_{ik}\in\mathbb{R}$. Identifying then parameters $\kappa$ and $\Delta$ of the multivariate distribution (4) with (26), it comes straightforwardly that, under the condition (8),

[TABLE]

where $\propto$ denotes here equality up to a constant. This means that we can use indifferently the multivariate Von-Mises distribution (4) or the Mahalanobis distance as a cost function if we add to the latter the constraint $|\phi_{i}|=1$, $\forall i\in\{1,\ldots,M\}$. $\square$