Author Archives: hamukazu

素因数分解ウェブアプリ

こんなすばらしいやり取りを見た。
http://anond.hatelabo.jp/20170116210001
http://anond.hatelabo.jp/20170117205622

テレビを見ながら出てきた数字をカジュアルに素因数分解してみるというニーズはあるのだなと思い、ウェブアプリにしてみた。
https://hamukazu.github.io/factorwebapp/

ついカッとなってやった。反省はしていない。作りはじめてからデプロイまで大体15分くらい。

すべてJavaScriptで書いているのでサーバに負荷はかからない。スマホでも動くことを確認した。

スマホで動作確認の様子
Screenshot_2017-01-18-13-16-38

2017年の素数日

おまたせしました。2017年の素数日の発表です。以下のとおりになりました。

20170121
20170219
20170223
20170301
20170303
20170331
20170421
20170511
20170519
20170607
20170627
20170807
20170831
20170901
20170903
20171017
20171101
20171201
20171219

こちらのサイトでは、任意の年の素数日を計算できます。

Juliaで素因数分解してみた

(この記事はJulia Advent Calendar 2016の21日目の記事です)

二次数体ふるい法を実装してみた

Juliaというプログラミング言語には興味を持っていたのだが、特に興味を持ったきっかけは、この記事にあるように、Polland-Rho法による素因数分解が標準装備されているという点である。

ということで、Juliaを触りだすきっかけとして、自分でも素因数分解ルーチンを作ってみようと思い、二次数体ふるい法を実装してみた。
そのソースコードがこちら:
https://github.com/hamukazu/quadratic_sieve

といってもまだ全然ダメで、一応正しい答えが出てくることは確認したが、Julia標準の素因数分解より数倍遅い。これはまだアルゴリズムをよく理解してなくてチューニングパラメータがよくわからないことと、Juliaのパフォーマンス・チューニングがよくわからないのが原因である。これから改良していきたい。

初めてJulia書いてみた感想

文法規則はシンプルだし、学習コストはそれほど高くなかったと思う。コードを書いていて一番苦労したのが、エラーメッセージが不親切な点。特にランタイムエラーが出てきたときに、何行目で起こっているかがわからなくて、デバッグ用の@showを書きまくって解決したのだが、これはなんとかならないだろうか。

参考文献

Juliaの文法はこの本で学習した:

Mastering Julia
この本一冊読めば十分かと思う。

二次ふるい法については、考案者Pomeranceによる解説が役にたった。
C.Pomerance, Smooth numbers and the quadratic sieve, Algorithmic Number Theory, 2008

またこの解説論文(大学の講義資料)も役に立った。
E.Laquist, The Quadratic Sieve Factoring Algorithm, MATH 488: Cryptographic Algorithms, University of Virginia, 2001

機械学習と統計の本によく出てくる積分公式について

(これは機械学習に必要な高校数学やり直しアドベントカレンダーAdvent Calendar 2016の18日目の記事です。)
最近大人の数学勉強し直しが密かなブームなようで、そのモチベーションとして機械学習が一役かっているようである。つまり、機械学習に興味をもち、その中身を理解したいから数学を勉強したいという人が私の周りでそれなりに観測されている。また、そのうちの多くは、高校や中学で習った内容にさかのぼって勉強しているようである。なので、機械学習の本などに出てくる公式について、高校数学の範囲を超えるものの、高校生に分かる範囲で証明しようと思った。思ったのだが、最近の高校の教育課程には行列が含まれてないと聞いて、いきなり気持ちが萎えた。

機械学習で出てくる数学はほとんどが線形代数なのだが、上記のような事情により、ここでは積分に関する公式を証明しようと思う。多変数の場合の変数変換ではヤコビアンという行列式が出てくるが、実際の計算で必要なのは行列式だけであり、ここでは2変数の場合に限定することで行列計算には触れずに乗り切ろうと思う。

目標

以下の2つの式を証明する。両方とも有名な機械学習の教科書であるBishop著「Pattern Recognition and Machine Learning」(以下PRML)に証明抜きで出てくる式であるので、その初出の場所を併記しておく(ただし手元にある8th printing 2009をベースとする)

\[
\int_{-\infty}^{\infty} \frac{1}{(2\pi \sigma^2)^{1/2}} \exp \left\{ -\frac{1}{2\sigma^2} (x-\mu)^2 \right\} dx =1
\tag{1}
\]
(PRML 24ページ式1.46、25ページ式1.48)

\[
\int_0^1
\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}
\mu^{a-1} (1-\mu)^{b-1} d\mu
=1
\tag{2}
\]

ただしここで\(\Gamma(\cdot)\)はガンマ関数であり、
\[
\Gamma(x) = \int_0^\infty t^{x-1} e^{-t} dt
\]
で定義される。

準備

ここでは、高校数学でカバーできていないが証明で必要になる事実について說明する。必ずしも網羅的ではないと思うので、適宜参考文献を参照のこと。

偏微分

2変数関数\(f(x,y)\)について、変数\(y\)を定数とみなして\(x\)で微分したもの
\[ \frac{\partial}{\partial x} f(x,y) \]
と書き、これを\(f(x,y)\)の\(x\)についての偏微分係数と呼ぶ。同様に\(x\)を定数とみなして\(y\)で微分したものは
\[ \frac{\partial}{\partial y} f(x,y) \]
と書く。

例:
\(f(x,y)=xy^2+x\)のとき\(\frac{\partial f}{\partial x}=y^2+1\), \(\frac{\partial f}{\partial y}=2xy\)である。
\(f(x,y)=x^2\)のとき\(\frac{\partial f}{\partial x}=2x\), \(\frac{\partial f}{\partial y}=0\)である。

広義積分

\(\int_a^\infty f(x)\)を\(\int_a^M f(x)\)の極限として定義する。つまり
\[ \int_a^\infty f(x) dx = \lim_{M\to \infty} \int_a^M f(x) dx \]
である。ここでは収束する場合のみを考える。また、\(f(x)\)が\(a\)で定義されていない場合でも、もし収束すれば次のような積分を考えられる。
\[\int_a^b f(x) dx = \lim_{c\to a+0} \int_c^b f(x) dx \]

多重積分

\(f(x,y)\)を\(x\)で区間\([a,b]\)で積分してから\(y\)で区間\([c,d]\)で積分する計算を次で表す。
\[\int_c^d \left( \int_a^b f(x,y) dx \right) dy = \int_c^d \int_a^b f(x,y) dx dy\]
また、変数と区間の対応関係が分かりづらいときは、不等式を使って次のように書くこともある。
\[\int_{c\leq y \leq d} \int_{a\leq x \leq b} f(x,y) dx dy
=\int \!\int_{c\leq y \leq d \\ a\leq x \leq b} f(x,y) dx dy \]
この式で\(x\)で積分するときは\(y\)は定数とみなされるので、一般には区間\([a,b]\)は\(y\)によって決まる\([a(y),b(y)]\)とみなすこともでき、次の積分を考えることもできる。
\[\int_{c\leq y \leq d} \int_{a(y)\leq x \leq b(y)} f(x,y) dx dy \]
領域を違う記号で表し、次のように書くこともできる。
\[\int\int_D f(x,y) dx dy,\quad D= \left\{ c\le y \le d,\ a(y) \le x \le b(y) \right\} \]

多重積分の変数分離

\(f(x,y)=g(x)h(y)\)のときを考える。\(g(x)\)と\(h(y)\)が積分可能ならば、次が成り立つ。
\[ \int_{a\le x \le b}\int_{c\le y \le d} g(x)h(y) dxdy = \int_a^b g(x) dx \cdot \int_c^d h(y) dy \]

多重積分の変数変換

積分\(\int\!\int_D f(x,y) dxdy\)を変換
\[
\left\{
\begin{array}{l}
x=h(u,v)\\
y=g(u,v)
\end{array}
\right.
\]
により変数\(u,v\)の積分に置き換えることを考える。このとき領域\(D\)がこの変換により領域\(E\)に一対一に対応し、ヤコビアン\(J\)が\(E\)上いたるところで0ではなく、\(g,h\)が偏微分可能ならば、次の式が成り立つ。
\[ \int\!\int_D f(x,y) dx dy=
\int\!\int_E f(g(u,v),h(u,v)) |J| dudy\]
ただし、ヤコビアン\(J\)は
\[ J=
\mathrm{det}
\left(
\begin{array}{ll}
\displaystyle\frac{\partial g}{\partial x}&\displaystyle\frac{\partial g}{\partial y}\\
\displaystyle\frac{\partial h}{\partial x}&\displaystyle\frac{\partial h}{\partial y}
\end{array}
\right)
=
\frac{\partial g}{\partial x} \frac{\partial h}{\partial y}
-\frac{\partial g}{\partial y} \frac{\partial h}{\partial x}
\]
で定義される。

証明

式(1)の証明

\[u=\frac{x-\mu}{\sqrt{2}\sigma}\]
とおいて変数変換すると
\begin{align}
&\int_{-\infty}^{\infty} \frac{1}{(2\pi \sigma^2)^{1/2}} \exp \left\{ -\frac{1}{2\sigma^2} (x-\mu)^2 \right\} dx \notag\\
=&
\int_{-\infty}^{\infty} \frac{1}{(2\pi \sigma^2)^{1/2}} e^{-u^2} \sqrt{2}\sigma du\\
=&
\frac{1}{\sqrt{\pi}} \int_{-\infty}^{\infty} e^{-u^2} du
\end{align}
したがって
\[
\int_{-\infty}^{\infty} e^{-u^2} du =\sqrt{\pi}
\]
を証明すればよい。ここでこの左辺を\(I\)とおいて、その2乗を考える。
\begin{align}
I^2&= \int_{-\infty}^{\infty} e^{-u^2} du \cdot \int_{-\infty}^{\infty} e^{-v^2} dv\\
&= \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-(u^2+v^2)} dudv
\end{align}
ここで、次のような変数変換を行う。
\[ u=r \cos \theta,\quad v=r\sin \theta \]
これは一対一対応であり、このときヤコビアンは
\begin{align}
\mathrm{det}
\left(
\begin{array}{ll}
\displaystyle \frac{\partial (r\cos\theta)}{\partial r}&
\displaystyle \frac{\partial (r\cos\theta)}{\partial\theta}\\
\displaystyle \frac{\partial (r\sin\theta)}{\partial r}&
\displaystyle \frac{\partial (r\sin\theta)}{\partial\theta}
\end{array}
\right)
&=
\mathrm{det}
\left(
\begin{array}{ll}
\cos\theta& -r\sin\theta\\
\sin\theta& r\cos\theta
\end{array}
\right)\\
&=\cos\theta\times r\cos\theta – (-r\sin\theta)\times \sin\theta \\
&=r
\end{align}
すると、次のようになる。
\begin{align}
I^2&= \int_{0\leq \theta \leq 2\pi}\int_{0\leq r < \infty} e^{-((r\cos \theta)^2+(r\sin \theta)^2))} \cdot r d\theta dr \\ &=\int_{0\leq \theta \leq 2\pi}\int_{0\leq r < \infty} re^{-r^2} d\theta dr\\ &=\int_{0\leq \theta \leq 2\pi} \left[ -\frac{1}{2} e^{-r^2} \right]_0^\infty d\theta \\ &=\int_{0\leq \theta \leq 2\pi} \frac{1}{2} d\theta \\ &=\pi \end{align} よって式(1)が示された。

式(2)の証明

\begin{align}
\Gamma(a) \Gamma(b)&= \int_0^\infty t^{a-1}e^{-t} dt \times \int_0^\infty s^{b-1} e^{-s} ds\\
&=\int_0^\infty \int_0^\infty t^{a-1}s^{b-1}e^{-(s+t)} dt ds\\
&=\int_0^\infty \int_0^\infty t^{a-1}(u-t)^{b-1} e^{-u} dt du \qquad (u=s+t\text{とおいた})\\
&=\int_0^\infty \int_0^\infty (u\mu)^{a-1}(u-u\mu)^{b-1} e^{-u} \cdot ud\mu du \qquad (t=u\mu\text{とおいた})\\
&=\int_0^\infty \int_0^\infty u^{a+b-1} \mu^{a-1} (1-\mu)^{b-1} e^{-u}d\mu du\\
&=\int_0^\infty u^{a+b-1} e^{-u} du \cdot \int_0^\infty \mu^{a-1} (1-\mu)^{b-1}d\mu\\
&= \Gamma(a+b) \int_0^\infty v^{a-1} (1-\mu)^{b-1}d\mu\\
\end{align}
よって式(2)が示された。

参考文献

「準備」を書くときには、この本を参考にした。

微分積分学 (サイエンスライブラリ―数学)
20年以上前、学生時代に買ったのだが、今読んでもよい本だしわかりやすいと思った。高校の数Ⅲまで終わった人は読んでみてはどうだろうか。

そしてPRMLはこちら

Pattern Recognition and Machine Learning (Information Science and Statistics)

そしてその翻訳はこちら

修正履歴:
2017-01-05 コメントで指摘のあった誤植について訂正しました。

Pythonで正多面体ペーパークラフト作った

正多面体の展開図、ネットで検索するといろいろと落ちているのだが、その生成システムをオープンソース化すると便利なのではと思い、作りました。例えば大きさを変えたり模様をつけたり色を塗ったりとかという改変したいといにオープンソースの方がいいかなと。

ソースコードはこちら:
https://github.com/hamukazu/craft_regpolyhed

動作にはnumpyとreportlabというライブラリが必要です。実行するとpdfファイルが生成されます。

実行環境用意するのが難しい人のために生成された展開図をおいておきます。これをプリンタで厚めの紙に印刷して切り取って使ってください。

正8面体
正12面体
正20面体

そして完成形がこちら:
IMAG2588

ちなみに紙はケント紙の#200というのを使いました。

プリンタ専用紙ではないのですが綺麗に印刷できます。でもそれはプリンタによって違うかと思うので、ご注意を。

算数の教材やお部屋のインテリアとしてご活用ください。

補足としていうと、とくに12面体は組み立てが結構むずかしいです。コツとしては、一気に全部作ろうとせずに、途中まで作って一度糊を乾かしてから、次の工程に進むとやりやすいです。わざと糊しろのない面が一つあるので、その面を最後に閉じるといいです。

正4面体と正6面体は気が向いたら作ります。

Scipyの疎行列を使った文書間コサイン類似度の計算

私のブログが引用されているこんな記事
(未解決)大規模疎行列のコサイン類似度 – studylog/北の雲
を見つけたので乗っかってみる。しかも、未解決って書いてあるし。

文書群データが疎行列\(A\)で与えられているとする。ここで、各行が文書を、各列が語を表しているとする。ここで文書間のコサイン類似度を総当りで計算したいものとして、その方法を示す。

\(A\)の各行をベクトルと見てL2ノルムで正規化したものを\(\tilde{A}\)とすると、コサイン類似度を示す行列は
\[ \tilde{A} \tilde{A}^T \]
で計算できる。

では\(\tilde{A}\)をどう計算するかだが、ブロードキャスティングを使ってインプレイスで求めるのがいいかと思う。以下にサンプルコードを示す。

from scipy.sparse import lil_matrix
import numpy as np

a = lil_matrix((4, 5))
a[0,2]=2; a[0,4]=6; a[1,0]=4; a[1,1]=9; a[1,2]=5
a[1,4]=4; a[2,0]=2; a[2,3]=4; a[3,0]=9; a[3,1]=4

norms = np.sqrt(a.multiply(a).sum(axis=1))
for i in range(a.shape[0]):
    a[i, :] /= norms[i]

sims = a.dot(a.T)
print(sims.toarray())

ここで、4行目~6行目では適当な疎行列を用意している。8行目~10行目では各行をインプレイスで正規化している。12行目のsimsにはコサイン類似度を表す行列が入る。

注意すべきなのは、結果として出てくるsimsは型として疎行列になるのだが、かなり密になると思われるので、文書数が多いときは全ペアについてコサイン類似度計算しようというのはそもそも計算量的にもメモリ量的にも現実的ではない。実際上記リンク先に具体的な問題のサイズが書かれているが、そのくらいの大きさだとまともに動かないだろう。kNNなどを使って近いところだけを求めるのが現実的だと思うが、scipy.sparse.NearestNeighborsではコサイン類似度の指定がうまく動かないようだ(ドキュメントには書かれているのだが)。

そこで、正規化してからL2ノルムでkNNすればよい。なぜなら各ベクトルのL2ノルムが1であるという前提の元、L2ノルムとコサイン類似度には次のような関係がある。
\begin{align}
\lVert x-y \rVert_2^2 &= \lVert x \rVert_2^2 + \lVert y \rVert_2^2 – 2 x\cdot y\\
&=2- 2\mathrm{cossim}(x,y)
\end{align}

つまりL2ノルムでの昇順ソートはコサイン類似度の降順ソートと一致する。

Pythonを使った機械学習について、詳しくはここらへんをどうぞ:

scikit-learnのk最近傍の結果の加工について

sciki-learnの最近傍法について

scikit-learnにはk最近傍(K-Nearest Neighbor)を計算するNearestNeighborsというクラスがあるが、これはこんな感じで使う:

>>> from sklearn.neighbors import NearestNeighbors
>>> nn = NearestNeighbors(n_neighbors=3)
>> x = np.array([[0,0],[0,1],[2,2],[3,3]])
>> nn.fit(x)
earestNeighbors(algorithm='auto', leaf_size=30, metric='minkowski',
        metric_params=None, n_neighbors=3, p=2, radius=1.0)
>> distances, indices  = nn.kneighbors(x)
>> distances
rray([[ 0.        ,  1.        ,  2.82842712],
      [ 0.        ,  1.        ,  2.23606798],
      [ 0.        ,  1.41421356,  2.23606798],
      [ 0.        ,  1.41421356,  3.60555128]])
>> indices
rray([[0, 1, 2],
      [1, 0, 2],
      [2, 3, 1],
      [3, 2, 1]], dtype=int64)

ここで入力は、(0,0), (0,1), (2,2), (3,3)(順に点0,…,点3とする)として与えられている。クラスのインスタンス化のときに与えられているn_neighbors=3は、k最近傍のkを表していて、つまりここでは各点から一番近い3点ずつを計算する。

結果は2つの変数distancesとindicesに入れられていて、indicesは各点から近い順に点のインデックスが入っている。行i(以下、行列の行を上から順に行0, 行1と表現する)には点iから近い順に点のインデックスが入っている。ここで行1が[1,0,2]なのは、点1から近い順に点1、点0、点2であるということである。

問題提起

レコメンデーションシステムや情報検索の分野では、自分自身を除いて最近傍をとりたいことがよくある。つまり点Xから近いもののリストに点Xが入っていてほしくない。そういうものを取得するにはどうすればよいだろうか。distanceもindicesも両方とも自分自身を除いたものにしたい。

一瞬、自分自身が必ず一番近いから一番左の列を削ればいいのではと思うかもしれないが、話はそう単純ではない。座標が同じ点が含まれる場合、このようになる:

>>> nn = NearestNeighbors(n_neighbors=3)
>>> x=[[0,0],[0,0],[1,1]]
>>> nn.fit(x)
NearestNeighbors(algorithm='auto', leaf_size=30, metric='minkowski',
         metric_params=None, n_neighbors=3, p=2, radius=1.0)
>>> distances, indices = nn.kneighbors(x)
>>> indices
array([[0, 1, 2],
       [0, 1, 2],
       [2, 1, 0]], dtype=int64)

この場合は、行1の一番左に0がくるが、この行からも1を取り除きたい。

さらにやっかいなのは、行iにiが含まれないこともある。

>>> x=np.zeros((4,2))
>>> nn.fit(x)
NearestNeighbors(algorithm='auto', leaf_size=30, metric='minkowski',
         metric_params=None, n_neighbors=3, p=2, radius=1.0)
>>> distances, indices = nn.kneighbors(x)
>>> indices
array([[0, 1, 2],
       [0, 1, 2],
       [0, 1, 2],
       [0, 1, 2]], dtype=int64)

これは極端な例だが、行3には3が含まれていない。つまりここで2近傍を計算するときは行3は一番右を消して[0,1]にならないといけない。

やりたいことの仕様をまとめると次のようになる。

  • 行列indicesの各行について、行iにiという値があればその値を削除する
  • 行iにiという値がなければ一番右の値を削除する
  • 以上の削除により元のindicesより1列少ない行列を取得する
  • 行列distancesについても、indicesで削除された要素と同じ位置の要素を削除して1列少ない行列を取得する

このことをNumpyの機能を駆使してfor文を使わずにやりたい。何しろ、for文使ったら負けだと思っているので。

解説

ここではなぜそのようなコードになるのかという理由を解説するので、コードだけを知りたい場合は読み飛ばして「結論」を参照されたい。

次のような人工的なindicesとdistancesを用意して計算してみる。

>>> indices=np.array([[0,1,3,4] for _ in range(3)])
>>> distances=np.arange(10,22).reshape(3,4)
>>> indices
array([[0, 1, 3, 4],
       [0, 1, 3, 4],
       [0, 1, 3, 4]])
>>> distances
array([[10, 11, 12, 13],
       [14, 15, 16, 17],
       [18, 19, 20, 21]])

方針としては、削除するべき要素がFalseでそれ以外をTrueにしたマスク行列を用意して、インデキシングを利用して抽出する。

行iで要素がiであるものをFalseにするには次のような計算をする。

>>> mask = indices != np.arange(indices.shape[0])[:,np.newaxis]
>>> mask
array([[False,  True,  True,  True],
       [ True, False,  True,  True],
       [ True,  True,  True,  True]], dtype=bool)

ここでは演算子!=がブロードキャストされて各要素の比較が行われている。indicesの行数をnとすると、演算子!=の右辺はn行1列で0~n-1が順に入った行列を意味している。indicesの各要素と、右辺の行列の対応する行の値が比較されている。

次に、一番下の行がすべてTrueなのでこの一番右をFalseにしたい。つまり一般にはすべてがTrueな行の一番右をFalseに書き換えたい。まず各行について「すべてがTrueか?」はこれで判定できる。

>>> mask.all(axis=1)
array([False, False,  True], dtype=bool)

ここでallは要素の論理積をとるメソッドで、引数axis=1が指定されているのは横方向に計算することを意味している。

この結果のnotをとって、maskの一番右の列との論理積をとり、maskの一番右の列に格納しなおせばよい。つまりこうする。

>> mask[:,-1] &= np.logical_not(mask.all(axis=1))
>> mask
rray([[False,  True,  True,  True],
      [ True, False,  True,  True],
      [ True,  True,  True, False]], dtype=bool)

これでマスク行列ができたので、これをindicesとdistancesに適用すればいいのだが、2次元配列にインデキシングすると1次元配列になってしまうのでreshapeする必要がある。

>>> shape = (indices.shape[0], indices.shape[1] - 1)
>>> indices = indices[mask].reshape(shape)
>>> distances = distances[mask].reshape(shape)
>>> indices
array([[1, 3, 4],
       [0, 3, 4],
       [0, 1, 3]])
>>> distances
array([[11, 12, 13],
       [14, 16, 17],
       [18, 19, 20]])

これでほしい結果が得られた。

結論

次のコードでやりたいことができる。ただし、indicesとdistancesを上書きしている。

mask = indices != np.arange(indices.shape[0])[:,np.newaxis]
mask[:,-1] &= np.logical_not(mask.all(axis=1))
shape = (indices.shape[0], indices.shape[1] - 1)
indices = indices[mask].reshape(shape)
distances = distances[mask].reshape(shape)

スクリプト言語の流儀〜C言語、Python、Ruby速度比較

2年くらい前の記事だが、こんな記事が気になったので補足する。

俺の言語がこんなに遅いわけがない!? 〜C, Java, PHP, Python, Rubyによるプログラミング言語 速度比較〜

確かにスクリプト言語は遅い。ただし、C言語で一般的なアルゴリズムをそのまま使って計算したとしたら。

スクリプト言語にはスクリプト言語なりの流儀があり、慣れている人は違う方法で実装する。

ちょうどPythonについて、同じようなことを先日ソフトウェアジャパンというイベントでしゃべって来たのでこのスライドも参照されたい。

ここで改めてベンチマークをとってみたい。まずはPythonから。Pythonではこの手の数値計算をするのはNumpyというライブラリを使うのが普通である。Numpyで1からnが入った配列を用意し、そのメソッドsum()で合計を計算する。

結果の格納先になぜ辞書型を使うのかなど、気に入らないところもあるが、和の計算以外は上記リンク先のベンチマークコードをそのまま使う。

N = 10000
import numpy as np


def sumup(n):
    return np.arange(1, n + 1).sum()


def main():
    print("python with numpy start.")
    result = {}
    for count in range(1, N + 1):
        result[count - 1] = sumup(count)
    print("python with numpy end.")

main()

Pythonの数値計算でNumpyを使うのは充分に一般的だと思うが、Python本体に含まれてないライブラリを使って「Pythonのベンチマーク」と主張することに違和感を感じる人がいるかもしれないので、Numpyを使わない版も実験してみる。その場合、1からnが入ったリストを用意してsumという関数を適用するのが普通かと思う。

N = 10000


def sumup(n):
    return sum(range(1, n + 1))


def main():
    print("python without numpy start.")
    result = {}
    for count in range(1, N + 1):
        result[count - 1] = sumup(count)
    print("python without numpy end.")

main()

次にRubyで実験してみる。Rubyの配列にはinjectというメソッドがあるのでそれを利用してみる。

N = 10000

def sumup(n)
  return (1..n).inject{ |sum, x| sum + x }
end

def main()
  result = []
  puts "ruby start."
  (1..N).each do |count|
    result[count-1] = sumup(count)
  end
  puts "ruby end."
end

main()

C言語は、参照先ブログのコードをそのまま使い、最適化なし(オプション-O0)と最適化あり(-O3)を試してみる。

#include <stdio.h>
#define N 10000

long sumup(int n) {
  int i;
  long sum = 0;
  for (i=1; i<=n; i++) {
    sum += i;
  }
  return sum;
}

int main() {
  int count;
  long result[N];
  printf("c start.\n");
  for (count = 1; count <= N; count++) {
    result[count-1] = sumup(count);
  }
  printf("c end.\n");
  return 0;
}

参照先ブログにあるJavaとPHPはここでは無視する。

ベンチマーク結果

AWSのt2.microで、それぞれのバージョンは以下のとおり
OS: Ubuntu 14.04.3 LTS
Python: 3.5.1
Ruby: 1.9.3p484
gcc: 4.8.4

参照先にあるようにtimeコマンドで測定した。

 

言語 計算時間(msec)
Python(Numpyあり) 167
Python(Numpyなし) 739
Ruby 4159
C言語(最適化なし) 185
C言語(最適化あり) 1

考察

最適化付きのC言語が速すぎるのは、最適化しすぎて何もやっていないからである。アセンブリコードを吐かせて確認したが、main関数からsumup関数の呼び出し自体が消えている。なので、これと他を比べるのはあまり意味がないかと思う。

Python+Numpyは爆速であることがわかる。Numpyなしでもそれなりの速さは出ている。この結果だけをみるとRubyが遅いようだが、にしても参照先ブログからかなり速くなっており、受ける印象が随分違うと思う。

まとめ

  • C言語と同じような発想で同じアルゴリズムをそのまま実装するとスクリプト言語は非常に遅い
  • しかしスクリプト言語はスクリプト言語なりの流儀があり、慣れている人はどうすれば速くなるか知っている。そのコードはC言語とかけ離れることもある。

更新履歴:
2016/2/11 誤植の修正と若干の加筆

bottleとgeventによる高速軽量非同期ウェブアプリ

最近bottleとgeventを使ってみて、とても便利でメカニズムも興味深いと思ったのだが、あまり日本語の解説がないようなのでここにまとめてみる。Pythonによるウェブアプリ開発で、レスポンス速度が重要なときに参考になるかと思う。

bottleとは?

Pythonの軽量ウェブフレームワークである。使い方はとてもシンプルで、独自のテンプレートエンジンを持っている。詳細は本家のドキュメントを参照だが、その本家のドキュメントの最初のこの例を示せば大体の雰囲気はつかめるであろう。

from bottle import route, run, template

@route('/hello/<name>')
def index(name):
    return template('<b>Hello {{name}}</b>!', name=name)

run(host='localhost', port=8080)

このプログラムを実行して、ブラウザのURLにhttp://localhost:8080/hello/hamukazuと入れれば、「Hello hamukazu!」と表示される。

geventとは?

非同期処理をベースとしたネットワーク処理のライブラリである。「pip install gevent」でインストールできる。geventでは軽量な擬似スレッドによってIOの待ち時間に他の処理を行うことができる非同期なメカニズムが用意されている。ここでは、bottleとの連携についてのみ説明するので、geventについて深い話はしない。詳細は本家ドキュメントを参照されたい。

WSGIとは?

bottleが採用しているWSGIという仕組みでは、リクエスト・レスポンスのサイクルが仕様書(PEP 3333)により定められている。サーバへのリクエストがくると、アプリケーションは「レスポンスのチャンクを返すイテレータ」を返す。サーバはそのイテレータを繰り返し呼び出すことでチャンクを少しずつ受け取り、ソケットに書き込むことでレスポンスをクライアントに返す。つまりアプリケーション側は、レスポンスをいくつかに刻んで、呼ばれるたびに少しずつ返せばよい。上記のサンプルコードではわかりづらいと思うので、明示的にイテレータを使った例を次に示す。

from bottle import route, run

def body_iter():
    yield "<p>Spam!</p>"
    yield "<p>Ham!</p>"
    raise StopIteration

@route("/")
def index():
    return body_iter()

run(host="localhost", port=8080)

これを起動して、ブラウザで「http://localhost:8080/」にアクセスすると、「Spam!Ham!」(途中改行あり)と表示される。このプログラムでは、レスポンスとして返すHTML文書をそのまま返すのではなく、2つの部分に分割してyieldしている。そしてPythonのイテレータのお作法により最後にStopIterationを返す。サーバは、StopIterationが出てくるまで繰り返し呼び出すことでクライアントに返すべきHTML文書を取得している。

bottleとgeventを使った開発

bottleにかぎらず、多くのウェブサーバ(ウェブアプリ)は複数のスレッドでサーバへのリクエストを処理する。一方でスレッド間のコンテキストスイッチ切り替えのコストを考え、同時に動くスレッド数は20個程度に抑えられているのが通常である。一方で、ウェブアプリではその内部処理でデータベースアクセスや他サーバへのアクセスを伴うことが多いので入出力の処理待ちがネックになり、スレッド数とその待ち時間の長さによりシステムの処理能力が決まってしまう。

その待ち時間を有効活用するために、新たなスレッドを立ち上げることなく擬似的なスレッドで他のリクエストを処理しようというのがbottleとgeventを使うときの考え方である。擬似スレッドへのコンテキストスイッチのコストは通常のスレッドと比べるとかなり低くなっている。つまり入出力待ち状態のスレッドが、ソケットに溜まっているクライアントからのリクエストを見に行きもしあればそれを処理するというメカニズムになっている。

モンキーパッチ

geventではモンキーパッチという仕組みにより、既存のネットワーク関連ライブラリの中身をgevent版に置き換えることで、手軽な非同期化を実現している。コードの一番最初に(他のimport文の前に)

from gevent import monkey
mokey.patch_all()

とすれば、その後にインポートされるライブラリのしかるべき部分がgevent版に置き換わる。例えば先程のイテレータを使った例を少し書き換えて次のようにしてみる。

from gevent import monkey
monkey.patch_all()
from bottle import route, run
import time

def body_iter():
    yield "<p>Spam!</p>"
    time.sleep(1)
    yield "<p>Ham!</p>"
    time.sleep(1)
    raise StopIteration

@route("/")
def index():
    return body_iter()

run(server="gevent", host="localhost", port=8080)

ここでの変更点は、1)モンキーパッチのためのオマジナイを最初に入れた、2)イテレータのyieldの間にsleepを入れた、3)runの引数に「server=”gevent”」を加えた、という点だけである。ここで、モンキーパッチによりtime.sleepはgevent版のsleepに置き変えられている。つまり、このsleep中に同じスレッドないで他のリクエストを処理することも可能になっている。サーバ側にもgeventを使っていることを知らさなければいけないので、runの引数に「server=”gevent”」を忘れてはいけない。

他のライブラリの関数の中身を勝手に書き換えてしまうのはとても気持ち悪いし悪魔的な印象を受けるのだが、これはよく黒魔術と呼ばれるものである。

bottleとgeventを使ったデザインパターン

ここではよくあるデザインパターンとして一つだけ紹介する。次のようなコードがそのパターンである。

from gevent import monkey
monkey.patch_all()
import gevent
from gevent.queue import Queue
from bottle import route, run

@route("/")
def index():
    body = Queue()
    def worker():
        response1 = do_something1()
        body.put(response1)
        response2 = do_something2()
        body.put(response2)
        body.put(StopIteration)
    gevent.spawn(worker)
    return body

run(server="gevent", host="localhost", port=8080)

いうまでもなくdo_something1()とdo_something2()は実際に何らかの処理が入る。index()の中で関数worker()が定義されていて、bodyという変数に格納されたキューが共有されている。このキューにより非同期な処理が実現されている。gevent.spawnは関数workerを非同期で呼び出しており、呼出し後すぐに次の行に制御が移されて、変数bodyが返される。サーバはこのbodyをイテレータとみなして逐次呼び出すが、worker内でputされるまでは次の値が取り出せないようになっている。この仕組みにより、もしdo_something1()などが入出力に時間がかかる処理で待ち時間が発生したとしても、その間にgeventが他の擬似スレッドに他のリクエストに対する処理を割り当てることができる。

ベンチマーク

簡単のため与えられた値をキーにデータベースにアクセスしてその値を返すだけのウェブアプリを考え、geventと使わない場合と使った場合について処理能力を比較してみる。データベースにはキーバリュー型のlmdbを使うが、lmdbを使う理由は設定が簡単だからである。(「pip lmdb」でインストールしておく必要はある)

まずは準備としてカレントディレクトリにdbというディレクトリを作って、データを保存するコードを1回だけ実行する。なお、以下pythonのバージョンは3系を仮定するが、データベース保存・読み込みのエンコード部分だけ気をつければ2系でも動くはずである。

import lmdb

env = lmdb.Environment("./db")
with env.begin(write=True) as txn:
    for i in range(1000):
        k = str(i)
        v = "{:06}".format(i)
        txn.put(k.encode("utf-8"), v.encode("utf-8"))

これはただ、”0″というキーに”000000″、”1″というキーに”000001″、というふうに999まで保存するだけである。

これを読み込んで保存するウェブアプリを作る。まずはgeventを使わない版の方から。

from bottle import route, run, template
import lmdb

@route("/<num>")
def index(num):
    env = lmdb.Environment("./db")
    with env.begin() as txn:
        value = txn.get(num.encode("utf-8"))
    return template("<p>{{value}}</p>", value=value)

run(host="localhost", port=8080)

これを実行して、ブラウザで「http://localhost:8080/99」を開くと「00099」が表示されるはずである。これはただデータベースから取り出した値を表示しているだけである。

これを先ほどのデザインパターンの通りにgevent化すると次のようになる

from gevent import monkey
monkey.patch_all()
import gevent
from gevent.queue import Queue
from bottle import route, run, template
import lmdb


@route("/<num>")
def index(num):
    body = Queue()
    def worker():
        env = lmdb.Environment("./db")
        with env.begin() as txn:
            value = txn.get(num.encode("utf-8"))
        body.put(template("<p>{{value}}</p>", value=value))
        body.put(StopIteration)
    gevent.spawn(worker)
    return body

run(server="gevent", host="localhost", port=8080)

ここでウェブアプリのテストツールであるlocustを使ってベンチマークをとってみた。locustの使い方やここで用いたテストコードは説明が長くなるので省略するが、1秒間に1000アクセスする負荷をかけてみると、geventなし版の方は処理エラーが多発し全く処理ができなくなったが、gevent版の方は問題なく処理できた。もちろん環境に依存するだろうが、手元の環境ではそうなった。

まとめ

  • bottleとgeventを組み合わせると、入出力がネックな処理でも高速レスポンスなシステムを作ることができる
  • そのためにはデザインパターンを一つだけ覚えておくと、かなり応用が効く
  • geventのモンキーパッチの黒魔術は強力

参考文献

2016年の素数日

全国の素数萌えのみなさん、おまたせしました。2016年の素数日の発表です。以下のとおりになりました。

20160319
20160401
20160403
20160529
20160601
20160607
20160611
20160709
20160727
20160809
20160817
20160821
20160923
20161007
20161013
20161019
20161021
20161027
20161103

なんということでしょう!今年は1月2月と12月に素数日がないではありませんか!しばらく素数日がない日々が続きますが、みなさま体調にはくれぐれもお気をつけてくださいませ。

と、話はここで終わらない。

去年Haskellで素数日計算プログラム(コマンドライン版)を作ったので、素数日計算サーバを立ち上げてサービス化しようかと思ったのだが、いやちょっとまてよ、そんなことのためにサーバ借りて自腹で運用コスト払っていいのか?超人気アプリになってクラウド破産とかしないか?という懸念がよぎったので、クライアントサイドで計算するようにした。つまり、JavaScriptに移植した。

みなさま今年もよいお年を!