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

# 量子アニーリングとは?

今年の[量子コンピューター Advent Calendar 2020](https://qiita.com/advent-calendar/2020/quantum)には執筆時点でまだアニーリングの平易な解説がのっていないようなので、簡単に説明します。

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

## 組合せ最適化問題

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

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

このような問題は[**組合せ最適化問題**](https://ja.wikipedia.org/wiki/%E7%B5%84%E5%90%88%E3%81%9B%E6%9C%80%E9%81%A9%E5%8C%96)とも呼ばれます。また、ここで問題を定義している$f(\mathbf{x})$を**目的関数**と呼び、これを以下のように数式で表現したものをBinary Quadric Model(BQM)やQuadratic Unconstrained Binary Optimization(QUBO)と呼びます。

```math
\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})$を加えたものを目的関数とすることで実現されます。

```math
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`](https://github.com/SpoonQ/rustqubo)[^1]を用いて[巡回セールスマン問題](https://ja.wikipedia.org/wiki/%E5%B7%A1%E5%9B%9E%E3%82%BB%E3%83%BC%E3%83%AB%E3%82%B9%E3%83%9E%E3%83%B3%E5%95%8F%E9%A1%8C)を解きます[^2]。

## Rustのインストール

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

[Rust をインストール](https://www.rust-lang.org/ja/tools/install)を参照

## プロジェクトの作成

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

```bash
$ 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`を編集します。

```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);
}
```

## 実行

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

```bash
$ 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](https://github.com/recruit-communications/pyqubo)がありますが、`RustQUBO`とは若干特徴が異なっています。
[^2]: 量子アニーリングではなく、`RustQUBO`に内蔵されているSAを用いています。