進化的アルゴリズム（EA）の概略とpythonでのEA計算用パッケージDEAPのチュートリアル
この記事について
進化的アルゴリズムを実装するためのフレームワークである、DEAPについて説明します。
この記事では、進化的アルゴリズムの概要について触れつつ、DEAPの公式のドキュメンテーション(DEAP1.3.0 documentation)のFirst stepsとBasic tutorialsを、要約しつつ日本語訳します。
進化的アルゴリズムについて
進化的アルゴリズム(evolutionary algorithm、EA)は、進化的計算の一分野で、個体群ベースのメタヒューリスティックな最適化アルゴリズムの総称です。解の探索において、解空間構造が不明で、その問題に固有な知識が使えず、全探索が不可能であるような広大な解空間を持つようなタイプの問題に対して、主に威力を発揮します。
EAは主に次の4つの手法からなります。
- 遺伝的アルゴリズム
- 遺伝的プログラミング
- 進化戦略
- 進化的プログラミング
この章ではそれぞれのアルゴリズムについて、簡単に説明します。
遺伝的アルゴリズム(GA)について
遺伝的アルゴリズムとは、生物進化における遺伝と適者生存の仕組みをソフトウェア的に真似することで、複雑な問題に対する最適解（の近似値）を模索する（メタヒューリスティクスな）方法です。
概要
遺伝的アルゴリズムはデータ（解の候補）を遺伝子で表現した「個体」を複数用意し、適応度の高い個体を優先的に選択して交叉・突然変異などの操作を繰り返しながら解を探索する手法です。適応度は適応度関数によって定義します。
この手法の利点は、評価関数の可微分性や単峰性などの知識がない場合であっても適用可能なことです。 必要とされる条件は評価関数の全順序性と、探索空間が位相（トポロジー）を持っていることです。
また、遺伝子の表現の仕方によっては組合せ最適化問題やNP困難な問題などのさまざまな問題に適用可能です。
具体的には、ある命題の解の候補を遺伝子(gene)とその集合体である、染色体(chromosome)で表現した個体(individual)を複数用意し、適応度(fitness)の高い個体を優先して選択し、交叉(crossover)、突然変異(mutation)などの操作を繰り返しながら、最適解の探索を行います。
流れ
遺伝的アルゴリズムは一般に以下の流れで実装されます。なお、下記では個体数を N, 最大世代数を G と置きます。
- あらかじめ N 個の個体が入る集合を二つ用意します。以下、この二つの集合を「現世代」、「次世代」と呼ぶことにします。
- 現世代に N 個の個体をランダムに生成します。
- 評価関数により、現世代の各個体の適応度をそれぞれ計算します。
- ある確率で次の3つの動作のどれかを行い、その結果を次世代に保存します。
- 個体を二つ選択（選択方法は後述）して交叉（後述）を行います。
- 個体を一つ選択して突然変異（後述）を行います。
- 個体を一つ選択してそのままコピーします。
- 次世代の個体数が N 個になるまで上記の動作を繰り返します。
- 次世代の個体数が N 個になったら次世代の内容を全て現世代に移します。
- 3.以降の動作を最大世代数 G 回まで繰り返し、最終的に「現世代」の中で最も適応度の高い個体を「解」として出力します。
遺伝的アルゴリズムの問題点
遺伝的アルゴリズムには、下に挙げるようないくつか注意しなければならない問題点があります。
初期収束
初期収束とは、最初の方の世代で偶然他の個体より、圧倒的に適応度の高い個体が生まれてしまったとき、集団の中でその個体の遺伝子が爆発的に増えてしまい、探索がかなり早い段階で収束してしまう減少です。ルーレット選択の設定が甘いときや、突然変異も効果がうまく表れないときに起こりやすいです。
対策としては、ルーレット選択を行う際に適切な設定をする、適応する問題に対して効果的になるように突然変異の操作を変更する、突然変異率を増やす、または、集団の数を増やすなどがあります。しかし、明確にこれで大丈夫という方法はなく、結局各パラメーターを何度も設定し直すしかありません。ヒッチハイキング
詳細は省きますが、交叉の手法として、二点交叉を用いた場合、最適解の近くで解がなかなか収束しないという問題です。一様交叉を使うとこの問題は発生しません。
進化戦略について
概要
進化戦略では、実数のベクトルで解を表現し、探索を行うと同時に自己変異用のパラメーターも更新していきます。
First steps
Overview（概観）
まずは、DEAP上で遺伝的プログラミングを実装する際の概要を説明します。これでDEAPのすばらしさが解ると思います。
Types
最初にやるべきことは、適切な問題のタイプを考えることです。これは、Creatingメソッドで行います。
Initilization
Typeを作ったらそこにランダムな初期値を入れます。これには、toolsというメソッドを使います。
Operators
Operatorの定義もInitializationと似ています。いくつかの演算はtoolsの中に予め入っていて、また、ユーザーが自分で演算を定義することも可能です。
Algorithms
上の3つで前段階は終わりです。いよいよアルゴリズムの実装に入ります。これは通常main functionの中で行われます。
Installation
Requirements
DEAPを使用するためには、Python3.4かPython2.7 以上のバージョンである必要があります。この記事では、Python3系の利用を想定しています。また、SCOOPとNumpyが導入されている必要があります。
Install DEAP
DEAPは次のpipコマンドでインストールできます。
pip install deap
Porting Guide
後方互換性についての話が乗っています。特に必要がない場合は読み飛ばしていいでしょう。
Creating Type
このチュートリアルでは、どのようにして、creatorメソッドでtypesが作られるのかとtooboxで初期化がされるのかを見せます。
Fitness（適応度）
Fitnessはweightアトリビュートを必要とする抽象クラス（抽象メソッドを持つクラス）です。適応度を最小化するためには、weightのパラメータを負にします。逆に適応度を最大化する場合は、正にします。下のコマンドはFitnessMinという名前の最小化適応度を作っています。
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
create()関数は少なくとも２つの引数を必要とします。新しく作るクラスの名前と、ベースとなるクラスです。それ以外の引数は全て、クラスのアトリビュートとなります。
Indivisual（個体）
個体を定義するクラスについての説明です。個体とは、遺伝子の集合として表されます。このクラスはNumpyから、引き継ぐことが可能なので、Numpyから引き継ぎたい人は、下の”Numpyからの引き継ぎチュートリアル”を参照してください。
https://deap.readthedocs.io/en/master/examples/ga_onemax_numpy.html
ここでは、いくつかの個体のタイプについて説明します。
List of Floats（実数値エンコーディング）
最初は遺伝子を、小数を含む単純なリストとして表現する実数値エンコーディングです。このような個体のリストを定義するためには、creatorメソッドを使って、まずIndivisualクラスを作る必要があります。
import random from deap import base from deap import creator from deap import tools creator.create("FitnessMax", base.Fitness, weights=(1.0,)) creator.create("Individual", list, fitness=creator.FitnessMax) IND_SIZE=10 toolbox = base.Toolbox() toolbox.register("attr_float", random.random) toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=IND_SIZE)
ここで新しく出てきたregister()メソッドの説明をします。registerメソッドは、少なくとも２つの引数を必要とします。エイリアスとそのエイリアスに渡す関数です。それ以外の引数はすべて、functools.partial()のコマンドが呼びだされたときに、関数に渡されます。
permutation（順列エンコーディング）
各遺伝子を順序を示す数字として表現する順列エンコーディングについての説明です。やり方は上でやったランダムな小数列によるエンコーディングとほとんど同じです。違いはtoolbox.registerの部分だけです。
import random from deap import base from deap import creator from deap import tools creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) creator.create("Individual", list, fitness=creator.FitnessMin) IND_SIZE=10 toolbox = base.Toolbox() toolbox.register("indices", random.sample, range(IND_SIZE), IND_SIZE) toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.indices)
Arismetic Expression(木構造エンコーディング)
最後に算術表現によく使われる木構造エンコーディングについて説明します。この場合は、PrimitiveSetという算術を定義の集合を作っておく必要があります。下のサンプルコードでは、PrimitiveSetは、"MAIN"という名前で呼び出されています。演算の取る引数の個数をarityといいます。また、個体は、木構造を生成するgenHalfAndHalf()関数によって生成されています。
import operator from deap import base from deap import creator from deap import gp from deap import tools pset = gp.PrimitiveSet("MAIN", arity=1) pset.addPrimitive(operator.add, 2) pset.addPrimitive(operator.sub, 2) pset.addPrimitive(operator.mul, 2) creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) creator.create("Individual", gp.PrimitiveTree, fitness=creator.FitnessMin, pset=pset) toolbox = base.Toolbox() toolbox.register("expr", gp.genHalfAndHalf, pset=pset, min_=1, max_=2) toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.expr)
Evolution Strategy（進化戦略）
進化戦略のための個体を作成する手順は、２つのリストを含むという点で、今までのものとは少しばかり異なります。リストのうちの一つは、実際の個体を表し、もう一つは、変異用のパラメータを表しています。今回のサンプルコードでは、listをベースクラスに使う代わりに、array.arrayを継承します。1つのオブジェクトの中に２つの異なるベクトルを生成する関数は用意されていないので、自前でそれを定義してやる必要があります。initES()関数が２つのクラスを受け取り、与えらた範囲内でランダムなパラメータを生成することで、それらのインスタンスを作成します。
import array import random from deap import base from deap import creator from deap import tools creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) creator.create("Individual", array.array, typecode="d", fitness=creator.FitnessMin, strategy=None) creator.create("Strategy", array.array, typecode="d") def initES(icls, scls, size, imin, imax, smin, smax): ind = icls(random.uniform(imin, imax) for _ in range(size)) ind.strategy = scls(random.uniform(smin, smax) for _ in range(size)) return ind IND_SIZE = 10 MIN_VALUE, MAX_VALUE = -5., 5. MIN_STRAT, MAX_STRAT = -1., 1. toolbox = base.Toolbox() toolbox.register("individual", initES, creator.Individual, creator.Strategy, IND_SIZE, MIN_VALUE, MAX_VALUE, MIN_STRAT, MAX_STRAT)
particle
粒子もまた、速度といままで通った経路の中で一番良い場所についての情報を持っているという点で特別なタイプの個体です。このタイプの個体は、今までやってきたようにcreateメソッドによって個体を定義することに加えて、speedとbestそして速度の上限と下限である(smin, smax)を定めてやる必要があります。ここでも、関数を自分で定義してあげる必要があります。
import random from deap import base from deap import creator from deap import tools creator.create("FitnessMax", base.Fitness, weights=(1.0, 1.0)) creator.create("Particle", list, fitness=creator.FitnessMax, speed=None, smin=None, smax=None, best=None) def initParticle(pcls, size, pmin, pmax, smin, smax): part = pcls(random.uniform(pmin, pmax) for _ in xrange(size)) part.speed = [random.uniform(smin, smax) for _ in xrange(size)] part.smin = smin part.smax = smax return part toolbox = base.Toolbox() toolbox.register("particle", initParticle, creator.Particle, size=2, pmin=-6, pmax=6, smin=-3, smax=3)
A Funky One
もし、取り扱っている問題が上で扱った以外の特別なパラメータを必要とする場合でもDEAPでは、簡単に実装することができます。 下のサンプルコードでは、ininCycle() functionによって、整数と実数のlistを作っています。
import random
from deap import base
from deap import creator
from deap import tools
creator.create("FitnessMax", base.Fitness, weights=(1.0, 1.0))
creator.create("Individual", list, fitness=creator.FitnessMax)
toolbox = base.Toolbox()
INT_MIN, INT_MAX = 5, 10
FLT_MIN, FLT_MAX = -0.2, 0.8
N_CYCLES = 4
toolbox.register("attr_int", random.randint, INT_MIN, INT_MAX)
toolbox.register("attr_flt", random.uniform, FLT_MIN, FLT_MAX)
toolbox.register("individual", tools.initCycle, creator.Individual,
(toolbox.attr_int, toolbox.attr_flt), n=N_CYCLES)
population（集団）
集団のクラスの作成方法も個体のクラスの作成方法と似ています。個体のクラスではアトリビュートで初期化していたところを、個体や戦略、粒子等で埋めるだけです。ここから、いくつかのタイプの集団についてのクラスの作成方法を説明します。
Bag
bag population(日本語名不明)は、最もよく使われる集団のタイプです。一般的にlistを用いて実装しますが、この集団の集合は、順序を持ちません。bagは特別なアトリビュートを持たないので、特別なクラスを必要としません。この集団は、toolboxのinitRepeat() functionを用いて、直接的に作成できます。
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
Grid
grid populationは近接する個体同士が相互作用し合うような場合に用いる特別な集団のクラスです。各個体は、格子の中に割り当てられます。それぞれの格子には1個体までしか入ることができません。この集団の実装は、複数のlistからなるという点以外はbagの場合とそれほど変わりません。
toolbox.register("row", tools.initRepeat, list, toolbox.individual, n=N_COL) toolbox.register("population", tools.initRepeat, list, toolbox.row, n=N_ROW)
Swarm(群)
swarm populationは粒子群最適化をするときに使われる手法です。クラスの内部に通信用のネットワークを持っているという点で上の２つとは異なります。使うことがしばらくなさそうなので、今回はこれの詳細についての説明は省かさせていただきます。
Demes
demeは副次的な集団で、集団の中に含まれます。使うことがしばらくなさそうなので、今回はこれの詳細についての説明は省かさせていただきます。
populationの実装例
import json from deap import base from deap import creator creator.create("FitnessMax", base.Fitness, weights=(1.0, 1.0)) creator.create("Individual", list, fitness=creator.FitnessMax) def initIndividual(icls, content): return icls(content) def initPopulation(pcls, ind_init, filename): with open(filename, "r") as pop_file: contents = json.load(pop_file) return pcls(ind_init(c) for c in contents) toolbox = base.Toolbox() toolbox.register("individual_guess", initIndividual, creator.Individual) toolbox.register("population_guess", initPopulation, list, toolbox.individual_guess, "my_guess.json") population = toolbox.population_guess()
operators and Algorithms
複雑なアルゴリズムを始める前に、DEAPについている基本的なアルゴリズムについて見てみましょう。最初は、（前節を参考にして）単純な個体クラスを作成することから始めます。そしていくつかの異なる演算子によって、その個体同士を相互作用させます。その後、アルゴリズムやその他のツールをどのようにして利用するか学びます。
A First Indivisual
最初に必要なモジュールをインポートして、個体クラスを作るのに要求される関数を登録します。
import random from deap import base from deap import creator from deap import tools IND_SIZE = 5 creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0)) creator.create("Individual", list, fitness=creator.FitnessMin) toolbox = base.Toolbox() toolbox.register("attr_float", random.random) toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, n=IND_SIZE)
最初の個体は次のようにして作成します。
ind1 = toolbox.individual() print(ind1) # [0.86..., 0.27..., 0.70..., 0.03..., 0.87...] print(ind1.fitness.valid) # False
Evaluation
評価は、進化的アルゴリズムの中で最も個別的な部分です。この部分だけは、ライブラリを使わないで、自分で定義を書かなければなりません。典型的な評価方法は、個体一つを引数としてとり、その個体の適応度をタプルにして返り値にしています。例えば、下のサンプルコードは、さっき作ったind1に対して適合度を計算しています。
def evaluate(individual):
# Do some hard computing on the individual
a = sum(individual)
b = len(individual)
return a, 1. / b
ind1.fitness.values = evaluate(ind1)
print ind1.fitness.valid # True
print ind1.fitness # (2.73, 0.2)
注意点として、評価関数の返り値は必ずタプル形式である必要があります。
Mutation（突然変異）
次に紹介する演算子は突然変異演算子です。deap.toolsモジュールの中にたくさんの突然変異演算子が入っています。それぞれの突然変異には、固有の特徴があり、個体のタイプに合わせて選択する必要があります。望まない挙動を避けるためにも、自分が使用するオペレーターのドキュメンテーションはしっかり読みましょう。
突然変異に関する一般的なルールとしては、それらは変異させることだけするということです。つまり、突然変異演算子を作用させる前に、もし変異する前の個体を保存しておきたいなら、その個体のコピーを作っておいてあげる必要があるということです。
mutant = toolbox.clone(ind1) ind2, = tools.mutGaussian(mutant, mu=0.0, sigma=0.2, indpb=0.2) del mutant.fitness.values
一番下の行で適合度の値を削除したのは、そこには変異させる前の個体の適合度が格納されているからです。また、上で突然変異演算子を作用させたことによって、mutantとind2が同じであることに注意してください。
print (ind2 is mutant) # True print (mutant is ind1) # False
Crossover（交叉）
次の紹介する演算子は、交叉演算子です。突然変異演算子同様deap.toolsというところに交叉演算子がたくさん入っていて、個体のタイプに合わせて、選んでやる必要があります。
交叉演算子に関する一般的なルールとしては、この演算子は、個体同士を配偶させるだけであるということです。つまり、交叉演算子を作用させる前に、もし交叉させる前の個体を保存しておきたいなら、その個体のコピーを作っておいてあげる必要があるということです。 では実際に予めクローンをつくておいた個体に対して、交叉演算子を作用させてみます。
child1, child2 = [toolbox.clone(ind) for ind in (ind1, ind2)] tools.cxBlend(child1, child2, 0.5) del child1.fitness.values del child2.fitness.values
Selection（選択）
選択は集団クラスに対して作用させる演算子で、deap.toolsに格納されています。選択演算子の第一引数には、個体のiterableな変数を入れる必要があり、第二引数には選択する個体の数を入れます。そして返り値は、選択された個体のリストとなっています。
selected = tools.selBest([child1, child2], 2) print child1 in selected # True
そして、たいていは選択のあとか、個体を変化させる前に、すべての集団の複製を作ります。
selected = toolbox.select(population, LAMBDA) offspring = [toolbox.clone(ind) for ind in selected]
Using Toolbox
toolboxの使い方について書いてあります。要約すると、deap.toolsの中にある諸々の関数はToolboxに登録して使うと管理しやすいというようなことが書いてあります。Toolboxには、register()とunregister()の２つのメソッドを主に使って、toolの登録と解除を行います。
from deap import base from deap import tools toolbox = base.Toolbox() def evaluateInd(individual): # Do some computation return result, toolbox.register("mate", tools.cxTwoPoint) toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.2) toolbox.register("select", tools.selTournament, tournsize=3) toolbox.register("evaluate", evaluateInd)
Computing Statistics
しばしば、最適化の過程で何が行われているのかを調べるために、統計量を計算したくなることがある。Statisticsは任意のデザインされたオブジェクトのアトリビュートに対して、統計量を取得することができます。そのためには、望む統計オブジェクトが内包されている統計関数をtoolboxに登録する必要があります。
stats = tools.Statistics(key=lambda ind: ind.fitness.values)
statisticsオブジェクトは最初の引数にkeyを指定する必要があります。
stats.register("avg", numpy.mean) stats.register("std", numpy.std) stats.register("min", numpy.min) stats.register("max", numpy.max)
うえのコードによって、統計関数が登録されました。register関数は、第一引数にエイリアスを要求して、ベクトルに作用する関数を第二引数に要求します。
Predefined Algorithms
eaSimple(), eaMuPlusLambda(), eaMuCommaLambda(), eaGenerateUpdate()というような予め用意されているアルゴリズムを使う場合、統計量オブジェクトは、アルゴリズムに引数として渡せば良いです。
pop, logbook = algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=0, stats=stats, verbose=True)
こうすると、集団のすべての世代ごとに統計量が毎回計算されます。verboseオプションをTrueにしてると、最適化が行われている間にも、統計量が表示されるようになります。アルゴリズムが終了すると、終状態における集団とLogbackが出力されます。
Writing Your Own Algorythm
自分でアルゴリズムを書くときは、統計量を含めることは、とても簡単です。compile()メソッドを呼び出すだけで良いです。
record = stats.compile(pop)
そしてこの統計量の値は、print()関数で呼び出すことができます。
>>> print(record)
{'std': 4.96, 'max': 63.0, 'avg': 50.2, 'min': 39.0}
Multi-objective Statistics
numpyの関数を使って、統計量は直接計算されているので、すべてのオブジェクトはnumpyのデフォルトの振る舞いでまとめることができます。それゆえ、どれについて統計関数を作用させるかは、statsオブジェクトを宣言するときに、指定する必要があります。
stats = tools.Statistics(key=lambda ind: ind.fitness.values) stats.register("avg", numpy.mean, axis=0) stats.register("std", numpy.std, axis=0) stats.register("min", numpy.min, axis=0) stats.register("max", numpy.max, axis=0)
Multiple Statistics
Logging Data
統計関数で計算された統計量は、Logbookに保存され、あとから呼び出せるようになっています。
logbook = tools.Logbook() logbook.record(gen=0, evals=30, **record)
Logbookから特定の統計量だけを抜き出してくることもできます。
gen, avg = logbook.select("gen", "avg")
Logbookはピッケル化可能なオブジェクトです。
import pickle
pickle.dump(logbook, lb_file)
Printing to Screen
logbookはスクリーンやファイルに出力することができます。
logbook.header = "gen", "avg", "spam"
>>> print(logbook) gen avg spam 0 [ 50.2]
Dealing with Multi-statistics
logbookはMultiStatisticsオブジェクトの返り値である辞書の辞書型の値にも対応しています。それぞれの辞書型変数は、chaptersの中に格納されています。なので、複数のオブジェクトの統計も一つの場合と同じように扱うことができます。
logbook = tools.Logbook() logbook.record(gen=0, evals=30, **record)
ひとつの違いは、カラムの順序で、したのようにします。
logbook.header = "gen", "evals", "fitness", "size" logbook.chapters["fitness"].header = "min", "avg", "max" logbook.chapters["size"].header = "min", "avg", "max"
出力は以下のようになります。
>>> print(logbook)
fitness size
------------------------- ---------------
gen evals min avg max min avg max
0 30 0.165572 1.71136 6.85956 3 4.54 7
データを引き出すには、下のようにchaptersを介して行います。
gen = logbook.select("gen") fit_mins = logbook.chapters["fitness"].select("min") size_avgs = logbook.chapters["size"].select("avg")
Some Plotting Suger
最適化が終わったときに最も一般的によくする操作の一つは、進化の過程のデータをプロットすることです。Logbookによって、この作業はとても効率的に行なえます。これは、matplotlibを使って行います。
gen = logbook.select("gen") fit_mins = logbook.chapters["fitness"].select("min") size_avgs = logbook.chapters["size"].select("avg") import matplotlib.pyplot as plt fig, ax1 = plt.subplots() line1 = ax1.plot(gen, fit_mins, "b-", label="Minimum Fitness") ax1.set_xlabel("Generation") ax1.set_ylabel("Fitness", color="b") for tl in ax1.get_yticklabels(): tl.set_color("b") ax2 = ax1.twinx() line2 = ax2.plot(gen, size_avgs, "r-", label="Average Size") ax2.set_ylabel("Size", color="r") for tl in ax2.get_yticklabels(): tl.set_color("r") lns = line1 + line2 labs = [l.get_label() for l in lns] ax1.legend(lns, labs, loc="center right") plt.show()
Using Multiple Processors
この節には、分散処理のやり方についての説明が書いてありますが、今回は省略します。
