Sampling and Estimation

6.1 Execution Primitives

The rotation prepares \(R_y(\theta)\ket0=\cos(\theta/2)\ket0+\sin(\theta/2)\ket1\). Two natural questions about it have different answers. The probability of measuring one is \(p(1)=\sin^2(\theta/2)\). The average of \(Z\) is \(\cos^2(\theta/2)-\sin^2(\theta/2)=\cos\theta\). Which interface should answer each question, and where in its result do we find the number? We take the interfaces first and their results second.

6.1.1 V2 Sampler and Estimator

We already have a sampler, and from its counts we could compute \(\langle Z\rangle\) ourselves. Why would we want a separate interface for averages? For one qubit and one observable, we would not. But an observable such as \(-ZZ-XX\) needs measurements in more than one basis, and larger observables need many. Choosing those bases and combining their results is work that the software, or a remote service, can do for us, and can do with techniques of its own.

So Qiskit offers two primitives, the standard interfaces we met as reference primitives. A Sampler takes measured circuits and returns their outcomes. An Estimator takes circuits without measurements, together with observables, and returns expectation values. Their current versions are called V2.

Both accept many requests in one call. Each request is written as a tuple, a fixed-order group of values in parentheses, such as (qc, SparsePauliOp("Z"), values). For an Estimator the order is circuit, observable, parameter values. Qiskit calls one such request a primitive unified bloc, or PUB. A PUB is a grouping of inputs, nothing more: it is not a gate, and it does not promise that its parts run together on a device. The companion uses the local StatevectorSampler and StatevectorEstimator, so no account is involved.

6.1.2 Local implementations and result structures

Once a run finishes, where is our number? The result holds one entry per PUB, in the order we sent them, so result[0] belongs to the first PUB. For an Estimator, result[0].data.evs is an array of expectation values, one for each set of parameter values in that PUB.

For a Sampler, result[0].data is organized by classical register name. We have always used measure_all(), which creates a register named meas, so we read data.meas. A circuit with its own register named c would be read as data.c. So an access expression copied from another program can fail, or read the wrong register, unless we check the register names first.

6.2 Statistical Readout

The local Estimator computes \(\cos\theta\) exactly. A device cannot: it repeats the circuit a finite number of times and averages. How far from the true value should we expect such an estimate to land? And when we only have bit strings, how do we turn them into an average ourselves? We derive the uncertainty first, then the decoding.

6.2.1 Shots, precision, and uncertainty

Suppose we run the circuit \(R\) times at \(\theta=\pi/2\), where \(\langle Z\rangle=0\). How close to 0 will our estimate be?

Each shot gives a bit, and we record \(s=+1\) for outcome 0 and \(s=-1\) for outcome 1, the two eigenvalues of \(Z\). The mean of \(s\) is \(p(0)-p(1)\), and since \(p(0)=1-p(1)\), \[\langle Z\rangle=1-2p(1).\] Because \(s^2=1\), the variance of one shot is \(\mathbb E[s^2]-\langle Z\rangle^2=1-\langle Z\rangle^2\). As with the coin, averaging \(R\) independent shots divides the variance by \(R\): \[\operatorname{Var}\bar s=\frac{1-\langle Z\rangle^2}{R},\qquad \operatorname{Var}\hat p(1)=\frac{p(1)\bigl(1-p(1)\bigr)}{R}.\] The two agree: \(\bar s=1-2\hat p(1)\), so its variance is four times that of \(\hat p(1)\), and \(4p(1-p)=1-(1-2p)^2\).

At \(\theta=\pi/2\) with \(R=4096\), the standard deviation of \(\bar s\) is \(1/64\approx0.016\). This standard deviation of an estimate is called its standard error. Doubling the shots divides it by \(\sqrt2\), so each extra digit of accuracy costs a hundred times more shots.

An Estimator lets us ask for a target standard error, its precision. StatevectorEstimator uses precision 0 by default and then returns the exact value. With a positive precision it adds random noise of about that size, to imitate a finite-shot estimate. Either way it is not a model of a noisy device.

6.2.2 Bit strings and classical postprocessing

Now suppose we have two-qubit counts and want \(\langle ZZ\rangle\). How do we compute it? The Z outcome of qubit \(j\) is \(+1\) when its bit is 0 and \(-1\) when it is 1. \(ZZ\) is the product of the two, so it is \(+1\) when the bits are equal and \(-1\) when they differ. Hence \[\widehat{\langle ZZ\rangle}=\frac{N_{\text{equal}}-N_{\text{unequal}}}{R}.\]

For \(ZI\), by contrast, only the left character, qubit 1, matters. So we must decode each string with the measurement mapping before computing. A mapping that swaps the bits leaves \(ZZ\) unchanged but exchanges \(ZI\) and \(IZ\).

What about small entries of a count table? Dividing counts by \(R\) gives an empirical distribution. An outcome with probability \(p\) has standard error \(\sqrt{p(1-p)/R}\), which relative to \(p\) is about \(1/\sqrt{pR}\). At \(p=0.001\) and \(R=4096\) we expect about 4 counts, and the relative uncertainty is about 50%. Small entries are the least reliable.

6.3 Repeated Workloads

We usually want \(\langle Z\rangle\) not at one angle but at many. Must we build and send one circuit per angle? A symbolic parameter lets one circuit serve every angle. What remains is to describe the list of angles so that the primitive runs the experiments we mean, and to keep their results in order.

6.3.1 Parameter sweeps

How do we ask for \(\theta=0,\pi/2,\pi\) in one PUB? The companion writes values = [[0], [np.pi / 2], [np.pi]], a list of lists. Each inner list is one row, holding one value for each parameter of the circuit. Our circuit has one parameter, so each row has one entry, and three rows ask for three experiments. The array has shape \(3\times1\), and the predicted results are \(\cos\theta=1,0,-1\).

Shapes decide which experiments run, so a mistake can change the workload without any warning. Suppose we also want \(\langle X\rangle\) at the same three angles. Passing the two observables as a flat list, of shape 2, beside three rows, of shape 3, raises an error, because the two shapes cannot be matched. Passing them as [[Z], [X]], of shape \(2\times1\), gives a \(2\times3\) result: every observable at every angle. The Estimator combines the shapes by NumPy’s broadcasting rule, which stretches a dimension of length 1 to match the other array. So we check the shape of evs before trusting the order of its entries.

6.3.2 Batching inputs and interpreting results

When we have several PUBs, why send them in one call rather than one at a time? Each call has overhead, and on a remote service each job may wait in a queue. Sending several PUBs in one call, called batching, saves that overhead. It does not reduce the work: every circuit, parameter row, observable and shot still costs what it costs. Since results come back in PUB order, we keep the parameter order and the result order together in the saved record.

The companion keeps its two sources of output apart. First the Estimator computes the exact three-point sweep and checks it with np.allclose. The middle value prints as about \(2\times10^{-16}\) rather than 0, a rounding error of the kind tolerances exist for. Then the program binds \(\theta=\pi/2\), adds measurements and samples 4096 shots, so the counts should be close to an even split.

"""Runnable companion for Quantum Programming, Chapter 6."""
import numpy as np
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")
qc = QuantumCircuit(1)
qc.ry(theta, 0)
values = [[0], [np.pi / 2], [np.pi]]
result = StatevectorEstimator().run([(qc, SparsePauliOp("Z"), values)]).result()
assert np.allclose(result[0].data.evs, [1, 0, -1])
measured = qc.assign_parameters({theta: np.pi / 2})
measured.measure_all()
counts = StatevectorSampler(seed=9).run([measured], shots=4096).result()[0].data.meas.get_counts()
print("Exact means:", result[0].data.evs, "sampled counts:", counts)

Companion script: chapter_06.py.

With seed 9 the sampled counts were 2051 zeros and 2045 ones, which gives \(\bar s=6/4096\approx0.0015\), well within one standard error of the exact value 0.

Every circuit so far acted on inputs we wrote in by hand. Many algorithms instead assume an oracle, a circuit that evaluates a function on whatever input it receives. How do we build one, rather than assume it exists?

Sources and further study: IBM, StatevectorEstimator; IBM, StatevectorSampler.