Quantum simulation

Control readings

A control decides whether the gate fires by what it reads, and all ones is only the default. control_values says what each control has to hold, in wire order:

tq.apply_mcx(tq.zeros_state(3), [0, 1, 2])                      # fires on |11.>
tq.apply_mcx(tq.zeros_state(3), [0, 1, 2], control_values=[0, 0])  # fires when both are 0

apply_controlled(state, mat, wires, control_values=...) is the same idea for an arbitrary operator, with the controls leading the wire list.

Batched parameters

A gate matrix may carry one entry per batch item, which gives every state in the batch its own parameter:

theta = tp.tensor([0.0, 3.14159])          # per batch item
states = tq.zeros_state(1, batch=2)
tq.apply_ry(states, 0, theta)               # first qubit stays 0, second flips

Ansatz templates

Layer templates lay out a trainable circuit instead of hand by hand. nlocal_circuit is the general form: reps blocks of “rotate every wire, then entangle”, a closing rotation layer, and one trainable angle per rotation gate.

tq.real_amplitudes(3, reps=2)          # RY rotations and CNOT pairs
tq.efficient_su2(3, reps=2)            # RY and RZ rotations, CNOT pairs
tq.nlocal_circuit(3, ("ry", "p"), reps=1, entanglement="circular")

The entanglement pattern is one of "linear", "reverse_linear", "circular" or "full"; entanglement_pairs(n, pattern) returns the pairs directly. Angles arrive as module parameters in (-pi, pi), so the template can go straight to an optimizer:

qc = tq.real_amplitudes(3, reps=2)
loss = tq.expectation(qc.state(), "ZI")
loss.backward()
len(list(qc.parameters()))            # 9: three wires, two blocks, one closing layer

Reading values back out of shots

A Pauli word that touches an X or a Y wire cannot be read off a plain histogram until the wires are turned: apply_measurement_basis puts the state into the frame where measuring in Z reports each factor’s eigenvalue. The histogram then averages those signs, which is expectation_from_counts.

readable = tq.apply_measurement_basis(state, "XZ")
counts = tq.measure(readable, shots=4096, wires=[0, 2])
tq.expectation_from_counts(counts, "XZ", wires=[0, 2])   # [batch]

Words that agree on every wire where both of them act can share one such rotation: each is a sub-word of the group’s product, so a single histogram reads all of them. commuting_groups finds the groups, joint_expectation runs them and returns one estimate per word — fewer settings means fewer shots for the same noise.

tq.commuting_groups(["XZ", "XI", "ZZ"])   # [['XZ', 'XI'], ['ZZ']]
tq.joint_expectation(state, ["XZ", "XI", "ZZ"], shots=4096)   # [batch, 3]

A readout that lies

A device reports a bit that was prepared wrong now and then, and the two directions usually differ. A readout matrix records that: row p is what gets reported for a prepared p, so it is a distribution over the reported bit.

noisy = tq.readout_error_matrix(0.02)               # 2% in both directions
worse = tq.readout_error_matrix(0.02, 0.05)         # 2% of 1s, 5% of 0s

measure takes one and decides every shot’s reading on its own, which is how a real readout behaves:

counts = tq.measure(state, shots=4096, readout=noisy)

Give one matrix per wire as [n, 2, 2] when the wires were calibrated apart. apply_readout_error applies the same map to a histogram or an exact distribution, so an offline histogram can be corrected — or biased, which is what it is usually wanted for.

Reading a state back out of a histogram

A histogram in the computational basis is the diagonal of a density matrix and nothing else. density_from_counts hands that diagonal back, which is the state the measurements cannot rule out and no more: every off-diagonal term is a phase, and a basis measurement never sees a phase.

rho = tq.density_from_counts(counts)                # populations, no phases
state = tq.state_from_counts(counts)                # the square roots of them

state_from_counts names the member of that family a histogram leaves free whose amplitudes are real and non-negative, which is the one overlapping the measured state the most. Comparing it with a theory compares populations; the phases come from histograms taken in rotated bases, which is what joint_expectation reads.

An outcome that never came out of a finite number of shots stays out. Pass pseudocount= to lift every outcome off the floor before normalizing, which makes the density matrix full rank for a tomography fit — and makes the answer wrong by that much, which is what it says on the tin.

A system that does not close

A channel applies all at once. Noise that acts over time is the generator d rho / dt = -i [H, rho] + sum_r (L_r rho L_r^dagger - 1/2 {L_r^dagger L_r, rho}), and Lindbladian is that generator, held as an object and applied as often as an integration needs.

generator = tq.Lindbladian(h2, jump_ops=["XI", ("IZ", 0.05)])
final = generator.evolve(rho, time=2.0)

The generator itself is never built. Written out densely it is 2**n times the size of rho squared – on ten qubits, a matrix of a trillion entries – so each term is applied on its own and the integral is formed by composition alone.

A jump operator is a Pauli word, a (wires, matrix) pair, or a full-register matrix, and a bare Pauli word means a rate of one. A rate can be written next to it, which is the short way to say how fast the error happens:

tq.Lindbladian(h2, [("XI", 0.02), ("ZI", 0.01)])        # Pauli errors at rates
tq.Lindbladian(h2, [([0], 0.05 ** 0.5 * excitation), n_qubits=4])

A word of identities and Zs is a sign pattern rather than a general matrix, so its whole dissipator is one factor per entry of rho and costs a single pass. A stack of jump operators is read one at a time, the same way.

A jump carries its rate in its own size — sqrt(gamma) * L is an error happening at rate gamma, which is the master equation’s own form of the statement — so amplitude_damping_kraus(gamma)[1] is the sqrt(gamma) |0><1| this wants.

A channel is not a generator, and handing one in as if it were is the one way to get a plausible answer that is quietly wrong. A channel is the exponential of a generator: summing dissipators over its Kraus operators applies the channel once and subtracts the state once per operator, which is not the exponential of anything. The population decay comes out right and the coherence decay does not, by a few percent, and nothing looks amiss. So a stack handed in here is read as a list of jumps, and a channel is applied to a density matrix with apply_channel instead — which is where a channel belongs.

Two routes integrate the generator, and both keep the trace:

generator.evolve(rho, 2.0, method="rk4")     # fixed fourth-order steps
generator.evolve(rho, 2.0, method="taylor")  # a truncated series of the span

"rk4" is the one that can follow a Hamiltonian that moves with time, given as a callable of t, and the one that reports intermediate times. "taylor" sums a series of the whole span, which costs several times less per unit of time covered and lands closer to the exact answer — the default wherever it applies, short spans and long ones alike.

No span is out of reach for the series. Its terms climb to about exp(span) while the sum they add to stays the size of the state, so one sum can only carry so far before it cancels itself away in the rounding; past that point the span is covered by several sums of a size one can carry, in a row. Advancing by a then by b is advancing by a + b, each half applied to the matrix the one before left, and the propagator is still never built — a propagator can only be squared as an operator, and this generator is never built to square.

What an experiment wants is the observables along the way:

states, values = generator.trajectory(
    rho, 2.0, n_sample=101, observables=[h2, "IZ"]
)

values is [time, observable, batch], real. Without a sample count, mesolve answers with the final matrix on its own; sesolve is the closed-system counterpart, which advances a state rather than its square and so costs 2**n per step instead of 4**n.

tq.sesolve(state, h2, time=2.0, order=2, n_sample=51, observables=[h2])
tq.mesolve(rho, h2, None, time=2.0, n_sample=51, observables=[h2])   # same answer

Gradients without a graph

Two situations call for something else. One is a value that is not built out of differentiable primitives at all: an energy estimated from a finite number of measurements has nothing to backpropagate through. The other is a circuit deep enough that backpropagating through every rotation costs more than evaluating the whole thing twice.

sample_energy answers the first — the Hamiltonian measured with as few measurement settings as its terms allow, optionally through a readout matrix so the estimate carries the error a device would add:

value = tq.sample_energy(state, h2, shots=4000)

param_shift answers the second, and answers the first too. For a value that depends on one angle through exp(-i θ G), the derivative is a difference of the value at two displaced angles, and it is exact — independent of how deep the circuit is and of what sits in the middle of it, a sampler included:

def energy(angles):
    return tq.energy(tq.apply_cnot(
        tq.apply_ry(tq.apply_ry(tq.zeros_state(2, batch=angles.shape[0]),
                                0, angles[:, 0]),
                    1, angles[:, 1]), 0, 1), h2)

angles = tp.tensor([[0.2, -0.4]])
gradients = tq.param_shift(energy, angles)          # no graph anywhere
sampled = lambda a: tq.sample_energy(build(a), h2, 2000)
noisy = tq.param_shift(sampled, angles)             # a value with no graph either

Every angle gate here is a rotation about half a Pauli word, whose eigenvalues are ±1/2, so the default displacement of π/2 and weight of one half are the rule for all of them. A generator whose eigenvalues differ by ω instead wants π / (2 ω) and a weight of ω / (2 sin(ω shift)); both are arguments, so the rule follows the gate rather than assuming it.

The rule reaches an angle that enters through one rotation and stops there. A controlled rotation is not one exponential — where the control is unset its generator reads zero — so the value carries half the frequency as well as the whole one, and no single displacement reaches both. Those angles want the differentiated gradient, which covers them.

One run, three ways

The three routes are the same run to a variational search: the same circuit, the same Hamiltonian, the same optimizer, and a different way of getting the direction to walk in.

h2 = tq.Hamiltonian([(-1.0523732, "II"), (0.3979374, "IZ"),
                     (-0.3979374, "ZI"), (-0.1809312, "XX")])
exact = float(tp.linalg.eigvalsh(h2.matrix())[0])        # -1.868555

def circuit(angles):                                     # 3 layers, 12 angles
    state = tq.zeros_state(2)
    for layer in range(3):
        for qubit in range(2):
            state = tq.apply_ry(state, qubit, angles[0, 4 * layer + qubit])
            state = tq.apply_rz(state, qubit, angles[0, 4 * layer + qubit + 1])
        state = tq.apply_cnot(state, 0, 1)
    return state

energy = lambda a: tq.energy(circuit(a), h2)
theta = nn.Parameter(tp.full((1, 12), 0.4))
optimizer = tp.optim.Adam([theta], lr=0.1)
for step in range(80):
    optimizer.zero_grad()
    energy(theta).backward()                             # through the circuit
    # theta.backward(tq.param_shift(energy, theta.detach()))          # by displacing
    # theta.backward(tq.param_shift(
    #     lambda a: tq.sample_energy(circuit(a), h2, 4000),
    #     theta.detach()))                                              # and in shots
    optimizer.step()

Backpropagating lands about two thousandths above the exact ground energy, the displaced rule about three millionths — it is exact, so what is left is the optimizer’s own step size — and four thousand shots per evaluation about a hundred-thousandth, which is where the shot noise puts it.

Precision

A register’s own dtype is what a whole circuit runs at: a double-precision state stays double through every gate, and a single-precision one stays single even when a matrix arrives wider. The amplitude a gate produces is only as precise as the register holding it, so the place to decide is where the state is built. An angle written as a plain number is read in double precision on the way to the matrix, which costs nothing and keeps a double-precision register from being fed a single-precision matrix.

wide = tq.zeros_state(2, dtype=tp.complex128)   # this is what decides
tq.apply_ry(wide, 0, 0.3)                        # double-precision amplitudes
tq.apply_ry(tq.zeros_state(2), 0, 0.3)          # single, whatever the angle says

One qubit as a point in the ball

A single qubit is a point of the unit ball, and reading a simulation off in those coordinates is what makes a comparison readable: the distance between two states is the angle between their vectors, a mixed state is inside rather than on, and a fit lands somewhere you can look at.

vector = tq.sv2bloch(state)          # (x, y, z), real, per batch entry
tq.dens2bloch(density)               # the same reading off a mixed state
tq.bloch2dens(vector)                # back to the matrix it names
tq.bloch2sv(vector)                  # and the state, when there is one

A normalized state is on the surface; how far inside a mixed one sits says how mixed it is, being twice the larger eigenvalue less one. The reading is real whatever the state’s own precision, and the reverse is exact for a mixed vector.