Rustで差分進化を実装

# はじめに
差分進化(DE: Differential Evolution)とはStornとPriceが1997年に提案した最適化アルゴリズムである[^1]。
簡単に言えば、入力`X`に対して何らかの出力`Y`を返す評価関数`Y=F(X)`を使って`X`を最適化する手法である。
このとき、評価関数`F`はなんでもよく、**ブラックボックス**[^2]であっても良い[^3]。
(そのため、しばしばBlack Box Optimizationとも呼ばれる）
評価関数`F`の目的は、入力された`X`に対して結果`Y`を出力し、得られた結果`Y`の良し悪しを定量的に評価できれば良い。
差分進化の欠点としては、得られた結果が大域的最適解となることを保証していないことが挙げられる。

[^1]: R. Storn and K. Price、 "Differential Evolution-A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces," J. Global Optim., vol. 11, no. 4, pp. 341-359, Dec. 1997
[^2]: 内部の処理がわからなくとも、ある入力`X`に対しては常に同じ結果`Y`が返ってくる性質（参照透過性?）を満たしていれば
[^3]: すなわち、勾配法などとは異なり、対象とする関数の微分可能性を考えずに使える。

差分進化の構成要素としては、

- `N`個の候補`X`を作成し、それらを評価する(**初期化**)。
- 全ての候補`x`∈`X`に対して：
    - 候補`x`に対して他の候補と組み合わせた変異体`m`を作成し、評価する(**評価**）。
    - 生成した変異体`m`が`x`よりも優れた特性を示すのであれば`x`の中身を`m`で更新する(**更新**)。
- 設定した停止条件を満たすまで(**評価**)→(**更新**)を繰り返していく。

本稿では、差分進化を構成する関数として `initialization()`、`evaluate()`、そして`mutation()`を実装する。

なお、ソースコード中のコメントはこの記事を執筆するにあたって追記したものであり、
実際には最低限の英語のコメントに留めている。

# 問題設定
候補`x`は多項式形式の確率質量関数(PMF: Probability Mass Function)で表現された分布であるとする。

```rust
// 今回用いる理論解析の都合上、ELEMとKMAXをそれぞれ定義する
pub const ELEM: usize = 20; //配列の要素数 (> KMAX)
pub const KMAX: usize = 16; //多項式の最大次数

#[derive(Copy、 Clone)]
pub struct Pmf {
    v: [f64; ELEM], // PMFを定義する多項式
    g: f64,         // tを達成するときのg
    t: f64,         // 評価する値(特性)
}
```

候補`x`を作成するにはarrayを与える[^4]。

```rust
impl Pmf {
    fn new(pmf: [f64;ELEM]) -> Pmf {
        // PMFなので全部足すと1になるかチェック
        let sum: f64 = pmf.iter().sum();
        if sum > 0.0 {
            assert_approx_eq!(sum, 1.0);
        }

        // `struct Pmf`を返す
        Pmf { v: pmf, g: 0.0, t: 0.0 }
    }
}
```

ただ、ほぼ`Pmf::new([0.0;ELEM]);`とすることが多かったので、
結果的に下記のようになった。

```rust
impl Pmf {
    fn new() -> Pmf {
        // `struct Pmf`を返す
        Pmf { v: [0.0;ELEM], g: 0.0, t: 0.0 }
    }
}
```


[^4]: `Vec<>`にしなかったのはマルチスレッド化を見越して。サイズがわかってるものしか受け付けてくれない`Arc`を使いたいから。今回の場合だと別に`Vec<>`でも良い。
評価関数`F`に`x`を入力すると結果`t`が返ってくる。
今回は、結果`t`は高いほど良いとし、
もし`x.t`より生成した候補の`m.t`が高かったら`x`の中身を`m`で書き換える。

# 差分進化の実装

全候補数を`NUMBER_OF_CANDS`とする。
また、停止条件を全ての候補が更新されなかったときとする。

全体の処理としてはこんな感じになる。

```rust
pub fn run() {
    // Mutant作成時に自分以外に3つの候補が必要なため
    assert!(NUMBER_OF_CANDS > 4);
    assert!(ELEM > KMAX);

    // 乱数生成器の初期化。同じプログラム内で使い回す
    let mut rng = rand::thread_rng();
    
    // 初期化
    let mut candidates = initialization(&mut rng);

    // 変異体生成に使うVec
    let range: Vec<usize> = (0..NUMBER_OF_CANDS).collect();
    // 更新
    while update(&range, &mut rng, &mut candidates) > 0 {};
    
    // 最も良い特性を示す候補を表示
    search_best(&candidates).print();
}
```

## 初期化

`Pmf`を`NUMBER_OF_CANDS`個格納するarrayを作成し返す。

```rust
// 初期化
fn initialization(rng: &mut rand::rngs::ThreadRng) -> [Pmf; NUMBER_OF_CANDS] {
    let mut ret = [Pmf::new(); NUMBER_OF_CANDS];
    for idx in 0..NUMBER_OF_CANDS {
        // 次数分布の作成
        for d in 2..KMAX + 1 {
            ret[idx].v[d] = rng.gen();
        }
        // 分布なので全ての係数を足して1となるように正規化する必要がある
        ret[idx].normalize();
        
        // 生成した候補の特性を評価
        evaluate(&mut ret[idx]);
    }
    return ret;
}

// 評価
fn evaluate(mut pmf: &mut Pmf) {
    // 内部でPmf.gとPmf.tをPmf.vが達成する特性に更新
    // 与えられた候補に対して評価する何らかの関数
    some_theoretical_analysis::run_pmf(&mut pmf);
}
```

全ての係数を足して`1.0`となるように`Pmf.v`を正規化するメソッド。

```rust
impl Pmf {
    ...
    fn normalize(&mut self) {
        let sum: f64 = self.v.iter().sum();
        for x in 0..ELEM {
            self.v[x] /= sum;
        }
    }
}
```

## 更新

```rust
fn update(
    range: &Vec<usize>、
    mut rng: &mut rand::rngs::ThreadRng、
    candidates: &mut [Pmf; NUMBER_OF_CANDS]、
) -> u64 {
    // 更新回数をカウントする変数
    let mut swap = 0;

    for idx in 0..NUMBER_OF_CANDS {
        // 変異体の生成
        let mutant = generate_mutant(&idx, &range, &mut rng, &candidates);

        // 変異体の特性が元の候補より優れていたら更新
        if mutant.t > candidates[idx].t {
            swap += 1;
            candidates[idx].copy_from(&mutant);
        }
    }
    // 更新された候補数を返す
    return swap;
}
```

`update()`が更新された候補の数を返すようにすることで、

```rust
while update(&range, &mut rng, &mut candidates) > 0 {};
```
とし、繰り返しの停止条件を「全ての候補が更新がされなくなるまで」とすることができる。
もちろん、

```rust
while update(&range, &mut rng, &mut candidates) < 5 {};
```
とすることで更新が5未満だったら停止するようにしたり、

```rust
for _ in 0..NUMBER_OF_ITERATION {
    // 更新がないのに続けても無意味なのでこの停止条件を入れる
    if update(&range, &mut rng, &mut candidates) == 0 {
         break;
    }
}
```
とすることで規定の回数だけ更新するなど、
停止条件は用途に合わせて設定することができる。
`.copy_from()`は自分の中身を別の`Pmf`の中身で更新するメソッド。

```rust
impl Pmf {
    ...
    fn copy_from(&mut self, p: &Pmf) {
        for x in 0..ELEM {
            self.v[x] = p.v[x];
        }
        self.g = p.g;
        self.t = p.t;
    }
}
```

### 変異体の生成
候補`x`の変異体`m`の生成には、まず対象となる候補以外の3つの候補の係数を用いて中間体`s`と呼ばれるものを作る。
$$
s = \sum_i \{a_i - \mu(b_i + c_i)\}x^i
$$
ここで、$a_i$は候補`a`のPMFの`i`番目の次数の係数を表す。
$\mu$は `[0、1)`の値をとる**スケーリング係数(SF: Scaling Factor)**と呼称されるパラメータ。
次に生成した`s`と**交叉率(CR: Crossover Rate)**を用いて変異体`m`を作成する。
変異体`m`の各係数$m_i$は、
$$
m_i =
\begin{cases}
s_i&(r > CR)\newline
x_i&(\mathrm{otherwise})
\end{cases}
$$
となる。ここで、$r$は`[0、1)`の値をとる乱数である。

下記のソースコードを見た方がわかりやすいかもしれない。

```rust

const CR: f64 = 0.5;  // Crossover Rate
const SF: f64 = 0.8;  // Scaling Factor

fn generate_mutant(
    idx: &usize、
    range: &Vec<usize>、
    mut rng: &mut rand::rngs::ThreadRng、
    candidates: &[Pmf; NUMBER_OF_CANDS]、
) -> Pmf {
    // 自分(idx)以外のインデックスを3つ用意
    let (a, b, c) = get_three_indices(idx, &range, &mut rng);
    // 変異体のPmfを用意
    let mut mutant = Pmf::new();

    // 変異体の各係数を決める
    for x in 0..KMAX + 1 {
        if rng.gen::<f64>() > CR {
            mutant.v[x] =
                (candidates[a].v[x] + SF * (candidates[b].v[x] - candidates[c].v[x])).abs();
        } else {
            mutant.v[x] = candidates[*idx].v[x];
        }
    }
    // 分布なので正規化
    mutant.normalize();
    // 作成した変異体を評価
    evaluate(&mut mutant);
    return mutant;
}
```

### 補助関数
自分以外のインデックスを3つ返す関数。
`0`から`NUMBER_OF_CANDS-1`まで値を取る`Vec`からランダムに3つ選んで、
その中に自分のインデックスがなかったら返すようにする。
事前に`.filter()`で`idx`を`range`から弾くよりも`.clone()`が一回少なく済む。

```rust
fn get_three_indices(
    idx: &usize,
    range: &Vec<usize>,
    mut rng: &mut rand::rngs::ThreadRng,
) -> (usize, usize, usize) {
    loop {
　　　　　// (0..NUMBER_OF_CANDS)で初期化されたVecから要素を3つ取り出す
        let indices: Vec<usize> = range.choose_multiple(&mut rng, 3).cloned().collect();

        // そこに元のインデックスが含まれていなければloopを抜けてreturn
        if !indices.contains(&idx) {
            return (indices[0], indices[1], indices[2]);
        }
    }
}
```

## 最後に

特性が最も良かった候補を返す関数を定義する。

```rust
fn search_best(candidates: &[Pmf; NUMBER_OF_CANDS]) -> Pmf{
    let mut best = Pmf::new();
    for x in 0..NUMBER_OF_CANDS {
        if candidates[x].t > best.t {
            best.copy_from(&candidates[x]);
        }
    }
    return Pmf;
}

fn run() {
   ...
   search_best(&candidates).print();
}
```

`.print()`は出力用のメソッド。

```rust
impl Pmf {
    ...
    fn print(&self) {
        print!("[ ");
        for x in self.v.iter() {
            if x > &0.0 {
                print!("{:.8} ", x);
            }
        }
        println!("] {:.5} {:.5}", self.g, self.t);
    }
}
```

# まとめ
今回は`Rust`で差分進化の実装を行った。
CRやSFを可変とする様々な亜種[^5]が提案されているのでそちらも試してみるのも良いかも。

[^5]: [参考](https://www.jstage.jst.go.jp/article/jceeek/2014/0/2014_541/_pdf/-char/ja)
# 参考

- [原著](https://link.springer.com/article/10.1023%2FA%3A1008202821328?LI=true&utm_medium=affiliate&utm_source=commission_junction&utm_campaign=3_nsn6445_brand_PID4003003&utm_content=de_textlink)
- [A tutorial on Differential Evolution with Python](https://pablormier.github.io/2017/09/05/a-tutorial-on-differential-evolution-with-python/)
- [差分進化 (Wikipedia)](https://ja.wikipedia.org/wiki/差分進化)