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]))
../../_images/f27d68de65c7984fae3dd6904834a8bd3a6d6e8387d7e79429bb9dab22acb9e9.png

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')
../../_images/8da103114e7da1ac71c1d5d62785484f0a8f7f5ce7facde603ab6394f711652a.png
createCircuit(np.array((np.pi,np.pi,1,2)),G,2).draw(output='mpl')
../../_images/f6a04d8f4e7bdc0f116d617675aa8d532cdc6a30730f4fa6ebcb9d0cd6a54efc.png

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)
../../_images/7b9a2bf88434da7efa5606e0c4674bddf997b230de35dd5f10513729810cea35.png

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)
../../_images/550ee1e86c35895d5243b0a25152dc63e727c76ba97cf6fd71b657fad00f4d6d.png

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())
../../_images/8ef6fe6229175ef69eafa58f22a23c2185a1784a1e5d06ec880a75d1fdf5f246.png

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())
../../_images/ce215e9c5739b08e9422e19e5c5c034db3a4b516afcbf510ccc5251f1a993eb8.png

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$')
../../_images/2e53de1f4e3fdce5a0ba6bb63786dce6f0dd9eea96a3034ddc8210540fdaff7b.png
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$')
../../_images/c4ccedc8d35537228e003c27779e48b8d155725c41469e51a909388aa53c768f.png

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())
../../_images/7ca5375e5c871f4eefc345755b80328bd0e560d00ec61d4c56da79acf0781f90.png

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)
../../_images/8b7d0f35b6fc21afadb838bcd92f145d3f15e9deaac02d93fee7290b6d2ef465.png

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')
../../_images/6076d2b592107ea1f3a62cca29ec4f7bf516e6a637e3b685960f4230c9318f70.png

how does the circuit look on a real device?

plot_gate_map(FakeVigo())
../../_images/9cdfefdcdfa3fc0dc22288a97d02eb91a8ddf40f26c68cc9512e09d6d5913833.png
transpile(createCircuit(np.array((0.4,0.1)),G,1),backend=sim_vigo).draw(output='mpl')
../../_images/1c38c712d5332813c3eebd08c50f3f4e30165cf48aed49a811eb2a23b50013b5.png

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$')
../../_images/e0816804e632bdf0010d8e1c49001efee0b3468886080991a1f767be9dd72d94.png

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
../../_images/6fdde514a77ccbb3ace6e6ca0a4cd29a3803293fab74b1e0519e171addaa07a8.png
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.
../../_images/f0833f38c0fe75a4cd2ab7af005671c25a1f150ffc44df523d1ef003d481a776.png

Where to go next? Try increasing the depth…