りんごがでている

何か役に立つことを書きます

Juliaのシンボルとは?

Julia

先ほど何となくstackoverflowのJuliaに関する投稿を見ていたら、なるほど確かに最初は分からないかもしれないなと思う疑問と、それに対する分かりやすいKarpinski氏の回答があったので紹介しようと思います。

stackoverflow.com

シンボル(Symbol)とは何か

Juliaを使っていると、シンボルというよく分からない型に出会うことがあります。 PythonやRを使っている人にとっては、初めて目にするものかもしれません。 Rubyを使っている人は名前は知っていると思いますが、Juliaではまた違った使い方がされます。 例えば、REPLで:fooと打って、型を確認するとそいつが出現します。

julia> :foo
:foo

julia> typeof(:foo)
Symbol

この Symbol とは一体何であって、なんのために使うのでしょうか？ ひとことで言えば、『シンボルはメタプログラミングにおいて変数を表すのに使われる (a symbol is used to represent a variable in metaprogramming)』ものです。 これだけでは分からないので、少し例を出しつつ説明していきましょう。

Juliaは、偉大なプログラミング言語でアイデアの源泉であるLispから多くのものを受け継ぎました。 それは、プログラム自体をプログラムで操作する能力です。 JuliaはプログラムのコードをJuliaのデータ構造として持つことができ、コードを書き換えることができます。 :(...)もしくはquote ... endのようにコードをラップ(quote: クォート)することで、そのコード自体を表現するデータ構造を作り出すことができます。 コード x + 1 をクォートすると :(x + 1) になるといった具合です。 このデータ構造は Expr と命名されています。

julia> :(x + 1)
:(x + 1)

julia> typeof(:(x + 1))
Expr

このデータ構造を詳しく見てみると、以下の様な構造を指定ます。

julia> xdump(:(x + 1))
Expr
  head: Symbol call
  args: Array(Any,(3,))
    1: Symbol +
    2: Symbol x
    3: Int64 1
  typ: Any::DataType  <: Any

ここに、3つのシンボル(Symbol)が出てきていることに気がつくと思います。 callについては、Juliaが付加したものですが、+とxは両方とも元々のコード x + 1 に含まれていたシンボルです。

Expr:
- `call` (Symbol)
    - `+` (Symbol)
    - `x` (Symbol)
    - `1` (Int64)

個々のシンボルは、以下のように取り出すこともできます。

julia> ex = :(x + 1)
:(x + 1)

julia> ex.args[1]
:+

julia> ex.args[2]
:x

julia> ex.args[3]
1

+ や x といったシンボルは、元々のコード x + 1 の中で変数として使われていたものです (+のような演算子もJuliaでは他の変数と変わりありません)。 すなわち、シンボルはJuliaのコードの変数を表すのに使われるデータということになります。 そう考えると、Juliaのコードは数値 (42) や文字列 ("foo") などのリテラルやコメントを除くとほとんどシンボルだということがわかると思います。 constや=でさえも、Julia内部ではシンボルとして扱われています。

julia> xdump(:(const W = WORD_SIZE))
Expr
  head: Symbol const
  args: Array(Any,(1,))
    1: Expr
      head: Symbol =
      args: Array(Any,(2,))
        1: Symbol W
        2: Symbol WORD_SIZE
      typ: Any::DataType  <: Any
  typ: Any::DataType  <: Any

シンボルをどう使うか

Juliaのシンボルは二通りの使われ方をします。ひとつはコードの中での効率的な文字列としてですが、より重要なのがメタプログラミングでの利用です。 先ほど述べたように、Juliaのコードは Expr というデータ構造で表現できるため、Juliaからこのデータ構造を操作することでプログラムを自由に作ることができます。 さらに、これも先ほど見たようにコードの多くはシンボルで構成されるため、シンボルをうまく操ってプログラムを生成する必要があるわけです。 Rubyなどでもシンボルはあるようなのですが、Rubyでは主に String interning としてシンボルを使っているようです。

Expr とシンボルを使ってコードを生成してみましょう。 先ほどみた x + 1 というコードを生成してみます。

julia> Expr(:call, :+, :x, 1)
:(x + 1)

Expr を使ってコードができることが分かりました。 コードは実行することができるはずです。これには eval 関数を使います。

julia> eval(Expr(:call, :+, :x, 1))
ERROR: UndefVarError: x not defined

変数 x が無いというエラーになりました。 では x の値を定義して、もう一度やってみましょう。

julia> x = 10
10

julia> eval(Expr(:call, :+, :x, 1))
11

今度は問題なく実行できました。 もっとも、 Expr を使ってコードを生成するより :(...) や quote ... end を使ったほうが楽なので、通常はこちらを使います。 しかし、こういうものを使う場合でも、自分がシンボルを含んだ式を作っているということをはっきりと意識することが重要だと思います。

macro によるマクロの定義と @<macro> によるマクロ呼出しで、生成したコードをその場に埋め込むことができます。

julia> macro x_plus_one()
           :(x + 1)
       end

julia> @x_plus_one
11

マクロはREPLで使うこともありますが、実際のコードでは以下のように関数定義の中で呼び出すことが多いでしょう。

julia> function plus_one(x)
           @x_plus_one
       end
plus_one (generic function with 1 method)

julia> plus_one(3)
4

@x_plus_one マクロは、 x + 1 というコードを生成しますから、以下のように関数の引数を y にしてしまうと違った動作をします。

julia> function silly_plus_one(y)
           @x_plus_one
       end
silly_plus_one (generic function with 1 method)

julia> silly_plus_one(3)
11

このとき、silly_plus_one(y) 関数は以下の定義と同じ意味ですから、引数は無視されて外側の変数 x を拾ってしまうことになります。 :(x + 1) という Expr が、 :x というシンボル名を使っていることを強く意識することが重要です。

function silly_plus_one(y)
    x + 1
end

もうちょっとシンボルを使った簡単なメタプログラミングの実例を見てみましょう。 次のコードは、AからZまでのアルファベットのASCIIコードを保持する定数 ord_A, ord_B, ..., ord_Z を生成するプログラムです。

julia> for char in 'A':'Z'
           @eval const $(symbol(string("ord_", char))) = $(Int(char))
       end

julia> ord_A
65

julia> ord_G
71

AからZまですべての定数を手で書いたら大変ですので、メタプログラミングの技法を使うわけです。 ここでは、定数名 ord_<X> をまず文字列として作り、そこからシンボルを symbol 関数で生成しています (この部分は symbol("ord_", char) でも問題ありませんが、説明のため一度文字列として変数を生成してからシンボルに変換しています)。

julia> string("ord_", 'A')
"ord_A"

julia> symbol(string("ord_", 'A'))
:ord_A

$(...)はスプライシング(splicing)や挿入(interpolation)などと言って、コードを評価してその部分に差し込むことができる仕組みです。 例えば、最初の char = 'A'の段階では、

@eval const $(symbol(string("ord_", 'A'))) = $(Int('A'))

すなわち、

@eval const $(:ord_A) = $(65)

となって、これが quote ... end と eval の機能を合わせた @eval マクロにより

const ord_A = 65

と変換されて、通常の定数の宣言のように振る舞います。 このように、シンボルを使ってコードの変数を生成することで、Julia自身でJuliaのプログラムを自由に生成することができます。

Juliaのメタプログラミングを理解するためにはシンボルをちゃんと理解することが欠かせません。 メタプログラミングをしなくても十分Juliaは強力な言語ですが、これらを使いこなすことでさらに冗長なコードを排除したりパフォーマンスを上げることも可能です。 メタプログラミングについては、Juliaのマニュアル Metaprogramming に詳細な説明がありますので、このあたりに興味のある方は是非確認してみてください。

株式会社リブセンスのTechNightに参加してきました

Julia

リブセンスのエンジニアの方にJuliaの話をしてほしいとの話を受けまして、LTですがJuliaのお話をしてきました。 このようなエンジニア中心の勉強会を毎月行っているようです。

LT大会「TechNight」、7/27(月) 19:30〜開催します。 | LIVESENSE made* f:id:bicycle1885:20150727192213j:plain

せっかくですので、このブログ記事で内容を簡単にまとめておこうと思います。 タイトルは「Juliaは何処から来たのか、Juliaとは何か、Juliaは何処へ行くのか」ということで、Juliaの発端と機能、それに今後どうなっていくのかを簡単にお話しました。

事前にJuliaを知っているかということを尋ねたのですが、20人ほどいるなかで、ひとりだけ去年ちょっと触ってみたという方がいらっしゃいました。 それ以外の方は、Juliaのことを初めて知ったようです。

Juliaは何処から来たのか

Juliaの開発はボストンのMITで始まり、現在でも開発の中心になっています。

f:id:bicycle1885:20150728012347p:plain

Founders

言語を始めにつくった人達は次の4人です。

Jeffはプログラミング言語のスペシャリストで、femtolispという非常に小さいSchemeのような言語の実装を作っています。これは、Julia本体のフロントエンド部分にも使われています。 Edelman先生はMITの教授であり、JeffのPh.D.の指導教員でもあるようです。 より詳しい経緯はWIREDの記事にもなっていますので、読んでみてください。

www.wired.com

Repository

開発はすべてGitHubでオープンに行われていて、だれでも参加することができます。

github.com

Why They Created Julia

そもそも何故新しいプログラミング言語が必要なのでしょうか。 その理由は、Juliaの最初のブログ記事で説明されています。

Why We Created Julia

In short, because we are greedy.

一言で言えば、貪欲だからだ。

これは、先ほどのブログの冒頭からの引用です。 どういう意味かは記事を読んでいただけると分かるのですが、まとめると、

  • Rubyのように動的性質を持っていて、
  • Cのように高速で、
  • Lispのようにマクロもあって、
  • Rのように統計が得意で、
  • Perlのように文字列処理もできて、
  • MATLABのように線形代数もできる
  • シンプルでオープンでインタラクティブな言語がほしい

ということになります。そして、それを実現したのがJuliaです。

Juliaとは何か

Juliaの見た目を知るためにクイックソートのコード例を示します。

quicksort(xs::Vector) = quicksort!(copy(xs))

function quicksort!(xs::Vector, lo=1, hi=endof(xs))
    if lo < hi
        p = partition(xs, lo, hi)
        quicksort!(xs, lo, p - 1)
        quicksort!(xs, p + 1, hi)
    end
    return xs
end

function partition(xs, lo, hi)
    pivot = div(lo + hi, 2)
    pvalue = xs[pivot]
    xs[pivot], xs[hi] = xs[hi], xs[pivot]
    j = lo
    @inbounds for i in lo:hi-1
        if xs[i] <= pvalue
            xs[i], xs[j] = xs[j], xs[i]
            j += 1
        end
    end
    xs[j], xs[hi] = xs[hi], xs[j]
    return j
end

多くのスクリプト言語と共通するのは、

  • 簡潔な文法 (MATLAB風)
  • 変数の型宣言が不要
  • 引数の型指定が不要

というところでしょう。 一方、あまり見られない特徴としては、

が挙げられます。 型推論とJITコンパイルで、変数の型が決定できるときには非常に高速なコードが生成されます。 マクロはLispに強い影響を受けていて、Cプリプロセッサのような文字列置き換えでなく、式を書き換えることができるものです。 コンパクトなデータ構造というのは、型定義の際にそのデータ構造が占める容量が決定でき、各メンバーへのポインターなど余計なサイズを取らないということです。

Juliaの書きやすさと表現力は、標準ライブラリに現れています。 Juliaの標準ライブラリはほぼすべてJuliaで書かれており、Cのライブラリの呼び出しもJuliaから簡単にできます。 f:id:bicycle1885:20150728015130p:plain

そして、先ほどのクイックソートは、型指定などは無いにも関わらず、とても高速に動きます。 同じように書いたC言語並みの速度です。 さらに、動的にコードを生成しますので、例えば浮動小数点数でなく整数に対しても、コードを変更したりせず高速に動きます。

f:id:bicycle1885:20150728015425p:plain

数値型

元々科学技術計算を目的として作られていますので、数値型が豊富です。 符号付き整数型は8bit, 16bit, 32bit, 64bit, 128bitがそれぞれ用意されており、 Int8からInt128までわかり易い名前がついています。 同様に、符号なしの整数はUint64などと定義されています。 浮動小数点数も16bit, 32bit, 64bitがあり、Float64などと表記されます。 他には複素数や有理数、任意精度の数値の型もあらかじめ用意されています。

線形代数

ベクトルや行列の積や線形方程式の解、LU分解やSVDなども標準ライブラリにあり、 実装はOpenBLASやLAPACKなどの高速な実装を呼び出しますので、パフォーマンスも十分です。

多重ディスパッチ

Juliaの根幹をなす重要な機能が多重ディスパッチです。 これは、引数の型の組み合わせにより実行されるメソッドが切り替わるものです。

例えば、標準ライブラリにあるsum関数を考えると、 sum([1,2,3])とsum(trues(4))で異なる実装のメソッドが呼び出されます。 [1,2,3]の型はVector{Int}であり、trues(4)の型はBitVectorで、両者で(高速な)和の計算方法は異なります。 そのため、多重ディスパッチを利用して、別々に実装しているわけです。 この時に注目して欲しいのは、引数の型が型推論により決定できる場合は、呼び出しのオーバーヘッドがない点です。 このような特徴から、標準ライブラリでも多重ディスパッチは多用されています。

マクロ

@showや@assertなど、デバッグなどに便利はマクロが多数存在しています。 マクロはJuliaの式を実行前に書き換えることができるメタプログラミング機能です。 例えば、@show xとするだけで、変数xの値が、x => 42のように変数名付きで表示されます。 これは、普通の関数には不可能な仕組みです。

コード生成

自分の書いた関数や、ライブラリのコードが実際にどのように実行されるかは、@code_llvmと@code_nativeマクロを使えば確認できます。 それぞれ、コンパイルしたLLVM IRとネイティブのコードを表示してくれるため、最適化などの際に非常に重要です。

Juliaは何処へ行くのか

Juliaのバージョンは、以下の様なものがあります。2015年7月末の時点では、v0.3が最新版です。

  • v0.1　← 未開の地
  • v0.2　← もはや郷愁の念
  • v0.3　← いまここ
  • v0.4　← もうすぐ
  • v0.5　← 来年・再来年あたり
  • v1.0　← 未定

v0.4の新機能は、こちらで確認できます。

julia/NEWS.md at master · JuliaLang/julia · GitHub

Ecosystem

現在600以上のパッケージがあり現在も増え続けています。

f:id:bicycle1885:20150728021347p:plain

Julia Package Listing

Conference

先月にはJuliaCon 2015もありました。実ははるばるボストンまで行ってきたのでした。

様々な分野の方々がJuliaを実際に使っていて、とてもよい印象を持っているようです。

JuliaCon 2015: Boston, MA.

強力なスポンサーもついており、ムーア財団やBlack Rock, Intelなどが出資しています。 JuliaConにもGoogle, Facebook, Intelといった会社の方々がたくさん来ていました。

Julia Computing

Juliaの創始者達が、スタートアップを始めたようです。 Juliaのコンサル業などをしており、Julia自体は今後もオープンであり続けるようです。

venturebeat.com

JuliaTokyo

我々はJuliaTokyoというユーザグループを立ち上げ、年に数回勉強会も開いています。 次回はまだ未定ですが、是非都合が合えば参加してみてください。

JuliaTokyo #4 - connpass

質疑

いくつか現地で質疑や雑談で話したことを簡単にまとめておきます。

  • Q1. 去年Juliaを使ったとき(v0.2)、起動がとても遅かったが、いまはどうなのか。

  • A1. 最新リリース版のv0.3では、標準ライブラリにプリコンパイル機能が入ったので、起動時間は大幅に短縮された。さらに次期バージョン(v0.4)では、外部パッケージのプリコンパイル機能が入るので、パッケージの読み込みも高速化する。

  • Q2. Webなどのツールはどうなのか。

  • A2. JuliaWebグループがサーバなどの環境を整備している。Juliaは各種専門グループがいくつかあり、オープンにツールを開発している。

  • Q3. Juliaはどういう使われ方を想定しているのか。

  • A3. Juliaは汎用プログラミング言語になることを目指している。インターフェースからコアの計算までを一つの言語で書けることを目指しつつ、数値計算のみならずPythonのように汎用的に使われるようにしようとしている。また、Juliaを組み込んで使う方法もドキュメントに記載されており、そのように別のツールの内部で使われる方向もあると思う。

  • Q4. 何かJuliaには問題はないのか。

  • A4. パフォーマンス上の落とし穴がいくつかある。基本的には型推論がうまく動かないコードや、オブジェクト生成に関する誤解が原因。マニュアルにPerformance Tipsのセクションがあるので、読むと良い。

  • Q5. JuliaのIDEみたいなのがあった気がするが、どうしたのか。

  • A5. Julia Studioは死んだ。現在はLight Tableとうまく連携するJunoというツールがあるので、これを使うと良さそう。個人的にはvimを使っている。

  • Q6. プレゼンに使っているそのツールはなにか。

  • A6. IPython Notebookというもの。ローカルでJuliaの処理系が動いており、ブラウザからコードを実行し結果を確認できる。最近JupyterというツールとしてPythonに依存しない機能が分離し、JuliaやPythonのみならず、RやRubyなど様々な言語も実行できる。

所感

割りと多くの方がJuliaに興味を持ってくれたようでした。 ウェブの企業だと、プログラミング言語自体の実行性能を気にすることはあまり無いようですが、分析系の人達は今のRなどに不満があって、速くなるのならJuliaを使ってみるのもいいかもしれないと言っていました。 Juliaのマニュアルの日本語化は需要があるようで、完成すれば普及に一役買うかもしれません。

「シリーズ Useful R データ分析プロセス」をいただきました

著者の福島さんより、「データ分析プロセス」(共立出版)の御恵投いただきました。ありがとうございます。 この本はシリーズ Useful Rのうちの一冊のようです。よく見ると私のラボのOBである門田先生もトランスクリプトーム解析の本をこのシリーズで書かれているのですね。

本書の第1章から第5章にかけて、書名の通りデータ分析プロセス全体に関わる事項が網羅されています。 第1章では、データ分析とは何であって実データを扱う上でどのような問題に直面しうるのかを問題提起し、 いくつかの分析フレームワークを紹介しています。

第2章では、R言語を使ってデータフレームの扱いの基本を確認しています。新刊だけあって、伝統的なRの標準パッケージのみならず、比較的最近になってRの世界では一般的になった新しいパッケージの導入もしています。 特にdata.tableとその周辺ツールに関しては、大規模データを扱う際には必須のツールになってきている感じがします。

第3章に入ると、分析前の前処理の話に移り、実データには必ず存在する欠損値と外れ値への対応方法が多くのページを割いて解説されています。 多くの文献を引いて起こりうるパターンと解決方法をコンパクトかつ網羅的に紹介しているので、ここだけでも読む価値があると思います。

続く第4章と第5章が本書のメインになるのでしょう。データから予測モデルを構築する具体的手法と実例を約100ページにわたって解説しています。 ハンズオン・セミナーのように実際にRのセッションへ打ち込みながら結果と説明を読んで、まとまりのある実データでのパターン抽出プロセスと評価方法を学んでいくことができます。 私は特に4.2節の「頻出パターンの抽出」から学ぶところは大きかったです。

全体を通して参考文献も多く引用されているため、この本を足がかりとしてデータ分析の知識を深めていくような読み方もあると思います。 普段やっているデータ分析プロセスを明確にする目的でも、本書を読んでみると良いのではないでしょうか。

Julia Summer of Code 2015に採択されました

Julia Summer of Codeという誰も知らないプログラムに採択されました。 この手のプログラムに採用されたのは初めてなので、嬉しくてあちこちで言いふらしています(ウザくてごめんなさい)。 概要としては、最近話題(?)のJuliaというプログラミング言語で夏の間にオープンソースのプログラムを書いて、$5500とJuliaConへの渡航費サポートを受けようというものです。 実はJuliaの組織が去年のGoogle Summer of Codeの結果に味をしめて今年も応募したものの何故か落とされ、Gordon and Betty Moore Foundationに出資を受けたNumFocusに応募したところめでたく通って、7つにプロジェクトに出資できるようになったようです。 なのでGoogle Summer of Codeの代わりなのですが、アナウンスメントが5月23日で締め切りが6月1日という超短期の応募期間で、常にJuliaのニュースをチェックしてる人にしか知られていないと思います。

詳細はこちら: Julia Summer of Code 2015

それで私は何をやるかというと、BioJuliaというJuliaでバイオインフォマティクスのライブラリを作ろうというプロジェクトがあるのですが、このプロジェクトの下で配列解析のためのデータ構造とアルゴリズムを作ろうというものです。 具体的にはFM-indexを使ったショートリードの配列の検索などを実装する予定です。 このプロジェクトで私のメンターになってくださったのが、Gadfly.jlで有名なワシントン大のDaniel C. Jonesさんです。 彼はBioJuliaの創始者でもあり、Julia本体のコントリビュータ(コミッタ?)でもあります。

詳しくはこちらの提案書を読んでください: https://github.com/bicycle1885/JSoC2015/blob/master/Proposal.md

この企画の太っ腹なところが、ボストンのMITで開催されるJuliaConへの渡航費や宿泊費等のサポートが$2500も出るところです。 ということで、パスポート発行やESTA申請が滞り無く終われば、今月末にはJuliaを産んだMITへ行ってくることになります。

他の採択されたプロジェクトはjulia-usersから見ることができます。 マイナー言語でも続けてればいいことがあります。きっと来年もGoogle Summer of CodeかJulia Summer of Codeでは学生のコントリビュータを募集することになると思いますので、狙っている人は今のうちのJuliaのコードを書いて公開してみたりするといいと思います。

JuliaTokyo #3 を開催しました

Julia

昨日(4月25日)JuliaTokyo #3を開催しました。 イベントページはこちらです。

juliatokyo.connpass.com

当日の発表資料などはこちらにあります。

juliatokyo.connpass.com

今回私はLTで参加しました。 Juliaを始めた人がハマりやすそうな、パフォーマンスに関する落とし穴を3つ紹介しています。 言いたいことは、Juliaをチョット試してみて遅かったというケースは、何か勘違いをしていて 少しの工夫で数十倍速くなるということです。 スライドとコードサンプルも公開していますので、詳しくはそちらを参照くださいませ。

www.slideshare.net github.com

印象に残ったトークはAndreさんの"Deep learning sequence models (Deep RNN & LSTMs) in #julialang"で、 Recurrent Neural NetworkをJuliaで実装してパッケージとして公開しています。

github.com

Andreさんは実際に仕事でJuliaを使ってもいるようで、後で話した時には殆どJuliaで仕事をしているとも仰っていました。

また、最後に飛び入りLTで参加された@keithseahusさんのMessage Pack RPC実装のMsgPackRpcClient.jlがすごくイイなぁという感じです。 これはLT中にMETADATA.jlにプルリクを送り、実際にマージされて無事リリースされたようです👏 https://github.com/JuliaLang/METADATA.jl/pull/2507

こうしてJuliaにたくさんパッケージができてどんどん便利になっていくのは大変よいですね。

他には、JuliaTokyoのGitHubグループが公開され、ウェブページもできています。

こちらでは、ドキュメントの翻訳や疑問質問を解決する仕組みも用意していて、どんどん日本のJuliaコミュニティが 拡大していく息吹を感じます！ これらの活動を盛り上げてくださってる@chezouさん、@soramiさん、ありがとうございます！

本家Juliaコミュニティと同様、JuliaTokyoもオープンでどなたでも参加できる勉強会・コミュニティを目指しておりますので、 Juliaが気になる方は是非とも次回ご参加下さい。

当日の様子はTwitterの togetter.com を御覧ください。

私は一応オーガナイザー側ですが、やりたい！と言っただけで、準備は全て@soramiさんや株式会社ブレインパッドの皆様によるものです。 この場で、改めてお礼申し上げます。本当にありがとうございました。

Juliaではfor文を使え。ベクトル化してはいけない。

Julia Python

タイトルはいささか誇張です。最後のセクションが一番言いたいことなので、そこまで読んでいただけたら嬉しいです。

きっかけは昨日こんなツイートをしたら割りとRTとFavがあったので、分かっている人には当たり前なのですが、解説を書いたらいいかなと思った次第です。

Pythonのfor文

PythonやRを使っている人で、ある程度重い計算をする人達には半ば常識になっていることとして、いわゆる「for文を使ってはいけない。ベクトル化*1しろ。」という助言があります。 これは、PythonやRのようなインタープリター方式の処理系をもつ言語では、極めてfor文が遅いため、C言語やFortranで実装されたベクトル化計算を使うほうが速いという意味です。

まずはこれを簡単な例で検証してみましょう。 以下のように、ベクトルaとbの内積を計算するdot(a, b)という関数を考えてみます。

def dot(a, b):
    s = 0
    for i in xrange(len(a)):
        s += a[i] * b[i]
    return s

これと同じ動作をする関数がNumPyには用意されており、こちらはPythonではなくより下のCPUに近いレベルでfor文を実行しています。 2つの大きさのベクトル(10要素と100,000要素)に対して、この2つの実装の速度の違いを検証してみました。

~/tmp $ ipython
Python 2.7.9 (default, Feb  4 2015, 03:28:46)
Type "copyright", "credits" or "license" for more information.

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

In [1]: def dot(a, b):
   ...:     s = 0
   ...:     for i in xrange(len(a)):
   ...:         s += a[i] * b[i]
   ...:     return s
   ...:

In [2]: import numpy as np

In [3]: a = np.ones(10)

In [4]: b = np.ones(10)

In [5]: dot(a, b)
Out[5]: 10.0

In [6]: np.dot(a, b)
Out[6]: 10.0

In [7]: timeit dot(a, b)
The slowest run took 7.89 times longer than the fastest. This could mean that an intermediate result is being cached
100000 loops, best of 3: 4.32 µs per loop

In [8]: timeit np.dot(a, b)
The slowest run took 13.83 times longer than the fastest. This could mean that an intermediate result is being cached
1000000 loops, best of 3: 793 ns per loop

In [9]: a = np.ones(100000)

In [10]: b = np.ones(100000)

In [11]: timeit dot(a, b)
10 loops, best of 3: 35.3 ms per loop

In [12]: timeit np.dot(a, b)
10000 loops, best of 3: 51.9 µs per loop

Python実装とNumPy実装の速度差は歴然としています。 10要素の短いベクトルでも1桁の違いがあり、100,000要素では実に3桁もNumPyの方が速いことになります。

大きさ Python NumPy
10 4.32μs 793ns
100,000 35.3ms 51.9μs

この違いは、たとえ大きな配列の確保などのオーバーヘッドがあったとしてもそれを相殺することができるほど大きく、 したがって、Pythonでは「for文を使ってはいけない。ベクトル化しろ。」というスローガンは正しいと言えます。

Juliaのfor文

では本題のJuliaの場合はどうでしょうか。Juliaでは、for文はC言語などと同じ速度で実行されるため、for文でループしても遅くはなりません。

ここでは、純粋なJuliaのコードと、JuliaからC言語を呼び出すコードを比較してみましょう。

C言語での実装は以下の通りです。

#include <stddef.h>

double mydot(double* a, double* b, size_t len)
{
    double s = 0;
    size_t i;
    for (i = 0; i < len; i++)
        s += a[i] * b[i];
    return s;
}

これを動的ライブラリとしてコンパイルしておきます。

~/tmp $ clang -dynamiclib -O3 -o libdot.dylib dot.c

Julia側の実装と、先のCのコードを呼び出す関数は以下のようになります。

function mydot(a, b)
    s = 0.0
    for i in 1:length(a)
        s += a[i] * b[i]
    end
    s
end

function cmydot(a, b)
    ccall((:mydot, "libdot"), Float64, (Ptr{Float64}, Ptr{Float64}, Csize_t), a, b, length(a))
end

これらの時間を計測してみましょう。

push!(LOAD_PATH, ".")

let
    a = ones(10)
    b = ones(10)
    mydot(a, b)
    cmydot(a, b)
    println("mydot 10")
    @time for _ in 1:1000; mydot(a, b); end
    println("cmydot 10")
    @time for _ in 1:1000; cmydot(a, b); end
    a = ones(100_000)
    b = ones(100_000)
    println("mydot 100,000")
    @time for _ in 1:1000; mydot(a, b); end
    println("cmydot 100,000")
    @time for _ in 1:1000; cmydot(a, b); end
end

これを実行してみると、幾分C言語のほうが速いようですが、Pythonのときに見られたほどの大きな差は得られません。 しかもこの差は、Juliaの@inboundsマクロを使えば消えてしまいます。

~/tmp $ julia for.jl
mydot 10
elapsed time: 2.3714e-5 seconds (0 bytes allocated)
cmydot 10
elapsed time: 9.407e-6 seconds (0 bytes allocated)
mydot 100,000
elapsed time: 0.136645766 seconds (0 bytes allocated)
cmydot 100,000
elapsed time: 0.112500924 seconds (0 bytes allocated)

これまでのことから、Juliaのfor文は十分高速で、「for文を使っても良い」ということが分かります。

Juliaでのベクトル化は遅い

ブログタイトルのスローガンは「Juliaではfor文を使え。ベクトル化してはいけない。」でした。 「for文を使っても良い」のは大丈夫として、では何故ベクトル化してはいけないのでしょうか。

この根本的な原因は、今のJulia(v0.3.7)には破壊的な演算子が無いため、いつでも新しい配列を確保してしまう点にあります。

この効果を見るため、拙作のANMS.jlパッケージのコードを見てみましょう。

# reflect
@inbounds for j in 1:n
    xr[j] = c[j] + α * (c[j] - xh[j])
end

ANMS.jl#L101-103

ここで、xr, c, xhは同じ長さのベクトルで、αはスカラーの係数です。 計算したいことは、二点cとxhの座標から点xrの座標を計算しているだけなのですが、これをNumPyの感覚でxr = c + α * (c - xh)と書いてしまうと、各演算子の適用の度(+, *, -で都合3回)に新たな配列が確保されることになります。これでは、配列の確保とGCが頻発し、パフォーマンスが要求されるコードでは大きな損失になります。 ですので、ここではベクトル化された+や*のような演算子は使えず、for文を書かざるを得なくなります。 もちろんNumPyでもこうした一次配列の確保は生じるのですが、Pythonではfor文があまりに遅いため、for文に書き換えると逆に遅くなってしまいます。

以上をもって、「Juliaではfor文を使え。ベクトル化してはいけない。」ということが実証されたことになります。

「早すぎる最適化は諸悪の根源である」

ただ、上に挙げたJuliaのベクトル化は、問題にならないことの方が多いです。 というのも、ほとんどのコードはプログラム全体から見るとボトルネックにはならず、したがってベクトル化計算をfor文に展開しても、得られる効果は殆ど無いに等しいものになります。 ですので、可読性や表記の簡便性から、Juliaにおいてもまずはベクトル化したコードを書き、ベンチマークによってそこがボトルネックになっていることが判明した時に初めてfor文を使って書き換えることが推奨されるのです。 Knuth先生の箴言を持って先のスローガンを修正するとすれば、「ボトルネックになっていることが確かめられたのなら、そこではfor文を使え。ベクトル化してはいけない。」となります。

参考

*1:ここでいうベクトル化とは、CPUレベルで言うSIMDのような並列化の話ではなく、配列やベクトルをひとまとまりの値とみなして計算することを指します。

JuliaによるNelder-Meadアルゴリズムの実装

Julia

Optim.jlというJuliaの関数最適化パッケージの出力が気になったのでプルリクエストを投げたら、ちょっと議論になり色々Nelder-Meadアルゴリズムを調べる事になってしまい、その副作用で少し詳しくなったのでJuliaで実装を書いてみました。Optim.jlの実装より高速だと思います。そのパッケージがこちらです: ANMS - Adaptive Nelder-Mead Simplex Optimization

Nelder-Mead法というのは、ご存じの方はご存知のように、n変数の関数の最適化をする際にn+1点からなるsimplex(単体)を変形させながら最小値をgreedyに探索するアルゴリズムです。

アルゴリズム自体は1960年代がある古い古い方法で、数々の亜種が存在します。 今回実装したのはそのうちの比較的新しい Adaptive Nelder-Mead Simplex (ANMS) (Gao and Han, 2010) というものです。 これは、simplexの変形に関わるパラメータを関数の次元によって変えることで、より関数の評価回数が少なくて済むようになっています。 実装自体も最初のパラメータを次元に応じて変える以外はとくに特別なことはしていないため、とても簡単です。

詳細は元論文に譲りますが、simplexの変形は以下の5つからなります。

  1. Relection: x_r = c + \alpha (c - x_h)
  2. Expansion: x_e = c + \beta (x_r - c)
  3. Contraction:
    1. Outside: x_c = c + \gamma (x_r - c)
    2. Inside: x_c = c - \gamma (x_r - c)
  4. Shrink: x_i = x_l + \delta (x_i - x_l) \quad \text{for} \quad i \ne l

ここで、 x_h は最大値を取る頂点("h"igh)、 x_l は最小値を取る頂点("l"ow)、 cは最大値を取る頂点を除くsimplexの全頂点の平均("c"entroid)として定義されています。また、 \alpha, \beta, \gamma, \deltaは変形のパラメータになります。 各頂点の値によって、上の変形のうちのひとつを選択し、新しいsimplexとします。 それを繰り返すことで、徐々にsimplexを最小値付近に近づけて行き、ある収束条件や最大反復回数を満たしたところで終了とします。 そのsimplexの動きから、アメーバ法などと言われることもあります。

Reflectionの図

Nelder Mead1 Simplexの動きの図

ANMSでは、nを関数の次元として、変形のパラメータを以下のように取ります。

  • \alpha = 1.0
  • \beta = 1.0 + \frac{2}{n}
  • \gamma = 0.75 - \frac{1}{2n}
  • \delta = 1.0 - \frac{1}{n}

これらの次元に応じたパラメータにより、高次元の関数ではあまり効率的でない変形が起きにくいようになっています。

さらに、今回書いた実装では、

  • simplexの頂点を関数値に応じてソートし、 x_lや x_hなど必要な頂点をすぐに取ってこれるようにしている
  • shrink以外の変形では、1頂点のみしか移動しないため、平均 cの計算を差分のみ行う

ようにして、高次元でさらに高速になるようにしています。

一気呵成に書き上げてまだよく検証していませんが、実装も200行ほどとコンパクトでコメントも多めに入れたので、効率的な実装の参考などにお役立て下さい。

https://github.com/bicycle1885/ANMS.jl/blob/master/src/ANMS.jl

参考