# Fast and Robust Phase Retrieval via Deep Expectation-Consistent
Approximation

Saurav K. Shastri, and Philip Schniter S. K. Shastri and P. Schniter are with
the Dept. of Electrical and Computer Engineering, The Ohio State University,
Columbus, OH 43201, USA, Email<EMAIL_ADDRESS>
###### Abstract

Accurately recovering images from phaseless measurements is a challenging and
long-standing problem. In this work, we present “deepECpr,” which combines
expectation-consistent (EC) approximation with deep denoising networks to
surpass state-of-the-art phase-retrieval methods in both speed and accuracy.
In addition to applying EC in a non-traditional manner, deepECpr includes a
novel stochastic damping scheme that is inspired by recent diffusion methods.
Like existing phase-retrieval methods based on plug-and-play priors,
regularization by denoising, or diffusion, deepECpr iterates a denoising stage
with a measurement-exploitation stage. But unlike existing methods, deepECpr
requires far fewer denoiser calls. We compare deepECpr to the state-of-the-art
prDeep (Metzler et al., 2018), Deep-ITA (Wang et al., 2020), and Diffusion
Posterior Sampling (Chung et al., 2023) methods for noisy phase-retrieval of
color, natural, and unnatural grayscale images on oversampled-Fourier and
coded-diffraction-pattern measurements and find improvements in both PSNR and
SSIM with 5x fewer denoiser calls.

##  I Introduction

In phase retrieval (PR), one aims to recover an unknown signal or image
$\bm{x}_{\mathsf{true}}\in{\mathbb{R}}^{d}$ from phaseless measurements

$\displaystyle\bm{y}=|\bm{Ax}_{\mathsf{true}}|+\bm{w},$ (1)

where $\bm{A}\in{\mathbb{C}}^{m\times d}$ is a known linear operator,
$|\cdot|$ denotes element-wise squared absolute value, and
$\bm{w}\in{\mathbb{R}}^{m}$ is noise. PR is needed when it is infeasible to
measure phase information, as often occurs in optics [1], holography [2],
X-ray crystallography [3, 4], diffraction imaging [5], non-line-of-sight
imaging [6], astronomy [7], computational biology [8], speech recognition [9],
and other fields. PR is challenging because, even with full-rank $\bm{A}$ and
no noise, there exist many hypotheses of $\bm{x}$ that explain the phaseless
measurements $\bm{y}$ [10].

Various approaches to PR have been proposed. Classical methods, like the
Gerchberg-Saxton (GS) [11] and Hybrid Input-Output (HIO) [12] algorithms, are
based on iterative projection. Although these algorithms and their variants
[13, 14, 15] are simple to implement and fast to execute, their output
qualities don’t compete with those of contemporary methods, as we show in the
sequel.

A more modern approach is to formulate PR as negative log-likelihood (NLL)
minimization and solve it using gradient-based iterative methods. The PR NLL
is non-convex, however, and so spectral initialization strategies have been
proposed in an attempt to avoid bad local minima [16, 17, 18]. Spectral
initialization tends to work well with randomized $\bm{A}$ [16, 17, 18] but
often fails for the deterministic Fourier $\bm{A}$ (see, e.g., [19]). As a
more direct attack on non-convexity, convex relaxations such as PhaseLift [20]
and PhaseMax [21] have been proposed (see [22] for an overview). However,
these relaxations manifest as semidefinite programs in $d^{2}$-dimensional
space, which are computationally impractical for imaging applications with
$d$-pixel images. Approximate Message Passing (AMP) algorithms have also been
proposed for PR under sparsity-based signal models [23, 24]. Although near-
optimal for high-dimensional i.i.d. or rotationally invariant random $\bm{A}$
[25, 26], they tend to diverge for deterministic Fourier $\bm{A}$.

For image PR, various approaches have been proposed to exploit the prior
knowledge that $\bm{x}$ is an image. For example, the plug-and-play (PnP) [27]
and RED [28, 29] frameworks have been adapted to PR in [30, 19, 31]. These
methods iterate between NLL reduction and denoising, allowing them to harness
the power of deep-network-based image denoisers such as DnCNN [32], and they
tend to work well with both random and deterministic $\bm{A}$. AMP-based PnP
PR approaches have also been proposed [33], but they tend to struggle with
deterministic Fourier $\bm{A}$ [19]. The Compressed-Sensing Generative Model
(GSGM) framework [34] was applied to PR in [35], where—given an image
generator $g_{\bm{\theta}}(\cdot)$ with fixed parameters $\bm{\theta}$—one
searches for a code vector $\bm{z}$ for which
$|\bm{A}g_{\bm{\theta}}(\bm{z})|$ matches the phaseless measurements $\bm{y}$.
A variation [36, 37] inspired by Deep Image Prior (DIP) [38] instead fixes
$\bm{z}$ and optimizes the generator parameters $\bm{\theta}$. End-to-end deep
networks have also been proposed to learn the inverse mapping from $\bm{y}$ to
$\bm{x}$ in PR [39, 40], but they often fail with deterministic $\bm{A}$ [31].
More recently, diffusion methods have been proposed to sample from the
posterior distribution [41, 42] in PR. As we show in Sec. IV, however, DPS
[41] PR suffers from instability and image artifacts.

In this work, we propose a novel approach to PR based on the expectation-
consistent (EC) approximation framework [43, 44], which we refer to as
“deepECpr.” Different from existing EC-based methods like VAMP [45] and GVAMP
[24], deepECpr avoids the need for large random $\bm{A}$ by employing EC in a
different way. In addition, deepECpr includes a novel “stochastic damping”
scheme to further mitigate potential deviations from EC modeling assumptions.
Like PnP, RED, and diffusion-based approaches to PR, deepECpr iteratively
calls a deep-net denoiser, but it converges in significantly fewer iterations
to more accurate solutions. We attribute deepECpr’s excellent performance to
its ability to track error-variances, thereby allowing efficient use of the
denoiser. This paper builds on our conference publication [46] by adding
stochastic damping, detailed explanations and derivations, and greatly
expanded experimental results.

The paper is organized as follows. Section II reviews generalized linear
models and EC, Sec. III details our proposed approach, Sec. IV presents
numerical experiments, and Sec. V concludes.

##  II Background

###  II-A Generalized Linear Model (GLM)

A wide variety of linear and nonlinear inverse problems fall under the
Generalized Linear Model (GLM) framework [47], which models the relationship
between the measurements $\bm{y}\in\mathcal{Y}^{m}$ and the signal/image
hypothesis $\bm{x}\in{\mathbb{R}}^{d}$ or ${\mathbb{C}}^{d}$ using

$\displaystyle
p_{\mathsf{y}|\mathsf{x}}(\bm{y}|\bm{x})=\prod_{i=1}^{m}p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i})\text{~{}for~{}}\bm{z}=\bm{Ax},$
(2)

for scalar “measurement channel” $p_{\mathsf{y}|\mathsf{z}}$ and forward
transform $\bm{A}\in{\mathbb{R}}^{m\times d}$ or ${\mathbb{C}}^{m\times d}$.
By choosing appropriate $p_{\mathsf{y}|\mathsf{z}}$, one can model, e.g.,
additive noise of an arbitrary distribution, logistic regression [48], Poisson
regression [49], dequantization [50], and phase retrieval [10, 51].

$p_{\mathsf{x}}(\cdot)$$\bm{A}$$\prod_{i}p_{\mathsf{y}|\mathsf{z}}(\cdot|z_{i})$
$\bm{y}$ $\bm{x}$$\bm{z}$ Figure 1: Generalized linear model relating signal
$\bm{x}$ to measurements $\bm{y}$.

A simple instance of the GLM is the Standard Linear Model (SLM), where
$p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i})=\mathcal{N}(y_{i};z_{i},v)$ and so

$\bm{y}=\bm{z}_{\mathsf{true}}+\bm{w}=\bm{Ax}_{\mathsf{true}}+\bm{w},\quad\bm{w}\sim\mathcal{N}(\bm{0},v\bm{I}).$
(3)

Here, $v$ is the variance of the additive white Gaussian noise (AWGN) vector
$\bm{w}$. Through appropriate choice of $\bm{A}$, the SLM covers, e.g.,
denoising, deblurring, compressed sensing, super-resolution, inpainting,
magnetic resonance imaging, and computed tomography.

For phase retrieval, several variants of $p_{\mathsf{y}|\mathsf{z}}$ have been
employed, such as
$p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i})=\mathcal{N}(y_{i},|z_{i}|^{2},v)$
[16],
$p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i})=\mathrm{Poisson}(y_{i};|z_{i}|^{2})$
[52], the Rician model resulting from $y_{i}\sim|z_{i}+w_{i}|$ with
$w_{i}\sim\mathcal{N}(0,v)$ [23], and

$\displaystyle
p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i})=\mathcal{N}(y_{i};|z_{i}|,v).$ (4)

For algorithm design, we will focus on (4), as it has proven to be effective
in practice [53] and is adopted by state-of-the-art techniques like prDeep
[19] and Deep-ITA [31].

###  II-B Expectation-Consistent Approximation

We now consider statistical estimation of $\bm{x}\sim p_{\mathsf{x}}$ from
$\bm{y}$ under the GLM (2). When both $p_{\mathsf{x}}$ and
$p_{\mathsf{y}|\mathsf{z}}$ are Gaussian, the maximum a posteriori (MAP) and
minimum mean-squared error (MMSE) estimates of $\bm{x}$ coincide and are
analytically computable. But Gaussian $p_{\mathsf{y}|\mathsf{z}}$ is not
useful for phase retrieval. When both $p_{\mathsf{x}}$ and
$p_{\mathsf{y}|\mathsf{z}}$ are log-concave, the posterior
$p_{\mathsf{x}|\mathsf{y}}$ is also log-concave and thus the MAP estimate of
$\bm{x}$ can be computed using standard convex-optimization algorithms [54,
55]. Although the $p_{\mathsf{y}|\mathsf{z}}$ in (4) is indeed log-concave,
high-fidelity image priors $p_{\mathsf{x}}$ are not. Furthermore, the
widespread use of PSNR as a recovery-performance metric suggests the use of
MMSE estimation over MAP.

The expectation-consistent (EC) approximation framework [43, 44] is a well
established method to approximate the MMSE estimate of $\bm{x}$ from $\bm{y}$,
i.e., the conditional mean $\operatorname{E}\\{\bm{x}|\bm{y}\\}$. The
application of EC to the SLM (3) or GLM (2) is sometimes referred to as vector
AMP (VAMP) [56] and generalized VAMP [24], respectively. Rigorous analyses of
VAMP have established its MSE optimality under high-dimensional rotationally
invariant $\bm{A}$ and various assumptions on $p_{\mathsf{x}}$ and
$p_{\mathsf{y}|\mathsf{z}}$ [45, 57, 58, 59].

In the above works, EC is employed to MMSE-estimate $\bm{x}$ from $\bm{y}$
using the prior $p_{\mathsf{x}}$ and the likelihood
$p_{\mathsf{y}|\mathsf{x}}$ from (2). This approach, however, fails for many
deterministic choices of $\bm{A}$. We instead employ EC to MMSE-estimate
$\bm{z}$ from $\bm{y}$, using the prior $p_{\mathsf{z}}$ and likelihood
$p_{\mathsf{y}|\mathsf{z}}$. Separately, we estimate $\bm{x}$ from
$\bm{z}=\bm{Ax}$. Thus, when reviewing EC below, we do so in the context of
estimating $\bm{z}$ using $p_{\mathsf{z}}$ and $p_{\mathsf{y}|\mathsf{z}}$.
The more “traditional” application of EC (estimating $\bm{x}$ using
$p_{\mathsf{x}}$ and $p_{\mathsf{y}|\mathsf{x}}$) is identical up to variable
substitutions.

To begin our review of EC, notice that one can write the true posterior
distribution $p_{\mathsf{z}|\mathsf{y}}(\cdot|\bm{y})$ without approximation
as

$\displaystyle
p_{\mathsf{z}|\mathsf{y}}(\cdot|\bm{y})=\operatorname*{arg\,min}_{q}D\big{(}q\|p_{\mathsf{z}|\mathsf{y}}(\cdot|\bm{y})\big{)}$
(5)
$\displaystyle=\operatorname*{arg\,min}_{q}D\big{(}q\|p_{\mathsf{y}|\mathsf{z}}(\bm{y}|\cdot)\big{)}+D(q\|p_{\mathsf{z}})+H(q)$
(6)
$\displaystyle=\operatorname*{arg\,min}_{q_{1}=q_{2}=q_{3}}\underbrace{D\big{(}q_{1}\|p_{\mathsf{y}|\mathsf{z}}(\bm{y}|\cdot)\big{)}+D(q_{2}\|p_{\mathsf{z}})+H(q_{3})}_{\displaystyle\triangleq
J_{\text{\sf Gibbs}}(q_{1},q_{2},q_{3})},$ (7)

where $D(q\|p)\triangleq\int
q(\bm{z})\ln\frac{q(\bm{z})}{p(\bm{z})}\mathop{}\\!\mathrm{d}\bm{z}$ denotes
Kullback-Leibler (KL) divergence, $H(q)\triangleq-\int q(\bm{z})\ln
q(\bm{z})\mathop{}\\!\mathrm{d}\bm{z}$ denotes differential entropy, and Bayes
rule
$p_{\mathsf{z}|\mathsf{y}}(\bm{z}|\bm{y})=p_{\mathsf{y}|\mathsf{z}}(\bm{y}|\bm{z})p_{\mathsf{z}}(\bm{z})/p_{\mathsf{y}}(\bm{y})$
was used for (6). Because the optimization in (7) is intractable, Opper and
Winther [60] proposed to relax the equality constraints to moment-matching
constraints, i.e.,

$\displaystyle\operatorname*{arg\,min}_{q_{1},q_{2},q_{3}}J_{\text{\sf
Gibbs}}(q_{1},q_{2},q_{3})\text{~{}such that~{}}$ (8)
$\displaystyle~{}\begin{cases}\operatorname{E}\\{\bm{z}|q_{1}\\}=\operatorname{E}\\{\bm{z}|q_{2}\\}=\operatorname{E}\\{\bm{z}|q_{3}\\}\triangleq\hat{\bm{z}}\\\
\operatorname{tr}(\operatorname{Cov}\\{\bm{z}|q_{1}\\})=\operatorname{tr}(\operatorname{Cov}\\{\bm{z}|q_{2}\\})=\operatorname{tr}(\operatorname{Cov}\\{\bm{z}|q_{3}\\})\triangleq
m\widehat{v},\end{cases}$

where $\operatorname{E}\\{\bm{z}|q_{i}\\}$ and
$\operatorname{Cov}\\{\bm{z}|q_{i}\\}$ denote the mean and covariance of
$\bm{z}$ under $\bm{z}\sim q_{i}$ and $\operatorname{tr}(\cdot)$ denotes the
trace, and then approximate $\operatorname{E}\\{\bm{z}|\bm{y}\\}$ by
$\hat{\bm{z}}$. The solution to (8) takes the form

$\displaystyle
q_{1}(\bm{z};\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)})$
$\displaystyle\propto
p_{\mathsf{y}|\mathsf{z}}(\bm{y}|\bm{z})\mathcal{N}(\bm{z};\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)}\bm{I})~{}$
(9) $\displaystyle
q_{2}(\bm{z};\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})$
$\displaystyle\propto
p_{\mathsf{z}}(\bm{z})\mathcal{N}(\bm{\bar{z}}^{\scriptscriptstyle(2)};\bm{z},\overline{v}^{\scriptscriptstyle(2)}\bm{I})~{}$
(10) $\displaystyle q_{3}(\bm{z};\hat{\bm{z}},\widehat{v})$
$\displaystyle=\mathcal{N}(\bm{z};\hat{\bm{z}},\widehat{v}\bm{I})$ (11)

for some
$\bm{\bar{z}}^{\scriptscriptstyle(1)},\bm{\bar{z}}^{\scriptscriptstyle(2)}$
and
$\overline{v}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(2)}$
known as the “extrinsic” means and variances, respectively. Furthermore,

$\displaystyle\widehat{v}$
$\displaystyle=\big{(}1/\overline{v}^{\scriptscriptstyle(1)}+1/\overline{v}^{\scriptscriptstyle(2)}\big{)}^{-1}$
(12) $\displaystyle\hat{\bm{z}}$
$\displaystyle=\bigg{(}\frac{\bm{\bar{z}}^{\scriptscriptstyle(1)}}{\overline{v}^{\scriptscriptstyle(1)}}+\frac{\bm{\bar{z}}^{\scriptscriptstyle(2)}}{\overline{v}^{\scriptscriptstyle(2)}}\bigg{)}\widehat{v}.$
(13)

Note that $q_{1}$ combines the true likelihood $p_{\mathsf{y}|\mathsf{z}}$
with the Gaussian “pseudo-prior”
$\bm{z}\sim\mathcal{N}(\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)}\bm{I})$,
while $q_{2}$ combines the true prior $p_{\mathsf{z}}$ with the Gaussian
“pseudo-likelihood”
$\mathcal{N}(\bm{\bar{z}}^{\scriptscriptstyle(2)};\bm{z},\overline{v}^{\scriptscriptstyle(2)}\bm{I})$,
i.e., the model

$\displaystyle\bm{\bar{z}}^{\scriptscriptstyle(2)}=\bm{z}_{\mathsf{true}}+\bm{e}\text{~{}with~{}}\bm{e}\sim\mathcal{N}(\bm{0},\overline{v}^{\scriptscriptstyle(2)}\bm{I}).$
(14)

To solve for
$(\bm{\bar{z}}^{\scriptscriptstyle(1)},\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(2)})$,
Opper and Winther [60] proposed Alg. 1, which iterates the following four
steps: compute the posterior mean and trace-covariance of
$q_{2}(\cdot;\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})$,
update the extrinsic quantities
$(\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)})$
to obey (12)-(13), compute the posterior mean and trace-covariance of
$q_{1}(\cdot;\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)})$,
and update the extrinsic quantities
$(\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})$
to obey (12)-(13). The non-informative initialization of
$\overline{v}^{\scriptscriptstyle(2)}$ and
$\bm{\bar{z}}^{\scriptscriptstyle(2)}$ in line 4 ensures that
$q_{2}=p_{\mathsf{z}}$ initially, so that the initial
$\widehat{v}^{\scriptscriptstyle(2)}$ and
$\hat{\bm{z}}^{\scriptscriptstyle(2)}$ are the prior variance and mean (as are
the initial $\overline{v}^{\scriptscriptstyle(1)}$ and
$\bm{\bar{z}}^{\scriptscriptstyle(1)}$). Algorithm 1 can be interpreted [60]
as an instance of expectation propagation (EP) from [61].

Algorithm 1 EC to infer $\bm{z}\sim p_{\mathsf{z}}$ from $\bm{y}\sim
p_{\mathsf{y}|\mathsf{z}}(\cdot|\bm{z})$

1:$p_{\mathsf{z}}(\cdot),p_{\mathsf{y}|\mathsf{z}}(\bm{y}|\cdot)$

2:$q_{1}(\bm{z};\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)})\propto
p_{\mathsf{y}|\mathsf{z}}(\bm{y}|\bm{z})\mathcal{N}(\bm{z};\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)}\bm{I})$

3:$q_{2}(\bm{z};\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})\propto
p_{\mathsf{z}}(\bm{z})\mathcal{N}(\bm{z};\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)}\bm{I})$

4:$\overline{v}^{\scriptscriptstyle(2)}\leftarrow\infty$ and
$\bm{\bar{z}}^{\scriptscriptstyle(2)}\leftarrow\text{arbitrary finite vector
in }{\mathbb{R}}^{m}$

5:repeat

6: // Exploit prior $p_{\mathsf{z}}(\cdot)$ and pseudo-likelihood

7:
$\widehat{v}^{\scriptscriptstyle(2)}\leftarrow\tfrac{1}{m}\operatorname{tr}(\operatorname{Cov}\\{\bm{z}|q_{2}(\cdot;\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})\\})$

8:
$\hat{\bm{z}}^{\scriptscriptstyle(2)}\leftarrow\operatorname{E}\\{\bm{z}|q_{2}(\cdot;\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})\\}$

9:
$\overline{v}^{\scriptscriptstyle(1)}\leftarrow(1/\widehat{v}^{\scriptscriptstyle(2)}-1/\overline{v}^{\scriptscriptstyle(2)})^{-1}$

10:
$\bm{\bar{z}}^{\scriptscriptstyle(1)}\leftarrow(\hat{\bm{z}}^{\scriptscriptstyle(2)}/\widehat{v}^{\scriptscriptstyle(2)}-\bm{\bar{z}}^{\scriptscriptstyle(2)}/\overline{v}^{\scriptscriptstyle(2)})\overline{v}^{\scriptscriptstyle(1)}$

11: // Exploit likelihood $p_{\mathsf{y}|\mathsf{z}}(\bm{y}|\cdot)$ and
pseudo-prior

12:
$\widehat{v}^{\scriptscriptstyle(1)}\leftarrow\tfrac{1}{m}\operatorname{tr}(\operatorname{Cov}\\{\bm{z}|q_{1}(\cdot;\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)})\\})$

13:
$\hat{\bm{z}}^{\scriptscriptstyle(1)}\leftarrow\operatorname{E}\\{\bm{z}|q_{1}(\cdot;\bm{\bar{z}}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)})\\}$

14:
$\overline{v}^{\scriptscriptstyle(2)}\leftarrow(1/\widehat{v}^{\scriptscriptstyle(1)}-1/\overline{v}^{\scriptscriptstyle(1)})^{-1}$

15:
$\bm{\bar{z}}^{\scriptscriptstyle(2)}\leftarrow(\hat{\bm{z}}^{\scriptscriptstyle(1)}/\widehat{v}^{\scriptscriptstyle(1)}-\bm{\bar{z}}^{\scriptscriptstyle(1)}/\overline{v}^{\scriptscriptstyle(1)})\overline{v}^{\scriptscriptstyle(2)}$

16:return $\hat{\bm{z}}^{\scriptscriptstyle(2)}$ as the EC approximation of
$\operatorname{E}\\{\bm{z}|\bm{y}\\}$

##  III Proposed Method

In the following subsections we propose several modifications of the standard
EC approach from Alg. 1 that involve stochastic damping, deep networks, and
simplifications/approximations specific to phase retrieval.

###  III-A Stochastic Damping

The EC Alg. 1 is not guaranteed to converge for general $p_{\mathsf{z}}$ and
$p_{\mathsf{y}|\mathsf{z}}$. One technique that helps promote convergence is
damping [62, 63], which slows the updates in a way that preserves the EC fixed
points. For Gaussian $p_{\mathsf{z}}$ and $p_{\mathsf{y}|\mathsf{z}}$,
appropriate damping guarantees convergence [44].

The damping technique from [63] considers the
$\overline{v}^{\scriptscriptstyle(i)}$ and
$\bm{\bar{z}}^{\scriptscriptstyle(i)}$ computed in lines 9-10 and 14-15 of
Alg. 1 as “raw” quantities and damps them using the additional steps

$\displaystyle\overline{v}^{\scriptscriptstyle(i)}$
$\displaystyle\leftarrow\big{(}\mu^{\scriptscriptstyle(i)}\sqrt{\smash{\overline{v}^{\scriptscriptstyle(i)}_{\mathsf{raw}}}\vphantom{\overline{v}^{\scriptscriptstyle(i)}}}+(1-\mu^{\scriptscriptstyle(i)})\sqrt{\smash{\overline{v}^{\scriptscriptstyle(i)}_{\mathsf{old}}}\vphantom{\overline{v}^{\scriptscriptstyle(i)}}}\big{)}^{2}$
(15a) $\displaystyle\bm{\bar{z}}^{\scriptscriptstyle(i)}$
$\displaystyle\leftarrow\mu^{\scriptscriptstyle(i)}\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(i)}+(1-\mu^{\scriptscriptstyle(i)})\bm{\bar{z}}_{\mathsf{old}}^{\scriptscriptstyle(i)},$
(15b)

where $\overline{v}_{\mathsf{old}}^{\scriptscriptstyle(i)}$ and
$\bm{\bar{z}}_{\mathsf{old}}^{\scriptscriptstyle(i)}$ denote the values of
$\overline{v}^{\scriptscriptstyle(i)}$ and
$\bm{\bar{z}}^{\scriptscriptstyle(i)}$ from the previous iteration and
$\mu^{\scriptscriptstyle(i)}\in(0,1]$ is a fixed damping constant. When
$\mu^{\scriptscriptstyle(1)}=1=\mu^{\scriptscriptstyle(2)}$, the original EC
is recovered.

We now propose a “stochastic damping” procedure inspired by diffusion methods
like [64]. It damps $\overline{v}^{\scriptscriptstyle(2)}$ as in (15a) but
constructs $\bm{\bar{z}}^{\scriptscriptstyle(2)}$ by adding AWGN of a
prescribed variance:

$\displaystyle\overline{v}^{\scriptscriptstyle(2)}$
$\displaystyle\leftarrow\big{(}\mu^{\scriptscriptstyle(2)}\sqrt{\smash{\overline{v}^{\scriptscriptstyle(2)}_{\mathsf{raw}}}\vphantom{\overline{v}^{\scriptscriptstyle(2)}}}+(1-\mu^{\scriptscriptstyle(2)})\sqrt{\smash{\overline{v}^{\scriptscriptstyle(2)}_{\mathsf{old}}}\vphantom{\overline{v}^{\scriptscriptstyle(2)}}}\big{)}^{2}$
(16a) $\displaystyle\bm{\bar{z}}^{\scriptscriptstyle(2)}$
$\displaystyle\leftarrow\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(2)}+\bm{\epsilon},\quad\bm{\epsilon}\sim\mathcal{N}(\bm{0},[\overline{v}^{\scriptscriptstyle(2)}-\overline{v}^{\scriptscriptstyle(2)}_{\mathsf{raw}}]\bm{I}).$
(16b)

To better understand (16b), recall from (14) that, under ideal conditions,
$\bm{\bar{z}}^{\scriptscriptstyle(2)}=\bm{z}_{\mathsf{true}}+\bm{e}$ for
$\bm{e}\sim\mathcal{N}(\bm{0},\overline{v}^{\scriptscriptstyle(2)}\bm{I})$ and
$\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(2)}=\bm{z}_{\mathsf{true}}+\bm{\varepsilon}$
for
$\bm{\varepsilon}\sim\mathcal{N}(\bm{0},\overline{v}_{\mathsf{raw}}^{\scriptscriptstyle(2)}\bm{I})$.
This implies that a valid $\bm{\bar{z}}^{\scriptscriptstyle(2)}$ can be
constructed from $\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(2)}$ by
adding additional AWGN of variance
$\overline{v}^{\scriptscriptstyle(2)}-\overline{v}_{\mathsf{raw}}^{\scriptscriptstyle(2)}$,
as done in (16b). Here, it is expected that
$\overline{v}^{\scriptscriptstyle(2)}>\overline{v}^{\scriptscriptstyle(2)}_{\mathsf{raw}}$
because the variances decrease over the iterations and damping aims to slow
down that decrease. In Sec. IV, we show that (16) has advantages over (15),
perhaps because (16) attempts to enforce the EC model (14).

###  III-B deepEC for general $\bm{A}$

As formulated in Alg. 1, EC estimates $\bm{z}\sim p_{\mathsf{z}}$ from
observations $\bm{y}\sim p_{\mathsf{y}|\mathsf{z}}(\cdot|\bm{z})$. To solve
inverse problems using the GLM model (2), we instead want to estimate
$\bm{x}\sim p_{\mathsf{x}}$ from $\bm{y}\sim
p_{\mathsf{y}|\mathsf{z}}(\cdot|\bm{z})$ where $\bm{z}=\bm{Ax}$. In principle,
this can be done using a minor modification of Alg. 1 that rewrites
$\hat{\bm{z}}^{\scriptscriptstyle(2)}$ in line 8 as a function of an MMSE
estimate of $\bm{x}$:

$\displaystyle\hat{\bm{z}}^{\scriptscriptstyle(2)}=\operatorname{E}\\{\bm{z}|q_{2}(\cdot;\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})\\}$
(17) $\displaystyle=$
$\displaystyle\operatorname{E}\\{\bm{z}|\bm{\bar{z}}^{\scriptscriptstyle(2)}=\bm{z}+\bm{e}\\},~{}\bm{z}\sim
p_{\mathsf{z}},~{}\bm{e}\sim\mathcal{N}(0,\overline{v}^{\scriptscriptstyle(2)}\bm{I})$
(18) $\displaystyle=$
$\displaystyle\bm{A}\underbrace{\operatorname{E}\\{\bm{x}|\bm{\bar{z}}^{\scriptscriptstyle(2)}=\bm{Ax}+\bm{e}\\},~{}\bm{x}\sim
p_{\mathsf{x}},~{}\bm{e}\sim\mathcal{N}(0,\overline{v}^{\scriptscriptstyle(2)}\bm{I})}_{\displaystyle\triangleq\hat{\bm{x}}^{\scriptscriptstyle(2)}},\qquad$
(19)

and then returns $\hat{\bm{x}}^{\scriptscriptstyle(2)}$ as the EC output. The
practical challenge, however, is that computing
$\hat{\bm{x}}^{\scriptscriptstyle(2)}$ via (19) involves MMSE estimation of
$\bm{x}$ under the SLM (3), which is itself non-trivial.

Fortunately, it is now commonplace to train a deep networks to solve SLMs
[65]. Writing
$\hat{\bm{x}}^{\scriptscriptstyle(2)}\approx\bm{f}_{\bm{\theta}}(\bm{\bar{z}}^{\scriptscriptstyle(2)};\overline{v}^{\scriptscriptstyle(2)})$
for deep network $\bm{f}_{\bm{\theta}}$, we could train the network parameters
$\bm{\theta}$ via

$\displaystyle\arg\min_{\bm{\theta}}\sum_{t=1}^{T}\operatorname{E}\big{\\{}\|\bm{x}_{t}-\bm{f}_{\bm{\theta}}(\bm{Ax}_{t}+\bm{e};\overline{v}^{\scriptscriptstyle(2)})\|^{2}\big{\\}}$
(20)

with training data $\\{\bm{x}_{t}\\}_{t=1}^{T}$, random noise
$\bm{e}\sim\mathcal{N}(\bm{0},\overline{v}^{\scriptscriptstyle(2)}\bm{I})$,
and random noise variance
$\overline{v}^{\scriptscriptstyle(2)}\sim\mathrm{Unif}[0,v_{\max}]$ for some
$v_{\max}$. Likewise, we could train another deep network $h_{\bm{\phi}}$ to
approximate the posterior variance $\widehat{v}^{\scriptscriptstyle(2)}$.
Writing $\widehat{v}^{\scriptscriptstyle(2)}\approx
h_{\bm{\phi}}(\bm{\bar{z}}^{\scriptscriptstyle(2)};\overline{v}^{\scriptscriptstyle(2)})$,
its network parameters $\bm{\phi}$ might be trained via

$\arg\min_{\bm{\phi}}\sum_{t=1}^{T}\operatorname{E}\big{\\{}\big{|}\tfrac{1}{m}\|\bm{z}_{t}-\bm{A}\hat{\bm{x}}_{t}^{\scriptscriptstyle(2)}\|^{2}-h_{\bm{\phi}}(\bm{Ax}_{t}+\bm{e};\overline{v}^{\scriptscriptstyle(2)})\big{|}\big{\\}},$
(21)

where
$\hat{\bm{x}}_{t}^{\scriptscriptstyle(2)}\triangleq\bm{f}_{\bm{\theta}}(\bm{Ax}_{t}+\bm{e};\overline{v}^{\scriptscriptstyle(2)})$
and where $\bm{x}_{t}$, $\bm{e}$, and $\overline{v}^{\scriptscriptstyle(2)}$
are constructed as they were for (20).

We will refer to EC with these deep-network approximations as “deepEC.” To
summarize, deepEC is EC Algorithm 1 but with
$\hat{\bm{z}}^{\scriptscriptstyle(2)}$ computed via
$\bm{Af}_{\bm{\theta}}(\bm{\bar{z}}^{\scriptscriptstyle(2)};\overline{v}^{\scriptscriptstyle(2)})$
and $\widehat{v}^{\scriptscriptstyle(2)}$ computed via
$h_{\bm{\phi}}(\bm{\bar{z}}^{\scriptscriptstyle(2)};\overline{v}^{\scriptscriptstyle(2)})$.
As such, deepEC solves GLM problems by plugging deep SLM networks into EC Alg.
1.

Algorithm 2 deepECpr

1:transform $\bm{A}$, channel $p_{\mathsf{y|z}}(\bm{y},\cdot)$, denoiser
$\bm{f}_{\bm{\theta}}(\cdot,\cdot)$, denoising factor $\beta\in(0,1)$, image
initialization $\hat{\bm{x}}_{\mathsf{init}}$, variance initialization
$\overline{v}_{\mathsf{init}}$, variance initialization factor $\zeta\approx
1.2$, damping factors $\mu^{\scriptscriptstyle(1)}\in(0,1]$ and
$\mu^{\scriptscriptstyle(2)}\in(0,1]$

2:initialize:
$\bm{\bar{z}}^{\scriptscriptstyle(2)}\leftarrow\bm{A}\hat{\bm{x}}_{\mathsf{init}}+\bm{n}$
with $\bm{n}\sim\mathcal{CN}(\bm{0},\overline{v}_{\mathsf{init}}\bm{I})$, and
$\overline{v}^{\scriptscriptstyle(2)}\leftarrow\zeta\overline{v}_{\mathsf{init}}$

3:for $j=1,\dots,J$ do

4: // Exploit prior $\bm{x}\sim p_{\mathsf{x}}$ (via denoiser
$\bm{f}_{\bm{\theta}}$) and pseudo-likelihood
$p(\bm{\bar{z}}^{\scriptscriptstyle(2)}|\bm{x})=\mathcal{CN}(\bm{\bar{z}}^{\scriptscriptstyle(2)};\bm{Ax},\overline{v}^{\scriptscriptstyle(2)}\bm{I})$

5: denoise:
$\hat{\bm{x}}^{\scriptscriptstyle(2)}\leftarrow\bm{f}_{\bm{\theta}}(\Re\\{\bm{A}^{\textsf{H}}\bm{\bar{z}}^{\scriptscriptstyle(2)}\\},0.5\overline{v}^{\scriptscriptstyle(2)})$

6: deep-network approximated posterior mean and variance:
$\hat{\bm{z}}^{\scriptscriptstyle(2)}\leftarrow\bm{A}\hat{\bm{x}}^{\scriptscriptstyle(2)}$
and
$\widehat{v}^{\scriptscriptstyle(2)}\leftarrow\beta\overline{v}^{\scriptscriptstyle(2)}$

7: extrinsic variance:
$\overline{v}_{\mathsf{raw}}^{\scriptscriptstyle(1)}\leftarrow(1/\widehat{v}^{\scriptscriptstyle(2)}-1/\overline{v}^{\scriptscriptstyle(2)})^{-1}$
and
$\overline{v}_{\mathsf{old}}^{\scriptscriptstyle(1)}\leftarrow\overline{v}^{\scriptscriptstyle(1)}$

8: extrinsic mean:
$\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(1)}\leftarrow(\hat{\bm{z}}^{\scriptscriptstyle(2)}/\widehat{v}^{\scriptscriptstyle(2)}-\bm{\bar{z}}^{\scriptscriptstyle(2)}/\overline{v}^{\scriptscriptstyle(2)})\overline{v}^{\scriptscriptstyle(1)}$
and
$\bm{\bar{z}}_{\mathsf{old}}^{\scriptscriptstyle(1)}\leftarrow\bm{\bar{z}}^{\scriptscriptstyle(1)}$

9: damping:
${\overline{v}^{\scriptscriptstyle(1)}}\leftarrow\big{(}\mu^{\scriptscriptstyle(1)}\sqrt{{\overline{v}_{\mathsf{raw}}^{\scriptscriptstyle(1)}}}+(1-\mu^{\scriptscriptstyle(1)})\sqrt{\smash{\overline{v}^{\scriptscriptstyle(1)}_{\mathsf{old}}}\vphantom{\overline{v}^{\scriptscriptstyle(1)}}}\big{)}^{2}$
and
${\bm{\bar{z}}^{\scriptscriptstyle(1)}}\leftarrow\mu^{\scriptscriptstyle(1)}{\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(1)}}+(1-\mu^{\scriptscriptstyle(1)}){\bm{\bar{z}}_{\mathsf{old}}^{\scriptscriptstyle(1)}}$

10: // Exploit likelihood $p_{\mathsf{y|z}}(y_{i}|z_{i})$ and pseudo-prior
$z_{i}\sim\mathcal{CN}(\overline{z}_{i}^{\scriptscriptstyle(1)},\overline{v}^{\scriptscriptstyle(1)})$
for $i=1,\dots,m$

11: Laplace-approximated posterior mean and variance:
$(\widehat{z}_{i}^{\scriptscriptstyle(1)},\widehat{v}_{i}^{\scriptscriptstyle(1)})~{}\forall
i$ and set
$\widehat{v}^{\scriptscriptstyle(1)}\leftarrow\tfrac{1}{m}\sum_{i=1}^{m}\widehat{v}_{i}^{\scriptscriptstyle(1)}$

12: extrinsic variance:
$\overline{v}_{\mathsf{raw}}^{\scriptscriptstyle(2)}\leftarrow(1/\widehat{v}^{\scriptscriptstyle(1)}-1/\overline{v}^{\scriptscriptstyle(1)})^{-1}$
and
$\overline{v}_{\mathsf{old}}^{\scriptscriptstyle(2)}\leftarrow\overline{v}^{\scriptscriptstyle(2)}$

13: extrinsic mean:
$\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(2)}\leftarrow(\hat{\bm{z}}^{\scriptscriptstyle(1)}/\widehat{v}^{\scriptscriptstyle(1)}-\bm{\bar{z}}^{\scriptscriptstyle(1)}/\overline{v}^{\scriptscriptstyle(1)})\overline{v}^{\scriptscriptstyle(2)}$

14: stochastic damping:
${\overline{v}^{\scriptscriptstyle(2)}}\leftarrow\big{(}\mu^{\scriptscriptstyle(2)}\sqrt{\smash{{\overline{v}}^{\scriptscriptstyle(2)}_{\mathsf{raw}}}\vphantom{{\overline{v}}^{\scriptscriptstyle(2)}}}+(1-\mu^{\scriptscriptstyle(2)})\sqrt{\smash{\overline{v}^{\scriptscriptstyle(2)}_{\mathsf{old}}}\vphantom{\overline{v}^{\scriptscriptstyle(2)}}}\big{)}^{2}$
and
${\bm{\bar{z}}^{\scriptscriptstyle(2)}}\leftarrow\bm{\bar{z}}_{\mathsf{raw}}^{\scriptscriptstyle(2)}+\bm{\epsilon}$
with
$\bm{\epsilon}\sim\mathcal{CN}(\bm{0},[{\overline{v}^{\scriptscriptstyle(2)}}-\overline{v}_{\mathsf{raw}}^{\scriptscriptstyle(2)}]\bm{I})$

15:end for

16:return $\hat{\bm{x}}^{\scriptscriptstyle(2)}$

###  III-C deepEC when $\bm{A}^{\textsf{H}}\bm{A}=\bm{I}$

One practical drawback to deepEC is that the deep networks
$\bm{f}_{\bm{\theta}}$ and $h_{\bm{\phi}}$ are dependent on the forward
operator $\bm{A}$, and so the networks must be retrained for each new choice
of $\bm{A}$. There is, however, an important exception to this rule, which is
when $\bm{A}^{\textsf{H}}\bm{A}=\bm{I}$. For such column-orthogonal $\bm{A}$,
a sufficient statistic [66] for the estimation of
$\hat{\bm{x}}^{\scriptscriptstyle(2)}$ in (19) is

$\displaystyle\bm{\bar{r}}^{\scriptscriptstyle(2)}\triangleq\bm{A}^{\textsf{H}}\bm{\bar{z}}^{\scriptscriptstyle(2)}=\bm{x}+\bm{\epsilon},~{}~{}\bm{\epsilon}\sim\mathcal{N}(\bm{0},\overline{v}^{\scriptscriptstyle(2)}\bm{I}),$
(22)

in which case $\hat{\bm{x}}^{\scriptscriptstyle(2)}$ solves the MMSE denoising
problem

$\displaystyle\hat{\bm{x}}^{\scriptscriptstyle(2)}=\operatorname{E}\\{\bm{x}|\bm{\bar{r}}^{\scriptscriptstyle(2)}=\bm{x}+\bm{\epsilon}\\},~{}\bm{x}\sim
p_{\mathsf{x}},~{}\bm{\epsilon}\sim\mathcal{N}(\bm{0},\overline{v}^{\scriptscriptstyle(2)}\bm{I}),$
(23)

and we can approximate $\hat{\bm{x}}^{\scriptscriptstyle(2)}$ using a deep
_denoising_ network
$\bm{f}_{\bm{\theta}}(\bm{\bar{r}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})$
that is invariant to $\bm{A}$. Likewise, we have that

$\displaystyle\widehat{v}^{\scriptscriptstyle(2)}$
$\displaystyle=\tfrac{1}{m}\operatorname{tr}(\operatorname{Cov}\\{\bm{z}|\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)}\\})$
(24)
$\displaystyle=\tfrac{1}{m}\operatorname{tr}(\bm{A}\operatorname{Cov}\\{\bm{x}|\bm{\bar{z}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)}\\}\bm{A}^{\textsf{H}})$
(25)
$\displaystyle=\tfrac{1}{m}\operatorname{tr}(\bm{A}\operatorname{Cov}\\{\bm{x}|\bm{\bar{r}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)}\\}\bm{A}^{\textsf{H}})$
(26)
$\displaystyle=\tfrac{1}{m}\operatorname{tr}(\bm{A}^{\textsf{H}}\bm{A}\operatorname{Cov}\\{\bm{x}|\bm{\bar{r}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)}\\})$
(27)
$\displaystyle=\tfrac{1}{m}\operatorname{tr}(\operatorname{Cov}\\{\bm{x}|\bm{\bar{r}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)}\\}),$
(28)

and so the network $h_{\bm{\phi}}$ that deepEC uses to approximate
$\widehat{v}^{\scriptscriptstyle(2)}$ can also be invariant to $\bm{A}$. One
way to construct $h_{\bm{\phi}}$ would be to use the Monte-Carlo scheme from
[67], which employs random $\\{\bm{p}_{c}\\}_{c=1}^{C}$ such that
$\operatorname{E}\\{\bm{p}_{c}\bm{p}_{c}^{\top}\\}=\bm{I}$ and small
$\delta>0$ in

$\displaystyle\operatorname{tr}(\operatorname{Cov}\\{\bm{x}|\bm{\bar{r}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)}\\})$
$\displaystyle\approx\frac{1}{C}\sum_{c=1}^{C}\frac{\bm{p}_{c}^{\top}[\bm{f}_{\bm{\theta}}(\bm{\bar{r}}^{\scriptscriptstyle(2)}\\!+\\!\delta\bm{p}_{c},\overline{v}^{\scriptscriptstyle(2)})-\bm{f}_{\bm{\theta}}(\bm{\bar{r}}^{\scriptscriptstyle(2)},\overline{v}^{\scriptscriptstyle(2)})]}{\delta}$
(29)

at the expense of $C\geq 1$ additional calls of the denoiser
$\bm{f}_{\bm{\theta}}$.

In summary, for any forward operator $\bm{A}$ that obeys
$\bm{A}^{\textsf{H}}\bm{A}=\bm{I}$, the proposed deepEC scheme computes
approximate-MMSE solutions to GLM problems by iteratively calling a deep
denoiser, similar to the PnP [27] or RED [28, 29] schemes. But it requires far
fewer iterations, as we demonstrate in Sec. IV.

For phase retrieval, $\bm{A}^{\textsf{H}}\bm{A}=\bm{I}$ holds for the two
classes of $\bm{A}$ that dominate the literature: i) oversampled-Fourier (OSF)
and ii) coded diffraction pattern (CDP) [68]. For OSF,

$\displaystyle\bm{A}=\bm{F}_{m}\bm{O}\text{~{}~{}for~{}~{}}\bm{O}\triangleq\begin{bmatrix}\bm{I}_{d}\\\
\bm{0}_{m-d\times d}\end{bmatrix}\in{\mathbb{R}}^{m\times d},$ (30)

where $\bm{F}_{m}\in{\mathbb{C}}^{m\times m}$ is a unitary 2D Fourier
transform and $\bm{O}$ pads the vectorized image with zeros. For CDP,

$\displaystyle\bm{A}=\frac{1}{\sqrt{K}}\begin{bmatrix}\bm{F}_{d}\operatorname{Diag}(\bm{c}_{1})\\\\[-2.84526pt]
\vdots\\\ \bm{F}_{d}\operatorname{Diag}(\bm{c}_{K})\end{bmatrix},$ (31)

where $\bm{F}_{d}\in{\mathbb{C}}^{d\times d}$ is a unitary 2D Fourier
transform and $\bm{c}_{k}\in{\mathbb{C}}^{d}$ are random code vectors with
entries drawn independently and uniformly from the unit circle in
${\mathbb{C}}$. By inspection, $\bm{A}^{\textsf{H}}\bm{A}=\bm{I}$ under (30)
with $m\geq d$ and (31) with $K\geq 1$.

###  III-D deepEC for Phase Retrieval

For PR, we propose to apply the deepEC algorithm from Sec. III-C with the
stochastic damping scheme from Sec. III-A. In addition, we use the minor
modifications described below. The overall scheme, “deepECpr,” is detailed in
Alg. 2.

_Initialization:_ Unlike the non-informative
$\bm{\bar{z}}^{\scriptscriptstyle(2)}$ initialization used in Alg. 1, deepECpr
uses the HIO algorithm [12] to obtain an initial estimate of $\bm{x}$, just
like the state-of-the-art prDeep [19] and Deep-ITA [31] methods. Then, to
initialize $\bm{\bar{z}}^{\scriptscriptstyle(2)}$, deepECpr transforms the HIO
recovery $\hat{\bm{x}}_{\mathsf{init}}$ to the $\bm{z}$-domain and adds
circular-complex AWGN:

$\bm{\bar{z}}^{\scriptscriptstyle(2)}\leftarrow\bm{A}\hat{\bm{x}}_{\mathsf{init}}+\bm{n}\text{~{}with~{}}\bm{n}\sim\mathcal{CN}(\bm{0},\overline{v}_{\mathsf{init}}\bm{I}).$
(32)

By using a large noise variance (nominally
$\overline{v}_{\mathsf{init}}\approx 100^{2}$ for pixel amplitudes in
$[0,255]$) the structured error artifacts in
$\bm{A}\hat{\bm{x}}_{\mathsf{init}}$ are suppressed. Because the error
variance of $\bm{\bar{z}}^{\scriptscriptstyle(2)}$ is larger than
$\overline{v}_{\mathsf{init}}$ due to the error in
$\bm{A}\hat{\bm{x}}_{\mathsf{init}}$, the initial
$\overline{v}^{\scriptscriptstyle(2)}$ is set to
$\zeta\overline{v}_{\mathsf{init}}$ for some $\zeta>1$ (nominally
$\zeta=1.2$).

_Real-valued $\bm{x}$ and complex-valued $\bm{A}$:_ The sufficient statistic
$\bm{\bar{r}}^{\scriptscriptstyle(2)}=\bm{A}^{\textsf{H}}\bm{\bar{z}}^{\scriptscriptstyle(2)}$
in (22) would be appropriate if both $\bm{x}$ and $\bm{A}$ were real-valued,
or if both were complex-valued. However, for the PR problems that we consider,
$\bm{x}$ is real-valued and $\bm{A}$ is complex-valued, and so we form the
sufficient statistic as

$\displaystyle\bm{\bar{r}}^{\scriptscriptstyle(2)}\triangleq\Re\\{\bm{A}^{\textsf{H}}\bm{\bar{z}}^{\scriptscriptstyle(2)}\\}=\bm{x}+\bm{\epsilon},\quad\bm{\epsilon}\sim\mathcal{N}(\bm{0},0.5\overline{v}^{\scriptscriptstyle(2)}\bm{I}),$
(33)

and call the deepECpr denoiser as
$\hat{\bm{z}}^{\scriptscriptstyle(2)}=\bm{A}\bm{f}_{\bm{\theta}}(\bm{\bar{r}}^{\scriptscriptstyle(2)},0.5\overline{v}^{\scriptscriptstyle(2)})$.

_Simplified variance-prediction:_ As described in Sec. III-B, the deep
variance predictor $h_{\bm{\phi}}$ aims to estimate the error variance
$\widehat{v}^{\scriptscriptstyle(2)}$ of the transformed denoiser output
$\hat{\bm{z}}^{\scriptscriptstyle(2)}$. Although the Monte-Carlo scheme (29)
is an option, it involves additional denoiser calls, which would increase
runtime. We thus propose to use the simple approximation

$\displaystyle\widehat{v}^{\scriptscriptstyle(2)}\approx\beta\overline{v}^{\scriptscriptstyle(2)},$
(34)

where $\beta\in(0,1)$ is a denoiser-specific constant. In words, (34) models
the denoiser’s output error variance $\widehat{v}^{\scriptscriptstyle(2)}$ as
a fixed fraction of its input error variance
$\overline{v}^{\scriptscriptstyle(2)}$.

_Laplace approximation:_ For PR, we adopt the likelihood

$\displaystyle
p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i})=\mathcal{N}(y_{i};|z_{i}|,v),~{}~{}i=1,\dots,m,$
(35)

due to its excellent performance in practice [53]. However, the resulting EC
posterior

$\displaystyle
q_{2}(z_{i};\overline{z}^{\scriptscriptstyle(1)}_{i},\overline{v}^{\scriptscriptstyle(1)})$
$\displaystyle=\frac{p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i})\mathcal{CN}(z_{i};\overline{z}^{\scriptscriptstyle(1)}_{i},\overline{v}^{\scriptscriptstyle(1)})}{\int
p_{\mathsf{y}|\mathsf{z}}(y_{i}|z_{i}^{\prime})\mathcal{CN}(z_{i}^{\prime};\overline{z}^{\scriptscriptstyle(1)}_{i},\overline{v}^{\scriptscriptstyle(1)})\mathop{}\\!\mathrm{d}z_{i}^{\prime}}$
(36) $\displaystyle\triangleq
p_{\mathsf{z}|\mathsf{y}}(z_{i}|y_{i};\overline{z}^{\scriptscriptstyle(1)}_{i},\overline{v}^{\scriptscriptstyle(1)})$
(37)

does not have analytically tractable mean or variance. Consequently, we employ
the Laplace approximation [69], which assigns
$\widehat{z}_{i}^{\scriptscriptstyle(1)}$ to the posterior mode (i.e., the MAP
estimate) and $\widehat{v}_{i}^{\scriptscriptstyle(1)}$ to the trace of the
inverse Hessian of $-\ln
p_{\mathsf{z}|\mathsf{y}}(z_{i}|y_{i};\overline{z}^{\scriptscriptstyle(1)}_{i},\overline{v}^{\scriptscriptstyle(1)})$
at $z_{i}=\widehat{z}_{i}^{\scriptscriptstyle(1)}$, both of which have closed-
form expressions. In particular, the MAP estimate of $z_{i}$ is given by [31]

$\displaystyle\widehat{z}^{\scriptscriptstyle(1)}_{i}=\frac{\overline{v}^{\scriptscriptstyle(1)}y_{i}+2v|\overline{z}^{\scriptscriptstyle(1)}_{i}|}{\overline{v}^{\scriptscriptstyle(1)}+2v}e^{j\angle\overline{z}^{\scriptscriptstyle(1)}_{i}},$
(38)

and the trace-inverse-Hessian at
$z_{i}=\widehat{z}_{i}^{\scriptscriptstyle(1)}$ is

$\displaystyle\widehat{v}_{i}^{\scriptscriptstyle(1)}=\frac{\overline{v}^{\scriptscriptstyle(1)}(\overline{v}^{\scriptscriptstyle(1)}y_{i}+4v|\overline{z}^{\scriptscriptstyle(1)}_{i}|)}{2|\overline{z}^{\scriptscriptstyle(1)}_{i}|(\overline{v}^{\scriptscriptstyle(1)}+2v)}$
(39)

as derived in App. B. The overall posterior variance
$\widehat{v}^{\scriptscriptstyle(1)}$ is then set to the average value over
all $i$, i.e.,
$\widehat{v}^{\scriptscriptstyle(1)}=\tfrac{1}{m}\sum_{i=1}^{m}\widehat{v}_{i}^{\scriptscriptstyle(1)}$.

Figure 2: The thirty $256\times 256$ FFHQ test images.

Figure 3: The natural (left) and unnatural (right) grayscale test images.

##  IV Numerical Experiments

In this section, we present numerical experiments with Poisson-shot-noise-
corrupted OSF and CDP phaseless measurements, comparing the performance of the
proposed deepECpr algorithm to the classical HIO [12] and state-of-the-art
prDeep [19], Deep-ITA (both “F” and “S” variants) [31], and DPS [41]
approaches.

###  IV-A Test Data and Denoisers

We test on three image datasets. As in the DPS paper [41], we use $256\times
256$ color images from the FFHQ test partition [70]. In particular, we use the
first thirty images, shown in Fig. 2. As in the Deep-ITA paper [31], we use
the six “natural” $128\times 128$ grayscale images shown in Fig. 3, and as in
the prDeep paper [19], we use the six “unnatural” $128\times 128$ grayscale
images shown in Fig. 3.

When recovering the FFHQ test images, we use a denoiser that was constructed
by rescaling (see [41, eq. (9)]) an unconditional DDPM diffusion model [64]
trained on 49 000 FFHQ training images [71]. The resulting diffusion denoiser
“$\bm{f}_{\mathsf{diff}}(\bm{r}_{\mathsf{noisy}},v_{\mathsf{in}})$” accepts a
noisy image $\bm{r}_{\mathsf{noisy}}$ and an input-noise variance
$v_{\mathsf{in}}$, and thus is directly compatible with deepECpr and DPS. The
prDeep and Deep-ITA methods, however, do not have a way to estimate
$v_{\mathsf{in}}$ and instead sequence through a bank of four “blind”
denoisers designed to work at input-noise standard deviations (SD) of 60, 40,
20, 10 (relative to pixel values in $[0,255]$), respectively. Thus, for prDeep
and Deep-ITA, we use
$\bm{f}_{\mathsf{diff}}(\bm{r}_{\mathsf{noisy}},v_{\mathsf{in}})$ with
$v_{\mathsf{in}}\in\\{60^{2},40^{2},20^{2},10^{2}\\}$ for FFHQ image recovery.

When recovering the grayscale images, we follow prDeep and Deep-ITA in using a
bank of four blind DnCNN [32] denoisers
$\\{\bm{f}_{\mathsf{dncnn\text{-}60}}(\cdot),\bm{f}_{\mathsf{dncnn\text{-}40}}(\cdot),\bm{f}_{\mathsf{dncnn\text{-}20}}(\cdot),\bm{f}_{\mathsf{dncnn\text{-}10}}(\cdot)\\}$.
Although deepECpr would likely work better with a non-blind denoiser, we
wanted a fair comparison with prDeep and Deep-ITA. These DnCNN denoisers were
trained over the SD intervals $\\{[40,60],[20,40],[10,20],[0,10]\\}$,
respectively, on the BSD400 [72] dataset using the bias-free approach from
[73]. Note that our natural grayscale test images are similar to (but distinct
from) the BSD400 training data, while our unnatural grayscale test images
could be considered out-of-distribution.

###  IV-B Measurement Generation

Like the prDeep and Deep-ITA papers, we employ OSF and CDP forward operators
$\bm{A}$ (recall (30)-(31)) with 4$\times$ oversampling and the Poisson shot-
noise mechanism [19]

$\displaystyle
y_{i}^{2}=|z_{i}|^{2}+w_{i}\text{~{}with~{}}w_{i}\sim\mathcal{N}(0,\alpha^{2}|z_{i}|^{2}),$
(40)

where
$y_{i}^{2}/\alpha^{2}\sim\mathcal{N}(|z_{i}|^{2}/\alpha^{2},|z_{i}|^{2}/\alpha^{2})$
is approximately $\mathrm{Poisson}((|z_{i}|/\alpha)^{2})$ for sufficiently
small noise levels $\alpha$.

###  IV-C Performance Evaluation

We quantify performance using PSNR and SSIM [74] after resolving certain
ambiguities fundamental to PR.

Phaseless Fourier measurements are unaffected by spatial translations,
conjugate flips, and global phase [75]. The phase ambiguity is circumvented
when the image pixels are non-negative, as with all of our test images.
Furthermore, the translation ambiguity is avoided when oversampling is used
with images that have a sufficient number of non-zero edge pixels, as with our
color FFHQ test images, but not our grayscale test images. Thus, when
recovering the grayscale test images from OSF measurements, we correct for
both translations and flips, but when recovering the color FFHQ test images,
we correct for flips only. With color images, however, one further
complication arises with HIO. Because HIO separately recovers each color
channel, the flip ambiguity can manifest differently in each channel. To
address this latter problem, we fix the first channel of the HIO recovery and
flip of each remaining channel to maximize its correlation with the first
channel.

Phaseless CDP measurements experience only a global-phase ambiguity [76],
which is inconsequential with our test images, and so we perform no ambiguity
resolution.

###  IV-D Algorithm Setup

For HIO, we ported the MATLAB implementation from [77] to Python and set the
tunable step-size parameter to $0.9$. For prDeep, we use the Python
implementation from [78]. For Deep-ITA, we implemented the algorithm in Python
since no public information is available. For DPS, we use the Python
implementation from [71]. As in the prDeep and Deep-ITA papers, we use
$p_{\mathsf{y}|\mathsf{z}}$ from (4) for all algorithms (except HIO, which
does not use a likelihood). For deepECpr, we use damping factors
$\mu^{\scriptscriptstyle(1)}=0.3$ and $\mu^{\scriptscriptstyle(2)}=0.075$ and
initialization factor $\zeta=1.2$. EM-based auto-tuning [79] is used to re-
evaluate $\overline{v}^{\scriptscriptstyle(2)}$ for the first 3 and 10
iterations in the FFHQ and grayscale experiments, respectively. The initial
noise level $\overline{v}_{\mathsf{init}}$ is set to $120^{2}$ for the FFHQ
experiments and $70^{2}$ for the grayscale experiments. Experiment-specific
settings are described later.

_Initialization/Reporting for OSF:_ For the OSF experiments, we apply the
initialization/reporting protocol used in both the prDeep and Deep-ITA papers
to those algorithms, as well as to HIO and deepECpr. First, HIO is run 50
times, for 50 iterations each, from a random initialization. The candidate
$\hat{\bm{x}}$ with the lowest measurement residual $\|\bm{y}-|\bm{Ax}|\|$ is
then selected as the HIO initialization, after which HIO is run for 1000
iterations. In the case of color images, the second and third channels are
flipped to best match the first. Finally, the resulting HIO output is used to
initialize prDeep, Deep-ITA, and deepECpr. The entire procedure is repeated
three times, and the reconstruction with the lowest measurement residual is
reported as the final output for each algorithm. For DPS, we follow the
authors’ recommendation to run the algorithm four times and report the
reconstruction with the lowest measurement residual.

_Initialization/Reporting for CDP:_ For the CDP experiments, we use the
initialization
$\hat{\bm{x}}_{\mathsf{init}}=e^{j\angle{\bm{A}^{\textsf{H}}\bm{1}}}\bm{A}^{\textsf{H}}\bm{y}$
for HIO, prDeep, Deep-ITA, and deepECpr, where
$e^{j\angle{\bm{A}^{\textsf{H}}\bm{1}}}$ is computed and applied element-wise.
Each algorithm was run only once.

###  IV-E OSF Phase Retrieval of FFHQ Images

TABLE I: Average PSNR and SSIM for FFHQ phase retrieval with OSF and CDP
operators and noise level $\alpha$.

method | OSF | CDP
---|---|---
$\alpha=4$ | $\alpha=6$ | $\alpha=8$ | $\alpha=5$ | $\alpha=15$ | $\alpha=45$
PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
HIO | 27.37 | 0.6759 | 26.08 | 0.6163 | 25.03 | 0.5664 | 36.48 | 0.9140 | 26.79 | 0.6064 | 17.58 | 0.2330
Deep-ITA-F | 35.05 | 0.9420 | 34.94 | 0.9374 | 34.50 | 0.9321 | 41.69 | 0.9786 | 37.20 | 0.9556 | 29.93 | 0.8402
Deep-ITA-S | 34.01 | 0.9365 | 34.54 | 0.9367 | 35.15 | 0.9412 | 41.87 | 0.9795 | 37.15 | 0.9531 | 27.28 | 0.7817
prDeep | 37.69 | 0.9654 | 35.28 | 0.9523 | 33.68 | 0.9410 | 42.43 | 0.9816 | 37.48 | 0.9584 | 23.14 | 0.4477
DPS | 27.22 | 0.7674 | 25.84 | 0.7530 | 24.57 | 0.7408 | 41.52 | 0.9786 | 35.68 | 0.9381 | 30.12 | 0.8428
deepECpr (proposed) | 39.75 | 0.9720 | 37.01 | 0.9567 | 34.86 | 0.9404 | 43.12 | 0.9846 | 37.55 | 0.9589 | 32.15 | 0.8941

Figure 4: Top: Example FFHQ image recoveries from phaseless OSF measurements
at noise level $\alpha=8$, with PSNR indicated in the top right corner of each
image. Bottom: Zoomed versions of the cyan regions in the top row. Figure 5:
Average PSNR versus iteration for FFHQ phase retrieval from OSF measurements
at noise level $\alpha=6$. Figure 6: Evolution of the true and estimated
$\bm{z}$ errors in deepECpr for OSF phase retrieval at noise level $\alpha=8$.

We first investigate the recovery of FFHQ test images from noisy, phaseless
OSF measurements. For prDeep and Deep-ITA, we sequence through the denoisers
$\bm{f}_{\mathsf{diff}}(\cdot,60^{2}),\bm{f}_{\mathsf{diff}}(\cdot,40^{2}),\bm{f}_{\mathsf{diff}}(\cdot,20^{2})$,
and $\bm{f}_{\mathsf{diff}}(\cdot,10^{2})$ for $300$ iterations each, for a
total of $1200$ iterations, as recommended by the authors. Also, we set the
$\lambda$ parameter to $0.1\sqrt{v}$ for prDeep, and $0.05\sqrt{v}$ for Deep-
ITA, to maximize validation PSNR. For DPS, we use $1000$ steps as recommended
by the authors and set the scale parameter to $0.0075$ to maximize validation
PSNR. For deepECpr, we use $\bm{f}_{\mathsf{diff}}$ as specified in Alg. 2 for
$200$ iterations with denoising factor $\beta=0.15$.

Table I shows the average PSNR and SSIM at noise levels
$\alpha\in\\{4,6,8\\}$. There, deepECpr outperforms the other methods in both
PSNR and SSIM at $\alpha\in\\{4,6\\}$, where the PSNR advantage is significant
($\approx 2$ dB). At $\alpha=8$, deepECpr trails the best-performing
technique, Deep-ITA-S, by a small margin: $0.29$ dB in PSNR and $0.0008$ in
SSIM.

Figure 4 shows example recoveries at $\alpha=8$. There we see severe artifacts
in the HIO and DPS recoveries. Although the prDeep, Deep-ITA-F, Deep-ITA-S,
and deepECpr approaches avoid unwanted artifacts, the proposed deepECpr does
the best job of recovering fine details in the ground-truth image.

Figure 5 plots PSNR (averaged over the FFHQ test images) versus iteration at
$\alpha=6$. Note that the prDeep, Deep-ITA, DPS, and deepECpr algorithms all
call the denoiser once per iteration and thus have similar per-iteration
complexities. From the figure, we see that deepECpr converges an order-of-
magnitude faster than its competitors.

Figure 6 plots the standard-deviations (SDs) of the deepECpr-estimated and
true errors in the $\bm{z}$ estimates for an typical OSF recovery at noise
level $\alpha=8$. In particular, it plots
$\sqrt{\overline{v}^{\scriptscriptstyle(i)}}$ in comparison to the true error
$m^{-1/2}\|\bm{\bar{z}}^{\scriptscriptstyle(i)}-\bm{z}\|$, and
$\sqrt{\widehat{v}^{\scriptscriptstyle(i)}}$ in comparison to the true error
$m^{-1/2}\|\hat{\bm{z}}^{\scriptscriptstyle(i)}-\bm{z}\|$, for
$i\in\\{1,2\\}$. The figure shows the accuracy of
$\overline{v}^{\scriptscriptstyle(2)}$ over all iterations and the accuracy of
the other estimates after iteration $100$. Furthermore, it shows the EC fixed-
point condition
$\widehat{v}^{\scriptscriptstyle(1)}=\widehat{v}^{\scriptscriptstyle(2)}$
being satisfied.

Table II presents an ablation study comparing the proposed stochastic damping
(16) to the deterministic damping (15) when used in line 14 of Alg. 2. With
stochastic damping, deepECpr yields higher PSNR and SSIM values at all tested
noise levels $\alpha$. Figure 8 plots the SD of the AWGN noise injected by
stochastic damping as a function of the deepECpr iteration. The plot shows
that the SD starts high but decreases to zero with the iterations.

###  IV-F CDP Phase Retrieval of FFHQ Images

Figure 7: Top: Example FFHQ image recoveries from phaseless CDP measurements
at noise level $\alpha=45$, with PSNR indicated in the top right corner of
each image. Bottom: Zoomed versions of the cyan regions in the top row. TABLE
II: PSNR and SSIM of deterministic versus stochastic damping for FFHQ OSF
phase retrieval under noise level $\alpha$

type of damping used | $\alpha=4$ | $\alpha=6$ | $\alpha=8$
---|---|---|---
in line 14 of Alg. 2 | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
deterministic damping (15) | 38.05 | 0.9565 | 36.3 | 0.9468 | 34.20 | 0.9267
stochastic damping (16) | 39.75 | 0.9720 | 37.01 | 0.9567 | 34.86 | 0.9404

Figure 8: Standard deviation of the noise added by deepECpr’s stochastic
damping scheme versus iteration for an example of FFHQ phase retrieval with
OSF measurements at noise level $\alpha=8$. Figure 9: Average PSNR versus
iteration for FFHQ phase retrieval from CDP measurements at noise level
$\alpha=5$. TABLE III: Average PSNR and SSIM for natural and unnatural
grayscale image phase retrieval with CDP measurements at noise level $\alpha$.

method | natural | unnatural
---|---|---
$\alpha=5$ | $\alpha=15$ | $\alpha=45$ | $\alpha=5$ | $\alpha=15$ | $\alpha=45$
PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
HIO | 36.39 | 0.9541 | 26.65 | 0.7663 | 17.13 | 0.4155 | 36.45 | 0.8891 | 26.94 | 0.6133 | 18.55 | 0.3175
Deep-ITA-F | 38.67 | 0.9783 | 29.57 | 0.8595 | 22.77 | 0.6301 | 39.42 | 0.9607 | 29.80 | 0.7526 | 21.40 | 0.4352
Deep-ITA-S | 38.80 | 0.9797 | 28.56 | 0.8357 | 17.25 | 0.4260 | 39.80 | 0.9679 | 28.56 | 0.7077 | 18.74 | 0.3211
prDeep | 38.73 | 0.9785 | 32.49 | 0.9388 | 26.38 | 0.8167 | 39.63 | 0.9660 | 33.84 | 0.9267 | 27.57 | 0.8220
deepECpr (proposed) | 39.46 | 0.9845 | 32.75 | 0.9439 | 27.06 | 0.8428 | 40.92 | 0.9808 | 34.17 | 0.9356 | 28.14 | 0.8292

Figure 10: Top: Example unnatural image recoveries from phaseless CDP
measurements at noise level $\alpha=15$, with PSNR indicated in the top right
corner of each image. Bottom: Zoomed versions of the cyan regions in the top
row.

We now investigate the recovery of FFHQ test images from noisy, phaseless CDP
measurements. Deep-ITA, prDeep, DPS, and deepECpr are configured as described
in Sec. IV-E, except that the Deep-ITA uses $\lambda=0.01\sqrt{v}$ and DPS
uses scale parameter $0.05$, in both cases to maximize validation PSNR.
deepECpr uses the denoising factor $\beta=0.002$.

Table I shows the average PSNR and SSIM values at noise levels
$\alpha\in\\{5,15,45\\}$. There we see that deepECpr outperforms the other
methods in both PSNR and SSIM at all $\alpha$ and gives a PSNR advantage of
$>2$ dB at $\alpha=45$.

Figure 7 shows example recoveries at $\alpha=45$. There we see that both HIO
and prDeep fall prey to noise-like artifacts. Meanwhile, Deep-ITA-F and Deep-
ITA-S over-smooth the image and introduce artifacts. DPS more accurately
recovers textures, but at the expense of hallucinations. Finally, deepECpr
does the best job of faithfully recovering the true image.

Figure 9 plots the PSNR (averaged over the FFHQ test images) versus iteration
at $\alpha=5$. As before, deepECpr converges an order-of-magnitude faster than
its competitors.

###  IV-G CDP Phase Retrieval of Grayscale Images

We now investigate the recovery of natural and unnatural grayscale test images
from noisy, phaseless CDF measurements. We omit DPS from the comparison due to
the unavailability of a pretrained unconditional diffusion model that is
compatible with our natural and unnatural grayscale test images, and the
impossibility of training a good one with a small dataset such as BSD400. The
algorithms under test are configured as described in Sec. IV-D. In addition,
prDeep uses the author-recommended $\lambda=0.1\sqrt{v}$, while Deep-ITA uses
$\lambda=0.25\sqrt{v}$ to maximize validation PSNR. For deepECpr, we use
$\beta=0.045,0.035,0.025,0.020$ for
$\bm{f}_{\mathsf{dncnn\text{-}60}},\bm{f}_{\mathsf{dncnn\text{-}40}},\bm{f}_{\mathsf{dncnn\text{-}20}},\bm{f}_{\mathsf{dncnn\text{-}10}}$,
respectively.

Table III shows the average PSNR and SSIM values at noise levels
$\alpha\in\\{5,15,45\\}$. There deepECpr achieves the highest PSNR and SSIM
values at all tested noise levels for both natural and unnatural grayscale
images.

Figure 10 shows example recoveries at $\alpha=15$. There we see that
deepECpr’s recovery looks cleaner and sharper than those of the competing
methods.

TABLE IV: Average PSNR and SSIM for natural and unnatural grayscale test
images with OSF phase retrieval at noise level $\alpha$.

method | natural | unnatural
---|---|---
$\alpha=4$ | $\alpha=6$ | $\alpha=8$ | $\alpha=4$ | $\alpha=6$ | $\alpha=8$
PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
HIO | 24.32 | 0.6932 | 22.75 | 0.6267 | 21.63 | 0.5702 | 26.66 | 0.6221 | 25.10 | 0.5653 | 24.04 | 0.5231
Deep-ITA-F | 31.19 | 0.8996 | 29.64 | 0.8642 | 28.35 | 0.8278 | 28.30 | 0.6904 | 26.76 | 0.6466 | 26.58 | 0.6450
Deep-ITA-S | 30.98 | 0.8933 | 29.71 | 0.8603 | 28.78 | 0.8343 | 28.20 | 0.6894 | 26.96 | 0.6472 | 26.75 | 0.6390
prDeep | 35.71 | 0.9632 | 32.56 | 0.9242 | 29.57 | 0.8670 | 30.22 | 0.7533 | 26.66 | 0.6630 | 25.96 | 0.6297
deepECpr (proposed) | 37.30 | 0.9766 | 34.06 | 0.9546 | 31.85 | 0.9302 | 30.98 | 0.7675 | 27.45 | 0.6897 | 27.09 | 0.6780

Figure 11: Top: Example natural image recoveries from phaseless OSF
measurements at noise level $\alpha=8$, with PSNR indicated in the top right
corner of each image. Bottom: Zoomed versions of the cyan regions in the top
row.

###  IV-H OSF Phase Retrieval of Grayscale Images

Finally, we investigate the recovery of natural and unnatural grayscale test
images from noisy, phaseless OSF measurements. The algorithms are set up as
described in Sec. IV-G, except that prDeep and Deep-ITA use the author-
recommended choices $\lambda=\sqrt{v}$ and $\lambda=0.025\sqrt{v}$,
respectively, and deepECpr uses $\beta=0.45,0.35,0.25,0.20$ for
$\bm{f}_{\mathsf{dncnn\text{-}60}},\bm{f}_{\mathsf{dncnn\text{-}40}},\bm{f}_{\mathsf{dncnn\text{-}20}},\bm{f}_{\mathsf{dncnn\text{-}10}}$
over $150,100,75,75$ iterations, respectively.

Table IV shows the average PSNR and SSIM values at noise levels
$\alpha\in\\{4,6,8\\}$. There, for both natural and unnatural test images,
deepECpr achieves the highest PSNR and SSIM at all tested noise levels. For
natural images, deepECpr’s PSNR exceeds the second-best PSNR by $\geq 1.5$ dB.

Figure 11 shows example recoveries at $\alpha=8$. There we see that deepECpr
recovers the fine details and texture of the original image, while the other
approaches show noise-like artifacts and overly rough textures.

##  V Conclusion

In this paper, we proposed “deepECpr,” which combines expectation-consistent
(EC) approximation with deep denoising networks to advance the state-of-the-
art in phase retrieval. In addition to a non-traditional application of EC,
deepECpr employs a novel stochastic damping scheme. Experimental results with
oversampled-Fourier and coded-diffraction-pattern operators, and with color,
natural, and unnatural grayscale images, demonstrate deepECpr’s advantages in
image recovery from noise-corrupted phaseless measurements. In almost all
cases, deepECpr yielded improved PSNR and SSIM over the state-of-the-art
prDeep, Deep-ITA, and DPS methods with 5x fewer denoiser calls. Although not
explicitly tested, the proposed deepEC approach should be directly applicable
to other generalized-linear-model problems like dequantization, logistic
regression, and image recovery in non-Gaussian noise.

## References

  * [1] A. Walther, “The question of phase retrieval in optics,” Optica Acta, vol. 10, no. 1, pp. 41–49, 1963.
  * [2] J. Zhang, N. Pégard, J. Zhong, H. Adesnik, and L. Waller, “3D computer-generated holography by non-convex optimization,” Optica, vol. 4, no. 10, pp. 1306–1313, 2017.
  * [3] R. P. Millane, “Phase retrieval in crystallography and optics,” J. Opt. Soc. Amer. A, vol. 7, pp. 394–411, Mar. 1990.
  * [4] R. W. Harrison, “Phase problem in crystallography,” J. Opt. Soc. Amer. A, vol. 10, no. 5, pp. 1046–1055, 1993.
  * [5] J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of X-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature, vol. 400, no. 6742, pp. 342–344, 1999.
  * [6] C. A. Metzler, F. Heide, P. Rangarajan, M. M. Balaji, A. Viswanath, A. Veeraraghavan, and R. G. Baraniuk, “Deep-inverse correlography: Towards real-time high-resolution non-line-of-sight imaging,” Optica, vol. 7, no. 1, pp. 63–71, 2020.
  * [7] C. Fienup and J. Dainty, “Phase retrieval and image reconstruction for astronomy,” Image recovery: Theory and Application, vol. 231, p. 275, 1987\.
  * [8] M. Stefik, “Inferring DNA structures from segmentation data,” Artificial Intelligence, vol. 11, no. 1-2, pp. 85–114, 1978.
  * [9] R. Balan, P. G. Casazza, and D. Edidin, “On signal reconstruction without phase,” Appl. Comput. Harmonic Anal., vol. 20, pp. 345–356, May 2006.
  * [10] Y. Shechtman, Y. C. Eldar, O. Cohen, H. N. Chapman, J. Miao, and M. Segev, “Phase retrieval with application to optical imaging: A contemporary overview,” IEEE Signal Process. Mag., vol. 32, no. 3, pp. 87–109, 2015\.
  * [11] R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of the phase from image and diffraction plane pictures,” Optika, vol. 35, no. 2, pp. 237–246, 1972.
  * [12] J. R. Fienup, “Phase retrieval algorithms: A comparison,” Appl. Optics, vol. 21, pp. 2758–2769, Aug. 1982.
  * [13] H. H. Bauschke, P. L. Combettes, and D. R. Luke, “Hybrid projection-reflection method for phase retrieval,” J. Opt. Soc. Amer. A, vol. 20, no. 6, pp. 1025–1034, 2003.
  * [14] V. Elser, “Phase retrieval by iterated projections,” J. Opt. Soc. Amer. A, vol. 20, no. 1, pp. 40–55, 2003.
  * [15] C.-C. Chen, J. Miao, C. Wang, and T. Lee, “Application of optimization technique to noncrystalline X-ray diffraction microscopy: Guided hybrid input-output method,” Physical Rev. B, vol. 76, no. 6, p. 064113, 2007\.
  * [16] E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval via Wirtinger flow: Theory and algorithms,” IEEE Trans. Info. Theory, vol. 61, no. 4, pp. 1985–2007, 2015.
  * [17] M. Mondelli and A. Montanari, “Fundamental limits of weak recovery with applications to phase retrieval,” in Proc. Conf. Learn. Theory, pp. 1445–1450, 2018.
  * [18] W. Luo, W. Alghamdi, and Y. M. Lu, “Optimal spectral initialization for signal recovery with applications to phase retrieval,” IEEE Trans. Signal Process., vol. 67, no. 9, pp. 2347–2356, 2019.
  * [19] C. A. Metzler, P. Schniter, A. Veeraraghavan, and R. G. Baraniuk, “prDeep: Robust phase retrieval with flexible deep neural networks,” in Proc. Intl. Conf. Mach. Learn., pp. 3501–3510, 2018.
  * [20] E. J. Candès, T. Strohmer, and V. Voroninski, “PhaseLift: Exact and stable signal recovery from magnitude measurements via convex programming,” Commun. Pure & Appl. Math., vol. 66, no. 8, pp. 1241–1274, 2013.
  * [21] T. Goldstein and C. Studer, “Phasemax: Convex phase retrieval via basis pursuit,” IEEE Trans. Info. Theory, vol. 64, no. 4, pp. 2675–2689, 2018\.
  * [22] A. Fannjiang and T. Strohmer, “The numerics of phase retrieval,” Acta Numerica, vol. 29, pp. 125–228, 2020.
  * [23] P. Schniter and S. Rangan, “Compressive phase retrieval via generalized approximate message passing,” IEEE Trans. Signal Process., vol. 63, pp. 1043–1055, Feb. 2015.
  * [24] P. Schniter, S. Rangan, and A. K. Fletcher, “Vector approximate message passing for the generalized linear model,” in Proc. Asilomar Conf. Signals Syst. Comput., pp. 1525–1529, 2016.
  * [25] A. Maillard, B. Loureiro, F. Krzakala, and L. Zdeborová, “Phase retrieval in high dimensions: Statistical and computational phase transitions,” in Proc. Neural Info. Process. Syst. Conf., vol. 33, pp. 11071–11082, 2020\.
  * [26] M. Mondelli and R. Venkataramanan, “Approximate message passing with spectral initialization for generalized linear models,” in Proc. Intl. Conf. Artificial Intell. Statist., pp. 397–405, 2021.
  * [27] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg, “Plug-and-play priors for model based reconstruction,” in Proc. IEEE Global Conf. Signal Info. Process., pp. 945–948, 2013.
  * [28] Y. Romano, M. Elad, and P. Milanfar, “The little engine that could: Regularization by denoising (RED),” SIAM J. Imag. Sci., vol. 10, no. 4, pp. 1804–1844, 2017.
  * [29] E. T. Reehorst and P. Schniter, “Regularization by denoising: Clarifications and new interpretations,” IEEE Trans. Comput. Imag., vol. 5, pp. 52–67, Mar. 2019.
  * [30] F. Heide, S. Diamond, M. Nießner, J. Ragan-Kelley, W. Heidrich, and G. Wetzstein, “Proximal: Efficient image optimization using proximal algorithms,” ACM Trans. Graph., vol. 35, no. 4, pp. 1–15, 2016.
  * [31] Y. Wang, X. Sun, and J. Fleischer, “When deep denoising meets iterative phase retrieval,” in Proc. Intl. Conf. Mach. Learn., pp. 10007–10017, 2020.
  * [32] K. Zhang, W. Zuo, Y. Chen, D. Meng, and L. Zhang, “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,” IEEE Trans. Image Process., vol. 26, pp. 3142–3155, July 2017.
  * [33] C. A. Metzler, A. Maleki, and R. G. Baraniuk, “BM3D-PRGAMP: Compressive phase retrieval based on BM3D denoising,” in Proc. IEEE Intl. Conf. Image Process., pp. 2504–2508, 2016.
  * [34] A. Bora, A. Jalal, E. Price, and A. G. Dimakis, “Compressed sensing using generative models,” in Proc. Intl. Conf. Mach. Learn., pp. 537–546, 2017\.
  * [35] P. Hand, O. Leong, and V. Voroninski, “Phase retrieval under a generative prior,” in Proc. Neural Info. Process. Syst. Conf., vol. 31, 2018.
  * [36] F. Wang, Y. Bian, H. Wang, M. Lyu, G. Pedrini, W. Osten, G. Barbastathis, and G. Situ, “Phase imaging with an untrained neural network,” Light: Science & Applications, vol. 9, no. 1, p. 77, 2020.
  * [37] E. Bostan, R. Heckel, M. Chen, M. Kellman, and L. Waller, “Deep phase decoder: Self-calibrating phase microscopy with an untrained deep neural network,” Optica, vol. 7, no. 6, pp. 559–562, 2020.
  * [38] D. Ulyanov, A. Vedaldi, and V. Lempitsky, “Deep image prior,” in Proc. IEEE Conf. Comp. Vision Pattern Recog., 2018.
  * [39] A. Sinha, J. Lee, S. Li, and G. Barbastathis, “Lensless computational imaging through deep learning,” Optica, vol. 4, no. 9, pp. 1117–1125, 2017.
  * [40] Y. Rivenson, Y. Zhang, H. Günaydın, D. Teng, and A. Ozcan, “Phase recovery and holographic image reconstruction using deep learning in neural networks,” Light: Science & Applications, vol. 7, no. 2, p. 17141, 2018\.
  * [41] H. Chung, J. Kim, M. T. McCann, M. L. Klasky, and J. C. Ye, “Diffusion posterior sampling for general noisy inverse problems,” in Proc. Intl. Conf. Learn. Rep., 2023.
  * [42] N. Torem, R. Ronen, Y. Y. Schechner, and M. Elad, “Complex-valued retrievals from noisy images using diffusion models,” in Proc. IEEE Intl. Conf. Comput. Vis., pp. 3810–3820, 2023.
  * [43] M. Opper and O. Winther, “Expectation consistent approximate inference,” J. Mach. Learn. Res., vol. 1, pp. 2177–2204, 2005.
  * [44] A. K. Fletcher, M. Sahraee-Ardakan, S. Rangan, and P. Schniter, “Expectation consistent approximate inference: Generalizations and convergence,” in Proc. IEEE Intl. Symp. Info. Theory, pp. 190–194, 2016.
  * [45] S. Rangan, P. Schniter, and A. K. Fletcher, “Vector approximate message passing,” IEEE Trans. Info. Theory, vol. 65, pp. 6664–6684, Oct. 2019\.
  * [46] S. K. Shastri, R. Ahmad, and P. Schniter, “Deep expectation-consistent approximation for phase retrieval,” in Proc. Asilomar Conf. Signals Syst. Comput., (Pacific Grove, CA), pp. 910–914, Oct. 2023.
  * [47] P. McCullagh and J. A. Nelder, Generalized Linear Models.  London: Chapman & Hall/CRC, 2nd ed., 1989.
  * [48] T. Hastie, R. Tibshirani, J. H. Friedman, and J. H. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction.  Springer, 2nd ed., 2009.
  * [49] R. M. Willett, R. F. Marcia, and J. M. Nichols, “Compressed sensing for practical optical imaging systems: A tutorial,” Optical Eng., vol. 50, July 2011.
  * [50] N. S. Jayant and P. Noll, Digital Coding of Waveforms: Principles and Applications to Speech and Video.  Prentice-Hall, 1984.
  * [51] J. Dong, L. Valzania, A. Maillard, T.-A. Pham, S. Gigan, and M. Unser, “Phase retrieval: From computational imaging to machine learning: A tutorial,” IEEE Signal Process. Mag., vol. 40, no. 1, pp. 45–57, 2023.
  * [52] V. Katkovnik and J. Astola, “Phase retrieval via spatial light modulator phase modulation in 4f optical setup: Numerical inverse imaging with sparse regularization for phase and amplitude,” J. Opt. Soc. Amer. A, vol. 29, no. 1, pp. 105–116, 2012.
  * [53] L.-H. Yeh, J. Dong, J. Zhong, L. Tian, M. Chen, G. Tang, M. Soltanolkotabi, and L. Waller, “Experimental robustness of Fourier ptychography phase retrieval algorithms,” Opt. Express, vol. 23, no. 26, pp. 33214–33240, 2015.
  * [54] S. Boyd and L. Vandenberghe, Convex Optimization.  Providence, RI: Cambridge, 2004.
  * [55] P. L. Combettes and J.-C. Pesquet, “Proximal splitting methods in signal processing,” in Fixed-Point Algorithms for Inverse Problems in Science and Engineering (H. Bauschke, R. Burachik, P. Combettes, V. Elser, D. Luke, and H. Wolkowicz, eds.), pp. 185–212, Springer, 2011.
  * [56] S. Rangan, P. Schniter, and A. K. Fletcher, “Vector approximate message passing,” in Proc. IEEE Intl. Symp. Info. Theory, pp. 1588–1592, 2017\.
  * [57] A. K. Fletcher, P. Pandit, S. Rangan, S. Sarkar, and P. Schniter, “Plug-in estimation in high-dimensional linear inverse problems: A rigorous analysis,” in Proc. Neural Info. Process. Syst. Conf., pp. 7440–7449, 2018\.
  * [58] P. Pandit, M. Sahraee-Ardakan, S. Rangan, P. Schniter, and A. K. Fletcher, “Inference with deep generative priors in high dimensions,” IEEE J. Sel. Areas Info. Theory, vol. 1, no. 1, pp. 336–347, 2020.
  * [59] K. Takeuchi, “Rigorous dynamics of expectation-propagation-based signal recovery from unitarily invariant measurements,” IEEE Trans. Info. Theory, vol. 66, pp. 368–386, Jan. 2020.
  * [60] M. Opper and O. Winther, “Expectation consistent free energies for approximate inference,” in Proc. Neural Info. Process. Syst. Conf., pp. 1001–1008, 2005.
  * [61] T. Minka, A Family of Approximate Algorithms for Bayesian Inference.  PhD thesis, Dept. Comp. Sci. Eng., MIT, Cambridge, MA, USA, 2001.
  * [62] S. Sarkar, S. Rangan, A. K. Fletcher, and P. Schniter, “Bilinear recovery using adaptive vector-AMP,” IEEE Trans. Signal Process., vol. 67, no. 13, pp. 3383–3396, 2019.
  * [63] S. Sarkar, R. Ahmad, and P. Schniter, “MRI image recovery using damped denoising vector AMP,” in Proc. IEEE Intl. Conf. Acoust. Speech Signal Process., pp. 8108–8112, 2021.
  * [64] J. Ho, A. Jain, and P. Abbeel, “Denoising diffusion probabilistic models,” in Proc. Neural Info. Process. Syst. Conf., vol. 33, pp. 6840–6851, 2020.
  * [65] S. Arridge, P. Maass, O. Öktem, and C.-B. Schönlieb, “Solving inverse problems using data-driven models,” Acta Numerica, vol. 28, pp. 1–174, June 2019.
  * [66] H. V. Poor, An Introduction to Signal Detection and Estimation.  New York, NY, USA: Springer, 2nd ed., 1994.
  * [67] S. Ramani, T. Blu, and M. Unser, “Monte-Carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms,” IEEE Trans. Image Process., vol. 17, pp. 1540–1554, Sept. 2008.
  * [68] E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval from coded diffraction patterns,” Appl. Comput. Harmonic Anal., vol. 39, no. 2, pp. 277–299, 2015.
  * [69] C. M. Bishop, Pattern Recognition and Machine Learning.  New York: Springer, 2007.
  * [70] T. Karras, S. Laine, and T. Aila, “A style-based generator architecture for generative adversarial networks,” in Proc. IEEE Conf. Comp. Vision Pattern Recog., pp. 4396–4405, 2019.
  * [71] H. Chung, J. Kim, M. T. McCann, M. L. Klasky, and J. C. Ye, “diffusion-posterior-sampling.” https://github.com/DPS2022/diffusion-posterior-sampling, Mar. 2023.
  * [72] D. Martin, C. Fowlkes, D. Tal, and J. Malik, “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,” in Proc. IEEE Intl. Conf. Comput. Vis., vol. 2, pp. 416–423, July 2001.
  * [73] S. Mohan, Z. Kadkhodaie, E. P. Simoncelli, and C. Fernandez-Granda, “Robust and interpretable blind image denoising via bias-free convolutional neural networks,” in Proc. Intl. Conf. Learn. Rep., 2020.
  * [74] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” IEEE Trans. Image Process., vol. 13, pp. 600–612, Apr. 2004.
  * [75] T. Bendory, R. Beinert, and Y. C. Eldar, “Fourier phase retrieval: Uniqueness and algorithms,” in MATHEON Conf. Compressed Sensing and its Applications, pp. 55–91, 2015.
  * [76] A. Fannjiang, “Absolute uniqueness of phase retrieval with random illumination,” Inverse Problems, vol. 28, no. 7, p. 075008, 2012.
  * [77] C. A. Metzler, “prdeep.” https://github.com/ricedsp/prDeep/tree/master, 2018\.
  * [78] D. Hekstra, I. Hunt-Isaak, J. Greisman, and J. Russell, “phase-retrieval.” https://github.com/Hekstra-Lab/phase-retrieval, 2018.
  * [79] A. K. Fletcher, M. Sahraee-Ardakan, S. Rangan, and P. Schniter, “Rigorous dynamics and consistent estimation in arbitrarily conditioned linear systems,” in Proc. Neural Info. Process. Syst. Conf., pp. 2542–2551, 2017\.

##  Appendix A Variance of
$p(z_{i}|y_{i};\overline{z}^{\scriptscriptstyle(1)}_{i},\overline{v}^{\scriptscriptstyle(1)})$

For brevity, we omit the variable index “$i$” since the development is
identical for all indices $i$, as well as the superscript
$(\cdot)^{\scriptscriptstyle(1)}$. We set $\widehat{v}$ at the approximation
of the posterior variance given by the Laplace approximation [69], which is

$\displaystyle\operatorname{tr}\\{\bm{H}(\widehat{z})^{-1}\\}\text{~{}~{}for~{}~{}}\bm{H}(\widehat{z})\triangleq\nabla_{z}^{2}(-\ln
p(z|y;\overline{z},\overline{v}))\big{|}_{z=\widehat{z}},$ (41)

where $\widehat{z}$ is the MAP estimate given in (38), i.e.,

$\displaystyle\widehat{z}=\frac{\overline{v}y+2v|\overline{z}|}{\overline{v}+2v}e^{j\angle\overline{z}},$
(42)

and $\nabla_{z}^{2}$ denotes the Hessian with respect to the real and
imaginary parts of $z$.

In App. B, $\operatorname{tr}\\{\bm{H}(\widehat{z})^{-1}\\}$ is derived to be

$\displaystyle\operatorname{tr}\big{\\{}\bm{H}(\widehat{z})^{-1}\big{\\}}$
$\displaystyle=v\frac{2\frac{|\widehat{z}|}{y}-\frac{\overline{v}}{\overline{v}+2v}}{\frac{\overline{v}+2v}{\overline{v}}\frac{|\widehat{z}|}{y}-1}$
(43)

Assuming that $y\geq-2v|\overline{z}|/\overline{v}$, we have

$\displaystyle 2\frac{|\widehat{z}|}{y}-\frac{\overline{v}}{\overline{v}+2v}$
$\displaystyle=\frac{2\overline{v}+4v|\overline{z}|/y}{\overline{v}+2v}-\frac{\overline{v}}{\overline{v}+2v}=\frac{\overline{v}+4v|\overline{z}|/y}{\overline{v}+2v}$
(44)

and

$\displaystyle\frac{\overline{v}+2v}{\overline{v}}\frac{|\widehat{z}|}{y}-1$
$\displaystyle=\frac{\overline{v}+2v}{\overline{v}}\cdot\frac{\overline{v}+2v|\overline{z}|/y}{\overline{v}+2v}-1=\frac{2v|\overline{z}|/y}{\overline{v}}$
(45)

so that

$\displaystyle\operatorname{tr}\big{\\{}\bm{H}(\widehat{z})^{-1}\big{\\}}$
$\displaystyle=v\frac{\overline{v}+4v|\overline{z}|/y}{\overline{v}+2v}\cdot\frac{\overline{v}}{2v|\overline{z}|/y}=\frac{\overline{v}(\overline{v}y+4v|\overline{z}|)}{2|\overline{z}|(\overline{v}+2v)}.$
(46)

##  Appendix B Derivation of the Trace Inverse Hessian

Here we derive $\nabla_{z}^{2}(-\ln p(z|y;\overline{z},\overline{v}))$, the
Hessian of $-\ln p(z|y;\overline{z},\overline{v})$ with respect to the real
and imaginary components of $z\in{\mathbb{C}}$. We will use $z_{r}$ and
$z_{j}$ to denote the real and imaginary components of $z$, and
$\overline{z}_{r}$ and $\overline{z}_{j}$ to denote the real and imaginary
components of $\overline{z}$, respectively. We have

$\displaystyle-\ln
p(z|y;\overline{z},\overline{v})=\tfrac{1}{2v}(y-|z|)^{2}+\tfrac{1}{\overline{v}}|{z}-\overline{z}|^{2}+\text{const}$
(47)
$\displaystyle=\tfrac{1}{2v}\Big{(}y-\sqrt{z_{r}^{2}+z_{j}^{2}}\Big{)}^{2}+\tfrac{1}{\overline{v}}\big{[}({z}_{r}-\overline{z}_{r})^{2}+({z}_{j}-\overline{z}_{j})^{2}\big{]}$
$\displaystyle\quad+\text{const}.$ (48)

This implies that

$\displaystyle-\frac{\partial}{\partial z_{r}}\ln
p(z|y;\overline{z},\overline{v})$ (49) $\displaystyle=$
$\displaystyle\frac{1}{v}\Big{(}\sqrt{{z}_{r}^{2}+{z}_{j}^{2}}-y\Big{)}\frac{1}{2}({z}_{r}^{2}+{z}_{j}^{2})^{-1/2}2{z}_{r}+\frac{2}{\overline{v}}({z}_{r}-\overline{z}_{r})\qquad$
$\displaystyle=$
$\displaystyle\frac{1}{v}\bigg{(}1-\frac{y}{\sqrt{\smash{{z}_{r}^{2}+{z}_{j}^{2}}\vphantom{z_{r}^{2}}}}\bigg{)}{z}_{r}+\frac{2}{\overline{v}}({z}_{r}-\overline{z}_{r}),$
(50)

and similarly that

$\displaystyle-\frac{\partial}{\partial z_{j}}\ln
p(z|y;\overline{z},\overline{v})$
$\displaystyle=\frac{1}{v}\bigg{(}1-\frac{y}{\sqrt{\smash{{z}_{r}^{2}+{z}_{j}^{2}}\vphantom{z_{r}^{2}}}}\bigg{)}{z}_{j}+\frac{2}{\overline{v}}({z}_{j}-\overline{z}_{j}).$
(51)

The second derivatives are

$\displaystyle-\frac{\partial^{2}}{\partial z_{r}^{2}}\ln
p(z|y;\overline{z},\overline{v})$ $\displaystyle=\frac{\partial}{\partial
z_{r}}\bigg{(}\frac{1}{v}\bigg{(}1-\frac{y}{\sqrt{\smash{{z}_{r}^{2}+{z}_{j}^{2}}\vphantom{z_{r}^{2}}}}\bigg{)}{z}_{r}+\frac{2}{\overline{v}}({z}_{r}-\overline{z}_{r})\bigg{)}$
(52)
$\displaystyle=\frac{1}{v}\bigg{(}1-\frac{y}{\sqrt{\smash{{z}_{r}^{2}+{z}_{j}^{2}}\vphantom{z_{r}^{2}}}}+\frac{yz_{r}}{2[{{z}_{r}^{2}+{z}_{j}^{2}}]^{3/2}}2z_{r}\bigg{)}+\frac{2}{\overline{v}}$
(53)
$\displaystyle=\frac{1}{v}\bigg{(}1-\frac{y({z}_{r}^{2}+{z}_{j}^{2})}{[{{z}_{r}^{2}+{z}_{j}^{2}}]^{3/2}}+\frac{yz_{r}^{2}}{[{{z}_{r}^{2}+{z}_{j}^{2}}]^{3/2}}\bigg{)}+\frac{2}{\overline{v}}$
(54)
$\displaystyle=\frac{1}{v}\bigg{(}1-\frac{y{z}_{j}^{2}}{|z|^{3}}\bigg{)}+\frac{2}{\overline{v}},$
(55)

and

$\displaystyle-\frac{\partial^{2}}{\partial z_{j}^{2}}\ln
p(z|y;\overline{z},\overline{v})$
$\displaystyle=\frac{1}{v}\bigg{(}1-\frac{y{z}_{r}^{2}}{|z|^{3}}\bigg{)}+\frac{2}{\overline{v}},$
(56)

and

$\displaystyle-\frac{\partial^{2}}{\partial z_{j}\partial z_{r}}\ln
p(z|y;\overline{z},\overline{v})$ $\displaystyle=\frac{\partial}{\partial
z_{j}}\bigg{(}\frac{1}{v}\bigg{(}1-\frac{y}{\sqrt{\smash{{z}_{r}^{2}+{z}_{j}^{2}}\vphantom{z_{r}^{2}}}}\bigg{)}{z}_{r}+\frac{2}{\overline{v}}({z}_{r}-\overline{z}_{r})\bigg{)}$
(57)
$\displaystyle=\frac{1}{v}\frac{yz_{r}}{2[{{z}_{r}^{2}+{z}_{j}^{2}}]^{3/2}}2z_{j}=\frac{yz_{r}z_{j}}{v|z|^{3}},$
(58)

and so the Hessian matrix is

$\displaystyle\bm{H}(z)$
$\displaystyle=\begin{bmatrix}-\frac{\partial^{2}}{\partial z_{r}^{2}}\ln
p(z|y;\overline{z},\overline{v})&-\frac{\partial^{2}}{\partial z_{j}\partial
z_{r}}\ln p(z|y;\overline{z},\overline{v})\\\ -\frac{\partial^{2}}{\partial
z_{j}\partial z_{r}}\ln
p(z|y;\overline{z},\overline{v})&-\frac{\partial^{2}}{\partial z_{j}^{2}}\ln
p(z|y;\overline{z},\overline{v})\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\frac{1}{v}\big{(}1-\frac{y{z}_{j}^{2}}{|z|^{3}}\big{)}+\frac{2}{\overline{v}}&\frac{yz_{r}z_{j}}{v|z|^{3}}\\\
\frac{yz_{r}z_{j}}{v|z|^{3}}&\frac{1}{v}\big{(}1-\frac{y{z}_{r}^{2}}{|z|^{3}}\big{)}+\frac{2}{\overline{v}}\end{bmatrix}$
(59)
$\displaystyle=\bigg{(}\frac{1}{v}+\frac{2}{\overline{v}}\bigg{)}\bm{I}_{2}-\frac{y}{v|z|^{3}}\begin{bmatrix}z_{j}\\\
-z_{r}\end{bmatrix}\begin{bmatrix}z_{j}&-z_{r}\end{bmatrix}.$ (60)

Using the matrix inversion lemma,

$\displaystyle\operatorname{tr}\big{\\{}(a\bm{I}_{2}+b\bm{cc}^{\top})^{-1}\big{\\}}=\operatorname{tr}\Big{\\{}\frac{1}{a}\Big{(}\bm{I}_{2}-\frac{b\bm{cc}^{\top}}{a+b\bm{c}^{\top}\bm{c}}\Big{)}\Big{\\}}$
(61)
$\displaystyle=\frac{1}{a}\Big{(}2-\frac{b\bm{c}^{\top}\bm{c}}{a+b\bm{c}^{\top}\bm{c}}\Big{)}=\frac{2+a^{-1}b\bm{c}^{\top}\bm{c}}{a+b\bm{c}^{\top}\bm{c}}.$
(62)

Applying this to (60) using

$\displaystyle a$
$\displaystyle=\frac{1}{v}+\frac{2}{\overline{v}}=\frac{\overline{v}+2v}{\overline{v}v}$
(63) $\displaystyle b\bm{c}^{\top}\bm{c}$
$\displaystyle=-\frac{y}{v|z|^{3}}\begin{bmatrix}z_{j}&-z_{r}\end{bmatrix}\begin{bmatrix}z_{j}\\\
-z_{r}\end{bmatrix}=-\frac{y}{v|z|},$ (64)

we get

$\displaystyle\operatorname{tr}\big{\\{}\bm{H}(z)^{-1}\big{\\}}$
$\displaystyle=\frac{2-\frac{\overline{v}v}{\overline{v}+2v}\frac{y}{v|z|}}{\frac{\overline{v}+2v}{\overline{v}v}-\frac{y}{v|z|}}=v\frac{2\frac{|z|}{y}-\frac{\overline{v}}{\overline{v}+2v}}{\frac{\overline{v}+2v}{\overline{v}}\frac{|z|}{y}-1}.$
(65)