背景

皆さんはどのような用途にRustを使用していますか?私は現職では物流オペレーションの最適化アルゴリズムを開発しており、人のスケジューリング、トラックのルート最適化、商品の梱包など離散的な構造の組合せ最適化を扱うことが多く、この分野にRustを使用しています。

機械学習ではPythonを利用することが多いですが組合せ最適化ではPythonはあまり活躍できません。Pythonが機械学習で多く利用されるのはNumPyやPyTorchなどを使用して行列計算が中心のアルゴリズムを実装するのが圧倒的に得意だからですが、木構造や順列などは行列演算とは異なるタイプの演算だからです。

以前はC++でアルゴリズムを実装し、PyBind11を用いてPythonから呼べる形にして全体のアプリケーションをPythonとC++で開発していましたが、開発環境の構築やデバック/プロファイリングがしにくいなどの問題がありました。RustはC++とほぼ同等の実行速度を誇り、洗練された文法とcargoという優れたビルドシステムによってまさにこう言った組合せ最適化のアルゴリズムとアプリケーションを通しで実装するのに最も適している言語だと思います。

今回は組合せ最適化問題を解く方法の中でも最も汎用的な局所探索(Local Search)とそれを行うための自作の汎用フレームワークであるlocalsearchについて紹介したいと思います。

組合せ最適化と局所探索の基礎

最適化問題は関数$f(x)$を最小(最大)にするような$x$を探すという問題です。以下の3つの要素を特定することで最適化問題を特定できます。

独立変数, $x$
目的関数, $f$、スコア、コストとも呼ぶ
制約条件, $x$の取りうる範囲や各次元間の制約

通常は制約条件を満たす$x$のことを実行可能解といいますが、この記事内では単純に解と呼びます。

連続最適化問題

例えば機械学習の訓練フェーズも最適化問題です。まず$y=\mathrm{model}(x|w)$のようにあるパラメトリックなモデル関数を用いて$y$を予測し、実際の観測値$y_{obs}$との差分$\mathrm{loss}(y|y_{obs})=\rm{loss}(w|x,y_{obs})$を最小にするようなパラメータ$w$を求める問題です。

独立変数 -> $w$、任意の実数ベクトル
目的関数 -> $\mathrm{loss}$
制約条件 -> 通常は無し

このような問題は制約なしの連続最適化と呼ばれ、勾配降下法を用いて解くことが一般的です。

離散最適化問題, 組合せ最適化問題

一方で組合せ最適化問題では独立変数$x$が連続な実数ベクトルではなく、順列や組合せなど複雑な制約条件を持った離散的な構造を扱います。

例

巡回セールスマン問題, TSP

独立変数: $x \in \mathrm{Perm}(k)$, 要素がkの順列の集合
目的関数: 合計移動距離、時間、コストなど
制約条件: 各都市の到達と出発の時間窓など(TSPTW)

ナップサック問題

$x_{ki} \in [0, 1]$: ナップサック$k$にアイテム$i$を入れるかどうかの指示変数
目的関数: 必要なナップサックの合計数量、ナップサックに入れられたアイテムの合計価値など
制約条件: ナップサックの最大数量、各ナップサックの合計容積、必要なアイテムの合計価値など

このような問題に対しては目的関数が微分できませんので勾配降下法を用いることができません。代わりに問題をグラフ構造で表して効率よく探索したり、現在の解を少し動かした近傍の中で最もいいものを逐次探索していくようなアルゴリズムがよく用いられます。

局所探索法

局所探索法(Local Search)は組合せ最適化問題を解く最も汎用的なフレームワークで以下のステップからなります。

初期解 $x_0$ を生成
目的関数で評価 $s_0 = f(x_0)$
現在の値を設定$(x, s) = (x_0, s_0)$
$x$からわずかに変化させた複数の試行解を生成し、その中で最もスコアが良いものを選択 -> $(x', s')$
一定のルールに従って$(x, s)$を$(x', s')$に更新
終了条件を満たすまで4,5を繰り返す
探索中に見つかった最もスコアが良い解を出力

4,5のステップを繰り返すことによって少しずつ解を改善していくという流れになっていることがわかります。特に4のステップは現在の解$x$から試行解を作成する手順のみあればいいわけですので特に$x$が実数値ベクトルである必要もありませんし、複雑な制約条件を満たすような生成法が構成できれば効率的に探索することが可能です。そういった意味で初めにも書いたように最も汎用的なフレームワークというわけです。

探索する変数の表現も最終的な解の表現$x$である必要はありません。解をある中間表現$z$から決定論的な処理で構成できるようにし、この中間表現を局所探索し、評価する前にまず中間表現$z$から解$x$決定論的な処理でデコードしてから評価するという手順も考えられます。特に前者の中間表現の空間のことを探索空間、最終的な解の空間を解空間と呼びます。例えば梅谷俊治先生の著書¹では長方形積み込み問題(Bin Packing問題)に対してこのアプローチを紹介しています。Bin Pakcing問題にはbottom-left-fill(BLF)など様々な決定論的なヒューリスティックアルゴリズムがあり、複数の長方形をどの順番で入力するかで最終的な結果が変わるのでこの順番を中間表現として局所探索で確率的に探索し、その時の順番でBLF法を使って最終的な配置とその評価値を求めるという手順になります。このアプローチはとても強力で私も実務で様々な問題に適用しています。²

5の更新ルールを調整することによって探索と活用(いわゆるExploration & Exploitation)のバランスを調整して収束の速さや局所最適解からの脱出を実現することができます。代表的なルールをいくつかあげます。

Hill Clibming, 山登り法: s'がsより改善した場合のみ更新する
Epsilon Greedy: 一定確率εでs'がsより改悪した場合も更新する
Metropolis: 改悪の場合も exp(-β(s'=s))の確率で更新する。分子シミュレーションやベイズの事後分布のサンプリングでおなじみ
Simulated Annealing: Metropolisにおいて(逆)温度パラメータを高いところから始めてだんだんと冷ましていく方法
Tabu Search: 以前の更新を繰り返さないように、その更新の逆の操作を禁止リストとして保存しておき、スコアは見ずに禁止リストに抵触するかどうかで更新を決める

ライブラリ化の方針

局所探索を適用していくうえで自分が解きたい問題依存の個所と汎用的な個所を分けることが重要です。前節のステップでいえば1,2,4つまり初期解の生成、試行解の生成、目的関数が問題依存のステップで、それ以外は共通で使える処理になります。つまり1,2,4をインターフェース(Rustではtrait)として抽象化しておき、主に5の部分をライブラリ内に直接実装すれば様々な問題に対して統一的に利用できるようになります。

自作の localsearch においては以下のようにtraitを定義しています。

pub trait OptModel: Sync + Send {
    /// Type of the Score
    type ScoreType: Ord + Copy + Sync + Send;
    /// Type of the Solution
    type SolutionType: Clone + Sync + Send;
    /// Type of the Transition
    type TransitionType: Clone + Sync + Send;

    /// Randomly generate a solution
    fn generate_random_solution<R: rand::Rng>(
        &self,
        rng: &mut R,
    ) -> AnyResult<(Self::SolutionType, Self::ScoreType)>;

    /// Generate a new trial solution from current solution
    fn generate_trial_solution<R: rand::Rng>(
        &self,
        current_solution: Self::SolutionType,
        current_score: Self::ScoreType,
        rng: &mut R,
    ) -> (Self::SolutionType, Self::TransitionType, Self::ScoreType);
}

ちょっとネーミングがいまいちですが、generate_random_solutionが初期解の生成で³ generate_trial_solutionが試行解の生成です。問題によっては解の生成と同時にスコアの評価も行うと高速化できるケースもあるのであえて解とスコアを同時に返すインターフェースにしています。trialのほうはスコアの差分を効率よく計算できるケースも考えてcurrent_scoreも引数として渡す設計です。TransitionはTabu Searchの場合のみ使う試行解の生成の逆操作に相当します。基本的にはこの２つのmethodを実装すれば前節で紹介したような様々なアルゴリズムで最適化を行うことができます。

なお、v0.22.0の段階で以下の複数の探索アルゴリズムを実装しています。

Hill Climbing.
Epsilon Greedy Search
Tabu Search
Metropolis
Simulated Annealing
Parallel Tempering
Population Annealing
Adaptive Annealing
Relative Annealing
Logistic Annealing
Tsallis Relative Annealing

特に9以降はオリジナルの方法です。何か面白い探索法のアイディアがありましたらぜひPull Requestやメッセージをいただければと思います。

使用例

例えば localsearch のexamples内ではTSPの問題モデルを実装してあります。

type Edge = (usize, usize);
type SolutionType = Vec<usize>;
// remvoed edges and inserted edges
type TransitionType = ([Edge; 2], [Edge; 2]);
type ScoreType = NotNan<f64>;

struct TSPModel {
    start: usize,
    distance_matrix: HashMap<Edge, f64>,
}

...

impl OptModel for TSPModel {
    type SolutionType = SolutionType;
    type TransitionType = TransitionType;
    type ScoreType = ScoreType;
    fn generate_random_solution<R: rand::Rng>(
        &self,
        rng: &mut R,
    ) -> AnyResult<(Self::SolutionType, Self::ScoreType)> {
        ...
    }

    fn generate_trial_solution<R: rand::Rng>(
        &self,
        current_solution: Self::SolutionType,
        current_score: Self::ScoreType,
        rng: &mut R,
    ) -> (Self::SolutionType, Self::TransitionType, Self::ScoreType) {
        ...
    }

解は訪れる都市のインデックスの順列 Vec<usize> でスコアは合計移動距離 NotNan<f64> です。generate_trial_solution内ではいわゆる2-optと呼ばれる操作で試行解を作成し、追加された辺と削除された辺を特定して距離の差を使ってスコアを更新するという方法で効率的に計算しています。

そして例えばこのモデルを以下のようにしてSimulated Annealingで最適化できます⁴。

# 全部で何ステップ行うか
let n_iter: usize = 100000;
# 計算時間の上限
let time_limit = Duration::from_secs(60);
# 収束判定、この回数たっても新しい最善解が見つからなかったら終了
let patience = n_iter / 2;
# 近傍を作る際にいくつ試行解を作るか
let n_trials = 16;
# この回数だけ新しい最善解が見つからなかったら探索を今の最善解からやり直し
let return_iter = n_iter / 50;
# 初期逆温度
let initial_beta = 1e-2;
# 冷却速度、逆温度で扱っているので実際には増やす方向
let cooling_rate = 1.01;
# 何ステップおきに逆温度を更新するか
let update_freq = 100;

# Optimizerインスタンスを生成
let optimizer = SimulatedAnnealingOptimizer::new(
    patience,
    n_trials,
    return_iter,
    initial_beta,
    cooling_rate,
    update_freq,
);

# Optimizerにtsp_modelと他パラメータを渡して実行
# callbackにはプログレスバーや最適解の定期ファイル保存などを実装できる(詳細はexampleのコードを見てください)
let (final_solution, final_score) = optimizer
    .run_with_callback(&tsp_model, None, n_iter, time_limit, &mut callback)
    .unwrap();

# 最終的な答えを表示
println!(
    "final score = {}, num of cities {}",
    final_score,
    final_solution.len()
);

以下はTSPLIBというTSPの問題集の中のberlin52という52都市の問題を解いた結果です。最適解の距離は7544.3659でそれを見つけることに成功しています。

このようにして実際に現場の業務に対して最適なオペレーションを求めるという際に、業務を如何に上記のOptModelとしてモデリングするかに集中すればいいわけです。

まとめ

RustはPythonやC++に代わって組合せ最適化問題のアプリケーションの開発にとても適した言語です
局所探索は様々な組合せ最適化問題に利用できるとても汎用的なフレームワークです
localsearch crateを使用することで実務の際には対象の問題を OptModel のtraitとして実装することに集中すればいいです
組合せ最適化問題は楽しいよ!

しっかり学ぶ数理最適化モデルからアルゴリズムまで ↩
実はこの考え方フレームワークは組合せ最適化問題を強化学習で解く際のフレームワークでもあります。 https://github.com/ai4co/rl4co ↩
初期解はランダムに決めることが多いのでこのような関数名にしていますが、将来的に変えるかもしれません。 ↩
この例では初期逆温度と冷却速度を定数で与えていますが、事前にサンプリングして目標とする初期受容確率から逆算するためのメソッドも実装してあります。] ↩