Advanced examples#

Here are a few short advanced how to examples.

How to select hardware devices?#

Qibo supports execution on different hardware configurations including CPU with multi-threading, single GPU and multiple GPUs. Here we provide some useful information on how to control the devices that Qibo uses for circuit execution in order to maximize performance for the available hardware configuration.

Switching between CPU and GPU#

If a GPU with CUDA support is available in the system and tensorflow or qibojit (cupy) are installed then circuits will be executed on the GPU automatically unless the user specifies otherwise. One can change the default simulation device using qibo.set_device:

import qibo
qibo.set_device("/CPU:0")
final_state = c() # circuit will now be executed on CPU

The syntax of device names follows the pattern '/{device type}:{device number}' where device type can be CPU or GPU and the device number is an integer that distinguishes multiple devices of the same type starting from 0. For more details we refer to Tensorflow’s tutorial on manual device placement. Alternatively, running the command CUDA_VISIBLE_DEVICES="" in a terminal hides CUDA GPUs from this terminal session.

In most cases the GPU accelerates execution compared to CPU, however the following limitations should be noted:

  • For small circuits (less than 10 qubits) the overhead from casting tensors on GPU may be larger than executing the circuit on CPU, making CPU execution preferable. In such cases disabling CPU multi-threading may also increase performance (see next subsection).

  • A standard GPU has 12-16GB of memory and thus can simulate up to 30 qubits on single-precision or 29 qubits with double-precision when Qibo’s default gates are used. For larger circuits one should either use the CPU (assuming it has more memory) or a distributed circuit configuration. The latter allows splitting the state vector on multiple devices and is useful both when multiple GPUs are available in the system or even for re-using a single GPU (see relevant subsection bellow).

Note that if the used device runs out of memory during a circuit execution an error will be raised prompting the user to switch the default device using qibo.set_device.

Setting the number of CPU threads#

Qibo by default uses the qibojit backend which is based on custom operators. This backend uses OpenMP instructions for parallelization. In most cases, utilizing all available CPU threads provides better performance. However, for small circuits the parallelization overhead may decrease performance making single thread execution preferrable.

You can restrict the number of threads by exporting the OMP_NUM_THREADS environment variable with the requested number of threads before launching Qibo, or programmatically, during runtime, as follows:

import qibo
# set the number of threads to 1
qibo.set_threads(1)
# retrieve the current number of threads
current_threads = qibo.get_threads()

On the other hand, when using the tensorflow backend Qibo inherits Tensorflow’s defaults for CPU thread configuration. Tensorflow allows restricting the number of threads as follows:

import tensorflow as tf
tf.config.threading.set_inter_op_parallelism_threads(1)
tf.config.threading.set_intra_op_parallelism_threads(1)
import qibo

Note that this should be run during Tensorflow initialization in the beginning of the script and before creating the qibo backend.

Using multiple GPUs#

Qibo supports distributed circuit execution on multiple GPUs. This feature can be used as follows:

from qibo import Circuit, gates

# Define GPU configuration
accelerators = {"/GPU:0": 3, "/GPU:1": 1}
# this will use the first GPU three times and the second one time
# leading to four total logical devices
# construct the distributed circuit for 32 qubits
c = Circuit(32, accelerators)

Gates can then be added normally using c.add and the circuit can be executed using c(). Note that a memory_device is passed in the distributed circuit (if this is not passed the CPU will be used by default). This device does not perform any gate calculations but is used to store the full state. Therefore the distributed simulation is limited by the amount of CPU memory.

Also, note that it is possible to reuse a single GPU multiple times increasing the number of “logical” devices in the distributed calculation. This allows users to execute circuits with more than 30 qubits on a single GPU by reusing several times using accelerators = {"/GPU:0": ndevices}. Such a simulation will be limited by CPU memory only.

For systems without GPUs, the distributed implementation can be used with any type of device. For example if multiple CPUs, the user can pass these CPUs in the accelerator dictionary.

Distributed circuits are generally slower than using a single GPU due to communication bottleneck. However for more than 30 qubits (which do not fit in single GPU) and specific applications (such as the QFT) the multi-GPU scheme can be faster than using only CPU.

Note that simulating a circuit using multiple GPUs partitions the state in multiple pieces which are distributed to the different devices. Creating the full state as a single tensor would require merging these pieces and using twice as much memory. This is disabled by default, however the user may create the full state as follows:

# Create distributed circuits for two GPUs
c = Circuit(32, {"/GPU:0": 1, "/GPU:1": 1})
# Add gates
c.add(...)
# Execute (``result`` will be a ``DistributedState``)
result = c()

# ``DistributedState`` supports indexing and slicing
print(result[40])
# will print the 40th component of the final state vector
print(result[20:25])
# will print the components from 20 to 24 (inclusive)

# Access the full state (will double memory usage)
final_state = result.state()
# ``final_state`` is a ``tf.Tensor``

How to use callbacks?#

Callbacks allow the user to apply additional functions on the state vector during circuit execution. An example use case of this is the calculation of entanglement entropy as the state propagates through a circuit. This can be implemented easily using qibo.callbacks.EntanglementEntropy and the qibo.gates.CallbackGate gate. For example:

from qibo import models, gates, callbacks

# create entropy callback where qubit 0 is the first subsystem
entropy = callbacks.EntanglementEntropy([0])

# initialize circuit with 2 qubits and add gates
c = models.Circuit(2) # state is |00> (entropy = 0)
c.add(gates.CallbackGate(entropy)) # performs entropy calculation in the initial state
c.add(gates.H(0)) # state is |+0> (entropy = 0)
c.add(gates.CallbackGate(entropy)) # performs entropy calculation after H
c.add(gates.CNOT(0, 1)) # state is |00> + |11> (entropy = 1))
c.add(gates.CallbackGate(entropy)) # performs entropy calculation after CNOT

# execute the circuit using the callback
final_state = c()

The results can be accessed using indexing on the callback objects. In this example entropy[:] will return [0, 0, 1] which are the values of entropy after every gate in the circuit.

The same callback object can be used in a second execution of this or a different circuit. For example

# c is the same circuit as above
# execute the circuit
final_state = c()
# execute the circuit a second time
final_state = c()

# print result
print(entropy[:]) # [0, 0, 1, 0, 0, 1]

The callback for entanglement entropy can also be used on state vectors directly. For example

How to use parametrized gates?#

Some Qibo gates such as rotations accept values for their free parameter. Once such gates are added in a circuit their parameters can be updated using the qibo.models.circuit.Circuit.set_parameters() method. For example:

from qibo import Circuit, gates
# create a circuit with all parameters set to 0.
c = Circuit(3)
c.add(gates.RX(0, theta=0))
c.add(gates.RY(1, theta=0))
c.add(gates.CZ(1, 2))
c.add(gates.fSim(0, 2, theta=0, phi=0))
c.add(gates.H(2))

# set new values to the circuit's parameters
params = [0.123, 0.456, (0.789, 0.321)]
c.set_parameters(params)

initializes a circuit with all gate parameters set to 0 and then updates the values of these parameters according to the params list. Alternatively the user can use circuit.set_parameters() with a dictionary or a flat list. The keys of the dictionary should be references to the gate objects of the circuit. For example:

c = Circuit(3)
g0 = gates.RX(0, theta=0)
g1 = gates.RY(1, theta=0)
g2 = gates.fSim(0, 2, theta=0, phi=0)
c.add([g0, g1, gates.CZ(1, 2), g2, gates.H(2)])

# set new values to the circuit's parameters using a dictionary
params = {g0: 0.123, g1: 0.456, g2: (0.789, 0.321)}
c.set_parameters(params)
# equivalently the parameter's can be update with a list as
params = [0.123, 0.456, (0.789, 0.321)]
c.set_parameters(params)
# or with a flat list as
params = [0.123, 0.456, 0.789, 0.321]
c.set_parameters(params)

If a list is given then its length and elements should be compatible with the parametrized gates contained in the circuit. If a dictionary is given then its keys should be all the parametrized gates in the circuit.

The following gates support parameter setting:

  • RX, RY, RZ, U1, CU1: Accept a single theta parameter.

  • qibo.gates.fSim: Accepts a tuple of two parameters (theta, phi).

  • qibo.gates.GeneralizedfSim: Accepts a tuple of two parameters (unitary, phi). Here unitary should be a unitary matrix given as an array or tf.Tensor of shape (2, 2).

  • qibo.gates.Unitary: Accepts a single unitary parameter. This should be an array or tf.Tensor of shape (2, 2).

Note that a np.ndarray or a tf.Tensor may also be used in the place of a flat list. Using qibo.models.circuit.Circuit.set_parameters() is more efficient than recreating a new circuit with new parameter values. The inverse method qibo.models.circuit.Circuit.get_parameters() is also available and returns a list, dictionary or flat list with the current parameter values of all parametrized gates in the circuit.

It is possible to hide a parametrized gate from the action of qibo.models.circuit.Circuit.get_parameters() and qibo.models.circuit.Circuit.set_parameters() by setting the trainable=False during gate creation. For example:

c = Circuit(3)
c.add(gates.RX(0, theta=0.123))
c.add(gates.RY(1, theta=0.456, trainable=False))
c.add(gates.fSim(0, 2, theta=0.789, phi=0.567))

print(c.get_parameters())
# prints [(0.123,), (0.789, 0.567)] ignoring the parameters of the RY gate

[(0.123,), (0.789, 0.567)]

This is useful when the user wants to freeze the parameters of specific gates during optimization. Note that trainable defaults to True for all parametrized gates.

How to collapse state during measurements?#

As mentioned in the How to perform measurements? measurement can by default be used only in the end of the circuit and they do not have any effect on the state. In this section we describe how to collapse the state during measurements and re-use measured qubits in the circuit. Collapsing the state means projecting to the |0> or |1> subspace according to the sampled result for each measured qubit.

The state is collapsed when the collapse=True is used during instantiation of the qibo.gates.M gate. For example

from qibo import Circuit, gates

c = Circuit(1, density_matrix=True)
c.add(gates.H(0))
output = c.add(gates.M(0, collapse=True))
c.add(gates.H(0))
result = c(nshots=1)
print(result)
# prints |+><+| if 0 is measured
# or |-><-| if 1 is measured

In this example the single qubit is measured while in the state (|0> + |1>) and is collapsed to either |0> or |1>. The qubit can then be re-used by adding more gates that act to this. The outcomes of collapse=True measurements is not contained in the final result object but is accessible from the output object returned when adding the gate to the circuit. output supports the output.samples() and output.frequencies() functionality as described in How to perform measurements?.

Collapse gates are single-shot by default because the state collapse is not well-defined for more than one shots. If the user passes the nshots arguments during the circuit execution (eg. result = c(nshots=100) in the above example), then the circuit execution will be repeated nshots times using a loop:

for _ in range(nshots):
    result = c()

Note that this will be more time-consuming compared to multi-shot simulation of standard (non-collapse) measurements where the circuit is simulated once and the final state vector is sampled nshots times. For multi-shot simulation the outcomes are still accessible using output.samples() and output.frequencies().

Using normal measurements and collapse measurements in the same circuit is also possible:

from qibo import Circuit, gates

c = Circuit(2)
c.add(gates.H(0))
c.add(gates.H(1))
output = c.add(gates.M(0, collapse=True))
c.add(gates.H(0))
c.add(gates.M(0, 1))
result = c(nshots=100)

In this case output will contain the results of the first collapse=True measurement while result will contain the results of the standard measurement.

Conditioning gates on measurement outcomes#

The output of collapse=True measurements can be used as a parameter in any parametrized gate as follows:

import numpy as np
from qibo import Circuit, gates

c = Circuit(2, density_matrix=True)
c.add(gates.H(0))
output = c.add(gates.M(0, collapse=True))
c.add(gates.RX(1, theta=np.pi * output.symbols[0] / 4))
result = c()

In this case the first qubit will be measured and if 1 is found a pi/4 X-rotation will be applied to the second qubit, otherwise no rotation. Qibo allows to use output as a parameter during circuit creation through the use of sympy.Symbol objects. These symbols can be accessed through the output.symbols list and they acquire a numerical value during execution when the measurement is performed. As explained above, if nshots > 1 is given during circuit execution the execution is repeated using a loop.

If more than one qubits are used in a collapse=True measurement gate the output.symbols list can be indexed accordingly:

import numpy as np
from qibo import Circuit, gates

c = Circuit(3, density_matrix=True)
c.add(gates.H(0))
output = c.add(gates.M(0, 1, collapse=True))
c.add(gates.RX(1, theta=np.pi * output.symbols[0] / 4))
c.add(gates.RY(2, theta=np.pi * (output.symbols[0] + output.symbols[1]) / 5))
result = c()

How to invert a circuit?#

Many quantum algorithms require using a specific subroutine and its inverse in the same circuit. Qibo simplifies this implementation via the qibo.models.circuit.Circuit.invert() method. This method produces the inverse of a circuit by taking the dagger of all gates in reverse order. It can be used with circuit addition to simplify the construction of algorithms, for example:

from qibo import Circuit, gates

# Create a subroutine
subroutine = Circuit(6)
subroutine.add([gates.RX(i, theta=0.1) for i in range(5)])
subroutine.add([gates.CZ(i, i + 1) for i in range(0, 5, 2)])

# Create the middle part of the circuit
middle = Circuit(6)
middle.add([gates.CU2(i, i + 1, phi=0.1, lam=0.2) for i in range(0, 5, 2)])

# Create the total circuit as subroutine + middle + subroutine^{-1}
circuit = subroutine + middle + subroutine.invert()

Note that circuit addition works only between circuits that act on the same number of qubits. It is often useful to add subroutines only on a subset of qubits of the large circuit. This is possible using the qibo.models.circuit.Circuit.on_qubits() method. For example:

from qibo import models, gates

# Create a small circuit of 4 qubits
smallc = models.Circuit(4)
smallc.add((gates.RX(i, theta=0.1) for i in range(4)))
smallc.add((gates.CNOT(0, 1), gates.CNOT(2, 3)))

# Create a large circuit on 8 qubits
largec = models.Circuit(8)
# Add the small circuit on even qubits
largec.add(smallc.on_qubits(*range(0, 8, 2)))
# Add a QFT on odd qubits
largec.add(models.QFT(4).on_qubits(*range(1, 8, 2)))
# Add an inverse QFT on first 6 qubits
largec.add(models.QFT(6).invert().on_qubits(*range(6)))

How to write a VQE?#

The VQE requires an ansatz function and a Hamiltonian object. There are examples of VQE optimization in examples/benchmarks:

  • vqe.py: a simple example with the XXZ model.

Here is a simple example using the Heisenberg XXZ model Hamiltonian:

import numpy as np
from qibo import models, gates, hamiltonians

nqubits = 6
nlayers  = 4

# Create variational circuit
circuit = models.Circuit(nqubits)
for l in range(nlayers):
    circuit.add((gates.RY(q, theta=0) for q in range(nqubits)))
    circuit.add((gates.CZ(q, q+1) for q in range(0, nqubits-1, 2)))
    circuit.add((gates.RY(q, theta=0) for q in range(nqubits)))
    circuit.add((gates.CZ(q, q+1) for q in range(1, nqubits-2, 2)))
    circuit.add(gates.CZ(0, nqubits-1))
circuit.add((gates.RY(q, theta=0) for q in range(nqubits)))

# Create XXZ Hamiltonian
hamiltonian = hamiltonians.XXZ(nqubits=nqubits)
# Create VQE model
vqe = models.VQE(circuit, hamiltonian)

# Optimize starting from a random guess for the variational parameters
initial_parameters = np.random.uniform(0, 2*np.pi,
                                        2*nqubits*nlayers + nqubits)
best, params, extra = vqe.minimize(initial_parameters, method='BFGS', compile=False)

For more information on the available options of the vqe.minimize call we refer to the Optimizers section of the documentation. Note that if the Stochastic Gradient Descent optimizer is used then the user has to use a backend based on tensorflow primitives and not the default custom backend, as custom operators currently do not support automatic differentiation. To switch the backend one can do qibo.set_backend("tensorflow"). Check the How to use automatic differentiation? section for more details.

When using a VQE with more than 12 qubits, it may be useful to fuse the circit implementing the ansatz using qibo.models.Circuit.fuse(). This optimizes performance by fusing the layer of one-qubit parametrized gates with the layer of two-qubit entangling gates and applying both as a single layer of general two-qubit gates (as 4x4 matrices).

circuit = models.Circuit(nqubits)
for l in range(nlayers):
    circuit.add((gates.RY(q, theta=0) for q in range(nqubits)))
    circuit.add((gates.CZ(q, q+1) for q in range(0, nqubits-1, 2)))
    circuit.add((gates.RY(q, theta=0) for q in range(nqubits)))
    circuit.add((gates.CZ(q, q+1) for q in range(1, nqubits-2, 2)))
    circuit.add(gates.CZ(0, nqubits-1))
circuit.add((gates.RY(q, theta=0) for q in range(nqubits)))
circuit = circuit.fuse()

How to write a custom variational circuit optimization?#

Similarly to the VQE, a custom implementation of a Variational Quantum Circuit (VQC) model can be achieved by defining a custom loss function and calling the qibo.optimizers.optimize() method.

Here is a simple example using a custom loss function:

import numpy as np
from qibo import models, gates
from qibo.optimizers import optimize

# custom loss function, computes fidelity
def myloss(parameters, circuit, target):
    circuit.set_parameters(parameters)
    final_state = circuit().state()
    return 1 - np.abs(np.conj(target).dot(final_state))

nqubits = 6
nlayers  = 2

# Create variational circuit
c = models.Circuit(nqubits)
for l in range(nlayers):
    c.add((gates.RY(q, theta=0) for q in range(nqubits)))
    c.add((gates.CZ(q, q+1) for q in range(0, nqubits-1, 2)))
    c.add((gates.RY(q, theta=0) for q in range(nqubits)))
    c.add((gates.CZ(q, q+1) for q in range(1, nqubits-2, 2)))
    c.add(gates.CZ(0, nqubits-1))
c.add((gates.RY(q, theta=0) for q in range(nqubits)))

# Optimize starting from a random guess for the variational parameters
x0 = np.random.uniform(0, 2*np.pi, 2*nqubits*nlayers + nqubits)
data = np.random.normal(0, 1, size=2**nqubits)

# perform optimization
best, params, extra = optimize(myloss, x0, args=(c, data), method='BFGS')

# set final solution to circuit instance
c.set_parameters(params)

How to use the QAOA?#

The quantum approximate optimization algorithm (QAOA) was introduced in arXiv:1411.4028 and is a prominent algorithm for solving hard optimization problems using the circuit-based model of quantum computation. Qibo provides an implementation of the QAOA as a model that can be defined using a qibo.hamiltonians.Hamiltonian. When properly optimized, the QAOA ansatz will approximate the ground state of this Hamiltonian. Here is a simple example using the Heisenberg XXZ Hamiltonian:

import numpy as np
from qibo import models, hamiltonians

# Create XXZ Hamiltonian for six qubits
hamiltonian = hamiltonians.XXZ(6)
# Create QAOA model
qaoa = models.QAOA(hamiltonian)

# Optimize starting from a random guess for the variational parameters
initial_parameters = 0.01 * np.random.uniform(0,1,4)
best_energy, final_parameters, extra = qaoa.minimize(initial_parameters, method="BFGS")

In the above example the initial guess for parameters has length four and therefore the QAOA ansatz consists of four operators, two using the hamiltonian and two using the mixer Hamiltonian. The user may specify the mixer Hamiltonian when defining the QAOA model, otherwise qibo.hamiltonians.X will be used by default. Note that the user may set the values of the variational parameters explicitly using qibo.models.QAOA.set_parameters(). Similarly to the VQE, we refer to Optimizers for more information on the available options of the qaoa.minimize.

QAOA uses the |++...+> as the default initial state on which the variational operators are applied. The user may specify a different initial state when executing or optimizing by passing the initial_state argument.

The QAOA model uses Solvers to apply the exponential operators to the state vector. For more information on how solvers work we refer to the How to simulate time evolution? section. When a qibo.hamiltonians.Hamiltonian is used then solvers will exponentiate it using its full matrix. Alternatively, if a qibo.hamiltonians.SymbolicHamiltonian is used then solvers will fall back to traditional Qibo circuits that perform Trotter steps. For more information on how the Trotter decomposition is implemented in Qibo we refer to the Using Trotter decomposition example.

When Trotter decomposition is used, it is possible to execute the QAOA circuit on multiple devices, by passing an accelerators dictionary when defining the model. For example the previous example would have to be modified as:

from qibo import models, hamiltonians

hamiltonian = hamiltonians.XXZ(6, dense=False)
qaoa = models.QAOA(hamiltonian, accelerators={"/GPU:0": 1, "/GPU:1": 1})

How to use automatic differentiation?#

As a deep learning framework, Tensorflow supports automatic differentiation. This can be used to optimize the parameters of variational circuits. For example the following script optimizes the parameters of two rotations so that the circuit output matches a target state using the fidelity as the corresponding loss function.

import qibo
qibo.set_backend("tensorflow")
import tensorflow as tf
from qibo import gates, models

# Optimization parameters
nepochs = 1000
optimizer = tf.keras.optimizers.Adam()
target_state = tf.ones(4, dtype=tf.complex128) / 2.0

# Define circuit ansatz
params = tf.Variable(
    tf.random.uniform((2,), dtype=tf.float64).astype(tf.complex128)
)
c = models.Circuit(2)
c.add(gates.RX(0, params[0]))
c.add(gates.RY(1, params[1]))

for _ in range(nepochs):
    with tf.GradientTape() as tape:
        c.set_parameters(params)
        final_state = c().state()
        fidelity = tf.math.abs(tf.reduce_sum(tf.math.conj(target_state) * final_state))
        loss = 1 - fidelity
    grads = tape.gradient(loss, params)
    grads = tf.math.real(grads)
    optimizer.apply_gradients(zip([grads], [params]))

Note that the "tensorflow" backend has to be used here because other custom backends do not support automatic differentiation.

The optimization procedure may also be compiled, however in this case it is not possible to use qibo.circuit.Circuit.set_parameters() as the circuit needs to be defined inside the compiled tf.GradientTape(). For example:

import qibo
qibo.set_backend("tensorflow")
import tensorflow as tf
from qibo import gates, models

nepochs = 1000
optimizer = tf.keras.optimizers.Adam()
target_state = tf.ones(4, dtype=tf.complex128) / 2.0
params = tf.Variable(tf.random.uniform((2,), dtype=tf.float64))

@tf.function
def optimize(params):
    with tf.GradientTape() as tape:
        c = models.Circuit(2)
        c.add(gates.RX(0, theta=params[0]))
        c.add(gates.RY(1, theta=params[1]))
        final_state = c().state()
        fidelity = tf.math.abs(tf.reduce_sum(tf.math.conj(target_state) * final_state))
        loss = 1 - fidelity
    grads = tape.gradient(loss, params)
    grads = tf.math.real(grads)
    optimizer.apply_gradients(zip([grads], [params]))

for _ in range(nepochs):
    optimize(params)

The user may also use tf.Variable and parametrized gates in any other way that is supported by Tensorflow, such as defining custom Keras layers and using the Sequential model API to train them.

How to perform noisy simulation?#

Qibo can perform noisy simulation with two different methods: by repeating the circuit execution multiple times and applying noise gates probabilistically or by using density matrices and applying noise channels. The two methods are analyzed in the following sections.

Moreover, Qibo provides functionality to add bit-flip errors to measurements after the simulation is completed. This is analyzed in Measurement errors.

Using density matrices#

Qibo circuits can evolve density matrices if they are initialized using the density_matrix=True flag, for example:

import qibo
qibo.set_backend("qibojit")

from qibo import models, gates

# Define circuit
c = models.Circuit(2, density_matrix=True)
c.add(gates.H(0))
c.add(gates.H(1))
# execute using the default initial state |00><00|
result = c() # will be |++><++|

will perform the transformation

\[|00 \rangle \langle 00| \rightarrow (H_1 \otimes H_2)|00 \rangle \langle 00|(H_1 \otimes H_2)^\dagger = |++ \rangle \langle ++|\]

Similarly to state vector circuit simulation, the user may specify a custom initial density matrix by passing the corresponding array when executing the circuit. If a state vector is passed as an initial state in a density matrix circuit, it will be transformed automatically to the equivalent density matrix.

Additionally, Qibo provides several gates that represent channels which can be used during a density matrix simulation. We refer to the Channels section of the documentation for a complete list of the available channels. Noise can be simulated using these channels, for example:

from qibo import models, gates

c = models.Circuit(2, density_matrix=True) # starts with state |00><00|
c.add(gates.X(1))
# transforms |00><00| -> |01><01|
c.add(gates.PauliNoiseChannel(0, [("X", 0.3)]))
# transforms |01><01| -> (1 - px)|01><01| + px |11><11|
result = c()
# result.state() will be tf.Tensor(diag([0, 0.7, 0, 0.3]))

will perform the transformation

\[\begin{split}|00\rangle \langle 00|& \rightarrow (I \otimes X)|00\rangle \langle 00|(I \otimes X) = |01\rangle \langle 01| \\& \rightarrow 0.7|01\rangle \langle 01| + 0.3(X\otimes I)|01\rangle \langle 01|(X\otimes I)^\dagger \\& = 0.7|01\rangle \langle 01| + 0.3|11\rangle \langle 11|\end{split}\]

Measurements and callbacks can be used with density matrices exactly as in the case of state vector simulation.

Using repeated execution#

Simulating noise with density matrices is memory intensive as it effectively doubles the number of qubits. Qibo provides an alternative way of simulating the effect of channels without using density matrices, which relies on state vectors and repeated circuit execution with sampling. Noise can be simulated by creating a normal (non-density matrix) circuit and repeating its execution as follows:

import numpy as np
from qibo import models, gates

# Define circuit
c = models.Circuit(5)
thetas = np.random.random(5)
c.add((gates.RX(i, theta=t) for i, t in enumerate(thetas)))
# Add noise channels to all qubits
c.add((gates.PauliNoiseChannel(i, [("X", 0.2), ("Y", 0.0), ("Z", 0.3)])
       for i in range(5)))
# Add measurement of all qubits
c.add(gates.M(*range(5)))

# Repeat execution 1000 times
result = c(nshots=1000)

In this example the simulation is repeated 1000 times and the action of the qibo.gates.PauliNoiseChannel gate differs each time, because the error X, Y and Z gates are sampled according to the given probabilities. Note that when a channel is used, the command c(nshots=1000) has a different behavior than what is described in How to perform measurements?. Normally c(nshots=1000) would execute the circuit once and would then sample 1000 bit-strings from the final state. When channels are used, the full is executed 1000 times because the behavior of channels is probabilistic and different in each execution. Note that now the simulation time required will increase linearly with the number of repetitions (nshots).

Note that executing a circuit with channels only once is possible, however, since the channel acts probabilistically, the results of a single execution are random and usually not useful on their own. It is possible also to use repeated execution with noise channels even without the presence of measurements. If c(nshots=1000) is called for a circuit that contains channels but no measurements measurements then the circuit will be executed 1000 times and the final 1000 state vectors will be returned as a tensor of shape (nshots, 2 ^ nqubits). Note that this tensor is usually large and may lead to memory errors, therefore this usage is not advised.

Unlike the density matrix approach, it is not possible to use every channel with sampling and repeated execution. Specifically, qibo.gates.UnitaryChannel and qibo.gates.PauliNoiseChannel can be used with sampling, while qibo.gates.KrausChannel requires density matrices.

Adding noise after every gate#

In practical applications noise typically occurs after every gate. Qibo provides the qibo.models.circuit.Circuit.with_pauli_noise() method which automatically creates a new circuit that contains a qibo.gates.PauliNoiseChannel after every gate. The user can control the probabilities of the noise channel using a noise map, which is a dictionary that maps qubits to the corresponding probability triplets. For example, the following script

from qibo import models, gates

c = models.Circuit(2)
c.add([gates.H(0), gates.H(1), gates.CNOT(0, 1)])

# Define a noise map that maps qubit IDs to noise probabilities
noise_map = {0: list(zip(["X", "Z"], [0.1, 0.2])), 1: list(zip(["Y", "Z"], [0.2, 0.1]))}
noisy_c = c.with_pauli_noise(noise_map)

will create a new circuit noisy_c that is equivalent to:

noisy_c2 = models.Circuit(2)
noisy_c2.add(gates.H(0))
noisy_c2.add(gates.PauliNoiseChannel(0, [("X", 0.1), ("Y", 0.0), ("Z", 0.2)]))
noisy_c2.add(gates.H(1))
noisy_c2.add(gates.PauliNoiseChannel(1, [("X", 0.0), ("Y", 0.2), ("Z", 0.1)]))
noisy_c2.add(gates.CNOT(0, 1))
noisy_c2.add(gates.PauliNoiseChannel(0, [("X", 0.1), ("Y", 0.0), ("Z", 0.2)]))
noisy_c2.add(gates.PauliNoiseChannel(1, [("X", 0.0), ("Y", 0.2), ("Z", 0.1)]))

Note that noisy_c uses the gate objects of the original circuit c (it is not a deep copy), while in noisy_c2 each gate was created as a new object.

The user may use a single tuple instead of a dictionary as the noise map In this case the same probabilities will be applied to all qubits. That is noise_map = list(zip(["X", "Z"], [0.1, 0.1])) is equivalent to noise_map = {0: list(zip(["X", "Z"], [0.1, 0.1])), 1: list(zip(["X", "Z"], [0.1, 0.1])), ...}.

As described in the previous sections, if qibo.models.circuit.Circuit.with_pauli_noise() is used in a circuit that uses state vectors then noise will be simulated with repeated execution. If the user wishes to use density matrices instead, this is possible by passing the density_matrix=True flag during the circuit initialization and call .with_pauli_noise on the new circuit.

Using a noise model#

In a real quantum circuit some gates can be highly faulty and introduce errors. In order to simulate this behavior Qibo provides the qibo.noise.NoiseModel class which can store errors that are gate-dependent using the qibo.noise.NoiseModel.add() method and generate the corresponding noisy circuit with qibo.noise.NoiseModel.apply(). The corresponding noise is applied after every instance of the gate in the circuit. It is also possible to specify on which qubits the noise will be added.

The current quantum errors available to build a custom noise model are: qibo.noise.PauliError, qibo.noise.ThermalRelaxationError and qibo.noise.ResetError.

Here is an example on how to use a noise model:

import numpy as np
from qibo import models, gates
from qibo.noise import NoiseModel, PauliError

# Build specific noise model with 3 quantum errors:
# - Pauli error on H only for qubit 1.
# - Pauli error on CNOT for all the qubits.
# - Pauli error on RX(pi/2) for qubit 0.
noise = NoiseModel()
noise.add(PauliError([("X", 0.5)]), gates.H, 1)
noise.add(PauliError([("Y", 0.5)]), gates.CNOT)
is_sqrt_x = (lambda g: np.pi/2 in g.parameters)
noise.add(PauliError([("X", 0.5)]), gates.RX, qubits=0, conditions=is_sqrt_x)

# Generate noiseless circuit.
c = models.Circuit(2)
c.add([gates.H(0), gates.H(1), gates.CNOT(0, 1), gates.RX(0, np.pi/2),  gates.RX(0, 3*np.pi/2), gates.RX(1, np.pi/2)])

# Apply noise to the circuit according to the noise model.
noisy_c = noise.apply(c)

The noisy circuit defined above will be equivalent to the following circuit:

noisy_c2 = models.Circuit(2)
noisy_c2.add(gates.H(0))
noisy_c2.add(gates.H(1))
noisy_c2.add(gates.PauliNoiseChannel(1, [("X", 0.5)]))
noisy_c2.add(gates.CNOT(0, 1))
noisy_c2.add(gates.PauliNoiseChannel(0, [("Y", 0.5)]))
noisy_c2.add(gates.PauliNoiseChannel(1, [("Y", 0.5)]))
noisy_c2.add(gates.RX(0, np.pi/2))
noisy_c2.add(gates.PauliNoiseChannel(0, [("X", 0.5)]))
noisy_c2.add(gates.RX(0, 3*np.pi/2))
noisy_c2.add(gates.RX(1, np.pi/2))

The qibo.noise.NoiseModel class supports also density matrices, it is sufficient to pass a circuit which was initialized with density_matrix=True.

Measurement errors#

qibo.measurements.CircuitResult provides qibo.measurements.CircuitResult.apply_bitflips() which allows adding bit-flip errors to the sampled bit-strings without having to re-execute the simulation. For example:

import numpy as np
from qibo import models, gates

thetas = np.random.random(4)
c = models.Circuit(4)
c.add((gates.RX(i, theta=t) for i, t in enumerate(thetas)))
c.add([gates.M(0, 1), gates.M(2, 3)])
result = c(nshots=100)
# add bit-flip errors with probability 0.2 for all qubits
result.apply_bitflips(0.2)
# add bit-flip errors with different probabilities for each qubit
error_map = {0: 0.2, 1: 0.1, 2: 0.3, 3: 0.1}
result.apply_bitflips(error_map)

The corresponding noisy samples and frequencies can then be obtained as described in the How to perform measurements? example.

Note that qibo.measurements.CircuitResult.apply_bitflips() modifies the measurement samples contained in the corresponding state and therefore the original noiseless measurement samples are no longer accessible. It is possible to keep the original samples by creating a copy of the states before applying the bitflips:

import numpy as np
from qibo import models, gates

thetas = np.random.random(4)
c = models.Circuit(4)
c.add((gates.RX(i, theta=t) for i, t in enumerate(thetas)))
c.add([gates.M(0, 1), gates.M(2, 3)])
result = c(nshots=100)
# add bit-flip errors with probability 0.2 for all qubits
result.apply_bitflips(0.2)
# add bit-flip errors with different probabilities for each qubit
error_map = {0: 0.2, 1: 0.1, 2: 0.3, 3: 0.1}
result.apply_bitflips(error_map)

Alternatively, the user may specify a bit-flip error map when defining measurement gates:

import numpy as np
from qibo import models, gates

thetas = np.random.random(6)
c = models.Circuit(6)
c.add((gates.RX(i, theta=t) for i, t in enumerate(thetas)))
c.add(gates.M(0, 1, p0=0.2))
c.add(gates.M(2, 3, p0={2: 0.1, 3: 0.0}))
c.add(gates.M(4, 5, p0=[0.4, 0.3]))
result = c(nshots=100)

In this case result will contain noisy samples according to the given bit-flip probabilities. The probabilities can be given as a dictionary (must contain all measured qubits as keys), a list (must have the sample as the measured qubits) or a single float number (to be used on all measured qubits). Note that, unlike the previous code example, when bit-flip errors are incorporated as part of measurement gates it is not possible to access the noiseless samples.

Moreover, it is possible to simulate asymmetric bit-flips using the p1 argument as result.apply_bitflips(p0=0.2, p1=0.1). In this case a probability of 0.2 will be used for 0->1 errors but 0.1 for 1->0 errors. Similarly to p0, p1 can be a single float number or a dictionary and can be used both in qibo.measurements.CircuitResult.apply_bitflips() and the measurement gate. If p1 is not specified the value of p0 will be used for both errors.

Simulating IBMQ’s quantum hardware#

Qibo can perform a simulation of a real quantum computer using the qibo.noise.IBMQNoiseModel class. It is possible by passing the circuit instance that we want to simulate and the noise channels’ parameters as a dictionary. In this model, the user must set the relaxation times t1 and t2 for each qubit, an approximated gate times, and depolarizing errors for each one-qubit (depolarizing_one_qubit) and two-qubit (depolarizing_two_qubit) gates. Additionally, one can also pass single-qubit readout error probabilities (readout_one_qubit).

from qibo import Circuit, gates
from qibo.noise import IBMQNoiseModel

nqubits = 2
circuit = Circuit(2, density_matrix=True)
circuit.add(
    [
        gates.H(0),
        gates.X(1),
        gates.Z(0),
        gates.X(0),
        gates.CNOT(0,1),
        gates.CNOT(1, 0),
        gates.X(1),
        gates.Z(1),
        gates.M(0),
        gates.M(1),
    ]
)

print("raw circuit:")
print(circuit.draw())

parameters = {
    "t1": {"0": 250*1e-06, "1": 240*1e-06},
    "t2": {"0": 150*1e-06, "1": 160*1e-06},
    "gate_times" : (200*1e-9, 400*1e-9),
    "excited_population" : 0,
    "depolarizing_one_qubit" : 4.000e-4,
    "depolarizing_two_qubit": 1.500e-4,
    "readout_one_qubit" : {"0": (0.022, 0.034), "1": (0.015, 0.041)},
    }

noise_model = IBMQNoiseModel()
noise_model.from_dict(parameters)
noisy_circuit = noise_model.apply(circuit)

print("noisy circuit:")
print(noisy_circuit.draw())

noisy_circuit is the new circuit containing the error gate channels.

How to perform error mitigation?#

Noise and errors in circuits are one of the biggest obstacles to face in quantum computing. Say that you have a circuit \(C\) and you want to measure an observable \(A\) at the end of it, in general you are going to obtain an expected value \(\langle A \rangle_{noisy}\) that can lie quiet far from the true one \(\langle A \rangle_{exact}\). In Qibo, different methods are implemented for mitigating errors in circuits and obtaining a better estimate of the noise-free expected value \(\langle A \rangle_{exact}\).

Let’s see how to use them. For starters, let’s define a dummy circuit with some RZ, RX and CNOT gates:

import numpy as np

from qibo import Circuit, gates

# Define the circuit
nqubits = 3
hz = 0.5
hx = 0.5
dt = 0.25
circ = Circuit(nqubits, density_matrix=True)
circ.add(gates.RZ(q, theta=-2 * hz * dt - np.pi / 2) for q in range(nqubits))
circ.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits))
circ.add(gates.RZ(q, theta=-2 * hx * dt + np.pi) for q in range(nqubits))
circ.add(gates.RX(q, theta=np.pi / 2) for q in range(nqubits))
circ.add(gates.RZ(q, theta=-np.pi / 2) for q in range(nqubits))
circ.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2))
circ.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(0, nqubits - 1, 2))
circ.add(gates.CNOT(q, q + 1) for q in range(0, nqubits - 1, 2))
circ.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
circ.add(gates.RZ(q + 1, theta=-2 * dt) for q in range(1, nqubits, 2))
circ.add(gates.CNOT(q, q + 1) for q in range(1, nqubits, 2))
# Include the measurements
circ.add(gates.M(*range(nqubits)))

# visualize the circuit
print(circ.draw())

#  q0: ─RZ─RX─RZ─RX─RZ─o────o────────M─
#  q1: ─RZ─RX─RZ─RX─RZ─X─RZ─X─o────o─M─
#  q2: ─RZ─RX─RZ─RX─RZ────────X─RZ─X─M─

remember to initialize the circuit with density_matrix=True and to include the measuerement gates at the end for expectation value calculation.

As observable we can simply take \(Z_0 Z_1 Z_2\) :

from qibo.symbols import Z
from qibo.hamiltonians import SymbolicHamiltonian
from qibo.backends import GlobalBackend

backend = GlobalBackend()

# Define the observable
obs = np.prod([Z(i) for i in range(nqubits)])
obs = SymbolicHamiltonian(obs, backend=backend)

We can obtain the exact expected value by running the circuit on any simulation backend. To mimic the execution on the real quantum hardware, instead, we can use a noise model:

# Noise-free expected value
exact = obs.expectation(backend.execute_circuit(circ).state())
print(exact)
# 0.9096065335014379

from qibo.noise import DepolarizingError, ReadoutError, NoiseModel
from qibo.quantum_info import random_stochastic_matrix

# Define the noise model
noise =  NoiseModel()
# depolarizing error after each CNOT
noise.add(DepolarizingError(0.1), gates.CNOT)
# readout error
# randomly initialize the bitflip probabilities
prob = random_stochastic_matrix(
    2**nqubits, diagonally_dominant=True, seed=2, backend=backend
)
noise.add(ReadoutError(probabilities=prob), gate=gates.M)
# Noisy expected value without mitigation
noisy = obs.expectation(backend.execute_circuit(noise.apply(circ)).state())
print(noisy)
# 0.5647937721701448

Note that when running on the quantum hardware, you won’t need to use a noise model anymore, you will just have to change the backend to the appropriate one.

Now let’s check that error mitigation produces better estimates of the exact expected value.

Readout Mitigation#

Firstly, let’s try to mitigate the readout errors. To do this, we can either compute the response matrix and use it modify the final state after the circuit execution:

from qibo.models.error_mitigation import get_expectation_val_with_readout_mitigation, get_response_matrix

nshots = 10000
# compute the response matrix
response_matrix = get_response_matrix(
    nqubits, backend=backend, noise_model=noise, nshots=nshots
)
# define mitigation options
readout = {"response_matrix": response_matrix}
# mitigate the readout errors
mit_val = get_expectation_val_with_readout_mitigation(circ, obs, noise, readout=readout)
print(mit_val)
# 0.5945794816381054

Or use the randomized readout mitigation:

from qibo.models.error_mitigation import apply_randomized_readout_mitigation

# define mitigation options
readout = {"ncircuits": 10}
# mitigate the readout errors
mit_val = get_expectation_val_with_readout_mitigation(circ, obs, noise, readout=readout)
print(mit_val)
# 0.5860884499785314

Alright, the expected value is improving, but we are still far from the ideal one. Readout mitigation alone is not enough, let’s try to use some more advanced methods to get rid of the depolarizing error we introduced in the CNOT gates.

Zero Noise Extrapolation (ZNE)#

To run ZNE, we just need to define the noise levels to use. Each level corresponds to the number of CNOT or RX pairs (depending on the value of insertion_gate) inserted in the circuit in correspondence to the original ones. Since we decided to simulate noisy CNOTs:

Level 1
q0: ─X─  -->  q0: ─X───X──X─
q1: ─o─  -->  q1: ─o───o──o─

Level 2
q0: ─X─  -->  q0: ─X───X──X───X──X─
q1: ─o─  -->  q1: ─o───o──o───o──o─

.
.
.

For example if we use the five levels [0,1,2,3,4] :

from qibo.models.error_mitigation import ZNE

# Mitigated expected value
estimate = ZNE(
    circuit=circ,
    observable=obs,
    noise_levels=np.arange(5),
    noise_model=noise,
    nshots=10000,
    insertion_gate='CNOT',
    backend=backend,
)
print(estimate)
# 0.8332843749999996

we get an expected value closer to the exact one. We can further improve by using ZNE combined with the readout mitigation:

# we can either use
# the response matrix computed earlier
readout = {'response_matrix': response_matrix}
# or the randomized readout
readout = {'ncircuits': 10}

# Mitigated expected value
estimate = ZNE(
    circuit=circ,
    observable=obs,
    backend=backend,
    noise_levels=np.arange(5),
    noise_model=noise,
    nshots=10000,
    insertion_gate='CNOT',
    readout=readout,
)
print(estimate)
# 0.8979124892467807

Clifford Data Regression (CDR)#

For CDR instead, you don’t need to define anything additional. However, keep in mind that the input circuit is expected to be decomposed in the set of primitive gates \(RX(\frac{\pi}{2}), CNOT, X\) and \(RZ(\theta)\).

from qibo.models.error_mitigation import CDR

# Mitigated expected value
estimate = CDR(
    circuit=circ,
    observable=obs,
    n_training_samples=10,
    backend=backend,
    noise_model=noise,
    nshots=10000,
    readout=readout,
)
print(estimate)
# 0.8983676333969615

Again, the mitigated expected value improves over the noisy one and is also slightly better compared to ZNE.

Variable Noise CDR (vnCDR)#

Being a combination of ZNE and CDR, vnCDR requires you to define the noise levels as done in ZNE, and the same caveat about the input circuit for CDR is valid here as well.

from qibo.models.error_mitigation import vnCDR

# Mitigated expected value
estimate = vnCDR(
    circuit=circ,
    observable=obs,
    n_training_samples=10,
    backend=backend,
    noise_levels=np.arange(3),
    noise_model=noise,
    nshots=10000,
    insertion_gate='CNOT',
    readout=readout,
)
print(estimate)
# 0.8998376314644383

The result is similar to the one obtained by CDR. Usually, one would expect slightly better results for vnCDR, however, this can substantially vary depending on the circuit and the observable considered and, therefore, it is hard to tell a priori.

Importance Clifford Sampling (ICS)#

The use of iCS is straightforward, analogous to CDR and vnCDR.

from qibo.models.error_mitigation import ICS

# Mitigated expected value
estimate = ICS(
    circuit=circ,
    observable=obs,
    n_training_samples=10,
    backend=backend,
    noise_model=noise,
    nshots=10000,
    readout=readout,
)
print(estimate)
# 0.9183495097398502

Again, the mitigated expected value improves over the noisy one and is also slightly better compared to ZNE. This was just a basic example usage of the three methods, for all the details about them you should check the API-reference page Error Mitigation.

How to simulate time evolution?#

Simulating the unitary time evolution of quantum states is useful in many physics applications including the simulation of adiabatic quantum computation. Qibo provides the qibo.models.StateEvolution model that simulates unitary evolution using the full state vector. For example:

import numpy as np
from qibo import hamiltonians, models

# Define evolution model under the non-interacting sum(Z) Hamiltonian
# with time step dt=1e-1
nqubits = 4
evolve = models.StateEvolution(hamiltonians.Z(nqubits), dt=1e-1)
# Define initial state as |++++>
initial_state = np.ones(2 ** nqubits) / np.sqrt(2 ** nqubits)
# Get the final state after time t=2
final_state = evolve(final_time=2, initial_state=initial_state)

When studying dynamics people are usually interested not only in the final state vector but also in observing how physical quantities change during the time evolution. This is possible using callbacks. For example, in the above case we can track how <X> changes as follows:

import numpy as np
from qibo import hamiltonians, models, callbacks

nqubits = 4
# Define a callback that calculates the energy (expectation value) of the X Hamiltonian
observable = callbacks.Energy(hamiltonians.X(nqubits))
# Create evolution object using the above callback and a time step of dt=1e-3
evolve = models.StateEvolution(hamiltonians.Z(nqubits), dt=1e-3,
                               callbacks=[observable])
# Evolve for total time t=1
initial_state = np.ones(2 ** nqubits) / np.sqrt(2 ** nqubits)
final_state = evolve(final_time=1, initial_state=initial_state)

print(observable[:])
# will print an array of shape ``(1001,)`` that holds <X>(t) values

Note that the time step dt=1e-3 defines how often we calculate <X> during the evolution.

In the above cases the exact time evolution operator (exponential of the Hamiltonian) was used to evolve the state vector. Because the evolution Hamiltonian is time-independent, the matrix exponentiation happens only once. It is possible to simulate time-dependent Hamiltonians by passing a function of time instead of a qibo.hamiltonians.Hamiltonian in the qibo.models.StateEvolution model. For example:

import numpy as np
from qibo import hamiltonians, models

# Defina a time dependent Hamiltonian
nqubits = 4
ham = lambda t: np.cos(t) * hamiltonians.Z(nqubits)
# and pass it to the evolution model
evolve = models.StateEvolution(ham, dt=1e-3)
initial_state = np.ones(2 ** nqubits) / np.sqrt(2 ** nqubits)
final_state = evolve(final_time=1, initial_state=initial_state)

The above script will still use the exact time evolution operator with the exponentiation repeated for each time step. The integration method can be changed using the solver argument when executing. The solvers that are currently implemented are the default exponential solver ("exp") and two Runge-Kutta solvers: fourth-order ("rk4") and fifth-order ("rk45"). For more information we refer to the Solvers section.

Using Trotter decomposition#

Trotter decomposition provides a way to represent the unitary evolution of quantum states as a sequence of local unitaries. This allows to represent the physical process of time evolution as a quantum circuit. Qibo provides functionality to perform this transformation automatically, if the underlying Hamiltonian object is defined as a sum of commuting parts that consist of terms that can be exponentiated efficiently. Such Hamiltonian can be implemented in Qibo using qibo.hamiltonians.SymbolicHamiltonian. The implementation of Trotter decomposition is based on Sec. 4.1 of arXiv:1901.05824. Below is an example of how to use this object in practice:

from qibo import hamiltonians

# Define TFIM model as a non-dense ``SymbolicHamiltonian``
ham = hamiltonians.TFIM(nqubits=5, dense=False)
# This object can be used to create the circuit that
# implements a single Trotter time step ``dt``
circuit = ham.circuit(dt=1e-2)

This is a standard qibo.core.circuit.Circuit that contains qibo.gates.Unitary gates corresponding to the exponentials of the Trotter decomposition and can be executed on any state.

Note that in the transverse field Ising model (TFIM) that was used in this example is among the pre-coded Hamiltonians in Qibo and could be created as a qibo.hamiltonians.SymbolicHamiltonian simply using the dense=False flag. For more information on the difference between dense and non-dense Hamiltonians we refer to the Hamiltonians section. Note that only non-dense Hamiltonians created using dense=False or through the qibo.hamiltonians.SymbolicHamiltonian object can be used for evolution using Trotter decomposition. If a dense Hamiltonian is used then evolution will be done by exponentiating the full Hamiltonian matrix.

Defining custom Hamiltonians from terms can be more complicated, however Qibo simplifies this process by providing the option to define Hamiltonians symbolically through the use of sympy. For more information on this we refer to the How to define custom Hamiltonians using symbols? example.

A qibo.hamiltonians.SymbolicHamiltonian can also be used to simulate time evolution. This can be done by passing the Hamiltonian to a qibo.models.StateEvolution model and using the exponential solver. For example:

import numpy as np
from qibo import models, hamiltonians

nqubits = 5
# Create a critical TFIM Hamiltonian as ``SymbolicHamiltonian``
ham = hamiltonians.TFIM(nqubits=nqubits, h=1.0, dense=False)
# Define the |+++++> initial state
initial_state = np.ones(2 ** nqubits) / np.sqrt(2 ** nqubits)
# Define the evolution model
evolve = models.StateEvolution(ham, dt=1e-3)
# Evolve for total time T=1
final_state = evolve(final_time=1, initial_state=initial_state)

This script creates the Trotter circuit for dt=1e-3 and applies it repeatedly to the given initial state T / dt = 1000 times to obtain the final state of the evolution.

Since Trotter evolution is based on Qibo circuits, it also supports distributed execution on multiple devices (GPUs). This can be enabled by passing an accelerators dictionary when defining the qibo.models.StateEvolution model. We refer to the How to select hardware devices? example for more details on how the accelerators dictionary can be used.

How to simulate adiabatic time evolution?#

Qibo provides the qibo.models.AdiabaticEvolution model to simulate adiabatic time evolution. This is a special case of the qibo.models.StateEvolution model analyzed in the previous example where the evolution Hamiltonian is interpolated between an initial “easy” Hamiltonian and a “hard” Hamiltonian that usually solves an optimization problem. Here is an example of adiabatic evolution simulation:

import numpy as np
from qibo import hamiltonians, models

nqubits = 4
T = 1 # total evolution time
# Define the easy and hard Hamiltonians
h0 = hamiltonians.X(nqubits)
h1 = hamiltonians.TFIM(nqubits, h=0)
# Define the interpolation scheduling
s = lambda t: t
# Define evolution model
evolve = models.AdiabaticEvolution(h0, h1, s, dt=1e-2)
# Get the final state of the evolution
final_state = evolve(final_time=T)

According to the adiabatic theorem, for proper scheduling and total evolution time the final_state should approximate the ground state of the “hard” Hamiltonian.

If the initial state is not specified, the ground state of the easy Hamiltonian will be used, which is common for adiabatic evolution. A distributed execution can be used by passing an accelerators dictionary during the initialization of the AdiabaticEvolution model. In this case the default initial state is |++...+> (full superposition in the computational basis).

Callbacks may also be used as in the previous example. An additional callback (qibo.callbacks.Gap) is available for calculating the energies and the gap of the adiabatic evolution Hamiltonian. Its usage is similar to other callbacks:

import numpy as np
from qibo import hamiltonians, models, callbacks

nqubits = 4
h0 = hamiltonians.X(nqubits)
h1 = hamiltonians.TFIM(nqubits, h=0)

ground = callbacks.Gap(mode=0)
# define a callback for calculating the gap
gap = callbacks.Gap()
# define and execute the ``AdiabaticEvolution`` model
evolution = models.AdiabaticEvolution(h0, h1, lambda t: t, dt=1e-1,
                                      callbacks=[gap, ground])

final_state = evolution(final_time=1.0)
# print the values of the gap at each evolution time step
print(gap[:])

The scheduling function s should be a callable that accepts one (s(t)) or two (s(t, p)) arguments. The first argument accepts values in [0, 1] and corresponds to the ratio t / final_time during evolution. The second optional argument is a vector of free parameters that can be optimized. The function should, by definition, satisfy the properties s(0, p) = 0 and s(1, p) = 1 for any p, otherwise errors will be raised.

All state evolution functionality described in the previous example can also be used for simulating adiabatic evolution. The solver can be specified during the initialization of the qibo.models.AdiabaticEvolution model and a Trotter decomposition may be used with the exponential solver. The Trotter decomposition will be used automatically if h0 and h1 are defined using as qibo.hamiltonians.SymbolicHamiltonian objects. For pre-coded Hamiltonians this can be done simply as:

from qibo import hamiltonians, models

nqubits = 4
# Define ``SymolicHamiltonian``s
h0 = hamiltonians.X(nqubits, dense=False)
h1 = hamiltonians.TFIM(nqubits, h=0, dense=False)
# Perform adiabatic evolution using the Trotter decomposition
evolution = models.AdiabaticEvolution(h0, h1, lambda t: t, dt=1e-1)
final_state = evolution(final_time=1.0)

When Trotter evolution is used, it is also possible to execute on multiple devices by passing an accelerators dictionary in the creation of the qibo.models.AdiabaticEvolution model.

Note that h0 and h1 should have the same type, either both qibo.hamiltonians.Hamiltonian or both qibo.hamiltonians.SymbolicHamiltonian.

Optimizing the scheduling function#

The free parameters p of the scheduling function can be optimized using the qibo.models.AdiabaticEvolution.minimize() method. The parameters are optimized so that the final state of the adiabatic evolution approximates the ground state of the “hard” Hamiltonian. Optimization is similar to what is described in the How to write a VQE? example and can be done as follows:

import numpy as np
from qibo import hamiltonians, models

# Define Hamiltonians
h0 = hamiltonians.X(3)
h1 = hamiltonians.TFIM(3)
# Define scheduling function with a free variational parameter ``p``
sp = lambda t, p: (1 - p) * np.sqrt(t) + p * t
# Define an evolution model with dt=1e-2
evolution = models.AdiabaticEvolution(h0, h1, sp, dt=1e-2)
# Find the optimal value for ``p`` starting from ``p = 0.5`` and ``T=1``.
initial_guess = [0.5, 1]
# best, params, extra = evolution.minimize(initial_guess, method="BFGS", options={'disp': True})
print(best) # prints the best energy <H1> found from the final state
print(params) # prints the optimal values for the parameters.

Note that the minimize method optimizes both the free parameters p of the scheduling function as well as the total evolution time. The initial guess for the total evolution time is the last value of the given initial_guess array. For a list of the available optimizers we refer to Optimizers.

How to define custom Hamiltonians using symbols?#

In order to use the VQE, QAOA and time evolution models in Qibo the user has to define Hamiltonians based on qibo.hamiltonians.Hamiltonian which uses the full matrix representation of the corresponding operator or qibo.hamiltonians.SymbolicHamiltonian which uses a more efficient term representation. Qibo provides pre-coded Hamiltonians for some common models, such as the transverse field Ising model (TFIM) and the Heisenberg model (see Hamiltonians for a complete list of the pre-coded models). In order to explore other problems the user needs to define the Hamiltonian objects from scratch.

A standard way to define Hamiltonians is through their full matrix representation. For example the following code generates the TFIM Hamiltonian with periodic boundary conditions for four qubits by constructing the corresponding 16x16 matrix:

import numpy as np
from qibo import hamiltonians, matrices

# ZZ terms
matrix = np.kron(np.kron(matrices.Z, matrices.Z), np.kron(matrices.I, matrices.I))
matrix += np.kron(np.kron(matrices.I, matrices.Z), np.kron(matrices.Z, matrices.I))
matrix += np.kron(np.kron(matrices.I, matrices.I), np.kron(matrices.Z, matrices.Z))
matrix += np.kron(np.kron(matrices.Z, matrices.I), np.kron(matrices.I, matrices.Z))
# X terms
matrix += np.kron(np.kron(matrices.X, matrices.I), np.kron(matrices.I, matrices.I))
matrix += np.kron(np.kron(matrices.I, matrices.X), np.kron(matrices.I, matrices.I))
matrix += np.kron(np.kron(matrices.I, matrices.I), np.kron(matrices.X, matrices.I))
matrix += np.kron(np.kron(matrices.I, matrices.I), np.kron(matrices.I, matrices.X))
# Create Hamiltonian object
ham = hamiltonians.Hamiltonian(4, matrix)

Although it is possible to generalize the above construction to arbitrary number of qubits this procedure may be more complex for other Hamiltonians. Moreover constructing the full matrix does not scale well with increasing the number of qubits. This makes the use of qibo.hamiltonians.SymbolicHamiltonian preferrable as the qubit number increases, as this Hamiltonians is not based in the full matrix representation.

To simplify the construction of Hamiltonians, Qibo provides the qibo.hamiltonians.SymbolicHamiltonian object which allows the user to construct Hamiltonian objects by writing their symbolic form using sympy symbols. Moreover Qibo provides quantum-computation specific symbols (qibo.symbols.Symbol) such as the Pauli operators. For example, the TFIM on four qubits could be constructed as:

import numpy as np
from qibo import hamiltonians
from qibo.symbols import X, Z

# Define Hamiltonian using Qibo symbols
# ZZ terms
symbolic_ham = sum(Z(i) * Z(i + 1) for i in range(3))
# periodic boundary condition term
symbolic_ham += Z(0) * Z(3)
# X terms
symbolic_ham += sum(X(i) for i in range(4))

# Define a Hamiltonian using the above form
ham = hamiltonians.SymbolicHamiltonian(symbolic_ham)
# This Hamiltonian is memory efficient as it does not construct the full matrix

# The corresponding dense Hamiltonian which contains the full matrix can
# be constructed easily as
dense_ham = ham.dense
# and the matrix is accessed as ``dense_ham.matrix`` or ``ham.matrix``.

Defining Hamiltonians from symbols is usually a simple process as the symbolic form is very close to the form of the Hamiltonian on paper. Note that when a qibo.hamiltonians.SymbolicHamiltonian is used for time evolution, Qibo handles automatically automatically the Trotter decomposition by splitting to the appropriate terms.

Qibo symbols support an additional commutative argument which is set to False by default since quantum operators are non-commuting objects. When the user knows that the Hamiltonian consists of commuting terms only, such as products of Z operators, switching commutative to True may speed-up some symbolic calculations, such as the sympy.expand used when calculating the Trotter decomposition for the Hamiltonian. This option can be used when constructing each symbol:

from qibo import hamiltonians
from qibo.symbols import Z

form = Z(0, commutative=True) * Z(1, commutative=True) + Z(1, commutative=True) * Z(2, commutative=True)
ham = hamiltonians.SymbolicHamiltonian(form)

How to calculate expectation values using samples?#

It is possible to calculate the expectation value of a qibo.hamiltonians.Hamiltonian on a given state using the qibo.hamiltonians.Hamiltonian.expectation() method, which can be called on a state or density matrix. For example

from qibo import Circuit, gates
from qibo.hamiltonians import XXZ

circuit = Circuit(4)
circuit.add(gates.H(i) for i in range(4))
circuit.add(gates.CNOT(0, 1))
circuit.add(gates.CNOT(1, 2))
circuit.add(gates.CNOT(2, 3))

hamiltonian = XXZ(4)

result = circuit()
expectation_value = hamiltonian.expectation(result.state())

In this example, the circuit will be simulated to obtain the final state vector and the corresponding expectation value will be calculated through exact matrix multiplication with the Hamiltonian matrix. If a qibo.hamiltonians.SymbolicHamiltonian is used instead, the expectation value will be calculated as a sum of expectation values of local terms, allowing calculations of more qubits with lower memory consumption. The calculation of each local term still requires the state vector.

When executing a circuit on real hardware, usually only measurements of the state are available, not the state vector. Qibo provides qibo.hamiltonians.Hamiltonian.expectation_from_samples() to allow calculation of expectation values directly from such samples:

from qibo import Circuit, gates
from qibo.hamiltonians import Z

circuit = Circuit(4)
circuit.add(gates.H(i) for i in range(4))
circuit.add(gates.CNOT(0, 1))
circuit.add(gates.CNOT(1, 2))
circuit.add(gates.CNOT(2, 3))
circuit.add(gates.M(*range(4)))

hamiltonian = Z(4)

result = circuit(nshots=1024)
expectation_value = hamiltonian.expectation_from_samples(result.frequencies())

This example simulates the circuit similarly to the previous one but calculates the expectation value using the frequencies of shots, instead of the exact state vector. This can also be invoked directly from the result object:

expectation_value = result.expectation_from_samples(hamiltonian)

The expectation from samples currently works only for Hamiltonians that are diagonal in the computational basis.

How to modify the transpiler?#

Logical quantum circuits for quantum algorithms are hardware agnostic. Usually an all-to-all qubit connectivity is assumed while most current hardware only allows the execution of two-qubit gates on a restricted subset of qubit pairs. Moreover, quantum devices are restricted to executing a subset of gates, referred to as native. This means that, in order to execute circuits on a real quantum chip, they must be transformed into an equivalent, hardware specific, circuit. The transformation of the circuit is carried out by the transpiler through the resolution of two key steps: connectivity matching and native gates decomposition. In order to execute a gate between two qubits that are not directly connected SWAP gates are required. This procedure is called routing. As on NISQ devices two-qubit gates are a large source of noise, this procedure generates an overall noisier circuit. Therefore, the goal of an efficient routing algorithm is to minimize the number of SWAP gates introduced. An important step to ease the connectivity problem, is finding anoptimal initial mapping between logical and physical qubits. This step is called placement. The native gates decomposition in the transpiling procedure is performed by the unroller. An optimal decomposition uses the least amount of two-qubit native gates. It is also possible to reduce the number of gates of the resulting circuit by exploiting commutation relations, KAK decomposition or machine learning techniques. Qibo implements a built-in transpiler with customizable options for each step. The main algorithms that can be used at each transpiler step are reported below with a short description.

The initial placement can be found with one of the following procedures: - Trivial: logical-physical qubit mapping is an identity. - Custom: custom logical-physical qubit mapping. - Random greedy: the best mapping is found within a set of random layouts based on a greedy policy. - Subgraph isomorphism: the initial mapping is the one that guarantees the execution of most gates at the beginning of the circuit without introducing any SWAP. - Reverse traversal: this technique uses one or more reverse routing passes to find an optimal mapping by starting from a trivial layout.

The routing problem can be solved with the following algorithms: - Shortest paths: when unconnected logical qubits have to interact, they are moved on the chip on the shortest path connecting them. When multiple shortest paths are present, the one that also matches the largest number of the following two-qubit gates is chosen. - Sabre: this heuristic routing technique uses a customizable cost function to add SWAP gates that reduce the distance between unconnected qubits involved in two-qubit gates.

Qibolab unroller applies recursively a set of hard-coded gates decompositions in order to translate any gate into single and two-qubit native gates. Single qubit gates are translated into U3, RX, RZ, X and Z gates. It is possible to fuse multiple single qubit gates acting on the same qubit into a single U3 gate. For the two-qubit native gates it is possible to use CZ and/or iSWAP. When both CZ and iSWAP gates are available the chosen decomposition is the one that minimizes the use of two-qubit gates.

Multiple transpilation steps can be implemented using the qibo.transpiler.pipeline.Pipeline:

import networkx as nx

from qibo import gates
from qibo.models import Circuit
from qibo.transpiler.pipeline import Passes, assert_transpiling
from qibo.transpiler.optimizer import Preprocessing
from qibo.transpiler.router import ShortestPaths
from qibo.transpiler.unroller import Unroller, NativeGates
from qibo.transpiler.placer import Random

# Define connectivity as nx.Graph
def star_connectivity():
    Q = [i for i in range(5)]
    chip = nx.Graph()
    chip.add_nodes_from(Q)
    graph_list = [(Q[i], Q[2]) for i in range(5) if i != 2]
    chip.add_edges_from(graph_list)
    return chip

# Define the circuit
circuit = Circuit(2)
circuit.add(gates.H(0))
circuit.add(gates.CZ(0, 1))

# Define custom passes as a list
custom_passes = []
# Preprocessing adds qubits in the original circuit to match the number of qubits in the chip
custom_passes.append(Preprocessing(connectivity=star_connectivity()))
# Placement step
custom_passes.append(Random(connectivity=star_connectivity()))
# Routing step
custom_passes.append(ShortestPaths(connectivity=star_connectivity()))
# Gate decomposition step
custom_passes.append(Unroller(native_gates=NativeGates.default()))

# Define the general pipeline
custom_pipeline = Passes(custom_passes, connectivity=star_connectivity(), native_gates=NativeGates.default())

# Call the transpiler pipeline on the circuit
transpiled_circ, final_layout = custom_pipeline(circuit)

# Optinally call assert_transpiling to check that the final circuit can be executed on hardware
# For this test it is necessary to get the initial layout
initial_layout = custom_pipeline.get_initial_layout()
assert_transpiling(
    original_circuit=circuit,
    transpiled_circuit=transpiled_circ,
    connectivity=star_connectivity(),
    initial_layout=initial_layout,
    final_layout=final_layout,
    native_gates=NativeGates.default()
)

In this case circuits will first be transpiled to respect the 5-qubit star connectivity, with qubit 2 as the middle qubit. This will potentially add some SWAP gates. Then all gates will be converted to native. The qibo.transpiler.unroller.Unroller transpiler used in this example assumes Z, RZ, GPI2 or U3 as the single-qubit native gates, and supports CZ and iSWAP as two-qubit natives. In this case we restricted the two-qubit gate set to CZ only. The final_layout contains the final logical-physical qubit mapping.