Import Libraries
from numpy import *
from scipy.linalg import hadamard
import matplotlib.pyplot as plt
import time
from tabulate import tabulate
from scipy.optimize import curve_fit
Handle and arbitrary number of qubits¶
Program Grover's search for an arbitrary number of qubits N. Check that the optimum number of Grover iterations is $\frac{\pi}{4}\sqrt{2^{N}}$, for $N>3$. Then increase $N$ as much as possible.
The Hadamard Gate
A Hadamard gate is described by the matrix
For example, with $N=4$ qubits, the operation that applies a Hadamard gate to qubit 1 is
nn = 4
N = 2**nn
I = eye(2)
H = hadamard(2)
factor = 1/sqrt(2)
def had(i):
if i == nn:
product = H
else:
product = I
for j in range(nn-1,0,-1):
if j!=i:
product = kron(I,product)
else:
product = kron(H,product)
return array(product).astype(int)
H1, H2, H3 = had(1), had(2), had(3)
print(f"The Hadamard gate that applied to qubit 1 is ~{round(factor,4)} multiplied by: \n{H1}")
Grover's Algorithm
This algorthm requires the oracle $\hat{O}$ and the $\hat{J}$ matricies. For $N=4$, the product $\hat{H}^{(4)}\cdot\hat{H}^{(3)}\cdot \hat{H}^{(2)}\cdot \hat{H}^{(1)}$ is an operator that applies a Hadamard gate to qubits 1, 2, 3, and 4. This operator is equvalent to the 'hadamard' function from 'scipy.linalg' for $2^{N}$, where $N=4$, multiplied by the factor $\frac{1}{4}$. Let's program Grover's search for an arbitrary number of qubits $N$.
def opO(state,N):
O = identity(N,dtype=int)
O[state,state] = -1
return O
def opJ(N):
J = identity(N,dtype=int)
J[0,0] = -1
return J
def groverN(nn,N,state,rep,question):
qb = zeros(N)
qb[state] = 1
qb = array([qb]).T
oracle = opO(question,N)
J = opJ(N)
factor = (1/(sqrt(2)))**nn
H = factor*hadamard(N)
qb = H @ qb
for i in range(rep):
qb = H @ J @ H @ oracle @ qb
return qb
Let's test Grover's Algorithm for $N=4$ qubits and $3$ iterations (approximately the optimum number) on the state $\lvert0000\rangle$, with the correct question $\lvert0000\rangle$.
print(f"The optimum # of repititions is {round((pi/4)*sqrt(N),2)}"
f" for ~90% return on your state\n")
state = 0
rep = 3
g = groverN(nn,N,state,rep,state)
g = g**2
x = linspace(0,N-1,N)
y = g.T
y = y[0]
plt.bar(x,y,color='purple')
plt.title("Result Probabilities for each Quantum State")
plt.ylim(0,1)
plt.show()
reps = 19
arr = zeros(reps)
rarr = linspace(1,reps,reps)
for i in range(1,reps):
g = groverN(nn,N,state,i,state)
arr[i] = g[0]**2
def g(x,A,w,s,C):
return A*sin(w*x+s)+C
guesses = (1,1,1,1)
(A,w,s,C), cc = curve_fit(g,rarr,arr,p0=guesses)
xmod = linspace(rarr[0],rarr[-1],100)
ymod = g(xmod,A,w,s,C)
plt.plot(rarr,arr,'ko')
plt.plot(xmod,ymod,color='purple')
plt.xlim(1,reps)
plt.title("Quantum State Probabilities vs. Number of Iterations")
plt.show()
The measured result is $\lvert0000\rangle$ more than 90% of the time, as expected, and the return rate on the quantum state is cyclic. Let's see how long it takes to compute Grover's Algorithm.
arr = list()
m = 12
qarr = linspace(2,m,m-1).astype(int)
times = list()
for r in range(len(qarr)):
time_start = time.perf_counter()
g = groverN(qarr[r],2**qarr[r],state,rep,state)
time_elapsed = (time.perf_counter() - time_start)
arr.append([qarr[r],time_elapsed])
times.append(time_elapsed)
print(tabulate(arr,headers=["N-qubits","Time (s)"]))
Now we can plot this data and create a best fit!
def f(x,A,p):
return A*exp(x*p)
guesses = (1,1)
(A,p), cc = curve_fit(f,qarr,array(times),p0=guesses)
xmod = linspace(qarr[0],qarr[-1],100)
ymod = f(xmod,A,p)
plt.plot(qarr,array(times),'ko')
plt.plot(xmod,ymod,color='purple')
plt.xlabel("N-qubits")
plt.ylabel("Grover Execution time (s)")
plt.show()
The expoential fit can be used to predict how long it should take the algorithm to do N-qubit computations (given it has enough RAM). For 4GB RAM and 16 bytes per complex amplitude, the largest possible state is 27-qubit.
c = 27
v = 21
predict = list()
prx = linspace(v,c,c-v-1).astype(int)
pry = f(prx,A,p)/(60**2*24*365*100)
for g in range(len(prx)): predict.append([prx[g],pry[g]])
print(tabulate(predict,headers=["N-qubits","Time (centuries)"]))