Implementing Deutsch’s Algorithm in qiskit and cirq
Check out the code for this post here: github.com/przecze/deutschs_algorithm_in_qiskit_vs_cirq
Deutsch’s algorithm is the simplest example that shows the subtle power of Quantum Computing. I decided to try to implement the simplest version, one qubit version of it using two popular quantum computing libraries: qiskit and cirq. This turned out to be a very interesting hands-on exercise to get started with quantum computing in Python.
Deutsch’s algorithm
Deustch’s algorithm, is also known as Deutsch-Jozsa algorithm. Given a binary function f, assuming f is either unbalanced (constant) or balanced (returning 0 and 1 equally often), the algorithm enables us to distinguish between the two options in a single evaluation. The quantum circuit looks as follows:
The idea behind the circuit is: we input the superposition of all possible inputs (prepared with the Hadamard gate) and we receive the superposition of all possible outputs. Then we measure the output state in a special basis which lets us extract only the information if f is balanced or not without knowing the specific output values.
Let’s try implementing the algorithm in qiskit
Implementing Deutsch’s algorithm in qiskit
We start with creating the “preparation” circuit:
import qiskit as qs
from qiskit import QuantumCircuit, QuantumRegister,ClassicalRegisterregs = [QuantumRegister(2, 'q'), ClassicalRegister(1, 'c')]
# state initialization: apply x on output bit, and H on all
init = QuantumCircuit(*regs)
init.x(1)
init.h(0)
init.h(1)
init.barrier() # only for separation in printing
qiskit has built-in printing of circuits in ascii art:
In [2]: print(init)
┌───┐ ░
q_0: ┤ H ├──────░─
├───┤┌───┐ ░
q_1: ┤ X ├┤ H ├─░─
└───┘└───┘ ░
c: 1/═════════════Which is very nice and useful for debugging. The “measure” circuit will look like this:
# at the end, apply h to input bit and measure it
end = QuantumCircuit(*regs)
end.h(0)
end.measure(0, 0)Another fun feature: circuits can be combined with simple + operator:
In [3]: print(init + end)
┌───┐ ░ ┌───┐┌─┐
q_0: ┤ H ├──────░─┤ H ├┤M├
├───┤┌───┐ ░ └───┘└╥┘
q_1: ┤ X ├┤ H ├─░───────╫─
└───┘└───┘ ░ ║
c: 1/═══════════════════╩═
0Of course, we are still missing here the most interesting part: function f itself. Let’s implement the two cases separately
Unbalanced function
An unbalanced function means the output is always the same. In other words, if we consider the standard reversible representation of a quantum function: f(|x, y⟩) = |x, f(x)⊕ y⟩ we want f(x) to always be 0 or always be 1. That means we have two possible cases:
unbalanced0: f(|x, y⟩) = |x, f(x)⊕ y⟩ = |x, 0⊕ y⟩ = |x, y⟩
This means our circuit for f will be simply an identity:
unbalanced0 = QuantumCircuit(*regs)
unbalanced0.barrier()unbalanced1: f(|x, y⟩) = |x, f(x)⊕ y⟩ = |x, 1⊕ y⟩ = |x,~y⟩
So we always flip y regardless of x. A quantum circuit will simply apply a NOT to the y register:
unbalanced1 = QuantumCircuit(*regs)
unbalanced1.x(1)
unbalanced1.barrier()Note: NOT gate is called .x is qiskit. This is because applying a NOT can also be represented as “rotating the qubit 180° around the x-axis in the Bloch sphere”. This has nothing to do with our notation where x represents the input qubit.
Let’s see the two resulting circuits:
In [4]: print(init + unbalanced0 + end)
┌───┐ ░ ░ ┌───┐┌─┐
q_0: ┤ H ├──────░──░─┤ H ├┤M├
├───┤┌───┐ ░ ░ └───┘└╥┘
q_1: ┤ X ├┤ H ├─░──░───────╫─
└───┘└───┘ ░ ░ ║
c: 1/══════════════════════╩═
0
In [5]: print(init + unbalanced1 + end)
┌───┐ ░ ░ ┌───┐┌─┐
q_0: ┤ H ├──────░───────░─┤ H ├┤M├
├───┤┌───┐ ░ ┌───┐ ░ └───┘└╥┘
q_1: ┤ X ├┤ H ├─░─┤ X ├─░───────╫─
└───┘└───┘ ░ └───┘ ░ ║
c: 1/═══════════════════════════╩═
0Balanced function
A balanced function is a tiny bit trickier: f(x) actually depends on x now. Let’s consider the simplest case of f(x)=x and see what it means to our y register: f(∣x, y⟩) = |x, f(x)⊕ y⟩ = |x, x⊕ y⟩
Which is no less and no more than the well-known quantum CNOT gate. A NOT operation is applied on the second register depending on the state of the first register. If |x⟩ is |1⟩ we flip the |y⟩, otherwise we leave |y⟩ unchanged. Let’s create the balanced subcircuit:
balanced = QuantumCircuit(*regs)
balanced.cx(0,1)
balanced.barrier()Note: as NOT is also called x , CNOT is applied with cx . See the explanation above.
Let’s see the full circuit for the balanced function:
In [6]: print(init + balanced + end)
┌───┐ ░ ░ ┌───┐┌─┐
q_0: ┤ H ├──────░───■───░─┤ H ├┤M├
├───┤┌───┐ ░ ┌─┴─┐ ░ └───┘└╥┘
q_1: ┤ X ├┤ H ├─░─┤ X ├─░───────╫─
└───┘└───┘ ░ └───┘ ░ ║
c: 1/═══════════════════════════╩═
0Evaluate Deutsch’s algorithm
The premise of Deutsch’s algorithm is that we can distinguish balanced from unbalanced function in a single run. In our code we will run 10 shots, just to make sure there’s no probability involved and that for the given function the result is always the same.
To evaluate the circuit in qiskit, we have to transpile it, feed it to a simulator and extract the results:
sim = QasmSimulator()
# expected: only 1s for balanced, only 0s for unbalanced
for kind, oracle in (('b', balanced),
('u', unbalanced0),
('u', unbalanced1)):
circuit = init + oracle + end
transpiled = qs.transpile(circuit)
counts = sim.run(transpiled, shots=10).result().get_counts()
print(kind, counts)Before we test our code, let’s take a look once again at the full script:
import qiskit as qs
from qiskit import QuantumCircuit, QuantumRegister,ClassicalRegister
from qiskit.providers.aer import QasmSimulator
from qiskit import transpileregs = [QuantumRegister(2, 'q'), ClassicalRegister(1, 'c')]# state initialization: apply x on output bit, and H on all
init = QuantumCircuit(*regs)
init.x(1)
init.h(0)
init.h(1)
init.barrier()# balanced is a c-not - output depends on input
balanced = QuantumCircuit(*regs)
balanced.cx(0,1)
balanced.barrier()# unbalanced0 is the unity
unbalanced0 = QuantumCircuit(*regs)
unbalanced0.barrier()# unbalanced1 is a 'not' - output is flipped independent of input
unbalanced1 = QuantumCircuit(*regs)
unbalanced1.x(1)
unbalanced1.barrier()# at the end, apply h to input bit and measure it
end = QuantumCircuit(*regs)
end.h(0)
end.measure(0, 0)sim = QasmSimulator()
# expected: only 1s for balanced, only 0s for unbalanced
for kind, oracle in (('b', balanced),
('u', unbalanced0),
('u', unbalanced1)):
circuit = init + oracle + end
transpiled = qs.transpile(circuit)
counts = sim.run(transpiled, shots=10).result().get_counts()
print(kind, counts)
And running the script results in… (drumroll)
In [9]: %run deutsch_in_qiskit.py
b {'1': 10}
u {'0': 10}
u {'0': 10}Full success! Having this simple script as a baseline, let’s try to do the same in the cirq library.
Implementing Deutsch’s algorithm in cirq
Building circuits is a bit different here than it was in qiskit: you have to create each gate separately and pass it to the circuit constructor:
import cirq# Create the qubits
qubits = [cirq.LineQubit(i) for i in (0, 1)]# Create a circuit
init = cirq.Circuit(
cirq.H(qubits[0]),
cirq.X(qubits[1]),
cirq.H(qubits[1])
)end = cirq.Circuit(
cirq.H(qubits[0]),
cirq.measure(qubits[0], key='m')
)
The good news is, the trick with adding two circuits works just as well. Also, cirq just like qiskit has some build-in ascii visualization:
In [17]: print(init)
0: ───H───────
1: ───X───H───
In [18]: print(end)
0: ───H───M('m')───
In [19]: print(init + end)
0: ───H───────H───M('m')───
1: ───X───H────────────────Arguably more primitive than the qiskit block diagrams. But, as they say, simple is beautiful.
We already know the story behind our balanced and unbalanced subcircuits, so we can write down the cirq versions without further ado:
balanced = cirq.Circuit(
cirq.CNOT(*qubits)
)
unbalanced0 = cirq.Circuit()unbalanced1 = cirq.Circuit(
cirq.X(qubits[1])
)
And the evaluation code is much simpler: no need to transpile the circuit, and you can print the result right away to get what you need:
# Simulate the circuit several times.
simulator = cirq.Simulator()
# expected: only 1s for balanced, only 0s for unbalanced
for kind, oracle in (('b', balanced),
('u', unbalanced0),
('u', unbalanced1)):
result = simulator.run(init + oracle + end, repetitions=10)
print(kind, result) Before running the script let’s look at it in its fullness:
import cirq# Create the qubits
qubits = [cirq.LineQubit(0), cirq.LineQubit(1)]# Create a circuit
init = cirq.Circuit(
cirq.H(qubits[0]),
cirq.X(qubits[1]),
cirq.H(qubits[1])
)balanced = cirq.Circuit(
cirq.CNOT(*qubits)
)unbalanced0 = cirq.Circuit()unbalanced1 = cirq.Circuit(
cirq.X(qubits[1])
)end = cirq.Circuit(
cirq.H(qubits[0]),
cirq.measure(qubits[0], key='m')
)# Simulate the circuit several times.
simulator = cirq.Simulator()
# expected: only 1s for balanced, only 0s for unbalanced
for kind, oracle in (('b', balanced),
('u', unbalanced0),
('u', unbalanced1)):
result = simulator.run(init + oracle + end, repetitions=10)
print(kind, result)
Executing the code above results in:
In [23]: %run deutsch_in_cirq.py
b m=1111111111
u m=0000000000
u m=0000000000Debugging the circuit — getting the state vector
When preparing the circuits above, I didn’t get it right the first time. I made some mistakes and the measurements didn’t come out as I expected them to. At this point, I thought it would be very helpful to check the full quantum state after each subcircuit (something that you’re fundamentally not allowed to do on a real quantum machine — but there’s no reason not to do it on a simulator). This turned out to be pretty complicated in qiskit. You have to do the following:
- add a special “snapshot” gate to your circuit
- transpile the circuit and run it
- extract the snapshot state vector as an array from the result
- pretty print the complex array
To this end I created some helper functions:
from qiskit.providers.aer.extensions.snapshot import Snapshotdef to_braket(array):
""" helper for pretty printing """
state = []
basis = ('|00>', '|10>', '|01>', '|11>')
for im, base_state in zip(array, basis):
if im:
if not im.imag:
state.append(f'{im.real:.1f}{base_state}')
else:
state.append(f'({im:.1f}){base_state}')
return " + ".join(state)
def get_output_state(circuit, quantum_registers):
""" returns the braket-formatted output state of a circuit """
c = circuit.copy()
c.append(Snapshot('snap', 'statevector', 2), *quantum_registers)
sim = QasmSimulator()
result = sim.run(qs.transpile(c, sim), shots=1).result()
snap = result.data()['snapshots']['statevector']['snap'][0]
return to_braket(snap)
Which then can be used for debugging. For example, to see that the init subcircuit properly initializes the state as the superposition of all possible |00⟩, …,|11⟩ inputs:
In [36]: get_output_state(init, regs[0])
Out[36]: '0.5|00> + 0.5|10> + (-0.5+0.0j)|01> + (-0.5+0.0j)|11>'Or that after running through a full circuit with a balanced oracle, x qubit is always in |1⟩ state:
In [37]: get_output_state(init + balanced + end, regs[0])
Out[37]: '(0.7-0.0j)|10> + (-0.7+0.0j)|11>'After needing to hack around the qiskit to get the very simple debug information I wanted, I tried to do the same thing in cirq. And I was almost blown away:
def get_output_state(circuit):
""" That's all it takes! """
return to_braket(circuit.final_state_vector())It works just as well. And there’s no need to pass around the registers:
In [47]: get_output_state(init)
Out[47]: '0.5|00> + -0.5|10> + 0.5|01> + -0.5|11>'
In [48]: get_output_state(init + balanced + end)
Out[48]: '0.7|01> + -0.7|11>'Conclusions
We took the simplest, one-qubit version of Deutsch’s algorithm and tried to see it in action using both qiskit and cirq libraries. Based on this simple example, in comparison to more popular qiskit, cirq seems more primitive at first glance, but after playing with it for a bit it turns out to have a more intuitive interface and much easier access to the simulated quantum state.
Based on that, for simple projects and learning quantum computing I recommend starting with cirq right away, without bothering with qiskit.
