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…