Frequency-Resolved Optical Gating Recovery via Smoothing Gradient

TL;DR

This paper introduces a new algorithm for reconstructing ultrashort laser pulses from FROG data, which is more effective with incomplete data and is simple, scalable, and theoretically sound.

Contribution

The paper proposes a novel smoothing gradient-based algorithm for FROG pulse recovery, improving performance with incomplete data and providing convergence guarantees.

Findings

Outperforms state-of-the-art with incomplete FROG traces

Computational cost comparable to existing methods

Converges to critical points near the true solution

Abstract

Frequency-resolved optical gating (FROG) is a popular technique for complete characterization of ultrashort laser pulses. The acquired data in FROG, called FROG trace, is the Fourier magnitude of the product of the unknown pulse with a time-shifted version of itself, for several different shifts. To estimate the pulse from the FROG trace, we propose an algorithm that minimizes a smoothed non-convex least-squares objective function. The method consists of two steps. First, we approximate the pulse by an iterative spectral algorithm. Then, the attained initialization is refined based upon a sequence of block stochastic gradient iterations. The algorithm is theoretically simple, numerically scalable, and easy-to-implement. Empirically, our approach outperforms the state-of-the-art when the FROG trace is incomplete, that is, when only few shifts are recorded. Simulations also suggest that the proposed algorithm exhibits similar computational cost compared to a state-of-the-art technique for both complete and incomplete data. In addition, we prove that in the vicinity of the true solution, the algorithm converges to a critical point. A Matlab implementation is publicly available at https://github.com/samuelpinilla/FROG.

Figures

Tables

Table 1. TABLE I : Comparison of iteration count and time cost

	Algorithms	Iterations	Time (s)
$L=1$	BSGA	60	1.451
	Ptych	36	1.325
$L=2$	BSGA	111	1.567
	Ptych	125	1.954
$L=4$	BSGA	265	1.772
	Ptych	300	2.013

Equations

$$\mathbf{Z}[p,k]$$

$$k$$

$$\text{dist}(\mathbf{x},\mathbf{w}):=\frac{{\Vert \sqrt{\mathbf{Z}}-\sqrt{\mathbf{W}}\Vert }_{\text{F}}}{{\Vert \sqrt{\mathbf{Z}}\Vert }_{\text{F}}},$$

$$\underset{\mathbf{z}\in {\mathbb{C}}^n}{\min }h(\mathbf{z},\mu )=\underset{\mathbf{z}\in {\mathbb{C}}^n}{\min }\frac{1}{NR}\sum _{k=0}^{N-1}\sum _{p=0}^{R-1}{\ell }_{k,p}(\mathbf{z},\mu ),$$

$${\ell }_{k,p}(\mathbf{z},\mu ):={\left[{\phi }_{\mu }\left(\left|\sum _{n=0}^{N-1}\mathbf{z}[n]\mathbf{z}[n+pL]e^{-2\pi ink/N}\right|\right)-\sqrt{\mathbf{Z}[p,k]}\right]}^2.$$

$${\phi }_{\mu }(w):=\sqrt{w^2+{\mu }^2},$$

$${\mathbf{f}}_k^H:=\left[{\omega }^{-0(k-1)},{\omega }^{-1(k-1)},\cdots ,{\omega }^{-(n-1)(k-1)}\right],$$

$$\frac{\partial h(\mathbf{z},\mu )}{\partial \overline{\mathbf{z}}[\ell ]}:=$$

$${\overline{q}}_{\ell ,p}:=$$

$${\mathbf{g}}_p:=$$

$$\frac{\partial h(\mathbf{z},\mu )}{\partial \overline{\mathbf{z}}}:={\left[\frac{\partial h(\mathbf{z},\mu )}{\partial \overline{\mathbf{z}}[0]},\cdots ,\frac{\partial h(\mathbf{z},\mu )}{\partial \overline{\mathbf{z}}[N-1]}\right]}^H.$$

$${\mathbf{x}}^{(t+1)}:={\mathbf{x}}^{(t)}-\alpha \frac{\partial h({\mathbf{x}}^{(t)},{\mu }^{(t)})}{\partial \overline{\mathbf{z}}},$$

$${\mathbf{d}}_{{\Gamma }_{(t)}}[\ell ]=\sum _{p,k\in {\Gamma }_{(t)}}\left({\mathbf{f}}_k^H{\mathbf{g}}_p^{(t)}-{\upsilon }_{k,p,t}\right){\overline{q}}_{\ell ,p}^{(t)}{e}^{2\pi i\ell k/N},$$

$$\underset{t\to \infty }{\lim }{\mu }^{(t)}=0,\text{ and }\underset{t\to \infty }{\lim }{\Vert \frac{\partial h({\mathbf{x}}^{(t)},{\mu }^{(t)})}{\partial \overline{\mathbf{z}}}\Vert }_2=0,$$

$$\mathbf{Y}[p,\ell ]=\frac{1}{N}\sum _{k=0}^{N-1}\mathbf{Z}[p,k]e^{-2\pi ik\ell /N}$$

$$=$$

$$\mathbf{Y}[p,\ell ]$$

$$=tr(\mathbf{X}{\mathbf{A}}_{p,\ell }),$$

$${\mathbf{y}}_{\ell }={\mathbf{G}}_{\ell }{\mathbf{x}}_{\ell },$$

$${\mathbf{G}}_{\ell }[p,n]:={\mathbf{x}}_{pL}[n]{\overline{\mathbf{x}}}_{pL}[n+\ell ].$$

$${\mathbf{x}}_{ini\_pty}[r]:=\mathbf{v}[r]\exp (i\mathbf{\theta }[r]),$$

$$\mathbf{v}[r]$$

$$:=\sum _{n=0}^{N-1}|\mathbf{x}[n]{|}^2|\mathbf{x}[n+rL]{|}^2.$$

$${\mathbf{x}}_{\ell }^{(0)}=\text{diag}({\mathbf{X}}_0^{(0)},\ell ),$$

$${\mathbf{X}}_0^{(0)}={\mathbf{x}}_{ini\_pty}{\mathbf{x}}_{ini\_pty}^H.$$

$$\text{diag}({\mathbf{X}}_0^{(t)},\ell )={\mathbf{x}}_{\ell }^{(t)}.$$

$${\mathbf{G}}_{\ell }^{(t)}[p,n]={\mathbf{x}}_{pL}^{(t)}[n]{\overline{\mathbf{x}}}_{pL}^{(t)}[n+\ell ],$$

$$\underset{{\mathbf{p}}_{\ell }\in {\mathbb{C}}^N}{\min }$$

$${\mathbf{x}}_{\ell }^{(t)}={\mathbf{B}}_{\ell ,t}^{-1}{\mathbf{e}}_{\ell ,t},$$

$${\mathbf{B}}_{\ell ,t}$$

$${\mathbf{e}}_{\ell ,t}$$

$$\text{diag}({\mathbf{X}}_0^{(t)},\ell )={\mathbf{B}}_{\ell ,t}^{-1}{\mathbf{e}}_{\ell ,t},\ell =0,\cdots ,N-1.$$

Full text

I Introduction

Phase retrieval (PR) is the inverse problem of recovering a signal from its Fourier magnitude [1]. PR arises in many fields in science and engineering, such as optics [2, 3], astronomical imaging [4], microscopy [5], and X-ray crystallography [6]. In this work, we focus on a popular technique for full characterization of ultrashort pulses called frequency-resolved optical gating (FROG) [7, 8]. The acquired data in FROG corresponds to the squared Fourier magnitude of the product of the unknown pulse with its delayed replica, for several different time shifts. The product of the signal with itself is usually performed using a second harmonic generation crystal. This measured data is called the FROG trace. In this paper, we focus on the inverse problem of recovering a pulse from its second-harmonic generation FROG trace.

Recent works have studied conditions under which a pulse can be uniquely identified, up to trivial ambiguities, from its FROG trace [9, 10]. In Section II we present and discuss these results. In particular, it has been shown that in theory not all the delay steps are needed to recover the pulse. The most commonly used algorithm to estimate a pulse from its FROG trace is the principal component generalized projections (PCGP), originally introduced in [11]. PCGP follows classical algorithms in PR based on alternating projections. Specifically, PCGP is initialized by a Gaussian pulse with random phases. It then builds an auxiliary matrix by rearranging the columns in the time-delay plane which is Fourier transformed. The pulse is next updated as the leading eigenvector of this matrix. The FROG trace of the updated pulse is then constructed and its magnitude is replaced by the acquired phaseless measurements. This procedure is repeated until convergence.

A recent algorithm, called Retrieved-Amplitude N-grid Algorithmic (RANA) [12], exploits the expected continuity of the signal to construct several initial estimates of the power spectrum of the pulse. These estimates are obtained from a set of smaller grids of the FROG trace, a strategy called multi-grid. The initial guesses are then refined using the generalized projections method in [13], in parallel, where the one that best fits the acquired measurements is chosen as the reconstructed pulse. The multi-grid initialization procedure allows RANA to converge faster than PCGP and leads to more accurate reconstructions. Importantly, both RANA and PCGP were developed to retrieve the pulse when all the delay steps of the FROG trace are acquired.

In [14], the authors suggest an alternative recovery strategy, inspired by blind ptychography [1]. This method starts with the integrated measured FROG trace over the frequency dimension as an initial guess. Then, the initialization is refined using a stochastic descent strategy, which involves a single time shift of the FROG trace per iteration to update the estimated pulse. This paper claims to attain better estimation of the pulse compared to PCGP. However, the non-smooth cost function may increase the amount of required measurements (sample complexity) to recover the pulse, which affects the success rate as will be shown in the numerical results. This limitation appears since the non-smoothness of the objective may lead to unbounded gradients [15, 16].

In this paper we propose a block stochastic gradient algorithm (BSGA) for FROG recovery that minimizes a smoothed amplitude-based least-squares empirical loss. Amplitude-based objectives have shown improved results in standard PR [17]. BSGA is initialized by a spectral method that requires only few iterations. This procedure can be seen as a modification of the strategy proposed in [18] that approximates the signal $\mathbf{x}$ from the FROG trace as the leading eigenvector of a carefully designed matrix. The two stages (initialization and gradient iterations) differ from prior contributions by the new initialization technique and the inclusion of a smoothing function, following [15]. Specifically, in contrast to [12], the proposed initialization aims to estimate the pulse, rather than its power spectrum. Numerical results show that our initialization returns a more accurate estimation of the pulse compared to the starting point of the ptychography (Ptych) method in [14]. Furthermore, BSGA shows improvements in recovering both the magnitude and the phase of the pulse from incomplete data ($L>1$) compared to Ptych. As aforementioned, alternative methods such as PCGP and RANA are designed for $L=1$, and therefore do not work well in this regime. In Theorem 1 we provide partial theoretical justification for the success of the algorithm by showing that in the vicinity of the true solution, BSGA converges to a critical point.

The paper is organized as follows. We begin in Section II by introducing necessary background on FROG. Section III presents the proposed initialization technique and introduces an iterative procedure to refine the solution by minimizing a smooth least-squares objective. Section V presents numerical results and compares our approach with competitive algorithms. Finally, Section VI concludes the paper.

We denote by $\mathbb{R}_{+}:=\{w\in\mathbb{R}:w\geq 0\}$ and $\mathbb{R}_{++}:=\{w\in\mathbb{R}:w>0\}$ the sets of positive and strictly positive real numbers, respectively. The conjugate and the conjugate transpose of the vector $\mathbf{w}\in\mathbb{C}^{N}$ are denoted as $\overline{\mathbf{w}}\in\mathbb{C}^{N}$ and $\mathbf{w}^{H}\in\mathbb{C}^{N}$, respectively. The $n$th entry of a vector $\mathbf{w}$, which is assumed to be periodic, is written as $\mathbf{w}[n]$. We denote by $\tilde{\mathbf{w}}$ and $\hat{\mathbf{w}}$ the Fourier transform of a vector and its conjugate reflected version (that is, $\hat{\mathbf{w}}[n]:=\overline{\mathbf{w}}[-n]$). The notation $\text{diag}(\mathbf{W},\ell)$ refers to a column vector with entries $\mathbf{W}[j,(j+\ell)\mod N]$ for $j=0,\cdots,N-1$. For vectors, $\|\mathbf{w}\|_{p}$ is the $\ell_{p}$ norm. Additionally, we use $\odot,$ and $*$ for the Hadamard (point-wise) product, and convolution, respectively. Finally, $\mathbb{E}[\cdot]$ represents the expected value.

II Problem formulation

FROG is probably the most commonly-used approach for full characterization of ultrashort optical pulses due to its simplicity and good experimental performance [7]. Experimentally, a FROG apparatus produces a two-dimensional intensity diagram, also known as FROG trace, of an input pulse by interacting the pulse with delayed versions of itself in a nonlinear-optical medium, usually using a second harmonic generation (SHG) crystal [7]. Mathematically, the FROG trace of a signal $\mathbf{x}\in\mathbb{C}^{N}$ is defined as

[TABLE]

with $R=\lceil N/L\rceil$ where $L<N$ and $i:=\sqrt{-1}$. This work assumes that the signal $\mathbf{x}$ is periodic, that is, $\mathbf{x}[n]=\mathbf{x}[n+lN]$ for any $l\in\mathbb{Z}$.

The FROG trace defined in (1) can be considered as a map $\mathbb{C}^{N}\rightarrow\mathbb{R}_{+}^{\lceil N/L\rceil}$ that has three types of symmetry, usually called trivial ambiguities in the PR literature. These ambiguities are summarized in Proposition 1, using the following definition of a bandlimited signal.

Definition 1**.**

We say that $\mathbf{x}\in\mathbb{C}^{N}$ is a $\hskip 100.00015ptB-bandlimited$ signal if its Fourier transform $\tilde{\mathbf{x}}\in\mathbb{C}^{N}$ contains $N-B$ consecutive zeros. That is, there exists $k$ such that $\tilde{\mathbf{x}}[k]=\cdots=\tilde{\mathbf{x}}[N+k+B-1]=0$.

Proposition 1**.**

([10]) Let $\mathbf{x}\in\mathbb{C}^{N}$ be the underlying signal and let $\tilde{\mathbf{x}}\in\mathbb{C}^{N}$ be its Fourier transform. Let $\mathbf{Z}[p,k]$ be the FROG trace of $\mathbf{x}$ defined as in (1) for some fixed $L$. Then, the following signals have the same FROG trace as $\mathbf{x}$:

the rotated signal $\mathbf{x}e^{i\phi}$ for some $\phi\in\mathbb{R}$; 2. 2.

the translated signal $\mathbf{x}^{\ell}$ obeying $\mathbf{x}^{\ell}[n]=\mathbf{x}[n-\ell]$ for some $\ell\in\mathbb{Z}$ (equivalently, a signal with Fourier transform $\tilde{\mathbf{x}}^{\ell}$ obeying $\tilde{\mathbf{x}}^{\ell}[k]=\tilde{\mathbf{x}}[k]e^{-2\pi i\ell k/N}$ for some $\ell\in\mathbb{Z}$); 3. 3.

the reflected signal $\hat{\mathbf{x}}$ obeying $\hat{\mathbf{x}}[n]:=\overline{\mathbf{x}}[-n]$.

If $\mathbf{x}$ is a $B$-bandlimited signal for some $B\leq N/2$, then the translation ambiguity is continuous. Namely, any signal with a Fourier transform such that $\tilde{\mathbf{x}}^{\psi}[k]:=\tilde{\mathbf{x}}[k]e^{i\psi k}$ for some $\psi\in\mathbb{R}$, has the same FROG trace as $\mathbf{x}$.

Our goal is to estimate the signal $\mathbf{x}$, up to trivial ambiguities, from the FROG trace $\mathbf{Z}$. The work [10] established that the pulse $\mathbf{x}$ can be uniquely identified (up to trivial ambiguities) from the FROG trace under rather mild conditions as summarized in the following proposition.

Proposition 2**.**

([10]) Let $\mathbf{x}\in\mathbb{C}^{N}$ be a $B$-bandlimited signal as in Definition 1 for some $B\leq N/2$. If $N/L\geq 4$, then almost all signals are determined uniquely from their FROG trace $\mathbf{Z}[p,k]$, up to trivial ambiguities, from $m\geq 3B$ measurements. If in addition we have access to the signal’s power spectrum and $N/L\geq 3$, then $m\geq 2B$ measurements suffice.

Proposition 2 has been recently extended to the case of blind ptychography, or blind FROG, in which the goal is to estimate two signals simultaneously [19]. Evidently, Proposition 2 allows choices of $L>1$ meaning that not all the delay steps are needed to recover the pulse, and therefore a method that works in this regime as well is desired.

To take the ambiguities into account, we measure the relative error between the true signal $\mathbf{x}$ and any $\mathbf{w}\in\mathbb{C}^{N}$ as

[TABLE]

where $\mathbf{Z}$ is the FROG trace of $\mathbf{x}$ according to (1), $\sqrt{\cdot}$ is the point-wise square root, $\mathbf{W}$ is the FROG trace of $\mathbf{w}$, and $\lVert\cdot\rVert_{\text{F}}$ denotes the Frobenius norm. Note that if $\text{dist}(\mathbf{x},\mathbf{w})=0$, and the uniqueness conditions of Proposition 2 are met, then for almost all signals $\mathbf{w}$ is equal to $\mathbf{x}$ up to trivial ambiguities.

In recent years, many papers have examined the problem of recovering a signal from phaseless quadratic random measurements. A popular approach is to minimize the intensity least-squares objective; see for instance [20]. Recent works have shown that minimizing the amplitude least-squares objective leads to better reconstruction under noisy scenarios [17, 21, 22]. However, the latter cost function is non-smooth and thus may lead to a biased descent direction [15]. To overcome the non-smoothness of the objective function, we follow the smoothing strategy proposed in [15].

The smooth objective to recover the underlying pulse considered in this work is

[TABLE]

where

[TABLE]

The function $\varphi_{\mu}:\mathbb{R}\rightarrow\mathbb{R}_{++}$ in (4) is defined as

[TABLE]

with $\mu\in\mathbb{R}_{++}$ (a tunable parameter). Notice that if $\mu=0$, then (4) reduces to the non-smooth formulation. In [21], the authors addressed the non-smoothness by introducing truncation parameters into the gradient step in order to eliminate the errors in the estimated descent direction. However, this procedure can modify the search direction and increase the sample complexity of the phase retrieval problem [15].

In this work we propose a block stochastic gradient algorithm (BSGA) to solve (3), that is initialized by a spectral procedure which requires only a few iterations. Section III explains in detail the proposed algorithm.

III Reconstruction Algorithm

In order to solve the optimization problem in (3), we develop a gradient-based algorithm, called BSGA. The algorithm is initialized by the outcome of a spectral method approximating the signal $\mathbf{x}$ which will be explained in Section IV.

To refine the initial estimate we use the Wirtinger derivatives as introduced in [23]. Let us define the vector $\mathbf{f}_{k}^{H}$ as

[TABLE]

with $\omega=e^{\frac{2\pi i}{n}}$ the $n$th root of unity. Then, the Wirtinger derivative of $h(\mathbf{z},\mu)$ in (3) with respect to $\overline{\mathbf{z}}[\ell]$ is given by

[TABLE]

where $\upsilon_{k,p}:=\sqrt{\mathbf{Z}[p,k]}\frac{\mathbf{f}_{k}^{H}\mathbf{g}_{p}}{\varphi_{\mu}\left(\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}\right\rvert\right)}$, and

[TABLE]

The gradient of $h(\mathbf{z},\mu)$ is then

[TABLE]

Using (8), we define a standard gradient algorithm, taking the form of

[TABLE]

where $\alpha$ is the step size.

To alleviate the memory requirements and computational complexity required for large $N$, we suggest a block stochastic gradient descent strategy. Instead of calculating (6), we choose only a random subset of the sum for each iteration $t$, that is,

[TABLE]

where the set $\Gamma_{(t)}$ is chosen uniformly and independently at random at each iteration $t$ from subsets of $\hskip 50.00008pt\{1,\cdots,N\}\times\{1\cdots,R\}$ with cardinality $Q$. Specifically, the gradient in (8) is uniformly sampled using a minibatch of data, in this case of size $Q$ for each step update, such that in expectation is (6) [24, page 130].

As mentioned in Section III, choosing $\mu>0$ prevents bias in the update direction. Since the function $h$ is smooth, we are able to construct a descent rule for $\mu$ (Line 13 of Algorithm 1) in order to guarantee convergence to a first-order optimal point, that is, a point with zero gradient, in the vicinity of the solution.

Theorem 1**.**

Let $\mathbf{x}$ be $B$-bandlimited for some $B\leq N/2$, satisfying $\text{dist}(\mathbf{x},\mathbf{x}^{(t)})\leq\rho$ for some sufficiently small constant $\rho>0$. Suppose that $L=1$ and $\Gamma_{(t)}$ is sampled uniformly at random from all subsets of $\{1,\cdots,N\}\times\{1\cdots,R\}$ with cardinality $Q$, independently for each iteration. Then for almost all signals, Algorithm 1 with step size $\alpha\in(0,\frac{2}{U}]$ satisfies

[TABLE]

for some constant $U>0$ depending on $\rho$.

Proof.

See Appendix A. ∎

IV Initialization Algorithm

In this section we devise a method to initialize the gradient iterations. This strategy approximates the signal $\mathbf{x}$ from the FROG trace as the leading eigenvector of a carefully designed matrix. We divide the exposition of the initialization procedure into two cases, $L=1$ and $L>1$, explained in Sections IV-A and IV-B, respectively.

IV-A Initialization for $L=1$

Instead of directly dealing with the FROG trace in (1), we consider the acquired data in a transformed domain by taking its 1D DFT with respect to the frequency variable (normalized by $1/N$). Our measurement model is then

[TABLE]

where $p,\ell=0,\cdots,N-1$. Observe that for fixed $p$, $\mathbf{Y}[p,\ell]$ is the autocorrelation of $\mathbf{x}\odot\mathbf{x}_{pL}$, where $\mathbf{x}_{pL}[n]=\mathbf{x}[n+pL]$.

Let $\mathbf{D}_{pL}\in\mathbb{C}^{N\times N}$ be a diagonal matrix composed of the entries of $\mathbf{x}_{pL}$, and let $\mathbf{C}_{\ell}$ be a circulant matrix that shifts the entries of a vector by $\ell$ locations, namely, $(\mathbf{C}_{\ell}\mathbf{x})[n]=\mathbf{x}[n+\ell]$. Then, the matrix $\mathbf{X}:=\mathbf{x}\mathbf{x}^{H}$ is linearly mapped to $\mathbf{Y}[p,\ell]$ as follows:

[TABLE]

where $\mathbf{A}_{p,\ell}=\mathbf{C}_{-\ell}\mathbf{D}_{pL}\overline{\mathbf{D}}_{pL+\ell}$, and $tr(\cdot)$ denotes the trace function. Observe that $\mathbf{C}_{\ell}^{T}=\mathbf{C}_{-\ell}$. Thus, we have that

[TABLE]

for a fixed $\ell\in\{0,\cdots,N-1\}$, where $\mathbf{y}_{\ell}[n]=\mathbf{Y}[n,\ell]$ and $\mathbf{x}_{\ell}=\text{diag}(\mathbf{X},\ell)$. The $(p,n)$th entry of the matrix $\mathbf{G}_{\ell}\in\mathbb{C}^{\lceil\frac{N}{L}\rceil\times N}$ is given by

[TABLE]

Since $L=1$, it follows from (15) that $\mathbf{G}_{\ell}$ is a circulant matrix. Therefore, $\mathbf{G}_{\ell}$ is invertible if and only if the DFT of its first column, in this case $\mathbf{x}\odot(\mathbf{C}_{\ell}\overline{\mathbf{x}})$, is non-vanishing.

Using (14), we propose a method to estimate the signal $\mathbf{x}$ from measurements (1) using an alternating scheme: fixing $\mathbf{G}_{\ell}$, solving for $\mathbf{x}_{\ell}$, updating $\mathbf{G}_{\ell}$ and so forth. The new methodology proposed in [18] cannot be directly employed since here the matrices $\mathbf{G}_{\ell}$ are also unknown. Thus, our approach estimates the matrices $\mathbf{G}_{\ell}$ together with $\mathbf{x}_{\ell}$.

We start the alternating scheme by the initialization suggested in [14]

[TABLE]

where $\boldsymbol{\theta}[r]\in[0,2\pi)$ is chosen uniformly at random for all $r\in\{0,\cdots,N-1\}$. The $r$th entry of $\mathbf{v}$ corresponds to the summation of the measured FROG trace over the frequency axis:

[TABLE]

Once the vector $\mathbf{x}_{ini\_pty}$ is constructed, the vectors $\mathbf{x}_{\ell}^{(t)}$ at $t=0$ can be built as

[TABLE]

where

[TABLE]

Then, from (18) we proceed with an alternating procedure between estimating the matrix $\mathbf{G}_{\ell}$, and updating the vector $\mathbf{x}_{\ell}$ as follows.

• •

*Update rule for $\mathbf{G}_{\ell}$: *In order to update $\mathbf{G}_{\ell}$, we update the matrix $\mathbf{X}_{0}^{(t)}$ as

[TABLE]

Observe that if $\mathbf{x}_{\ell}^{(t)}$ is close to $\mathbf{x}_{\ell}$ for all $\ell$, then $\mathbf{X}_{0}^{(t)}$ is close to $\mathbf{x}\mathbf{x}^{H}$. Letting $\mathbf{w}^{(t)}$ be the leading (unit-norm) eigenvector of the matrix $\mathbf{X}_{0}^{(t)}$ constructed in (20), from (15) each matrix $\mathbf{G}^{(t)}_{\ell}$ at iteration $t$ is given by

[TABLE]

where $\mathbf{x}^{(t)}_{pL}[n]=\mathbf{w}^{(t)}[n+pL]$.

• •

*Optimization with respect to $\mathbf{x}_{\ell}$: *Fixing $\mathbf{G}_{\ell}^{(t-1)}$, one can estimate $\mathbf{x}^{(t)}_{\ell}$ at iteration $t$ by solving the linear least-squares (LS) problem

[TABLE]

The relationship between the vectors $\mathbf{x}^{(t)}_{\ell}$ is ignored at this stage. If $\mathbf{G}_{\ell}^{(t-1)}$ is invertible, then the solution to this problem is given by $(\mathbf{G}_{\ell}^{(t-1)})^{-1}\mathbf{y}_{\ell}$. Since $\mathbf{G}_{\ell}^{(t-1)}$ is a circulant matrix, it is invertible if and only if the DFT of $\mathbf{x}^{(t-1)}\odot(\mathbf{C}_{\ell}\overline{\mathbf{x}}^{(t-1)})$ is non-vanishing. This condition cannot be ensured in general. Thus, we propose a surrogate proximal optimization problem to estimate $\mathbf{x}_{\ell}^{(t)}$ by

[TABLE]

where $\lambda>0$ is a regularization parameter. In practice $\lambda$ is a tunable parameter [25]. In particular, for this work the value of $\lambda$ was determined using a cross-validation strategy such that each simulation uses the value that results in the smallest relative error according to (2). The surrogate optimization problem in (23) is strongly convex [25], and admits the following closed form solution

[TABLE]

where

[TABLE]

with $\mathbf{I}\in\mathbb{R}^{N\times N}$ the identity matrix. Clearly $\mathbf{B}_{\ell,t}$ in (25) is always invertible. The update step for each $\mathbf{x}_{\ell}^{(t)}$ is computed in Line 9 of Algorithm 2.

Finally, in order to estimate $\mathbf{x}$, the (unit-norm) principal eigenvector of $\mathbf{X}_{0}^{(T)}$ is normalized by

[TABLE]

where $\mathcal{S}:=\left\{n:\left(\mathbf{B}_{0,T}^{-1}\mathbf{e}_{0,T}\right)[n]>0\right\}$. Observe that (26) results from the fact that $\sum_{n=0}^{N-1}\text{diag}(\mathbf{X},0)[n]=\lVert\mathbf{x}\rVert^{2}_{2}$.

After a few iterations of this two-step procedure, the output is used to initialize the gradient algorithm described in Section III. This alternating scheme is summarized in Algorithm 2.

IV-B FROG initialization step for $L>1$

Until now we focused on the case $L=1$. If $L>1$, then the linear system in (14) is underdetermined and $\mathbf{y}_{\ell}$ can be viewed as a subsampled version of (12) by a factor $L$. Therefore, in order to increase the number of equations when $L>1$, we up-sample $\mathbf{y}_{\ell}$ by a factor $L$. Specifically, we follow the proposed scheme in [18] that expands the measurement vector $\mathbf{y}_{\ell}$ by low-pass interpolation. Once the measurements are upsampled, we proceed as for $L=1$. This initialization, for $L>1$, is summarized in Algorithm 3. From Line 3 to Line 5 the low-pass interpolation by a factor $L$ is computed, and then in Line 6, Algorithm 2 generates the initial estimation of the underlying signal.

V Numerical Results

This section evaluates the numerical performance of BSGA and compares the results with the stochastic gradient algorithm Ptych proposed in [14]. We used the following parameters for Algorithm 1: $\gamma_{1}=0.1$, $\gamma=0.1$, $\alpha=0.6$, $\mu_{0}=65$, and $\epsilon=1\times 10^{-10}$. The number of indices that are chosen uniformly at random is fixed as $Q=N$. A cubic interpolation was used in Algorithm 3 (see Line 4), and the regularization parameter was fixed to $\lambda=0.5$.

Five tests were conducted to evaluate the performance of the proposed method under noisy and noiseless scenarios at different values of signal-to-noise-ratio (SNR), defined as SNR$=10\log_{10}(\lVert\mathbf{Z}\rVert^{2}_{\text{F}}/\lVert\boldsymbol{\sigma}\rVert^{2}_{\text{2}})$, where $\boldsymbol{\sigma}$ is the variance of the noise. First, we examine the empirical success rate of BSGA for different values of $L$. The second experiment assesses the performance of the initialization technique and its impact on the reconstruction quality. Third, we show several examples of reconstructed pulses attained with BSGA and Ptych under noisy and noiseless scenarios, when the complete FROG trace is used. The fourth experiment investigates the performance of the proposed method and Ptych in reconstructing the pulses when $L>1$ and the FROG trace is corrupted by noise. The last test compares the computational complexity between the reconstruction methods in terms of their running time to reach a given relative error.

The signals used in the simulations were constructed as follows. For all tests, we built a set of $\left\lceil\frac{N-1}{2}\right\rceil$-bandlimited pulses that conform to a Gaussian power spectrum centered at 800 nm. Specifically, each pulse ($N=128$ grid points) is produced via the Fourier transform of a complex vector with a Gaussian-shaped amplitude with a cutoff frequency of $150$ femtoseconds*-1* (fsec*-1*). Next, we multiply the obtained power spectrum by a uniformly distributed random phase. In the experiments we used the inverse Fourier of this signal as the underlying pulse.

All simulations were implemented in Matlab R2019a on an Intel Core i7 3.41Ghz CPU with 32 GB RAM. The code for BSGA is publicly available at https://github.com/samuelpinilla/FROG. The code of Ptych was downloaded from the authors’ website111https://oren.net.technion.ac.il/homepage/.

V-A Empirical Probability of Success

This section numerically evaluates the success rate of BSGA. To this end, BSGA and Ptych are initialized at $\mathbf{x}^{(0)}=\mathbf{x}+\delta\zeta$, where $\delta$ is a fixed constant and $\zeta$ takes values on $\{-1,1\}$ with equal probability, while $L$ ranges from 1 to 6. A trial is declared successful when the returned estimate attains a relative error as in (2) that is smaller than $10^{-6}$. We numerically determine the empirical success rate among 100 trials. Fig. 1 summarizes these results, and shows that BSGA performs better than Ptych, since it is able to retrieve the signal for larger values of $L$.

V-B Relative Error of the Initialization Procedure

This section examines the impact of the designed initialization described in Algorithms 2 and 3, under noisy and noiseless scenarios. We compare the relative error between the starting vector in (16), and the returned solution $\mathbf{x}^{(0)}$ of the proposed initialization procedure. The number of iterations to attain the vector $\mathbf{x}^{(0)}$ using the designed initialization was fixed as $T=2$, and we numerically determine the relative error averaged over 100 trials. These numerical results are summarized in Fig. 2, and indicate that the proposed initialization algorithm outperforms $\mathbf{x}_{ini\_pty}$.

In order to illustrate the effect of the initial guesses, we ran Algorithm 1 initialized by $\mathbf{x}_{ini\_pty}$ and $\mathbf{x}^{(0)}$ with $L=4$. Fig. 3 shows the attained reconstructions. Notice that the proposed reconstruction algorithm fails in estimating the input pulse when it was initialized by $\mathbf{x}_{ini\_pty}$.

We numerically determine the performance of the proposed initialization at different SNR levels, with $L$ ranging from 1 to 8. Specifically, we added white noise to the FROG measurements at different SNR levels: SNR = 8dB, 12dB, 16dB and 20dB. Fig. 4 displays the relative error attained by the proposed initialization for different SNR and $L$ values.

From Fig. 4 it can be seen that the returned initialization at levels of SNR $\leq 16$dB is, approximately, independent of the value of $L$ when $L\leq 6$. Combining these numerical results with Fig. 1, we conclude that BSGA is able to better estimate the underlying pulse (up to trivial ambiguities) if $L\leq 4$ for both noiseless and noisy scenarios compared to Ptych.

Finally, we numerically determine the empirical success rate of BSGA with increasing $L$, in the absence of noise, when Algorithm 1 is initialized with $\mathbf{x}_{ini\_pty}$, a random vector and $\mathbf{x}^{(0)}$. A trial is declared successful when the returned estimate attains a relative error as in (2) that is smaller than $1\times 10^{-6}$. The results are summarized in Fig. 5, where the number of iterations that BSGA requires to reach the given relative error for $L=1$ is also presented. The success rate and the number of iterations are averaged over 100 pulses. The reported results show the effectiveness of Algorithm 1 when it is initialized by $\mathbf{x}^{(0)}$ for $L>1$.

V-C Pulse Reconstruction Examples for $L=1$

In this section we show the performance of BSGA in recovering two pulses under noiseless and noisy scenarios for $L=1$. The results are presented in Figs. 6, and 7, respectively, where the attained relative errors by BSGA and Ptych are included. For the second scenario, the FROG trace is corrupted by Gaussian noise with SNR = 20dB.

From the results in Fig. 6 it can be observed that both methods, BSGA and Ptych, are able to estimate the pulses and provide similar results for the noiseless case.

On the other hand, in Fig. 7, the attained reconstructions, for the noisy scenario, indicate that BSGA is able to better estimate the pulse compared to Ptych. This advantage is obtained because of the effectiveness of the proposed smoothing update step and initialization strategy from complete data as reported in Fig. 1, and Figs. 2, 4, respectively.

V-D Pulse Reconstruction Examples for $L>1$

Next, we examine the recovery performance of BSGA from noisy incomplete data by adding Gaussian noise with SNR = 20dB, for $L\in\{2,4,8\}$. Figs. 8 and 9 illustrate the attained reconstructions for BSGA and Ptych; their attained relative errors are also reported in Fig. 9. These figures suggest that BSGA better estimates the pulse and its FROG trace compared to Ptych over a range of values of $L$. This advantage is obtained because of the effectiveness of the proposed smoothing update step and initialization strategy from incomplete data as reported in Fig. 1, and Figs. 2, 4, respectively.

V-E Computational Complexity

Simulations were conducted to compare the speed of convergence of the algorithms in the absence of noise, for $L=1,2$ and $4$. Table I reports the number of iterations and running time required by BSGA and Ptych to achieve a relative error of $1\times 10^{-6}$, averaged over 100 pulses. The experiment shows that BSGA is similar in time and number of iterations compared to Ptych for a range of values of $L$.

VI Conclusion

This paper presented a new method called BSGA, to fully characterize a pulse from its FROG trace. Our algorithm consists of two steps: a spectral initialization followed by successive refinements based upon a sequence of block stochastic gradient iterations. The principles of the proposed method were adopted from algorithms developed for the problems of recovering a signal from random quadratic measurements, and from the short-time Fourier phaseless data.

Numerical experiments were conducted to evaluate the performance of BSGA in comparison to the ptychography technique developed in [14]. The results show improvements in recovering the pulse for both magnitude and phase, from noisy incomplete data. Additionally, the numerical results suggest the effectiveness of the proposed initialization under both noiseless and noisy scenarios with incomplete data. Future work should include implementing BSGA on real data to further validate its performance. Another interesting research direction is to examine similar strategies for blind FROG in which two signals are estimated simultaneously [7].

Appendix A Proof Theorem 1

Let us define the search set as

[TABLE]

for some small constant $\rho>0$. Recall that $\mathbf{z}$ is a $B$-bandlimited signal if there exists $k$ such that $\tilde{\mathbf{z}}[k]=\cdots=\tilde{\mathbf{z}}[N+k+B-1]=0$, where $\tilde{\mathbf{z}}$ is the Fourier transform of $\mathbf{z}$. The bandlimit condition guarantees that we have unique solution, according to Proposition 2.

In order to prove Theorem 1, the function $h(\mathbf{z},\mu)$ in (3) needs to satisfy the four requirements stated in the following lemma, which are used in the analysis of convergence for stochastic gradient methods [26].

Lemma 1**.**

The function $h(\mathbf{z},\mu)$ in (3) and its Wirtinger derivative in (8) satisfy the following properties.

The cost function $h(\mathbf{z},\mu)$ in (3) is bounded below. 2. 2.

The set $\mathcal{J}$ as defined in (27) is closed and bounded. 3. 3.

There exists a constant $U>0$, such that

[TABLE]

holds for all $\mathbf{z}_{1},\mathbf{z}_{2}\in\mathcal{J}$. 4. 4.

For all $\mathbf{z}\in\mathcal{J}$

[TABLE]

for some $\zeta>0$, where $\mathbf{d}_{\Gamma_{(t)}}$ is as in Line 9 of Algorithm 1.

Proof.

See Appendix B. ∎

To prove Theorem 1, denote the set $\hskip 80.00012pt\mathcal{K}_{1}:=\{t|\mu^{(t+1)}=\gamma_{1}\mu^{(t)}\}$ with $\gamma_{1}\in(0,1)$, which is a tunable parameter [27]. If the set $\mathcal{K}_{1}$ is finite, then according to Lines 13-16 in Algorithm 1 there exists an integer $\grave{t}$, such that, for all $t>\grave{t}$

[TABLE]

with $\gamma\in(0,1)$. Taking $\grave{\mu}=\mu^{(\grave{t})}$, the optimization problem (3) reduces to

[TABLE]

Now, considering the properties stated in Lemma 1, from [26, Theorem 2.1] we get

[TABLE]

It can be readily seen that (32) contradicts the assumption $\left\lVert\mathbf{d}_{\Gamma_{(t)}}\right\rVert_{2}\geq\gamma\mu^{(\grave{t})}$, for all $t>\grave{t}$. This shows that $\mathcal{K}_{1}$ must be infinite and $\displaystyle\lim_{t\rightarrow\infty}\mu^{(t)}=0$.

Given that $\mathcal{K}_{1}$ is infinite, we deduce that

[TABLE]

where the second line follows from the Jensen inequality. Therefore, from (33) the result of Theorem 1 holds.

Appendix B Proof of Lemma 1

The proof of Lemma 1 is obtained by individually proving the following four requirements.

1. Following from the definition of $h(\mathbf{z},\mu)$ in (3) it is clear that $h(\mathbf{z},\mu)\geq 0$ and thus bounded below.

2. This holds by definition.

3. From (6) it follows that the $\ell$-th entry of $\frac{\partial h(\mathbf{z},\mu)}{\partial\overline{\mathbf{z}}}$ is given by

[TABLE]

where $\upsilon_{k,p}=\sqrt{\mathbf{Z}[p,k]}\frac{\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})}{\varphi_{\mu}\left(\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})\right\rvert\right)}$, and

[TABLE]

Let $\mathbf{D}_{p}(\mathbf{z})$ be a diagonal matrix composed of the entries of $\overline{\mathbf{z}}_{pL}[n]=\overline{\mathbf{z}}[n+pL]$. Using (5), the term $\overline{q}_{\ell,p}e^{2\pi i\ell k/N}$ can be rewritten as

[TABLE]

Thus,

[TABLE]

where

[TABLE]

and

[TABLE]

To prove 3) we establish that any $f_{k,p}(\mathbf{z})$ and $g_{k,p}(\mathbf{z})$ satisfy

[TABLE]

and

[TABLE]

for all $\mathbf{z}_{1},\mathbf{z}_{2}\in\mathcal{J}$ and some constants $r_{k,p},s_{k,p}>0$. In fact, once we prove (39), it can be performed a similar analysis for $g_{k,p}(\mathbf{z})$, and thus the result of this third part holds.

From the definition of $f_{k,p}(\mathbf{z})$, for any $\mathbf{z}_{1},\mathbf{z}_{2}\in\mathcal{J}$ we have that

[TABLE]

considering that $\mathbf{D}_{p}(\mathbf{z}_{1})$ and $\mathbf{D}_{p}(\mathbf{z}_{2})$ are diagonal matrices, and $\lVert\mathbf{f}_{k}\rVert_{2}=\sqrt{N}$. Observe that from (38) and (41) it can be obtained that

[TABLE]

where the second inequality comes from the fact that $\varphi_{\mu}\left(\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z}_{1})\right\rvert\right)\geq\mu$. The term $p_{1}$ in (42) can be upper bounded as

[TABLE]

Recall that $\mathcal{J}$ is a closed bounded set, and thus compact. Since $\varphi_{\mu}(\cdot)$ is a continuous function, there exists a constant $M_{\varphi_{\mu}}$ such that $\varphi_{\mu}\left(\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})\right\rvert\right)\leq M_{\varphi_{\mu}}$ for all $\mathbf{z}\in\mathcal{J}$. Also, from Lemma 2 in [15] we have that $\varphi_{\mu}(\cdot)$ is a 1-Lipschitz function. Combining this with (43) we get

[TABLE]

and thus

[TABLE]

where (45) results from applying the triangular inequality. Putting together (42) and (45) we obtain that

[TABLE]

Observe that the upper bound in (46) directly depends on a term of the form $\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})$ for some $\mathbf{z}\in\mathcal{J}$, which might be zero. However, Lemma 2 proves that $\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})\right\rvert>0$ or equivalently $\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})\not=0$, for almost all $\mathbf{z}\in\mathcal{J}$.

Lemma 2**.**

Let $\mathbf{z}\in\mathcal{J}$ where $\mathcal{J}$ as defined in (27). Then, for almost all $\mathbf{z}\in\mathcal{J}$ the following holds

[TABLE]

for all $k,p\in\{0,\cdots,N-1\}$, with $\mathbf{g}_{p}(\mathbf{z})$ as in (B).

Proof.

We prove this lemma by contradiction. Suppose that $\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})\right\rvert=0$. Then, from (1) we have that

[TABLE]

Observe that (48) is a quartic polynomial equation with respect to the entries of $\mathbf{z}$. However, for almost all signals $\mathbf{z}\in\mathcal{J}$ the left hand side of (48) will not be equal to zero which leads to a contradiction [10]. ∎

Then, proceeding to bound the term $\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z})\right\rvert$, notice that from (1) we have that

[TABLE]

in which the second inequality arises from $\lVert\mathbf{z}\rVert_{2}\leq\sqrt{N}\lVert\mathbf{z}\rVert_{\infty}$ and $\lVert\mathbf{z}\rVert_{1}\leq\sqrt{N}\lVert\mathbf{z}\rVert_{2}$. Combining (46) and (49) we get

[TABLE]

Now, we have to analyze the term $\left\lvert\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z}_{1})-\mathbf{f}_{k}^{H}\mathbf{g}_{p}(\mathbf{z}_{2})\right\rvert$ in (50). Specifically, from (1) it can be obtained that

[TABLE]

where the second inequality results from $\lVert\mathbf{z}\rVert_{2}\leq\sqrt{N}\lVert\mathbf{z}\rVert_{\infty}$ and $\lVert\mathbf{z}\rVert_{1}\leq\sqrt{N}\lVert\mathbf{z}\rVert_{2}$. Combining (50) and (51) we obtain that

[TABLE]

where $r_{k,p}$ is given by

[TABLE]

Since the set $\mathcal{J}$ is bounded, then $\lVert\mathbf{z}\rVert_{2}<\infty$ for all $\mathbf{z}\in\mathcal{J}$. Therefore, $0<r_{k,p}<\infty$, and from (52) the result holds.

4. We proceed to prove (29). Observe that

[TABLE]

in which the inequality comes from the fact that $\left\lVert\mathbf{w}_{1}+\mathbf{w}_{2}\right\rVert^{2}_{2}\leq 2\left(\left\lVert\mathbf{w}_{1}\right\rVert^{2}_{2}+\left\lVert\mathbf{w}_{2}\right\rVert^{2}_{2}\right)$ for any $\hskip 50.00008pt\mathbf{w}_{1},\mathbf{w}_{2}\in\mathbb{C}^{N}$. Combining (28) and (54) we have that

[TABLE]

for some $U>0$. Recall that $\Gamma_{(t)}$ is sampled uniformly at random from all subsets of $\{1,\cdots,N\}\times\{1\cdots,R\}$ with cardinality $Q$. From the definition of $\mathbf{d}_{\Gamma_{(t)}}$ in Line 9 of Algorithm 1, it can be concluded that

[TABLE]

using the fact that $\left\lVert\mathbf{w}_{1}+\mathbf{w}_{2}\right\rVert^{2}_{2}\leq 2\left(\left\lVert\mathbf{w}_{1}\right\rVert^{2}_{2}+\left\lVert\mathbf{w}_{2}\right\rVert^{2}_{2}\right)$ for any $\mathbf{w}_{1},\mathbf{w}_{2}\in\mathbb{C}^{N}$. Furthermore, since $f_{k,p}(\mathbf{z})$ and $g_{k,p}(\mathbf{z})$ satisfy (39) and (40), respectively, we conclude that

[TABLE]

for some constants $r_{k,p},s_{k,p}>0$. Thus, combining (55) and (57) we obtain that

[TABLE]

where $\zeta$ is defined as

[TABLE]

Notice $\zeta<\infty$ because the set $\mathcal{J}$ is bounded. Thus, from (58) the result holds.