Skip to content

Multi-qubit: GHZ phase estimation

In this example, we define a circuit to prepare a GHZ state. This state, ubiquitous in quantum information processing tasks, is highly entangled; and in the context of quantum metrology, saturates the Heisenberg limit.

import itertools

import equinox as eqx
import jax
import jax.numpy as jnp
import seaborn as sns
import ultraplot as uplt
from rich.pretty import pprint

from squint import Circuit
from squint.interface.base import SharedGate, Wire
from squint.interface.dv import Conditional, DiscreteVariableState, HGate, RZGate, x
from squint.backends.tensornetwork.simulator import Simulator
from squint.utils import partition_op
from squint.math.information_matrices import quantum_fisher_information_matrix, classical_fisher_information_matrix

n = 3  # number of qubits
wires = [Wire(dim=2, idx=i) for i in range(n)]

circuit = Circuit()
for w in wires:
    circuit.add(DiscreteVariableState(wires=(w,), n=(0,)))

circuit.add(HGate(wires=(wires[0],)))
for i in range(n - 1):
    circuit.add(Conditional(ufunc=x, wires=(wires[i], wires[i + 1])))

circuit.add(
    SharedGate(op=RZGate(wires=(wires[0],), phi=0.0 * jnp.pi), wires=tuple(wires[1:])),
    "phase",
)

for w in wires:
    circuit.add(HGate(wires=(w,)))

pprint(circuit)

params, static = partition_op(circuit, "phase")
sim = Simulator(static=static, params=params)

def forward_probs(p):
    return jnp.abs(sim.forward(p)) ** 2

grad_probs = jax.jacfwd(forward_probs)

phis = jnp.linspace(-jnp.pi, jnp.pi, 100)
params = eqx.tree_at(lambda pytree: pytree.ops["phase"].op.phi, params, phis)

probs = jax.vmap(forward_probs)(params)
grads = jax.vmap(grad_probs)(params).ops["phase"].op.phi
qfims = jax.vmap(lambda p: quantum_fisher_information_matrix(sim.forward, sim.grad, p))(params)
cfims = jax.vmap(lambda p: classical_fisher_information_matrix(forward_probs, grad_probs, p))(params)

colors = sns.color_palette("Set2", n_colors=jnp.prod(jnp.array(probs.shape[1:])))
fig, axs = uplt.subplots(nrows=3, figsize=(6, 4), sharey=False)

for i, idx in enumerate(
    itertools.product(*[list(range(ell)) for ell in probs.shape[1:]])
):
    axs[0].plot(phis, probs[:, *idx], label=f"{idx}", color=colors[i])

    axs[1].plot(phis, grads[:, *idx], label=f"{idx}", color=colors[i])
axs[0].legend()
axs[0].set(ylabel=r"$p(\mathbf{x} | \varphi)$")
axs[1].set(ylabel=r"$\partial_{\varphi} p(\mathbf{x} | \varphi)$")

axs[2].plot(phis, qfims.squeeze(), color=colors[0], label=r"$\mathcal{I}_\varphi^Q$")
axs[2].plot(phis, cfims.squeeze(), color=colors[-1], label=r"$\mathcal{I}_\varphi^C$")
axs[2].set(
    xlabel=r"Phase, $\varphi$",
    ylabel=r"$\mathcal{I}_\varphi^C$",
    ylim=[0, 1.05 * jnp.max(qfims)],
)
axs[2].legend();

img