GEOPE – Geodesic Pulse Engineering¤

Overview¤

geope finds piecewise-constant control pulses that implement a target quantum gate on an $n$-qubit system. Given a target unitary $U_T$, a set of available control generators (the projected basis), and optionally fixed drift generators, the optimiser searches for real-valued parameters $\phi$ such that

\[ U(\phi) \;=\; \prod_{g=1}^{N_g} \exp\!\Bigl(i \sum_{k}\phi_{g,k}\,G_k\Bigr) \;\approx\; U_T, \]

where each $H_g = \sum_k \phi_{g,k}\,G_k$ is a linear combination of basis generators on segment $g$.

The core algorithm is the geodesic method: at each step it computes the shortest path on $U(d)$ from the current unitary to the target, projects that direction onto the controllable subspace, then solves a convex least-squares problem and a one-dimensional line search to take a parameter step. This is distinct from gradient-based methods like GRAPE that follow the fidelity gradient directly.

The entry point is Parameters — a state object that bundles every input the optimiser needs (basis, control, drift, target, constraints, pulse constraints, param_transform, bounds, init values, seed, projective flag). Pass it to Geope and call .optimize(max_steps=...). The returned Parameters carries the live/final parameters and fidelity (and to_dict()); the full run trajectory and best_* helpers live on an opt-in History logger (geope.history).

The optimisation functions themselves (unitary computation, fidelity, Jacobian, geodesic Hamiltonian, gammas/omegas, Hessian) are pure function factories in engine.py. They are built lazily and cached on the Parameters object on first use, so there is no separate engine object — Geope, Grape and Gecko all read the functions straight off the (shared) Parameters. The factories return un-jitted callables; JIT compilation happens once, when the optimiser's update_step is first traced inside optimize().

Class hierarchy¤

Basis, Hamiltonian, Unitary      (lie.py)

engine.py  — pure function factories (get_jacobian_fn, get_geodesic_hamiltonian_fn,
             get_gammas_and_omegas_fn, get_hessian_fn, fidelity helpers, …)
                ↓ built lazily & cached on
Parameters                       (parameters.py)
                ↘
                  Geope / Grape  (geope.py / grape.py)
                  Gecko          (gecko.py)

Lie group classes (`lie.py`)¤

`Basis`¤

Represents a set of Lie algebra generators (e.g. Pauli strings) as a rank-3 tensor of shape $(K, d, d)$.

Basis(basis, labels=None, local_dim=2, n_qubits=None,
      interaction_graph=None, interaction_map=None)

Parameter	Description
`basis`	`np.ndarray` of shape `(K, 2ⁿ, 2ⁿ)` — Hermitian generators
`labels`	list of Pauli-string labels, e.g. `["XI", "IX", "ZZ"]`
`local_dim`	local Hilbert-space dimension, default 2
`n_qubits`	override for qubit count when $d \neq 2^n$
`interaction_graph`	list of qubit tuples to keep, e.g. `[(1,2), (2,3)]`
`interaction_map`	dict of qubit-tuple → allowed interaction labels

Key properties:

Property	Description
`basis`	the `(K, d, d)` tensor
`lie_algebra_dim`	$K$ — number of generators
`dim`	$d$ — matrix dimension
`n`	number of qubits
`labels`	string labels
`plot_labels`	LaTeX strings, e.g. `"$X_{1}Z_{2}$"`
`interaction_qubits`	tuple of qubit indices for each generator
`interaction_graph`, `interaction_map`	as above

Key methods:

overlap(other) — boolean mask over other's basis, true where there is nonzero trace overlap with self. Used by Parameters to build its index masks.
verify() — orthogonality check under the trace inner product.
linear_span(parameters) — $\sum_i \phi_i G_i$.
generate_parameter_list(parameter_map) — converts a dict like {1: {"x": 0.5}, (1,2): {"zz": 0.3}} to a flat parameter array.
generate_bounds(bounds_map, piecewise_steps) — converts {"x": (-1, 1)} to (lower, upper) arrays.
apply_interaction_graph(graph) / apply_interaction_map(map) — prune to hardware connectivity.

`Hamiltonian`¤

Represents $H = \sum_i \phi_i G_i$ and its unitary $U = e^{iH}$.

Hamiltonian(basis, parameters)

Attribute	Description
`basis`	the `Basis` object
`parameters`	coefficient vector $\phi$
`matrix`	$\sum_i \phi_i G_i$
`unitary`	`Unitary(expm(i·matrix))`

Methods:

geodesic_hamiltonian(target_unitary) — a Hamiltonian whose parameters are the geodesic direction $-i\log(U^\dagger U_T)$ decomposed in the basis.
fidelity(unitary_matrix) — $|\mathrm{Tr}(U^\dagger V)|/d$.
parameters_from_hamiltonian(H, basis) (static) — coefficients via $\mathrm{Re}\,\mathrm{Tr}(G_i H)/d$.

`Unitary`¤

Wraps a unitary matrix with validation. Validates $UU^\dagger = I$ on construction.

parameters(basis) — Lie-algebra coefficients via the principal logm.
fidelity(other).
geodesic_hamiltonian(basis, target).

Basis construction utilities (`utils.py`)¤

Function	Description
`construct_full_pauli_basis(n)`	all $4^n - 1$ non-identity Pauli strings
`construct_two_body_pauli_basis(n)`	1-body and 2-body terms only
`construct_Heisenberg_pauli_basis(n)`	1-body + same-type 2-body (XX, YY, ZZ)
`construct_restricted_pauli_basis(n, restriction)`	custom restriction (list or dict)
`construct_full_spin_boson_basis(n_spins, n_bosons, truncation)`	spin-boson hybrid
`construct_restricted_spin_boson_basis(...)`	restricted spin-boson
`filter_basis_by_control(basis, control)`	filter an existing `Basis` by a control dict (handy when $d \neq 2^n$)

Restriction formats¤

construct_restricted_pauli_basis accepts two formats.

List — allowed interaction types as lower-case strings:

control = geope.construct_restricted_pauli_basis(2, ['x', 'z'])
control = geope.construct_restricted_pauli_basis(3, ['x', 'y', 'z'])
drift   = geope.construct_restricted_pauli_basis(3, ['zz'])

Dict — allowed interactions per qubit or qubit pair (1-indexed):

control = geope.construct_restricted_pauli_basis(2, {1: ['x'], 2: ['x'], (1,2): ['zz']})

Drift parameter values¤

Drift coefficients are specified via generate_parameter_list on the drift basis (or passed directly through Parameters(drift_values=...)):

drift_basis = geope.construct_restricted_pauli_basis(3, ['zz'])
drift_values = drift_basis.generate_parameter_list({
    (1, 2): {"zz": 1.0},
    (2, 3): {"zz": 1.0},
    (1, 3): {"zz": 1.0},
})
# → [1.0, 1.0, 1.0] matching basis order ["ZZI", "ZIZ", "IZZ"]

Linear equality constraints (global controls)¤

Constraints enforce that selected projected parameters maintain fixed ratios. Use generate_parameter_list to build constraint vectors:

control = geope.construct_restricted_pauli_basis(3, ['x', 'z'])

global_x = control.generate_parameter_list({1: {"x": 1}, 2: {"x": 1}, 3: {"x": 1}})
# → ties X₁ = X₂ = X₃

global_z = control.generate_parameter_list({1: {"z": 1}, 2: {"z": 1}, 3: {"z": 1}})
# → ties Z₁ = Z₂ = Z₃

Pass via Parameters(constraints=[global_x, global_z], ...).

Pulse-shape constraints¤

Pulse constraints fix the relative values of specified parameters across piecewise steps — the temporal shape is frozen while the overall scale is optimised. This is an alternative to the drift_basis + drift_values route when you want a drift-like term whose amplitude is still tuned by the optimiser.

projected = geope.construct_restricted_pauli_basis(3, ['x', 'z', 'zz'])
params = geope.Parameters(
    basis=geope.construct_full_pauli_basis(3),
    control={1: ['x', 'z'], 2: ['x', 'z'], 3: ['x', 'z'],
             (1, 2): ['zz'], (2, 3): ['zz'], (1, 3): ['zz']},
    target=U_T,
    piecewise_steps=10,
    pulse_constraints={(1, 2): ['zz'], (2, 3): ['zz'], (1, 3): ['zz']},
)

pulse_constraints uses the same {qubit_index_or_tuple: [interaction]} dict format as control — here it freezes the temporal shape of the three zz terms.

Approach	Use case
`drift_basis` + `drift_values`	drift is truly fixed and not optimised
`pulse_constraints` on projected params	drift-like terms whose amplitude is optimised but whose temporal profile is fixed

Other utilities:

prepare_random_parameters(proj_indices, expander, spread, seed) — random initial parameters respecting constraints.
golden_section_search(f, a, b, tol) — JIT-compatible 1-D minimiser (used internally by the line search).
adam_line_search(f, a, b, lr, num_steps, finite_difference) — JIT-compatible 1-D Adam minimiser (finite-difference or exact-gradient; used by the "adam*" line-search methods).
merge_constraints(constraints) — merges overlapping linear constraints.
qft_unitary(n), multicontrol_unitary(U, n_controls) — common target unitaries.
make_per_element_transform(transforms) — helper to build a param_transform from per-element callables.

The `Parameters` object¤

Parameters is the recommended entry point.

Parameters(basis=None, control=None, drift=None,
           init_values=None, drift_values=None,
           target=None, piecewise_steps=1, fixed_drift=True,
           constraints=None, pulse_constraints=None, bounds=None,
           init_spread=0.1, seed=None,
           param_transform=None, n_experimental_params=None,
           projective=True)

Parameter	Description
`basis`	the full `Basis`; defaults to 2-qubit full Pauli basis if `None`
`control`	dict picking the projected (controllable) subset
`drift`	dict picking the drift subset
`init_values`	dict in `control` format, or `ndarray` of full-basis shape, or `None` (random)
`drift_values`	dict, `ndarray`, or `None` (ones)
`target`	target unitary as `ndarray`
`piecewise_steps`	number of gate segments $N_g$
`fixed_drift`	whether drift is held fixed during optimisation
`constraints`	list of constraint vectors / dicts
`pulse_constraints`	control-format dict `{site: [ops]}` (same format as `control`) of projected terms whose time-shape is fixed
`bounds`	dict `{label: (lo, hi)}` — consumed by `Geope.bound(...)`, not by the main loop
`init_spread`	half-width of uniform random init, in units of $\pi$
`seed`	random seed
`param_transform`	callable mapping experimental params to basis coefficients
`n_experimental_params`	length of the experimental input; defaults to `projected_basis.lie_algebra_dim`
`projective`	`True` (default) for projective fidelity, `False` for phase-sensitive

Attributes populated after construction:

Attribute	Description
`basis`, `projected_basis`, `drift_basis`	the three `Basis` objects
`target`	the target
`drift_parameters`	drift coefficients (or `None`)
`constraint_arrays`, `constraint_expander`	merged constraints and reduced-space mapping
`bounds`	pre-built bounds tuple (or `None`)

Live optimisation state — seeded at construction and updated in place by Geope:

Attribute	Description
`parameters`	current parameter array, shape `(N_g, K_{full})`; seeded to the initial guess, holds the final result after `optimize()`
`fidelity`	current fidelity (`None` before a run)
`infidelity`	`1 - fidelity` (`None` before a run)
`basis_coefficients`	current parameters mapped through `param_transform` if set
`to_dict()`	current solution as a control-style dict

The full run trajectory and the best_* helpers live on the opt-in History logger, not on Parameters.

Metadata, functions, and optimiser¤

Algebraic metadata (cached on `Parameters`)¤

Parameters derives the index masks relating the three bases as cached properties (computed once from basis / projected_basis / drift_basis):

projected_indices — shape $(K_{\text{full}},)$, which full-basis elements are controllable (projected_basis.overlap(basis)).
drift_indices — shape $(K_{\text{full}},)$, which are fixed drift.
proj_drift_indices = projected_indices | drift_indices.
proj_indices_projdrift_basis — projected mask within the proj+drift subspace, shape $(K_{\text{pd}},)$.
drift_indices_projdrift_basis — drift mask within the proj+drift subspace.
proj_drift_basis — Basis containing only the projected + drift elements; used for all JIT computations.

Optimisation functions (`engine.py`, cached on `Parameters`)¤

The optimisation primitives are pure function factories in engine.py. Each is built lazily and cached on the Parameters object on first access, so the optimisers read them as params.<name> (there is no engine object). They return un-jitted callables that fuse into the optimiser's @jax.jit update_step — so compilation happens once, on the first optimize() call, not at construction.

params.compute_U_fn(params_list) — scans over piecewise steps via jax.lax.scan: $\,U = \prod_g \exp\!\bigl(i \sum_i \phi_{g,i}\,G_i\bigr).$ Input shape: $(N_g, K_{\text{pd}})$.
params.fid_U_fn(U) — $|\mathrm{Tr}(U_T^\dagger U)|/d$ when projective=True, $\mathrm{Re}\,\mathrm{Tr}(U_T^\dagger U)/d$ when projective=False.
params.project_omegas_fn — projects matrices onto the full Pauli basis via trace inner products $\mathrm{Tr}(G_i M)$. For $n > 5$ uses on-the-fly batched projection to manage memory.
params.jac_fn — Jacobian $\partial U/\partial\phi_{g,k}$ (JAX autodiff, holomorphic=True; a real/imag-split Jacobian under param_transform).
params.geo_fn — geodesic tangent at $U$: $U \cdot \bigl(-i\log(U^\dagger U_T)\bigr)$, with the global-phase generator subtracted when projective=True.
params.gammas_and_omegas — the combined geodesic-coefficients-and-Jacobian-projection used by the GEOPE update step.
params.infid_U_fn, params.infid_fn, params.grad_fn, params.hess_fn — infidelity helpers; grad_fn/hess_fn are used by GRAPE.

Because the functions are cached on Parameters, sharing one Parameters between a Geope and a Gecko reuses the identical callables — so JAX reuses the compiled traces instead of recompiling.

`Geope` (`geope.py`)¤

Geope(params, verbose=False, history=None)

Geope requires a Parameters object as its single positional argument. The optimisation functions, initial parameters, drift, constraints, pulse constraints, seed, initialisation spread, projective flag and param_transform are all read from params. Passing anything other than a Parameters raises TypeError.

Parameter	Description
`params`	a `Parameters` instance bundling all inputs
`verbose`	print per-step progress
`history`	optional `History` logger (`None` = no logging)

The iteration cap, the line search, and the three run-control knobs are arguments of optimize, not constructor fields:

from geope import adam, GoldenSection

optimize(max_steps=1000,
         line_search=GoldenSection(),        # default; or adam(1e-2)
         precision=0.9999999,
         max_step_size=0.9, gram_schmidt_step_size=1.3)

`optimize` argument	Description
`max_steps`	maximum number of optimisation steps
`line_search`	a `LineSearch` object tuning the geodesic step size; defaults to `GoldenSection()`
`precision`	target fidelity threshold (host-side; no recompile)
`max_step_size`	maximum line-search step (baked into the JIT; changing it recompiles)
`gram_schmidt_step_size`	step size for the Gram–Schmidt fallback (host-side; a falsy value disables it)

The line searches are immutable config objects (frozen dataclasses):

GoldenSection(tol=1e-5) — golden-section search (the default), stateless.
adam(lr=0.05, num_steps=30, finite_difference=True, warm_start=False, ...) — 1-D Adam line search. finite_difference=False uses an exact autodiff gradient; warm_start=True seeds each step from the previous step's t.

The line-search object and max_step_size bake into JIT-compiled functions that optimize builds on first use and reuses across calls; the frozen-dataclass value equality means two equal line searches (e.g. the per-call default GoldenSection()) reuse the compiled functions, while changing the object or max_step_size triggers a one-off recompile. precision and gram_schmidt_step_size are host-side only — changing them never recompiles.

Live state and logging:

The current parameters and fidelity live on params (params.parameters, params.fidelity); Geope updates them in place each step, and optimize(max_steps=...) returns the Parameters instance itself — so the user has a single handle for both inputs and the final result.
step_size — the transient last line-search step size.
history — an optional History logger (None unless one was passed). When supplied, the full run trajectory and best_* helpers are available on it (see below).

`History` (`history.py`)¤

History(logging_fn=None)

An opt-in, configurable run log. Pass one to Geope (history=History()) and the full trajectory is recorded into geope.history; leave it None and no history is kept (the final answer still lives on params).

By default each step records five columns — parameters (a full-basis snapshot), fidelities, infidelities, step_sizes, and an integer steps counter derived from the log length. Pass logging_fn to record arbitrary per-step values instead: it receives the running Geope and returns a dict of column -> value (e.g. History(logging_fn=lambda g: {"fid": float(g.params.fidelity)})).

Member	Description
`record(geope)`	append one row via `logging_fn`; called by `Geope` each step
`reset()`	drop all rows
`len(history)`	number of recorded rows
`history.<col>` / `history["<col>"]`	a logged column (the same list)
`keys()`	the logged column names
`to_dataframe()`	the logs as a `pandas.DataFrame`
`best_fidelity`	`max(fidelities)` (or `None`)
`best_parameters`	parameters at the highest-fidelity step (or `None`)
`best_basis_coefficients`	best parameters mapped through `param_transform` if set
`to_dict()`	best solution as a control-style dict (`{}` if unavailable)

The best-over-trajectory helpers need the default fidelities/parameters columns; under a custom logging_fn that omits them they degrade to None/{} rather than raising. Note params.parameters is the single current array while history.parameters is the list of per-step snapshots. A History is meant for a single run.

Core algorithm: `optimize()`¤

for each step:
    1. Extract free_params = parameters[:, proj_drift_indices]

    2. Compute the geodesic direction:
       U  = compute_U_fn(free_params)
       g  = -i · logm(U† U_T)                       # generator in u(d)
       g  = g - Tr(g)/d · I        if projective    # drop global-phase generator
       Γ  = U · g                                   # geodesic tangent
       γ  = project(Γ) / d                          # coefficients in basis

    3. Compute the Jacobian projections:
       ω[g, k] = project(i · ∂U/∂φ_{g,k})

    4. Solve the constrained least-squares problem:
       sol = argmin ||ω^T · sol - γ||
       (optionally through a constraint+pulse expander E)

    5. Normalise and line-search:
       coeffs = sol · sqrt(N_g) / ||sol||
       dt     = argmin infid(φ + t · coeffs)        # over t ∈ [-t_max, 0]
       φ_new  = φ + dt · coeffs

    6. If fidelity decreased, Gram–Schmidt fallback:
       proj_c = random_direction ⊥ coeffs
       try ±proj_c, keep the side with higher fidelity

The line search interval $[-t_{\max}, 0]$ is the toward-target half-line under the algorithm's sign convention: solving $\omega^\top \cdot \mathrm{sol} = \gamma$ orients coeffs such that negative dt reduces infidelity. The minimiser operates on infid_U_fn, which is always non-negative, so the search is well-defined in both projective=True and projective=False modes. Geope reports fidelity = 1 - infid at the chosen step.

Key functions¤

gammas_and_omegas(free_params) — per-iteration core. Computes the unitary, the geodesic Hamiltonian, the projection $\gamma$, the full Jacobian $\partial U/\partial\phi$, and the per-parameter projections $\omega$. Returns $(\gamma, \omega)$.
linear_comb_projected_coeffs_multigate(ω, γ, E) — least-squares solve, optionally through a constraint expander $E$.
update_linesearch(params, coeffs, piecewise_steps) — golden-section minimisation of $\mathrm{infid}(\phi + t \cdot \mathrm{coeffs})$ over $t \in [-t_{\max}, 0]$.

Constraints¤

Linear equality constraints¤

constraints (or Parameters.constraints) takes a list of vectors $c$ of length $K_{\text{proj}}$ enforcing $c \cdot \phi^{\text{proj}} = 0$, or dicts in control format that are converted into such vectors. Internally, overlapping constraints are merged via an expander matrix $C$ that maps free parameters to the full projected space, and the least-squares solve becomes $\min \|\omega^\top C \tilde c - \gamma\|$ with $c = C \tilde c$.

Pulse-shape constraints¤

pulse_constraints fixes the relative shape of selected parameters across piecewise steps. For each constrained label $k$ the time profile $\phi_k(g)$ is constrained to a one-dimensional subspace:

\[ \phi_k(g) \;=\; \alpha_k\, t_k(g), \qquad \|t_k\| = 1, \quad \alpha_k \in \mathbb{R}. \]

The template $t_k$ is read off the current solution at the moment optimize() is called (or the flat template $\mathbf{1}/\sqrt{N_g}$ if the column is empty), and after every iteration $\phi_k$ is re-projected:

\[ \phi_k \;\leftarrow\; \bigl(\phi_k \cdot t_k\bigr) t_k. \]

Formally, the flat parameter vector $\Phi \in \mathbb{R}^{N_g K_{\text{proj}}}$ is replaced by free parameters $\psi$ via $\Phi = E\,\psi$, where $E$ has $N_g$ identity columns for each unconstrained $k$ and a single template column for each constrained $k$. When combined with a linear-equality expander $C$, the combined expander is $E_{\text{comb}} = (I_{N_g} \otimes C)\,(I_{N_g} \otimes C)^{+}\,E$.

With param_transform, pulse constraints reference parameter indices in $\phi^{\text{exp}}$ rather than projected-basis labels.

Experimental parameters (`param_transform`)¤

GEOPE's native parameters are basis coefficients $\phi^{\text{proj}}_{g,k}$. In practice the experimentally controllable quantities are often different — an amplitude–phase pair driving two basis elements through $\cos/\sin$, a small set of pulse-shape coefficients, a calibration map $\phi^{\text{proj}} = f(\text{voltage}, \text{frequency})$. param_transform lets you optimise directly over those experimental knobs $\phi^{\text{exp}}$:

\[ \phi^{\text{proj}}_{g,\cdot} \;=\; \tau\bigl(\phi^{\text{exp}}_{g,\cdot}\bigr) \;\;\;\text{or}\;\;\; \tau\bigl(\phi^{\text{exp}}_{g,\cdot},\, g\bigr). \]

Contract¤

param_transform must be a JAX-traceable callable. Accepted signatures:

Step-independent: tau(phi) with phi.shape == (n_experimental_params,).
Step-dependent: tau(phi, step_index) with a scalar int32 step index.

The output is a 1-D array whose length is either:

projected_basis.lie_algebra_dim — taken as projected-basis coefficients;
basis.lie_algebra_dim — relevant projected entries extracted automatically via projected_basis.overlap(basis).

Parameters.n_experimental_params sets the input dimension. When param_transform is set, the engine's compute_U_fn is wrapped to apply vmap(τ) over the gate axis, embed the result into the proj+drift slots, broadcast drift coefficients, and delegate to the unitary-product code. The Jacobian is replaced by a split-real-imaginary version (real intermediates in τ would otherwise drop the imaginary part under holomorphic autodiff).

Helper: `make_per_element_transform`¤

For element-wise transforms:

import jax.numpy as jnp

tau = geope.make_per_element_transform([
    jnp.cos,                 # phi[0] → cos(phi[0])
    jnp.sin,                 # phi[1] → sin(phi[1])
    lambda x: 0.5 * x,       # phi[2] → 0.5 * phi[2]
    None,                    # phi[3] passes through
])

Worked example: Rabi rotation in $(A, \varphi)$¤

import numpy as np
import jax.numpy as jnp
import geope

basis = geope.construct_full_pauli_basis(1)

def rabi(phi):                                   # phi = (A, varphi)
    A, varphi = phi[0], phi[1]
    return jnp.array([A * jnp.cos(varphi),       # X coefficient
                      A * jnp.sin(varphi),       # Y coefficient
                      0.0])                      # Z coefficient

theta = np.pi / 3
RX = np.array([[np.cos(theta/2), -1j*np.sin(theta/2)],
               [-1j*np.sin(theta/2),  np.cos(theta/2)]], dtype=complex)

params = geope.Parameters(
    basis=basis, control={1: ['x', 'y', 'z']}, target=RX,
    piecewise_steps=4,
    param_transform=rabi, n_experimental_params=2,
    init_spread=0.3, seed=0,
)
g = geope.Geope(params, history=geope.History())
g.optimize(max_steps=300, precision=1 - 1e-7)
print(float(g.params.fidelity))        # final fidelity (lives on Parameters)
print(g.params.basis_coefficients)     # current params mapped through param_transform
print(g.history.best_fidelity)         # best fidelity over the trajectory

Practical implications¤

Gecko's null-space methods (speed, length, robust) must use parameter_indices, not parameter_labels, when param_transform is set — labels no longer correspond to optimised parameters. Gecko raises ValueError otherwise. (Gecko supports experimental parameters in every construction mode: when reusing a Geope the engine is already wrapped; when built from params it re-wraps a fresh engine.)
Internally param_transform mode uses float64; basis-coefficient mode uses complex128. Tolerances and bounds you supply should match.

Phase-sensitive vs projective¤

The two fidelities differ in how the trace is taken:

\[ F_{\text{proj}}(U, U_T) = \frac{|\mathrm{Tr}(U_T^\dagger U)|}{d}, \qquad F_{\text{full}}(U, U_T) = \frac{\mathrm{Re}\,\mathrm{Tr}(U_T^\dagger U)}{d}. \]

$F_{\text{proj}} \in [0,1]$ is invariant under $U \mapsto e^{i\theta}U$ (the global phase is unobservable). $F_{\text{full}} \in [-1,1]$ is not. Use projective=False only when the absolute phase matters — for example, when the gate is a sub-block of a larger coherent unitary, or when stitching multiple gates whose relative phase enters the composite fidelity.

Two pathologies to keep in mind for phase-sensitive mode:

Traceless targets (Hadamard, single-qubit $X/Y/Z$, etc.) make the gradient of $F_{\text{full}}$ vanish at $U = I$ in every direction; a random init near identity has no descent direction. Use larger init_spread or non-zero init_values.
Stopping criterion. precision = 0.9999999 is meaningful for $F_{\text{proj}}$. For $F_{\text{full}}$ the same threshold is valid near the optimum (both fidelities agree as $U \to U_T$), but the optimiser may transit negative-fidelity regions on its way — that's geometry, not a bug.

Null-space optimisation (`Gecko`)¤

After the main GEOPE loop has converged, the null space of the Jacobian $\omega$ represents directions in parameter space that don't change the unitary to first order. Stepping along these lets you optimise secondary objectives while preserving fidelity.

These passes live on a separate optimiser, Gecko, which post-processes a solution. A Gecko is constructed from a Parameters object — the same object a Geope uses:

Gecko(p) — operate on whatever solution p.parameters holds.
Gecko(g.params) — post-process a Geope result straight after g.optimize(...). Because the optimisation functions are cached on params, this reuses Geope's already-compiled functions rather than recompiling.

The solution does not have to come from Geope. Gecko operates on the current params.parameters — that array can be a Geope result, but it can equally be a solution found by any other method (a different optimiser, an analytic/hand-crafted pulse, an imported result, …). Just put the parameters into a Parameters object describing the same system (basis, projected_basis/drift_basis, target, piecewise_steps, and any param_transform) and call Gecko(p); it refines the imported solution while preserving its fidelity. (When params has never been evaluated, Gecko computes the baseline fidelity itself on construction.)

When you pass a Geope's params, the Parameters object is shared with that Geope, so a pass with piecewise_steps_multiplier > 1 advances the shared state forward (params.parameters and params.piecewise_steps move to the new count together).

Available objectives (methods on `Gecko`)¤

Method	Cost minimised	Purpose
`smooth(...)`	$\sum_g \\|\phi_{g+1} - \phi_g\\|^2$	reduce variation across segments
`smooth_frequency(...)`	$\sum_{m \ge 1, k}	\widehat{\phi_k}(m)
`filter_frequency(filter_fn, ...)`	$\\|\widehat\phi - \mathcal{F}(\widehat\phi)\\|^2$	drive $\phi$ toward a user-defined filtered version (= $L^2$ distance by Parseval)
`speed(parameter_*, ...)`	$\max_{g, k \in P}	\phi_{g,k}
`length(parameter_*, ...)`	$\sum_g \sqrt{\sum_{k \in P}\phi_{g,k}^2 + \\|d_g\\|^2}$	reduce total pulse length (drift contribution included)
`robust(parameter_*, delta, num_samples, ...)`	$1 - \min_{\delta \in [-\Delta,+\Delta]^{	P
`bound(bounds, method, ...)`	$\max(\phi - u_b, l_b - \phi)$	enforce a box constraint via `'projected_gradient'` / `'pg'` or `'mid_point'` / `'mp'`

Each returns (success, iters). Pass piecewise_steps_multiplier > 1 to subdivide existing segments before the pass (linear interpolation), giving more null-space degrees of freedom.

Null-space algorithm: `Gecko._null_space_optimisation()`¤

1. Optionally subdivide piecewise steps (piecewise_steps_multiplier)
2. Build the combined expander (pulse × linear-equality)
3. For each iteration:
   a. Compute Jacobian projections ω
   b. SVD of ω → null-space basis N (right-singular vectors below the rank)
   c. Compute the cost gradient ∇C(φ)
   d. Project the negative gradient onto the null space:
        x = lstsq(N, -∇C)
   e. Step: φ ← φ + rate · N·x / ||x||
   f. Enforce pulse templates if applicable
   g. Recompute fidelity (preserved to first order)

Parameter spaces and index mappings¤

The codebase uses three basis spaces with boolean masks mapping between them:

full_basis  (dim K_full)         — all generators
  ├── projected_indices          — shape (K_full,)
  ├── drift_indices              — shape (K_full,)
  └── proj_drift_indices         — projected | drift

proj_drift_basis  (dim K_pd)     — only projected + drift elements
  ├── proj_indices_projdrift_basis    — shape (K_pd,)
  └── drift_indices_projdrift_basis   — shape (K_pd,)

projected_basis  (dim K_proj)    — only the controllable elements

Parameters are stored in full-basis space $(N_g, K_{\text{full}})$. JIT functions operate on the proj+drift subspace $(N_g, K_{\text{pd}})$. With param_transform, the engine indices are overridden so the optimisation runs uniformly on $\phi^{\text{exp}} \in \mathbb{R}^{N_g \times n_{\text{exp}}}$.

Usage¤

import numpy as np
import geope

# Bases
full    = geope.construct_full_pauli_basis(3)
control = {1: ['x', 'z'], 2: ['x', 'z'], 3: ['x', 'z']}
drift   = {(1, 2): ['zz'], (2, 3): ['zz'], (1, 3): ['zz']}

# Target: Toffoli
target = geope.multicontrol_unitary(np.array([[0, 1], [1, 0]]), 2)

# Bundle everything in a Parameters object
params = geope.Parameters(
    basis=full,
    control=control,
    drift=drift,
    drift_values={(1, 2): {"zz": 1.0},
                  (2, 3): {"zz": 1.0},
                  (1, 3): {"zz": 1.0}},
    target=target,
    piecewise_steps=20,
    seed=0,
)

# Run — updates params (parameters/fidelity) in place; returns the same Parameters
g = geope.Geope(params, history=geope.History())
result = g.optimize(max_steps=1000, precision=0.9999)
print(float(result.fidelity))        # final fidelity (lives on Parameters)
print(result.to_dict())              # current solution as a control dict
print(g.history.best_fidelity)       # best over the trajectory

# Null-space passes — fidelity preserved — live on Gecko, which
# shares the converged optimiser's Parameters (and its cached functions).
gk = geope.Gecko(g.params)
gk.smooth(piecewise_steps_multiplier=2, smoothing_rate=0.05, diff_tol=1e-3)
gk.smooth_frequency(smoothing_rate=0.05, diff_tol=1e-3)
gk.bound({"x": (-1, 1), "z": (-1, 1)}, method='projected_gradient')
gk.robust(parameter_labels=["XII", "IXI", "IIX"], delta=0.01)
gk.speed(parameter_labels=["XII", "IXI", "IIX"])
gk.length()

Refining a solution from another method¤

Gecko does not require the solution to have been produced by Geope. Drop any fidelity-achieving solution — from a different optimiser, an analytic construction, or an imported result — into a Parameters describing the same system, then build a Gecko directly from it:

# `phi` is a (piecewise_steps, K_full) parameter array obtained elsewhere.
params = geope.Parameters(
    basis=full,
    control=control,
    drift=drift,
    target=target,
    piecewise_steps=phi.shape[0],
    init_values=phi,          # the externally-found solution
)

gk = geope.Gecko(params)          # no Geope needed
gk.smooth(piecewise_steps_multiplier=2, smoothing_rate=0.05, diff_tol=1e-3)
print(float(gk.params.fidelity))  # baseline computed on construction, preserved by the pass

Building a `Parameters` from pre-built bases¤

If you've already constructed Basis objects (e.g. via construct_restricted_pauli_basis) and don't want to re-express them as control / drift dicts, pass them directly via the projected_basis and drift_basis kwargs:

projected = geope.construct_restricted_pauli_basis(3, ['x', 'z'])
drift_b   = geope.construct_restricted_pauli_basis(3, ['zz'])

params = geope.Parameters(
    basis=full,
    projected_basis=projected,
    drift_basis=drift_b,
    drift_values=drift_b.generate_parameter_list({
        (1, 2): {"zz": 1.0},
        (2, 3): {"zz": 1.0},
        (1, 3): {"zz": 1.0},
    }),
    target=target,
    piecewise_steps=20,
    seed=0,
)

This is the escape hatch for cases where the projected subset can't be expressed as a control dict. projected_basis and control are mutually exclusive; same for drift_basis and drift.

Function	Description
`construct_full_pauli_basis(n)`	all \(4^n - 1\) non-identity Pauli strings
`construct_two_body_pauli_basis(n)`	1-body and 2-body terms only
`construct_Heisenberg_pauli_basis(n)`	1-body + same-type 2-body (XX, YY, ZZ)
`construct_restricted_pauli_basis(n, restriction)`	custom restriction (list or dict)
`construct_full_spin_boson_basis(n_spins, n_bosons, truncation)`	spin-boson hybrid
`construct_restricted_spin_boson_basis(...)`	restricted spin-boson
`filter_basis_by_control(basis, control)`	filter an existing `Basis` by a control dict (handy when \(d \neq 2^n\))

Method	Cost minimised	Purpose
`smooth(...)`	\(\sum_g \\|\phi_{g+1} - \phi_g\\|^2\)	reduce variation across segments
`smooth_frequency(...)`	$\sum_{m \ge 1, k}	\widehat{\phi_k}(m)
`filter_frequency(filter_fn, ...)`	\(\\|\widehat\phi - \mathcal{F}(\widehat\phi)\\|^2\)	drive \(\phi\) toward a user-defined filtered version (= \(L^2\) distance by Parseval)
`speed(parameter_*, ...)`	$\max_{g, k \in P}	\phi_{g,k}
`length(parameter_*, ...)`	\(\sum_g \sqrt{\sum_{k \in P}\phi_{g,k}^2 + \\|d_g\\|^2}\)	reduce total pulse length (drift contribution included)
`robust(parameter_*, delta, num_samples, ...)`	$1 - \min_{\delta \in [-\Delta,+\Delta]^{	P
`bound(bounds, method, ...)`	\(\max(\phi - u_b, l_b - \phi)\)	enforce a box constraint via `'projected_gradient'` / `'pg'` or `'mid_point'` / `'mp'`

GEOPE – Geodesic Pulse Engineering¤

Overview¤

Class hierarchy¤

Lie group classes (lie.py)¤

Basis¤

Hamiltonian¤

Unitary¤

Basis construction utilities (utils.py)¤

Restriction formats¤

Drift parameter values¤

Linear equality constraints (global controls)¤

Pulse-shape constraints¤

The Parameters object¤

Metadata, functions, and optimiser¤

Algebraic metadata (cached on Parameters)¤

Optimisation functions (engine.py, cached on Parameters)¤

Geope (geope.py)¤

History (history.py)¤

Core algorithm: optimize()¤

Key functions¤

Constraints¤

Linear equality constraints¤

Pulse-shape constraints¤

Experimental parameters (param_transform)¤

Contract¤

Helper: make_per_element_transform¤

Worked example: Rabi rotation in \((A, \varphi)\)¤

Practical implications¤

Phase-sensitive vs projective¤

Null-space optimisation (Gecko)¤

Available objectives (methods on Gecko)¤

Null-space algorithm: Gecko._null_space_optimisation()¤

Parameter spaces and index mappings¤

Usage¤

Refining a solution from another method¤

Building a Parameters from pre-built bases¤

Lie group classes (`lie.py`)¤

`Basis`¤

`Hamiltonian`¤

`Unitary`¤

Basis construction utilities (`utils.py`)¤

The `Parameters` object¤

Algebraic metadata (cached on `Parameters`)¤

Optimisation functions (`engine.py`, cached on `Parameters`)¤

`Geope` (`geope.py`)¤

`History` (`history.py`)¤

Core algorithm: `optimize()`¤

Experimental parameters (`param_transform`)¤

Helper: `make_per_element_transform`¤

Null-space optimisation (`Gecko`)¤

Available objectives (methods on `Gecko`)¤

Null-space algorithm: `Gecko._null_space_optimisation()`¤

Building a `Parameters` from pre-built bases¤