What is the "sum-over-paths" picture of quantum computing?

Question

In chapter 9 of Scott Aaronson's book Quantum Computing Since Democritus (see online progenitor lecture notes here), he introduces a strange (to me) conceptualization of applying two Hadamard-like operations to a qbit, summed up in this diagram:

I had the opportunity to ask him about this and he made three interesting statements:

1. This is called the "Feynman" or "sum-over-paths" picture of quantum mechanics
2. This picture is interesting because it reveals the effect of destructive interference
3. If a QC is simulated using the sum-over-paths method, it takes exponential time instead of space

How does the sum-over-paths picture work conceptually for quantum computation? For example, given a series of quantum logic gates applied to some qbits, how would I calculate the resulting amplitude values of those qbits using the sum-over-paths method? If this is too complicated to explain in an answer on this site, are there any good introductions to this topic elsewhere?

Also, do any of the main quantum languages/frameworks include a sum-over-paths simulator?

Answer 1

Initially perhaps it is best to compare and contrast a couple of different approaches to simulating the output of a random quantum circuit $U$ of about $m=1000$ gates acting on, say, $n=53$ qubits initialized to $\vert 0^{53}\rangle$, as in Sycamore. Suppose we execute our circuit $U$ on the $53$ qubits of Sycamore, and we sample to get an output string $b\in\{0,1\}^{53}$. We wish to determine the probability of getting an output $b$, e.g. determine $\Vert\langle b\vert U\vert 0^{53}\rangle\Vert^2$

1. In a first "Schrödinger" approach, we can store all $2^{53}$ amplitudes in memory, and update each amplitude based on the $m=1000$ gates, one for each step. This requires memory on the order of $2^n$ (e.g. a lookup table for each amplitude), but time on the order of $m 2^n$ (just keep updating the wavefunction for each gate), followed by a quick look-up of the amplitude/probability for $b$.
2. In a second "Feynman" approach, one instead puts primacy on the actual output $b$ that is sampled, and performs the sum-over-paths calculations to determine the corresponding probability of obtaining $b$. Because only the measured output string $b$ is stored, memory usage is much improved to about $m+n$, but time takes roughly $4^m$ (doubling for each branch).
3. One can combine the two approaches as considered by Aaronson and Chen - into what they termed the "Schrödinger-Feynman" hybrid algorithm, that uses $\mathrm{poly}(m,n)$ space and $m^{O(n)}$ time.

Given a decision problem, we can put the output of the decision into a register $b\in\{0,1\}$; say $b=0$ if the decision is NO but $b=1$ if the decision is YES. If we have a quantum circuit to determine $b$, we can see that we can run a Schrödinger algorithm or a Feynman algorithm to determine whether $b=1$.

There are some implications for complexity classes of Feynman's path integral approach that I am not versed in, but I think Aaronson is fond of saying something like "Feynman won his Nobel by proving $\mathrm{BQP}\subseteq\mathrm{P}^\mathrm{\#P}$", well before $\mathrm{BQP}$ was even defined.

Answer 2

Expanding on Mark S's answer with some math:

Let's take Aaronson's example, where you are applying the following unitary $U$ twice to a qbit initialized to $|0\rangle$:

$$ U = \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{-1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \end{bmatrix} $$

The conventional (Schrödinger) picture would represent this computation as follows:

$$ UU|0\rangle = \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{-1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \end{bmatrix} \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{-1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \end{bmatrix} \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} \frac{1}{\sqrt{2}} & \frac{-1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\ \end{bmatrix} \begin{bmatrix} \frac{1}{\sqrt{2}} \\ \frac{1}{\sqrt{2}} \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \end{bmatrix} = |1\rangle $$

Suppose you didn't care about the value of the entire state vector, but rather just the value of a single amplitude. Could you calculate this amplitude without keeping the entire $2^n$-sized state vector around in memory? Turns out yes, you can, with the the Feynman (or sum-over-paths) picture. Consider the following equation:

$$ |x\rangle = U_m \ldots U_1|\psi\rangle $$

How do we calculate the value of the final amplitude $i$, $|x\rangle[i]$? It's a simple recursive algorithm. Let's define the recursive function $f(|\psi\rangle, U, i)$, which returns the value of amplitude $i$ after the (1-qbit (2x2) or 2-qbit (4x4)) unitary operations in $U$ are applied to initial state $|\psi\rangle$:

1. Base case: $U = [\ ]$, return $|\psi\rangle[i]$
2. General case: $U = [U_m, \ldots, U_{1}]$, $U_m$ is a 1-qbit 2x2 matrix: in the $i$th row of the full matrix $M = \mathbb{I}_2 \otimes \ldots \otimes U_m \otimes \ldots \otimes \mathbb{I_2}$ (which would be multiplied against the state vector) there will be at most two nonzero entries at some indices $a$ and $b$; this means we need the amplitudes at index $a$ and $b$ for the state prior to $U_m$ being applied, so return: $$ M[i,a]\cdot f(|\psi\rangle, Tail(U), a) + M[i,b] \cdot f(|\psi\rangle, Tail(U), b) $$
3. General case: $U = [U_m, \ldots, U_{1}]$, $U_m$ is a 2-qbit 4x4 matrix: in the $i$th row of the full matrix $M = \mathbb{I}_2 \otimes \ldots \otimes U_m \otimes \ldots \otimes \mathbb{I_2}$ (which would be multiplied against the state vector) there will be at most four nonzero entries at some indices $a$, $b$, $c$, and $d$; this means we need the amplitudes at index $a$, $b$, $c$, and $d$ for the state prior to $U_m$ being applied, so return: $$ M[i,a]\cdot f(|\psi\rangle, Tail(U), a) + M[i,b] \cdot f(|\psi\rangle, Tail(U), b) + M[i,c] \cdot f(|\psi\rangle, Tail(U), c) + M[i,d] \cdot f(|\psi\rangle, Tail(U), d) $$

For our example above, this looks as follows; calculating the value of $f(|0\rangle, [U, U], 2)$ (aka the final amplitude for state $|1\rangle$):

$$ f(|0\rangle, [U, U], 2) = U[2,1] f(|0\rangle, [U], 1) + U[2,2] f(|0\rangle, [U], 2) $$ $$ f(|0\rangle, [U, U], 2) = \frac{1}{\sqrt{2}} f(|0\rangle, [U], 1) + \frac{1}{\sqrt{2}} f(|0\rangle, [U], 2) $$ $$ f(|0\rangle, [U, U], 2) = \frac{1}{\sqrt{2}} \left( \frac{1}{\sqrt{2}} f(|0\rangle, [\ ], 1) + \frac{-1}{\sqrt{2}} f(|0\rangle, [\ ], 2) \right) + \frac{1}{\sqrt{2}} f(|0\rangle, [U], 2) $$ $$ f(|0\rangle, [U, U], 2) = \frac{1}{\sqrt{2}} \left( \frac{1}{\sqrt{2}} \cdot 1 + \frac{-1}{\sqrt{2}} \cdot 0 \right) + \frac{1}{\sqrt{2}} f(|0\rangle, [U], 2) $$ $$ f(|0\rangle, [U, U], 2) = \frac{1}{2}+ \frac{1}{\sqrt{2}} \left( \frac{1}{\sqrt{2}} f(|0\rangle, [\ ], 1) + \frac{1}{\sqrt{2}} f(|0\rangle, [\ ], 2) \right) $$ $$ f(|0\rangle, [U, U], 2) = \frac{1}{2}+ \frac{1}{\sqrt{2}} \left( \frac{1}{\sqrt{2}} \cdot 1 + \frac{1}{\sqrt{2}} \cdot 0 \right) $$ $$ f(|0\rangle, [U, U], 2) = 1 $$

Given $m$ unitary operations on $n$ qbits, this will clearly take $O(4^m)$ time if every operation is a two-qbit 4x4 matrix, since it will lead to four additional recursive calls at every level. It will take at least $O(m)$ space, since the recursive call stack will grow to length proportional to $m$ (other sources say $O(m + n)$ but I'm not sure where the $n$ term comes in).

It's somewhat unclear how you could efficiently determine the indices of amplitudes affected by the operator at each recursive level, but I'm sure if I thought about it hard enough there'd be some obvious trick. I'll also have to think some more about how this recursive algorithm relates to the branching diagram up above.