Use case: Optimization
We are going to implement QAOA for the MaxCut problem
we start by importing qiskit to simulate circuits and some other basic libraries
from qiskit import *
import numpy as np
import pylab as pl
from mpl_toolkits.axes_grid1 import make_axes_locatable
import networkx as nx
from qiskit.visualization import *
networkx can be used to handle graphs
V = np.arange(0,2,1)
E =[(0,1,1.0)]
G = nx.Graph()
G.add_nodes_from(V)
G.add_weighted_edges_from(E)
nx.draw_networkx(G)
np.array(list(G.nodes)), V
(array([0, 1]), array([0, 1]))

We define a function that creates our circuit
def createCircuit(x,G,depth):
V = list(G.nodes)
num_V = len(V)
q = QuantumRegister(num_V)
c = ClassicalRegister(num_V)
circ = QuantumCircuit(q,c)
#uniform superposition
circ.h(range(num_V))
circ.barrier()
for d in range(depth):
gamma=x[2*d]
beta=x[2*d+1]
#go through all edges and add Rzz gate
for edge in G.edges():
i=int(edge[0])
j=int(edge[1])
w = G[i][j]['weight']
circ.cx(q[i],q[j])
circ.rz(w*gamma,q[j])
circ.cx(q[i],q[j])
circ.barrier()
#add the mixer
circ.rx(2*beta,range(num_V))
circ.barrier()
circ.measure(q,c)
return circ
draw example circuits
createCircuit(np.array((np.pi,np.pi)),G,1).draw(output='mpl')

createCircuit(np.array((np.pi,np.pi,1,2)),G,2).draw(output='mpl')

in order to evaluate a solution we define a function that gives us the cost
def cost(x,G):
C=0
for edge in G.edges():
i = int(edge[0])
j = int(edge[1])
w = G[i][j]['weight']
C = C + w/2*(1-(2*x[i]-1)*(2*x[j]-1))
return C
brute force function that lists all 2^n possibilities and finds the best solutions (use with caution)
def listcosts(G):
costs={}
maximum=0
solutions=[]
V = list(G.nodes)
num_V = len(V)
for i in range(2**num_V):
binstring="{0:b}".format(i).zfill(num_V)
y=[int(i) for i in binstring]
costs[binstring]=cost(y,G)
maximum = max(maximum,costs[binstring])
for key in costs:
if costs[key]==maximum:
solutions.append(key)
return costs, maximum, solutions
listcosts(G)
({'00': 0.0, '01': 1.0, '10': 1.0, '11': 0.0}, 1.0, ['01', '10'])
the result of a circuit contains a dictionary containing bitstrings together with how many times they have occured
we define a function that gives us (an approximation of) the expectationvalue based on this
def expectationValue(data,G):
res=data.result().results
E=[]
V = list(G.nodes)
num_V = G.number_of_nodes()
for result in res:
n_shots = result.shots
counts = result.data.counts
e = 0
for hexkey in list(counts.keys()):
count = counts[hexkey]
binstring = "{0:b}".format(int(hexkey,0)).zfill(num_V)
y=[int(i) for i in binstring]
e += cost(y,G)*count/n_shots
E.append(-e)
return E
we will use an ideal simulator
ideal_sim = Aer.get_backend('qasm_simulator')
to get the energy/training/cost landscape for depth p=1, we sample the region \([0,\frac{\pi}{2}]^2\)
circuits=[]
n=16
for gamma in np.linspace(0,np.pi,n):
for beta in np.linspace(0,np.pi,n):
circuits.append(createCircuit(np.array((gamma,beta)),G,1))
job_sim = execute(circuits, ideal_sim, shots=1024*2*2*2)
val=expectationValue(job_sim,G)
E=np.array(val).reshape(n,n)
f = pl.figure(figsize=(6, 6), dpi= 80, facecolor='w', edgecolor='k');
_=pl.xlabel(r'$\gamma$')
_=pl.ylabel(r'$\beta$')
ax = pl.gca()
im = ax.imshow(E,interpolation='bicubic',origin='lower',extent=[0,np.pi,0,np.pi])
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
_=pl.colorbar(im, cax=cax)

we import minimizers from scipy to do local minimization
from scipy import optimize as opt
the minimize funciton needs a function to evalution, which we create next
def getval(x, backend):
circ=createCircuit(x,G,1)
tcirc = transpile(circ, backend)
j = execute(tcirc, backend, shots=1024*2*2*2)
val=expectationValue(j,G)
return val[0]
out=opt.minimize(getval, x0=(1, 1), method='Nelder-Mead',\
args=(ideal_sim),\
options={'xatol': 1e-2, 'fatol': 1e-2, 'disp': True})
Optimization terminated successfully.
Current function value: -1.000000
Iterations: 18
Function evaluations: 36
f = pl.figure(figsize=(6, 6), dpi= 80, facecolor='w', edgecolor='k');
_=pl.xlabel(r'$\beta$')
_=pl.ylabel(r'$\gamma$')
ax = pl.gca()
im = ax.imshow(E,interpolation='bicubic',origin='lower',extent=[0,np.pi,0,np.pi])
_=pl.plot(out.x[1],out.x[0],'xr')
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)
_=pl.colorbar(im, cax=cax)

let’s plot the histogram
j = execute(createCircuit(out.x,G,1), ideal_sim, shots=1024*2*2*2)
plot_histogram(j.result().get_counts())

next, let’s see the effect of noise
we use a backend that emulates a real device “Vigo”
from qiskit.providers.aer import AerSimulator
from qiskit.test.mock import FakeVigo
device_backend = FakeVigo()
sim_vigo = AerSimulator.from_backend(device_backend)
plot_error_map(FakeVigo())

again, we sample the energy/training/cost landscape for depth p=1
noisy_circuits=[]
n=16
for gamma in np.linspace(0,np.pi,n):
for beta in np.linspace(0,np.pi,n):
circ=createCircuit(np.array((gamma,beta)),G,1)
#we need to transpile the circuit for the backend
tcirc=transpile(circ,sim_vigo)
noisy_circuits.append(tcirc)
job_sim_noise = execute(noisy_circuits, sim_vigo, shots=1024*2*2*2)
val_noise=expectationValue(job_sim_noise,G)
E_noise=np.array(val_noise).reshape(n,n)
we can now compare the landscape from the ideal and noisy simulation
(the circuit consist of only 2 qubits and 2 cnot gates)
f, axarr = pl.subplots(1,2,figsize=(10, 10), dpi= 80, facecolor='w', edgecolor='k')
vmin=min(np.min(E),np.min(E_noise))
vmax=min(np.max(E),np.max(E_noise))
_=axarr[0].imshow(E,vmin=vmin,vmax=vmax,origin='lower',extent=[0,np.pi,0,np.pi])
_=axarr[1].imshow(E_noise,vmin=vmin,vmax=vmax,origin='lower',extent=[0,np.pi,0,np.pi])
for i in range(2):
_=axarr[i].set_xlabel(r'$\beta$')
_=axarr[i].set_ylabel(r'$\gamma$')

out_noise=opt.minimize(getval, x0=(1,1),\
args=(sim_vigo), method='Nelder-Mead',\
options={'xatol': 1e-2, 'fatol': 1e-2, 'disp': True})
Optimization terminated successfully.
Current function value: -0.894897
Iterations: 20
Function evaluations: 48
f, axarr = pl.subplots(1,2,figsize=(10, 10), dpi= 80, facecolor='w', edgecolor='k')
vmin=min(np.min(E),np.min(E_noise))
vmax=min(np.max(E),np.max(E_noise))
_=axarr[0].imshow(E,vmin=vmin,vmax=vmax,origin='lower',extent=[0,np.pi,0,np.pi],interpolation='bicubic')
_=axarr[1].imshow(E_noise,vmin=vmin,vmax=vmax,origin='lower',extent=[0,np.pi,0,np.pi],interpolation='bicubic')
_=axarr[0].plot(out.x[1],out.x[0],'xr')
_=axarr[1].plot(out_noise.x[1],out_noise.x[0],'ok')
for i in range(2):
_=axarr[i].set_xlabel(r'$\beta$')
_=axarr[i].set_ylabel(r'$\gamma$')

let’s plot the histogram. How did the noise effect it?
circ=createCircuit(out_noise.x,G,1)
tcirc=transpile(circ,sim_vigo)
j = execute(tcirc, sim_vigo, shots=1024*2*2*2)
plot_histogram(j.result().get_counts())

finally, let’s create a larger circuit
V = np.arange(0,5,1)
E =[(0,1,1.0),(0,2,1.0),(2,3,1.0),(3,1,1.0),(3,4,1.0),(4,2,1.0)]
G = nx.Graph()
G.add_nodes_from(V)
G.add_weighted_edges_from(E)
nx.draw_networkx(G)

what are the solutions?
l,m,maxcost=listcosts(G)
print(m,maxcost)
5.0 ['01100', '01101', '10010', '10011']
how does the circuit look?
createCircuit(np.array((0.4,0.1)),G,1).draw(output='mpl')

how does the circuit look on a real device?
plot_gate_map(FakeVigo())

transpile(createCircuit(np.array((0.4,0.1)),G,1),backend=sim_vigo).draw(output='mpl')

create energy landscapes for ideal and noisy simulations
circuits=[]
n=16
for gamma in np.linspace(0,np.pi,n):
for beta in np.linspace(0,np.pi,n):
circ=createCircuit(np.array((gamma,beta)),G,1)
circuits.append(circ)
job_sim = execute(circuits, ideal_sim, shots=1024*2*2*2)
val=expectationValue(job_sim,G)
E=np.array(val).reshape(n,n)
noisy_circuits=[]
n=16
for gamma in np.linspace(0,np.pi,n):
for beta in np.linspace(0,np.pi,n):
circ=createCircuit(np.array((gamma,beta)),G,1)
#we need to transpile the circuit for the backend
tcirc=transpile(circ,sim_vigo)
noisy_circuits.append(tcirc)
job_sim_noisy = execute(noisy_circuits, sim_vigo, shots=1024*2*2*2)
val_noisy=expectationValue(job_sim_noisy,G)
E_noise=np.array(val_noisy).reshape(n,n)
a comparison of the energy landscapes
f, axarr = pl.subplots(1,2,figsize=(10, 10), dpi= 80, facecolor='w', edgecolor='k')
vmin=min(np.min(E),np.min(E_noise))
vmax=min(np.max(E),np.max(E_noise))
_=axarr[0].imshow(E,vmin=vmin,vmax=vmax,origin='lower',extent=[0,np.pi,0,np.pi],interpolation='bicubic')
_=axarr[1].imshow(E_noise,vmin=vmin,vmax=vmax,origin='lower',extent=[0,np.pi,0,np.pi],interpolation='bicubic')
# _=axarr[0].plot(out.x[1],out.x[0],'xr')
# _=axarr[1].plot(out_noise.x[1],out_noise.x[0],'ok')
for i in range(2):
_=axarr[i].set_xlabel(r'$\beta$')
_=axarr[i].set_ylabel(r'$\gamma$')

histograms of the best solution at depth p=1
out=opt.minimize(getval, x0=(1,1),\
args=(ideal_sim), method='Nelder-Mead',\
options={'xatol': 1e-2, 'fatol': 1e-2, 'disp': True})
circ=createCircuit(out.x,G,1)
j = execute(circ, ideal_sim, shots=1024*2*2*2)
plot_histogram(j.result().get_counts())
Optimization terminated successfully.
Current function value: -3.890503
Iterations: 14
Function evaluations: 28

out_noise=opt.minimize(getval, x0=(1,1),\
args=(sim_vigo), method='Nelder-Mead',\
options={'xatol': 1e-2, 'fatol': 1e-2, 'disp': True})
circ=createCircuit(out_noise.x,G,1)
tcirc=transpile(circ,sim_vigo)
j = execute(tcirc, sim_vigo, shots=1024*2*2*2)
plot_histogram(j.result().get_counts())
Warning: Maximum number of function evaluations has been exceeded.

Where to go next? Try increasing the depth…