Alternating Projection, Ptychographic Imaging and Phase Synchronization

Introduction

The reconstruction of a scattering potential from measurements of scattered intensity in the far-field has occupied scientists and applied mathematicians for over a century, and arises in fields as varied as optics , astronomy , X-ray crystallography , tomographic imaging , holography , electron microscopy and particle scattering generally. Although phase-less diffraction measurements using short wavelength (such as X-ray, neutron, or electron wave packets) have been at the foundation of some of the most dramatic breakthrough in science - such as the first direct confirmation of the existence of atoms , the structure of DNA , RNA and over $100,000$ proteins or drugs involved in human life - the solution to the scattering problem for a general object was generally thought to be impossible for many years. Nevertheless, numerous experimental techniques that employ forms of interferometric/holographic measurements, gratings , and other phase mechanisms like random phase masks, sparsity structure, etc to help overcome the problem of phase-less measurements have been proposed over the years .

More recently an experimental technique has emerged that enables to image what no-one was able to see before: macroscopic specimens in 3D at wavelength (i.e. potentially atomic) resolution, with chemical state specificity. Ptychography was proposed in 1969 to improve the resolution in electron or x-ray microscopy by combining microscopy with scattering measurements. This technique enables one to build up very large images at wavelength resolution by combining the large field of view of a high precision scanning microscope system with the resolution enabled by diffraction measurements. In other words, the diffractive imaging and the scanning microscope techniques are combined together.

Initially, technological problems made ptychography impractical. Now, thanks to advances in source brightness and detector speed , research institutions around the world are rushing to develop hundreds of ptychographic microscopes to help scientists understand ever more complex nano-materials, self-assembled devices, or to study different length-scales involved in life, from macro-molecular machines to bones , and whenever observing the whole picture is as important as recovering local atomic arrangement of the components.

Experimentally, ptychography works by retrofitting a scanning microscope with a parallel detector. In a scanning microscope, a small beam is focused onto the sample via a lens, and the transmission is measured in a single-element detector. The image is built up by plotting the transmission as a function of the sample position as it is rastered across the beam. In such microscope, the resolution of the image is given by the beam size. In ptychography, one replaces the single element detector with a two-dimensional array detector such as a CCD and measures the intensity distribution at many scattering angles, much like a radar detector system for the microscopic world. Each recorded diffraction pattern contains short spatial Fourier frequency information about features that are smaller than the beam-size, enabling higher resolution. At short wavelengths however it is only possible to measure the intensity of the diffracted light. To reconstruct an image of the object, one needs to retrieve the phase. The phase retrieval problem is made tractable in ptychography by recording multiple diffraction patterns from the same region of the object, compensating phase-less information with a redundant set of measurements.

While reconstruction methods often work well in practice, fundamental mathematical questions concerning their convergence remain unresolved. The reader of an experimental paper is often left to wonder if the image and the resulting claims are valid, or one possibility among many solutions. Retractions of experimental results do happen (see for a discussion of controversial results in the optical community), and the problem is exacerbated because reproducing an image a nanoscale object is often not practical. What are often referred to as convergence results for projection algorithms are far from what we need for global convergence .

A popular algorithm for solving the phase retrieval problem was proposed in 1972. In their famous paper, Gerchberg and Saxton , independently of previous mathematical results for projections onto convex sets, proposed a simple algorithm for solving phase retrieval problems in two dimensions. In the algorithm was recognized as a projection algorithm that involves alternating projections between measurement space and object space. In 1982 Fienup generalized the Gerchberg-Saxton algorithm and analyzed many of its properties, showing, in particular, that the directions of the projections in the generalized Gerchberg-Saxton algorithm are formally similar to directions of steepest descent for a distance metric. One particular algorithm we focus on this paper is the alternating projection (AP) algorithm, which iteratively alternates between enforcing two pieces of information about the phase retrieval problem: the solution has known measured amplitude, and the illumination geometry is known. The main purpose of the AP algorithm is finding the solution that satisfies both conditions simultaneously.

Projection algorithms for convex sets have been well understood since 1960s. The phase retrieval problem, however, involves nonconvex sets. For this reason, the convergence properties of the Gerchberg-Saxton algorithm and its variants is still an open question except in very special cases .

There are two main results reported in this paper. We survey the relation between the AP algorithm and the uniqueness result shown in , and based on show that locally the stagnation set of the AP algorithm coincides with the unique solution up to a global phase factor in Theorem 3.16. With the help of the above results, in Theorem 3.18 we demonstrate the necessary and sufficient conditions of the local convergence of the AP algorithm to the unique solution up to a global phase factor. We show that the AP algorithm can fail to converge, in which case the step size can become arbitrarily small even though the limit is not a stagnation point. This issue has led to some confusion throughout the literature.

Second, we survey the intimate relationship between the ptychography imaging problem and the notion of phase synchronization. We form the connection graph and study the synchronization function of the ptychography imaging problem, which motivates the application of the recently developed technique graph connection Laplacian (GCL). In particular, in the ptychography imaging problem, phase synchronization based on GCL is applied to quickly construct an accurate initial guess for the AP algorithm to accelerate convergence speed for large scale diffraction data problems. With the help of the above results, in Section 5 we show some numerical results using different new algorithms. We also propose a new lens design and synchronization strategies that achieve over $80\times$ convergence rate and exhibit linear convergence. Numerical tests with noise exhibit linear relationship between the norm of the noise and and the final reconstruction error. While these numerical results are encouraging, they raise several questions and have practical implications, which we discuss in the conclusion.

The paper is organized as following. In Section 2 we introduce the ptychography experimental setup and notation. In Section 3 we show the necessary and sufficient conditions of the local convergence of AP. In addition, we discuss the relationship between the AP algorithm and optimization and show that the second derivative of the associated objective function is positive close to the solution. In Section 4 we discuss the relationship between the AP algorithm and the notion of phase synchronization, and propose methods based on GCL to obtain an accurate initial guess. In Section 5 we show numerical results of proposed methods and propose a new lens design and synchronization strategies that achieve over $40\times$ faster convergence than the AP algorithm and $10\times$ faster than the relaxed averaged alternating reflection (RAAR) algorithm.

Background and notations

2. The mathematical framework of the ptychography experiment

In a ptychography experiment, an object of interest is illuminated by a coherent beam, and the resulting diffraction pattern intensity is discretized by a pixellated camera. Numerically, the illuminated portion of the object is discretized to enable fast numerical methods. Such approximation is a valid representation of the physical experiment when the illumination function is smaller than the maximum bandwidth allowed by detector. We refer to to situations when these conditions are not strictly satisfied.

and similarly the support of $\psi$ is denoted as $\text{supp}(\psi)$.

With these notations, the relationship between the diffraction measurements collected in a ptychography experiment and $\psi$ can be represented compactly as

or $\bm{\mathsf{a}}=|\bm{\mathsf{FQ}}\psi^{\vee}|$, where

where $\alpha,\beta=0,\ldots,m-1$. The objective of the ptychographic reconstruction problem is to find $\psi$ given $\bm{\mathsf{a}}$ and the form (1).

The alternating projection algorithm and its convergence result

In this section, we describe the general phase retrieval problem and study the convergence of the alternating projection (AP) algorithm.

is it possible to recover $\psi_{0}$ from $\bm{\mathsf{a}}$? From now on, we assume that $\bm{\mathsf{a}}(i)\neq 0$ for all $i=1,\ldots,M$. Indeed, if there is any zero entry, we could remove the $i$-th vector $\bm{f}_{i}$ from the frame, as the phase information of the $i$-th component is not meaningful and we do not need to recover anything.

2. The alternating projection algorithm

that is, $P_{\bm{\mathsf{S}}}$ projects a complex vector to $R_{\bm{\mathsf{S}}}$.

Note that $(\bm{\mathsf{S}}^{*}\bm{\mathsf{S}})^{-1}$ exists since $\{\bm{f}_{i}\}_{i=1}^{M}$ is assumed to be a frame.

that is, $P_{\bm{\mathsf{a}}}$ substitutes the amplitude of $\bm{\mathsf{z}}(j)$ by $\bm{\mathsf{a}}(j)$ and preserve the phase informationNote that there are infinite different ways to define $P_{\bm{\mathsf{a}}}$ when $\bm{\mathsf{z}}$ has at least zero entries. Indeed, when the $i$-th entry of $\bm{\mathsf{z}}$ is zero, we could define the $i$-th entry of $P_{\bm{\mathsf{a}}}\bm{\mathsf{z}}$ to be $\bm{\mathsf{a}}(i)e^{i\theta}$, where $\theta\neq 0$. Here we focus on our definition for the sake of its simple appearance. Thus, we could view an entry with value as having the amplitude and phase and clearly $P_{\bm{\mathsf{a}}}$ is discontinuous at $\bm{\mathsf{z}}$ when there is at least one zero entry..

are both satisfied. Once we find the solution, the object of interest $\psi_{0}$ is estimated by

In the AP algorithm, the problem (5) is tackled by the following iterative scheme

3. Fundamental results

The main theorem in we count on is the following.

4. Some quantities and basic properties

To study the convergence behavior of the AP algorithm, we need the following definition.

The stagnation set (or the fixed points) of the AP algorithm when the given data is $\bm{\mathsf{a}}$ is defined as

Now, we can compare the definition of the stagnation set with the solution set of the phase retrieval problem. Clearly the solution set $S_{\bm{\mathsf{a}}}\subset\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}$. The stagnation set reflects the fact that $P_{\bm{\mathsf{a}}}\zeta-\zeta\neq 0$ does not imply $\zeta\neq P_{\bm{\mathsf{S}}}P_{\bm{\mathsf{a}}}\zeta$; that is, when $\zeta=P_{\bm{\mathsf{S}}}P_{\bm{\mathsf{a}}}\zeta$, $\zeta$ may or may not be the solution.

5. Some properties of the stagnation set

We take a closer look at the $\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}$ set. An immediate observation is the following co-dimension $1$ quantification of the stagnation set.

The stagnation set $\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}$ is of co-dimension $1$.

Suppose $\eta\in\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}$. By definition we have $\eta=P_{\bm{\mathsf{S}}}P_{\bm{\mathsf{a}}}\eta$. Thus we know $P_{\bm{\mathsf{S}}}(I-P_{\bm{\mathsf{a}}})\eta=0$ and hence $\eta^{*}(I-P_{\bm{\mathsf{a}}})\eta=0$. A direct expansion leads to $\sum_{k=1}^{M}(|\eta(k)|-\bm{\mathsf{a}}(k))|\eta(k)|=0$. Note that this equality is equivalent to the following

This equality leads to the co-dimension one conclusion. ∎

The equation (8) indicates that the non-negative real vector associated with $\eta$ in the stagnation set is on the sphere with the center $\bm{\mathsf{a}}/2$ and the radius $\|\bm{\mathsf{a}}\|_{2}/2$. Define

$\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}\cap Z$ is a closed subset on $Z\cap R_{\bm{\mathsf{S}}}$.

is a closed subset on $Z\cap R_{\bm{\mathsf{S}}}$. ∎

We know that the solution set $S_{\bm{\mathsf{a}}}$ is a closed $S^{1}$ set. We now show that the same geometric feature holds for a vector in the stagnation point when its all entries are non-zero.

If $\eta\in\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}\cap Z$, $e^{i\theta}\eta\in\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}\cap Z$ for all $\theta\in[0,2\pi)$.

When all entries of $\eta$ are non-zero, it is clear that $P_{\bm{\mathsf{a}}}e^{i\theta}\eta=e^{i\theta}P_{\bm{\mathsf{a}}}\eta$. Since $P_{\bm{\mathsf{S}}}$ is linear, we further conclude that $P_{\bm{\mathsf{S}}}P_{\bm{\mathsf{a}}}e^{i\theta}\eta=e^{i\theta}\eta$, which concludes the proof.

In this subsection, we take a closer look at the $I-P_{\bm{\mathsf{a}}}$ operator, which is related to the optimization approach discussed in Section 3.8.

On the other hand, we know that $I-P_{\bm{\mathsf{a}}}$ is not one-to-one. A quick observation of (11) is that when there is an entry in $\zeta$, we could find more than one $\eta$ so that $\eta-P_{\bm{\mathsf{a}}}\eta=\zeta$, and the more entries of $\zeta$ are zero, the more $\eta$ we could find. We now take a closer look at this one-to-one issue. Clearly by our definition, when $M=1$ so that $\bm{\mathsf{a}}=a>0$, we have $(I-P_{\bm{\mathsf{a}}})(0)=-a$. For the non-zero input to $I-P_{\bm{\mathsf{a}}}$, we have the following Lemma.

when $|\zeta|=0$, $\eta=ae^{i\theta}$, where $\theta\in[0,2\pi)$, are all solutions to $(I-P_{\bm{\mathsf{a}}})\eta=\zeta$.

When $\zeta(k)\geq\bm{\mathsf{a}}(k)$ for all $k$, there is a unique point in $\eta\in Z$ so that $(I-P_{\bm{\mathsf{a}}})\eta=\zeta$; that is, $I-P_{\bm{\mathsf{a}}}$ is one-to-one only on the set

Moreover, we have that $I-P_{\bm{\mathsf{a}}}$ is $2^{k}$-to-one on the set

and $I-P_{\bm{\mathsf{a}}}$ is infinite-to-one on the set

Clearly $Z=Y_{0}\cup\left(\cup_{k=1}^{M}Y_{k}\right)\cup Y_{\infty}$. Note the difference between $I-P_{\bm{\mathsf{a}}}$ and $P_{\bm{\mathsf{a}}}$ – $P_{\bm{\mathsf{a}}}$ is an infinity to one map. The results of Lemma 3.11, Lemma 3.12 and Corollary 3.1 are summarized in Figure 4, which illustrates the complicated behavior of the operator $I-P_{\bm{\mathsf{a}}}$.

where $b_{k}=|\zeta(k)|$. By the inner products $\eta_{1}^{*}(\eta_{1}-P_{\bm{\mathsf{a}}}\eta_{1})=\eta_{2}^{*}(\eta_{2}-P_{\bm{\mathsf{a}}}\eta_{2})=0$, we have

where $B:=\sum^{t_{3}}_{k=t_{2}+1}|\eta_{1}(k)|(|\eta_{1}(k)|-\bm{\mathsf{a}}(k))=\sum^{t_{3}}_{k=t_{2}+1}|\eta_{2}(k)|(|\eta_{2}(k)|-\bm{\mathsf{a}}(k))>0$ and $C:=\sum^{M}_{k=t_{3}+1}|\eta_{1}(k)|(\bm{\mathsf{a}}(k)-|\eta_{1}(k)|)=\sum^{M}_{k=t_{3}+1}|\eta_{2}(k)|(\bm{\mathsf{a}}(k)-|\eta_{2}(k)|)>0$. Here, the relationships in $B$ and $C$ come from (17). First, assume that $B-C\geq 0$. Then, by the relationship in (17), we have the inequalities

which is absurd. Similarly, if we have $B-C<0$, we use the inner products $\eta_{1}^{*}(\eta_{2}-P_{\bm{\mathsf{a}}}\eta_{2})=\eta_{2}^{*}(\eta_{1}-P_{\bm{\mathsf{a}}}\eta_{1})=0$ and get

Note that Corollary 3.1 and Lemma 3.13 do not imply that $\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}$ is in $Y_{0}$. It is possible that $\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}\subset Y_{k}$ such that $I-P_{\bm{\mathsf{a}}}$ is one-to-one. We have the following property restricting the stagnation set.

We could find $\epsilon>0$ small enough so that $(Z\cap B_{\epsilon}(0))\cap\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}=\emptyset$.

Take $\eta\in\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}\cap Z$ so that $0<|\eta(i)|<\epsilon\ll 1$ for all $i=1,\ldots,M$. By Lemma 3.12, we know that $\eta^{*}(\eta-P_{\bm{\mathsf{a}}}\eta)=0$ is equivalent to

where we denote $\zeta:=\eta-P_{\bm{\mathsf{a}}}\eta$ and $\pm$ depends on the possible $\eta$ associated with $\zeta$. Clearly $|\zeta(k)|=\bm{\mathsf{a}}(k)-|\eta(k)|<\bm{\mathsf{a}}(k)$, so $|\zeta(k)|-\bm{\mathsf{a}}(k)<0$ and $||\zeta(k)|-\bm{\mathsf{a}}(k)|<\epsilon$ for all $k$. Thus, we claim that (18) could not hold. If (18) holds, we should have $1\leq i_{1}<i_{2}<\ldots<i_{k}\leq M$ for some $1\leq k<M$ so that

Note that $\sum_{i=1}^{k}(|\zeta(i_{k})|+\bm{\mathsf{a}}(i_{k}))|\zeta(i_{k})|>0$ and $\sum_{j\neq i_{1},\ldots,i_{k}}(|\zeta(k)|-\bm{\mathsf{a}}(k))|\zeta(k)|<0$. While there are only finite possibilities of $1\leq i_{1}<i_{2}<\ldots<i_{k}\leq M$ for (19), we know that when $\epsilon$ is small enough, (19) does not hold. To be more precise, take $\sum_{i=1}^{M-1}(|\zeta(i)|+\bm{\mathsf{a}}(i))|\zeta(i)|=(\bm{\mathsf{a}}(M)-|\zeta(M)|)|\zeta(M)|$ as an example. Since $(\bm{\mathsf{a}}(M)-|\zeta(M)|)|\zeta(M)|<\epsilon\bm{\mathsf{a}}(M)$, when $\epsilon$ is small enough, $\sum_{i=1}^{M-1}(|\zeta(i)|+\bm{\mathsf{a}}(i))|\zeta(i)|=(\bm{\mathsf{a}}(M)-|\zeta(M)|)|\zeta(M)|$ fails. ∎

7. Convergence of the AP algorithm

In this subsection, we show an if and only if condition for the local convergence of the AP algorithm. Recall that we assume without loss of generality that $S_{\bm{\mathsf{a}}}\subset Z$.

where the equality holds when $w=P_{\bm{\mathsf{a}}}\zeta$.

where the equality holds when $z=P_{\bm{\mathsf{S}}}w$.

For all nonzero $\zeta\in R_{\bm{\mathsf{S}}}$, $P_{\bm{\mathsf{a}}}\zeta$ is not perpendicular to $R_{\bm{\mathsf{S}}}$.

When $M\geq 4N-2$ and $R_{\bm{\mathsf{S}}}$ generic, given $\zeta\in R_{\bm{\mathsf{S}}}$, $P_{\bm{\mathsf{a}}}\zeta\in R_{\bm{\mathsf{S}}}$ holds if and only if $\zeta\in S_{\bm{\mathsf{a}}}$.

The proof of (b) is directly from the fact the $P_{\bm{\mathsf{S}}}$ is a projection operator.

For (c), denote $\zeta=(b_{i}e^{i\theta_{i}})_{i=1}^{M}\in R_{\bm{\mathsf{S}}}\backslash\{0\}$, where $b_{i}\geq 0$ and $\theta_{i}\in[0,2\pi)$. Suppose $b_{i}>0$ for all $i$. Then by definition $P_{\bm{\mathsf{a}}}\zeta=(a_{i}e^{i\theta_{i}})_{i=1}^{M}$. Then it is clear that $\langle P_{\bm{\mathsf{a}}}\zeta,\zeta\rangle>0$, which shows the claim. When $b_{i}=0$ for some $i$, the $i$-th term does not contribute to $\langle P_{\bm{\mathsf{a}}}\zeta,\zeta\rangle$ and hence the argument holds.

The statement (d) is direct from Theorem 3.4.

The following theorem states the local convergence of the AP algorithm.

When $M\geq 4N-2$ and $R_{\bm{\mathsf{S}}}$ is generic, we could find an open neighborhood $U_{\bm{\mathsf{S}}}$ of $S_{\bm{\mathsf{a}}}$ so that $U_{\bm{\mathsf{S}}}\cap\Theta^{\textup{AP}}_{\bm{\mathsf{a}},0}=S_{\bm{\mathsf{a}}}$.

The AP algorithm can be studied in the non-convex optimization framework . Given a set of subsets $S_{i}$, $i=1,\ldots,L$ of a metric space $X$ so that $S:=\cap_{i=1}^{L}S_{i}\neq\emptyset$. To find $S$, we may consider the proposed sequence of successive projections (SOSP) scheme, which successively project the estimator to $S_{i}$. When the initial value $x_{0}$ of the SOSP $\{x_{n}\}_{n\geq 0}$ is a point of attraction [23, Definition 4.4] of an ordered collection of proximal sets in a metric space whose intersection $S$ is not empty, then either $\{x_{n}\}_{n\geq 0}$ converges to a point in $S$ or the set of the cluster points of $\{x_{n}\}_{n\geq 0}$ is a nontrivial continuum in $S$ [23, Theorem 4.3].

Here $\{\alpha_{l},\beta_{l}\}$ depend on $\bm{\mathsf{S}}$ and $\bm{\mathsf{a}}$. In particular, when $M\geq 4N-2$, $R_{\bm{\mathsf{S}}}$ generic and $\zeta^{(0)}\in U_{\bm{\mathsf{S}}}$, $\alpha_{l}<1$ and $\beta_{l}<1$. Moreover, if we denote $\zeta^{(l)}=(b_{k}^{(l)}e^{i\bm{\phi}^{(l)}_{k}})_{k=1}^{M}$, where $\bm{\phi}^{(l)}_{k}\in[0,2\pi)$ when $b_{k}^{(l)}>0$ and $\bm{\phi}^{(l)}_{k}=0$ when $b_{k}^{(l)}=0$, the following inequality holds:

Based on Lemma 3.15(a), we have the following inequalities. First,

due to Lemma 3.15 (a); by Lemma 3.15 (b), we have

When $M\geq 4N-2$ and $\zeta^{(0)}\in U_{\bm{\mathsf{S}}}$, the equality can not hold since $\zeta^{(l-1)}\notin\Theta_{\bm{\mathsf{a}},0}^{\text{AP}}$ due to Theorem 3.16. Similarly, by Lemma 3.15, we have (20). Now, since $\zeta^{(l)}=(b_{k}^{(l)}e^{i\bm{\phi}^{(l)}_{k}})_{k=1}^{M}$, we have

The equations (20) and (21) imply monotonic decrease and the equation (22) relates the phase step with the decrease in equation (20). We mention that (20) and (21), which are also shown in , do not imply convergence to the solution nor to a stagnation point. Also note that (20) and (21) do not imply

Indeed, note that $P_{\bm{\mathsf{a}}}\zeta^{(l)}-\zeta^{(l+1)}$ is perpendicular to $\zeta^{(l+1)}-\zeta^{(l)}$. Thus we have

where when (20) and (21) hold, it is still possible that $\|(P_{\bm{\mathsf{S}}}P_{\bm{\mathsf{a}}}-I)\zeta^{(l)}\|>\|(P_{\bm{\mathsf{S}}}P_{\bm{\mathsf{a}}}-I)\zeta^{(l-1)}\|$. See Figure 12 in the numerical section for an example. We finally come to our main Theorem regarding the if and only if condition of the local convergence of the AP algorithm.

When $M\geq 4N-2$ and $R_{\bm{\mathsf{S}}}$ generic, the following three conditions are equivalent when the initial point is inside $U_{\bm{\mathsf{S}}}$:

AP algorithm converges to the solution set ;

$\|(P_{\bm{\mathsf{a}}}-I)\zeta^{(l)}\|\to 0$ ;

$\|(P_{\bm{\mathsf{S}}}-I)\zeta^{(l+1/2)}\|\to 0$ .

so (1) implies (2). Similarly, we have (1) implies (3) since

due to the fact that $(P_{\bm{\mathsf{S}}}-I)$ is continuous. In addition, since $P_{\bm{\mathsf{S}}}$ is a projection operator, we have

Next, we show (2) implies (3) and (4). Note that when $\|(P_{\bm{\mathsf{a}}}-I)\zeta^{(l)}\|\to 0$, we have $\|\zeta^{(l+1)}-\zeta^{(l)}\|\to 0$ and $\|(P_{\bm{\mathsf{S}}}-I)P_{\bm{\mathsf{a}}}\zeta^{(l)}\|\to 0$ by (23).

Finally, we show that that (3) implies (1). Since $P_{\bm{\mathsf{S}}}\zeta^{(l+1/2)}\in R_{\bm{\mathsf{S}}}$, (3) means $\zeta^{(l+1/2)}\in P_{\bm{\mathsf{a}}}$ converges to a point located on $R_{\bm{\mathsf{S}}}\cap T_{\bm{\mathsf{a}}}$; that is, $\zeta^{(l+1/2)}$ converges to the solution set when the initial point is inside $U_{\bm{\mathsf{S}}}$. Thus we have finished the claim that (1), (2) and (3) are equivalent.

and hence the convergence. Second, suppose $\liminf_{l\to\infty}\alpha_{l}=1$, that is, $\lim_{l\to\infty}\alpha_{l}=1$. Clearly the series $p_{n}:=\Pi_{l=1}^{n}\alpha_{l}$ converges as $n\to\infty$ since $\alpha_{l}<1$. If the infinite product $\Pi_{l=1}^{\infty}\alpha_{l}$ diverges to , the AP algorithm converges to the solution, but at a slow rate, which might be as slow as possible. Note that $\Pi_{l=1}^{\infty}\alpha_{l}$ converges if and only if the series $\sum_{l=1}^{\infty}(1-\alpha_{l})$ converges.

8. The Relationship between the AP Algorithm and Optimization

To better understand the AP algorithm, we assume $M\geq 4N-2$ in this section. Define an objective function

Note that we take the transpose since $r(\bm{\mathsf{z}})$ is a real vector. The objective function $\rho$, when restricted on $R_{\bm{\mathsf{S}}}$, gauges how far we are to the solution. Recall that the solution is located on $R_{\bm{\mathsf{S}}}\cap Z$ by assumption. To evaluate the gradient and Hessian of $\rho$, we prepare the following calculations . First, we evaluate the derivative of $r(\bm{\mathsf{z}})$ with respect to $\bm{\mathsf{z}}$ at $\bm{\mathsf{z}}\in Z$:

Thus, by the chain rule we obtain the derivative of $\rho(\bm{\mathsf{z}})$ with respect to $\bm{\mathsf{z}}$ and $\overline{\bm{\mathsf{z}}}$ at $\bm{\mathsf{z}}\in Z$:

Next we evaluate the following quantities evaluated at $\bm{\mathsf{z}}$:

The Hessian of $\rho$ at $\bm{\mathsf{z}}$, denoted by $\nabla^{2}\rho|_{\bm{\mathsf{z}}}$, by a direct calculation is given by

which leads to the following evaluation of the curvature of the $\rho$. Take $\bm{\mathsf{w}}\in Z$. Denote $\bm{\mathsf{w}}=(b_{i}e^{i\theta_{i}})_{i=1}^{M}$, where $\theta_{i}\in[0,2\pi)$ when $b_{i}>0$ and $\theta_{i}=0$ when $b_{i}=0$. Then by a direct expansion, the second derivative of $\rho$ in the direction $\bm{\mathsf{w}}$ at $\bm{\mathsf{z}}$ is

We have the following observations about the gradient and Hessian:

Note that we can view the AP algorithm as the projected gradient descent algorithm related to the objective function $\rho$ . Indeed, we have

when $\zeta^{(l)}\in U_{\bm{\mathsf{F}}}\backslash S_{\bm{\mathsf{a}}}$. By (43), for $\zeta^{(l)}=\bm{\mathsf{b}}^{(l)}e^{i\bm{\phi}^{(l)}}$ we have

By Lemma 3.6, for a generic $R_{\bm{\mathsf{S}}}$, the gradient of $\rho$ on $R_{\bm{\mathsf{S}}}$ is zero only at $S_{\bm{\mathsf{a}}}$ since the only points on $R_{\bm{\mathsf{S}}}$ that have modulations $\bm{\mathsf{a}}$ are the points in the solution set. Also, by Theorem 3.16 when $\zeta^{(l)}\in U_{\bm{\mathsf{S}}}\backslash S_{\bm{\mathsf{a}}}$, $\nabla\rho|_{\zeta^{(l)}}$ is not perpendicular to $R_{\bm{\mathsf{S}}}$, since $P_{\bm{\mathsf{S}}}\nabla\rho|_{\zeta^{(l)}}=P_{\bm{\mathsf{S}}}(I-P_{\bm{\mathsf{a}}})\zeta^{(l)}\neq 0$ on $U_{\bm{\mathsf{S}}}\backslash S_{\bm{\mathsf{a}}}$. Furthermore, when $\zeta^{(l)}\notin S_{\bm{\mathsf{a}}}$, $\nabla\rho|_{\zeta^{(l)}}$ does not locate on $R_{\bm{\mathsf{S}}}$. Indeed, if $\nabla\rho|_{\zeta^{(l)}}\in R_{\bm{\mathsf{S}}}$, then $\zeta^{(l+1)}=\zeta^{(l)}-(I-P_{\bm{\mathsf{a}}})\zeta^{(l)}=P_{\bm{\mathsf{a}}}\zeta^{(l)}$; that is, $P_{\bm{\mathsf{a}}}\zeta^{(l)}\in R_{\bm{\mathsf{S}}}$ and hence $\zeta^{(l)}\in S_{\bm{\mathsf{a}}}$.

For $\bm{\mathsf{z}}=e^{it}\bm{\mathsf{S}}\psi_{0}=\bm{\mathsf{a}}e^{i(\bm{\phi}^{\bm{\mathsf{a}}}+t)}\in S_{\bm{\mathsf{a}}}$, for some $t\in[0,2\pi)$, and $\bm{\mathsf{w}}=\bm{\mathsf{b}}e^{i\bm{\theta}}\neq 0$, by (53) we know

which is always non-negative since $\sin^{2}\leq 1$. When $\bm{\theta}=\bm{\phi}^{\bm{\mathsf{a}}}+t+\pi/2$, $\nabla^{2}\rho|_{\bm{\mathsf{z}}}(\bm{\mathsf{w}})=0$.

The ptychography imaging problem and phase synchronization

In this section, we focus ourselves on the ptychography problem – how to find a good initial value for the iterative algorithm like AP, so that we could have a convergence result and speed up the algorithm. To simplify the discussion, we assume that $\text{supp}(\omega)=D_{r}^{m}$. A general setup can be easily adapted to $\text{supp}(\omega)\subsetneqq D_{r}^{m}$ We make the following assumption about the illumination scheme:

The chosen illumination scheme $\mathcal{X}_{K}$ satisfies the following two conditions

$\bm{\mathsf{x}}_{i}\neq\bm{\mathsf{x}}_{j}$ for all $i\neq j$;

$\mathcal{X}_{K}$ is ordered so that $\cup_{i=1}^{l}\iota_{\bm{\mathsf{x}}_{i}}(\text{supp}(\omega))\subsetneqq\cup_{i=1}^{l+1}\iota_{\bm{\mathsf{x}}_{i}}(\text{supp}(\omega))$, where $l=1,\ldots K-1$, $\cup_{i=1}^{K-1}\iota_{\bm{\mathsf{x}}_{i}}(\text{supp}(\omega))\subsetneqq D_{r}^{n}$ and $\cup_{i=1}^{K}\iota_{x_{i}}(\text{supp}(\omega))=D_{r}^{n}$;

For each $i$, there exists $j$ so that $\iota_{\bm{\mathsf{x}}_{i}}(\text{supp}(\omega))\cap\iota_{\bm{\mathsf{x}}_{j}}(\text{supp}(\omega))\neq\emptyset$.

The third assumption essentially says that each subregion is overlapped by at least one other subregion so that there is a channel for these subregions to “exchange information”.

Given $\mathcal{X}_{K}$, the object of interest $\psi$ is connected with respect to $\mathcal{X}_{K}$.

Given $\bm{\mathsf{a}}=|\bm{\mathsf{FQ}}\psi|$, we combine the essences of the AP algorithm and consider the following optimization problem:

which is a 1 to 1 map providing the index of the entry $\bm{\mathsf{r}}_{k}$ of the $k$-th illumination window $\iota_{\bm{\mathsf{x}}_{k}}(D_{r}^{m})$ in the long stack vector. Recall that $\texttt{r}_{m}$ is defined in (2) and $\bm{\mathsf{r}}_{k}$ and $D_{r}^{m}$ are defined in Section 2.2. For $j=1,\ldots,K$ and $\bm{\mathsf{s}}\in D_{r}^{m}$, define a set

which contains the indices of all illumination windows covering $\bm{\mathsf{x}}_{j}+\bm{\mathsf{s}}$. Also define a subset of $D_{r}^{m}$

which collects the indices of the pixels in all illumination windows which cover $\bm{\mathsf{x}}_{j}+\bm{\mathsf{s}}$. We choose to use this seeming complicated index since we would like to make clear the relationship between the illumination windows and their pixels. By Assumption 4.1 and a direct calculation, we know that $\bm{\mathsf{Q}}^{*}\bm{\mathsf{Q}}$ is a $n^{2}\times n^{2}$ non-degenerate diagonal matrix describing how many illumination windows cover a given pixel of the object of interest, where the $\texttt{r}_{n}(\bm{\mathsf{x}}_{j}+\bm{\mathsf{r}}_{j})$-th diagonal entry is $\sum_{\bm{\mathsf{r}}\in J_{\bm{\mathsf{x}}_{j},\bm{\mathsf{r}}_{j}}}|\omega(\bm{\mathsf{r}})|^{2}$. So, the matrix $P_{\bm{\mathsf{Q}}}:=\bm{\mathsf{Q}}(\bm{\mathsf{Q}}^{*}\bm{\mathsf{Q}})^{-1}\bm{\mathsf{Q}}^{*}$ satisfies

where $(i,\bm{\mathsf{r}}_{i})\sim(j,\bm{\mathsf{r}}_{j})$ means all illumination windows covering the pixel $\iota_{\bm{\mathsf{x}}_{i}}(\bm{\mathsf{r}}_{i})$. Geometrically, $P_{\bm{\mathsf{Q}}}$ describes how two illumination windows in the spatial domain are intersected and how the overlapped pixels are related via the illuminating function $\omega$. Note that when $\zeta_{(i)}$ contains the right amplitude and phase, $F^{*}\zeta_{(i)}$ is the correct image on $\iota_{\bm{\mathsf{x}}_{i}}(D_{r}^{m})$. Thus, maximizing $\zeta^{*}P_{\bm{\mathsf{FQ}}}\zeta$ is equivalent to requiring that the images on a pair of overlapping illumination windows match in the overlapping region. In particular, by Assumption 4.1, phases on one illumination window will be synchronized with at least one different illumination window if we maximize $\zeta^{*}P_{\bm{\mathsf{FQ}}}\zeta$. Also, by Assumption 4.3, the phases in different disconnected regions of $\psi$ associated with $\mathcal{X}_{K}$ are guaranteed to interact with each other so that the phase can be synchronized in the end.

To better understand (55), we further consider the relationship between the phases when the illumination windows overlap. We start from studying the Hermitian matrix $P_{\bm{\mathsf{FQ}}}$ in (55). The amplitude information, $\bm{\mathsf{a}}$, will be taken into account later. Consider the following phase synchronization problem:

Denote $\mathsf{O}_{ij}:=\iota_{\bm{\mathsf{x}}_{i}}(D_{r}^{m})\cap\iota_{\bm{\mathsf{x}}_{j}}(D_{r}^{m})$ to be the overlap of two illumination windows. A direct expansion of (58) leads to

where $\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}:=\bm{\mathsf{x}}_{i}-\bm{\mathsf{x}}_{j}$ and the last equality comes from the fact that $\bm{\mathsf{T}}_{\bm{\mathsf{x}}_{i}}\bm{\mathsf{T}}^{*}_{\bm{\mathsf{x}}_{j}}=\bm{\mathsf{T}}_{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}$ and $\bm{\mathsf{T}}^{*}_{\bm{\mathsf{x}}_{i}}\bm{\mathsf{T}}_{\bm{\mathsf{x}}_{i}}=I$. Clearly if $\mathsf{O}_{ij}=\emptyset$, $\bm{\mathsf{T}}_{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}\bm{\mathsf{R}}^{*}$ is a zero matrix. Note that $\bm{\mathsf{T}}_{\bm{\mathsf{x}}_{i}}(\bm{\mathsf{Q}}^{*}\bm{\mathsf{Q}})^{-1}\bm{\mathsf{T}}^{*}_{\bm{\mathsf{x}}_{i}}$, as the conjugation of $(\bm{\mathsf{Q}}^{*}\bm{\mathsf{Q}})^{-1}$ by $\bm{\mathsf{T}}_{\bm{\mathsf{x}}_{i}}$, is diagonal. It actually translates the $\texttt{r}_{n}(\bm{\mathsf{x}}_{i})$-th diagonal entry to the $1$-st diagonal entry. Also note that the overlapping information about the $i$-th and $j$-th illumination windows is preserved in $\bm{\mathsf{T}}_{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}\bm{\mathsf{R}}^{*}$.

Now we move $\bm{\mathsf{T}}_{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}$ out of $F{\text{diag}}(w)\bm{\mathsf{R}}\bm{\mathsf{T}}_{\bm{\mathsf{x}}_{i}}(\bm{\mathsf{Q}}^{*}\bm{\mathsf{Q}})^{-1}\bm{\mathsf{T}}^{*}_{\bm{\mathsf{x}}_{i}}\bm{\mathsf{T}}_{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}\bm{\mathsf{R}}^{*}{\text{diag}}(w^{*})F^{*}$ by a direct expansion:

where $\bm{\mathsf{M}}^{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}$ is a $m^{2}\times m^{2}$ masking matrix which is diagonal and depends on ${\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}$:

and $\mathsf{D}_{ij}:=\iota_{(0,0)}D_{r}^{m}\cap[\bm{\mathsf{T}}_{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}\iota_{(0,0)}D_{r}^{m}]$. This equality indicates the influence of the restriction matrix $\bm{\mathsf{R}}$ – the non-overlapped parts of the two overlapping subregions cannot be eliminated. Next, for $\bm{\mathsf{r}}_{i},\bm{\mathsf{r}}_{j}\in D_{r}^{m}$, when $\mathsf{O}_{ij}\neq\emptyset$, the $m^{2}\times m^{2}$ matrix $F\bm{\mathsf{M}}^{\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}F^{*}{\text{diag}}([e^{-i\texttt{q}_{m}^{-1}(1)\cdot\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}},\ldots,e^{-i\texttt{q}_{m}^{-1}(m^{2})\cdot\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}}])$ satisfies

where $\Phi_{\bm{\mathsf{r}}_{i}\bm{\mathsf{r}}_{j}}:=\bm{\mathsf{r}}_{i}-\bm{\mathsf{r}}_{j}$,

Recall that the Fourier-Wigner transform of $\omega_{ij}$ is also called the ambiguity function of $\omega_{ij}$, which measures the spatial lag $\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}}$ and frequency shift $\Phi_{\bm{\mathsf{r}}_{i}\bm{\mathsf{r}}_{j}}$ between the two diffraction images when $\iota_{\bm{\mathsf{x}}_{i}}(D_{r}^{m})\cap\iota_{\bm{\mathsf{x}}_{j}}(D_{r}^{m})\neq\emptyset$. It is well-known that the absolute value of the ambiguity function gauges how difficult we can distinguish two objects, that is, how similar two objects are [36, p.33]. Thus $V_{\omega_{ij}}(\Delta_{\bm{\mathsf{x}}_{i}\bm{\mathsf{x}}_{j}},\Phi_{\bm{\mathsf{r}}_{i}\bm{\mathsf{r}}_{j}})$ can be viewed as a sort of affinity measuring the relationship between two illumination windows. Also, from (60) we know that the phase information of $F\text{diag}(\omega)$ gets involved in $V_{\omega_{ij}}$, in particular when $i=j$. Indeed, when we are working with the same patch, $\bm{\mathsf{M}}^{0}$ is a diagonal matrix with real entries $|\omega|^{2}$, so $F\bm{\mathsf{M}}^{0}F^{*}$ contains only the phase information of $F\text{diag}(\omega)$, which influences the phase estimation.

Another intuition behind the ptychgraphy is the following. If two illumination windows overlap, they have common information in the Fourier space up to some phase difference determined by the relative position of the illuminations, while this information is contaminated by the non-overlapping parts of the two illuminations.

2. Spectral relaxation and phase synchronization

Based on the above understanding regarding the $P_{\bm{\mathsf{FQ}}}$ and the amplitude information, in this section we propose two relaxations of the non-convex optimization problems discussed above to estimate the phase, which lead to a better initial value of the AP algorithm.

The first algorithm is directly motivated by (62) where we take the affinity information among vertices and phase relationship into account. We have the following observations.

the phase between vertices $(i,\bm{\mathsf{r}}_{i})$ and $(j,\bm{\mathsf{r}}_{j})$ are related by a non-unitary transform $\Omega((i,\bm{\mathsf{r}}_{i}),(j,\bm{\mathsf{r}}_{j}))$, which modulation indicating the affinity;

the larger the amplitude $\bm{\mathsf{a}}_{(i)}(\bm{\mathsf{r}}_{i})$ is, the more effort we should put in recovering the phase;

In addition, the phase ramping effect, denoted as

and a real $Km^{2}\times Km^{2}$ diagonal matrix $\bm{\mathsf{D}}$ so that

The synchronization property of GCL has been studied in . While noise is inevitable in real data, the robustness of GCL to different kinds of noises have been studied in the framework of block random matrix and reported in . In addition, under the manifold setup , it asymptotically converges to the heat kernel of the associated connection Laplacian, which top eigenvector-field is the most parallel vector field branded in the manifold structure. We refer the reader to the appendix of for a summary of the above results.

The second algorithm we propose has the same flavor, but we consider the amplitude information in a different way compared with (62). Indeed, the amplitude is taken into consideration as a truncation threshold leading to the following relaxation of (58) to estimate the phase. Based on the amplitude, we define a thresholding matrix

where $\epsilon_{a}\geq 0$ is the threshold chosen by the user, and evaluate the following functional

which is equivalent to finding the top eigenvector of the Hermitian matrix $T_{\bm{\mathsf{a}}}P_{\bm{\mathsf{FQ}}}T_{\bm{\mathsf{a}}}$. Our second proposed estimator of the phase to the ptychography problem is then the phase of the top eigenvector of $T_{\bm{\mathsf{a}}}P_{\bm{\mathsf{FQ}}}T_{\bm{\mathsf{a}}}$. We call this approach to the truncation phase synchronization (t-PS) algorithm. See Section 5 for its numerical performance. This optimization problem is essentially different from (56) due to the thresholding, and this difference plays an essential role in the optimization. Its theoretical property is beyond the scope of this paper and will be reported in another paper.

Numerical results

We begin with describing the two lens we use. The first one is a typical illumination probe in an experimental system. The illuminating beam is formed by a small lens, with a dark “beam-stop” to sort-out harmonic contaminations formed by diffractive Fresnel lenses, represented by a circular aperture in the Fourier domain. The lens is denoted as $\omega_{\text{s}}$ and is illustrated in the top row of Figure 6. The second is a band-limited random (BLR) lens, denoted as $\omega_{\text{BLR}}$ which we describe now. Note that a small lens can only “connect” Fourier frequencies that are close together, while a wide lens produces a small illumination and the illumination scheme can only connect frames that are near each other. The intuition behind the synchronization analysis of the ptychographic problem leads us to suggest a different lens that enables to connect pixels across the data space. Experimental observations confirm that diffuse probes , and wide apertures produce better results in ptychography. We design our second lens by setting the amplitude and a random phase of an annular aperture in the Fourier domain, then iteratively adjust the amplitude in real and Fourier domains to determine a lens with a circular focus and given amplitude. The motivation for the limited size of the focus is to reduce the requirements of the experimental detector response function (such as pixel size). Such lens can be fabricated using lithographic techniques . The second lens is described in the bottom row of Figure 6.

We begin with a small problem – an object of size $256\times 256$ pixels, that is $n=256$, shown in Figure 7, using the lens $\omega_{\text{s}}$. We collect $k=32\times 32$ frames, with $128\times 128$ pixels, that is $m=128$. The frames are distributed uniformly to cover the object: we start by setting the positions $\bm{\mathsf{x}}_{i}=(x_{i},y_{j})$ on a square grid lattice, with $x_{i}-x_{i+1}=\Delta x$ and $y_{i}-y_{i+1}=\Delta y$. In this first experiment, we take $\Delta x=\Delta y=8$. Then we shear odd rows, that is, $x_{i}$, by $\Delta x/2$ and perturb the position by a random perturbation randomly sampled uniformly from $[-1.5,+1.5]$ in both $x_{i}$ and $y_{i}$. Fractional pixel shifts are accounted by interpolation of the illumination matrix. We use the following algorithms, where PS is the abbreviation of phase synchronization.

start with random object: $\zeta^{(0)}=\bm{\mathsf{FQ}}(\text{random})$ ;

find the largest eigenvalue $v_{0}$ of the GCL matrix $\bm{\mathsf{D}}^{-1}\bm{\mathsf{S}}$;

$\psi_{\text{GCL-PS}}=(\bm{\mathsf{Q}}^{*}\bm{\mathsf{Q}})^{-1}\bm{\mathsf{Q}}^{*}\bm{\mathsf{F}}^{*}P_{\bm{\mathsf{a}}}v_{0}$.

find the largest eigenvalue $v_{0}$ of the phase synchronization matrix $T_{\bm{\mathsf{a}}}P_{\bm{\mathsf{FQ}}}T_{\bm{\mathsf{a}}}$, where $T_{\bm{\mathsf{a}}}=\text{diag}(\chi_{\bm{\mathsf{a}}>\epsilon_{a}})$;

$\psi_{\text{t-PS}}=(\bm{\mathsf{Q}}^{*}\bm{\mathsf{Q}})^{-1}\bm{\mathsf{Q}}^{*}\bm{\mathsf{F}}^{*}P_{\bm{\mathsf{a}}}v_{0}$.

find the largest eigenvalue $v_{0}$ of $\bm{\mathsf{D}}^{-1}\bm{\mathsf{S}}$;

find the largest eigenvalue $v_{0}$ of $T_{\bm{\mathsf{a}}}P_{\bm{\mathsf{FQ}}}T_{\bm{\mathsf{a}}}$;

The result of the first experiment is shown in Figure 7.

We repeat the same experiment with an image of a self-assembled cluster of $50$ nm colloidal gold nanoparticles obtained by Scanning Electron Microscopy. To produce a complex image, the gray-scale value are projected onto a circle in the complex plane. The size is $256\times 256$ pixels and we use the lens $\omega_{\text{s}}$. The result of the second experiment is shown in Figure 8.

We compare these two illumination functions, $\omega_{\text{s}}$ and $\omega_{\text{BLR}}$, with the same two objects with the same parameters as before. The results are shown in Figure 9 and Figure 10. Clearly t-PS produces a better start with the new illumination. In this example, such better start also leads to higher rate of convergence.

Yet next, we test the algorithm in a larger problem, an object of $512\times 512$ pixels, that is $n=512$, with the same lens size ($128\times 128$). We increase the field of view of the illumination scheme with increased spacing among frames $\Delta x=16$ and $\Delta y=16$. One of the issues of projection algorithms such as AP is that frames that are far apart communicate very weakly with each other, this leads to slower rate of convergence. This is an issue when we are limited by the number of iterations, due to high data rate and finite computational resources. In Figure 11 we show the result of $101$ iterations of AP with holes in the scarf, while t-PS gives a good initial start that leads to improved SNR. Notice that the hole in the scarf and other defects are produced by AP alone.

In our next numerical experiment, we introduce new algorithms that lead to over $80\times$ acceleration in the rate of convergence. First, we use the RAAR algorithm described below which is popular among the optical community (using RAAR in combination with a shrink-wrap algorithm to enforce sparsity) because it often leads to improved convergence rate. Second, we introduce a frame-wise synchronization technique to adjust the phase of every frame at every iteration based on existing frame-wide local information. Finally, we combine frame-wise synchronization with projected conjugate gradient (CG).

start with random object $\zeta^{(0)}=\bm{\mathsf{FQ}}(\text{random})$

find the largest eigenvalue $v_{0}$ of the kernel $T_{\bm{\mathsf{a}}}P_{\bm{\mathsf{FQ}}}T_{\bm{\mathsf{a}}}$

find the largest eigenvalue $v_{0}$ of the kernel $T_{\bm{\mathsf{a}}}P_{\bm{\mathsf{FQ}}}T_{\bm{\mathsf{a}}}$. Start

where $\bm{\mathsf{B}}$ is a $K\times K$ diagonal block matrix with its diagonal the $m^{2}\times 1$ row vector $\bm{1}^{T}$ that distributes the frame-wise phase to all the pixels;

t-PS: see above to initialize $\zeta^{(0)}$

repeat (2)-(3) until convergences or maximum iterations

The frame-wise synchronization, step (2), is motivated by the augmented approach . We estimate a phase factor for each frame based on the existing phase estimator of each frames, which leads to long-range phase synchronization across the image. Indeed, we consider

where $\bm{\mathsf{B}}$ is a $K\times K$ diagonal block matrix with its diagonal the $m^{2}\times 1$ row vector $\bm{1}^{T}$ that distributes the phase over the frame. We can re-write as:

The scaling factor $\left(\bm{\mathsf{Q}}^{\ast}\bm{\mathsf{Q}}\right)$ in $P_{\bm{\mathsf{FQ}}}$ can be weighted out by considering the pairwise relationship:

by swapping the diagonal matrix $T^{\ast}_{(i)}R^{\ast}Q_{(i)}$. We optimize the frame-wise phase vector $\bm{\xi}$ based on the existing estimator

We tested these algorithms, as well as the AP and t-PS+AP algorithms, on the same data setup in Figure 11, and the convergence results of different algorithms are shown in Figure 12 for comparison. Notice the change of scale in the last plot, where convergence is over 80$\times$ faster than the AP algorithm.

In our final test, we test the AP algorithm with noise. Noisy data is simulated using a proxy for Poisson statistics. We define $\bm{\sigma}$ a randomly distributed gaussian noise, and simulate noisy data and define the measurement error $\varepsilon_{\sigma}$:

Conclusions

In this paper, we demonstrate the the necessary and sufficient conditions of the local convergence of the alternating projection (AP) algorithm to the unique solution up to a global phase factor, and apply it to the ptychography imaging problem. To be more precise, we have conditions so that the user can check if the AP algorithm gives the inverse transform of the phase retrieval problem when the frame is generic. We also survey the intimate relationship between the AP algorithm and the notion of phase synchronization and propose two algorithm, GCL-PS and t-PS, to quickly construct an accurate initial guess for the AP algorithm for large scale diffraction data problems. In addition, by combining the RAAR algorithm or conjugate gradient method with the frame-wise synchronization, the convergence is over $40-80\times$ faster than the AP algorithm and is about $10\times$ faster than the RAAR algorithm.

There are several problems left unanswered in this paper. We mention at least the following four directions. First, in addition to the global convergence issue of the AP algorithm, how to design the best lens and illumination scheme so that we can obtain an accurate reconstruction for the real samples; given a detector, with a limited rate, dynamic range and response function, what is the best scheme to encode more information per detector channel. Second, the noise influence on the convergence behavior needs further investigation. Experimental uncertainties include not only photon-counting statistics but also perturbations of the lens , illumination scheme (positions), incoherent measurements, detector response and discretization, time dependent fluctuations, etc. Third, spectral methods such as the proposed algorithms in this paper (GCL-PS and t-PS) have the potential to be scaled up on high-performance computing architectures to handle the big imaging data in the coming new light source era . Last, although RAAR, synchro-RAAR and other iterative schemes perform well in practice, their convergence behavior needs to be further studied. Can we design better iterative methods based on our findings that exploit phase synchronization schemes more efficiently?

Acknowledgements

This work is partially supported by the Center for Applied Mathematics for Energy Research Applications (CAMERA), which is a partnership between Basic Energy Sciences (BES) and Advanced Scientific Computing Research (ASRC) at the U.S. Department of Energy (SM) and by AFOSR grant FA9550-09-1-0643 (HT). The authors would like to thank Professor Arthur Szlam, Dr. Jeffrey J. Donatelli and Dr. Wenjing Liao for their inputs to improve the paper. H.-T. Wu thanks Professor Ingrid Daubechies and Professor Albert Fannajing for the discussion. We acknowledge NVIDIA for providing us with a Tesla K40 GPU for our tests.

References