Algorithms
Internally one can divide the implemented algorithms into two general types. Classical-Algorithms and General-Optimization-Algorithms, where the former refers to algorithms that have been specifically introduced to solve the phase-retrieval problem in ultrafast optics, while the latter refers to generic optimization algorithms.
The algorithms described here belong to either type. In addition several method-specific algorithms have been published and are discussed in Method-Specific Algorithms.
Sidenote: Even though almost none of the algorithms require the delay and frequency axis to be FFT-conform, it seems that this increases the performance.
Classical Algorithms
Generalized Projection [2, 12]
Starting from an initial guess the algorithm works by constructing the nonlinear signal field \(S(\theta, \omega)\) (\(\theta\) refers to an arbitrary scanning parameter). The amplitude of this guessed signal field is replaced by the measured amplitude of the trace to yield \(S'(\theta, \omega)\). Thus the signal field is projected into the space of signal fields that could produce the measured trace. The updated guess of the pulse is obtained through an iterative optimization of the Z-error.
\(Z = \sum |S'(\theta, t) - S(\theta, t)|^2\)
\(E'(\omega) = E(\omega) - \frac{dZ}{dE^*(\omega)}\)
Below is a small working example of the Generalized Projection Algorithm to an SHG-FROG. Apart from the nonlinear optimization methods described in the Nonlinear Optimization section, the main tunable parameters are the number of gradient descent steps (no_steps_descent) and the step size (global_gamma).
from pulsedjax.frog import GeneralizedProjection
gp = GeneralizedProjection(delay, frequency, trace, "shg")
population = gp.create_initial_population(5, "random")
gp.global_gamma = 1
gp.no_steps_descent = 15
final_result = gp.run(population, 50)
Ptychographic Iterative Engine [7, 13, 14]
The general idea of Ptychography is to move a probe \(P(\theta)\) over an object \(O\). The resulting signal is simply the product of the two. Thus the Ptychographic Iterative Engine (PIE) demands the nonlinear signal to be of the form \(S(\theta, t) = O(t)\cdot P(\theta,t)\). This makes it suitable for most pulse characterization methods, with the exception inteferometry based methods such as i-FROG and doubleblind retrieval of 2D-SI and VAMPIRE.
The PIE can be described in two equivalent ways. Either as an iterative optimization of the regularized Z-error:
\(Z = \sum |S'(\theta, t) - O'(t)\cdot P(\theta, t)|^2 + \sum U(t)\cdot |O'(t) - O(t)|^2\)
with the heuristic weighting \(U(t)\).
Alternatively it can be described by the weighted iterative optimization of:
\(L = \sum |\sqrt{T(\theta, \omega)} - |S(\theta, \omega)||^2\)
via
\(O'(t) = O(t) - U'(t)\cdot \frac{dL(\theta)}{dO^*(t)}\)
where \(U'\neq U\) and determines the subtype of the PIE.
Usually the PIE works in local iterations, where optimization steps are taken subsequently over random permutations of the measured trace along \(\theta\). In this implementation an additional version performs global steps by optimizing for all measured data along \(\theta\) simultaneously.
Below is a small working example of the PIE to an SHG-FROG. Apart from the nonlinear optimization methods described in the Nonlinear Optimization section, the main tunable parameters are the PIE version, the weighting parameter \(\alpha\), as well as the local and global stepsizes.
from pulsedjax.frog import PtychographicIterativeEngine
pie = PtychographicIterativeEngine(delay, frequency, trace, "shg", pie_method="PIE")
population = pie.create_initial_population(5, "random")
pie.alpha = 0.15
pie.local_gamma = 1e-1
pie.global_gamma = 1e-3
final_result = pie.run(population, 50, 50)
COPRA [15, 16]
The original implementation of the Common-Pulse-Retrieval-Algorithm (COPRA) is pypret.
Roughly speaking, COPRA is a combination of the Generalized-Projection (GP) and the PIE. The algorithm works in two phases. In the first phase, similar to the PIE, local iterations are performed on random permutations along \(\theta\). However instead of using an PIE update, COPRA performs single steps to optimize the Z-error. After a sufficient number of local iterations COPRA switches to global iterations. A global iteration optimizes two different error quantities in an alternating fashion. These are the r-error:
\(r = \sum |T(\theta,\omega) - \mu\cdot |S(\theta,\omega)|^2|^2\)
and the Z-error. Using single gradient-descent iterations the minimization of r-error is used to obtain updated signal fields \(S'(\theta,t)\) which are subsequently used to update the pulse guess via the minimization of the Z-error using single gradient-descent iterations. Thus in a global iteration smooth and gradual transitions are employed over plain projections.
\(S'(\theta, \omega) = S(\theta, \omega) - \frac{dr}{dS^*(\theta,\omega)}\)
\(E'(\omega) = E(\omega) - \frac{dZ}{dE^*(\omega)}\)
In addition for all gradient-descent iterations COPRA estimates an optimal step size via the Polyak-Stepsize.
Below is a small working example of the COPRA to an SHG-FROG. Apart from the nonlinear optimization methods described in the Nonlinear Optimization section, the main tunable parameters are the local and global stepsizes.
from pulsedjax.frog import COPRA
copra = COPRA(delay, frequency, trace, "shg")
population = copra.create_initial_population(5, "random")
copra.local_gamma = 1
copra.global_gamma = 0.25
final_result = copra.run(population, 50, 50)
LSF [17]
The Linesearch-Frog-Algorithm (LSF) is a gradient-free algorithm. Despite its name it works with all measurement techniques. Simply speaking the algorithm performs a high-dimensional randomized bisection search. This works by selecting a random complex-valued direction. By setting the maximum pulse magnitude to one, a lower and an upper boundry can be defined. By performing bisection search along the random direction inside the boundary, a slightly improoved pulse guess can be found.
The original algorithm only uses bisection, however it is straight forward to disect the search direction into N slices. This should in theory reduce the required number of iterations in the sectional searches.
In addition on top of the fully randomized version, in this implementation one can restrict the random directions to exhibit some degree of smoothness. However, this may not necesserily improve the algorithm performance.
Below is a small working example of the LSF to an SHG-FROG. The main tunable parameters are the number of iterations per bisection-search as well as a tuning parameter for the smoothness in case of a continuous direction search.
from pulsedjax.frog import LSF
lsf = LSF(delay, frequency, trace, "shg")
population = lsf.create_initial_population(5, "random")
lsf.no_sections = 2
lsf.number_of_disection_iterations = 12
lsf.direction_mode = "random"
# lsf.ratio_points_for_continuous = 0.25 # only used when random_direction_mode="continuous"
final_result = lsf.run(population, 50)
Additional Constraint and Momentum [14]
All of the Classical Algorithms, are able to incorporate the pulse-spectrum as an additional constraint. This may improve the retrieval result, but it also may possibly cause issues if the provided spectrum and measured trace are somehow inconsistent.
For GeneralizedProjection, COPRA and LSF the spectrum is used as a true constraint, meaning that the phase factor is optimized while the amplitude is obtained from the measured spectrum.
A spectrum can be incorporated in any algorithm as shown below.
from pulsedjax.frog import GeneralizedProjection
gp = GeneralizedProjection(delay, frequency, trace, "shg")
gp.use_measured_spectrum(frequency_spectrum, pulse_spectrum, pulse_or_gate="pulse")
population = gp.create_initial_population(5, "random")
gp.global_gamma = 1
gp.no_steps_descent = 15
final_result = gp.run(population, 50)
Additionally the gradient-descent like optimizations can be accelerated via the incroporation of momentum. This means that a velocity-map of the pulses is updated each iteration and used to additionally update the current guess. This can be used with any Classical Algorithm as shown below.
from pulsedjax.frog import GeneralizedProjection
gp = GeneralizedProjection(delay, frequency, trace, "shg")
gp.momentum(eta=0.9)
population = gp.create_initial_population(5, "random")
gp.global_gamma = 1
gp.no_steps_descent = 15
final_result = gp.run(population, 50)
General Optimization Algorithms
Differential Evolution [18, 19]
This implementation of the Differential Evolution algorithm (DE), follows the same paper as the scipy-implementation. It works by taking a parent population and creating a child population through mutations and crossover operations. A fitness evaluation then returns the population back to its original size.
In addition to the standard implementation an additional crossover strategy (smooth) and an additional selection mechanism (global) have been implemented. The smooth-crossover is motivated by the fact that the population can be parametrized on a grid. Conventional crossovers (bin, exp) will introduce sharp jumps and cutoffs into the population. The smooth crossover works similar to the exp-crossover but introduces a smooth transition via a tanh() function.
The additional selection mechanism is motivated by the fact that the greedy-selection performs a pair-based comparison of parent and child population. The global-selection performs a global comparison of all individuals and randomly selects the fittest individuals. This selection can be tuned via a temperature of a Fermi-Dirac-like probability distribution.
Below is a small working example of the DE to an SHG-FROG. The main tunable parameters are the strategy, selection mechanism, the mutation rate, the crossover rate and in the temperature in case of a global selection.
from pulsedjax.frog import DifferentialEvolution
de = DifferentialEvolution(delay, frequency, trace, "shg",
strategy="best1_bin", selection_mechanism="greedy",
mutation_rate=0.5, crossover_rate=0.7)
population = de.create_initial_population(100,
amp_type="bsplines_5", phase_type="bsplines_5",
no_funcs_amp=20, no_funcs_phase=20)
final_result = de.run(population, 50)
Evosax [20, 21]
This “Algorithm” uses the python/jax package evosax to employ advanced evolutionary solvers for pulse retrieval. The availablilty of a multitude of such solvers comes at a drawback. Evosax does not respect the shape of pytrees. Internally all pytrees are flattened before being subjected to an evolutionary algorithm. This causes mixing of variables that should not necessarily be mixed. For example the position and the magnitude of a gaussian. In this implementation precautions have been taken to prevent the mixing of amplitude and phase variables. But within the phases and amplitudes evosax will happily mixing variables arbitrarily. Essentially this heavily promotes exploration over exploitation.
Below is a small working example of Evosax with an SHG-FROG. One should be able to choose any evosax solver. These solvers can be tuned via their input parameters.
from pulsedjax.frog import Evosax
from evosax.algorithms import DifferentialEvolution as DifferentialEvolutionEvosax
# the evosax-solver should not be initialized
evo = Evosax(delay, frequency, trace, "shg", solver=DifferentialEvolutionEvosax)
population = evo.create_initial_population(100,
amp_type="gaussian", phase_type="sigmoidal",
no_funcs_amp=5, no_funcs_phase=15)
evo.solver_params = None # None, causes the default params to be used.
evo.solver_kwargs = {} # Some solvers need non-standard input args/kwargs at initialization
final_result = evo.run(population, 50)
AutoDiff [22, 23, 24]
This “Algorithm” uses the python/jax packages optax and optimistix to employ their implemented optimization algorithms.
Below is a small working example of AutoDiff with an SHG-FROG. One should be able to choose any optax/optimistix solver. These solvers can be tuned via their input parameters.
from pulsedjax.frog import AutoDiff
import optimistix
import optax
solver = optax.adam(learning_rate=1e-1)
# solver = optimistix.LBFGS(1, 1) # the tolerances dont matter, since a fixed number of iterations is done
ad = AutoDiff(delay, frequency, trace, "shg", solver=solver)
population = ad.create_initial_population(5,
amp_type="gaussian", phase_type="sigmoidal",
no_funcs_amp=5, no_funcs_phase=15)
final_result = ad.run(population, 50)
Additional Constraint and Loss-Function
As in the Classical Algorithms, the General Algorithms are able to incorporate the pulse-spectrum as an additional constraint. All implementations incorporate this constraint as a true constraint optimization. Meaning that if a spectrum is provided only the phase will be optimized. For doubleblind retrievals the provision of spectra for pulse and gate-pulse is necessary to avoid retrieval ambiguities.
A spectrum can be incorporated in any algorithm as shown below.
from pulsedjax.frog import AutoDiff
import optimistix
import optax
solver = optax.adam(learning_rate=1e-1)
# solver = optimistix.LBFGS(1, 1) # the tolerances dont matter, since a fixed number of iterations is done
ad = AutoDiff(delay, frequency, trace, "shg", solver=solver)
ad.use_measured_spectrum(frequency_spectrum, pulse_spectrum, pulse_or_gate="pulse")
population = ad.create_initial_population(5,
amp_type="gaussian", phase_type="sigmoidal",
no_funcs_amp=5, no_funcs_phase=15)
final_result = ad.run(population, 50)
Since the General Algorithms dont rely on analytic derivatives one is not restricted to a specific Error/Loss-function. By default the G-error is used:
\(G = \frac{1}{N_\theta\cdot N_\omega}\sum |T(\theta, \omega) - |S(\theta,\omega)|^2|^2\)
However one can specify an arbitrary user-defined loss-function [25].
from pulsedjax.frog import AutoDiff
import optax
import jax.numpy as jnp
def my_loss_function(trace, measured_trace, exponent_trace=1, L_norm=1):
"""
An example of a modified loss function for general algorithms.
Inspired by ELFPIE: S. Zhang, T.T.J.M. Berendschot and J. Zhou, Signal Processing 210 109088, 10.1016/j.sigpro.2023.109088, (2023)
Args:
trace (jnp.array): the current trace
measured_trace (jnp.array): the measured trace
kwargs:
exponent_trace (int, float): the exponent applied to trace before building the residual
L_norm (int, float): the norm exponent to be used
"""
measured_trace = measured_trace/jnp.max(jnp.abs(measured_trace))
trace = trace/jnp.max(jnp.abs(trace))
residual = jnp.sign(trace)*jnp.abs(trace)**(exponent_trace) - jnp.sign(measured_trace)*jnp.abs(measured_trace)**(exponent_trace)
residual = jnp.gradient(residual)
error = jnp.sum(jnp.abs(residual)**L_norm)
return error
solver = optax.adam(learning_rate=1e-1)
ad = AutoDiff(delay, frequency, trace, "shg", solver=solver)
population = ad.create_initial_population(5,
amp_type="gaussian", phase_type="sigmoidal",
no_funcs_amp=5, no_funcs_phase=15)
ad.error_metric = my_loss_function
final_result = ad.run(population, 50)
Cross-correlation and Doubleblind Retrievals
All algorithms are capable of retrieving the pulse from a cross-correlation measurement, given that the gate pulse has been provided. In addition the keyword cross_correlation=True needs to be set during algorithm instantiation.
myalgorithm = MyAlgorithm(theta, frequency, trace, nonlinear_method, cross_correlation=True, ... )
myalgorithm.get_gate_pulse(frequency_gate, gate_pulse)
If the gate-pulse is not known a doubleblind-retrieval of both the pulse and the gate-pulse can be attempted with most algorithms. It can be enabled by setting cross_correlation="doubleblind". Those algorithms who are not meant for doubleblind retrievals will/should raise an assertion error if this is attempted.
The simultaneous retrieval of pulse and gate-pulse is non-unique, thus additional constraints in the form of the pulse-spectrum and gate-pulse spectrum should be provided. Examples for cross-correlation and doubleblind retrievals can be found here.
Some algorithms like CPCGPA retrieve the actual gate and not the gate pulse. In such cases the provision of a spectrum for the gate-pulse would introduce inconsistencies. However this algorithm incorporates the required additional constraints through the “principal-component procdure”.
Retrieval of the Calibration Curve
Usually, the calibration curve of a spectrometer or an optical setup is known. In such cases it can be provided via algorithm.get_calibration_curve(frequency, calibration_curve). However in rare cases the calibration curve may not easily be accessible. Luckily the calibration curve \(\mu\) can be recovered via linear least-squares from the G-error. [26]
\(G = \sum\limits_{mn} (T_{mn}^{\mathrm{exp}} - \mu_n\cdot T_{mn})^2\)
\(\quad \Rightarrow \quad \mu_n = \dfrac{\sum_m T_{mn}^{\mathrm{exp}} \cdot T_{mn}}{\sum_m T_{mn}^2}\)
Thus for a given pulse guess the trace and accompanying calibration curve can be generated. The calibration curve itself is used in the iterative update of the pulse guess. This feedback/additional degree of freedom requires additional constraints on the optimization in order to yield useful results. These can be provided through a calibrated spectrum of the pulse.
# provision of known calibration curve, can also be directly multiplied to meausured_trace before algorithm initialization
myalgorithm.get_calibration_curve(frequency, calibration_curve)
# optimize calibration curve in local or global stages
myalgorithm.local_optimize_calibration_curve = True
myalgorithm.global_optimize_calibration_curve = True
In practice it seems that the calibration curve retrieval tends to be only successfuly with iterative updates to the signal field \(S'(\theta,t)\). Thus, one should enable this approach via:
myalgorithm.r_local_method = "iteration"
myalgorithm.r_global_method = "iteration"
Possible reasons for this are that a plain projection likely leads to a more eratic movement on the error-surface. Alternatively, during the projection \(\mu\) is incorporated via division, which amplifies noise.