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.