はじめに

Python、DEAPを使って遺伝的アルゴリズムを用いた特徴選択について調査し、コマンドを実装したのでメモしておく。

遺伝的アルゴリズムとは

他にも用途はあるが、今回の特徴選択では、組み合わせ最適化手法の一つとして使うため、その前提で説明する。
遺伝的アルゴリズムのとっつきにくいところは、「個体」「遺伝子」「選択」「交差」「突然変異」「適応度」といった生物学的な用語で説明されるため、その概念が理解しにくいということである。

組み合わせ最適化の文脈で、自分の理解でアルゴリズムをざっくりと説明すると以下のとおりである。

組み合わせ解をいくつか適当に生成。（組み合わせ解の1つ1つを「個体」いう）
組み合わせ解1つ1つを何らかの基準(「適応度」という）で評価する
評価に基づき「選択」「交差」「突然変異」等の操作を行い新たな組み合わせ解を何個か生成。
1～3を何回か(これを世代という)実行し、適応度がこれ以上改善しないぐらいのところでやめて、よさそうな組み合わせ解を採用する

特徴選択における遺伝的アルゴリズム適用の目的

特徴選択とは、大量の説明変数を、学習に寄与しそうな説明変数に絞るということである。特徴選択をすることで、ノイズが除去される、オーバーフィットの防止になる、学習時間が速くなる等のメリットがある。特徴選択をする場合、一つ一つの特徴と目的変数との関係によって選択する方法があるが、自分の経験上、精度についてはあまり効果がない場合もある。これは説明変数同士に関係があり、その組み合わも考慮しなければならないからだと考えられる。
（現実にはそのケースの方が圧倒的に多いと思われる）

そこで遺伝的アルゴリズムの登場である。

特徴同士の最適な組み合わせを見つけるということを、遺伝的アルゴリズムを使って行うわけである。
この場合、「適応度」は何に基づいて行えばいいだろうか。特徴を減らして精度をよくしたいわけだから、

特徴数の少なさ
その特徴をつかって予測モデルを作成した場合の精度の高さ

で評価すればよいわけである。

前置きが長くなったが、これをDEAPというフレームワークを使ってやってみようと思う。

実施環境

実施環境は以下の通り。
- Windows10
- Anaconda
- Python3.7
- DEAP1.3
- scikit-learn 0.21

DEAPとは

Pythonによる進化的計算フレームワークである。ライセンスは GNU Lesser General Public License v3.0 である。
遺伝的アルゴリズムの他にも、以下URLのように様々なことができる。（PSO以外は詳しくは知らないが。。。)
https://deap.readthedocs.io/en/master/examples/index.html

DEAPのインストール

Anaconda 環境では以下一発でインストールできる。

conda install deap

ソースコード

遺伝的アルゴリズムによる特徴選択を行うコマンドを作成。ほぼ、こちらの流用になっている。
変更したところは後で記載する。

GAFeatureSelection.py

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import RidgeCV
import random
from deap import creator, base, tools, algorithms
import argparse


# 評価関数
def get_fitness(individual, X, y):

    n_features = sum(individual)

    if n_features == 0:
        return 9999, -9999

    x = X.iloc[:, [bool(val) for val in individual]].values

    scores = []

    for _ in range(30):
        X_train, X_test, y_train, y_test = train_test_split(x, y, test_size=0.4)

        # normalizeされていないデータをinputと想定しているため、normarized=Trueとした。
        model = RidgeCV(scoring="r2", normalize=True)
        model.fit(X_train, y_train)
        scores.append(model.score(X_test, y_test))

    score = np.array(scores).mean()

    print("n_feature={0}, score={1}".format(n_features, score))

    return n_features, score


def main():

    parser = argparse.ArgumentParser()
    parser.add_argument("-input", type=str, required=True)
    parser.add_argument("-output", type=str, required=True)

    args = parser.parse_args()
    print("Input File : " + args.input)
    print("Output File : " + args.output)

    # データの読み込み
    df = pd.read_csv(args.input, index_col=0)

    # 目的変数の取得
    y = df.iloc[:, 0].values

    # 説明変数の取得
    df_data = df.iloc[:, 1:]

    # 遺伝的アルゴリズム開始
    print("Starting Genetic Algorizm By Random Forest With {0} column".format(len(df_data.columns)))

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

    toolbox = base.Toolbox()

    toolbox.register("attr_bool", random.randint, 0, 1)
    toolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_bool, len(df_data.columns))
    toolbox.register("population", tools.initRepeat, list, toolbox.individual)

    toolbox.register("evaluate", get_fitness, X=df_data, y=y)

    toolbox.register("mate", tools.cxOnePoint)

    # 下は正解率が低く特徴数が減らないため参考サイトから変更
    # toolbox.register("mate", tools.cxTwoPoint)

    toolbox.register("mutate", tools.mutFlipBit, indpb=0.05)
    toolbox.register("select", tools.selTournament, tournsize=3)

    # 以下は正解率が低く特徴数が減らないため参考サイトから変更
    # toolbox.register("select", tools.selNSGA2, k=10)

    # 世代数の設定
    pop = toolbox.population(n=30)
    hof = tools.ParetoFront()
    stats = tools.Statistics(lambda ind: ind.fitness.values)
    stats.register("avg", np.mean, axis=0)
    stats.register("std", np.std, axis=0)
    stats.register("min", np.min, axis=0)
    stats.register("max", np.max, axis=0)

    # 実行
    pop, log = algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2,
                                   ngen=50, stats=stats, halloffame=hof,
                                   verbose=True)

    # 殿堂入り(Hall Of Fame)の解を全てファイルに保存
    for individual in hof:
        remained_column = []
        print("{0},{1}".format(individual.fitness.values[0], individual.fitness.values[1]))

        for column, is_support in zip(df_data.columns, individual):
            if is_support == 1:
                remained_column.append(column)

        # 複数の解があるため、全てをファイルに保存する
        print("Writing Output File With {0} remaind columns".format(len(remained_column)))
        drop_columns = list(set(df_data.columns) - set(remained_column))
        df.drop(drop_columns, axis=1).to_csv(args.output +"_{0}_{1}.csv".format(individual.fitness.values[0], individual.fitness.values[1]))


if __name__ == "__main__":
    main()

ソース解説

ほぼ元ネタの流用なので、要点だけ説明する。

評価関数の登録

get_fitnessの中で評価を行っている。問題はこの関数をどうDEAPに認識させ、その関数のINとOUTはどうかけばよいのかとうことである。上のソースでDEAPに関数を登録するところは、toolbox.register("evaluate", get_fitness, X=x, y=y)の箇所である。ここでは、get_fitnessという関数であり、その引数がXとyをとるというものであるということを示している。(Xとyは、本記事では機械学習の予測が目的であるため、学習データと目的変数に対応して設定したもの）

評価関数の実装

def get_fitness(individual, X, y):がそれになる。
第一引数には、DEAPが生成した個体の情報（組み合わせ解）、第二、第三引数には登録時にしたXとyを受け取るものととして定義する。

individualの形として、今回の場合は、特徴の数、各特徴についてその特徴を使うか(=1)、使わないか(0)、つまりビットの並びになるわけで、それをtoolbox.register("individual", tools.initRepeat, creator.Individual, toolbox.attr_bool, len(df_data.columns))で指定している。

評価関数の中は、Xの中からindividualで指定された特徴および、yをつかってクロスバリデーションを行い、戻り値として特徴数と精度を返却している。

評価についてはcreator.create("FitnessMax", base.Fitness, weights=(-1.0, 1.0))で指定している。ソース中の肝ともいえる部分である。これは評価関数 get_fitnessからの戻り値(特徴数、スコア）の評価方法を設定しており、-1は低い方をよく評価し、1は高い方を評価するという意味である。つまり特徴数が少なく精度が高いほどを優秀な遺伝子として残していくわけである（ここでようやく遺伝的アルゴリズムっぽい単語が登場）。このことが理解できれば、このアルゴリズムを別の問題にも適用できるだろう。

結果の出力

algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2,ngen=50, stats=stats, halloffame=hof,verbose=True)でアルゴリズムが実行される。終了したらhof(Hall Of Fame=殿堂入り)に、よさげな最適解がいくつか格納されているので、今回は全ての最適解についてその最適解に特徴が絞り込まれたファイルを出力するようにした。ファイル名はoutputオプションで指定したで指定したパスに、特徴数とスコアを刻印したファイル名としている。

元ネタからの変更箇所

元ネタからの変更箇所は以下の通り。
-RidgeDVであるが、normalizeされていないデータをinputと想定しているため、normarized=Trueのオプションを付与した。
-"mate", "select"の指定は元ネタの通りtools.cxTwoPoint, tools.selNSGA2だと正解率が低く特徴数が減らなかったため、参考サイトから変更した。

その他の箇所

自分の勉強不足のため解説できる状態ではないが、「選択」「交差」「突然変異」等の具体的なやり方を指示しているのだろうと思われる。

実行方法

以下の通り実行

python GAFeatureSelection.py -input test\train.csv -output test\output

オプションの説明
- input 入力ファイルを指定。csv形式を前提。1列目はID、2列目は目的変数、3列目以降を説明変数と仮定。
- output 出力ファイルを指定。実査には、これに特徴数とスコアを刻印されたcsvファイルが最終出力となる。

実行例

化合物ライブラリRDKitに付属のデータを使って生成した1562件の特徴をもつ1200件のデータに対し、30個体、500世代の計算を実行したところ、以下の9件の特徴数とスコアの組み合わせのものが殿堂入りしたものとなった。計算時間は、約6～7時間を要した。

525.0,0.8213247480199452
Writing Output File With 525 remaind columns
527.0,0.8325450328562906
Writing Output File With 527 remaind columns
528.0,0.8356185010777725
Writing Output File With 528 remaind columns
537.0,0.8361035885472596
Writing Output File With 537 remaind columns
538.0,0.8438213291722809
Writing Output File With 538 remaind columns
558.0,0.8448773229193409
Writing Output File With 558 remaind columns
740.0,0.8464607450447291
Writing Output File With 740 remaind columns
755.0,0.8471958542235395
Writing Output File With 755 remaind columns
765.0,0.853716386001891
Writing Output File With 765 remaind columns

初期個体の特徴数は、乱数処理によって実質約半分(=今回781件)になるため、そこら約200件程度しか特徴数が減っていない。
しかも精度がさがっている。。。（涙）

所感

他のデータでも色々やってはみたものの、中々特徴選択のメリットを享受した体験がなく、自分の中で評価手法が確立していない。サンプル数が少なく、線形回帰手法など大量に特徴を与えるとオーバーフィッティングしやすい手法を使う場合、効果はあるとは思われる。あと元から大量のノイズのあるデータにも効果はあるだろう。自分がよく取り組んでいる化合物データでは、そもそも化学構造から特徴を生成するケースが大半であり、それ以外に得られるデータも少ないため、使える特徴は全て使いたいことが多く、ノイズが多いケースは少ないだろう。逆に大量のデータがあって、勾配ブースティングやDeepLearnigを使う場合などは、ほぼ全ての特徴を使うとういこともありうるだろう（あくまで個人の経験則だが)。いずれよいベンチマークを探して評価手法を確立したい。また説明できなかったところも含め、DEAPをもっと調査し、別の進化的アルゴリズムであるPSOによる特徴抽出も試してみたい。

DEAPで遺伝的アルゴリズムによる特徴選択をしてみた