PennylaneでQAOA w/勾配降下

そういえばPennylaneでQAOAを勾配降下で解く実装を記事にしていませんでした。

# 全結合MAX CUT をQAOAで解く

例えば10ノードの全結合グラフのMAX CUT を考えます。

```python
import pennylane as qml
from pennylane import qaoa
from pennylane import numpy as np
from networkx import Graph, erdos_renyi_graph
import matplotlib.pyplot as plt

n_qubits = 10 # number of qubits
n_layers = 1
graph = erdos_renyi_graph(n_qubits, p=1, directed=False)
# Defines the wires and the graph on which MaxCut is being performed
wires = range(n_qubits)

import networkx as nx
nx.draw(graph, with_labels = True)
plt.show()
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/60d2a2f9-cef9-5ccf-9369-e04c0a0b0963.png)

QAOA Ansatzを組み立てるのに必要なパーツを先に定義します。

```python
# Defines the QAOA cost and mixer Hamiltonians
cost_h, mixer_h = qaoa.maxcut(graph)

# Defines a layer of the QAOA ansatz from the cost and mixer Hamiltonians
def qaoa_layer(gamma, alpha):
    qaoa.cost_layer(gamma, cost_h)
    qaoa.mixer_layer(alpha, mixer_h)
```

## 勾配降下

```python
# Repeatedly applies layers of the QAOA ansatz
dev = qml.device('default.qubit', wires=len(wires))
@qml.qnode(dev,diff_method='backprop')
def cost_function(params, **kwargs):
    for w in wires:
        qml.Hadamard(wires=w)
    qml.layer(qaoa_layer, n_layers, params[0:n_layers,0], params[0:n_layers,1])
    return qml.expval(cost_h)

opt = qml.GradientDescentOptimizer(0.01)
hist_cost = []
var_init = np.random.randn(n_layers,2)
var = var_init

import time
# 処理前の時刻
t1 = time.time() 
for it in range(10):
    var, _cost = opt.step_and_cost(cost_function, var)
    print("Iter: {:5d} | Cost: {:0.7f} ".format(it, _cost))
    hist_cost.append(_cost)
    
# 処理後の時刻
t2 = time.time() 
# 経過時間を表示
elapsed_time = t2-t1
print(f"経過時間：{elapsed_time}")
```
`Iter:     0 | Cost: -15.4306719 
Iter:     1 | Cost: -19.8801702 
Iter:     2 | Cost: -22.1785671 
Iter:     3 | Cost: -22.4678397 
Iter:     4 | Cost: -22.4965657 
Iter:     5 | Cost: -22.4996299 
Iter:     6 | Cost: -22.4999601 
Iter:     7 | Cost: -22.4999958 
Iter:     8 | Cost: -22.4999996 
Iter:     9 | Cost: -22.5000000 
経過時間：2.244999885559082`

```python
plt.plot(hist_cost,'o-')
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/e645513d-6ed2-740f-faa4-7122b806f00a.png)

##　非勾配
scipy.minimizeで非勾配な最適化をしてみます。

```python
dev = qml.device('lightning.qubit', wires=len(wires))
# Repeatedly applies layers of the QAOA ansatz
@qml.qnode(dev)
def cost_function(var, **kwargs):
    var_2d_array = np.reshape(var,(n_layers,2))
    for w in wires:
        qml.Hadamard(wires=w)
    qml.layer(qaoa_layer, n_layers, var_2d_array[0:n_layers,0], var_2d_array[0:n_layers,1])
    return qml.expval(cost_h)

from scipy.optimize import minimize
import time

var_init = np.random.uniform(low=-1, high=1, size=(2*n_layers)) # one-dimensional array
hist_cost = []
var = var_init

count = 0
def cbf(Xi):
    global count
    global hist_cost
    cost_now = cost_function(Xi)
    hist_cost.append(cost_now)
    print('iter = '+str(count)+' | cost = '+str(cost_now))
    count += 1


t1 = time.time()
result = minimize(fun=cost_function, x0=var_init, method='Nelder-Mead', callback=cbf, options={"maxiter":11})
# 処理後の時刻
t2 = time.time() 
# 経過時間を表示
elapsed_time = t2-t1
print(f"経過時間：{elapsed_time}")
```

`iter = 0 | cost = -24.05773509630845
iter = 1 | cost = -24.119313285482026
iter = 2 | cost = -24.24372585868572
iter = 3 | cost = -24.37502879518174
iter = 4 | cost = -24.40971631000042
iter = 5 | cost = -24.472231940032042
iter = 6 | cost = -24.472231940032042
iter = 7 | cost = -24.472231940032042
iter = 8 | cost = -24.50816706605304
iter = 9 | cost = -24.50816706605304
経過時間：3.0799999237060547`

```python
plt.plot(hist_cost,'o-')
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/878664/7aa73afb-4705-976d-1845-9bf76a8381d6.png)

この例では、所要時間はあまり変わりません。

しかしQAOAのレイヤー数、つまりパラメータの数を増やすと劇的に違います。
例えばレイヤー数を10、scipy.minimizeのmethodをmethod='Powell'とした場合では、
非勾配は数分かかりますが勾配法は30秒程度で終了します。
パラメータ数が数十、数百となってくると、勾配法のメリットが出てきます。