Quantum Volume¶

Quantum volume is a measurement of the error amount characterizing a chosen quantum hardware. The quantum volume is a result of running a circuit based on principles of randomness and statistical analysis, which provides a single number to compare between different hardwares.

The scheme of the quantum volume is as following [1]:

1. For a given number of qubits $n$, a circuit is made of $n$ quantum layer.
2. Each layer consist of unitary operation between pairing of the $n$ qubits. The pairs are chosen at random, and in the case $n$ is odd, one will not have operation.
3. The unitary operation between each pair is the Haar random matrix, i.e SU(4) operation containing random complex number in such manner that the probability to measure a quantum state is kept with uniform distribution.
4. For a single circuit made for $n$ qubits, a measurements are made, and the heavy output probability (i.e. the probability to measure the states above the median value) is calculated. Due to the nature of the distribution of random complex number, one can evaluate that for an ideal case (no noises), the heavy output probability should be ~0.85. For assessment of the quantum volume the demand is subsided to the following inequality: (1) P{heavy_outputs} >= 2/3 .
5. For a given output, in order to get the quantum volume, one repeats steps 1-4 for increasing number of qubits until the inequality described in step 4 does not hold. In order to assure that, the circuits are created many times and the average and standard deviation are taken into account.
6. The quantum volume will be the 2 to the power the number of qubits, such that pass inequality (1), under the procedure described in 1-5.

The heavy output probability is a good measurement of the quality of the circuit, as noise will reduce the probabilities to the uniform distribution. While so, one needs to consider that there are many components to the results of the procedure - including the hardware noises, but also the connectivity map, the quality of the transpilation, and even the quantum software which translated the circuit into basis gates for the hardware, and thus contributes to the circuit's depth.

In this demonstration we show the code to implement the steps to calculate the quantum volume using the Classiq platform, and show example to such calculation made for several quantum simulators and hardwares.

Step 1: create Haar random unitary matrix¶

We hereby create a function with generate a (n,n) sized Haar random unitary matrix [2]. This matrix contains random complex number which distribute evenly in the $2^n$ space of quantum states. The Haar distribution tells us how to weight the elements of $𝑈(𝑁)$ such that such uniform distribution occurs in the parameter space.

In [ ]:

import numpy as np
from numpy.linalg import qr
from scipy.stats import unitary_group


def haar(n):
    # n=4
    u1 = unitary_group.rvs(n)
    u2 = unitary_group.rvs(n)
    Z = u1 + 1j * u2

    Q, R = qr(Z)

    Lambda = np.diag([R[i, i] / np.abs(R[i, i]) for i in range(n)])

    return np.dot(Q, Lambda)

import numpy as np from numpy.linalg import qr from scipy.stats import unitary_group def haar(n): # n=4 u1 = unitary_group.rvs(n) u2 = unitary_group.rvs(n) Z = u1 + 1j * u2 Q, R = qr(Z) Lambda = np.diag([R[i, i] / np.abs(R[i, i]) for i in range(n)]) return np.dot(Q, Lambda)

step 2: create quantum volume circuit¶

The function qv_model creates for given $N$ number of qubit the quantum volume model. For $N$ qubits, the circuit must include $N$ quantum volume layers. The layers are built using qv_layer function which create a random pairing between the $N$ qubits (for odd number, a random chosen qubit will not be operated). Between each pair, a unitary gate will be operated which consist of Haar random unitary matrix of size 4.

In [ ]:

import math
import random

from classiq import Model, QReg, QUInt, RegisterUserInput
from classiq.builtin_functions import UnitaryGate


def qv_layer(N, model, qubit_to_wires, qubit_list):
    for idx in range(math.floor(N / 2)):
        # step 1: isolate the qubit pairs for the layers
        a = qubit_list[idx]
        b = qubit_list[math.floor(N / 2) + idx]

        # step 2: generate the random matix : this need to change for the random matrix when its possible
        gate_matrix = haar(4).tolist()
        unitary_gate_params = UnitaryGate(data=gate_matrix)
        uni_out = model.UnitaryGate(
            unitary_gate_params,
            in_wires=QReg.concat(qubit_to_wires[a], qubit_to_wires[b]),
        )

        # step 3: refresh the wiring map
        qubit_to_wires[a], qubit_to_wires[b] = (
            uni_out["TARGET"][0],
            uni_out["TARGET"][1],
        )

    return (model, qubit_to_wires)

import math import random from classiq import Model, QReg, QUInt, RegisterUserInput from classiq.builtin_functions import UnitaryGate def qv_layer(N, model, qubit_to_wires, qubit_list): for idx in range(math.floor(N / 2)): # step 1: isolate the qubit pairs for the layers a = qubit_list[idx] b = qubit_list[math.floor(N / 2) + idx] # step 2: generate the random matix : this need to change for the random matrix when its possible gate_matrix = haar(4).tolist() unitary_gate_params = UnitaryGate(data=gate_matrix) uni_out = model.UnitaryGate( unitary_gate_params, in_wires=QReg.concat(qubit_to_wires[a], qubit_to_wires[b]), ) # step 3: refresh the wiring map qubit_to_wires[a], qubit_to_wires[b] = ( uni_out["TARGET"][0], uni_out["TARGET"][1], ) return (model, qubit_to_wires)

In [ ]:

def qv_model(N):
    # step 1: create the qubit list
    qubit_list = list(range(N))

    # step 2: create the model, qubits and wiring map
    inputs = RegisterUserInput(size=N)
    model = Model()
    input_dict = model.create_inputs({"TARGET": QReg[N]})
    qubit_to_wires = [input_dict["TARGET"][k] for k in range(N)]

    # step 3: send the model to add layers
    for idx in range(N):
        model, qubit_to_wires = qv_layer(N, model, qubit_to_wires, qubit_list)
        random.shuffle(qubit_list)

    model.sample()
    return model.get_model()

def qv_model(N): # step 1: create the qubit list qubit_list = list(range(N)) # step 2: create the model, qubits and wiring map inputs = RegisterUserInput(size=N) model = Model() input_dict = model.create_inputs({"TARGET": QReg[N]}) qubit_to_wires = [input_dict["TARGET"][k] for k in range(N)] # step 3: send the model to add layers for idx in range(N): model, qubit_to_wires = qv_layer(N, model, qubit_to_wires, qubit_list) random.shuffle(qubit_list) model.sample() return model.get_model()

Step 3: execution and analysis¶

The execution and analysis part consist of 2 functions:

• execute_qv which send a given number of shots, circuit and preferences sends the circuit to execution in the quantum hardware specified in the preferences. This function returns the results of the execution from the hardware.
• heavy_outputs_prob which analyze the results from execution, and returns the heavy output probability, i.e. the probability for a single state in the space to be greater then the median value (median = "middle" of a sorted list of numbers).

The functionround_significant rounds a number for one significant figure.

In [ ]:

from classiq import execute, set_quantum_program_execution_preferences
from classiq.execution import (
    AzureBackendPreferences,
    ClassiqBackendPreferences,
    ExecutionDetails,
    ExecutionPreferences,
)


def execute_qv(qprog, num_shots, preferences):
    qprog = set_quantum_program_execution_preferences(
        qprog,
        ExecutionPreferences(num_shots=num_shots, backend_preferences=preferences),
    )
    results = execute(qprog).result()
    return results[0].value

from classiq import execute, set_quantum_program_execution_preferences from classiq.execution import ( AzureBackendPreferences, ClassiqBackendPreferences, ExecutionDetails, ExecutionPreferences, ) def execute_qv(qprog, num_shots, preferences): qprog = set_quantum_program_execution_preferences( qprog, ExecutionPreferences(num_shots=num_shots, backend_preferences=preferences), ) results = execute(qprog).result() return results[0].value

In [ ]:

def heavy_outputs_prob(results):
    d = list(results.counts.values())
    med = np.median(d)
    heavy_outputs_prob = 0
    # print(med)
    for count, item in enumerate(d):
        if item >= med:
            heavy_outputs_prob = heavy_outputs_prob + item
    return heavy_outputs_prob

def heavy_outputs_prob(results): d = list(results.counts.values()) med = np.median(d) heavy_outputs_prob = 0 # print(med) for count, item in enumerate(d): if item >= med: heavy_outputs_prob = heavy_outputs_prob + item return heavy_outputs_prob

In [ ]:

from math import floor, log10


def round_significant(x):
    return round(x, -int(floor(log10(abs(x)))))

from math import floor, log10 def round_significant(x): return round(x, -int(floor(log10(abs(x)))))

Step 4: find quantum volume algorithm¶

Using the previously defined functions, find_qv finds the quantum volume value, for defined parameters including hardware definitions. The find_qv function, send for each number of qubits defined (between min_qubit and max_qubits) the value of heavy output probability. This is repeated num_trials times. Then, the heavy output probability is averaged, and the standard deviation is calculated. If the number of qubits chosen for the circuit is less than the number of qubits in the chosen hardware, the qubits will be randomly picked for run according to the rules of the hardware provider.

The quantum volume qubits number is defined as the larger number of qubits for which the heavy output probability, decrease by two sigma (2 times the standard deviation), is more or equal to 2/3. The quantum volume will be 2 to the power the number of quantum volume qubits.

One must note, that if the result given for the log2 of the quantum volume is the same as the chosen max_qubits, there is a possibility that the quantum volume is greater then found by the function, and we recommend to run the program for a greater span.

In [ ]:

from tqdm import tqdm

from classiq import synthesize


def find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences):
    ### initialization
    qubit_num = range(min_qubits, max_qubits + 1)
    heavy_list = np.zeros(max_qubits - min_qubits + 1)
    std_list = np.zeros(max_qubits - min_qubits + 1)
    qubit_v = 0

    ### calculate the heavy outputs for each number of qubits
    for num in tqdm(qubit_num):
        heavy_outputs = 0
        std = 0
        heavytostd = np.zeros(num_trials)
        for idx in tqdm(range(num_trials)):
            model = qv_model(num)
            qprog = synthesize(model)
            results = execute_qv(qprog, num_shots, preferences)
            heavy_temp = heavy_outputs_prob(results)
            heavy_outputs = heavy_outputs + heavy_temp
            heavytostd[idx] = heavy_temp
        s = num - min_qubits
        heavy_list[s] = heavy_outputs / (num_trials * num_shots)
        temp_hl = heavy_outputs / (num_trials * num_shots)
        std = np.std(heavytostd) / (num_trials * num_shots)
        std_list[s] = std
        temp_std = round_significant(std)
        print(
            "for",
            num,
            "qubits the heavy outputs probability is:",
            temp_hl,
            "with",
            temp_std,
            "standard deviation",
        )

    ###determine the quantum volume
    for num in qubit_num:
        s = num - min_qubits
        heavy_is = heavy_list[s] - 2 * (std_list[s])
        if heavy_is >= 2 / 3:
            qubit_v = num
        else:
            break

    qv = 2**qubit_v
    print("##### The quantum volume is", qv, "#####")

    return qv

from tqdm import tqdm from classiq import synthesize def find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences): ### initialization qubit_num = range(min_qubits, max_qubits + 1) heavy_list = np.zeros(max_qubits - min_qubits + 1) std_list = np.zeros(max_qubits - min_qubits + 1) qubit_v = 0 ### calculate the heavy outputs for each number of qubits for num in tqdm(qubit_num): heavy_outputs = 0 std = 0 heavytostd = np.zeros(num_trials) for idx in tqdm(range(num_trials)): model = qv_model(num) qprog = synthesize(model) results = execute_qv(qprog, num_shots, preferences) heavy_temp = heavy_outputs_prob(results) heavy_outputs = heavy_outputs + heavy_temp heavytostd[idx] = heavy_temp s = num - min_qubits heavy_list[s] = heavy_outputs / (num_trials * num_shots) temp_hl = heavy_outputs / (num_trials * num_shots) std = np.std(heavytostd) / (num_trials * num_shots) std_list[s] = std temp_std = round_significant(std) print( "for", num, "qubits the heavy outputs probability is:", temp_hl, "with", temp_std, "standard deviation", ) ###determine the quantum volume for num in qubit_num: s = num - min_qubits heavy_is = heavy_list[s] - 2 * (std_list[s]) if heavy_is >= 2 / 3: qubit_v = num else: break qv = 2**qubit_v print("##### The quantum volume is", qv, "#####") return qv

Running examples¶

We will hereby run the code to find the quantum volume of several quantum simulators and hardwares.

1. Running with IBM simulator¶

In [ ]:

num_trials = 10  # number of times to run qv circuit for each number of qubits. best - 200 or more
num_shots = 100  # number of runs for each execution. best - 1000 or more
preferences = ClassiqBackendPreferences(backend_name="aer_simulator")
min_qubits = 3
max_qubits = 6

qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

num_trials = 10 # number of times to run qv circuit for each number of qubits. best - 200 or more num_shots = 100 # number of runs for each execution. best - 1000 or more preferences = ClassiqBackendPreferences(backend_name="aer_simulator") min_qubits = 3 max_qubits = 6 qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

Since this is a simulator with no errors, we expact the heavy output probability for any number of qubits to be approximatly 0.85

2. Running with Rigetti Aspen M-3¶

In [ ]:

num_trials = 10  # number of times to run qv circuit for each number of qubits
num_shots = 3  # number of runs for each execution
preferences = AzureBackendPreferences(backend_name="Rigetti.Qpu.Aspen-M-3")
min_qubits = 2
max_qubits = 3

# qv = find_qv_trials, num_(numshots, min_qubits, max_qubits, preferences)

num_trials = 10 # number of times to run qv circuit for each number of qubits num_shots = 3 # number of runs for each execution preferences = AzureBackendPreferences(backend_name="Rigetti.Qpu.Aspen-M-3") min_qubits = 2 max_qubits = 3 # qv = find_qv_trials, num_(numshots, min_qubits, max_qubits, preferences)

3. Running with IBM machines¶

We will attempt to run a few IBM machines: ibmq_lima with reported quantum volume of 8; ibmq_quito with reported quantum volume of 16; ibmq_manila with reported quantum volume of 32 [3].

In [ ]:

from classiq.execution import IBMBackendPreferences, IBMBackendProvider

ibm_provider = IBMBackendProvider()

preferences = IBMBackendPreferences(
    backend_name="ibmq_lima",
    access_token="insert_token_number",
    provider=ibm_provider,
)

num_trials = 5  # number of times to run qv circuit for each number of qubits
num_shots = 10  # number of runs for each execution
min_qubits = 2
max_qubits = 4

# qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

from classiq.execution import IBMBackendPreferences, IBMBackendProvider ibm_provider = IBMBackendProvider() preferences = IBMBackendPreferences( backend_name="ibmq_lima", access_token="insert_token_number", provider=ibm_provider, ) num_trials = 5 # number of times to run qv circuit for each number of qubits num_shots = 10 # number of runs for each execution min_qubits = 2 max_qubits = 4 # qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

In [ ]:

from classiq.execution import IBMBackendPreferences, IBMBackendProvider

ibm_provider = IBMBackendProvider()

preferences = IBMBackendPreferences(
    backend_name="ibmq_quito",
    access_token="insert_token_number",
    provider=ibm_provider,
)

num_trials = 1  # number of times to run qv circuit for each number of qubits
num_shots = 10  # number of runs for each execution
min_qubits = 2
max_qubits = 3

# qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

from classiq.execution import IBMBackendPreferences, IBMBackendProvider ibm_provider = IBMBackendProvider() preferences = IBMBackendPreferences( backend_name="ibmq_quito", access_token="insert_token_number", provider=ibm_provider, ) num_trials = 1 # number of times to run qv circuit for each number of qubits num_shots = 10 # number of runs for each execution min_qubits = 2 max_qubits = 3 # qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

In [ ]:

from classiq.execution import IBMBackendPreferences, IBMBackendProvider

ibm_provider = IBMBackendProvider()

preferences = IBMBackendPreferences(
    backend_name="ibmq_manila",
    access_token="insert_token_number",
    provider=ibm_provider,
)

num_trials = 1  # number of times to run qv circuit for each number of qubits
num_shots = 10  # number of runs for each execution
min_qubits = 2
max_qubits = 3

# qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

from classiq.execution import IBMBackendPreferences, IBMBackendProvider ibm_provider = IBMBackendProvider() preferences = IBMBackendPreferences( backend_name="ibmq_manila", access_token="insert_token_number", provider=ibm_provider, ) num_trials = 1 # number of times to run qv circuit for each number of qubits num_shots = 10 # number of runs for each execution min_qubits = 2 max_qubits = 3 # qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)

References¶

[1] https://arxiv.org/pdf/1811.12926.pdf

[2] How to generate a random unitary matrix by Maris Ozols - http://home.lu.lv/~sd20008/papers/essays/Random%20unitary%20[paper].pdf

[3] Computer resources from IBM official website - https://quantum-computing.ibm.com/services/resources?tab=yours