列生成法をQUBO求解に応用する論文 (Hirama and M. Ohzeki (2023)) の紹介

この記事はJij Inc. Advent Calendar 2023の20日目の記事です。
はじめまして、株式会社Jijの西村光嗣です。

概要

QUBO形式の問題に制約条件が加わった問題を数理最適化の手法である列生成法を用いて解く論文である
S. Hirama and M. Ohzeki, Efficient Algorithm for Binary Quadratic Problem by Column Generation and Quantum Annealing, J. Phys. Soc. Jpn. 92, (2023).
の紹介をする。

はじめに

目的関数が2次式で、かつ変数$x_i$が0か1しかとらない次のような問題

$$H = \sum_{i < j} Q_{ij} x_i x_j$$

はQUBO (Quadratic Unconstrained Binary Optimization)と呼ばれる種類であり、組合せ最適化問題をこの形式に帰着させた後、量子アニーリング等を含むメタヒューリスティクスを用いて求解を行う手法が注目を浴びている。

一方、実用的な問題には、例えば変数の総和がある値$C$にならなければならない：

$$\sum_{i} x_i = C$$

といったような制約条件が多く含まれ、この制約条件をQUBO形式で扱うために

$$H = \sum_{i < j} Q_{ij} x_i x_j + \lambda \left( \sum_{i} x_i - C \right)^2$$

のように元の制約条件に対応する量にペナルティをかけて対処する必要があった。

本論文は列生成法と呼ばれる数理最適化の手法を用いることにより、制約ありのQUBO問題を解く手法を提案している。

列生成法

列生成法のわかりやすい具体例はブレインパッドのブログ記事が参考になる。

https://blog.brainpad.co.jp/entry/2020/10/23/000003

列生成法のエッセンスを上記のブログに沿って説明すると以下の通りとなる。

今、複数の配送ルートと、それぞれの配送ルートに関連する制約が存在したとして、最適な配送ルートの組を選択したい状況を考える。

第3回：はじめての配送計画の列生成法【ブレインパッドの数理最適化ブログ】(https://blog.brainpad.co.jp/entry/2020/10/23/000003) より引用

数式で記述すると以下のようになる。

$$\begin{align} \min_{x} &\sum_{r \in R} d_r x_r \\ \mathrm{s.t. }\ & \sum_{r \in R}^{n} a_{ri}x_r \geq 1 \ \forall i \in S\\ x_r &= \{0, 1\} \end{align}$$

ここで、$x_r$が配送ルート$r$を選択するか否かを0, 1 で表現した変数であり、上記の問題は係数 $d_r$の重みにより定義されたコスト関数を最小化する$x_r$ の組合せを決定する問題である。ただし、制約式 (2)を守る必要がある。

しかしながら。もし配送ルート候補の数 (つまり集合$R$の要素数)が途方もない数存在する場合、上記の数理最適化問題を直接解くことは極めて難しくなる問題がある。

考えるうる配送ルートの中には、明らかに最適解からかけ離れていて実用的には考慮する必要がないものがあると考えると、考慮すべき配送ルートの候補の数を削減したいというモチベーションがある。

ここで、上記の問題にて一旦$x_r$を$0 \leq x_r \leq 1$ の値を取る連続変数とひとまず置くことで、元の問題を緩和した問題を定義してみる。

$$\begin{align} \min_{x} &\sum_{r \in R} d_r x_r \\ \mathrm{s.t. }\ & \sum_{r \in R}^{n} a_{ri}x_r \geq 1 \ \forall i \in S\\ 0 &\leq x_r \leq 1 \end{align}$$

こうすることで、元の問題の双対問題が定義できる。

$$\begin{align} \max_{y} &\sum_{i \in S} y_i \\ \mathrm{s.t. }\ & \sum_{i \in S} a_{ri}y_i \leq d_r \ \forall r \in R\\ \end{align}$$

配送ルート候補の数 (集合$R$の要素数)が途方もない数存在する状況は、双対問題において、(8)で定義された形状の制約が大量に存在する場合に対応する、すなわち、多くの制約条件が意味をなさず、一部の制約条件のみが寄与していることが予想できる。イメージ図は以下のようなものとなる。

はじめての列生成法, 宮本　裕一郎 (https://orsj.org/wp-content/corsj/or57-4/or57_4_198.pdf) より引用

これにより、以下のアルゴリズムが構築できる。

列生成法アルゴリズム

1. 配送ルート候補の数を適当に決め、その集合を$R_0$ とする。

2. ルート候補を選択する問題の緩和問題である

   $$\begin{align} \min_{x} &\sum_{r \in R_0} d_r x_r \\ \mathrm{s.t. }\ & \sum_{r \in R_0} a_{ri}x_r \geq 1 \ \forall i \in S\\ 0 &\leq x_r \leq 1 \end{align}$$

   を解く。これにより、双対問題である

   $$\begin{align} \max_{y} &\sum_{i \in S} y_i \\ \mathrm{s.t. }\ & \sum_{i \in S} a_{ri}y_i \leq d_r \ \forall r \in R_0\\ \end{align}$$

   の解 $y_i^{*}$ も同時に解けたことになる。

3. 双対問題の解 $y_i^{*}$ のもとで、制約式 (13) を破るような$d_r$、$a_{ri}$ を探す。すなわち、

$d_r - \sum_{i \in S} a_{ri} y_i^{*} < 0$ を満たす $d_r$、$a_{ri}$ を探す。これを達成するには$d_r - \sum_{i \in S} a_{ri} y_i^{*}$ を最小化する $d_i, a_{ri}$ を求める最適化問題を解けば良い。

1. $d_r - \sum_{i \in S} a_{ri} y_i^{*}$ の最小値が0以上となってしまう場合、制約 (13)を破るような$d_r$、$a_{ri}$ が存在しないことを意味するので、手順2を行い$x_r$ を求めてアルゴリズムを終了する。
2. $d_r - \sum_{i \in S} a_{ri} y_i^{*}$ の最小値が0以下となる場合、制約 (13)を破るような$d_r$、$a_{ri}$ が存在したことを意味するので、この$d_r$、$a_{ri}$ を持つ新たなルート候補$r^*$を追加し、$R_0 \leftarrow R_0 \cup \{r^*\}$ とアップデートする。その後手順2へと戻る。

本論文の手法について

配送ルートの組を選択する上記の問題を、QUBOの変数の組合せを求める問題に読み替えると本記事で紹介する論文と同義となる。用語、変数の具体的な対応関係は

ブレインパッド記事	本論文
配送ルート$r$	QUBO変数$x^p = (x^p_1, x^p_2, …)$ の0-1 数列を指し示すインデックス$p$
配送ルート候補を選択する問題	様々なQUBO変数$x^p = (x^p_1, x^p_2, …)$ の0-1 数列の集合から、制約を満たしつつ最もコスト関数を小さくする数列$x^p$を選択する問題
$x_r$	$\lambda_p$
$y_i^{*}$	$\rho_k^{*}, \pi^{*}_0$
$d_r$、$a_{ri}$	QUBO変数$x^p = (x^p_1, x^p_2, …)$ の0-1 数列

といった関係になる。

問題設定

次のような、QUBOに制約条件を付加した問題を解きたい。

$$\begin{align} \min_{x} &\sum_{ij} Q_{ij} x_i x_j \\ \mathrm{s.t. }\ & \sum_{ij} A_{kij} x_i x_j \leq b_k \ \forall k \\ \end{align}$$

しかしながら、考慮するべきQUBO変数の数列の集合は膨大な数存在する。 (2^(バイナリ変数の個数))

ここで、上記の問題を次のように読み替える。

$$\begin{align} \min_{\lambda} &\sum_{p \in \mathcal{P}}\sum_{ij} Q_{ij} x^p_i x^p_j \lambda^p \\ \mathrm{s.t. }\ & \sum_{p\in \mathcal{P}}\sum_{ij} A_{kij} x^p_i x^p_j \lambda_p \leq b_k \ \forall k \\ \sum_{p} \lambda^p &= 1 \\ \lambda^p &= \{0, 1\} \ \forall p \in \mathcal{P} \end{align}$$

上式は、予めQUBO数列 $x^p = (x^p_1, x^p_2, \ldots)$ を $p$ 個用意したのち、この数列の中で (14)式のコスト関数を最も小さくし、かつ(15)式の制約を満たす数列を1つだけ選んでくる問題を指しているため、QUBO数列 $x^p_i = (x^p_1, x^p_2, \ldots)$ が全てのバイナリ変数の組合せを網羅するように用意されていれば (14), (15)式と等価な問題となる。

列生成法アルゴリズム

このセクションでは、上述の列生成法アルゴリズムに対応させる形で本論文の提案手法を記述する。

1. QUBO変数$x^p = (x^p_1, x^p_2, \ldots)$ の数列を適当な数用意し、その集合を$\mathcal{P}_0$ とする。

2. QUBO変数$x^p = (x^p_1, x^p_2, \ldots)$ の数列を1つだけ選んでくる問題の緩和問題である

   $$\begin{align} \min_{\lambda} &\sum_{p \in \mathcal{P}_0}\sum_{ij} Q_{ij} x^p_i x^p_j \lambda^p \\ \mathrm{s.t. }\ & \sum_{p\in \mathcal{P}_0}\sum_{ij} A_{kij} x^p_i x^p_j \lambda_p \leq b_k \ \forall k \\ \sum_{p} \lambda^p &= 1 \\ 0 &\leq \lambda^p\leq 1 \ \forall p \in \mathcal{P}_0 \end{align}$$

   を解く。これにより、双対問題である

   $$\begin{align} \max_{\rho, \pi_0} &\sum_{k} b_k \rho_k + \pi_0 \\ \mathrm{s.t. }\ & \sum_{k} \sum_{ij} A_{kij} x^p_i x^p_j \rho_k + \pi_0 \leq \sum_{ij} Q_{ij} x^p_i x^p_j \ \forall p \in \mathcal{P}_0\\ \end{align}$$

   の解 $\rho_k^{*}, \pi^{*}_0$ も同時に解けたことになる。

3. 双対問題の解 $\rho_k^{*}, \pi^{*}_0$ のもとで、制約式 (25) を破るようなQUBO変数の数列$x^p$を探す。すなわち、

$$\sum_{k} \sum_{ij} A_{kij} x^p_i x^p_j \rho_k + \pi_0 > \sum_{ij} Q_{ij} x^p_i x^p_j$$

を満たすQUBO変数の数列$x^p$を探す。これを達成するには論文中の式 (10)の最適化問題

$$\min_x \sum_{ij} \left\{ \left( Q_{ij} - \sum_{k=1}^m \rho_{k}^* A_{kij} \right) x_i x_j \right\} - \pi_0^*$$

を解けば良い。この式(10)に制約条件が存在しないことがミソとなる (QUBOソルバーで扱いやすい)

1. 論文中の式 (10)の最小値が0以上となってしまう場合、制約 (13)を破るようなQUBO変数の数列$x^p$が存在しないことを意味するので、手順2を行い$\lambda_p$ を求めてアルゴリズムを終了する。
2. 最小値が0以下となる場合、制約 (25)を破るようQUBO変数の数列$x^p$がが存在したことを意味するので、このQUBO変数の数列を$x^{p+1}$として定義し、$R_0 \leftarrow R_0 \cup \{r^*\}$ とアップデートする。その後手順2へと戻る。

まとめ

このブログではQUBO形式の問題に制約条件が含まれるような問題に対し、数理最適化の手法の一つである列生成法を応用する論文を紹介した。

今後も一層、ただQUBO問題を直接求解するだけではなく、数理最適化の技法をハイブリッドした手法が出てくるのではと予想している。

募集

最後に
＼Rustエンジニア・数理最適化エンジニア募集中！／
株式会社Jijでは、数学や物理学のバックグラウンドを活かし、量子計算と数理最適化のフロンティアで活躍するRustエンジニア、数理最適化エンジニアを募集しています！
詳細は下記のリンクからご覧ください。皆さんのご応募をお待ちしております！
Rustエンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/51062
数理最適化エンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/75132