Quantum key distribution using BB84

Suppose Alice has a secret message she wants to send to Bob. They have two information channels: a quantum channel (such as a stream of photons) and an unsecure classical channel (such as unencrypted email). Let’s assume that the message and key are both encoded using binary and the key is a one-time pad. They will eventually XOR message and key to encrypt and decrypt the message.

BB84 (Bennett & Brassard, 1984) is an approach to quantum key distribution that makes it possible for Alice to send an encryption key to Bob and for them to verify with a high level of confidence that the key hasn’t been intercepted. Once convinced that the key was sent securely, Alice can use it to encrypt the secret message and send it over the unsecure classical channel. Bob safely decrypts at the other side.

BB84 takes advantages of three related features of quantum mechanics.

  1. If a qubit is in superposition, e.g., \(\frac{|0\rangle + |1\rangle}{\sqrt{2}}\), then it’s impossible to infer its state from one observation. In particular, suppose an eavesdropper, Eve, measured the qubit and observed \(|1\rangle\). She couldn’t tell whether it was in the state \(\frac{|0\rangle + |1\rangle}{\sqrt{2}}\) or \(|1\rangle\).
  2. Once you measure a qubit in superposition, it collapses to a basis state, so Eve couldn’t intercept it and forward it onto Bob.
  3. You also can’t clone a qubit that is in an unknown state so, e.g., Eve couldn’t intercept a qubit from Alice, clone it, measure one of the pair and pass the other unmeasured qubit onto Bob.

Here’s how BB84 works.

Alice begins by generating a stream of random bits:

110110100110000101110…

Some of these bits will be used for the key.

Next, she generates a random stream of bases, classical basis (\(+\)) or the Hadamard basis (\(\times\)):

\(+\)\(\times\)\(+\)\(+\)\(\times\)\(+\)\(\times\)\(+\)\(\times\)\(\times\)\(+\)\(\times\)\(\times\)\(+\)\(\times\)\(+\)\(\times\)\(\times\)\(+\)\(\times\)\(+\)…

These will determine which basis will be used to send each bit of the key. The figure below shows how Alice will encode the 1s and 0s on these bases (Mayers, 2001).

For a photon, the classical basis corresponds to vertical/horizontal polarisation and the Hadamard basis to the diagonal polarisations. Note in this case the state space corresponds nicely with the physical realisation of the qubit; that won’t always be the case, e.g., consider spin up and down states for electrons.

Here are the qubit states Alice sends:

\(\Psi(0, +) = |0\rangle\)

\(\Psi(1, +) = |1\rangle\)

\(\displaystyle \Psi(0, \times) = H|0\rangle = \frac{|0\rangle + |1\rangle}{\sqrt{2}}\)

\(\displaystyle \Psi(1, \times) = H|1\rangle = \frac{|0\rangle -|1\rangle}{\sqrt{2}}\)

Where \(\displaystyle H = \frac{1}{\sqrt2} \begin{bmatrix}1 & 1\\1 & -1\end{bmatrix}\).

Alice sends the potential key bits using the randomly chosen bases.

Bob measures them, randomly deciding whether to assume a classical or Hadamard basis – for the latter applying the Hadamard operator before measuring. The possible outcomes are as follows:

Alice sends Bob assumes classical Bob assumes Hadamard
\(\Psi(0, +) = |0\rangle\) \(P(|0\rangle) = 1\) \(P(|0\rangle) = \frac{1}{2}\)
\(\Psi(1, +) = |1\rangle\) \(P(|1\rangle) = 1\) \(P(|1\rangle) = \frac{1}{2}\)
\(\Psi(0, \times) = H|0\rangle\) \(P(|0\rangle) = \frac{1}{2}\) \(P(|0\rangle) = 1\)
\(\Psi(1, \times) = H|1\rangle\) \(P(|1\rangle) = \frac{1}{2}\) \(P(|1\rangle) = 1\)

So if Bob randomly chooses the correct basis, he gets the correct bit; if not, he gets a random bit, 50-50 1 or 0.

Bob tells Alice what bases he used on the unsecure channel. At this point, doing so doesn’t help any eavesdroppers since the qubits have already been measured.

Alice tells Bob which bases were correct, again on the unsecure channel, so Bob can discard measurements for the rest. She also sacrifices a random sample of the potential key bits by sharing what the right answer is for those. If Eve the eavesdropper intercepted any, then they won’t all have made it to Bob.

If Alice and Bob are satisfied that the test qubits made it across okay, then they use the remainder for the key, sending the encrypted message through the unsecure channel.

That’s the gist. The devil’s in the detail, e.g., working out how many bits to send to ensure enough remain after checking for intercepts; how many to sacrifice to check that Eve wasn’t eavesdropping; and taking account of errors in transmission.

The security of BB84 depends on the quality of the implementation. The field of quantum hacking explores a variety of ways to exploit imperfections, e.g., classical side channels that leak information about what was sent, which can be intercepted without interfering with the qubit. Researchers devise clever ways to defend against these attacks. See, e.g., Dixon et al. (2017).

References

Bennett, C. H. & Brassard, G. (1984). Quantum cryptography: Public key distribution and coin tossing. Proceedings of IEEE International Conference on Computers, Systems, and Signal Processing (pp. 175-179), India.

Dixon, A. R., Dynes, J. F., Lucamarini, M., Fröhlich, B., Sharpe, A. W., Plews, A., Tam, W., Yuan, Z. L., Tanizawa, Y., Sato, H., Kawamura, S., Fujiwara, M., Sasaki, M., & Shields, A. J. (2017). Quantum key distribution with hacking countermeasures and long term field trial. Scientific Reports, 7, 1978.

Mayers, D. (2001). Unconditional security in quantum cryptography. Journal of the ACM, 48, 351–406.

Quantum entanglement

This brief post shows how to establish an entangled state on IBM Quantum’s computers – here using ibmq_manila.

The circuit to establish \(\displaystyle \frac{|00\rangle + |11\rangle}{\sqrt{2}}\) is as follows:

A few clicks later, and the circuit is added to the queue…

The total wait to get onto ibmq_manila was about an hour and the circuit ran 1,000 times in 5 seconds.

Here are the results:

It worked across 93.9% of runs: both qubits had the same measurement outcome, roughly 50-50 split between \(|00\rangle\) and \(|11\rangle\).

 

Why do quantum circuits often end with an exclusive-OR (XOR)?

In a previous post I worked through the arithmetic of the Deutsch quantum algorithm to show that it works, but I avoided attempting to explain why any of it works. In this post I’ll explain one corner of the circuit: why the gate implementing the function tested by the algorithm outputs the exclusive-OR (XOR) of \(y\) and \(f(x)\): \(y \oplus f(x)\). See the picture below (Mermin, 2007, p. 42). The question is, why is that output not just \(f(x)\), since that’s what the circuit is supposed to compute?

There are two parts to the explanation. Firstly, quantum gates have to be reversible, because that’s how the physics works. The reason for that is far bigger than this blog post and I’m not a physicist. Try Richard Feynman’s (1985) explanation. If we believe Feynman, then the reason XOR appears in many algorithms is easy to grasp. Below I’ll expand the arithmetic provided by Mermin (2007, p. 37, the paragraph around equation 2.3).

Let’s start with the punchline.

Double application of \(U_f\) leads to the bottom register of the circuit evaluating to \([y \oplus f(x)] \oplus f(x)\), which is equivalent to \(y \oplus [f(x) \oplus f(x)]\) by the associativity of \(\oplus\). This simplifies to \(y\), as the table below shows for the four combinations of \(y\) and \(f(x)\):

\(y\) \(f(x)\) \(\overbrace{f(x) \oplus f(x)}^{\textstyle z}\) \(y \oplus z\)
0 0 0 0
0 1 0 0
1 0 0 1
1 1 0 1

Note how \(f(x) \oplus f(x)\) cancels out \(f(x)\), leaving \(y\), so the first and last column are equal. The XOR ensures that double application of \(U_f\) gives us \(y\) again!

Another way to write this without the table:

\([y \oplus f(x)] \oplus f(x)\)

= { associativity of \(\oplus\) }

\(y \oplus [f(x) \oplus f(x)]\)

= { since \(1 \oplus 1 = 0\oplus 0 = 0\) }

\(y \oplus 0\)

= { definition of \(\oplus\) }

\(y\).

Now we just have to show that

\(U_f U_f (|x\rangle \otimes |y\rangle)\)

sends us to

\(| x\rangle \otimes | y \oplus f(x) \oplus f(x)\rangle\)

so that the sums above are relevant:

\(U_f U_f (|x\rangle \otimes |y\rangle)\)

= { by application of \(U_f\) }

\(U_f (|x\rangle \otimes |y \oplus f(x)\rangle)\)

= { by application of \(U_f\) again }

\(|x\rangle \otimes |[y \oplus f(x)] \oplus f(x)\rangle\).

We’re done.

More on reversibility

The requirement that gates are reversible is often described by saying that quantum algorithms change the state of qubits using unitary transformations. A matrix \(U\) representing such a transformation is unitary if and only if \(U U^\dagger = I\), where \(I\) is the identity matrix and \(U^\dagger\) is the conjugate transpose of \(U\) (see Adedoyin et al., 2022). If the elements of \(U\) don’t have any imaginary parts then \(U^\dagger\) is just the transpose of \(U\).

Let’s try it for the most complicated (relatively speaking!) function from the Deutsch problem, \(U_{f_2}\):

As we saw, this is represented by the following matrix:

\(U_{f_2} = \mbox{CNOT} (I \otimes X) = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}\)

Is it unitary?

\(U_{f_2} U_{f_2}^\dagger\)

=

\(\begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}\)

=

\(\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}\)

=

\(I\).

Yes.

References

Feynman, R. P. (1985). Quantum mechanical computers. Optics News, 11(2), 11–20.

Mermin, N. D. (2007). Quantum computer science: an introduction. Cambridge University Press.

J., A., Adedoyin, A., Ambrosiano, J., Anisimov, P., Casper, W., Chennupati, G., Coffrin, C., Djidjev, H., Gunter, D., Karra, S., Lemons, N., Lin, S., Malyzhenkov, A., Mascarenas, D., Mniszewski, S., Nadiga, B., O’Malley, D., Oyen, D., Pakin, S., … Lokhov, A. Y. (2022). Quantum Algorithm Implementations for Beginners. ACM Transactions on Quantum Computing, 3(4), 1–92.

IBM Quantum ibmq_quito runs of the Deutsch algorithm

 

The Deutsch algorithm decides whether a function \(f : \{0, 1\} \rightarrow \{0, 1\}\) is constant (\(f(0) = f(1)\)) or balanced (for this problem, \(f(0) \ne f(1)\), but see also the Deutsch–Jozsa algorithm).

There are four possible functions, \(f_i\), as follows:

Function \(f_i(0)\) \(f_i(1)\) Type
\(f_0\) 0 0 Constant
\(f_1\) 0 1 Balanced
\(f_2\) 1 0 Balanced
\(f_3\) 1 1 Constant

Classically, we would need to query \(f\) twice to decide whether an unknown function is constant or balanced. The quantum algorithm uses superposition to get the answer in one query.

I worked through the sums and circuit notation in this post. Below, circuits for running the algorithm for each of the four functions on one of IBM Quantum‘s quantum computers, ibmq_quito. I ran each 1,000 times to see how often an actual quantum computer(!) gets the right answer. I’ve also pasted in the Python code, which was generated automatically, and can be used to setup the circuits again.

The measurement outcome should be 1 for a constant function and 0 for a balanced function. As you can see, the modal answer across runs is correct, but there is a bit of noise.

f0 (constant)

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from numpy import pi
qreg_q = QuantumRegister(2, 'q')
creg_c = ClassicalRegister(1, 'c')
circuit = QuantumCircuit(qreg_q, creg_c)
circuit.reset(qreg_q[0])
circuit.reset(qreg_q[1])
circuit.x(qreg_q[0])
circuit.x(qreg_q[1])
circuit.h(qreg_q[0])
circuit.h(qreg_q[1])
circuit.h(qreg_q[0])
circuit.measure(qreg_q[0], creg_c[0])
# @columns [0,0,1,1,2,2,3,4]

f1 (balanced)

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from numpy import pi
qreg_q = QuantumRegister(2, 'q')
creg_c = ClassicalRegister(1, 'c')
circuit = QuantumCircuit(qreg_q, creg_c)
circuit.reset(qreg_q[0])
circuit.reset(qreg_q[1])
circuit.x(qreg_q[0])
circuit.x(qreg_q[1])
circuit.h(qreg_q[0])
circuit.h(qreg_q[1])
circuit.cx(qreg_q[0], qreg_q[1])
circuit.h(qreg_q[0])
circuit.measure(qreg_q[0], creg_c[0])
# @columns [0,0,1,1,2,2,3,4,5]

f2 (balanced)

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from numpy import pi
qreg_q = QuantumRegister(2, 'q')
creg_c = ClassicalRegister(1, 'c')
circuit = QuantumCircuit(qreg_q, creg_c)
circuit.reset(qreg_q[0])
circuit.reset(qreg_q[1])
circuit.x(qreg_q[0])
circuit.x(qreg_q[1])
circuit.h(qreg_q[0])
circuit.h(qreg_q[1])
circuit.x(qreg_q[1])
circuit.cx(qreg_q[0], qreg_q[1])
circuit.h(qreg_q[0])
circuit.measure(qreg_q[0], creg_c[0])
# @columns [0,0,1,1,2,2,3,4,5,6]

f3 (balanced)

from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from numpy import pi
qreg_q = QuantumRegister(2, 'q')
creg_c = ClassicalRegister(1, 'c')
circuit = QuantumCircuit(qreg_q, creg_c)
circuit.reset(qreg_q[0])
circuit.reset(qreg_q[1])
circuit.x(qreg_q[0])
circuit.x(qreg_q[1])
circuit.h(qreg_q[0])
circuit.h(qreg_q[1])
circuit.x(qreg_q[1])
circuit.h(qreg_q[0])
circuit.measure(qreg_q[0], creg_c[0])
# @columns [0,0,1,1,2,2,3,4,5]

First time running a quantum circuit

This eve I clicked a few things, pieced together a simple quantum computing circuit and ran it on a quantum computer – an actual quantum computer(!), IBM Quantum’s ibm_nairobi.

Here’s the circuit, the Deutsch algorithm with \(U_{f_2}\) slotted in as the function to test. Pretend you can’t see what function it is – the result of the algorithm will tell us (something about) which it is.

In hindsight, I could have simplified it…

Wait a bit…

Then finally, a result!

The zeroth bit (least significant bit) is (almost always) zero, so the function is balanced.

 

The arithmetic of the Deutsch algorithm

The Deutsch algorithm (named after its originator David Deutsch) is a vivid and simple illustration of the power of quantum computing. Beyond achieving this admirable goal, it has no practical purpose, though it does seem to have inspired more useful algorithms. This blog post works through the arithmetic to show that the algorithm works (I did the sums to convince myself as I couldn’t grasp what was going on in Dirac notation). I don’t attempt to explain why it works.

The problem

Here’s the problem the Deutsch algorithm solves.

Suppose you are presented with a function \(f : \{0, 1\} \rightarrow \{0, 1\}\). You don’t know how \(f\) is defined, but you can call it with a \(0\) or a \(1\) and see what result it gives. The result will always be a \(0\) or a \(1\).

Your task is to decide whether or not \(f\) is constant, that is, \(f(0) = f(1)\). If \(f\) isn’t constant, it is said to be balanced since for one input the answer is \(0\) and for the other it is \(1\).

There are four possible functions, \(f_i\), as follows:

Function \(f_i(0)\) \(f_i(1)\) Type
\(f_0\) 0 0 Constant
\(f_1\) 0 1 Balanced
\(f_2\) 1 0 Balanced
\(f_3\) 1 1 Constant

Classically, to distinguish between the constant and balanced functions we would simply call \(f(0)\) and \(f(1)\), look at the results, and see which function it is. The Deutsch quantum algorithm can distinguish between a constant and balanced function by calling the function only once. It can’t tell us which function it is, just whether it is constant or balanced.

The circuit

Here is the circuit that solves it (diagram mildly simplified from Alice Flarend and Bob Hilborn, 2022, Fig 11.3, p. 163):

The gate \(U_f\) is a representation of the unknown function, \(H\) is the Hadamard gate (named after Jacques Hadamard), and \(M\) is measurement. Let’s work through the arithmetic. I’ll start with the easy bits at the input and output before tackling \(U_f\) in the middle.

The inputs

To get it going, first we need to apply \(H\) to the basis state \(|1\rangle\).

\(\displaystyle H = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}\),

\(\displaystyle |1\rangle = \begin{bmatrix} 0 \\ 1  \end{bmatrix}\), so

\(\displaystyle H|1\rangle = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1  \end{bmatrix}\).

Next, since the circuit starts with two qubits both in the state \(H|1\rangle\), we need to calculate the tensor product \(H|1\rangle \otimes H|1\rangle\):

\(\displaystyle H|1\rangle \otimes H|1\rangle = \frac{1}{2} \begin{bmatrix} 1 \\ -1 \\ -1 \\ 1  \end{bmatrix}\).

In ket notation, and reading top to the bottom of the column matrix, this is equivalent to:

\(\displaystyle \frac{1}{2}|00\rangle-\frac{1}{2}|01\rangle-\frac{1}{2}|10\rangle+\frac{1}{2}|11\rangle\).

\(|00\rangle\) means the inputs (reading the circuit top to bottom) were \(0\) and \(0\), \(|01\rangle\) means \(0\) and \(1\), and so on.

This pair of qubits is in superposition. If we measured after this step, then each of the four input possibilities would be observed with equal probability. (Use Born’s rule to square the state magnitudes; each is \(1/4\).)

The output

After the \(U_f\) magic happens (next section), we undo the Hadamard on the top path (register) of the circuit. Another tensor product gives us the appropriate matrix, this time of Hadamard and the identity matrix, \(I\), defined:

\(\displaystyle I = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\).

So, \(\displaystyle H \otimes I = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \\ 1 & 0 & -1 & 0 \\ 0 & 1 & 0 & -1 \end{bmatrix}\).

The middle

Now we are left with the four implementations of \(U_f\), the overall structure of which is as follows (Mermin, 2007, p. 42):

Note that the first input, \(|x\rangle\) to the gate remains unchanged by the output and the second out is the exclusive-OR (XOR) of \(y\) and \(f(x)\). (In a future blog post I explain why this helps.) Here are the sums for each \(U_f\) (the pictures are from Mermin, 2007, p. 42; I did the ‘rithmetic):

\(U_{f_0}\):

This is the easiest – it’s just the identity matrix. In the spirit of how we typically create a product of inputs, this is the following tensor product:

\(\displaystyle U_{f_0} = I \otimes I = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \otimes \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}\),

the structure of which pleases me for some reason (how the \(4 \times 4\) identity matrix has a couple of \(2 \times 2\) identity matrices in top-left and bottom-right corners, and how \(\otimes\) makes it happen!).

\(U_{f_1}\):

This is the controlled NOT gate:

\(U_{f_1} = \mbox{CNOT} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix}\).

\(U_{f_2}\):

This takes the NOT (a.k.a. \(X\)) of the second input and then applies CNOT to both inputs. Let’s take the first step first:

\(I \otimes X = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \otimes \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix}\).

(This is mildly pleasing: note the \(2 \times 2\) NOT gate at the top left and bottom right.)

Now multiply CNOT by it:

\(U_{f_2} = \mbox{CNOT} (I \otimes X) = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}\)

\(U_{f_3}\):

Another easier one: just NOT the second input, so we need \(I \otimes X\) and we already calculated it above:

\(U_{f_3} = I \otimes X = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \otimes \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix}\).

Gluing it all together

Now the fun bit in a big table. Let

\(\mbox{Inp} = H|1\rangle \otimes H|1\rangle\) and
\(\mbox{Undo} = H \otimes I\).

That is, \(\mbox{Inp}\) is the input and \(\mbox{Undo}\) is the step in the Deutsch circuit before measurement which applies the Hadamard gate.

Function \(\mbox{Inp}\ U_{f_i}\) \(\mbox{Inp}\ U_{f_i}\ \mbox{Undo}\)
\(f_0\)
(constant)
\(\displaystyle \frac{1}{2} \begin{bmatrix} 1 \\ -1 \\ -1 \\ 1  \end{bmatrix}\) \(\displaystyle \frac{1}{\sqrt{2}} \begin{bmatrix} 0 \\ 0 \\ 1 \\ -1  \end{bmatrix}\)
\(f_1\)
(balanced)
\(\displaystyle \frac{1}{2} \begin{bmatrix} 1 \\ -1 \\ 1 \\ -1  \end{bmatrix}\) \(\displaystyle \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1 \\ 0 \\ 0  \end{bmatrix}\)
\(f_2\)
(balanced)
\(\displaystyle \frac{1}{2} \begin{bmatrix} -1 \\ 1 \\ -1 \\ 1  \end{bmatrix}\) \(\displaystyle \frac{1}{\sqrt{2}} \begin{bmatrix} -1 \\ 1 \\ 0 \\ 0  \end{bmatrix}\)
\(f_3\)
(constant)
\(\displaystyle \frac{1}{2} \begin{bmatrix} -1 \\ 1 \\ 1 \\ -1  \end{bmatrix}\) \(\displaystyle \frac{1}{\sqrt{2}} \begin{bmatrix} 0 \\ 0 \\ -1 \\ 1  \end{bmatrix}\)

To see what we would measure, use Born’s rule and square the amplitudes. This gives the following probabilities for each qubit pair:

\(P(|00\rangle)\) \(P(|01\rangle)\) \(P(|10\rangle)\) \(P(|11\rangle)\)
\(f_0\)
(constant)
\(0\) \(0\) \(\frac{1}{2}\) \(\frac{1}{2}\)
\(f_1\)
(balanced)
\(\frac{1}{2}\) \(\frac{1}{2}\) \(0\) \(0\)
\(f_2\)
(balanced)
\(\frac{1}{2}\) \(\frac{1}{2}\) \(0\) \(0\)
\(f_3\)
(constant)
\(0\) \(0\) \(\frac{1}{2}\) \(\frac{1}{2}\)

If the function is constant, the probability that the first qubit is a \(1\), which I’ll write \(P(|1\_\rangle)\), is:

\(\displaystyle P(|1\_\rangle) = P(|10\rangle) + P(|11\rangle) = \frac{1}{2} + \frac{1}{2} = 1\).

We observe either \(|10\rangle\) or \(|11\rangle\) with equal probability but the first qubit is always \(1\).

If the function is balanced, the probability that the first qubit is a \(0\) is:

\(\displaystyle P(|0\_\rangle) = P(|00\rangle) + P(|01\rangle) = \frac{1}{2} + \frac{1}{2} = 1\).

For this possibility, we observe either \(|00\rangle\) or \(|01\rangle\) with equal probability but the first qubit is always a \(0\).

It worked! Looking at this first qubit tells us whether the function was constant or balanced. The second qubit has a 50-50 chance of being a \(1\) or \(0\) so doesn’t tell us anything helpful.

See this post for the output when I ran the algorithm on an actual quantum computer, ibmq_quito.

But why does it work?!

As I mentioned at the outset, I am not going to attempt to explain why it works, but Mermin’s (2007, p. 41) potted history of the algorithm may be of interest:

“Deutsch’s problem is the simplest example of a quantum tradeoff that sacrifices particular information to acquire relational information. A crude version of it appeared in a 1985 paper by David Deutsch that, together with a 1982 paper by Richard Feynman, launched the whole field. In that early version the trick could be executed successfully only half the time. It took a while for people to realize that the trick could be accomplished every single time.”

References

Flarend, A., & Hilborn, B. (2022). Quantum computing: from Alice to Bob. Oxford University Press.

Mermin, N. D. (2007). Quantum computer science: an introduction. Cambridge University Press.

Quantum computing in R

Quantum computing is developing rapidly. You can now write and run programs on real quantum computers(!) via the web(!!) for free(!!!), e.g., on IBM Quantum. The most powerful quantum computer open to the general public only has 7 qubits. But that’s enough to play with quantum goodies like superposition and entangled states.

All the cool quantum kids use Qiskit, an open-source Python package that simulates programs and talks to IBM’s actual quantum computers so you can run them for real. I quite like R, so was wondering whether it has any quantum computing packages. qsimulatR caught my eye; it simulates programs with up to 24 qubits and automagically translates R descriptions of programs to Qiskit so they can run on IBM’s machines.

I have started reading Alice Flarend and Bob Hilborn’s (2022) intro to quantum computing, so thought I’d give qsimulatR a go to see how easy it is to use it to answer problems on the Bell circuit (pp. 140-143). The Bell circuit is interesting because it demonstrates both superposition and  entanglement.

Here’s a picture of the Bell circuit (drawn automatically by qsimulatR – I’ll get to that shortly):

The H is the Hadamard gate which converts a computational basis state into a superposition state.

\(\displaystyle H = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}\)

Let’s start with that. First, here’s the qsimulatR code for \(|0\rangle\):

> qstate(nbits = 1, coef = c(1,0))
( 1 ) * |0>

And for \(|1\rangle\):

> qstate(nbits = 1, coef = c(0,1))
( 1 ) * |1>

To apply Hadamard to these, we just multiply. Here is \(H|0\rangle = \frac{1}{\sqrt{2}} |0\rangle + \frac{1}{\sqrt{2}} |1\rangle\):

> H(1) * qstate(nbits = 1, coef = c(1,0))
  ( 0.7071068 ) * |0> 
+ ( 0.7071068 ) * |1>

And here is \(H|1\rangle = \frac{1}{\sqrt{2}} |0\rangle- \frac{1}{\sqrt{2}} |1\rangle\):

> H(1) * qstate(nbits = 1, coef = c(0,1))
   ( 0.7071068 ) * |0> 
+ ( -0.7071068 ) * |1>

The second stage of the Bell circuit is the controlled NOT gate (CNOT):

\(\displaystyle \mbox{CNOT} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix}\)

We will need two qubits for CNOT. Here’s how to make the four combinations of computational basis states, by default labelled \(|00\rangle\), \(|01\rangle\), \(|10\rangle\), and \(|11\rangle\):

> qstate(nbits = 2, coefs = c(1,0,0,0))
( 1 ) * |00> 
> qstate(nbits = 2, coefs = c(0,1,0,0))
( 1 ) * |01> 
> qstate(nbits = 2, coefs = c(0,0,1,0))
( 1 ) * |10> 
> qstate(nbits = 2, coefs = c(0,0,0,1))
( 1 ) * |11>

Gluing it all together is easy. The H function has a parameter that tells it which qubit it should be applied to (we want the first, as per the circuit diagram). Suppose the whole circuit is called \(\mbox{BC}\), then \(\mbox{BC}|00\rangle\) is setup using:

> CNOT() * (H(1) *
+           qstate(nbits = 2, coefs = c(1,0,0,0)))
  ( 0.7071068 ) * |00> 
+ ( 0.7071068 ) * |11>

It worked!

\(\displaystyle \mbox{BC}|00\rangle = \frac{1}{\sqrt{2}} |00\rangle+ \frac{1}{\sqrt{2}} |11\rangle\)

Given a single qubit in superposition, we can’t peek inside to see its state since measuring a qubit collapses its state to one of its basis states. Given the entangled pair \(\frac{1}{\sqrt{2}} |00\rangle+ \frac{1}{\sqrt{2}} |11\rangle\), there’s a 50% chance we will observe \(|00\rangle\), 50% chance we will observe \(|11\rangle\), and a 0% chance we will observe the other basis states \(|01\rangle\) or \(|10\rangle\). The probabilities come from Born’s rule – square the amplitudes (actually the absolute value of the amplitudes since in general they’re complex numbers):

\(\left( \frac{1}{\sqrt{2}}\right)^2 = \frac{1}{2}\).

However, we can measure the states of a sequence of qubits that have been prepared in the same way and use the frequencies to estimate superposition states. This is quantum tomography, and qsimulatR allows us to simulate it. Here I have simulated 10,000 trials.

> measure(bc00, rep = 10000) |> summary()
All bits have been measured 10000 times with the outcome:
  ( 4996 ) * |00> 
+ ( 5004 ) * |11>

As expected (with some sampling error), \(P(|00\rangle) \approx P(|11\rangle) \approx \frac{1}{2}\).

We can wrap the Bell circuit in a function like so:

bell <- function(coefs) {
  inputs <- qstate(nbits = 2,
                   coefs = coefs,
                   basis = c("|00>","|10>","|01>","|11>"))
  
  cat("Input: ")
  print(inputs)
  cat("\n")

  CNOT() * (H(1) * inputs)
}

Note, I have explicitly named the basis states to match the bit order used in the book, which differs from the default order in the package (the package numbers bits right to left, the book uses left to right).

Here are all four combinations:

> bell(c(1,0,0,0))
Input:    ( 1 )	* |00> 

   ( 0.7071068 )	* |00> 
 + ( 0.7071068 )	* |11> 
 
> bell(c(0,1,0,0))
Input:    ( 1 )	* |10> 

   ( 0.7071068 )	* |00> 
 + ( -0.7071068 )	* |11> 
 
> bell(c(0,0,1,0))
Input:    ( 1 )	* |01> 

   ( 0.7071068 )	* |10> 
 + ( 0.7071068 )	* |01> 
 
> bell(c(0,0,0,1))
Input:    ( 1 )	* |11> 

   ( -0.7071068 )	* |10> 
 + ( 0.7071068 )	* |01>

Exercise (“Try It”) 10.6 asks us to reverse \(\mbox{BC}|00\rangle\), i.e., the reverse Bell circuit. That’s easy to do:

> RBC <- H(1) * (CNOT() * bell(c(1,0,0,0)))
> RBC
   ( 1 )	* |00>

So we’re back to \(|00\rangle = |0\rangle \otimes |0\rangle\).

The very cool thing about qsimulatR is that it draws circuit diagrams for us. Here’s the Bell circuit piped into the reverse Bell circuit:

RBC |> plot()

If that was very cool, then this is surely exceedingly cool – a one-liner to export the R object to IBM’s Qiskit Python format, so you can run it on one of their quantum computers:

export2qiskit(RBC, filename = "rbc.py")

This gives:

# automatically generated by qsimulatR
qc = QuantumCircuit(2)
qc.h(0)
qc.cx(0, 1)
qc.cx(0, 1)
qc.h(0)

A bit of trial and error and copy and paste and ta-da, it worked. IBM Quantum Lab parsed the code and drew the circuit diagram. A small tweak was needed to tell it the name of the circuit:

Now, I’m trying to wrap my head around the Deutsch–Jozsa algorithm. Watch this space…

Edited to add: here’s what happened next.

References

Flarend, A., & Hilborn, B. (2022). Quantum computing: from Alice to Bob. Oxford University Press.