進化計算ライブラリDEAPで高次方程式を解く

はじめに

こんにちは。株式会社コズムの開発部でインターンをしています、鈴木です。

前回に引き続き、遺伝的アルゴリズム第２弾！

ということで今回は、強力な進化計算ライブラリであるDEAPを使って数学の難問にチャレンジします。

遺伝的アルゴリズムの基礎については、前回の記事で「０から実装」を通して説明していますので、ぜひご一読ください。

前回記事はこちら

https://zenn.dev/cosmcorp/articles/2e01cd279f5e63

本題

今回は、以下の五次方程式の解を遺伝的アルゴリズム（GA）で求めます。あえて適当に作ったため、私は左辺の関数の振る舞いについて何も知りません。

$$x^5 + 2x^4 + 3x^3 + 4x^2 + 5x + 6 = 0$$

一般に、5次以上の方程式は解の公式が存在せず、代数的に解くのが不可能とされています。こういった直接的に求めることの難しい解を探索的に求めるのには、GAはぴったり。

5次方程式は数学的にも興味深い話題なので、ぜひ調べてみてください。

使用外部ライブラリ

• DEAP（Distributed Evolutionary Algorithms in Python）

  進化的アルゴリズム、特に遺伝的アルゴリズム、遺伝的プログラミング、および関連する手法を簡単に実装できるように設計されたPythonライブラリ。

https://github.com/DEAP/deap

実装

それでは実装に移ります！

ライブラリの使い方については概要の説明に留めているので、詳細は公式ドキュメントを確認してみてください。

ステップ１：下準備

まずは必要なライブラリをインポートし、最適化に使う目的関数と評価関数を定義。

左辺の多項式が0になれば良いので、func(x)の値の絶対値を評価関数としてそれを最小化するようにします。

from deap import base, creator, tools, algorithms
import random
import sys

def func(x):  #　目的関数
    return x**5 + 2*x**4 + 3*x**3 + 4*x**2 + 5*x +6

def evaluate(individual):  # 評価関数
    x = individual[0]
    return abs(func(x)),

ステップ２：DEAPの設定

ステップ２では、DEAP内のクラスに関数を定義します。

toolboxに、使用する関数がすべて登録できたら完了。DEAPに用意されていないものを使用したい場合は、creator.createを用いてあらかじめ作成してからtoolboxに登録します。今回は個体生成の関数と評価関数のみ自前で作成し、交叉、突然変異、選択の関数はDEAPに用意されているものを使いました。

"""
creator.create()関数で新しいクラスを動的に生成する。
"""
creator.create("FitnessMin", base.Fitness, weights=(-1.0,)) # 適応度クラスを作成。最大化問題の場合、weightsを正の値に設定。
creator.create("Individual", list, fitness=creator.FitnessMin) # 個体クラスを作成。リストを基にし、FitnessMaxを適応度として使用。

"""
世代生成の関数をtoolboxに登録する
"""
toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, -sys.maxsize, sys.maxsize)  # 属性生成関数を登録（すべての実数からランダムに選択）
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, 1)　　 # 個体生成関数を登録（すでに定義したcreator.Individualとtoolbox.attr_floatを参照）
toolbox.register("population", tools.initRepeat, list, toolbox.individual)　 # 個体群（世代）生成関数を登録(直前で定義したtoolbox.individualを参照）

"""
評価、交差、突然変異、選択の関数をtoolboxに登録
"""
toolbox.register("evaluate", evaluate)　　　# ステップ１で作成したevaluate関数を登録
toolbox.register("mate", tools.cxBlend, alpha=0.5)　 #　あらかじめ用意されているcxblend（（ブレンド交叉）を登録
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.2)　　#　あらかじめ用意されているmutGaussian(ガウス突然変異)を登録
toolbox.register("select", tools.selTournament, tournsize=3) # あらかじめ用意されているselTournament(トーナメント選択）を登録

6行目で定義しているtoolbox.populationは、toolbox.individualを繰り返し実行させリストに格納します。さらにこのtoolbox.individualは「toolbox.attr_floatで決定した要素をcreator.Individualに渡す」という動作を繰り返すものです。

すこし構造が複雑に思えるかもしれませんが、小さい構成要素から順に作っていって大きい枠組みに組み込んでいく、という様子をイメージすると分かりやすいです。

＜補足＞

DEAPでは、toolsというクラスにデフォルトの関数が登録されています。
そこで、toolsに用意されている交叉関数についていくつかご紹介します。

今回は６番のtools.cxBlendを使用しています。（ステップ２の一番下のブロックの上から２行目）

1. 一点交叉 (One Point Crossover)

• 関数名: cxOnePoint

  一つの点をランダムに選び、その点を境にして二つの親の遺伝子を交換する。

2. 二点交叉 (Two Point Crossover)

• 関数名: cxTwoPoint

  二つの点をランダムに選び、その二点の間の遺伝子を交換する。

3. 一様交叉 (Uniform Crossover)

• 関数名: cxUniform

  各遺伝子座について一定の確率で親を交換する。

5. 順序交叉 (Order Crossover, OX)

• 関数名: cxOrdered

  順序問題に適した交叉方法で、親の一部の遺伝子を選び、その順序を保持しながら残りの遺伝子を他の親から受け継ぐ。

6. ブレンド交叉 (Blend Crossover, BLX-alpha)

• 関数名: cxBlend

  実数値の遺伝子を持つ個体に対して使用され、親の遺伝子値を基にして新たな遺伝子値を生成します。alpha パラメータでブレンドの度合いを調整する。（今回は標準値の0.5に設定）

そのほかにも、部分的交叉、サイクル交叉、シンプレックス交叉などがあります。toolsには、交叉意外のステップにおいてもGAに有効な関数が多数登録されているので、ぜひ確認してみてください。

ステップ３：実行

いよいよ実行です。高次方程式は複数の実数解を持つことがあるので、アルゴリズムを繰り返し実行し新しい解が見つかるたびに保存していきます。

DEAPでは、ステップ２でtoolboxに登録した関数をalgorithmsというクラスで簡単に実行することができます。まさに至れり尽くせりですね。

NGEN = 200  　　　　　　　#　世代数
TIMES = 10　　　　　　　　　　　#　アルゴリズムの実行回数
answers = set()　#　解を格納する集合

"""
TIMES回実行して、すべての実数解を探す
"""

for i in range (TIMES):
    print(f"times:{i}")  # 何回目の実行か表示 
    population = toolbox.population(n=1000)
    for gen in range(NGEN):
        offspring = algorithms.varAnd(population, toolbox, cxpb=0.3, mutpb=0.3)  # DEAPのアルゴリズム実行ツール（交叉確率と変異確率を設定）
        fits = toolbox.map(toolbox.evaluate, offspring)  # それぞれの個体に評価関数を適応し評価を出す
        for fit, ind in zip(fits, offspring):
            ind.fitness.values = fit

        fitnesses = [ind.fitness.values[0] for ind in offspring]
        best_fitness = min(fitnesses)
        mean_fitness = sum(fitnesses) / len(fitnesses)
        population = toolbox.select(offspring, k=len(population))
        if gen % 20 == 0:  # 実行中の状況を確認
            print(f"Generation {gen}: Best {best_fitness}, Mean {mean_fitness}")
            
    best_ind = tools.selBest(population, 1)[0]  # 最良個体の表示
    best_solution = best_ind[0]
    if func(best_solution)==0 and best_solution not in answers:  # 新しい解が見つかったら保存
        answers.add(best_solution)
        print("answers is added")
        
for answer in answers:  # 最後に見つかったすべての解を表示
    print(f"answer: {answer}")
   

結果

200世代×10回の実行結果が以下です。1番下に表示されているのが見つかった方程式の解。

times:0
Generation 0: Best 4.6440861097513555e+76, Mean 1.0754159363438404e+94
Generation 50: Best 1.3578197416208582e-05, Mean 2.6787667526115246
Generation 100: Best 0.0, Mean 2.6088655495421524
Generation 150: Best 0.0, Mean 1.9740695328248732
answers is added
times:1
Generation 0: Best 2.524635447533104e+80, Mean 9.433914469237196e+93
Generation 50: Best 3.146895106986847e-05, Mean 1.6126328972580957
Generation 100: Best 0.0, Mean 3.423718851412908
Generation 150: Best 0.0, Mean 2.4557431563103234
times:2
Generation 0: Best 2.8160492436575244e+80, Mean 1.0446997116737236e+94
Generation 50: Best 1.4597320994624852e-05, Mean 2.705770482005267
Generation 100: Best 0.0, Mean 1.7825787093178675
Generation 150: Best 0.0, Mean 3.2847364316431964
times:3
Generation 0: Best 1.2483796394022157e+75, Mean 1.019867444758366e+94
Generation 50: Best 5.816265231217699e-06, Mean 1.7750558801759515
Generation 100: Best 0.0, Mean 4.857211533537652
Generation 150: Best 0.0, Mean 0.8952336215028022
times:4
Generation 0: Best 8.216398627417412e+80, Mean 9.528776804292449e+93
Generation 50: Best 8.090321568232639e-08, Mean 2.1255774876329716
Generation 100: Best 0.0, Mean 1.440400321605238
Generation 150: Best 0.0, Mean 2.4343796949422027
times:5
Generation 0: Best 7.43206908386961e+79, Mean 1.0157365974981632e+94
Generation 50: Best 1.630879285841047e-06, Mean 1.9867055110379508
Generation 100: Best 0.0, Mean 0.9283067129741072
Generation 150: Best 0.0, Mean 5.206093297610345
times:6
Generation 0: Best 1.197063202726468e+80, Mean 1.0073418090102469e+94
Generation 50: Best 0.000649166893611941, Mean 2.3151620389981824
Generation 100: Best 0.0, Mean 2.4452749120272985
Generation 150: Best 0.0, Mean 1.4484867765675278
times:7
Generation 0: Best 1.1058331955362586e+73, Mean 1.0708417794951476e+94
Generation 50: Best 5.672377598742173e-07, Mean 3.303427614302371
Generation 100: Best 0.0, Mean 1.9295717596476036
Generation 150: Best 0.0, Mean 2.119754511564904
times:8
Generation 0: Best 1.6826087253367982e+78, Mean 9.12611836962651e+93
Generation 50: Best 0.00018675920866861873, Mean 1.8177097821485593
Generation 100: Best 0.0, Mean 1.4240113016013958
Generation 150: Best 0.0, Mean 1.6984712106743147
times:9
Generation 0: Best 7.076660584969771e+81, Mean 1.03657571858495e+94
Generation 50: Best 8.537185470558484e-06, Mean 2.4502106011232514
Generation 100: Best 0.0, Mean 1.4444652635092747
Generation 150: Best 0.0, Mean 3.275122879950779
answer: -1.4917979881399008

なんと、解が１つしか見つかりませんでした、、、

試行回数を増やして、再度実行してみる

TIMES(試行回数)を100回に増やして再度実行してみます。見やすくするために途中経過のプリントを省きました。

answers is added
answer: -1.4917979881399008

やはり、解が１つしか見つかりません。

ここで考えられるのは、局所解に陥ってしまい探索が十分にできていないか、そもそも適当に作ったこの方程式が実数解を１つしか持たないかの２パターンです。

文明の利器、使います

ずるいような気がしますが、気にせず使います。グラフ描写ツールのGeoGebraです。

https://www.geogebra.org/graphing?lang=ja
なんと私が適当に作った関数は、実数解を１つしか持たないものでした。お題ミス。

高次の関数はパラメータをほんの少し変えるだけで振る舞いが大きく変化するので、こういうこともあります。
x軸との交点（方程式の解）は先ほどの結果とピッタリですね。

コード全文

コメントを省いて全文載せておくので、お好きな関数を入力したりして遊んでみてください。

from deap import base, creator, tools, algorithms
import random
import sys

def func(x):
    return x**5 + 2*x**4 + 3*x**3 + 4*x**2 + 5*x +6

def evaluate(individual):
    x = individual[0]
    return abs(func(x)),

creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)

toolbox = base.Toolbox()
toolbox.register("attr_float", random.uniform, -sys.maxsize, sys.maxsize)
toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_float, 1)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("evaluate", evaluate)
toolbox.register("mate", tools.cxBlend, alpha=0.5)
toolbox.register("mutate", tools.mutGaussian, mu=0, sigma=1, indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)

NGEN = 200
TIMES = 100
answers = set()

for i in range (TIMES):
    print(f"times:{i}")
    population = toolbox.population(n=1000)
    for gen in range(NGEN):
        offspring = algorithms.varAnd(population, toolbox, cxpb=0.3, mutpb=0.3)
        fits = toolbox.map(toolbox.evaluate, offspring)
        for fit, ind in zip(fits, offspring):
            ind.fitness.values = fit

        fitnesses = [ind.fitness.values[0] for ind in offspring]
        best_fitness = min(fitnesses)
        mean_fitness = sum(fitnesses) / len(fitnesses)
        population = toolbox.select(offspring, k=len(population))
        if gen % 50 == 0:
            print(f"Generation {gen}: Best {best_fitness}, Mean {mean_fitness}")

    best_ind = tools.selBest(population, 1)[0]
    best_solution = best_ind[0]
    if func(best_solution)==0 and best_solution not in answers:
        answers.add(best_solution)
        print("answers is added")
        
for answer in answers:
    print(f"answer: {answer}")

まとめ

今回は遺伝的アルゴリズムのライブラリであるDEAPを用いて、高次方程式の解を求めました。

ライブラリを用いることで、前回記事のように０から実装したときと比べはるかに短いコードで実装できた上に、交叉や選択などのデフォルトで用意されている関数が非常に強力です。

ここまで完璧にサポートされている進化計算のライブラリを私は知らないので、スピード感を持って実装したい場合の選択肢はDEAP一択でしょう。

正直なところ、この問題に関しては関数解析のソフトを用いた方が簡単かつ高速に答えにありつけるのはご存知かと思いますが、今回の趣旨はライブラリの紹介だったので、その点はご容赦いただければと思います。

最後までお読みいただきありがとうございました！