強化学習とニューラルネットワークを組合せてみた。（その2）

昨日は強化学習の関数近似としてニューラルネットワークを使うときの勾配計算について書いた。

ニューラルネットワークの仕様

まず、ざっとした仕様を。

構造
- 3層ニューラルネットワーク
- 入力層のユニット数、中間層のユニット数は、イニシャライザの引数で指定
- 出力層のユニットは1つ
初期化
- 各ユニットの重みは、平均0、分散 $\frac{1}{(\mathrm{ユニットへの入力のサイズ}) + 1}$ のガウス分布に従って初期化する
活性化関数
- 中間層の活性化関数 $f^{(2)}$ には、以下の正規化線形関数をちょっと変えたものを使う：
  $f^{(2)}(u) = \left\{ \begin{array}{l} u \qquad (u \ge 0) \\ 0.1u \quad (u \lt 0) \end{array} \right.$
- 出力層の活性化関数 $f^{(3)}$ には、以下のシグモイド関数を近似した関数をちょっと変えたものを使う：
  $f^{(3)}(u) = \left\{ \begin{array}{l} 0.1u + 0.9o_{\mathrm{min}} \quad (u \lt o_{\mathrm{min}}) \\ u \qquad \qquad \qquad (o_{\mathrm{min}} \le u \le o_{\mathrm{max}}) \\ 0.1u + 0.9o_{\mathrm{max}} \quad (o_{\mathrm{max}} \lt u) \end{array} \right.$
  ただし、 $o_{\mathrm{min}}$ と $o_{\mathrm{max}}$ はそれぞれ出力層で望まれる出力の最小値と最大値で、イニシャライザの引数で指定
機能
- 入力を与えると、出力と、出力の重みに関する勾配を返す
- 入力と出力の重みに関する勾配を与えると、勾配を重みに足し込んだ状態での出力を返す
- 重みを指定された重みの差分で更新する

なお、上記で正規化線形関数やシグモイド関数を近似した関数をちょっと変えているのは、微分したときの値が0になってしまうことを防ぐため。

デルタの計算を見ると分かる通り、微分したときの値が0になってしまうと、そのユニットの重みすべてが更新されないということが起こる。
それが特定の入力に対してだけ微分が0になるのなら、まだ別の入力で重みが変化する可能性があるけど、すべての入力に対して微分が0になってしまうということが仮に起きた場合、それ以上学習が出来なくなってしまう。
もちろん、そうなる可能性は低いし、そうなってしまった場合、ちょっと傾きがついてるくらいではあまり意味がないかもしれないけど、念のため。

NormalDistRandomクラス

まずはガウス分布（正規分布）に従った乱数を得るためのクラス。
といっても、強化学習について学んでみた。（その5） - いものやま。のものを再利用するだけなんだけど。

#====================
# normal_dist_random.rb
#--------------------
# 正規分布に従う乱数生成器
#====================

class NormalDistRandom
  include Math

  @@random = Random.new

  # 期待値exp, 分散varの正規分布に従った乱数を生成する
  # 乱数生成器を作成する。
  def initialize(exp=0.0, var=1.0)
    @exp = exp
    @var = var
    @values = Array.new(0)
  end

  def get_random
    if @values.size == 0
      # ボックス＝ミュラー法で乱数を生成する
      
      # NOTE
      # Random#randは[0, 1)で値を返すので、
      # (0, 1]に変換する。
      a = 1.0 - @@random.rand
      b = 1.0 - @@random.rand

      z1 = sqrt(-2.0 * log(a)) * cos(2 * PI * b)
      z2 = sqrt(-2.0 * log(a)) * sin(2 * PI * b)

      rand1 = z1 * sqrt(@var) + @exp
      rand2 = z2 * sqrt(@var) + @exp

      @values.push(rand1, rand2)
    end

    return @values.shift
  end
end

ValueNNクラス

さて、肝心のニューラルネットワークの実装。
ValueNNクラスとして実装していく。

#====================
# value_nn.rb
#--------------------
# 価値ベクトルを関数近似するためのニューラルネットワーク
#
# 中間層では、以下の活性化関数を使う：
#   f(u) = u     (u >= 0)
#          0.1*u (otherwise)
# なお、この導関数は
#   f'(u) = 1   (u >= 0)
#           0.1 (otherwise)
#
# 出力層では、以下の活性化関数を使う：
#   f(u) = u   (output_min <= u <= output_max)
#          0.1u + 0.9output_max (output_max < u)
#          0.1u + 0.9output_min (u < output_min)
# なお、この導関数は
#   f'(u) = 1   (output_min <= u <= output_max)
#           0.1 (otherwise)
#
# また、ドロップアウトも使える。（未実装）
#====================

require_relative "normal_dist_random"

class ValueNN

# 続く

ちゃっかり「ドロップアウトも使える」とか書いてるけど、未実装w

イニシャライザ

まずは初期化から。

# 続き

  def initialize(input_size, hidden_unit_size, output_min, output_max)
    @input_size = input_size
    @hidden_unit_size = hidden_unit_size
    @output_min = output_min
    @output_max = output_max

    hidden_unit_weight_variance = 1.0 / (@input_size + 1.0)
    hidden_unit_weight_generator = NormalDistRandom.new(0.0, hidden_unit_weight_variance)
    @hidden_units_weights = Array.new
    @hidden_units_bias = Array.new
    @hidden_unit_size.times do
      weights = Array.new(@input_size) do
        hidden_unit_weight_generator.get_random
      end
      bias = hidden_unit_weight_generator.get_random
      @hidden_units_weights.push(weights)
      @hidden_units_bias.push(bias)
    end

    output_unit_weight_variance = 1.0 / (@hidden_unit_size + 1.0)
    output_unit_weight_generator = NormalDistRandom.new(0.0, output_unit_weight_variance)
    @output_unit_weights = Array.new(@hidden_unit_size) do
      output_unit_weight_generator.get_random
    end
    @output_unit_bias = output_unit_weight_generator.get_random
  end

# 続く

イニシャライザは引数として、入力層のユニット数、中間層（隠れ層）のユニット数、出力層で望まれる出力の最小値と最大値をとる。
これらはインスタンス変数として保持。

そのあと、隠れ層の重みとバイアスの初期化と、出力層の重みとバイアスの初期化を行っている。

なお、重みの初期化のときに指定している分散については、ニューラルネットワークについて学んでみた。（その5） - いものやま。を参照。
ここではボードの各マスの表現が、平均0、分散1になっていることを前提としている。
けど、実際は・・・
このあたりはオンライン学習だと難しいところ。
（といっても、出力層の重みの初期化は中間層の出力の平均と分散に依存するわけで、元々の理屈がどれほど有効なのかも怪しいところではある）

出力と勾配の計算

次に、出力と勾配の計算。

# 続き

  def get_value_and_weights_gradient(input, drop_rate=0.0)
    # calculate output by forward propagation
    
    hidden_units_output = Array.new
    hidden_units_activation_gradient = Array.new
    @hidden_unit_size.times do |hidden_unit_index|
      input_sum = @hidden_units_bias[hidden_unit_index]
      @input_size.times do |input_index|
        input_sum += @hidden_units_weights[hidden_unit_index][input_index] * input[input_index]
      end
      output, activation_gradient = hidden_unit_activation_and_gradient(input_sum)
      hidden_units_output.push output
      hidden_units_activation_gradient.push activation_gradient
    end

    hidden_units_output_sum = @output_unit_bias
    @hidden_unit_size.times do |hidden_unit_index|
      hidden_units_output_sum += @output_unit_weights[hidden_unit_index] * hidden_units_output[hidden_unit_index]
    end
    output_unit_output, output_unit_activation_gradient = output_unit_activation_and_gradient(hidden_units_output_sum)

    # calculate delta by back propagation

    output_unit_delta = output_unit_activation_gradient

    hidden_units_delta = Array.new
    @hidden_unit_size.times do |hidden_unit_index|
      delta = output_unit_delta * @output_unit_weights[hidden_unit_index] * hidden_units_activation_gradient[hidden_unit_index]
      hidden_units_delta.push delta
    end

    # calculate weights gradient

    weights_gradient = Array.new
    @hidden_unit_size.times do |hidden_unit_index|
      hidden_unit_delta = hidden_units_delta[hidden_unit_index]
      weights_gradient.push hidden_unit_delta # hidden unit bias gradient
      @input_size.times do |input_index|
        hidden_unit_weight_gradient = hidden_unit_delta * input[input_index]
        weights_gradient.push hidden_unit_weight_gradient
      end
    end
    weights_gradient.push output_unit_delta # output unit bias gradient
    @hidden_unit_size.times do |hidden_unit_index|
      output_unit_weight_gradient = output_unit_delta * hidden_units_output[hidden_unit_index]
      weights_gradient.push output_unit_weight_gradient
    end

    [output_unit_output, weights_gradient]
  end

# 続く

変数名が凄いことになってるけど、適切な名前が難しくて。

やっている処理を大雑把に説明すると以下の通り：

中間層での出力 hidden_units_output と活性化関数の微分の値 hidden_units_activation_gradient の計算
出力層での出力 output_unit_output と活性化関数の微分の値 output_unit_activation_gradient の計算
出力層でのデルタ output_unit_delta の計算（これは output_unit_activation_gradient と同じ値になる）
中間層でのデルタ hidden_units_delta の計算
出力の重みに関する勾配 weights_gradient の計算とベクトル化

最後、ベクトル化と書いているのは、一言で「重み」といっても、中間層の各ユニットごとのバイアスと重み、出力層のバイアスと重みは、それぞれ別のインスタンス変数として保持されているので、それぞれの重みに対する勾配を一本のベクトルにまとめる作業が必要となるから。

感度確認

次は、感度確認の計算。
ここでいう感度確認というのは、重みを仮に変化させたときに、出力がどのように変わるのかを確認すること。

# 続き

  def get_value_with_weights_gradient(input, weights_gradient, alpha=1.0)
    weights_gradient_copy = weights_gradient.dup

    # calculate output by forward propagation, with weights_gradient
    
    hidden_units_output = Array.new
    @hidden_unit_size.times do |hidden_unit_index|
      input_sum = @hidden_units_bias[hidden_unit_index] + alpha * weights_gradient_copy.shift
      @input_size.times do |input_index|
        input_sum += (@hidden_units_weights[hidden_unit_index][input_index] + alpha * weights_gradient_copy.shift) * input[input_index]
      end
      output, _ = hidden_unit_activation_and_gradient(input_sum)
      hidden_units_output.push output
    end

    hidden_units_output_sum = @output_unit_bias + alpha * weights_gradient_copy.shift
    @hidden_unit_size.times do |hidden_unit_index|
      hidden_units_output_sum += (@output_unit_weights[hidden_unit_index] + alpha * weights_gradient_copy.shift) * hidden_units_output[hidden_unit_index]
    end
    output_unit_output, _ = output_unit_activation_and_gradient(hidden_units_output_sum)

    output_unit_output
  end

# 続く

やっていることは、基本的には順伝播による出力の計算。
ただし、引数として指定された出力に関する重みの勾配を、オプションの引数として指定されたステップサイズ（デフォルトは1.0）の分だけ仮に足し込んで、計算を行うようにしている。
これによって、指定された出力に関する重みの勾配の分だけ仮に重みを変化させたら、出力がどのように変わるのかをチェックすることが出来る。
基本的には勾配の方向に重みを変化させているので、出力される値は増えることになる。

重みの取得と更新

そして、重みの取得と更新。

# 続き

  def get_weights
    weights = Array.new
    @hidden_unit_size.times do |hidden_unit_index|
      weights.push @hidden_units_bias[hidden_unit_index]
      @input_size.times do |input_index|
        weights.push @hidden_units_weights[hidden_unit_index][input_index]
      end
    end
    weights.push @output_unit_bias
    @hidden_unit_size.times do |hidden_unit_index|
      weights.push @output_unit_weights[hidden_unit_index]
    end
    weights
  end

  def add_weights(weights_diff)
    weights_diff_copy = weights_diff.dup
    @hidden_unit_size.times do |hidden_unit_index|
      @hidden_units_bias[hidden_unit_index] += weights_diff_copy.shift
      @input_size.times do |input_index|
        @hidden_units_weights[hidden_unit_index][input_index] += weights_diff_copy.shift
      end
    end
    @output_unit_bias += weights_diff_copy.shift
    @hidden_unit_size.times do |hidden_unit_index|
      @output_unit_weights[hidden_unit_index] += weights_diff_copy.shift
    end
  end

# 続く

重みを取得した場合には、出力と勾配の計算のときと同じようにベクトル化された重みを返して、重みを更新する場合も、ベクトル化された重みの差分を使って更新を行うようにしている。
（これによって、ベクトル化された重みの内部的な順番はカプセル化されていることになる）

活性化関数と微分の値

最後に、活性化関数と、その微分の値の計算。
これらは外部からは参照されないので、privateメソッド。

# 続き

  private

  def hidden_unit_activation_and_gradient(u)
    if u >= 0.0
      [u, 1.0]
    else
      [0.1 * u, 0.1]
    end
  end

  def output_unit_activation_and_gradient(u)
    if u < @output_min
      [0.1 * u + 0.9 * @output_min, 0.1]
    elsif u <= @output_max
      [u, 1.0]
    else
      [0.1 * u + 0.9 * @output_max, 0.1]
    end
  end
end

# 続く

特に難しいことはなく。
正規化線形関数やシグモイド関数を近似した関数（をちょっと変えた関数）を使う場合、計算がシンプルなのが、実装的にも速度的にも嬉しい。

動作確認

動作確認として、次のようなコードを書いてみた。

# 続き

if __FILE__ == $PROGRAM_NAME
  require "pp"

  value_nn = ValueNN.new(3, 10, -1.0, 1.0)
  pp value_nn

  [
    [1.0, 1.0, 1.0],
    [1.0, 1.0, 0.0],
    [1.0, 0.0, 1.0],
    [0.0, 1.0, 1.0],
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
    [0.0, 0.0, 1.0],
    [0.0, 0.0, 0.0],
  ].each do |input|
    output, weights_gradient = value_nn.get_value_and_weights_gradient(input)
    output_with_weights_gradient = value_nn.get_value_with_weights_gradient(input, weights_gradient)
    diff = output_with_weights_gradient - output

    alpha = 1.0
    upper_bound = nil
    lower_bound = nil
    last = 100
    100.times do |t|
      if diff < 0.0
        upper_bound = alpha
        if lower_bound.nil?
          alpha /= 2.0
        else
          alpha = (upper_bound + lower_bound) / 2.0
        end
      elsif diff < 0.9
        lower_bound = alpha
        if upper_bound.nil?
          alpha /= diff
        else
          alpha = (upper_bound + lower_bound) / 2.0
        end
      elsif diff > 1.1
        upper_bound = alpha
        if lower_bound.nil?
          alpha /= diff
        else
          alpha = (upper_bound + lower_bound) / 2.0
        end
      else
        last = t
        break
      end

      output_with_weights_gradient_and_alpha = value_nn.get_value_with_weights_gradient(input, weights_gradient, alpha)
      diff = output_with_weights_gradient_and_alpha - output
    end

    output_with_01alpha = value_nn.get_value_with_weights_gradient(input, weights_gradient, 0.1 * alpha)
    diff_01alpha = output_with_01alpha - output

    puts "input: #{input}, output: #{output}"
    puts "  alpha: #{alpha}, iterations: #{last}"
    puts "  diff (alpha): #{diff}, diff (0.1*alpha): #{diff_01alpha}"

    weights_diff = weights_gradient.map{|i| i * 0.1 * alpha}
    value_nn.add_weights(weights_diff)
    new_output, _ = value_nn.get_value_and_weights_gradient(input)
    new_diff = new_output - output
    puts "  new output: #{new_output}, diff: #{new_diff}"
  end
end

途中やっているのは、勾配のサイズ調整のテスト。

勾配をそのまま足しこんだ場合、出力がどれくらい変わるのかは一定ではない。
そこで、バイナリサーチを行って、出力の変化が0.9〜1.1に収まる勾配のステップサイズを求めている。
こうすると、ステップサイズが短ければ変化はほぼ線形だろうという前提の元で、出力の変化が安定すると考えられる。

実際に実行してみると、次のような出力になる：

$ ruby value_nn.rb
#<ValueNN:0x007f924204c2a8
 @hidden_unit_size=10,
 @hidden_units_bias=
  [0.6388043194142576,
   0.8111472412930117,
   -0.9083398002751207,
   0.9012593799813717,
   0.4632177509972459,
   0.16069458236820633,
   -0.05268819785257466,
   0.05983555291194059,
   0.12025108968752526,
   0.08906106292892609],
 @hidden_units_weights=
  [[1.5061408191590566, 0.5559511579546162, 1.3795212066557985],
   [-1.0897041450730245, 0.9499441685536663, -0.1092271362299939],
   [0.549813751326248, -1.1526277079218483, 0.5945344045037418],
   [0.07042470560313681, 0.38558842594277487, 0.038954478972395325],
   [0.6343409035894474, 0.9496561329804604, -0.19788872777872757],
   [-0.594958117370738, -0.4022702212659905, -0.46741823267255517],
   [-0.0636861557815281, -0.16573514134489922, -0.296306156575204],
   [-0.513782423127337, -0.36264664009108954, -0.23523227542021724],
   [0.7327466434937593, 0.5200046703370264, 0.5001916419587498],
   [0.8616650920348585, 0.07843810512608271, 0.975236254617662]],
 @input_size=3,
 @output_max=1.0,
 @output_min=-1.0,
 @output_unit_bias=0.5194442005797152,
 @output_unit_weights=
  [0.11791334520427822,
   0.21087471230022242,
   0.13394362366766135,
   0.1548651002475027,
   -0.06531708340852825,
   -0.19569150418273915,
   -0.04761938363767847,
   -0.055579817882346776,
   0.3885456237749503,
   0.3543651721264999]>
input: [1.0, 1.0, 1.0], output: 1.1674510673289527
  alpha: 3.011134431135842, iterations: 1
  diff (alpha): 1.0516788262092518, diff (0.1*alpha): 0.09820411622643666
  new output: 1.2656551835553893, diff: 0.09820411622643666
input: [1.0, 1.0, 0.0], output: 1.1679823025003728
  alpha: 4.452649521234082, iterations: 4
  diff (alpha): 0.9983796457619765, diff (0.1*alpha): 0.08845108868548857
  new output: 1.2564333911858614, diff: 0.08845108868548857
input: [1.0, 0.0, 1.0], output: 1.291316835256338
  alpha: 3.6689947286853926, iterations: 4
  diff (alpha): 1.0029093301274044, diff (0.1*alpha): 0.08828419106250429
  new output: 1.3796010263188423, diff: 0.08828419106250429
input: [0.0, 1.0, 1.0], output: 1.3277405519021526
  alpha: 3.747875495975956, iterations: 3
  diff (alpha): 0.9162226758139969, diff (0.1*alpha): 0.07779803091728899
  new output: 1.4055385828194416, diff: 0.07779803091728899
input: [1.0, 0.0, 0.0], output: 1.2723117168693887
  alpha: 5.8891672682359015, iterations: 3
  diff (alpha): 0.9591138387358009, diff (0.1*alpha): 0.07684628488900391
  new output: 1.3491580017583926, diff: 0.07684628488900391
input: [0.0, 1.0, 0.0], output: 1.2847721638335425
  alpha: 5.344092313185152, iterations: 3
  diff (alpha): 0.9324888142646408, diff (0.1*alpha): 0.07711935055919072
  new output: 1.3618915143927333, diff: 0.07711935055919072
input: [0.0, 0.0, 1.0], output: 1.3815914295833625
  alpha: 5.377858892455183, iterations: 3
  diff (alpha): 0.9762747025728995, diff (0.1*alpha): 0.07665354250340606
  new output: 1.4582449720867685, diff: 0.07665354250340606
input: [0.0, 0.0, 0.0], output: 1.2030895602113556
  alpha: 11.511167259879109, iterations: 3
  diff (alpha): 1.0332763089611887, diff (0.1*alpha): 0.07704880717926832
  new output: 1.280138367390624, diff: 0.07704880717926832

各入力に対する出力と、出力の変化を0.9〜1.1に収める適切なステップサイズが求められている。
そして、そのステップサイズを1/10にすると、出力の変化もだいたい1/10になっていることが分かる。
さらに、実際に勾配を求めたステップサイズの1/10で変化させたときに新しく出力を求めると、ちゃんと期待された差分だけ出力が増えているのが分かる。
（ただ、出力が-1.0から1.0に収まっていない・・・これは入力の各成分の平均が0でなくて0.5になってしまっているので、その分だけ上方向にズレてしまっているのだと思う）

今日はここまで！