模擬焼き鈍し法用Qubo生成pythonパッケージquboassistの使い方解説

著者：榎本 觀
連絡先：[e-mail]

概要

本記事では、筆者が開発したQubo生成pythonパッケージquboassistの使用方法を解説する。

導入

まず、Quboとは二次制約なし二値最適化問題（Quadratic unconstrainted binary optimization）の略で、即ち

$${\rm minimize\ \ \ }f(x) = x^t A x,\ \ s.t.\ \ x \in \{0,1\}^N$$

という問題である。ここで、$A$は$N$次実対称行列である。

quboassistは、各制約条件へのある重み付けの下で

$${\rm minimize}\ \ \ f(x) =x^t Ax\\ s.t. \ \ \ \forall i, a_i \leq x_i\leq b_i,x_i\in\mathbb Z,\ \ \forall j, \sum_ kc_{jk} x_k \geq d_j\ \ (c_{jk}, d_j \in \mathbb Z)$$

という問題をQuboに変換する。
この問題は、目的函数が二次であり、各変数は有界な範囲の中で整数値をとり、かつ各制約条件は整数係数且つ線型である。
巡回セールスマン問題など、様々な計算量理論的に難しい問題がこの形式に定式化出来る。

本パッケージは、内部アルゴリズムのコア部分を、効率的なデータ構造を用いてRustで並列計算を駆使して実装し、また直感的な問題入力インタフェースをpythonで実装している為、高度な計算を簡便に行うことが可能になる。変換の理論的基礎については、github（https://github.com/enomotokan/quboassist）やPyPI（https://pypi.org/project/quboassist/）にあるREADME.mdファイルの最後の辺りを参照されたい。

使用方法

以下、本パッケージquboassistの使用方法を解説する。

まずはpythonがインストールされている環境で
pip install quboassist
してローカル環境にquboassistのビルド済みファイルをインストールする。記事執筆時点（2024/11/15）でpythonバージョンは3.13以下3.8以上、OSはubuntu、MacOS、Windowsに対応している。

早速使ってみよう。コマンドラインでpythonを起動し、
from quboassist import *
で内部で定義されている三つのクラスVariable,Formula,Problemをインポートする。quboassistで変数を定義するには、Variableクラスのインスタンスを作る。

x = Variable("x", 2, 5)

これで、$2,3,4,5$の値をとる$x$という名前の整数変数を一つ定義出来たことになる。一度に大量の変数を定義したい場合は、リスト内包表記で

x = [Variable("x{i}".format(i), 2, 5) for i in range(10000)]

とすれば良い。次は、定義した変数を用いて式を作りたい。例えば、

$$f = - x_0^2 -3 x_1$$

に対応する式をquboassistで作るには、

f = - x[0]**2 -  3 * x[1]

とする。内部では、Variableクラスのインスタンスである変数たちがFormulaクラスのインスタンスfとなり、その係数情報がRustベクタでメモリ上に載る。
そもそも解きたい問題が二次なので、式も二次式までのみ取り扱え、それ以上を取り扱おうとするとエラーになる。

不等式や等式を作るには、例えば

$$x_0 > x_1$$

に対応する式を

g = x[0] > x[1]

としてgに格納する。これらを用いて、問題を定義することが出来る。

P = Problem() 
P.add_objective(f)
P.add_constraint(g)

PにProblemクラスのインスタンスを格納し、その目的函数を上で定義したfとし、制約条件をgとする。これで、最適化問題

$${\rm minimize}\ \ \ f(x) = - x_0^2 -3 x_1 \ \ \ s.t.\ \ x \in \{2,3,4,5\}^2, x_0 > x_1$$

をPに格納できたことになる。実際に逐次代入することで簡単に確かめることが出来るが、この問題の最適解は

$$(x_0, x_1) = (5, 4)$$

である。
この問題をQubo形式に変換するには、例えば制約条件gの重み付けを$10$として

qubo = P.compile([10])

とする。今回は一つのみであるが、各制約条件の重みを表すベクトルをcompileメソッドに代入すると、Quboが生成されてquboに格納される。内部では、まず全ての不等式制約に対して補助変数を加えて等式制約に変換し、次に補助変数を含め全ての変数を$0,1$二値変数の線型結合で表現し、それら等式制約の二乗を重みを掛けて目的函数に加える。即ち、ラグランジュ緩和である。

後は、Quboをソルバーに代入して問題を解いてみよう。
pip install dwave-neal
で模擬焼き鈍し法（シミュレーテッドアニーリング）ソルバーをインストールしておく。

sampler = neal.SimulatedAnnealingSampler()

でソルバーをsamplerに格納する。以下、既に生成したquboをtodictメソッドでsamplerの受け付ける辞書型に変換して解かせる。

result = sampler.sample_qubo(qubo.todict()).first.sample

resultには、変換した問題の二値変数の値が辞書式で格納されている。solutionメソッドで、resultから元の問題の解を構成できる。

solution = P.solution(result, "neal")
print(solution)

出力される解は、ほとんど$1$に近い確率で以下のようになるはずである。

({'x0': 5, 'x1': 4}, [True])

注意してほしいのは、模擬焼き鈍し法は発見的手法（ヒューリスティクス）であり、必ずしも最適解を出力するわけではないことである。さらに、本パッケージの場合ラグランジュ緩和で元の問題を無理やり制約のない問題に変換しているので、必ずしも出力される解が全ての制約条件を充たすとは限らぬのである。
その情報はsolutionの二番目にあり、これはそれぞれの制約条件が充たされているか否かを表す真理値の格納されたリストである。

重み付けを如何にすべきか

重みは、$0$の場合は完全に制約条件を無視し、巨大になれば目的函数を無視することになるので、その中間のいい塩梅を、重みを調節しながら何回も最適化を実行することで見つける必要がある。実際、Quboはこのような使われ方をすることが多いため、quboassistの内部では、問題Pの前回のcompile時に既に追加されていた制約条件に対応するQuboは保存され、その繰り返しの中で同じ計算を反復することを避けている。

前章の例題では、制約条件の重みを$10$として解けば高い確率で最適解が得られるが、一般の、大量の変数と制約条件を持つ問題に対し、全ての制約条件の重み付けを手動で調節するのは大変である。そこで、ベイズ最適化によるハイパーパラメータ調節パッケージであるoptunaを用いることを薦めたい。optunaそのものの解説については他の文献に譲るとして、ここではoptunaを併用することで上の目的が達成されることを見る。

コマンドラインで
pip install optuna
しておく。

最適化問題を焼き鈍し法で解く時のハイパーパラメータである重み$w$は、全ての制約条件を充たす程度に大きく、且つ目的函数をなるべく最小化できる程度に小さくしたい。そこで例えば、問題を解いたときの出力の二番目の真理値を並べたベクトル$b$に対し、目的函数を定数$M>0$を用いて

$$f(w,b) = \sum_i \left\{ w_i + M (1 - b_i) \right \}$$

と定義すれば上の目標を実現する目的函数となる。
ここで、それぞれ真理値$\rm True$と整数$1$、$\rm False$と$0$を同一視している。実際、pythonではこのように暗黙に型変換が行われる。

この目的函数を用いて、最適化問題

$${\rm minimize}\ \ f(x) = -x_0^2 - x_1$$

$$s.t.\ \ x \in \{0,1,2,3\}^3, x_0 > x_1, x_0 + x_1 = x_2$$

を解くためのサンプルコードを以下に示す。このサンプルコードでは$M = 10$としているが、これはそれぞれの問題に対して適宜調節されたい。

import neal
import quboassist
import optuna
import numpy as np
from copy import copy

x = [quboassist.Variable("x{}".format(i), 0, 3) for i in range(10)]

best_solution = []
best_val = np.inf

P = quboassist.Problem()

f = - x[0]**2 - x[1]
P.add_objective(f)

g = x[0] > x[1]
P.add_constraint(g)

h = x[0] + x[1] == x[2]
P.add_constraint(h)

sampler = neal.SimulatedAnnealingSampler()

def objective(trial):
    
    w = [trial.suggest_float("w{}".format(i), 0, 5) for i in range(2)  
    qubo = P.compile(w)

    result = sampler.sample_qubo(qubo.todict()).first
    solution = P.solution(result.sample)

    print("\n")
    print(solution)

    obj = w[0] + w[1] + 10 * sum(np.logical_not(solution[1]))
    val = result.energy
    
    # Note that result.energy and the value of objective function may differ by a constant which appears when expanding the product of variables!

    global best_solution, best_val
    
    if np.all(solution[1]) and val < best_val:
        best_val = val
        best_solution = copy(solution)
    
    return obj

study = optuna.create_study(direction="minimize")
study.optimize(objective, n_trials=20)

print("\n\nBest Solutuion")
print(best_solution)

追記

筆者は、普段は数理最適化専門のフリーランスエンジニアとして活動しています。詳しい自己紹介は以のwantedlyページを参考にして下さい。
https://www.wantedly.com/id/enomotokan

もし本記事やコードに関して疑問点などあれば以下のメールアドレスにご連絡下さい。
また、数理最適化及びヒューリスティクスに関する技術の相談、開発の依頼なども同メールアドレス宛に連絡下されば幸いです。

[e-mail]