common.title
Cloud support

Nobisuke

Dekisugi

RAG

autoQAOA
RAG for dev
Fortune telling app
Annealing
DEEPSCORE
Translation
Overview
Service overview
Terms of service

Privacy policy

Contact
Research
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