Simulation of quantum circuits on the QLM: introduction

This notebook aims at introducing how simulations are run in the QLM. It is valid for both noisy and ideal circuit simulators.

• When simulating a quantum circuit, one can aim for different kinds of results:
  
  

  • You might want to get the full amplitude vector back, e.g for detailed analysis and debugging of your circuit implementation.
    
    
  • On the contrary, you may be interested in strictly emulating the behavior of a quantum computer, and getting a list of measurement results as a result of your computation.
    
    
  • Finally, perhaps you are only interested in the average value of an observable at the output of your circuit. It might be useful when dealing, for instance, when dealing with hybrid variational approaches (QAOA, VQE...).
    
    

• Table of contents of this notebook

  • Introduction of the process, with default options
  • Working with a subset of qubits
  • Getting a list of measurement results
  • Grouping measurement results by value: aggregate_data
  • Getting the average value of an observable
  • A little hack: directly getting the table of amplitudes as a numpy array

Within this notebook, we only work with a simple Bell-state-creation quantum circuit, and simulate it with our generic simulator PyLinalg (based on linear algebra).

For more details about how to write observables see this notebook and the Sphinx Documentation

Overall process

A simulation is started by sending a job to a qpu (quantum processing unit) via its submit method.

In our case, a qpu is a simulator. The job is created from a quantum circuit. Simulation options are specified at the creation of this job.

The following snippet example of the process:

from qat.lang.AQASM import Program, H, CNOT
from qat.qpus import PyLinalg

# qpu creation
qpu = PyLinalg() 

# program creation and gate applications
my_prog = Program()
qbits = my_prog.qalloc(2)
my_prog.apply(H,qbits[0])
my_prog.apply(CNOT,qbits)

# converting into a circuit
circ = my_prog.to_circ()

job = circ.to_job() # specify simulation options. Here: default.

result = qpu.submit(job)

for state in result:
    print(state)

Sample(_amplitude=ComplexNumber(re=0.7071067811865475, im=0.0), probability=0.4999999999999999, _state=b'\x00', err=None, intermediate_measurements=None, qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=ComplexNumber(re=0.7071067811865475, im=0.0), probability=0.4999999999999999, _state=b'\x03', err=None, intermediate_measurements=None, qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])

What happened here ?

"result" contains all states with non-zero amplitude.

The job was created with default arguments. Which in particular that the "nbshots" argument was equal to 0 (see docstring below).

• Concerning nbshots, the convention regarding the behavior depending on its value is:
  
  
  • if nbshots = 0 then the qpu returns the best it can do. In the case of LinAlg the best it can do is just to return the full state distribution. Because states are, by default, filtered with an amplitude threshold (see docstring below), the result here consists in the two possible outputs, for a Bell state. It might be the case, for instance with an actual quantum chip, that "the best the qpu can do" is to return a list of measurement result, as long as reasonably possible.
    
    
  • if nbshots > 0 then the qpu returns a list of samples, obtained by measuring the output probability distribution of the circuit. As we will see later, in that case, the format of the output also depends on the "aggregate_data" argument.
    
    

As you can see, it is also at job creation that one can specify the subset of qubits of interest. Apart from the number of qubits in the returned samples, as we will see, the behavior is strictly identical to the default case where all qubits are measured.

help(circ.to_job)

Help on method to_job in module qat.core.wrappers.circuit:

to_job(job_type='SAMPLE', qubits=None, nbshots=0, aggregate_data=True, amp_threshold=9.094947017729282e-13, **kwargs) method of qat.core.wrappers.circuit.Circuit instance
    Generates a Job containing the circuit and some post
    processing information.

    Args:
        job_type (str): possible values are "SAMPLE" for computational basis sampling
            of some qubits, or "OBS" for observable evaluation
            (see :py:mod:`qat.core.Observable`
            for more information about this mode).
        qubits (optional, list<int>, list<QRegister>): the list of qubits to measure
            (in "SAMPLE" mode).
            If some quantum register is passed instead, all the qubits of the register will be
            measured. Moreover, if the register was "typed"
            (see :meth:`qat.lang.Program.qalloc`
            for more information), the results will be cast into the register type.
            Defaults to None,
            meaning all qubits are to be measured.
        nbshots (optional, int): The number of shots to perform. Defaults to zero.
            If set to zero,
            the convention is that the QPU should do its best: a quantum processor will
            use the largest
            amount of shot authorized by its configuration, while a simulator
            will try to output
            all the possible states constained in the final distribution, together with their
            probabilities (and possible amplitudes).
        aggregate_data (optional, bool): if set to True, and `nbshots` is not zero,
            the samples will be
            aggregated and their probability field will be used to store their
            observed frequencies of apparition. Defaults to True.
        amp_threshold (optional, double): amplitude threshold under which states are
            not returned in the result structure. Useful to prune states that
            are unlikely to be measured out of the
            returned samples. Defaults to 1/2^40.

    Keyword Args:
        observable (:class:`~qat.core.Observable`): see
            :class:`~qat.core.Observable`. Used for the
            "OBS" mode only).

We are now going to play around with all these arguments, in order to illustrate the behaviors they lead to.

Focusing on a subset of qubits

job = circ.to_job(qubits=[0]) # focusing on the first qubit. 

result = qpu.submit(job)

for state in result:
    print(state) # 1-qubit states, because we focus on one qubit.

Sample(_amplitude=None, probability=0.4999999999999999, _state=b'\x00', err=None, intermediate_measurements=None, qregs=[DefaultRegister(key=None, length=1, start=0, msb=None, _subtype_metadata=None)])
Sample(_amplitude=None, probability=0.4999999999999999, _state=b'\x01', err=None, intermediate_measurements=None, qregs=[DefaultRegister(key=None, length=1, start=0, msb=None, _subtype_metadata=None)])

Or, equivalently, but more general (e.g when working with several registers of qubits):

job = circ.to_job(qubits=[qbits[0].index]) # focusing on the first qubit. the index contains the 
                                           # index of the qubit within the entire set of qubits,    
                                           # composed of all registers.
result = qpu.submit(job)

for state in result:
    print(state)

Sample(_amplitude=None, probability=0.4999999999999999, _state=b'\x00', err=None, intermediate_measurements=None, qregs=[DefaultRegister(key=None, length=1, start=0, msb=None, _subtype_metadata=None)])
Sample(_amplitude=None, probability=0.4999999999999999, _state=b'\x01', err=None, intermediate_measurements=None, qregs=[DefaultRegister(key=None, length=1, start=0, msb=None, _subtype_metadata=None)])

Strict emulation: getting a list of samples from the output probability distribution.

By "strict emulation", we mean that, even though here we are using a classical simulator to process our quantum circuits, the result looks exactly like the raw output of an actual quantum computer.

job = circ.to_job(nbshots=10, aggregate_data=False) # see below for aggregate_data.

result = qpu.submit(job)

for state in result: # 10 results are printed, as required.
    print(state) 

Sample(_amplitude=None, probability=None, _state=b'\x00', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x00', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x00', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x03', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x03', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x03', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x00', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x03', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x03', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=None, _state=b'\x00', err=None, intermediate_measurements=[], qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])

Grouping results by value: aggregate_data option (defaults to True)

When this option is active (which is the default), all samples containing a given state are "aggregated" together. The result object then contains one unique sample per possible output. It comes with an empirical estimation of the probability of that output, and an "error" associated to that estimation.

job = circ.to_job(nbshots=10) # see below for aggregate_data.

result = qpu.submit(job)

for state in result: # 10 results are printed, as required.
    print(state) 

Sample(_amplitude=None, probability=0.4, _state=b'\x00', err=0.1632993161855452, intermediate_measurements=None, qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=0.6, _state=b'\x03', err=0.1632993161855452, intermediate_measurements=None, qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])

Notice how the error decreases when taking a greater number of samples: TODECIDE: to we remove this or do we implement the error estimation for linalg ?

job = circ.to_job(nbshots=1000) # see below for aggregate_data.

result = qpu.submit(job)

for state in result: # 10 results are printed, as required.
    print(state) 

Sample(_amplitude=None, probability=0.479, _state=b'\x03', err=0.015805341148131185, intermediate_measurements=None, qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])
Sample(_amplitude=None, probability=0.521, _state=b'\x00', err=0.015805341148131185, intermediate_measurements=None, qregs=[QRegister(scope=<qat.lang.AQASM.program.Program object at 0x7f68a81353d0>, length=2, start=0, msb=None, _subtype_metadata=None, qbits=[<qat.lang.AQASM.bits.Qbit object at 0x7f68a8134fb0>, <qat.lang.AQASM.bits.Qbit object at 0x7f68a81365a0>])])

Average value of an observable

Often, for instance when dealing with hybrid variational algorithms, we are actually interested in the average value of an observable, and not measurement values.

Here, we only give a simple example of the "ZZ" observable, which is always equal to 1 for a bell pair. For more advanced used of observables, and details regarding their construction, we refer the reader to Sphinx Documentation and this notebook.

from qat.core import Observable, Term

obs = Observable(2., pauli_terms=[Term(1., "ZZ", [0, 1])])
job = circ.to_job("OBS", observable=obs)

result = qpu.submit(job)

print("Should be 1, as we are working with a Bell pair: ", result.value)

Should be 1, as we are working with a Bell pair:  0.9999999999999998

Getting the tables of amplitude directly, as a numpy array: wavefunction

from qat.core.simutil import wavefunction

qpu = PyLinalg() 

wf = wavefunction(circ, qpu)

print(wf) # simple and efficient.

[0.70710678+0.j 0.        +0.j 0.        +0.j 0.70710678+0.j]