QAOA variational parameter estimation in machine learning

I would like to look at QAOA, a parameterization of quantum adiabatic computation, which is a theory for dealing with discrete optimization problems on a quantum computer. Usually, the parameters of QAOA start randomly and converge, but I found a reference that estimates the initial parameter positions using LSTM, so I would like to do something similar.

Reference
Learning to learn with quantum neural networks

Check the QAOA circuit

Although blueqat automatically calculates QAOA, it is almost a black box.

from blueqat import vqe
from blueqat.pauli import *

#hamiltonian
hamiltonian = -1.0*Z[0] + 1.0*Z[1] + 0.5*Z[0]*Z[1]

#qaoa steps
step = 1

#Inside it the variational algorithm is running
result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step)).run()
print(result.most_common(4))

(((0, 1), 0.8616821697760932), ((1, 0), 0.09374669398046599), ((1, 1), 0.022285568121720205), ((0, 0), 0.0222855681217202))

You do not need to know much about QAOA to run it, as it automatically parameterizes and even optimizes, but this time we want to optimize the parameters, so we do not want to black box it.

QAOA is time evolved using the time evolution operator.

from blueqat import Circuit
a = [term.get_time_evolution() for term in hamiltonian]

time_evolution = Circuit()
angle = np.random.rand()*2*np.pi
for evo in a:
    evo(time_evolution, angle)
    
time_evolution.run_with_ibmq(returns='draw', output='mpl')

<Figure size 327.252x204.68 with 1 Axes>

We will also try to add time evolution of a transverse magnetic field Hamiltonian to this.

transverse = X[0] + X[1]
b = [term.get_time_evolution() for term in transverse]

time_evolution_x = Circuit()
angle2 = np.random.rand()*2*np.pi
for evo in b:
    evo(time_evolution_x, angle2)
    
time_evolution_x.run_with_ibmq(returns='draw', output='mpl')

<Figure size 267.052x204.68 with 1 Axes>

Putting these two together will bring us closer to a QAOA with p=1.

(time_evolution + time_evolution_x).run_with_ibmq(returns='draw', output='mpl')

<Figure size 507.852x204.68 with 1 Axes>

The initial state of the quantum adiabatic calculation takes the eigenstates corresponding to the transverse magnetic field Hamiltonian, so we create a |+> state.

(Circuit(2).h[:] + time_evolution + time_evolution_x).run_with_ibmq(returns='draw', output='mpl')

<Figure size 568.052x204.68 with 1 Axes>

In this case, the angle is a random value, but by optimizing this angle, the quantum adiabatic calculation can be used as a variational calculation to derive a near-optimal solution.

To look at the circuit we ran

result.circuit

Circuit(2).h[:].cx[0, 1].rz(-16.66276103353925)[1].cx[0, 1].rz(33.3255220670785)[0].rz(-33.3255220670785)[1].rx(4.517772793364023)[:]

Thus, a parameter corresponding to the Hamiltonian of the problem setup and a parameter corresponding to the transverse magnetic field Hamiltonian are found. These parameters are obtained by parameterizing the angles and solving the optimization problem with scipy minimize. The minimized values are the expected values obtained from the Hamiltonian and the state vector.

We would like to plot the values during the optimization process, starting from random and seeing how they are optimized.

The parameters to be optimizerd are set to beta and gamma, starting at random.

time_evolution = Circuit(2).h[:]
beta = np.random.rand()*2*np.pi
for evo in a:
    evo(time_evolution, beta)

gamma = np.random.rand()*2*np.pi
for evo in b:
    evo(time_evolution, gamma)

time_evolution.run_with_ibmq(returns='draw', output='mpl')

<Figure size 568.052x204.68 with 1 Axes>

The state vector after running this circuit is

sv = time_evolution.run()
sv

array([ 0.30809098-0.44281209j,  0.21083871-0.43441041j,
       -0.4292161 -0.02455541j,  0.30809098-0.44281209j])

We use this and the Hamiltonian to find the expectation value.

#function to get expectation value from state vector
def sv_expt(sv_tmp, h_param):
    z0 = np.array([[1,0],[0,-1]])
    i0 = np.eye(2)
    h0 = h_param[0] * np.kron(z0, i0) + h_param[1] * np.kron(i0, z0) + h_param[2] * np.kron(z0, z0)
    expt = sv_tmp.conj()@h0@sv_tmp
    return expt.real

#try a problem
h_init = [-1, 1, 0.5]
expt = sv_expt(sv, h_init)
print(expt)

-0.014666664541780838

Optimize parameters.

import copy

arr_expt = [expt] #history of expectation value
arr_beta = [beta] #hitstory of angle beta
arr_gamma = [gamma] #history of angle gamma

step = 1000

#return the expectation value from the quantum circuit
def te_circuit(te_a, te_b, angle_beta, angle_gamma, h_param):
    te_tmp = Circuit(2).h[:]

    for evo in te_a:
        evo(te_tmp, angle_beta)

    for evo in te_b:
        evo(te_tmp, angle_gamma)

    sv_tmp = te_tmp.run()
    return sv_expt(sv_tmp, h_param)

#loop of optimization
for _ in range(step):
    
    #for differentiation
    h = 0.01
    
    #learning rate
    e = 0.1
    
    #get the expectation value
    expt = te_circuit(a, b, beta, gamma, h_init)

    #expectation value of beta + h
    beta_tmp = copy.copy(beta) + h
    expt_delta_beta = te_circuit(a, b, beta_tmp, gamma, h_init)

    #expectation value of gamma + h
    gamma_tmp = copy.copy(gamma) + h
    expt_delta_gamma = te_circuit(a, b, beta, gamma_tmp, h_init)

    #update the angle using differentiation
    beta -= (expt_delta_beta - expt) * e
    gamma -= (expt_delta_gamma - expt) * e
    
    #check the updated expectation value
    expt = te_circuit(a, b, beta, gamma, h_init)

    arr_beta.append(beta)
    arr_gamma.append(gamma)
    arr_expt.append(expt)

Let's see the history of expectation value.

import matplotlib.pyplot as plt

plt.plot(arr_expt)
plt.show()

<Figure size 432x288 with 1 Axes>

And then look at the optimized history of angles.

plt.plot(arr_beta, label="beta")
plt.plot(arr_gamma, label="gamma")

plt.legend()
plt.show()

<Figure size 432x288 with 1 Axes>

Next, since p=1 in this case, I illustrated how the angle affects the minimum value.

reference : https://pennylane.ai/qml/demos/learning2learn.html

grid = 30

beta_grid = np.linspace(0, 2*np.pi, grid)
gamma_grid = np.linspace(0, 2*np.pi, grid)
Beta, Gamma = np.meshgrid(beta_grid, gamma_grid)

data = np.zeros((grid, grid))

for i in range(grid):
    for j in range(grid):
        data[i][j] = te_circuit(a, b, beta_grid[i], gamma_grid[j], h_init)
        
#draw the contourf
plt.contourf(Beta, Gamma, data)
plt.colorbar()

#draw the optimized trajectory of angles
plt.plot(arr_beta, arr_gamma, linestyle="--", color="red")
plt.plot(arr_beta[-1], arr_gamma[-1], linestyle="--", color="red", marker="x")

plt.xlabel(r"$\beta$")
plt.ylabel(r"$\gamma$")

plt.show()

<Figure size 432x288 with 2 Axes>

You can see that it is gradually moving from the initial state to a lower energy state with darker color.

The graph itself is easier to understand when viewed in three dimensions.

from mpl_toolkits.mplot3d import axes3d

ax = plt.figure().add_subplot(projection='3d')
ax.contour3D(Beta, Gamma, data, 200, cmap='rainbow')
plt.show()

<Figure size 432x288 with 1 Axes>

We have seen how the angle parameter is optimized using mathematical optimization to find a solution to QAOA. Next, we would like to estimate the initial value of the angle of this cost function and train the classical NN to start as close to the solution as possible. First, we create the data. This time, we will arbitrarily create a 2-qubit Hamiltonian and learn its converged expectation value as a label.

def make_data():
    rand1 = np.random.rand()*2-1
    rand2 = np.random.rand()*2-1
    rand3 = np.random.rand()*2-1
    
    h_make = [rand1, rand2, rand3]
    
    hami1 = rand1*Z[0] + rand2*Z[1] + rand3*Z[0]*Z[1]
    transverse2 = X[0] + X[1]
    
    aa = [term.get_time_evolution() for term in hami1]
    bb = [term.get_time_evolution() for term in transverse2]
    
    beta = np.random.rand()*2*np.pi
    gamma = np.random.rand()*2*np.pi

    min1 = 0
    min2 = rand1
    min3 = rand2
    min4 = rand1 + rand2 + rand3
    minimum = min(min1, min2, min3, min4)
    
    step = 1000
    h = 0.01
    e = 0.1
    
    for _ in range(step):
        expt = te_circuit(aa, bb, beta, gamma, h_make)

        beta_tmp = copy.copy(beta) + h
        expt_delta_beta = te_circuit(aa, bb, beta_tmp, gamma, h_make)

        gamma_tmp = copy.copy(gamma) + h
        expt_delta_gamma = te_circuit(aa, bb, beta, gamma_tmp, h_make)

        beta -= (expt_delta_beta - expt) * e
        gamma -= (expt_delta_gamma - expt) * e
    
    #更新された値をチェック
    expt = te_circuit(a, b, beta, gamma, h_init)
        
    if expt < minimum+0.5:
        return [rand1, rand2, rand3], [beta, gamma]
    else:
        return 0

Create training data and separate it into training and evaluation data.

from sklearn.neural_network import MLPRegressor
from sklearn.model_selection import train_test_split

X_data = []
y_data = []

for _ in range(1000):
    get_data = make_data()
    if get_data != 0:
        X_data.append(get_data[0])
        y_data.append(get_data[1])
        print("*", end="")
    else:
        print("x", end="")

print("")
print(len(X_data))

X_train, X_test, y_train, y_test = train_test_split(X_data, y_data, test_size=0.2, random_state=0)

xx******xxxxxx******x*x*x*xx***x**xx***x*****x**xx**x**x*xx**xxx***x****x*x****x*x*xxxx**x*xx*x****x*xxx***xxx**xx*x*xx*xxx*x****x**xx*x**xxx*xx*x*x*x*xxx*x*xx**xxxxxxxxx*xx*x*x*x**x**x*x*xx**x***x***xxxxxx**xxxxxx*xxxx**x*xx*xx*x*xx**xx*x*x*xxxx**xxxxx***x*x*x*****xx***x***xx*x*xxxx*x*xxx**xxx*****x**xxxx***x******xxxx****x*xxxxx*xxxx*x***xx*x*xxx**xx**xx***x*xxx*xx*xxxxxx*x****x*xx***x*x***xx**x*x***x*x***x***xx**xxxx*x*xx****xxx*x****xx*xx**xx*xx*x*x*x**x*xx*xxx**xxx***xx****xx***xxxx*xxx*x**x*xxx*xx*xxxxx**xxx*xx****x**x*x**x*xxx*****x*xxxxxxx*xx**x**xxx*xx*x**x*xx*x*xx*x**xx**x*x*x*x**x**xxx**xx*******x*x*xx*xxxxx*x*x*xxxxxxx**x*xx*xxx**xx*xx*xx**x**x*xx*x*x*xxx*x***x*xxx***x**x***x***xxx*****xx*x**x***xx*xxx*x*x**xxxxx**xx*xx*xxx*xx*****x**x**x*x******xxx*xxxx*x*x*x*xxxxxxxxx*xx********xxxx*x*x*x*x**x*xx**x***x*xx*xx****x**xx**xxx*x*xxxxxxxxxxxxx***xx***x**x*x*x**xx**x**x*x**x*xxxx***xxx***x*xxxx**xx****x***xx*x**x**x*xx*xx*x**x*x*xxxx**x***xxxxx*xxx**x**xx*xx**x***x*xxxx***xxxxx
483

Evaluate the trained model.

model_nn = MLPRegressor(random_state=42)
model_nn.fit(X_train, y_train)

score = model_nn.score(X_train, y_train)
print("Training score:", score)

score2 = model_nn.score(X_test, y_test)
print("Valuation score:", score2)

Training score: 0.01067089406467775
Valuation score: -0.020423018316892194

/opt/conda/lib/python3.8/site-packages/sklearn/neural_network/_multilayer_perceptron.py:614: ConvergenceWarning: Stochastic Optimizer: Maximum iterations (200) reached and the optimization hasn't converged yet.
  warnings.warn(

We can actually estimate the variational parameters from the new Hamiltonian.

#Create Hamiltonian from random coefficients
rand1_pred = np.random.rand()*2-1
rand2_pred = np.random.rand()*2-1
rand3_pred = np.random.rand()*2-1
h_pred = [rand1_pred, rand2_pred, rand3_pred]

#Estimating angles
x_pred = model_nn.predict([h_pred])
x_pred = x_pred[0]
print(x_pred)

beta = x_pred[0]
gamma = x_pred[1]

#Create time development from Hamiltonian
hami_pred = rand1_pred*Z[0] + rand2_pred*Z[1] + rand3_pred*Z[0]*Z[1]
transverse = X[0] + X[1]

a_pred = [term.get_time_evolution() for term in hami_pred]
b_pred = [term.get_time_evolution() for term in transverse]

expt_pred = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
print(expt_pred)

arr_expt = [expt_pred]
arr_beta = [beta]
arr_gamma = [gamma]

step = 1000
h = 0.01
e = 0.1
   
for _ in range(step):
    expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)

    beta_tmp = copy.copy(beta) + h
    expt_delta_beta = te_circuit(a_pred, b_pred, beta_tmp, gamma, h_pred)

    gamma_tmp = copy.copy(gamma) + h
    expt_delta_gamma = te_circuit(a_pred, b_pred, beta, gamma_tmp, h_pred)

    beta -= (expt_delta_beta - expt) * e
    gamma -= (expt_delta_gamma - expt) * e 
    
    expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)

    arr_beta.append(beta)
    arr_gamma.append(gamma)
    arr_expt.append(expt)

[3.19008567 3.48471883]
-0.1294131295898297

Next, let's generate it at random as well.

beta = np.random.rand()*2*np.pi
gamma = np.random.rand()*2*np.pi

expt_pred = te_circuit(a_pred, b_pred, beta, gamma, h_pred)
print(expt_pred)

arr_expt2 = [expt_pred]
arr_beta2 = [beta]
arr_gamma2 = [gamma]
   
for _ in range(step):
    expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)

    beta_tmp = copy.copy(beta) + h
    expt_delta_beta = te_circuit(a_pred, b_pred, beta_tmp, gamma, h_pred)

    gamma_tmp = copy.copy(gamma) + h
    expt_delta_gamma = te_circuit(a_pred, b_pred, beta, gamma_tmp, h_pred)

    beta -= (expt_delta_beta - expt) * e
    gamma -= (expt_delta_gamma - expt) * e 
    
    expt = te_circuit(a_pred, b_pred, beta, gamma, h_pred)

    arr_beta2.append(beta)
    arr_gamma2.append(gamma)
    arr_expt2.append(expt)

0.013404047639859196

Compare convergence with random generation

import matplotlib.pyplot as plt
plt.plot(arr_expt, label="from NN")
plt.plot(arr_expt2, label="random")
plt.legend()
plt.show()

<Figure size 432x288 with 1 Axes>

Compare convergence with random generation

grid = 30

beta_grid = np.linspace(0, 2*np.pi, grid)
gamma_grid = np.linspace(0, 2*np.pi, grid)
Beta, Gamma = np.meshgrid(beta_grid, gamma_grid)

data = np.zeros((grid, grid))

for i in range(grid):
    for j in range(grid):
        data[i][j] = te_circuit(a_pred, b_pred, beta_grid[i], gamma_grid[j], h_pred)
        
plt.contourf(Beta, Gamma, data)
plt.colorbar()

plt.plot(arr_beta, arr_gamma, linestyle="--", color="red")
plt.plot(arr_beta[-1], arr_gamma[-1], linestyle="--", color="red", marker="x")

plt.plot(arr_beta2, arr_gamma2, linestyle="--", color="blue")
plt.plot(arr_beta2[-1], arr_gamma2[-1], linestyle="--", color="blue", marker="x")

plt.xlabel(r"$\beta$")
plt.ylabel(r"$\gamma$")

plt.show()

<Figure size 432x288 with 2 Axes>

from mpl_toolkits.mplot3d import axes3d

ax = plt.figure().add_subplot(projection='3d')
ax.contour3D(Beta, Gamma, data, 200, cmap='rainbow')
plt.show()

<Figure size 432x288 with 1 Axes>

I tried it several times, and although it may be somewhat comforting, I think that the parameters of the one that was estimated to some extent were more likely to converge overall. Sometimes random happens to get better values. That's all.

QAOA variational parameter estimation in machine learning

Yuichiro Minato

QAOA variational parameter estimation in machine learning

Check the QAOA circuit