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

自然数分割問題とは、ある自然数の集合を２つのグループA, Bに分割し、それぞれのグループ内の自然数の和が同じになるような分割方法を探す問題です。
解きたい問題のQUBOを作成します。
N個の自然数の $i$ 番目の自然数を $n_i$ とし、その自然数がどちらのグループに属するかを $q_i$ で表します。
$n_i$ がグループAに属する時には
$q_i=1$ 、グループBに属する時には $q_i=0$ とします。
ここで、２つのグループ内のそれぞれの和が等しい時に最小となるようなコスト関数 $E$ を考えます。
この場合、

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

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

展開すると、

$E=(\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=(\sum_{i=1}^{N}n_iq_i)^2 - \sum_{i=1}^{N}n_iq_i\sum_{j=1}^{N}n_j$

また、 $q_i=1$ または $q_i=0$ のとき、 $q_i^2=q_i$ です。また、 $\sum_{j=1}^N{n_j}$ はnの総和で定数ですので、
$n_{sum}$ とします。さらに展開して整理すると</br>

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

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

$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<j:
      qubo[i][j]=2*numbers[i] * numbers[j]
    
print(qubo)

[[-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.]
 [   0.    0.    0.    0.    0. -325.   70.   30.   30.   60.   70.   30.
    50.   30.   20.   20.   20.]
 [   0.    0.    0.    0.    0.    0. -441.   42.   42.   84.   98.   42.
    70.   42.   28.   28.   28.]
 [   0.    0.    0.    0.    0.    0.    0. -201.   18.   36.   42.   18.
    30.   18.   12.   12.   12.]
 [   0.    0.    0.    0.    0.    0.    0.    0. -201.   36.   42.   18.
    30.   18.   12.   12.   12.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0. -384.   84.   36.
    60.   36.   24.   24.   24.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0. -441.   42.
    70.   42.   28.   28.   28.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0. -201.
    30.   18.   12.   12.   12.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
  -325.   30.   20.   20.   20.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
     0. -201.   12.   12.   12.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
     0.    0. -136.    8.    8.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
     0.    0.    0. -136.    8.]
 [   0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.    0.
     0.    0.    0.    0. -136.]]

ここで、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

QAOAで自然数分割問題

Yuichiro Minato

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