# Shor's algorithm - why doesn't the final collapse of the auxiliary qubits cripple the computation?

Let's look at quantum subroutine of Shor's algorithm (image source):

1. Hadamard gates create superposition of all (exponential number) values for input qubits.
2. Then we perform a classical function on them, which is here: $f(a) = y^a \textrm{ mod } N$, where $N$ is the number we would like to factorize, $y$ is some chosen (usually random) number smaller than $N$.
3. Then we perform measurement of value of this function $f(a)=m$ (random). This measurement restricts the original ensemble to only input values $a$, such that $f(a)=m$.
4. Mathematics says that this restricted ensemble has to be periodic, this period can be concluded from value of (quantum) Fourier transform, and allows to conclude the factors.

However, quantum computers require reversible/unitary operations, e.g. we cannot just use OR gate: $(x,y) \to x\ \mathrm{or}\ y$, and instead we need to use e.g. $(x, y, z) \to (x, y, z\ \mathrm{xor}\ (x\ \mathrm{or}\ y))$, which is own inverse $-$ but which requires one additional auxiliary qubit initialized as $z=0$. So their required number is comparable with the number of gates of the classical function - can be quite large.

But what's happening with all these auxiliary qubits at the end of computation, and after the computation is over?

Measuring the classical function has led to the crucial restriction of the original ensemble - can we ensure that auxiliary qubits aren't finally also measured/collapsed, also restricting the ensemble?

Is there some time interval when such measurement restricts the ensemble? (Analogously: required time order between QFT and measurement of classical function?)

If not, can we ensure that restriction from (inevitable?) collapse of auxiliary qubits does not cripple our computation?

Peter Shor has confirmed (below) the problem with auxiliary qubits, requiring to "uncompute" them to fixed values for proper computation process.

In the factoring algorithm, there are three kinds of qubits. In the OP's notation, there are "input qubits", which start in a superposition of all possible values, and which you eventually take the Fourier transform of. There are "value qubits", in which you compute the function $y^a \pmod{N}$, where $a$ is the value in the input qubits. And there are "auxiliary qubits", which you use as workspace to help do this computation.

In order to make the factoring algorithm work properly, you need to reset all the auxiliary qubits, which started as $|0\rangle$ at the beginning of the computation, to $|0\rangle$ at the end of the computation. This is called "uncomputing" these qubits. (Actually, you can set them to anything you please as long as it is a constant independent of the workings of the algorithm.) Theorems about reversible classical computation ensure that it is possible to do this.

If you reset the auxiliary qubits to $|0\rangle$, then if the environment, or somebody, measures them, nothing is revealed about the computation, and the computation is not "crippled". If you forget to reset them to $|0\rangle$, you probably won't get the right answer, whether or not anybody measures them.

• I think that when the asker writes "auxiliary qubit" that they are referring to the qubits in the second register (the one storing $B^k \text{ mod } R$). Whereas your answer (and the diagram in included in the question) are referring to the extra workspace qubits used when computing the second register when they say "auxiliary". Could you make sure to clarify this distinction, since reseting the second register is not necessary and uncomputing it would be very bad? – Craig Gidney Nov 25 '17 at 19:42
• @Craig: Read the original post carefully. He's referring to the extra workspace. "e.g. we cannot just use OR gate: (x,y)→x or y(x,y)→x or y, and instead we need to use [the Toffoli gate], which is own inverse — but which requires one additional auxiliary qubit." – Peter Shor Nov 25 '17 at 19:45
• Thanks for clarifying the distinction. I agree that the question appears clear, but based on commenting/chatting with the OP I think they were confused about needing to measure the value qubits (e.g. in chat they said "I would say that not measuring the [second register], we should have the original ensemble on input: of all 2^n values, hence period returned by QFT would be 1.") – Craig Gidney Nov 25 '17 at 20:25
• Measuring the qubits in the "orthogonal basis" isn't going to work. You need to "uncompute" them. This is discussed briefly on page 9 of my paper, and in much more detail in C. H. Bennett (1973), Logical reversibility of computation. – Peter Shor Nov 25 '17 at 23:20
• +1 Great answer. By the way, I wanted to inform you about the new Stack Exchange site on Quantum Computing, which is in private beta at the moment, but will soon go public. In case you're interested, consider joining. We need more expert users like you on that site. :) – user139621 Mar 15 '18 at 17:02

Update: I originally thought the question was referring to the "value" qubits when the asker said "auxilliary". This answer explains why you don't need to measure the value qubits. For the actual auxilliary qubits, which are used as workspace while computing the value qubits, it should also be okay to measure them afterwards but only because a proper circuit uncomputes them back to 0.

After the value qubits are computed (the ones storing $B^k \text{ mod } R$, not the ones you used as helpers while computing that value!), you can just throw them out. You don't need to measure them or protect them or care for them. Just drop them on the floor. Nothing anyone does to them can hurt the remaining part of the computation. See this walkthrough of Shor's algorithm.

Let's do a simple example via my simulator Quirk. We will initialize a uniform superposition of qubits, and then compute their parity onto an auxiliary qubit (click the image to manipulate the circuit in the simulator):

The two green boxes are showing a representation of the density matrix of the top three qubits. We can show that information without disturbing the system because this is a simulator.

Before the parity computation, the qubits are entirely coherent. Afterwards, some of the off-diagonal indicators have disappeared (become zero). This indicates a partial loss of coherence. The states with an even number of ones have been decohered from the states that have an odd number of ones.

Now let's try to use the auxiliary qubit to "mess with" the top three qubits. If we succeed, the density matrix display will show something different. First thing to try is measurement:

Nothing's different.

Maybe we measured along the wrong axis? Let's rotate the qubit before measuring it:

Still no change!

In fact, no matter what we do to the bottom qubit, we can't change the density matrix of the top three qubits. Not without some kind of operation that crosses between them, or some kind of conditioning (e.g. only consider the subset of states where the measurement of the bottom qubit returned a particular result).

If you find this hard to believe, I recommend simply messing around in Quirk for a while, trying to make the densities of the top three qubits change by operating only on the bottom qubit.

Another way to confirm that it doesn't matter if you measure the auxilliary qubits is to just do the algebra and check.

The initial state is:

$$|\psi_0\rangle = |0\rangle_{\text{main}} \otimes |0\rangle_{\text{aux}} = |0\rangle_{\text{all}}$$

Then we Hadamard transform the main register:

$$|\psi_1\rangle = H_{\text{main}} |\psi_0\rangle = \sum_{k=0}^{2^n-1} |k\rangle_{\text{main}} \otimes |0\rangle_{\text{aux}}$$

Note that I am ignoring normalizing factors. In the end my argument is going to be based on the proportional size of various cases, instead of absolute size, so this is fine.

Then we pick a random base $B$ apply the modular exponentiation operation, which adds B^k mod R into the auxiliary register where k is the computational basis value of the main register. On an actual machine we'd use some temporary workspace to implement this operation, but that all gets cleaned up so here we only care about the effect on the main and aux registers:

$$M = \Big[ \text{aux} \text{ += } B^{\text{main}} \text{ mod } R \Big]$$

$$|\psi_2\rangle = M \cdot |\psi_1\rangle = \sum_{k=0}^{2^n-1} |k\rangle_{\text{main}} \otimes |B^k \text{ mod } R\rangle_{\text{aux}}$$

Now we can rewrite $k$ to be in terms of the unknown period of $B^k \text{ mod } R$. We'll use $k = l \cdot m + s$ where $l$ is the period, $s$ is an iteration variable for the offset between 0 and $l$, and $m$ is an iteration variable. With that in mind, we rewrite $|\psi_2\rangle$ as:

$$|\psi_2\rangle = \sum_{m=0}^{\;\;\lceil 2^n / l \rceil-1\;\;} \sum_{s=0}^{\text{min}(l, 2^n-lm)-1} |lm+s\rangle_{\text{main}} \otimes |B^{lm+s} \text{ mod } R\rangle_{\text{aux}}$$

Note that $B^{lm+s} = B^{s} \pmod{R}$. Also note that the complicated boundary conditions on $s$ can be simplified by approximating our actual sum with a sum that goes up to the first multiple of $l$ after $2^n$. This is a good approximation as long as $2^n >> l$, which is true since $n$ is chosen such that $2^n > R^2$ and we know that $R > l$. Anyways, after applying that simplification and approximation we get:

$$|\psi_2\rangle \approx \sum_{m=0}^{\;\;\lceil 2^n / l \rceil-1\;\;} \sum_{s=0}^{l-1} |lm+s\rangle_{\text{main}} \otimes |B^{s} \text{ mod } R\rangle_{\text{aux}}$$

Because the boundary condition of $s$ doesn't depend on $m$ anymore, we can rearrange the order of the sums. Which gives us:

$$|\psi_2\rangle \approx \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} |lm+s\rangle_{\text{main}} \right)$$

Now we apply the inverse Fourier transform operation to the main register. Notice that it can be moved from the outside of the sum to the inside:

\begin{align} |\psi_3\rangle &= \text{QFT}^{\dagger}_{\text{main}} \cdot |\psi_2\rangle \\ &\approx \text{QFT}^{\dagger}_{\text{main}} \cdot \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} |lm+s\rangle_{\text{main}} \right) \\ &= \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{QFT}^{\dagger}_{\text{main}} \cdot |lm+s\rangle_{\text{main}} \right) \end{align}

Then expand the definition of the QFT into a sum over a variable $j$, and move that sum outward:

\begin{align} |\psi_3\rangle &\approx \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{m=0}^{\lceil 2^n / l \rceil-1} \;\; \sum_{j=0}^{2^n-1} |j\rangle_{\text{main}} \cdot \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right) \\ &= \sum_{s=0}^{l-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes \sum_{j=0}^{2^n-1} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right) \\ &= \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} \left( |B^{s} \text{ mod } R\rangle_{\text{aux}} \otimes |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right) \end{align}

Now we're going to measure the main register. The probability of getting the result $r$ is the total squared magnitude of states where the first register is $r$. Algebraically:

\begin{align} P(r) &= \sum_{a,b | a=r} \Big| (\langle a |_{\text{main}} \otimes \langle b |_{\text{aux}}) \cdot | \psi_3 \rangle \Big|^2 \\ &= \sum_{b} \Big| (\langle r |_{\text{main}} \otimes \langle b |_{\text{aux}}) \cdot | \psi_3 \rangle \Big|^2 \\ &\approx \sum_{b} \left| \langle r |_{\text{main}} \langle b |_{\text{aux}} \cdot \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} |B^{s} \text{ mod } R\rangle_{\text{aux}} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \end{align}

Because all of our basis kets are perpendicular, any summands that fail to satisfy $b=B^s \pmod{R}$ and $r=lm+s$ will be zero'd out. The remaining terms have bras and kets that exactly match, giving an inner product of 1. I'll do this over a few steps because it simplifies the sum considerably:

\begin{align} P(r) &\approx \sum_{b} \left| \langle r |_{\text{main}} \langle b |_{\text{aux}} \cdot \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} |B^{s} \text{ mod } R\rangle_{\text{aux}} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{b} \left| \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} \langle r |_{\text{main}} \langle b |_{\text{aux}} \cdot |B^{s} \text{ mod } R\rangle_{\text{aux}} |j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \text{exp}(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{b} \left| \sum_{s=0}^{l-1} \sum_{j=0}^{2^n-1} \langle r | j\rangle_{\text{main}} \langle b | B^{s} \text{ mod } R\rangle_{\text{aux}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{j=0}^{2^n-1} \langle r | j\rangle_{\text{main}} \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot j \cdot (lm+s)) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot (lm+s)) \right|^2 \end{align}

Now we're getting somewhere. The next thing to do is get rid of that pesky $s$. Factor the $s$ component out of the inner sum, which lets you factor it out of the squared-magnitude, at which point you realize it contributes nothing and the sum can turn into a multiplication by $l$:

\begin{align} P(r) &\approx \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot (lm+s)) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \cdot \exp(i\tau \cdot 2^{-n} \cdot r \cdot s) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \exp(i\tau \cdot 2^{-n} \cdot r \cdot s) \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \right|^2 \\ &= \sum_{s=0}^{l-1} \big| \exp(i\tau \cdot 2^{-n} \cdot r \cdot s) \big|^2 \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \right|^2 \\ &= \sum_{s=0}^{l-1} \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot r \cdot lm) \right|^2 \\ &= l \cdot \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot rl \cdot m) \right|^2 \end{align}

Nearly there. To make the structure of the sum blatant, we extract a variable $\omega = \exp(i\tau rl / 2^{n})$:

\begin{align} P(r) &\approx l \cdot \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \exp(i\tau \cdot 2^{-n} \cdot rl \cdot m) \right|^2 \\ &= l \cdot \left| \sum_{m=0}^{\lceil 2^n / l \rceil-1} \omega^{m} \right|^2 \text{ where } \omega = \exp(i\tau r l / 2^{n}) \end{align}

The inner sum will be largest when all of its terms point in the same direction, i.e. when $\omega \approx 1$. Which means $\exp(i\tau rl / 2^{n}) \approx 1$, which in turn means that $2^{-n} r l$ is nearly an integer $d$. Rewrite $2^{-n} r l \approx d$ and you get:

$$r \approx 2^n \cdot \frac{d}{l}$$

In other words, if the period is $l$ then the values you are most likely to measure are placed near multiples of $2^n / l$. In practice you recover $l$ by solving "which possible multiple was my measurement nearest to?".

I leave it as an exercise for the reader to work out exactly how much more likely you are to measure values of $r$ that give values of $\omega$ close to 1.

• A small self contained example might help. – DanielSank Nov 19 '17 at 5:19
• Sure there is no FTL (or back in time) communication - like in EPR or Delayed Choice QE, there is no way to actually send information this way. However, Shor uses the fact that measurement of classical function affects the computation (by restricting ensemble, delay seems unimportant?), so analogous effect should be true for auxiliary qbits (?), as physically situation is practically the same. Just releasing them into environment should collapse them - what seems analogous to measurement (?) Maybe we should e.g. shoot such photons into space, such that they could fly uncollapsed infinitely? – Jarek Duda Nov 19 '17 at 6:54
• @DanielSank Done. – Craig Gidney Nov 19 '17 at 19:48
• @JarekDuda It's fine if they collapse. All of the decoherence damage they can do was already done by computing them in the first place. Collapsing them won't do any more damage to the state of the remaining qubits. I have rewritten my answer. – Craig Gidney Nov 19 '17 at 19:49
• @CraigGidney, so why measuring/collapsing the qbits of classical function influences the original ensemble, but it is not true for auxiliary qbits? From physical point of view they seem nearly the same: used in gates on the way - what is the difference between them? Maybe collapsing these auxiliary qbits is one of the reasons of the problem with decoherence ? Can we ensure that qbit will never collapse, e.g. by shooting such photon into space? – Jarek Duda Nov 20 '17 at 7:46