Thoughts on Control Flow

[moodlezoup]

Perhaps the biggest open question in our project thus far is how to deal with classical control flow. When branches are introduced to a program via e.g. if_then statements in pyQuil, we would like to (1) characterize the probability of each branch and (2) accurately estimate the density matrix at an arbitrary breakpoint in one of the branches.

Our current paradigm centers around the following picture, describing the most basic case of classical control in a Quil program:

      A
      |   # measure(B)
     / \
    /   \
Q> X     Y

A, X, and Y in this picture correspond to "basic blocks", i.e. instances of the QuilBlock class in control_flow_graph.py. The measure(B) indicates that the branch is conditioned on a classical register which depends on measurements of qubits in some subset B of the entire system. We imagine a breakpoint in X, at which the user queries for the state of some other qubit subset Q.

What we would like to compute, to draw a Bayesian analogy, is the density function of Q after running A + X, conditioned on the measure(B) satisfying the branch condition. In (informal) mathematical notation P(state(Q) | A+X, measure(B) in M_x), where M_x denotes the subset of {0, 1}^|B| which would lead us down the X branch.

There are three cases of B and Q that we need to consider: two of which have simple solutions, and the third of which has required a lot more thought. In what follows, we will use the notation ent(B) to denote the entanglement set of B induced by the multi-qubit gates in A + X (see the function entanglement_set(B)). Similarly, let ent_A(B) denote the entanglement set of B induced by the multi-qubit gates in A alone.

1. ent(Q) and ent(B) are disjoint. This is the simplest case. We can characterize the branch probabilities by introducing an ancilla bit which gets flipped in one branch but not the other. We can then run A many times and see how many times the ancilla bit gets flipped. We can then trim A+X to only include gates involving ent(Q) and run state tomography to estimate the desired density matrix.
2. ent_X(Q) and (ent_A(B) - B) are disjoint, ent_A(Q) and ent_A(B) are disjoint. In other words, Q and B are not entangled at any point in A, but after B is measured we know what state B will be in going into. In this scenario, we can run A many times to determine the probability of each possible measurement in M_x, i.e. post-measurement states of B that satisfy the branch condition. Then we can "clamp" B to one such post-measurement state, trim gates involving ent_A(B) from A, then run state tomography on the resulting A'+X. Repeat for each state in M_x.
3. The general case, i.e. no assumptions on the entanglement of Q and B. This presents a problem: on one hand, we don't know how measuring B affects the state of Q; on the other, we need to be in the right branch, or at least know which branch a particular execution branch would have gone down.

The trickiness of the general case led us to explore several options. I have compiled here a list of four such options:

1. Enforce case 2: In other words, don't handle the general case at all. Easy, but not exactly desirable.
2. Branch rewiring: (related) In the above example, we'd replace Y with RESET; JUMP @A. This would effectively force the execution down the X branch through trial and error. This a nice and clean solution, but has two problems: (a) forest-benchmarking's state tomography doesn't support classical control, so we'd have to either make a PR or rely on an update on their end; (b) this is computationally expensive if we care about a low-probability branch.
3. Process tomography: Run state tomography on A, and process tomography on X. Then by combining the two results, we arrive at an estimate of the density matrix for A+X. This is probably feasible, but process tomography is computationally expensive.
4. Make a guess and run with it: Run state tomography on A and compute a particular ensemble of wavefunctions from the resulting density matrix, e.g. via eigendecomposition. For each wavefunction psi_i in this ensemble, let X_i = [initialize qubits to psi_i] + X. Run state tomography on each X_i. The problem with this option is that the wavefunction ensemble can only be correct up to a unitary transform, which could completely change the state of Q.
5. QVM hybrid: We observe that if we had an oracle for the wavefunction of ent_A(B), we'd basically reduce the general case to the second case. So, we'll let the QVM play the role of this oracle. Trim A to only include gates involving ent_A(B), and run the QVM on this trimmed program to obtain the wavefunction ensemble for ent_A(B) at the branching point. Then we proceed as we did in case 2: clamp to one of the wavefunctions in the ensemble, and trim gates involving ent_A(B) from A, and run state tomography for Q on A'+X. The only problem we anticipate from this approach is if |ent_A(B)| is too large to tractably simulate on the QVM. But we hope that this will be a good compromise for most practical applications.

We currently plan to proceed with the QVM hybrid as our main approach. We can also fall back to the simpler approaches for cases 1 and 2 if we determine we are in one of those two scenarios.

[moodlezoup]

Open question: how do we generalize to arbitrary DAGs?