LoginSignup
3
3

More than 3 years have passed since last update.

@Scstechr

Rustで差分進化を実装

Last updated at Posted at 2020-02-09

はじめに

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

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

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

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

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

問題設定

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

// 今回用いる理論解析の都合上、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。

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]);とすることが多かったので、
結果的に下記のようになった。

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

評価関数Fにxを入力すると結果tが返ってくる。
今回は、結果tは高いほど良いとし、
もしx.tより生成した候補のm.tが高かったらxの中身をmで書き換える。

差分進化の実装

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

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

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を作成し返す。

// 初期化
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を正規化するメソッド。

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

更新

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()が更新された候補の数を返すようにすることで、

while update(&range, &mut rng, &mut candidates) > 0 {};

とし、繰り返しの停止条件を「全ての候補が更新がされなくなるまで」とすることができる。
もちろん、

while update(&range, &mut rng, &mut candidates) < 5 {};

とすることで更新が5未満だったら停止するようにしたり、

for _ in 0..NUMBER_OF_ITERATION {
    // 更新がないのに続けても無意味なのでこの停止条件を入れる
    if update(&range, &mut rng, &mut candidates) == 0 {
         break;
    }
}

とすることで規定の回数だけ更新するなど、
停止条件は用途に合わせて設定することができる。
.copy_from()は自分の中身を別のPmfの中身で更新するメソッド。

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)の値をとる乱数である。

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


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()が一回少なく済む。

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

最後に

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

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()は出力用のメソッド。

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が提案されているのでそちらも試してみるのも良いかも。

参考

  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. すなわち、勾配法などとは異なり、対象とする関数の微分可能性を考えずに使える。 ↩

  4. Vec<>にしなかったのはマルチスレッド化を見越して。サイズがわかってるものしか受け付けてくれないArcを使いたいから。今回の場合だと別にVec<>でも良い。 ↩

  5. 参考 ↩

3
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
3