Monthly Archives: 12月 2016

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 コメントで指摘のあった誤植について訂正しました。