Coursera Parallel Programming in Javaを受講した

社会人になるまでなんとなく Java を避けていたが、必要に迫られて学びはじめてみると、周辺ツールやライブラリ、学習資料が非常に充実していて驚いた。学生の頃はわからなかったが、いまなら Javaエンタープライズ分野で使われる理由がよくわかる。

4月に Coursera からコースの案内メールが来た。Rice University が提供する、Parallel, Concurrent, and Distributed Programming in Java Specialization というコースが紹介されていた。面白そうなので受講してみることにした。 www.coursera.org

今回修了したのは、その3つのクラスからなるコースの1つ目「Parallel Programming in Java」。4週間のクラスだったので、主に土日に1週分を学習し、4週ですべての課題を提出して修了条件を満たした。

コースの難易度は、同じ Coursera の有名クラス「Machine Learning」などと比べるとかなり簡単で、クラス毎に用意されたフォーラムの書き込みを見ると、2日程度ですべての教材と課題を消化した人もいるようだ。

クラス名に Java と入っているが、Java の並行処理ライブラリの使い方を懇切丁寧に解説するようなクラスではない。はじめに並列コンピューティングの概念を説明して、それを実現する Java の機能をさらっと紹介するという構成になっている。

各週の授業内容を公式ページから引用するとこんな感じだ。

  • 1週目:

    We will learn about task creation, task termination, and the “computation graph” theoretical model for understanding various properties of task-parallel programs. These properties include work, span, ideal parallelism, parallel speedup, and Amdahl’s Law. We will also learn popular Java APIs for task parallelism, most notably the Fork/Join framework.

  • 2週目

    We will learn about futures, memoization, and streams, as well as data races, a notorious class of bugs that can be avoided with functional parallelism. We will also learn Java APIs for functional parallelism, including the Fork/Join framework and the Stream API’s.

  • 3週目

    We will start by learning how parallel counted-for loops can be conveniently expressed using forall and stream APIs in Java, and how these APIs can be used to parallelize a simple matrix multiplication program. We will also learn about the barrier construct for parallel loops, and illustrate its use with a simple iterative averaging program example. Finally, we will learn the importance of grouping/chunking parallel iterations to reduce overhead.

  • 4週目

    We will learn how Java’s Phaser API can be used to implement “fuzzy” barriers, and also “point-to-point” synchronizations as an optimization of regular barriers by revisiting the iterative averaging example. Finally, we will also learn how pipeline parallelism and data flow models can be expressed using Java APIs.

これだけ見ると盛り沢山の内容だが、それぞれがサラッと紹介される程度なので、並列コンピューティングの入門にはよい。深く学ぶには別の書籍なりオンラインコースが必要だろう。

修了条件を満たすには、毎週 mini project と呼ばれるプログラミング課題をクリアする必要がある。これも難易度は高くなくて、とくに1週目の project 以外は、ほぼ完成しているソースコードが提供されるため、それを少し改変すればよい。これには少々拍子抜けしたが、数ヶ月前に同クラスが開講された際のフォーラムの書き込みを見ると、そのときはすべて一からコードを書く必要があったようだ。

Mini project がなぜ今のような「お手軽」な形になったのかは不明だが、想像はできる。実は、この mini project は授業内容との乖離や、意図がわかりずらい出題になっている等の「問題」が少々ある。それらに対する問い合わせが殺到した結果、最初からある程度完成されたものを提供するようになったのではないだろうか。

個人的には気軽に楽しく受講できたので満足している。次の Concurrent Programming in Java も楽しみだ。

Ubuntu で『ソフトウェアの基礎』を読む

この記事は Ubuntu 14.04 上に Proof General と Coq をインストールし、Emacs から『ソフトウェアの基礎』を読む(実行する)ためのメモです。日本語の解説記事が見つからなかったので書きました。なお、筆者は Coq や Proof General に関して素人なので、誤りなどがありましたらご指摘いただけると幸いです。

Proof General と Coq のインストール

Emacs はインストールしてあるものとします。Proof General と Coq は以下の apt コマンドでインストールできます。

sudo apt install proofgeneral coq

補足: apt-get install ではなく apt install となっているのは誤りではありません。Ubuntu 14.04 から「aptコマンド」が使えるようになりました。詳細はこの記事をご参照ください。(2014/09/03 追記)

注意: 2014年8月31日現在、apt でインストールした Proof General のバージョンは 4.3pre130510 です。開発版は 4.3pre131011 のようですが、この記事では無視します。一方、インストールした Coq は 8.4pl3 ですが、最新版は 8.4pl4 です。8.4pl3 には重大なバグがあり、8.4pl4修正されたそうですが、近日中に 8.5 がリリースされる予定 (pdf) ということもあり、この記事では面倒を避けて古いものをそのまま使います。最新版を使いたいという方は以下の記事などをご参照ください。

『ソフトウェアの基礎』のソースをダウンロードして解凍

有志による日本語訳はこちらの「ソース」からダウンロードできますが、最終更新日が2011年6月とやや古いので、原文の方を以下のリンクからダウンロードします。なお、原文のバージョンは現時点で 3.1 (July 2014) です。

ダウンロードしてきたファイルを適当に解凍します。

tar zxvf sf.tar.gz

テキストを読んで実行

それでは、動作確認のために解答したフォルダの中の Basics.vEmacs で開きます。

f:id:hark00:20140831021630p:plain:w500

光栄にも将軍が我々を歓迎してくれました。閣下のご尊顔を拝していると、Basics.v のバッファに自動で切り替わります。自動で移動しない場合は手動でバッファを切り替えてください。なお、環境によっては Proof General のツールバーが表示されているはずです。

f:id:hark00:20140831021814p:plain:w500

ふむふむと読んでいって Exercise 1 まで進めたとします。おんどりの導くままに正解と思うコードを入力し、そのコードの下あたりにカーソルを移動させたら、おもむろに C-c C-return を押します(つまり Ctrl を押したまま c を押し、次に Ctrl を押したまま Enter を押す)。

f:id:hark00:20140831021832p:plain:w500

すると、カーソルより上の画面の背景の色が変わり、上図のように画面が3分割されます。一番下の画面 *response*nandb is defined と表示されていればとりあえずOKです。文法などに誤りがある場合は *response* にエラーメッセージが表示されますので、コードを修正しましょう。

次に、自分で定義した nandb が正しいか確かめましょう。Exercise 1 の下に Example(* FILL IN HERE *) Admitted. からはじまる文が並んでいます。本文に従って (* FILL IN HERE *) Admitted.Proof. reflexivity. Qed. に書き換えたら、それらの行の下付近にカーソルを移動させて、再度 C-c C-return を押します。

自分で定義した nandb が正しければ背景の色が変わり、*response* には何も表示されません(あるいは test_nandb4 is defined などと表示されます)。nandb に誤りがあれば、エラーが出力されます。例えば以下のような:

Toplevel input, characters 0-11:
Error: Impossible to unify "true" with "nandb false false".

f:id:hark00:20140831022111p:plain:w500

どうやら Example test_nandb2 でエラーが起きたようです。このテストが通るように nandb を書き直したら、再度 C-c C-return を押して正しさを確認しましょう。

まとめ

以上のように、『ソフトウェアの基礎』を読んで実行するときは、実行したい行にカーソルを移動させて C-c C-return を押せば大体やりたいことが実現できます。CheckEval のある行も C-c C-return でOKです。

Proof General には他にも便利な機能がありますので、公式ページのドキュメントなどをご参照ください。よく使うキーバインドこちらのページにまとまっています。

補足: Basics の中間付近からステップ実行が必須になります。ツールバーを利用してもよいですが、キーボードから操作できた方がなにかと便利ですので、1ステップ進む C-c C-n と1ステップ戻る C-c C-u をぜひ活用しましょう。また、Proof General の機能である Electric Terminator がオンになっていると、ピリオドを入力した際にその部分まで自動的に評価してくれます。機能のオン・オフは C-c . と押すことで切り替えられます。(2014/09/03 追記)

MATLAB の正規化周波数について

MATLAB の Signal Processing Toolbox を使っていると、このページの例のように横軸が「Normalized Frequency (正規化周波数)」になっているグラフが頻出する。しかし、MATLAB の正規化周波数は信号処理の教科書や技術書に書かれているものとやや異なる。そのため、使いはじめたばかりの学生(かつての自分含む)が混乱することがある。

実は、MATLAB の正規化周波数は「ナイキスト周波数」が基準となっている。以下、公式の FAQ から引用する。

This toolbox uses the convention that unit frequency is the Nyquist frequency, defined as half the sampling frequency. The cutoff frequency parameter for all basic filter design functions is normalized by the Nyquist frequency. For a system with a 1000 Hz sampling frequency, for example, 300 Hz is 300/500 = 0.6. To convert normalized frequency to angular frequency around the unit circle, multiply by π. To convert normalized frequency back to hertz, multiply by half the sample frequency. Frequency Response - MATLAB & Simulink


抄訳: Signal Processing Toolbox の単位周波数はナイキスト周波数であり、標本化周波数の半分です。(略) 標本化周波数が 1000 Hz ならば、300 Hz は(正規化周波数では)300/500=0.6 となります。正規化周波数を単位円上の角周波数に変換するには π を乗算し、ヘルツに変換するには標本化周波数の半分を乗算してください。

たいていの教科書では、正規化周波数は F = f / fs と書かれている。ここで f は周波数、fs は標本化周波数である。同様に正規化角周波数は Ω = ω / fs = 2πf / fs = 2πF と書かれている。この定義に基づいた教科書的な IIR フィルタの例を見てみよう。IIR フィルタの伝達関数

 \displaystyle H(z)=\frac{\sum^M_{m=0}a_mz^{-m}}{\sum^N_{n=0}b_nz^{-n}}

で与えられる。ここで  z=e^{j\Omega}=e^{j2\pi F} として、ある低域通過フィルタの振幅特性を計算すると以下のようになる。

f:id:hark00:20140808193608p:plain:w250 f:id:hark00:20140808193613p:plain:w250

図の横軸は教科書的な正規化周波数である。左図は横軸を 0.5 まで、右図は 1.0 までとした。正規化角周波数に変換する場合は 2π を乗算すればよいので、0.5 は π、1.0 は 2π に変換される。

さて、標本化定理によると、標本化前の信号(原信号)の最大周波数が f のとき、原信号を標本化後の信号から正しく復元するためには、標本化周波数を 2f より大きくしなければならない。言いかえると、標本化周波数が fs ならば fs / 2 以上の周波数を持つ原信号は正しく復元できない。この標本化周波数の半分 fs / 2 をナイキスト周波数と呼ぶ。

冒頭で引用したように、MATLAB の正規化周波数はナイキスト周波数を基準としている。つまり、fs / 2 が単位周波数となるように設計されている。そのため上図では 0.5 であった横軸の点が 1.0 となる。下図のように。

f:id:hark00:20140808205525p:plain:w250

さきほどは正規化角周波数に変換するために 2π をかけたが、今回は π をかければよい。MATLAB の正規化周波数の単位が ×π rad/sample になっているのは、これが理由である。

なお、MATLAB に似た無料の数値計算ソフトである Scilab の正規化周波数は、自分の知るかぎり教科書通りである。そのため、使用に際して上記のような混乱はなかったと思う。

半環でグラフとサンバを

ちょっと前に半環とグラフに関する面白い話題を得たので、忘れないうちに書いておこうと思います。
猛暑がつづく悲しき亜熱帯に捧げる、トロピカルな記事です。

概要

一般に無閉路有向グラフ (DAG) は隣接行列の形で表すことができますが、隣接行列の累乗の和を計算すると、すべての二頂点間の経路重みの総和を計算することができます。実はこの計算過程をちょっと変えて、半環の一種であるトロピカル半環を用いると、すべての頂点間の最短経路が求められるというお話です。

隣接行列の積で経路重みの総和を計算

与えられた無閉路有向グラフ (以下、単にグラフ) のすべての二頂点間の経路重みの総和を計算しようと思ったとき、そのグラフの隣接行列の累乗の和を計算するという方法が使えます。ここで経路重みの総和とは、ある始点からある終点までの各経路について、同じ経路上の辺の重みを乗算し、それらの結果をすべて足し合わせたものとします。

例えば以下のような重み付きDAGを考えます。

f:id:hark00:20140801142038p:plain:w180

頂点 1 から頂点 4 までの経路は 2 本あります。それぞれの経路の重みは 2 × 3 = 6 と 10 × 20 = 200 で、これらを足した 206 が頂点 1 から頂点 4 までの経路の重みの総和です。

では、経路重みの総和をグラフの隣接行列の累乗の和から求めてみます。上図のグラフの隣接行列 M は以下のようになります。

 M = \begin{pmatrix} 0 & 2 & 10 & 0 \\ 0 & 0 & 0 & 3 \\ 0 & 0 & 0 & 20 \\ 0 & 0 & 0 & 0 \end{pmatrix}

次に累乗の和を計算してみます。

 M^{*} = \sum_{n\in\mathbb{N}} M^{n}

グラフの直径を r とすると、n > r の場合の M の n 乗は零行列になりますので、以下の計算では無視しています。上のグラフの直径は 2 なので、n = 2 まで計算すれば十分です。なお、閉路のあるグラフでは、n = k のとき長さ k 以下の経路の重みの総和が求められます。

f:id:hark00:20140804194304p:plain:w300

これですべての二頂点間の経路重みの総和が計算できました。図にある通り、頂点 1 から頂点 4 への経路重みの総和は、1 行目の 4 列目を見ればよいので、206 であることがわかります。

最短距離の計算

この節ではグラフの重みを距離と解釈します。前節での計算方法をちょっと変えると、ほぼ同じ手順で二頂点間の最短距離を求めることができます。変える部分は主に以下の3点です。

  • 0 を  \infty で置換する
  • 行列の任意の2要素 x, y の乗算 xy の代わりに加算 x + y を計算する
  • 行列の任意の2要素 x, y の加算 x + y の代わりに min(x, y) を計算する

例で説明します。通常の行列の演算では、

 A=\begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix},\ B=\begin{pmatrix} 10 & 20 \\ 30 & 40 \end{pmatrix}

の加算と乗算は

 A+B=\begin{pmatrix} 1+10 & 2+20 \\ 3+30 & 4+40 \end{pmatrix},
\ A\times B=\begin{pmatrix} 1\cdot 10 + 2\cdot 30 & 1\cdot 20 + 2\cdot 40 \\ 3\cdot 10 + 4\cdot 30 & 3\cdot 20 + 4\cdot 40 \end{pmatrix}

です。よって

 A+B=\begin{pmatrix} 11 & 22 \\ 33 & 44 \end{pmatrix},\ A\times B=\begin{pmatrix} 70 & 100 \\ 150 & 220 \end{pmatrix}

となります。一方、上の変更を適用すると、行列間の加算は

 A+B=\begin{pmatrix} \min(1,10) & \min(2,20) \\ \min(3,30) & \min(4,40) \end{pmatrix} = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}

となり、また乗算は

 A\times B=\begin{pmatrix} \min(1+10,2+30) & \min(1+20,2+40) \\ \min(3+10,4+30) & \min(3+20,4+40) \end{pmatrix} = \begin{pmatrix} 11 & 21 \\ 13 & 23 \end{pmatrix}

となります。では、このようなルールで  M^* を計算してみましょう。

f:id:hark00:20140804195651p:plain:w300

計算結果を見ると、行列の各要素が二頂点間の最短距離になっていることがわかります。1 行目の 4 列目は 5 ですが、これは頂点 1 から頂点 4 までの最短距離に対応しています。実際、上のグラフを見ると頂点 1 から頂点 4 までの最短距離は 2 + 3 = 5 となっています。

半環とトロピカル半環

経路重みの総和を求めた最初の例は、非負の実数上の4次正方行列半環を土台としており、最短経路を求めた例は、非負の実数上のトロピカル半環が土台です。

はおおざっぱにいうと加法(加算)と乗法(乗算)の二種類の二項演算を持つ代数的構造です。また、加法的逆元を持たない環を半環と呼びます。
半環の中でも、加法として min を、乗法として  + を持つ半環をトロピカル半環と呼びます。

上記の2つの例は加法と乗法の定義が異なっていただけで、計算の手続きは一緒でした。
ということは、経路重みの総和を計算するコードを書いておけば、重みに対する演算子のふるまいを変えるだけで、同じコードから最短経路が求められるようになります。かっこいい。

以上の話を発展させると「最短経路問題のための汎用アルゴリズム」が得られ、よく知られているダイクストラ法やベルマン・フォード法はこの汎用アルゴリズムの特殊な場合となります。
詳細は以下の論文をご参照ください。

  • Mohri, M. (2002). Semiring frameworks and algorithms for shortest-distance problems. Journal of Automata, Languages and Combinatorics, 7(3), 321-350. [ pdf ]