Circuits and exact gradients

This tutorial builds circuits, samples them, and computes gradients of expectation values in three ways. It also shows why the adjoint method is the default: its cost does not grow with the number of parameters.

import time

import matplotlib.pyplot as plt
import numpy as np

import metatron_qso as q

print("QSO", q.__version__)

QSO 1.2.0

A GHZ state

Qubit 0 is the least significant bit, so bit strings print qubit 2 first.

ghz = q.Circuit(3).h(0).cx(0, 1).cx(1, 2)
state = q.simulate(ghz)
print(np.round(state.amplitudes, 4))
counts = state.sample(2000, seed=7)
print(counts)
print("entanglement entropy of qubit 0:", state.entanglement_entropy([0]), "bit")

[0.7071+0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j
 0.    +0.j 0.7071+0.j]
{'000': 1044, '111': 956}
entanglement entropy of qubit 0: 0.9999999999999999 bit

fig, ax = plt.subplots(figsize=(5, 2.5))
labels = [q.bitstring(k, 3) for k in range(8)]
ax.bar(labels, [counts.get(b, 0) for b in labels])
ax.set_xlabel("outcome")
ax.set_ylabel("counts (2000 shots)")
plt.tight_layout()
plt.show()

An energy landscape and its derivative

One parameter θ, observable H = Z₀Z₁ − 0.5·X₀. The adjoint gradient is exact; finite differences agree to their truncation error.

c = q.Circuit(2)
theta = c.new_parameter()
c.ry(0, theta).cx(0, 1).rx(1, 0.5 * theta + 0.2)
h = q.PauliSum("Z0 Z1 - 0.5 X0")

xs = np.linspace(-np.pi, np.pi, 121)
values, adjoint, fd = [], [], []
for x in xs:
    v, g = q.expectation_and_gradient(c, h, [x])
    _, g_fd = q.expectation_and_gradient(c, h, [x], method="finite_difference")
    values.append(v)
    adjoint.append(g[0])
    fd.append(g_fd[0])
print("largest |adjoint − finite difference|:", np.max(np.abs(np.array(adjoint) - fd)))

fig, ax = plt.subplots(figsize=(6, 3))
ax.plot(xs, values, label="⟨H⟩(θ)")
ax.plot(xs, adjoint, label="d⟨H⟩/dθ (adjoint)")
ax.set_xlabel("θ")
ax.legend()
plt.tight_layout()
plt.show()

largest |adjoint − finite difference|: 2.3853558017705723e-11

Cost of a full gradient

For an ansatz with P parameters, parameter shift needs about 2P circuit evaluations; the adjoint method needs about three, whatever P is.

h = q.transverse_field_ising(8, 1.0, 0.7)

def seconds(f, repeats=5):
    f()
    t = time.perf_counter()
    for _ in range(repeats):
        f()
    return (time.perf_counter() - t) / repeats

layers = [1, 2, 4, 8, 12]
params, t_adj, t_shift = [], [], []
for L in layers:
    ansatz = q.efficient_su2(8, L)
    theta = np.linspace(0, 1, ansatz.num_parameters)
    params.append(ansatz.num_parameters)
    t_adj.append(seconds(lambda: q.expectation_and_gradient(ansatz, h, theta)))
    t_shift.append(seconds(lambda: q.expectation_and_gradient(ansatz, h, theta, method="parameter_shift")))

fig, ax = plt.subplots(figsize=(5, 3))
ax.loglog(params, t_adj, "o-", label="adjoint")
ax.loglog(params, t_shift, "s-", label="parameter shift")
ax.set_xlabel("number of parameters")
ax.set_ylabel("seconds per gradient")
ax.legend()
plt.tight_layout()
plt.show()
print(f"{params[-1]} parameters: parameter shift is {t_shift[-1] / t_adj[-1]:.0f}× slower")

208 parameters: parameter shift is 122× slower

OpenQASM round trip

Export binds the parameters; import gives an equivalent circuit.

text = c.to_qasm([0.4], measure=True)
print(text)
again = q.Circuit.from_qasm(text)
print("fidelity:", q.simulate(c, [0.4]).fidelity(q.simulate(again)))

OPENQASM 2.0;
include "qelib1.inc";
qreg q[2];
creg c[2];
ry(0.4) q[0];
cx q[0],q[1];
rx(0.4) q[1];
measure q -> c;

fidelity: 1.0