東大数学を遺伝的アルゴリズムで解く

はじめに

こんにちは、コズムの開発部でインターンをしています、鈴木です。
私は機械学習やアルゴリズムを用いた開発をメインに行っており、その中でも特に思い入れの深いアルゴリズムについて執筆する運びとなりました。

ということで今回は、プログラムを駆使して東京大学の入試数学を倒します。
東大生エンジニアも多数在籍する弊社ですが、非東大生である私が日本一の入試問題を解くにはどうしたらよいでしょうか、、、。

それはやはり、アルゴリズム。とりわけ遺伝的アルゴリズム一択でしょう！

遺伝的アルゴリズムってなんだ

進化論的考え方

遺伝的アルゴリズム(genetic algorithm、以下GA)とは、生物の遺伝様式を模した、最適化アルゴリズムの一種です。自然界では、環境に適応した強者が生き残り、遺伝を繰り返していくうちに先の世代では強者の要素をもった個体が支配的となります。

最初、キリンの首の長さには個体差があったのに、高い木になった果実を食べられる首の長いキリンが有利に生き残った結果、後の世代では首の長いキリンだらけになった、という例は有名ですね。進化論の詳細については本記事の主題ではないので割愛します。

どうコードに落とし込むか

コーディングにおけるGAで、キリンの例で言うところの「高い木」があるような世界を用意するのは、我々エンジニアです。何をもって優秀とするのかは、どの世界で生きるのかによって変わってきますよね。生き残らせたい特性から逆算して、エンジニアは世界を設定します。

いよいよ本題

ここで本題に入っていきます。今回取り組む問題と解答は以下です。

問題

（引用元：https://www5a.biglobe.ne.jp/~t-konno/math/mathematics_tokyo.htm)

解答

（引用元：https://www5a.biglobe.ne.jp/~t-konno/math/mathematics_tokyo.htm)

東京大学から1973年の出題です。比較的簡単に解けますが、場合分けが少し面倒。

ここで、受験生の皆さんはまずGAの実装を検討しましょう！

実装

今回は、ライブラリの使用は必要最小限にとどめて、GAを０から実装していきます。

最初の準備

まず、計算に必要なライブラリのインポートをします。その上で多少の数学的考察を挟み、問題文で与えられたV(a)を定義しています。

import random
import sys

# g(x) = f(x) - ax　とする。
# g(x) = 1-ax(1<=x<=2のとき),　(1-a)x-1(2<=x<=3のとき)　であり、gのグラフはすべての区間において/
# 直線であるから、gは区間の端点で最大最小をとる。

def V(a):
    #最大値、最小値の候補は以下のg1,g2,g3
    g1 = 1-a   
    g2 = 1-2*a
    g3 = 2-3*a
    V = max(g1,g2,g3) - min(g1,g2,g3)
    return V

初期世代生成

GAの最初のステップである、初期世代の生成をここで行います。

解の可能性を限定しないために、多様性をもった初期世代を生成することが重要です。ここではランダム関数を用いてこれを実現しました。

先ほどの例で言うところの、「色々な首の長さをもったキリン集団」がここで生み出されます。

def makeIndividual():
		"""
		初期世代を構成する個体を１つ生み出す関数
		"""
    a = random.uniform(-sys.maxsize, sys.maxsize)  # aはすべての実数にわたり動く
    return {'a': a, 'V(a)': V(a)}

def makePopulation(populationSize):
		"""
		makeIndividualを繰り返し実行して世代を作る関数
		"""
    return [makeIndividual() for _ in range(populationSize)]

選択

２つめのステップである選択では、前の世代から優秀な個体を選別し次の世代を生み出す親に任命します。この選択の手法は複数ありますが、一番シンプルなエリート選択を採用しました。

今回の場合、最終行のラムダ式の部分がどの個体を残すかの判断基準になっており、これを一般に評価関数と呼びます。良い親を選ばなければ優秀な特性が消滅してしまうので、大変重要なステップです。

ここでは「V(a)の値が大きいほど優秀である」という世界を設定しました。

def selectParents(population):
		"""
		次の世代を生み出す親を決定する関数
		"""
    return sorted(population, key=lambda x: x['V(a)'])[:2]

交叉

交叉では、選択した個体を「親」として、次の世代を構成する「子」を作ります。

親の良い特性が子に引き継がれるようにしなければアルゴリズムが破綻してしまうので、評価関数と同様に重要なステップです。

ここでは、２つの親個体のもつaの値の平均を子に引き継いでいます。

def crossover(parent1, parent2):
		"""
		親個体２つから子を作る関数
		"""
    child_a = (parent1['a'] + parent2['a']) / 2  # 平均をとる
    return {'a': child_a, 'V(a)': V(child_a)}

突然変異

突然変異により親にはない遺伝情報を偶発的に与えることで、局所解に陥ることを防ぎます。
局所最適化を防ぐ方法が多様性を保つことであるのは、現実世界においても示唆的ですね。

def mutate(individual, mutationRate=0.01):
		"""
		突然変異を行う関数（変異確率はデフォルトで0.01に設定）
		"""
    if random.random() < mutationRate:
        individual['a'] += random.uniform(-1, 1)  
    individual['V(a)'] = V(individual['a'])

実行用の関数

最後に実行用の関数を作り、これを呼び出すだけでアルゴリズムが実行されるようにします。

１世代の個体数と、生み出す世代数を引数にとり、V(a)の最小値とそれを与えるaの値を返します。実行時に解の収束の様子を見えるように、50世代ごとにその時点での状況をプリントするようにしました。

def geneticAlgorithm(populationSize, generations):
		"""
		ここまで作ってきた関数を組み合わせて、実行用の大きな関数を作る
		"""
    population = makePopulation(populationSize)
    for generation in range(generations):
        parents = selectParents(population)
        child = crossover(parents[0], parents[1])
        mutate(child)
        population = sorted(population + [child], key=lambda x: x['V(a)'])[:populationSize]
        if generation%50 == 0:
            print('世代数{}:  min = {}, a = {}'.format(generation,population[0]['V(a)'],population[0]['a']))
    return population[0]

実行してみましょう。

実行

1000個の個体を持つ集団を、遺伝を繰り返し1000世代生成します。

best_individual = geneticAlgorithm(1000,1000)
print(f"min: {best_individual['V(a)']} at a = {best_individual['a']}")

結果

解けました！100世代目でほとんど正解に到達し、150世代目以降では完全に収束しています。

これで東大数学１完です。

世代数0:  min = 307661156085759.0, a = 153830578042880.0
世代数50:  min = 0.9678494095794417, a = 0.9678494095794417
世代数100:  min = 0.5000000000000018, a = 0.5000000000000018
世代数150:  min = 0.5, a = 0.5
世代数200:  min = 0.5, a = 0.5
世代数250:  min = 0.5, a = 0.5
世代数300:  min = 0.5, a = 0.5
世代数350:  min = 0.5, a = 0.5
世代数400:  min = 0.5, a = 0.5
世代数450:  min = 0.5, a = 0.5
世代数500:  min = 0.5, a = 0.5
世代数550:  min = 0.5, a = 0.5
世代数600:  min = 0.5, a = 0.5
世代数650:  min = 0.5, a = 0.5
世代数700:  min = 0.5, a = 0.5
世代数750:  min = 0.5, a = 0.5
世代数800:  min = 0.5, a = 0.5
世代数850:  min = 0.5, a = 0.5
世代数900:  min = 0.5, a = 0.5
世代数950:  min = 0.5, a = 0.5
min: 0.5 at a = 0.5

もう少し細かく見てみよう

優秀でない個体が淘汰されていく様子を見ていくために、世代数を絞ってもう一度実行してみます。

ここでは、世代数を200に変更した上、細かく追跡するため10世代毎に表示するよう書き換えました。GAでは、初期世代生成のランダム性によって収束に必要な世代数は実行ごとに変動するので、前回の収束結果の150から少し余裕を持たせて200に世代数を設定しています。

最初とてつもなく大きかったaの値が徐々に収束しています。
優秀でない個体が淘汰されていく様子が分かりますね。

世代数0:  min = 5.043892365767885e+16, a = 2.5219461828839424e+16
世代数10:  min = 27779222839967.0, a = 13889611419984.0
世代数20:  min = 17100688490.34375, a = -8550344244.671875
世代数30:  min = 94278199.79462111, a = -47139099.397310555
世代数40:  min = 60744.947265480354, a = -30371.973632740177
世代数50:  min = 73.60448820509683, a = 37.302244102548414
世代数60:  min = 0.5547715335606661, a = 0.44522846643933384
世代数70:  min = 0.5000381038324566, a = 0.4999618961675434
世代数80:  min = 0.5000000493185384, a = 0.4999999506814615
世代数90:  min = 0.5000000000424871, a = 0.5000000000424871
世代数100:  min = 0.5000000000001764, a = 0.49999999999982364
世代数110:  min = 0.5, a = 0.5
世代数120:  min = 0.5, a = 0.5
世代数130:  min = 0.5, a = 0.5
世代数140:  min = 0.5, a = 0.5
世代数150:  min = 0.5, a = 0.5
世代数160:  min = 0.5, a = 0.5
世代数170:  min = 0.5, a = 0.5
世代数180:  min = 0.5, a = 0.5
世代数190:  min = 0.5, a = 0.5
min: 0.5 at a = 0.5

コード全文

ここまでご説明してきたコードの全文を添付しておきます。

煮るなり焼くなりご自由にお役立てください。

import random
import sys

# g(x) = f(x) - ax　とする。
# g(x) = 1-ax(1<=x<=2のとき),　(1-a)x-1(2<=x<=3のとき)　であり、gのグラフはすべての区間において直線であるから、gは区間の端点で最大最小をとる。

def V(a):
    #最大値、最小値の候補は以下のg1,g2,g3
    g1 = 1-a   
    g2 = 1-2*a
    g3 = 2-3*a
    V = max(g1,g2,g3) - min(g1,g2,g3)
    return V

def makeIndividual():
    a = random.uniform(-sys.maxsize, sys.maxsize)  
    return {'a': a, 'V(a)': V(a)}

def makePopulation(populationSize):
    return [makeIndividual() for _ in range(populationSize)]

def selectParents(population):
    return sorted(population, key=lambda x: x['V(a)'])[:2]

def crossover(parent1, parent2):
    child_a = (parent1['a'] + parent2['a']) / 2 
    return {'a': child_a, 'V(a)': V(child_a)}

def mutate(individual, mutationRate=0.01):
    if random.random() < mutationRate:
        individual['a'] += random.uniform(-1, 1)  
    individual['V(a)'] = V(individual['a'])

def geneticAlgorithm(populationSize, generations):
    population = makePopulation(populationSize)
    for generation in range(generations):
        parents = selectParents(population)
        child = crossover(parents[0], parents[1])
        mutate(child)
        population = sorted(population + [child], key=lambda x: x['V(a)'])[:populationSize]
        if generation % 10 == 0:
            print('世代数{}:  min = {}, a = {}'.format(generation, population[0]['V(a)'], population[0]['a']))
    return population[0]

best_individual = geneticAlgorithm(1000,200)
print(f"min: {best_individual['V(a)']} at a = {best_individual['a']}")

まとめ

今回は、東京大学の入試数学を題材として遺伝的アルゴリズムを０から実装しました。

シンプルな1次関数だったのでスムーズに解にたどりつけましたが、極値が複数あるものや高次の関数の場合にはもう少し工夫が必要です。

多数ある最適化アルゴリズムの中で、面倒ごとを抜きにしてスピーディーに構築できるのが遺伝的アルゴリズムの良さですね。次回はもっと難度を上げた問題でさらにその威力をご紹介できたらと思います！

読んでいただきありがとうございました。