# Energy dependent forward B $\rightarrow$ $J/\psi$ measurements in p+p
collisions at PHENIX

P-25, MS H846, Physics Divison, LANL, Los Alamos, NM, 87545, USA
E-mail

###### Abstract:

The heavy flavor studies at RHIC help improve the knowledge of the
bottom/charm quark production and can test Quantum Chromodynamics (QCD).
Compared to the LHC/Tevatron, the RHIC heavy flavor production originates from
different partonic sub-processes and has a complementary kinematic coverage.
The PHENIX forward rapidity silicon vertex detector (FVTX) provides precise
determination of the event vertex, tracking and the Distance of Closest
Approach (DCA) of charged tracks. This detector allows direct access to the
$B$ meson production via measurements of non-prompt $J/\psi\to\mu^{+}+\mu^{-}$
within $1.2<|y|<2.2$ rapidity in $p$+$p$ collisions at $\sqrt{s}=$ 510 and 200
GeV. Comparison among PHENIX measurements of the $B\rightarrow J/\psi$
fraction with integrated $J/\psi$ $p_{T}$ up to 5 GeV$/c$ and higher energy
results at the Tevatron and the LHC presents a smooth center of mass energy
dependence from 0.2 to 13 TeV in $p$+$p$ ($p$+$\bar{p}$) collisions. The Next-
To-Leading order Perturbative QCD (NLO pQCD) calculations are in reasonable
agreement with the extracted total $b\bar{b}$ cross section based on the
$B\rightarrow J/\psi$ fraction measurements at PHENIX.

##  1 Introduction

Heavy quarks are produced in the early stage of high energy hadronic
collisions due to their high mass ($m_{\rm{c,b}}\gg\Lambda_{\rm{QCD}}$). Since
the heavy flavor quarks do not vanish or change into other flavors during the
hard scattering and the sequential fragmentation processes, they can be
treated as hard probes to study the initial state partons. At RHIC energies,
the bottom production is dominated by the gluon-gluon fusion (pair creation)
process in contrast to the Tevatron and the LHC measurements which are
dominated by the flavor excitation process [1]. In additional to this, the
kinematic range of the nucleon (nuclear) gluon distribution function (PDF)
accessed by the RHIC heavy flavor measurements is different from what has been
measured at the LHC. Therefore, measuring the heavy flavor production at RHIC
provides a unique test of Quantum Chromodynamics (QCD).

Yields of non-prompt $J/\psi$ from $B$ meson decay are proportional to the
total $B$ meson production. Extracting $B$ meson production via this channel
provides a better signal to background ratio compared to the bottom semi-
leptonic decay measurements as backgrounds are suppressed by the $J/\psi$
identification. Decayed particles with a non-zero displaced vertex can be
identified with the help of the Forward Silicon Vertex Detector (FVTX) at
PHENIX, which can precisely measure the Distance of Closet Approach (DCA) of
final state particles. Using this technique, $B$ mesons in the
forward/backward rapidities have been studied at PHENIX in $p$+$p$ collisions.
The $B\rightarrow J/\psi$ fractions have been measured and the $b\bar{b}$
cross sections are studied in 510 and 200 GeV $p$+$p$ collisions. The NLO pQCD
calculations are in good agreements with these results.

Figure 1:  The $B\to J/\psi$ fraction fit to muon $\textrm{DCA}_{\rm R}$ of
$J/\psi$ sample in the (a) $1.2<y<2.2$ and (b) $-2.2<y<-1.2$ regions. The red
solid curve stands for the total fit which includes the prompt $J/\psi$ (solid
blue), the $B$-meson $\rightarrow J/\psi$ (green filled region), the
combinatorial background (magenta dashed curve), the $c\bar{c}+b\bar{b}$
background (brown long-dashed curve) and the detector mismatching background
(purple short-dashed curve). Figures from [2].

##  2 The fraction of $J/\psi$ from $B$ meson decay in $\sqrt{s}$ = 510 and
200 GeV $p$+$p$ collisions

The fraction of $J/\psi$ from $B$ meson decay at forward/backward rapidities
($1.2<|y|<2.2$) has been first studied at PHENIX in the 510 GeV $p$+$p$
collisions during the 2012 RHIC run. Identification of $J/\psi$ from $B$ meson
decay is based on measurements of the DCA radial projection
($\textrm{DCA}_{\rm R}$) of muons from $J/\psi$ decay as the FVTX has a better
spatial resolution along the radial direction compared to its azimuthal angle
component. The $\textrm{DCA}_{\rm R}$ of muons from prompt $J/\psi$ decay
follows a symmetric distribution. Because of the decay kinematics, muons from
$B$ meson decayed $J/\psi$ produce an asymmetric tail in their
$\textrm{DCA}_{\rm R}$ distribution.

Figure 2:  The $B\to J/\psi$ fraction with integrated $J/\psi$ $p_{T}$ and
$1.2<|y|<2.2$ rapidity in $p$+$p$ collisions at $\sqrt{s}=$ 510 (left) GeV and
$\sqrt{s}=$ 200 (right) GeV. In the left panel, the red shaded band stands for
the FONLL+CEM calculated $p_{T}$ dependent $B\to J/\psi$ fraction within
$1.2<|y|<2.2$ in 500 GeV $p$+$p$ collisions and the blue shaded curve
represents the $p_{T}$ integrated ($p_{T}<$ 5 GeV$/c$)$B\to J/\psi$ fraction
within the same kinematic region. The yellow shaded curve in the right panel
stands for the FONLL+CEM predicted rapidity dependent $B\to J/\psi$ fractions
in 200 GeV $p$+$p$ collisions.

This analysis requires muon tracks pass through quality cuts and have good
matching between the FVTX and the Muon Tracker (MuTr). $J/\psi$ candidates are
selected from unlike-sign dimuon pairs within the mass region of 2.7-3.5
GeV/$c^{2}$. After verifications of consistent $\textrm{DCA}_{\rm R}$
resolutions in both data and the full simulation of
PYTHIA+GEANT+reconstruction with realistic vertex and dead maps etc, the
$\textrm{DCA}_{\rm R}$ distribution shapes of muons from prompt $J/\psi$ and
$B$ meson decayed $J/\psi$ are determined in the full simulation. The detector
mis-alignments have been corrected before and after the data production. Any
remaining $\textrm{DCA}_{\rm R}$ offset is determined by identified prompt
hadrons in data. Besides the prompt $J/\psi$ and $B\rightarrow J/\psi$
components, the dimuon sample in data also contains the combinatorial
background, the detector mis-matching background and the heavy flavor
($c\bar{c}+b\bar{b}$) continuum background. The combinatorial background which
represents the mis-identified muon and/or hadron pairs is determined by
unlike-sign track pairs in normalized mixed event samples. Since a nearly 1
$m$ long hadron absorber is located between the FVTX and the MuTr, there is a
certain probability of matching an incorrect FVTX track to a MuTr track after
the multiple scattering within the absorber. The FVTX-MuTr mis-matching
background is determined in mixed events with topologically rotated FVTX. The
heavy flavor continuum background is determined using a full simulation.
Details of the analysis procedure have been discussed in [2]. A maximum log-
likelihood fit of the muon $\textrm{DCA}_{\rm R}$ distribution is applied to
data to simultaneously determine the raw yields of prompt $J/\psi$ and non-
prompt $J/\psi$ from $B$ meson decay. Figure 1 shows the fit results within
$-2.2<y<-1.2$ and $1.2<y<2.2$ rapidity regions in 510 GeV $p$+$p$ collisions.

Figure 3:  Comparison of PHENIX $B$ $\rightarrow$ $J/\psi$ fraction measured
in 510 GeV [2] and 200 GeV $p$+$p$ collisions [3] with the global data from
CDF, ALICE, CMS and LHCb experiments as a function of center of mass energy
integrated in the $J/\psi$ $0<p_{T}<5$ GeV/$c$ interval. The uncertainty is
statistical and systematic combined.

Figure 4:  Comparison of the $b\bar{b}$ cross section extracted from the
PHENIX $B$ $\rightarrow$ $J/\psi$ fraction measured in 510 GeV [2] and 200 GeV
[3] $p$+$p$ collisions with previous measurements at Fermilab, RHIC, LHC and
NLO pQCD calculations. The uncertainty of the NLO pQCD calculations is shown
in green dashed line, and the systematic uncertainty of experimental results
is shown in filled boxes.

The detector acceptance$\times$efficiency corrected fraction of $J/\psi$ from
$B$ meson decay with integrated $J/\psi$ $p_{T}$ ($p_{T}<5$ GeV/c) and
$1.2<|y|<2.2$ rapidity coverage measured in 510 GeV $p$+$p$ collisions is
shown in the left panel of Figure 2. Reasonable agreement between this result
and the Fixed-Order Next-To-Leading Logarithm and Color-Evaporation-Model
(FONLL+CEM) [4, 5, 6] calculations has been achieved. Following the same
analysis procedure, the $B\rightarrow J/\psi$ fraction is measured in 200 GeV
$p$+$p$ collisions with the 2015 RHIC data. The right panel of Figure 2 shows
the measured results in 200 GeV $p$+$p$ collisions with integrated $J/\psi$
$p_{T}$ in the forward and backward rapidity regions, which are also in good
agreement with the derived rapidity dependent $B\rightarrow J/\psi$ fractions
based on the FONLL [7] and the CEM [8] calculations.

##  3 Energy dependent $B\to J/\psi$ fraction and extrapolated $b\bar{b}$
cross section in $p$+$p$ collisions

In order to understand the energy dependent $B$ hadron production, the PHENIX
measured $B\rightarrow J/\psi$ fractions in 510 and 200 GeV $p$+$p$ collisions
are compared with other Tevatron and LHC measurements [9, 10, 11, 12, 13] with
integrated $J/\psi$ $p_{T}$ up to 5 GeV/c at different rapidities and higher
center of mass energies. As shown in Figure 3, the center of mass energy
dependent $B\rightarrow J/\psi$ fraction results indicate a smooth transition
of $B$ production from the RHIC energies which are dominated by the gluon-
gluon fusion process to the LHC energies which are dominated by the gluon
splitting process. To further confirm this, the total $b\bar{b}$ cross
sections in 200 and 510 GeV $p$+$p$ collisions are extracted from the forward
$B\to J/\psi$ fraction results at PHENIX with the measured inclusive $J/\psi$
cross section [14], the branching ratio of $B$ meson decay to $J/\psi$ and the
phase space scaling factor calculated by the FONLL. The preliminary results of
the extracted $b\bar{b}$ cross sections in $p$+$p$ collisions, as shown in
Figure 4, are in reasonable agreement with the NLO pQCD calculations [15]. The
PHENIX extracted $b\bar{b}$ cross sections and other global measurements
follow the NLO pQCD predicted energy dependence of the bottom production.

##  4 Summary and Outlook

First measurements of $J/\psi$ from $B$ meson decay in forward/backward
rapidities with integrated $p_{T}$ starting from zero are achieved in both 510
and 200 GeV $p$+$p$ collisions at PHENIX. A smooth transition is found for the
$B\rightarrow J/\psi$ fractions with integrated $J/\psi$ $p_{T}$ up to 5
GeV$/c$ measured in $p$+$p$ or $p$+$\bar{p}$ collisions from 0.2 to 13 TeV.
The NLO pQCD calculated $b\bar{b}$ cross sections are consistent with the
extracted values based on the $B\rightarrow J/\psi$ fraction measurements at
PHENIX in 200 and 510 GeV $p$+$p$ collisions. Global data on the bottom cross
section follow the energy dependence predicted by the NLO pQCD calculations.

Large data sets in $p$+$p$ and $p$+Au collisions collected by PHENIX provide
opportunities to further study the open heavy flavor production in
forward/backward rapidities. This allows us to study the $p_{T}$ dependent
$B\rightarrow J/\psi$ yields and explore the cold nuclear matter effects on
the bottom production. The method used for the $B$ meson decay to $J/\psi$
fraction analysis through analyzing the muon $DCA_{R}$ can be extended to the
study of $B$ and $D$ meson semi-leptonic decays to muons in the
forward/backward rapidities to study the charm and bottom production in the
low $p_{T}$ region.

## References

  * [1] Norrbin E and Sjostrand T 2000 Eur. Phys. J. C 17 137 (Preprint arXiv: hep-ph/0005110)
  * [2] Aidala C et al. (PHENIX Collaboration) 2017 Phys. Rev. D 95 092002 (Preprint arXiv:1701.01342)
  * [3] Aidala C et al. 2017 (Preprint arXiv:1702.01085)
  * [4] Bedjidian M et al. 2004 Hard probes in heavy ion collisions at the LHC: Heavy flavor physics (Preprint arXiv: hep-ph/0311048)
  * [5] Cacciari M, Nason P and Vogt R 2005 Phys. Rev. Lett. 95(12) 122001
  * [6] Nelson R E, Vogt R and Frawley A D 2013 Phys. Rev. C 87 014908 (Preprint arXiv: 1210.4610)
  * [7] M Cacciari M G and Nason P 1998 JHEP 9805 007 (Preprint arXiv: hep-ph/9803400)
  * [8] Frawley A D, Ullrich T and Vogt R 2008 Phys. Rept. 462 125–175 (Preprint arXiv: 0806.1013)
  * [9] Acosta D et al. 2005 Phys. Rev. D 71 032001
  * [10] Abelev B et al. 2012 JHEP 11 065
  * [11] Khachatryan V et al. 2011 Eur. Phys. J. C 71 1575
  * [12] Aaij R et al. (LHCb Collaboration) 2013 JHEP 06 064 (Preprint arXiv: 1304.6977)
  * [13] Aaij R et al. 2015 JHEP 10 172
  * [14] Adare A et al. (PHENIX Collaboration) 2012 Phys. Rev. D 85 092004 (Preprint arXiv: 1105.1966)
  * [15] Vogt R 2008 The European Physical Journal Special Topics 155 213–222


# Phase retrieval by power iterations

Stefano Marchesini  Advanced Light Source, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720

###### Abstract

I show that the power iteration method applied to the phase retrieval problem
converges under special conditions. One is given the relative phases between
small non-overlapping groups of pixels of a recorded intensity pattern, but no
information on the phase between the groups of pixels. Numerical tests show
that the inverse block iteration recovers the solution in 1 iteration.

††preprint: LBNL-5939E

##  I Introduction

Given a set of intensity measurements, $a^{2}\in\mathbb{R}^{M}$, an unknown
object $\psi\in\mathbb{C}^{N}$ represented by a complex $n\times n$ image
($N=n^{2}$), a known “illumination matrix” or support matrix ${\bf Q}$
($\mathbb{C}^{M\times N}$ matrix), a known propagation operator (typically
one, or a stack of 2D FFT operators) $\bm{F}$ of dimension $M\times M$ and set
of frames ${\bf z}\in\mathbb{C}^{M}$, which are related by:

${\bf z}=\bm{F}{\bf Q}\psi,\qquad|{\bf z}|=a.$

Our goal is to find $\psi$ or the intermediate variable ${\bf z}$, given
$\bm{F}$, ${\bf Q}$ and $a$. To do so, we need to find a phase $\phi$ such
that ${\bf z}=a\phi$ is in the range of $\bm{F}{\bf Q}$.

We can eliminate $\psi$ by using the operator $P_{\bf Q}$ to project a vector
${\bf z}$ onto the range of $\bm{F}{\bf Q}$:

$\displaystyle P_{\bf Q}$ $\displaystyle=$ $\displaystyle\bm{F}{\bf Q}({\bf
Q}^{\ast}{\bf Q})^{-1}{\bf Q}^{\ast}\bm{F}^{\ast}.$ (1)

Which ensures that the unknown vector $\psi$ can be obtained from the frame
${\bf z}$ by $\psi=({\bf Q}^{\ast}{\bf Q})^{-1}{\bf Q}^{\ast}\bm{F}^{\ast}{\bf
z}$.

A popular approach is to find a vector $z$ such that:

$\displaystyle\left\|[I-P_{Q}]{\bf z}\right\|^{2}$ $\displaystyle=$
$\displaystyle 0,$ (2) $\displaystyle\left\|[I-P_{a}]{\bf z}\right\|^{2}$
$\displaystyle=$ $\displaystyle\||{\bf z}|-a\|^{2}=0,$ (3)

are satisfied simultaneously, and where the Fourier magnitude projection
$P_{a}$ when applied to a vector ${\bf z}$, yields:

$\displaystyle P_{{\bf a}}{\bf z}=\bm{a}\frac{{\bf z}}{|{\bf
z}|},\quad\bm{a}=\mbox{Diag}\left(a\right)\ \ $ (4)

where division are intended as element-wise operations Marchesini  _et al._
(2012); Marchesini (2007).

##  II Phase optimization

Here we want to minimize Eq. 2 w.r.t. a phase vector $\phi$
($\phi_{i}^{\ast}\phi_{i}=1,\forall i$). That is, we want to find:

$\displaystyle\arg\min_{\phi}$ $\displaystyle\left\|[I-P_{\bf
Q}]\mathrm{Diag}(a)\phi\right\|^{2},$ $\displaystyle\arg\min_{\phi}$
$\displaystyle\phi^{\ast}\bm{a}\left[I-P_{\bf Q}\right]\bm{a}\phi,$ (5)

I discuss three approaches that relax the phase modulus condition
($\phi_{i}^{\ast}\phi_{i}=1\quad\forall i$) to synchronize the relative
phases.

### Power iteration

By changing variable ${\bf z}=\bm{a}\phi$, we write:

$\displaystyle\arg\min_{\bf z}\quad{\bf z}^{\ast}\left(I-P_{\bf Q}\right){\bf
z}$ (6)

By relaxing $(a_{i}^{2}|\phi_{i}|^{2}=a_{i}^{2},\forall i)$ and using $\|{\bf
z}\|=\|\bm{a}\phi\|^{2}=\|\bm{a}\|^{2}$, we can re-write Eq. (6) as finding
the eigenvector with largest eigenvalueSinger (2011). Since $\|{\bf
z}\|=\|\bm{a}\phi\|=\|a\|$ is constant, we rewrite Eq, (5) as:

$\displaystyle\arg\max_{\bf z}\quad{\bf z}^{\ast}P_{\bf Q}{\bf z},\quad$ (7)

we apply one step of power iteration:

$\displaystyle\nu^{\ell+1}=P_{\bf Q}{\bf z}$ (8)

We then form a projection on the unit torus to ensure that
$|\nu^{\ell+1}_{i}|=1$ (or $|{\bf z}_{i}|=a_{i}$) by element-wise
normalization:

$\displaystyle z^{(\ell+1)}=\bm{a}\frac{\nu}{|\nu|}$ $\displaystyle=$
$\displaystyle P_{a}P_{\bf Q}z^{(\ell)}$

Here we have obtained the classical alternating projection method, which is
known to stagnate with classical CDI but to converge (slowly) in ptychographic
imaging.

### Greedy phase optimization

Since the diagonal term ${\bf z}^{\ast}\mbox{Diag}\left(P_{\bf Q}\right){\bf
z}$ is also independent on the choice of $\phi$ (for $|\phi_{i}|=1$), one can
remove it when computing the power iteration:

$\displaystyle\nu^{(\ell+1)}=[P_{Q}-\mbox{Diag}\left(P_{Q}\right)]{\bf
z}^{(\ell)}$ (9)

After we apply the projection of $\phi^{\ell+1}$ to the unit torus, we obtain
the following update Waldspurger  _et al._ (2012):

$z^{(\ell+1)}=P_{a}\left(P_{\bf
Q}-\mbox{Diag}\left(P_{Q}\right)\right)z^{(\ell)}$

In classical CDI, $P_{Q_{ii}}={\|{\bf Q}\|^{2}\over\|\bm{1}\|^{2}}$ is simply
the sum of the support volume (or area) normalized by the oversampled volume,
in ptychographic imaging $\mbox{Diag}\left(P_{\bf Q}\right)$ is the ratio of
intensities $P_{Q_{ii}}={\|Q_{i}\|^{2}\over{\|{\bf Q}\|^{2}}}$ for every pixel
of a frame $i$ generate by a submatrix $Q_{i}$. At the first iteration, using
data generated from the object in Fig. 6 with a random phase as a starting
guess, Eq. (9) appears to out-perform Eq. (6), however the two methods
converge to similar local minimum within ten iterations. The relaxations in
Eqs. (6,9) are similar. By removing diagonal components we change the
relaxation. In Eq. (6) we have $\|\bm{a}^{2}\phi\|$ constant, in Eq. (9)
$\|[1-\mbox{Diag}\left(P_{\bf Q}\right)]\bm{a}^{2}\phi\|$ is constant. However
$\mbox{Diag}\left(P_{Q}\right)$ is often constant and the two relaxations are
equivalent, giving more weight to high intensity values.

### Inverse iterationMarchesini  _et al._ (2012).

If we solve the minimization problem (Eq. (5) with a different relaxation,
setting $\|\phi\|=k$ to a constant, we re-write the problem as

$\displaystyle\arg\max\phi^{\ast}H^{-1}\phi,\quad H={\bf a}[I-P]{\bf a}$ (10)

and apply the power iteration:

$\displaystyle H\nu^{(\ell+1)}=\phi^{\ell}$ (11)

This method is commonly referred to as inverse iteration and it is used to
find the smallest eigenvector of a matrix. We note however that any $\nu$
written in the following way:

$\nu=\left({\bf a}^{-1}P_{\bf Q}{\bf a}\right){\bf x}$ (12)

is an eigenvector with 0 eigenvalue of ${\bf a}(I-P_{\bf Q}){\bf a}$,
therefore the inverse iteration method cannot be applied directly.

When $H$ is singular, then instead of power iteration we may want to find the
smallest modification of the phase that is in the null space of $H$, which we
can write it as a re-weighted LSQ problem of the form:

$\displaystyle\arg\min_{\bf z}\left\|\frac{1}{{\bf a}}({\bf z}-{\bf
z}^{\ell})\right\|,\quad\text{s.t. ${\bf z}=P_{\bf Q}{\bf z}$}$ (13)

which provides a search direction toward the solution that differs from
standard projection algorithms. Another approach is to include additional
restrictions on $\nu$ before applying the inverse iteration as described in
the following.

### Inverse block iteration

In Marchesini  _et al._ (2012) it was observed that computing the exact
solution to Eq. (11) after “binning”, or fixing the relative phase between
groups of pixels, improved convergence rate in large scale ptychographic
imaging. The use eigensolvers for the interferometric case was also suggested
in Alexeev  _et al._ (2012), for the connection Laplacian of a graph.

Let us introduce a binning matrix $\bm{T}^{\ast}$ composed of a series of
masks $T_{i}$ that integrate over a region of dimension $M/k$ of the data (in
Fourier domain). For example, we can partition our data in 3, creating a tall
matrix of dimension $M\times 3$:

$\bm{T}=\left(\begin{array}[]{c}{\bm{1}_{M/3},\bm{0}_{M/3},\bm{0}_{M/3}}\\\
{\bm{0}_{M/3},\bm{1}_{M/3},\bm{0}_{M/3}}\\\
{\bm{0}_{M/3},\bm{0}_{M/3},\bm{1}_{M/3}}\end{array}\right)$

where
$T_{1}=(\bm{1}^{\ast}_{M/3},\bm{0}^{\ast}_{M/3},\bm{0}^{\ast}_{M/3})^{\ast}$
is a vector of length $M$.

We restrict our search of the solution to Eq. (11) by restricting $\nu$ to be:

$\displaystyle\nu=\mbox{Diag}\left(\phi^{\ell}\right)T\omega,$ (14)

If we multiply from the left by
$T^{\ast}\mbox{Diag}\left(\phi^{\ell}\right)^{\ast}$ in Eq. (11) we obtain the
inverse iteration step with initial 0-phase vector as first guess:

$\displaystyle\hat{H}^{(\ell)}\omega^{(\ell+1)}=\lambda_{1}\bm{1},\quad$ (15)

Where

$\hat{H}^{(\ell)}=T^{\ast}\mbox{Diag}\left(z^{(\ell)}\right)^{\ast}[I-P_{\bf
Q}]\mbox{Diag}\left(z^{(\ell)}\right)T,$

$\lambda_{1}$ is a scalar multiplicative factor, and $\bm{1}$ is a vector of
appropriate length ($=3$ in this example). By computing $\omega^{(\ell+1)}$
from Eq. (15), and $\nu^{\ell+1}$ from Eq. (14), and projecting on the unit
torus we obtain the update $\phi^{(\ell+1)}$:

$\displaystyle\phi^{(\ell+1)}=\frac{\nu^{(\ell+1)}}{|\nu^{(\ell+1)}|}=\mbox{Diag}\left(\phi^{(\ell)}\right)\frac{T\omega^{(\ell+1)}}{|T\omega^{(\ell+1)}|}.$
(16)

In the following section we’ll show an example of the inverse iteration
method.

##  III Numerical example

Here $\psi$ consists of the cameraman image of $32\times 32$ pixels, embedded
in a matrix of $64\times 64$ pixels (Fig. 6). The “illumination matrix” is the
support of the object, ${\bf Q}=\mbox{Diag}\left(S\right)$. The support is 1
inside the $32\times 32$ box containing the image, and 0 otherwise. $({\bf
Q}^{\ast}{\bf Q})^{-1}$ is replaced by the pseudoinverse ${\bf Q}={\bf
Q}^{\ast}$. The Fourier transform of $\psi$ was perturbed by $32\times 32$
randomly distributed phases (Fig. 6), each multiplying a bin of $2\times 2$
pixels (Fig. 6). Upon perturbation, the image in real space (Fig. 6) cannot be
distinguished. Many iterations of Eq. 6 or Eq. 9 cannot converge (Fig. 6
showing Eq. (6 ) updates), while 1 iteration of Eqs. (15,16) converges to the
solution (Fig. 6).

##  IV Conclusions

I have shown that power iteration methods can recover phase perturbations
under special circumstances. If one is given the relative phases between a
small group of pixels (binned) and a random perturbation of the phase between
all the groups of pixels (the bins), then the inverse block iteration can
recover the solution in 1 iteration. In Marchesini  _et al._ (2012) it was
observed that the inverse block iteration improved convergence rate in large
scale ptychographic imaging. The inverse block iteration was also shown to
recover perturbations in the experimental geometry such as position errors and
intensity fluctuations. More work is needed to determine the optimal
combination of Eqs. (6,9,13,15,16), and the properties of $T$, in large scale
phase retrieval problems.

I acknowledge usefull discussions with Jeff Donatelli of UC Berkeley. This
work was stimulated by the Phase Retrieval workshop at the Erwin Schroedinger
International Institute for Mathematical Physics (ESI) organized by Karlheinz
Gröchenig and Thomas Strohmer. This work was supported by the Laboratory
Directed Research and Development Program of Lawrence Berkeley National
Laboratory under the U.S. Department of Energy contract number DE-
AC02-05CH11231.

## Disclaimers

This document was prepared as an account of work sponsored by the United
States Government. While this document is believed to contain correct
information, neither the United States Government nor any agency thereof, nor
the Regents of the University of California, nor any of their employees, makes
any warranty, express or implied, or assumes any legal responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product,
or process disclosed, or represents that its use would not infringe privately
owned rights. Reference herein to any specific commercial product, process, or
service by its trade name, trademark, manufacturer, or otherwise, does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof, or the Regents of the
University of California. The views and opinions of authors expressed herein
do not necessarily state or reflect those of the United States Government or
any agency thereof or the Regents of the University of California.

## References

  * Marchesini  _et al._ (2012) S. Marchesini, A. Schirotzek, F. Maia, and C. Yang, (2012), arXiv:1209.4924 [physics.optics] .
  * Marchesini (2007) S. Marchesini, Rev Sci Instrum 78, 011301 (2007), arXiv:physics/0603201 .
  * Singer (2011) A. Singer, Applied and Computational Harmonic Analysis 30, 20 (2011), arXiv:0905.3174 .
  * Waldspurger  _et al._ (2012) I. Waldspurger, A. d’Aspremont, and S. Mallat, ArXiv (2012), arXiv:1206.0102 [math.OC] .
  * Alexeev  _et al._ (2012) B. Alexeev, A. S. Bandeira, M. Fickus, and D. G. Mixon, (2012), arXiv:1210.7752 [cs.IT] .

Figure 1:  Object $\psi$ ($64\times 64$) used to simulate diffraction data

Figure 2:  Each column of the matrix $T$ extracts an area of $2\times 2$
pixels out of an image of ($64\times 64$) pixels.

Figure 3:  Random perturbation: ($64\times 64$) phases generated by $32\times
32$ random phases each spread over a bin of ($2\times 2$) pixels.

Figure 4:  Image in real space ($\left|\bm{F}^{\ast}{\bf z}\right|$) after
random phase perturbation (Fig. 6)

Figure 5:  Image in real space after random phase perturbation (Fig. 6), using
Eq. 6 updates ( $\left|\bm{F}^{\ast}(P_{\bm{a}}P_{S})^{1000}{\bf z}\right|$).

Figure 6:  Image in real space after one step of Eqs. (15,16) update.