How to do measurements in PennyLane

PennyLane offers great resources to design your own gate-based quantum circuits and explore the capabilities of quantum technologies. One of its noteworthy features is that we can either run our circuit on real quantum machines or simulate it classically using our own computer. Even though the future goal is, of course, to use actual quantum computers, the second option allows us to “cheat”. Since we are computing the wavefunction analytically, it gives us access to quantities that would otherwise be out of reach. As we will see, using classical simulations can make our lives much easier if our goal is to figure out if the circuit we are trying to build actually works!

In this how-to, we will use a simple quantum circuit as an example to explain when and how to implement the different types of measurements that are available in PennyLane.

Suppose that we are interested in building a quantum circuit with a specific functionality. There are two essential questions that one needs to ask. We first need to figure out the sequence of gates that compose our algorithm. But this is not the end of the story. It is also very important to define the output of our circuit! In particular, we need to figure out the observable of interest that we want to measure and what kind of information we want to get from this observable.

We will start by building the following simple quantum circuit with two qubits:

In order to do this, we first need to import PennyLane

import pennylane as qml

and then initialize a device with two qubits.

dev = qml.device("default.qubit", wires=2)

Now we are ready to build the circuit and explore its possible outcomes.

Expectation value and variance of an observable

Our circuit in PennyLane is defined as a function that takes a classical input \theta , implements a sequence of quantum operations (one operation per line) and, as highlighted in this how-to, returns a measurement function as an output.

@qml.qnode(dev)
def circuit_expval(theta):
    qml.RX(theta, wires=0)
    qml.Hadamard(wires=1)
    qml.CNOT(wires=[0, 1])
    return qml.expval(qml.PauliY(0))

Here we have chosen to return the expectation value qml.expval() of the operator Y applied to the first qubit. More specifically, this means that we are computing \langle \psi \vert Y \vert \psi \rangle , where \vert \psi \rangle is the 2^n -dimensional vector generated by the circuit (n being the number of wires) and Y is an operator that acts on the first qubit only. When we run the circuit for \theta = 0.5 , we get a number which corresponds exactly to this expectation value as an output.

>>> circuit_expval(0.5)
-0.4794255386042029

A similar but much less common output is the variance of an observable, qml.var().

@qml.qnode(dev)
def circuit_var(theta):
    qml.RX(theta, wires=0)
    qml.Hadamard(wires=1)
    qml.CNOT(wires=[0, 1])
    return qml.var(qml.PauliY(0))

Obtaining again a number as an output.

>>> circuit_var(0.5)
0.7701511529340698

In most cases, we will be interested in obtaining expectation values and variances of observables. Therefore these two measurements will do the job. Nevertheless, this is not the only information we can extract from the circuit.

Samples and probabilities

In a real experiment, we usually obtain a set of actual measurements sampled from the wave function probabilities. In that case, the actual wave function is not accessible, but again, we can “cheat” because we are simulating our circuit classically.

Therefore, if we want to mimic a real device, we can output a set of measurements from the resulting state. To do this, we use the function qml.sample(), for which we need to specify the shots (number of samples) when defining the device.

dev = qml.device("default.qubit", wires=2, shots=4)

Then we can build the circuit

@qml.qnode(dev)
def circuit_sample(theta):
    qml.RX(theta, wires=0)
    qml.Hadamard(wires=1)
    qml.CNOT(wires=[0, 1])
    return qml.sample(qml.PauliY(0))

and obtain the measurement values as output.

>>> circuit_sample(0.5)
array([ 1.,  1.,  1., -1.])

Here we could go a bit further and ask how we obtained these values. What PennyLane is doing is simply taking the eigenvalues \{\lambda_i\} of the observable and sampling them with their associated probabilities p(\lambda_i) = \vert \langle \xi_i \vert \psi \rangle \vert^2 . Here, \vert \xi_i \rangle is the corresponding basis state in the basis of the observable. It is important to note that the sampling measurement is the only one that is not auto-differentiable, since it uses a stochastic process!

This measurement has a bit of a different functionality if we do not provide an observable (return qml.sample()). In this case, the wavefunction will be sampled in the computational basis, with \{0, 1\} as possible eigenvalues. Additionally, since we did not specify which wires to measure, the output will contain one sample for each qubit.

>>> circuit_sample(0.5)
tensor([[0, 1],
        [0, 0],
        [1, 1],
        [0, 0]], requires_grad=True)

We can also specify a set of wires and then the wave function will be sampled only in the computational basis of the qubits that we have selected.

Instead of sampling from the probability density of the wave function, we might be interested in directly obtaining the probabilities of each computational basis state. Well, then PennyLane allows us to do this using the qml.probs() function.

@qml.qnode(dev)
def circuit_probs(theta):
    qml.RX(theta, wires=0)
    qml.Hadamard(wires=1)
    qml.CNOT(wires=[0, 1])
    return qml.probs(wires=[0, 1])

If we reset the shots argument to None (for analytical results),

dev = qml.device("default.qubit", wires=2, shots=None)

the function will return the following array with 2^n numbers, where n is the amount of wires that we input.

>>> circuit_probs(0.5)
tensor([0.46939564, 0.46939564, 0.03060436, 0.03060436], requires_grad=True)

Note that the order of the basis states is lexicographic, which in this case corresponds to \{\vert 00 \rangle, \vert 01 \rangle, \vert 10 \rangle, \vert 11 \rangle\} . For example, we see that the probability of obtaining the state \vert 10 \rangle is approximately 3%.

If we instead specify the number of shots when defining the device, the measurement will result in “noisy” probabilities given by the sampling process, in contrast to the case above in which we analytically compute the exact probabilities from the amplitudes of the final state.

Moreover, if we input an observable instead, we obtain the probabilities of the computational basis states rotated by that observable. For example, the circuit

@qml.qnode(dev)
def circuit_probs_H(theta):
    qml.RX(theta, wires=0)
    qml.Hadamard(wires=1)
    qml.CNOT(wires=[0, 1])
    return qml.probs(op=qml.Hadamard(wires=0))

will return

>>> circuit_probs_H(0.5)
tensor([0.81027229, 0.18972771], requires_grad=True)

which corresponds to a probability of 81% of obtaining the rotated \vert 0 \rangle state when measuring the first qubit and 19% of obtaining the rotated \vert 0 \rangle state.

State and density matrix

Finally, the qml.state() measurement gives us the most complete output possible, this is, the full wave function in the computational basis. It accepts no observables and returns a pure state. The returned array is, as before, in lexicographic order.

@qml.qnode(dev)
def circuit_state(theta):
    qml.RX(theta, wires=0)
    qml.Hadamard(wires=1)
    qml.CNOT(wires=[0, 1])
    return qml.state()
>>> circuit_state(0.5)
tensor([0.68512454+0.j, 0.68512454+0.j,
        0.-0.17494102j, 0., -0.17494102j],
        requires_grad=True)

Sometimes, instead of obtaining the state as a wave function, we might be interested in obtaining the quantum density matrix in the computational basis. Then, we can use qml.density_matrix(). In particular, this is interesting because by specifying the wires argument, we can trace out a part of the system (as in qml.density_matrix([0])). This can result in a mixed state, which can be only represented by the reduced density matrix.

@qml.qnode(dev)
def circuit_density_matrix(theta):
    qml.RX(theta, wires=0)
    qml.Hadamard(wires=1)
    qml.CNOT(wires=[0, 1])
    return qml.density_matrix(wires=0)
>>> circuit_density_matrix(0.5)
tensor([[0.93879128+0.j, 0.+0.23971277j],
        [0.-0.23971277j, 0.06120872+0.j]], requires_grad=True)

---

In this how-to, we have seen that PennyLane offers much more than expectation values of observables when it comes to returning the output of a function and that sometimes choosing the right measurement function can make our lives much easier!