Linq: noisy simulation

This notebook assumes that you are familiar with the basics of this package (Gate, Circuit, Translator, Simulator), and focuses on how users can perform noisy simulation on simulators such as cirq, qiskit or qulacs, using a common interface. Users are free to pick between different compute backends, favoring the one with the best performance of most adequate features, and switch between them with minimal changes to their code.

In particular, we also take the time here to briefly introduce the various quantum error channels supported for the sake of clarity, occasionally sprinkling some (gentle) mathematical equations on top.

Finally, the last section of this notebook delves a bit deeper into the strengths and weaknesses of the different backends supported, in order to give users some insight about what simulator may work best for their experiments, and yield the best performance.

Table of contents

• 1. NoiseModel class
• 2. Noisy Simulation
• 2.1 Pauli error
• 2.2 Depolarization error
• 3. About the supported backends

---

Requirements

In order to run the contents of this notebook, I simply recommend that you install tangelo as per the installation instructions, as well as one or several backends supporting noisy simulation (cirq, qiskit, qulacs…). You can install them with pip. Executing the cell below provides a quick way of installing Tangelo, as well as the minimal set of backends to execute code cells in this notebook.

# Installation of tangelo if not already installed.
try:
    import tangelo
except ModuleNotFoundError:
    !pip install git+https://github.com/goodchemistryco/Tangelo.git@develop  --quiet
    !pip install qiskit qulacs amazon-braket-sdk --quiet

1. NoiseModel class

The noisy_simulation module is home to the NoiseModel class, which provides a generic way to define a noise model that can later on be passed to the Simulator to perform the noisy simulation of a Circuit object on backends supporting that feature.

A NoiseModel object contains a dictionary that maps abstract gates to a list of noises and their parameters, which are usually probabilities. The type of noise currently available can be reached using the SUPPORTED_NOISE_MODELS variable, and the documentation of the NoiseModel class explains what parameters are valid for each of them. The noisy_gates method allows users to quickly check the set of gates included in the noise model.

This representation tells the simulation backend that whenever one of the gates listed in this dictionary is encountered at runtime, the corresponding additional noise(s) should be applied to the quantum state afterwards.

Below, we show the available noise models as well as how users can instantiate a NoiseModel object and then specify noise for any gate relevant to them (here, a Pauli noise with uniform probability distribution over X, Y and Z) using the add_quantum_error method.

from tangelo.linq.noisy_simulation import NoiseModel, SUPPORTED_NOISE_MODELS

print(f'Supported noise models: {SUPPORTED_NOISE_MODELS}\n')

nmp = NoiseModel()
nmp.add_quantum_error("X", 'pauli', [1/3, 1/3, 1/3])
nmp.add_quantum_error("CNOT", 'pauli', [1/3, 1/3, 1/3])

print(f'Noisy gates: {nmp.noisy_gates}\n{nmp}')

Supported noise models: {'pauli', 'depol'}

Noisy gates: {'CNOT', 'X'}
{'X': [('pauli', [0.3333333333333333, 0.3333333333333333, 0.3333333333333333])], 'CNOT': [('pauli', [0.3333333333333333, 0.3333333333333333, 0.3333333333333333])]}

As you can see above, attempting to print a NoiseModel object prints the underlying dictionary. Each noisy gate is associated with a list of 2-tuples, corresponding to the noise model and its parameters.

Using add_quantum_error several times on the same gate allows users to specify potentially complex noise models, cumulating different types of error channels. Combined with variables spanning sets of gates and the Python syntax, we can quickly express rather complete noise models, such as the arbitrary one below:

from tangelo.linq import ONE_QUBIT_GATES, TWO_QUBIT_GATES

nm = NoiseModel()

for g in ONE_QUBIT_GATES:
    nm.add_quantum_error(g, 'pauli', [1/4]*3)
nm.add_quantum_error('CNOT', 'pauli', [1., 0., 0.])
nm.add_quantum_error('CNOT', 'depol', 0.72)
print(nm)

{'PHASE': [('pauli', [0.25, 0.25, 0.25])], 'T': [('pauli', [0.25, 0.25, 0.25])], 'RX': [('pauli', [0.25, 0.25, 0.25])], 'S': [('pauli', [0.25, 0.25, 0.25])], 'Z': [('pauli', [0.25, 0.25, 0.25])], 'Y': [('pauli', [0.25, 0.25, 0.25])], 'RY': [('pauli', [0.25, 0.25, 0.25])], 'H': [('pauli', [0.25, 0.25, 0.25])], 'X': [('pauli', [0.25, 0.25, 0.25])], 'RZ': [('pauli', [0.25, 0.25, 0.25])], 'CNOT': [('pauli', [1.0, 0.0, 0.0]), ('depol', 0.72)]}

Note: You can directly access the underlying dictionary of the NoiseModel object through its private member _quantum_errors if you wish to perform some “surgery” on it directly, in order to modify an existing object. We however recommend to use the add_quantum_error whenever possible to avoid bad surprises and subtle bugs.

Note: Attempting to cumulate quantum error channels of same type (ex: two depol error channels) will return an error.

2. Noisy Simulation

Once a noise model has been created, it can be passed to a simulator (preferably at instantiation). Not all simulators can perform noisy simulation: check out the backend_info dictionary to see which targets shows noisy_simulation set to True. Qiskit and qulacs should be among the available options.

from tangelo.linq import backend_info
from pprint import pprint 

pprint(backend_info)

{'cirq': {'noisy_simulation': True,
          'statevector_available': True,
          'statevector_order': 'lsq_first'},
 'qdk': {'noisy_simulation': False,
         'statevector_available': False,
         'statevector_order': None},
 'qiskit': {'noisy_simulation': True,
            'statevector_available': True,
            'statevector_order': 'msq_first'},
 'qulacs': {'noisy_simulation': True,
            'statevector_available': True,
            'statevector_order': 'msq_first'}}

You can then instantiate a Simulator object for the target of your choice, and provide both the number of shots and noise model required. We recommend to set these parameters at instantiation and avoid modifying them later.

Assuming you have defined an abstract circuit (e.g an object of the Circuit class), the simulate method can then be used to perform the noisy simulation for the desired target and number of shots. It returns a 2-tuple whose second entry is None, which we will disregard as a noisy simulation does not return a statevector.

from tangelo.linq import get_backend, Gate, Circuit

# Define a noise model
nmp = NoiseModel()
nmp.add_quantum_error("X", 'pauli', [1/3, 1/3, 1/3])

# Instantiate a simulator supporting noisy simulation
s_nmp = get_backend(target='qulacs', n_shots=10**6, noise_model=nmp)

# Define a Circuit
c1 = Circuit([Gate("X",0)])

# Simulate and look at the results
res_pauli1, _ = s_nmp.simulate(c1)
print(res_pauli1)

{'0': 0.665706, '1': 0.334294}

2.1 Pauli error

A Pauli error channel can be defined by an arbitrary probability distribution \([p_X, p_Y, p_Z]\) over the 3 Pauli operators. This channel randomly applies one of the X, Y, Z Pauli operators after a noisy gate has been encountered, by sampling from this probability distribution. Implicitly, the identity is applied with probability \(1- (p_X + p_Y + p_Z)\), and no noise is added in that case.

In the example above, our circuit merely consists of a single X gate applied to a single qubit, initially in state \(|0\rangle\). Our noise model always randomly applies, with equal probability, an extra X, Y or Z gate to the resulting quantum state.

As the number of shots grows, the results of the above should converge towards \(|1\rangle\) (resp \(|0\rangle\)) observed with probability \(1/3\) (resp \(2/3\)). The theoretical verification is pretty straightforward:

• After X has been applied to our qubit, our system is in state \(|1\rangle\)
• With probability \(2/3\), an extra X or Y gate is applied, and our system is sure to be measured in state \(|0\rangle\)
• With probability \(1/3\), an extra Z gate is applied, and our system is sure to be measured in state \(|1\rangle\)

2.2 Depolarization error

Depolarizing quantum error channels are also available. The channel can be defined as:

\[E_{\lambda, n}(\rho) = (1 - \lambda)\rho + \lambda Tr[\rho] \frac{I}{2^{n}}\]

where:

• \(n\) is the number of qubits in the system
• \(\rho\) is the density matrix representing the quantum state of the system
• \(\lambda\) is the parameter of the depolarizing channel, and is such that \(0 \le \lambda \le \frac{4^n}{4^n-1}\)

In particular:

• If \(\lambda=0\), \(E\) is the identity channel (no noise)
• If \(\lambda=1\), \(E(\rho)=\frac{I}{2^{n}}\)
• If \(\lambda=\frac{4^n}{4^n-1}\), \(E\) is a uniform Pauli error channel

The cell below shows that calculations for our “X-gate circuit” align with this formula. In particular, the case where \(\lambda=\frac{4}{3}\) yields the same error channel than the Pauli error channel introduced in the previous section.

for l in [0., 1., 4/3]:
    nmd = NoiseModel()
    nmd.add_quantum_error("X", 'depol', l)
    s_nmd = get_backend(target='qulacs', n_shots=10**6, noise_model=nmd)
    res_depol1, _ = s_nmd.simulate(c1) # This circuit only contains a single X gate
    print(f'lambda = {l:.2f} :: {res_depol1}')

lambda = 0.00 :: {'1': 1.0}
lambda = 1.00 :: {'0': 0.500183, '1': 0.499817}
lambda = 1.33 :: {'0': 0.666488, '1': 0.333512}

Similarly, we verify that in the example below, the results coincide with the mathematical formula above, and yield \(\frac{1}{2}\rho + \frac{1}{8}I\), \(\rho\) being the density matrix of the system, still left in state \(|00\rangle\) after the CNOT gate has been applied.

# A circuit with a single CNOT gate
c2 = Circuit([Gate("CNOT", target=1, control=0)])

nmd = NoiseModel()
nmd.add_quantum_error("CNOT", 'depol', 0.5)
s_nmd = get_backend(target='qulacs', n_shots=10**6, noise_model=nmd)
res_depol1, _ = s_nmd.simulate(c2)
print(res_depol1)

{'11': 0.124993, '00': 0.624951, '01': 0.125071, '10': 0.124985}

3. About the supported backends

get_backend provides a uniform interface to several noisy simulators, such as qiskit and qulacs, hiding the complexity of the noisy simulation process from the user. The inner workings are a bit different depending on the simulator selected, and the simulators themselves have their own strengths and weaknesses, which we gloss over in order to provide you with some insight about what tool may help you the most in your developments.

Qiskit [This section is outdated and needs to be updated]

qiskit already has their own NoiseModel class. Our package simply translates our NoiseModel object into the qiskit one at runtime, which is then passed to their shot-based QASM simulator, in charge of transpiling and applying this noise model to the circuit. The package exposes this function as get_qiskit_noise_model, but you should ideally let the backend object take care of it. I recommend you import qiskit’s NoiseModel as QiskitNoiseModel to avoid insidious bugs, if you ever need both in the same script.

Note: The rest of this subsection may not hold for long, as a pull-request introducing support for rx, ry and rz gates has been merged into the Qiskit repository, and may be avilable through the pip package soon (Oct 2020?)

Qiskit can be restrictive, as it enforces transpilation of your circuit into the gate set supported by their qasm simulator. For instance, the qasm simulator gate set does not include the rx,ry,rz rotation gates, but instead converts them at runtime into u1,u2 oru3 gates (reference), which are IBM’s gates. This has several implications:

• Users may not “get what they ask for”: there is a loss of transparency due to what’s going on under the hood, with the transpiling step.

• In particular all rotations gates may map into u3 gates ar runtime, since it can apply any single-qubit unitary with the right parameters. This means that specifying different noise channel for the different rotation gates (which may make sense according to some researchers) is ill-defined: what noise should then be applied to u3 ? Our package limits cumulating quantum error channels of different natures only to avoid stacking redundant noise to mitigate the issue, but it is better to probably switch to another backend to stay true to your research.

from tangelo.linq.noisy_simulation import get_qiskit_noise_model

nm = NoiseModel()
nm.add_quantum_error('RX', 'depol', 4/3)
nm.add_quantum_error('RX', 'pauli', [1.,0.,0.])
qnm = get_qiskit_noise_model(nm)

# Any rotation gate will have its quantum error channel applied to u1, u2, u3 as we do not know
# which one the transpiler will generate and need to apply noise to. 
# Qiskit always shows a warning when several noise channels are composed onto the same basis gate
print(qnm)

WARNING: all-qubit error already exists for instruction "u1", composing with additional error.
WARNING: all-qubit error already exists for instruction "u2", composing with additional error.
WARNING: all-qubit error already exists for instruction "u3", composing with additional error.

NoiseModel:
  Basis gates: ['cx', 'id', 'rz', 'sx', 'u1', 'u2', 'u3']
  Instructions with noise: ['u2', 'u1', 'u3']
  All-qubits errors: ['u1', 'u2', 'u3']

Qulacs

qulacs does not have built-in noise model objects, but has noisy gates instead, which means users can build their own noisy circuit by adding these noisy gates when relevant. That noisy circuit is then simulated by the qulacs C++ engine. This very flexible approach, coupled to the fact the qulacs simulator does not enforce any specific basis gate set (it even supports the IBM u gates), lets you do pretty much whatever you want when it comes to noise.

The translate_circuit function, taking in an abstract circuit and returning a qulacs circuit thus supports an extra arguments compared to the other one: it can take in a NoiseModel object, and will return the noisy circuit. The Simulator is tasked with calling this translation function with the noise model, under the hood. The noisy circuit is however of limited interest to the user to have access to this circuit, as the noisy gates show up as “generic gates” and very little information is available through Python (the underlying qulacs code being written in C++): all you know is basically where noisy gates have been added.

If you’d like to know more about what these noisy gates are, please have a look at the following: - Brief look at Probabilistic gates - Code listing: Noise functions

cirq

Much like qulacs, cirq does not have built-in noise model objects, but has noisy gates instead, which means users can build their own noisy circuit by adding these noisy gates when relevant. That noisy circuit is then simulated by cirq. The translate_circuit function takes a NoiseModel object to implement this approach, returning the noisy circuit.

One further note, the DensityMatrixSimulator is used for cirq simulation of noisy or mixed state circuits.

If you’d like to know more about what these noisy gates are, please have a look at the following: - Guide: Noise - Code listing: All Gates

Performance

When it comes to noisy simulation, we noticed that the performance of the different backends we support was not as we expected from publicly available benchmarks, which usually focus on noiseless simulation. This may be explained by both the difference in the approaches taken by these backends, as well as the quality of our wrappers. Performance also depends on:

• the size of your circuit (number of qubits, depth)
• your hardware and your environment (OS, number of OpenMP threads, virtual machine or not…)

We recommend proceeding on a case-by-case basis. Below, a quick example on our both of them perform on a circuit repeating CNOT ladders a few times. Feel free to change some parameters.

Note: Jupyter notebooks may be considerably slower than running Python scripts directly, depending on your environment. They may not be the most reliable when it comes to making a choice. In the future, there may be a benchmark folder available in this package, which could help you figure out quickly what is best in your environment and how fast things can go.

import time
import numpy as np

theta = 2.1
n_qubits = 5
n_repetitions = 3

def ladder_gates(n_qubits, theta, n_repetitions):
    gates = [Gate("H", i) for i in range(n_qubits)]
    gates += [Gate("CNOT", target=i+1, control=i) for i in range (n_qubits-1)]
    gates += [Gate("RY", n_qubits-1, parameter=theta)]
    gates += [Gate("CNOT", target=i+1, control=i) for i in range (n_qubits-1)][::-1]
    return gates * n_repetitions

# CNOT ladder circuit, a common pattern in quantum chemistry ansatz circuits
c_base = Circuit(ladder_gates(n_qubits, theta, n_repetitions))
print(c_base)

# Define noise model
nmb = NoiseModel()
nmb.add_quantum_error('RY', 'pauli', [.5, 0.5, 0.])
nmb.add_quantum_error('CNOT', 'depol', 0.2)
print(nmb,'\n')

# Elapsed time may differ greatly depending on your machine and development environment
for backend in {'qulacs', 'qiskit', 'cirq'}:
    sim = get_backend(backend, n_shots=10**6, noise_model=nmb)
    tstart = time.time()
    res, _ = sim.simulate(c_base)
    print(f'{backend}:: {res} \t (elapsed: {time.time()-tstart:.1f} s)\n')

Circuit object. Size 42 

H         target : [0]   
H         target : [1]   
H         target : [2]   
H         target : [3]   
H         target : [4]   
CNOT      target : [1]   control : [0]   
CNOT      target : [2]   control : [1]   
CNOT      target : [3]   control : [2]   
CNOT      target : [4]   control : [3]   
RY        target : [4]   parameter : 2.1
CNOT      target : [4]   control : [3]   
CNOT      target : [3]   control : [2]   
CNOT      target : [2]   control : [1]   
CNOT      target : [1]   control : [0]   
H         target : [0]   
H         target : [1]   
H         target : [2]   
H         target : [3]   
H         target : [4]   
CNOT      target : [1]   control : [0]   
CNOT      target : [2]   control : [1]   
CNOT      target : [3]   control : [2]   
CNOT      target : [4]   control : [3]   
RY        target : [4]   parameter : 2.1
CNOT      target : [4]   control : [3]   
CNOT      target : [3]   control : [2]   
CNOT      target : [2]   control : [1]   
CNOT      target : [1]   control : [0]   
H         target : [0]   
H         target : [1]   
H         target : [2]   
H         target : [3]   
H         target : [4]   
CNOT      target : [1]   control : [0]   
CNOT      target : [2]   control : [1]   
CNOT      target : [3]   control : [2]   
CNOT      target : [4]   control : [3]   
RY        target : [4]   parameter : 2.1
CNOT      target : [4]   control : [3]   
CNOT      target : [3]   control : [2]   
CNOT      target : [2]   control : [1]   
CNOT      target : [1]   control : [0]   

{'RY': [('pauli', [0.5, 0.5, 0.0])], 'CNOT': [('depol', 0.2)]} 

cirq:: {'10110': 0.031467, '11001': 0.03113, '00000': 0.031401, '01100': 0.031785, '10100': 0.031793, '11111': 0.030839, '01111': 0.031061, '10101': 0.03082, '00001': 0.031185, '11100': 0.03136, '11010': 0.031607, '01010': 0.031367, '01001': 0.031188, '00011': 0.030849, '00101': 0.03113, '11110': 0.031378, '00010': 0.031566, '10000': 0.03142, '11000': 0.031563, '00111': 0.030893, '10001': 0.031271, '00100': 0.031254, '01101': 0.030801, '11011': 0.03114, '01011': 0.03132, '01000': 0.031301, '10011': 0.031099, '10111': 0.031005, '01110': 0.031418, '11101': 0.030834, '10010': 0.031388, '00110': 0.031367}     (elapsed: 5.5 s)

qiskit:: {'10111': 0.029946, '01100': 0.030131, '00000': 0.030416, '11110': 0.030231, '01011': 0.032338, '10000': 0.032232, '11000': 0.030366, '01001': 0.030344, '11010': 0.032555, '10011': 0.032224, '00001': 0.031977, '10101': 0.032195, '10100': 0.030172, '11001': 0.032078, '01110': 0.032224, '00011': 0.030167, '10110': 0.032283, '11101': 0.030029, '00111': 0.032251, '01010': 0.030193, '01101': 0.032227, '11111': 0.032376, '01000': 0.032486, '11011': 0.030211, '10010': 0.030649, '00100': 0.032344, '00110': 0.030454, '10001': 0.030305, '00010': 0.032109, '01111': 0.03028, '00101': 0.030018, '11100': 0.032189}     (elapsed: 2.1 s)

qulacs:: {'10110': 0.031455, '00100': 0.031153, '11111': 0.030975, '10001': 0.031072, '11100': 0.031301, '01110': 0.03146, '00111': 0.031039, '01001': 0.030899, '01011': 0.030996, '11101': 0.03128, '11010': 0.03132, '00110': 0.031647, '11110': 0.031263, '10010': 0.03098, '10101': 0.031471, '01010': 0.031262, '11011': 0.031197, '11000': 0.031448, '00011': 0.03094, '10011': 0.031287, '00010': 0.031363, '01111': 0.031272, '00001': 0.03129, '10100': 0.03135, '01100': 0.031506, '01101': 0.031216, '01000': 0.031404, '00101': 0.030979, '00000': 0.031368, '10111': 0.03109, '11001': 0.031172, '10000': 0.031545}    (elapsed: 7.5 s)

---

Valentin Senicourt