{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Quantum circuits in Python using nothing but Numpy\n", "Joris Kattemölle\n", "\n", "There are fancy packages that can (classically) run small quantum circuits in Python, such as Cirq, Qiskit, ProjectQ and even QuTiP. It is, of course, possible to do this without using any of these packages. I think this is much more enlightening. Here, I show how it can be done in a minimalistic, but high-performance way. Though simple, the code below is actually competitive with aforementioned packages in terms of speed. \n", "\n", "In the following, I'll assume you have a good understanding of quantum mechanics, and a basic working knowledge of Python and Numpy. I'll use expressions like $a$ =a, where $a$ is in standard mathematical notation, and a is a representation of that same thing in Python. The following is going to involve a lot of explanation, but in the end only a few lines of code, and overview of which can be found at the end of this page. This page can be downloaded as a [Jupyter notebook](http://www.kattemolle.com/other/QCinPY.ipynb) or as a [pdf](http://www.kattemolle.com/other/QCinPY.pdf).\n", "\n", "\n", "Contents:
\n", "[States](#States)
\n", "[Single qubit operators](#Single-qubit-operators)
\n", "[The register class](#The-register-class)
\n", "[Two-qubit operators](#Two-qubit-operators)
\n", "[Measurement](#Measurement)
\n", "[Conclusion](#Conclusion)
\n", "[Code](#Code)
\n", "\n", "## States\n", "For simplicity, let us consider $n=4$ spin-1/2 particles, or qubits. Any state $|\\psi \\rangle$ can be expanded as\n", "\n", "$|\\psi\\rangle = \\sum_{ijkl}\\psi_{ijkl}|i\\rangle|j\\rangle|k\\rangle|l\\rangle.$\n", "\n", "To store this state in Python we just store the values $\\{\\psi_{ijkl}\\}$ in a Numpy array psi in such a way that $\\psi_{ijkl}=$ psi[i,j,k,l].\n", "\n", "Note that the shape of psi is (2,2,2,2). Or in other words, it has four indices, each of which can take two values. \n", "\n", "### Example 1\n", "Let $|\\psi \\rangle = (|0000\\rangle + |1111\\rangle)/\\sqrt{2}$, or in component notation: $\\psi_{0000}=\\psi_{1111}=1/\\sqrt{2}$ with all other components vanishing. This state is stored as" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "psi=np.zeros((2,2,2,2)) # Create an array of zeros with the right shape. \n", "psi[0,0,0,0]=psi[1,1,1,1]=1/np.sqrt(2) # Set the right entries to 1/sqrt(2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Single qubit operators\n", "The Hadamard operator is given by $H=\\frac{1}{\\sqrt{2}}\\left(\\begin{array}{rr}1&1\\\\ 1&-1\\end{array}\\right)$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "code_folding": [] }, "outputs": [], "source": [ "H_matrix=1/np.sqrt(2)*np.array([[1, 1],\n", " [1,-1]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Say we want to act with this operator on the 0th qubit, $|\\psi'\\rangle=X\\otimes I \\otimes I \\otimes I \\ |\\psi\\rangle$. (We count qubits from left to right, starting from 0. The same holds for indices.) In component notation, this means we want to create an array with components\n", "\n", "psi_prime[i,j,k,l] $=\\psi'_{ijkl}=\\sum_{i'}X_{ii'}\\psi_{i'jkl}.$\n", "\n", "The can be achieved in a high-performance way by using [np.tensordot](https://numpy.org/doc/stable/reference/generated/numpy.tensordot.html?highlight=tensordot#numpy.tensordot), \n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "psi_prime=np.tensordot(H_matrix,psi,(1,0)) # Contract the 1st index of H_matrix with the 0th index of psi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> Remark 1. Of course, we could also represent states as long, one-dimensional arrays, in the way you were probably thought in your quantum mechanics 101 course. In this case, the operator $X\\otimes I \\otimes I \\otimes I \\ |\\psi\\rangle$ would be constructed by taking four Kronecker products. However, this approach is not very efficient, nor very elegant, for systems with $n\\gg4$ qubits. This is because merely constructing the operator $X\\otimes I \\otimes \\ldots \\otimes I$ requires you to create an (albeit sparse) $2^n\\times2^n$ matrix. \n", "\n", "There is an important caveat if we want to apply $H$ to the 1st qubit (or 2nd or 3rd). Namely, by construction, np.tensordot(H_matrix,psi,(1,1)), which contracts the 1st index of H_matrix with the 1st index of psi, produces an array with entries psi_prime[j,i,k,l] $=\\sum_{j'}H_{jj'}\\psi_{ij'kl}$. Note the order of the indices: np.tensordot contracts indices, but leaves the free indices in the order reading from left to right. This is not what we're used to in quantum mechanics. We don't want an operator to change the order of the qubits. To fix this, we have to move the indices (or axes) back in the right places with [np.moveaxis](https://numpy.org/doc/stable/reference/generated/numpy.moveaxis.html?highlight=moveaxis). So, the correct way to construct psi_prime is by" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "psi_prime=np.tensordot(H_matrix,psi,(1,1)) # Contract the 1st index of H_matrix with the 1st index of psi \n", "psi_prime=np.moveaxis(psi_prime,0,1) # Put axes in the right place by putting axis 0 at position 1, keeping the other axis in order" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The register class\n", "\n", "To make the process a bit smoother, let's define the register as a class, and define the Hadamard gate as a function that acts on this class. The class only contains the state of the register, psi, and the number of qubits, n. The state of the register is initialized to the $n$-qubit state $|0,0,\\ldots,0\\rangle$." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [], "source": [ "class Reg: \n", " def __init__(self,n):\n", " self.n=n\n", " self.psi=np.zeros((2,)*n) # make array of zeros with right shape\n", " self.psi[(0,)*n]=1 # put psi[0,0,...,0] to 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now act with the Hadamard operator on the $i$th qubit by invoking the function" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def H(i,reg): # Alter the array psi into the array where H has acted on the ith qubit\n", " reg.psi=np.tensordot(H_matrix,reg.psi,(1,i)) \n", " reg.psi=np.moveaxis(reg.psi,0,i)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Examples" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n" ] } ], "source": [ "reg=Reg(4)\n", "print(reg.psi.flatten())" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.70710678 0. 0. 0. 0. 0.\n", " 0. 0. 0.70710678 0. 0. 0.\n", " 0. 0. 0. 0. ]\n" ] } ], "source": [ "reg=Reg(4)\n", "H(0,reg)\n", "print(reg.psi.flatten())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-qubit operators\n", "We now extend the previous to two-qubit operations. As an example, take the $\\mathrm{CNOT}$ operator, defined by \n", "\n", "$\\mathrm{CNOT}|00\\rangle=|00\\rangle,\\quad \\mathrm{CNOT}|01\\rangle=|01\\rangle, \\quad \\mathrm{CNOT}|10\\rangle=|11\\rangle, \\quad \\mathrm{CNOT}|11\\rangle=|10\\rangle$.\n", "\n", "In component notation, \n", "\n", "$\\mathrm{CNOT}_{0000}=\\mathrm{CNOT}_{0101}=\\mathrm{CNOT}_{1011}=\\mathrm{CNOT}_{1110}=1$, with other components vanishing. \n", "\n", "Here the 0th qubit is called the control qubit and the 1st qubit the target qubit. Say we want to implement $\\mathrm{CNOT}\\otimes I \\otimes I |\\psi \\rangle$, or, in component notation, construct the array with components psi_prime[i,j,k,l] $=\\psi'_{ijkl}=\\sum_{k'l'}\\mathrm{CNOT}_{ij k' l'}\\psi_{k'l'kl}$.\n", "\n", "Now, before we proceed, I want to note that more often than not, the CNOT operator is depicted as a $2\\times2$ array, \n", "\n", "$\\mathrm{CNOT}=\\left(\\begin{array}{cccc}1&0&0&0\\\\ 0&1&0&0\\\\ 0&0&0&1 \\\\ 0&0&1&0 \\end{array}\\right).$ " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "CNOT_matrix=np.array([[1,0,0,0],\n", " [0,1,0,0],\n", " [0,0,0,1],\n", " [0,0,1,0]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is the form you obtain following the method of Remark 1. However, for our purposes, this tensor does not have the right shape. It has only two axes, or indices, whereas we need four in the component formulas above. Fortunately, we can easily get it in the right shape by using [np.reshape](https://numpy.org/doc/stable/reference/generated/numpy.reshape.html?highlight=reshape#numpy.reshape)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "CNOT_tensor=np.reshape(CNOT_matrix, (2,2,2,2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again taking into account the right order of the indices, we can implement this operator by applying the function" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def CNOT(control, target, reg):\n", " # Contract 2nd index of CNOT_tensor with control index, and 3rd index of CNOT_tensor with target index.\n", " reg.psi=np.tensordot(CNOT_tensor, reg.psi, ((2,3),(control, target))) \n", " # Put axes back in the right place\n", " reg.psi=np.moveaxis(reg.psi,(0,1),(control,target)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example 2\n", "When acting on psi $=|00\\ldots\\rangle$, the following function constructs the (generalization of) the state from Example 1." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[0.70710678 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0. 0. 0.\n", " 0. 0. 0. 0.70710678]\n" ] } ], "source": [ "def generate_GHZ(reg): \n", " H(0,reg)\n", " for i in range(reg.n-1):\n", " CNOT(i,i+1,reg)\n", " \n", "# Usage\n", "reg=Reg(4)\n", "generate_GHZ(reg)\n", "print(reg.psi.flatten())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Measurement" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A measurement of the $i$th qubit in the computational basis changes the state $|\\psi\\rangle$ into $|\\psi'\\rangle= \\frac{|j\\rangle_i\\langle j|_i |\\psi\\rangle}{\\lVert|j\\rangle_i\\langle j|_i |\\psi\\rangle \\rVert}$ with probability $p(j)=\\lVert\\,|j\\rangle_i\\langle j|_i |\\psi \\rangle\\rVert^2$. (Here e.g. $\\langle j |_0=\\langle j| \\otimes I \\otimes I \\otimes I\\ \\$). \n", "\n", "First, let us implement $|j\\rangle_i \\langle j|_i$ for the two different vales of $j$. Given a value for $j$, the operator $|j\\rangle_i \\langle j|_i$ can be written as a matrix, and we may implement this matrix in the same way we've implemented H_matrix." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "projectors=[np.array([[1,0],[0,0]]), np.array([[0,0],[0,1]]) ] # list containing the projectors |0><0| and |1><1|\n", "\n", "def project(i,j,reg): # RETURN state with ith qubit of reg projected onto |j>\n", " projected=np.tensordot(projectors[j],reg.psi,(1,i))\n", " return np.moveaxis(projected,0,i)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now to emulate measurement, we need to act with the right projector with the right probability:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "from scipy.linalg import norm \n", "\n", "def measure(i,reg): # Change reg.psi into post-measurement state w/ correct probability. Return measurement value as int\n", " projected=project(i,0,reg) \n", " norm_projected=norm(projected.flatten()) \n", " if np.random.random() 