common.title

Generate Image

Overview
Service overview
Terms of service

Privacy policy

Contact
Sign in
Sign up
common.title

QAOAで自然数分割問題

Yuichiro Minato

2021/07/17 10:18

#qaoa

1

自然数分割問題を計算する

自然数分割問題とは、ある自然数の集合を２つのグループA, Bに分割し、それぞれのグループ内の自然数の和が同じになるような分割方法を探す問題です。

解きたい問題のQUBOを作成します。 N個の自然数のii番目の自然数をnin_iとし、その自然数がどちらのグループに属するかをqiq_iで表します。 nin_iがグループAに属する時には qi=1q_i=1、グループBに属する時にはqi=0q_i=0とします。 ここで、２つのグループ内のそれぞれの和が等しい時に最小となるようなコスト関数EEを考えます。

この場合、

　E={∑i=1Nni∗(2qi−1)}2E=\{\sum_{i=1}^{N}n_i*(2q_i-1)\}^2

とすれば、自然数nin_iがグループＡに属するとき2qi−1=12q_i-1=1、グループBに属するとき2qi−1=−12q_i-1=-1 になりますので、グループＡとグループＢに属する自然数の和が等しいときに E=0E=0になり、異なるとE>0E>0になります。

展開すると、

　E=(∑i=1N2niqi)2−2(∑i=1N2niqi)(∑j=1Nnj)+(∑i=1Nni)2E=(\sum_{i=1}^{N}2n_iq_i)^2 - 2(\sum_{i=1}^{N}2n_iq_i)(\sum_{j=1}^{N}n_j) + (\sum_{i=1}^{N}n_i)^2

コスト関数Eは最小化すれば良いので、最後の定数項は要らなくなります。またコスト関数は大きさのみ関係あるので、全体を４で割って

　E=(∑i=1Nniqi)2−∑i=1Nniqi∑j=1NnjE=(\sum_{i=1}^{N}n_iq_i)^2 - \sum_{i=1}^{N}n_iq_i\sum_{j=1}^{N}n_j

また、qi=1q_i=1またはqi=0q_i=0のとき、qi2=qiq_i^2=q_iです。また、∑j=1Nnj\sum_{j=1}^N{n_j} はnの総和で定数ですので、 nsumn_{sum}とします。さらに展開して整理すると

　E=∑i=1N(ni2−nsumni)qi+2∑i<jninjqiqjE=\sum_{i=1}^{N}(n_i^2 - n_{sum}n_i)q_i +2 \sum_{i < j}n_in_jq_iq_j

これを行列表記すると下記のようになります。

　qubo=[n12−nsumn12n1n22n1n32n1n4...0n22−nsumn22n2n32n2n4...00n32−nsumn32n3n4...000n42−nsumn4..................]qubo = \left[\begin{array}{rrrrr}n_1^2 - n_{sum}n_1 & 2n_1n_2 & 2n_1n_3 & 2n_1n_4 & ...\\ 0 & n_2^2 - n_{sum}n_2 & 2n_2n_3& 2n_2n_4 &...\\ 0 & 0 & n_3^2 - n_{sum}n_3 & 2n_3n_4 & ...\\ 0 & 0 & 0 & n_4^2 - n_{sum}n_4 & ...\\ ... & ... & ... & ... &... \end{array} \right]

これをpythonプログラムで書き、シミュレータを実行して結果を得ます。

import numpy as np
from blueqat.wq import *
from blueqat import vqe

numbers = np.array([3,2,6,9,2,5,7,3,3,6,7,3,5,3,2,2,2])
qubo = np.zeros((numbers.size,numbers.size))
for i in range(numbers.size):
  for j in range(numbers.size):
    if i == j:
      qubo[i][i]=numbers[i]**2-numbers.sum()*numbers[i]
    elif i
[[-201.   12.   36.   54.   12.   30.   42.   18.   18.   36.   42.   18.

    30.   18.   12.   12.   12.]

 [   0. -136.   24.   36.    8.   20.   28.   12.   12.   24.   28.   12.

    20.   12.    8.    8.    8.]

 [   0.    0. -384.  108.   24.   60.   84.   36.   36.   72.   84.   36.

    60.   36.   24.   24.   24.]

 [   0.    0.    0. -549.   36.   90.  126.   54.   54.  108.  126.   54.

    90.   54.   36.   36.   36.]

 [   0.    0.    0.    0. -136.   20.   28.   12.   12.   24.   28.   12.

    20.   12.    8.    8.    8.]

ここで、blueqatにはqubo matrixをハミルトニアンに直接直す機能があります。

hamiltonian = pauli(qubo)
print(hamiltonian)
3.0*Z[0]*Z[1] + 9.0*Z[0]*Z[2] + 13.5*Z[0]*Z[3] + 3.0*Z[0]*Z[4] + 7.5*Z[0]*Z[5] + 10.5*Z[0]*Z[6] + 4.5*Z[0]*Z[7] + 4.5*Z[0]*Z[8] + 9.0*Z[0]*Z[9] + 10.5*Z[0]*Z[10] + 4.5*Z[0]*Z[11] + 7.5*Z[0]*Z[12] + 4.5*Z[0]*Z[13] + 3.0*Z[0]*Z[14] + 3.0*Z[0]*Z[15] + 3.0*Z[0]*Z[16] - 1133.5*I + 6.0*Z[1]*Z[2] + 9.0*Z[1]*Z[3] + 2.0*Z[1]*Z[4] + 5.0*Z[1]*Z[5] + 7.0*Z[1]*Z[6] + 3.0*Z[1]*Z[7] + 3.0*Z[1]*Z[8] + 6.0*Z[1]*Z[9] + 7.0*Z[1]*Z[10] + 3.0*Z[1]*Z[11] + 5.0*Z[1]*Z[12] + 3.0*Z[1]*Z[13] + 2.0*Z[1]*Z[14] + 2.0*Z[1]*Z[15] + 2.0*Z[1]*Z[16] + 27.0*Z[2]*Z[3] + 6.0*Z[2]*Z[4] + 15.0*Z[2]*Z[5] + 21.0*Z[2]*Z[6] + 9.0*Z[2]*Z[7] + 9.0*Z[2]*Z[8] + 18.0*Z[2]*Z[9] + 21.0*Z[2]*Z[10] + 9.0*Z[2]*Z[11] + 15.0*Z[2]*Z[12] + 9.0*Z[2]*Z[13] + 6.0*Z[2]*Z[14] + 6.0*Z[2]*Z[15] + 6.0*Z[2]*Z[16] + 9.0*Z[3]*Z[4] + 22.5*Z[3]*Z[5] + 31.5*Z[3]*Z[6] + 13.5*Z[3]*Z[7] + 13.5*Z[3]*Z[8] + 27.0*Z[3]*Z[9] + 31.5*Z[3]*Z[10] + 13.5*Z[3]*Z[11] + 22.5*Z[3]*Z[12] + 13.5*Z[3]*Z[13] + 9.0*Z[3]*Z[14] + 9.0*Z[3]*Z[15] + 9.0*Z[3]*Z[16] + 5.0*Z[4]*Z[5] + 7.0*Z[4]*Z[6] + 3.0*Z[4]*Z[7] + 3.0*Z[4]*Z[8] + 6.0*Z[4]*Z[9] + 7.0*Z[4]*Z[10] + 3.0*Z[4]*Z[11] + 5.0*Z[4]*Z[12] + 3.0*Z[4]*Z[13] + 2.0*Z[4]*Z[14] + 2.0*Z[4]*Z[15] + 2.0*Z[4]*Z[16] + 17.5*Z[5]*Z[6] + 7.5*Z[5]*Z[7] + 7.5*Z[5]*Z[8] + 15.0*Z[5]*Z[9] + 17.5*Z[5]*Z[10] + 7.5*Z[5]*Z[11] + 12.5*Z[5]*Z[12] + 7.5*Z[5]*Z[13] + 5.0*Z[5]*Z[14] + 5.0*Z[5]*Z[15] + 5.0*Z[5]*Z[16] + 10.5*Z[6]*Z[7] + 10.5*Z[6]*Z[8] + 21.0*Z[6]*Z[9] + 24.5*Z[6]*Z[10] + 10.5*Z[6]*Z[11] + 17.5*Z[6]*Z[12] + 10.5*Z[6]*Z[13] + 7.0*Z[6]*Z[14] + 7.0*Z[6]*Z[15] + 7.0*Z[6]*Z[16] + 4.5*Z[7]*Z[8] + 9.0*Z[7]*Z[9] + 10.5*Z[7]*Z[10] + 4.5*Z[7]*Z[11] + 7.5*Z[7]*Z[12] + 4.5*Z[7]*Z[13] + 3.0*Z[7]*Z[14] + 3.0*Z[7]*Z[15] + 3.0*Z[7]*Z[16] + 9.0*Z[8]*Z[9] + 10.5*Z[8]*Z[10] + 4.5*Z[8]*Z[11] + 7.5*Z[8]*Z[12] + 4.5*Z[8]*Z[13] + 3.0*Z[8]*Z[14] + 3.0*Z[8]*Z[15] + 3.0*Z[8]*Z[16] + 21.0*Z[9]*Z[10] + 9.0*Z[9]*Z[11] + 15.0*Z[9]*Z[12] + 9.0*Z[9]*Z[13] + 6.0*Z[9]*Z[14] + 6.0*Z[9]*Z[15] + 6.0*Z[9]*Z[16] + 10.5*Z[10]*Z[11] + 17.5*Z[10]*Z[12] + 10.5*Z[10]*Z[13] + 7.0*Z[10]*Z[14] + 7.0*Z[10]*Z[15] + 7.0*Z[10]*Z[16] + 7.5*Z[11]*Z[12] + 4.5*Z[11]*Z[13] + 3.0*Z[11]*Z[14] + 3.0*Z[11]*Z[15] + 3.0*Z[11]*Z[16] + 7.5*Z[12]*Z[13] + 5.0*Z[12]*Z[14] + 5.0*Z[12]*Z[15] + 5.0*Z[12]*Z[16] + 3.0*Z[13]*Z[14] + 3.0*Z[13]*Z[15] + 3.0*Z[13]*Z[16] + 2.0*Z[14]*Z[15] + 2.0*Z[14]*Z[16] + 2.0*Z[15]*Z[16]

こちらの式をQAOAに入れます。式が長いので多少計算に時間がかかります。

result = vqe.Vqe(vqe.QaoaAnsatz(hamiltonian, step=1)).run()
answer = result.most_common()[0][0]
print(answer)
(0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0)

自然数が２つのグループに分けられ、和が等しくなっています。

group1_string = ""
group2_string = ""
group1_sum = 0
group2_sum = 0
for i in range(numbers.size):
  if answer[i] == 0:
    group1_string+= '+' + str(numbers[i])
    group1_sum+=numbers[i]
  else:
    group2_string+= '+' + str(numbers[i])
    group2_sum+=numbers[i]

print(group1_string[1:],"=",group1_sum)
print(group2_string[1:],"=",group1_sum)
3+2+2+3+3+3+3+2+2+2 = 25

6+9+5+7+6+7+5 = 25

© 2024, blueqat Inc. All rights reserved