Raspberry Pi + SORACOM Air セットアップ

　前からRaspberry PiやArduino等に興味はあったものの自分では試せていなかったのですが、先日 IoT Technology Conference if-up 2017 で 3G SIM の USBドングルをいただいたので、これを機に自分でもRaspberry Piを購入して色々と試してみることにしました。同じカンファレンスで紹介されていてとても良さそうだったIoTエンジニア養成読本も買って、まずはそのハンズオンの内容をベースに動くものを作ってみようと思います。

gihyo.jp

Raspberry Piのセットアップ

　まずは下記URLからRaspbian OSのイメージをダウンロードします。HDMIモニタやUSBキーボードを接続せずにセットアップしたかったので、デフォルトでsshサーバが起動する2016-09-28のイメージを選択しました。これより後のバージョンだとデフォルトではsshサーバは起動しないようです。（結局LANケーブルが見つからずHDMIモニタとUSBキーボードを繋いでセットアップしたので、最新のを選択してもよかったのですが。。）

http://ftp.jaist.ac.jp/pub/raspberrypi/raspbian_lite/images/raspbian_lite-2016-09-28/

　ダウンロードしたファイルは解凍してimgファイルを取り出しておきます。

　続いてイメージファイルをSDカードに書き込んでいきます。今回は16GBのマイクロSDにアダプタをつけてMacBook ProのSDカードスロットにさします。

　diskutilコマンドでストレージを確認します。

pi  $ diskutil list
/dev/disk0 (internal, physical):
   #:                       TYPE NAME                    SIZE       IDENTIFIER
   0:      GUID_partition_scheme                        *500.3 GB   disk0
   1:                        EFI EFI                     209.7 MB   disk0s1
   2:          Apple_CoreStorage Macintosh HD            499.4 GB   disk0s2
   3:                 Apple_Boot Recovery HD             650.0 MB   disk0s3

/dev/disk1 (internal, virtual):
   #:                       TYPE NAME                    SIZE       IDENTIFIER
   0:                            Macintosh HD           +499.1 GB   disk1
                                 Logical Volume on disk0s2
                                 D522C0D0-F775-4496-8BDA-640948662DCD
                                 Unlocked Encrypted

/dev/disk2 (internal, physical):
   #:                       TYPE NAME                    SIZE       IDENTIFIER
   0:     FDisk_partition_scheme                        *15.5 GB    disk2
   1:             Windows_FAT_32 NO NAME                 15.5 GB    disk2s1

　サイズから判断して /dev/disk2 がSDカードなので、ddコマンドでSDカードにOSのイメージを書き込みます。マウントしていると書き込めないので、diskutil unmountDisk コマンドでアンマウントしてから実行します。

pi  $ sudo dd if=2016-09-23-raspbian-jessie-lite.img of=/dev/rdisk2 bs=1m
Password:
dd: /dev/rdisk2: Resource busy
pi  $ 
pi  $ diskutil unmountDisk /dev/disk2
Unmount of all volumes on disk2 was successful
pi  $ 
pi  $ sudo dd if=2016-09-23-raspbian-jessie-lite.img of=/dev/rdisk2 bs=1m
1325+0 records in
1325+0 records out
1389363200 bytes transferred in 94.477593 secs (14705743 bytes/sec)
pi  $ 

　MacBook ProからSDカードを取り出してアダプタを外し、Raspberry PiのSDカードスロットにさします。今回買ったのはRaspberry Pi 3 Model Bです。

　HDMIモニタとUSBキーボードを接続して最後に電源ケーブルを接続すると、OSが起動します。

　ログインプロンプトが表示されたらデフォルトのログインIDとパスワードでログインし、無線LANへの接続の設定を行います。wpa_passphraseコマンドを利用するのですが、wpa_supplicant.confはrootユーザしか書き込み権限を持っていないので、書籍でも紹介されているようにパイプでつないで tee コマンドをsudoで使って書き込むのが良いのですが、手元のキーボードだとパイプが入力できず、キーボードの設定を変更すれば良いと思うのですが面倒だったので、一時的にroot以外にも書き込み権限を与えてリダイレクトでファイルに追記し、そのあとでまた権限を戻しました。MY_AP_SSID と MY_AP_PASSWORD は接続するAPのSSIDとパスワードに置き換えてください。設定後はOSを再起動します。

$ sudo chmod 606 /etc/wpa_supplicant/wpa_supplicant.conf
$ wpa_passphrase MY_AP_SSID MY_AP_PASSWORD >> /etc/wpa_supplicant/wpa_supplicant.conf
$ cat /etc/wpa_supplicant/wpa_supplicant.conf
$ sudo chmod 600 /etc/wpa_supplicant/wpa_supplicant.conf
$ sudo reboot

　再起動後にipコマンドで無線LANへ接続できていてIPアドレスが割り当てられていることを確認します。

pi@raspberrypi:~ $ ip a show dev wlan0
3: wlan0: <BROADCAST,MULTICAST,UP,LOWER_UP> mtu 1500 qdisc pfifo_fast state UP group default qlen 1000
    link/ether b8:27:eb:e6:89:f8 brd ff:ff:ff:ff:ff:ff
    inet 192.168.10.10/24 brd 192.168.10.255 scope global wlan0
       valid_lft forever preferred_lft forever
    inet6 2408:212:2862:5c00:b342:9793:7abe:c897/64 scope global noprefixroute dynamic 
       valid_lft 2591786sec preferred_lft 604586sec
    inet6 fe80::bc23:73fb:4394:7972/64 scope link 
       valid_lft forever preferred_lft forever

　同じLAN内のRaspberry Pi端末には raspberrypi.local でアクセスできるので、下記のようにMacBook Proからsshでログインします。

pi  $ ssh pi@raspberrypi.local
The authenticity of host 'raspberrypi.local (2408:212:2862:5c00:b342:9793:7abe:c897)' can't be established.
ECDSA key fingerprint is SHA256:B98CBwQsKcPnKMeIGBNQ065GNnMXZvBm1pKoHc7+0Zw.
Are you sure you want to continue connecting (yes/no)? yes
Warning: Permanently added 'raspberrypi.local,2408:212:2862:5c00:b342:9793:7abe:c897' (ECDSA) to the list of known hosts.
pi@raspberrypi.local's password: 

The programs included with the Debian GNU/Linux system are free software;
the exact distribution terms for each program are described in the
individual files in /usr/share/doc/*/copyright.

Debian GNU/Linux comes with ABSOLUTELY NO WARRANTY, to the extent
permitted by applicable law.
Last login: Sat Apr 29 16:26:21 2017

　下記コマンドでRaspbianを最新の状態にアップデートします。

% sudo apt-get update
% sudo apt-get upgrade
% sudo apt-get dist-upgrade

　最後にpasswdコマンドでデフォルトのログインパスワードをオリジナルのものに変更して、基本的なセットアップは終了です。

SORACOM Air での接続セットアップ

　カンファレンスでいただいた3G Sim の USBドングルを使って、SORACOM Airでネットワーク接続できるようにセットアップします。事前にSORACOM AirのSimを購入して、SORACOMのユーザアカウントの作成とSimの登録を済ませておきます。

　いただいたUSBドングルはABIT AK-020です。

　3G接続に必要なパッケージをインストールします。

pi@raspberrypi:~ $ sudo apt-get install -y usb-modeswitch wvdial
Reading package lists... Done
Building dependency tree       
Reading state information... Done
usb-modeswitch is already the newest version.
The following extra packages will be installed:
  libpcap0.8 libuniconf4.6 libwvstreams4.6-base libwvstreams4.6-extras ppp
The following NEW packages will be installed:
  libpcap0.8 libuniconf4.6 libwvstreams4.6-base libwvstreams4.6-extras ppp wvdial
0 upgraded, 6 newly installed, 0 to remove and 0 not upgraded.
Need to get 1390 kB of archives.
After this operation, 3127 kB of additional disk space will be used.
Get:1 http://mirrordirector.raspbian.org/raspbian/ jessie/main libpcap0.8 armhf 1.6.2-2 [121 kB]
Get:2 http://mirrordirector.raspbian.org/raspbian/ jessie/main libwvstreams4.6-base armhf 4.6.1-7 [235 kB]
Get:3 http://mirrordirector.raspbian.org/raspbian/ jessie/main libwvstreams4.6-extras armhf 4.6.1-7 [448 kB]
Get:4 http://mirrordirector.raspbian.org/raspbian/ jessie/main libuniconf4.6 armhf 4.6.1-7 [173 kB]
Get:5 http://mirrordirector.raspbian.org/raspbian/ jessie/main ppp armhf 2.4.6-3.1 [306 kB]
Get:6 http://mirrordirector.raspbian.org/raspbian/ jessie/main wvdial armhf 1.61-4.1 [107 kB]
Fetched 1390 kB in 2s (493 kB/s)
Can't set locale; make sure $LC_* and $LANG are correct!
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LC_CTYPE = "UTF-8",
        LANG = "en_GB.UTF-8"
    are supported and installed on your system.
perl: warning: Falling back to a fallback locale ("en_GB.UTF-8").
locale: Cannot set LC_CTYPE to default locale: No such file or directory
locale: Cannot set LC_ALL to default locale: No such file or directory
Preconfiguring packages ...
Selecting previously unselected package libpcap0.8:armhf.
(Reading database ... 31414 files and directories currently installed.)
Preparing to unpack .../libpcap0.8_1.6.2-2_armhf.deb ...
Unpacking libpcap0.8:armhf (1.6.2-2) ...
Selecting previously unselected package libwvstreams4.6-base.
Preparing to unpack .../libwvstreams4.6-base_4.6.1-7_armhf.deb ...
Unpacking libwvstreams4.6-base (4.6.1-7) ...
Selecting previously unselected package libwvstreams4.6-extras.
Preparing to unpack .../libwvstreams4.6-extras_4.6.1-7_armhf.deb ...
Unpacking libwvstreams4.6-extras (4.6.1-7) ...
Selecting previously unselected package libuniconf4.6.
Preparing to unpack .../libuniconf4.6_4.6.1-7_armhf.deb ...
Unpacking libuniconf4.6 (4.6.1-7) ...
Selecting previously unselected package ppp.
Preparing to unpack .../ppp_2.4.6-3.1_armhf.deb ...
Unpacking ppp (2.4.6-3.1) ...
Selecting previously unselected package wvdial.
Preparing to unpack .../wvdial_1.61-4.1_armhf.deb ...
Unpacking wvdial (1.61-4.1) ...
Processing triggers for man-db (2.7.0.2-5) ...
Processing triggers for systemd (215-17+deb8u6) ...
Setting up libpcap0.8:armhf (1.6.2-2) ...
Setting up libwvstreams4.6-base (4.6.1-7) ...
Setting up libwvstreams4.6-extras (4.6.1-7) ...
Setting up libuniconf4.6 (4.6.1-7) ...
Setting up ppp (2.4.6-3.1) ...
Setting up wvdial (1.61-4.1) ...
locale: Cannot set LC_CTYPE to default locale: No such file or directory
locale: Cannot set LC_ALL to default locale: No such file or directory

Sorry.  You can retry the autodetection at any time by running "wvdialconf".
   (Or you can create /etc/wvdial.conf yourself.)
Processing triggers for libc-bin (2.19-18+deb8u7) ...
Processing triggers for systemd (215-17+deb8u6) ...

　SORACOM Airでネットワーク接続するためのスクリプトをダウンロードして実行権限を付与します。

pi@raspberrypi:~ $ curl -O http://soracom-files.s3.amazonaws.com/connect_air.sh
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100  2843  100  2843    0     0  23012      0 --:--:-- --:--:-- --:--:-- 23113
pi@raspberrypi:~ $ 
pi@raspberrypi:~ $ ls -l
total 4
-rw-r--r-- 1 pi pi 2843 Apr 30 04:08 connect_air.sh
pi@raspberrypi:~ $ 
pi@raspberrypi:~ $ chmod +x connect_air.sh 
pi@raspberrypi:~ $ 
pi@raspberrypi:~ $ ls -l
total 4
-rwxr-xr-x 1 pi pi 2843 Apr 30 04:08 connect_air.sh
pi@raspberrypi:~ $ 
pi@raspberrypi:~ $ sudo mv connect_air.sh /usr/local/sbin/
pi@raspberrypi:~ $ 
pi@raspberrypi:~ $ ls -l /usr/local/sbin/connect_air.sh 
-rwxr-xr-x 1 pi pi 2843 Apr 30 04:08 /usr/local/sbin/connect_air.sh

　そしてUSBドングルにSimを入れて、Raspberry PiのUSBポートに挿し、スクリプトを実行します。

pi@raspberrypi:~ $ sudo /usr/local/sbin/connect_air.sh 
Found AK-020
Configuring modem ... done.
waiting for modem device..done.
Resetting modem ...done
could not initialize AK-020
waiting for modem device
--> WvDial: Internet dialer version 1.61
--> Cannot get information for serial port.
--> Initializing modem.
--> Sending: ATZ
ATZ
OK
--> Sending: ATQ0 V1 E1 S0=0 &C1 &D2 +FCLASS=0
ATQ0 V1 E1 S0=0 &C1 &D2 +FCLASS=0
OK
--> Sending: AT+CGDCONT=1,"IP","soracom.io"
AT+CGDCONT=1,"IP","soracom.io"
OK
--> Modem initialized.
--> Sending: ATD*99***1#
--> Waiting for carrier.
ATD*99***1#
CONNECT 21000000
--> Carrier detected.  Starting PPP immediately.
--> Starting pppd at Sun Apr 30 04:19:32 2017
--> Pid of pppd: 1873
--> Using interface ppp0
--> local  IP address 10.247.81.162
--> remote IP address 10.64.64.64
--> primary   DNS address 100.127.0.53
--> secondary DNS address 100.127.1.53

　SORACOMのユーザーコンソールから確認すると、SessionがOnlineになっていて、接続できていることが確認できます。

ifup されてきました（if-up2017に参加してきました）

　先日ソラコムさん主催の IoT Technology Conference if-up 2017 に行ってきました。ifup というのはLinuxでネットワークインタフェースを有効にするためのコマンドで、参加者それぞれの IoT Technology に関するインタフェースをUPしてほしいということでカンファレンスのタイトルに使われたそうです。セッションの内容はどれも興味深く、参考になるものばかりで、まんまと私も ifup されてきたわけですが、特にソラコムCTOの安川さんによる、「SORACOM Inside」という、SORACOMの設計思想やアーキテクチャを紹介されたキーノートの内容がとても参考になったので、その中で特に印象深かった点を紹介させていただきます。

　安川さんのセッション資料はこちらで公開されています。

if-up 2017 | SORACOM Inside from SORACOM,INC

www.slideshare.net

Polaris, Dipper, Hubble

　SORACOMプラットフォームは、パケット転送等の主要機能を担う Polaris（北極星）、認証や課金などの周辺機能を担う Dipper（北斗七星）、監視やデプロイを担う Hubble（宇宙望遠鏡）という要素で構成されているそう。プロダクトの役割と思想がマッチしたこういうネーミングはかっこいいなーと思いました。内部のメンバーとしても、それぞれのプロダクトの役割について共通の認識が持てるのではないかと思います。

ローンチ時からグローバル対応

　SORACOMユーザーコンソールはローンチ時からグローバル対応していて、HTML/JS/CSSによるSPAとして構成し、APIでSORACOMプラットフォームと連携し、S3にアップロードすることでCloudFrontで世界に配信されるようにしているとのこと。また、最初から多言語に対応し、タイムゾーンはUTCに統一しているとのことでした。一度普通のWebアプリとして構成してしまったものを後からSPAの構成に変更するのは大変そうですが、最初からこの構成を意識していれば、確かにS3にアップロードするだけでユーザーコンソールアプリがデプロイできるのは良いなと思いました。タイムゾーンもJSTをベースにしようとすると、色々な箇所でミスなくJSTに統一するのも大変ですし、UTCに統一することができれば後から世界展開するときにも確かに面倒なことが少なそうです。

疎結合化と非同期化

　SORACOMでは各サービス間の連携はs3を経由するなどして疎結合化されていて、一方に障害が起こっても、他方に直接的な影響がないようにしているとのこと。また、データ形式さえ決めておけば、それぞれのサービスの開発も非同期で行えるので、開発速度の向上にも有効なようです。いわゆるマイクロサービスの構成かと思いますが、前述のSPAの構成と同様で、一度モノリシックな構成で作ってしまうと後からマイクロサービスに分けるのは現実的ではなかったりします。スタートアップのサービス立ち上げ当初は特にモノリシックな一つのWebアプリとして一気に作ってしまうケースが多いかと思いますが、マイクロサービスの構成を最初から意識して作られている点はすごいと思いました。

DevOpsとOpsDev

　ソラコムの開発者は全員DevOpsを実践されているそうですが、運用の守りが手薄にならないよう、運用作業省力化のための開発を専門的に行うOpsDevエンジニアを導入されているそう。Hubbleの中でも障害を検知した時は、インスタンスの入れ替えで復旧できるようなものは、自動的に復旧させるようにしているとのことです。検知から通知まではどこの会社でも普通に実施していると思いますが、自動復旧までやれているケースはあまりないのではないでしょうか。また、OpsDevという考え方は聞いたことがありませんでしたが、運用省力化のための開発が好きなエンジニアがいるケースであれば、有効な選択肢だと思いました。うちの会社でも取り入れてみようかと思っています。

当たり前をちゃんとやれてるのがすごい

　講演後に安川さんとお話しさせていただいたときに、「内容としては当たり前のことばかりなので、講演前は聞いている方の反応がどうなるか不安があった」ということでした。マイクロサービスの考え方など、確かに一つ一つは広まっているものも多いですが、それぞれをしっかりとやり切れているというのはなかなかないと思いますし、とても参考になりました。また、システムの構成以外にも、ソラコムさんではリーダーシップ・ステートメントや、毎日行われているSyncというミーティングなど、チームアップの面でも参考にさせていただきたい点が多いので、またお話させていただきたいなぁと思っています。

　それと、今回参加者には参加特典としてIoTデバイスがプレゼントされていて、私は3G SIMのUSBドングルをもらいました。そこでこれを機にRaspberry Piも買ったので、色々やってみようと思います。

　今回一つ心残りだったのは、休憩時間中にまつもとゆきひろさんと玉川さんと3人でお話させていただく機会があったのですが、なかなかこのお二人と一緒に話をさせてもらう機会はないと思うので、3人で写真を撮らせてもらえば良かったと後になって思いました。次に機会があれば撮らせていただこうと思っています。

各種パラメータ最適化手法の実装（SGD, Momentum, AdaGrad, Adam）

　今回は「ゼロから作るDeepLearning」で紹介されている各種パラメータ最適化手法を、書籍のPythonのサンプルコードをベースに、Rubyで実装してみました。

www.oreilly.co.jp

　各手法のロジックについては書籍で説明されていますので割愛します。また、前回の記事で書いたように、Rubyでは値の受け渡しが参照の値渡しになるので、パラメータのハッシュの各値は配列として保持する前提です。

SGD（確率的勾配降下法）

　SGDは前回の記事でもすでに使っていたのと同じで、別クラスとして分けただけのものです。

class SGD
  def initialize(lr: 0.01)
    @lr = lr
  end

  def update(params:, grads:)
    params.keys.each do |key|
      params[key][0] -= @lr * grads[key]
    end
  end
end

Momentum

　paramsのハッシュの各値を配列として扱っている以外は、PythonのコードをそのままRubyに置き換えています。インスタンス変数 v にはparamsと同じ構造の値を持ちますが、インスタンス内でしか使わないため、ハッシュの値は配列にはせずにそのままパラメータを保持しています。

　ゼロ行列の生成は Numo::NArray.zeros メソッドを使っています。

Numo::NArray.zeros
http://ruby-numo.github.io/narray/narray/Numo/NArray.html#zeros-class_method

class Momentum
  def initialize(lr: 0.01, momentum: 0.9)
    @lr = lr
    @momentum = momentum
    @v = nil
  end

  def update(params:, grads:)
    if @v.nil?
      @v = {}
      params.each do |key, value|
        @v[key] = Numo::DFloat.zeros(value.first.shape)
      end
    end

    params.keys.each do |key|
      @v[key] = @momentum * @v[key] - @lr * grads[key]
      params[key][0] += @v[key]
    end
  end
end

AdaGrad

　こちらもMomentumと同様にPythonのコードを置き換えています。

　行列に対しての平方根の計算は、Numo::DFloat::Math.sqrt メソッドを使っています。

Numo::DFloat::Math.sqrt
http://ruby-numo.github.io/narray/narray/Numo/DFloat/Math.html#sqrt-class_method

class AdaGrad
  def initialize(lr: 0.01)
    @lr = lr
    @h = nil
  end

  def update(params:, grads:)
    if @h.nil?
      @h = {}
      params.each do |key, value|
        @h[key] = Numo::DFloat.zeros(value.first.shape)
      end
    end

    params.keys.each do |key|
      @h[key] += grads[key] * grads[key]
      params[key][0] -= @lr * grads[key] / (Numo::DFloat::Math.sqrt(@h[key]) + 1e-7)
    end
  end
end

Adam

　こちらも同様の置き換えです。

class Adam
  def initialize(lr: 0.001, beta1: 0.9, beta2: 0.999)
    @lr = lr
    @beta1 = beta1
    @beta2 = beta2
    @iter = 0
    @m = nil
    @v = nil
  end

  def update(params:, grads:)
    if @m.nil?
      @m = {}
      @v = {}
      params.each do |key, value|
        @m[key] = Numo::DFloat.zeros(value.first.shape)
        @v[key] = Numo::DFloat.zeros(value.first.shape)
      end
    end

    @iter += 1
    lr_t = @lr * Numo::DFloat::Math.sqrt(1.0 - @beta2 ** @iter) / (1.0 - @beta1 ** @iter)

    params.keys.each do |key|
      @m[key] += (1 - @beta1) * (grads[key] - @m[key])
      @v[key] += (1 - @beta2) * (grads[key] ** 2 - @v[key])

      params[key][0] -= lr_t * @m[key] / (Numo::DFloat::Math.sqrt(@v[key]) + 1e-7)
    end
  end
end

MNISTデータセットによる最適化手法の比較

　上記の最適化手法の実装について、MNISTデータセットを用いた学習の比較を行います。こちらも基本的には書籍のPython実装をベースにしていて、5 層のニューラルネットワークで、各層 100 個のニューロンを持つ ネットワークという構成です。活性化関数にはReLUを用いています。

　まずは複数レイヤのネットワークの処理を行うためのMultiLayerNetクラスの実装です。以前のTwoLayerNetを、3層以上のネットワークにも対応させた形です。

require 'numo/narray'
require './layers.rb'

class MultiLayerNet
  def initialize(input_size:, hidden_size_list:, output_size:, activation: :relu, weight_init_std: :relu, weight_decay_lambda: 0)
    @input_size          = input_size
    @output_size         = output_size
    @hidden_size_list    = hidden_size_list
    @hidden_layer_num    = hidden_size_list.size
    @weight_decay_lambda = weight_decay_lambda
    @params              = {}

    # 重みの初期化
    init_weight(weight_init_std)

    # レイヤの生成
    activation_layer = {
      sigmoid: Sigmoid,
      relu:    Relu
    }
    @layers = {}
    (1..@hidden_layer_num).each do |idx|
      @layers["Affine#{idx}"] = Affine.new(w: @params["w#{idx}"], b: @params["b#{idx}"])
      @layers["Activation_function#{idx}"] = activation_layer[activation].new
    end

    idx = @hidden_layer_num + 1
    @layers["Affine#{idx}"] = Affine.new(w: @params["w#{idx}"], b: @params["b#{idx}"])

    @last_layer = SoftmaxWithLoss.new
  end

  def params
    @params
  end

  def init_weight(weight_init_std)
    all_size_list = [@input_size] + @hidden_size_list + [@output_size]
    (1..(all_size_list.size - 1)).each do |idx|
      scale = weight_init_std
      if %i(relu he).include?(weight_init_std)
        scale = Numo::DFloat::Math.sqrt(2.0 / all_size_list[idx - 1])
      elsif %i(sigmoid xavier).include?(weight_init_std)
        scale = Numo::DFloat::Math.sqrt(1.0 / all_size_list[idx - 1])
      end

      Numo::NArray.srand
      @params["w#{idx}"] = [scale * Numo::DFloat.new(all_size_list[idx - 1], all_size_list[idx]).rand_norm]
      @params["b#{idx}"] = [Numo::DFloat.zeros(all_size_list[idx])]
    end
  end

  def predict(x:)
    @layers.values.inject(x) do |x, layer|
      x = layer.forward(x: x)
    end
  end

  # x: 入力データ, t: 教師データ
  def loss(x:, t:)
    y = predict(x: x)

    weight_decay = 0
    (1..(@hidden_layer_num + 1)).each do |idx|
      w = @params["w#{idx}"].first
      weight_decay += 0.5 * @weight_decay_lambda * (w ** 2).sum
    end
    @last_layer.forward(x: y, t: t) + weight_decay
  end

  def accuracy(x:, t:)
    y = predict(x: x)
    y = y.max_index(1) % 10
    if t.ndim != 1
      t = t.max_index(1) % 10
    end

    y.eq(t).cast_to(Numo::UInt16).sum / x.shape[0].to_f
  end

  def gradient(x:, t:)
    # forward
    loss(x: x, t: t)

    # backward
    dout = 1
    dout = @last_layer.backward(dout: dout)

    layers = @layers.values.reverse
    layers.inject(dout) do |dout, layer|
      dout = layer.backward(dout: dout)
    end

    grads = {}
    (1..(@hidden_layer_num + 1)).each do |idx|
      grads["w#{idx}"] = @layers["Affine#{idx}"].dw + @weight_decay_lambda * @layers["Affine#{idx}"].w.first
      grads["b#{idx}"] = @layers["Affine#{idx}"].db
    end

    grads
  end
end

　そして各手法での学習と、結果のグラフ描画を行うためのスクリプトの実装です。基本的な処理は前回の誤差逆伝播法での学習処理と同じで、SGD以外にもMomentum、AdaGrad、Adamでの学習を行い、結果を比較しています。

require 'numo/gnuplot'
require './mnist.rb'
require './optimizers.rb'
require './multi_layer_net.rb'

# 0: MNISTデータの読み込み
x_train, t_train, x_test, t_test = load_mnist(normalize: true)

train_size = x_train.shape[0]
batch_size = 128
max_iterations = 1500

# 1: 実験の設定
optimizers = {
  sgd:      SGD.new,
  momentum: Momentum.new,
  adagrad:  AdaGrad.new,
  adam:     Adam.new
}

networks = {}
train_loss = {}
optimizers.each do |key, optimizer|
  networks[key]   = MultiLayerNet.new(input_size: 784, hidden_size_list: [100, 100, 100, 100], output_size: 10)
  train_loss[key] = []
end

# 2: 訓練の開始
max_iterations.times do |i|
  Numo::NArray.srand
  batch_mask = Numo::Int32.new(batch_size).rand(0, train_size)
  x_batch = x_train[batch_mask, true]
  t_batch = t_train[batch_mask]

  optimizers.each do |key, optimizer|
    grads = networks[key].gradient(x: x_batch, t: t_batch)
    optimizers[key].update(params: networks[key].params, grads: grads)

    loss = networks[key].loss(x: x_batch, t: t_batch)
    train_loss[key] << loss
  end

  next unless i % 100 == 0

  puts "========== iteration: #{i} =========="
  optimizers.keys.each do |key|
    loss = networks[key].loss(x: x_batch, t: t_batch)
    puts "#{key}: #{loss}"
  end
end

# 3: グラフの描画
x = (0..(max_iterations - 1)).to_a
Numo.gnuplot do
  set xlabel: 'iterations'
  set ylabel: 'loss'
  set yrange: 0...1
  plot x, train_loss[:sgd],      { w: :lines, t: 'SGD',      lc_rgb: 'green',  lw: 1 },
       x, train_loss[:momentum], { w: :lines, t: 'Momentum', lc_rgb: 'orange', lw: 1 },
       x, train_loss[:adagrad],  { w: :lines, t: 'AdaGrad',  lc_rgb: 'red',    lw: 1 },
       x, train_loss[:adam],     { w: :lines, t: 'Adam',     lc_rgb: 'blue',   lw: 1 }
end

　これをirbから実行すると下記のようにコンソールに結果が100イテレーションごとに表示され、最後にグラフが描画されます。

irb(main):001:0> load './optimizer_compare_mnist.rb'
========== iteration: 0 ==========
sgd: 2.490228492804417
momentum: 2.492025460726973
adagrad: 2.0825493921580134
adam: 2.276848398313694
========== iteration: 100 ==========
sgd: 1.5001470235084333
momentum: 0.42706984858496944
adagrad: 0.2004084345838115
adam: 0.34096154614776963
========== iteration: 200 ==========
sgd: 0.7664658855425125
momentum: 0.2733076953685949
adagrad: 0.11048856792711875
adam: 0.24828428524427817
========== iteration: 300 ==========
sgd: 0.5159466996794285
momentum: 0.25594092908625543
adagrad: 0.13587365073017388
adam: 0.20545236418926946
========== iteration: 400 ==========
sgd: 0.5042479159281286
momentum: 0.24538385847033825
adagrad: 0.11124732792005954
adam: 0.1797581753729203
========== iteration: 500 ==========
sgd: 0.32290967978019125
momentum: 0.1599522422679423
adagrad: 0.05731788233265379
adam: 0.0888823836035264
========== iteration: 600 ==========
sgd: 0.44467997494741673
momentum: 0.2578398459161452
adagrad: 0.1316129116477675
adam: 0.22439066383913017
========== iteration: 700 ==========
sgd: 0.28407622085704987
momentum: 0.10056311655844065
adagrad: 0.07989693533502204
adam: 0.0821317626167635
========== iteration: 800 ==========
sgd: 0.2706466278682429
momentum: 0.15550352523100197
adagrad: 0.06489312962717893
adam: 0.08528336870483003
========== iteration: 900 ==========
sgd: 0.2240422184352822
momentum: 0.11062658792202897
adagrad: 0.059913720263603615
adam: 0.03302573552710864
========== iteration: 1000 ==========
sgd: 0.3832020077768542
momentum: 0.13726781722583942
adagrad: 0.03417701415203686
adam: 0.053558781255080776
========== iteration: 1100 ==========
sgd: 0.38619949224379774
momentum: 0.15175966760909282
adagrad: 0.04222220798211423
adam: 0.06822940475295906
========== iteration: 1200 ==========
sgd: 0.2998755819916694
momentum: 0.07572742725924923
adagrad: 0.075962654676941
adam: 0.02748595912322749
========== iteration: 1300 ==========
sgd: 0.25322815781566416
momentum: 0.06003606774412698
adagrad: 0.03236788855975958
adam: 0.053987864918752786
========== iteration: 1400 ==========
sgd: 0.2720482764348912
momentum: 0.09105631835160209
adagrad: 0.044338756972504875
adam: 0.0668196066196452
=> true

　イテレーションの回数を1,500回としていますが、手元のVM環境ではこれ以上回数を増やすと途中で強制終了されてしまいました。Pythonコードをそのままの構成で移植しているので、もっとRubyに最適化した実装を考慮する必要がありそうです。

　Numo.gnuplot でのグラフ描画時は、setメソッドによる設定はplotより前に実行しておかないとグラフに反映されないようでした。

　コードは下記リポジトリにも公開しています。

github.com

ニューラルネットワークの誤差逆伝播法による学習アルゴリズムの実装

　今回も引き続き「ゼロから作るDeepLearning」をベースに、前回数値微分で実装した学習アルゴリズムの誤差逆伝播法版をRubyで実装してみました。計算の内容等は書籍を参照いただくとして、Rubyで実装した際のポイントを説明していきたいと思います。

www.oreilly.co.jp

活性化関数レイヤ実装

　今回はニューラルネットワークの各レイヤをそれぞれ一つのクラスとして実装しています。活性化関数レイヤはReLUレイヤとSigmoidレイヤです。

ReLUレイヤ

　ReLUレイヤは下記のように実装しました。

class Relu
  def initialize
    @mask = nil
  end

  def forward(x:)
    @mask = (x <= 0)
    out = x.copy
    out[@mask] = 0
    out
  end

  def backward(dout:)
    dout[@mask] = 0
    dout
  end
end

　順伝播時のforwardメソッドの引数としてはNArray配列を期待しています。NArray配列 x に対して (x <= 0) をすると0以下の要素が1、それ以外が0のBit配列を返しますので、それをマスクに使い、入力に対して0以下の要素を0に置き換えた結果を返しています。また、マスクは逆伝播時にも使うのでインスタンス変数に保持します。

　逆伝播時のbackwardメソッドでは順伝播時に保持したマスクを元に、上流からの入力に対してマスクの要素が1になっている要素に0を設定しています。

Sigmoidレイヤ

　Sigmoidレイヤの実装は下記のように行いました。

require './sigmoid.rb'

class Sigmoid
  def initialize
    @out = nil
  end

  def forward(x:)
    @out = sigmoid(x)
    @out
  end

  def backword(dout:)
    dout * (1.0 - @out) * @out
  end
end

　順伝播時は以前実装したsigmoidメソッドを呼んでいるだけで、その結果をインスタンス変数に保持しておきます。

　逆伝播時は順伝播時の出力を元に計算を行なった結果を返します。

AffineレイヤとSoftmaxレイヤの実装

　ニューラルネットワークの順伝播において重みの計算とバイアスの加算を行なっていた層をAffineレイヤとして実装し、出力層ではソフトマックス関数を用いて出力を正規化するSoftmaxレイヤを実装し、それぞれの順伝播、逆伝播の処理を実装します。

Affineレイヤ

　Affineレイヤの実装は下記の通りです。

class Affine
  def initialize(w:, b:)
    @w = w
    @b = b

    @x = nil
    @original_x_shape = nil

    # 重み・バイアスパラメータの微分
    @dw = nil
    @db = nil
  end

  def dw
    @dw
  end

  def db
    @db
  end

  def forward(x:)
    # テンソル対応
    @original_x_shape = x.shape
    x = x.reshape(x.shape[0], nil)
    @x = x

    @x.dot(@w.first) + @b.first
  end

  def backward(dout:)
    dx = dout.dot(@w.first.transpose)
    @dw = @x.transpose.dot(dout)
    @db = dout.sum(0)

    dx.reshape(*@original_x_shape)
  end
end

　initializeメソッドでは重み、バイアスパラメータをインスタンス変数に格納し、それ以外にも処理に必要になるインスタンス変数を定義しています。

　順伝播時は入力の行列の形と入力行列を保持し、入力行列と重みパラメータの内積にバイアスを加算した結果を返します。

　逆伝播時はまず内積の逆伝播の計算として重みパラメータと順伝播時の入力値の転置行列を使った計算を行なっています。NArray行列の転置行列はtransposeメソッドで取得できます。

Numo::NArray#transpose
http://ruby-numo.github.io/narray/narray/Numo/NArray.html#transpose-instance_method

　バイアスの逆伝播の計算では行列の0番目の軸（データ単位）に対しての合計を求めるため、sumメソッドのパラメータに0を指定しています。

Numo::UInt32#sum
http://ruby-numo.github.io/narray/narray/Numo/Int32.html#sum-instance_method

　最後に入力値の逆伝播の計算結果を入力値と同じ行列の形にreshapeして返しています。入力値の形状は変数に配列データとして格納していますが、reshapeメソッドのパラメータは配列ではないので、* で配列を展開して渡しています。

Softmaxレイヤ

　Softmaxレイヤは損失関数である交差エントロピー誤差の計算も含めて、SoftmaxWithLossクラスとして下記のように実装しました。

require './softmax.rb'
require './cross_entropy_error.rb'

class SoftmaxWithLoss
  def initialize
    @loss = nil
    @y = nil # softmaxの出力
    @t = nil # 教師データ
  end

  def forward(x:, t:)
    @t = t
    @y = softmax(x)
    @loss = cross_entropy_error(@y, @t)

    @loss
  end

  def backward(dout: 1)
    batch_size = @t.shape[0]
    if @t.size == @y.size # 教師データがon-hot-vectorの場合
      return (@y - @t) / batch_size
    end

    dx = @y.copy

    (0..(batch_size - 1)).to_a.zip(@t).each do |index_array|
      dx[*index_array] -= 1
    end

    dx / batch_size
  end
end

　順伝播時は以前実装したsoftmaxメソッドとcross_entropy_errorメソッドを使用しています。入力値をsoftmaxメソッドで正規化し、その結果と教師データをcross_entropy_errorメソッドに渡して交差エントロピー誤差を計算しています。

　逆伝播時はデータ一個あたりの誤差を伝播するために誤差をデータ数で割っています。

　行列の要素の参照方法については、Pythonのndarrayの場合は配列を二つ渡すと、それぞれを行・列のインデックスとして要素を参照してくれますが、NArray行列で同じような指定をすると、一つ目の配列で指定した全ての行で、二つ目の配列に指定した全ての要素が取得されてしまいます。

>>> array
array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12]])
>>> array[[0, 2], [1, 3]]
array([ 2, 12])

irb(main):021:0> array                
=> Numo::UInt32#shape=[3,4]           
[[1, 2, 3, 4],                        
 [5, 6, 7, 8],                        
 [9, 10, 11, 12]]                     
irb(main):022:0> array[[0, 2], [1, 3]]
=> Numo::UInt32(view)#shape=[2,2]     
[[2, 4],                              
 [10, 12]]                            

　そこでzipメソッドを使ってそれぞれの配列から一つずつデータを取得し、直接一つの要素を参照するようにしました。

irb(main):024:0* [0, 2].zip([1, 3]).each do |array_index| 
irb(main):025:1*   puts array[*array_index]               
irb(main):026:1> end                                      
2                                                         
12                                                        
=> [[0, 1], [2, 3]]                                       

誤差逆伝播法でのニューラルネットワーク実装

　上記で実装したレイヤを組み合わせて、誤差逆伝播法でのニューラルネットワークを実装します。下記のように一つのクラスとして実装します。

require 'numo/narray'
require './numerical_gradient.rb'
require './layers.rb'

class TwoLayerNet
  def initialize(input_size:, hidden_size:, output_size:, weight_init_std: 0.01)
    # 重みの初期化
    Numo::NArray.srand
    @params = {
      w1: [weight_init_std * Numo::DFloat.new(input_size, hidden_size).rand_norm],
      b1: [Numo::DFloat.zeros(hidden_size)],
      w2: [weight_init_std * Numo::DFloat.new(hidden_size, output_size).rand_norm],
      b2: [Numo::DFloat.zeros(output_size)]
    }

    # レイヤの生成
    @layers = {
      affine1: Affine.new(w: @params[:w1], b: @params[:b1]),
      relu1:   Relu.new,
      affine2: Affine.new(w: @params[:w2], b: @params[:b2])
    }
    @last_layer = SoftmaxWithLoss.new
  end

  def params
    @params
  end

  def predict(x:)
    @layers.values.inject(x) do |x, layer|
      x = layer.forward(x: x)
    end
  end

  # x: 入力データ, t: 教師データ
  def loss(x:, t:)
    y = predict(x: x)
    @last_layer.forward(x: y, t: t)
  end

  def accuracy(x:, t:)
    y = predict(x: x)
    y = y.max_index(1) % 10
    if t.ndim != 1
      t = t.max_index(1) % 10
    end

    y.eq(t).cast_to(Numo::UInt16).sum / x.shape[0].to_f
  end

  def numerical_gradients(x:, t:)
    loss_w = lambda { loss(x: x, t: t) }

    {
      w1: numerical_gradient(loss_w, @params[:w1].first),
      b1: numerical_gradient(loss_w, @params[:b1].first),
      w2: numerical_gradient(loss_w, @params[:w2].first),
      b2: numerical_gradient(loss_w, @params[:b2].first)
    }
  end

  def gradient(x:, t:)
    # forward
    loss(x: x, t: t)

    # backward
    dout = 1
    dout = @last_layer.backward(dout: dout)

    layers = @layers.values.reverse
    layers.inject(dout) do |dout, layer|
      dout = layer.backward(dout: dout)
    end

    {
      w1: @layers[:affine1].dw,
      b1: @layers[:affine1].db,
      w2: @layers[:affine2].dw,
      b2: @layers[:affine2].db
    }
  end
end

　initializeでは重みとバイアスパラメータの初期化と、各レイヤのインスタンスの生成を行なっています。重みとバイアスのパラメータを配列として保持しているのは、学習結果をTwoLayerNetクラス内のパラメータに反映したものを各レイヤで参照できるようにするためです。Affineクラスのインスタンス生成時にパラメータを w: @params[:w1] という形で渡していますが、Rubyでは参照の値渡しになるため、TwoLayerNet側で @params[:w1] に計算結果を代入しても、 Affineインスタンス側では初期化時に渡されたパラメータを参照し続けます。そこで参照先を配列にして、計算結果はその配列の中身を更新する形にしています。書籍のPythonコードでは参照渡しになるため、配列として保持しなくてもTwoLayerNet側での変更がAffineインスタンス側で参照されています。

　それと、パラメータをランダムに生成する前に、Numo::NArray.srandメソッドでseedが変わるようにしています。これをしないとスクリプト実行時に毎回同じ内容のパラメータが作成されてしまいます。

　また、ニューラルネットワークでは各レイヤが実行される順番が重要なので、Pythonの場合は通常のDictionaryではなくOrderedDictを使っていますが、RubyのHashでは順番が保持されるため、通常のHashをそのまま使っています。

　勾配の計算処理のgradientメソッドでは、まず順伝播の処理を行うためにlossメソッドを実行します。lossメソッドではpredictメソッドを呼び出し、injectメソッドで各レイヤのforward処理を実行し、順伝播の処理を行なっています。そして最後に last_layer変数に保持しているSoftmaxWithLossインスタンスの順伝播処理で交差エントロピー誤差を計算しています。

　続いて逆伝播処理ではまずSoftmaxWithLossインスタンスの逆伝播処理を行なったあと、各レイヤを逆順に逆伝播処理を行ない、計算結果を返しています。

　精度確認用のaccuracyメソッドでは、まず各レイヤの順伝播処理を行い、その結果一番確度の高い要素のインデックスと教師データを付き合わせて正解率を計算しています。max_indexメソッドでは各データの最大値を判定したいので、引数で1軸目を指定しています（0だと全ての要素の中からの最大値を判定する）。結果は行列全ての要素の中でのインデックス値を返すので、10で割ることで各データのインデックス値に変換しています。

　eqメソッドでは一致する要素は1、一致しない要素は0のNumo::Bit配列を返します。正解数としてBitが1の要素の合計を計算するため、Numo::Bit配列をcast_toメソッドでInt配列に変換しています。この時、各ビット値は1か0なので、Int8でも正しく変換できますが、合計を計算するときにInt8の範囲を超えると正しく合計値が取得できません。今回の入力データの件数は60,000件あるので、Int8だと範囲を超えてしまうため、UInt16に変換した上で合計値を取得しています。

誤差逆伝播法の勾配確認

　誤差逆伝播法で求めた勾配が正しいかどうかを確認するため、数値微分で求めた勾配と比較して確認します。

require './mnist.rb'
require './two_layer_net.rb'

x_train, t_train, x_test, t_test = load_mnist(normalize: true, one_hot_label: true)

network = TwoLayerNet.new(input_size: 784, hidden_size: 50, output_size: 10)

x_batch = x_train[0..2, true]
t_batch = t_train[0..2, true]

grad_numerical = network.numerical_gradients(x: x_batch, t: t_batch)
grad_backprop = network.gradient(x: x_batch, t: t_batch)

grad_numerical.keys.each do |key|
  diff = (grad_backprop[key] - grad_numerical[key]).abs.mean
  puts "#{key}: #{diff}"
end

　MNISTデータの先頭3件を使い、数値微分での計算と誤差逆伝播法での計算を一度行い、それぞれの結果の差を計算しています。実行結果は下記のようになり、ほぼ差がないことが確認できます。

[vagrant@localhost vagrant]$ ruby gradient_check.rb
w1: 2.616843054226442e-13                          
b1: 8.347379796928845e-13                          
w2: 1.0232621862468414e-12                         
b2: 1.2012612987666316e-10                         

誤差逆伝播法を使った学習

　前回の記事で実装したミニバッチ学習と評価を、誤差逆伝播法を使うように変更します。前回との違いはnumerical_gradientsメソッドではなくgradientメソッドを使うようにした点です。

require 'numo/narray'
require 'numo/gnuplot'
require './mnist.rb'
require './two_layer_net.rb'

# データの読み込み
x_train, t_train, x_test, t_test = load_mnist(normalize: true, one_hot_label: true)

network = TwoLayerNet.new(input_size: 784, hidden_size: 50, output_size: 10)

iters_num = 10_000 # 繰り返し回数
train_size = x_train.shape[0]
batch_size = 100
learning_rate = 0.1

train_loss_list = []
train_acc_list = []
test_acc_list = []

iter_per_epoch = [train_size / batch_size, 1].max

iters_num.times do |i|
  Numo::NArray.srand
  batch_mask = Numo::Int32.new(batch_size).rand(0, train_size)
  x_batch = x_train[batch_mask, true]
  t_batch = t_train[batch_mask, true]

  # 勾配の計算
  #grad = network.numerical_gradients(x_batch, t_batch)
  grad = network.gradient(x: x_batch, t: t_batch)

  # パラメータの更新
  %i(w1 b1 w2 b2).each do |key|
    network.params[key][0] -= learning_rate * grad[key]
  end

  loss = network.loss(x: x_batch, t: t_batch)
  train_loss_list << loss

  next if i % iter_per_epoch != 0

  train_acc = network.accuracy(x: x_train, t: t_train)
  test_acc = network.accuracy(x: x_test, t: t_test)
  train_acc_list << train_acc
  test_acc_list << test_acc
  puts "train acc, test acc | #{train_acc}, #{test_acc}"
end

# グラフの描画
x = (0..(train_acc_list.size - 1)).to_a
Numo.gnuplot do
  plot x, train_acc_list, { w: :lines, t: 'train acc', lc_rgb: 'blue' },
       x, test_acc_list, { w: :lines, t: 'test acc', lc_rgb: 'green' }
  set xlabel: 'epochs'
  set ylabel: 'accuracy'
  set yrange: 0..1
end

　実行結果は下記のようになります。シェルからスクリプトファイルを実行してもグラフの描画が行われなかったので、irbからロードすることで実行し、グラフが表示されるようにしています。

irb(main):001:0> load './train_neuralnet.rb'
train acc, test acc | 0.12578333333333333, 0.1225
train acc, test acc | 0.9026833333333333, 0.9071
train acc, test acc | 0.92225, 0.923
train acc, test acc | 0.9348833333333333, 0.9331
train acc, test acc | 0.9436, 0.9403
train acc, test acc | 0.9513666666666667, 0.9503
train acc, test acc | 0.9568833333333333, 0.9565
train acc, test acc | 0.9614666666666667, 0.9594
train acc, test acc | 0.96415, 0.9611
train acc, test acc | 0.9659, 0.9616
train acc, test acc | 0.9677833333333333, 0.9622
train acc, test acc | 0.97175, 0.9655
train acc, test acc | 0.97275, 0.9666
train acc, test acc | 0.9745666666666667, 0.968
train acc, test acc | 0.9764166666666667, 0.9682
train acc, test acc | 0.9780666666666666, 0.9696
train acc, test acc | 0.9788166666666667, 0.9704
=> true

実行速度

　手元のVagrant環境で実行した結果、Ruby版とPython版のそれぞれの実行速度は下記のようになりました。

Ruby: 4分41秒 Python: 39秒

　Python版の方がかなり早いです。今回はPython版をベースにRuby版を実装したので、パフォーマンスチューニングが可能かやってみたいところです。

　コードは下記リポジトリでも公開しています。

github.com

ニューラルネットワークの数値微分による学習アルゴリズムの実装

　今回も引き続き「ゼロから作るDeepLearning」をベースに、数値微分による学習アルゴリズムをRubyで実装してみました。

www.oreilly.co.jp

　最初に書いておくと、今回の数値微分での実装は、実装はシンプルなもののその分処理に時間がかかり、手元の環境では繰り返し処理を終わらせることができませんでした。次回以降で書籍の次の章で解説されている、誤差逆伝播法での実装に変更してみたいと思います。

2層ニューラルネットワークのクラス

　書籍のPython実装に基づき、2層のニューラルネットワークをTwoLayerNetという一つのクラスとして実装します。まずはコード全体です。

require 'numo/narray'
require './softmax.rb'
require './sigmoid.rb'
require './cross_entropy_error.rb'
require './numerical_gradient.rb'

class TwoLayerNet
  def initialize(input_size, hidden_size, output_size, weight_init_std = 0.01)
    @params = {}
    @params[:w1] = weight_init_std * Numo::DFloat.new(input_size, hidden_size).rand_norm
    @params[:b1] = Numo::DFloat.zeros(hidden_size)
    @params[:w2] = weight_init_std * Numo::DFloat.new(hidden_size, output_size).rand_norm
    @params[:b2] = Numo::DFloat.zeros(output_size)
  end

  def params
    @params
  end

  def predict(x)
    w1 = @params[:w1]
    w2 = @params[:w2]
    b1 = @params[:b1]
    b2 = @params[:b2]

    a1 = x.dot(w1) + b1
    z1 = sigmoid(a1)
    a2 = z1.dot(w2) + b2
    softmax(a2)
  end

  def loss(x, t)
    y = predict(x)
    cross_entropy_error(y, t)
  end

  def accuracy(x, t)
    y = predict(x)
    y = y.max_index(1) % 10
    t = t.max_index(1) % 10

    y.eq(t).cast_to(Numo::UInt32).sum / x.shape[0].to_f
  end

  def numerical_gradients(x, t)
    loss_w = lambda {|w| loss(x, t) }

    grads = {}
    grads[:w1] = numerical_gradient(loss_w, @params[:w1])
    grads[:b1] = numerical_gradient(loss_w, @params[:b1])
    grads[:w2] = numerical_gradient(loss_w, @params[:w2])
    grads[:b2] = numerical_gradient(loss_w, @params[:b2])

    grads
  end
end

初期化

　まずはinitializeメソッドで重みとバイアスのパラメータを初期化し、インスタンス変数にHashとして保持します。合わせて外部からパラメータを参照するためのメソッドも用意しておきます。

def initialize(input_size, hidden_size, output_size, weight_init_std = 0.01)
  @params = {}
  @params[:w1] = weight_init_std * Numo::DFloat.new(input_size, hidden_size).rand_norm
  @params[:b1] = Numo::DFloat.zeros(hidden_size)
  @params[:w2] = weight_init_std * Numo::DFloat.new(hidden_size, output_size).rand_norm
  @params[:b2] = Numo::DFloat.zeros(output_size)
end

def params
  @params
end

　initializeメソッドの引数には、入力層のニューロン数、隠れ層のニューロン数、出力層のニューロン数を渡します。

　重みの初期値の生成は、まず Numo::DFloat.new で “入力層ニューロン数 x 隠れ層ニューロン数” もしくは “隠れ層ニューロン数 x 出力層ニューロン数” の行列を用意し、rand_normメソッドで標準正規分布に従う乱数を設定します。

Class: Numo::DFloat — Documentation by YARD 0.9.8

　バイアスの初期値は Numo::DFloat.zeros メソッドで全てゼロの配列を用意します。

Class: Numo::NArray — Documentation by YARD 0.9.8

推論処理

　画像データを引数にとり、推論処理を行います。

def predict(x)
  w1 = @params[:w1]
  w2 = @params[:w2]
  b1 = @params[:b1]
  b2 = @params[:b2]

  a1 = x.dot(w1) + b1
  z1 = sigmoid(a1)
  a2 = z1.dot(w2) + b2
  softmax(a2)
end

　処理内容としては前回の記事で実装したものと同じで、隠れ層の活性化関数にシグモイド関数、出力層の活性化関数にソフトマックス関数を使用しています。

def sigmoid(x)
  1 / (1 + Numo::DFloat::Math.exp(-x))
end

def softmax(a)
  if a.ndim == 2
    a = a.transpose
    a = a - a.max(0)
    y = Numo::DFloat::Math.exp(a) / Numo::DFloat::Math.exp(a).sum(0)
    return y.transpose
  end

  c = a.max
  exp_a = Numo::DFloat::Math.exp(a - c)
  sum_exp_a = exp_a.sum
  exp_a / sum_exp_a
end

損失関数

　入力データ（画像データ）と教師データ（正解ラベル）を引数にとり、損失関数の値を求めます。

def loss(x, t)
  y = predict(x)
  cross_entropy_error(y, t)
end

　画像データからpredictメソッドで推論を行なった結果と正解ラベルから交差エントロピー誤差を求めています。

def cross_entropy_error(y, t)
  if y.ndim == 1
    t = t.reshape(1, t.size)
    y = y.reshape(1, y.size)
  end

  batch_size = y.shape[0]
  -(t * (Numo::DFloat::Math.log(y))).sum / batch_size # one-hot表現用
end

認識精度の計算

　入力データ（画像データ）と教師データ（正解ラベル）を引数にとり、認識精度を計算します。

def accuracy(x, t)
  y = predict(x)
  y = y.max_index(1) % 10
  t = t.max_index(1) % 10

  y.eq(t).cast_to(Numo::UInt32).sum / x.shape[0].to_f
end

　画像データからpredictメソッドで推論を行なった結果は0〜9の数字である確度（確率）の配列として返されるため、どの数字である確率が最も高いかを max_index メソッドを用いて取得します。推論はバッチ処理で行うため、複数の画像データに対する結果として二次元配列として返ってくるので、二次元目を基準に最大値を求めるため、max_indexの引数には1を渡しています（0次元目、1次元目と考えるため）。ただし、max_index では多次元配列の場合、全ての要素数に対してのインデックス値（10 x 10の配列だったら0〜99の値）として戻ってくるため、10で割ったあまりを求めることで、0〜9のラベルを表すようにしています。

y = y.max_index(1) % 10

　これは正解ラベルについても同様で、 one-hot表現で渡されている複数データについての正解ラベルは二次元配列なので、10で割ることで0〜9のラベルを表すようにしています。

t = t.max_index(1) % 10

　そして最後に推論結果のラベル配列と正解ラベル配列をeqメソッドで比較し、正解数（配列の値が1になっている数）を入力データの数で割ることで、正解率を求めています。この辺りのeqメソッドの使い方等は前回記事と同様です。

重みパラメータに対する勾配の計算

　入力データ（画像データ）と教師データ（正解ラベル）を引数にとり、各パラメータに対する勾配を計算します。

def numerical_gradients(x, t)
  loss_w = lambda {|w| loss(x, t) }

  grads = {}
  grads[:w1] = numerical_gradient(loss_w, @params[:w1])
  grads[:b1] = numerical_gradient(loss_w, @params[:b1])
  grads[:w2] = numerical_gradient(loss_w, @params[:w2])
  grads[:b2] = numerical_gradient(loss_w, @params[:b2])

  grads
end

　損失関数を計算結果を取得するlambdaを用意し、勾配計算用のnumerical_gradientメソッドに各パラメータ共に渡し、結果をHashに格納して返します。

def numerical_gradient(f, x)
  h = 1e-4
  grad = Numo::DFloat.zeros(x.shape)

  x.size.times do |i|
    tmp_val = x[i]

    x[i] = tmp_val + h
    fxh1 = f.call(x)

    x[i] = tmp_val - h
    fxh2 = f.call(x)

    grad[i] = (fxh1 - fxh2) / (2 * h)
    x[i] = tmp_val
  end

  grad
end

　numerical_gradientメソッドでは中心差分での数値微分によって各パラメータの勾配を計算します。まず結果の格納用に入力パラメータと同じ形のゼロ配列を用意します。そして入力パラメータの各要素について (f(x + h) - f(x - h)) / 2h を計算して結果格納用の配列に格納し、最後にその配列を返しています。Numo::NArrayの多次元配列では単一のインデックスで配列の内容を参照する場合、全要素をフラットに並べた時のインデックス値を意味するので、x.size で全要素数を取得して、timesでその要素数分ループを回して、順次配列の内容を参照するようにしています。

ミニバッチ学習と評価の実装

　TwoLayerNetクラスを使ってMNISTデータセットを用いた学習と評価を行います。学習の実装は、訓練データから無作為に一部のデータを取り出して入力データとする、ミニバッチ学習で行います。まずはコード全体です。

require 'numo/narray'
require 'numo/gnuplot'
require './mnist.rb'
require './two_layer_net.rb'

# データの読み込み
x_train, t_train, x_test, t_test = load_mnist(true, true, true)

network = TwoLayerNet.new(784, 50, 10)

iters_num = 10_000 # 繰り返し回数
train_size = x_train.shape[0]
batch_size = 100
learning_rate = 0.1

train_loss_list = []
train_acc_list = []
test_acc_list = []

iter_per_epoch = [train_size / batch_size, 1].max

iters_num.times do |i|
  batch_mask = Numo::Int32.new(batch_size).rand(0, train_size)
  x_batch = x_train[batch_mask, true]
  t_batch = t_train[batch_mask, true]

  # 勾配の計算
  grad = network.numerical_gradients(x_batch, t_batch)

  # パラメータの更新
  %i(w1 b1 w2 b2).each do |key|
    network.params[key] -= learning_rate * grad[key]
  end

  loss = network.loss(x_batch, t_batch)
  train_loss_list << loss

  next if i % iter_per_epoch != 0

  train_acc = network.accuracy(x_train, t_train)
  test_acc = network.accuracy(x_test, t_test)
  train_acc_list << train_acc
  test_acc_list << test_acc
  puts "train acc, test acc | #{train_acc}, #{test_acc}"
end

# グラフの描画
x = (0..(train_acc_list.size - 1)).to_a
Numo.gnuplot do
  plot x, train_acc_list, { w: :lines, t: 'train acc', lc_rgb: 'blue' },
       x, test_acc_list, { w: :lines, t: 'test acc', lc_rgb: 'green' }
  set xlabel: 'epochs'
  set ylabel: 'accuracy'
  set yrange: 0..1
end

MNISTデータの読み込みとパラメータ初期化

　まずは以前実装したMNISTデータのロード処理を使ってMNISTデータを読み込んだ後、TwoLayerNetの初期化と、繰り返し数等の設定を行います。

x_train, t_train, x_test, t_test = load_mnist(true, true, true)

network = TwoLayerNet.new(784, 50, 10)

iters_num = 10_000 # 繰り返し回数
train_size = x_train.shape[0]
batch_size = 100
learning_rate = 0.1

　ニューロン数の構成は、入力層は28 x 28の画像データなので 784、出力層は0-9の数字を表すので 10、隠れ層は50としています。また、確率勾配降下法によるパラメータ更新の繰り返し数は10,000回で、ミニバッチのサイズは100としています。

ミニバッチデータの取得

　繰り返し処理の中では、まず60,000件の画像データの中からミニバッチデータを取得しています。

batch_mask = Numo::Int32.new(batch_size).rand(0, train_size)
x_batch = x_train[batch_mask, true]
t_batch = t_train[batch_mask, true]

　Numo::Int32.new でバッチサイズ分の要素数の配列を用意し、randメソッドで0以上60,000未満のランダム値を生成しています。これによって、60,000件の画像データのどのインデックス値のデータを取得するかを決めています。そしてそのインデックス値に該当する画像データと正解ラベルを取得していますが、それぞれ 60,000 x 784、60,000 x 10 の二次元配列になっているため、ランダムに生成したインデックス値の配列が一次元目を表すことを意味するように、二次元目はtrueを指定し、該当する二次元目のデータは全て取得するようにしています。

勾配の計算とパラメータの更新

　ミニバッチデータをTwoLayerNetの勾配計算メソッドに渡して勾配データを取得し、学習率をかけたものを各パラメータから引くことで、各パラメータを更新します。

# 勾配の計算
grad = network.numerical_gradients(x_batch, t_batch)

# パラメータの更新
%i(w1 b1 w2 b2).each do |key|
  network.params[key] -= learning_rate * grad[key]
end

損失関数の計算と記録

　ミニバッチデータから損失関数を計算し、繰り返しごとの結果を格納し、経過を記録します。

loss = network.loss(x_batch, t_batch)
train_loss_list << loss

認識精度の計算

　トレーニングデータとテストデータに対する認識精度を計算します。

next if i % iter_per_epoch != 0

train_acc = network.accuracy(x_train, t_train)
test_acc = network.accuracy(x_test, t_test)
train_acc_list << train_acc
test_acc_list << test_acc
puts "train acc, test acc | #{train_acc}, #{test_acc}"

　繰り返しごとに毎回認識精度を計算すると時間がかかってしまうため、1エポック（学習において全ての訓練データを使い切った時の回数。今回は60,000件の訓練データに対してミニバッチサイズが100なので、 60,000 / 100 = 600回）ごとに認識精度を計算してその値を結果配列に格納し、コンソールにも表示します。

認識精度推移のグラフ描画

　最後に認識精度の推移をグラフに描画します。

x = (0..(train_acc_list.size - 1)).to_a
Numo.gnuplot do
  plot x, train_acc_list, { w: :lines, t: 'train acc', lc_rgb: 'blue' },
       x, test_acc_list, { w: :lines, t: 'test acc', lc_rgb: 'green' }
  set xlabel: 'epochs'
  set ylabel: 'accuracy'
  set yrange: 0..1
end

　グラフの描画にはgnuplotを使い、Rubyでgnuplotを扱うためのgemとして、Numo::Gnuplot を使わせていただきました。

github.com

　こちらも参考にさせていただきました。

MF / 【Ruby】numo-gnuplotで遊ぶ

　以前に Jupyter Notebook 上で Nyaplot を使わせてもらっていたことはあるのですが、ターミナル上で動かして手軽にグラフを表示したいというときは Numo::Gnuplot が手軽に使えて良さそうです。

次は誤差逆伝播法

　冒頭でも書きましたが、今回実装はしたものの、繰り返しの一回の処理でもかなり時間がかかり、10,000回の繰り返しは現実的な時間では終わりそうもありませんでした。書籍でも言われていますが、数値微分による実装は実装はシンプルなものの処理にはかなり時間がかかるので、次は誤差逆伝播法での実装で試してみたいと思います。

　今回のコードは下記リポジトリで公開してあります。

github.com

ニューラルネットワークの推論処理

　前回に引き続き書籍「ゼロから作るDeepLearning」をベースに、前回NArray配列として扱うようにしたMNISTデータに対して推論処理を行うニューラルネットワークを実装してみます。

www.oreilly.co.jp

　ニューロンの構成は下記の通りです。

• 入力層: 784 # 28 x 28 の画像データのピクセル数
• 隠れ層1: 50
• 隠れ層2: 100
• 出力層: 10 # 結果として10クラス（0から9の数字）に分類する

サンプルコード全体

　まずはサンプルコードの全体を掲載しておきます。

require 'numo/narray'
require 'json'
require './sigmoid.rb'
require './softmax.rb'
require './mnist.rb'

def get_data
  x_train, t_train, x_test, t_test = load_mnist(true, true, false)
  [x_test, t_test]
end

def init_network
  nw = JSON.load(File.read('sample_weight.json'))
  network = {}
  nw.each do |k, v|
    network[k.to_sym] = Numo::DFloat[*v]
  end
  network
end

def predict(network, x)
  w1 = network[:w1]
  w2 = network[:w2]
  w3 = network[:w3]
  b1 = network[:b1]
  b2 = network[:b2]
  b3 = network[:b3]

  a1 = x.dot(w1) + b1
  z1 = sigmoid(a1)
  a2 = z1.dot(w2) + b2
  z2 = sigmoid(a2)
  a3 = z2.dot(w3) + b3
  softmax(a3)
end

x, t = get_data
network = init_network

batch_size = 100
accuracy_cnt = 0
x.to_a.each_slice(batch_size).with_index do |x_batch, idx|
  y_batch = predict(network, Numo::DFloat[*x_batch])
  p = y_batch.max_index(1) % 10
  accuracy_cnt += p.eq(t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]).cast_to(Numo::UInt8).sum
end

puts "Accuracy: #{accuracy_cnt.to_f / x.shape[0]}"

MNISTデータの取得

　前回実装したMNISTデータのロード処理（mnist.rb）を使ってMNISTデータのテストデータを取得します。

def get_data
  x_train, t_train, x_test, t_test = load_mnist(true, true, false)
  [x_test, t_test]
end

重みパラメータのロード

　学習済みの重みとバイアスのパラメータは、書籍のサンプルコードで pickle ファイルとして提供されているものをあらかじめJSONに変換してファイルに保存しておき（sample_weight.json）、それを読み込んでいます。

def init_network
  nw = JSON.load(File.read('sample_weight.json'))
  network = {}
  nw.each do |k, v|
    network[k.to_sym] = Numo::DFloat[*v]
  end
  network
end

推論処理

　上記でロードしたテストデータと重みパラメータに対して、隠れ層での活性化関数にはシグモイド関数、出力層での活性化関数にはソフトマックス関数を使って推論処理を行います。シグモイド関数は sigmoid.rb、ソフトマックス関数は softmax.rb として保存して読み込んでおきます。

def sigmoid(x)
  1 / (1 + Numo::DFloat::Math.exp(-x))
end

def softmax(a)
  c = a.max
  exp_a = Numo::DFloat::Math.exp(a - c)
  sum_exp_a = exp_a.sum
  exp_a / sum_exp_a
end

def predict(network, x)
  w1 = network[:w1]
  w2 = network[:w2]
  w3 = network[:w3]
  b1 = network[:b1]
  b2 = network[:b2]
  b3 = network[:b3]

  a1 = x.dot(w1) + b1
  z1 = sigmoid(a1)
  a2 = z1.dot(w2) + b2
  z2 = sigmoid(a2)
  a3 = z2.dot(w3) + b3
  softmax(a3)
end

処理の実行

　上記メソッドを使って処理を実行します。

x, t = get_data
network = init_network

batch_size = 100
accuracy_cnt = 0
x.to_a.each_slice(batch_size).with_index do |x_batch, idx|
  y_batch = predict(network, Numo::DFloat[*x_batch])
  p = y_batch.max_index(1) % 10
  accuracy_cnt += p.eq(t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]).cast_to(Numo::UInt8).sum
end

puts "Accuracy: #{accuracy_cnt.to_f / x.shape[0]}"

　get_data で取得した画像データを100件ずつバッチ処理します。NArrayの二次元配列データはそのままループすると一次元配列として並べて各要素が参照されてしまうので、to_a で通常の配列データに変換した上で、 each_slice で100件ずつのまとまりにし、with_index でインデックスを取得します。

x.to_a.each_slice(batch_size).with_index do |x_batch, idx|
  ...
end

　ループの中では100件分の画像データを再度NArray配列に変換してpredictメソッドに渡し、推論処理を実行した結果を y_batch として受け取っています。predict では100件分の各画像データについて、0-9の各数字に対しての確度が確率として返されるので、100 x 10 の二次元配列になります。

y_batch = predict(network, Numo::DFloat[*x_batch])

　その y_batch に対して max_index メソッドで二次元目の配列データについて最も数値が大きい要素のインデックスを取得しています。ただし取得されるインデックスは 100 x 10 を一次元配列として並べた場合のインデックスになるので、10で割った余りを取得することで、0-9の分類に変換しています。

Class: Numo::Int32 — Documentation by YARD 0.9.8

p = y_batch.max_index(1) % 10

　処理している画像データに対応するラベルデータを取得します。each_sliceとwith_indexを使った場合、インデックスとしてはeach_sliceのまとまりごとに0から順番に振られるので、インデックスにバッチサイズ（100）をかけることで配列データの始点を特定し、それにさらにバッチサイズを加算したものから1を引くことで、終点を特定して、ラベルデータ配列 t から該当するラベルデータを取得します。

t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]

　それを分類結果データ p と eq メソッドで比較します。

Class: Numo::Int32 — Documentation by YARD 0.9.8

　eqメソッドでは配列の各要素を比較し、一致する場合は1を、一致しない場合は0を配列として返します。

irb(main):043:0* a
=> Numo::Int32#shape=[5]
[1, 3, 5, 7, 9]
irb(main):044:0> b
=> Numo::Int32#shape=[5]
[1, 2, 5, 7, 8]
irb(main):045:0> a.eq(b)
=> Numo::Bit#shape=[5]
[1, 0, 1, 1, 0]

　返される配列のデータ型はNumo::BitなのでNumo::UInt8に変換し、合計値を取得することで推論が正しかった要素の数が取得できます。

accuracy_cnt += p.eq(t[(idx * batch_size)..(idx * batch_size + (batch_size - 1))]).cast_to(Numo::UInt8).sum

　そして最後に正解数を要素数で割ることで、正解率を計算しています。

puts "Accuracy: #{accuracy_cnt.to_f / x.shape[0]}"

　今回実装したコードはこちらにも公開しました。

github.com

MNISTデータセットをNArray配列として扱う

　RubyからMNISTのデータセットを扱えるように、書籍「ゼロから作るDeepLearning」のコードをベースにMNISTのデータセットをNArray配列としてロードするコードを書いてみました。

www.oreilly.co.jp

コード全体

　まずはコード全体を掲載します。

require 'open-uri'
require 'zlib'
require 'fileutils'
require 'numo/narray'

URL_BASE = 'http://yann.lecun.com/exdb/mnist/'
KEY_FILES = {
  train_img:   'train-images-idx3-ubyte.gz',
  train_label: 'train-labels-idx1-ubyte.gz',
  test_img:    't10k-images-idx3-ubyte.gz',
  test_label:   't10k-labels-idx1-ubyte.gz'
}

DATASET_DIR = "#{File.dirname(__FILE__)}/dataset"
SAVE_FILE = "#{DATASET_DIR}/mnist.dump"

IMG_SIZE = 784

def download(file_name)
  puts "Downloading #{file_name} ..."
  open(URL_BASE + file_name) do |file|
    open("#{DATASET_DIR}/#{file_name}", 'w+b') do |out|
      out.write(file.read)
    end
  end
  puts "Done."
end

def download_mnist
  FileUtils.mkdir_p(DATASET_DIR)
  KEY_FILES.each do |k, file|
    download(file)
  end
end

def load_img(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = gz.each_byte.to_a[16..-1].each_slice(IMG_SIZE).to_a
    data = Numo::UInt8[*data]
  end
  puts "Done"
  data
end

def load_label(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = Numo::UInt8[*gz.each_byte.to_a[8..-1]]
  end
  puts "Done"
  data
end

def convert_narray
  dataset = {}
  dataset[:train_img] = load_img(KEY_FILES[:train_img])
  dataset[:train_label] = load_label(KEY_FILES[:train_label])
  dataset[:test_img] = load_img(KEY_FILES[:test_img])
  dataset[:test_label] = load_label(KEY_FILES[:test_label])
  dataset
end

def init_mnist
  download_mnist
  dataset = convert_narray
  puts "Creating dump file ..."
  File.write(SAVE_FILE, Marshal.dump(dataset))
  puts "Done!"
end

def change_one_hot_label(x)
  one_hot_arrays = x.to_a.map do |v|
    one_hot_array = Array.new(10, 0)
    one_hot_array[v] = 1
    one_hot_array
  end
  Numo::UInt8[*one_hot_arrays]
end

def load_mnist(normalize = true, flatten = true, one_hot_label = false)
  unless File.exist?(SAVE_FILE)
    init_mnist
  end

  dataset = Marshal.load(File.read(SAVE_FILE))

  if normalize
    %i(train_img test_img).each do |key|
      dataset[key] = dataset[key].cast_to(Numo::DFloat)
      dataset[key] /= 255.0
    end
  end

  if one_hot_label
    dataset[:train_label] = change_one_hot_label(dataset[:train_label])
    dataset[:test_label] = change_one_hot_label(dataset[:test_label])
  end

  unless flatten
    %i(train_img test_img).each do |key|
      dataset[key] = dataset[key].reshape(dataset[key].shape[0], 28, 28)
    end
  end

  [dataset[:train_img], dataset[:train_label], dataset[:test_img], dataset[:test_label]]
end

データのダウンロード

　初回実行時はMNISTデータがホスティングされている Yann LeCun’s MNIST からデータをダウンロードします。対象ファイルは下記のように定義しています。

URL_BASE = 'http://yann.lecun.com/exdb/mnist/'
KEY_FILES = {
  train_img:   'train-images-idx3-ubyte.gz',
  train_label: 'train-labels-idx1-ubyte.gz',
  test_img:    't10k-images-idx3-ubyte.gz',
  test_label:   't10k-labels-idx1-ubyte.gz'
}

　これをデータセット用のディレクトリにダウンロードします。

def download(file_name)
  puts "Downloading #{file_name} ..."
  open(URL_BASE + file_name) do |file|
    open("#{DATASET_DIR}/#{file_name}", 'w+b') do |out|
      out.write(file.read)
    end
  end
  puts "Done."
end

def download_mnist
  FileUtils.mkdir_p(DATASET_DIR)
  KEY_FILES.each do |k, file|
    download(file)
  end
end

画像データのロード

　ダウンロードしたMNISTデータセットの画像データをロードします。データ形式は一般的な画像フォーマットではなく、Yann LeCun’s MNISTに下記のように説明があります。

TRAINING SET IMAGE FILE (train-images-idx3-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000803(2051) magic number 
0004     32 bit integer  60000            number of images 
0008     32 bit integer  28               number of rows 
0012     32 bit integer  28               number of columns 
0016     unsigned byte   ??               pixel 
0017     unsigned byte   ??               pixel 
........ 
xxxx     unsigned byte   ??               pixel
Pixels are organized row-wise. Pixel values are 0 to 255. 0 means background (white), 255 means foreground (black).

TEST SET IMAGE FILE (t10k-images-idx3-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000803(2051) magic number 
0004     32 bit integer  10000            number of images 
0008     32 bit integer  28               number of rows 
0012     32 bit integer  28               number of columns 
0016     unsigned byte   ??               pixel 
0017     unsigned byte   ??               pixel 
........ 
xxxx     unsigned byte   ??               pixel

　また、下記サイトも参考にさせていただきました。

TensorFlow : MNIST データ・ダウンロード (コード解説) – TensorFlow

　上記内容と書籍のコードを元に、下記のように実装しました。

def load_img(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = gz.each_byte.to_a[16..-1].each_slice(IMG_SIZE).to_a
    data = Numo::UInt8[*data]
  end
  puts "Done"
  data
end

　ダウンロードしたgzファイルをopenし、byte単位で読み込みます。先頭16byteは画像以外の情報なので、16byte以降から末尾までを読み込み、画像サイズ（28 x 28 = 784）ごとの配列で分割して二次元配列にして、それをNArrayの配列データに変換します。

ラベルデータのロード

　ラベルデータのフォーマットは下記の通りです。

TRAINING SET LABEL FILE (train-labels-idx1-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000801(2049) magic number (MSB first) 
0004     32 bit integer  60000            number of items 
0008     unsigned byte   ??               label 
0009     unsigned byte   ??               label 
........ 
xxxx     unsigned byte   ??               label
The labels values are 0 to 9.

TEST SET LABEL FILE (t10k-labels-idx1-ubyte):

[offset] [type]          [value]          [description] 
0000     32 bit integer  0x00000801(2049) magic number (MSB first) 
0004     32 bit integer  10000            number of items 
0008     unsigned byte   ??               label 
0009     unsigned byte   ??               label 
........ 
xxxx     unsigned byte   ??               label
The labels values are 0 to 9.

　実装は下記の通りです。基本的には画像データの場合と同じですが、画像データと違って正解ラベルの一次配列データなので、画像サイズで分割するような処理はありません。

def load_label(file_name)
  puts "Converting #{file_name} to NArray ..."
  data = nil
  Zlib::GzipReader.open("#{DATASET_DIR}/#{file_name}") do |gz|
    data = Numo::UInt8[*gz.each_byte.to_a[8..-1]]
  end
  puts "Done"
  data
end

MNISTデータのセーブとロード

　毎回MNISTデータをサイトからダウンロードしたらgzファイルからロードするのは時間がかかるので、初回実行時にデータをローカルに保持するようにします。今回はひとまず汎用性はあまり考慮せず、このスクリプトからのみ扱う前提で、 Marshal.dump したものをファイルに保存します。

File.write(SAVE_FILE, Marshal.dump(dataset))

　次回実行時はdumpファイルが存在すればダウンロードや変換処理はスキップして、dumpファイルをロードします。

dataset = Marshal.load(File.read(SAVE_FILE))

オプション

　基本的な処理はここまでですが、書籍の内容に沿って、3つのオプションを用意します。

　まず normalize オプションは画像データの各ピクセルデータを正規化するかどうかの指定で、元のデータは 0 - 255 の範囲の数値ですが、true を指定するとこれを 0.0 - 1.0 の値に正規化します。

if normalize
  %i(train_img test_img).each do |key|
    dataset[key] = dataset[key].cast_to(Numo::DFloat)
    dataset[key] /= 255.0
  end
end

　次に one_hot_label オプションは、正解ラベルデータを one_hot 表現の配列データとして返すかどうかの指定です。元データは画像データがどの数字なのかを表す一次元配列で、例えば画像データが 3, 5, 2 であれば、[3, 5, 2] という形ですが、 one_hot 表現だと、それぞれについて要素数10の配列を用意し、正解のインデックスのみ1で、それ以外は0の二次元配列になります。

[
 [0, 0, 0, 1, 0, 0, 0, 0, 0, 0], # 3が正解なのでインデックス3の要素だけ1
 [0, 0, 0, 0, 0, 1, 0, 0, 0, 0], # 5が正解なのでインデックス5の要素だけ1
 [0, 0, 1, 0, 0, 0, 0, 0, 0, 0]  # 2が正解なのでインデックス2の要素だけ1
]

def change_one_hot_label(x)
  one_hot_arrays = x.to_a.map do |v|
    one_hot_array = Array.new(10, 0)
    one_hot_array[v] = 1
    one_hot_array
  end
  Numo::UInt8[*one_hot_arrays]
end

if one_hot_label
  dataset[:train_label] = change_one_hot_label(dataset[:train_label])
  dataset[:test_label] = change_one_hot_label(dataset[:test_label])
end

　最後に flatten オプションは各画像のデータを 28 x 28 の二次元配列として返すか、784要素の一次元配列として返すかの指定になります。元データは784要素の一次元配列なので、flatten オプションに false が指定されれば 28 x 28 の二次元配列に変換します。

unless flatten
  %i(train_img test_img).each do |key|
    dataset[key] = dataset[key].reshape(dataset[key].shape[0], 28, 28)
  end
end

ロードしたMNISTデータを参照するサンプル

　上記コードによって読み込んだMNISTデータを参照してみます。上記のロード用コードは mnist.rb として保存しています。

require 'rmagick'
require './mnist.rb'

x_train, t_train, x_test, t_test = load_mnist

img = x_train[0, true]
label = t_train[0]
puts img.shape
puts img.max
puts img.min
puts label

image = Magick::Image.new(28, 28)
image.import_pixels(0, 0, 28, 28, 'I', img, Magick::FloatPixel)
image.display

　画像データとラベルデータからそれぞれ一つ目を取り出し、ラベルはターミナルに表示しています。画像データは ImageMagick(RMagick)で 28 x 28 の画像オブジェクトを作成した上で、 import_pixelsメソッドでインポートして、displayメソッドで別ウィンドウに表示します。先頭データは5なので、5の画像が表示されます。

　今回実装したコードは下記にも公開しています。

github.com