都内でデータサイエンティストしている者です。はじめてのアドベントカレンダーの投稿となります。今年もはやいものですね。
よければ、暇つぶしにどうぞ。
逆解析の何がウレシイの！？
皆さん、日々機械学習を用いて、データからあらゆるラベルであったり、数値であったりを予測していると思います。私も業務で、スコアリングだったり、売上だったり、様々指標を教師あり学習器を用いて、予測しています。しかし、予測値を出してプロジェクトが終わりということはなかなかありません。特に多いのは、なぜその予測値が出るのか？という説明性を求められるケースです。さらには、ほしい予測値を出すことはできるのか？といったケースもあります。この記事では後者に絞った話をします。
まず初めに、哲学的な問を投げさせてください。そもそも、モデル理解するとは何でしょう？何がどうなれば、理解したと言えるのでしょう？これにはいくつか主張はあるかと思いますが、ここでは「モデルが出す出力をコントロールできるようになること＝モデルを理解していること」と考えます。例えば、あるユーザーが将来優良顧客になる確率を予測するモデルがあったとします。そのモデルはユーザーの属性や過去の行動を元に、「0.3」と予測したとします。では、このユーザーのスコアを「0.8」に上げたい場合は、どうすれば良いのでしょう？ユーザーの属性や行動には現実的に変えられないもの（性別、登録時期）と変えられるもの（１日あたりのPV数）があります。ユーザーの、変えられる特徴量を上手く変更して、0.8まで上げたいとすると、どの特徴量がどうなれば、このユーザーのスコアは0.8になるのでしょうか？
上ではスコアリングの例を挙げましたが、似たような問は例えば売上の予測のような数値の予測でも起きます。例えば、売上の予測値が100と出たとします。しかし売上目標は120となっており、このままだと未達の可能性があります（こういったプロジェクトを着地点予測といったりします）。さて、目標の120が必要な場合、何をどうすれば、120に改善できるのでしょうか？
上は、機械学習モデルの例ですが、例えば、数理最適化等で物流の効率化などに取り組んでいる例でも、適用可能です。数理最適化でトラックの配送ルートを決めたとしましょう。ですが、現実では、よりコストインパクトがほしくロジックで出した答えよりも、良い答えが欲しい場合もあります。こういった場合でも逆解析が威力を発揮します。
いかがだったでしょうか？こういったシーンは、現実のプロジェクトでは良く起きると思います。特にビジネスサイドがデータサイエンティストに意見を求めたいケースはこういったケースではないでしょうか？しかし、意外にもこのような問にダイレクトに答える方法はあまり世の中に示されていないように思います。各変数の貢献度だけでは、上のような問にダイレクトに答えることは難しいように思います。
逆解析はデータサイエンティスト的には自分が作成したモデルを理解する良い手助けにもなりますし、ビジネスサイドからしても、打ち手につながります。私は、逆解析こそが「攻めのIT」に相応しい技術かなと思ったりしています。
逆解析とブラックボックス最適化について
前置きが長くなりましたが、ここからは技術的な話をします。勘の良い方は、この話は、モデルの入力と出力を逆転させたいという話だということに気づかれたのではないでしょうか？X→yを予測する機械学習モデルがあったとして、今の状況はy→Xを出したいのです。y = f(X)というモデルがあったとして、そのモデルのパラメータをいじることはせずに、欲しいyからそのyをできるだけ正しく出力するXを出したいという話です。
面白いことに、これは機械学習の世界を離れ、数理最適化問題の枠組みで定式化することができます。
左側のf^{-1}は関数fの逆関数を指しています。f(x)=x^2とかだとf^{-1}=√xですね。改めてですが、欲しい答え（y）からそれを再現するXを出したいのです。次に右側です、数理計画法として記述されています。こちらではyは固定値（定数）で\hat{y}および、Xが変数です。要するに、Xをいろいろ試して、\hat{y}をいろいろ出力し、欲しいyとの誤差が0になれば、そのときのXが逆関数の出力だと理解できます。
どうしてわざわざ、数理計画法の形式で書くのかと思った方もいるかもしれません。\min: |y - f(X)|でいいじゃん！というツッコミも来そうなところです。数理計画法として書くことでXに対する制約条件を追加することができます。↑に上げた例でいえば、ユーザーの性別、年齢といった属性値は変えることはできません。そういった変えられない変数に対しては、x_{age}=30のように変数を固定してしまえば良いのです。
さらに、f(X)の関数系は特別指定していません。fはランダムフォレストでもLightGBMでもNNでもなんでも良いのです。Xを入れたらyが出てくれば関数系は問いません。このように関数系を問わない最適化をブラックボックス最適化と呼んだりするようです（詳しいことは専門家ではないので知りません。私はあくまで実務家です。）
お試し
ではでは、簡単に↑のブラックボックス最適化問題を解くコードをご紹介します。ブラックボックス最適化ができるライブラリはいくつかあるようですが、私は使い慣れているLocalsolverを使いました。ちなみに、こちらは有償ツールとなります。無料ツールであれば、hyperopt, Kurobako, といったいくつかあるようです。（これらは使ったことがないのでわかりません。）
ここでは、簡単のため、制約条件は入れずに解きます。以下のような、整数値、バイナリ値、連続値が入り混じった変数を10列用意し、連続値yを予測するモデルをランダムフォレストを用いて作成しました。4行しか表示していませんが、実際には10000行程度、用意しました。
逆解析はここからが、本番です。この学習済みのランダムフォレストを固定したまま、今度は欲しい出力の一覧を用意し、ランダムフォレストにこの出力を出させてみましょう。与えるインプットは以下のデータです。
なかなか見慣れないインプットではないでしょうか？普通は逆を思い浮かべませんか？Xが埋まっていて、yが空欄かと思いきやその逆で、yだけが埋まっています。このyを出力するXを埋めるというのがやりたきことです。入出力が逆転しているので、逆解析と呼ばれます。逆解析をしているpythonコードを簡単に紹介します。
import localsolver
import numpy as np
import pandas as pd
import logging, os
from sklearn.ensemble import RandomForestRegressor
def solve_inverse_machine_learning(regressor,
y_target,
list_idx,
list_cname,
cname_map_range,
cname_map_type,
test_size,
ls_time_limit=10):
"""ブラックボックス最適化を用いて、逆関数値を出力する
:param regressor: 機械学習モデル
:param y_target: 欲しい出力のリスト
:param list_idx: Xの行のインデックス
:param list_cname: Xの列のインデックス
:param cname_map_range: 列ごとの取り得る値の範囲
:param cname_map_type: 列ごとの型（整数orバイナリor連続）
:param test_size: y_targetの長さ
:param ls_time_limit: 計算時間
"""
def _predict(args):
"""args must be 1-d list
"""
X = np.array(args).reshape(test_size, len(list_cname))
res = np.abs(y_target - regressor.predict(X))
return res.sum()
result = []
with localsolver.LocalSolver() as ls:
model = ls.model
X = []
idx_cname_map_var = {}
for idx in list_idx:
for cname in list_cname:
type_ = cname_map_type[cname]
if type_=='B':
v = model.bool()
elif type_=='I':
v = model.int(cname_map_range[cname][0],
cname_map_range[cname][1])
elif type_=='C':
v = model.float(cname_map_range[cname][0],
cname_map_range[cname][1])
else:
raise Exception('invalid column type')
v.name = cname + '_%s'%idx
X.append(v)
idx_cname_map_var[idx, cname] = v
objective = model.double_external_function(_predict)
objective.external_context.lower_bound = 0.
model.minimize(objective(X))
model.close()
ls.param.set_time_limit(ls_time_limit)
ls.param.set_log_file('./log_%d.log'%len(list_cname))
ls.solve()
sol = ls.solution
for x in X:
result.append(sol.get_value(x))
return result
if __name__=='__main__':
TIMELIMIT = 60 # * 60
NB_COLS = [10] # [10**1, 10**2, 10**3]
INTEGER_COL_RATE = 0.3
BINARY_COL_RATE = 0.1
FLOAT_COL_RANGE = (0., 10.)
INTEGER_COL_RANGE = (0, 5)
TRAIN_SIZE = 10**4
TEST_SIZE = 10**2
Y_SCALE = 100
DATA_DIR = './data'
OUTPUT_PATH = os.path.join(DATA_DIR, 'inverse_ml.xlsx')
nbcols_df_map = {}
np.random.seed(3655)
for nb_cols in NB_COLS:
logger.info('start make datamart %s'%nb_cols)
integer_cols = int(nb_cols * INTEGER_COL_RATE)
binary_cols = int(nb_cols * BINARY_COL_RATE)
float_cols = nb_cols - integer_cols - binary_cols
assert integer_cols+binary_cols+float_cols==nb_cols
list_float_cname = ['x_%s'%i for i in range(1, float_cols+1)]
df_float = pd.DataFrame(np.random.random_sample((TRAIN_SIZE, float_cols)) * FLOAT_COL_RANGE[1],
columns=list_float_cname)
list_int_cname = ['x_%s'%i for i in range(1+float_cols, 1+float_cols+integer_cols)]
df_int = pd.DataFrame(np.random.randint(INTEGER_COL_RANGE[0], INTEGER_COL_RANGE[1] + 1,
(TRAIN_SIZE, integer_cols)),
columns=list_int_cname)
list_binary_cname = ['x_%s'%i for i in range(1+float_cols+integer_cols, 1+float_cols+integer_cols+binary_cols)]
df_binary = pd.DataFrame(np.random.randint(0, 2, (TRAIN_SIZE, binary_cols)),
columns=list_binary_cname)
list_cname = list_float_cname + list_int_cname + list_binary_cname
datamart = df_float.join(df_int).join(df_binary)
# datamart['y'] = np.random.randint(0, 2, TRAIN_SIZE) # set binary target
datamart['y'] = np.random.random_sample(TRAIN_SIZE) * Y_SCALE # set score target
logger.info('end make datamart %s'%nb_cols)
cname_map_type = {}
cname_map_range = {}
for c in list_cname:
if c in list_float_cname:
v = 'C'
range_ = FLOAT_COL_RANGE
elif c in list_int_cname:
v = 'I'
range_ = INTEGER_COL_RANGE
elif c in list_binary_cname:
v = 'B'
range_ = None
else:
raise Exception('invalid cname')
cname_map_type[c] = v
cname_map_range[c] = range_
X, y = datamart.iloc[:, :nb_cols], datamart['y']
logger.info('start fit regressor %s'%nb_cols)
# regressor = LogisticRegression()
regressor = RandomForestRegressor(n_estimators=10)
regressor.fit(X=X, y=y)
logger.info('end fit regressor %s'%nb_cols)
logger.info('start optimize %s'%nb_cols)
y_target = np.random.random_sample(TEST_SIZE) * Y_SCALE # set target score
list_idx = list(range(1, y_target.shape[0] + 1))
res_X = solve_inverse_machine_learning(
regressor=regressor,
y_target=y_target,
list_idx=list_idx,
list_cname=list_cname,
cname_map_type=cname_map_type,
cname_map_range=cname_map_range,
test_size=TEST_SIZE,
ls_time_limit=TIMELIMIT
)
logger.info('end optimize %s'%nb_cols)
logger.info('start draw errors %s'%nb_cols)
res_X = np.array(res_X).reshape(TEST_SIZE, nb_cols)
# y_pred_array = regressor.predict_proba(res_X)
y_pred_array = regressor.predict(res_X)
df_result = pd.DataFrame(res_X, columns=list_cname)
df_result['y_target'] = y_target
df_result['y_pred'] = y_pred_array
df_result['y_error'] = np.abs(y_target - y_pred_array)
nbcols_df_map[nb_cols] = df_result
logger.info('end draw errors %s'%nb_cols)
writer = pd.ExcelWriter(OUTPUT_PATH, engine='xlsxwriter')
for nb_cols, df_result in sorted(nbcols_df_map.items(), key=lambda x:x[0]):
df_result.to_excel(excel_writer=writer, sheet_name='result_%d'%nb_cols, index=False)
writer.save()
def _predict(args): となっているところが目的関数を表しています。resには、欲しい出力（y_target）とモデル出力f(X)の絶対誤差が格納されています。この絶対誤差をできるだけ小さくしようとしています。model.minimize(objective(X))で実際に目的関数に設定しています。
また、Xには連続値とバイナリと整数が入り混じっているので、その型を定義しているのが、v = model.bool()といった箇所になります。
「cname_map_type」で、列名に対して、その列の型が返ってきます。
さてさて、このコードを実行すると、以下のように表が見事に埋まってくれました。
y_targetは元々欲しかったyで、y_predはf(X)の返り値となります。y_errorは|y_target-y_pred|です。完全一致していませんが、十分に近い値になっていることが見て取れます。これにより、めでたく、y=f(X)の逆関数値が完全ではないですが、求まったことになります。
ビジネスシーンで適用するならば、y_targetがユーザーの何かしらのスコアを表しているとし、x_iが何かしらの行動を指しているすれば、1行目のユーザーで言えば、欲しいスコア71.6点のスコアが欲しいとするならば、x_8（何かしらの回数）を4回にすればよいといった示唆が出せます。
今回は全ての値をフリーで解きましたが、x_10が性別を表しているとすれば、x_10=1と値を固定した上で解くことも可能です。（属性は変えずに、行動量を変化させた上で欲しいスコアを導き出すことができます。）
まとめと注意点
この記事では、機械学習プロジェクトのビジネスシーンで求められる例を挙げ、それを解決する一つの方法として逆解析を紹介しました。逆解析は、数理計画法として自然に定式化することができ、その数理計画問題を解くことで、近似的に逆関数値を出せます。もしかすると、機械学習の専門家から、いろいろとツッコミが来るのかもしれませんが、私は論文を書いたこともなければ、学会発表もしたことがない身ですので、専門的な背景は深くは知りません。本業はあくまでクライアント業ですので、クライアント業をしていく中で、自然と、こういったのあれば良いやん！と思い立った次第です。ぜひとも、逆解析が盛り上がってほしいなと思います。
それでは、良い年末を！^^
参考
以下の記事を参考にしました。
機械学習モデルを逆解析する