CUBE SUGAR CONTAINER

技術系のこと書きます。

Python: Pykka でアクターモデルについて学ぶ

アクターモデルというのは、並行処理のプログラミングモデルの一つだ。 並行処理という言葉からは、まずマルチスレッドとかをイメージすると思うけど、それよりも抽象度の高い概念となっている。 つまり、アクターモデルというのはマルチスレッドなどを用いて構築することになる。 どちらかといえばプロセス間通信 (IPC) の技法であって、共有メモリやロック、RPC と比較するものかもしれない。

そんなアクターモデルは、概念とか使ったときの嬉しさを理解・実感するのがなかなか難しいモデルだとも思う。 理由としては、使い始めるまでに必要なコード量が多かったり、それなりの規模のアプリケーションで使わないとメリットが分かりづらい点が挙げられる。 ただ、これはあくまで主観的なものだけど、アクターモデルをベースに組まれたアプリケーションは規模が大きくなっても並行処理をしているコードが読みやすい。 共有メモリやロックを使う場合だと、平行処理できるパートの量と可読性がトレードオフになりがちな気がするのとは対照的な感じがしている。

上記の理由としては、それぞれの手法のアプローチが次のように根本的に異なるためだと思う。

  • 共有メモリやロックを使うやり方
    • 考え方のベースはシングルスレッド
    • ボトルネックになっている末端の処理を局所的に並行にしていく
  • アクターモデル
    • 考え方のベースはマルチスレッド
    • 末端の処理はシングルスレッドで処理される
      • ただし、アクターを横に並べることで並行度を上げられる

今回は、そんなアクターモデルを Python で実装するためのフレームワークである Pykka を使ってみることにする。 ざっと調べた感じ、この Pykka が Python で汎用的なアクターモデルを実装するためのパッケージとしては最も有名そうだった。 Pykka は Scala/Java でアクターモデルを実装するためのフレームワークである Akka の影響を受けている。 使われている用語の多くも Akka から取られているらしいけど、まあ使う分には特に意識する必要はない。

ただし、Pykka は単純に Akka を Python で再実装したものではない。 次のコンポーネントは明示的にサポートしないことが宣言されている。 もし欲しければ自分で作り込まなきゃいけない。

  • スーパーバイザ
    • セルフヒーリングとかを実現するやつ
  • メッセージルータ
    • 並行処理を書きやすくしたりアプリケーションの骨子を作るためのやつ
  • ネットワーク越しの分散処理
    • 並行処理のタスクをネットワーク越しの各ホストに分配するやつ

今回使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.12.3
BuildVersion:   16D32
$ python --version
Python 3.6.0

インストール

インストールは Python のパッケージマネージャである pip を使ってさくっとできる。

$ pip install pykka

もし pip が入っていないときはインストールしよう。

$ curl https://bootstrap.pypa.io/get-pip.py | sudo python

あと、システムの実行環境を汚さないように Python 仮想環境を使えるようにしておくのがおすすめ。 これには例えば virtualenv なんかを使う。

$ sudo pip install virtualenv

そこらへんはお好みで。

はじめてのアクターモデル

アクターモデルでは、その名の通りアクターというコンポーネントを組み合わせてアプリケーションを構築していく。

アクターという概念を構成する要素は、実のところたったの三つしかない。 まず、他のアクターから何らかのデータを受け取るために用意されたメッセージキューが一つ目。 そして、その受け取ったデータを渡す先となるメッセージハンドラが二つ目。 最後に、メッセージキューからデータを取り出して、それをメッセージハンドラにディスパッチする役目のスレッドが三つ目だ。 この三つの要素が、それぞれのアクターインスタンスに必ず備わっている。 メッセージキューとスレッドはデフォルトで一つずつ、メッセージハンドラについては自分で何をするか記述する。

重要な点としては、アクター同士はメモリなどを共有することがない。 することがないというよりも、する必要がない・するとアクターモデルを使う意味が薄れるという方が正しい。 アクターモデルでは、何かを共有することなくお互いのメッセージキューにメッセージを放り込み合うことで情報をやり取りして処理を進めていく。 いわゆるメッセージパッシングと呼ばれるものだ。

また、アクターの中と外の世界は非同期の壁で断絶されている。 アクターの中身を、その外側から同期的に直接触ることはできない。 触ることができてしまうと、処理の競合が起こってしまうからだ。 アクターの中を触るには、いわゆる Future を必ず介すことになる。 処理はキューイングされていつかは実行されるけど、具体的にそれがいつかは分からない状態でアクセスする。

説明が長くなったけど、そろそろ具体的なサンプルコードを見ていくことにしよう。 以下は Pykka を使って作る、最も単純なアプリケーションの一つだ。 このアプリケーションではメッセージを受け取って、それを出力するだけのアクターを定義している。 説明はコードの中にコメントとして記述している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class GreetingActor(pykka.ThreadingActor):
    """メッセージを受け取って挨拶をするアクター"""

    def on_receive(self, message):
        """メッセージが届いたときに実行されるハンドラ"""
        # メッセージの中から名前を取り出す
        name = message['name']
        # 内容をプリントする
        print('Hello, {0}!'.format(name))


def main():
    # アクターを起動する
    actor_ref = GreetingActor.start()

    # アクターにメッセージを送る (一方向)
    actor_ref.tell({'name': 'World'})

    # アクターを停止する (本当はちょっと待った方が良いかも)
    actor_ref.stop()


if __name__ == '__main__':
    main()

アクターは pykka.ThreadingActor というクラスを継承して作る。 on_receive() というメソッドはアクターがメッセージを受け取って、それがディスパッチされたときに発火するハンドラだ。 サンプルコードでは、そのアクターを起動してメッセージを送りつけている。 アクターは start() メソッドが呼ばれた瞬間からメッセージキューを読み出すスレッドが起動して動き出す。 あとは、アクターに tell() メソッドでメッセージを送ることができる。

また、ポイントとしてはアクターを起動したときに返ってくる内容が ActorRef というアクターへの間接的な参照になっている点だ。 つまり、起動したアクターそのものを直接触ることはできず、必ずこの参照を経由してアクセスすることになる。 直接触ることを防ぐことによって、アクターの中で競合を起こさないようにしているのだろう。

ちなみに、上記のサンプルコードではメッセージを送りつけた直後にアクターを停止しているけど、本当はちょっと待った方がお行儀が良いのかもしれない。 ちょっと待つというのは、具体的にはメインスレッドをスリープさせるということだ。 何故なら、キューからのメッセージの取り出しおよびハンドラの発火は別スレッドで進むため、すぐに実行される保証がないから。 とはいえ、他にやっていることもないからすぐに処理されるはずで、現実的には問題ないんだけど。

上記のサンプルコードを実行すると、次のようになる。

$ python hellopykka.py
Hello, World!

ただのハローワールドだけど、これはアクターモデルで実装されている!

アクターモデルのコンストラクタに引数を渡す

ここからは、アクターモデルというより Pykka の使い方の説明に入っていく。

先ほどのサンプルコードでは、アクターにメッセージの形で引数を渡した。 とはいえ、起動時の初期パラメータを設定したい場合も多いと思う。 そんなときはコンストラクタに引数を用意すれば良い。 親クラスのコンストラクタは必ず呼ぶこと。 というより、呼ばないと例外になる。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class GreetingActor(pykka.ThreadingActor):

    def __init__(self, greeting='Hello'):
        # 親クラスのコンストラクタは必ず呼ぶ
        super().__init__()

        # 挨拶の内容をカスタマイズできるようにする
        self.greeting = greeting

    def on_receive(self, message):
        name = message['name']
        print('{0}, {1}!'.format(self.greeting, name))


def main():
    # アクターを起動するときに使ったパラメータがコンストラクタに渡される
    actor_ref = GreetingActor.start('Hajimemashite')

    actor_ref.tell({'name': 'Sekai'})

    actor_ref.stop()


if __name__ == '__main__':
    main()

コンストラクタに渡す引数はアクターの起動時に渡す。

上記のサンプルコードを実行した結果は次の通り。

$ python initpykka.py
Hajimemashite, Sekai!

コンストラクタで渡したメッセージが使われている。

返り値のあるメッセージパッシング

これまでの例では、アクターにメッセージを送っておしまいという、一方的なものだった。 アクターに対して一方的にメッセージを送りつけるには tell() メソッドを使う。 とはいえ、何かメッセージを渡して、その結果として返り値を戻してほしいというのはよくある。

そんなときはメッセージを送りつけるときに ask() メソッドを使う。 また、メッセージハンドラでは返り値を戻すように作る。 返り値を受け取る側では、いくつか受け取り方を選ぶことができる。 次のサンプルコードでは、戻り値のあるパターンを実装している。 そして、返り値の受け取り側でもいくつかのやり方を試している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class GreetingActor(pykka.ThreadingActor):

    def on_receive(self, message):
        name = message['name']
        # 返り値があると ask() メソッドでその値が得られるようになる
        return 'Hello, {}!'.format(name)


def main():
    actor_ref = GreetingActor.start()

    # ActorRef#ask() を使うとアクターから返り値が得られる (返事があるまでブロックする)
    answer = actor_ref.ask({'name': 'World'})
    print(answer)

    # ブロックする時間にタイムアウトを設定したいときは timeout を指定する
    answer = actor_ref.ask({'name': 'Timeout'}, timeout=1)
    print(answer)

    # 処理をブロックさせたくないときは block=False を指定すると
    # 返り値が Future オブジェクトになる
    future = actor_ref.ask({'name': 'Future'}, block=False)
    # Future オブジェクトから結果を取り出す (ここでブロックする)
    answer = future.get()
    print(answer)

    future = actor_ref.ask({'name': 'Future-Timeout'}, block=False)
    # Future から結果を取り出すときのタイムアウトも指定できる
    answer = future.get(timeout=1)
    print(answer)

    actor_ref.stop()


if __name__ == '__main__':
    main()

上記で返り値を受け取るやり方がいくつかあることについて補足しておく。 前述した通り、アクターの内外は非同期の壁で断絶されている。 つまり、メッセージに対する返事が返ってくるタイミングは、将来であることは分かるものの具体的にいつかは分からないのだ。

そのため、返事の受け取り方として「返ってくるまですぐにブロックする」か「とりあえず受理された旨の返事だけ受け取る」やり方を選べるわけだ。 後者では、いわゆる Future を受け取っている。

先ほどのサンプルコードを実行した結果は次の通り。

$ python askpykka.py
Hello, World!
Hello, Timeout!
Hello, Future!
Hello, Future-Timeout!

アクターの中で発生した例外を処理する

もし、アクターの返り値のある処理で例外が発生したことに気づくことができないと大変だ。 Pykka では ActorRef#ask() メソッド経由で呼び出された処理の中で例外が起こったときは、呼び出し元にそれが伝搬するようになっている。

次のサンプルコードではメッセージハンドラの中で意図的に例外を上げている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class ExceptionActor(pykka.ThreadingActor):

    def on_receive(self, message):
        raise Exception('Oops!')


def main():
    actor_ref = ExceptionActor.start()

    try:
        # ちなみに ActorRef.tell() を使うときは、投げっぱなしなので
        # 呼び出し先で例外が起こっても呼び出し元でキャッチできない
        actor_ref.ask({'name': 'World'})
    except Exception as e:
        print(e)

    actor_ref.stop()


if __name__ == '__main__':
    main()

上記のサンプルコードを実行した結果は次の通り。

$ python trypykka.py
Oops!

親スレッドが死んだときにアクターのスレッドも停止させたい

ちなみに、先ほどのサンプルコードで例外をキャッチしないでいるとどうなるか試してみただろうか。 実はアプリケーションの実行はそのまま継続することになる。 どうしてそんなことが起こるかというと、子スレッドがデーモンスレッドになっていないためだ。 アプリケーションは全てのスレッドが停止するまで処理を継続することになる。

もし、親スレッドが停止したときに子スレッドも停止させたい、つまり一家心中を図りたいときは次のようにする。 アクターには use_daemon_thread というアトリビュートが用意されているので、そのフラグを True にしよう。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class ExceptionActor(pykka.ThreadingActor):
    # 親のスレッドが終了したら、同時に停止させたいときはこうする
    # 親が死んだとき一家心中になるスレッドをデーモンスレッドという
    use_daemon_thread = True

    def on_receive(self, message):
        raise Exception('Oops!')


def main():
    actor_ref = ExceptionActor.start()

    # デーモンスレッドになっていないと、ここで例外になって親スレッドが死ぬ
    actor_ref.ask({'name': 'World'})

    # すると、このコードが実行されないのでアクターのスレッドが停止しない
    actor_ref.stop()


if __name__ == '__main__':
    main()

こうすれば、例外を取りこぼしたことでメインスレッドが停止したときに、そこから起動されたアクターのスレッドも同時に停止する。

$ python daemonpykka.py
Traceback (most recent call last):
  File "daemonpykka.py", line 27, in <module>
    main()
  File "daemonpykka.py", line 20, in main
    actor_ref.ask({'name': 'World'})
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/actor.py", line 435, in ask
    return future.get(timeout=timeout)
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/threading.py", line 52, in get
    compat.reraise(*self._data['exc_info'])
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/compat.py", line 24, in reraise
    raise value
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/actor.py", line 201, in _actor_loop
    response = self._handle_receive(message)
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/actor.py", line 304, in _handle_receive
    return self.on_receive(message)
  File "daemonpykka.py", line 13, in on_receive
    raise Exception('Oops!')
Exception: Oops!

アクターのフラグを切り替えて、挙動の違いを確認してみよう。

アクタープロキシを使って処理を呼び出す

これまでの例では pykka.ThreadingActor の持つ on_receive() というメソッドをオーバーライドすることでアクターを実装してきた。 しかし、このやり方はだいぶプリミティブで Pykka の低レイヤーな API を使っている。 実際には、アクタープロキシという機構を使って、アクターのメソッドをあたかも直接呼び出しているような風に処理できる。

次のサンプルコードではアクターで独自のメソッド greeting を定義している。 その呼び出し側では ActorRef から ActorProxy を生成して、あたかもメソッドを直接呼び出しているように扱っている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class GreetingActor(pykka.ThreadingActor):

    def greeting(self):
        print('Hello, World!')


def main():
    actor_ref = GreetingActor.start()

    # ActorProxy を取得する
    actor_proxy = actor_ref.proxy()

    # ActorProxy では定義されているメソッドを直接呼べる (ように見える)
    actor_proxy.greeting()

    actor_ref.stop()


if __name__ == '__main__':
    main()

上記のサンプルコードを実行した結果は次の通り。

$ python proxypykka.py
Hello, World!

アクタープロキシを使って返り値のある場合

先ほどのサンプルコードでは、アクタープロキシを使って返り値のない場合を扱っていた。 今度は返り値のあるときを試してみる。 アクタープロキシを使って返り値を受け取るときは、実は必ず Future を介すことになる。

次のサンプルコードではアクタープロキシを使った上で返り値を受け取っている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class GreetingActor(pykka.ThreadingActor):

    def greeting(self):
        return 'Hello, World!'


def main():
    actor_ref = GreetingActor.start()
    actor_proxy = actor_ref.proxy()

    # 返り値のある呼び出しでは Future オブジェクトが返る
    future = actor_proxy.greeting()
    answer = future.get()
    print(answer)

    actor_ref.stop()


if __name__ == '__main__':
    main()

アクタープロキシ経由で実行したメソッドは、必ず Future 経由で返り値を受け取ることになる。

上記の実行結果は次の通り。

$ python futurepykka.py
Hello, World!

アクターのアトリビュートにアクタープロキシ経由でアクセスする

前述した通り、アクターの中と外は非同期の壁で断絶されていることから、直接同期的に触ることはできない。 それをすると、アクターの処理が競合する恐れがあるためだ。

では、アクターのアトリビュートを触るのにも、必ずアクセサに相当するメソッドを定義する必要があるのだろうか? だとすると相当に面倒くさいと思うんだけど、これは定義する必要がない。 アクタープロキシ経由であればアトリビュートも Future 経由にはなるが、触ることができる。

次のサンプルコードではアクターに message というアトリビュートを定義している。 メインスレッドからは、それをアクタープロキシ経由で読み出している。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class MyActor(pykka.ThreadingActor):

    def __init__(self):
        super().__init__()
        # アクターでないオブジェクトをアトリビュートとして持つ
        self.message = 'Hello, World!'


def main():
    actor_ref = MyActor.start()
    actor_proxy = actor_ref.proxy()

    # アクターでないアトリビュートにアクセスすると Future が返る
    future = actor_proxy.message
    # 値を取り出す
    answer = future.get()
    print(answer)

    actor_ref.stop()


if __name__ == '__main__':
    main()

ポイントとしては、アトリビュートへのアクセスであっても Future が返る点だ。 こうすればアクセスが非同期に処理されるため、処理が競合することはない。

上記の処理結果は次の通り。

$ python attrpykka.py
Hello, World!

アクターのメンバが持つメソッドをアクタープロキシ経由で呼び出す

先ほどのサンプルコードではアクターのメンバがただの文字列だった。 次は、これがメソッドを持ったインスタンスだったときの話。

次のサンプルコードでは、アクター MyActor がアクターでないアトリビュート NonActor のインスタンスをメンバに持っている。 それを、アクタープロキシ経由で外から呼び出している。 この場合、メンバとなるアクターでないオブジェクトには pykka_traversable というフラグが必要になる。 このフラグを立てることで、アクタープロキシから呼び出すことができるようになる。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class NonActor(object):
    """アクターのアトリビュートになる非アクタークラス"""
    # pykka_traversable アトリビュートを True にするとプロキシ経由でアクセスできるようになる
    pykka_traversable = True

    def greeting(self):
        """プロキシ経由で外部から呼びたいメソッド"""
        return 'Hello, World!'


class MyActor(pykka.ThreadingActor):

    def __init__(self):
        super().__init__()
        # アクターでないオブジェクトをアトリビュートとして持つ
        self.non_actor = NonActor()


def main():
    actor_ref = MyActor.start()
    actor_proxy = actor_ref.proxy()

    # アクターのアトリビュートのオブジェクトが持つメソッドを呼び出せるようになる
    future = actor_proxy.non_actor.greeting()
    # Future が返ってくるので値を取り出す
    answer = future.get()
    print(answer)

    actor_ref.stop()


if __name__ == '__main__':
    main()

このパターンでも返り値は Future になる。

上記のサンプルコードの実行結果は次の通り。

$ python methodpykka.py
Hello, World!

アクターモデルを使った並行処理

これまでのサンプルコードでは、アクターは定義しているものの平行処理をしていなかった。 なので、説明の都合上、仕方なかったとはいえ一体何のためにこんなことを、という印象だったかもしれない。 ここからは、ようやくアクターモデルを使って並行処理をしてみる。

アクターモデルの並行処理では、基本的にアクターをたくさん用意してタスクをそれぞれのアクターに振り分けていく。 アクターには、それぞれメッセージを処理するスレッドが起動しているので、処理をバラバラに並列で進めることができるわけ。

次のサンプルコードでは、 アクター MyActor に擬似的に処理に時間のかかるメソッド process() を用意している。 時間がかかる、というところはスレッドをスリープすることで表現した。 ここでは、時間のかかる処理を 4 回呼び出している。 それぞれの処理には 1 秒かかる。 一つのアクターで直列に処理すれば、単純に 4 秒かかるはずの処理を、ここでは二つのアクターに振り分けている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division

import time

import pykka


class MyActor(pykka.ThreadingActor):

    def process(self, processing_id, sleep_sec):
        """時間がかかる処理を模したメソッド"""
        # time.sleep() を使って時間がかかる処理をエミュレートする
        time.sleep(sleep_sec)
        # 処理が終わったら内容を出力する
        print('Completed: ID {0} in {1} s'.format(processing_id, sleep_sec))


def main():
    # アクターを二つ起動する
    actor_proxy_a = MyActor.start().proxy()
    actor_proxy_b = MyActor.start().proxy()

    # 処理を開始した時刻を記録しておく
    start_time = time.time()

    # それぞれのアクターに仕事を二つずつ割り当てる
    # 処理はアクターのスレッドがキューから取り出して開始する
    future1 = actor_proxy_a.process(1, 1)
    future2 = actor_proxy_a.process(2, 1)

    future3 = actor_proxy_b.process(3, 1)
    future4 = actor_proxy_b.process(4, 1)

    # 計算が終わるのを待つ (threading#join() のようなものと考えると分かりやすい)
    future1.get()
    future2.get()
    future3.get()
    future4.get()

    # 処理が終わったら時刻を記録する
    end_time = time.time()

    # 一連の処理にかかった時間を出力する
    elapsed_time = end_time - start_time
    print('Elapsed Time: {0} s'.format(elapsed_time))

    # すべてのアクターを停止させる
    pykka.ActorRegistry.stop_all()


if __name__ == '__main__':
    main()

さて、処理が終わるのに実際は何秒かかるだろうか?

上記のサンプルコードの実行結果は次の通り。 二つのアクターで処理した結果として本来一つのアクターなら 4 秒かかるところ 2 秒で完了している。

$ python cuncurrentpykka.py
Completed: ID 3 in 1 s
Completed: ID 1 in 1 s
Completed: ID 2 in 1 s
Completed: ID 4 in 1 s
Elapsed Time: 2.0071158409118652 s

アクターモデルでも設計によっては競合が起こる

アクターモデルを使えば、競合が起きないというような説明をされることがたまにあるようだ。 しかし、その「競合が起きない」というのは、より低レイヤーな側面からアクターモデルを見たときの話だと思う。 あくまで、アクターのメッセージキューを処理するスレッドが一つだけなので、アクターの中で処理が競合しないというのに過ぎない。 より高レイヤーに、アプリケーションという側面からアクターモデルを見たとき、競合は容易に起こる。 それぞれの競合は全く別の次元の話なんだけど、もし混ぜて考えてしまうと危険だ。 次は、その競合について見てみよう。

サンプルコードでは、銀行口座を模したモデルを用意した。 AccountActor という銀行口座を模したアクターを UserActor というユーザを模したアクターが操作する。 UserActor はまず銀行口座から預金残高を読み取って、その値を元に預金残高を更新する。 ここでポイントなのは AccountActor の API がアトミックになっていないという点だ。 つまり、銀行口座の読み取りから更新までの間に別のアクターからの処理が割り込む余地がある。 サンプルコードでは、口座の読み取りから更新までの間に意図的にスリープを入れることで競合が起こりやすいようにしている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import time

import pykka


class AccountActor(pykka.ThreadingActor):
    """銀行口座を模したアクター"""

    def __init__(self, initial_money=0):
        super(AccountActor, self).__init__()

        self._money = initial_money

    def get_balance(self):
        # 口座の残高を得る
        return self._money

    def set_balance(self, money):
        # 口座の残高を設定する
        self._money = money


class UserActor(pykka.ThreadingActor):
    """銀行口座を使うユーザを模したアクター"""

    def __init__(self, username, account_actor):
        super(UserActor, self).__init__()

        self._username = username
        self._account_actor = account_actor

    def save(self, amount):
        # 残高を読み出す
        future = self._account_actor.get_balance()
        balance = future.get()
        print('{0} read balance: {1}'.format(self._username, balance))
        # 残高を読んでから次の処理に移るのに時間がかかった!
        time.sleep(1)
        print('{0} sleep...'.format(self._username))
        # 口座残高を更新する
        new_balance = balance + amount
        future = self._account_actor.set_balance(new_balance)
        future.get()
        print('{0} set balance: {1}'.format(self._username, new_balance))


def main():
    # アクターを起動する
    account_actor_ref = AccountActor.start(100)
    account_actor_proxy = account_actor_ref.proxy()

    user_a_ref = UserActor.start('UserA', account_actor_proxy)
    user_a_proxy = user_a_ref.proxy()

    user_b_ref = UserActor.start('UserB', account_actor_proxy)
    user_b_proxy = user_b_ref.proxy()

    # ユーザ A と B が同時に口座に 50 を入金する
    future_a = user_a_proxy.save(50)
    future_b = user_b_proxy.save(50)

    future_a.get()
    future_b.get()

    # 最終的な口座の残高を出力する
    print('Final Balance: {0}'.format(account_actor_proxy.get_balance().get()))

    pykka.ActorRegistry.stop_all()


if __name__ == '__main__':
    main()

二つのユーザを模したアクターは銀行残高をそれぞれ 50 だけ増やすべく処理を実施する。 初期の預金残高は 100 なので、最終的には 200 になっていないとおかしい。

では、上記のサンプルコードを実行してみよう。

$ python conflictpykka.py
UserA read balance: 100
UserB read balance: 100
UserA sleep...
UserB sleep...
UserA set balance: 150
UserB set balance: 150
Final Balance: 150

最終的な口座残高は 150 となってしまった。 見事に競合が発生している。

問題は、アクターの API が更新に対してアトミックになっていなかった点だ。 次の改良したサンプルコードでは、預金残高を更新するメソッドを update_balance() という形でアトミックにしている。 現在の預金残高を元に、口座のアクター自体が更新をかけるので、間に別のアクターの処理が入り込む余地はない。 アクターの中はメッセージキューを処理する一つのスレッドで処理されているので、同時に同じメソッドが走る心配はない。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import pykka


class AccountActor(pykka.ThreadingActor):

    def __init__(self, initial_money=0):
        super(AccountActor, self).__init__()

        self._money = initial_money

    def get_balance(self):
        return self._money

    def update_balance(self, money):
        # アトミックに口座残高を更新する
        self._money = self._money + money


class UserActor(pykka.ThreadingActor):

    def __init__(self, username, account_actor):
        super(UserActor, self).__init__()

        self._username = username
        self._account_actor = account_actor

    def save(self, amount):
        # 口座残高を更新する
        future = self._account_actor.update_balance(amount)
        future.get()
        print('{0} update balance: {1}'.format(self._username, amount))


def main():
    # アクターを起動する
    account_actor_ref = AccountActor.start(100)
    account_actor_proxy = account_actor_ref.proxy()

    user_a_ref = UserActor.start('UserA', account_actor_proxy)
    user_a_proxy = user_a_ref.proxy()

    user_b_ref = UserActor.start('UserB', account_actor_proxy)
    user_b_proxy = user_b_ref.proxy()

    # ユーザ A と B が同時に口座に 50 を入金する
    future_a = user_a_proxy.save(50)
    future_b = user_b_proxy.save(50)

    future_a.get()
    future_b.get()

    # 最終的な口座の残高を出力する
    print('Final Balance: {0}'.format(account_actor_proxy.get_balance().get()))

    pykka.ActorRegistry.stop_all()


if __name__ == '__main__':
    main()

上記のサンプルコードを実行した結果は次の通り。 今度はちゃんと預金残高が 200 に更新されている。

$ python conflictpykka2.py
UserA update balance: 50
UserB update balance: 50
Final Balance: 200

アクターモデルでも設計によってはデッドロックが起こる

典型的なデッドロックは、マルチスレッドでリソースをロックする順番が不定になっていたときに起こる。 アクターモデルはロックを使わないので、それならデッドロックが起こらないかというと、そんなことはない。

次のサンプルコードではデッドロックを発生させている。 ActorAActorB のメソッドを呼んで、ActorBActorA のメソッドを呼ぶ、というように操作が交差している。 それぞれのメソッドは各アクターのスレッドで非同期に処理されるので、返り値で得られる Future を get() すると完了するまでブロックする。 すなわち、双方がお互いの処理の完了を待ってデッドロックを起こしてしまう。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import logging
import signal
import os

import pykka
from pykka import debug


class ActorA(pykka.ThreadingActor):

    def method_a(self, actor_b_proxy):
        # 2. ActorB の method_b を呼ぶ
        future = actor_b_proxy.method_b()
        # 3. Future の中身を取り出そうとすると ActorA のスレッドは取り出せるまでブロックする
        future.get()


class ActorB(pykka.ThreadingActor):

    def __init__(self, actor_a_proxy):
        super(ActorB, self).__init__()
        self.actor_a_proxy = actor_a_proxy

    def method_b(self):
        # 4. ActorA の method_a を呼ぶ (本当は自身のプロキシを渡す必要があるけど)
        future = self.actor_a_proxy.method_a()
        # 5. Future の中身を取り出そうとすると ActorB のスレッドは取り出せるまでブロックする
        future.get()
        # XXX: ActorA のスレッドは既にブロックしているのでデッドロックを起こす


def main():
    """意図的にデッドロックを起こしてみる"""

    # デバッグレベルのログが出力されるようにする
    logging.basicConfig(level=logging.DEBUG)

    # シグナル SIGUSR1 でスレッドのトレースバックを出力する
    signal.signal(signal.SIGUSR1, debug.log_thread_tracebacks)

    # アクターを起動する
    actor_a_ref = ActorA.start()
    actor_a_proxy = actor_a_ref.proxy()
    # ActorB には ActorA のプロキシを渡しておく
    actor_b_ref = ActorB.start(actor_a_proxy)
    actor_b_proxy = actor_b_ref.proxy()

    # 1. ActorA の method_a を呼ぶ (引数として ActorB を渡す)
    actor_a_proxy.method_a(actor_b_proxy)

    # メインスレッドでは PID を出力する
    pid = os.getpid()
    print(f'PID: {pid}')

    # アクターがデーモンスレッドではないので、ここまできてもプログラムは終了しない


if __name__ == '__main__':
    main()

上記を実行すると、何も表示されずに処理が止まる。

$ python deadlockpykka.py
DEBUG:pykka:Registered ActorA (urn:uuid:5b470837-65b9-408d-9843-23982f2efcc5)
DEBUG:pykka:Starting ActorA (urn:uuid:5b470837-65b9-408d-9843-23982f2efcc5)
DEBUG:pykka:Registered ActorB (urn:uuid:02ab87fd-705b-4ec4-9718-78870b99eb90)
DEBUG:pykka:Starting ActorB (urn:uuid:02ab87fd-705b-4ec4-9718-78870b99eb90)
PID: 9199

サンプルコードでプロセス ID を表示するようにしているところに気づいただろうか? また Pykka の debug.log_thread_tracebacks という関数がシグナル SIGUSR1 に対して登録されている。 そこで、別のターミナルからこのハンドラに対してシグナルを送ってみよう。

$ kill -SIGUSR1 9199

すると、今のスレッドが何処を処理しているのかトレースバックが表示される。 この表示からは、二つのアクターのスレッドが future.get() でブロックしていることが分かる。

$ python deadlockpykka.py
DEBUG:pykka:Registered ActorA (urn:uuid:5b470837-65b9-408d-9843-23982f2efcc5)
DEBUG:pykka:Starting ActorA (urn:uuid:5b470837-65b9-408d-9843-23982f2efcc5)
DEBUG:pykka:Registered ActorB (urn:uuid:02ab87fd-705b-4ec4-9718-78870b99eb90)
DEBUG:pykka:Starting ActorB (urn:uuid:02ab87fd-705b-4ec4-9718-78870b99eb90)
PID: 9199
CRITICAL:pykka:Current state of ActorB-2 (ident: 123145451573248):
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 884, in _bootstrap
    self._bootstrap_inner()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 916, in _bootstrap_inner
    self.run()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 864, in run
    self._target(*self._args, **self._kwargs)
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/actor.py", line 201, in _actor_loop
    response = self._handle_receive(message)
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/actor.py", line 295, in _handle_receive
    return callee(*message['args'], **message['kwargs'])
  File "deadlockpykka.py", line 31, in method_b
    future.get()
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/threading.py", line 50, in get
    self._data = self._queue.get(True, timeout)
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/queue.py", line 164, in get
    self.not_empty.wait()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 295, in wait
    waiter.acquire()

CRITICAL:pykka:Current state of ActorA-1 (ident: 123145446318080):
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 884, in _bootstrap
    self._bootstrap_inner()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 916, in _bootstrap_inner
    self.run()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 864, in run
    self._target(*self._args, **self._kwargs)
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/actor.py", line 201, in _actor_loop
    response = self._handle_receive(message)
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/actor.py", line 295, in _handle_receive
    return callee(*message['args'], **message['kwargs'])
  File "deadlockpykka.py", line 18, in method_a
    future.get()
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/threading.py", line 50, in get
    self._data = self._queue.get(True, timeout)
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/queue.py", line 164, in get
    self.not_empty.wait()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 295, in wait
    waiter.acquire()

CRITICAL:pykka:Current state of MainThread (ident: 140736596534208):
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 1290, in _shutdown
    t.join()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 1056, in join
    self._wait_for_tstate_lock()
  File "/Users/amedama/.pyenv/versions/3.6.0/lib/python3.6/threading.py", line 1072, in _wait_for_tstate_lock
    elif lock.acquire(block, timeout):
  File "/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages/pykka/debug.py", line 68, in log_thread_tracebacks
    stack = ''.join(traceback.format_stack(frame))

Pykka の並行処理に関しては、上記のようにしてデバッグができる。

スレッドプールっぽいものを作って並行処理をしてみる

次はもうちょっと実践的に、スレッドプールっぽく並行処理をしてみよう。 とはいえ、あまり汎用的には作らずベタ書きだけど。

次のサンプルコードでは、先ほどの例と同じように MyActor に処理に時間のかかるメソッド process() を用意している。 アクターを配列で 10 個生成した上で、タスクを 20 回呼び出してみよう。 ただし、通常のスレッドプールのような形で終わったものから新しいタスクを割り当てる、というものではない。 あらかじめ、各アクターに決まった番号のタスクを割り当てるような形になっている。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division

import time
import random

import pykka


class MyActor(pykka.ThreadingActor):

    def process(self, processing_id):
        """時間がかかる処理を模したメソッド"""
        # 処理時間をランダムに用意する
        sleep_sec = random.randint(1, 10) / 10
        # time.sleep() を使って時間がかかる処理をエミュレートする
        time.sleep(sleep_sec)
        # 処理が終わったら内容を出力する
        print('Completed: ID {0} in {1} s'.format(processing_id, sleep_sec))
        # 処理の ID を返す
        return processing_id


def main():
    # スレッドプールのようにアクターを用意する
    actor_proxy_pool = [MyActor.start().proxy() for _ in range(10)]

    # それぞれのアクターに仕事を割り当てていく
    # かかる時間はランダムにする
    futures = [
        actor_proxy_pool[i % len(actor_proxy_pool)].process(i)
        for i in range(20)
    ]

    # すべての処理が完了するのを待ち合わせる
    answers = pykka.get_all(futures)

    print('*** Calculation Completed ***')

    # 完了した結果を出力する
    for answer in answers:
        print('ID: {0}'.format(answer))

    # すべてのアクターを停止させる
    pykka.ActorRegistry.stop_all()


if __name__ == '__main__':
    main()

実行結果は次の通り。 各タスクが順番バラバラで実行されているのが分かる。

$ python tpoolpykka.py
Completed: ID 3 in 0.1 s
Completed: ID 0 in 0.3 s
Completed: ID 7 in 0.3 s
Completed: ID 13 in 0.2 s
Completed: ID 1 in 0.5 s
Completed: ID 4 in 0.5 s
Completed: ID 5 in 0.6 s
Completed: ID 2 in 0.7 s
Completed: ID 8 in 0.7 s
Completed: ID 10 in 0.4 s
Completed: ID 17 in 0.4 s
Completed: ID 6 in 0.8 s
Completed: ID 9 in 1.0 s
Completed: ID 15 in 0.5 s
Completed: ID 11 in 0.6 s
Completed: ID 12 in 0.5 s
Completed: ID 16 in 0.5 s
Completed: ID 14 in 0.9 s
Completed: ID 19 in 0.5 s
Completed: ID 18 in 0.9 s
*** Calculation Completed ***
ID: 0
ID: 1
ID: 2
ID: 3
ID: 4
ID: 5
ID: 6
ID: 7
ID: 8
ID: 9
ID: 10
ID: 11
ID: 12
ID: 13
ID: 14
ID: 15
ID: 16
ID: 17
ID: 18
ID: 19

出力の最後を見て分かる通り、実行結果として得られる配列の順番は変わっていない。 これは、単純に内部でやっているのが順番に Future#get() をしていく、という処理だからだ。 とはいえ、最終的に得られる結果の順番が入れ替わらないという特性は並行処理をする上で優れた特徴だと思う。 場合によっては終わったものから返ってきて、後からソートし直さないといけないことなんかもあるし。

ラウンドロビンで処理するメッセージルータを実装してみる

次は、お試しで作ってみたメッセージルータについて。 これはほんとに試しに作ってみただけなのでプルーフオブコンセプトと思ってほしい。

次のサンプルコードでは MyActor をラップするメッセージルータの RoundRobinRouter を定義している。 RoundRobinRouter では、内部的に MyActor のプロキシを複数生成して、処理を振り分ける。 先ほどのべた書きのスレッドプールを、もうちょっと洗練させた感じ。 RoundRobinRouter は、まず処理をフォワードしたいアクターと同時並行数 (アクターをいくつ生成するか) を指定して起動する。 その上で RoundRobinRouter#forward() メソッドに呼び出したいメソッド名と引数を指定するだけ。 注意点としては、返り値が二重の Future になっているところ。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from __future__ import division

import time
import random

import pykka


class MyActor(pykka.ThreadingActor):

    def process(self, processing_id):
        """時間がかかる処理を模したメソッド"""
        # 処理時間をランダムに用意する
        sleep_sec = random.randint(1, 10) / 10
        # time.sleep() を使って時間がかかる処理をエミュレートする
        time.sleep(sleep_sec)
        # 処理が終わったら内容を出力する
        print('Completed: ID {0} in {1}s by {2}'.format(
            processing_id,
            sleep_sec,
            self.actor_urn,
        ))
        # 処理の ID を返す
        return processing_id


class RoundRobinRouter(pykka.ThreadingActor):
    """ラップしたアクターにラウンドロビンで処理を割り当てていくメッセージルータ"""

    def __init__(self, actor_type, size):
        super(RoundRobinRouter, self).__init__()

        # アクターのプロキシを生成しておく
        self._proxy_pool = [actor_type.start().proxy()
                            for _ in range(size)]
        # 次に処理させるアクターを決めるジェネレータ
        self._index_gen = self._next_actor_index()

    def forward(self, name, *args, **kwargs):
        # 次の処理させるアクターのインデックス
        actor_index = next(self._index_gen)
        # 次の処理させるアクターのプロキシ
        actor_proxy = self._proxy_pool[actor_index]
        # フォワード対象のアトリビュートをプロキシから取り出す
        attribute = getattr(actor_proxy, name)
        # 引数を渡して呼び出す
        return attribute(*args, **kwargs)

    def _next_actor_index(self):
        """次の処理対象のアクターのインデックスを生成するジェネレータ"""
        i = 0
        while True:
            yield i
            # 0, 1, 2, 3 ... 0, 1, 2, 3 と繰り返す
            i = (i + 1) % len(self._proxy_pool)


def main():
    # 五つの MyActor をラップしたルータを生成する
    router_ref = RoundRobinRouter.start(MyActor, 5)
    router_proxy = router_ref.proxy()

    futures = [router_proxy.forward('process', i)
               for i in range(10)]

    # このルータを経由させると Future の二段重ねになる
    inner_futures = pykka.get_all(futures)
    answers = pykka.get_all(inner_futures)

    print('*** Calculation Completed ***')

    for answer in answers:
        print('ID: {0}'.format(answer))

    # すべてのアクターを停止させる
    pykka.ActorRegistry.stop_all()


if __name__ == '__main__':
    main()

上記の実行結果は次の通り。 タスクを処理したアクターの識別子も同時に出力するようになっている。 見ると、タスクの識別子が 5 つの周期で同じアクターによって処理されていることが分かる。 上手く動作しているようだ。

$ python routerpykka.py
Completed: ID 1 in 0.3s by urn:uuid:a5d3f2e6-6d59-4351-9be8-e47023428e99
Completed: ID 2 in 0.6s by urn:uuid:ee6f93cf-66df-4500-b416-84defff8b195
Completed: ID 3 in 0.6s by urn:uuid:c6003463-87fc-4988-b8bb-34186e867630
Completed: ID 0 in 0.7s by urn:uuid:9a2fe8bc-49fe-4b19-96d5-6658ecae3677
Completed: ID 7 in 0.1s by urn:uuid:ee6f93cf-66df-4500-b416-84defff8b195
Completed: ID 4 in 0.8s by urn:uuid:0e3b5668-6825-43f5-bd95-9cf2c706a03d
Completed: ID 6 in 0.5s by urn:uuid:a5d3f2e6-6d59-4351-9be8-e47023428e99
Completed: ID 9 in 0.3s by urn:uuid:0e3b5668-6825-43f5-bd95-9cf2c706a03d
Completed: ID 5 in 0.8s by urn:uuid:9a2fe8bc-49fe-4b19-96d5-6658ecae3677
Completed: ID 8 in 1.0s by urn:uuid:c6003463-87fc-4988-b8bb-34186e867630
*** Calculation Completed ***
ID: 0
ID: 1
ID: 2
ID: 3
ID: 4
ID: 5
ID: 6
ID: 7
ID: 8
ID: 9

まとめ

今回は Python でアクターモデルを実装するためのフレームワークである Pykka を試してみた。

アクターモデルというのは、並行処理のプログラミングモデルの一つ。 ただし、マルチスレッドやマルチプロセスよりも、抽象度の高い概念になっている。 対比するとしたら共有メモリやロック、RPC などがそれに当たると思う。

アクターモデルの概念を理解する上で重要なのは、アクターを構成する三つの要素を理解すること。 一つ目がメッセージキュー、二つ目がメッセージハンドラ、三つ目がそれらをディスパッチするスレッド。 生成されたアクターには、それぞれにこの三つの要素が必ず備わっている。 メッセージキューとスレッドはあらかじめアクターごとに一つずつ生成されるし、メッセージハンドラは自分で定義する。

アクターモデルでは、内部で競合を起こさないようにアクターの内外が非同期の壁で分断されている。 アクターの内部に、外側から同期的に直接触れることはできない。 必ず、触るとしたら Future を経由することになる。 同期的に触っているように見えたとしても、内部的には必ず経由している。

アクターモデルで並行処理を実現するには、アクターを横にたくさん並べてタスクを割り当てていく。 ただし、それぞれのタスクはアトミックに作られていないといけない。 アクターモデルで競合が起こらないというのは、あくまで「低レイヤーで内部的には」という但し書きがつく。 また、ロックがないからといって、デッドロックが起こらないというわけでもない。

Pykka の場合、Akka の機能を全て揃えているわけではないので、割りとよく使いそうな機能でも自分で書く必要がある。 また、正直アプリケーションで直接 Pykka を使うというシチュエーションは、あまり無いのかもしれないと思った。 なぜなら Pykka の提供する API は、直接使うにはあまりに低レイヤーすぎるからだ。 なので、Pykka をベースにした何らかのフレームワークをまず書いて、その上でアプリケーションを組むことになりそうな印象を受けた。

いじょう。

Python: Keras/TensorFlow の学習を GPU で高速化する (Ubuntu 16.04 LTS)

以前、このブログで Keras/TensorFlow の学習スピードを GPU を使って速くする記事を書いた。 ただし、このとき使った OS は Mac OS X (macOS Sierra) だった。

blog.amedama.jp

とはいえ NVIDIA の dGPU を積んだ Mac がどれだけあるんだというと、正直なかなか無いと思う。 実際にやってみるとしたら Linux だよねということで、今回は Ubuntu 16.04 LTS を使う場合について書く。

インストールの手順については次の公式ドキュメントをベースに進める。

Installing TensorFlow on Ubuntu  |  TensorFlow

環境について

今回使った OS のバージョンなどは次の通り。

$ cat /etc/lsb-release
DISTRIB_ID=Ubuntu
DISTRIB_RELEASE=16.04
DISTRIB_CODENAME=xenial
DISTRIB_DESCRIPTION="Ubuntu 16.04.2 LTS"
$ uname -r
4.4.0-66-generic

また、学習に GPU を使う場合はハードウェアの環境にも制約があるので、ここで記述しておく。 まず、当たり前だけど大前提として CUDA をサポートしている NVIDIA のグラフィックカードが必ず必要になる。 ちなみに TensorFlow のドキュメントを見ると一部 Radeon も対応しているらしいことが書かれていたけど、試してはいない。

さらに、NVIDIA のグラフィックカードなら何でも良いのかというと、そういうわけでもない。 次の NVIDIA 公式ページを見てほしい。 ここには NVIDIA の各 GPU ごとの Compute Capability という数値が載っている。

CUDA GPUs | NVIDIA Developer

Compute Capability という字面だけを見ると、性能を表す数値のような印象を受けるけど、そうではない。 数値はパフォーマンスとは関係がなく、あくまで GPU の世代というか載っている機能のバージョンらしい。 問題は、この数値が少なくとも 3.0 以上ないと現行の CUDA 8.0 がサポートしていない、という点だ。

また、PyPI で配布されている TensorFlow の GPU 版バイナリパッケージは CUDA 8.0 でビルドされている。 なので、もちろん動かす環境にも CUDA 8.0 が必要になる。 以上のことから Compute Capability が 3.0 以上の GPU を用意すべきだ。

ちなみに、ドキュメントを見ると一応 TensorFlow は今のところ CUDA 7.5 もサポートしているらしい。 CUDA 7.5 であれば Compute Capability は 2.0 以上をサポートしている。 ただし、そのときは自分でソースコードからビルドしなければいけない。 それについて、この記事では扱わない。

今回 GPU には GeForce GTX 1050 Ti を用意した。 これの Compute Capability は 6.1 となっている。 具体的なグラフィックカードは次のもの。

MSI GeForce GTX 1050 Ti 4G OC グラフィックスボード VD6191

MSI GeForce GTX 1050 Ti 4G OC グラフィックスボード VD6191

GPU のチョイスについては、現時点 (2017/3) で最もコストパフォーマンスが高くなるものを選んだ。 PassMark というベンチマークのスコアを値段で割った値を考えると GeForce GTX 1050 Ti が一番高くなる。 あとは、CPU に比べると GPU の学習速度はそれこそ一桁か、下手すると二桁は違ってくる。 つまり、GPU が載っていること自体がまずは重要で、性能はその次かなという感じ。 ようするに 60 秒かかった処理が 6 秒になるのは大きく違うけど、6 秒が 3 秒になってどれだけ嬉しいか?っていう。 それもあってハイエンドモデルではなくエントリーモデルを選ぶことにした。

グラフィックカードを刺して GPU を Ubuntu 16.04 LTS が認識することを確認しておこう。

$ lspci | grep -i nvidia
01:00.0 VGA compatible controller: NVIDIA Corporation Device 1c82 (rev a1)
01:00.1 Audio device: NVIDIA Corporation Device 0fb9 (rev a1)

CUDA Toolkit をインストールする

インストールの第一歩として CUDA Toolkit をインストールする。 そのために、まずは NVIDIA 公式サイトを参照する。

developer.nvidia.com

インストール方法は色々とあるけど、ここでは CUDA のリモートリポジトリを APT に登録するやり方を取る。 これなら、後からマイナーアップデートがあったときにも自動で更新できる。 上記ウェブサイトで、次のように操作して deb ファイルのダウンロード URL を取得しよう。

Linux > x86_64 > Ubuntu > 16.04 > deb (network)

上記で得られた URL から deb ファイルをダウンロードしてくる。

$ wget http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1604/x86_64/cuda-repo-ubuntu1604_8.0.61-1_amd64.deb

deb ファイルをインストールする。 これで CUDA のリポジトリが APT パッケージシステムに登録される。

$ sudo dpkg -i cuda-repo-ubuntu1604_8.0.61-1_amd64.deb
$ sudo apt-get update

これで APT 経由で CUDA をインストールできるようになる。 ちなみに、インストールにはとても長い時間がかかるので覚悟しておこう。 先ほどのインストール方法でローカルインストールを選んだ場合には 1.9GB のファイルを落とすことになる。 かかる時間は、推して知るべし。

$ sudo apt-get -y install cuda

上記が終わると、依存パッケージとして一緒に NVIDIA のグラフィックドライバなんかも入る。

$ dpkg -l | grep nvidia-[0-9]
ii  nvidia-375                            375.26-0ubuntu1                          amd64        NVIDIA binary driver - version 375.26
ii  nvidia-375-dev                        375.26-0ubuntu1                          amd64        NVIDIA binary Xorg driver development files

cuDNN をインストールする

次に cuDNN をインストールする。

cuDNN のインストールには、あらかじめ NVIDIA のサイトでデベロッパー登録が必要になる。 登録を済ませると cuDNN をダウンロードできるようになる。

https://developer.nvidia.com/rdp/cudnn-download

cuDNN の現行のバージョンは 5.1 となっている。 ダウンロードに使うパッケージは CUDA Toolkit のバージョンと組み合わせになっている。 今回であれば v5.1 + CUDA 8.0 のパッケージを選ぼう。

Download cuDNN v5.1 (Jan 20, 2017), for CUDA 8.0 > cuDNN v5.1 Library for Linux

ダウンロードしてきたら、それを CUDA のインストールされたディレクトリに放り込む。

$ sudo tar xvf cudnn-8.0-linux-x64-v5.1.tgz -C /usr/local
cuda/include/cudnn.h
cuda/lib64/libcudnn.so
cuda/lib64/libcudnn.so.5
cuda/lib64/libcudnn.so.5.1.10
cuda/lib64/libcudnn_static.a

libcupti-dev をインストールする

次に The CUDA Profiling Tools Interface もインストールしておく。

$ sudo apt-get -y install libcupti-dev

ここまで終わったら、一旦再起動しておこう。

$ sudo shutdown -r now

Keras/TensorFlow (GPU)

さて、次はいよいよ Keras/TensorFlow をインストールする…と言いたいところだけど、その前に。 システムの Python 実行環境を汚したくないので Python の仮想環境を作れるようにしよう。

ここでは virtualenvwrapper を使う。

$ sudo apt-get -y install virtualenvwrapper

シェルを起動し直す。

$ exec $SHELL

これで Python 仮想環境が作れるようになった。 tensorflow-with-gpu という名前で Python3 の仮想環境を作る。

$ mkvirtualenv tensorflow-with-gpu -p $(which python3)

これで最低限のパッケージが入った、システムから独立した仮想環境ができた。

(tensorflow-with-gpu) $ pip list --format=columns
Package       Version
------------- --------
appdirs       1.4.3
packaging     16.8
pip           9.0.1
pkg-resources 0.0.0
pyparsing     2.2.0
setuptools    34.3.2
six           1.10.0
wheel         0.30.0a0

ここに Keras と GPU 版 TensorFlow をインストールする。

(tensorflow-with-gpu) $ pip install keras tensorflow-gpu

インストールすると、こんな感じ。

(tensorflow-with-gpu) $ pip list --format=columns | grep -i -e keras -e tensorflow
Keras          1.2.2
tensorflow-gpu 1.0.1

ちなみに、以降はシェルのプレフィックスを表記しないけど Python 仮想環境上で実行し続けている。

GPU で学習させる

まずはベンチマーク用のアプリケーションで使うデータセット (MNIST) をダウンロードしよう。 Python の REPL を起動する。

$ python

Keras の mnist パッケージをインポートする。

>>> from keras.datasets import mnist
Using TensorFlow backend.
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcublas.so.8.0 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcudnn.so.5 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcufft.so.8.0 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcuda.so.1 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcurand.so.8.0 locally

load_data() 関数でデータセットのダウンロードが始まる。 これには少し時間がかかる。

>>> mnist.load_data()
Downloading data from https://s3.amazonaws.com/img-datasets/mnist.pkl.gz
...
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]]], dtype=uint8), array([7, 2, 1, ..., 4, 5, 6], dtype=uint8)))

終わったら REPL から抜ける。

>>> exit()

続いて学習のベンチマークに使うアプリケーションをダウンロードしてくる。

$ curl -O https://raw.githubusercontent.com/fchollet/keras/master/examples/mnist_cnn.py
$ echo 'K.clear_session()' >> mnist_cnn.py

ベンチマーク用のアプリケーションを実行する。 これは MNIST データセットを CNN で認識するものになっている。 初めての実行だと GPU の環境を整えるので学習を始まるまで時間がかかるかも。

$ time python mnist_cnn.py
Using TensorFlow backend.
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcublas.so.8.0 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcudnn.so.5 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcufft.so.8.0 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcuda.so.1 locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcurand.so.8.0 locally
X_train shape: (60000, 28, 28, 1)
60000 train samples
10000 test samples
Train on 60000 samples, validate on 10000 samples
Epoch 1/12
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE3 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.1 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.
I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:910] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero
I tensorflow/core/common_runtime/gpu/gpu_device.cc:885] Found device 0 with properties:
name: GeForce GTX 1050 Ti
major: 6 minor: 1 memoryClockRate (GHz) 1.455
pciBusID 0000:01:00.0
Total memory: 3.94GiB
Free memory: 3.84GiB
I tensorflow/core/common_runtime/gpu/gpu_device.cc:906] DMA: 0
I tensorflow/core/common_runtime/gpu/gpu_device.cc:916] 0:   Y
I tensorflow/core/common_runtime/gpu/gpu_device.cc:975] Creating TensorFlow device (/gpu:0) -> (device: 0, name: GeForce GTX 1050 Ti, pci bus id: 0000:01:00.0)
60000/60000 [==============================] - 8s - loss: 0.3661 - acc: 0.8871 - val_loss: 0.0864 - val_acc: 0.9726
Epoch 2/12
60000/60000 [==============================] - 6s - loss: 0.1337 - acc: 0.9604 - val_loss: 0.0613 - val_acc: 0.9806
Epoch 3/12
60000/60000 [==============================] - 6s - loss: 0.1046 - acc: 0.9694 - val_loss: 0.0524 - val_acc: 0.9838
Epoch 4/12
60000/60000 [==============================] - 6s - loss: 0.0879 - acc: 0.9737 - val_loss: 0.0428 - val_acc: 0.9858
Epoch 5/12
60000/60000 [==============================] - 6s - loss: 0.0755 - acc: 0.9773 - val_loss: 0.0393 - val_acc: 0.9870
Epoch 6/12
60000/60000 [==============================] - 6s - loss: 0.0683 - acc: 0.9801 - val_loss: 0.0368 - val_acc: 0.9875
Epoch 7/12
60000/60000 [==============================] - 6s - loss: 0.0648 - acc: 0.9806 - val_loss: 0.0367 - val_acc: 0.9888
Epoch 8/12
60000/60000 [==============================] - 6s - loss: 0.0593 - acc: 0.9820 - val_loss: 0.0353 - val_acc: 0.9889
Epoch 9/12
60000/60000 [==============================] - 6s - loss: 0.0548 - acc: 0.9833 - val_loss: 0.0328 - val_acc: 0.9890
Epoch 10/12
60000/60000 [==============================] - 6s - loss: 0.0529 - acc: 0.9845 - val_loss: 0.0316 - val_acc: 0.9897
Epoch 11/12
60000/60000 [==============================] - 6s - loss: 0.0508 - acc: 0.9848 - val_loss: 0.0308 - val_acc: 0.9902
Epoch 12/12
60000/60000 [==============================] - 6s - loss: 0.0484 - acc: 0.9852 - val_loss: 0.0298 - val_acc: 0.9899
Test score: 0.0297728348608
Test accuracy: 0.9899

real  1m28.273s
user  1m39.684s
sys   0m8.740s

全体の実行が一分半ほどで終わった。 一回のエポックは大体 6 秒ほどで終わっている。

ちなみに GPU の仕事っぷりは nvidia-smi コマンドから見られる。

$ nvidia-smi
Sun Mar 12 20:24:56 2017
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 375.26                 Driver Version: 375.26                    |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|===============================+======================+======================|
|   0  GeForce GTX 105...  Off  | 0000:01:00.0      On |                  N/A |
| 35%   38C    P0    58W /  75W |   3862MiB /  4037MiB |     77%      Default |
+-------------------------------+----------------------+----------------------+

+-----------------------------------------------------------------------------+
| Processes:                                                       GPU Memory |
|  GPU       PID  Type  Process name                               Usage      |
|=============================================================================|
|    0      1510    G   /usr/lib/xorg/Xorg                              59MiB |
|    0      3374    C   python                                        3799MiB |
+-----------------------------------------------------------------------------+

毎秒ごとに表示を更新したいときは、こんな感じで。

$ nvidia-smi -l 1

CPU で学習させる

ちなみに対比として CPU を使ったときの例も載せておく。

使った CPU は次の通り。 第一世代 Core i7 なので、なかなか古い。 現行世代なら Core i3 くらいの性能らしい。

$ cat /proc/cpuinfo | grep -i name | head -n 1
model name      : Intel(R) Core(TM) i7 CPU         860  @ 2.80GHz
$ cat /proc/cpuinfo | grep -i name | wc -l
8

余談ながら、当然チップセットやマザーボードも古いので今回使ったグラフィックカードを認識するかドキドキした。

先ほどインストールした GPU 版 TensorFlow をアンインストールして CPU 版 TensorFlow をインストールする。

$ pip uninstall -y tensorflow-gpu
$ pip install tensorflow

先ほどと同じようにベンチマーク用のアプリケーションを実行する。

$ time python mnist_cnn.py
Using TensorFlow backend.
X_train shape: (60000, 28, 28, 1)
60000 train samples
10000 test samples
Train on 60000 samples, validate on 10000 samples
Epoch 1/12
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE3 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.1 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.
60000/60000 [==============================] - 117s - loss: 0.3733 - acc: 0.8858 - val_loss: 0.0885 - val_acc: 0.9723
Epoch 2/12
60000/60000 [==============================] - 117s - loss: 0.1347 - acc: 0.9596 - val_loss: 0.0635 - val_acc: 0.9799
Epoch 3/12
60000/60000 [==============================] - 115s - loss: 0.1031 - acc: 0.9693 - val_loss: 0.0565 - val_acc: 0.9818
Epoch 4/12
60000/60000 [==============================] - 115s - loss: 0.0887 - acc: 0.9740 - val_loss: 0.0448 - val_acc: 0.9852
Epoch 5/12
60000/60000 [==============================] - 114s - loss: 0.0780 - acc: 0.9773 - val_loss: 0.0415 - val_acc: 0.9868
Epoch 6/12
60000/60000 [==============================] - 113s - loss: 0.0707 - acc: 0.9788 - val_loss: 0.0376 - val_acc: 0.9870
Epoch 7/12
60000/60000 [==============================] - 112s - loss: 0.0651 - acc: 0.9808 - val_loss: 0.0352 - val_acc: 0.9893
Epoch 8/12
60000/60000 [==============================] - 112s - loss: 0.0604 - acc: 0.9818 - val_loss: 0.0352 - val_acc: 0.9888
Epoch 9/12
60000/60000 [==============================] - 112s - loss: 0.0555 - acc: 0.9840 - val_loss: 0.0332 - val_acc: 0.9884
Epoch 10/12
60000/60000 [==============================] - 112s - loss: 0.0542 - acc: 0.9835 - val_loss: 0.0322 - val_acc: 0.9895
Epoch 11/12
60000/60000 [==============================] - 111s - loss: 0.0495 - acc: 0.9853 - val_loss: 0.0319 - val_acc: 0.9897
Epoch 12/12
60000/60000 [==============================] - 110s - loss: 0.0477 - acc: 0.9857 - val_loss: 0.0326 - val_acc: 0.9893
Test score: 0.032642004689
Test accuracy: 0.9893

real  22m54.715s
user  123m6.744s
sys   46m53.296s

今度は 23 分ほどかかった。

GPU が 1 分 28 秒で終わったのに対して CPU では 22 分 54 秒かかっている。 つまり、GPU を使うことで学習が 15.6 倍高速化できた。 一回のエポック当たりでは 19.1 倍速くなっている。

まとめ

今回は OS として Ubuntu 16.04 LTS を使って Keras/TensorFlow の学習を GPU で高速化してみた。 CUDA がサポートしている GPU は Compute Capability という数値で表現されている。 この数値が大きいものほど、新しい世代のものになっている。 使う GPU については NVIDIA で、少なくとも Compute Capability が 3.0 以上のものを選ぼう。 もちろん、将来の CUDA で最小要件の Compute Capability は上がっていくはずなので、なるべく大きい方が良い。

いじょう。

MSI GeForce GTX 1050 Ti 4G OC グラフィックスボード VD6191

MSI GeForce GTX 1050 Ti 4G OC グラフィックスボード VD6191

Python: Keras/TensorFlow の学習を CPU の拡張命令で高速化する (Mac OS X)

今回のネタは TensorFlow を使っていると、いつも目にしていた警告について。

それは、次のようなもの。

W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.1 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX instructions, but these are available on your machine and could speed up CPU computations.

どうやら TensorFlow は CPU の拡張命令を使うことができるらしい。 しかし、PyPI 経由で一般に配布されているバイナリパッケージでは、それらは有効になっていない。 もちろん、同じ x86 系 CPU であってもどんな拡張命令セットを持っているかは製品によって異なるので、これは当然の判断だろう。

そこで、今回は CPU の拡張命令を有効にして TensorFlow をソースコードからコンパイルしてみることにした。 そうすることで、どれくらい学習を高速化できるだろうか?というのが主題となる。

先に結論から書いてしまうと、今回使った環境では大体 30% くらい学習が速くなった。 ぶっちゃけ、これくらいであれば素直に GPU を使える環境を用意した方が良いと思う。

blog.amedama.jp

ただし、自分で使う環境向けにチューンするという面において、自前でコンパイルするのは良い選択肢かもしれないと感じた。

今回使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.12.3
BuildVersion:   16D32
$ python --version
Python 3.6.0

あらかじめ、計算速度のベンチマークに使うアプリケーションをダウンロードしておこう。

$ curl -O https://raw.githubusercontent.com/fchollet/keras/master/examples/mnist_cnn.py
$ echo 'K.clear_session()' >> mnist_cnn.py

拡張命令を使わない場合

まずは、CPU の拡張命令を使わない場合にかかる時間を調べておこう。

PyPI から汎用のバイナリパッケージで TensorFlow をインストールする。 ベンチマークにするアプリケーションを実行するために Keras も入れておく。

$ pip install keras tensorflow

インストールできたら Python の REPL を起動する。

$ python

REPL が起動したら MNIST のデータセットをダウンロードしておこう。 これには少し時間がかかる。

>>> from keras.datasets import mnist
Using TensorFlow backend.
>>> mnist.load_data()
Downloading data from https://s3.amazonaws.com/img-datasets/mnist.pkl.gz

終わったら REPL から抜けておく。

>>> exit()

それでは、まずは CPU の拡張命令を使わない場合にかかる時間を測ろう。 time コマンド経由でベンチマーク用のアプリケーションを実行する。

$ time python mnist_cnn.py
Using TensorFlow backend.
X_train shape: (60000, 28, 28, 1)
60000 train samples
10000 test samples
Train on 60000 samples, validate on 10000 samples
Epoch 1/12
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.1 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX instructions, but these are available on your machine and could speed up CPU computations.
60000/60000 [==============================] - 92s - loss: 0.3741 - acc: 0.8859 - val_loss: 0.0880 - val_acc: 0.9718
...(省略)...
python mnist_cnn.py  2452.87s user 859.80s system 302% cpu 18:15.11 total

CPU の拡張命令を使わない場合には 2450 秒かかった。

拡張命令を使う場合

次は CPU の拡張命令を使う場合について。 警告メッセージを見る限り、今回使うマシンでは SSE4.1 SSE4.2 AVX という拡張命令が使える。

それらの拡張命令を有効にするためにソースコードから TensorFlow をインストールする。 ソースコードからのインストールには次のドキュメントを参照すると良い。 Download and Setup  |  TensorFlow

まずは、Homebrew を使って依存パッケージをインストールする。

$ brew install bazel swig

次に TensorFlow のソースコードをクローンする。

$ git clone https://github.com/tensorflow/tensorflow.git
$ cd tensorflow

先ほどベンチマークしたのと同じバージョンのタグをチェックアウトしよう。

$ git checkout v1.0.1

次にインストールするための環境を configure で設定していく。 今回は CUDA を使わないので、さほど難しくはない。

$ ./configure
Please specify the location of python. [Default is /Users/amedama/.virtualenvs/py36/bin/python]:
Please specify optimization flags to use during compilation [Default is -march=native]:
Do you wish to use jemalloc as the malloc implementation? (Linux only) [Y/n] n
jemalloc disabled on Linux
Do you wish to build TensorFlow with Google Cloud Platform support? [y/N]
No Google Cloud Platform support will be enabled for TensorFlow
Do you wish to build TensorFlow with Hadoop File System support? [y/N]
No Hadoop File System support will be enabled for TensorFlow
Do you wish to build TensorFlow with the XLA just-in-time compiler (experimental)? [y/N]
No XLA support will be enabled for TensorFlow
Found possible Python library paths:
  /Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages
Please input the desired Python library path to use.  Default is [/Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages]

Using python library path: /Users/amedama/.virtualenvs/py36/lib/python3.6/site-packages
Do you wish to build TensorFlow with OpenCL support? [y/N]
No OpenCL support will be enabled for TensorFlow
Do you wish to build TensorFlow with CUDA support? [y/N]
No CUDA support will be enabled for TensorFlow
Configuration finished
Extracting Bazel installation...
..........................................................
INFO: Starting clean (this may take a while). Consider using --expunge_async if the clean takes more than several minutes.
.........................................................
INFO: All external dependencies fetched successfully.

設定が終わったら次に bazel を使ってビルドする。 ここで使用する CPU の拡張命令を指定する。 --copt=-m の後ろに拡張命令の名前を入力しよう。 この処理には時間がかかるので気長に待つ。

$ bazel build -c opt --copt=-mavx --copt=-msse4.1 --copt=-msse4.2 //tensorflow/tools/pip_package:build_pip_package

もし、ここで指定しなければ CPU 拡張命令を使わない汎用なパッケージになる。

続けて Wheel パッケージをビルドする。 二番目に渡している /tmp/tensorflow_pkg が Wheel パッケージを保存する場所になる。

$ bazel-bin/tensorflow/tools/pip_package/build_pip_package /tmp/tensorflow_pkg

これで CPU の拡張命令が有効になった Wheel パッケージができた。

$ ls /tmp/tensorflow_pkg
tensorflow-1.0.1-cp36-cp36m-macosx_10_12_x86_64.whl

既にインストールされている汎用なパッケージをアンインストールした上で、インストールしよう。

$ pip uninstall -y tensorflow
$ pip install /tmp/tensorflow_pkg/tensorflow-1.0.1-cp36-cp36m-macosx_10_12_x86_64.whl
$ pip list --format=columns | grep -i tensorflow
tensorflow       1.0.1

さて、先ほどと同じようにベンチマークとなるアプリケーションを実行してみよう。 どれくらい速くなるだろうか。

$ time python mnist_cnn.py
Using TensorFlow backend.
X_train shape: (60000, 28, 28, 1)
60000 train samples
10000 test samples
Train on 60000 samples, validate on 10000 samples
Epoch 1/12
60000/60000 [==============================] - 75s - loss: 0.3735 - acc: 0.8859 - val_loss: 0.0888 - val_acc: 0.9724
...(省略)...
python mnist_cnn.py  1871.67s user 901.15s system 296% cpu 15:34.90 total

今度は 1870 秒で終わった。

結果は CPU の拡張命令を使わない場合が 2450 秒で、使った場合が 1870 秒となった。 つまり、処理速度が 30% ほど速くなった。 また、一回のエポックでいうと、だいたい 20% ほど速くなっているようだ。

まとめ

今回は、普段 TensorFlow を使っているときに表示される警告が気になってソースコードからビルドを試してみた。 PyPI で配布されている TensorFlow は CPU に依存しない汎用なバイナリになっている。 もし、自分の CPU に合わせて拡張命令を使ったバイナリを作りたいときは、ソースコードからコンパイルする必要がある。 CPU の拡張命令を有効にすると、今回のケースでは 30% ほど学習が速くなった。 GPU を使った学習の高速化に比べると、これは微々たるもののように感じる。 しかし、自分の環境に向けてチューンされたバイナリを作るという意味では、ソースコードからのコンパイルは有用かもしれない。

いじょう。

Python: python-fire の CLI 自動生成を試す

今回は Google が公開した python-fire というパッケージを試してみた。 python-fire では、クラスやモジュールを渡すことで、定義されている関数やメソッドを元に CLI を自動で生成してくれる。

ただし、一つ注意すべきなのは、できあがる CLI はそこまで親切な作りではない、という点だ。 実際にユーザに提供するような CLI を実装するときは、従来通り Click のようなフレームワークを使うことになるだろう。 では python-fire はどういったときに活躍するかというと、これは開発時のテストだと思う。 実装した内容をトライアンドエラーするための CLI という用途であれば python-fire は非常に強力なパッケージだと感じた。

今回使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.12.3
BuildVersion:   16D32
$ python --version
Python 3.6.0

インストール

インストールは Python のパッケージツールの pip を使ってできる。

$ pip install fire

もし pip がインストールされていないときは、あらかじめ入れておこう。

$ curl https://bootstrap.pypa.io/get-pip.py | sudo python

あと、これは完全に蛇足だけどシステムの Python 実行環境にそのまま入れるのはおすすめしない。 色んなパッケージを試すときは、システムからは独立した Python 仮想環境を作って入れた方が良い。 以下は、例えば virtualenv を使ったやり方。

$ sudo pip install virtualenv
$ mkdir -p venv
$ virtualenv venv

これで、最低限のパッケージの入った Python 仮想環境ができる。

(venv) $ pip list --format=columns
Package    Version
---------- -------
appdirs    1.4.3
packaging  16.8
pip        9.0.1
pyparsing  2.2.0
setuptools 34.3.1
six        1.10.0
wheel      0.29.0

python-fire をインストールすると、次のようなパッケージが入る。 意外と依存パッケージが多い。

(venv) $ pip install fire
(venv) $ pip list --format=columns
Package          Version
---------------- -------
appdirs          1.4.2
appnope          0.1.0
decorator        4.0.11
fire             0.1.0
ipython          5.3.0
ipython-genutils 0.1.0
packaging        16.8
pexpect          4.2.1
pickleshare      0.7.4
pip              9.0.1
prompt-toolkit   1.0.13
ptyprocess       0.5.1
Pygments         2.2.0
pyparsing        2.2.0
setuptools       34.3.1
simplegeneric    0.8.1
six              1.10.0
traitlets        4.3.2
wcwidth          0.1.7
wheel            0.29.0

基本的な使い方

まずは python-fire の基本的な使い方から見ていく。 とはいっても、その使い方は至ってシンプル。 例えば CLI を自動生成したいクラスがあるなら、それを fire.Fire() コマンドに渡すだけ。

次のコマンドを実行すると helloworld.py というファイルでサンプルコードが保存される。

$ cat << 'EOF' > helloworld.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-


class MyClass(object):

    def greet(self):
        print('Hello, World!')


if __name__ == '__main__':
    # fire をインポートする
    import fire
    # コマンドラインで呼び出したいクラスを指定する
    fire.Fire(MyClass)
EOF

上記のサンプルコードでは MyClass というクラスを python-fire に渡している。 このクラスには greet() というメソッドがあって、文字列をプリントする内容になっている。

上記を何の引数も渡さずに実行してみよう。

$ python helloworld.py
Type:        MyClass
String form: <__main__.MyClass object at 0x10bb64f28>
File:        ~/Documents/temporary/helloworld.py

Usage:       helloworld.py
             helloworld.py greet

すると Usage が出力されることがわかる。

上記の出力に沿ってオプションを渡してみよう。 具体的には、MyClass の持つメソッド名である greet を指定する。

$ python helloworld.py greet
Hello, World!

これだけで MyClass#greet() の処理が呼び出された。

関数・メソッドに引数を渡す

先ほどの例ではメソッドの呼び出しを試した。 ただ、そのメソッドには引数がなかった。 次は引数のあるメソッドを呼び出すときについて見てみる。

以下のコマンドを実行すると、引数のあるメソッドを定義したサンプルコードができる。

$ cat << 'EOF' > calculate.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-


class MyClass(object):

    def add(self, a, b=1):
        return a + b


if __name__ == '__main__':
    import fire
    fire.Fire(MyClass)
EOF

上記では MyClass#add() メソッドが引数として ab を受け取る。 そのうち b に関してはデフォルト引数が指定されている。

上記で用意したサンプルコードを先ほどと同じようにメソッド名だけ入れて実行してみよう。 この場合、必要な引数 a が入力されていないためトレースが表示される。

$ python calculate.py add
Fire trace:
1. Initial component
2. Instantiated class "MyClass" (calculate.py:5)
3. Accessed property "add" (calculate.py:7)
4. ('The function received no value for the required argument:', 'a')

Type:        method
String form: <bound method MyClass.add of <__main__.MyClass object at 0x109fba518>>
File:        ~/Documents/temporary/calculate.py
Line:        7

Usage:       calculate.py add A [B]
             calculate.py add --a A [--b B]

Usage を見ると、メソッド名である add に続いて引数を渡せば良いことがわかる。

上記にもとづいてメソッド名に続いて引数を渡してみよう。 それぞれ引数の ab に対応する。

$ python calculate.py add 1 2
3

引数の b はデフォルト引数が指定されているので省略できる。

$ python calculate.py add 1
2

引数に渡す内容の指定を順不定にしたいときは -- を使って引数名を指定してやる。

$ python calculate.py add --a=1 --b=2
3

コンストラクタに引数を渡す

先ほどまでの例では CLI を自動生成したいクラスのコンストラクタに引数がなかった。 現実には、何も引数を取らないコンストラクタの方が珍しいと思う。 そこで、次はコンストラクタに引数を取る場合を試してみる。

次のコマンドを実行するとコンストラクタに引数を取るサンプルコードができる。

$ cat << 'EOF' > constructor.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-


class MyClass(object):

    def __init__(self, msg):
        """コンストラクタに引数 msg を取る"""
        self.msg = msg

    def greet(self):
        print(self.msg)


if __name__ == '__main__':
    import fire
    fire.Fire(MyClass)
EOF

先ほどの例と同じように引数を何も指定せず実行してみよう。 するとコンストラクタの引数が指定されていないという内容でトレースが表示される。

$ python constructor.py
Fire trace:
1. Initial component
2. ('The function received no value for the required argument:', 'msg')

Type:        type
String form: <class '__main__.MyClass'>
File:        ~/Documents/temporary/constructor.py
Line:        5

Usage:       constructor.py MSG
             constructor.py --msg MSG

どうやらコンストラクタの引数をまずは指定する必要があるらしい。

上記の表示に沿って、コンストラクタの --msg 引数を指定してみよう。

$ python constructor.py --msg="Hello World"
Type:           MyClass
String form:    <__main__.MyClass object at 0x102b303c8>
File:           ~/Documents/temporary/constructor.py
Init docstring: コンストラクタに引数 msg を取る

Usage:          constructor.py --msg='Hello World' 
                constructor.py --msg='Hello World' greet
                constructor.py --msg='Hello World' msg

すると、今度は実行できるメソッド名が表示されるようになった。 メソッド greet() の下にある msg はインスタンスのメンバとして引数を保存しているため表示されているんだろう。

上記の表示に沿ってコンストラクタの引数を入力しながらメソッドを指定する。

$ python constructor.py --msg='Hello World' greet
Hello World

これでコンストラクタに引数を取るクラスであっても実行できるようになった。

関数から自動生成する

これまでの例ではクラスに対して CLI を自動生成する例だった。 次は関数に対して試してみることにしよう。

次のコマンドを実行すると関数を使うパターンのサンプルコードができる。

$ cat << 'EOF' > function.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-


def greet(msg='Hello, World!'):
    print(msg)


if __name__ == '__main__':
    import fire
    fire.Fire(greet)
EOF

とはいえ、やっていることはこれまでと変わらない。 渡すものがクラスの代わりに関数になっただけ。

今度も、特に何も指定せずまずは実行してみよう。

$ python function.py
Hello, World!

今回に関しては、これだけで事足りてしまった。

ちなみに、上記の場合はデフォルト引数が指定されていたので動いた。 もちろん、次のようにして引数を指定することもできる。

$ python function.py --msg="Hajimemashite Sekai"
Hajimemashite Sekai

ちなみに、ここで面白い挙動に気づいた。 引数を指定するときに内容にカンマが含まれていると、型が自動的にタプルになるようだ。

$ python function.py --msg="Hajimemashite, Sekai"
('Hajimemashite', 'Sekai')

型の自動判別

先ほどの例を元に、どういった内容を引数に指定すると、どの型と判別するのか調べてみた。

次のサンプルコードでは受け取った引数の型を出力する。

$ cat << 'EOF' > argtype.py 
#!/usr/bin/env python
# -*- coding: utf-8 -*-


def print_type(obj):
    print(type(obj))


if __name__ == '__main__':
    import fire
    fire.Fire(print_type)
EOF

色んな内容を渡してみた結果が次の通り。 いきなりタプルになったのは驚いたけど、文字列を指定したいなら '' で囲むのがセオリーなようだ。

$ python argtype.py --obj="foo"
<class 'str'>

$ python argtype.py --obj="1"  
<class 'int'>
$ python argtype.py --obj="1.0"
<class 'float'>
$ python argtype.py --obj="'1'"
<class 'str'>

$ python argtype.py --obj="1,2"
<class 'tuple'>
$ python argtype.py --obj="'1,2'"
<class 'str'>

$ python argtype.py --obj="[1, 2]"
<class 'list'>
$ python argtype.py --obj="{1, 2}"
<class 'set'>
$ python argtype.py --obj="{1: 2}"
<class 'dict'>

モジュールから自動生成する

次は、おそらく一番使う場面が多そうなやり方。 モジュールをそのまま指定して CLI を自動生成してしまうやり方。

次のコマンドを実行するとモジュールをそのまま指定するサンプルコードができる。 この場合は、特に何も指定せず fire.Fire() をインスタンス化すれば良い。

$ cat << 'EOF' > module.py
#!/usr/bin/env python
# -*- coding: utf-8 -*-

class MyClass(object):

    def greet(self):
        print('Hello, World!')


def greet(msg='Hello, World!'):
    print(msg)


def add(a, b):
    return a + b


if __name__ == '__main__':
    import fire
    fire.Fire()
EOF

モジュールには MyClass クラスや greet() および add() 関数が定義されている。

次も、特に何も指定せずまずは実行してみよう。 ただ、これだけだと何が指定できるのかよく分からない出力になる。

$ python module.py
MyClass: <class '__main__.MyClass'>
greet:   <function greet at 0x10e53ee18>
add:     <function add at 0x10e67d378>
fire:    <module 'fire' from '/Users/amedama/.virtualenvs/fire/lib/python3.6/site-packages/fire/__init__.py'>

そこで --help オプションを指定してみよう。 間にある -- は python-fire 自体に渡す引数と、自動生成した CLI に渡す引数を区別するために使う。

$ python module.py -- --help
Type:        dict
String form: {'__name__': '__main__', '__doc__': None, '__package__': None, '__loader__': <_frozen_importlib_e <...> ule 'fire' from '/Users/amedama/.virtualenvs/fire/lib/python3.6/site-packages/fire/__init__.py'>}
Length:      13

Usage:       module.py 
             module.py MyClass
             module.py greet
             module.py add
             module.py fire

すると、実行できる内容が出力されるようになった。

関数に関しては、これまでと特に変わらないやり方で実行できる。

$ python module.py greet
Hello, World!
$ python module.py add 1 2
3

インスタンスメソッドについては、クラス名とメソッド名を続けて入力すれば良い。

$ python module.py MyClass greet
Hello, World!

インタラクティブシェルに入る

最初に依存パッケージの中に ipython があったことに気づいていたかもしれない。 ipython は高機能な Python の REPL で python-fire の引数に --interactive を指定すると、それが起動するようになっている。 起動した ipython では特にインポートなどをすることなくモジュールの中で定義されているクラスなどを使うことができる。

$ python helloworld.py -- --interactive
Fire is starting a Python REPL with the following objects:
Modules: fire
Objects: MyClass, component, helloworld.py, result, trace

Python 3.6.0 (default, Feb 25 2017, 20:17:10)
Type "copyright", "credits" or "license" for more information.

IPython 5.3.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]: MyClass
Out[1]: __main__.MyClass

In [2]: MyClass().greet()
Hello, World!

In [3]: exit()

通常通り ipython を指定しただけでは、こうはならない。 インポートしていないオブジェクトを使おうとするとエラーになる。

$ ipython
Python 3.6.0 (default, Feb 25 2017, 20:17:10)
Type "copyright", "credits" or "license" for more information.

IPython 5.3.0 -- An enhanced Interactive Python.
?         -> Introduction and overview of IPython's features.
%quickref -> Quick reference.
help      -> Python's own help system.
object?   -> Details about 'object', use 'object??' for extra details.

In [1]: MyClass
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-1-d1eb688265f8> in <module>()
----> 1 MyClass

NameError: name 'MyClass' is not defined

本来であれば次のようにしてインポートしてやる必要がある。 python-fire 経由で起動すると、この手間をはぶくことができるようだ。

In [2]: from helloworld import MyClass

In [3]: MyClass
Out[3]: helloworld.MyClass

呼び出しをトレースする

これまでの例でもエラーになったときに何度か表示されていたけど、呼び出し内容をトレースすることもできる。 これには python-fire の引数として --trace オプションを指定する。

$ python helloworld.py greet -- --trace
Fire trace:
1. Initial component
2. Instantiated class "MyClass" (helloworld.py:5)
3. Accessed property "greet" (helloworld.py:7)

上記のように各ステップでどういったことが実行されているのかが表示される。

補完スクリプトを出力する

python-fire では bash 用の補完スクリプトを出力することもできる。 それには、次のようにして --completion を python-fire の引数に渡す。

$ python helloworld.py -- --completion
# bash completion support for helloworld.py
# DO NOT EDIT.
# This script is autogenerated by fire/completion.py.

_complete-helloworldpy()
{
  local start cur opts
  COMPREPLY=()
  start="${COMP_WORDS[@]:0:COMP_CWORD}"
  cur="${COMP_WORDS[COMP_CWORD]}"

  opts=""


  if [[ "$start" == "helloworld.py" ]] ; then
    opts="greet"
  fi

  if [[ "$start" == "helloworld.py greet" ]] ; then
    opts="--self"
  fi

  COMPREPLY=( $(compgen -W "${opts}" -- ${cur}) )
  return 0
}

complete -F _complete-helloworldpy helloworld.py

ここで出力される内容を bash の設定ファイルに入れて読み込んでやるとシェルで補完が効くようになる。

$ python helloworld.py -- --completion >> ~/.bashrc
$ source ~/.bashrc

ただし、出力される補完スクリプトはパスが通っている場所で、かつスクリプトに実行権限がついていることを前提になっている。 そのため、上記のように python コマンド経由で実行する限りでは補完が効かない。 もし使うとしたら、こんな感じにしないとダメそう。

$ chmod +x helloworld.py
$ sudo mv helloworld.py /usr/local/bin/

ちなみに bashcompinit を有効にして試してみたけど zsh では動かなかった。

詳細を出力する

これまで見てきたように --help オプションを指定すると詳しい使い方が出力された。 次のように --verbose オプションを指定すると、特殊メソッドやアトリビュートを含むさらに詳しい内容が出力される。

$ python helloworld.py -- --verbose
Type:        MyClass
String form: <__main__.MyClass object at 0x1055cc4a8>
File:        ~/Documents/temporary/helloworld.py

Usage:       helloworld.py
             helloworld.py __class__
             helloworld.py __delattr__
             helloworld.py __dict__
             helloworld.py __dir__
             helloworld.py __doc__
             helloworld.py __eq__
             helloworld.py __format__
             helloworld.py __ge__
             helloworld.py __getattribute__
             helloworld.py __gt__
             helloworld.py __hash__
             helloworld.py __init__
             helloworld.py __init_subclass__
             helloworld.py __le__
             helloworld.py __lt__
             helloworld.py __module__
             helloworld.py __ne__
             helloworld.py __new__
             helloworld.py __reduce__
             helloworld.py __reduce_ex__
             helloworld.py __repr__
             helloworld.py __setattr__
             helloworld.py __sizeof__
             helloworld.py __str__
             helloworld.py __subclasshook__
             helloworld.py __weakref__
             helloworld.py greet

まとめ

今回は Google の作った CLI を自動生成するパッケージである python-fire を試してみた。 生成される CLI は割りと雑だけど、開発者が使う分には全く問題ないレベルだと思う。 これまでは挙動を軽く試すのにもスクリプトを毎回書き換えたり ipython を起動したりと面倒が多かった。 一つ一つの作業にかかる時間は短くとも、長い目で見れば多くの時間を浪費していることだろう。 そんなとき python-fire を使えば、スクリプトに数行を追加するだけでその手間をはぶくことができるのは大きいと感じた。

いじょう。

Python: Keras/TensorFlow の学習を GPU で高速化する (Mac OS X)

Keras というのは Python を使ってニューラルネットワークを組むためのフレームワーク。 Python でニューラルネットワークのフレームワークというと、他にも TensorFlow とか Chainer なんかが有名どころ。 Keras はそれらに比べると、より高い抽象度の API を提供しているところが特徴みたい。 実のところ Keras はデフォルトで TensorFlow をバックエンドとして動作する。 バックエンドとしては、他にも Theano が選べるらしい。

今回は Keras で組んだニューラルネットワークを GPU で学習させてみることにした。 そのとき CPU と比べて、どれくらい速くなるかを試してみたい。

使った環境は次の通り。

$ sw_vers
ProductName:    Mac OS X
ProductVersion: 10.12.3
BuildVersion:   16D32
$ python --version
Python 3.5.3

当然だけどマシンには dGPU のグラフィックカードが載っている必要がある。 今回は GeForce GTX 675MX を使った。

まずは CPU で学習させるための下準備

まずは Python のパッケージ管理システムである pip を使って Keras をインストールする。

$ pip install keras

もし pip がそもそも入っていないときは、インストールしておこう。

$ curl -O https://bootstrap.pypa.io/get-pip.py | sudo python

ついでに言うと、システムにサードパーティ製のパッケージをどんどん突っ込んでいくと混乱の元になる。 なので virtualenv などを使って仮想環境を作れるようにした方が良い。 ここらへんに何を使うかは好みによる。

続いてバックエンドとして動作する TensorFlow (CPU 版) もインストールしておこう。

$ pip install tensorflow

インストールが終わったら python コマンドを実行して REPL を立ち上げよう。

$ python

REPL が起動したら Keras の API を使って MNIST データセットをダウンロードしよう。 これにはちょっと時間がかかる。

>>> from keras.datasets import mnist
Using TensorFlow backend.
>>> mnist.load_data()
Downloading data from https://s3.amazonaws.com/img-datasets/mnist.pkl.gz

ダウンロードが終わったら exit() 関数で REPL から抜けよう。

>>> exit()

続いて Keras のサンプルプログラムをダウンロードする。 これは Convolutional Neural Network (CNN) を使って MNIST データセットの画像を分類するものになっている。

$ curl -O https://raw.githubusercontent.com/fchollet/keras/master/examples/mnist_cnn.py

ちなみに、このプログラムにはリソースの開放関係でちょっとした問題があるっぽいのでパッチを当てておく。 これはそのうちいらなくなるかもしれないけど、余分にあったとしても特に問題はないコードなので。

$ echo 'K.clear_session()' >> mnist_cnn.py

CPU を使って学習させる

さて、上記でまずは Keras/TensorFlow を CPU を使って学習させる下準備が整った。

早速、先ほどダウンロードしてきたサンプルプログラムを実行してみよう。 time コマンドを先頭に入れることで学習にかかった時間を測ることにする。

$ time python mnist_cnn.py
Using TensorFlow backend.
X_train shape: (60000, 28, 28, 1)
60000 train samples
10000 test samples
Train on 60000 samples, validate on 10000 samples
Epoch 1/12
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.1 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX instructions, but these are available on your machine and could speed up CPU computations.
60000/60000 [==============================] - 78s - loss: 0.3744 - acc: 0.8854 - val_loss: 0.0887 - val_acc: 0.9716
Epoch 2/12
60000/60000 [==============================] - 78s - loss: 0.1357 - acc: 0.9595 - val_loss: 0.0634 - val_acc: 0.9800
Epoch 3/12
60000/60000 [==============================] - 79s - loss: 0.1033 - acc: 0.9693 - val_loss: 0.0572 - val_acc: 0.9818
Epoch 4/12
60000/60000 [==============================] - 78s - loss: 0.0890 - acc: 0.9740 - val_loss: 0.0448 - val_acc: 0.9853
Epoch 5/12
60000/60000 [==============================] - 76s - loss: 0.0784 - acc: 0.9767 - val_loss: 0.0423 - val_acc: 0.9865
Epoch 6/12
60000/60000 [==============================] - 80s - loss: 0.0710 - acc: 0.9788 - val_loss: 0.0384 - val_acc: 0.9871
Epoch 7/12
60000/60000 [==============================] - 84s - loss: 0.0656 - acc: 0.9808 - val_loss: 0.0364 - val_acc: 0.9891
Epoch 8/12
60000/60000 [==============================] - 84s - loss: 0.0606 - acc: 0.9818 - val_loss: 0.0354 - val_acc: 0.9888
Epoch 9/12
60000/60000 [==============================] - 79s - loss: 0.0555 - acc: 0.9839 - val_loss: 0.0336 - val_acc: 0.9883
Epoch 10/12
60000/60000 [==============================] - 76s - loss: 0.0545 - acc: 0.9837 - val_loss: 0.0325 - val_acc: 0.9895
Epoch 11/12
60000/60000 [==============================] - 76s - loss: 0.0500 - acc: 0.9856 - val_loss: 0.0325 - val_acc: 0.9900
Epoch 12/12
60000/60000 [==============================] - 78s - loss: 0.0476 - acc: 0.9860 - val_loss: 0.0328 - val_acc: 0.9893
Test score: 0.0328271338152
Test accuracy: 0.9893
python mnist_cnn.py  2470.93s user 664.19s system 326% cpu 16:00.66 total

MNIST を分類したときの正解率が 98.93% となっている。 かかった時間は 2471 秒で、ようするに 40 分もかかったことになる。

もし仮に一からニューラルネットワークを組むとしたら、もちろんこんな風に一発では上手くいかない。 たくさんのトライアンドエラーを繰り返して、その毎回にこんな時間がかかるとしたら気が遠くなってくる。

TensorFlow (GPU 版) をインストールする

CPU での学習にかかる時間に絶望したところで、次は代わりに GPU を使って学習させてみよう。 GPU は並列化しやすい単純な計算が得意で、ニューラルネットワークで内部的に用いられる行列演算もそれに当たる。

ニューラルネットワークがどうして行列演算を使うのかを知りたいときは、次の本を読むのがおすすめ。 この本はディープラーニングについて理論と実装の両面から理解を深められる貴重な本だと思う。

ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装

ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装

話が脱線したけど Keras/TensorFlow で組むニューラルネットワークを GPU で学習させるには CUDA が必要になる。 また、バックエンドとして動作する TensorFlow についても GPU 対応版のものをインストールする必要がある。 代わりに、サンプルコードに関しては、学習に CPU を使うときと GPU を使うときで同じものが使える。

まずは Homebrew を使って coreutils をインストールしておく。

$ brew install coreutils

もし Homebrew が入っていないという場合には公式サイトの記述をもとにインストールしておこう。

続いて CUDA をインストールする。 これも Homebrew Cask を使って入れられる。

$ brew tap caskroom/cask
$ brew cask install cuda

続いて cuDNN をインストールする。 これだけはコマンドラインでちょちょいとインストールすることができない。 なぜかというと、サイトでユーザ登録をしないとダウンロードできないようになっているから。

なので、まずは以下の公式サイトでユーザ登録をしよう。

NVIDIA cuDNN | NVIDIA Developer

その上で、今なら “cuDNN v5.1 Library for OSX” というバージョンをダウンロードしてくる。

ダウンロードできたら、それを CUDA のディレクトリに解凍して放り込む。

$ sudo tar -xvf cudnn-8.0-osx-x64-v5.1.tgz -C /usr/local

ここで念のためマシンを再起動しておこう。

$ sudo shutdown -r now

再起動が終わったら、まずは CPU 版の TensorFlow をアンインストールする。 そして、改めて GPU 版の TensorFlow をインストールしよう。 これには名前のサフィックスとして -gpu がついている。

$ pip uninstall -y tensorflow
$ pip install tensorflow-gpu

GPU を使って学習させる

さて、これで GPU を使ってニューラルネットワークを学習させるための準備が整った。 先ほどと同じサンプルコードを実行して GPU を使った学習を試してみよう。

$ time python mnist_cnn.py
Using TensorFlow backend.
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcublas.8.0.dylib locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcudnn.5.dylib locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcufft.8.0.dylib locally
I tensorflow/stream_executor/dso_loader.cc:126] Couldn't open CUDA library libcuda.1.dylib. LD_LIBRARY_PATH: :/usr/local/cuda/lib
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcuda.dylib locally
I tensorflow/stream_executor/dso_loader.cc:135] successfully opened CUDA library libcurand.8.0.dylib locally
X_train shape: (60000, 28, 28, 1)
60000 train samples
10000 test samples
Train on 60000 samples, validate on 10000 samples
Epoch 1/12
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.1 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use SSE4.2 instructions, but these are available on your machine and could speed up CPU computations.
W tensorflow/core/platform/cpu_feature_guard.cc:45] The TensorFlow library wasn't compiled to use AVX instructions, but these are available on your machine and could speed up CPU computations.
I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:874] OS X does not support NUMA - returning NUMA node zero
I tensorflow/core/common_runtime/gpu/gpu_device.cc:885] Found device 0 with properties:
name: GeForce GTX 675MX
major: 3 minor: 0 memoryClockRate (GHz) 0.719
pciBusID 0000:01:00.0
Total memory: 1023.69MiB
Free memory: 639.50MiB
I tensorflow/core/common_runtime/gpu/gpu_device.cc:906] DMA: 0
I tensorflow/core/common_runtime/gpu/gpu_device.cc:916] 0:   Y
I tensorflow/core/common_runtime/gpu/gpu_device.cc:975] Creating TensorFlow device (/gpu:0) -> (device: 0, name: GeForce GTX 675MX, pci bus id: 0000:01:00.0)
60000/60000 [==============================] - 16s - loss: 0.3717 - acc: 0.8865 - val_loss: 0.0904 - val_acc: 0.9714
Epoch 2/12
60000/60000 [==============================] - 13s - loss: 0.1365 - acc: 0.9596 - val_loss: 0.0616 - val_acc: 0.9809
Epoch 3/12
60000/60000 [==============================] - 13s - loss: 0.1048 - acc: 0.9694 - val_loss: 0.0515 - val_acc: 0.9841
Epoch 4/12
60000/60000 [==============================] - 13s - loss: 0.0894 - acc: 0.9736 - val_loss: 0.0439 - val_acc: 0.9857
Epoch 5/12
60000/60000 [==============================] - 13s - loss: 0.0769 - acc: 0.9771 - val_loss: 0.0401 - val_acc: 0.9869
Epoch 6/12
60000/60000 [==============================] - 13s - loss: 0.0700 - acc: 0.9790 - val_loss: 0.0360 - val_acc: 0.9885
Epoch 7/12
60000/60000 [==============================] - 13s - loss: 0.0646 - acc: 0.9806 - val_loss: 0.0364 - val_acc: 0.9881
Epoch 8/12
60000/60000 [==============================] - 13s - loss: 0.0611 - acc: 0.9819 - val_loss: 0.0345 - val_acc: 0.9882
Epoch 9/12
60000/60000 [==============================] - 13s - loss: 0.0553 - acc: 0.9838 - val_loss: 0.0323 - val_acc: 0.9889
Epoch 10/12
60000/60000 [==============================] - 13s - loss: 0.0534 - acc: 0.9841 - val_loss: 0.0307 - val_acc: 0.9898
Epoch 11/12
60000/60000 [==============================] - 13s - loss: 0.0504 - acc: 0.9851 - val_loss: 0.0297 - val_acc: 0.9895
Epoch 12/12
60000/60000 [==============================] - 13s - loss: 0.0474 - acc: 0.9860 - val_loss: 0.0298 - val_acc: 0.9899
Test score: 0.0298051692682
Test accuracy: 0.9899
python mnist_cnn.py  138.85s user 22.06s system 91% cpu 2:55.87 total

今度は 139 秒で学習が終わった。

CPU を使ったときが 2471 秒なので、およそ 17.8 倍も学習が高速化できた。 学習時間が約 5.6% に短縮されているともいえる。 これだけ違うとトライアンドエラーを繰り返すときにかかる時間は全く変わってくる。

まとめ

  • Keras/TensorFlow で学習に GPU を使いたいときは、以下のものを入れる
    • CUDA
    • cuDNN
    • GPU 対応版 TensorFlow
  • ソースコードは CPU/GPU どちらを学習に使うときであっても共用できる
  • 学習にかかる時間を CPU と GPU で比較してみた
    • 今回試したパターンでは GPU を使うと CPU 比で 17.8 倍も高速化できた
  • 結論: ディープラーニングをやるなら GPU は絶対に使おう

めでたしめでたし。

統計: Python と R で重回帰分析してみる

今回は R と Python の両方を使って重回帰分析をしてみる。 モチベーションとしては、できるだけ手に慣れた Python を使って分析をしていきたいという気持ちがある。 ただ、計算結果が意図通りのものになっているのかを R の結果と見比べて確かめておきたい。

また、分析にはボストンデータセットを用いる。 このデータセットはボストンの各地区ごとの不動産の平均価格と、それに付随する情報が入っている。

今回使った環境は次の通り。 R と Python は、あらかじめインストール済みであることを想定している。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G1212
$ python --version
Python 3.5.2
$ R --version
R version 3.3.2 (2016-10-31) -- "Sincere Pumpkin Patch"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under the terms of the
GNU General Public License versions 2 or 3.
For more information about these matters see
http://www.gnu.org/licenses/.

もし、まだ R をインストールしていないのであれば、こちらの記事が参考になるかも。

blog.amedama.jp

R

まずは R のやり方から。

最初に R を起動する。

$ R

続いてボストンデータセットの入っているライブラリ MASS を読み込む。

> library(MASS)

データセットの中身は次のようになっている。

> head(Boston)
     crim zn indus chas   nox    rm  age    dis rad tax ptratio  black lstat
1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21
  medv
1 24.0
2 21.6
3 34.7
4 33.4
5 36.2
6 28.7

各変数は次のような意味を持っている。 ただ、今回はボストンデータセットについて詳しく見ることが目的ではないので詳しくは立ち入らない。

1. CRIM: per capita crime rate by town
2. ZN: proportion of residential land zoned for lots over 25,000 sq.ft.
3. INDUS: proportion of non-retail business acres per town
4. CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
5. NOX: nitric oxides concentration (parts per 10 million)
6. RM: average number of rooms per dwelling
7. AGE: proportion of owner-occupied units built prior to 1940
8. DIS: weighted distances to five Boston employment centres
9. RAD: index of accessibility to radial highways
10. TAX: full-value property-tax rate per $10,000
11. PTRATIO: pupil-teacher ratio by town
12. B: 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
13. LSTAT: % lower status of the population
14. MEDV: Median value of owner-occupied homes in $1000's
(https://archive.ics.uci.edu/ml/datasets/Housing より)

今回はシンプルに線形回帰モデルを使うので lm() 関数を使う。 medv は目的変数である「その地区の不動産の平均価格」になっている。 次の処理では、目的関数をその他の全ての変数を使って重回帰している。

> lmfit = lm(medv ~ ., data=Boston)

重回帰した結果の詳細は summary() 関数で取得できる。

> summary(lmfit)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max
-15.595  -2.730  -0.518   1.777  26.199

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 **
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288
chas         2.687e+00  8.616e-01   3.118 0.001925 **
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 **
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

自由度調整済み決定係数が 0.7338 ということで、あまり上手く説明できていなさそうな感じ。

Python

続いて Python でのやり方について。 今回は statsmodels というライブラリを使うことにする。 ちなみに、線形回帰モデルは別のライブラリにも色々な実装がある。 ただ、フィッティングしたときに得られるパラメータが statsmodels は細かく取れて R に近いので、今回はこちらを使ってみた。

まずは必要なライブラリをインストールする。 今回は回帰に使わない scikit-learn を入れているのはボストンデータセットを読み込むため。

$ pip install statsmodels scikit-learn

最初に Python の REPL を起動しよう。

$ python

起動したら scikit-learn からボストンデータセットを読み込む。

>>> from sklearn import datasets
>>> dataset = datasets.load_boston()

このデータセットには次のようなメンバがある。

>>> dir(dataset)
['DESCR', 'data', 'feature_names', 'target']

target メンバが目的変数となる各地区ごとの不動産物件の平均価格となっている。

>>> dataset.target[:10]
array([ 24. ,  21.6,  34.7,  33.4,  36.2,  28.7,  22.9,  27.1,  16.5,  18.9])

そして説明変数が data に入っている。 ボストンを 506 地区に分割して 13 の変数を与えたデータになっている。

>>> dataset.data.shape
(506, 13)

変数名は feature_names から取得できる。

>>> dataset.feature_names
array(['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD',
       'TAX', 'PTRATIO', 'B', 'LSTAT'],
      dtype='<U7')

データセットを読み込めたところで statsmodels をインポートしよう。

>>> from statsmodels import api as sm

線形回帰モデルを使うときは OLS クラスを使う。 第一引数に目的変数を、第二変数に説明変数を渡す。

>>> model = sm.OLS(dataset.target, dataset.data)

準備ができたら fit() メソッドを使って分析する。

>>> result = model.fit()

分析結果からは summary() メソッドを使って R の summary() 関数に近い内容が取得できる。

>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.959
Model:                            OLS   Adj. R-squared:                  0.958
Method:                 Least Squares   F-statistic:                     891.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):               0.00
Time:                        22:59:04   Log-Likelihood:                -1523.8
No. Observations:                 506   AIC:                             3074.
Df Residuals:                     493   BIC:                             3129.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
x1            -0.0916      0.034     -2.675      0.008        -0.159    -0.024
x2             0.0487      0.014      3.379      0.001         0.020     0.077
x3            -0.0038      0.064     -0.059      0.953        -0.130     0.123
x4             2.8564      0.904      3.160      0.002         1.080     4.633
x5            -2.8808      3.359     -0.858      0.392        -9.481     3.720
x6             5.9252      0.309     19.168      0.000         5.318     6.533
x7            -0.0072      0.014     -0.523      0.601        -0.034     0.020
x8            -0.9680      0.196     -4.947      0.000        -1.352    -0.584
x9             0.1704      0.067      2.554      0.011         0.039     0.302
x10           -0.0094      0.004     -2.393      0.017        -0.017    -0.002
x11           -0.3924      0.110     -3.571      0.000        -0.608    -0.177
x12            0.0150      0.003      5.561      0.000         0.010     0.020
x13           -0.4170      0.051     -8.214      0.000        -0.517    -0.317
==============================================================================
Omnibus:                      204.050   Durbin-Watson:                   0.999
Prob(Omnibus):                  0.000   Jarque-Bera (JB):             1372.527
Skew:                           1.609   Prob(JB):                    9.11e-299
Kurtosis:                      10.399   Cond. No.                     8.50e+03
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.5e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

めでたしめでたし、と行きたいところだけどちょっと待ってほしい。 先ほど取得した R の結果と内容が異なっている。

実は R の分析では説明変数にバイアス項という変数が自動で挿入されていた。 これは、不動産物件のベースとなる価格と捉えることができる。

statsmodels を使っているときにバイアス項を説明変数に挿入するときは statsmodels.add_constant() 関数を使う。 バイアス項を挿入した説明変数を X という名前で保存しておこう。 ついでに目的変数も Y という変数で保存しておく。

>>> X = sm.add_constant(dataset.data)
>>> Y = dataset.target

バイアス項を追加した説明変数を用いて、先ほどと同じように回帰してみよう。

>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          6.95e-135
Time:                        23:00:34   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         36.4911      5.104      7.149      0.000        26.462    46.520
x1            -0.1072      0.033     -3.276      0.001        -0.171    -0.043
x2             0.0464      0.014      3.380      0.001         0.019     0.073
x3             0.0209      0.061      0.339      0.735        -0.100     0.142
x4             2.6886      0.862      3.120      0.002         0.996     4.381
x5           -17.7958      3.821     -4.658      0.000       -25.302   -10.289
x6             3.8048      0.418      9.102      0.000         2.983     4.626
x7             0.0008      0.013      0.057      0.955        -0.025     0.027
x8            -1.4758      0.199     -7.398      0.000        -1.868    -1.084
x9             0.3057      0.066      4.608      0.000         0.175     0.436
x10           -0.0123      0.004     -3.278      0.001        -0.020    -0.005
x11           -0.9535      0.131     -7.287      0.000        -1.211    -0.696
x12            0.0094      0.003      3.500      0.001         0.004     0.015
x13           -0.5255      0.051    -10.366      0.000        -0.625    -0.426
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                     1.51e+04
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.51e+04. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

今度は R で得られた内容と一致している!

AIC

重回帰分析では、モデリングに用いる説明変数の選択基準として AIC という数値を小さくするやり方があるので、その計算方法について。

R

R では AIC() 関数が用意されているので、それに回帰結果を渡す。

> AIC(lmfit)
[1] 3027.609

Python

Python では statsmodels の回帰結果のインスタンスに aic というメンバがあって、そこに格納されている。 ちなみに AIC はサマリーの中にも表示されていた。

>>> result.aic
3025.6767200074596

標準化

標準化というのはデータセットを平均が 0 で標準偏差が 1 になるように加工すること。 説明変数が標準化されていると、それぞれの説明変数が目的変数に対してどれくらい影響を与えているかを見たりするのに分かりやすい。

R

R には標準化のための関数として scale() が用意されている。

> X <- scale(Boston[, -14])
> Y <- Boston[, 14]

説明変数を標準化して回帰してみよう。

> lmfit <- lm(Y ~ X)
> summary(lmfit)

Call:
lm(formula = Y ~ X)

Residuals:
    Min      1Q  Median      3Q     Max
-15.595  -2.730  -0.518   1.777  26.199

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.53281    0.21095 106.814  < 2e-16 ***
Xcrim       -0.92906    0.28269  -3.287 0.001087 **
Xzn          1.08264    0.32016   3.382 0.000778 ***
Xindus       0.14104    0.42188   0.334 0.738288
Xchas        0.68241    0.21884   3.118 0.001925 **
Xnox        -2.05875    0.44262  -4.651 4.25e-06 ***
Xrm          2.67688    0.29364   9.116  < 2e-16 ***
Xage         0.01949    0.37184   0.052 0.958229
Xdis        -3.10712    0.41999  -7.398 6.01e-13 ***
Xrad         2.66485    0.57770   4.613 5.07e-06 ***
Xtax        -2.07884    0.63379  -3.280 0.001112 **
Xptratio    -2.06265    0.28323  -7.283 1.31e-12 ***
Xblack       0.85011    0.24521   3.467 0.000573 ***
Xlstat      -3.74733    0.36216 -10.347  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

説明変数が標準化されたことで、どの変数がどれくらい影響を与えているのかが分かりやすくなった。 例えば lstat などは影響が大きそうだ。

Python

Python の場合は scikit-learn の preprocessing() 関数を使うと簡単に標準化できる。

>>> X = dataset.data
>>> from sklearn import preprocessing
>>> X_scaled = preprocessing.scale(X)

標準化した上でバイアス項を追加して回帰してみよう。

>>> X_scaled_with_constant = sm.add_constant(X_scaled)
>>> model = sm.OLS(Y, X_scaled_with_constant)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.741
Model:                            OLS   Adj. R-squared:                  0.734
Method:                 Least Squares   F-statistic:                     108.1
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          6.95e-135
Time:                        23:50:36   Log-Likelihood:                -1498.8
No. Observations:                 506   AIC:                             3026.
Df Residuals:                     492   BIC:                             3085.
Df Model:                          13
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const         22.5328      0.211    106.807      0.000        22.118    22.947
x1            -0.9204      0.281     -3.276      0.001        -1.472    -0.368
x2             1.0810      0.320      3.380      0.001         0.453     1.709
x3             0.1430      0.421      0.339      0.735        -0.685     0.971
x4             0.6822      0.219      3.120      0.002         0.253     1.112
x5            -2.0601      0.442     -4.658      0.000        -2.929    -1.191
x6             2.6706      0.293      9.102      0.000         2.094     3.247
x7             0.0211      0.371      0.057      0.955        -0.709     0.751
x8            -3.1044      0.420     -7.398      0.000        -3.929    -2.280
x9             2.6588      0.577      4.608      0.000         1.525     3.792
x10           -2.0759      0.633     -3.278      0.001        -3.320    -0.832
x11           -2.0622      0.283     -7.287      0.000        -2.618    -1.506
x12            0.8566      0.245      3.500      0.001         0.376     1.338
x13           -3.7487      0.362    -10.366      0.000        -4.459    -3.038
==============================================================================
Omnibus:                      178.029   Durbin-Watson:                   1.078
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              782.015
Skew:                           1.521   Prob(JB):                    1.54e-170
Kurtosis:                       8.276   Cond. No.                         9.82
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
"""

結果が R の内容と一致していることが分かる。

ちなみに標準化の操作自体は単純なので、自前でやることもできる。 数式に直すと、こんな感じ。 各要素から平均 ( \mu) を引いて標準偏差 ( \sigma) で割る。

 Z = \frac{X - \mu}{\sigma}

コードにするとこんな感じ。

>>> X_scaled = (X - X.mean(axis=0)) / X.std(axis=0)

交互作用モデル

交互作用モデルというのは、複数の変数をまとめて説明変数として扱うやり方。 ある説明変数と別の説明変数の相乗作用があるときに使うと良い。

R

R では lm() 関数を実行するときに説明変数を ^2 するだけで交互作用モデルを使える。

> lmfit = lm(medv ~ .^2, data=Boston)
> summary(lmfit)

Call:
lm(formula = medv ~ .^2, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max
-7.9374 -1.5344 -0.1068  1.2973 17.8500

Coefficients:
                Estimate Std. Error t value Pr(>|t|)
(Intercept)   -1.579e+02  6.800e+01  -2.323 0.020683 *
crim          -1.707e+01  6.554e+00  -2.605 0.009526 **
zn            -7.529e-02  4.580e-01  -0.164 0.869508
indus         -2.819e+00  1.696e+00  -1.663 0.097111 .
chas           4.451e+01  1.952e+01   2.280 0.023123 *
nox            2.006e+01  7.516e+01   0.267 0.789717
rm             2.527e+01  5.699e+00   4.435 1.18e-05 ***
age            1.263e+00  2.728e-01   4.630 4.90e-06 ***
dis           -1.698e+00  4.604e+00  -0.369 0.712395
rad            1.861e+00  2.464e+00   0.755 0.450532
tax            3.670e-02  1.440e-01   0.255 0.798978
ptratio        2.725e+00  2.850e+00   0.956 0.339567
black          9.942e-02  7.468e-02   1.331 0.183833
lstat          1.656e+00  8.533e-01   1.940 0.053032 .
crim:zn        4.144e-01  1.804e-01   2.297 0.022128 *
crim:indus    -4.693e-02  4.480e-01  -0.105 0.916621
crim:chas      2.428e+00  5.710e-01   4.251 2.63e-05 ***
crim:nox      -1.108e+00  9.285e-01  -1.193 0.233425
crim:rm        2.163e-01  4.907e-02   4.409 1.33e-05 ***
crim:age      -3.083e-03  3.781e-03  -0.815 0.415315
crim:dis      -1.903e-01  1.060e-01  -1.795 0.073307 .
crim:rad      -6.584e-01  5.815e-01  -1.132 0.258198
crim:tax       3.479e-02  4.287e-02   0.812 0.417453
crim:ptratio   4.915e-01  3.328e-01   1.477 0.140476
crim:black    -4.612e-04  1.793e-04  -2.572 0.010451 *
crim:lstat     2.964e-02  6.544e-03   4.530 7.72e-06 ***
...(省略)...
rad:tax        3.131e-05  1.446e-03   0.022 0.982729
rad:ptratio   -4.379e-02  8.392e-02  -0.522 0.602121
rad:black     -4.362e-04  2.518e-03  -0.173 0.862561
rad:lstat     -2.529e-02  1.816e-02  -1.392 0.164530
tax:ptratio    7.854e-03  2.504e-03   3.137 0.001830 **
tax:black     -4.785e-07  1.999e-04  -0.002 0.998091
tax:lstat     -1.403e-03  1.208e-03  -1.162 0.245940
ptratio:black  1.203e-03  3.361e-03   0.358 0.720508
ptratio:lstat  3.901e-03  2.985e-02   0.131 0.896068
black:lstat   -6.118e-04  4.157e-04  -1.472 0.141837
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.852 on 414 degrees of freedom
Multiple R-squared:  0.9212, Adjusted R-squared:  0.9039
F-statistic: 53.18 on 91 and 414 DF,  p-value: < 2.2e-16

結果から説明変数の数が一気に増えてることが分かる。 これは 13 種類の説明変数から 2 種類を組み合わせで取り出して、それらを交互作用モデルの説明変数として追加しているため。

また、自由度調整済み決定係数が 0.9039 となっており、ただ 13 種類の変数を使うときよりも上手くモデリングできていることが分かる。 まあ、この値だけを見ていても過適合の恐れがあるけど。

Python

statsmodels には回帰するときのやり方を R のような書き方 (R-style らしい) で指定できる。 このやり方を使うとさくっと相互作用モデルも使えそうな感じ。

Fitting models using R-style formulas — statsmodels 0.7.0 documentation

ただまあ、これを使うよりも自前で用意する方が楽しそう。 なので、今回は自分で計算してみることにする。

Python で組み合わせを計算するときは itertools モジュールを使うと良い。 この中にある combinations() 関数を使えば簡単に組み合わせが取れる。 例えば説明変数の組み合わせを計算したいなら、こう。

>>> import itertools
>>> list(itertools.combinations(dataset.feature_names, 2))
[('CRIM', 'ZN'), ('CRIM', 'INDUS'), ('CRIM', 'CHAS'), ('CRIM', 'NOX'), ('CRIM', 'RM'), ('CRIM', 'AGE'), ('CRIM', 'DIS'), ('CRIM', 'RAD'), ('CRIM', 'TAX'), ('CRIM', 'PTRATIO'), ('CRIM', 'B'), ('CRIM', 'LSTAT'), ('ZN', 'INDUS'), ('ZN', 'CHAS'), ('ZN', 'NOX'), ('ZN', 'RM'), ('ZN', 'AGE'), ('ZN', 'DIS'), ('ZN', 'RAD'), ('ZN', 'TAX'), ('ZN', 'PTRATIO'), ('ZN', 'B'), ('ZN', 'LSTAT'), ('INDUS', 'CHAS'), ('INDUS', 'NOX'), ('INDUS', 'RM'), ('INDUS', 'AGE'), ('INDUS', 'DIS'), ('INDUS', 'RAD'), ('INDUS', 'TAX'), ('INDUS', 'PTRATIO'), ('INDUS', 'B'), ('INDUS', 'LSTAT'), ('CHAS', 'NOX'), ('CHAS', 'RM'), ('CHAS', 'AGE'), ('CHAS', 'DIS'), ('CHAS', 'RAD'), ('CHAS', 'TAX'), ('CHAS', 'PTRATIO'), ('CHAS', 'B'), ('CHAS', 'LSTAT'), ('NOX', 'RM'), ('NOX', 'AGE'), ('NOX', 'DIS'), ('NOX', 'RAD'), ('NOX', 'TAX'), ('NOX', 'PTRATIO'), ('NOX', 'B'), ('NOX', 'LSTAT'), ('RM', 'AGE'), ('RM', 'DIS'), ('RM', 'RAD'), ('RM', 'TAX'), ('RM', 'PTRATIO'), ('RM', 'B'), ('RM', 'LSTAT'), ('AGE', 'DIS'), ('AGE', 'RAD'), ('AGE', 'TAX'), ('AGE', 'PTRATIO'), ('AGE', 'B'), ('AGE', 'LSTAT'), ('DIS', 'RAD'), ('DIS', 'TAX'), ('DIS', 'PTRATIO'), ('DIS', 'B'), ('DIS', 'LSTAT'), ('RAD', 'TAX'), ('RAD', 'PTRATIO'), ('RAD', 'B'), ('RAD', 'LSTAT'), ('TAX', 'PTRATIO'), ('TAX', 'B'), ('TAX', 'LSTAT'), ('PTRATIO', 'B'), ('PTRATIO', 'LSTAT'), ('B', 'LSTAT')]

組み合わせは  {}_{13} C_2 = \frac{13!}{2!(13 - 2)!} = 78 通り。

>>> len(list(itertools.combinations(dataset.feature_names, 2)))
78

それぞれの説明変数を取り出すための添字にするとこんな感じ。

>>> list(itertools.combinations(range(dataset.data.shape[1]), 2))
[(0, 1), (0, 2), (0, 3), (0, 4), (0, 5), (0, 6), (0, 7), (0, 8), (0, 9), (0, 10), (0, 11), (0, 12), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 8), (1, 9), (1, 10), (1, 11), (1, 12), (2, 3), (2, 4), (2, 5), (2, 6), (2, 7), (2, 8), (2, 9), (2, 10), (2, 11), (2, 12), (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (3, 10), (3, 11), (3, 12), (4, 5), (4, 6), (4, 7), (4, 8), (4, 9), (4, 10), (4, 11), (4, 12), (5, 6), (5, 7), (5, 8), (5, 9), (5, 10), (5, 11), (5, 12), (6, 7), (6, 8), (6, 9), (6, 10), (6, 11), (6, 12), (7, 8), (7, 9), (7, 10), (7, 11), (7, 12), (8, 9), (8, 10), (8, 11), (8, 12), (9, 10), (9, 11), (9, 12), (10, 11), (10, 12), (11, 12)]

交互作用モデルで使うそれぞれの説明変数は、説明変数同士を乗算することで計算できる。

>>> b = np.array([X[:, 0] * X[:, 1]])

計算した説明変数を、既存の説明変数にカラムとして追加するには次のようにすれば良い。

>>> a = X
>>> np.concatenate((a, b.T), axis=1)
array([[  6.32000000e-03,   1.80000000e+01,   2.31000000e+00, ...,
          3.96900000e+02,   4.98000000e+00,   1.13760000e-01],
       [  2.73100000e-02,   0.00000000e+00,   7.07000000e+00, ...,
          3.96900000e+02,   9.14000000e+00,   0.00000000e+00],
       [  2.72900000e-02,   0.00000000e+00,   7.07000000e+00, ...,
          3.92830000e+02,   4.03000000e+00,   0.00000000e+00],
       ...,
       [  6.07600000e-02,   0.00000000e+00,   1.19300000e+01, ...,
          3.96900000e+02,   5.64000000e+00,   0.00000000e+00],
       [  1.09590000e-01,   0.00000000e+00,   1.19300000e+01, ...,
          3.93450000e+02,   6.48000000e+00,   0.00000000e+00],
       [  4.74100000e-02,   0.00000000e+00,   1.19300000e+01, ...,
          3.96900000e+02,   7.88000000e+00,   0.00000000e+00]])

これで説明変数が一つ増える。

>>> np.concatenate((a, b.T), axis=1).shape
(506, 14)

上記にもとづいて相互作用モデルで使う全ての変数を揃える。

>>> X = dataset.data
>>> for i, j in itertools.combinations(range(X.shape[1]), 2):
...     interaction_column = np.array([X[:, i] * X[:, j]])
...     X = np.concatenate((X, interaction_column.T), axis=1)
...
>>> X.shape
(506, 91)

そして、バイアス項を追加する。

>>> X = sm.add_constant(X)
>>> X.shape
(506, 92)

上記の説明変数を使って回帰する。

>>> model = sm.OLS(Y, X)
>>> result = model.fit()
>>> result.summary()
<class 'statsmodels.iolib.summary.Summary'>
"""
                            OLS Regression Results
==============================================================================
Dep. Variable:                      y   R-squared:                       0.921
Model:                            OLS   Adj. R-squared:                  0.904
Method:                 Least Squares   F-statistic:                     53.21
Date:                Thu, 22 Dec 2016   Prob (F-statistic):          5.85e-181
Time:                        19:16:55   Log-Likelihood:                -1197.3
No. Observations:                 506   AIC:                             2579.
Df Residuals:                     414   BIC:                             2967.
Df Model:                          91
Covariance Type:            nonrobust
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
const       -159.4473     67.980     -2.346      0.019      -293.076   -25.819
x1           -16.9621      6.549     -2.590      0.010       -29.835    -4.089
x2            -0.0769      0.458     -0.168      0.867        -0.977     0.823
x3            -2.8163      1.695     -1.662      0.097        -6.148     0.516
x4            44.6397     19.519      2.287      0.023         6.271    83.008
x5            22.1561     75.150      0.295      0.768      -125.567   169.879
x6            25.4191      5.700      4.459      0.000        14.214    36.624
x7             1.2569      0.273      4.612      0.000         0.721     1.793
x8            -1.6398      4.604     -0.356      0.722       -10.690     7.410
x9             1.8056      2.462      0.733      0.464        -3.034     6.646
x10            0.0374      0.144      0.260      0.795        -0.246     0.320
x11            2.7520      2.849      0.966      0.335        -2.849     8.353
x12            0.1008      0.074      1.354      0.177        -0.046     0.247
x13            1.6587      0.853      1.944      0.053        -0.018     3.336
x14            0.4124      0.180      2.287      0.023         0.058     0.767
x15           -0.0467      0.448     -0.104      0.917        -0.927     0.834
x16            2.4377      0.571      4.271      0.000         1.316     3.560
x17           -1.2206      0.920     -1.327      0.185        -3.029     0.588
x18            0.2117      0.049      4.345      0.000         0.116     0.307
x19           -0.0031      0.004     -0.828      0.408        -0.011     0.004
x20           -0.1956      0.106     -1.848      0.065        -0.404     0.012
x21           -0.6599      0.581     -1.135      0.257        -1.803     0.483
x22            0.0349      0.043      0.815      0.416        -0.049     0.119
x23            0.4895      0.333      1.471      0.142        -0.164     1.143
x24           -0.0004      0.000     -2.513      0.012        -0.001 -9.71e-05
x25            0.0295      0.007      4.512      0.000         0.017     0.042
...(省略)...
x82         3.151e-05      0.001      0.022      0.983        -0.003     0.003
x83           -0.0439      0.084     -0.523      0.601        -0.209     0.121
x84           -0.0004      0.003     -0.165      0.869        -0.005     0.005
x85           -0.0251      0.018     -1.384      0.167        -0.061     0.011
x86            0.0078      0.003      3.135      0.002         0.003     0.013
x87        -2.303e-06      0.000     -0.012      0.991        -0.000     0.000
x88           -0.0014      0.001     -1.167      0.244        -0.004     0.001
x89            0.0012      0.003      0.349      0.727        -0.005     0.008
x90            0.0038      0.030      0.128      0.898        -0.055     0.062
x91           -0.0006      0.000     -1.510      0.132        -0.001     0.000
==============================================================================
Omnibus:                      121.344   Durbin-Watson:                   1.524
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              628.079
Skew:                           0.942   Prob(JB):                    4.11e-137
Kurtosis:                       8.123   Cond. No.                     1.18e+08
==============================================================================

Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 1.18e+08. This might indicate that there are
strong multicollinearity or other numerical problems.
"""

上記を見ると R で相互作用モデルを使って計算したときの内容と一致していることが分かる。

まとめ

今回は R と Python (statsmodels) で線形モデルを使った重回帰分析をしてみた。

統計: 二つのグループの平均と分散を合成する

例えば、あるグループ A と B が別々にテストを受けたとする。 それぞれのグループの人数と平均点、そして分散は分かっているとしよう。 このとき、グループ A と B を合わせた全体での平均や分散は計算することができるだろうか?

結論から言うと、これはできる。 今回は、その手順について書いてみることにする。 また、同時に Python を使ってサンプルデータを生成したり計算の裏付けをしてみよう。

使った環境は次の通り。

$ sw_vers 
ProductName:    Mac OS X
ProductVersion: 10.11.6
BuildVersion:   15G1212
$ python --version
Python 3.5.2

題材とするデータを生成する

平均や分散の計算を楽にするために、まずは NumPy をインストールしておこう。

$ pip install numpy

そして Python の REPL を起動する。

$ python

起動したら NumPy をインポートしよう。

>>> import numpy as np

そして、まずはグループ A を生成する。 パラメータは適当だけど、点数は 40 から 90 の範囲で 10 人分作っておく。

>>> A = np.random.randint(40, 90, 10)
>>> A
array([57, 67, 74, 79, 82, 53, 80, 46, 74, 47])

上記の平均と分散はそれぞれ 65.9176.09 になった。

>>> A.mean()
65.900000000000006
>>> A.var()
176.09

同じようにグループ B も生成する。 こちらは平均が 72.466... で分散が 60.9155... になった。

>>> B = np.random.randint(55, 85, 15)
>>> B
array([83, 77, 77, 59, 67, 57, 77, 71, 68, 77, 72, 76, 84, 64, 78])
>>> B.mean()
72.466666666666669
>>> B.var()
60.915555555555549

そして、上記を合わせたグループについても計算しておこう。 便宜的に、ここではこれをグループ C と呼ぶことにする。

>>> C = np.r_[A, B]
>>> C
array([57, 67, 74, 79, 82, 53, 80, 46, 74, 47, 83, 77, 77, 59, 67, 57, 77,
       71, 68, 77, 72, 76, 84, 64, 78])

グループ C の平均と分散は次のようになった。 今回は、この値をグループ A と B の平均と分散とデータ点数からだけで求めるのが最終的な目的となる。

>>> C.mean()
69.840000000000003
>>> C.var()
117.3344

平均を合成する

平均の合成はすごく簡単。 まず、それぞれのグループの平均値に人数をかけると各グループの点数の合計が得られる。 あとは、それをまとめてから全体の人数で割れば合成した平均が得られる。

数式にすると、こんな感じ。

 \bar{x_C} = \frac { \bar{x_A} \times n_A + \bar{x_B} + n_B }{ n_A + n_B }

ここで  \bar{x_A} \bar{x_B} は各グループの平均値を表している。 そして  n_A n_B は各グループの人数を表している。

実際に Python を使って計算してみよう。

>>> C_average = (A.mean() * 10 + B.mean() * 15) / (10 + 15)
>>> C_average
69.840000000000003

上記の値は、実際にグループをまとめて計算した平均値と一致している。

>>> C.mean()
69.840000000000003

分散を合成する

分散の合成については、もうちょっと厄介な計算が必要になる。 また、分散の計算方法についても詳しく知っておかないといけない。

分散について

まず、最も基本的な分散の計算式は次のようになる。 分散というのは、各要素から平均値を引いたもの (偏差) を二乗して足し合わせて、それをデータ点数で割ることで得られる。

 s^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2}

上記で  \bar{x} x の平均値となっているので、計算式で書くとこう。

 \bar{x} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i

また、分散の公式は次のように式変形できる。 この式は、ようするに各要素を二乗して足し合わせてデータ点数で割ったものから、平均値の二乗を引いて計算することができる、ということ。

 s^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} (x_i - \bar{x})^{2} = \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i^{2} - \bar{x}^{2}

また、上記の「二乗したものを足し合わせてデータ点数で割る」という処理に着目してみよう。 これは言い換えると二乗したものの平均を計算していることになる。 バー記号は平均を表現する。

 \frac{1}{n} \displaystyle \sum_{i = 1}^{n} x_i^{2} = \bar{x^{2}}

つまり、分散は次のような式でも表現できることが分かった。 ようするに、分散は各要素を二乗した平均値 ( \bar{x^{2}})から各要素の平均値の二乗 ( \bar{x}^2) を引くことで計算できる。

 s^{2} = \bar{x^{2}} - \bar{x}^2

分散を合成する手順

さて、ここまでの説明で分散の計算方法は分かった。 改めて、求めなければいけないものの式を書いてみよう。

 s_C^{2} = \bar{x_C^{2}} - \bar{x_C}^2

上記の中で  \bar{x_C} については、既に平均の合成をしたところで計算が済んでいる。 問題は  \bar{x_C^{2}} の方で、これはまだ未知の値なので計算しなきゃいけない。

しかし  \bar{x_C^{2}} は、どうすれば求まるだろうか? 実は、これも先ほどの分散の公式から導くことができる。

 s^{2} = \bar{x^{2}} - \bar{x}^2

上記の中で、各グループ毎の  s^{2} \bar{x} は既に分かっている。 つまり、式を移行して  \bar{x^{2}} を求めるようにできるということ。

 \bar{x^{2}} = s^{2} + \bar{x}^2

次のように、各グループごとに二乗した平均値を求めることができる。

 \bar{x_A^{2}} = s_A^{2} + \bar{x_A}^2

 \bar{x_B^{2}} = s_B^{2} + \bar{x_B}^2

各グループごとの二乗した平均値を求めることができれば、あとは簡単で最初にやった平均の合成をすれば良い。

 \bar{x_C^{2}} = \frac{ \bar{x_A^{2}} \times n_A + \bar{x_B^{2}} \times n_B }{ n_A + n_B }

上記が計算できれば合成した分散を計算するのに必要な要素が揃うので、あとは定義通りにやるだけ。

 s_C^{2} = \bar{x_C^{2}} - \bar{x_C}^2

確かめてみる

それでは、上記の式を Python で実行して確かめてみよう。

まずはグループ A の二乗した平均値から求める。

>>> A_square_average = A.var() + A.mean() ** 2
>>> A_square_average
4518.9000000000005

上記は、次の式に対応している。

 \bar{x_A^{2}} = s_A^{2} + \bar{x_A}^2

そして、同じようにグループ B の二乗した平均値を求める。

>>> B_square_average = B.var() + B.mean() ** 2
>>> B_square_average
5312.333333333333

対応する式は次の通り。

 \bar{x_B^{2}} = s_B^{2} + \bar{x_B}^2

そして、グループをまとめたときの二乗した平均値を求める。

>>> C_square_average = (A_square_average * 10 + B_square_average * 15) / (10 + 15)
>>> C_square_average
4994.96

対応する式はこちら。

 \bar{x_C^{2}} = \frac{ \bar{x_A^{2}} \times n_A + \bar{x_B^{2}} \times n_B }{ n_A + n_B }

上記で必要な材料は揃ったので、あとは定義通りに計算するだけ。

>>> C_var = C_square_average - C_average ** 2
>>> C_var
117.33439999999973

実際に計算した値と比較してみよう。

>>> C.var()
117.3344
>>> C_var
117.33439999999973

微妙に誤差が出てるけど、バッチリ合ってるね!

まとめ

  • 二つのグループの平均と分散とデータ点数だけが分かっている状況であってもグループ全体の平均と分散は計算できる