remarkRemark hypothesisHypothesis claimClaim assumptionAssumption Stochastic
ADMM for PtychographyKevin Bui and Zichao (Wendy) Di ex_supplement

# A Stochastic ADMM Algorithm for Large-Scale Ptychography with Weighted
Difference of Anisotropic and Isotropic Total Variation††thanks: Submitted to
the editors DATE. This material is based upon work supported by the U.S.
Department of Energy, Office of Science, under contract number DE-
AC02-06CH11357.

Kevin Bui Department of Mathematics, University of California at Irvine,
Irvine, CA 92697 USA<EMAIL_ADDRESS>Zichao (Wendy) Di Mathematics and
Computer Science Division, Argonne National Laboratory, Lemont, IL ().
<EMAIL_ADDRESS>
###### Abstract

Ptychography is an imaging technique that has various scientific applications,
ranging from biology to optics. The method scans the object of interest in a
series of overlapping positions, thereby generating a set of multiple Fourier
magnitude measurements that are potentially corrupted by noise. From these
measurements, an image of the object can be reconstructed depending on how the
related inverse problem is formulated and solved. In this paper, we propose a
class of variational models that incorporate the weighted
anisotropic–isotropic total variation (AITV), an effective regularizer for
image recovery. This class of models is applicable to measurements corrupted
by either Gaussian or Poisson noise. In order to have the models applicable
for large number of ptychographic scans, we design an efficient stochastic
alternating direction method of multipliers algorithm and establish its
convergence. Numerical experiments demonstrate that from a large set of highly
corrupted Fourier measurements, the proposed stochastic algorithm with AITV
regularization can reconstruct complex-valued images with satisfactory
quality, especially for the phase components.

###### keywords:

phase retrieval, total variation, ptychography, ADMM, Poisson/Gaussian noise,
nonconvex optimization, stochastic optimization

65F22, 65K10, 68U10, 90C06, 90C15, 90C26

##  1 Introduction

Ptychography is a popular imaging technique that combines both coherent
diffractive imaging and scanning transmission microscopy. It has been used in
various industrial and scientific applications, including biology [34, 45,
58], crystallography [14], and optics [44, 48]. To perform a ptychographic
experiment (see Figure 1), a coherent beam is scanned across the object of
interest, where each scan may have overlapping positions with another. The
scanning procedure provides a set of phaseless measurements that can be used
to reconstruct an image of the object of interest.

We describe the 2D ptychography in the discrete setting. Let
$z\in\mathbb{C}^{n^{2}}$ be the object of interest with $n\times n$ pixels and
$\omega\in\mathbb{C}^{m^{2}}$ be the localized 2D probe with $m\times m$
pixels, where $m<n$. Both the object $z$ and the probe $\omega$ are expressed
as vectors in lexiographical order. We denote the set of $N$ masks by
$\\{S_{j}\\}_{j=1}^{N}$, where each $S_{j}\in\mathbb{R}^{m^{2}\times n^{2}}$
is a binary matrix that represents a $(m\times m)$-size window over the image
$z$. The set of phaseless measurements $\\{d_{j}\\}_{j=1}^{N}$ is obtained by
$d_{j}=|\mathcal{F}(P_{j}z)|^{2}=|\mathcal{F}(\omega\circ S_{j}z)|^{2}$, where
$P_{j}\coloneqq\omega\circ S_{j}$ is the $j$th probe,
$\mathcal{F}\in\mathbb{C}^{m^{2}\times m^{2}}$ is the 2D discrete Fourier
operator, the operation $\circ$ is elementwise multiplication, and the
operation $|\cdot|$ is the elementwise absolute value of a vector. We aim to
solve the following ptychographic phase retrieval problems. When the probe is
unknown, the blind ptychographic phase retrieval problem is expressed by

(1) $\displaystyle\text{BP-PR: \qquad To find
}\omega\in\mathbb{C}^{m^{2}}\text{ and }z\in\mathbb{C}^{n^{2}}\text{ such that
}|\mathcal{F}(\omega\circ S_{j}z)|^{2}=d_{j},\;j=1,\ldots,N.$

When the probe $\omega$ is known, (1) reduces to the non-blind case where we
only find $z\in\mathbb{C}^{n^{2}}$.

Figure 1: Schematic of a ptychography experiment.

Multiple algorithms have been developed to solve the non-blind and blind phase
retrieval problems. One of the most popular methods is the ptychographical
iterative engine (PIE) [42], where later refinements led to ePIE [33] and rPIE
[32]. The PIE methods are based on gradient descent applied to each
measurement sequentially. Other gradient-based methods for phase retrieval
include Wirtinger flow [6] and its variants [13, 56, 57], which use careful
initialization by a spectral method and adaptive step sizes. PIE is also one
class of projection algorithms for phase retrieval. Other projection-based
algorithms are hybrid projection-reflection [1], Douglas-Rachford splitting
[39, 46], and the relaxed averaged alternating reflections [31]. The phase
retrieval problem can be formulated as a semidefinite optimization problem.
For example, PhaseLift [7] solves the phase retrieval problem as a trace
(nuclear) norm minimization problem. A nonconvex variant called PhaseLiftOff
subtracts the trace norm by the Frobenius norm in the objective function [55].
PhaseCut proposes a different semidefinite formulation of the phase retrieval
problem by explicitly separating the amplitude and phase variables and
optimize only the values of the phase variables [47]. The phase retrieval
problem can alternatively be written as a saddle point problem [51], solved by
alternating direction method of multipliers (ADMM) [4]. A globally convergent
ADMM algorithm has recently been developed to solve the BP-PR problem [8].
Another globally convergent algorithm is proximal alternating linearized
minimization (PALM) [2], which has also been adapted to solve the BP-PR
problem [12, 21]. For a detailed survey of numerical algorithms for phase
retrieval, please refer to [17].

For large-scale ptychography, when a huge number of scans are collected, many
of the aforementioned algorithms may be inapplicable or may need to be adapted
because of the demanding memory footprint and computational cost. Various
parallel algorithms have been developed. For example, an asynchronous parallel
version of ePIE has been implemented on GPUs, where each partition of a
measurement set is asynchronously processed to obtain a sub-image and the sub-
images are later fused together to form the entire image [35]. A parallel
version of relaxed averaged alternating reflections has recently been
developed for GPU implmentation [15]. Unfortunately, some of these parallel
algorithms require a GPU, which many computers do not have. However, there are
efficient algorithms for large-scale ptychography without the need for a GPU,
although having one could speed up the processing time. A multigrid
optimization framework has been proposed to accelerate large-scale gradient-
based methods for phase retrieval [52]. An overlapping domain decomposition
method combined with ADMM leads to a highly parallel algorithm with good load
balance [9]. To the best of our knowledge, there does not yet exist a
stochastic optimization algorithm for large-scale ptychography that
iteratively processes a batch of measurements. Such an algorithm would be
useful for practitioners who do not have access to multiple cores to perform
parallel computing.

To improve the image reconstruction quality in phase retrieval, total
variation (TV) [43] has been incorporated for the cases when the measurements
are corrupted with Gaussian noise [11] or with Poisson noise [10]. Both cases
consider the isotropic TV approximation:

(2) $\displaystyle\|\nabla
z\|_{2,1}=\sum_{i=1}^{n^{2}}\sqrt{|(\nabla_{x}z)_{i}|^{2}+|(\nabla_{y}z)_{i}|^{2}},$

where $\nabla_{x}$ and $\nabla_{y}$ are the horizontal and vertical difference
operators, respectively, and $(\nabla_{x}z)_{i}$ and $(\nabla_{y}z)_{i}$ are
the $i$th entries of the vectors $\nabla_{x}z$ and $\nabla_{y}z$,
respectively. However, it has been known that isotropic TV tends to blur
oblique edges. An alternative approximation that preserves sharper edges is
the anisotropic TV [16]:

(3) $\displaystyle\|\nabla
z\|_{1}=\sum_{i=1}^{n^{2}}\left(|(\nabla_{x}z)_{i}|+|(\nabla_{y}z)_{i}|\right).$

Overall, TV is meant to approximate the $\ell_{0}$ “norm” of the image
gradient, i.e., $\|\nabla z\|_{0}$, because TV is based on the $\ell_{1}$
norm, a convex relaxation of $\ell_{0}$. A nonconvex alternative to $\ell_{1}$
is $\ell_{1}-\alpha\ell_{2},0<\alpha\leq 1$, which performs well in recovering
sparse solutions in various compressed sensing problems [27, 28, 29, 54]. The
superior performance of $\ell_{1}-\alpha\ell_{2}$ in sparse recovery has
motivated the development of the weighted difference of anisotropic and
isotropic total variation (AITV) [30], which applies $\ell_{1}-\alpha\ell_{2}$
on each gradient vector of an image. Mathematically, AITV is formulated by

(4) $\displaystyle\|\nabla z\|_{1}-\alpha\|\nabla
z\|_{2,1}=\sum_{i=1}^{n^{2}}\left[|(\nabla_{x}z)_{i}|+|(\nabla_{y}z)_{i}|-\alpha\sqrt{|(\nabla_{x}z)_{i}|^{2}+|(\nabla_{y}z)_{i}|^{2}}\right].$

AITV has demonstrated better performance than TV in image denoising, image
deconvolution, image segmentation, and MRI reconstruction [5, 30, 38],
especially in preserving sharper edges.

In this work, we consider the large-scale ptychography problem where the
measurements are corrupted by either Gaussian or Poisson noise. To improve the
image reconstruction quality, the AITV regularization is incorporated. The
overall problem is formulated as a general variational problem, where we
develop an ADMM algorithm to solve it. The ADMM algorithm has subproblems that
can be approximately solved by stochastic gradient descent (SGD) [40].
Although SGD is a generic and popular algorithm for unconstrained optimization
problems whose objective functions have a finite-sum structure, it may not be
directly applicable to the subproblem being solved in the case of
ptychography. Hence, we show how to appropriately apply SGD in order to
develop our specialized stochastic ADMM algorithm. To further modify the
algorithm, we incorporate adaptive step size based on the PIE algorithms [32,
33, 42]. Instead of using all measurements per iteration, this stochastic ADMM
algorithm can iteratively process a batch of measurements to accurately
perform image reconstruction.

The paper is organized as follows. In Section 2, we review notations and
definitions that will be used throughout the paper. Next in Section 3 we
describe the AITV-regularized variational models to solve the image
ptychography problem. Within this section, we design the stochastic ADMM
algorithms to solve these models. Convergence analysis of the stochatsic ADMM
algorithm follows in Section 4. In Section 5, we illustrate the performance of
our proposed stochastic ADMM algorithms and compare them with other competing
algorithms. Lastly, in Section 6, we conclude the paper with summary and
future works.

##  2 Preliminaries

In this section, we describe basic notations used throughout the paper.

Let $z\in\mathbb{C}^{n^{2}}$. The $i$th entry of $z$ is denoted by $(z)_{i}$.
The vector $\mathbf{1}$ is a vector whose entries are all ones. The vector
$\mathbf{0}$ is also defined similarly. The real transpose and conjugate
transpose of $z$ are denoted by $z^{\top}$ and $z^{*}$, respectively. The same
superscript notations are used for the real transpose and conjugate transpose
of matrices. The sign of a complex value $z^{\prime}\in\mathbb{C}$ is given by

$\displaystyle\operatorname{sgn}(z^{\prime})=\begin{cases}\displaystyle\frac{z^{\prime}}{|z^{\prime}|}&\text{
if }z^{\prime}\neq 0,\\\
c\in\\{c^{\prime}\in\mathbb{C}:|c^{\prime}|=1\\}&\text{ if
}z^{\prime}=0.\end{cases}$

The sign of a vector $z\in\mathbb{C}^{n^{2}}$ is denoted by
$\operatorname{sgn}(z)$ and is defined elementwise by
$\operatorname{sgn}(z)_{i}=\operatorname{sgn}(z_{i}),\;i=1,\ldots,n^{2}.$ The
standard basis vectors of $\mathbb{C}^{n^{2}}$ are denoted by
$\\{\mathbf{e}_{i}\\}_{i=1}^{n^{2}}$, where $\mathbf{e}_{i}$ is a vector whose
$i$th component is 1 while all other components are zeros. The diagonal matrix
of a vector $z\in\mathbb{C}^{n^{2}}$ is denoted by
$D_{z}=\text{Diag}(z)=z\mathbf{1}^{\top}\circ I_{n^{2}\times n^{2}}.$

For $p=(p_{x},p_{y})\in\mathbb{C}^{n^{2}}\times\mathbb{C}^{n^{2}}$, its $i$th
entry is $p_{i}=\begin{bmatrix}(p_{x})_{i}\\\
(p_{y})_{i}\end{bmatrix}\in\mathbb{C}^{2}$. We define the following norms on
$\mathbb{C}^{n^{2}}\times\mathbb{C}^{n^{2}}$:

$\displaystyle\|p\|_{1}$
$\displaystyle=\sum_{i=1}^{n^{2}}|(p_{x})_{i}|+|(p_{y})_{i}|,\quad\|p\|_{2}=\sqrt{\sum_{i=1}^{n^{2}}|(p_{x})_{i}|^{2}+|(p_{y})_{i}|^{2}},$
$\displaystyle\|p\|_{2,1}$
$\displaystyle=\sum_{i=1}^{n^{2}}\sqrt{|(p_{x})_{i}|^{2}+|(p_{y})_{i}|^{2}}=\sum_{i=1}^{n^{2}}\|p_{i}\|_{2}.$

The discrete gradient operator
$\nabla:\mathbb{C}^{n^{2}}\rightarrow\mathbb{C}^{n^{2}}\times\mathbb{C}^{n^{2}}$
when specifically applied to the image $z$ is given by $\nabla
z=(\nabla_{x}z,\nabla_{y}z)$, where $\nabla_{x}$ and $\nabla_{y}$ are the
forward horizontal and vertical difference operators.

We also define the proximal operator of a function
$f:\mathbb{C}^{n^{2}}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ as

$\displaystyle\text{prox}_{f(\cdot)}(z^{\prime})=\operatorname*{arg\,min}_{z}f(z)+\frac{1}{2}\|z-z^{\prime}\|_{2}^{2},\,\forall
z^{\prime}\in\mathbb{C}^{n^{2}}.$

##  3 Proposed Model

Throughout the paper, we assume that among the mask set
$\\{S_{j}\\}_{j=1}^{N}$, there exists $j^{\prime}$ such that
$\|S_{j^{\prime}}\mathbf{e}_{i}\|_{1}=1$ for each $i=1,\ldots,n^{2}$. This
assumption ensures that each pixel of an image $z\in\mathbb{C}^{n^{2}}$ is
sampled at least once.

Suppose that the measurements $\\{d_{j}\\}_{j=1}^{N}$ are corrupted by
independent and identically distributed (iid) noise, i.e.,
$d_{j}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{Noise}(|\mathcal{F}(\omega\circ
S_{j}z)|^{2})\text{ for }j=1,\ldots,N.$ We assume that the noise is either
Gaussian or Poisson, both of which are common in phase retrieval. Given an
unknown probe $\omega$, the blind variational model [8] is

(5)
$\displaystyle\min_{\omega\in\mathbb{C}^{m^{2}},z\in\mathbb{C}^{n^{2}}}\sum_{j=1}^{N}\mathcal{B}(|\mathcal{F}(\omega\circ
S_{j}z)|^{2},d_{j}),$

where

(6)
$\displaystyle\mathcal{B}(g,f)=\displaystyle\begin{cases}\displaystyle\frac{1}{2}\|\sqrt{g}-\sqrt{f}\|_{2}^{2},&\text{
amplitude Gaussian metric (AGM)
\cite[cite]{[\@@bibref{}{wen2012alternating}{}{}]}},\\\ \\\
\displaystyle\frac{1}{2}\langle g-f\circ\log(g),\mathbf{1}\rangle,&\text{
intensity Poisson metric (IPM)
\cite[cite]{[\@@bibref{}{chang2018total}{}{}]}}.\end{cases}$

Note that $\sqrt{\cdot}$ is elementwise square root. When the probe $\omega$
is known, (5) simplifies to the non-blind case as a special case where we only
need to find $z\in\mathbb{C}^{n^{2}}$. Hence, throughout the rest of this
section, we will focus on the blind case. To improve image recovery, we
propose a class of AITV-regularized variants of (5).

###  3.1 AITV model

For image ptychography, we propose the following AITV-regularized model:

(7)
$\displaystyle\min_{\omega\in\mathbb{C}^{m^{2}},z\in\mathbb{C}^{n^{2}}}\sum_{j=1}^{N}\mathcal{B}(|\mathcal{F}(\omega\circ
S_{j}z)|^{2},d_{j})+\lambda\left(\|\nabla z\|_{1}-\alpha\|\nabla
z\|_{2,1}\right),\;\lambda>0,\;\alpha\in[0,1].$

To develop an ADMM algorithm of (7), we introduce auxiliary variables
$u=(u_{1},\ldots,u_{N})\in\mathbb{C}^{m^{2}\times N}$ and
$v=(v_{x},v_{y})\in\mathbb{C}^{n^{2}}\times\mathbb{C}^{n^{2}}$ so that we
obtain an equivalent constrained optimization problem

(8) $\displaystyle\begin{aligned}
\min_{u,v,z}\quad&\sum_{j=1}^{N}\mathcal{B}(|u_{j}|^{2},d_{j})+\lambda\left(\|v\|_{1}-\alpha\|v\|_{2,1}\right)\;\textrm{s.t.}&u_{j}=\mathcal{F}(\omega\circ
S_{j}z),\;j=1,\ldots,N,\text{ and }v=\nabla z.\end{aligned}$

The augmented Lagrangian of (8) is

(9) $\displaystyle\begin{aligned}
\mathcal{L}(u,\omega,v,z,\Lambda,y)=&\sum_{j=1}^{N}\left[\mathcal{B}(|u_{j}|^{2},d_{j})+\mathbb{R}\left(\langle\Lambda_{j},u_{j}-\mathcal{F}(\omega\circ
S_{j}z)\rangle\right)+\frac{\beta_{1}}{2}\|u_{j}-\mathcal{F}(\omega\circ
S_{j}z)\|_{2}^{2}\right]\\\
&+\lambda\left(\|v\|_{1}-\alpha\|v\|_{2,1}\right)+\mathbb{R}\left(\langle
y,v-\nabla z\rangle\right)+\frac{\beta_{2}}{2}\|v-\nabla
z\|_{2}^{2},\end{aligned}$

where $\mathbb{R}(\cdot)$ denotes the real component of a complex number;
$\langle\cdot,\cdot\rangle$ denotes the complex inner product between two
vectors; $\Lambda=(\Lambda_{1},\ldots,\Lambda_{N})\in\mathbb{C}^{m^{2}\times
N}$ and $y=(y_{x},y_{y})\in\mathbb{C}^{n^{2}\times n^{2}}$ are Lagrange
multipliers; and $\beta_{1},\beta_{2}>0$ are penalty parameters. The ADMM
algorithm iterates as follows:

(10a) $\displaystyle u^{t+1}$
$\displaystyle\in\operatorname*{arg\,min}_{u}\mathcal{L}(u,\omega^{t},v^{t},z^{t},\Lambda^{t},y^{t}),$
(10b) $\displaystyle\omega^{t+1}$
$\displaystyle\in\operatorname*{arg\,min}_{\omega}\mathcal{L}(u^{t+1},\omega,v^{t},z^{t},\Lambda^{t},y^{t}),$
(10c) $\displaystyle v^{t+1}$
$\displaystyle\in\operatorname*{arg\,min}_{v}\mathcal{L}(u^{t+1},\omega^{t+1},v,z^{t},\Lambda^{t},y^{t}),$
(10d) $\displaystyle z^{t+1}$
$\displaystyle\in\operatorname*{arg\,min}_{z}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z,\Lambda^{t},y^{t}),$
(10e) $\displaystyle\Lambda_{j}^{t+1}$
$\displaystyle=\Lambda_{j}^{t}+\beta_{1}\left(u_{j}^{t+1}-\mathcal{F}(\omega^{t+1}\circ
S_{j}z^{t+1})\right),\quad j=1,\ldots,N,$ (10f) $\displaystyle y^{t+1}$
$\displaystyle=y^{t}+\beta_{2}\left(v^{t+1}-\nabla z^{t+1}\right).$

Next we explain how to solve each subproblem.

####  3.1.1 $u$-subproblem

In (10a), we solve $u_{j}$ independently of each other. For each
$j=1,\ldots,N$, we have

(11) $\displaystyle\begin{aligned}
&u_{j}^{t+1}\in\operatorname*{arg\,min}_{u_{j}}\mathcal{B}(|u_{j}|^{2},d_{j})+\mathbb{R}\left(\langle\Lambda_{j}^{t},u_{j}-\mathcal{F}(P_{j}^{t}z^{t})\rangle\right)+\frac{\beta_{1}}{2}\|u_{j}-\mathcal{F}(P_{j}^{t}z^{t})\|_{2}^{2}\\\
&=\operatorname*{arg\,min}_{u_{j}}\frac{1}{\beta_{1}}\mathcal{B}(|u_{j}|^{2},d_{j})+\frac{1}{2}\left\|u_{j}-\mathcal{F}(P_{j}^{t}z^{t})+\frac{1}{\beta_{1}}\Lambda_{j}^{t}\right\|_{2}^{2}=\text{prox}_{\frac{1}{\beta_{1}}\mathcal{B}(|\cdot|^{2},d_{j})}\left(\mathcal{F}(P_{j}^{t}z^{t})-\frac{1}{\beta_{1}}\Lambda_{j}^{t}\right),\end{aligned}$

where $P_{j}^{t}=\omega^{t}\circ S_{j}$. The proximal operator for each
fidelity term in (6) has a closed-form solution provided in [10, 11], so we
have

(12) $\displaystyle\begin{aligned}
&u_{j}^{t+1}=\begin{cases}\frac{\sqrt{d_{j}}+\beta_{1}\left|\mathcal{F}(P_{j}^{t}z^{t})-\frac{1}{\beta_{1}}\Lambda_{j}^{t}\right|}{1+\beta_{1}}\circ\operatorname{sgn}\left(\mathcal{F}(P_{j}^{t}z^{t})-\frac{1}{\beta_{1}}\Lambda_{j}^{t}\right),&\text{AGM,}\\\
\frac{\beta_{1}|\mathcal{F}(P_{j}^{t}z^{t})-\frac{1}{\beta_{1}}\Lambda_{j}^{t}|+\sqrt{\left(\beta_{1}|\mathcal{F}(P_{j}^{t}z^{t})-\frac{1}{\beta_{1}}\Lambda_{j}^{t}|\right)^{2}+4(1+\beta_{1})d_{j}}}{2(1+\beta_{1})}\circ\operatorname{sgn}\left(\mathcal{F}(P_{j}^{t}z^{t})-\frac{1}{\beta_{1}}\Lambda_{j}^{t}\right),&\text{IPM}.\end{cases}\end{aligned}$

####  3.1.2 $\omega$-subproblem

The $\omega$-subproblem (10b) can be rewritten as

(13) $\displaystyle\begin{aligned}
\omega^{t+1}&\in\operatorname*{arg\,min}_{\omega}\sum_{j=1}^{N}\left[\frac{\beta_{1}}{2}\left\|\mathcal{F}^{-1}\left(u_{j}^{t+1}+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right)-\omega\circ
S_{j}z^{t}\right\|_{2}^{2}\right],\end{aligned}$

which shows that updating $\omega$ requires access to all $N$ probes. Hence,
we develop an alternative update scheme that uses only $b\leq N$ probes.
Instead of solving (10b) exactly, we linearize it as done in [26, 37] to
obtain

(14)
$\displaystyle\omega^{t+1}\in\operatorname*{arg\,min}_{\omega}\langle\nabla_{\omega}\mathcal{L}(u^{t+1},\omega^{t},z^{t},\Lambda^{t},y^{t}),\omega-\omega^{t}\rangle+\frac{1}{2\delta_{\omega}^{t}}\|\omega-\omega^{t}\|_{2}^{2}$

for some constant $\delta_{\omega}^{t}>0$ at iteration $t$. Then (14) is
equivalent to performing gradient descent with step size
$\delta_{\omega}^{t}$:

(15)
$\displaystyle\omega^{t+1}=\omega^{t}-\delta_{\omega}^{t}\nabla_{\omega}\mathcal{L}(u^{t+1},\omega^{t},v^{t},z^{t},\Lambda^{t},y^{t}).$

Next we approximate $\nabla\mathcal{L}_{\omega}$ by its stochastic estimator
$\tilde{\nabla}_{\omega}\mathcal{L}$. thereby updating $\omega^{t+1}$ by SGD
with step size $\delta_{\omega}^{t}>0$:

(16)
$\displaystyle\omega^{t+1}=\omega^{t}-\delta_{\omega}^{t}\tilde{\nabla}_{\omega}\mathcal{L}(u^{t+1},\omega^{t},v^{t},z^{t},\Lambda^{t},y^{t}).$

We derive some candidates for $\tilde{\nabla}_{\omega}\mathcal{L}$. Let
$\mathcal{G}^{t}_{j}(\omega)=\frac{\beta_{1}}{2}\left\|\mathcal{F}^{-1}\left(u_{j}^{t+1}+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right)-\omega\circ
S_{j}z^{t}\right\|_{2}^{2}$, which means that
$\nabla\mathcal{G}_{j}^{t}(\omega)=-\beta_{1}(S_{j}z^{t})^{*}\circ\left[\mathcal{F}^{-1}\left(u_{j}^{t+1}+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right)-\omega\circ
S_{j}z^{t}\right]$. (15) can be rewritten as a gradient descent step with
$N\delta_{\omega}^{t}$:

(17)
$\displaystyle\omega^{t+1}=\omega^{t}-N\delta_{\omega}^{t}\left(\frac{1}{N}\sum_{j=1}^{N}\nabla
G_{j}^{t}(\omega^{t})\right).$

The SGD estimator [3, 40] of $\frac{1}{N}\sum_{j=1}^{N}\nabla
G_{j}^{t}(\omega^{t})$ is $\frac{1}{b}\sum_{j\in n_{t}}\nabla
G_{j}^{t}(\omega^{t})$, where $n^{t}\subset\\{1,\ldots,N\\}$ is a sub-batch of
$N$ masks such that $|n_{t}|=b$. Then the SGD variant (after scaling
$\delta_{\omega}^{t}$) of (15) is

(18)
$\displaystyle\omega^{t+1}=\omega^{t}-\delta_{\omega}^{t}\left(\frac{1}{b}\sum_{j\in
n^{t}}\nabla G_{j}^{t}(\omega^{t})\right).$

This implies that one candidate stochastic estimator for
$\nabla_{\omega}\mathcal{L}$ is

(19) $\displaystyle\begin{aligned}
\tilde{\nabla}_{\omega}^{SGD}\mathcal{L}(u^{t+1},\omega^{t},v^{t},z^{t},\Lambda^{t},y^{t})&=\frac{1}{b}\sum_{j\in
n^{t}}\nabla G_{j}^{t}(\omega^{t})\\\ &=-\frac{\beta_{1}}{b}\sum_{j\in
n^{t}}(S_{j}z^{t})^{*}\circ\left[\mathcal{F}^{-1}\left(u_{j}^{t+1}+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right)-\omega^{t}\circ
S_{j}z^{t}\right].\end{aligned}$

We can further modify (19) by incorporating spatially varying step sizes
inspired by the PIE algorithms [32, 33, 42]. Let

(20)
$\displaystyle\Phi_{j}^{t}=\begin{cases}\displaystyle\frac{\mathbf{1}}{\|S_{j}z^{t}\|_{\infty}^{2}},&\text{ePIE
\cite[cite]{[\@@bibref{}{maiden2009improved}{}{}]}},\\\
\displaystyle\frac{\|S_{j}z^{t}\|_{1}\mathbf{1}}{\|S_{j}z^{t}\|_{\infty}\left(|S_{j}z^{t}|^{2}+\gamma_{\omega}\|S_{j}z^{t}\|_{\infty}^{2}\mathbf{1}\right)},&\text{PIE
\cite[cite]{[\@@bibref{}{rodenburg2004phase}{}{}]}},\\\
\displaystyle\frac{\mathbf{1}}{(1-\gamma_{\omega})|S_{j}z^{t}|^{2}+\gamma_{\omega}\|S_{j}z^{t}\|_{\infty}^{2}\mathbf{1}},&\text{rPIE
\cite[cite]{[\@@bibref{}{maiden2017further}{}{}]}},\end{cases}$

where $\gamma_{\omega}\in[0,1]$ for PIE and rPIE and division is elementwise.
Incorporating $\Phi_{j}^{t}$ into (19), we have another class of stochastic
estimators

(21) $\displaystyle\begin{aligned}
\tilde{\nabla}_{\omega}^{PIE}\mathcal{L}(u^{t+1},\omega^{t},v^{t},z^{t},\Lambda^{t},y^{t})=-\frac{\beta_{1}}{b}\sum_{j\in
n^{t}}\Phi_{j}^{t}\circ(S_{j}z^{t})^{*}\circ\left[\mathcal{F}^{-1}\left(u_{j}^{t+1}+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right)-\omega^{t}\circ
S_{j}z^{t}\right].\end{aligned}$

####  3.1.3 $v$-subproblem

Expanding (10c) gives

(22) $\displaystyle\begin{aligned}
v^{t+1}&\in\operatorname*{arg\,min}_{v}\frac{\lambda}{\beta_{2}}\left(\|v\|_{1}-\alpha\|v\|_{2,1}\right)+\frac{1}{2}\left\|v-\nabla
z^{t}+\frac{y^{t}}{\beta_{2}}\right\|_{2}^{2}\\\
&=\operatorname*{arg\,min}_{v}\sum_{i=1}^{n^{2}}\frac{\lambda}{\beta_{2}}\left(\|v_{i}\|_{1}-\alpha\|v_{i}\|_{2}\right)+\frac{1}{2}\left\|v_{i}-(\nabla
z^{t})_{i}+\frac{(y^{t})_{i}}{\beta_{2}}\right\|_{2}^{2},\end{aligned}$

which means that the solution $v^{t+1}$ can be solved elementwise. As a
result, the subproblem simplifies to

(23)
$\displaystyle(v^{t+1})_{i}=\text{prox}_{\frac{\lambda}{\beta_{2}}(\|\cdot\|_{1}-\alpha\|\cdot\|_{2})}\left((\nabla
z^{t})_{i}-\frac{(y^{t})_{i}}{\beta_{2}}\right).$

A closed-form solution for the proximal operator of $\ell_{1}-\alpha\ell_{2}$
is provided in [27] but only for real-valued vectors. We generalize it to the
complex case in Lemma 3.1, whose proof is delayed to Appendix A.

######  Lemma 3.1.

Given $x^{\prime}\in\mathbb{C}^{n}$, $\lambda>0$. and $\alpha\geq 0$, we have
the following cases:

  1. 1.

When $\|x^{\prime}\|_{\infty}>\lambda$, we have

$\displaystyle x^{*}=(\|\xi\|_{2}+\alpha\lambda)\frac{\xi}{\|\xi\|_{2}},\text{
where }\xi=\operatorname{sgn}(x^{\prime})\circ\max(|x^{\prime}|-\lambda,0).$

  2. 2.

When $(1-\alpha)\lambda<\|x^{\prime}\|_{\infty}\leq\lambda$, we have $x^{*}$
as a 1-sparse vector such that one chooses an index
$i\in\operatorname*{arg\,max}_{j}(|(x^{\prime})_{j}|)$ and have

$\displaystyle(x^{*})_{j}=\begin{cases}(|(x^{\prime})_{j}|+(\alpha-1)\lambda)\operatorname{sgn}((x^{\prime})_{j})&\text{
if }$j = i$,\\\ 0&\text{ if }j\neq i.\end{cases}$

  3. 3.

When $\|x^{\prime}\|_{\infty}\leq(1-\alpha)\lambda$, we have $x^{*}=0$.

Then $x^{*}$ is an optimal solution to

(24)
$\displaystyle\mathrm{prox}_{\lambda(\|\cdot\|_{1}-\alpha\|\cdot\|_{2})}(x^{\prime})=\operatorname*{arg\,min}_{x}\|x\|_{1}-\alpha\|x\|_{2}+\frac{1}{2\lambda}\|x-x^{\prime}\|_{2}^{2}.$

####  3.1.4 $z$-subproblem

(10d) can be rewritten as

(25) $\displaystyle\begin{aligned}
z^{t+1}&\in\operatorname*{arg\,min}_{z}\sum_{j=1}^{N}\left[\frac{\beta_{1}}{2}\left\|u_{j}^{t+1}-\mathcal{F}(P_{j}^{t+1}z)+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right\|_{2}^{2}\right]+\frac{\beta_{2}}{2}\left\|v^{t+1}-\nabla
z+\frac{y^{t}}{\beta_{2}}\right\|_{2}^{2},\end{aligned}$

which implies that $z^{t+1}$ must satisfy the first-order optimality condition

(26) $\displaystyle\begin{aligned}
\left(\beta_{1}\sum_{j=1}^{N}(P_{j}^{t+1})^{*}P_{j}^{t+1}-\beta_{2}\Delta\right)z^{t+1}=\sum_{j=1}^{N}\beta_{1}(P_{j}^{t+1})^{*}\mathcal{F}^{-1}\left(u_{j}^{t+1}+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right)+\beta_{2}\nabla^{\top}\left(v^{t+1}+\frac{y^{t}}{\beta_{2}}\right),\end{aligned}$

where the Laplacian $\Delta=-\nabla^{\top}\nabla$. Since the coefficient
matrix of $z^{t+1}$ is invertible, solving (26) can be performed exactly, but
it could be computationally expensive if the matrix system is extremely large
because of the image size of $z$. Since the coefficient matrix tends to be
sparse, conjugate gradient [22] can be used to solve (26) like in [10, 11],
but it needs access to all $N$ probes and requires at most $n^{2}$ iterations
to attain an exact solution, assuming exact arithmetic. Moreover, it is
sensitive to roundoff error [19].

Alternatively, like in Section 3.1.2 we linearize (25) to obtain the gradient
descent step with step size $\delta_{z}^{t}>0$:

(27) $\displaystyle
z^{t+1}=z^{t}-\delta_{z}^{t}\nabla_{z}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t}).$

We approximate $\nabla_{z}\mathcal{L}$ by its stochastic estimator
$\tilde{\nabla}_{z}\mathcal{L}$ that only has access to $b\leq N$ probes.
Replacing $\nabla_{z}\mathcal{L}$ with $\tilde{\nabla}_{z}\mathcal{L}$ in (27)
gives

(28) $\displaystyle
z^{t+1}=z^{t}-\delta_{z}^{t}\tilde{\nabla}_{z}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t}).$

To design candidates for $\tilde{\nabla}_{z}\mathcal{L}$, we will use the
following lemma:

######  Lemma 3.2.

Let $S\in\mathbb{R}^{m^{2}\times n^{2}}$. If $\mathbf{e}_{i}\in\text{ker}(S)$
for some index $i$, then for any $x\in\mathbb{C}^{m^{2}}$, we have
$(S^{\top}x)_{i}=0$.

######  Proof 3.3.

We have $(S^{\top}x)_{i}=\langle S^{\top}x,\mathbf{e}_{i}\rangle=\langle
x,S\mathbf{e}_{i}\rangle=\langle x,0\rangle=0.$

For brevity, we denote the vectors

(29) $\displaystyle\begin{aligned}
A_{j}^{t}&=-\beta_{1}\left[(P_{j}^{t+1})^{*}\mathcal{F}^{-1}\left(u_{j}^{t+1}+\frac{\Lambda_{j}^{t}}{\beta_{1}}\right)-(P_{j}^{t+1})^{*}P_{j}^{t+1}z^{t}\right],\\\
B^{t}&=-\beta_{2}\left[\nabla^{\top}\left(v^{t+1}+\frac{y^{t}}{\beta_{2}}\right)+\Delta
z^{t}\right].\end{aligned}$

At each element $i=1,\ldots,n^{2}$, (27) becomes

(30) $\displaystyle\begin{aligned}
(z^{t+1})_{i}&=(z^{t})_{i}-\delta_{z}^{t}(\nabla_{z}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t}))_{i}=(z^{t})_{i}-\delta_{z}^{t}\left[\sum_{j=1}^{N}(A_{j}^{t})_{i}+(B^{t})_{i}\right].\end{aligned}$

By Lemma 3.2, since $(P_{j}^{t+1})^{*}=(\omega^{t+1}\circ
S_{j})^{*}=S_{j}^{\top}D_{(\omega^{t+1})^{*}}$, we have $(A_{j}^{t})_{i}=0$ if
$\mathbf{e}_{i}\in\text{ker}(S_{j})$, which means that element $i$ is not
scanned by the mask matrix $S_{j}$. For each $i=1,\ldots,n^{2}$, we define
$N_{i}=\\{j:\mathbf{e}_{i}\not\in\text{ker}(S_{j})\\}$ to be the set of
indices corresponding to the mask matrices that scan element $i$. As a result,
(30) reduces to and can be rewritten as

(31) $\displaystyle\begin{aligned}
(z^{t+1})_{i}&=(z^{t})_{i}-\delta_{z}^{t}\left[\sum_{j\in
N_{i}}(A_{j}^{t})_{i}+(B^{t})_{i}\right]=(z^{t})_{i}-|N_{i}|\delta_{z}^{t}\left[\frac{1}{|N_{i}|}\sum_{j\in
N_{i}}\left((A_{j}^{t})_{i}+\frac{1}{|N_{i}|}(B^{t})_{i}\right)\right].\end{aligned}$

Comparing (30) and (31), we observe that

$\displaystyle\frac{1}{|N_{i}|}\sum_{j\in
N_{i}}\left((A_{j}^{t})_{i}+\frac{1}{|N_{i}|}(B^{t})_{i}\right)\propto(\nabla_{z}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t}))_{i}.$

Thus, a candidate for the stochastic estimator $\tilde{\nabla}_{z}\mathcal{L}$
is the SGD estimator $\tilde{\nabla}_{z}^{SGD}\mathcal{L}$ given by

(32)
$\displaystyle\left(\tilde{\nabla}_{z}^{SGD}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t})\right)_{i}=\frac{1}{|n_{i}^{t}|}\sum_{j\in
n_{i}^{t}}\left((A_{j}^{t})_{i}+\frac{1}{|N_{i}|}(B^{t})_{i}\right),$

where $n_{i}^{t}\subset N_{i}$ is a mini-batch sampled from $N_{i}$ at
iteration $t$ [3, 40].

We can further modify (32) by incorporating spatially varying step sizes
inspired by the PIE algorithms [32, 33, 42]. We define

(33)
$\displaystyle\Psi_{i,j}=\begin{cases}\displaystyle\frac{1}{\|\omega^{t+1}\|_{\infty}^{2}}&\text{ePIE
\cite[cite]{[\@@bibref{}{maiden2009improved}{}{}]}},\\\
\displaystyle\frac{\|P_{j}^{t+1}\mathbf{e}_{i}\|_{1}}{\|\omega^{t+1}\|_{\infty}\left(\|P_{j}^{t+1}\mathbf{e}_{i}\|_{1}^{2}+\gamma_{z}\|\omega^{t+1}\|_{\infty}^{2}\right)}&\text{PIE
\cite[cite]{[\@@bibref{}{rodenburg2004phase}{}{}]}},\\\
\displaystyle\frac{1}{(1-\gamma_{z})\|P_{j}^{t+1}\mathbf{e}_{i}\|_{1}^{2}+\gamma_{z}\|\omega^{t+1}\|_{\infty}^{2}}&\text{rPIE
\cite[cite]{[\@@bibref{}{maiden2017further}{}{}]}},\end{cases}$

with $\gamma_{z}\in[0,1]$ for PIE and rPIE. A class of PIE candidates for the
stochastic estimator is

(34)
$\displaystyle\left(\tilde{\nabla}_{z}^{PIE}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t})\right)_{i}=\frac{1}{|n_{i}^{t}|}\sum_{j\in
n_{i}^{t}}\Psi_{i,j}\left((A_{j}^{t})_{i}+\frac{1}{|N_{i}|}(B^{t})_{i}\right).$

The overall stochastic ADMM algorithm that solves (7) is provided by Algorithm
1. Notice that the non-blind problem is just a special case with the probe
$\omega$ fixed.

Algorithm 1 Stochastic ADMM to solve (7)

Input:  set of masks $\\{S_{j}\\}_{j=1}^{N}$; model parameters $\lambda>0$,
$\alpha\in[0,1]$; penalty parameters $\beta_{1},\beta_{2}>0$; sequence of step
sizes $\\{(\delta_{\omega}^{t},\delta_{z}^{t})\\}_{t=1}^{\infty}$; batch size
$b\leq N$; PIE factors $\gamma_{z},\gamma_{\omega}\in[0,1]$.

1: Initialize $\omega^{0}$, $z^{0}$,
$\\{u_{j}^{0}\\}_{j=1}^{N}=\\{\Lambda_{j}^{0}\\}_{j=1}^{N}$, $y^{0}=\nabla
z^{0}$.

2: for $t=0$ to $T-1$ do

3: Uniformly sample without replacement the subset
$n^{t}\subset\\{1,\ldots,N\\}$ of batch size $b$.

4: Compute $n_{i}^{t}$ from $n^{t}$, i.e., $n_{i}^{t}=\displaystyle\sum_{j\in
n^{t}}\|S_{j}\mathbf{e}_{i}\|_{1}$.

5: Update $u_{j}^{t+1}$ according to (LABEL:eq:u_closed_form) for each $j\in
n^{t}$.

6: if $\omega$ is unknown then

7: Update
$\omega^{t+1}=\omega^{t}-\delta_{\omega}^{t}\tilde{\nabla}_{\omega}\mathcal{L}(u^{t+1},\omega^{t},v^{t},z^{t},\Lambda^{t},y^{t})$.
See (19) and (21) for a candidate $\tilde{\nabla}_{\omega}\mathcal{L}$.

8: else

9: $\omega^{t+1}=\omega^{t}$.

10: end if

11: Compute
$\displaystyle(v^{t+1})_{i}=\text{prox}_{\frac{\lambda}{\beta_{2}}(\|\cdot\|_{1}-\alpha\|\cdot\|_{2})}\left((\nabla
z^{t})_{i}-\frac{(y^{t})_{i}}{\beta_{2}}\right)$ for all $i=1,\ldots,n^{2}$;
see Lemma 3.1.

12: Update
$(z^{t+1})_{i}=(z^{t})_{i}-\delta_{z}^{t}\left(\tilde{\nabla}_{z}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t})\right)_{i}$
for all $i$ such that $n_{i}^{t}\neq 0$. See (32) and (34) for a candidate
$\tilde{\nabla}_{z}\mathcal{L}$.

13: Compute  $\displaystyle\Lambda_{j}^{t+1}$
$\displaystyle=\Lambda_{j}^{t}+\beta_{1}\left(u_{j}^{t+1}-\mathcal{F}(\omega^{t+1}\circ
S_{j}z^{t+1})\right),\;\forall j\in n^{t},$ $\displaystyle y^{t+1}$
$\displaystyle=y^{t}+\beta_{2}\left(v^{t+1}-\nabla z^{t+1}\right).$

14: end for

Output: $\omega^{*}=\omega^{T}$, $z^{*}=z^{T}$

##  4 Convergence Analysis

We discuss the convergence of Algorithm 1. Although global convergence for
ADMM can be established using Kurdyka-Łojasiewicz assumptions [49], the result
does not apply for our models because our models contain the gradient
operator, which does not satisfy the necessary surjectivity assumption. Hence,
we will prove up to subsequential convergence. The convergence analysis is
based on the analyses done in [10, 11, 51], where under certain assumptions,
they showed that the iterate subsequences of the ADMM algorithms converge to
Karush-Kuhn-Tucker (KKT) points. To simplify notation, let $Z=(u,\omega,v,z)$
and $\Omega=(\Lambda,y)$. We also write $\mathcal{L}(\omega)$, for example, to
represent the Lagrangian with respect to $\omega$ with all other variables
fixed at their most recent values. A KKT point $(Z^{\star},\Omega^{\star})$ of
the Lagrangian (9) satisfies the KKT conditions given by

(35a) $\displaystyle 0$
$\displaystyle\in\begin{cases}\partial|u_{j}^{\star}|(|u_{j}^{\star}|-\sqrt{d_{j}})+\Lambda_{j}^{\star},&\text{if
AGM},\\\
\partial|u_{j}^{\star}|\left(|u_{j}^{\star}|-\displaystyle\frac{\sqrt{d_{j}}}{|u_{j}^{\star}|}\right)+\Lambda_{j}^{\star},&\text{if
IPM,}\end{cases}\quad\text{for }j=1,\ldots,N,$ (35b)
$\displaystyle-\frac{y^{\star}}{\lambda}$
$\displaystyle\in\partial(\|v^{\star}\|_{1}-\alpha\|v^{\star}\|_{2,1}),$ (35c)
$\displaystyle u_{j}^{\star}$ $\displaystyle=\mathcal{F}(\omega^{\star}\circ
S_{j}z^{\star})\quad\text{for }j=1,\ldots,N,$ (35d) $\displaystyle v^{\star}$
$\displaystyle=\nabla z^{\star},$ (35e)
$\displaystyle\nabla_{\omega}\mathcal{L}(Z^{\star},\Omega^{\star})$
$\displaystyle=0,$ (35f)
$\displaystyle\nabla_{z}\mathcal{L}(Z^{\star},\Omega^{\star})$
$\displaystyle=0.$

Because we implement SGD to solve for the probe $\omega$ and the image $z$ in
the ADMM algorithm, we replace (35e) and (35f) with the following conditions,
respectively:

(35e′)
$\displaystyle\mathbb{E}\left[\left\|\nabla_{\omega}\mathcal{L}(Z^{\star},\Omega^{\star})\right\|_{2}^{2}\right]$
$\displaystyle=0$ (35f′)
$\displaystyle\mathbb{E}\left[\left\|\nabla_{z}\mathcal{L}(Z^{\star},\Omega^{\star})\right\|_{2}^{2}\right]$
$\displaystyle=0.$

We say a point $(Z^{\star},\Omega^{\star})$ is a stochastic KKT point if it
satisfies (35a)-(35d) and (35e′)-(35f′).

Where $\mathbb{E}_{t}$ denotes the expectation conditioned on the first $t$
iterations of the stochastic ADMM algorithm, we impose the following
assumption adapted from [3, Assumption 4.3] relating to the stochastic
gradient estimators $\tilde{\nabla}_{\omega}\mathcal{L}$ and
$\tilde{\nabla}_{z}\mathcal{L}$. Let $\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$
be a sequence of iterates generated by Algorithm 1. Suppose that at each
iteration $t$, the stochastic gradient estimators
$\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\coloneqq\tilde{\nabla}_{\omega}\mathcal{L}(u^{t+1},\omega^{t},v^{t},z^{t},\Lambda^{t},y^{t})$
and
$\tilde{\nabla}_{z}\mathcal{L}(z^{t})\coloneqq\tilde{\nabla}_{\omega}\mathcal{L}(u^{t+1},\omega^{t+1},v^{t+1},z^{t},\Lambda^{t},y^{t})$
satisfy the following:

  1. (a)

There exist constants $K_{U}\geq K_{L}>0$ such that

(37)
$\displaystyle\mathbb{R}\left(\mathbb{E}_{t}\left[\langle\nabla_{\omega}\mathcal{L}(\omega^{t}),\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\rangle\right]\right)$
$\displaystyle\geq
K_{L}\mathbb{E}_{t}\left[\|\nabla_{\omega}\mathcal{L}(\omega^{t})\|_{2}^{2}\right]$
(38)
$\displaystyle\left\|\mathbb{E}_{t}[\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})]\right\|_{2}^{2}$
$\displaystyle\leq
K_{U}\mathbb{E}_{t}\left[\|\nabla_{\omega}\mathcal{L}(\omega^{t})\|_{2}^{2}\right]$
(39)
$\displaystyle\mathbb{R}\left(\mathbb{E}_{t}\left[\langle\nabla_{z}\mathcal{L}(z^{t})),\tilde{\nabla}_{z}\mathcal{L}(z^{t})\rangle\right]\right)$
$\displaystyle\geq
K_{L}\mathbb{E}_{t}\left[\|\nabla_{z}\mathcal{L}(z^{t})\|_{2}^{2}\right]$ (40)
$\displaystyle\left\|\mathbb{E}_{t}[\tilde{\nabla}_{z}\mathcal{L}(z^{t})]\right\|_{2}^{2}$
$\displaystyle\leq
K_{U}\mathbb{E}_{t}\left[\|\nabla_{z}\mathcal{L}(z^{t})\|_{2}^{2}\right].$

  2. (b)

There exists constant $M,M_{V}\geq 0$ such that

(41)
$\displaystyle\mathbb{E}_{t}\left[\|\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\|_{2}^{2}\right]-\left\|\mathbb{E}_{t}[\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})]\right\|_{2}^{2}$
$\displaystyle\leq
M+M_{V}\mathbb{E}_{t}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]$
(42)
$\displaystyle\mathbb{E}_{t}\left[\|\tilde{\nabla}_{z}\mathcal{L}(z^{t})\|_{2}^{2}\right]-\left\|\mathbb{E}_{t}[\tilde{\nabla}_{z}\mathcal{L}(z^{t})]\right\|_{2}^{2}$
$\displaystyle\leq
M+M_{V}\mathbb{E}_{t}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right].$

To prove the convergence of Algorithm 1, we require the following preliminary
results. Lemma 4.1 provides useful inequalities while Proposition 4.2 bounds
the iterates $\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ and establishes some
bounded property of the gradients
$\\{(\nabla_{\omega}\mathcal{L}(\omega^{t}),\nabla_{z}\mathcal{L}(z^{t}))\\}_{t=1}^{\infty}$.

######  Lemma 4.1.

Let $\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ be a sequence of iterates
generated by Algorithm 1 that satisfies Assumption 4. Suppose that
$\\{(\omega^{t},z^{t})\\}_{t=1}^{\infty}$ is bounded. For each iteration $t$,
we have

(43)
$\displaystyle\mathbb{E}_{t}[\mathcal{L}(\omega^{t+1})]-\mathbb{E}_{t}[\mathcal{L}(\omega^{t})]$
$\displaystyle\leq-\left(K_{L}-\frac{\delta^{t}_{\omega}L_{\omega}(M_{V}+K_{U})}{2}\right)\delta^{t}_{\omega}\mathbb{E}_{t}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]+\frac{(\delta_{\omega}^{t})^{2}L_{\omega}M}{2}$
(44)
$\displaystyle\mathbb{E}_{t}[\mathcal{L}(z^{t+1})]-\mathbb{E}_{t}[\mathcal{L}(z^{t})]$
$\displaystyle\leq-\left(K_{L}-\frac{\delta^{t}_{z}L_{z}(M_{V}+K_{U})}{2}\right)\delta^{t}_{z}\mathbb{E}_{t}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]+\frac{(\delta_{z}^{t})^{2}L_{z}M}{2}$

for some constants $L_{\omega},L_{z}>0$.

######  Proposition 4.2.

Let $\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ be a sequence of iterates
generated by Algorithm 1 that satisfies Assumption 4. Suppose
$\\{(\omega^{t},z^{t})\\}_{t=1}^{\infty}$ is bounded,
$\sum_{t=1}^{\infty}\|\Omega^{t+1}-\Omega^{t}\|_{2}^{2}<\infty$, and

(45) $\displaystyle\begin{aligned}
\sum_{t=1}^{\infty}\delta_{\omega}^{t}=\infty,\quad\sum_{t=1}^{\infty}(\delta_{\omega}^{t})^{2}<\infty,\\\
\sum_{t=1}^{\infty}\delta_{z}^{t}=\infty,\quad\sum_{t=1}^{\infty}(\delta_{z}^{t})^{2}<\infty.\end{aligned}$

Then $\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ is bounded and

(46)
$\displaystyle\sum_{t=1}^{\infty}\mathbb{E}\left[\delta_{\omega}^{t}\|\nabla_{\omega}\mathcal{L}(\omega^{t})\|_{2}^{2}+\delta_{z}^{t}\|\nabla_{z}\mathcal{L}(z^{t})\|_{2}^{2}\right]$
$\displaystyle<\infty.$

The convergence of Algorithm 1 is finally established below (see proofs in
Appendix B).

######  Theorem 4.3.

Let $\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ be generated by Algorithm 1.
Under the same assumption as Proposition 4.2, there exists a subsequence of
$\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ whose accumulation point
$(Z^{\star},\Omega^{\star})$ is a stochastic KKT point almost surely (a.s.) of
(9).

We note that the requirement
$\sum_{t=1}^{\infty}\|\Omega^{t+1}-\Omega^{t}\|_{2}^{2}<\infty$ is rather
strong, but similar assumption was made in other nonconvex ADMM algorithms
[24, 25, 53] that do not satisfy the necessary assumptions for global
convergence [49].

##  5 Numerical Results

In this section, we evaluate the performance of Algorithm 1 on two complex
images presented in Figure 2. The probe size used for both images is
$256\times 256$, and the probe is scanned across an image from left to right
and top to bottom, giving us $N=100$ measurements. The measurements
$\\{d_{j}\\}_{j=1}^{N}$ are either corrupted by Gaussian noise or Poisson
noise. More specifically, when the measurements are corrupted by Gaussian
noise, we have

(47) $\displaystyle d_{j}=(|\mathcal{F}(P_{j}z)|+\epsilon)^{2},$

where $\epsilon$ is an i.i.d. Gaussian random vector. When the measurements
are corrupted by Poisson noise, we have

(48) $\displaystyle d_{j}=\text{Poisson}(|\mathcal{F}(P_{j}z_{\zeta})|),$

where $z_{\zeta}=\zeta z$ for some constant $\zeta>0$. Note that Poisson noise
is stronger when $\zeta$ is smaller.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Two complex sample images and their probes examined in the experiments. Left column: sample magnitude; middle column: sample phase with inserted proportionally sized probe magnitude; left column: the magnitude differences between the ground-truth probe and the initial probe $\omega^{0}$. Table 1: Parameter settings for each method. Note that $b$ refers to the batch size. |  Total Epochs | $\beta_{1}=\beta_{2}$ | $\delta_{z}^{t}$ | $\Psi_{i,j}$ | $\delta_{\omega}^{t}$ | $\Phi_{i,j}$
---|---|---|---|---|---|---
AGM | 600 | 0.25 | $\begin{cases}2\sqrt{b}&1\leq t\leq 300\\\ \frac{1}{5}\sqrt{b}&300<t\leq 450\\\ \frac{1}{50}\sqrt{b}&450<t\leq 600\end{cases}$ |  rPIE ($\gamma_{z}=0.1$) | $\begin{cases}\sqrt{b}\times 10^{-3}&1\leq t\leq 300\\\ \sqrt{b}\times 10^{-4}&300<t\leq 450\\\ \sqrt{b}\times 10^{-5}&450\leq t\leq 600\end{cases}$ |  rPIE ($\gamma_{\omega}=0.025$)
IPM | 300 | 0.25 | $\begin{cases}15\sqrt{b}&1\leq t\leq 150\\\ \frac{3}{2}\sqrt{b}&150<t\leq 225\\\ \frac{3}{20}\sqrt{b}&225<t\leq 300\end{cases}$ |  ePIE | $\begin{cases}2\sqrt{b}\times 10^{-3}&1\leq t\leq 300\\\ 2\sqrt{b}\times 10^{-4}&300<t\leq 450\\\ 2\sqrt{b}\times 10^{-5}&450\leq t\leq 600\end{cases}$ |  ePIE

For numerical evaluation, we compute the SSIMs [50] of the magnitudes and
phases between the reconstructed image $z^{**}$ and the ground-truth image
$z^{g}$, where $z^{**}_{i}=\zeta^{*}z^{*}_{i+t^{*}}$ is adjusted for scaling
by $\zeta^{*}$ and translation by $t^{*}$ and
$(\zeta^{*},t^{*})=\displaystyle\operatorname*{arg\,min}_{\zeta\in\mathbb{C},t\in\mathbb{Z}}\sum_{i=1}^{n^{2}}|\zeta
z^{*}_{i+t}-z^{g}_{i}|^{2}.$ We compare the proposed stochastic ADMM
algorithms with its deterministic, full-batch counterparts (i.e., (26) and
(13) are solved exactly) and its isotropic TV counterparts based on [10, 11].
The results are also compared with Douglas-Rachford splitting [46], rPIE [33],
and PHeBIE [21]. We follow the implementation of Douglas-Rachford from [8].

Table 2: SSIM results of the algorithms applied to the Gaussian corrupted measurements with SNR = 40. The stochastic algorithms (e.g., AITV and isoTV, $b\in\\{5,10,20,50\\}$) are ran three times to obtain the average SSIM values. Bold indicates best value; underline indicates second best value. | Non-blind | Blind
---|---|---
| Chip | Cameraman/Baboon | Chip | Cameraman/Baboon
|  mag. SSIM |  phase SSIM |  mag. SSIM |  phase SSIM |  mag. SSIM |  phase SSIM |  mag. SSIM |  phase SSIM
DR | 0.8130 | 0.8089 | 0.8701 | 0.5191 | 0.8008 | 0.7642 | 0.8009 | 0.3207
rPIE | 0.8886 | 0.9073 | 0.8930 | 0.6055 | 0.9070 | 0.9120 | 0.8890 | 0.6145
PHeBIE | 0.8004 | 0.8019 | 0.8725 | 0.5718 | 0.8612 | 0.8438 | 0.8846 | 0.5756
isoTV $(b=5)$ | 0.9501 | 0.9027 | 0.9393 | 0.7578 | 0.9426 | 0.8919 | 0.9324 | 0.7547
isoTV $(b=10)$ | 0.9498 | 0.9004 | 0.9387 | 0.7475 | 0.9429 | 0.8891 | 0.9326 | 0.7477
isoTV $(b=20)$ | 0.9514 | 0.8981 | 0.9385 | 0.7302 | 0.9447 | 0.8850 | 0.9298 | 0.7289
isoTV $(b=50)$ | 0.9355 | 0.9193 | 0.9294 | 0.7050 | 0.9322 | 0.9047 | 0.9153 | 0.7025
isoTV (full batch) | 0.9578 | 0.9145 | 0.9769 | 0.7338 | 0.9527 | 0.8698 | 0.9589 | 0.5774
AITV $(b=5)$ | 0.9585 | 0.9556 | 0.9438 | 0.7720 | 0.9490 | 0.9477 | 0.9373 | 0.7775
AITV $(b=10)$ | 0.9620 | 0.9579 | 0.9515 | 0.7747 | 0.9534 | 0.9481 | 0.9450 | 0.7772
AITV $(b=20)$ | 0.9629 | 0.9583 | 0.9538 | 0.7707 | 0.9547 | 0.9470 | 0.9468 | 0.7690
AITV $(b=50)$ | 0.9585 | 0.9550 | 0.9490 | 0.7358 | 0.9514 | 0.9432 | 0.9391 | 0.7342
AITV (full batch) | 0.9674 | 0.9513 | 0.9814 | 0.7463 | 0.9676 | 0.9296 | 0.9725 | 0.5956

|  |  |  |  |
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |

(a) isoTV ($b=50$) (b) AITV ($b=20$)

(c) AITV (full) (d) DR

(e) rPIE (f) PHeBIE

(g) isoTV $(b=50)$ (h) AITV ($b=20$)

(i) AITV (full) (j) DR

(k) rPIE (l) PHeBIE

(m) isoTV ($b=5$) (n) AITV ($b=10$)

(o) AITV (full) (p) DR

(q) rPIE (r) PHeBIE

(s) isoTV ($b=5$) (t) AITV ($b=10$)

(u) AITV (full) (v) DR

(w) rPIE (x) PHeBIE

Figure 3: Reconstructions of the non-blind case for the Gaussian noise. Top
two rows: reconstructions of Figures 2(a)-2(b); bottom two rows:
reconstructions of Figs. 2(d)-2(e).

We initialize $z^{0}=\frac{1}{\sqrt{2}}(\mathbf{1}+i\mathbf{1})$ when using
AGM for Gaussian-corrupted measurements and
$z^{0}=\frac{\zeta}{\sqrt{2}}(\mathbf{1}+i\mathbf{1})$ when using IPM for
Poisson-corrupted measurements. When performing the blind experiments using
Algorithm 1, $\omega^{0}$ is initialized as the perturbation of the ground-
truth probe. The magnitude differences between the initial probe and the
ground-truth probe are shown in Figure 2. The selected parameters, except for
$\lambda$, are summarized in Table 1. The initial step sizes for
$\delta_{z}^{t}$ and $\delta_{\omega}^{t}$ are determined empirically, and
motivated by (45), we decrease them by a factor of 10 at the $1/2$ and $3/4$
of the total epochs. Decreasing the step size in this way is a popular
technique, especially in the deep learning community [18, 20]. Inspired from
[18], the step sizes are multiplied by a factor of $\sqrt{b}$ so that they are
scaled accordingly to the batch size $b$. For AITV regularization, we examine
$\alpha\in\\{0.2,0.4,0.6,0.8\\}$ and determine that $\alpha=0.8$ yields the
best results across all of our numerical examples. The batch sizes we examine
are $b\in\\{5,10,20,50\\}$ for Gaussian noise and $b\in\\{5,10,20,25\\}$ for
Poisson noise. For each parameter setting and image, we run three trials to
obtain the mean SSIM values.

The code for the experiments is available at
https://github.com/kbui1993/Stochastic_ADMM_Ptycho.

###  5.1 Gaussian noise

The SNR of the noisy measurements [10] is given by

$\displaystyle\text{SNR}\left(\\{\sqrt{d_{j}}\\}_{j=1}^{N},\\{|\mathcal{F}(P_{j}z)|\\}_{j=1}^{N}\right)=-10\log_{10}\left(\frac{\displaystyle\sum_{j=1}^{N}\|\sqrt{d_{j}}-|\mathcal{F}(P_{j}z)|\|_{2}^{2}}{\displaystyle\sum_{j=1}^{N}\|\mathcal{F}(P_{j}z)\|_{2}^{2}}\right),$

so determined by the SNR value, the noise level $\epsilon$ in (47) can be
calculated by

$\displaystyle\epsilon=\sqrt{\frac{10^{-\text{SNR}/10}\displaystyle\sum_{j=1}^{N}\|\mathcal{F}(P_{j}z)\|_{2}^{2}}{Nm^{2}}}.$

|  |  |  |  |
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |
|  |  |  |  |

(a) isoTV ($b=5$) (b) AITV ($b=20$)

(c) AITV (full) (d) DR

(e) rPIE (f) PHeBIE

(g) isoTV $(b=5)$ (h) AITV $(b=20)$

(i) AITV (full) (j) DR

(k) rPIE (l) PHeBIE

(m) isoTV ($b=10$) (n) AITV ($b=10$)

(o) AITV (full) (p) DR

(q) rPIE (r) PHeBIE

(s) isoTV ($b=10$) (t) AITV ($b=10$)

(u) AITV (full) (v) DR

(w) rPIE (x) PHeBIE

Figure 4: Reconstructions of the blind case for Gaussian noise.

Figure 5: Amplitude Gaussian metric plotted across 600 epochs for the blind
algorithms on the complex image given by Figure 2(d)-2(e).

For both the non-blind and blind case, we examine the case when
$\text{SNR}=40$ for the noisy measurements. We set $\lambda=10.0$. The
numerical results are recorded in Table 2. For both the non-blind and blind
cases, DR, rPIE, and PHeBIE yield magnitude images with the worst SSIM values,
and AITV outperforms its corresponding isotropic TV counterpart by having
better SSIM values for both the magnitude and phase images. The stochastic
AITV, particularly $b=10$ or $b=20$, has slightly lower magnitude SSIM values
by at most $0.04$ than the best results obtained from the deterministic, full-
batch AITV. In fact, stochastic AITV attains the second best magnitude SSIM
values in three out of the four cases considered. On the other hand,
stochastic AITV with either $b=10$ or $b=20$ has the best phase SSIM values,
outperforming its deterministic version by up to 0.19. For the blind case of
the cameraman/baboon image, the deterministic AITV and isotropic TV have worse
phase SSIM values than their stochastic counterparts. In general, the
stochastic algorithm does best in recovering phase images with superior SSIM
values while recovering the magnitude images with comparable SSIM values as
the deterministic algorithm.

Figure 6: Magnitude and phase SSIMs over different Gaussian noise level for
the complex image given by Figure 2(d)-2(e) for the blind case.

The reconstructed images for the non-blind experiments are presented in Figure
3. DR, rPIE, PHeBIE, and the stochastic algorithms have artifacts in all four
corners of the magnitude images because the corners are scanned significantly
less than in the middle of the image. However, the deterministic AITV has no
artifacts because (26) is solved exactly for the image solution $z$. As a
result, it has higher magnitude SSIM values than their stochastic
counterparts. Nevertheless, the stochastic algorithms yield better phase
images with less noise artifacts than any other algorithms. For example, the
phase images of Figure 2(e) reconstructed from stochastic isoTV and AITV have
the least amount of cameraman remnants.

Figure 5 shows the results of the blind algorithms. The phase images
reconstructed by the deterministic AITV, DR, ePIE and PHeBIE are significantly
worse than the stochastic algorithms. For example in Figure 2(b), the
contrasts of the reconstructed images by the deterministic AITV, DR, and
PHeBIE are inconsistent as they become darker from left to right while the
contrasts are more consistent with the stochastic algorithms. For Figure 2(e),
the stochastic algorithms perform the best in recovering the phase image while
deterministic AITV is unable to recover the left half of the image and DR,
rPIE, and PHeBIE have strong remnants of the cameraman present. Like in the
non-blind case, stochastic AITV reconstructs the phase image the best.

In Figure 5, we examine the convergence of the blind algorithms applied to the
Cameraman/Baboon image by recording their AGM values for each epoch. We omit
the convergence curves for isoTV since their curves are similar to their AITV
counterparts. Overall, the curves for our proposed stochastic algorithms are
decreasing, validating the numerical convergence of Algorithm 1 with AGM
fidelity. However, their curves are slightly above the deterministic ADMM
algorithm and rPIE. The reason why rPIE outperforms the AITV algorithms is
because it seeks to only minimize AGM while the AITV algorithms minimize a
larger objective function given by (9). Overall, after several hundred epochs,
our proposed stochastic algorithms can give comparable AGM values as the
deterministic AITV and rPIE algorithms.

Lastly, we analyze the robustness of the blind algorithms applied to Figures
2(d)-2(e) corrupted by different levels of Gaussian noise, from SNR 25 to 65.
The fidelity parameter $\lambda$ varies for different noise level of the
image: $\lambda=100$ for SNR = 25; $\lambda=50$ for SNR = 30, 35; $\lambda=10$
for SNR =40, 45; $\lambda=5$ for SNR = 50, 55, 60; and $\lambda=3$ for SNR =
65. The SSIMs for the magnitude and phase images across different SNRs are
plotted in Figure 6. For $\text{SNR}\geq 40$, the deterministic algorithms
have the best magnitude SSIMs than the other algorithms while their stochastic
counterparts have slightly lower SSIMs. When $\text{SNR}<40$, the stochastic
algorithms perform the best. In fact, stochastic AITV has magnitude SSIM at
least 0.90 across different noise levels. For the phase image, the stochastic
algorithms have the highest SSIMs up to SNR = 55. For $\text{SNR}\geq 60$, the
rPIE algorithm has the best phase SSIM while stochastic AITV has the second
best. Overall, stochastic AITV is the most stable across different levels of
Gaussian noise.

###  5.2 Poisson noise

Table 3: SSIM results of the algorithms applied to the Poisson corrupted measurements with $\eta=0.01$. The stochastic algorithms (e.g., AITV and isoTV, $b\in\\{5,10,20,25\\}$) are ran three times to obtain the average SSIM values. Bold indicates best value; underline indicates second best value. | Non-blind | Blind
---|---|---
| Chip | Cameraman/Baboon | Chip | Cameraman/Baboon
|  mag. SSIM |  phase SSIM |  mag. SSIM |  phase SSIM |  mag. SSIM |  phase SSIM |  mag. SSIM |  phase SSIM
DR | 0.8523 | 0.8455 | 0.8704 | 0.5043 | 0.8431 | 0.7387 | 0.7630 | 0.2529
PHeBIE | 0.9404 | 0.9398 | 0.9271 | 0.6791 | 0.9280 | 0.9082 | 0.8678 | 0.5470
isoTV $(b=5)$ | 0.9460 | 0.9151 | 0.9301 | 0.7193 | 0.9394 | 0.9001 | 0.9238 | 0.7105
isoTV $(b=10)$ | 0.9365 | 0.9212 | 0.9250 | 0.7085 | 0.9342 | 0.9064 | 0.9188 | 0.6979
isoTV $(b=20)$ | 0.9335 | 0.9397 | 0.9224 | 0.7008 | 0.9355 | 0.9248 | 0.9153 | 0.6905
isoTV $(b=25)$ | 0.9353 | 0.9465 | 0.9220 | 0.6992 | 0.9376 | 0.9315 | 0.9150 | 0.6907
isoTV (full batch) | 0.9767 | 0.9590 | 0.9773 | 0.7093 | 0.9655 | 0.9192 | 0.9588 | 0.4920
AITV $(b=5)$ | 0.9526 | 0.9680 | 0.9319 | 0.7477 | 0.9409 | 0.9530 | 0.9242 | 0.7418
AITV $(b=10)$ | 0.9590 | 0.9685 | 0.9375 | 0.7366 | 0.9472 | 0.9533 | 0.9321 | 0.7318
AITV $(b=20)$ | 0.9598 | 0.9682 | 0.9383 | 0.7171 | 0.9494 | 0.9526 | 0.9322 | 0.7114
AITV $(b=25)$ | 0.9585 | 0.9676 | 0.9373 | 0.7112 | 0.9492 | 0.9525 | 0.9307 | 0.7055
AITV (full batch) | 0.9803 | 0.9644 | 0.9782 | 0.7084 | 0.9741 | 0.9354 | 0.9671 | 0.4975

For both the non-blind and blind case, we examine the measurements corrupted
with Poisson noise with $\eta=0.01$ according to (48). We set $\lambda=0.15$.
The numerical results are recorded in Table 3. Note that rPIE results are
excluded because the algorithm is tailored towards measurements corrupted with
Gaussian noise [32]. Across all cases, deterministic AITV attains the highest
magnitude SSIM values and stochastic AITV attains the highest phase SSIM
values while DR performs the worst in reconstructing images from Poisson-
corrupted measurements. We observe general improvement in SSIM values for both
magnitude and phase images by using AITV over isoTV. Although the stochastic
algorithms have lower SSIM values than their deterministic counterparts for
the magnitude images, the difference is at most 0.047 for AITV and at most
0.056 for isoTV. Moreover, the SSIM values of the magnitude images from the
stochastic algorithms are at least 0.91. Similar to the Gaussian noise case,
stochastic AITV reconstructs the phase image well while recovering the
magnitude image with satisfactory quality.

We examine the robustness of the blind algorithms on Figures 2(d)-2(e) with
different level of Poisson noise. The noise levels we examine are
$\eta\in\\{0.005k\\}_{k=1}^{9}$. We set the fidelity parameter to be
$\lambda=15\times\eta$. The SSIMs for the magnitude and phase images across
different Poisson noise levels are plotted in Figure 7. We observe that the
deterministic algorithms yield the best magnitude SSIMs while the stochastic
algorithms yield the best phase SSIMs. DR yields the worst results for both
magnitude and phase components. Although stochastic AITV yields the third best
SSIMs for the magnitude image, its SSIMs are at least 0.90. Moreover, it has
the best phase SSIMs, significantly more than its deterministic counterpart by
at about 0.20. In summary, stochastic AITV is a robust method across different
levels of Poisson noise.

##  6 Conclusion

In this work, we propose AITV-regularized variational models for image
ptychography, where the measurements are corrupted by either Gaussian or
Poisson noise. To adapt the algorithm for large number of measurements, we
design a stochastic ADMM algorithm that incorporates adaptive step sizes based
on the PIE algorithms. Overall, using both AITV regularization and stochastic
ADMM, we are able to reconstruct an image of satisfactory quality from heavily
corrupted measurements as demonstrated in our numerical experiments. In fact,
the phase component of the image is best recovered through our proposed
algorithms. Lastly, we prove theoretical convergence for the proposed
stochastic ADMM algorithm under certain conditions and demonstrate numerical
convergence in our experiments.

Future directions include the design of a globally convergent algorithm for
the AITV-regularized ptychography model, and incorporation of variance-reduced
stochastic gradient estimators, such as SVRG [23] and SARAH [36], to
accelerate convergence and improve reconstruction quality.

Figure 7: Magnitude and phase SSIMs over different Poisson noise level for the
complex image given by Figure 2(d)-2(e) for the blind case.

##  Appendix A Proof of Lemma 3.1

######  Proof A.1.

If $x^{\prime}=0$, then it is trivial, so for the rest of the proof, we assume
that $x^{\prime}\neq 0$.

Suppose $x^{*}$ is the optimal solution to (24). We show that
$\operatorname{sgn}(x^{*})=\operatorname{sgn}(x^{\prime})$. If $x^{*}=0$, then
we can choose $c_{i}\in\\{c^{\prime}\in\mathbb{C}:|c^{\prime}|=1\\}$ such that
$\operatorname{sgn}(x^{*})_{i}=c_{i}=\operatorname{sgn}(x^{\prime})_{i}$ for
each $i=1,\ldots,n^{2}$, giving the desired result. Suppose that $x^{*}\neq
0$. Because $\|\cdot\|_{1}-\alpha\|\cdot\|_{2,1}$ is rotation invariant, we
only need to examine and expand the quadratic term in (24). We see that

$\displaystyle\|x^{*}-x^{\prime}\|_{2}^{2}$
$\displaystyle=\sum_{i=1}^{n^{2}}|(x^{*})_{i}-(x^{\prime})_{i}|^{2}=\sum_{i=1}^{n^{2}}\left(|(x^{*})_{i}|^{2}+|(x^{\prime})_{i}|^{2}-2|(x^{*})_{i}||(x^{\prime})_{i}|\cos\theta_{i}\right),$

where $\theta_{i}$ is the angle between the components $(x^{*})_{i}$ and
$(x^{\prime})_{i}$. This term is minimized when $\theta_{i}=0$ for all $i$.
This means that $\text{sgn}(x^{*})_{i}=\text{sgn}(x^{\prime})_{i}$ for all
$i$, or otherwise $x^{*}$ would not be an optimal solution to (24). Hence,
$\operatorname{sgn}(x^{*})=\operatorname{sgn}(x^{\prime})$.

After establishing that
$\operatorname{sgn}(x^{*})=\operatorname{sgn}(x^{\prime})$, we simplify (24)
to an optimization problem with respect to $|x|$ given by

$\displaystyle|x^{*}|=\operatorname*{arg\,min}_{\rho\in\mathbb{R}^{n},\;(\rho)_{i}\geq
0}\|\rho\|_{1}-\alpha\|\rho\|_{2}+\frac{1}{2\lambda}\|\rho-|x^{\prime}|\|_{2}^{2}.$

Therefore, by applying [27, Lemma 1] to the optimization problem followed by
multiplying the solution $|x^{*}|$ by $\operatorname{sgn}(x^{*})$, we obtain
the desired results.

##  Appendix B Proofs of Section 4

Before proving our main results, we present preliminary tools necessary for
the convergence analysis.

######  Definition B.1 ([41]).

Let $h:\mathbb{R}^{n^{2}}\rightarrow(-\infty,+\infty]$ be a proper and lower
semicontinuouous function and $\text{dom
h}\coloneqq\\{x\in\mathbb{R}^{n^{2}}:h(x)<\infty\\}$.

  1. (a)

The Fréchet subdifferential of $h$ at the point $x\in\text{dom }h$ is the set

$\displaystyle\hat{\partial}h(x)=\left\\{v\in\mathbb{R}^{n^{2}}:\liminf_{y\neq
x,y\rightarrow x}\frac{h(y)-h(x)-\langle v,y-x\rangle}{\|y-x\|}\geq
0\right\\}.$

  2. (b)

The limiting subdifferential of $h$ at the point $x\in\text{dom }h$ is the set

$\displaystyle\partial
h(x)=\left\\{v\in\mathbb{R}^{n^{2}}:\exists\\{(x^{t},y^{t})\\}_{t=1}^{\infty}\text{
s.t. }x^{t}\rightarrow x,h(x^{t})\rightarrow h(x),\;\hat{\partial}h(x^{t})\ni
y^{t}\rightarrow y\right\\}.$

We note that the limiting subdifferential is closed [41]:

$\displaystyle(x^{t},y^{t})\rightarrow(x,y),h(x^{t})\rightarrow
h(x),y^{t}\in\partial h(x^{t})\implies y\in\partial h(x).$

After establishing the definitions of subdifferentials in the real case, we
extend them to the complex case. If $z=z_{1}+z_{2}\mathbf{i}$, where
$z_{1},z_{2}\in\mathbb{R}^{n^{2}}$, then the limiting subdifferential of a
proper and lower semicontinuous function
$f:\mathbb{C}^{n^{2}}\rightarrow(-\infty,+\infty]$ is defined by

(49) $\displaystyle\partial
f(z)\coloneqq\partial_{z_{1}}f(z)+\partial_{z_{2}}f(z)\mathbf{i}.$

If the function $f$ is continuously differentiable at the point $z$, then with
slight abuse of notation, we denote its gradient by $\nabla f(z)$, and
$\partial f(z)=\\{\nabla f(z)\\}$ [41].

###  B.1 Proof of Lemma 4.1

######  Proof B.2.

If $\\{(\omega^{t},z^{t})\\}_{t=1}^{\infty}$ is bounded, then there exists a
constant $C$ such that $\|\omega^{t}\|_{\infty},\|z^{t}\|_{\infty}\leq C$ for
all $t\in\mathbb{N}$. We establish that $\mathcal{L}(\omega)$ has a Lipschitz
continuous gradient with respect to $\omega$. For any
$\omega_{1},\omega_{2}\in\mathbb{C}^{m^{2}}$ at iteration $t$, we have

$\displaystyle\left\|\nabla_{\omega}\mathcal{L}(\omega_{2})-\nabla_{\omega}\mathcal{L}(\omega_{1})\right\|_{2}\leq\left\|\sum_{j=1}^{N}\beta_{1}(S_{j}z^{t})^{*}\circ(\omega_{2}-\omega_{1})\circ(S_{j}z^{t})\right\|_{2}$
$\displaystyle\leq\left[\sum_{j=1}^{N}\beta_{1}\left\|(S_{j}z^{t})^{*}\circ(S_{j}z^{t})\right\|_{\infty}\right]\|\omega_{2}-\omega_{1}\|_{2}\leq\beta_{1}NC^{2}\|\omega_{2}-\omega_{1}\|_{2}.$

Hence, we observe that $\mathcal{L}(\omega)$ has a Lipschitz continuous
gradient with Lipschitz constant $L_{\omega}\coloneqq\beta_{1}NC^{2}$. By the
descent property [8, Definition 1], at iteration $t$ we have

$\displaystyle\mathcal{L}(\omega^{t+1})-\mathcal{L}(\omega^{t})$
$\displaystyle\leq\mathbb{R}(\langle\nabla_{\omega}\mathcal{L}(\omega^{t}),\omega^{t+1}-\omega^{t}\rangle)+\frac{L_{\omega}}{2}\|\omega^{t+1}-\omega^{t}\|_{2}^{2}$
$\displaystyle=-\delta_{\omega}^{t}\mathbb{R}(\langle\nabla_{\omega}\mathcal{L}(\omega^{t}),\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\rangle)+\frac{L_{\omega}(\delta^{t}_{\omega})^{2}}{2}\|\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\|_{2}^{2},$

where the last equality is due to (16). Taking the expectation with respect to
the first $t$ iterations, we obtain

$\displaystyle\mathbb{E}_{t}[\mathcal{L}(\omega^{t+1})]-\mathbb{E}_{t}[\mathcal{L}(\omega^{t})]$
$\displaystyle=-\delta_{\omega}^{t}\mathbb{R}\left(\mathbb{E}_{t}\left[\langle\nabla_{\omega}\mathcal{L}(\omega^{t}),\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\rangle\right]\right)+\frac{L_{\omega}(\delta^{t}_{\omega})^{2}}{2}\mathbb{E}_{t}\left[\|\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\|_{2}^{2}\right]$
$\displaystyle\leq-\left(K_{L}-\frac{\delta^{t}_{\omega}L_{\omega}(M_{V}+K_{U})}{2}\right)\delta^{t}_{\omega}\mathbb{E}_{t}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]+\frac{(\delta_{\omega}^{t})^{2}L_{\omega}M}{2},$

where the last inequality is due to combining (37), (38). and (41). Similarly,
we can estimate (44) because we can compute that $\mathcal{L}(z)$ has a
Lipschitz continuous gradient with Lipschitz constant
$L_{z}\coloneqq\beta_{1}NC^{2}+\beta_{2}\|\Delta\|$ and follow the same steps
as above.

###  B.2 Proof of Proposition 4.2

######  Proof B.3.

Because $\sum_{t=1}^{\infty}\|\Omega^{t+1}-\Omega^{t}\|_{2}^{2}<\infty$, we
have
$\displaystyle\lim_{t\rightarrow\infty}\Lambda_{j}^{t+1}-\Lambda_{j}^{t}=0$
for each $j=1,\ldots,N$ and
$\displaystyle\lim_{t\rightarrow\infty}y^{t+1}-y^{t}=0$, which implies from
(10e)-(10f) that

(50)
$\displaystyle\lim_{t\rightarrow\infty}u_{j}^{t}-\mathcal{F}(\omega^{t}\circ
S_{j}z^{t})$ $\displaystyle=0\;\quad\forall j=1,\ldots,N,$ (51)
$\displaystyle\lim_{t\rightarrow\infty}v^{t}-\nabla z^{t}$
$\displaystyle=0.\;$

It follows that $\\{(u^{t},v^{t})\\}_{t=1}^{\infty}$ is bounded since
$\\{(\omega^{t},z^{t})\\}_{t=1}^{\infty}$ is bounded. By
(LABEL:eq:u_closed_form), when $\mathcal{B}(\cdot,\cdot)$ is AGM, we have

$\displaystyle\|u_{j}^{t+1}\|_{2}=\left\|\frac{\sqrt{d_{j}}+\beta_{1}\left|\mathcal{F}(\omega^{t}\circ
S_{j}z^{t})-\frac{1}{\beta_{1}}\Lambda_{j}^{t}\right|}{1+\beta_{1}}\right\|_{2}\geq\frac{\beta_{1}}{1+\beta_{1}}\left(\frac{1}{\beta_{1}}\left\|\Lambda_{j}^{t}\right\|_{2}-\|\mathcal{F}(\omega^{t}\circ
S_{j}z^{t})\|_{2}\right),$

or equivalently,

(52)
$\displaystyle(1+\beta_{1})\|u_{j}^{t+1}\|_{2}+\beta_{1}\|\mathcal{F}(\omega^{t}\circ
S_{j}z^{t})\|_{2}\geq\|\Lambda_{j}^{t}\|_{2}.$

Similarly, when $\mathcal{B}(\cdot,\cdot)$ is IPM, we have the same inequality
as (52). As a result, $\\{\Lambda^{t}\\}_{t=1}^{\infty}$ is bounded. Finally,
we show that $\\{y^{t}\\}_{t=1}^{\infty}$ is bounded. By Lemma 3.1, we have
two cases. When

$\displaystyle\left\|(\nabla
z^{t})_{i}-\frac{(y^{t})_{i}}{\beta_{2}}\right\|_{\infty}\leq\frac{\lambda}{\beta_{2}},$

we have

$\displaystyle\frac{\lambda}{\beta_{2}}\geq\left\|(\nabla
z^{t})_{i}-\frac{(y^{t})_{i}}{\beta_{2}}\right\|_{\infty}\geq\frac{1}{\beta_{2}}\left\|(y^{t})_{i}\right\|_{\infty}-\left\|(\nabla
z^{t})_{i}\right\|_{\infty},$

or $\lambda+\beta_{2}\|(\nabla
z^{t})_{i}\|_{\infty}\geq\|(y^{t})_{i}\|_{\infty}.$ Otherwise, we have

$\displaystyle\|(v^{t+1})_{i}\|_{\infty}\geq\left\|\left|(\nabla
z^{t})_{i}-\frac{(y^{t})_{i}}{\beta_{2}}\right|-\frac{\lambda}{\beta_{2}}\right\|_{\infty}$
$\displaystyle\geq\left\|(\nabla
z^{t})_{i}-\frac{(y^{t})_{i}}{\beta_{2}}\right\|_{\infty}-\frac{\lambda}{\beta_{2}}$
$\displaystyle\geq\frac{1}{\beta_{2}}\left\|(y^{t})_{i}\right\|_{\infty}-\left\|(\nabla
z^{t})_{i}\right\|_{\infty}-\frac{\lambda}{\beta_{2}},$

or $\beta_{2}\left(\|(v^{t+1})_{i}\|_{\infty}+\left\|(\nabla
z^{t})_{i}\right\|_{\infty}\right)+\lambda\geq\left\|(y^{t})_{i}\right\|_{\infty}$.
Altogether, $\\{y^{t}\\}_{t=1}^{\infty}$ is bounded since
$\\{(v^{t},z^{t})\\}_{t=1}^{\infty}$ is bounded. Therefore, we establish that
$\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ is bounded.

We see that

$\displaystyle\mathcal{L}(Z,\Omega)=$
$\displaystyle\sum_{j=1}^{N}\left[\mathcal{B}(|u_{j}|^{2},d_{j})+\frac{\beta_{1}}{2}\left\|u_{j}-\mathcal{F}(\omega\circ
S_{j}z)+\frac{\Lambda_{j}}{\beta_{1}}\right\|_{2}^{2}-\frac{1}{2\beta_{1}}\|\Lambda_{j}\|_{2}^{2}\right]$
$\displaystyle+\lambda(\|v\|_{1}-\alpha\|v\|_{2,1})+\frac{\beta_{2}}{2}\left\|v-\nabla
z+\frac{y}{\beta_{2}}\right\|_{2}^{2}-\frac{1}{2\beta_{2}}\|y\|_{2}^{2}$
$\displaystyle\geq\sum_{j=1}^{N}\left[\mathcal{B}(|u_{j}|^{2},d_{j})-\frac{1}{2\beta_{1}}\|\Lambda_{j}\|_{2}^{2}\right]-\frac{1}{2\beta_{2}}\|y\|_{2}^{2}.$

Because $\mathcal{B}(\cdot,\cdot)$ is bounded below according to (6) and
$\\{\Omega^{t}\\}_{t=1}^{\infty}$ is bounded,
$\\{\mathcal{L}(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ is bounded below by some
constant $\mathcal{L}_{\inf}$. By (10a) and (10c), we have
$\mathcal{L}(u^{t+1})\leq\mathcal{L}(u^{t})$ and
$\mathcal{L}(v^{t+1})\leq\mathcal{L}(v^{t})$, respectively, so taking
expectation with respect to the first $t$ iterations, we obtain

(53)
$\displaystyle\mathbb{E}_{t}[\mathcal{L}(\omega^{t})]=\mathbb{E}_{t}[\mathcal{L}(u^{t+1})]$
$\displaystyle\leq\mathbb{E}_{t}[\mathcal{L}(u^{t})],$ (54)
$\displaystyle\mathbb{E}_{t}[\mathcal{L}(z^{t})]=\mathbb{E}_{t}[\mathcal{L}(v^{t+1})]$
$\displaystyle\leq\mathbb{E}_{t}[\mathcal{L}(v^{t})]=\mathbb{E}_{t}[\mathcal{L}(\omega^{t+1})].$

In addition, we have

$\displaystyle\mathcal{L}(\Lambda^{t+1})-\mathcal{L}(\Lambda^{t})$
$\displaystyle=\sum_{j=1}^{N}\mathbb{R}(\langle\Lambda_{j}^{t+1}-\Lambda_{j}^{t},u_{j}^{t+1}-\mathcal{F}(\omega^{t+1}\circ
S_{j}z^{t+1})\rangle)$
$\displaystyle=\frac{1}{\beta_{1}}\sum_{j=1}^{N}\left\|\Lambda_{j}^{t+1}-\Lambda_{j}^{t}\right\|_{2}^{2}=\frac{1}{\beta_{1}}\|\Lambda^{t+1}-\Lambda^{t}\|_{2}^{2},$

where the second to last equality is due to (10e). Taking expectation with
respect to the first $t$ iterations gives

(55)
$\displaystyle\mathbb{E}_{t}[\mathcal{L}(\Lambda^{t+1})]-\mathbb{E}_{t}[\mathcal{L}(\Lambda^{t})]=\frac{1}{\beta_{1}}\mathbb{E}_{t}\left[\left\|\Lambda^{t+1}-\Lambda^{t}\right\|_{2}^{2}\right].$

Similarly, we obtain

(56)
$\displaystyle\mathbb{E}_{t}[\mathcal{L}(y^{t+1})]-\mathbb{E}_{t}[\mathcal{L}(y^{t})]$
$\displaystyle=\frac{1}{\beta_{2}}\mathbb{E}_{t}[\|y^{t+1}-y^{t}\|_{2}^{2}].$

Summing up (43)-(44), (53)-(56) and taking total expectation, we have

(57) $\displaystyle\begin{aligned}
\mathbb{E}[\mathcal{L}(Z^{t+1},\Omega^{t+1})]&-\mathbb{E}[\mathcal{L}(Z^{t},\Omega^{t})]\leq\frac{1}{\beta_{1}}\mathbb{E}_{t}\left[\left\|\Lambda^{t+1}-\Lambda^{t}\right\|_{2}^{2}\right]+\frac{1}{\beta_{2}}\mathbb{E}[\|y^{t+1}-y^{t}\|_{2}^{2}]\\\
&-\left(K_{L}-\frac{\delta^{t}_{\omega}L_{\omega}(M_{V}+K_{U})}{2}\right)\delta^{t}_{\omega}\mathbb{E}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]+\frac{(\delta_{\omega}^{t})^{2}L_{\omega}M}{2}\\\
&-\left(K_{L}-\frac{\delta^{t}_{z}L_{z}(M_{V}+K_{U})}{2}\right)\delta^{t}_{z}\mathbb{E}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]+\frac{(\delta_{z}^{t})^{2}L_{z}M}{2}.\end{aligned}$

By (45), $\displaystyle\lim_{t\rightarrow\infty}\delta^{t}_{\omega}=0$ and
$\displaystyle\lim_{t\rightarrow\infty}\delta^{t}_{z}=0$, which means that we
can assume without generality that
$\delta^{t}_{\omega}L_{\omega},\delta^{t}_{z}L_{z}<\frac{K_{L}}{M_{V}+K_{U}}$
for all $t\in\mathbb{N}$. Hence, summing up $t=1,\ldots,T$, we obtain

$\displaystyle\mathcal{L}_{\inf}-\mathbb{E}[\mathcal{L}(Z^{0},\Omega^{0})]$
$\displaystyle\leq\mathbb{E}[\mathcal{L}(Z^{T+1},\Omega^{T+1})]-\mathbb{E}[\mathcal{L}(Z^{0},\Omega^{0})]$
$\displaystyle\leq\sum_{t=1}^{T}\Bigg{(}C\|\Omega^{t+1}-\Omega^{t}\|_{2}^{2}-\frac{K_{L}\delta_{\omega}^{t}}{2}\mathbb{E}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]-\frac{K_{L}\delta_{z}^{t}}{2}\mathbb{E}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]$
$\displaystyle+\frac{(\delta_{\omega}^{t})^{2}L_{\omega}M}{2}+\frac{(\delta_{z}^{t})^{2}L_{z}M}{2}\Bigg{)},$

where $C=\max\\{\frac{1}{\beta_{1}},\frac{1}{\beta_{2}}\\}$. Rearranging the
inequality and letting $T\rightarrow\infty$ give us

$\displaystyle 0\leq$
$\displaystyle\frac{K_{L}}{2}\sum_{t=1}^{\infty}\left(\delta_{\omega}^{t}\mathbb{E}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]+\delta_{z}^{t}\mathbb{E}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]\right)\leq$
$\displaystyle\mathbb{E}[\mathcal{L}(Z^{0},\Omega^{0})]-\mathcal{L}_{\inf}+\sum_{t=1}^{\infty}\Bigg{(}C\|\Omega^{t+1}-\Omega^{t}\|_{2}^{2}+\frac{(\delta_{\omega}^{t})^{2}L_{\omega}M}{2}+\frac{(\delta_{z}^{t})^{2}L_{z}M}{2}\Bigg{)}.$

By the assumption, the right-hand side is bounded, so therefore it implies
(46).

###  B.3 Proof of Theorem 4.3

######  Proof B.4.

If $\\{(w^{t},z^{t})\\}_{t=1}^{\infty}$ is bounded, then
$\|w^{t}\|_{2},\|z^{t}\|_{2}\leq C$ for all $t\in\mathbb{N}$ for some constant
$C>0$. At a given iteration $t$, we have
$\nabla_{\omega}\mathcal{L}(\omega^{t+1})=0$ by the first-order optimality
condition of (10b) and $\nabla_{\omega}\mathcal{L}(\omega)$ to be Lipschitz.
It follows that

$\displaystyle\|\nabla_{\omega}\mathcal{L}(\omega^{t})\|_{2}=\|\nabla_{\omega}\mathcal{L}(\omega^{t+1})-\nabla_{\omega}\mathcal{L}(\omega^{t})\|_{2}\leq
L_{\omega}\|\omega^{t+1}-\omega^{t}\|_{2}\leq 2CL_{\omega}.$

Similarly, we have $\|\nabla_{z}\mathcal{L}(z^{t})\|_{2}\leq 2CL_{z}.$ As a
result, by squaring the inequalities and taking their expectation,
$\left\\{\left(\mathbb{E}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right],\mathbb{E}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]\right)\right\\}_{t=1}^{\infty}$
is bounded. The step size condition (45) implies that
$\displaystyle\lim_{t\rightarrow\infty}\delta^{t}_{\omega}=0$ and
$\displaystyle\lim_{t\rightarrow\infty}\delta^{t}_{z}=0$. Since the SGD steps
for $\omega$ and $z$ are

$\displaystyle\omega^{t+1}=\omega^{t}-\delta_{\omega}^{t}\tilde{\nabla}\mathcal{L}(\omega^{t}),\;z^{t+1}=z^{t}-\delta_{z}^{t}\tilde{\nabla}\mathcal{L}(z^{t}),$

it follows that by taking expectation with respect to the first $t$ iterations
and using Assumption 4,

$\displaystyle\mathbb{E}_{t}[\left\|\omega^{t+1}-\omega^{t}\right\|_{2}^{2}]=(\delta_{\omega}^{t})^{2}\mathbb{E}_{t}[\|\tilde{\nabla}_{\omega}\mathcal{L}(\omega^{t})\|_{2}^{2}]\leq(\delta_{\omega}^{t})^{2}\left(M+(M_{V}+K_{U})\mathbb{E}_{t}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]\right),$
$\displaystyle\mathbb{E}_{t}[\left\|z^{t+1}-z^{t}\right\|_{2}^{2}]=(\delta_{z}^{t})^{2}\mathbb{E}_{t}[\|\tilde{\nabla}_{z}\mathcal{L}(z^{t})\|_{2}^{2}]\leq(\delta_{z}^{t})^{2}\left(M+(M_{V}+K_{U})\mathbb{E}_{t}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]\right).$

Because$\left\\{\left(\mathbb{E}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right],\mathbb{E}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]\right)\right\\}_{t=1}^{\infty}$
is bounded, we apply total expectation and take the limit to obtain

(58)
$\displaystyle\lim_{t\rightarrow\infty}\mathbb{E}[\left\|\omega^{t+1}-\omega^{t}\right\|_{2}^{2}]$
$\displaystyle\leq\lim_{t\rightarrow\infty}(\delta_{\omega}^{t})^{2}\left(M+(M_{V}+K_{U})\mathbb{E}\left[\left\|\nabla_{\omega}\mathcal{L}(\omega^{t})\right\|_{2}^{2}\right]\right)=0,$
(59)
$\displaystyle\lim_{t\rightarrow\infty}\mathbb{E}[\left\|z^{t+1}-z^{t}\right\|_{2}^{2}]$
$\displaystyle\leq\lim_{t\rightarrow\infty}(\delta_{z}^{t})^{2}\left(M+(M_{V}+K_{U})\mathbb{E}\left[\left\|\nabla_{z}\mathcal{L}(z^{t})\right\|_{2}^{2}\right]\right)=0.$

Earlier, in the proof of Proposition 4.2, we obtain

(60)
$\displaystyle\lim_{t\rightarrow\infty}u_{j}^{t}-\mathcal{F}(\omega^{t}\circ
S_{j}z^{t})$ $\displaystyle=0,\;\quad\forall j=1,\ldots,N,$ (61)
$\displaystyle\lim_{t\rightarrow\infty}v^{t}-\nabla z^{t}$ $\displaystyle=0.$

By Proposition 4.2, we have $\\{(Z^{t},\Omega^{t})\\}_{t=1}^{\infty}$ to be
bounded and (46) to be true.

Because $\\{Z^{t}\\}_{t=1}^{\infty}$ is bounded, there exists a convergent
subsequence $\\{(Z^{t_{k}},\Omega^{t_{k}})\\}_{k=1}^{\infty}$ and a point
$(Z^{\star},\Omega^{\star})$ such that
$\displaystyle\lim_{k\rightarrow\infty}(Z^{t_{k}},\Omega^{t_{k}})=(Z^{\star},\Omega^{\star})$.
In addition, (46) and (58)-(61) hold for the subsequence and its further
subsequences. (58)-(59) imply that there is a further subsequence
$\\{(Z^{t_{k_{\ell}}},\Omega^{t_{k_{\ell}}})\\}_{\ell=1}^{\infty}$ such that

(62)
$\displaystyle\lim_{\ell\rightarrow\infty}\omega^{t_{k_{\ell}}+1}-\omega^{t_{k_{\ell}}}$
$\displaystyle=0\;\text{a.s.},$ (63)
$\displaystyle\lim_{\ell\rightarrow\infty}z^{t_{k_{\ell}}+1}-z^{t_{k_{\ell}}}$
$\displaystyle=0\;\text{a.s.}$

From (46), we obtain
$\displaystyle\lim_{\ell\rightarrow\infty}\mathbb{E}\left[\delta_{\omega}^{t_{k_{\ell}}}\|\nabla_{\omega}\mathcal{L}(\omega^{t_{k_{\ell}}})\|_{2}^{2}+\delta_{z}^{t_{k_{\ell}}}\|\nabla_{z}\mathcal{L}(z^{t_{k_{\ell}}})\|_{2}^{2}\right]=0$,
so we obtain

(64)
$\displaystyle\liminf_{\ell\rightarrow\infty}\mathbb{E}[\|\nabla_{\omega}\mathcal{L}(\omega^{t_{k_{\ell}}})\|_{2}^{2}]$
$\displaystyle=0$ (65)
$\displaystyle\liminf_{\ell\rightarrow\infty}\mathbb{E}[\|\nabla_{z}\mathcal{L}(z^{t_{k_{\ell}}})\|_{2}^{2}]$
$\displaystyle=0.$

Now we show that $(Z^{\star},\Omega^{\star})$ is a stochastic KKT point a.s.
From (60)-(61), we have

(66) $\displaystyle u_{j}^{\star}$
$\displaystyle=\lim_{\ell\rightarrow\infty}u_{j}^{t_{k_{\ell}}}=\lim_{\ell\rightarrow\infty}\mathcal{F}(\omega^{{t_{k_{\ell}}}}\circ
S_{j}z^{t_{k_{\ell}}})=\mathcal{F}(\omega^{\star}\circ S_{j}z^{\star}),$ (67)
$\displaystyle v^{\star}$
$\displaystyle=\lim_{\ell\rightarrow\infty}v^{t_{k_{\ell}}}=\lim_{\ell\rightarrow\infty}\nabla
z^{t_{k_{\ell}}}=\nabla z^{\star}.$

By (60) and (62)-(63), we have

$\displaystyle\begin{aligned}
\lim_{\ell\rightarrow\infty}u^{t_{k_{\ell}}+1}&=\lim_{\ell\rightarrow\infty}\mathcal{F}\left(\omega^{t_{k_{\ell}}+1}\circ
S_{j}z^{t_{k_{\ell}}+1}\right)=\mathcal{F}\left(\lim_{\ell\rightarrow\infty}\omega^{t_{k_{\ell}}+1}\circ
S_{j}\left(\lim_{\ell\rightarrow\infty}z^{t_{k_{\ell}}+1}\right)\right)\\\
&=\mathcal{F}\left(\lim_{\ell\rightarrow\infty}\omega^{t_{k_{\ell}}}\circ
S_{j}\left(\lim_{\ell\rightarrow\infty}z^{t_{k_{\ell}}}\right)\right)=\lim_{\ell\rightarrow\infty}\mathcal{F}\left(\omega^{t_{k_{\ell}}}\circ
S_{j}z^{t_{k_{\ell}}}\right)=\lim_{\ell\rightarrow\infty}u^{t_{k_{\ell}}}\;\text{a.s.}\end{aligned}$

As a result, we have

(68)
$\displaystyle\lim_{\ell\rightarrow\infty}u^{t_{k_{\ell}}+1}-u^{t_{k_{\ell}}}=0\;\text{a.s.}$

Similarly, by (61) and (63), we have

(69)
$\displaystyle\lim_{\ell\rightarrow\infty}v^{t_{k_{\ell}}+1}-v^{t_{k_{\ell}}}=0\;\text{a.s.}$

For the sake of brevity, we will omit “a.s.” henceforth.

Next we prove (35a). Suppose that $\mathcal{B}(\cdot,\cdot)$ is AGM. At
iteration $t_{k_{\ell}}$, the first-order optimality condition of (10a) is

(70) $\displaystyle
0\in\partial\left|u_{j}^{t_{k_{\ell}}+1}\right|\left(|u_{j}^{t_{k_{\ell}}+1}|-\sqrt{d_{j}}\right)+\Lambda_{j}^{t_{k_{\ell}}}+\beta_{1}\left(u_{j}^{t_{k_{\ell}}+1}-\mathcal{F}(\omega^{t_{k_{\ell}}}\circ
S_{j}z^{k_{\ell}})\right),$

where

$\displaystyle\left(\partial\left|u_{j}^{t_{k_{\ell}}+1}\right|\right)_{i}=\begin{cases}\displaystyle\frac{(u_{j}^{t_{k_{\ell}}+1})_{i}}{|(u_{j}^{t_{k_{\ell}}+1})_{i}|},&\text{
if }(u_{j}^{t_{k_{\ell}}+1})_{i}\neq 0\\\ \\{u\in\mathbb{C}:|u|\leq
1\\}&\text{ if }(u_{j}^{t_{k_{\ell}}+1})_{i}=0.\end{cases}$

If $(u_{j}^{\star})_{i}\neq 0$ for some $i\in\\{1,\ldots,n^{2}\\}$, then there
exists a neighborhood
$B_{r}((u_{j}^{\star})_{i})=\\{u\in\mathbb{C}:|u-(u_{j}^{\star})_{i}|\leq
r\\}$ such that all points in $B_{r}((u_{j}^{\star})_{i})$ are nonzero. By
(68) and the fact that
$\displaystyle\lim_{\ell\rightarrow\infty}u^{t_{k_{\ell}}}=u^{\star}$, we have
$(u_{j}^{t_{k_{\ell}}+1})_{i}\neq 0$ for all $t_{k_{\ell}}$ sufficiently
large. As a result, (70) becomes

(71) $\displaystyle
0=\frac{(u_{j}^{t_{k_{\ell}}+1})_{i}}{|(u_{j}^{t_{k_{\ell}}+1})_{i}|}\left(|(u_{j}^{t_{k_{\ell}}+1})_{i}|-(\sqrt{d_{j}})_{i}\right)+(\Lambda_{j}^{t_{k_{\ell}}})_{i}+\beta_{1}\left((u_{j}^{t_{k_{\ell}}+1})_{i}-(\mathcal{F}(\omega^{t_{k_{\ell}}}\circ
S_{j}z^{k_{\ell}}))_{i}\right),$

By (60) and (68), we take the limit to obtain

(72) $\displaystyle
0=\frac{(u_{j}^{\star})_{i}}{|(u_{j}^{\star})_{i}|}\left(|(u_{j}^{\star})_{i}|-(\sqrt{d_{j}})_{i}\right)+(\Lambda_{j}^{\star})_{i}.$

On the other hand, when $(u_{j}^{\star})_{i}=0$ for some index $i$, we have
$(u_{j}^{t_{k_{\ell}}})_{i}\in
B_{r}((u_{j}^{\star})_{i})\setminus\\{(u_{j}^{\star})_{i}\\}$ to be nonzero.
Taking the absolute value of (71) gives us

$\displaystyle\left|(\Lambda_{j}^{t_{k_{\ell}}})_{i}+\beta_{1}\left((u_{j}^{t_{k_{\ell}}+1})_{i}-(\mathcal{F}(\omega^{t_{k_{\ell}}}\circ
S_{j}z^{k_{\ell}}))_{i}\right)\right|$
$\displaystyle=\left||(u_{j}^{t_{k_{\ell}}+1})_{i}|-(\sqrt{d_{j}})_{i}\right|.$

Again, by (60) and (68), taking the limit yields

$\displaystyle\left|(\Lambda_{j}^{\star})_{i}\right|=\left|(\sqrt{d_{j}})_{i}\right|,$

which implies that there exists $u^{\prime}\in\\{u\in\mathbb{C}:|u|\leq 1\\}$
such that $-u^{\prime}(\sqrt{d_{j}})_{i}+(\Lambda_{j}^{\star})_{i}=0$. This
result and (72) give us (35a) for the AGM case. As for the IPM case, because
$\displaystyle\lim_{x\rightarrow 0^{+}}x-d\log x=+\infty$, it follows that
$(u_{j}^{\star})_{i}\neq 0$ for all $i$, so we only need to worry about the
nonzero case. Hence, verifying (35a) for IPM requires computation similar to
the nonzero case for AGM, so it is omitted for brevity.

At iteration $t_{k_{\ell}}$, the first-order optimality condition of (10c) is

(73)
$\displaystyle-\frac{y^{t_{k_{\ell}}}}{\lambda}-\frac{\beta_{2}}{\lambda}\left(v^{t_{k_{\ell}}}-\nabla
z^{t_{k_{\ell}}}\right)\in\partial\left(\|v^{t_{k_{\ell}}+1}\|_{1}-\alpha\|v^{t_{k_{\ell}}+1}\|_{2,1}\right).$

By (61) and (69), we have

$\displaystyle\lim_{\ell\rightarrow\infty}-\frac{y^{t_{k_{\ell}}}}{\lambda}-\frac{\beta_{2}}{\lambda}\left(v^{t_{k_{\ell}}}-\nabla
z^{t_{k_{\ell}}}\right)=-\frac{y^{\star}}{\lambda},$
$\displaystyle\lim_{\ell\rightarrow\infty}v^{t_{k_{\ell}}+1}=\lim_{\ell\rightarrow\infty}v^{t_{k_{\ell}}}=v^{\star}.$

By continuity,
$\displaystyle\lim_{\ell\rightarrow\infty}\|v^{t_{k_{\ell}}+1}\|_{1}-\alpha\|v^{t_{k_{\ell}}+1}\|_{2,1}=\|v^{\star}\|_{1}-\alpha\|v^{\star}\|_{2,1}.$
Altogether, by closedness of limiting subdifferential, we establish (35b).

Lastly, we prove (35e′)-(35f′). At iteration $t_{k_{\ell}}$, we have

$\displaystyle\nabla_{\omega}\mathcal{L}(\omega^{t_{k_{\ell}}})=$
$\displaystyle-\sum_{j=1}^{N}\left\\{\beta_{1}(S_{j}z^{t_{k_{\ell}}})^{*}\circ\left[\mathcal{F}^{-1}\left(u_{j}^{{t_{k_{\ell}}}+1}+\frac{\Lambda_{j}^{t_{k_{\ell}}}}{\beta_{1}}\right)-\omega^{t_{k_{\ell}}}\circ
S_{j}z^{t_{k_{\ell}}}\right]\right\\}.$

Taking the limit, we see that

$\displaystyle\nabla_{\omega}\mathcal{L}(Z^{\star},\Omega^{\star})$
$\displaystyle=-\sum_{j=1}^{N}\left\\{\beta_{1}(S_{j}z^{\star})^{*}\circ\left[\mathcal{F}^{-1}\left(u_{j}^{\star}+\frac{\Lambda_{j}^{\star}}{\beta_{1}}\right)-\omega^{\star}\circ
S_{j}z^{\star}\right]\right\\}=\lim_{\ell\rightarrow\infty}\nabla_{\omega}\mathcal{L}(\omega^{t_{k_{\ell}}}).$

Applying Fatou’s Lemma to (64), we have
$\mathbb{E}\left[\|\nabla_{\omega}\mathcal{L}(Z^{\star},\Omega^{\star})\|_{2}^{2}\right]\leq\displaystyle\liminf_{\ell\rightarrow\infty}\mathbb{E}[\|\nabla_{\omega}\mathcal{L}(\omega^{t_{k_{\ell}}})\|_{2}^{2}]=0.$
Proving (35f′) is similar to proving (35e′), so we omit the proof. Altogether,
$(Z^{\star},\Omega^{\star})$ is a stochastic KKT point a.s.

## References

  * [1] H. H. Bauschke, P. L. Combettes, and D. R. Luke, Hybrid projection–reflection method for phase retrieval, Journal of the Optical Society of America A, 20 (2003), pp. 1025–1034.
  * [2] J. Bolte, S. Sabach, and M. Teboulle, Proximal alternating linearized minimization for nonconvex and nonsmooth problems, Mathematical Programming, 146 (2014), pp. 459–494.
  * [3] L. Bottou, F. E. Curtis, and J. Nocedal, Optimization methods for large-scale machine learning, SIAM Review, 60 (2018), pp. 223–311.
  * [4] S. Boyd, N. Parikh, and E. Chu, Distributed optimization and statistical learning via the alternating direction method of multipliers, Now Publishers Inc, 2011.
  * [5] K. Bui, F. Park, Y. Lou, and J. Xin, A weighted difference of anisotropic and isotropic total variation for relaxed Mumford–Shah color and multiphase image segmentation, SIAM Journal on Imaging Sciences, 14 (2021), pp. 1078–1113.
  * [6] E. J. Candes, X. Li, and M. Soltanolkotabi, Phase retrieval via Wirtinger flow: Theory and algorithms, IEEE Transactions on Information Theory, 61 (2015), pp. 1985–2007.
  * [7] E. J. Candes, T. Strohmer, and V. Voroninski, PhaseLift: Exact and stable signal recovery from magnitude measurements via convex programming, Communications on Pure and Applied Mathematics, 66 (2013), pp. 1241–1274.
  * [8] H. Chang, P. Enfedaque, and S. Marchesini, Blind ptychographic phase retrieval via convergent alternating direction method of multipliers, SIAM Journal on Imaging Sciences, 12 (2019), pp. 153–185.
  * [9] H. Chang, R. Glowinski, S. Marchesini, X.-C. Tai, Y. Wang, and T. Zeng, Overlapping domain decomposition methods for ptychographic imaging, SIAM Journal on Scientific Computing, 43 (2021), pp. B570–B597.
  * [10] H. Chang, Y. Lou, Y. Duan, and S. Marchesini, Total variation–based phase retrieval for Poisson noise removal, SIAM Journal on Imaging Sciences, 11 (2018), pp. 24–55.
  * [11] H. Chang, Y. Lou, M. K. Ng, and T. Zeng, Phase retrieval from incomplete magnitude information via total variation regularization, SIAM Journal on Scientific Computing, 38 (2016), pp. A3672–A3695.
  * [12] H. Chang, S. Marchesini, Y. Lou, and T. Zeng, Variational phase retrieval with globally convergent preconditioned proximal algorithm, SIAM Journal on Imaging Sciences, 11 (2018), pp. 56–93.
  * [13] Y. Chen and E. J. Candès, Solving random quadratic systems of equations is nearly as easy as solving linear systems, Communications on Pure and Applied Mathematics, 70 (2017), pp. 822–883.
  * [14] L. De Caro, D. Altamura, M. Arciniegas, D. Siliqi, M. R. Kim, T. Sibillano, L. Manna, and C. Giannini, Ptychographic imaging of branched colloidal nanocrystals embedded in free-standing thick polystyrene films, Scientific Reports, 6 (2016), pp. 1–8.
  * [15] P. Enfedaque, H. Chang, B. Enders, D. Shapiro, and S. Marchesini, High performance partial coherent x-ray ptychography, in International Conference on Computational Science, Springer, 2019, pp. 46–59.
  * [16] S. Esedoḡlu and S. J. Osher, Decomposition of images by the anisotropic Rudin-Osher-Fatemi model, Communications on Pure and Applied Mathematics, 57 (2004), pp. 1609–1626.
  * [17] A. Fannjiang and T. Strohmer, The numerics of phase retrieval, Acta Numerica, 29 (2020), pp. 125–228.
  * [18] P. Goyal, P. Dollár, R. Girshick, P. Noordhuis, L. Wesolowski, A. Kyrola, A. Tulloch, Y. Jia, and K. He, Accurate, large minibatch SGD: Training Imagenet in 1 hour, arXiv preprint arXiv:1706.02677, (2017).
  * [19] A. Greenbaum, Behavior of slightly perturbed lanczos and conjugate-gradient recurrences, Linear Algebra and its Applications, 113 (1989), pp. 7–63.
  * [20] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 770–778.
  * [21] R. Hesse, D. R. Luke, S. Sabach, and M. K. Tam, Proximal heterogeneous block implicit-explicit method and application to blind ptychographic diffraction imaging, SIAM Journal on Imaging Sciences, 8 (2015), pp. 426–457.
  * [22] M. R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving, Journal of research of the National Bureau of Standards, 49 (1952), p. 409.
  * [23] R. Johnson and T. Zhang, Accelerating stochastic gradient descent using predictive variance reduction, Advances in neural information processing systems, 26 (2013), pp. 315–323.
  * [24] Y.-F. Ke and C.-F. Ma, Alternating direction methods for solving a class of sylvester-like matrix equations, Linear and Multilinear Algebra, 65 (2017), pp. 2268–2292.
  * [25] H. Liu, K. Deng, H. Liu, and Z. Wen, An entropy-regularized admm for binary quadratic programming, Journal of Global Optimization, (2022), pp. 1–33.
  * [26] Q. Liu, X. Shen, and Y. Gu, Linearized admm for nonconvex nonsmooth optimization with convergence analysis, IEEE Access, 7 (2019), pp. 76131–76144.
  * [27] Y. Lou and M. Yan, Fast L1–L2 minimization via a proximal operator, Journal of Scientific Computing, 74 (2018), pp. 767–785.
  * [28] Y. Lou, P. Yin, Q. He, and J. Xin, Computing sparse representation in a highly coherent dictionary based on difference of ${L}_{1}$ and ${L}_{2}$, Journal of Scientific Computing, 64 (2015), pp. 178–196.
  * [29] Y. Lou, P. Yin, and J. Xin, Point source super-resolution via non-convex ${L}_{1}$ based methods, Journal of Scientific Computing, 68 (2016), pp. 1082–1100.
  * [30] Y. Lou, T. Zeng, S. Osher, and J. Xin, A weighted difference of anisotropic and isotropic total variation model for image processing, SIAM Journal on Imaging Sciences, 8 (2015), pp. 1798–1823.
  * [31] D. R. Luke, Relaxed averaged alternating reflections for diffraction imaging, Inverse Problems, 21 (2004), p. 37.
  * [32] A. Maiden, D. Johnson, and P. Li, Further improvements to the ptychographical iterative engine, Optica, 4 (2017), pp. 736–745.
  * [33] A. M. Maiden and J. M. Rodenburg, An improved ptychographical phase retrieval algorithm for diffractive imaging, Ultramicroscopy, 109 (2009), pp. 1256–1262.
  * [34] J. Marrison, L. Räty, P. Marriott, and P. O’Toole, Ptychography–a label free, high-contrast imaging technique for live cells using quantitative phase information, Scientific Reports, 3 (2013), pp. 1–7.
  * [35] Y. S. Nashed, D. J. Vine, T. Peterka, J. Deng, R. Ross, and C. Jacobsen, Parallel ptychographic reconstruction, Optics Express, 22 (2014), pp. 32082–32097.
  * [36] L. M. Nguyen, J. Liu, K. Scheinberg, and M. Takáč, Sarah: A novel method for machine learning problems using stochastic recursive gradient, in International Conference on Machine Learning, PMLR, 2017, pp. 2613–2621.
  * [37] Y. Ouyang, Y. Chen, G. Lan, and E. Pasiliao Jr, An accelerated linearized alternating direction method of multipliers, SIAM Journal on Imaging Sciences, 8 (2015), pp. 644–681.
  * [38] F. Park, Y. Lou, and J. Xin, A weighted difference of anisotropic and isotropic total variation for relaxed mumford-shah image segmentation, in 2016 IEEE International Conference on Image Processing (ICIP), IEEE, 2016, pp. 4314–4318.
  * [39] M. Pham, A. Rana, J. Miao, and S. Osher, Semi-implicit relaxed Douglas-Rachford algorithm (sDR) for ptychography, Optics Express, 27 (2019), pp. 31246–31260.
  * [40] H. Robbins and S. Monro, A stochastic approximation method, The annals of mathematical statistics, (1951), pp. 400–407.
  * [41] R. T. Rockafellar and R. J.-B. Wets, Variational analysis, vol. 317, Springer Science & Business Media, 2009.
  * [42] J. M. Rodenburg and H. M. Faulkner, A phase retrieval algorithm for shifting illumination, Applied Physics Letters, 85 (2004), pp. 4795–4797.
  * [43] L. I. Rudin, S. Osher, and E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D: nonlinear phenomena, 60 (1992), pp. 259–268.
  * [44] 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 processing magazine, 32 (2015), pp. 87–109.
  * [45] A. Suzuki, K. Shimomura, M. Hirose, N. Burdet, and Y. Takahashi, Dark-field x-ray ptychography: Towards high-resolution imaging of thick and unstained biological specimens, Scientific Reports, 6 (2016), pp. 1–9.
  * [46] P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer, Probe retrieval in ptychographic coherent diffractive imaging, Ultramicroscopy, 109 (2009), pp. 338–343.
  * [47] I. Waldspurger, A. d’Aspremont, and S. Mallat, Phase recovery, maxcut and complex semidefinite programming, Mathematical Programming, 149 (2015), pp. 47–81.
  * [48] A. Walther, The question of phase retrieval in optics, Optica Acta: International Journal of Optics, 10 (1963), pp. 41–49.
  * [49] Y. Wang, W. Yin, and J. Zeng, Global convergence of admm in nonconvex nonsmooth optimization, Journal of Scientific Computing, 78 (2019), pp. 29–63.
  * [50] Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, Image quality assessment: from error visibility to structural similarity, IEEE transactions on image processing, 13 (2004), pp. 600–612.
  * [51] Z. Wen, C. Yang, X. Liu, and S. Marchesini, Alternating direction methods for classical and ptychographic phase retrieval, Inverse Problems, 28 (2012), p. 115010.
  * [52] S. Wu Fung and Z. Di, Multigrid optimization for large-scale ptychographic phase retrieval, SIAM Journal on Imaging Sciences, 13 (2020), pp. 214–233.
  * [53] Y. Xu, W. Yin, Z. Wen, and Y. Zhang, An alternating direction algorithm for matrix completion with nonnegative factors, Frontiers of Mathematics in China, 7 (2012), pp. 365–384.
  * [54] P. Yin, Y. Lou, Q. He, and J. Xin, Minimization of $\ell_{1-2}$ for compressed sensing, SIAM Journal on Scientific Computing, 37 (2015), pp. A536–A563.
  * [55] P. Yin and J. Xin, PhaseLiftOff: an accurate and stable phase retrieval method based on difference of trace and frobenius norms, Communications in Mathematical Sciences, 13 (2014), pp. 1033–1049.
  * [56] Z. Yuan, H. Wang, and Q. Wang, Phase retrieval via sparse Wirtinger flow, Journal of Computational and Applied Mathematics, 355 (2019), pp. 162–173.
  * [57] H. Zhang, Y. Zhou, Y. Liang, and Y. Chi, A nonconvex approach for phase retrieval: Reshaped Wirtinger flow and incremental algorithms, Journal of Machine Learning Research, 18 (2017).
  * [58] L. Zhou, J. Song, J. S. Kim, X. Pei, C. Huang, M. Boyce, L. Mendonça, D. Clare, A. Siebert, C. S. Allen, et al., Low-dose phase retrieval of biological specimens using cryo-electron ptychography, Nature Communications, 11 (2020), pp. 1–9.