ネコ探索アルゴリズムについて
こんにちは。今日は強化学習というかメタヒューリスティックな方法にであるネコ探索を用いたナップサック問題の近似値を求める方法について、論文を読み実装したので書いていきたいと思います。このようなコウモリ探索、オオカミ探索などメタヒューリスティックなアルゴリズムについて日本語の記事が殆どなかったため、記しておこうという理由です。具体的には以下の論文を読みました。 また本ブログでは実装例は載せますが、個人で行なったものであり必ずしも正しいかはわかりません。あくまで参考にして詳しくは論文を読んでいただけると幸いです。
Ricardo Soto , Broderick Crawford , Angelo Aste Toledo, Hanns de la Fuente-Mella , Carlos Castro , Fernando Paredes , and Rodrigo Olivares (2019) "Solving the Manufacturing Cell Design Problem through Binary Cat Swarm Optimization with Dynamic Mixture Ratios
概要
ネコ探索とは何かということを説明しますが、正直名前の由来はいまいちしっくりこないです。ただネコ科の習性を活かした探索方法ということがわかります。具体的に以下の2つのmodeに従って探索されていきます。
seeking_mode
tracing_mode
順に説明します。
seeking_mode
ネコ科の動物(ライオンなど)は基本的にほとんどの時間はじっとしております。しかしその時間は必ずしも休んでいるわけではなく、敵がどこにいるか、どの方角にいるかなど、見張っている時間も多いです。 また遺伝的アルゴリズム同様、子孫を増やしていきます。これを元にした探索手法です。正直初めて論文を読んだ時そして今もよくわからないです。がアルゴリズムはしっかりしているのでその部分をナップサック問題に合わせて簡潔に説明します。 まずナップサックにどの荷物を入れるかという1つの解が存在します。 その解は配列として用意し、0101001のようにi番目の荷物を入れる場合は1を入れない場合は0を立てる配列です。 これを数十個ほどコピーします。そしてそれらのコピーされた配列1つ1つに対して、ある確率のもとで1->0,0->1に反転していきます。すると異なった配列が複数できるわけです。当然これらは価値が異なってきますので、容量オーバーとなるもの以外で、価値を計算していきます。このなかで良い価値が生じたものを 次の世代の解として更新し、残していくという流れです。(この部分では必ずしも最も価値が高いものを選ぶわけではなくそれも確率的に決めるのですが説明するより論文を読んでいただけたほうがよいため割愛します)
tracing_mode
ここでは実際にネコ科の動物が獲物を狙いに行くことを模倣したmodeです。 正直意味はわからないですが,ありがたいことに論文にフローチャートが掲載されているので、それどおり実装するだけです。ですので概要的な説明になりますが、速度ベクトルが2つあり、それぞれ0->1にする速度、1->0にする速度2つがあります。それをこれまでの最高スコアを出している配列をもとに決めていきます。そしてそのベクトルの大きさがランダム値以上の時はじめの初期配列を更新しそうでない時は更新しないという流れです。僕は最近知ったのですがPSO(粒子群最適化法)と呼ばれる手法に非常に近い形で配列を更新していくことになります PSOについては以下の記事が非常に参考になりました。
ざっくりいえば値を更新する際、 もともとの慣性速度とこれまでの最大のスコアを出している位置とのベクトル和により速度を更新し初期解を変更しようという流れです。また今回、更新された配列がナップサックの容量を超えた時はそもそも更新しないという方法を取りましたが、論文では別の問題ではありますが、よりよいアルゴリズムが記されていたのでこれをナップサック問題に適応できる形に直して適応させるべきと考えましたがよくわかりませんでしたのでわかる方がいましたら教えてくれると幸いです。
まとめ
以上の方法を元に実装したところ、ある程度の大きさまでは厳密解と一致しました。 またその探索の流れを見てみると少しずつ解に近づいていく様子が観察できたため、実装に大きなバグはないのかなと思います。ただnを大きくするとやや厳密解とずれるため、そのへんはパラメータ調整でうまく行くのかそれともバグがあるのかはわからないです。(AtCoder上のナップサック問題に適応してみると半分くらいACでした)今回英語の比較的新しい論文を読んで実装した経験はほぼ初だったのでブログにまとめました。何度も述べますが、実装に誤りがあるかもわからないので、論文を読んで各人が実装していただければと思います。
実装
標準入力は以下のような形です。
N Max_Capacity W_i,...,W_N (重さ) C_i,...,C_N (価値)
C++による実装
#include <bits/stdc++.h> #include<time.h> #include <unistd.h> using namespace std; //vectorの中身を表示する関数 template <class T> void show(vector<T> v) { for (int i = 0; i < v.size(); i++) { cerr << v[i] << " "; } cerr << "\n"; } //容量オーバーしていないかチェックする関数 bool check(vector<int> &x, vector<int> &weight, int capacity, int n) { int weight_sum = 0; for (int i = 0; i < n; i++) { if (x[i]) { weight_sum += weight[i]; } } if (weight_sum > capacity) { return false; } else return true; } //ナップサックに入った価値を計算する関数 int calc_sum_profit(vector<int> &x, vector<int> &profit,int n) { int profit_sum = 0; for (int i = 0; i < n; i++) { if (x[i]) { profit_sum += profit[i]; } } return profit_sum; } //価値最大を計算する関数 int calc_max(int C,int n,vector<int> &weight, vector<int>&profit,vector<int>&x_best ,vector<vector<int>>& cat,int capacity,int max_profit){ for (int i = 0; i < C; i++) { if (check(cat[i], weight, capacity, n)) { int sum1 = calc_sum_profit(cat[i], profit, n); if (sum1 > max_profit) { max_profit = sum1; for (int j = 0; j < n; j++) { x_best[j] = cat[i][j]; } } } } return max_profit; } //グローバル変数にしておく double MR = 0.75, SMP = 15, CDC = 0.2, PMO = 0.76; //seeking _mode void seeking_mode(vector<int> & x, vector<int> &profit,vector<int>& weight,int capacity,int n){ vector<vector<int>> copy_cat(SMP, vector<int>(n, 0)); //SMP個分コピーする for (int i = 0; i < SMP; i++) { for (int j = 0; j < n; j++) { copy_cat[i][j] = x[j]; } } for (int i = 0; i < SMP; i++) { for (int j = 0; j < n; j++) { if((double)rand()/RAND_MAX < CDC){ if ((double)rand() / RAND_MAX < PMO){ copy_cat[i][j] = rand() % 2; } } } } vector<int> ok_cat; vector<int> fs(SMP); int fs_max = 0, fs_min = 1e9; int flag=0; for (int i = 0; i < SMP; i++) { if(check(copy_cat[i],weight,capacity,n)){ int tmp = calc_sum_profit(copy_cat[i], profit, n); fs_max = max(fs_max, tmp); fs_min = min(fs_min, tmp); fs[i] = tmp; flag = 1; } } //もしすべて容量オーバーしていたら何もしない if(!flag) return; //確率の // 重み付け確率を求めたい vector<pair<double,int>> prob; for (int i = 0; i < SMP; i++) { if (check(copy_cat[i], weight, capacity, n)) prob.push_back({ (double)((fs[i] - fs_min) / (double)(fs_max - fs_min+0.001)),i}); } sort(prob.begin(), prob.end()); vector<double> rui(prob.size() + 1); rui[0] = -1; rui[1] = prob[0].first; for (int i = 1; i < prob.size(); i++) { rui[i+1] = rui[i] + prob[i].first; } //大きい数字 int num = rand() % (int)(rui[rui.size() - 1]+1); //確率的にどれを選ぶかを決める int tmp_index = --lower_bound(rui.begin(), rui.end(), num) - rui.begin(); int index = prob[tmp_index].second; int candidate_Score = fs[index]; if(candidate_Score > calc_sum_profit(x,profit,n)){ for (int i = 0; i < n; i++) { x[i] = copy_cat[index][i]; } } } //tracing_mode void tracing_mode(vector<int> & x, vector<int> & x_best,int n) { //論文通りに実装 double d1 = 0,d0 = 0; double V1 = 0,V0 = 0,V_prime = 0; vector<double> t(n, 0); double c1 = 1, r1 = 0.5, w = 1; for (int i = 0; i < n; i++) { if(x_best[i]){ d1 = r1 * c1; d0 = -r1 * c1; } else{ d1 = -r1 * c1; d0 = r1 * c1; } V1 = w * V1 + d1; V0 = w * V0 + d0; if(x[i]){ V_prime = V0; } else{ V_prime = V1; } t[i] = (double) (1.0 / (1.0 +(double) exp(-V_prime))); if((double) rand()/RAND_MAX >= t[i]){ x[i] = x_best[i]; } } } //初期化関数 void init(vector<vector<int>> & cat,int C,int n ){ for (int i = 0; i < C; i++) { for (int j = 0; j < n; j++) { cat[i][j] = rand() % 2; } } } signed main(){ srand((unsigned int)time(NULL)); int n; int capacity; cin >> n; cin >> capacity; vector<int> weight(n); vector<int> profit(n); for (int i = 0; i < n; i++) { cin >> weight[i]; } for (int i = 0; i < n; i++) { cin >> profit[i]; } int C = 20; int Iterations = 50; vector<vector<int>> cat(C, vector<int>(n, 0)); vector<bool> seek_or_trace(C,0); vector<int> x_best(n); //初期化 init(cat, C, n); //暫定的な最大値の配列を取得 int max_profit = 0; for (int i = 0; i < C; i++) { if(check(cat[i],weight,capacity,n)){ int sum1 = calc_sum_profit(cat[i], profit, n); if (sum1 > max_profit) { max_profit = sum1; for (int j = 0; j < n; j++) { x_best[j] = cat[i][j]; } } } } //ここからスタート int counter = 0; time_t pretime = time(NULL); while (counter < Iterations) { ///max_profit = calc_max(C, n, weight, profit, x_best, cat, capacity, max_profit); for (int i = 0; i < C; i++) { vector<int>tmp=cat[i]; if((double) rand()/RAND_MAX < MR){ seeking_mode(cat[i],profit,weight,capacity,n); } else{ tracing_mode(cat[i], x_best, n); } //更新されたもとの配列が制約を満たしていなかったらもとに戻す if(!check(cat[i],weight,capacity,n)){ cat[i]=tmp; } } max_profit = calc_max(C, n, weight, profit, x_best, cat, capacity, max_profit); counter++; //ここのコメントアウトを外すと毎回の更新の際の配列の変化がわかります。 // show(x_best); // sleep(1); } show(x_best); cout << max_profit << endl; }
参考文献
Ricardo Soto , Broderick Crawford , Angelo Aste Toledo, Hanns de la Fuente-Mella , Carlos Castro , Fernando Paredes , and Rodrigo Olivares (2019) "Solving the Manufacturing Cell Design Problem through Binary Cat Swarm Optimization with Dynamic Mixture Ratios
