This document explains how to solve the knapsack problem by QAOA. There are 3 steps; prepare set of solutions, evaluate objective value of set, check if weight of solution exceeds the capacity. The circuit proceeding these 3 steps is propsed.
`import numpy as np
import qiskit
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import execute, IBMQ, BasicAer, Aer
from qiskit.tools import visualization
from qiskit.quantum_info.operators import Operator, Pauli
from qiskit.tools.visualization import plot_histogram
from qiskit.tools.monitor import job_monitor
import qiskit.ignis.verification as ver
import matplotlib.pyplot as plt
from itertools import product, combinations
N = 5
C = 6
cost = [3,5,9,7,4]
weight = [2,3,5,4,1]
beta = 0.1 np.pi
gamma = 0.1 np.pi
parameter
B=1
A = B*(max(cost)+100)
function
def X(i):
x = np.eye(1)
for k in range(N):
if k == i:
x = np.kron(x, np.array([[0,1],[1,0]]))
else:
x = np.kron(x, np.eye(2))
return x
def Z(i, coeff):
z = np.eye(1)
for k in range(N):
if k == i:
z = np.kron(z, np.array([[coeff,0],[0,-1*coeff]]))
else:
z = np.kron(z, np.eye(2))
return z
def Unitary(A,k):
vec = np.linalg.eig(A)[1]
val = np.linalg.eig(A)[0]
d_list = np.array([np.e*(ki*1j) for i in val])
diag = np.diag(d_list)
a=np.dot(vec,diag)
b=np.dot(a,np.linalg.inv(vec))
return b
def Cost(x):
return sum([int(x[i])*cost[N-1-i] for i in range(N)])
def Weight(x):
return sum([int(x[i])*weight[N-1-i] for i in range(N)])
term = dict()
term["null"] = 0
for i in range(N):
term[qr_N[i]] = 0
for n in range(C):
term[qr_C[n]] = 0
for i in range(N):
for n in range(C):
term[qr_N[i], qr_C[n]] = 0
for i,j in combinations(list(range(N)),2):
term[qr_N[i],qr_N[j]] = 0
for m,n in combinations(list(range(C)),2):
term[qr_C[m],qr_C[n]] = 0
for i in range(N):
term[qr_N[i]] += -1/2(A(weight[i]*2)-Bcost[i])
term["null"] += 1/2(A(weight[i]*2)-Bcost[i])
for n in range(C):
term[qr_C[n]] += -1/2(A(n2)-1)
term["null"] += 1/2(A(n2)-1)
for (i,j) in combinations(list(range(N)),2):
term["null"] += 2Aweight[i]weight[j]1/4
term[qr_N[i]] += 2Aweight[i]weight[j](-1/4)
term[qr_N[j]] += 2Aweight[i]weight[j](-1/4)
term[qr_N[i], qr_N[j]] += 2Aweight[i]weight[j](1/4)
for (m,n) in combinations(list(range(C)),2):
term["null"] += 2Amn1/4
term[qr_C[m]] += 2Amn(-1/4)
term[qr_C[n]] += 2Amn(-1/4)
term[qr_C[m], qr_C[n]] += 2Amn(1/4)
for i in range(N):
for n in range(C):
term["null"] += -2Aweight[i]n1/4
term[qr_N[i]] += -2Aweight[i]n(-1/4)
term[qr_C[n]] += -2Aweight[i]n(-1/4)
term[qr_N[i], qr_C[n]] += -2Aweight[i]n(1/4)
for t in term:
if type(t) == qiskit.circuit.quantumregister.Qubit:
circuit.p(gammaterm[t],t)
if type(t) == tuple:
circuit.p(2gammaterm[t],t[0])
circuit.p(2gamma*term[t],t[1])
0sophy1
commented
3 years ago
Abstract
This document explains how to solve the knapsack problem by QAOA. There are 3 steps; prepare set of solutions, evaluate objective value of set, check if weight of solution exceeds the capacity. The circuit proceeding these 3 steps is propsed.
Description
Members
@slackhandleemail:example@example.comDeliverable
`import numpy as np import qiskit from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister from qiskit import execute, IBMQ, BasicAer, Aer from qiskit.tools import visualization from qiskit.quantum_info.operators import Operator, Pauli from qiskit.tools.visualization import plot_histogram from qiskit.tools.monitor import job_monitor import qiskit.ignis.verification as ver import matplotlib.pyplot as plt from itertools import product, combinations
IBMQ.enable_account('7318f773894831d2364969f994dae508e50e9fcc12e300030a47a0e4fb100c1739aae35a0c13dfd3a4575323f7ec88d7b9efbb5dd6741db3f79aea0dc567840b')
knapsack problem
N = 5 C = 6 cost = [3,5,9,7,4] weight = [2,3,5,4,1] beta = 0.1 np.pi gamma = 0.1 np.pi
parameter
B=1 A = B*(max(cost)+100)
function
def X(i): x = np.eye(1) for k in range(N): if k == i: x = np.kron(x, np.array([[0,1],[1,0]])) else: x = np.kron(x, np.eye(2)) return x
def Z(i, coeff): z = np.eye(1) for k in range(N): if k == i: z = np.kron(z, np.array([[coeff,0],[0,-1*coeff]])) else: z = np.kron(z, np.eye(2)) return z
def Unitary(A,k): vec = np.linalg.eig(A)[1] val = np.linalg.eig(A)[0] d_list = np.array([np.e*(ki*1j) for i in val]) diag = np.diag(d_list) a=np.dot(vec,diag) b=np.dot(a,np.linalg.inv(vec)) return b
def Cost(x): return sum([int(x[i])*cost[N-1-i] for i in range(N)])
def Weight(x): return sum([int(x[i])*weight[N-1-i] for i in range(N)])
ready for quantum circuit
qr_N = QuantumRegister(N) qr_C = QuantumRegister(C) cr = ClassicalRegister(N)
Create a Quantum Circuit acting on the q register
circuit = QuantumCircuit(qr_N,qr_C,cr,name="QAOA")
initialization
for i in range(N+C): circuit.h(i)
hamiltonian
term = dict() term["null"] = 0 for i in range(N): term[qr_N[i]] = 0 for n in range(C): term[qr_C[n]] = 0 for i in range(N): for n in range(C): term[qr_N[i], qr_C[n]] = 0 for i,j in combinations(list(range(N)),2): term[qr_N[i],qr_N[j]] = 0 for m,n in combinations(list(range(C)),2): term[qr_C[m],qr_C[n]] = 0
for i in range(N): term[qr_N[i]] += -1/2(A(weight[i]*2)-Bcost[i]) term["null"] += 1/2(A(weight[i]*2)-Bcost[i])
for n in range(C): term[qr_C[n]] += -1/2(A(n2)-1) term["null"] += 1/2(A(n2)-1)
for (i,j) in combinations(list(range(N)),2): term["null"] += 2Aweight[i]weight[j]1/4 term[qr_N[i]] += 2Aweight[i]weight[j](-1/4) term[qr_N[j]] += 2Aweight[i]weight[j](-1/4) term[qr_N[i], qr_N[j]] += 2Aweight[i]weight[j](1/4)
for (m,n) in combinations(list(range(C)),2): term["null"] += 2Amn1/4 term[qr_C[m]] += 2Amn(-1/4) term[qr_C[n]] += 2Amn(-1/4) term[qr_C[m], qr_C[n]] += 2Amn(1/4)
for i in range(N): for n in range(C): term["null"] += -2Aweight[i]n1/4 term[qr_N[i]] += -2Aweight[i]n(-1/4) term[qr_C[n]] += -2Aweight[i]n(-1/4) term[qr_N[i], qr_C[n]] += -2Aweight[i]n(1/4)
for t in term: if type(t) == qiskit.circuit.quantumregister.Qubit: circuit.p(gammaterm[t],t) if type(t) == tuple: circuit.p(2gammaterm[t],t[0]) circuit.p(2gamma*term[t],t[1])
for i in range(N): circuit.rx(2*beta, qr_N[i])
for n in range(C): circuit.rx(2*beta, qr_C[n])
Calculate statevector
backend = Aer.get_backend('statevector_simulator') job = execute(circuit, backend) result = job.result() outputstate = result.get_statevector(circuit) print("output")
experiment on simulator
circuit.measure(qr_N,cr) circuit.draw(output = 'mpl') plt.show() backend = Aer.get_backend('qasm_simulator') job = execute(circuit, backend,shots = 100) result = job.result() count = result.get_counts(circuit)
Execute on quantum Device
provider = IBMQ.get_provider() backends = provider.get_backend('ibmq_16_melbourne') job = execute(circuit, backend,shots = 100) result = job.result() counts = result.get_counts(circuit) plot_histogram(counts) plt.show()
opt_sol = "" opt_obj = 0 for sol in counts: if counts[sol] >= 5: if (Cost(sol) > opt_obj) and (Weight(sol) <= 6): opt_sol = sol opt_obj = Cost(sol)
print("optimal solution:", opt_sol) print("optimal cost:", opt_obj)`
GitHub repo