【機械学習を基本から丁寧に】TensorFlow Eager Executionで多項式回帰の実行

スポンサーリンク

はじめに

前回の記事では単回帰で直線をフィッティングする問題を、TensorFlow Eager Executionで実行しました。 割と低レベルなTensorFlowの機能を使ったため、数式とほとんど直結した形で機械学習の動作を確認できたのではないかと思います。

www.hellocybernetics.tech

実は先のような単回帰をしっかりマスターすると、基本的にどんなモデルを仮定するのかというのがちょっと変わるだけで、いろいろな手法を扱うことができることが分かってくると思います。

今回はそのことを見るために、前回の記事内容とコードを少し改変して多項式回帰を行ってみます。

実行環境

Google colabで実行しました。

tensorflow 1.9.0rc1 matplotlib 2.1.2

理屈編

問題設定

ここでは最も単純な線形回帰である単回帰を、少し発展させた多項式回帰について説明します。 これは、入力変数がスカラーで x$x$、出力変数もスカラーで y$y$ のときに、 x$x$ と y$y$ の関係を多項式で表そうというものです。 データセット D=(x1,y1),(x2,y2),⋯,(xN,yN)$D=(x_1,y_1),(x_2,y_2),\cdots ,(x_N,y_N)$がある時に（ xi,yi$x_i,y_i$ はそれぞれスカラー）、多項式回帰を用いるというのは、入力 x$x$ と出力 y$y$ の間に

y=a1x+a2x2+a3x3+⋯+aKxK+b=∑k=1Kakxk+b$$y=a_{1}x+a_2{x}^2+a_3{x}^3+\cdots +a_K{x}^K+b=\sum _{k=1}^K{a}_k{x}^k+b$$

の関係があるとしてa1,a2,a3,⋯,aK$a_1,a_2,a_3,\cdots ,a_K$とb$b$を求めようというものになります。ここでK$K$は設計者が決める自然数です。 言うまでもなく、このK$K$の設定次第で、表しうる多項式の複雑さが変わってきます（k$k$は飛び飛びの自然数でも構わない。とにかく多項式で表して回帰をするんだということ）。

この場合には、これらのデータは、どのデータを見ても概ね yi=∑Kk=1akxki+b$y_i=\sum _{k=1}^K{a}_k{x}_i^k+b$ となっていると考えるわけですね。これを全てのデータでそれ相応に成り立たせるようなa1,a2,a3,⋯,aK,b$a_1,a_2,a_3,\cdots ,a_K,b$ を求めることにします。 ここで「それ相応」という意味については、単回帰のときと同様の設定で行きます。

損失関数

さて、 (xi,yi)$(x_i,y_i)$ というデータのズレ具合は li={yi−(∑Kk=1akxki+b)}2$l_i=\{y_i-(\sum _{k=1}^K{a}_k{x}_i^k+b){\}}^2$ と表すことにします。仮にズレが全く無かったのであれば li=0$l_i=0$ となることに注意してください（なぜなら私達の頭の中では、 yi=∑Kk=1akxki+b$y_i=\sum _{k=1}^K{a}_k{x}_i^k+b$と期待しているのですから ）。 私達がこれから獲得したい多項式というのは、このようなズレがデータ全体で見た時に少ないことが望まれます。

そこで、これらのデータのズレの総和を取りましょう。

L(a1,⋯,aK,b)=∑iNli=∑iN{yi−(∑k=1Kakxki+b)}2$$L(a_1,\cdots ,a_K,b)=\sum _i^N{l}_i=\sum _i^N\{y_i-(\sum _{k=1}^K{a}_k{x}_i^k+b){\}}^2$$

これが小さくなれば、データ全体考慮した上でそれ相応の線を引けた（ a1,a2,a3,⋯,aK,b$a_1,a_2,a_3,\cdots ,a_K,b$ を決められた）ことになります。 これを損失関数と呼び、L(a1,⋯,aK,b)$L(a_1,\cdots ,a_K,b)$ が大きいほど損失が大きく、うまく一次式が求められていないという話になります。

損失を小さくする勾配法

損失関数 L(a1,⋯,aK,b)$L(a_1,\cdots ,a_K,b)$ を小さくするような a1,a2,a3,⋯,aK,b$a_1,a_2,a_3,\cdots ,a_K,b$ をどのように求めればいいのでしょうか？ 微分積分学を知っている人ならば、 「a1$a_1$ で微分して 0$0$ とおく、…aK$a_K$で微分して0$0$とおく、 b$b$ で微分して 0$0$ とおく」をやって連立方程式を解けば良さそうだとわかるかもしれません。

いずれのパラメータを使おうと、パラメータak$a_k$について統一的に、

ak←ak−ϵ∂L(a1,⋯,aK,b)∂ak$$a_k\leftarrow a_k-ϵ\frac{\mathrm{\partial }L(a_1,\cdots ,a_K,b)}{\mathrm{\partial }{a}_k}$$

と表すことができます。

実践編

必要なライブラリの準備

今回はEager Executiuonを使って説明していきます。 import tensorflow.contrib.eager as tfeを忘れないようにしましょう。また、tf.enable_eager_execution()によって、以降のコードをすべてEagerモードで実行するようになります。

import numpy as np
import tensorflow as tf
import tensorflow.contrib.eager as tfe
import matplotlib.pyplot as plt
import seaborn as sns
tf.enable_eager_execution()

問題設定

今回は

y=x3+2x2−4x−1+2ϵ$$y=x^3+2x^2-4x-1+2ϵ$$

ϵ∼N(0,1)$$ϵ\sim \mathcal{N}(0,1)$$

というデータを観測したことにします。 実際には(x,y)$(x,y)$にはy=x3+2x2−4x−1+2ϵ$y=x^3+2x^2-4x-1+2ϵ$の関係があるのだが、y$y$を観測するときに平均0$0$分散1$1$のガウスノイズϵ$ϵ$が乗って更にノイズが２倍になってしまっているようなケースを想定していることになります。

もちろん、今回は「このデータが本来はy=x3+2x2−4x−1+2$y=x^3+2x^2-4x-1+2$なのだ！と判っている」前提で話を進めていますが、現実は「データだけが手元にあり、本来はどのような関係なのかわからない」状況であることは注意してください。 わからないからこそ、機械学習や統計解析によって関係性を知ろうということ行うのです。

def toy_polynomial_data():
  x = np.linspace(-3, 3, 50)
  y = x**3 + 2*x**2 - 4*x + 2*np.random.randn(50) - 1
  return x, y

x, y = toy_polynomial_data()
plt.scatter(x, y)
plt.show()

モデルの設計と損失関数

さて、上記のデータの関係性を明らかにするためのモデルの設計を行います。 ここではEagerモードでのモデルの書き方を簡単に説明しながら進めていきましょう。

全体像としては以下のようになります。先に答えを見ておきましょう。次に１つ１つのメソッドを簡単に説明していきます。

class Model(tf.keras.Model):
  def __init__(self):
    super(Model, self).__init__()
    self.a1 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=1)
    self.a2 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=-1)
    self.a3 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=1)
    self.b = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=2)
  
  def call(self, x):
    return self.b + self.a1*x + self.a2*x**2 + self.a3*x**3
  
  def loss_fn(self, x, y):
    y_pre = self(x)
    mse = 0.5 * (y - y_pre) ** 2 
    return tf.reduce_sum(mse)
  
  def grads_fn(self, x, y):
    with tfe.GradientTape() as tape:
      loss = self.loss_fn(x, y)
      return tape.gradient(loss, [self.a1,
                                  self.a2,
                                  self.a3,
                                  self.b])
    
  def update(self, x, y, lr=0.001):
    grads = self.grads_fn(x, y)
    ## variable.assign_sub(value)
    ## variable -= value
    self.a1.assign_sub(lr * grads[0])
    self.a2.assign_sub(lr * grads[1])
    self.a3.assign_sub(lr * grads[2])
    self.b.assign_sub(lr * grads[3])

tf.keras.Modelクラスでモデルの雛形を作る

基本的に、tf.keras.Modelクラスを継承して使うことになります。 実際には今回の問題はかなり原始的な話題なので、このクラスを使わなくてもいいのですが、 基本的な雛形はこのような形式になるのではないかと思われるので、慣習に習っておきます。

まず、def __init__(self)には今回のモデルで用いるパラメータ（やニューラルネットワークの層）を準備しておきます。 今回必要は a1x1+a2x2+a3x3+b$a_1{x}^1+a_2{x}^2+a_3{x}^3+b$を表すために a1,a2,a3,b$a_1,a_2,a_3,b$ 4つのパラメータを準備しておきます。 tf.contrib.eager.Variable()を利用し、初期値initial_valueを適当に設定しておきます（乱数を使ってもいいですが、今回は敢えて、答えの数値と明らかに違う数値にしておきます）。

次に、call()メソッドは、モデルの実際の計算を書き下します。今回はa1x1+a2x2+a3x3+b$a_1{x}^1+a_2{x}^2+a_3{x}^3+b$を計算するようにしておくだけなので簡単です。

  def __init__(self):
    super(Model, self).__init__()
    self.a1 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=1)
    self.a2 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=-1)
    self.a3 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=1)
    self.b = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=2)
  
  def call(self, x):
    return self.b + self.a1*x + self.a2*x**2 + self.a3*x**3

損失関数

次は損失関数の部分です。

L(a,b)=12∑iN{yi−(∑k=1Kakxki+b)}2$$L(a,b)=\frac{1}{2}\sum _i^N\{y_i-(\sum _{k=1}^K{a}_k{x}_i^k+b){\}}^2$$

というのを書きましょう。y_pre = self(x)で予測であるa1x1+a2x2+a3x3+b$a_1{x}^1+a_2{x}^2+a_3{x}^3+b$を計算します。そして、手元にある実データのy$y$との誤差を図るためにmse = 0.5 * (y - y_pre) ** 2を計算します。 そして、すべてのデータの和を取るためにtf.reduce_sum()を使います。

数式と直結していますので、書くのはそんなに難しくないですね（今は一次元なので簡単ですが、次元が増えていくと、どの次元で和を取るのか等を間違えないようにしなければならないです）。

  def loss_fn(self, x, y):
    y_pre = self(x)
    mse = 0.5 * (y - y_pre) ** 2 
    return tf.reduce_sum(mse)

勾配の計算

パラメータの更新は、微分を計算して、現在の値から引くという操作でしたね。

a1←a1−ϵ∂L(a1,a2,a3,b)∂a1$$a_1\leftarrow a_1-ϵ\frac{\mathrm{\partial }L(a_1,a_2,a_3,b)}{\mathrm{\partial }{a}_1}$$

というものでした。これを実行するには勾配（微分）の計算が必要になります。 今回の問題ならば微分を手で解くことも難しくありませんが、ニューラルネットワークなどのパラメータの微分を、全部て計算で書き下しておくなどやっていられません。 TensorFlowなどの深層学習フレームワークには、指定したパラメータによる微分を自動で計算してくれる仕組みが備わっています。その理論的な背景は「バックプロパゲーション」と呼ばれています。

今回はその理論的な背景には踏み込まず、その機能を有りがたく使わせていただくことにしましょう。 そのコードは以下のようになります。

  def grads_fn(self, x, y):
    with tfe.GradientTape() as tape:
      loss = self.loss_fn(x, y)
      return tape.gradient(loss, [self.a1,
                                  self.a2,
                                  self.a3,
                                  self.b])

忘れずにすべてのパラメータのリストを渡してあげましょう。 返してくれる順番は、引数で与えた順番に一致するので、ちゃんと意識しておきましょう。

パラメータの更新

いよいよパラメータの更新です。

  def update(self, x, y, lr=0.001):
    grads = self.grads_fn(x, y)
    ## variable.assign_sub(value)
    ## variable -= value
    self.a1.assign_sub(lr * grads[0])
    self.a2.assign_sub(lr * grads[1])
    self.a3.assign_sub(lr * grads[2])
    self.b.assign_sub(lr * grads[3])

あとは、このgradsの中身（パラメータの順番）に注意して、更新のコードを書くだけです。 ただし、

self.a -= lr * grads[0]

のような書き方はできません。少し数式の書き方とは離れてしまいますが、assign_subメソッドを使ってください。

ここまでで一通りモデルを書き終えました。再度モデルの全体像を載せておきます。

class Model(tf.keras.Model):
  def __init__(self):
    super(Model, self).__init__()
    self.a1 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=1)
    self.a2 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=-1)
    self.a3 = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=1)
    self.b = tf.contrib.eager.Variable(dtype=tf.float32,
                                       initial_value=2)
  
  def call(self, x):
    return self.b + self.a1*x + self.a2*x**2 + self.a3*x**3
  
  def loss_fn(self, x, y):
    y_pre = self(x)
    mse = 0.5 * (y - y_pre) ** 2 
    return tf.reduce_sum(mse)
  
  def grads_fn(self, x, y):
    with tfe.GradientTape() as tape:
      loss = self.loss_fn(x, y)
      return tape.gradient(loss, [self.a1,
                                  self.a2,
                                  self.a3,
                                  self.b])
    
  def update(self, x, y, lr=0.001):
    grads = self.grads_fn(x, y)
    ## variable.assign_sub(value)
    ## variable -= value
    self.a1.assign_sub(lr * grads[0])
    self.a2.assign_sub(lr * grads[1])
    self.a3.assign_sub(lr * grads[2])
    self.b.assign_sub(lr * grads[3])
      

実験

初期状態のモデル

早速作ったモデルをmodel=Model()でインスタンス化して、実験してみましょう。 まずは学習を全く行わずにデータをぶち込んでみます。

model = Model()
y_init = model(x)

plt.scatter(x, t)
plt.plot(x, y_init, c='r')
plt.show()

学習後のモデル

400回ほどパラメータの更新を行ってみましょう。 損失の変化を記録するためにloss=[]と、モデルの予測の記録を取るためにreg=[]を準備しておきます。 あとはfor文の中でmodel.update(x, y)により学習を進めます。(y_pre=model(x)は予測の記録を取るためだけのコードです)

reg = []
loss = []
for _ in range(400):
    y_pre = model(x)
    reg.append((x, y_pre))
    model.update(x, y, 2e-4)
    loss.append(model.loss_fn(x, y))

図示すると、赤破線が初期のモデルの予測、緑が学習後の予測で、青の点が実データです。

本当に学習は上手く行ったのか

損失の変化に関しては以下のようになっており、ほとんど収束しているように見えます。

plt.plot(range(400), loss)
plt.show()

しかし、これだけでは学習が上手く行ったのかは本当はわかりません。 同様のデータをもう１セット準備しておいて、テストをする必要があります。 今回の手法は、手持ちのデータに対して損失関数を最小化するに過ぎません。 手持ちの計測データがもしも完全にその計測を代表できるものではないとすれば、 そのデータにフィットした今回の予測も使えるとは保証できません。

通常はデータをトレーニングデータとテストデータに分けなければならないということを覚えておきましょう。 （また、深層学習ではハイパーパラメータが大量に出てくるため、それらの調整に使うためのバリデーションデータも別途必要になります。）

更に今回大事なことがもう１つあります。 私達はデータが3次の多項式であることを知っていたため、モデルも同様に3次の多項式として表現しました。 しかし普段手に入るデータが何次の多項式で表せそうなのか、それをひと目で判断する術はありません。つまり、今回は単なる出来レースなのです。

もしも3次の多項式にノイズが加わったデータを、一目見て3次とわからなかった場合には、もっと高次の多項式でモデルを考えることになるかもしれません。 そうなると、果たして今回のようにパッと見で上手く学習できたということが言えるのでしょうか。

補足

今回のコードはパラメータもデータもすべてスカラーとして表現しました。

y=∑k=1Kakxk+b$$y=\sum _{k=1}^K{a}_k{x}^k+b$$

について、a=(b,a1,a2,a3,⋯,aK)T$\mathbf{a}=(b,a_1,a_2,a_3,\cdots ,a_K{)}^T$ と x=(1,x1,x2,x3,⋯,xK)T$\mathbf{x}=(1,x_1,x_2,x_3,\cdots ,x_K{)}^T$ と表しておけば、ベクトルの内積を使って

y=aTx$$y={\mathbf{a}}^T\mathbf{x}$$

と表せることに注意してください。これを使えば、もう少し楽にコードを書く方法があります。

コード全体

（最後のセルのコードは、学習時に予測の直線がどのようにフィットしていくかをアニメーションで表示するためのコードです）

Sorry, something went wrong. Reload?

Sorry, we cannot display this file.

Sorry, this file is invalid so it cannot be displayed.

view raw polynomial_fit_scalar.ipynb hosted with ❤ by GitHub