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 hardware backends.
The scheme of the quantum volume [1]: 1. For a number of qubits \(n\), a circuit is made of the \(n\) quantum layer. 2. Each layer consists of a unitary operation between pairs of \(n\) qubits. The pairs are chosen at random. If \(n\) is odd, one of them does not have an operation. 3. The unitary operation between each pair is the Haar random matrix; i.e., SU(4) operation containing a random complex number in such a manner that the probability of measuring a quantum state is kept with uniform distribution. 4. A single circuit of \(n\) qubits is measured and the heavy output probability (i.e., the probability of measuring 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 subsides to the following inequality: (1) P{heavy_outputs} >= 2/3 . 5. For a given output, to get the quantum volume, repeat Items 1-4 for an increasing number of qubits until the inequality described in Item 4 does not hold. To ensure it, the circuits are created many times and the average and standard deviation are taken into account. 6. The quantum volume is 2 to the power of the number of qubits, such that they pass inequality (1) as per the procedure described in Items 1-5.
The heavy output probability is a good measurement of the quality of the circuit, as noise reduces the probabilities of uniform distribution. While this is so, consider that there are many components to the results of the procedureβnot only the hardware noises, but also the connectivity map, the quality of the transpilation, and even the quantum software that translates the circuit into basis gates for the hardware, thus contributing to the circuit depth.
This demonstration shows the code for implementing the steps to calculate the quantum volume using the Classiq platform and an example of such calculations for several quantum simulators and hardware backends.
Step 1: Create a Haar Random Unitary Matrix
Create a function, generating a (n,n) sized Haar random unitary matrix [2]. This matrix contains a random complex number that is distributed evenly in the \(2^n\) space of quantum states. The Haar distribution indicates how to weight the elements of \(π(π)\) such that uniform distribution occurs in the parameter space.
import numpy as np
from numpy.linalg import qr
from scipy.stats import unitary_group
def haar(n):
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 a Quantum Volume Circuit
The qv_model function creates the quantum volume model for a given \(N\) number of qubits. For \(N\) qubits, the circuit must include \(N\) quantum volume layers. The layers are built using the qv_layer function, which creates random pairing between the \(N\) qubits. (For an odd number, a randomly chosen qubit is not operational). Between each pair, a unitary gate operates, consisting of a Haar random unitary matrix of size 4.
import math
import random
from classiq.qmod import (
CArray,
CReal,
Output,
QArray,
QBit,
allocate,
bind,
create_model,
qfunc,
unitary,
)
@qfunc
def apply_2qbit_unitary(unitary_params: CArray[CArray[CReal]], x: QBit, y: QBit):
joined = QArray[QBit]("joined")
bind([x, y], joined)
unitary(unitary_params, joined)
bind(joined, [x, y])
# We will use flat-python in order to generate the random unitaries.
# That is why you don't see the @qfunc decorator here
def qv_layer(target: QArray[QBit], N):
# Step 1: start with a shuffle of the qubits
qubit_list = list(range(N))
random.shuffle(qubit_list)
for idx in range(math.floor(N / 2)):
# Step 2: Isolate the qubit pairs for the layers
a = qubit_list[idx]
b = qubit_list[math.floor(N / 2) + idx]
# Step 3: Generate the random matrix (this needs to change for the random matrix when possible)
gate_matrix = haar(4).tolist()
apply_2qbit_unitary(gate_matrix, target[a], target[b])
def qv_model(N):
@qfunc
def main(target: Output[QArray[QBit]]):
allocate(N, target)
for idx in range(N):
qv_layer(target, N)
return create_model(main)
Step 3: Execute and Analyze
The execution and analysis part consists of these functions:
* execute_qv sends a quantum program for execution on a given quantum hardware with a specified number of shots. The 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.
from classiq import execute, set_quantum_program_execution_preferences
from classiq.execution import (
AzureBackendPreferences,
ClassiqBackendPreferences,
ClassiqSimulatorBackendNames,
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
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
from math import floor, log10
def round_significant(x):
return round(x, -int(floor(log10(abs(x)))))
Step 4: Find the 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.
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
Examples
Run the code to find the quantum volume of several quantum simulators and hardwares.
Running with Classiq's Simulator
num_trials = 10 # number of times to run the 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=ClassiqSimulatorBackendNames.SIMULATOR
)
min_qubits = 3
max_qubits = 6
qv = find_qv(num_trials, num_shots, min_qubits, max_qubits, preferences)
0%| | 0/4 [00:00<?, ?it/s]
0%| | 0/10 [00:00<?, ?it/s]
[A
10%|β | 1/10 [00:03<00:34, 3.81s/it]
[A
20%|ββ | 2/10 [00:07<00:29, 3.63s/it]
[A
30%|βββ | 3/10 [00:10<00:24, 3.57s/it]
[A
40%|ββββ | 4/10 [00:14<00:21, 3.58s/it]
[A
50%|βββββ | 5/10 [00:18<00:18, 3.62s/it]
[A
60%|ββββββ | 6/10 [00:21<00:14, 3.63s/it]
[A
70%|βββββββ | 7/10 [00:25<00:10, 3.63s/it]
[A
80%|ββββββββ | 8/10 [00:28<00:07, 3.62s/it]
[A
90%|βββββββββ | 9/10 [00:31<00:03, 3.29s/it]
[A
100%|ββββββββββ| 10/10 [00:35<00:00, 3.46s/it]
[A
100%|ββββββββββ| 10/10 [00:35<00:00, 3.54s/it]
25%|βββ | 1/4 [00:35<01:46, 35.40s/it]
for 3 qubits the heavy outputs probability is: 0.811 with 0.008 standard deviation
0%| | 0/10 [00:00<?, ?it/s]
[A
10%|β | 1/10 [00:03<00:32, 3.57s/it]
[A
20%|ββ | 2/10 [00:07<00:28, 3.55s/it]
[A
30%|βββ | 3/10 [00:10<00:25, 3.58s/it]
[A
40%|ββββ | 4/10 [00:15<00:23, 3.93s/it]
[A
50%|βββββ | 5/10 [00:19<00:20, 4.19s/it]
[A
60%|ββββββ | 6/10 [00:23<00:16, 4.02s/it]
[A
70%|βββββββ | 7/10 [00:27<00:11, 3.94s/it]
[A
80%|ββββββββ | 8/10 [00:32<00:08, 4.19s/it]
[A
90%|βββββββββ | 9/10 [00:35<00:03, 4.00s/it]
[A
100%|ββββββββββ| 10/10 [00:39<00:00, 3.85s/it]
[A
100%|ββββββββββ| 10/10 [00:39<00:00, 3.91s/it]
50%|βββββ | 2/4 [01:14<01:15, 37.60s/it]
for 4 qubits the heavy outputs probability is: 0.846 with 0.004 standard deviation
0%| | 0/10 [00:00<?, ?it/s]
[A
10%|β | 1/10 [00:04<00:42, 4.72s/it]
[A
20%|ββ | 2/10 [00:08<00:33, 4.13s/it]
[A
30%|βββ | 3/10 [00:12<00:27, 3.91s/it]
[A
40%|ββββ | 4/10 [00:15<00:22, 3.82s/it]
[A
50%|βββββ | 5/10 [00:20<00:20, 4.09s/it]
[A
60%|ββββββ | 6/10 [00:25<00:17, 4.31s/it]
[A
70%|βββββββ | 7/10 [00:28<00:12, 4.12s/it]
[A
80%|ββββββββ | 8/10 [00:33<00:08, 4.28s/it]
[A
90%|βββββββββ | 9/10 [00:38<00:04, 4.39s/it]
[A
100%|ββββββββββ| 10/10 [00:41<00:00, 4.21s/it]
[A
100%|ββββββββββ| 10/10 [00:41<00:00, 4.18s/it]
75%|ββββββββ | 3/4 [01:56<00:39, 39.54s/it]
for 5 qubits the heavy outputs probability is: 0.835 with 0.005 standard deviation
0%| | 0/10 [00:00<?, ?it/s]
[A
10%|β | 1/10 [00:04<00:44, 4.90s/it]
[A
20%|ββ | 2/10 [00:09<00:38, 4.84s/it]
[A
30%|βββ | 3/10 [00:14<00:33, 4.85s/it]
[A
40%|ββββ | 4/10 [00:19<00:28, 4.80s/it]
[A
50%|βββββ | 5/10 [00:23<00:23, 4.76s/it]
[A
60%|ββββββ | 6/10 [00:28<00:19, 4.77s/it]
[A
70%|βββββββ | 7/10 [00:33<00:14, 4.76s/it]
[A
80%|ββββββββ | 8/10 [00:38<00:09, 4.96s/it]
[A
90%|βββββββββ | 9/10 [00:43<00:04, 4.89s/it]
[A
100%|ββββββββββ| 10/10 [00:48<00:00, 4.78s/it]
[A
100%|ββββββββββ| 10/10 [00:48<00:00, 4.82s/it]
100%|ββββββββββ| 4/4 [02:44<00:00, 42.95s/it]
100%|ββββββββββ| 4/4 [02:44<00:00, 41.14s/it]
for 6 qubits the heavy outputs probability is: 0.861 with 0.006 standard deviation
##### The quantum volume is 64 #####
Since this is a simulator with no errors, we expect the heavy output probability for any number of qubits to be approximately 0.85
Running with Rigetti Aspen M-3
num_trials = 10 # number of times to run the 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)
Running with IBM Machines
Try to run a few IBM machines: ibmq_lima with reported quantum volume of 8; ibmq_quito with reported quantum volume of 16; and ibmq_manila with reported quantum volume of 32 [3].
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 the 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_quito",
access_token="insert_token_number",
provider=ibm_provider,
)
num_trials = 1 # number of times to run the 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 the 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 the official IBM website: https://quantum-computing.ibm.com/services/resources?tab=yours.