不確定特異点

広く深く、ところどころ超深く

STM32F7 Discovery の開発環境 (2)

LED チカチカプログラム

まずは、お決まりの LED を点滅させるプログラムを作ります。回路図を見ると、STM32F456 の PI1 端子は LD1 という緑色の LED に接続されてます。PI1 は GPIOI で制御可能なので、GPIOI を 0/1 制御することで LED を点滅させてみます。

#include "stm32f7xx.h"

#define BIT(x, n) ((x) & (1 << (n)))

#define COUNTER_BIT 17

void Reset_Handler()
{
    int i = 0;

    SystemInit();
    
    /* Enable clock for each peripherals */
    RCC->AHB1ENR |= RCC_AHB1ENR_GPIOIEN;
    
    /* PI1 pin is assigned for GPIO output */
    GPIOI->MODER = 0x00000004;
    
    while (1) {
        if (i++ == (1<<COUNTER_BIT) - 1)
            i = 0;
        GPIOI->ODR = BIT(i,COUNTER_BIT-1) >> (COUNTER_BIT-2);
    }
}

リセットハンドラに実行後のミニマムの設定は、

  • RCC の AHB1ENR レジスタにより GPIOI に対してクロックを供給する
  • GPIOI の MODER レジスタにより出力設定にする

の 2 点です。あとは、カウンタ(i)の上位ビットを GPIOI の ODR レジスタから出力させているだけです。カウンタを 0 クリアするタイミングは、点滅が目視できる値に適当に調整します。

ROM のプログラム方法

Flash が 0x08000000 番地にマッピングされているため、この番地を ROM 領域として使うようにリンカスクリプトを設定して ELF ファイルを作成します。ちなみに Flash は 0x00200000 番地からの ITCM 経由と 0x08000000 番地からの AXIM 経由の 2 通りのバスによるマッピングがなされており、リンカスクリプトで ROM の開始アドレスを設定することにより、どちらのバスを使ってアクセスするかを選択できます。両方試してみましたが、ITCM の方が高速でした。命令フェッチにバスを占有できるからでしょうか。(謎な ART とかいうのもあったけど)

参考

image.elf という ELF ファイルを OpenOCD で転送するには次のように実行します。このとき、ROM を書き込むのはあくまで 0x08000000 番地を指定する、ということのようです。

$ arm-linux-gnueabi-objcopy -O binary image.elf image
$ openocd -s /usr/local/share/openocd/scripts/ -f board/stm32f7discovery.cfg -c "program image verify reset exit 0x08000000"
Open On-Chip Debugger 0.10.0-dev-00026-g97a8b08-dirty (2015-10-23-18:37)
Licensed under GNU GPL v2
For bug reports, read
        http://openocd.org/doc/doxygen/bugs.html
Info : The selected transport took over low-level target control. The results might differ compared to plain JTAG/SWD
adapter speed: 2000 kHz
adapter_nsrst_delay: 100
srst_only separate srst_nogate srst_open_drain connect_deassert_srst
Info : Unable to match requested speed 2000 kHz, using 1800 kHz
Info : Unable to match requested speed 2000 kHz, using 1800 kHz
Info : clock speed 1800 kHz
Info : STLINK v2 JTAG v24 API v2 SWIM v11 VID 0x0483 PID 0x374B
Info : using stlink api v2
Info : Target voltage: 3.211002
Info : stm32f7x.cpu: hardware has 8 breakpoints, 4 watchpoints
target state: halted
target halted due to debug-request, current mode: Thread
xPSR: 0x01000000 pc: 0x0800025c msp: 0x20020000
** Programming Started **
auto erase enabled
Info : device id = 0x10016449
Info : flash size = 1024kbytes
wrote 32768 bytes from file image in 1.673457s (19.122 KiB/s)
** Programming Finished **
** Verify Started **
verified 740 bytes in 0.224078s (3.225 KiB/s)
** Verified OK **
** Resetting Target **
shutdown command invoked

LED がチカチカしました。

PLL を使って 216 MHz までクロックアップ

正直、STM32CubeF7 のドライバを使った方が手早くて確実なのですが、スタートアップに必要な処理を理解するために自力で設定してみます。STM32F4 Discovery の開発環境 (3) で実装した EnablePLL() のコードを応用して、STM32F7 対応させてみます。変更点はそれほど多くはなく、

  • PLL の逓倍設定
    • HSE=16MHz ⇒ 168MHz 設定から HSE=25MHz ⇒ 216MHz 設定に変更
  • オーバードライブモードへの移行を追加
    • STM32F4 Discovery の時は電源系は Voltage Scaling = 1 のみしか設定しなかった。
  • Flash レイテンシの設定を変更
    • 5WS ⇒ 7WS へ変更

を追加変更しました。(※以下のプログラムでは端折ってますが、レジスタ設定はレジスタへの Write が反映されていることを確認するため、Read チェックを入れた方が確実かもしれません。)

例によって、EnablePLL() により所望の周波数設定になっていることを確認するため、タイマ割り込みを使って GPIOI 端子をトグルさせて周期を測定してみます。

#include "stm32f7xx.h"

#define PLL_M 25
#define PLL_N 432
#define PLL_P 2
#define PLL_Q 9

void EnablePLL()
{
    RCC->CR |= RCC_CR_HSEON;
    while (!(RCC->CR & RCC_CR_HSERDY))
        continue;
    
    /* Power supply setup */
    RCC->APB1ENR |= RCC_APB1ENR_PWREN;
    PWR->CR1 |= PWR_CR1_VOS;
    
    /* Enable Over-drive mode */
    PWR->CR1 |= PWR_CR1_ODEN;
    while (!(PWR->CSR1 & PWR_CSR1_ODRDY))
        continue;

    /* Switch to Over-drive mode */
    PWR->CR1 |= PWR_CR1_ODSWEN;
    while (!(PWR->CSR1 & PWR_CSR1_ODSWRDY))
        continue;

    /*
     * PLL settings
     * AHB  = SYSCLK / 1 = 216MHz
     * APB1 = SYSCLK / 4 = 54MHz
     * APB2 = SYSCLK / 2 = 108MHz
     */
    RCC->CFGR |= RCC_CFGR_HPRE_DIV1  |
                 RCC_CFGR_PPRE2_DIV2 |
                 RCC_CFGR_PPRE1_DIV4;
    
    /* PLLCLK = HSI(16 MHz) / M * N / P */
    RCC->PLLCFGR = PLL_M | (PLL_N << 6) | (((PLL_P >> 1) -1) << 16) |
                   RCC_PLLCFGR_PLLSRC_HSE | (PLL_Q << 24);

    /* PLL ON */
    RCC->CR |= RCC_CR_PLLON;

    /* Wait until PLL is locked up */
    while (!(RCC->CR & RCC_CR_PLLRDY))
        continue;
    
    /* Increase flash latency */
    FLASH->ACR = FLASH_ACR_LATENCY_7WS;

    RCC->CFGR &= ~RCC_CFGR_SW;
    RCC->CFGR |= RCC_CFGR_SW_PLL;
}

void Reset_Handler()
{
    SystemInit();
    
    /* Enable clock for each peripherals */
    RCC->AHB1ENR |= RCC_AHB1ENR_GPIOIEN;
    RCC->APB1ENR |= RCC_APB1ENR_TIM2EN;
    
    /* PI1 pin is assigned for GPIO output */
    GPIOI->MODER = 0x00000004;
    
    EnablePLL();
    
    /* Timer interrupting at regular interval */
    TIM2->PSC  = 0;            /* Timer clock  108 MHz / PSC = 104 KHz */
    TIM2->ARR  = 54*2*1000-1;  /* Auto reload: 108 KHz / ARR = 1 Hz */
    TIM2->DIER = TIM_DIER_UIE; /* Enable update interrupt */
    TIM2->EGR  = TIM_EGR_UG;   /* Initialize counter */
    TIM2->SR   = 0;            /* Clear all interrupts */
    TIM2->CR1 |= TIM_CR1_CEN;  /* Enable counter */

    NVIC_SetPriority(TIM2_IRQn, 2);
    NVIC_EnableIRQ(TIM2_IRQn);

    while (1)
        continue;
}

void TIM2_IRQHandler()
{
    if (TIM2->SR & TIM_SR_UIF) {
        TIM2->SR = ~TIM_SR_UIF;   /* Clear update interrupt */
        GPIOI->ODR ^= 0x00000002; /* Toggle GPIO output */
    }
}

108 MHz(SYSCLK/2) をタイマのオートリロード機能で分周して 1 KHz まで落とし、その 1 KHz 周期でかかる割り込み処理の中で GPIOI をトグルさせるので、PI1 端子上に 500 Hz が観測されるはずです。というわけで、GPIOI の出力をオシロで測定してみます。

f:id:tanakahx:20151026003750j:plainf:id:tanakahx:20151026011858p:plain

500 Hz ピッタリです。

STM32F7 Discovery の開発環境 (1)

STM32F7 Discovery という Cortex-M7 コア搭載の基板を入手したので、早速使ってみました。SRAM 320KB, SDRAM 16MB, 静電容量式タッチパネル付き TFT 液晶, USB, Ethernet, Micro SD, Audio 果ては カメラ I/F まで搭載し、Arduino も被せられるようになっている何でもあり基板になってます。これで 8000 円なのでかなりコスパがいい方だと思います。

f:id:tanakahx:20151025070035j:plainf:id:tanakahx:20151025070036j:plain

ST-LINK/V2-1 ファームウェアのアップデート

ST-LINK/V2-1 のファームウェアをアップデートすることが推奨されているようなので、アップデータをダウンロードして実行しました。Windows では EXE ファイルが提供されてますが、Linux/Mac OS向けには Java のプログラムが提供されていました。Java 版は以下のように起動します。

$ cd stsw-link007/AllPlatforms
$ java -jar STLinkUpgrade.jar

f:id:tanakahx:20151025062339p:plain

"Refresh device list"、"Open in update mode"、"Upgrade" の順に実行するとファームウェアの更新が始まるので、成功を祈りながら心静かに見守ります。

仮想マシンに ST-LINK の USB 接続を認識させる

僕はホスト OS に Mac OS を使ってますが、ARM の開発環境は VirtualBox 上の Ubuntu に準備しているため、STMF7 Discovery を USB で接続した場合に、Mac OS ではなく Ubuntu 側に認識させるための設定を行います。/etc/udev/rules.d/49-stlinkv2.rules のルールファイルに以下の行を追加します。

SUBSYSTEMS=="usb", ATTRS{idVendor}=="0483", ATTRS{idProduct}=="374b", \
    MODE:="0666", \
    SYMLINK+="stlinkv2_%n"

ルールファイルを反映するには以下のコマンドを実行します。

$ sudo udevadm control --reload-rules

下図のように VirtualBox 側にもゲスト OS に USB を見せられるよう、デバイスフィルタの設定をしておきます。udev のルールファイルに設定したベンダ ID やプロダクト ID も VirtualBox のデバイスフィルタから確認することができます。

f:id:tanakahx:20151025063014p:plain

OpenOCD のインストール

Linux から STM 基板にプログラムを書き込むツールとして OpenOCD を使います。現時点では STM32F7 Discovery に正式対応されていないようなので、以下のように最新版に対してパッチを当てる必要があります。(なお、これを試した時点での最新版ハッシュ値は a859befa67c4abfc8aafb95f9c9c7bc3629f1ff8 です。)

$ git clone http://openocd.zylin.com/openocd
$ git fetch http://openocd.zylin.com/openocd refs/changes/54/2754/5 && git checkout FETCH_HEAD

インストール手順は以下の通りです。

$ ./bootstrap
$ ./configure
$ make
$ sudo make install

OpenOCD を実行して STM 基板の ST-LINK 用の LED が点滅すれば接続に成功してます。

$ openocd -s /usr/local/share/openocd/scripts/ -f board/stm32f7discovery.cfg
Open On-Chip Debugger 0.10.0-dev-00026-g97a8b08-dirty (2015-10-23-18:37)
Licensed under GNU GPL v2
For bug reports, read
        http://openocd.org/doc/doxygen/bugs.html
Info : The selected transport took over low-level target control. The results might differ compared to plain JTAG/SWD
adapter speed: 2000 kHz
adapter_nsrst_delay: 100
srst_only separate srst_nogate srst_open_drain connect_deassert_srst
Info : Unable to match requested speed 2000 kHz, using 1800 kHz
Info : Unable to match requested speed 2000 kHz, using 1800 kHz
Info : clock speed 1800 kHz
Info : STLINK v2 JTAG v24 API v2 SWIM v11 VID 0x0483 PID 0x374B
Info : using stlink api v2
Info : Target voltage: 3.211002
Info : stm32f7x.cpu: hardware has 8 breakpoints, 4 watchpoints

次はテスト用のプログラムを作成して Flash ROM に焼いてみます。

LTK でウィンドウを手前に表示する

LTK とは Tk ツールキットを Common Lisp から使えるようにしたものです。Tcl/Tk をベースにしているため、幅広い OS に対応していることや手軽に UI を作成できることが魅力です。Common Lisp でちょっとした GUI を作りたいときにも便利そうですね。

LTK を使うには Tcl/Tk 付属の wish コマンドが使えればよく、Common Lisp ライブラリのインストールは Quicklisp から LTK を ql:quickload するだけですみます。

CL-USER> (ql:quickload :ltk)
To load "ltk":
  Load 1 ASDF system:
    ltk
; Loading "ltk"

(:LTK)

LTK には ltktest というサンプルプログラムが付属しており、以下のように試すことができます。また、ltktest のように export はされてませんが、ソースコードの中には ltk-eyes という xeyes 的なサンプルプログラムもありました。

CL-USER> (ltk:ltktest)
CL-USER> (ltk::ltk-eyes)

ところで上記のサンプルを動かしてみるとわかるのですが、LTK を実行して起動したウィンドウが手前に表示されないため、インタラクティブに UI を組み立てていくことが面倒です。これは別に LTK に限った話ではなく、wish でも同様の動作になるのですが、なんとかウィンドウを手前に持ってきたいので調べてみると、wish では次のようにすれば良いことがわかりました。

wm attribute . -topmost 1

しかし、LTK のドキュメントを見る限り、これに相当する直接的な API は存在しないようでした。そこで、LTK のコードを眺めてみると、LTK は実際には内部的に wish を起動し、format-wish という関数により、ストリーム経由で Tcl コマンドを wish に送ることで GUI を作る、という仕組みのようでした(たぶん)。というわけで、format-wish を使って上記の Tcl コマンドを送ることで LTK でも同じ効果が得られないか、実験してみました。

(defun main ()
  (ltk:with-ltk ()
    (let* ((b (make-instance 'ltk:button
                             :master nil
                             :text "hello"
                             :command #'(lambda () (format t "hello, world!~%"))))
           (c (ltk:make-canvas nil :width 200 :height 150))
           (r (ltk:create-rectangle c 20 20 100 100)))
      (ltk:itemconfigure c r "fill" "black")
      (ltk:pack c)
      (ltk:pack b)
      (ltk:format-wish "wm attribute . -topmost 1"))))

単なるサンプルとして UI を並べているだけなのでプログラム的にはあまり意味はありません(キャンバスとボタンがあって、キャンバスに四角をお絵描きするだけです)が、最終行にあるように format-wish により wm attribute . -topmost 1 を実行してあげると、期待通りウィンドウが手前に表示されました!

f:id:tanakahx:20151020234953p:plain

format-wish を使えば任意の Tcl コマンドを wish に投げられるので、これ以外にも、LTK の API として提供されていないような細かな設定が色々できると思います。

以下の記事を参考にさせて頂きました。

Ubuntu 14.04.3 LTS に Chainer をインストールする

Deep learning のフレームワークである Chainer を試用してみたところ、サンプルの MNIST の学習に時間がかかりすぎて、これは CPU に計算させる代物ではないことを悟りました。そこで、GPU を使って計算するための環境を整えることにしました。

環境

Chainer のバージョン

本記事の執筆時のバージョンは v1.4.0 です。(※2015/11/28追記)

Ubuntu 14.04.3 LTS のインストール

UNetbootin で USB に Ubuntu の iso ファイルを書き、USB ブートした上でインストールしました。SSDパーティション構成は以下のようにしました。

パーティション マウントポイント サイズ タイプ
/dev/sda1 /boot/efi 128MB efi
/dev/sda2 / 残りすべて ext4
/dev/sda3 - 2GB swap

BIOS のブートオプションでこの SSD を指定すると、/boot/efi から GRUB が起動されます。

システムの更新

カーネル含めて一度すべて最新の状態にします。

$ sudo apt-get update
$ sudo apt-get dist-upgrade

更新後は問題なく再起動できることを確認します。

Nouveau ドライバの無効化

NVIDIAGPU を積んだグラフィックカードの場合、デフォルトで Nouveau というドライバが使われていると思います。

$ lsmod | grep -i nouveau

以降の手順でインストールされる NVIDIA 製のドライバと競合するらしいので、いったん無効化しておきます。 /etc/modprobe.d/blacklist-nouveau.conf というファイルを新規に作成し、以下の設定を記述します。

blacklist nouveau
options nouveau modeset=0

カーネルモジュールを blacklist に登録した際は以下を実行しておくらしいです。

$ sudo update-initramfs -u

ここで再起動し、ディスプレイ解像度が落ちることで Nouveau が無効化されたことを確認します。

CUDA のインストール

CUDA 7.5 Downloads のサイトから deb ファイルをダウンロードして、以下の手順でインストールします。CUDA と一緒に NVIDIA のドライバもインストールされます。

$ sudo dpkg -i cuda-repo-ubuntu1404_7.5-18_amd64.deb
$ sudo apt-get update
$ sudo apt-get install cuda

CUDA のインストールが完了したら、以下の設定を .bashrc 等に追加しておきます。CUDA 用のコンパイラ(nvcc)や CUDA ライブラリが使えるようになります。

export PATH=/usr/local/cuda/bin:$PATH
export LD_LIBRARY_PATH=/usr/local/cuda/lib64:$LD_LIBRARY_PATH

この時点で再起動すると、CUDA と共に NVIDIA のドライバがインストールされるため、ディスプレイの解像度が元に戻ります。

$ lsmod | grep nvidia
nvidia_uvm             77824  0
nvidia               8605696  35 nvidia_uvm
drm                   344064  3 nvidia

NVIDIA のドライバを確認できました。以上で CUDA のインストールは完了です。

サンプルプログラムの動作確認

CUDA が正しくインストールされたことを確認するために、サンプルプログラムを動かしてみます。下記のように適当なディレクトリにコピーして make します。

$ cp -r /usr/local/cuda/samples ~
$ cd ~/samples
$ make

deviceQuery というサンプルを実行した結果です。

$ bin/x86_64/linux/release/deviceQuery
/home/tanaka/samples/bin/x86_64/linux/release/deviceQuery Starting...

 CUDA Device Query (Runtime API) version (CUDART static linking)

Detected 1 CUDA Capable device(s)

Device 0: "GeForce GTX 960"
  CUDA Driver Version / Runtime Version          7.5 / 7.5
  CUDA Capability Major/Minor version number:    5.2
  Total amount of global memory:                 4094 MBytes (4292719616 bytes)
  ( 8) Multiprocessors, (128) CUDA Cores/MP:     1024 CUDA Cores
  GPU Max Clock rate:                            1278 MHz (1.28 GHz)
  Memory Clock rate:                             3505 Mhz
  Memory Bus Width:                              128-bit
  L2 Cache Size:                                 1048576 bytes
  Maximum Texture Dimension Size (x,y,z)         1D=(65536), 2D=(65536, 65536), 3D=(4096, 4096, 4096)
  Maximum Layered 1D Texture Size, (num) layers  1D=(16384), 2048 layers
  Maximum Layered 2D Texture Size, (num) layers  2D=(16384, 16384), 2048 layers
  Total amount of constant memory:               65536 bytes
  Total amount of shared memory per block:       49152 bytes
  Total number of registers available per block: 65536
  Warp size:                                     32
  Maximum number of threads per multiprocessor:  2048
  Maximum number of threads per block:           1024
  Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
  Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)
  Maximum memory pitch:                          2147483647 bytes
  Texture alignment:                             512 bytes
  Concurrent copy and kernel execution:          Yes with 2 copy engine(s)
  Run time limit on kernels:                     Yes
  Integrated GPU sharing Host Memory:            No
  Support host page-locked memory mapping:       Yes
  Alignment requirement for Surfaces:            Yes
  Device has ECC support:                        Disabled
  Device supports Unified Addressing (UVA):      Yes
  Device PCI Domain ID / Bus ID / location ID:   0 / 1 / 0
  Compute Mode:
     < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >

deviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 7.5, CUDA Runtime Version = 7.5, NumDevs = 1, Device0 = GeForce GTX 960
Result = PASS

python-dev, pip のインストール

まず、Python のパッケージから必要とされることが多い python-dev をインストールします。後述の NumPy のインストールにも必要です。

$ sudo apt-get install python-dev

続いて、他のパッケージと同様に Chainer は pip からインストール可能なので、pip をインストールします。

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

apt-get でも python-pip なるパッケージとしてインストール可能ですが、最新版にするため get-pip.py からインストールしました。Python 3.4, 2.7.9 以降では ensurepip という仕組みで最初から pip が使えるようになっているようです。面倒なので Ubuntu も早く ensurepip にして欲しいな。

NumPy, SciPy のインストール

SciPy で必要になるため、あらかじめ BLAS/LAPACK ライブラリと gfortran を入れておきます。

$ sudo apt-get install libblas-dev liblapack-dev gfortran

そして、pip を使って NumPy, SciPy をインストールします。ここで明示的に pip install しなくても、後述の「Chainer のインストール」の過程で自動的にインストールされます。(※2015/11/28追記)

$ pip install --user numpy scipy

なお、以下のように、apt-get でインストールするのも手っ取り早いとは思いますが、最新版とは限りません。

$ sudo apt-get install python-numpy python-scipy

Chainer のインストール

上記の NumPy や SciPy と同様に、pip で chainer-cuda-deps と chainer というパッケージをインストールします。chainer-cuda-deps は GPU 実行する場合に必要になるようです。

pip で chainer パッケージをインストールします。なお、本記事執筆時の Chainer v1.4.0 では chainer-cuda-deps は不要です。v1.3.0 から不要になったようです。(※2015/11/28修正。情報をご提供くださった奥田さん、ありがとうございます!)

$ # pip install --user chainer-cuda-deps # これは不要
$ pip install --user chainer

以上で、インストール作業は完了です。

では、早速 Chainer の動作確認として、MNIST の手書き文字認識の学習データを用いた学習のサンプルを実行してみます。

$ sudo apt-get install git
$ git clone https://github.com/pfnet/chainer.git
$ cd chainer/examples/mnist
$ python train_mnist.py -g 0
...
('epoch', 20)
train mean loss=0.0444635186407, accuracy=0.986683343351
test  mean loss=0.0679141613924, accuracy=0.983100009561

計算中に nvidia-smi を実行するとこんな感じ。GPU が使われていることが確認できます。

$ nvidia-smi
Mon Oct 12 23:36:04 2015
+------------------------------------------------------+
| NVIDIA-SMI 352.39     Driver Version: 352.39         |
|-------------------------------+----------------------+----------------------+
| 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 960     Off  | 0000:01:00.0      On |                  N/A |
|  0%   40C    P2    54W / 130W |    167MiB /  4093MiB |     48%      Default |
+-------------------------------+----------------------+----------------------+

+-----------------------------------------------------------------------------+
| Processes:                                                       GPU Memory |
|  GPU       PID  Type  Process name                               Usage      |
|=============================================================================|
|    0      1212    G   /usr/bin/X                                      39MiB |
|    0      5940    C   python                                         113MiB |
+-----------------------------------------------------------------------------+

GPU(GTX 960 4GB) 版で epoch 20 まで達するに要した時間は約 1 分 30 秒くらいでした。-g のオプションを省略すると CPU で計算するようになるのですが、単純計算で 60 倍ほど時間がかかるようでしたので、Deep learning で何かやろうと思ったら GPU は必須だと実感しました。

学習とニューラルネットワーク (5)

ボルツマンマシン

今回はボルツマンマシンのシミュレーションを実装してみました。

Simple Boltzmann Machine Simulation

エネルギー関数として { E(x_1, x_2, x_3) = (x_1 + x_2 + x_3 - 3)^{2} } を仮定すると、エネルギーの極小点は { x_1 = x_2 = x_3 = 1 } の 1 点になります。これを実際にシミュレーションにより再現できるかを見てみます。(数値は左から右にかけて状態 (0,0,0)〜(1,1,1) の出現確率に対応)

BM> (bm:bm (clnu.mx:mx t (0 -2 -2) (-2 0 -2) (-2 -2 0))
           (clnu.mx:mx t (-5) (-5) (-5))
           :update-count 1000
           :temperature 0.5)
理論値: .000 .000 .000 .096 .000 .096 .096 .711
実測値: .000 .000 .000 .088 .000 .094 .100 .718
NIL

bm という関数に与えている 2 つの行列は、エネルギー関数を展開したときの 2 次式と 1 次式に相当します。

x1, x2, x3 は確率的に 0 または 1 の値を取りますが、それらの 8 パタンの組み合わせに対する出現確率を理論値と実測値で求めてます。理論値 0.711 が最も高いですが、これは先の極小点の状態に対応しており、確率的な挙動の正しさが確認できます。

続いてエネルギー関数として { E(x_1, x_2, x_3) = (x_1 + x_2 + x_3 - 1)^{2} } を仮定すると、エネルギーの極小点は { (x_1,  x_2, x_3) = (1, 0, 0), (0, 1, 0), (0, 0, 1) } の 3 点になるので、出現確率も 3 つの極値が存在するはずです。

BM> (bm:bm (clnu.mx:mx t (0 -2 -2) (-2 0 -2) (-2 -2 0))
           (clnu.mx:mx t (-1) (-1) (-1))
           :update-count 1000
           :temperature 0.5)
理論値: .038 .282 .282 .038 .282 .038 .038 .000
実測値: .040 .214 .326 .029 .305 .039 .047 .000
NIL

というわけで、上記の通り理論値 0.282 の確率で生起する状態が 3 状態確認できました。

エネルギー関数が持つ温度パラメータについて考察してみます。温度が低いと状態遷移する頻度が低くなるため、エネルギー関数の極小点が平衡状態における出現確率として浮かび上がりやすいですが、平衡状態は初期状態に左右されやすいため局所解にとらわれる可能性があります。逆に、温度が高いと状態遷移する頻度も高くなるため、あまりに温度が高すぎると出現確率が均一になってしまい、エネルギー関数の極小点が見えづらくなる、ということがわかりました。

実装は以下の本を参考にしました。(ちなみにこの本のボルツマンマシンのサンプルコードは、未初期化の配列xにアクセスしているため正しく動作しません…)

学習とニューラルネットワーク (電子情報通信工学シリーズ)

学習とニューラルネットワーク (電子情報通信工学シリーズ)

Python (NumPy) と Common Lisp (LLA) で行列積の実行速度を比較する

近年、機械学習や Deep Learning などのデータサイエンス分野を筆頭に数値計算の需要が高まってますが、その中でよく使われているのが Python の NumPy というライブラリです。NumPy を使うことで、動的で柔軟なスクリプト言語上で比較的高速に数値計算が可能になり、NumPy を使った SciPy や matplotlib などのライブラリも充実していることから、近頃の Python 人気が高まっているのだと思います。

個人的には Common Lisp が好きなので、NumPy との速度差が気になるところです。そこで単精度浮動小数点数を要素とする100行100列の行列積を対象に、PythonCommon Lisp、C でベンチマークをとって比較してみました。(あ、ちなみに僕は Python も好きですよ。)

実行環境は以下の通りです。

$ uname -a
Linux vagrant-ubuntu-trusty-64 3.13.0-62-generic #102-Ubuntu SMP Tue Aug 11 14:29:36 UTC 2015 x86_64 x86_64 x86_64 GNU/Linux

Python (NumPy)

まずは Python (NumPy) です。バージョンは 2.7.6 です。

# coding: utf-8
import numpy as np

N = 100

def simple_gemm(ma, mb):
    """NumPyの行列計算を使わずに実装した版。とてつもなく遅い。"""
    mc = np.zeros((ma.shape[0], mb.shape[0]))
    for i in range(ma.shape[0]):
        for j in range(mb.shape[1]):
            for k in range(ma.shape[1]):
                mc[i][j] += ma[i][k] * mb[k][j]
    return mc

def numpy_gemm(ma, mb):
    """NumPyによる高速な行列計算(BLASライブラリを利用)"""
    return np.dot(ma, mb)

# N行N列の行列を生成して、要素を[0, 1.0]の範囲の乱数で初期化
ma = np.random.rand(N, N)
mb = np.random.rand(N, N)

def run_simple(n = 10000):
    for i in range(n):
        simple_gemm(ma, mb)

def run_gemm(n = 10000):
    for i in range(n):
        numpy_gemm(ma, mb)

simple_gemm は比較用に NumPy を使わずに実装したもので、numpy_gemm が NumPy の機能を使ったものです。ipython 上で実行時間を計測してみます。

In [1]: run test.py

In [2]: time run_simple (10)
CPU times: user 7.5 s, sys: 8.29 ms, total: 7.51 s
Wall time: 7.68 s

In [3]: time run_gemm (10000)
CPU times: user 10.3 s, sys: 12.3 ms, total: 10.3 s
Wall time: 10.6 s

わかりづらいですが、simple_gemm を 10 回、numpy_gemm を 10000 回実行してます。simple_gemm の方は遅すぎて 10 回しか計算しませんでした。一方、 numpy_gemm の方は 70 倍くらい高速化されてます。この高速化の秘密を追って NumPy のソースコードを調べてみると、BLAS ライブラリに行き着きました。

BLAS (Basic Linear Algebra Subprograms)

BLAS とはもともと Fortran で書かれた線形代数演算ルーチンライブラリのデファクトスタンダードであり、現在はオリジナルのリファレンス BLAS の他に Intel MKL, ATLAS, OpenBLAS(GotoBLAS) など色々あるようです。どのライブラリがインストールされていてリンクされるかは次のコマンドで確認できます。

$ update-alternatives --display libblas.so.3
libblas.so.3 - manual mode
  link currently points to /usr/lib/libblas/libblas.so.3
/usr/lib/libblas/libblas.so.3 - priority 10
  slave libblas.so.3gf: /usr/lib/libblas/libblas.so.3
/usr/lib/openblas-base/libblas.so.3 - priority 40
  slave libblas.so.3gf: /usr/lib/openblas-base/libblas.so.3
Current 'best' version is '/usr/lib/openblas-base/libblas.so.3'.

使うライブラリを変更するには、--display の代わりに --config を指定します。

$ update-alternatives --config libblas.so.3
There are 2 choices for the alternative libblas.so.3 (providing /usr/lib/libblas.so.3).

  Selection    Path                                 Priority   Status
------------------------------------------------------------
  0            /usr/lib/openblas-base/libblas.so.3   40        auto mode
* 1            /usr/lib/libblas/libblas.so.3         10        manual mode
  2            /usr/lib/openblas-base/libblas.so.3   40        manual mode

Press enter to keep the current choice[*], or type selection number:

この OS には OpenBLAS とリファレンス BLAS がインストールされていて、リファレンス BLAS が使われていることがわかりました。

BLASループアンローリングやキャッシュヒット率の向上など等 CPU の性能を最大限に発揮できるような最適化が施されているため、これに Common Lisp Only で立ち向かうのは大変そうです。というわけで Common Lisp でも BLAS を使えるライブラリがないか調べてみました。へたれですみません。

Common Lisp (LLA)

LLA というライブラリがありました。

github.com

BLASLAPACK といった有名な数値計算ライブラリを Common Lisp から使えるようにするためのラッパーライブラリであり、NumPy のようなリッチな機能は搭載されてませんが、シンプルな分、今回の目的には十分です。

LLA が使うライブラリはデフォルトで libblas.so が指定されてますが、別のパスを設定する場合は README にあるように、*lla-configuration* という CL-USER パッケージのシンボルにライブラリへのパスを .sbclrc に設定しておきます。

(defvar *lla-configuration*
        '(:libraries ("/usr/lib/libblas/libblas.so.3.0")))

通常は libblas.so から適当なライブラリをシンボリックリンクしているため、上記の設定は不要だと思いますが、システム非標準のパスにあるライブラリを指定する場合などに使うことになるのかもしれません。

以下のようなベンチマーク用プログラムを作りました。Common Lisp の処理系は SBCL 1.2.15 です。

(defparameter *N* 100)

(defun make-matrix (rows cols)
  (make-array (list rows cols) :initial-element 0.0 :element-type 'single-float))

(defun randomize (m)
  "行列の要素を[0, 1]の範囲の乱数で初期化する"
  (destructuring-bind (rows cols)
      (array-dimensions m)
    (dotimes (row rows)
      (dotimes (col cols)
        (setf (aref m row col) (random 1.0)))))
  m)

(defun simple-gemm (ma mb)
  "Common Lispのみを使った行列積の計算"
  (declare (optimize (speed 3) (debug 0) (safety 0)))
  (declare (type (simple-array single-float (* *)) ma mb))
  (let ((rows (array-dimension ma 0))
        (cols (array-dimension mb 1)))
    (declare (type fixnum rows cols))
    (let ((result (make-matrix rows cols)))
      (declare (type (simple-array single-float (* *)) result))
      (dotimes (row rows)
        (dotimes (col cols)
          (dotimes (k cols)
            (incf (aref result row col)
                  (* (aref ma row k) (aref mb k col))))))
      result)))

(defun run (gemm &optional (n 10000))
  "gemmをn回計算して実行時間を測定する。行列の初期化は計測しない。"
  (let ((ma (randomize (make-matrix *N* *N*)))
        (mb (randomize (make-matrix *N* *N*))))
    (time
     (loop repeat n
        do (funcall gemm ma mb)
        finally (return nil)))))

(defun run-simple-gemm (&optional (n 10000))
  "simple-gemmのベンチマーク"
  (run #'simple-gemm n))

(defun run-lla-mm (&optional (n 10000))
  "llaのmm関数のベンチマーク"
  (run #'lla:mm n))

Python版と同様に、BLAS を使わない simple-gemmBLAS を使った lla-mm との比較をします。単精度浮動小数点数を要素とする100行100列の行列積を 10000 回実行した場合の実行時間です。

CL-USER> (run-simple-gemm 10000)
Evaluation took:
  33.395 seconds of real time
  33.236077 seconds of total run time (33.148582 user, 0.087495 system)
  [ Run times consist of 0.032 seconds GC time, and 33.205 seconds non-GC time. ]
  99.52% CPU
  80,349,715,978 processor cycles
  401,280,592 bytes consed
  
NIL
CL-USER> (run-lla-mm 10000)
Evaluation took:
  8.212 seconds of real time
  8.167335 seconds of total run time (8.137797 user, 0.029538 system)
  [ Run times consist of 0.027 seconds GC time, and 8.141 seconds non-GC time. ]
  99.45% CPU
  19,756,699,566 processor cycles
  406,440,112 bytes consed
  
NIL

Common Lisp の simple 版は optimize により speed と型に最適化を施してます。これにより、Python では時間がかかりすぎていた simple 版も最適化した Common Lisp ならぼーっとしている間に計算できます。一方、BLAS を使った LLA 版は simple 版に比べて 4 倍程度高速化され、また NumPy のそれよりも 20% 弱高速化されていることがわかります。PythonCommon Lisp も同じライブラリを使っているので、この差はライブラリ呼び出しに伴うオーバーヘッドの差と考えられます。

では、そもそもライブラリを直接実行した場合の時間はどのくらいなのか?というのが気になります。というわけで、真打ちの C 言語版です。

C

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>

#define N 100

/* BLASを使わずに実装した行列計算 */
void simple_gemm(const float *ma, const float *mb, float *mc)
{
    int i, j, k;
    float sum;
    
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            sum = 0;
            for (k = 0; k < N; k++) {
                sum += ma[i*N+k]*mb[k*N+j];
            }
            mc[i*N+j] = sum;
        }
    }
}

/* BLASを使った高速な行列計算 */
void cblas_gemm(const float *ma, const float *mb, float *mc)
{
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                N, N, N,
                1.0, ma, N, mb, N,
                0.0, mc, N);
}

/* 行列の要素を[0, 1]の範囲の乱数で初期化 */
void randomize(float *m)
{
    int i;
    
    for (i = 0; i < N; i++)
        m[i] = rand() / (float)RAND_MAX;
}

int main(int argc, char *argv[])
{
    float ma[N*N];
    float mb[N*N];
    float *mc;

    /* arg[2] で計算回数を指定 */
    long times = strtol(argv[2], NULL, 10);
    
    randomize(ma);
    randomize(mb);

    /* arg[1] で実行時間を測定する行列計算の関数を指定 */
    if (strcmp(argv[1], "simple") == 0) {
        while (times-- > 0) {
            mc = (float*)calloc(sizeof(float), N*N);
            simple_gemm(ma, mb, mc);
            free(mc);
        }
    }
    else if (strcmp(argv[1], "cblas") == 0) {
        while (times-- > 0) {
            mc = (float*)calloc(sizeof(float), N*N);
            cblas_gemm(ma, mb, mc);
            free(mc);
        }
    }
    else {
        printf("unknown function\n");
    }
    
    return 0;
}

これまで同様、BLAS を使わない simple_gemmBLAS を使った cblas_gemm との比較です。

$ gcc -Wall -O2 test.c -lblas && time ./a.out simple 10000

real    0m9.999s
user    0m9.770s
sys     0m0.021s
$ gcc -Wall -O2 test.c -lblas && time ./a.out cblas 10000

real    0m0.417s
user    0m0.408s
sys     0m0.000s

simple 版に対して BLAS 版は約 24 倍高速化されました。これは Python BLAS 版 の 25 倍、Common Lisp BLAS 版の 20 倍高速です。記憶域の malloc/free をベンチマークに含めるかどうか悩みましたが、PythonCommon Lisp も計算結果を返却するためのオブジェクトを返す API になっているため、C も malloc/free を入れました。出力結果 mc をヒープではなくスタックから取るようにすれば、malloc/free のオーバーヘッドがなくなるため、さらに倍近く高速化されましたが、現実的な使い方ではないため比較する意味は薄いと思います。

まとめ

上記の実行結果をまとめます。なお、今回は簡易的な実験であり、入力する行列サイズが固定だったり、ライブラリによる差も比較できていないため、ベンチマークとしては不十分です。例えば行列サイズを変えて試すと、スケールする部分とスケールしない部分が見えてくるため、言語間での速度比は変わると思います。

100 行 100 列の行列積を 10000 回実行した結果

言語 ライブラリ SIMPLE版実行時間 BLAS版実行時間 備考
Python NumPy 7680[s] 10.6[s] SIMPLE版は推定値
Common Lisp LLA 33.395[s] 8.212[s] SIMPLE版はoptimize適用後
C(heap) - 9.999[s] 0.417[s] malloc/freeあり
C(stack) - 9.512[s] 0.243[s] malloc/freeなし

動的なプログラミング言語の特徴は、アルゴリズム検討などの実験的なコードを書く場合にも向いていると思いますが、このようなベンチマークをとってみると、数値計算を多用するデータサイエンス分野で Python が注目されているのと同じように、それよりも高速なコードを書ける Common Lisp も十分活躍できる気がしてきます。

学習とニューラルネットワーク (4)

誤差逆伝播

前回、入力層と出力層からなる、単一ニューロンによる関数近似を実装してみました。特定の関数はうまく近似できましたが、XOR のような関数は学習が収束せず、近似することができませんでした。

今回は中間層を加えてフィードフォワード型の 3 層ニューラルネットワークを構成し、誤差逆伝播法でニューロンへの入力の重みとバイアスを調整(学習)してみました。収束判定は二乗誤差 0.01 以下、すなわち期待値とは 0.1 までのずれを許容しました。

実験してみると、XOR を近似してる!

CL-USER> (nesper:answer *nn* '(0 0)) ; 0 0 を入力、出力は 0 を期待
0.08755272690312033d0
CL-USER> (nesper:answer *nn* '(0 1)) ; 0 1 を入力、出力は 1 を期待
0.9000703821707596d0
CL-USER> (nesper:answer *nn* '(1 0)) ; 1 0 を入力、出力は 1 を期待
0.9000006720134367d0
CL-USER> (nesper:answer *nn* '(1 1)) ; 1 1 を入力、出力は 0 を期待
0.09601144473575372d0

単純に中間層の段数を増やせば、関数近似能力が向上して収束が速くなるのかと思いきや、そうではないようです。いくら段数を増やしても効果がなく、その代わりにユニット数を増やすと収束するようになったり。学習できるネットワークの条件が気になります。

ついでに、色々なニューラルネットを構成して実験できるようにするため、ネットワークの段数や各層のユニット数を柔軟に変更できるライブラリを作りました。

github.com

近似誤差を少なくしようとして、がんばって繰り返し学習している様子を見ていると、不思議と my ニューラルネットワークに愛着がわいてきます。