Expectation values and parameter sweeps¶
Compute Pauli expectation values directly from the simulated state, without sampling, and sweep a circuit’s parameters with exact gradients.
import maestro
from maestro.circuits import QuantumCircuit
ghz = QuantumCircuit()
ghz.h(0)
ghz.cx(0, 1)
ghz.cx(1, 2)
result = ghz.estimate(["ZZI", "XXX", "IZI"])
print(result["expectation_values"])
[0.9999999999999998, 0.9999999999999998, 0.0]
Observables are Pauli strings, given as a list or as one ;-separated string.
Values come back in the order you gave them. The circuit needs no measurements.
For this GHZ state ⟨ZZI⟩ = ⟨XXX⟩ = 1, while a single ⟨Z⟩ averages to 0.
Warning
In a Pauli string, the character at position q acts on qubit q, so "ZII"
is Z on qubit 0. Qiskit’s SparsePauliOp uses the opposite order, where the
same operator is "IIZ".
From an OpenQASM string¶
maestro.simple_estimate takes an OpenQASM 2.0 program instead of a circuit:
ghz_qasm = """
OPENQASM 2.0;
include "qelib1.inc";
qreg q[3];
h q[0];
cx q[0], q[1];
cx q[1], q[2];
"""
values = maestro.simple_estimate(ghz_qasm, "ZZZ;ZZI;XXX")["expectation_values"]
print(values)
[0.0, 0.9999999999999998, 0.9999999999999998]
⟨ZZZ⟩ is 0 here: |000⟩ and |111⟩ have opposite parity under Z⊗Z⊗Z.
On a different simulation method¶
Pass a SimulatorConfig to estimate on another method, for example an MPS:
mps_config = maestro.SimulatorConfig(
simulation_type=maestro.SimulationType.MatrixProductState,
max_bond_dimension=16,
)
print(ghz.estimate(["ZZI", "XXX"], config=mps_config)["expectation_values"])
[1.0000000000000002, 1.0000000000000002]
Reductions before simulation¶
Before allocating a simulator, Maestro shrinks the problem:
Causal-cone reduction (
enable_causal_cone_reduction=True, default) removes gates and idle qubits outside the backward lightcone of the requested observables.Observable clustering (
enable_observable_clustering=True, default) splits the observables into groups whose lightcones are disjoint, and reduces and simulates each group on its own smaller register. Observables whose lightcones overlap, even indirectly through other observables, stay in one group.max_cluster_qubitsalso splits overlapping lightcones (see below).Parity-check reduction (
auto_reduce=True, default) detects invariant check subspaces, such as QAOA parity ladders with uncomputed CNOT frames, and simulates an effective register of dimension 2^rank(checks) on the backend you configured.
The result reports what was applied. Here two independent entangled pairs give two groups of two qubits each:
pairs = QuantumCircuit()
pairs.h(0)
pairs.cx(0, 1)
pairs.ry(2, 0.4)
pairs.cx(2, 3)
result = pairs.estimate(["XXII", "IIXX"])
for key in ["expectation_values", "auto_reduced", "qubits_before", "qubits_after",
"clusters_count", "largest_cluster_qubits", "cluster_methods"]:
print(f"{key}: {result[key]}")
expectation_values: [0.9999999999999998, 0.38941834230865047]
auto_reduced: False
qubits_before: 4
qubits_after: 2
clusters_count: 2
largest_cluster_qubits: 2
cluster_methods: [SimulationType.Statevector, SimulationType.Statevector]
qubits_after reports the largest group, and cluster_methods lists each
group’s SimulationType. The source circuit is unchanged. Turn the
reductions off with the matching config fields:
config = maestro.SimulatorConfig(
enable_causal_cone_reduction=False,
enable_observable_clustering=False,
auto_reduce=False,
)
result = pairs.estimate(["XXII", "IIXX"], config=config)
print(result["qubits_after"], result["clusters_count"])
4 1
Requires GPU access
GPU TensorNetwork estimates whose effective register is a single qubit run on the GPU Statevector of the same device instead, because cuTensorNet’s MPS factorization needs at least two qubits. See Get GPU access.
Bounding group size¶
When many observables are spread across a device, their lightcones often chain
together into one group covering every qubit, and clustering saves nothing.
max_cluster_qubits caps each group’s combined lightcone. Each group starts
from the largest remaining lightcone and adds the overlapping observable that
adds the fewest new qubits, while the union fits the budget. An observable whose
own lightcone is larger than the budget gets a group of its own.
Here every nearest-neighbour ⟨XX⟩ on a 12-qubit brickwork circuit overlaps the next, so without a budget they form one 12-qubit group:
n = 12
brickwork = QuantumCircuit()
for q in range(n):
brickwork.ry(q, 0.1 * (q + 1))
for q in range(0, n - 1, 2):
brickwork.cx(q, q + 1)
for q in range(1, n - 1, 2):
brickwork.cx(q, q + 1)
neighbours = ["I" * q + "XX" + "I" * (n - 2 - q) for q in range(n - 1)]
for budget in [None, 4]:
config = maestro.SimulatorConfig(
simulation_type=maestro.SimulationType.PauliPropagator,
max_cluster_qubits=budget,
)
result = brickwork.estimate(neighbours, config=config)
print(f"budget={budget}: {result['clusters_count']} groups, "
f"largest {result['largest_cluster_qubits']} qubits")
budget=None: 1 groups, largest 12 qubits
budget=4: 5 groups, largest 4 qubits
Smaller groups mean smaller registers but more simulator runs, so the best budget depends on the backend. A budget helps backends whose cost grows with register width (Pauli propagation, statevector); matrix product states pay per group and gain little. GPU Pauli propagation, which requires GPU access, ignores the budget and keeps all observables in one run.
Note
With PauliPropagator and pp_gates_between_trims above 1, trims land on
different gates in each group’s reduced circuit, so values can change slightly
with the budget, within the truncation error. Trimming and deduplicating after
every gate (pp_gates_between_trims=1, pp_gates_between_deduplications=1)
makes values independent of the grouping. Measurements, resets, quantum
channels, classical control and distributed GPU layouts fall back to full
execution.
Parameter sweeps and gradients¶
To evaluate one circuit at many parameter values, compile it once with
AdjointProgram and pass one row of parameters per point. The parameters are
the angles of the p, rx, ry, rz, cp, crx, cry, crz, u and
cu gates, in circuit order. jacobian also returns exact derivatives,
computed with the adjoint method on the CPU statevector.
import numpy as np
ansatz = QuantumCircuit()
ansatz.ry(0, 0.0) # parameter 0
ansatz.cx(0, 1)
ansatz.rz(1, 0.0) # parameter 1
program = maestro.AdjointProgram(ansatz)
print("parameters:", program.num_params)
thetas = np.linspace(0, np.pi, 5)
params = np.column_stack([thetas, np.zeros_like(thetas)])
result = program.jacobian(["ZI", "XX"], params)
values = np.asarray(result["expectation_values"]) # (points, observables)
grads = np.asarray(result["jacobian"]) # (points, observables, parameters)
for theta, (zi, xx), dzi in zip(thetas, values, grads[:, 0, 0]):
print(f"theta={theta:.2f} <ZI>={zi:+.3f} <XX>={xx:+.3f} d<ZI>/dtheta={dzi:+.3f}")
parameters: 2
theta=0.00 <ZI>=+1.000 <XX>=+0.000 d<ZI>/dtheta=+0.000
theta=0.79 <ZI>=+0.707 <XX>=+0.707 d<ZI>/dtheta=-0.707
theta=1.57 <ZI>=+0.000 <XX>=+1.000 d<ZI>/dtheta=-1.000
theta=2.36 <ZI>=-0.707 <XX>=+0.707 d<ZI>/dtheta=-0.707
theta=3.14 <ZI>=-1.000 <XX>=+0.000 d<ZI>/dtheta=-0.000
Call expectation_values instead of jacobian when you only need values, or
vjp to backpropagate a loss through the circuit. maestro.adjoint_jacobian
does the same in one call without keeping a compiled program.
When only the end of a circuit changes between runs,
maestro.PrefixCheckpointedSimulator simulates the shared prefix once
and runs each suffix from that saved state.
Mirror fidelity¶
Mirror fidelity runs a circuit, then the inverse of every gate in reverse order, and returns the probability of getting back to |0…0⟩. An exact simulation gives 1.0, so the shortfall measures simulation error, for example from a low MPS bond dimension.
mirror = QuantumCircuit()
mirror.h(0)
mirror.cx(0, 1)
mirror.rx(0, np.pi / 4)
print("sampled (1024 shots):", maestro.mirror_fidelity(mirror))
print("exact amplitude: ", maestro.mirror_fidelity(mirror, full_amplitude=True))
print("circuit method, MPS: ", mirror.mirror_fidelity(config=mps_config, shots=10000))
sampled (1024 shots): 1.0
exact amplitude: 0.9999999999999993
circuit method, MPS: 1.0
By default it samples 1024 shots. full_amplitude=True computes it exactly,
which only suits small circuits. Measurements in the circuit are skipped; only
unitary gates are inverted.
Observable mirror fidelity¶
For wide circuits, the full-register probability is intractable for Pauli
propagation and collapses towards 0 for a truncated MPS. Pass observables to
evaluate Pauli strings on the mirror circuit C C† instead. Z-type strings have
ideal value 1.0, so truncation error shows up as a value below 1.0 on exactly
the observables you study.
pp_config = maestro.SimulatorConfig(
simulation_type=maestro.SimulationType.PauliPropagator,
pp_max_pauli_weight=4,
pp_coefficient_threshold=1e-6,
)
print(mirror.mirror_fidelity(config=pp_config, observables="ZZ")) # one string: float
print(mirror.mirror_fidelity(config=pp_config, observables=["ZZ", "ZI"])) # list
1.0
[1.0, 1.0]
In this mode, circuit optimization is turned off so gate and inverse pairs are
not cancelled, and shots and full_amplitude are ignored.
Inner products¶
maestro.inner_product(c1, c2) returns the complex overlap ⟨ψ₁|ψ₂⟩ of two
circuits’ output states. It evaluates ⟨0|U₁†U₂|0⟩ directly, so with an MPS it
never builds a full statevector. abs(overlap) ** 2 is the state fidelity, and
the phase is the relative global phase.
bell_a = QuantumCircuit()
bell_a.h(0)
bell_a.cx(0, 1)
bell_b = QuantumCircuit()
bell_b.h(0)
bell_b.cx(0, 1)
plus = QuantumCircuit()
plus.h(0)
minus = QuantumCircuit()
minus.x(0)
minus.h(0)
print("<bell|bell> =", maestro.inner_product(bell_a, bell_b))
print("|<+|->| =", abs(maestro.inner_product(plus, minus)))
<bell|bell> = (1+0j)
|<+|->| = 0.0
Incremental time evolution¶
maestro.incremental_evolve applies a Trotter step repeatedly on one persistent
simulator and records expectation values at chosen steps, instead of
re-simulating from step 0 for each time point.
init = QuantumCircuit()
init.x(0)
step = QuantumCircuit()
step.rx(0, 0.1)
step.rx(1, 0.2)
step.cx(0, 1)
step.rz(1, 0.05)
step.cx(0, 1)
result = maestro.incremental_evolve(
init_circuit=init,
trotter_step=step,
measure_at_steps=[10, 20, 30, 50, 100],
observables=["ZI", "IZ", "ZZ"],
config=mps_config,
)
for s, evs in zip(result["steps"], result["expectation_values"]):
print(f"step {s:3d}: <ZI>={evs[0]:+.4f} <IZ>={evs[1]:+.4f} <ZZ>={evs[2]:+.4f}")
step 10: <ZI>=-0.5487 <IZ>=-0.3863 <ZZ>=+0.1960
step 20: <ZI>=+0.3491 <IZ>=-0.5084 <ZZ>=-0.3431
step 30: <ZI>=+0.9235 <IZ>=+0.8119 <ZZ>=+0.7427
step 50: <ZI>=-0.0807 <IZ>=-0.5004 <ZZ>=+0.0056
step 100: <ZI>=+0.9681 <IZ>=-0.4246 <ZZ>=-0.4471
Measuring at n points across N steps costs about n·N/2 steps when each point starts a fresh circuit, and N steps incrementally: 300 instead of about 3,000 for N = 300 and 20 points.
Parameter |
Type |
Default |
Description |
|---|---|---|---|
|
|
Required |
State preparation, run once |
|
|
Required |
One Trotter step, applied repeatedly |
|
|
Required |
Steps at which to compute expectation values |
|
|
Required |
Pauli strings to evaluate |
|
|
|
Backend configuration |