Example: quantum circuits

This example is adapted from a blogpost by Carlo Lepelaars and originally appeared on his personal website. It is reposted here with permission.

import numpy as np
from numpy.testing import assert_almost_equal
from fasttransform import Transform, Pipeline

Transform

Transform is a fundamental building block to manipulate data in Python. While extremely simple to use, it is also flexible. Transforms can be set up to change behavior based on the input type, also called type dispatch. Transforms can also be made reversible. Keep this reversibility in mind for the quantum part later! Here is a simple example of a Transform that can square a number and take square root as reverse:

class S(Transform):
    def encodes(self, x): return x ** 2
    def decodes(self, x): return x ** 0.5
            
assert S()(10) == 100
assert S().decode(100) == 10
assert S().decode(S()(10)) == 10

A Transform with only encode can be defined even simpler with a lambda function. By default, decode returns its input (i.e. does nothing / “no-op”):

square = Transform(lambda x: x ** 2)
assert square(10) == 100
assert square.decode(100) == 100

Transform can also be used as a decorator to turn a function into a Transform:

@Transform
def square(x): return x ** 2
square(10) # 100
type(square) # <class 'fastcore.transform.Transform'>

fasttransform.transform.Transform

A powerful feature of Transform is type dispatch:

class MultiS(Transform):
    def encodes(self, x: int | float | complex | tuple): return x**2
    def encodes(self, x: list): return [x**2 for x in x]
    def decodes(self, x: int | float | complex | tuple): return x**0.5
    def decodes(self, x: list): return [x**0.5 for x in x]

ms = MultiS()
ms

MultiS(enc:2,dec:2)

# Lists
# By default, Transform processes lists as a whole
# 2nd encodes method is called
assert ms([1, 2, 3]) == [1, 4, 9]

# 2nd decodes method is called
assert ms.decode([1, 4, 9]) == [1.0, 2.0, 3.0]

# Tuples
# By default, Transform processes tuples elementwise
# 1st decodes method is called
assert ms((1, 2, 3)) == (1, 4, 9)

# 1st decodes method is called
assert ms.decode((1, 4, 9)) == (1.0, 2.0, 3.0)

# Complex numbers
# 1st encodes method is called on complex number
assert ms(10.0j) == (-100+0j)

# 1st decodes method is called on complex number
assert ms.decode(ms(10.0j)) == (6.123233995736766e-16+10j)

Transform automatically routes input to the corresponding method based on the input type. Transform will operate on lists as a whole, while tuples are processed elementwise. Processing lists and arrays as one object is a common use case in data science, for example if we are processing batches of images (3D arrays) for machine learning. If you want a Transform that also transforms tuples as a whole, use ItemTransform. To support inplace transforms use InplaceTransform.

Pipeline

Now that we have a firm understanding of Transform, Pipeline can be used to chain them together:

class S(Transform):
    "Square a number. Reverse is square root."
    def encodes(self, x): return x ** 2
    def decodes(self, x): return x ** 0.5
    
class A(Transform):
    "Add 1. Reverse is subtract 1."
    def encodes(self, x): return x + 1
    def decodes(self, x): return x - 1

pipe = Pipeline([S(), A()])
assert pipe(10) == 101 # 10**2 + 1 = 101
assert pipe.decode(10) == 3.0 # (10 - 1)**0.5 = 3.0
assert pipe.decode(pipe(10)) == 10 # (10**2 + 1 - 1)**0.5 = 10

I hope you appreciate the simplicity of Transform and Pipeline. This year I’ve been working on an open-source library for quantum computing that connects several quantum computing frameworks and standards (For example, Qiskit, PennyLane and OpenQASM. At the base of it lies numpy. To generalize quantum circuits I tried using scikit-learn’s Pipeline functionality. Unfortunately this led to a lot of boilerplate code and unnecessary features. fastcore offers a more elegant and promising solution for the use case of constructing quantum circuits, which we will explore in the next section.

Quantum

Quantum computing processes can be simulated on classical computers by transforming a statevector (i.e. list of complex numbers) through a series of reversible quantum logic gates (i.e. matrix of complex numbers). The state and gates are subject to constraints, but the basic operation is (reversible) vector-matrix multiplication. This is a perfect use case for fastcore’s Transform and Pipeline. To illustrate this point we will start with manipulating a single qubit using Transform.

Single Qubit

A statevector contains 2 complex numbers which gives us the likelihood of obtaining a \(0\) or \(1\). This is called a superposition between \(0\) and \(1\). The numbers in the vector are called probability amplitudes and can be converted to probabilities by computing \(|x|^2\) (i.e. the absolute value squared), where \(x\) is the statevector. Probabilities must always sum to \(1\), so a valid qubit state \([\alpha, \beta]\) must have \(|\alpha|^2 + |\beta|^2 = 1\).

A shortcut used for writing quantum states is Dirac notation. For example \(|0\rangle=[1, 0]^T\) (i.e. always 0) and \(|1\rangle=[0, 1]^T\) (i.e. always 1). Other valid single qubit states include \(\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle) = \begin{bmatrix} \frac{1}{\sqrt{2}} \frac{1}{\sqrt{2}} \end{bmatrix}^T\) for a perfectly equal superposition and \(\frac{1+i}{2}|0\rangle + \frac{1-i}{2}|1\rangle = \begin{bmatrix} \frac{1+i}{2} \frac{1-i}{2} \end{bmatrix}^T\) for a superposition of both the real and complex parts.

Using Transform we can easily define quantum logic gates. We define a base Transform that can do vector-matrix multiplication in a reversible way and implement common quantum gates. The I (identity) gate does nothing, the X (NOT) gate flips the qubit from \(|0\rangle\) to \(|1\rangle\) and vice versa. The Hadamard gate turns a qubit into superposition (i.e. turn \(|0\rangle\) into \(\frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)\)).

class _Q(Transform):
    "Base transform for quantum gates"
    def encodes(self, x): return x @ self.gate
    def decodes(self, x): return x @ self.gate.conj().T
                     
class I(_Q):
    "Identity gate. Does nothing."
    gate = np.array([[1, 0], 
                     [0, 1]])
    
class X(_Q):
    "X (NOT) gate. Flips from |0> to |1> and vice versa."
    gate = np.array([[0, 1], 
                     [1, 0]])
    
class H(_Q):
    "Hadamard (Superposition) gate. Turns a qubit into a superposition."
    gate = np.array([[1, 1], 
                     [1, -1]]) / np.sqrt(2)

This allows us to easily play with quantum gates:

zero_state = [1+0j, 0+0j] # Basis state |0>
superposition_state = np.array([0.5+0.5j, 0.5-0.5j])  # Complex superposition

# Identity operation
i = I() 
assert_almost_equal(i(zero_state), np.array([1.+0.j, 0.+0.j])) # (|0>)
assert_almost_equal(i.decode(zero_state), np.array([1.+0.j, 0.+0.j]))  # (|0>)

# X (NOT) operation 
x = X()
assert_almost_equal(x(zero_state), np.array([0+0j, 1+0j]))  # (|1>)
assert_almost_equal(x.decode(x(zero_state)), np.array([1.+0.j, 0.+0.j]))  # (|0>)
assert_almost_equal(x(superposition_state), np.array([0.5-0.5j, 0.5+0.5j]))  # (flips sign of complex part)

# Hadamard gate tests
h = H()
assert_almost_equal(h(zero_state), np.array([0.707+0.j, 0.707+0.j]), decimal=3)  # (superposition)
assert_almost_equal(h(superposition_state), np.array([0.707+0.j, 0.+0.707j]), decimal=3)  # (phase state)
assert_almost_equal(h.decode(h(superposition_state)), np.array([0.5+0.5j, 0.5-0.5j]))  # (complex superposition)

I hope this gives you an appreciation for the simplicity of Transform and why it works for quantum, even if the details and meaning of these operations might not be clear yet. If you are interested in learning about quantum computing I highly recommend reading the amazing educational material on quantum.country. A great quantum textbook is “Quantum Computing and Quantum Information” by Michael Nielsen and Isaac Chuang.

Another essential component of quantum is measurement. This transforms the quantum state into a probability distribution we can sample from. Note that these operations are not reversible. This is because after measurement, a quantum state collapses:

class M(Transform):
    "Turn a quantum statevector into a probability distribution"
    def encodes(self, x): return np.abs(x)**2
    def decodes(self, x): return NotImplementedError("No inverse exists for absolute value.")

class Samp(Transform):
    "Sample from a probability distribution"
    def encodes(self, x): return format(np.random.choice(len(x), p=x), f'0{int(np.log2(len(x)))}b')
    def decodes(self, x): return NotImplementedError("Sampling is not reversible.")

m = M()
samp = Samp()

# Sampling from zero state (|0>)
zero_state = [1+0j, 0+0j]
mzs = m(zero_state) # Transforms [1+0j, 0+0j] -> [1, 0]

assert_almost_equal(mzs, np.array([1+0j,0+0j]))                     
assert samp(mzs) == '0'

# Sampling from equal superposition
equal_superposition = [0.707, 0.707]
mes = m(equal_superposition) # Transforms [0.707, 0.707] -> [0.5, 0.5] (A coin flip. i.e. Bernoulli distribution))

mes = mes / mes.sum()  # not in original blog but needed as otherwises numers dont add to 1 because of rounding
assert_almost_equal(mes, np.array([0.5,0.5]))
assert samp(mes) in '01' # 0 or 1 with equal probability

# Sampling from complex superposition
complex_superposition = [0.5+0.5j, 0.5-0.5j]
mcs = m(complex_superposition) # Transforms [0.5+0.5j, 0.5-0.5j] -> [0.5, 0.5] (A coin flip. i.e. Bernoulli distribution)
assert_almost_equal(mcs, np.array([0.5,0.5]))

assert samp(mcs) in '01'  # Result is 0 or 1 with equal probability

We can now build a full quantum circuit using Pipeline. Note that a quantum pipeline is only reversible if it does not include measurement or sampling:

qc = Pipeline([X(), H(), I(), M(), Samp()])
# X transforms [1, 0] -> [0, 1]
# H transforms [0, 1] -> [0.707+0j, -0.707+0j]
# I transforms [0.707+0j, -0.707+0j] -> [0.707+0j, -0.707+0j]
# M transforms [0.707+0j, -0.707+0j] -> [0.5, 0.5]
# Samp samples from random distribution [0.5, 0.5]
assert qc(zero_state) in '01' # 0 or 1 with equal probability

Multi Qubit

Even though we only discussed single qubit cases we can the potential of Transform for quantum. It gets even more powerful if we start working with multiple qubits. The representation of a quantum state on a classical computer grows exponentially with each qubit, because each qubit can be entangled with others. We therefore need a matrix of \(2^n\) x \(2^n\) to represent a transformation of \(n\) qubits. The statevector for \(n\) qubits contains \(2^n\) complex numbers. Single qubit gates can be combined through the Tensor (Kronecker) product, which we handle in Concat:

class Concat(Transform):
    "Combine single qubit gates into a multi-qubit gate"
    def __init__(self, gates): self.gates = gates
    # Concatenate 2 or more gates
    def encodes(self, x): return x @ np.kron(*[g.gate for g in self.gates])
    # Reverse propagation for all gates
    def decodes(self, x):
        for g in reversed(self.gates): x = x @ np.kron(g.gate.conj().T, np.eye(len(x) // g.gate.shape[0]))
        return x

By concatenating single qubit gates we can construct multi-qubit circuits, while keeping the code extremely simple. However, some gates are fundamentally multi-qubit and cannot be constructed from single qubits. One example is the Controlled NOT (CNOT) gate, which flips the 2nd qubit from \(|0\rangle\) to \(|1\rangle\) based on the value of the 1st qubit. When we combine Hadamard on the 1st qubit and a CNOT gate we get the well known Bell state. This is a classic example of fully entangling qubits where we obtain \(00\) or \(11\) with equal probability.

class CNOT(_Q):
    "Controlled NOT gate"
    def __init__(self): self.gate = np.array([[1, 0, 0, 0], 
                                              [0, 1, 0, 0], 
                                              [0, 0, 0, 1], 
                                              [0, 0, 1, 0]])
                     
two_qubit_zero_state = np.array([1+0j, 0+0j, 0+0j, 0+0j]) # |00>
qc = Pipeline([Concat([H(), I()]), CNOT(), M(), Samp()])
# Concat([H(), I()]) transforms [1, 0, 0, 0] -> [0.707, 0, 0.707, 0]
# CNOT() transforms [0.707, 0, 0.707, 0] -> [0.707, 0, 0, 0.707]
# M() transforms [0.707, 0, 0, 0.707] -> [0.5, 0, 0, 0.5]
# Samp() samples from [0.5, 0, 0, 0.5] (50% chance at 00 and 50% chance at 3, which is 11 in binary)
assert qc(two_qubit_zero_state) in ('00','11') # 00 or 11 with equal probability

These techniques can be used to simulate and analyze more complicated multi-qubit circuits. The nice thing about building quantum circuits like this is that we can analyze every step and get a good understanding of what is happening. It also allows us to precisely explore techniques like quantum error correction. However, for large scale quantum circuits the matrices are huge and real quantum computers are needed to do the computation. Real quantum computers directly leverage properties of quantum mechanics like entanglement, superposition and interference. These are things that a classical computer can simulate, but cannot natively perform like a real quantum computer. Exploiting quantum properties can result in potentially exponential speedups. For more information on where quantum computers excel check out Ronald de Wolf’s great paper on “The Potential Impact of Quantum Computers on Society (2017). This paper has stood the test of time even though new breakthroughs have been achieved.m

Closing

If you made it all the way to the end, congratulations! I salute you! 🫡 Having an understanding of both advanced Python and quantum computing is a rare skillset. If you are interested in this intersection you might be interested in exploring Quantum Machine Learning. This field combines the optimization we often see in data science and machine learning with quantum computing. If this piques your interest, one of the best textbooks around is “Machine Learning with Quantum Computers by Maria Schuld and Francesco Petruccione”.

Learning resources

If you are interested in learning more about quantum, I recommend the following resources: