nykergoto’s blog

機械学習とpythonをメインに

重みのスケールに依存しないSGD: Path Normalized Optimization in Deep Neural Network

表題の論文を読んだのでまとめます!

Path-SGD を考えたモチベーション

ニューラルネットワークがこの論文の主題です。

Rescaling

今、あるニューラルネットワークの $i, i+1$ 番目の隠れ層の重み $W_i, W_{i+1}$ を取り出して $i$ 番目の重みを $x$ 倍して $i+1$ 番目の重みを $1/x$ 倍する操作を考えてみます(バイアスの大きさは 0 とします)。$i$ 層の出力を $y_i$ 入力を $z_i$ と表すことにします。隠れ層での計算は入力に対し重みを行列積することで表現できますから上記の操作後の隠れ層での出力を $\hat{y}_i$ とすると

$$ \begin{aligned} y_{i+1} &= W_{i+1} \cdot y_{i} = W_{i+1} \cdot W_i \cdot z_{i} \\ \hat{y}_{i+1} &= \frac{1}{x} W_{i+1} \cdot \hat{y}_i = \frac{1}{x} W_{i+1} \cdot (Wx) \cdot z_{i} \\ &= \frac{1}{x} \times x W_{i+1} \cdot W_i \cdot z_i \\ &= y_{i+1} \end{aligned} $$

となり出力は変化しません。このような重みに対する操作を rescaling と呼ぶことにします。

勾配法の Scaling への弱さ

ニューラルネットワークの学習では、勾配法の近似である確率的勾配法 (Stochastic Gradient Descent) が使われ、勾配計算には Forward Backward が用いられます。 この時重みへの正則化としてL2正則化を考えて少し重みをゼロに近づけるような効果を持たせて SGD を使うことが多いです。 この正則化は WeightDecay と呼ばれることもあります。

さてこの確率的勾配法、ひいては勾配法はネットワークの rescaling に対して不変ではありません。 これは勾配法の重み更新で使う Forward/Backward の勾配が、関連するレイヤーの重みスケールに依存しているからです。

そのため勾配法をネットワークの各層の重みが違うスケールをもっている Unbalanced な場合に使うと、学習は上手く進みません。 以下のグラフはスケールがそろっているときと揃っていないときのロス関数値をプロットしたものです。

f:id:dette:20181103194415p:plain

以下の図はSGDでの重み更新の具体です。左のネットワークではすべての重みのスケールが 1~10 程度なので調度良く値が更新されています。一方で右の例では10~100 のものと 0.1 のものが混合しています。そのため勾配更新を行うと大きい値の物はほぼ変わらず、小さい値の重みが一気に更新されてしまっています。

f:id:dette:20181103194746p:plain

「じゃあ重みが歪んでいる時でもいい感じに更新する最適化方法無いの?」という疑問に応えるのが、この論文の趣旨になっています。それが表題にもある Path-SGD です。*1

準備

考えるのは $d$ 層のニューラルネットワークです。 活性化関数は RELU とします。 このネットワークをDAG $G(V,E)$ として表現します。 ここで $V$ はノードを表し $E$ はエッジを表しています。 $D, C \in \mathbb{N}$ を入力, 出力の次元数とします。即ち入力ノードは

$$ v_{in[1]}, \ldots, v_{in[D]} \in V $$

と表現できます。

一般的な正則化

はじめに、ネットワークの各レイヤーの重みごとにグルーピングされた正則化を考えていきます。 パラメータ $p \ge 1, q \le \infty$ を用いて一般的な正則化関数は以下のように表現できます。

$$ \mu_{p,q}(w) = \left( \sum_{v \in V} \left( \sum_{(u \to v) \in E} | w_{u \to v} |^p \right)^{q/p} \right)^{1/q} $$

各ノードごとの足し算の項 $\sum_{v \in V}$ が有るためこれはグルーピングされた正則化関数の一般化であることがわかります。簡単な例であると $p = q = 1$ の時は L1 正則化に $p = q = 2$ の時にL2正則化(Weight Decay)に対応します。

これら以外にも重みの最大値を正則化とする max-normalization というものがあります。 これは上記の一般形式で $q \to \infty$ を考えた時の値として表現できます。

$$ \mu_{p,\infty}(w) = \sup_{v \in V} \left( \sum_{(u \to v) \in E} | w_{u \to v} |^p \right)^{1/p} $$

重み不変な正則化

重み不変な正則化とは、重みに対して rescaling の操作を行った時に値が変わらない正則化関数を指しています。

max-normalization は重みの上限をとっているためこれは明らかに重み不変性を満たしていません。適当なノードの入力を小さくして出力を大きくすれば、いくらでも上限値を大きくすることができるからです。

ここで、もし多数のネットワークが重み不変、すなわち入力が同じなら出力も同じである、という性質を満たしていて、どれでも好きなものを選べると仮定しましょう。 どれでも出力が同じなら、そのなかでもっとも正則化関数の値が小さいものを選ぶことが望ましいはずです。

現実には重み不変を持つネットワーク全てに対して正則化を考えることは難しいですが、ラッキーなことに max-normalization の場合この最小値を計算することが出来ます。それが以下の path-normalization です。

$$ \phi_p(w) = | \pi(w) |_p = \left( \sum_{v_{in}[i] \to \cdots v_{out}[j]} \left| \prod_{k=1}^d w_{e_k} \right|^p \right)^{1/p} $$

ここで $\pi(w) \in \mathbb{R}^N$ はパスベクトルと呼び、次元数 $N$ は入力から出力までのパスの組み合わせ数です。 上記の式中では $v_{in}[i] \to \cdots v_{out}[j]$ と表されている部分が組み合わせ数に対応しています。 この式もやや込み入っていますが、入力 $i$ と出力 $j$ を色々と入れ替えた時にできるすべてのパス、を表しているのでパスの組み合わせに相当していることが解ると思います。

パスベクトルの各次元の値はそのパス上の重みをすべて積算したものになっています。 そしてパスベクトルの $p$ ノルムが $\phi_p$ です。ちょっとややこしいですね。

ややこしいものを持ちだしたのには理由があって、実はこの $\phi_p$ は max-normalization との間に以下のような面白い性質を持ちます。

$$ \phi_p(w) = \min_{\tilde{w} \sim w} \left( \mu_{p, \infty}(\tilde{{w}}) \right)^d $$

ここで $\tilde{w} \sim w$ で表される $\tilde{w}$ は $w$ を rescaling した中で任意の入力に対して出力が $w$ と同じになるという同値類です。

これは即ち path-normalization $\phi_p(w)$ の値は rescaling な重みの中での max-normalization の最小値であることを表しています。 ということは path-normalization を最小化すれば勝手に max-normalization の最小値も下がるということがわかります。しかもその値は rescaling を許したすべての max-normalization の値の中でもっとも小さいのです。これは嬉しい。

また明らかに任意の $w \sim \tilde{w}$ に対して $\phi_p(w) \sim \phi_p(\tilde{w})$ であるので、これより path-normalization は重み不変性を持っていることがわかります。

以上から path-normalization は

  1. max-normalization の最小値であること
  2. 重み不変であること

という良い正則化項であることがわかりました。

Path-Normalization を持つ SGD

では path-normalization を正則化項とするような目的関数を考え、これに対して勾配法を考えていきます。 この時すべての重みを厳密に最適化することが難しいので、特定のエッジの重み $w_e$ のみを更新することを考えます。

今第 $t$ ステップの重み $w^{(t)}$ を得ていて $t+1$ での $e$ の重み $\hat{w}_{e}^{t+1}$ を計算する場面を考えます。 ロス関数を $L$ として偏微分すると以下の式を得ます。

$$ \hat{w}_{e}^{t+1} = w_{e}^t - \frac{\eta}{\gamma_p(w^{(t)}, e)} \frac{\partial L}{\partial w}(w^{(t)}) $$

ここで $\gamma $ は以下の式で表される重みの積になります。

$$ \gamma_p (w, e) = \left( \sum_{v_{in}[i] \cdots \overset{e}{\to} \cdots v_{out}[j]} \prod_{e_k \ne e} | w_{e_k} |^p \right)^{2/p} $$

これは入力から出力への経路のうちで $e$ の重みだけを含まない $\pi(w)$ のノルムになっています。 このことから例えば $e$ に対応する重みがとても大きい値を持っていると仮定すると $\gamma$ の値はその値が除外されるため小さくなり、勾配に係る係数は大きくなります。反対に小さい時には $\gamma$ は大きくなり勾配に係る係数は小さくなります。

結果として自分の大きさに応じた勾配更新が行われるようになることがわかります。 これは通常のSGDでは得られなかった性質です。(前半のレイヤーごとに重みが違う場合の更新を考えてみるとわかりやすいと思います)

このことを数学的に行ったのが Theorem4.1 でこのPath-SGDの更新はスケールに対して不変であることが示されています。

じっけん

有名データセットを用いて数値実験を行っています。比較するのは Path-SGD を Unbalanced な条件で学習させたものと SGD/AdaGrad を Unbalanced/Balanced な条件で学習させたものです。

f:id:dette:20181103194205p:plain

まず SGD を Unbalanced で実行すると MNIST ですら学習できていないことがわかります。 AdaGrad はもうちょっとがんばっていますがやはり苦しそうです。 一方で Path-SGD では Unbalanced にもかかわらず学習が上手く進み、また Balanced の場合の AdaGrad と並ぶかむしろそれよりも学習が早いことがわかります。*2

まとめ

  • Unbalanced な重みでは SGD は機能しない。それは勾配法が重み不変性をもたないから
  • Path-Normalization は max-normalization の最小値であり、重み不変。
  • それを正則化項として持つ objective を考えた path-SGD も重み不変。故に重みが Unbalanced でも学習が進む。

感想

重み不変性という観点が純粋に面白かった。また max-normalization の重み不変空間上での最小値が path-normalization の値につながっているのも綺麗。

一方でそもそも重みのスケールが異なるようなネットワークって存在する時あるんだろうか (転移学習とかならありうる?) とか数値実験で目的関数違うやつを同じグラフにプロットするのはどうなの(せめて精度とかにして欲しかった)とか思ったり。

気になった点としては、やはり重みごとにその勾配に対する係数を保存しないとだめだという点 ($\gamma_e$ の所)。著者は効率的な実装方法があるよーと述べていて forward backward 方式で1回計算すればいいように済ませる実装方法を提案しているがそれだとしても毎回の更新時に重みの係数も更新するのは面倒だし、重みが増えてくるとメモリ的にもしんどそう(とはいえたかだか二倍程度だからしれているのか?)。自分で実装してみたいなと思った時もこの面倒さがボトルネックでまだ出来ていません。Adam とか他の適合的なアルゴリズムが楽すぎるってことなのかな…

*1:pathという言葉がでている理由も導出を見ていると解ると思います

*2:そもそも正則化項が違うため最適解も最適値も違っているので training loss にはあまり意味は無いかも知れませんが……

Boid Model による ALife の実装 with Javascript

計算機状で生命体を模したモデルを作り、その挙動を研究するという分野があり、一般に人工生命 (Alife) と呼ばれています。 その中でも「群れ」をモデル化したものに boid model と呼ばれるものがあります。

今回は Boid Model の簡単な説明とそれを Javascipt + vue.js を使ってブラウザ上でインタラクティブにあそべるところまでを実装していきます。

完成図

agitated-goldberg-260b08.netlify.com f:id:dette:20181020175739g:plain

Boid Model とは

Boid Model とは Craig Reynolds によって開発された群れを模した Alfie のモデルのことを指します。Boid とは 「bird-oid object」の省略系で、このモデルが鳥を模したものであることが名前からもわかります。

Boid は以下の3つの基準にしたがって行動します。

1. 分離

f:id:dette:20181020164507p:plain

自分の周りに沢山鳥がいると飛びづらいですよね。なので boid は自分の近くの boid から離れようとします。これが分離の動きです。

2. 平均化

f:id:dette:20181020164540p:plain

boid は周りの動きに敏感です。自分の近くにいる集団の速度と自分の速度を一致させようとします。これが平均化の動きです。

3. 結束

f:id:dette:20181020164600p:plain

boid は周りの boid から完全に離れることは望んでいません。ですので boid はまわりの boid の重心座標へ移動しようとします。これが結束の動きです。

基礎的な boid は以上の3つのルールにしたがって行動します。

一見単純な用に見えますが、これらの重み付けを変えたり、 boid が見える周りの範囲の大きさを変えたりすることで、多種多様な動きをするようになります。

実装

今回は boid モデルを javascript で実装します。

まず boid を扱う際には物理モデルを作成する必要があります。ここでいう物理モデルは座標 $x$ 速度 $v$ 加速度 $\alpha$ を物体は持っていて以下の式によって変化する、ということを指します。

$$ \dot{x}(t) = v(t) \\ \dot{v}(t) = \alpha(t) $$

実際にコンピューター上ではこれを離散化する必要があるため、ある短い時間 $\Delta t$ を用いて以下のように近似してやります。

$$ x(t+1) = x(t) + \Delta t v(t) \\ v(t+1) = v(t) + \Delta t \alpha(t) $$

必要な数式はこれだけです。(たぶん)

下準備

では実装に入ります。まず大量に座標を扱うことになるので、座標系クラス Coodinate を作ります。とりあえずは二次元で考えるので変数は x, y の2つにしておきます。 ここでは省略しましたが座標同士の足し算、引き算、掛け算、ノルム計算などを実装しました。

class Coodinate {
  constructor(x, y) {
    this.x = x
    this.y = y
  }

  /**
   * 座標に数値を掛け算した座標を返します
   *
   * @static
   * @param {Coodinate} a
   * @param {Number} x
   * @memberof Coodinate
   */
  static mult(a, x) {
    const retval = new Coodinate(a.x * x, a.y * x)
    return retval
  }
  // などなど
}

次に物理法則にしたがって動く物体の基本となるクラス AbstractObject を作ります。このオブジェクトは場所、速度、加速度、速度の最大値、自分の位置の履歴を管理します。

class AbstractObject {
  /**
   *Creates an instance of AbstractObject.
   * @param {Coodinate} initPosition
   * @param {Numbe} maxVerocity オブジェクトの最大速度
   * @memberof AbstractObject
   */
  constructor(initPosition, maxVerocity = 5) {
    this.position = initPosition
    this.verocity = new Coodinate(Math.random(), Math.random())
    this.acceleration = new Coodinate(0, 0)
    this.maxVerocity = maxVerocity
    this.history = []
    this.ratio = null
  }

  /**
   * 加速度を元に座標と速度を更新するメソッドです
   *
   * @param {Coodinate} newAcceleration
   * @memberof AbstractObject
   */
  update(newAcceleration) {
    this.acceleration = newAcceleration
    this.verocity.plus(this.acceleration)
    const vNorm = Coodinate.distance(this.verocity, new Coodinate(0, 0))
    this.ratio = vNorm

    if (this.maxVerocity < vNorm) {
      const ratio = this.maxVerocity / vNorm
      this.verocity = this.verocity.multBy(ratio)
    }
    this.position.plus(this.verocity)

    this.history.push(Coodinate.copy(this.position))
  }
}

物体の位置が更新されるときは加速度を元に速度、速度を元に位置、というふうに計算していくので、加速度を渡して位置情報を update する関数を用意しています。

次に実際に動く boid object を作ります。今回は Fish と名づけました(今になって bird のほうが適切だったと気づきましたがまあ魚も群れを作るので)

/**
 * 魚の boid model
 *
 * @class Fish
 * @extends {AbstractObject}
 */
class Fish extends AbstractObject {
  static maxAccelerationNorm = 0.1
  static meanForceRatio = 0.7
  static dislikeForceRatio = 5.0
  static counter = 0

  /**
   * 魚のインスタンスを作成します
   * @param {Coodinate} initPosition
   * @param {number} [sakuteki=100] 魚が見える範囲です。この範囲内の魚に対して boid の3原則を適用します
   * @param {number} [dislikeDistance=20] この範囲内の魚からは遠ざかる動きをします.
   * @memberof Fish
   */
  constructor(initPosition, sakuteki = 100, dislikeDistance = 20) {
    super(initPosition)
    this.sakuteki = sakuteki
    this.dislikeDistance = dislikeDistance
    Fish.counter += 1
    this.id = Fish.counter
  }
  /**
   * 次のフレームでの自分の加速度(意志)を決定します
   *
   * @param {[Fish]} nearlyFishes
   * @memberof Fish
   */
  nextAcceleration(nearlyFishes) {
    // 近くに魚が居ない時加速しません
    if (nearlyFishes.length === 0) return new Coodinate(0, 0)

    const vList = nearlyFishes.map(f => f.verocity)
    const vMean = Coodinate.mean(vList)
    const pMean = Coodinate.mean(nearlyFishes.map(f => f.position))

    // 速度の平均値にどれぐらい合わせるかを計算 (法則2)
    const vDirection = Coodinate.minus(vMean, this.verocity)
      .norm()
      .multBy(Fish.meanForceRatio)

    // 中心に移動するちからのベクトルを計算 (法則3)
    const pDirection = Coodinate.minus(pMean, this.position).norm()

    let direction = pDirection.plus(vDirection)

    // 近すぎるおさかな
    const tooNearFishPositions = filterByDistance(
      nearlyFishes,
      this.position,
      this.dislikeDistance
    )

    if (tooNearFishPositions.length > 0) {
      const tooNearMeanPos = Coodinate.mean(
        tooNearFishPositions.map(f => f.position)
      )
      // 近すぎる魚からどのぐらい離れるかのベクトルを計算(法則1)
      const tooNearDirection = Coodinate.minus(tooNearMeanPos, this.position)
        .norm()
        .multBy(Fish.dislikeForceRatio)
      direction.minus(tooNearDirection)
    }

    return Coodinate.normalize(direction).multBy(Fish.maxAccelerationNorm)
  }
}

constructor に sakuteki と dislikeDistance のパラメータが加わりました。

  • sakuteki
    • 自分の周りのこの範囲の魚を見ることができます。先の boid 3法則はこの範囲内の魚に対して適用されます。
  • dislikeDistance
    • boid の第一法則である分離に関するパラメータです。この範囲内にいる魚からは逃げようとする動きを取ります。(コメント中の近すぎるお魚がそれに相当します。)

距離で絞り込んでくる必要がでたので関数を一つ用意しました。 viewFrom からみて maxDistance 内にいる魚を返す関数です。

function filterByDistance(fishes, viewFrom, maxDistance) {
  return fishes.filter(fish => {
    const d = Coodinate.distance(fish.position, viewFrom)
    return d < maxDistance
  })
}

最後に魚を泳がせるフィールドを用意します。

/**
 * 魚を泳がせるフィールドクラス
 *
 * @class Field
 */
export class Field {
  /**
   *Creates an instance of Field.
   * @param {number} [width=1000] フィールドの横幅
   * @param {number} [height=500] フィールドの縦幅
   * @param {number} [sakutekiRange=200] 新しく作る魚がどれぐらいの範囲を見れるか
   * @param {number} [dislikeDistance=20] 新しく作る魚はこの範囲以下の魚から離れようとします。
   * @memberof Field
   */
  constructor(
    width = 1000,
    height = 500,
    sakutekiRange = 200,
    dislikeDistance = 20
  ) {
    this.width = width
    this.height = height
    this.fishes = []
    this.newFishes = []
    this.sakutekiRange = sakutekiRange
    this.dislikeDistance = dislikeDistance
    this.isUpdating = false
  }

  /**
   *ランダムな位置に魚を追加します
   *
   * @memberof Field
   */
  addFish() {
    const x = (this.width * (0.5 + Math.random())) / 2
    const y = (this.height * (0.5 + Math.random())) / 2
    const pos = new Coodinate(x, y)
    const newFish = new Fish(pos, this.sakutekiRange, this.dislikeDistance)
    this.newFishes.push(newFish)

    if (this.isUpdating) return
    this.fishes = this.fishes.concat(this.newFishes)
    this.newFishes = []
  }

  /**
   * フィールドの時間を一つ進めます
   *
   * @memberof Field
   */
  next() {
    this.fishes = this.fishes.concat(this.newFishes.slice())
    this.newFishes = []

    // 初めに魚全員の行動(加速度)を決定する
    const accelerations = this.fishes.map(fish => {
      const fishesCanView = filterByDistance(
        this.fishes.filter(f => f !== fish),
        fish.position,
        fish.sakuteki
      )
      return fish.nextAcceleration(fishesCanView)
    })

    // 決めた行動(加速度)にしたがって更新
    this.fishes.forEach((fish, idx) => {
      const acc = accelerations[idx]

      // フィールドから外に出ようとする魚に対しては強制的に元に戻るような加速度をつける
      if (fish.position.x > this.width) {
        acc.plus(new Coodinate(-1, 0))
      }
      if (fish.position.x < 0) {
        acc.plus(new Coodinate(1, 0))
      }

      if (fish.position.y > this.height) {
        acc.plus(new Coodinate(0, -1))
      }
      if (fish.position.y < 0) {
        acc.plus(new Coodinate(0, 1))
      }

      fish.update(accelerations[idx])
    })
  }
}

動き自体は Fish が担っているので Field では特定の Fish インスタンスがどの範囲が見えているのか(見えている魚の集合はなんなのか)の計算が主になっています。

あとはこれを vue.js 上でレンダリングすると完成です。要は

  1. FIeldインスタンス field を作成
  2. window.setInterval で定期的に field.next() を呼び出す

を繰り返せば OK です。

実装: https://github.com/nyk510/boid-model/blob/master/components/boid/BoidCanvas.vue *1

Deploy

今回は netlify を使って nuxt を静的ファイルとしてデプロイしました。

agitated-goldberg-260b08.netlify.com

パラメータ変えて遊んでみると、いろいろな周期と規則が見えてきて面白いです。

今回は本当にさわりだけでしたが、面白い!と思った人は Alife とか boid で検索して本とか論文とか読んで僕に教えてくれると嬉しいです;)

*1:vueファイルの解説も書きたかったですが、あまりにも長いので断念

Pytorch で DCGAN をつくったよ

DCGAN を二年ぶりに実装しました。

github.com

いつも MNIST も面白くないので、 MNIST にアルファベットが加わった EMNIST をデータセットとして用いました。

www.nist.gov

f:id:dette:20180502060329p:plain
Emnist の画像例

こんな感じの画像が10万オーダーで格納されています。大きさは Mnist と同じ (28, 28) の grayscale の画像です。

DCGAN とは

DCGAN はニューラルネットワークの生成モデルである GAN (Generative Adversarial Networks) の一種です。 DCGANでは画像を生成するモデル (generator) だけではなく、画像が偽物(generatorが生成したもの)かデータセット上にある本物の画像かを判定するモデル (dicriminator) も同時に訓練するという枠組みのモデルです。こうすることで生成モデルはより判別モデルを騙すような高度な画像を生成し、それに応じて判別モデルはより精緻に偽画像を判定するようになり、競い合いながら学習が進んでいきます。

アルゴリズムに関しては以前の記事とかぶるので DCGANをChainerで実装 - nykergoto’s blog や、とても良くまとまっているこちら など参考にしていただけると良いかなと思います。

じっそう

pytorh での実装について。 まずはモデル部分から。 generator / dicsriminator ともに似たような感じなので generator だけ

class Generator(nn.Module):
    """
    ランダムベクトルから shape = (1, 28, 28) の画像を生成するモジュール
    Note: はじめ Linear な部分には Batchnorm1d を入れていなかったが学習が Discriminator に全く追いつかず崩壊した。
    Generator のネットワーク構成はシビアっぽい
    """

    def __init__(self, z_dim=100):
        super().__init__()
        self.z_dim = z_dim

        self.flatten_dim = 5 ** 2 * 256
        self.fc = nn.Sequential(*[
            nn.Linear(z_dim, 1024),
            nn.BatchNorm1d(1024),
            nn.ReLU(True),
            nn.Linear(1024, self.flatten_dim),
            nn.BatchNorm1d(self.flatten_dim),
            nn.ReLU(True)
        ])

        self.upsample = nn.Sequential(
            make_layer(256, 128, kernel=3, stride=1, padding=0, deconv=True), # 5 -> 7
            make_layer(128, 64, kernel=4, stride=2, deconv=True), # 7 -> 14
            nn.ConvTranspose2d(64, 1, kernel_size=4, stride=2, padding=1) # 14 -> 28
        )

    def forward(self, input_tensor):
        h = self.fc(input_tensor)
        h = h.view(-1, 256, 5, 5)
        x = self.upsample(h)
        x = F.sigmoid(x)
        return x

工夫があるとすれば make_layer という形で nn.Module を返すヘルパ関数を用意した、という点でしょうか。

def make_layer(in_channels, out_channels, kernel=2, stride=2, padding=1, use_batchnorm=True, deconv=False):
    layers = []

    if deconv:
        conv = nn.ConvTranspose2d(in_channels, out_channels, kernel_size=kernel, stride=stride, padding=padding)
    else:
        conv = nn.Conv2d(in_channels, out_channels, kernel_size=kernel, stride=stride, padding=padding)

    layers.append(conv)
    if use_batchnorm:
        bn = nn.BatchNorm2d(out_channels)
        layers.append(bn)
    activate = nn.ReLU(True)
    layers.append(activate)

    return nn.Sequential(*layers)

今回は3層程度なので直接記述しても大したことは無いですが、 resnet50 のように大きなネットワークになるとこのようにある塊をひとつの関数で作るようにすることは必要になるのかなーと思っています。(このやり方で良いのかは微妙ですが……みんなどうしているのだろうか…)

つぎに学習を実行する Solver クラスを作ってみました。

class Solver(object):
    """
    abstract class for solver
    """

    def __init__(self, models, device=None):
        """
        :param tuple[torch.nn.Module] models: 学習中に重みが更新されるモデルの配列
        :param str device: train 時のデバイス. `"cuda"` or `"cpu"`
        """
        self.models = models

        if isinstance(device, str):
            if device == "cuda" and torch.cuda.is_available():
                self.device = torch.device("cuda")
            else:
                self.device = torch.device("cpu")
        else:
            self.device = device

        self.callbacks = None

    def start_epoch(self, epoch):
        for model in self.models:
            model.to(self.device)
            model.train()

        for c in self.callbacks:
            c.initialize(self.models, self.device)
            c.start_epoch(epoch)

    def end_epoch(self, epoch, logs):
        for c in self.callbacks:
            c.end_epoch(epoch, logs)

    def fit(self, train_iterator, valid_iterator=None, epochs=10, callbacks=None, initial_lr=0.01,
            optimizer_params=None):
        """
        訓練データを用いてモデルの fitting を行う
        :param train_iterator: Iterable Object. 各 Iteration で (x_train, label_train) の tuple を返す
        :param valid_iterator:
        :param int epochs:
        :param list[Callback] callbacks:
        :param float initial_lr:
        :param dict | None optimizer_params:
        :return:
        """
        if optimizer_params is None:
            optimizer_params = {}
        self.set_optimizers(lr=initial_lr, params=optimizer_params)
        self.callbacks = callbacks

        for epoch in range(1, epochs + 1):
            self.start_epoch(epoch)
            logs = self.train(train_iterator)
            self.end_epoch(epoch, logs)
        return self

    def set_optimizers(self, lr, params):
        raise NotImplementedError()

    def train(self, train_iterator):
        logs = {}
        count = 0
        total = len(train_iterator)
        for batch_idx, (x, _) in tqdm(enumerate(train_iterator), total=total):
            count += 1
            logs = self._forward_backward_core(x, _, logs)

        for k, v in logs.items():
            if "loss" in k:
                logs[k] /= count
        return logs

    def _forward_backward_core(self, x, t, logs):
        raise NotImplementedError()

このクラスを作ったのにはちょっと理由があって、 pytorch のチュートリアルでは大概のものが epoch を引数にとって学習を行うクラス def train(epoch) のみを実行して model 部分などは global スコープにある変数を参照する形式をとっていてとても気持ちが悪かった、というのがあります。(なんとなく関数の引数以外の変数を使うのってこわくないですか)

で前使っていた keras を参考に training ループを回す部分とエポックの始まり終わりの定期実行を実装したものを書いてみた、というのが先の Solver クラスです。

最適化モデルの設定や forward-backward の定義はモデルや学習スキームによって変化するのでその部分は継承クラスで定義、という形式をとっています。

で、実装してみて思ったことなのですが……
結局こうやって枠にはめてしまうと pytorch の持っている柔軟なネットワークの構成ができる利点が薄れてしまってあまり良くないなあという印象を持ちました。

そもそも、こういう型にはめて記述できるものは keras などもっと便利に隠匿してくれるパッケージで書くことができますしそちらのほうが楽なんだろうなーということです。 pytorch を使うような場合にはそれこそ def train(epoch), def test(epoch) のように関数ベースで定義してループもごりごり書くほうが結果としてわかりやすくなるのかもしれません。 この辺はまだ pytorch 初心者なのでもっと良い書き方があるのだろうなとは思っているのでいろんなコードを読んで勉強していきたいなという所存です。

実行結果

最後になりましたが実行結果を載せます。

f:id:dette:20180502055515g:plain
実行結果 0 - 100 epoch

今回はGPUがあったので無事現実的時間内に生成できました。楽しい。

まとめ

  • pytorch 書き方がよくわからない
  • GAN は楽しい

Pytorch ことはじめ

0.40 の update もできて話題の pytorch を触り始めたのでその時のメモ。

出てくるコードはすべて以下のリポジトリ/notebooks/start で見ることができます。

github.com

参考資料

インストール

公式を見ると conda を用いてインストールするように言われるのでおとなしく従います。

conda install -y pytorch torchvision cuda91 -c python

ubuntu16.04 + pyenv + miniconda3-4.1.11 の環境ではすんなりとインストールできました。
docker つかって環境分離できたほうがよいよなあと思い公式の DockerImage をビルドしようとしたもののどうにも上手くビルドができなかった。(どうやら pytorch 公式が Docker まわり適当に放置しているみたい?) conda でインストールするバージョンの Dockerfile を作ってみた のでこちらからイメージを作成しても良いかなと思います。

Pytorch のきほん

pytorch で扱うベクトルは torch.Tensor クラスのものを用います。 このクラスには numpy っぽいメソッドがぞろぞろと定義されていて、基本的に numpy と同じ操作ができると思ってOKです。 スライス操作とかもできるので numpy に慣れている人からするととてもありがたい。

>>> x = torch.randn(2, 3)
>>> x.__class__
<class 'torch.Tensor'>
>>> x.shape
torch.Size([2, 3])

>>> x.reshape(-1)
tensor([ 1.5101, -0.4492, -0.1451,  0.8803,  0.0047, -0.0675])

>>> x[None, :]
tensor([[[ 1.5101, -0.4492, -0.1451],
         [ 0.8803,  0.0047, -0.0675]]])

# 行列積の計算は tensor.mm
>>> y = torch.randn(2, 3)
>>> y.mm(x)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
RuntimeError: size mismatch, m1: [2 x 3], m2: [2 x 3] at /opt/conda/conda-bld/pytorch_1524585239153/work/aten/src/TH/generic/THTensorMath.c:2033
>>> y.mm(x.t())
tensor([[ 0.0598,  0.1263],
        [ 0.6540, -0.0636]])

ただ負の数のステップ幅のスライスには対応していないみたい。(逆順に取ってきたい時の x[::-1] が使えない)

>>> x[::-1, :]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: negative step not yet supported

もしやりたい時には

  1. 一回 numpy array に変換してからスライス操作をかけて
  2. そのオブジェクトをコピー
  3. コピーされた numpy array を torch.Tensor に変換

というステップを取る必要がある。 ちょっとめんどくさい。 コピーをしないと

>>> torch.from_numpy(x.numpy()[::-1, :].copy())
tensor([[ 0.8803,  0.0047, -0.0675],
        [ 1.5101, -0.4492, -0.1451]])

このステップ2のコピーをするところをサボると負の値でのスライスを検知してエラーになってしまうので注意。(将来のリリースで対応するよ!っていうメッセージがでる)

>>> torch.from_numpy(x.numpy()[::-1, :])
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: some of the strides of a given numpy array are negative. This is currently not supported, but will be added in future releases.

ニューラルネットワークを作ってみる

では実際にニューラルネットワークを構築してみましょう。 想定するニューラルネットワークは 100 次元の実数値 $x \in \mathbb{R}^{100}$ を入力にとり、1次元の実数予測値 $\hat{y} \in \mathbb{R}$ を返すようなものです。

$$ \begin{align} z = W_1^T x \\ h = {\rm relu}(z) \\ \hat{y} = W_2^T h \end{align} $$

ここで $W_1 \in \mathbb{R}^{100 \times 1000}, W_2 \in \mathbb{R}^{10 \times 100}$ はそれぞれネットワークの重みを指します。

最初から pytorch の便利モジュールを使うとあれなので

  1. 手計算で更新
  2. Autograd をつかって更新
  3. torch.nn.Module & optim つかって更新

というふうにステップを経て実装を強化して行きます。

手計算で更新

行列計算ができるので、重みに対する backward をチェーンルールをつかって手計算で出して、その更新ルールにしたがって重みをアップデートする、というパワープレイから。 今回のモデルに対する勾配はロス関数を RMSE として

$$ f(y, \hat{y}) = \| y - \hat{y} \|^2 $$

とすると重みに対する勾配は以下のようになります。

$$ \begin{align} \frac{\partial f}{\partial W_2} &= - 2 (y - \hat{y}) h \\ \frac{\partial f}{\partial W_1} &= (- 2 (y - \hat{y}) W_1 {\rm diag}(H_0(h)) x )^T \\ &= -2 x^T (y - \hat{y}) W_1 {\rm diag}(H_0(h)) \end{align} $$

なのでこれにしたがって更新するように自分で書けばOK。*1

dtype = torch.float
use_gpu = False

if use_gpu:
    device = torch.device("cuda:0")
else:
    device = torch.device("cpu")
variable_types = {
    "device": device,
    "dtype": dtype
}

# ネットワーク形状の定義
batch_size = 100
input_dim = 100
hidden_dim = 1000
outpu_dim = 10

x = torch.randn(batch_size, input_dim, **variable_types)
y = torch.randn(batch_size, outpu_dim, **variable_types)

# 重みの初期値はランダムに設定
w1 = torch.randn(input_dim, hidden_dim, **variable_types)
w2 = torch.randn(hidden_dim, outpu_dim, **variable_types)

epochs = 1000
losses = []
for i in range(epochs):
    z = x.mm(w1)

    # relu activetion
    a = z.clamp(min=0)
    pred = a.mm(w2)

    loss = ((pred - y) ** 2.).sum()
    losses.append(loss.item())

    grad_pred = 2 * (pred - y)
    grad_w2 = a.t().mm(grad_pred)
    grad_a = grad_pred.mm(w2.t())
    grad_a[z < 0] = 0
    grad_w1 = x.t().mm(grad_a)

    w1 -= lr * grad_w1 
    w2 -= lr * grad_w2
    
    if i % 200 == 0:
        lr *= .5

Autograd を使う

先の例では自分で勾配を計算していました。 pytorch には計算グラフから自動的に勾配を計算する機能 (autograd) が実装されています。 この機能を使うと予測値と正解ラベルとのロスを計算する順伝搬だけ計算すれば勾配の部分 (backward) を自動的に算出してくれます。

使うためにはテンソルの定義部分で require_grad=True にすればOKです。

autograd をつかった勾配法の実装は以下のようになりました。先程のように具体的に勾配を計算する必要が無いためシンプルに記述できています。

variable_types = {
    "device": device,
    "dtype": dtype
}

variable_types["requires_grad"] = True

w1 = torch.randn(input_dim, hidden_dim, **variable_types)
w2 = torch.randn(hidden_dim, outpu_dim, **variable_types)

lr = 1e-6

losses = []
for epoch in range(epochs):
    
    pred = x.mm(w1).clamp(min=0).mm(w2)
    loss = (y - pred).pow(2).sum()
    
    # loss に対して backward を呼び出すと grad が自動計算される
    loss.backward()
    losses.append(loss.item())

    # 以下で重みを update するがこの時の変更は autograd に追跡されたくない。
    # `with torch.no_grad()` の部分では追跡を行わないようになる。
    # 他の方法としては tensor の値を直接変更するという方法がある
    # ex). w1.data -= lr* w1.grad.data
    # しかしこの方法だと履歴を追跡できなくなる
    with torch.no_grad():
        w1 -= lr * w1.grad
        w2 -= lr * w2.grad
        
        # 破壊的に勾配の値を 0 にする
        # 破壊的な変更は suffix に _ がつく
        w1.grad.zero_()
        w2.grad.zero_()

torch.nn.Module & optim をつかって更新

最後に torch.nn.Module と optim を使って重みの更新部分を隠匿してみましょう。

torch.nn.Module の使い方はこのクラスを継承して super().__init__() を呼び出すだけです。 以下は今回考えているニューラルネットワークの定義を Module を使って定義したものです。

class SimpleNNModel(torch.nn.Module):
    """
    二層のニューラルネットワークモデル"""
    
    def __init__(self, input_dim=100, output_dim=10, hidden_dim=1000):
        super().__init__()
        
        self.dense1 = torch.nn.Linear(input_dim, hidden_dim, bias=False)
        self.dense2 = torch.nn.Linear(hidden_dim, outpu_dim, bias=False)
        
    def forward(self, x):
        """
        順伝搬の計算を行って予測値を返す"""
        
        h = self.dense1(x)
        
        # equal to relu activation
        h = h.clamp(min=0)
        
        pred = self.dense2(h)
        return pred

継承したクラスでは Module に定義されている様々なメソッドを利用可能です。例えば toインスタンス内のデータすべてを cpu or gpu のメモリに展開する手助けをします。

model = SimpleNNModel()
model = model.to(device)

あとはパラメータを列挙したり, プロパティ内の Module を列挙したりもできます

for c in model.children():
    print(c)

# Linear(in_features=100, out_features=1000, bias=False)
# Linear(in_features=1000, out_features=10, bias=False)

for p in model.parameters():
    print(p)

# Parameter containing:
# tensor(1.00000e-02 *
#        [[ 7.6395, -6.7308,  9.9475,  ...,  1.5201, -3.2150, -0.1836],
#         [-8.0257,  9.9207,  2.4222,  ...,  4.2445,  8.2807,  7.5381],
#         [ 3.3367,  0.5269, -0.0468,  ...,  2.7539,  9.4210,  7.7625],
#         ...,
#         [ 6.5938, -3.1973, -0.3842,  ..., -3.0842,  6.5831,  4.6253],
#         [ 0.4635, -2.5175, -9.1168,  ...,  5.4131, -5.4361, -6.6949],
#         [ 2.5444,  9.0936, -3.0305,  ...,  6.1671,  4.1751, -3.2943]])
# Parameter containing:
# tensor(1.00000e-02 *
#        [[-3.0642,  3.1571,  1.1157,  ...,  2.1682,  2.4327, -1.7542],
#         [ 0.6330, -3.1313, -1.7368,  ...,  2.4449, -2.3944,  2.3621],
#         [-0.6492,  1.6452,  0.9604,  ..., -2.4719, -0.7905,  2.3031],
#         ...,
#         [-1.0288, -3.1262, -1.0839,  ..., -0.7796, -2.2369,  2.7909],
#         [ 0.0137, -2.7335,  2.7928,  ..., -1.3900,  0.6610,  2.2154],
#         [-3.1476, -2.1980, -1.7415,  ..., -2.2474, -2.0570, -2.9524]])

optim モジュールには optimizer が多数定義されています。 基本的に第一引数にその optimizer で更新したいパラメータの配列をわたし、その他の optimizer 自体のパラメータをそれに続いて渡す、というふうにインスタンスを生成します。

Module には parameters メソッドでパラメータ列挙ができるので、組み合わせて使うと簡単に optimizer に更新パラメータを知らせることができます。

# 勾配の初期化, 重みの update などの方法(アルゴリズム) が torch.optima にいろいろ定義されている
# 第一引数にこの optimizer で update したいパラメータの配列を渡す. 
# 今回は単純な stochastic gradient descent を使う
optimizer = torch.optim.SGD(model.parameters(), lr=1e-5, momentum=.8, weight_decay=1e-8, nesterov=True)

optim には zero_grad step メソッドが定義されています。 これらは、登録されたパラメータに対してそれぞれ勾配の初期化と更新を行ってくれるので、いちいち列挙しなくてよくてらくちんです。

ロス関数に関しても RMSE や CrossEntropy のような普通のロスであれば torch.nn 配下に定義されているのでそれを用いるとらくちん

# `criterion`: 予測値と正解ラベルとの差分を計算するオブジェクト
# `size_averaging` を `False` にすると, 与えられたバッチ n それぞれの rmse を合計した loss を返す。
# デフォルトでは `True` (すなわちバッチそれぞれの値を配列として返す)
criterion = torch.nn.MSELoss(size_average=False)

これら全部入りで学習スクリプトを書くと以下のようになります。

losses_sgd = []
for epoch in range(epochs):
    pred = model.forward(x)
    loss = criterion(pred, y)
    losses_sgd.append(loss.item())

    # 登録されたパラメータの勾配を 0 にする
    optimizer.zero_grad()
    loss.backward()

    # 登録されたパラメータを grad を使って更新する
    optimizer.step()

ちなみに autograd のときと module を使った時のロスの値をプロットすると以下のようなグラフになります。

圧倒的に Module がよいように見えますがこれは Module で定義されるパラメータの初期値が良い感じに設定されているからです。*2ニューラルネットワークの初期値はとても大事。)

f:id:dette:20180429072452p:plain

plt.figure(figsize=(8, 6))
plt.plot(losses, label="naive update")
plt.plot(losses_sgd, label="SGD with Nesterov and momentum")
plt.legend(loc=1)
plt.xlabel("epochs")
plt.ylabel("loss (log scale)")
plt.yscale("log")
plt.savefig("loss_transition.png", dpi=150)

*1:$H_0$ はヘヴィサイドの階段関数を指します

*2:重みのパラメータ数の平方根の逆数を上下限に持つ一様分布からサンプルされるようになっています。 http://pytorch.org/docs/master/_modules/torch/nn/modules/linear.html#Linear