Zenn
🪬

数分割問題に対する差分法をRustで実装してみる

2023/12/16に公開
Rust
数理最適化
最適化
数分割
tech

この記事はOptimind Advent Calendar 2023の16日目の記事となります🪬

はじめに

先日社内の勉強会で数分割問題(partition problem)について学ぶ機会があり、そこでこの問題に対するシンプルなヒューリスティックである差分法(largest differencing method)について知ったのでRustの勉強がてら実装してみます。

数分割問題

数分割問題は与えられた正の整数からなるリストを、要素の総和が等しくなるような二つの集合に分割可能かを判定する問題であり、NP完全問題あることが知られています。

例1：\{6, 9, 19, 13, 17\}

\{19, 13\}と\{6, 9, 17\}と分割すれば両方の集合での要素の和が共に32となるため、この入力は分割可能です。

例2：\{3, 20, 13, 6, 30, 40, 73\}

この入力に対してはどのように分割してもそれぞれの集合の総和を一致させることは出来ないため（嘘だと思う方は是非試してみて下さい）、この入力は分割不可能と判定されます。

数分割問題は総和の差が最小になるような二つの集合への分割を見つける、という最適化問題バージョンを考えることも可能です。この最適化問題バージョンの場合、後者の例は\{73, 20\}と\{40, 6, 30, 3, 13\}と分割すれば総和の差は93 - 92 = 1となるためこれが最適解になります。（\{6, 13, 73\}と\{3, 20, 30, 40\}と分割しても差が1となるため、これも最適解になります。）
数分割問題の日常的な応用例としては、スポーツなどにおけるチーム分けやバッチ処理などを行う際の処理の負荷分散などがあります。（n\ge2の場合も含むより一般的な数分割問題はMultiway number partitioningと呼ばれます）

差分法

では本記事の本題の差分法です。差分法を擬似コードで書くと以下のようになります。

1 while リストに要素が2つ以上ある do
2     a <- リスト内の最大の要素
3     b <- リスト内の二番目に大きい要素
4     aとbをリストから削除し、a - bをリストに加える
5 return

要素の管理をheapで行えば、このアルゴリズムは要素数nに対してO(nlog(n))で動作します。最終的な分割においてaとbは異なる分割に入り、whileループ終了後にリストに残った要素が各分割の合計の差分を表します。
上で紹介した例1を使って差分法の流れを図示すると以下のようになります（分かりやすさのため要素は降順に並べ替えてあります）

最終的な解である分割は以下のように逆順に辿ることによって再構築出来ます。（子の内、大きい方の要素を親と同じ分割に入れる）

非常にシンプルなアルゴリズムですが、計算実験において（一部の比較的小さな入力を除き）アニーリング法を使ったヒューリスティックより優れた性能を示したことが報告されています。[1]（もちろん差分法もヒューリスティックであるため常に最適解を見つける保証はありません）

実装

まず愚直に上の擬似コードをRustで書いてみます。

use std::collections::BinaryHeap;

fn largest_differencing(numbers: Vec<u64>) -> u64 {
    assert!(!numbers.is_empty(), "numbers must be non-empty");

    let mut heap = BinaryHeap::from(numbers);
    while heap.len() > 1 {
        let a = heap.pop().unwrap();
        let b = heap.pop().unwrap();
        let diff = a - b;
        heap.push(diff);
    }

    heap.pop().unwrap()
}

fn main() {
    let x = vec![6, 9, 19, 13, 17];         // 最適値は0
    let y = vec![3, 20, 13, 6, 30, 40, 73]; // 最適値は1
    println!("{}", largest_differencing(x));
    println!("{}", largest_differencing(y));
}

実際にコードを実行してみると正しい答えが得られていることが分かります。

$ cargo run
0
1

しかしこれだと具体的な解（分割）を得ることが出来ないため、二つの要素を差分で置き換える操作をする際に、差分を親、二つの要素を子とするような二分木を作っておくことでアルゴリズム終了後に再起的に解を構築出来るようにします。ここでは書籍の実装[2]を参考に、以下のように列挙型で二分木を定義します。

enum BinaryTree {
    Empty,
    NonEmpty(Box<Node>),
}

struct Node {
    element: u64,
    left: BinaryTree,
    right: BinaryTree,
}

impl Node {
    fn new(element: u64, left: BinaryTree, right: BinaryTree) -> Self {
        Node {
            element,
            left,
            right,
        }
    }
}

これを使って最初に示した関数を書き換えると以下のようになります。（まだ解を構築する処理は入れてません）

fn largest_differencing(numbers: Vec<u64>) -> u64 {
    assert!(!numbers.is_empty(), "numbers must be non-empty");

    let numbers = numbers
        .into_iter()
        .map(|n| Node::new(n, BinaryTree::Empty, BinaryTree::Empty))
        .collect::<Vec<_>>();
    let mut heap = BinaryHeap::from(numbers);
    while heap.len() > 1 {
        let a = heap.pop().unwrap();
        let b = heap.pop().unwrap();
        let diff = Node::new(
            a.element - b.element,
            BinaryTree::NonEmpty(Box::new(a)), // 大きい方の要素を左の子にする
            BinaryTree::NonEmpty(Box::new(b)),
        );
        heap.push(diff);
    }

    heap.pop().unwrap().element
}

ここでひとまずcargo checkを実行すると、以下のように「NodeにOrd traitが実装されていない！」と怒られました。

error[E0277]: the trait bound `Node: Ord` is not satisfied
  --> src/main.rs:31:37
   |
31 |     let mut heap = BinaryHeap::from(numbers);
   |                    ---------------- ^^^^^^^ the trait `Ord` is not implemented for `Node`
   |                    |
   |                    required by a bound introduced by this call
   |

Nodeをheapで管理するためにはNode同士が比較可能でないとダメなのはぐう正論のため、NodeにOrd traitを実装します。またPartialOrd、PartialEq、Eqも必要になるためそちらもまとめて実装します。

impl PartialEq for Node {
    fn eq(&self, other: &Self) -> bool {
        self.element == other.element
    }
}

impl Eq for Node {}

impl PartialOrd for Node {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.element.cmp(&other.element))
    }
}

impl Ord for Node {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        self.element.cmp(&other.element)
    }
}

改めてcargo checkを実行するとエラーが出ず、コードがコンパイル出来る状態になりました！

次に、要素がどちらの分割に含まれているかを表すための列挙型を定義します。これは解を構築する際、大きい方の子を親と同じ分割にいれるための判定に使います。

enum Partition {
    Left,
    Right,
}

以上を踏まえて、差分法適用後の最終的なNodeから解を再構築する関数は以下のように書けます。

fn reconstruct_solution(node: Node, current_partition: Partition, left: &mut Vec<u64>, right: &mut Vec<u64>) {
    let large = node.left;
    let small = node.right;

    if let (BinaryTree::NonEmpty(large_node), BinaryTree::NonEmpty(small_node)) = (large, small) {
        match current_partition {
            Partition::Left => {
                let index = left.iter().position(|x| *x == node.element).unwrap();
                left.remove(index);
                left.push(large_node.element); // 大きい方の要素を同じ分割(left)に入れる
                right.push(small_node.element);
                reconstruct_solution(*large_node, Partition::Left, left, right);
                reconstruct_solution(*small_node, Partition::Right, left, right);
            }
            Partition::Right => {
                let index = right.iter().position(|x| *x == node.element).unwrap();
                right.remove(index);
                right.push(large_node.element); // 大きい方の要素を同じ分割(right)に入れる
                left.push(small_node.element);
                reconstruct_solution(*large_node, Partition::Right, left, right);
                reconstruct_solution(*small_node, Partition::Left, left, right);
            }
        }
    }
}

（本来は分割を表すleftとrightに対するデータ構造としてmultiset等を使うのが適していますが、Rustの標準ライブラリにそれに相当するものが無いためVectorで代用しています。）
最終的なコード全体は以下のようになります。

コード全体
main.rs
use std::collections::BinaryHeap;

enum Partition {
    Left,
    Right,
}

enum BinaryTree {
    Empty,
    NonEmpty(Box<Node>),
}

struct Node {
    element: u64,
    left: BinaryTree,
    right: BinaryTree,
}

impl Node {
    fn new(element: u64, left: BinaryTree, right: BinaryTree) -> Self {
        Node {
            element,
            left,
            right,
        }
    }
}

impl PartialEq for Node {
    fn eq(&self, other: &Self) -> bool {
        self.element == other.element
    }
}

impl Eq for Node {}

impl PartialOrd for Node {
    fn partial_cmp(&self, other: &Self) -> Option<std::cmp::Ordering> {
        Some(self.element.cmp(&other.element))
    }
}

impl Ord for Node {
    fn cmp(&self, other: &Self) -> std::cmp::Ordering {
        self.element.cmp(&other.element)
    }
}

type Solution = (Vec<u64>, Vec<u64>);

fn largest_differencing(numbers: Vec<u64>) -> Solution {
    assert!(!numbers.is_empty(), "numbers must be non-empty");

    let numbers = numbers
        .into_iter()
        .map(|n| Node::new(n, BinaryTree::Empty, BinaryTree::Empty))
        .collect::<Vec<_>>();
    let mut heap = BinaryHeap::from(numbers);
    while heap.len() > 1 {
        let a = heap.pop().unwrap();
        let b = heap.pop().unwrap();
        let diff = Node::new(
            a.element - b.element,
            BinaryTree::NonEmpty(Box::new(a)),　// 大きい方の要素を左の子にする
            BinaryTree::NonEmpty(Box::new(b)),
        );
        heap.push(diff);
    }

    let root = heap.pop().unwrap();
    let mut left = vec![];
    let mut right = vec![];
    left.push(root.element);
    reconstruct_solution(root, Partition::Left, &mut left, &mut right);

    (left, right)
}

fn reconstruct_solution(node: Node, current_partition: Partition, left: &mut Vec<u64>, right: &mut Vec<u64>) {
    let large = node.left;
    let small = node.right;

    if let (BinaryTree::NonEmpty(large_node), BinaryTree::NonEmpty(small_node)) = (large, small) {
        match current_partition {
            Partition::Left => {
                let index = left.iter().position(|x| *x == node.element).unwrap();
                left.remove(index);
                left.push(large_node.element); // 大きい方の要素を同じ分割(left)に入れる
                right.push(small_node.element);
                reconstruct_solution(*large_node, Partition::Left, left, right);
                reconstruct_solution(*small_node, Partition::Right, left, right);
            }
            Partition::Right => {
                let index = right.iter().position(|x| *x == node.element).unwrap();
                right.remove(index);
                right.push(large_node.element); // 大きい方の要素を同じ分割(right)に入れる
                left.push(small_node.element);
                reconstruct_solution(*large_node, Partition::Right, left, right);
                reconstruct_solution(*small_node, Partition::Left, left, right);
            }
        }
    }
}

fn main() {
    let x = vec![6, 9, 19, 13, 17]; // 最適値は0
    let y = vec![3, 20, 13, 6, 30, 40, 73]; // 最適値は1
    println!("{:?}", largest_differencing(x));
    println!("{:?}", largest_differencing(y));
}

実行すると実際の分割が得られていることが確認できました！

$ cargo run
([19, 13], [17, 6, 9])         # 32 - 32 = 0
([20, 73], [6, 13, 3, 30, 40]) # 93 - 92 = 1

参考文献

本記事を書くにあたって以下の論文や書籍を参考にさせて頂きました🙇‍♂️

脚注

  1. Optimization by simulated annealing: an experimental evaluation; part II, graph coloring and number partitioning ↩︎

  2. Programming Rust 2nd edition ↩︎

Discussion