Keras / TensorFlowで線形回帰・非線形回帰をやってみる

以前、与えられた標本点を補間するベジェ曲線を求めるためにTensorflowを使うという記事を書いた。

記事1. TensorFlowを使って、点列を補間するベジェ曲線を求める（その１）
記事2. TensorFlowを使って、点列を補間するベジェ曲線を求める（その２）
記事3. TensorFlowのOptimizerを比較する（ベジェ曲線編）

これらは、TensorFlowを身近な問題に適用してAPIの使い方に慣れてみようという趣旨だったので、TensorFlowのtf.Graph、tf.Tensor、tf.Variableなど（現在では低水準APIと呼ばれるクラス群）を使って課題を実装した。ところが。

TensorFlowサイトの方針が変わった

現在、TensorFlowのサイトは大きく様変わりしてしまっている。Tutorialのページは、Tensorflowで実装されたKerasツール（高水準APIと呼ばれる）を使ったコード例で占められ、低水準APIのコード例はどこかへ行ってしまった。

これは、TensorFlowの主張に違いない。Nueral NetworkのユーザはKerasなどの高水準APIを使うべきで、低水準APIは知らなくても良いということだろう。

Keras APIを使って実装してみよう

これからがKerasの時代ならば、KerasのTensorFlow実装に慣れなければならない。そこで、線形回帰と非線形回帰の問題を、Kerasで実装してみることにする。サンプルは記事3の中から取ることにしよう。

技術的には、次のような事柄が登場する。

• Sequential Modelの使い方（線形回帰）
• Layerのカスタマイズの方法（非線形回帰）

線形回帰

図1は記事3では図21にあたる。Bezier曲線を求めるための最適化法としてFtrl法を用いたとき、学習係数と収束率の関係を実際のデータ（赤丸）から取得したものである。目的は青線を求めることである。これは放物線（2次曲線）に見えるが、非対称なので2次ではなく3次曲線である。

従って問題は、3次式

y(x;\{a_k\}) = \sum_{k=0}^3{f_k(x)a_k} \tag{1}

f_k(x) = x^k \tag{2}

の係数$\{a_k\}$を標本点との距離の総和（を標本点の総数$N$で割ったもの）

D(\{a_k\}) = \frac{1}{N}\sum_{i=0}^{N-1}{|y(x_i;\{a_k\}) - y_i|^2} \tag{3}

を最小化するように係数を決めることになる。ここで(1)式は$\{a_k\}$に対する1次式になっている（$x$に対するのではなく）から、線形回帰問題である。

線形回帰のKerasによる表現

単純パーセプトロンから活性化関数とバイアス項を取り除くと、単なる線形回帰問題になる。そのため、Kerasで表わすと、とても単純なコードになる。

リスト1

0001: def main():
0002:   d = np.loadtxt("入力ファイル")
0003:   x = d[:, 0]
0004:   inputs = np.array([x * x * x, x * x, x, x / x], np.float32).transpose()
0005:   labels = d[:, 1].reshape(-1, 1)
0006:
0007:   model = tf.keras.Sequential()
0008:   model.add(tf.keras.layers.Dense(1, use_bias=False))
0009:   model.compile(optimizer=tf.train.AdamOptimizer(1), loss='mse')
0010:
0011:   history = model.fit(inputs, labels, epochs=5000000)
0012:   mses = history.history['loss']
0013:   param = model.get_weights()[0]
0014:   print(mses[-1], param)

なお、例によって次のimportが必要である。

リスト2

import numpy as np
import tensorflow as tf

また、入力ファイルは以下のようなデータである（y軸は10^5倍している）。

リスト3

0.95 1.44685
1.0 1.60368
1.05 1.73257
...
1.55 1.09186

線形回帰のコードの説明

• 0004: inputs = np.array([x * x * x, x * x, x, x / x], ...).transpose()**

リスト3の元入力の第1フィールド${x_i}$を、(2)式の4元ベクトル

\boldsymbol{f}_i = (x_i^3, x_i^2, x_i, 1)

に変換して、入力データとする。transpose（転置）した結果、inputs の shape は（N, 4) になる。最初の次元はバッチ次元であり、Nは標本点の個数である。なお、$1$ではなく$ x / x$としているのは、各ベクトル要素のshapeを揃えるテクニックである。

• 0005: labels = d[:, 1].reshape(-1, 1)

labels の shape を (N, ) から (N, 1) に変換している（バッチ次元とデータ次元に分離）。

• 0007: model = tf.keras.Sequential()

tf.keras.Model は、低次元APIでは Graph に相当する。Modelの派生クラスであるtf.keras.Sequentialは、入力層から出力層へ直列に並べるモデルである。構造の単純な Neural Network は、全て Sequential で書くことができる。

• 0008: tf.keras.layers.Dense(1, use_bias=False)

直列に並んだ Neural Network の各層は、tf.keras.layers.Layer の派生クラスとして表現する。このうち、tf.keras.layers.Dense は、全結合層を表現する派生クラスである。第1引数の1は出力側次元（データ次元。バッチ次元は除く）を表す。use_bias引数は、バイアス項が無いことを示す。活性化関数は activation引数で指定するが、ここでは未指定なので活性化関数は適用されない。

• 0008: model.add(...)

各層を前段から後段にむけて付加する。最後に加えた層が出力層である。線形回帰では、出力層だけしかない。

• 0009: model.compile(optimizer=..., loss='mse')

最適化方法（optimizer引数）とLoss関数（loss引数）を指定し、モデルの構造を確定する。'mse' は平均二乗誤差関数を表す。(3)式がまさに平均二乗誤差なので、これだけでOK。ここでは最適化関数として、学習係数を1としたAdamOptimizerを指定している。

• 0011: history = model.fit(inputs, labels, epochs=5000000)

model.fit関数で教師付き学習を行う。入力がinputs、教師データがlabelsである。epochsでバッチ処理の繰り返し回数を指定する。fit()の戻り値はHistoryオブジェクトである。

• 0012: mses = history.history['loss']

fit()の戻り値であるHistoryオブジェクトには、学習過程のepochごとのLoss関数の値が保存されている。

• 0013: param = model.get_weights()[0]

最適化パラメータの値を求める。get_weights()の戻り値は、各層毎に最適化パラメータのリストになっている。このModelでは、最初の要素だけに意味がある。

実用的価値はないが

以上のように線形回帰はKerasを使うと、とても簡単に書ける。といっても、実用的価値は全く無い。線形回帰はNumPyのpolyfit()を使うのが良い。しかも、Kerasの最適化では500万ステップ必要とする精度に対して、1ミリ秒で同じ精度を得ることができる。これは線形回帰が解析的に解けるためである。

リスト4

p = polyfit(d[:, 0], d[:, 1])  # dはリスト1参照

非線形回帰

図2は記事3の図13である。Bezier曲線を求めるための最適化法としてFtrl法を用いたとき、学習係数と収束率の関係を実際のデータ（青三角）から取得したものである（面白いことに、Adagrad法の結果と完全に一致する）。目的は黒い曲線を求めることである。これは、次の式で表される。

\mu = a (1 - \exp{(- b \eta^2)}) \tag{4}

$a$と$b$の2個のパラメータを決定したい。

非線形回帰のKerasによる表現

リスト5

0001: def main():
0002:   d = np.loadtxt(入力ファイル)
0003:   inputs = d[:, 0].reshape(-1, 1)
0004:   labels = d[:, 1].reshape(-1, 1)
0005:   param = initParam(inputs, labels) # 適当なパラメータの初期値を計算する
0006:
0007:   model = tf.keras.Sequential()
0008:   model.add(NonLinearLayer(param, d.shape[0]))
0009:   model.compile(optimizer=tf.train.AdamOptimizer(1e-4), loss='mse')
0010: 
0011:   history = model.fit(inputs, labels, verbose=0, epochs=2000)
0012:   mses = history.history['loss']
0013:   param = model.get_weights()[0]
0014:   print(mses[-1], param)

リスト5を線形回帰のリスト1と比較して頂きたい。ほとんど同じ形だが、8行目のLayerクラスが、DenseからNonLinearLayerに変わっている。NonLinearLayerクラスは、式4の非線形回帰のためにカスタマイズしたLayerである。5行目のinitParam()は、標本点の座標からパラメータの初期値（NumPy配列）を決める関数（内容は省略する）。初期値が適当でないと最適化が失敗してしまう。

Layerのカスタマイズ

リスト6

0031: class NonLinearLayer(tf.keras.layers.Layer):
0032:   def __init__(self, param, batch_size):
0033:     self.initializer = tf.keras.initializers.Constant(param)
0034:     self.shape = param.shape
0035:     super(NonLinearLayer, self).__init__(batch_size=batch_size)
0036:
0037:   def build(self, input_shape):
0038:     self.param = self.add_weight('param', self.shape, initializer=self.initializer)
0039:
0040:   def call(self, inputs):
0041:     return self.param[0] * (1 - tf.exp(- self.param[1] * inputs * inputs))
0042: 
0043:   def compute_output_shape(self, input_shape):
0044:     return input_shape

tf.keras.layers.Layerの派生クラスを作る際に必要なメソッドは以下である。

build(self, input_shape)

最適化パラメータの登録など。低次元APIのtf.Variableの登録に相当。一回目のfit()から呼ばれる

call(self, inputs)

低次元APIの言葉では、inputsを入力側tf.Tensorとして、戻り値を出力側tf.Tensorとする。実際の計算ではなくtf.Graphへの登録であるため、一度呼ばれるだけ。一回目のfit()から呼ばれる

compute_output_param(self, input_shape)

一回目のfit()から、build()の後、call()の前に呼ばれる。出力側のshapeを戻す

get_config(self)

Layerのインスタンスに必須であるプロパティをディクショナリにしたものを戻す。シリアライズのために必要。この例には必須ないので、定義していない

NonLinearLayerのコードの説明

• 0033: tf.keras.initializers.Constant(param)

最適化の最初にパラメータの初期値を決めなければならない。初期化を行うのが Initializerである。Initializer は、add_weight() メソッドの実行時に必要になる。

全結合層を表現するDenseを例にあげると、Denseで使用する Initializer には、線型結合部分の kernel_initializer と、バイアス部の bias_initializer の2種類がある。Dense の場合、2種類の Initializer を指定しないとき、kernel_initializer には決められた範囲の一様分布で値を決める関数がセットされ、bias_initializer にはゼロで初期化する関数がセットされる。

NonlinearLayerの場合は、2個の変分パラメータには、ランダムな値ではなくリスト5の5行目であらかじめ決めておいた値を初期値としてセットしたい。そうでないと最適化が失敗するからである。このように定数（の配列）をセットする Initializer が Constant() である。

Constant Initializer は、パラメータを決まった値に初期化したい場合には、どんな場合にも使える。

• 0038: self.param = self.add_weight('param', self.shape, initializer=self.initializer)

Modelに対して、最適化パラメータを登録する役目を担うのが、ベースクラスであるLayerクラスで定義された add_weight() メソッドである。内部では、低水準APIの tf.Variable を使っているのだろう。

最初の2引数、レイヤ名とパラメータのshapeは必須である。それ以外のオプション引数のうち、ここでは23行で作成した Constant Initializer のインスタンスを initializer引数にセットする。

これも、あまり実用的ではないが

非線形回帰についても、あまり実用的ではない。Keras を使わなくても、SciPyのcurve_fit() を使えば数行で実現することができる（curve_fit() においても、最適化パラメータに適切な初期値を与えないと最適化が失敗する）。

リスト7

def func(x, a, b):
  return a * (1 - np.exp(- b * x * x))

p0 = initParam(d[:, 0], d[:, 1])　# 適当な初期値を決める。dはリスト5を参照
param, pcov = scipy.optimize.curve_fit(func, d[:, 0], d[:, 1], p0=p0)
print(param)

とはいえ、Keras/TensorFlow の最適化アルゴリズムは curve_fit() よりも優れているように見えるので、パラメータの数が多い場合には curve_fit() よりも効果的になるだろう。

少し嘘をついていました

リスト1とリスト5では、コードの見通しよくするために、少し嘘をついていた。というのは、学習過程を表現する11行目の fit() メソッドは、必ずしもLoss関数の最小値を与えないからである。1 epochごとのLoss関数の値は図3のように激しく揺れ動く。おそらく、この現象は過学習を避けるアルゴリズムのためだろう。

実際に回帰を行うためには、リスト1とリスト5の両方とも、9行目および13行目以降をリスト8のように変更する必要がある。

リスト8

0009:   model.compile(..., callbacks=[mycallback]) # キーワード引数callbacksの付加
...
0013:   params = history.history['param']
0014:   im = None
0015:   mm = None
0016:   for i in range(len(mses)):
0017:     if mm is None or mm > mse:
0018:       mm = mses[i]
0019:       im = i
0020:   
0019:   print(i, mm, params[im])
0020:
0021: mycallback = ParamKeeper()
0022: class ParamKeeper(tf.keras.callbacks.Callback):
0023:   def on_epoch_end(self, epoch, logs=None):
0024:     param = self.model.get_weights()[0]
0025:     history = self.model.history.history
0026:     history.setdefault('param', []).append(param.tolist())

終わりに

記事1と記事2では、TensorFlowの低水準APIを使ってテーマを実装した。そこで気付いたのは、tf.Graphへのパラメータ、演算、最適化方法の登録やtf.Sessionを使っての初期化、学習、計算値の取得などが、あまりにも定型的に過ぎて、毎回、ほとんど同じコードになることだった。

その意味では、Kerasの表現方法はとても合理的で、不必要な部分を隠蔽してくれる。Kerasでは書けない部分も、そこだけをカスタマイズすればよいから、全体としての見通しがとても良い。Neural Networkのプログラミングの敷居が、一段階低くなったのは間違いない。

とはいえ、低水準APIの知識が全く必要ないかというとそうではない。あまり一般的でないネットワークを構成する際には、ベースクラス（tf.keras.layers.Layer や tf.keras.initializers.Initializer）や、既製の派生クラス（tf.keras.layers.Dense や tf.keras.initializers.Constant）のコードを読む必要がある。そのためには低水準APIの知識が必須になる。TensorFlowのサイトに低水準APIのチュートリアルが無くなってしまった（あるのかも知れないが探せていない）のは、少し残念だ。