Higher Order Unconstraint Binary Optimization: HUBO

通常のイジング模型に対して高次の項を導入したモデルを考えてみます。 具体的には２値のバイナリ変数$\sigma_{i} \in \{-1, +1\}$ または $\sigma_{i} \in \{0, 1\}$に対して次のようなエネルギー関数を考えてみましょう。

$$H= c+\sum_{i} h_{i} \sigma_{i}+\sum_{i<j} J_{i j} \sigma_{i} \sigma_{j}+\sum_{i<j<k} K_{i, j, k} \sigma_{i} \sigma_{j} \sigma_{k}+\cdots, \\ i=1,\ldots ,N$$

ここで、添え字$i,j,k,\ldots$はバイナリ変数を指定するインデックスで、$1\leq i,j,k,\dots \leq N$の値を取ります。$c$は定数であり、0次の項に対応します。 このエネルギー関数の最小値を与える変数の組$(\sigma_{1} ,\sigma_{2} ,\ldots,\sigma_{N} )$を求める問題は、higher-order unconstrained binary optimization (HUBO)やpolynomial unconstrained binary optimization (PUBO)などと呼ばれています。以下ではHUBOと呼称することにします。 この手の問題は通常のイジングモデルを自然に拡張したものとみなすことができ、例えば、量子化学の分野で現れます。 このチュートリアルでOpenJijを用いたHUBOの解法を紹介します。

OpenJijでは、HUBOを直接解くことができますが、HUBOをQUBOに変換して間接的に解く手法もよく用いられています。 このチュートリアルの最後では、HUBOに対する直接解法と間接解法によって、どれだけ結果が異なるかを調べます。

HUBOの例

本チュートリアルでは、簡単のため、以下のような3次以下の項のみが現れる$N=3$変数の問題を考えます。以下では変数は$\sigma_i \in \{-1,+1\}$のスピン変数とします。

$$H = -\sigma_1 -\sigma_1\sigma_2 +\sigma_1\sigma_2\sigma_3,\;\;\; \sigma_i \in \{-1, +1\}$$

OpenJijでこの問題を解くために、まずはこのエネルギー関数の相互作用$h_1=-1,J_{1,2}=-1,K_{1,2,3}=1$をpythonの辞書型で表現します。 通常のイジング模型に対する相互作用の指定の仕方と同様に、辞書型のkeyに相互作用に関わるスピン変数を指定するインデックス、valueに対応する相互作用の値を設定します。

polynomial = {(1,): -1, (1,2): -1, (1,2,3): 1}

HUBOの直接解法

通常HUBOを解く際は相互作用の次数を2次以下に落としたQUBOを生成してこれを解きますが、 OpenJijではHUBOをQUBOに変換することなく、以下のsample_huboメソッドにより直接シミュレーテッドアニーリングを用いて、解を求めることができます。

!pip install openjij

Requirement already satisfied: openjij in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (0.9.3.dev52+g1ebe26c6)
Requirement already satisfied: numpy<1.27.0,>=1.17.3 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from openjij) (1.26.4)
Requirement already satisfied: dimod<0.13.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from openjij) (0.12.15)
Requirement already satisfied: scipy<1.12.0,>=1.7.3 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from openjij) (1.11.4)
Requirement already satisfied: requests<2.32.0,>=2.28.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from openjij) (2.31.0)
Requirement already satisfied: jij-cimod<1.7.0,>=1.6.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from openjij) (1.6.2)
Requirement already satisfied: typing-extensions>=4.2.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from openjij) (4.12.2)
Requirement already satisfied: charset-normalizer<4,>=2 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from requests<2.32.0,>=2.28.0->openjij) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from requests<2.32.0,>=2.28.0->openjij) (3.7)
Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from requests<2.32.0,>=2.28.0->openjij) (2.2.2)
Requirement already satisfied: certifi>=2017.4.17 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from requests<2.32.0,>=2.28.0->openjij) (2024.6.2)

import openjij as oj

# sample_huboメソッドを利用するためには、SASamplerメソッドを用いる必要があります。
sampler = oj.SASampler()

# sample_huboメソッドに投げます。
# sample_huboメソッドの第２引数は変数のタイプ、"SPIN"もしくは"BINARY"を指定します。
# SPINでは{-1,1}、BINARYでは{0,1}が変数として指定されます。  
response = sampler.sample_hubo(polynomial, "SPIN")
# 結果を表示します。
print(response)

   1  2  3 energy num_oc.
0 +1 +1 -1   -3.0       1
['SPIN', 1 rows, 1 samples, 3 variables]

$\sigma_1=+1,\;\sigma_2=+1,\;\sigma_3=-1$が解として得られていることが分かります。 このときのエネルギーは-3であり、これは最適解です。

なお、sample_huboメソッドは、以下のように辞書のkeyとして数値以外も扱うことができます。

response = sampler.sample_hubo({('a',): -1, ('a', 'b'): -1, ('a', 'b', 'c'): 1}, "SPIN")
print(response)

   a  b  c energy num_oc.
0 +1 +1 -1   -3.0       1
['SPIN', 1 rows, 1 samples, 3 variables]

この場合は以下のような文字列インデックスを持つエネルギー関数を最適化したことになります。 $$H = -\sigma_{\rm a} -\sigma_{\rm a}\sigma_{\rm b} +\sigma_{\rm a}\sigma_{\rm b}\sigma_{\rm c},\;\;\; \sigma_i \in \{-1, +1\}$$ インデックスが違うだけで先のモデルと全く同じモデルになります。

QUBO変換による解法

HUBOを解く方法の一つとして、3次以上の高次の項を２次以下に変換して対応するQUBOを構成しこれを解くというものがあります。 この章ではその方法を説明します。

高次相互作用から対応するQUBOを生成するためにD-Waveのdimodというライブラリを使います。 ここで、5.0に指定されているペナルティの大きさ(strength)は、高次相互作用を２次以下に変換する際に生じた制約条件に対するペナルティで、小さすぎると生成したQUBOの最適解がもとのHUBOと一致しなくなります。逆に大きすぎるとそもそも最適解が得られなくなってしまいます。実際上、どれくらいの値を指定すればいいのかというのは難しい問題ですが、ここでは一旦5.0に指定しています。

import dimod

# HUBO、strengthの大きさ、変数のタイプを指定して対応するquadraticモデルを生成します。
bqm_dimod = dimod.make_quadratic(poly=polynomial, strength=5.0, vartype="SPIN")
print('0次の項:', bqm_dimod.offset)
print('1次の項:', dict(bqm_dimod.linear))    # bqm.linearはpythonのdictに変換して表示します。
print('2次の項:', dict(bqm_dimod.quadratic)) # bqm.quadraticもpythonのdictに変換して表示します。

0次の項: 10.0
1次の項: {1: -3.5, 3: -2.5, '1*3': -2.5, 'aux1,3': -5.0, 2: 0.0}
2次の項: {(3, 1): 2.5, ('1*3', 1): 2.5, ('1*3', 3): 2.5, ('aux1,3', 1): 5.0, ('aux1,3', 3): 5.0, ('aux1,3', '1*3'): 5.0, (2, 1): -1.0, (2, '1*3'): 1.0}

見ての通りですが、もとの変数 $\sigma_1,\sigma_2,\sigma_3$ に加えて $\sigma_{"1∗2"}$、$\sigma_{"\rm aux1,2"}$ という２つの文字列で表現された変数が現れています。一般にHUBOをQUBOに変換すると変数の数、相互作用の数が増えてしまいます。 今回のケースではHUBOの場合の変数が3個、相互作用の数も3個でしたが、QUBOに変換することで、変数が5個、相互作用の数が7個に増えています。

このQUBOをOpenJijで解きたいわけですが、OpenJijでは数値と文字列が混在した変数は扱えないため、文字列を全て整数に変換する必要があります。 そのためにdimodのrelabel_variablesというメソッドを使って”1∗2”と”aux1,2”を整数に変換します。

#文字列と整数の対応関係を作る関数を定義します。
def generate_mapping(variables, N):
    mapping = {}
    #もともと整数であったN個のインデックスは変化しないようにします。
    for i in range(1, N+1):
        mapping[i] = i 
    count = N+1
    
    #新たに現れた文字列を整数に変換します。
    for v in variables:
        if type(v) == str:
            mapping[v] = count
            count += 1
    return mapping

# 変換前と変換後の変数の対応関係を表した辞書を作ります。
mapping = generate_mapping(bqm_dimod.variables, 3)

# インデックスを1始まりの整数に変換します。
bqm_dimod.relabel_variables(mapping)

print('0次の項:', bqm_dimod.offset)
print('1次の項:', dict(bqm_dimod.linear))    # bqm.linearはpythonのdictに変換して表示します。
print('2次の項:', dict(bqm_dimod.quadratic)) # bqm.quadraticもpythonのdictに変換して表示します。
print('変数の対応関係:', mapping) # Relabelした後のインデックスと元のインデックスの対応関係を表示します。

0次の項: 10.0
1次の項: {1: -3.5, 3: -2.5, 4: -2.5, 5: -5.0, 2: 0.0}
2次の項: {(3, 1): 2.5, (4, 1): 2.5, (4, 3): 2.5, (5, 1): 5.0, (5, 3): 5.0, (5, 4): 5.0, (2, 1): -1.0, (2, 4): 1.0}
変数の対応関係: {1: 1, 2: 2, 3: 3, '1*3': 4, 'aux1,3': 5}

全てのインデックスが整数に変換されました。それではこのQUBOをOpenJijを利用して解いてみましょう。

# dimodのbqmをOpenJijのBinaryQuadraticModelに変換します。
bqm_oj = oj.BinaryQuadraticModel(dict(bqm_dimod.linear), dict(bqm_dimod.quadratic), bqm_dimod.offset, vartype="SPIN")

# sampleメソッドを使ってSAを行います。
response = sampler.sample(bqm_oj)
print(response) 

   1  2  3  4  5 energy num_oc.
0 +1 +1 -1 -1 +1   -3.0       1
['SPIN', 1 rows, 1 samples, 5 variables]

ここで得られたエネルギーは先程決めたstrengthの値によってはもとのHUBOのエネルギーと対応してない可能性があります。 したがって、改めてエネルギーを計算する必要があります。 今回の場合、もともと変数は$\sigma_1,\sigma_2,\sigma_3$だけだったので、これらのスピン配位だけをつかってエネルギーを計算します。 ここでは$\sigma_1=+1,\sigma_2=+1,\sigma_3=-1$となっています。

# 元のHUBOの解に焼き直します。
hubo_configuration = {i+1: response.record[0][0][i] for i in range(3)}
print('対応するHUBOの解:', hubo_configuration)
print('対応するHUBOの解のエネルギー:', dimod.BinaryPolynomial(polynomial, "SPIN").energy(hubo_configuration))

対応するHUBOの解: {1: 1, 2: 1, 3: -1}
対応するHUBOの解のエネルギー: -3.0

ここではエネルギーとして-3を与えるスピン配位が得られました。これは今回のエネルギー関数の最適解になっています。

このことはdimodのExactPolySolverという厳密な最適解を求めるソルバーを使って実際に確認することができます。

# 元のHUBOの厳密な最適解を確認します。
sampleset = dimod.ExactPolySolver().sample_hising(h = {}, J = polynomial)
print('最適解:',sampleset.first.sample)
print('対応するエネルギー:',sampleset.first.energy)

最適解: {1: 1, 2: 1, 3: -1}
対応するエネルギー: -3.0

元のHUBOに対する最適解のエネルギーは確かに-3.0であることがわかります。 今回は変数の数が3個と少なかったため簡単に厳密解を求めることができますが、通常は厳密解を求めることは困難であることに注意してください。

HUBOの直接解法とQUBO変換による解法の比較

最後に, HUBOの直接解法とQUBO変換による解法を比較してみましょう。シミュレーテッドアニーリングによるシミュレーションを100回行い、各シミュレーションで得られたエネルギーを比較してみます。 まずはsample_huboメソッドを用いた場合のエネルギーを求めます。

# SAによるシミュレーションを行う回数を指定します。
num_reads = 100

# num_readsというパラメータを設定すると、その回数分SAを独立に行います。
# 今回は100に指定します。デフォルトは1になっています。
response = sampler.sample_hubo(polynomial, "SPIN", num_reads=num_reads)

# 得られたエネルギーをenergy_huboに代入します
energy_hubo = response.energies

次にQUBO変換を用いた解法によるエネルギーを求めます。 前に設定したstrengthの大きさが適切である保証がないため、毎回HUBOのエネルギーを計算し直す必要があることに注意してください。

response = sampler.sample(bqm_oj, num_reads=num_reads)

# 得られたスピン配位からhuboのエネルギーを計算し直す関数を定義します。
def calculate_true_energy(polynomial, response, N):
    energy_quad = []
    for i in range(num_reads):
        hubo_configuration = {j+1: response.record[i][0][j] for j in range(N)}
        energy_quad.append(dimod.BinaryPolynomial(polynomial, "BINARY").energy(hubo_configuration))
    return energy_quad

# 得られたエネルギーをenergy_quadに代入します
energy_quad = calculate_true_energy(polynomial, response, 3)

得られた100回分の結果をヒストグラムにして比較してみます。

!pip install matplotlib

Collecting matplotlib
  Downloading matplotlib-3.9.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (11 kB)
Collecting contourpy>=1.0.1 (from matplotlib)
  Downloading contourpy-1.2.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (5.8 kB)
Collecting cycler>=0.10 (from matplotlib)
  Downloading cycler-0.12.1-py3-none-any.whl.metadata (3.8 kB)
Collecting fonttools>=4.22.0 (from matplotlib)
  Downloading fonttools-4.53.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl.metadata (162 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 162.2/162.2 kB 40.5 MB/s eta 0:00:00
?25hCollecting kiwisolver>=1.3.1 (from matplotlib)
  Downloading kiwisolver-1.4.5-cp39-cp39-manylinux_2_12_x86_64.manylinux2010_x86_64.whl.metadata (6.4 kB)
Requirement already satisfied: numpy>=1.23 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from matplotlib) (1.26.4)
Requirement already satisfied: packaging>=20.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from matplotlib) (24.1)
Collecting pillow>=8 (from matplotlib)
  Downloading pillow-10.3.0-cp39-cp39-manylinux_2_28_x86_64.whl.metadata (9.2 kB)
Collecting pyparsing>=2.3.1 (from matplotlib)
  Downloading pyparsing-3.1.2-py3-none-any.whl.metadata (5.1 kB)
Requirement already satisfied: python-dateutil>=2.7 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from matplotlib) (2.9.0.post0)
Collecting importlib-resources>=3.2.0 (from matplotlib)
  Downloading importlib_resources-6.4.0-py3-none-any.whl.metadata (3.9 kB)
Requirement already satisfied: zipp>=3.1.0 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from importlib-resources>=3.2.0->matplotlib) (3.19.2)
Requirement already satisfied: six>=1.5 in /opt/hostedtoolcache/Python/3.9.19/x64/lib/python3.9/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)
Downloading matplotlib-3.9.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (8.3 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 8.3/8.3 MB 114.6 MB/s eta 0:00:00
?25hDownloading contourpy-1.2.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (304 kB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 305.0/305.0 kB 65.8 MB/s eta 0:00:00
?25hDownloading cycler-0.12.1-py3-none-any.whl (8.3 kB)
Downloading fonttools-4.53.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (4.6 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4.6/4.6 MB 110.5 MB/s eta 0:00:00
?25hDownloading importlib_resources-6.4.0-py3-none-any.whl (38 kB)
Downloading kiwisolver-1.4.5-cp39-cp39-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (1.6 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.6/1.6 MB 102.7 MB/s eta 0:00:00
?25hDownloading pillow-10.3.0-cp39-cp39-manylinux_2_28_x86_64.whl (4.5 MB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 4.5/4.5 MB 112.3 MB/s eta 0:00:00
?25hDownloading pyparsing-3.1.2-py3-none-any.whl (103 kB)
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 103.2/103.2 kB 33.2 MB/s eta 0:00:00
?25hInstalling collected packages: pyparsing, pillow, kiwisolver, importlib-resources, fonttools, cycler, contourpy, matplotlib
Successfully installed contourpy-1.2.1 cycler-0.12.1 fonttools-4.53.0 importlib-resources-6.4.0 kiwisolver-1.4.5 matplotlib-3.9.0 pillow-10.3.0 pyparsing-3.1.2

import matplotlib.pyplot as plt
plt.hist(energy_hubo, label='HUBO', range=(-3, 1), bins=10, alpha=0.5)
plt.hist(energy_quad, label='Through QUBO', range=(-3, 1), bins=10, alpha=0.5)
plt.legend()
plt.xlabel('Energy')
plt.ylabel('Frequency')

Text(0, 0.5, 'Frequency')

HUBOを直接解いた方がQUBO変換を用いた解法に比べてわずかに多く最適解が得られています。 ただしこの問題では例えばstrengthを1に設定するとQUBO変換する解法でより多くの最適解が得られます。 つまりはじめに設定したstrength=5という値が大きすぎたということです。 実際に確認してみましょう。

# strengthを1に設定してQUBOに変換します。
bqm_dimod = dimod.make_quadratic(poly=polynomial, strength=1.0, vartype="SPIN")

#インデックスを整数に変換してからOpenJijで解きます。
bqm_dimod.relabel_variables(mapping)
bqm_oj = oj.BinaryQuadraticModel(dict(bqm_dimod.linear), dict(bqm_dimod.quadratic), bqm_dimod.offset, vartype="SPIN")
response = sampler.sample(bqm_oj, num_reads=num_reads)
energy_quad = calculate_true_energy(polynomial, response, 3)

# ヒストグラムを表示します。
plt.hist(energy_hubo, label='HUBO', range=(-3, 1), bins=10, alpha=0.5)
plt.hist(energy_quad, label='Through QUBO', range=(-3, 1), bins=10, alpha=0.5)
plt.legend()
plt.xlabel('Energy')
plt.ylabel('Frequency')

Text(0, 0.5, 'Frequency')

今の場合QUBO変換による解法でほとんど最適解が得られています。 ただし、最適なstrengthの値は事前に分からないことに注意してください。 今の場合strength=1として最適解が得られていますが、一般のHUBOに対して適切なstrengthを決定することは難しい問題です。 また、QUBO変換を行うと、変数の数や相互作用の数が増えてしまうことも問題です。これにより余分なメモリが必要になっています。

ここまで例として扱ってきたモデルは単純すぎたので、もう少し大きい問題で両者の解法を比較してみましょう。 変数の数を$N=10$、相互作用を3次の全結合にして値を-1から+1のランダムにしてみます。 まずは相互作用を定義します。

import random

N=10
polynomial = {}
for i in range(1, N+1):
    for j in range(i+1, N+1):
        for k in range(j+1, N+1):
            polynomial[(i,j,k)] = random.uniform(-1, +1)

今までと同様に100回シミュレーションを行い得られたエネルギーを比較してみます。QUBO変換の際のstrengthは2としました。

#HUBOソルバーで直接解きます。
response = sampler.sample_hubo(polynomial, "SPIN", num_reads=num_reads)
energy_hubo = response.energies

#QUBO変換を通して解きます。
bqm_dimod = dimod.make_quadratic(poly=polynomial, strength=2, vartype="SPIN")
mapping = generate_mapping(bqm_dimod.variables, N)
bqm_dimod.relabel_variables(mapping)
bqm_oj = oj.BinaryQuadraticModel(dict(bqm_dimod.linear), dict(bqm_dimod.quadratic), bqm_dimod.offset, vartype="SPIN")
response = sampler.sample(bqm_oj, num_reads=num_reads)
energy_quad = calculate_true_energy(polynomial, response, N)

# ヒストグラムを表示します。
max_e = max(max(energy_hubo), max(energy_quad))
min_e = min(min(energy_hubo), min(energy_quad))
plt.hist(energy_hubo, label='HUBO', range=(min_e, max_e), alpha=0.5)
plt.hist(energy_quad, label='Through QUBO', range=(min_e, max_e), alpha=0.5)
plt.legend()
plt.xlabel('Energy')
plt.ylabel('Frequency')

Text(0, 0.5, 'Frequency')

このモデルではかなり差が出ているのが分かります。HUBOを直接解いたほうがよりエネルギーの低い解が得られています。 もちろんstrengthの値をより適切なものにすればQUBO変換の解法による解も改善する可能性はありますが、それを行うのは簡単ではありません。

まとめ

QUBO変換を通した解法では、QUBOに変換するための前処理やstrengthというパラメータを決定する必要があります。 一方でHUBOソルバーを利用すればそのような処理は不要になり、得られる解も(少なくとも今回取り扱ったモデルに関しては)QUBO変換を行う解法と同程度以上の解が得られることが分かりました。