Getting Started with GEOPE

This notebook walks through a minimal example of using GEOPE (Geodesic Pulse Engineering) to synthesise a quantum gate.

We will: 1. Define a target unitary (the Hadamard gate on 1 qubit). 2. Build the Lie algebra basis. 3. Run the GEOPE optimiser. 4. Check the result.

import numpy as np
import matplotlib.pyplot as plt

from geope import (
    Geope,
    Parameters,
    History,
    Hamiltonian,
    construct_full_pauli_basis,
    construct_restricted_pauli_basis,
    multimatmul,
    fidelity
)

1. Define the target unitary

We target the Hadamard gate, a single-qubit gate that maps \(|0\rangle \to |+\rangle\) and \(|1\rangle \to |-\rangle\):

\[H = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\]

target = np.array([[1, 1], [1, -1]], dtype=complex) / np.sqrt(2)
print("Target unitary (Hadamard):")
print(target)

2. Build the Pauli bases

For a single qubit the full Pauli basis is \(\{\sigma_x, \sigma_y, \sigma_z\}\) (excluding identity). GEOPE uses this as the Lie algebra generators from which unitaries are constructed via the matrix exponential, a single time step expression of the unitary is:

\[U(\boldsymbol{\phi}) = \exp\!\bigl(i \sum_k \phi_k \sigma_k\bigr)\]

A restricted Pauli basis is constructed by the construct_restricted_pauli_basis(n, restriction) function. This function can takes the number of spins (n) and either a list of interaction strings (e.g. ['x', 'y'] or for 2-qubit terms ['xx', 'yy']) or a dictionary mapping qubit labels to interaction strings (e.g. {1: ['x', 'y'], 2: ['x', 'y'], (1,2): ['zz']} if there were two qubits with an Ising interaction term). The labels 'x', 'y', and 'z' refer to Pauli spin matrices \(\sigma^x\), \(\sigma^y\), and \(\sigma^z\).

The restricted Pauli basis is used to construct the control basis (projected_basis) and the drift basis (drift_basis).

In this case, we have control of local x and y fields, and there is a drift z term, giving Hamiltonian:

\[ H = \sigma^z + h_x(t) \sigma^x + h_y(t) \sigma^y.\]

The drift parameter is defined to be 1. This is also the default value for drift parameters if they are not specified.

The initial parameters are set to be 1. The default value are random parameters – see the documentation for details.

nqubits = 1
basis = construct_full_pauli_basis(nqubits)
projected_basis = construct_restricted_pauli_basis(nqubits, {1: ['x', 'y']})
drift_basis = construct_restricted_pauli_basis(nqubits, {1: ['z']})
init_parameters = [0.5, -0.5]
drift_parameters = [1.]

print(f"Number of qubits:       {nqubits}")
print(f"Hilbert space dim:      {basis.dim}")
print(f"Lie algebra dim:        {basis.lie_algebra_dim}")
print(f"Basis labels:           {basis.labels}")
print(f"Projected basis labels: {projected_basis.labels}")
print(f"Drift basis labels:     {drift_basis.labels}")

3. Create the parameters and optimiser

Parameters bundles the problem definition — the bases, target, initial/drift values, number of piecewise steps and seed. Pass it to Geope, which builds the JAX functions (cached on Parameters, compiled lazily on first optimize()) (unitary computation, Jacobian, geodesic Hamiltonian, fidelity) and wraps them in an optimisation loop.

Geope itself takes only verbose and an optional History. The run-control knobs (max_step_size, precision, gram_schmidt_step_size) are arguments of optimize(), not the constructor.

An optional History logger records the per-step trajectory (parameters, fidelities, …). Without it, only the live/final state on Parameters is kept.

The number of piecewise steps is set on Parameters (piecewise_steps).

params = Parameters(
    basis=basis,
    projected_basis=projected_basis,
    drift_basis=drift_basis,
    init_values=init_parameters,
    drift_values=drift_parameters,
    target=target,
    piecewise_steps=3,
    seed=42,
)

opt = Geope(
    params,
    verbose=True,
    history=History(),
)

4. Run the optimisation

optimize(max_steps=..., max_step_size=..., precision=...) runs the geodesic update loop. Each step: 1. Computes the geodesic Hamiltonian pointing from the current unitary towards the target. 2. Projects it onto the controllable basis to find the optimal parameter update direction. 3. Performs a line search along that direction.

It returns the bound Parameters object, whose parameters / fidelity hold the final result; convergence is reached when params.fidelity >= precision. With a History attached, the full trajectory is available on opt.history.

opt.optimize(max_steps=200, max_step_size=0.9, precision=0.9999999)

converged = opt.params.fidelity >= opt.precision
print(f"\nConverged:        {converged}")
print(f"Steps:            {opt.history.steps[-1]}")
print(f"Fidelity:         {opt.params.fidelity}")
print(f"Final parameters: {opt.params.parameters}")

5. Check the result

Convergence plot

The fidelity should climb rapidly towards 1.

fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(opt.history.steps, opt.history.fidelities, "o-", markersize=3)
ax.set_xlabel("Step")
ax.set_ylabel("Fidelity")
ax.set_title("GEOPE convergence \u2014 Hadamard gate")
ax.axhline(opt.precision, color="grey", linestyle="--", linewidth=0.8, label=f"Target = {opt.precision}")
ax.legend()
plt.tight_layout()
plt.show()

Recovered parameters

The optimised Lie-algebra coefficients \(\boldsymbol{\phi}\) define the Hamiltonian whose exponential gives the target gate.

final_params = opt.params.parameters
print("Optimised parameters:")
for t, param_step in enumerate(final_params):
    print(f"--\nTime step {t}")
    for label, phi in zip(basis.labels, param_step):
        print(f"\t{label}: {phi:.6f}")

Verify the synthesised unitary

Reconstruct the unitary from the optimised parameters and compare it to the target. We find the target unitary up to global phase.

Us = [Hamiltonian(basis, param_step).unitary.matrix for param_step in final_params]
unitary_matrix = multimatmul(Us[::-1]) # Note the order is reversed since the full unitary is U_2 U_1 U_0

print("Synthesised unitary:")
print(unitary_matrix)

print("\nTarget unitary:")
print(target)

fid = fidelity(unitary_matrix, target)
print(f"\nFidelity: {fid}")