Solving nonlinear integral equations for laser pulse retrieval with Newton's method

We present an algorithm based on numerical techniques that have become standard for solving nonlinear integral equations: Newton's method, homotopy continuation, the multilevel method, and random projection to solve the inversion problem that appears when retrieving the electric field of an ultrashort laser pulse from a two-dimensional intensity map measured with frequency-resolved optical gating (FROG), dispersion-scan, or amplitude-swing experiments. Here we apply the solver to FROG and specify the necessary modifications for similar integrals. Unlike other approaches we transform the integral and work in time domain where the integral can be discretized as an overdetermined polynomial system and evaluated through list autocorrelations. The solution curve is partially continues and partially stochastic, consisting of small linked path segments and enables the computation of optimal solutions in the presents of noise. Interestingly, this is an alternative method to find real solutions of polynomial systems, which are notoriously difficult to find. We show how to implement adaptive Tikhonov-type regularization to smooth the solution when dealing with noisy data, and we compare the results for noisy test data with a least-squares solver and propose the L-curve method to fine-tune the regularization parameter.

Figure 1

Illustrating example: $E\left(t\right)$ discretized with a piecewise constant function on $N=4$ intervals (top left). Then the associated products $F_{{\tau }_i}\left(t\right)$ and $F_{{\tau }_i}\left(t\right)F_{{\tau }_i}(t-\sigma )$ (bottom) are piecewise constant as well on small parallelograms such that the nonlinear integral $J\left[E\right](\tau ,\sigma )$ (top right) can be computed via list autocorrelations along all grid segments (red lines). We consider $J$ only in quadrand $(+,+), \tau ,\sigma \in [0,2]$ as the other quadrands are linked through discrete symmetries. $J$ is nonzero only below the diagonal (dashed line) as $E\left(t\right)$ is nonzero on a bounded domain $t\in [-1,1]$ by definition. $E\left(t\right)$ is generally complex and normalized such that $\text{max}\left|J\right[E\left]\right(\tau ,\sigma \left)\right|=1$.

Figure 2

Illustrating example: Synthetic measurement trace with $129\times 129$ data points on $\tau ,\sigma \in [0,2]$. Coarse-grained data (only imaginary part shown) on $21\times 21$ pixels (40 data points per pixel) enable fast computation of approximants to initialize refined retrievals. Every pixel of the lower triangular part ($21\times (21+1)/2$) is associated with one equation in (16). If $E^{+}$ and $E^{-}$ are real roots, then ${\langle J{[E^{+},E^{-}]}_{ij}\rangle }^{+}$ and ${\langle J{[E^{+},E^{-}]}_{ij}\rangle }^{-}$ are real and equivalent to $\text{Re}\left[\langle J{[E^{+},E^{-}]}_{ij}\rangle \right], \text{Im}\left[\langle J{[E^{+},E^{-}]}_{ij}\rangle \right]$. Note: The exponent $1/4$ is convenient for data examination as (12) constitutes a 4th order polynomial system in the components $E_k^{+},E_k^{-}$.

Figure 3

Solving the polynomial system (18) for a test case on $N=15$. Top left: Path tracking the solution curve consisting of small path segments each corresponding to a single reduced system with fixed random matrix $M$. Top right: Global distance to target trace of full (colored) and reduced (gray) system decreases approximately exponentially (on each segment linearly). Bottom left: Local residual after predictor step (colored) and after corrector step (Newton step) (gray). Bottom right: Relative number of up paths (increasing $\parallel \mathrm{\Delta }C\parallel$) and valid paths (successful Newton step) when doing trial steps at the end of each segment to find a new path (and new $M$) along which $\parallel \mathrm{\Delta }C\parallel$ decreases.

Figure 4

Solving the polynomial system including a Tikhonov-type penalty term. Top right: The size of $\lambda$ is constant along each path segment (see Fig. 3) and adaptively decreased at the end of each segment, if a threshold on the relative size of the regularization term is crossed ${\delta }_{\mathrm{\Delta }C}^{\text{regul}}<20\%$ [Eq. (29)] for three different noise levels 0% (green), 1% (yellow), 2% (blue). Bottom left: The relative size of the penalty term $K$ on each segment is growing and globally adaptively reduced. Top left: The distances to the target trace with (opaque colored) and without (full colored) penalty term are simultaneously decreasing along each path segment, while their difference is increasing because $F$ is deformed slightly towards an improved match with the target trace $C_1$ for a class of curves with similar mean curvature. Then in the noiseless case (green) [noisy case (yellow, blue)], ${\delta }_{\mathrm{\Delta }C}^{\text{regul}}$ vanishes [settles] where $\lambda$ is near optimal. Bottom right: Beginning at a zero phase Gaussian $E(s=0)$, connected by a series of smooth intermediate solutions the path tracker continues towards the target pulse (dashed gray, only first 200 iterations shown).

Figure 5

Performance overview: Retrieval time vs $N$ with $1\%$ additive noise and three different values for the termination criterion $p_{\text{up}}^{\text{stop}}$ and the resulting retrieval accuracy, number of iterations, and retrieval accuracy after applying additional refinement steps all averaged over 10 runs. Initial data for one level is the interpolated result from the next coarser level. The plot implies choosing as smaller $p_{\text{up}}^{\text{stop}}$ while cascading towards the final grid and do refined retrievals there.

Figure 6

Error convergence vs $N$: Average (full colored points) final pulse error $\epsilon$ (a, c, d) and trace error $\parallel \mathrm{\Delta }C\parallel /\parallel C\parallel$ (b) on 14 pixelizations for the test cases TBP1 and TBP2 with and without noise. Two examples of the retrieved electric field on a coarse (e) and medium (f) pixelized grid for TBP1. (g) vs (h): Comparison of our solver with the result of a least-squares solver without regularization. For more details see Sec. 6c. Note: Here more intermediate grids then necessary are used. Normally, a cascade like $N\to 20\to 40\to 65$ is sufficient.

Figure 7

Left: Applying the L-curve method for nine different values of $\lambda$, 10 retrievals each (opaque colored), mean value (full colored). Noisy oscillations begin to increase the length of the pulse $\parallel X\parallel$ when decreasing $\lambda$; the trace error $\parallel \mathrm{\Delta }C\parallel$ improves through overfitting. Too smooth solutions (too large $\lambda$) show deviations from the original solution and have larger trace errors. The optimal amount of $\lambda$ is within the corner of the L-curve where the pulse error $\epsilon$ is small. Right: 15 randomly offset retrievals of pulse TBP2 (gray lines) on $N=65$ intervals. As the position of the pulse is not fixed relative to the numerical grid; sampling on arbitrary intergrid locations is possible.