量子コンピュータで最適化を試してみる

# 量子コンピュータと最適化
量子コンピュータの応用として最適化を考えてみます。D-Waveをはじめとした量子アニーリング方式の適用分野としては最適化に特化していますが、量子ゲート方式でもQAOA(Quantum Approximate Optimization Algorithm)と呼ばれるアルゴリズムを利用して最適化を解くことができます。


- 量子アニーリングマシンであるD-waveの量子ビット

![D-Waveの量子ビット](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373940/43c76cd9-97ef-4b5d-469f-cfb5e6c261a2.png )

- 量子ゲートの量子回路の例

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373940/d43a2a7b-a73c-3b82-de9e-012c24eb5632.png)


そこで今回は最適化の典型的な問題を量子アニーリング方式と量子ゲート方式の両方で解いてみます。

# 自然数分割問題
最適化でよく使われる問題に自然数の分割問題があります。自然数分割問題はとても分りやすい問題です。自然数の集合を2つのグループに分割して、それぞれのグループの和が等しくするものです。

以下の例では1、3、4、6、2の数字を2つのグループに分けたものです。それぞれの和は8となります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373940/8efc269f-d378-fa43-18bc-6503e8e7d0eb.png)


# 量子アニーリング方式で解いてみる
アニーリング方式は量子アニーリングを利用できるD-Waveを用いて解いてみます。D-Waveに関してはこの3月から国内でもユーザー登録することによって無料の検証基盤であるLeapが利用できます。

Leapの登録はこちらから行うことができます。

https://cloud.dwavesys.com/leap/signup/

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373940/9127fee4-6f6f-30b2-3e38-bc011bad2714.png)

D-Waveで解くためにイジングモデルで表現してみましょう。
まずは自然数のグループを$q_i$で表してみます。1つめのグループを$q_1$、2つめのグループを$q_2$とすると$q_1$と$q_2$のグループの自然数の和が等しいときに最小になるような次のコスト関数Eを考ます。

$$E = ( \sum_{i=1}^N n_i \times (2q_i - 1)) ^2$$

上記のようにすると、
- 自然数$n_i$がグループAに属するときはq_i=1とします。そうすると$2q_i - 1 = 1$
- 自然数$n_i$がグループBに属するときはq_i=0とします、そうすると$2q_i - 1 = -1$
になります。
従って、双方のグループの和が等しい場合に$E = 1$となります。

例えば{1、3、4、6、2}がグループAの{1、3、4}とグループBの{2、6}に分割されたケースを考えてみましょう。
$$
(1*(2*1-1)+3*(2*1-1)+4*(2*1-1)+6*(2*0-1)+2*(2*0-1))^2 
   =(1 + 3 + 4 - 6 - 2)^2
 = 0

 $$

途中は省略しますが、エネルギー式Eは展開すると次のようになります。

　$$E=\sum_{i=1}^{N}(n_i^2 - n_{sum}n_i)q_i +2 \sum_{i < j}n_in_jq_iq_j$$


上記の式をD-Waveで計算するための入力形式で表します。

```python

from dwave.system.samplers import DWaveSampler
from dwave.system.composites import EmbeddingComposite
import numpy as np

# 自然数の集合を生成します

N = 5
np.random.seed(8123172)
arr = np.random.randint(1, N, N)

# qubo行列を作成します

Q = {}
numbers = arr.tolist()
print('get random')
sum_n = sum(numbers)
for i in range(len(numbers)):
  for j in range(len(numbers)):
    if i == j:
      Q[i,j]=numbers[i]**2-sum_n*numbers[i]
    elif i<j:
      Q[i,j]=2*numbers[i] * numbers[j]
    
```
D-Waveで実行して、結果を取得します。

```python
sampler = EmbeddingComposite(DWaveSampler(token='DEV-xxxxxxx', slover = 'DW_2000Q_5')) 

print('enbedding finished.....')
result = sampler.sample_qubo(Q, chain_strength=3, num_reads=1000)
print('get result from d-wave.....')
```

取得した結果から必要なデータを取り出して検証します。

その前にD-Waveから得られた結果(ここではresult)からエネルギー値、計算結果、発生頻度を取り出す関数をあらかじめ作成しておきます。

```python
def get_results(result, max):
    r = []
    R = iter(result)
    E = iter(result.data())
    print(len(result))
    i = 0
    for line in result:
        sample = next(R)
        data = next(E)
        energy = data.energy
        occurrences = data.num_occurrences
        r.append([energy, sample, occurrences])

        i = i + 1
        if i >= max:
            break
        #print(f'{energy},{sample}, {occurrences}')
    return r

def divide_groups(sample):
    keys = sample.keys()
    arr = [0] * len(keys)
    for key in keys:
        arr[key] = sample[key]
    return arr
```
結果を表示します。

```python

ret = get_results(result,5) # エネルギーの低いものから5個取り出す

for item in ret:
    ans = divide_groups(item[1])
    print('result = ', ans)
    r = range(len(ans))
    group_1 = []
    group_2 = []
    for i in r:
        if ans[i]==0:
            group_1.append(numbers[i])
        else:
            group_2.append(numbers[i])
    print('no of elements in group_1 = ',len(group_1),'sum in group_1 = ',sum(group_1))
    print('no of elements in group_2 = ',len(group_2),'sum in group_2 = ',sum(group_2))
```

## 結果
自然数の集合[1 2 3 4 2]に対して結果は以下の通りになります。
結果の配列は0がグループ1、1がグループ2に属していることを表わし、それぞれの和は
6となります。

```text
result = [0, 1, 0, 1, 0]
no of elements in group_1 = 3 sum in group_1 = 6
no of elements in group_2 = 2 sum in group_2 = 6
result = [0, 1, 0, 1, 0]
no of elements in group_1 = 3 sum in group_1 = 6
no of elements in group_2 = 2 sum in group_2 = 6
result = [1, 0, 1, 0, 1]
no of elements in group_1 = 2 sum in group_1 = 6
no of elements in group_2 = 3 sum in group_2 = 6
result = [0, 0, 0, 1, 1]
no of elements in group_1 = 3 sum in group_1 = 6
no of elements in group_2 = 2 sum in group_2 = 6
result = [1, 1, 1, 0, 0]
no of elements in group_1 = 2 sum in group_1 = 6
no of elements in group_2 = 3 sum in group_2 = 6
```

# 量子ゲート方式で解いてみる

量子ゲート方式のシミュレータおよびIBM Qのバックエンドが利用できるpythonの開発言語であつqiskitを利用して解いてみます。qiskitは量子ゲードをアセンブルして回路を作成することなく、ハイレベルなロジックが記述可能なqiskit aquaがあります。今回はそのqiskit aquaを利用して自然数の分割問題を解いてみます。

aquaのgithub上のソースコードは以下にあります。

https://github.com/Qiskit/qiskit-aqua

その中で最適化のアルゴリズムはこちらにあります。

https://github.com/Qiskit/qiskit-aqua/tree/master/qiskit/optimization/ising

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/373940/a4bacdb8-d2ee-4cc4-cde8-7c84d3959e17.png)

今回は自然数分割問題ということで、partition.pyを利用します。

## QAOAで解いてみる

IBMのqiskitを使ってQAOA(Quantum Approximate Optimization Algorithm)で解いてみます。

```python
from qiskit.aqua.translators.ising import partition as partition
import qiskit.aqua.translators.ising as ising
from qiskit.aqua.input import EnergyInput
import numpy as np

N = 5
np.random.seed(8123172)
l = np.random.randint(1, N, N)
print(l)
qubitOp, offset = partition.get_partition_qubitops(l)
algo_input = EnergyInput(qubitOp)
```

```python
from qiskit import BasicAer
from qiskit.aqua.algorithms import VQE, QAOA
from qiskit.aqua.components.optimizers import SPSA, COBYLA, POWELL
from qiskit.aqua.components.variational_forms import RY
from qiskit.aqua import QuantumInstance

# QAOAで解くための設定

#optimizer =  ()
optimizer = POWELL()
qaoa = QAOA(qubitOp, optimizer, 1, operator_mode='matrix')

# バックエンドに接続

seed = 10598
sim_backend = BasicAer.get_backend('statevector_simulator')
quantum_instance = QuantumInstance(sim_backend, seed_simulator=seed, seed_transpiler=seed)

# 実行する
result = qaoa.run(quantum_instance)
```

結果を取り出して、表示します。

```python
x = partition.sample_most_likely(result['eigvecs'][0])
print('energy:', result['energy'])
print('time:', result['eval_time'])
print('max-cut objective:', result['energy'] + offset)
#print('solution:', partition.get_graph_solution(x))
print('solution objective:', partition.partition_value(x, l))
print(l, x)
```
結果は以下の通りですが、どうもうまく行きません。

```python
result = [0. 0. 1. 1. 0.]
no of elements in group_1 = 3 sum in group_1 = 5
no of elements in group_2 = 2 sum in group_2 = 7
```
## VQEで解いてみる
気を取り直して量子変分法-VQE(Variational Quantum Eigensolver)にて解いてみます。

```python
# VQEで解くための設定

optimizer = SPSA(max_trials=3000)
ry = RY(qubitOp.num_qubits, depth=5, entanglement='linear')
vqe = VQE(qubitOp, ry, optimizer, 'matrix')
```

結果は以下の通りでうまく行きました。

```text
result = [0. 0. 0. 1. 1.]
no of elements in group_1 = 3 sum in group_1 = 6
no of elements in group_2 = 2 sum in group_2 = 6
```
# おわりに
今回は量子コンピュータのアプリケーションの中でも最適化問題について量子アニーリング方式と量子ゲート方式を利用して自然数分割問題を解いてみました。
QAOAではうまくいかなかったのはなぜという疑問点があるので、もう少し調べて試してみます。
また、最適化に関しては自然数分割問題をやりましたが、今後は両者で他のいくつかの問題にも挑戦してみます。