Random Number Generation

Download this notebook - random_numbers.ipynb

In this simple example we look at using random numbers with Guppy and Selene. The example randomly applies simple identities to a qubit. This style of program is often used in Randomized Benchmarking as the rate of unexpected outcomes can indicate error rates in the system.

There are two ways of using random numbers:

1. Generate them ahead of time in Python and store in the Guppy program.

2. Generate them on the fly.

We will look at both options.

import random

from guppylang import guppy
from guppylang.std.builtins import barrier, comptime, result
from guppylang.std.qsystem.random import RNG
from guppylang.std.qsystem.utils import get_current_shot
from guppylang.std.quantum import measure, qubit, x, y, z
from guppylang.defs import GuppyFunctionDefinition

Let’s decide on a seed and make sure we use it for both Python and Guppy randomness, to get deterministic results.

INIT_SEED = 91919
random.seed(INIT_SEED)

# number of shots we will run
N_SHOTS = 100

Let’s first write our simple experiment functions in Guppy.

@guppy
def rand_pauli(q: qubit, randval: int) -> None:
    """Applies a random pauli to a qubit given a number in {0, 1, 2}"""
    if randval == 0:
        x(q)
    elif randval == 1:
        y(q)
    else:
        z(q)


@guppy
def pauli_identity(q: qubit, randval: int) -> None:
    """Apply the same random pauli twice separated by a barrier.
    Acts as identity in noiseless case since Paulis are self-adjoint."""
    rand_pauli(q, randval)
    barrier(q)
    rand_pauli(q, randval)

Pre-compiled randomness

First we generate all the random numbers we need ahead of time in Python. This can be stored as an array of integers in a Guppy program using the comptime expression. At runtime we use the current shot value to index in to the array and get a fresh random integer for the shot.

We report out both the result of measuring the qubit, and the random number used.

# randint upper bound is inclusive
randints = [random.randint(0, 2) for _ in range(N_SHOTS)]


@guppy
def main() -> None:
    q = qubit()
    randval = comptime(randints)[get_current_shot()]
    pauli_identity(q, randval)
    result("rng", randval)
    result("c", measure(q))

We also define a helper function to run our Guppy program. It runs N_SHOTS repeats of the program using Stim simulation (since our program is Clifford).

No noise model is used so the measurement result “c” should always be 0. The function reports the mean of the random numbers over the shots - we expect this to be around 1.

def run(main_def: GuppyFunctionDefinition) -> float:
    """Run the rng example and return the mean of the random numbers generated.

    Checks measurement always returns 0 (identity operation on 0 state).
    """
    shots = main_def.emulator(n_qubits=1).stabilizer_sim().with_seed(INIT_SEED).with_shots(N_SHOTS).run()
    rands = []
    for out in shots.collated_shots():
        assert out["c"] == [0]
        rands.extend(out["rng"])
    return sum(rands) / len(rands)

run(main)

1.04

Runtime randomness

The Selene emulator supports random number generation at runtime. This has slightly different behaviour to Python:

• To ensure reproducibility you always have to provide a seed to initialise an RNG object. Methods on this object can return random values.

• The RNG type is linear (like a qubit) and has to be explicitly discarded. This is also to ensure determinism as it prevents random number generation calls being re-ordered by the compiler.

The benefit over generating ahead of time is that for large experiments you can avoid storing a huge amount of random numbers in the program.

We can use the RNG type to re-write our main function, using our INIT_SEED along with a shift by current shot number to get a different number in each shot:

@guppy
def main() -> None:
    q = qubit()
    rng = RNG(comptime(INIT_SEED) + get_current_shot())
    # random_int_bounded upper bound is exclusive,
    # chooses a number in {0, 1, 2}
    randval = rng.random_int_bounded(3)
    rng.discard()
    pauli_identity(q, randval)
    result("rng", randval)
    result("c", measure(q))

run(main)

0.95