Rust言語を使ってアニーリングで制約問題を解く

量子アニーリングとは?

今年の量子コンピューター Advent Calendar 2020には執筆時点でまだアニーリングの平易な解説がのっていないようなので、簡単に説明します。

なお、ここで紹介している用語の説明等は、文献によっては異なる意味で用いられている場合もあるのでご注意ください。

組合せ最適化問題

(現在のアニーリングマシンにおける)アニーリングは、二値変数の組$\mathbf{x}$を引数にとる高々2次の実数値関数$f(\mathbf{x})$の最適化問題、すなわち以下のように表せる解$\mathbf{x}$を求める問題を解く方法の1つです。

\mathbf{x} = \underset{\mathbf{x}}{\mathrm{argmin}}  f(\mathbf{x})

このような問題は組合せ最適化問題とも呼ばれます。また、ここで問題を定義している$f(\mathbf{x})$を目的関数と呼び、これを以下のように数式で表現したものをBinary Quadric Model(BQM)やQuadratic Unconstrained Binary Optimization(QUBO)と呼びます。

\begin{align}
f(\mathbf{x})&= \sum_{i, j=1}^N{Q_{ij}q_iq_j}\\
&= \sum_{i=1}^N{h_is_i}+\sum_{i\neq j}{J_{ij}s_is_j}
\end{align}

ここで、$q_i\in {1,0}$をバイナリ変数, $s_i\in{1,-1}$をスピン変数と呼びます。(以下、これらを合わせて$x_i$と表記します。)さらに、実数値係数$Q_{ij}$もまたQUBOと呼ばれます。$h_i$や$J_{ij}$はそれぞれバイアスやカップリングなどと呼ばれます。

制約の導入

このように、アニーリングを用いると組合せ最適問題を解くことができますが、それだけでなく、制約を含めた制約最適化問題(Constrained Optimization Problem, COP)を解くこともできます。これは、組合せ最適化問題を表す$f_{\mathrm{optim}}(\mathbf{x})$に対して、$n$番目の制約$C_n$を表す$f_{C_n}(\mathbf{x})$を加えたものを目的関数とすることで実現されます。

f(\mathbf{x})=f_{\mathrm{optim}}(\mathbf{x})+\sum_n{\lambda_nf_{C_n}}(\mathbf{x})

ここで、$\lambda_n>0$は制約$C_n$の強さを表すパラメータで、プレースホルダ等と呼ばれます。また関数$f_{C_n}(\mathbf{x})>0$は、$f_{C_n}(\mathbf{x})=0$を満たす二値変数が制約$\mathbf{x}$が制約$C_n$を満たすように定義されます。

ここで、$\lambda_n$を適切に設定することで、$f(\mathbf{x})$の最適解である$\underset{\mathbf{x}}{\mathrm{argmin}} f(\mathbf{x})$がすべての制約$C_n$を満たすようにすることができます。しかしながら、実際にこの組合せ最適化問題をアニーリングマシンで解くと、局所解という解が得られることがあります。局所解$\mathbf{x}$はすべての制約$C_n$を満たすことが保証されないため、アニーリングを行うたびに結果$\mathbf{x}$が制約${C_n}$を満たすことを確認し、満たされない場合は${\lambda_n}$を変更してアニーリングをやり直す、等をする必要があります。

次元削減

アニーリングによって直接扱える組合せ最適化問題は(今の所)高々2次までです。そのため、一般の$N$次の組合せ最適化問題(HUBOという)を解く場合は次元削減(グラフカット最適化ともいう)を行う必要があります。この方法にはいくつかのやり方がありますが、例として、(石川の方法)によればバイナリ変数の3次の項$tq_1q_2q_3$に対して、新たにアンシラ変数$q_0$を追加し

• $t>0$のとき、$\underset{q_0}{\mathrm{min}}\ t{q_0(q_1+q_2+q_3-1)+(q_1q_2+q_2q_3+q_1q_3)-(q_1+q_2+q_3)+1}$
• $t<0$のとき、$\underset{q_0}{\mathrm{min}}\ tq_0(q_1+q_2+q_3-2)$

に置き換えることで2次の式に変換することができます。

Rustを用いて巡回セールスマン問題を解く

ここでは、RustQUBO¹を用いて巡回セールスマン問題を解きます²。

Rustのインストール

Rustupをインストールし、rustのツールチェインを導入します。実際に使うコマンドはcargoです。

Rust をインストールを参照

プロジェクトの作成

適当なディレクトリを作成し、cargo initでプロジェクトを作成します。

$ mkdir tsp-test
$ cd tsp-test
$ cargo init

Cargo.tomlの編集

RustQUBOを使うため、以下の行をCargo.tomlに追加します。

Cargo.toml

# ... (省略) ...
[dependencies]
rustqubo={ git = "https://github.com/SpoonQ/rustqubo", rev = "e6a0e", branch = "main" }
# ↑この行を追加

cargo checkを実行して、正しくインストールされることを確認してください。

プログラムの作成

src/main.rsを編集します。

extern crate rustqubo;
use rustqubo::expr::Expr;
use rustqubo::solve::SimpleSolver;

#[derive(PartialEq, Eq, Copy, Clone, Debug, Hash, PartialOrd, Ord)]
struct TspQubit(usize, usize);

fn main() {
    let cities = 5usize;
    let distance = [
        [0.0, 5.0, 5.0, 3.0, 4.5],
        [5.0, 0.0, 3.5, 5.0, 7.0],
        [5.0, 3.5, 0.0, 3.0, 4.5],
        [3.0, 5.0, 3.0, 0.0, 2.5],
        [4.5, 7.0, 4.5, 2.5, 0.0],
    ];
    let hmlt_city = (0..cities).into_iter().fold(Expr::Number(0), |exp, c| {
        let inner = (0..cities)
            .into_iter()
            .fold(Expr::Number(-1), |e, o| e + Expr::Binary(TspQubit(c, o)));
        exp + Expr::Constraint {
            label: format!("city {:}", c),
            expr: Box::new(inner.clone() * inner),
        }
    });
    let hmlt_order = (0..cities).into_iter().fold(Expr::Number(0), |exp, o| {
        let inner = (0..cities)
            .into_iter()
            .fold(Expr::Number(-1), |e, c| e + Expr::Binary(TspQubit(c, o)));
        exp + Expr::Constraint {
            label: format!("order {:}", o),
            expr: Box::new(inner.clone() * inner),
        }
    });
    let mut hmlt_distance = Expr::Number(0);
    for i in (0..cities).into_iter() {
        for j in (0..cities).into_iter() {
            for k in (0..cities).into_iter() {
                let d_ij = Expr::Float(distance[i][j]);
                hmlt_distance = hmlt_distance
                    + d_ij
                        * Expr::Binary(TspQubit(i, k))
                        * Expr::Binary(TspQubit(j, (k + 1) % cities))
            }
        }
    }
    let hmlt: Expr<(), _, _> = Expr::Number(700) * (hmlt_city + hmlt_order) + hmlt_distance;
    let compiled = hmlt.compile();
    let mut solver = SimpleSolver::new(&compiled);
    solver.generations = 20;
    solver.beta_count = 50;
    solver.sweeps_per_beta = 1;
    solver.samples = 1;
    let (c, qubits, constraints) = solver.solve();
    assert!(constraints.len() == 0);
    println!("{:?}", &qubits);
}

実行

以下のように実行します。

$ cargo run
{TspQubit(1, 1): false, TspQubit(0, 2): false, TspQubit(1, 3): false, TspQubit(1, 2): false, TspQubit(0, 3): false, TspQubit(3, 0): false, TspQubit(2, 4): false, TspQubit(4, 1): true, TspQubit(2, 2): true, TspQubit(0, 0): true, TspQubit(2, 3): false, TspQubit(0, 1): false, TspQubit(4, 4): false, TspQubit(1, 4): true, TspQubit(4, 3): false, TspQubit(4, 2): false, TspQubit(3, 2): false, TspQubit(3, 4): false, TspQubit(2, 0): false, TspQubit(4, 0): false, TspQubit(2, 1): false, TspQubit(3, 1): false, TspQubit(0, 4): false, TspQubit(1, 0): false, TspQubit(3, 3): true}

この解から、$0, 4, 2, 3, 1, 0$の順でめぐるルートが最適であることがわかります。

1. 他に有名な実装としてPyQUBOがありますが、RustQUBOとは若干特徴が異なっています。 ↩

2. 量子アニーリングではなく、RustQUBOに内蔵されているSAを用いています。 ↩