Variational Quantum Algorithms with tket

Variational Algorithms in tket

The Variational Quantum Eigensolver (VQE) is a hybrid quantum-classical algorithm that finds the ground state energy of a molecular Hamiltonian. A quantum processor evaluates a parameterized circuit (the ansatz) and measures expectation values, while a classical optimizer adjusts the parameters to minimize the energy. This feedback loop runs until the energy converges.

VQE belongs to the broader family of variational quantum algorithms, which also includes QAOA for combinatorial optimization and variational quantum classifiers for machine learning. The core workflow is the same in every case: prepare a parameterized quantum state, measure a cost function, update parameters classically, repeat.

tket supports parameterized circuits through symbolic parameters built on top of sympy. You define a circuit with free parameters, compile and optimize it once, then substitute concrete values at execution time. This compile-once pattern is central to writing efficient VQE code with tket, and we cover it in detail below.

tket Symbolic Parameters in Depth

Pytket uses sympy Symbol objects to represent free parameters in quantum gates. Any gate that accepts an angle (such as Ry, Rx, Rz, XXPhase, ZZPhase) can take a Symbol instead of a float. The convention in tket is that angles are measured in half-turns (multiples of pi), so Ry(0.5, 0) applies a rotation of 0.5 * pi = pi/2 radians.

from pytket import Circuit
from sympy import Symbol

theta = Symbol("theta")

circ = Circuit(2)
circ.Ry(theta, 0)   # Symbolic Y-rotation: angle = theta * pi radians
circ.CX(0, 1)
print(circ.free_symbols())  # {theta}

A circuit can hold multiple symbolic parameters:

from pytket import Circuit
from sympy import Symbol

theta1 = Symbol("theta1")
theta2 = Symbol("theta2")
theta3 = Symbol("theta3")

circ = Circuit(2)
circ.Ry(theta1, 0)
circ.Ry(theta2, 1)
circ.CX(0, 1)
circ.Rz(theta3, 0)

print(circ.free_symbols())  # {theta1, theta2, theta3}

To evaluate the circuit, you call symbol_substitution with a dictionary mapping each Symbol to a float value. After substitution, every symbolic gate becomes a concrete gate and the circuit is ready for execution:

from copy import deepcopy

# Always copy before substitution to preserve the original symbolic circuit
bound = circ.copy()
bound.symbol_substitution({theta1: 0.3, theta2: -0.1, theta3: 0.7})
print(bound.free_symbols())  # set() -- no free symbols remain

Two important details about symbol_substitution:

1. It modifies the circuit in place. If you need the original symbolic circuit for the next optimizer iteration, copy it first. Forgetting this is a common source of bugs.
2. You can do partial substitution. If you substitute only some symbols, the circuit remains symbolic with the remaining free parameters.

The Hydrogen Molecule Problem

We estimate the ground state energy of H2 at the equilibrium bond length of 0.735 angstroms. At this geometry, the electronic Hamiltonian in the STO-3G basis set reduces (after qubit mapping and symmetry reduction) to a 2-qubit operator expressed as a sum of Pauli terms:

H = -1.0524 * II + 0.3979 * IZ - 0.3979 * ZI - 0.0112 * ZZ + 0.1809 * XX

The coefficients come from standard quantum chemistry calculations (for example, via PySCF and OpenFermion). We hard-code them here to focus on the tket workflow. The exact ground state energy for this Hamiltonian is approximately -1.137 Hartree.

Note: these Pauli terms use the parity mapping for fermion-to-qubit encoding. If you use the Jordan-Wigner or Bravyi-Kitaev mappings instead, you get different Pauli strings with different coefficients. The algorithm works the same way regardless of which mapping you choose.

Building the Parameterized Ansatz

Minimal Single-Parameter Ansatz

A minimal ansatz for H2 uses a single variational parameter controlling an entangling rotation:

from pytket import Circuit
from sympy import Symbol

theta = Symbol("theta")

ansatz = Circuit(2)
ansatz.Ry(theta, 0)    # Parameterized Y-rotation on qubit 0
ansatz.CX(0, 1)        # Entangle the two qubits

This creates a circuit where theta is a free symbolic parameter. The state produced is cos(theta*pi/2)|00> + sin(theta*pi/2)|11>, which can represent any state in the subspace spanned by |00> and |11>. Since the H2 ground state lies primarily in this subspace, a single parameter suffices.

UCCSD Ansatz for H2

For a more physically motivated ansatz, consider the Unitary Coupled Cluster Singles and Doubles (UCCSD) approach. For H2 with 2 qubits, the UCCSD ansatz has a single parameter and implements the unitary:

U(theta) = exp(i * theta * (Y0 X1 - X0 Y1) / 2)

This operator generates excitations between the reference state |01> (one electron in each spin-orbital) and the doubly-excited state |10>. The decomposition into elementary gates is:

from pytket import Circuit
from sympy import Symbol

theta = Symbol("theta")

uccsd = Circuit(2)
# Prepare the Hartree-Fock reference state |01>
uccsd.X(0)

# Decomposition of exp(i*theta*(Y0 X1 - X0 Y1)/2)
# This is a standard decomposition using CNOT and Ry gates
uccsd.CX(0, 1)
uccsd.Ry(theta, 0)    # theta in half-turns (multiply by pi for radians)
uccsd.CX(0, 1)

The UCCSD ansatz is chemically motivated: it preserves particle number and spin symmetry by construction. This means the optimizer searches only over physically meaningful states, which generally improves convergence compared to arbitrary parameterized circuits.

Hardware-Efficient Ansatz

An alternative is the hardware-efficient ansatz, which uses layers of single-qubit rotations and entangling gates without any chemistry-specific structure:

from pytket import Circuit
from sympy import symbols

# Two layers, 2 qubits: 8 total parameters
params = symbols("p0:8")

hea = Circuit(2)

# Layer 1: single-qubit rotations
hea.Ry(params[0], 0)
hea.Ry(params[1], 1)
hea.Rz(params[2], 0)
hea.Rz(params[3], 1)

# Entangling layer
hea.CX(0, 1)

# Layer 2: single-qubit rotations
hea.Ry(params[4], 0)
hea.Ry(params[5], 1)
hea.Rz(params[6], 0)
hea.Rz(params[7], 1)

Comparing the two ansatz styles:

Property	UCCSD	Hardware-Efficient
Parameters (H2)	1	4 to 8
Physical motivation	Yes, preserves symmetries	No, general unitary
Circuit depth	Moderate, may need transpilation	Shallow, hardware-native gates
Expressibility	Limited to excitation subspace	High, but may explore unphysical states
Convergence	Typically fast (fewer parameters)	Can be slow, barren plateaus possible
Best for	Small molecules, proof of concept	Large systems, specific hardware topologies

For H2, UCCSD is the better choice because a single parameter captures the relevant physics. For larger molecules or when running on specific hardware, hardware-efficient ansatze can be more practical.

Measuring Expectation Values

Each Pauli term in the Hamiltonian must be measured in a separate circuit (or grouped into compatible sets). For a term like XX, we need to rotate both qubits into the X basis before measurement. For ZZ, no rotation is needed since measurement is natively in the Z basis.

Two Pauli strings are compatible (measurable in the same circuit) only if they commute qubit-wise. For our Hamiltonian, ZZ and XX do not commute on qubits 0 and 1 (Z and X anticommute), so they require separate measurement circuits. The terms II, IZ, ZI, and ZZ are all diagonal in the computational basis and can be extracted from a single Z-basis measurement.

from pytket import Circuit
from copy import deepcopy


def measure_pauli_string(ansatz, pauli_string, backend, n_shots=4096):
    """
    Measure the expectation value of a Pauli string.

    Args:
        ansatz: A concrete (non-symbolic) pytket Circuit.
        pauli_string: dict mapping qubit index to 'I', 'X', 'Y', 'Z'.
        backend: A pytket Backend instance.
        n_shots: Number of measurement shots.

    Returns:
        Float expectation value in [-1, 1].
    """
    circ = ansatz.copy()

    # Rotate into the measurement basis for each qubit
    for qubit_idx, pauli in pauli_string.items():
        if pauli == "X":
            circ.H(qubit_idx)           # H rotates Z eigenstates to X eigenstates
        elif pauli == "Y":
            circ.Sdg(qubit_idx)         # S-dagger followed by H rotates Z to Y
            circ.H(qubit_idx)
        # Z and I require no rotation

    circ.measure_all()

    compiled = backend.get_compiled_circuit(circ)
    handle = backend.process_circuit(compiled, n_shots=n_shots)
    result = backend.get_result(handle)
    counts = result.get_counts()

    # Compute expectation value from measurement statistics
    # Each bitstring contributes +1 (even parity) or -1 (odd parity)
    total = sum(counts.values())
    expectation = 0.0
    for bitstring, count in counts.items():
        relevant_bits = [
            bitstring[i] for i, p in pauli_string.items() if p != "I"
        ]
        parity = (-1) ** sum(relevant_bits)
        expectation += parity * count / total

    return expectation

The parity calculation works because Pauli Z has eigenvalues +1 (for |0>) and -1 (for |1>). For a product of Pauli Z operators, the eigenvalue is the product of individual eigenvalues, which equals (-1) raised to the number of qubits measured as |1>. The basis rotations before measurement transform X and Y measurements into equivalent Z measurements.

Expectation Values from Statevector Simulation

During algorithm development, shot-based sampling introduces noise that makes debugging difficult. You can compute exact expectation values using statevector simulation instead. This removes all statistical fluctuation and runs much faster.

import numpy as np
from pytket import Circuit, Qubit
from pytket.extensions.qiskit import AerBackend


def expectation_from_statevector(circuit, pauli_string):
    """
    Compute the exact expectation value of a Pauli string using
    statevector simulation. No shots, no statistical noise.

    Args:
        circuit: A concrete (non-symbolic) pytket Circuit with no measurements.
        pauli_string: dict mapping qubit index to 'I', 'X', 'Y', 'Z'.

    Returns:
        Float expectation value (exact).
    """
    n_qubits = circuit.n_qubits

    # Build the Pauli operator as a matrix
    pauli_matrices = {
        "I": np.eye(2),
        "X": np.array([[0, 1], [1, 0]]),
        "Y": np.array([[0, -1j], [1j, 0]]),
        "Z": np.array([[1, 0], [0, -1]]),
    }

    # Construct the full operator via tensor product
    # tket uses big-endian qubit ordering for statevectors
    op = np.array([[1.0]])
    for q in range(n_qubits):
        p = pauli_string.get(q, "I")
        op = np.kron(op, pauli_matrices[p])

    # Get the statevector from the simulator
    backend = AerBackend()
    circ_copy = circuit.copy()
    compiled = backend.get_compiled_circuit(circ_copy)
    handle = backend.process_circuit(compiled)
    result = backend.get_result(handle)
    statevector = result.get_state()

    # Expectation value = <psi|O|psi>
    expectation = np.real(statevector.conj() @ op @ statevector)
    return float(expectation)

Use this function during development and testing. Switch to shot-based measurement when you are ready to run on hardware or need to study the effects of finite sampling.

The Compile-Once Pattern for VQE

A VQE optimization loop calls the energy function hundreds of times. Each call builds a circuit, compiles it, and executes it. Compilation (routing to hardware topology, gate optimization, rebasing to native gate sets) is computationally expensive. The naive approach compiles inside the energy function on every single iteration:

# BAD: compiles on every iteration
def energy_naive(params):
    bound = ansatz.copy()
    bound.symbol_substitution({theta: params[0]})
    compiled = backend.get_compiled_circuit(bound)  # expensive!
    # ... measure and return energy

The correct approach compiles the symbolic circuit once before the optimization loop starts. Tket’s compilation passes preserve symbolic parameters, so you can compile first, then substitute later:

from pytket import Circuit
from pytket.extensions.qiskit import AerBackend
from sympy import Symbol

theta = Symbol("theta")

ansatz = Circuit(2)
ansatz.Ry(theta, 0)
ansatz.CX(0, 1)

backend = AerBackend()

# Compile once with symbols still present
compiled_ansatz = backend.get_compiled_circuit(ansatz)
# compiled_ansatz still has theta as a free symbol

# GOOD: substitute into the already-compiled circuit each iteration
def energy_fast(params):
    bound = compiled_ansatz.copy()
    bound.symbol_substitution({theta: params[0]})
    # No compilation needed, just execute
    # ... measure and return energy

This pattern gives identical results but eliminates repeated compilation overhead. For backends with complex routing requirements (like superconducting hardware with limited connectivity), the speedup is substantial.

Note that get_compiled_circuit applies backend-specific compilation including rebasing to the native gate set and routing to the device topology. Symbolic parameters flow through these transformations because tket treats them algebraically. The compiled circuit may have different gates (for example, Ry(theta) might become Rz and SX gates with symbolic expressions), but the substitution step still works correctly.

The Full VQE Loop

Now we connect the ansatz, the Hamiltonian, and a classical optimizer:

from pytket import Circuit
from pytket.extensions.qiskit import AerBackend
from scipy.optimize import minimize
from sympy import Symbol
from copy import deepcopy

# Define symbolic parameter
theta = Symbol("theta")

# Build and compile ansatz once
ansatz = Circuit(2)
ansatz.Ry(theta, 0)
ansatz.CX(0, 1)

backend = AerBackend()
compiled_ansatz = backend.get_compiled_circuit(ansatz)

# H2 Hamiltonian coefficients and Pauli terms
hamiltonian = [
    (-1.0524, {}),                           # II (constant term)
    ( 0.3979, {1: "Z"}),                     # IZ
    (-0.3979, {0: "Z"}),                     # ZI
    (-0.0112, {0: "Z", 1: "Z"}),             # ZZ
    ( 0.1809, {0: "X", 1: "X"}),             # XX
]

def energy(params):
    """Evaluate <psi(theta)|H|psi(theta)> via shot-based sampling."""
    # Substitute parameter into compiled circuit
    bound_circuit = compiled_ansatz.copy()
    bound_circuit.symbol_substitution({theta: params[0]})

    total_energy = 0.0
    for coeff, pauli_term in hamiltonian:
        if not pauli_term:
            total_energy += coeff
        else:
            exp_val = measure_pauli_string(
                bound_circuit, pauli_term, backend, n_shots=4096
            )
            total_energy += coeff * exp_val

    return total_energy

# Run the classical optimizer
result = minimize(
    energy,
    x0=[0.5],              # Initial guess for theta (in half-turns)
    method="COBYLA",        # Gradient-free optimizer
    options={"maxiter": 50, "rhobeg": 0.5},
)

print(f"Optimal theta: {result.x[0]:.4f}")
print(f"Ground state energy: {result.fun:.4f} Ha")
# Expected: approximately -1.137 Hartree

Convergence Monitoring

To understand how the optimizer progresses, add a callback that records the energy at each iteration:

from pytket import Circuit
from pytket.extensions.qiskit import AerBackend
from scipy.optimize import minimize
from sympy import Symbol

theta = Symbol("theta")

ansatz = Circuit(2)
ansatz.Ry(theta, 0)
ansatz.CX(0, 1)

backend = AerBackend()
compiled_ansatz = backend.get_compiled_circuit(ansatz)

hamiltonian = [
    (-1.0524, {}),
    ( 0.3979, {1: "Z"}),
    (-0.3979, {0: "Z"}),
    (-0.0112, {0: "Z", 1: "Z"}),
    ( 0.1809, {0: "X", 1: "X"}),
]

# Track convergence history
energy_history = []

def energy_with_logging(params):
    """Energy function that records each evaluation."""
    bound_circuit = compiled_ansatz.copy()
    bound_circuit.symbol_substitution({theta: params[0]})

    total_energy = 0.0
    for coeff, pauli_term in hamiltonian:
        if not pauli_term:
            total_energy += coeff
        else:
            exp_val = measure_pauli_string(
                bound_circuit, pauli_term, backend, n_shots=4096
            )
            total_energy += coeff * exp_val

    energy_history.append(total_energy)
    print(f"Iteration {len(energy_history):3d}: "
          f"theta = {params[0]:+.4f}, energy = {total_energy:.6f} Ha")
    return total_energy

result = minimize(
    energy_with_logging,
    x0=[0.5],
    method="COBYLA",
    options={"maxiter": 50, "rhobeg": 0.5},
)

Expected output (values will vary slightly due to shot noise):

Iteration   1: theta = +0.5000, energy = -1.052437 Ha
Iteration   2: theta = +1.0000, energy = -1.081223 Ha
Iteration   3: theta = +0.7500, energy = -1.098742 Ha
...
Iteration  18: theta = +0.5912, energy = -1.136284 Ha
Iteration  19: theta = +0.5918, energy = -1.137012 Ha
Iteration  20: theta = +0.5915, energy = -1.136891 Ha

The energy decreases monotonically (on average) toward the true ground state energy of -1.137 Ha. Shot noise causes small fluctuations around the converged value.

To visualize the convergence:

import matplotlib.pyplot as plt

plt.figure(figsize=(8, 5))
plt.plot(energy_history, "o-", markersize=3)
plt.axhline(y=-1.137, color="r", linestyle="--", label="Exact ground state")
plt.xlabel("Iteration")
plt.ylabel("Energy (Ha)")
plt.title("VQE Convergence for H2")
plt.legend()
plt.tight_layout()
plt.savefig("vqe_convergence.png", dpi=150)

Why COBYLA (and When to Use Gradients Instead)

Gradient-based optimizers (L-BFGS-B, Adam) require computing partial derivatives of the cost function. On quantum hardware, this means running additional circuits for each parameter using the parameter-shift rule or finite differences. For a single-parameter problem this is manageable, but the overhead grows linearly with parameter count.

COBYLA and other gradient-free methods (Nelder-Mead, Powell) evaluate only the cost function itself. This makes them simpler to implement and often more practical on noisy hardware where gradient estimates are unreliable. The tradeoff: gradient-free methods typically need more function evaluations to converge, especially in high-dimensional parameter spaces.

Rules of thumb for optimizer selection:

• 1 to 5 parameters: COBYLA works well, gradient overhead is not justified
• 5 to 20 parameters: parameter-shift gradients with L-BFGS-B often converge faster
• 20+ parameters: consider stochastic gradient methods (SPSA, Adam) or simultaneous perturbation
• Noisy hardware: gradient-free methods are more robust to noise in the cost function

The Parameter-Shift Rule in tket

The parameter-shift rule computes exact analytical gradients of quantum expectation values. For a gate of the form exp(-i * theta * G / 2) where G has eigenvalues +/-1 (true for all standard Pauli rotation gates), the gradient with respect to theta is:

dE/dtheta = [E(theta + pi/2) - E(theta - pi/2)] / 2

This requires exactly two circuit evaluations per parameter, each with the parameter shifted by +/-pi/2. Here is how to implement it manually with tket:

import numpy as np
from pytket import Circuit
from pytket.extensions.qiskit import AerBackend
from scipy.optimize import minimize
from sympy import symbols

# Multi-parameter ansatz
p0, p1, p2, p3 = symbols("p0 p1 p2 p3")
param_symbols = [p0, p1, p2, p3]

ansatz = Circuit(2)
ansatz.Ry(p0, 0)
ansatz.Ry(p1, 1)
ansatz.CX(0, 1)
ansatz.Ry(p2, 0)
ansatz.Rz(p3, 1)

backend = AerBackend()
compiled_ansatz = backend.get_compiled_circuit(ansatz)

hamiltonian = [
    (-1.0524, {}),
    ( 0.3979, {1: "Z"}),
    (-0.3979, {0: "Z"}),
    (-0.0112, {0: "Z", 1: "Z"}),
    ( 0.1809, {0: "X", 1: "X"}),
]


def evaluate_energy(param_values, n_shots=4096):
    """Evaluate energy at a specific parameter vector."""
    bound = compiled_ansatz.copy()
    sub_map = {sym: val for sym, val in zip(param_symbols, param_values)}
    bound.symbol_substitution(sub_map)

    total_energy = 0.0
    for coeff, pauli_term in hamiltonian:
        if not pauli_term:
            total_energy += coeff
        else:
            total_energy += coeff * measure_pauli_string(
                bound, pauli_term, backend, n_shots=n_shots
            )
    return total_energy


def compute_gradient(param_values, n_shots=4096):
    """
    Compute the gradient using the parameter-shift rule.

    For each parameter, evaluate the energy at theta +/- pi/2
    (which is +/- 0.5 in tket's half-turn convention).
    """
    gradient = np.zeros(len(param_values))
    shift = 0.5  # pi/2 in half-turns

    for i in range(len(param_values)):
        # Shift parameter i forward
        params_plus = np.array(param_values, dtype=float)
        params_plus[i] += shift
        e_plus = evaluate_energy(params_plus, n_shots=n_shots)

        # Shift parameter i backward
        params_minus = np.array(param_values, dtype=float)
        params_minus[i] -= shift
        e_minus = evaluate_energy(params_minus, n_shots=n_shots)

        gradient[i] = (e_plus - e_minus) / 2.0

    return gradient


# Gradient-based optimization with L-BFGS-B
result = minimize(
    evaluate_energy,
    x0=[0.1, 0.2, 0.3, 0.1],
    jac=compute_gradient,          # Supply analytical gradient
    method="L-BFGS-B",
    options={"maxiter": 100},
)

print(f"Optimal parameters: {result.x}")
print(f"Ground state energy: {result.fun:.4f} Ha")

The parameter-shift rule costs 2 * N circuit evaluations per gradient computation, where N is the number of parameters. Each circuit evaluation itself requires one run per Pauli term in the Hamiltonian (or per group of compatible terms). For our 5-term H2 Hamiltonian with 4 parameters, each gradient step requires 2 * 4 * 4 = 32 circuit executions (4 non-identity terms, 4 parameters, 2 shifts each). Compare this to COBYLA, which uses only 4 circuit executions per step (one energy evaluation with 4 non-identity terms).

Optimizing Circuits with tket Compilation Passes

tket provides powerful compilation passes that reduce gate count and circuit depth. For VQE circuits, this translates directly into less noise on real hardware.

from pytket import Circuit
from pytket.passes import FullPeepholeOptimise, RemoveRedundancies
from sympy import Symbol

theta = Symbol("theta")

# Build a deliberately unoptimized circuit
circ = Circuit(2)
circ.H(0)
circ.H(0)           # Two consecutive H gates cancel
circ.Ry(theta, 0)
circ.CX(0, 1)
circ.CX(0, 1)       # Two consecutive CX gates cancel

print(f"Before optimization: {circ.n_gates} gates")

# RemoveRedundancies catches simple cancellations
RemoveRedundancies().apply(circ)
print(f"After RemoveRedundancies: {circ.n_gates} gates")

For concrete (non-symbolic) circuits, FullPeepholeOptimise applies aggressive optimization including gate commutation, rotation merging, and KAK decomposition:

from pytket import Circuit
from pytket.passes import FullPeepholeOptimise

concrete = Circuit(2)
concrete.Ry(0.3, 0)
concrete.Rz(0.2, 0)
concrete.Ry(0.1, 0)
concrete.CX(0, 1)
concrete.Rz(0.4, 1)
concrete.CX(0, 1)

print(f"Before: {concrete.n_gates} gates, depth {concrete.depth()}")
FullPeepholeOptimise().apply(concrete)
print(f"After: {concrete.n_gates} gates, depth {concrete.depth()}")

An important caveat: FullPeepholeOptimise has limited effect on symbolic circuits. Many of its optimizations (such as merging consecutive rotations) require knowing the actual angle values. It can still perform structural optimizations like cancelling adjacent inverse gates, but the heavy numerical optimizations only kick in after parameter substitution.

For VQE, the best strategy combines both approaches:

1. Compile the symbolic circuit with get_compiled_circuit (handles routing and rebasing).
2. In the energy function, after substituting parameters, optionally apply FullPeepholeOptimise to the concrete circuit for maximum gate reduction.

The cost of step 2 is small compared to circuit execution on real hardware, so it is usually worthwhile.

tket Backends and the BackendResult API

One of tket’s primary strengths is backend portability. The same VQE code runs on simulators and on hardware from different vendors. You only change the backend object.

Available backends (each requires its own pytket extension package):

Backend	Package	Hardware	Native Gates
AerBackend	pytket-qiskit	Qiskit Aer simulator	U3, CX
IBMQBackend	pytket-qiskit	IBM quantum hardware	SX, Rz, CX, X
QuantinuumBackend	pytket-quantinuum	Quantinuum H-series (trapped ion)	Rz, PhasedX, ZZPhase
ForestBackend	pytket-pyquil	Rigetti QPUs	Rz, RX(pi/2), CZ
BraketBackend	pytket-braket	AWS Braket (IonQ, Rigetti, OQC)	Varies by device

Switching between backends:

from pytket.extensions.qiskit import AerBackend

# Simulator for development
backend = AerBackend()

# To switch to IBM hardware (requires IBM Quantum account):
# from pytket.extensions.qiskit import IBMQBackend
# backend = IBMQBackend("ibm_brisbane")

# To switch to Quantinuum (requires Quantinuum account):
# from pytket.extensions.quantinuum import QuantinuumBackend
# backend = QuantinuumBackend(device_name="H1-1SC")  # H1-1SC is the emulator

# The rest of the VQE code remains identical
compiled_ansatz = backend.get_compiled_circuit(ansatz)

The get_compiled_circuit method handles all backend-specific compilation: rebasing to the native gate set, routing to the device topology, and inserting SWAP gates where needed. This means your ansatz definition stays hardware-agnostic.

Retrieving results from any backend follows the same API:

handle = backend.process_circuit(compiled_circuit, n_shots=4096)
result = backend.get_result(handle)

# Shot-based results
counts = result.get_counts()     # dict mapping tuples to counts
shots = result.get_shots()       # numpy array of shape (n_shots, n_qubits)

# Statevector (only available on simulators, circuit must have no measurements)
# state = result.get_state()     # numpy array of complex amplitudes

The get_counts() method returns a dictionary where keys are tuples of integers (0 or 1) and values are the number of times that bitstring was observed. The get_shots() method returns the raw shot data as a 2D numpy array.

Running on Quantinuum Hardware

Quantinuum’s H-series trapped-ion processors have all-to-all connectivity (any qubit can interact with any other qubit directly) and high-fidelity gates. This means no SWAP overhead for routing. The native gate set includes ZZPhase (a two-qubit entangling gate) instead of CX.

from pytket.extensions.quantinuum import QuantinuumBackend

# Use the H1-1 syntax checker for testing (no credits consumed)
backend = QuantinuumBackend(device_name="H1-1SC")

# Compile the ansatz for Quantinuum's native gate set
# tket automatically decomposes CX into ZZPhase gates
compiled_ansatz = backend.get_compiled_circuit(ansatz)

# Submit the circuit
handle = backend.process_circuit(compiled_ansatz, n_shots=4096)
result = backend.get_result(handle)
counts = result.get_counts()

Key differences when running on Quantinuum:

• CX gates are decomposed into ZZPhase gates automatically by tket’s compilation.
• No routing overhead because the trapped-ion architecture has full connectivity.
• Gate fidelities are typically higher than superconducting qubits (>99.5% for two-qubit gates).
• Circuit execution is slower (milliseconds per gate vs nanoseconds for superconducting).
• Use the syntax checker (H1-1SC) during development to validate circuits without spending credits.

Shot Noise Analysis

With finite measurement shots, the estimated expectation value of any Pauli operator fluctuates around the true value. Understanding this uncertainty is critical for knowing how many shots you need.

For a single Pauli operator P, the measurement outcome is always +1 or -1. The variance of a single measurement is:

Var(P) = 1 - <P>^2

This is at most 1 (when

= 0). The variance of the sample mean over N shots is:

Var(mean) = Var(P) / N <= 1 / N

So the standard deviation (uncertainty) of the expectation value estimate is at most 1 / sqrt(N).

For the full Hamiltonian H = sum(c_i * P_i), the variance of the energy estimate (assuming independent measurements for each term) is:

Var(E) = sum(c_i^2 * Var(<P_i>) / N_i)

where N_i is the number of shots allocated to term i.

Practical shot budgets for H2:

Shots per term	Max uncertainty per term	Energy uncertainty (approx)
1,024	0.031	~0.018 Ha
4,096	0.016	~0.009 Ha
16,384	0.008	~0.004 Ha
1,000,000	0.001	~0.0006 Ha

To achieve chemical accuracy (1.6 mHa = 0.0016 Ha), you need roughly 100,000 to 1,000,000 shots per Pauli term, depending on the coefficients. This is one reason why VQE on current hardware often uses thousands of shots per term as a practical compromise, accepting energy precision of a few mHa.

import numpy as np

def estimate_shots_needed(coefficients, target_precision):
    """
    Estimate the total shot budget for a given energy precision.

    Assumes uniform shot allocation and worst-case variance (Var = 1).
    Optimal allocation would weight shots by |coefficient|.
    """
    # Worst case: Var(<P_i>) = 1 for all terms
    # Var(E) = sum(c_i^2 / N) = sum(c_i^2) / N
    sum_c2 = sum(c ** 2 for c in coefficients)
    # We want sqrt(Var(E)) <= target_precision
    # sum_c2 / N <= target_precision^2
    # N >= sum_c2 / target_precision^2
    n_per_term = int(np.ceil(sum_c2 / target_precision ** 2))
    return n_per_term

# H2 non-identity coefficients
coeffs = [0.3979, 0.3979, 0.0112, 0.1809]
shots = estimate_shots_needed(coeffs, target_precision=0.001)
print(f"Shots per term for 1 mHa precision: {shots:,}")
# Output: roughly 356,000 shots per term

Comparison: tket vs Qiskit for VQE

Both frameworks can implement VQE, but they differ in philosophy and strengths:

Aspect	tket	Qiskit
Circuit optimization	Aggressive peephole optimization, KAK decomposition	Preset pass managers (optimization levels 0-3)
Backend support	IBM, Quantinuum, Rigetti, IonQ, OQC via extensions	Primarily IBM, other providers via Qiskit extensions
Symbolic parameters	sympy Symbols, compile-once pattern	Qiskit ParameterVector, bind_parameters
Built-in VQE	Manual implementation (more control)	EstimatorV2 primitive, built-in VQE class
Ansatz library	Manual construction	EfficientSU2, UCCSD via qiskit-nature
Gate reduction	Typically stronger for multi-qubit circuits	Good, especially at optimization level 3
Hardware integration	Uniform API across vendors	Deep integration with IBM runtime
Learning curve	Moderate, more explicit control	Lower for IBM-focused work, rich documentation

Choose tket when you need backend portability (running the same code on IBM, Quantinuum, and Rigetti), or when you want maximum control over compilation. Choose Qiskit when you are focused on IBM hardware and want to leverage the built-in primitives and application libraries.

In practice, many researchers use both: tket for compilation and Qiskit for high-level algorithm construction.

Common Mistakes and Pitfalls

1. Forgetting that symbol_substitution modifies in place

# WRONG: loses the symbolic circuit after first iteration
bound = compiled_ansatz
bound.symbol_substitution({theta: 0.5})
# compiled_ansatz is now also modified!

# CORRECT: copy first
bound = compiled_ansatz.copy()
bound.symbol_substitution({theta: 0.5})

Always call .copy() before substitution if you need the symbolic circuit again.

2. Measuring incompatible Pauli terms in one circuit

ZZ and XX cannot be measured in the same circuit because Z and X anticommute on the same qubit. Each requires a different measurement basis. If you try to measure both from a single Z-basis measurement, the XX expectation value will be incorrect.

The correct approach: group compatible Pauli terms (those that commute qubit-wise) and measure each group in a separate circuit. For our H2 Hamiltonian:

• Group 1 (Z-basis): IZ, ZI, ZZ (all diagonal, one circuit)
• Group 2 (X-basis on both qubits): XX (separate circuit)

3. Compiling inside the optimization loop

As discussed in the compile-once section, calling get_compiled_circuit on every iteration wastes time. Compile the symbolic circuit once, then substitute into the compiled version.

4. Applying FullPeepholeOptimise to symbolic circuits and expecting full optimization

FullPeepholeOptimise can simplify circuit structure (cancel adjacent inverse gates, commute gates past each other), but it cannot merge symbolic rotations like Ry(theta) followed by Ry(phi) into Ry(theta + phi) when theta and phi are symbols. The full numerical optimization only works after substitution.

If you need maximum gate reduction, apply FullPeepholeOptimise after symbol_substitution, not before.

5. Confusing qubit mapping with Hamiltonian encoding

The Pauli terms in this tutorial use the parity mapping for fermion-to-qubit encoding. If you generate the Hamiltonian with OpenFermion or PySCF using the Jordan-Wigner mapping instead, you get different Pauli strings. Make sure the mapping used to generate the Hamiltonian matches what your ansatz expects. Mixing mappings produces silently wrong results.

6. Ignoring tket’s angle convention

tket measures rotation angles in half-turns (multiples of pi). A call to Ry(0.5, 0) rotates by 0.5 * pi = pi/2 radians, not 0.5 radians. If you are porting code from Qiskit (which uses radians), divide all angles by pi. Getting this wrong shifts the energy landscape and causes the optimizer to converge to the wrong value.

Next Steps

• Replace the hard-coded Hamiltonian with one generated from PySCF and OpenFermion.
• Implement Pauli term grouping (qubit-wise commutativity) to reduce the number of measurement circuits.
• Try the UCCSD ansatz on larger molecules (LiH, BeH2) with more qubits and parameters.
• Add noise mitigation using tket’s SpamCorrecter for readout error correction on hardware runs.
• Experiment with different optimizers (SPSA, Adam with parameter-shift gradients) to compare convergence.
• Implement the simultaneous perturbation stochastic approximation (SPSA) gradient estimator, which uses only 2 circuit evaluations per gradient step regardless of parameter count.
• Explore operator grouping strategies (sorted insertion, graph coloring) to minimize the number of measurement circuits for larger Hamiltonians.