willzeng

from pyquil.quil import Program
from pyquil.api import QPUConnection
from pyquil.gates import H, CNOT

qpu = QPUConnection(device_name='19Q-Acorn')

bell_state = Program(H(0), CNOT(0, 1))
result = qpu.run_and_measure(bell_state, [0, 1])

willzeng

# the state of N qubits is a complex vector of dimension 2^N
n_qubits = 2
dimension = 2 ** n_qubits
state = [1j, 0, 0, 0]

def print_probs(state, n_qubits):
    # the elements of this state are squared to calculate outcome probabilities
    for bitstring in range(n_qubits ** 2):
        probability = np.abs(state[bitstring]) ** 2

willzeng

def random_unitary(dim: int) -> np.ndarray:
  # follows the algorithm in https://arxiv.org/pdf/math-ph/0609050.pdf
  # returns a unitary of size dim x dim
  Z = np.array([np.random.normal(0, 1) + np.random.normal(0, 1) * 1j for _ in range(dim ** 2)]).reshape(dim, dim)
  Q, R = np.linalg.qr(Z)
  diag = np.diagonal(R)
  lamb = np.diag(diag) / np.absolute(diag)
  unitary = np.matmul(Q, lamb)

  # this condition asserts that the matrix is unitary

willzeng

def simulate_probability(unitary: np.ndarray, bitstring:int) -> float:
  # simulates the probability of measuring bitstring when evolving from the ground state
  # according to the quantum program given unitary
  return np.abs(unitary[bitstring, 0])**2


def quantum_sample_probability(n_qubits: int, trials: int) -> List:
  # returns the probabilities of a randomly chosen bistring outcome over "trials" number of different
  # random quantum programs on n_qubits

willzeng

n_qubits = 4
porter_thomas = quantum_sample_probability(n_qubits, 10_000)

# theoretical Porter-Thomas distribution
dim = 2**n_qubits
xspace = np.linspace(0.0, 1.0, 100)
yspace = dim * np.exp(-dim*xspace)

# plot both empirical and theoretical calculations
plt.figure(figsize=(9, 6))

willzeng

def fidelity_xeb(n_qubits: int, trials: int, n_samples: int, sampler: Callable[[np.ndarray, int], float]) -> float:
  dim = 2**n_qubits
  # keep track of the ideal simulated probabilities
  ideal_probs = []
  # loop over the random programs
  for _ in range(trials):
      unitary = random_unitary(dim)
      sample_probs = [sampler(unitary, bb) for bb in range(dim)]
      samples = np.random.choice(dim, size=n_samples, p=sample_probs)
      for sample in samples:

willzeng

# sample f_xeb using the same parameters as in the Google paper
n_qubits = 6
f_xeb = fidelity_xeb(n_qubits=n_qubits, trials=10, n_samples=10**5, sampler=simulate_probability)
print("Empirical FXEB: ", f_xeb)

willzeng

dim = 2**n_qubits

def unif_dist(unitary, bitstring):
    return 1/dim # all bitstrings have the same probability
unif_xeb = fidelity_xeb(n_qubits=n_qubits, trials=10, n_samples=10**5, sampler=unif_dist)
print("Empirical FXEB of a uniform sampler: ", unif_xeb)

willzeng

def fidelity_xeb_noisy(n_qubits: int, trials: int, n_samples: int, prob_no_error: float):
    dim = 2**n_qubits

    # keep track of ideal output probabilities
    ideal_probs = []

    # identify depolarizing operators over n-qubit space
    depolarizing_ops = []
    for x in itertools.product(paulis, repeat=n_qubits):
        op = functools.reduce(lambda a, b: np.kron(a, b), x)

willzeng

# run the noisy experiment
noisy_xeb = fidelity_xeb_noisy(n_qubits=6, trials=10**3, n_samples=10, prob_no_error=0.7)
print("Empirical FXEB of a noisy simulation: ", noisy_xeb)