Build a Quantum Computing Simulator from Scratch — Part 2
From Hadamard gates to quantum circuits
Welcome back to the next part of building a quantum simulator from scratch (in Racket). This is part 2 of the series, here is Part 1. In this part, we will complete the functional core of the simulator. We will see:
- How to measure state vectors and count the results.
- How to create Quantum gates.
- How to combine state vectors and quantum gates via tensor product.
- How to combine all of these into a function that makes circuits for you.
Whew, let's get started!
Measurement of Qubits
In the previous part, we saw how to use a column vector to represent a qubit. Let's recap quickly. A column vector is a list of complex numbers arranged vertically. A single-column vector can represent any combination of qubits, if we have n qubits in the system, we need 2^n rows in the column vector. Thus, a qubit |q1> (n = 1), is represented by:
[a,
b]
Where |a|^2 + |b|^2 = 1The values a and b are known as probability amplitudes and the sum of squares of those amplitudes should be 1 (Born’s rule). For more than one qubit system the rule is the same. Here’s a column vector for n=3.
[a1,
a2,
a3,
a4,
a5,
a6,
a7,
a8]
Where Sum(magnitude(a(i))^2; i = 1 to 8) = 1 These column vectors are also known as state vectors as they represent the state of a qubit system before measurement.
Did you notice we used the term before measurement? All the probabilities of a state vector exist before the qubits are measured. Once they are measured, however, a qubit is either 0 or 1. So the result of measurement of an n qubit system is a string of 0s and 1s, just like a classical system. This phenomenon is one of the wonders of quantum physics which is still not completely understood. The popular term for this is known as Collapse of the Wave function.
Ok, so how do we implement the measurement of a quantum system? We will adopt a very simple mechanism:
- In a state vector, if a row is 0, then we know that binary string is unlikely ever to appear after measurement.
- For the non-zero ones, we square them to get their actual probabilities.
- Then we scale them by a factor of 100.
- After scaling, together the segments would cover the entire scale of 100 units. This is because, before scaling, they would together make 1 according to Born’s rule.
- Then we roll a dice and see under which segment the dice falls, that would be our result of measurement. The bigger the segment, the more likely the dice would fall in that segment.
Why 100? Just to make the segments lengthier and make it easy for the dice to fall there.
All right let us implement this part:
(define (range-map start probs obs roll i n)
(if (= i n)
(error "No observation from circuit - invalid state")
(if (= (list-ref probs i) 0)
(range-map start
probs
obs
roll
(+ i 1)
n)
(if (and (>= roll start)
(< roll (+ start (list-ref probs i))))
(list-ref obs i)
(range-map (+ start (list-ref probs i)) ; next start
probs
obs
roll
(+ i 1) ; next observation to match
n)))))
(define (measure-mat m)
(let* ([len (matrix-num-rows m)]
[blen (exact-round (log len 2))]
[obs (bits-of-len blen)]
[ampls (matrix->list m)]
[scaled-probs (map (λ (x) (* x x 100)) ampls)]
[roll (random 100)])
(range-map 0 scaled-probs obs roll 0 len)))The heavy-loading of the operation is done by a recursive function range-map and measure-mat drives it. Note that, a state vector is just a matrix of n rows and 1 column in Racket, see Part 1 for how we are constructing the matrix.
We are not quite done yet. Measuring a qubit only once would give us some result according to the probabilities, but since it is ultimately a rolling of dice, the dice may not fall on the bigger segment. Hence, we may not get the correct answer on measuring the qubit system only once. In practice, a quantum system is measured multiple times, and the resulting histogram indicates what are the results of the system. We need a function that calls measure-mat many times for a state vector and combine the results with their counts.
(define (counts q #:shots [shots 1024])
(let ([results (make-hash)])
(for ([i (range shots)])
(let ([obs (measure-mat q)])
(if (hash-has-key? results obs)
(hash-set! results obs (+ (hash-ref results obs) 1))
(hash-set! results obs 1))))
results))By default, we measure 1024 times, but we can vary this based on our requirements. Let us take this function for a spin.
(counts (qubit 1 0)) ; ==> '#hash(((0) . 1024))
(counts (qubit 0 1)) ; ==> '#hash(((1) . 1024))Unsurprisingly, measuring a |0> will always produce 0 and a |1> will produce 1. How about the special state we discussed last time?
(define h-factor (/ 1.0 (sqrt 2)))
(counts (qubit h-factor h-factor)) ; ==> '#hash(((0) . 499) ((1) . 525))The h-factor simply sets the probabilities 1/2 for each state (|0>, |1>) to appear. It is the superposition state where either result is possible. And that’s what we also see in the measurement — almost half the time we got 0 and half times we got 1.
We will use the counts function from now on to determine the results of a quantum system.
Quantum Gates
Just having qubits is no fun. We need to do some operations on them too. That’s where the quantum gates enter. They get the name gate from their classical computing counterparts like AND-gate, OR-gate, and NOT-gate. Just like classical computing, quantum computing is also made up of these quantum gates working on qubits and circuits made of quantum gates. But unlike a classical boolean gate, quantum gates have some properties that must be satisfied.
- Quantum gates can be represented by unitary matrices. If the operation of a quantum gate cannot be expressed via a corresponding unitary matrix, then it may not be considered as a valid gate at all.
- The operations of quantum gates are completely reversible. This is very different from boolean gates. In an AND gate, looking at the results we cannot determine the inputs. We can only say, if the result is 1, then both inputs are 1. But if the result is 0, we don’t know if the inputs are (1, 0), (0, 1), or (0, 0). But for a quantum gate, the input can be reconstructed by looking at the output alone.
Since quantum gates are just matrices, we can easily construct them in our code by defining the corresponding matrix. And applying a gate is just multiplying the gate with the qubit.
new_state_vector = U x old_state_vector
Where U = unitary matrix, quantum gateWe start by defining a function to apply a gate followed by four Pauli matrices and the venerable Hadamard gate.
(define ((apply-op op-matrix) q)
(matrix* op-matrix q))
(define pauli-x (matrix [[0 1]
[1 0]]))
(define gX (apply-op pauli-x))
(define pauli-z (matrix [[1 0]
[0 -1]]))
(define gZ (apply-op pauli-z))
(define pauli-y (matrix [[0 0-i]
[0+i 0]]))
(define gY (apply-op pauli-y))
(define hadamard (matrix [[h-factor h-factor]
[h-factor (- h-factor)]]))
(define gH (apply-op hadamard))
(define q0 (qubit 1 0)) ; just qubit |0>
;; shortcuts
(define H hadamard)
(define X pauli-x)
(define Y pauli-y)
(define Z pauli-z)apply-op is a higher-order function that takes a matrix gate and returns a function. This function can take any qubit and produce the results. Out of these gates, X and H (or Pauli X and Hadamard) are more interesting ones for now. X gate is also known as quantum NOT gate. So, it takes a qubit and reverses its state. The H or Hadamard gate takes a qubit and puts it into a superposition state. Let’s see examples of both.
;; X |0> == |1>
(gX q0) ; ==> (array #[#[0] #[1]])
;; X |1> == |0>
(gX (qubit 0 1)) ; ==> (array #[#[1] #[0]])
;; X X |0> == |0>
(gX (gX q0)) ; ==> (array #[#[1] #[0]])
;; H |0> == equally likely superposition
(gH q0) ; ==> (array #[#[0.7071067811865475] #[0.7071067811865475]])
;; X H |0> ==> same superposition
(gX (gH q0)) ; ==> (array #[#[0.7071067811865475] #[0.7071067811865475]])
;; Measure (H |0>)
(counts (gH q0)) ; ==> '#hash(((0) . 516) ((1) . 508))
;; Measure (X |0>)
(counts (gX q0)) ; ==> '#hash(((1) . 1024))⊗ — Combining state vectors and gates
We have seen how to use gates on single qubits, but can we use them on multiple qubits? Before we can do that, we need a way to combine multiple qubits and form their state vector. We have done this manually so far, but we need a way to do this automatically. Turns out, we can do this using a special matrix operation called tensor product represented by the symbol ⊗.
The details of tensor-product can be found in the Wiki, but here’s a simplified textual representation.
[a b [p q r
c d] tensor* x y z]
==> [ a * [p q r b * [p q r
x y z] x y z]
c * [p q r d * [p q r
x y z] x y z] ]
==> [ ap aq ar bp bq br
ax ay az bx by bz
cp cq cr dp dq dr
cx cy cz dx dy dz ]And here’s the implementation. Granted it is not the best one and not a Lisp-y one too, but it does the work for now.
(define (tensor* m1 m2)
(define (mat->list m1 m2)
(matrix->list*
(matrix-map (λ (x)
(matrix->list*
(matrix-map (λ (y)
(* x y))
m2)))
m1)))
(define (list->mat m)
(let ([rows '()])
(for ([row m])
(for ([i (range (length (first row)))])
(let ([row-line '()])
(for ([m row])
(set! row-line (append row-line (list-ref m i))))
(set! rows (append rows (list row-line))))))
rows))
(let* ([m (mat->list m1 m2)]
[l (list->mat m)])
(list*->matrix l)))
(define (t* . qbits)
(let ([n (length qbits)])
(cond
((< n 1) (error "No qubits given"))
((= n 1) (first qbits))
(else (tensor* (first qbits)
(apply t* (rest qbits)))))))The t* function works on any number of matrices and produces the tensor product of all of those. This will come in handy when we make a circuit.
Let's invoke it a few times.
;; |0> ⊗ |0> == |00>
(t* q0
q0) ; ==> (array #[#[1] #[0] #[0] #[0]])
;; |0> ⊗ |1> == |01>
(t* q0
(gX q0)) ; ==> (array #[#[0] #[1] #[0] #[0]])
;; |1> ⊗ |0> == |10>
(t* (gX q0)
q0) ; ==> (array #[#[0] #[0] #[1] #[0]])
;; |1> ⊗ |1> == |11>
(t* (gX q0)
(gX q0)) ; ==> (array #[#[0] #[0] #[0] #[1]])
;; (H q0) ⊗ |1>
(t* (gH q0)
q0) ; ==> (array #[#[0.7071067811865475] #[0] #[0.7071067811865475] #[0]])
(counts (t* (gH q0)
q0)) ; ==> '#hash(((0 0) . 518) ((1 0) . 506))The last one is interesting — if we pair up one qubit in H a superposition state with another that is just |0>, the measurements show that, half the time we will get |00> and the other half |10>. Now let's make it more interesting.
(counts (t* (gH q0)
(gH q0)))
; ==>'#hash(((0 1) . 260) ((1 1) . 261) ((0 0) . 245) ((1 0) . 258))When both the qubits are in a superposition state the measurements now show all possible combinations to be equally likely!
Multiple gates on multiple qubits
We saw how we can combine states of multiple qubits (which can be the result of applying single qubit gates to the qubits individually, like the above example where we apply H a gate on each of the qubits and then combine them). But how do we combine gates so that they can act as a single operator on multiple qubit state vectors? Turns out the answer is the same — just tensor product the gates! Let’s implement a function that creates a bigger gate by applying tensor products to multiple gates.
(define (G* . ms)
(apply-op (apply t* ms)))That was simple. One more thing, we assume that G* will get a list of gate matrices like H or X — but what if we just have one gate for one qubit, and the other qubits just need flow through without applying any transformation? In that case we apply the Identity gate. An identity gate is a square matrix whose diagonals are all 1. This is also a unitary matrix so it can be a quantum gate. Lets write a helper function to make us an identity gate.
(define (I n)
(identity-matrix n))Now we can create a gate that applies H gate to the first qubit and nothing to the second qubit (that is identity) and we can mesure the results.
┌───┐
q_0: ┤ H ├
└───┘
q_1: ─────(define new-gate (G* H
(I 2)))
(define initial-qubits (t* q0
q0))
(define state-vector (new-gate initial-qubits))
state-vector ; ==> (array #[#[0.7071067811865475] #[0] #[0.7071067811865475] #[0]])Little-endian vs Big-endian
What is the eqivalent code for the same thing above in Qiskit? It looks like this:
from qiskit import *
from qiskit_aer import StatevectorSimulator
circuit = QuantumCircuit(2)
circuit.h(0)
simulator = StatevectorSimulator()
result = simulator.run(circuit).result()
print(result.get_statevector())The output looks like this:
Statevector([0.70710678+0.j, 0.70710678+0.j, 0. +0.j,
0. +0.j],
dims=(2, 2))Ignoring the 0s in the vector, it looks like our state vector, but not quite the same. Seems the position of the non-zero elements are different. What’s going on?
Turns out, the order of the matrices matter! In our current order, i.e. the H above and (I 2) below, it makes it a big endian system, where the least significant bit is on the left. However, qiskit and others like Q# uses little endian systems - least significant bit is on the right. Although in the circuit diagram we see the least bit at the top, it is intepreted to be on the right! It is confusing but after a bit of back and forth, it becomes normal (well hopefully). So our code is right except our order needs to be opposite for it to match qiskit output. Lets try to modify our G* function to follow the same convention and see what comes out.
(define (G* . ms) ;; qiskit compatible
(apply-op (apply t* (reverse ms))))Now running the same Racket code, we get the below statevector
(array #[#[0.7071067811865475] #[0.7071067811865475] #[0] #[0]])Now it matches with qiskit output. In general, the little endian is preferable format and we will stick to it, but it is easy to switch to big endian by just removing the reverse from G*. In fact, the library in Github defines another function nG* to do the same.
Quantum Circuits
Now we have all the ingredients to write a function that creates a quantum circuit.
(define (make-circuit matrices #:assembler [assembler G*])
(let ([ops (map (λ (gs) (if (list? gs)
(apply assembler gs)
(apply-op gs)))
matrices)])
(λ (input-qbits)
(for/fold ([sv input-qbits])
([f ops])
(f sv)))))
(define (qubits n)
(apply t* (for/list ([i (range n)])
q0)))The helper function qubits creates an input state vector with all qubits as zero.
We will now conclude this part by trying to make a simple circuit like above.
(define c1 (make-circuit
(list (list (I 2)
X)
(list H
H))))The equivalent Qiskit code:
circuit = QuantumCircuit(2)
circuit.x(1)
circuit.barrier()
circuit.h(0)
circuit.h(1)
circuit.draw('mpl')And this is how we apply the circuit to a set of qubits in base condition and finally measure the output.
(counts (c1 (qubits 2)))
; ==> '#hash(((0 1) . 255) ((0 0) . 272) ((1 1) . 252) ((1 0) . 245))Simple isn’t it!
In the next part we will introduce CNOT gate, entanglement and more advanced concepts.
Be sure to check out Part 1 of the series.
Check out the code for the complete simulator.
If you like this, do give a few claps and follow for more. Thanks for reading!
