States and Operators

4.1 State Representations

A simulator can show the state hidden inside a circuit, which makes it a debugging instrument: we can compare the program with its predicted action before sampling or noise enter. What form should that state take? A list of amplitudes is the natural first answer. We will find two situations it cannot describe, a random preparation and one part of an entangled pair, and one representation that handles both.

4.1.1 State vectors and density matrices

What exactly does Statevector.from_instruction(qc) compute? It starts every qubit in \(\ket0\), applies the circuit’s gates, and returns the \(2^n\) amplitudes in index order, with index \(k=q_0+2q_1+\cdots\). That works only for a circuit with no measurement or reset, since those are not unitary steps.

Now recall the device that prepares \(\ket{00}\) or \(\ket{11}\) with probability one half each. Can some single vector \(\psi\) describe it? To decide, we need a form in which predictions are easy to average. The average of an observable \(A\) in a pure state is \(\psi^\dagger A\psi\), and we can rewrite this using the trace, the sum of a matrix’s diagonal entries: \[\psi^\dagger A\psi=\sum_{j,k}\psi_j^*A_{jk}\psi_k=\operatorname{tr}\!\left(A\,\psi\psi^\dagger\right).\] A random preparation averages its predictions, so it gives \(\operatorname{tr}(A\rho)\) with \[\rho=\tfrac12\ketbra{00}{00}+\tfrac12\ketbra{11}{11}.\]

The matrix \(\rho\) is called a density matrix. For a mixture of states \(\psi_i\) with probabilities \(p_i\), it is \(\rho=\sum_ip_i\,\psi_i\psi_i^\dagger\), and every prediction is \(\operatorname{tr}(A\rho)\).

Is our \(\rho\) secretly of the form \(\psi\psi^\dagger\)? Any unit vector satisfies \((\psi\psi^\dagger)^2=\psi(\psi^\dagger\psi)\psi^\dagger=\psi\psi^\dagger\). Our \(\rho\) instead gives \(\rho^2=\tfrac14\ketbra{00}{00}+\tfrac14\ketbra{11}{11}=\rho/2\). So no vector describes the mixture, and we need the more general object. Qiskit’s DensityMatrix class holds it, and DensityMatrix(state) turns a vector into \(\psi\psi^\dagger\).

4.1.2 Composite systems and subsystem order

Density matrices settle a question vectors could not. In a Bell pair, what can we say about qubit 0 alone? Any test on qubit 0 alone uses an observable of the form \(I\otimes A\), with \(A\) on the right because qubit 0 is the right factor. For the Bell state \(\Phi=(\ket{00}+\ket{11})/\sqrt2\), \[\Phi^\dagger(I\otimes A)\Phi=\tfrac12\left(\bra0A\ket0+\bra1A\ket1\right)=\operatorname{tr}\!\left(A\,\tfrac I2\right).\] So every test on qubit 0 behaves as if that qubit had the density matrix \(I/2\). Since \((I/2)^2=I/4\), no vector describes it either.

The operation that produces a subsystem’s density matrix is the partial trace. To discard qubit 1, we sum over its label: \[(\rho_0)_{ab}=\sum_{c}\bra{c\,a}\rho\ket{c\,b},\] where the left label \(c\) belongs to qubit 1. For the Bell pair only \(c=a=b=0\) and \(c=a=b=1\) contribute, each giving \(\tfrac12\), so \(\rho_0=I/2\). In Qiskit, partial_trace(state, [1]) takes the list of subsystems to discard, not those to keep.

Is \(I/2\) just another way of writing \(\ket+\)? Both give zero and one with probability one half. But \[\ketbra++=\begin{pmatrix}\tfrac12&\tfrac12\\\tfrac12&\tfrac12\end{pmatrix},\qquad \frac I2=\begin{pmatrix}\tfrac12&0\\0&\tfrac12\end{pmatrix}.\] The average of X is \(\operatorname{tr}(X\rho)=\rho_{01}+\rho_{10}\), which is 1 for \(\ket+\) and 0 for \(I/2\). So an X-basis test tells them apart.

4.2 Operators and Observables

A state tells us what the circuit produced from one input, all zeros. A wrong circuit could still produce the right state from that input. And the averages we compared were computed from states one test at a time. We therefore need operators of two kinds: the circuit’s whole matrix, which covers every input, and the observables whose averages we predict.

4.2.1 Unitary operators

How can we check a circuit on every input at once? Its matrix does that, because column \(k\) of the matrix is the output for basis input \(k\). Operator(qc) builds this matrix for a circuit with no measurement or reset. For it to be a valid quantum evolution, it must be unitary, \(U^\dagger U=I\), and the companion checks this.

The check uses NumPy, imported as np. NumPy stores matrices as arrays and supplies matrix arithmetic. Operator(qc).data is the matrix as such an array. The method .conj() takes complex conjugates and .T transposes, so .conj().T is the dagger. The operator @ multiplies matrices, and np.eye(4) is the \(4\times4\) identity.

Why not test equality with ==? Computers store \(1/\sqrt2\) only approximately. In Python, (2**-0.5)**2 gives 0.5000000000000001, not 0.5, so an exact comparison can fail for a correct circuit. np.allclose instead accepts entries that agree within a small tolerance. Its default tolerance is far tighter than any mistake we are looking for, and far looser than these rounding errors.

How far does this check scale? An \(n\)-qubit unitary is a \(2^n\times2^n\) matrix, so it has \(2^n\cdot2^n=4^n\) complex entries. At 20 qubits that is about \(1.1\times10^{12}\) entries, roughly 18 terabytes at 16 bytes each. A full matrix is a check for small circuits only.

4.2.2 Pauli observables and expectation values

When we tested the Bell pair with H on both qubits, we inferred the result from counts. Can we predict such averages directly? The observables we need are Pauli strings: tensor products of \(I,X,Y,Z\), written as letters. Following the bit-order convention, the rightmost letter acts on qubit 0, so ZI means \(Z\) on qubit 1 and nothing on qubit 0.

Each Pauli string has eigenvalues \(\pm1\), so a single measurement of it gives \(+1\) or \(-1\), and its average lies between \(-1\) and \(1\). For the Bell state, the two Z outcomes always agree, so their product is always \(+1\) and \(\langle ZZ\rangle=1\). Since \(XX\) swaps \(\ket{00}\) with \(\ket{11}\), it leaves the Bell state unchanged, so \(\langle XX\rangle=1\) as well. Qubit 1 alone gives each Z outcome half the time, so \(\langle ZI\rangle=0\). Measuring \(XX\) is the H-on-both test, since \(HXH=Z\): we apply H to each qubit and read the product of the Z outcomes.

Qiskit’s SparsePauliOp stores an observable as a list of Pauli strings with coefficients, such as \(-ZZ-XX\), without building the \(4^n\)-entry matrix. The method state.expectation_value(op) computes \(\langle\psi|\text{op}|\psi\rangle\) exactly. An experiment cannot do that. It estimates the same number by repeating measurements and averaging.

4.3 Simulation

We now hold three kinds of numbers for one circuit: exact amplitudes and density matrices, exact expectation values, and sampled counts. They come from different calculations and mean different things. When a result disagrees with a prediction, which of them should we doubt, and which could a laboratory even produce? We compare exact and sampled numbers first, then ask what physical measurement can recover.

4.3.1 Exact evolution and sampled results

Suppose a test reports \(\langle ZI\rangle=0.03\) instead of 0. Is the program wrong? That depends on where the number came from. An exact expectation from expectation_value checks the ideal mathematical action, so a 0.03 there would signal a bug. A mean of \(R\) sampled \(\pm1\) values fluctuates, with a standard deviation near \(1/\sqrt R\) when the true mean is 0, so at \(R=1024\) a value of 0.03 is ordinary.

Later we will add a third layer, a simulation of a noisy device, whose results differ from both. If we keep exact, sampled and noisy results separate, a discrepancy points to one cause: the code, the sampling, or the noise model.

4.3.2 Simulator access versus physical measurement

The companion prints the density matrix of qubit 0. Which of its numbers could an experiment recover, and how?

A single run of a device gives one bit string, never an amplitude. Repeated runs estimate probabilities, and these are the diagonal entries of a density matrix. The off-diagonal entries take more work. Writing a one-qubit density matrix as \(\rho=(I+xX+yY+zZ)/2\), with \(x=\langle X\rangle\) and so on, gives \(\rho_{01}=(x-iy)/2\). So averages in the X and Y bases recover them. This reconstruction from several measurement bases is called tomography, and for \(n\) qubits it needs \(4^n-1\) averages, so it is practical only for small systems. A global phase is never recoverable at all.

A simulator, by contrast, hands us every number for free. That is its value as a debugging tool, and also why it can mislead us about what an experiment will show. The companion checks the Bell correlations, the reduced state and unitarity, all with tolerances for floating-point rounding.

"""Runnable companion for Quantum Programming, Chapter 4."""
import numpy as np
from qiskit import QuantumCircuit
from qiskit.quantum_info import Statevector, DensityMatrix, Operator, SparsePauliOp, partial_trace

qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
state = Statevector.from_instruction(qc)
local = partial_trace(DensityMatrix(state), [1])
assert np.allclose(local.data, np.eye(2) / 2)
assert np.isclose(state.expectation_value(SparsePauliOp("XX")), 1)
assert np.isclose(state.expectation_value(SparsePauliOp("ZI")), 0)
assert np.allclose(Operator(qc).data.conj().T @ Operator(qc).data, np.eye(4))
print("Local state:", local.data)

Companion script: chapter_04.py.

To extend the check, replace H by a rotation \(R_y(\theta)\). The circuit then prepares \(\cos(\theta/2)\ket{00}+\sin(\theta/2)\ket{11}\). Predict qubit 0’s density matrix before running: \(\operatorname{diag}(\cos^2(\theta/2),\sin^2(\theta/2))\). Its off-diagonal entries vanish after the CNOT, even though the two-qubit state is pure.

Every circuit so far was written out gate by gate. When a construction is used many times, at different angles or on different wires, how do we write it once and trust each use?

Sources and further study: IBM, quantum information API.