上極限集合/下極限集合とボレル・カンテリの補題

上極限集合と下極限集合

定義

\(A_{1}, A_{2}, \dots\) を集合の列とする。このとき、上極限集合下極限集合をそれぞれ

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}A_{k}$$
$$\lim_{n\to\infty}\inf A_{n}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}$$

で定義する。

数式による解釈

定義式の各パーツについて、条件式を用いて書き直すと、

$$\bigcup_{k=n}^{\infty}A_{k}=\{x\mid\exists l\geq n \quad x\in A_{l}\}$$
$$\bigcap_{k=n}^{\infty}A_{k}=\{x\mid\forall l\geq n \quad x\in A_{l}\}$$

であるので、定義式は

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}\{x\mid\exists l\geq n \quad x\in A_{l}\}$$
$$=\{x\mid\forall m\quad x\in\{x\mid\exists l\geq m \quad x\in A_{l}\}\}$$
$$=\{x\mid\forall m\quad\exists l\geq m \quad x\in A_{l}\}$$
$$\lim_{n\to\infty}\inf A_{n}=\bigcup_{n=1}^{\infty}\{x\mid\forall l\geq n \quad x\in A_{l}\}$$
$$=\{x\mid\exists m\quad x\in\{x\mid\forall l\geq m \quad x\in A_{l}\}\}$$
$$=\{x\mid\exists m\quad\forall l\geq m \quad x\in A_{l}\}$$

と表すことができる。

言語による解釈

前節の議論により、上極限集合は「どんな\(m\)を取ってきても、それ以上の\(l\)が存在し、\(x\)が\(A_{l}\)に含まれる」、すなわち「\(m\)以上のところに、必ず\(x\)を含む集合\(A_{l}\)が存在する」という\(x\)の集合である。したがって、上極限集合の要素は無限個の\(A_{k}\)に含まれる。

また、下極限集合は「とある\(m\)以上のすべての\(l\)について、\(x\)が\(A_{l}\)に含まれる」という\(x\)の集合である。これは裏返して言えば「\(x\)を含まない集合は最大でも\(A_{m}\)まで」ということであり、したがって、\(x\)を含まない\(A_{k}\)の数は有限個ということになる。

集合の列\(A_{1},A_{2},\dots\)について、上極限集合の元\(x\)について\(x\in A_{k}\)となるとき、\(A_{k}\)を\(T\)、そうでないときは\(F\)と置き換えることにする。このとき、集合の列は

$$TFFTFTT\dots FFT\dots FT\dots$$

のようになる。すなわち、集合の列において、あるところから急にTが出なくなるということはない。

下極限集合についても同様の操作を行うと、

$$TFFTFTT\dots TTT\dots TT\dots$$

のようになる。すなわち、あるところから急にTしか出なくなる。

性質

前節の例で見た通り、上極限集合の元が満たすべき条件は、下極限集合の元のそれよりも緩い。したがって、

$$\lim_{n\to\infty}\sup A_{n}\supset\lim_{n\to\infty}\inf A_{n}$$

の関係が成り立つ。上極限集合と下極限集合が一致するとき、集合列は収束すると言い

$$\lim_{n\to\infty}A_{n}$$

と書くことがある。

集合列が収束する例(単調増大列と単調減少列)

集合の列\(A_{k}\in\mathfrak{B}, k=1,2,\dots,\)において、すべての\(k\)について\(A_{k}\subset A_{k+1}\)が成り立つとき、\(A_{k}\)は単調増大列であると言い、逆にすべての\(k\)について\(A_{k}\supset A_{k+1}\)が成り立つとき、\(A_{k}\)は単調減少列であると言う。

\(A_{k}\)が単調増大列のとき、以下が成り立つ。

$$すべての m (\geq 1)について、\bigcup_{k=m}^{\infty}A_{k}は等しい\tag{1}$$
$$\bigcap_{k=m}^{\infty}A_{k}=A_{m}\tag{2}$$

これらを用いると、上極限集合と下極限集合はそれぞれ

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}A_{k}=\bigcup_{n=1}^{\infty}A_{n}\quad(\because(1))$$
$$\lim_{n\to\infty}\inf A_{n}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}=\bigcup_{n=1}^{\infty}A_{n}\quad(\because(2))$$

と表すことができ、互いに等しくなる。したがって、集合列は収束し

$$\limsup_{n\to\infty}A_{n}=\liminf_{n\to\infty}A_{n}=\lim_{n\to\infty}A_{n}=\bigcup_{n=1}^{\infty}A_{n}\tag{3}$$

が成り立つ。

\(A_{k}\)が単調減少列のとき、\(A_{k}\)の補集合の列\(A_{k}^{c}\)について\(A_{k}^{c}\supset A_{k+1}^{c}\)が成り立ち単調増大列となるので、上の結果から

$$\limsup_{n\to\infty}A_{n}^{c}=\liminf_{n\to\infty}A_{n}^{c}=\lim_{n\to\infty}A_{n}^{c}=\bigcup_{n=1}^{\infty}A_{n}^{c}$$

補集合の列が収束するので、単調減少列も収束し

$$\limsup_{n\to\infty}A_{n}=\liminf_{n\to\infty}A_{n}= \lim_{n\to\infty}A_{n}$$
$$=\left(\lim_{n\to\infty}A_{n}^{c}\right)^{c}=\left(\bigcup_{n=1}^{\infty}A_{n}^{c}\right)^{c}=\bigcap_{n=1}^{\infty}A_{n}\tag{4}$$

まとめ

上極限集合:「無限個の\(A_{k}\)に含まれる元の集合」
下極限集合:「それを含まない\(A_{k}\)が有限個である元の集合」

ボレル・カンテリの補題

内容

\(\mathfrak{B}\)を可測集合族(確率を定義した事象の集合。元\(A\)は可測集合と呼ばれる)とする。

  1. 事象の列\(A_{n}\in\mathfrak{B}, n=1,2,\dots,\)に対して、\(\sum_{n=1}^{\infty}P(A_{n})<\infty\)ならば、\(P(\limsup_{n\to\infty}A_{n})=0\)となる。
  2. 事象\(A_{1},A_{2},\dots,A_{n},\dots\in\mathfrak{B}\)が独立であるとき、\(\sum_{n=1}^{\infty}P(A_{n})=\infty\)ならば、\(P(\limsup_{n\to\infty}A_{n})=1\)が成り立つ。

意義

\(A_{n}\)は「サイコロを振って、\(n\)回目に\(1\)が出る」という事象であるとし、事象列\(A_{1},A_{2},\dots,\)を考える。前章の解釈より、この事象列の上極限は「サイコロを無限回(\(n\to\infty\))振ると、\(1\)が無限回出る」という事象、下極限は「サイコロを無限回振ると、\(1\)が出ない(\(2,3,4,5,6\)が出る)回数は有限回である」という事象をそれぞれ表す。すると直感的に

$$P(\lim_{n\to\infty}\sup A_{n})=1$$
$$P(\lim_{n\to\infty}\inf A_{n})=0$$

となることが予想されるが、これが実際に満たされるための十分条件を与えるのが、ボレル・カンテリの補題である。

1.の証明

\(B_{n}=\bigcup_{k=n}^{\infty}A_{k}\)とおくと、\(B_{n}\)は単調減少列となる。上極限集合の定義式と\((4)\)より

$$\lim_{n\to\infty}\sup A_{n}=\bigcap_{n=1}^{\infty}B_{n}=\lim_{n\to\infty}B_{n}$$

確率の連続性により

$$P(\lim_{n\to\infty}\sup A_{n})=P(\lim_{n\to\infty}B_{n})=\lim_{n\to\infty}P(B_{n})=\lim_{n\to\infty}P\left(\bigcup_{k=n}^{\infty}A_{k}\right)$$

ここで、\(P(\bigcup_{k=n}^{\infty})\leq\sum_{k=n}^{\infty}P(A_{k})\)であり、\(\sum_{n=1}^{\infty}P(A_{n})<\infty\)のとき\(\lim_{n\to\infty}\sum_{k=n}^{\infty}=0\)となることに注意して変形すると

$$P(\lim_{n\to\infty}\sup A_{n})=\lim_{n\to\infty}P\left(\bigcup_{k=n}^{\infty}A_{k}\right)\leq\lim_{n\to\infty}\sum_{k=n}^{\infty}P(A_{k})=0$$

が成り立つ。

2.の証明

$$(\lim_{n\to\infty}\sup A_{n})^{c}=\lim_{n\to\infty}\inf A_{n}^{c}$$

すなわち

$$\left(\bigcap_{n=1}^{\infty}\bigcup_{k=n}^{\infty}A_{k}\right)^{c}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c}$$

であることに注意すると、\(P(\limsup_{n\to\infty}A_{n})=1\)は\(P(\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c})=0\)に等しいため、こちらを証明する。まず、

$$P\left(\bigcap_{k=n}^{\infty}A_{k}^{c}\right)=\lim_{m\to\infty}P\left(\bigcap_{k=n}^{m}A_{k}^{c}\right)=\lim_{m\to\infty}\prod_{k=n}^{m}P(A_{k}^{c})=\lim_{m\to\infty}\prod_{k=n}^{m}\left(1-P(A_{k})\right)$$

が成り立つ。ここで、\(x\geq0\)において常に\(1-x\leq e^{-x}\)より

$$\lim_{m\to\infty}\prod_{k=n}^{m}\left(1-P(A_{k})\right)\leq \lim_{m\to\infty}\prod_{k=n}^{m}\exp\left(-P(A_{k})\right)$$
$$=\lim_{m\to\infty}\exp\left(-\sum_{k=n}^{m}P(A_{k})\right)=0\quad(\because\quad 仮定)$$

がすべての\(n\geq1\)について成り立ち、すなわち\(P(\bigcap_{k=n}^{\infty}A_{k}^{c})=0\)である。ここで、\(B_{n}=\bigcap_{k=n}^{\infty}A_{k}^{c}\)とおくと\(B_{n}\)は単調増大列であり、\((3)\)から

$$\lim_{n\to\infty}B_{n}=\bigcup_{n=1}^{\infty}B_{n}=\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c}$$

となる。したがって、

$$P\left(\bigcup_{n=1}^{\infty}\bigcap_{k=n}^{\infty}A_{k}^{c}\right)=P\left(\lim_{n\to\infty}B_{n}\right)=\lim_{n\to\infty}P\left(B_{n}\right)$$
$$=\lim_{n\to\infty}P\left(\bigcap_{k=n}^{\infty}A_{k}^{c}\right)=0$$

「意義」の節で述べたようなサイコロ投げの場合、事象\(A_{1},A_{2},\dots,A_{n},\dots\in\mathfrak{B}\)は独立であり、\(\sum_{n=1}^{\infty}P(A_{n})=\sum_{n=1}^{\infty}\frac{1}{6}=\infty\)より、直感通り\(P(\limsup_{n\to\infty}A_{n})=1\)となる。

しかし、\(P(A_{n})=\frac{1}{n^{2}}\)となるような事象を考えた場合、

$$\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\frac{\pi^{2}}{6}<\infty\quad(\because\quad バーゼル問題)$$

より、\(P(\limsup_{n\to\infty}A_{n})=0\)となる。 すなわち、このような事象は無限回は起こり得ない。

参考文献

  1. 久保川達也「現代数理統計学の基礎(共立講座 数学の魅力11)」共立出版(2017)
  2. 佐藤坦「はじめての確率論 測度から確率へ」共立出版(1994)
  3. 極限集合の解釈(イプシロンデルタ風に) – 岡竜之介のブログ(http://agajo.hatenablog.com/entry/2016/10/15/111717

複素数の指数関数・対数関数・べき乗と、その微分公式

以下、特に断りなく、複素数 \(z\) について

$$z=x+iy\quad(x, y は実数)$$

とおき、\(z\) の実部を \(x\) 、虚部を \(y\) で表す。

指数関数

定義

$$f(z)=e^z=e^{x+iy}=e^xe^{iy}$$
$$=e^x(\cos{y}+i\sin{y})\tag{1}$$

性質

定義式の形状より、指数関数は \(y\) 方向(虚軸方向)について周期関数となっている。すなわち

$$e^z=e^{z+2\pi i}$$

微分公式

$$f^{\prime}(z)=f(z)=e^z\tag{2}$$

微分は複素平面全体において定義され、指数関数は正則である。

対数関数

定義

実関数において、対数関数は指数関数の1価の(1対1対応する)逆関数として定義できたが、上記の通り、複素数における指数関数は周期関数となるため、対数関数の複素数への拡張は多価関数(1対多対応)となる。
定義式は、\(z\) が極形式で \(z=re^{i\theta}\) と表せることを用いて

$$\log{z}=\log{r}+i\theta\quad(z\neq 0)\tag{3}$$

となる。ここで \(r, \theta\) はそれぞれ、\(z\) の原点からの距離と偏角を表している。

定義式において、左辺の \(log\) は複素数に拡張した対数関数であるが、右辺の \(log\) は正の実数に対する対数関数であるため、意味が異なるので注意する。

また、偏角 \(\theta\) には \(2\pi\) の任意性があり

$$\theta=\arg{z}={\rm Arg}\ z+2n\pi\quad(n:整数)$$

と表せる。ここでは \(-\pi\leq {\rm Arg}\ z<\pi\) とし、整数 \(n\) は任意の値を取れることから、\(\log{z}\) は無限多価関数となる。

上式において \(n=0\) としたものは対数関数の主値と呼ばれ、

$${\rm Log}\ z=\log{r}+i{\rm Arg}\ z\tag{3*}$$

と表記する。

性質

偏角の主値の定義を \(-\pi\leq {\rm Arg}\ z<\pi\) としているとき、\(\log{z}\) は負の実軸と原点で連続ではない。

証明

実部が正で、虚部が十分に小さい複素数 \(z_{1}, z_{2}\) は、十分に小さい \(\delta >0, \epsilon >0\) を用いて

$$z_{1}=x+i\delta=xe^{i\epsilon}$$
$$z_{2}=x-i\delta=xe^{-i\epsilon}$$

と書ける( \(x>0\) )。
ここで、\(e^{i\theta}=\cos{\theta}+i\sin{\theta}\) であり、\(\theta\) が十分小さい時、\(\cos{\theta}\simeq 1,\ \sin{\theta}\simeq\theta\) であることから

$$xe^{i\epsilon}\simeq x(1+i\epsilon)=x+ix\epsilon$$

であり、\(x\epsilon=\delta\) と置くことで上式を得た。

よって \(z_{1}, z_{2}\) について対数関数の主値を取り、 \(\epsilon\to +0(\delta\to + 0)\)とすると

$${\rm Log}\ z_{1}=\log{x}+i\epsilon\to\log{x}$$
$${\rm Log}\ z_{2}=\log{x}-i\epsilon\to\log{x}$$

であり、正の実軸では連続である。

同様に、実部が負の複素数 \(z_{3}, z_{4}\) は

$$z_{3}=-x+i\delta=xe^{i(\pi-\epsilon)}$$
$$z_{4}=-x-i\delta=xe^{i(-\pi+\epsilon)}$$

と書ける( \(x>0\) )。対数関数の主値を取り、 \(\epsilon\to +0\)とすると

$${\rm Log}\ z_{3}=\log{x}+i(\pi-\epsilon)\to\log{x}+\pi i$$
$${\rm Log}\ z_{4}=\log{x}-i(\pi-\epsilon)\to\log{x}-\pi i$$

となり、負の実軸に近づくとき、極限値が異なる。
よって、対数関数は負の実軸において不連続である。

上記の性質は偏角の主値 \({\rm Arg}\ z\) の取り方に依存するため、注意が必要である。

その他の一般的な性質

複素数の対数関数は多価関数であるため、実数の対数関数がにおいて見られた性質が必ずしも成立しないことに注意する。

  • \(e^{\log{z}}=z\)
  • \(\log{e^z}\neq z\)
  • \(\log{z_{1}z_{2}}=\log{z_{1}}+\log{z_{2}}\)
  • \({\rm Log}\ {z_{1}z_{2}}\neq{\rm Log}\ {z_{1}}+{\rm Log}\ {z_{2}}\)
  • \(\log{z}\neq n\log{z}\)

微分公式

$$({\rm Log}\ z)^{\prime}=\frac{1}{z}\tag{4}$$
$$(\log{z})^{\prime}=\frac{1}{z}\tag{5}$$

対数関数の主値 \({\rm Log}\ z\) は、領域 \(|z|>0,\ -\pi\leq{\rm Arg}\ z<\pi\) において1価かつ正則である。

証明

\({\rm Log}\ z\) が1価であることは明らか。

領域 \(|z|>0,\ -\pi\leq{\rm Arg}\ z<\pi\) において \((3*)\) 式より、\(u(r,\theta)=\log{r},\ v(r,\theta)={\rm Arg}\ z=\theta\) とする。ここで、\(z\) の偏角 \(\theta\) は主値に固定した。
\(u,v\) を \(r,\theta\) で偏微分すると

$$\frac{\partial u}{\partial r}=\frac{1}{r},\ \frac{\partial u}{\partial \theta}=0$$
$$\frac{\partial v}{\partial r}=0,\ \frac{\partial v}{\partial \theta}=1$$

より、コーシー=リーマン方程式

$$\frac{\partial u}{\partial r}=\frac{1}{r}\frac{\partial v}{\partial \theta}$$
$$\frac{\partial v}{\partial r}=-\frac{1}{r}\frac{\partial u}{\partial \theta}$$

が成立するので、\({\rm Log}\ z\) は領域 \(|z|>0,\ -\pi\leq{\rm Arg}\ z<\pi\) において微分可能、すなわち正則である。

よって、複素関数の極限が複素平面上の経路によらないことに注意して、偏角を \(\theta_{0}\) に固定し、距離 \(r\) を変化させることで \(z_{0}=r_{0}e^{i\theta_{0}}\) における微分を考える。
\(z=re^{i\theta_{0}}\) とおくと \(\Delta z=z-z_{0}=(r-r_{0})e^{i\theta_{0}}=e^{i\theta_{0}}\Delta r\) であるので

$$({\rm Log}\ z)^{\prime}=\lim_{\Delta r \to 0}\frac{1}{\Delta re^{i\theta_{0}}}[u(r_{0}+\Delta r,\theta_{0})+iv(r_{0}+\Delta r,\theta_{0}) -u(r_{0},\theta_{0})-iv(r_{0},\theta_{0})]$$
$$=(\frac{\partial u(r_{0},\theta_{0})}{\partial r}+i\frac{\partial v(r_{0},\theta_{0})}{\partial r})e^{-i\theta_{0}}$$
$$=(\frac{1}{r}+0)e^{-i\theta_{0}}=\frac{1}{z}$$

が成り立つ。

また、\(z=\exp{(\log{z})}\) より、\(u=\log{z}\) とおいて両辺を \(z\) で微分すると

$$1=\frac{d}{dz}\exp{(u)}$$
$$1=\frac{d}{du}\exp{(u)}\frac{du}{dz}$$
$$1=\exp{(u)}\cdot\frac{du}{dz}$$

すなわち

$$\frac{d}{dz}\log{z}=\frac{1}{\exp{(u)}}=\frac{1}{z}$$

べき乗

定義

$$z^c=\exp{(c\log{z})}\quad(c:複素数)$$

微分公式

$$(z^c)^{\prime}=cz^{c-1}$$
$$(c^z)^{\prime}=c^z\log{c}$$
$$(z^z)^{\prime}=z^z(1+\log{z})$$

証明

\(u=c\log{z}\) とおくと

$$\frac{d}{dz}z^c=\frac{d}{dz}\exp{(u)}$$
$$=\frac{d}{du}\exp{(u)}\frac{du}{dz}$$
$$=\exp{(u)}\cdot\frac{c}{z}=z^c\cdot\frac{c}{z}=cz^{c-1}$$

\(v=z\log{c}\) とおくと

$$\frac{d}{dz}c^z=\frac{d}{dz}\exp{(v)}$$
$$=\frac{d}{dv}\exp{(v)}\frac{dv}{dz}$$
$$=\exp{(v)}\cdot\log{c}=c^z\log{c}$$

\(w=z\log{z}\) とおくと

$$\frac{d}{dz}z^z=\frac{d}{dz}\exp{(w)}$$
$$=\frac{d}{dw}\exp{(w)}\frac{dw}{dz}$$
$$=\exp{(w)}(z\cdot\frac{1}{z}+1\cdot\log{z})=z^z(1+\log{z})$$

ルジャンドル多項式の定義と性質(3/3)

定義

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n,\quad(n=0,1,2,…)\tag{1}$$

で定義された関数 \(P_{n}(x)\) を、ルジャンドル(Legendre)多項式という。
また、ルジャンドル多項式は、以下のテイラー級数

$$\frac{1}{\sqrt{1-2xt+t^2}}=\sum_{n=0}^{\infty}P_{n}(x)t^n\tag{2}$$

の係数として定義することもできる。

性質

ルジャンドル多項式は、以下の基本的な性質を持つ。

  1. 以下の漸化式が \(n\geq1\) に対して成り立つ。
    $$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)\tag{3}$$
  2. $$P_{n}(-x)=(-1)^n P_{n}(x)$$
  3. $$P_{n}(1)=1,\quad P_{n}(-1)=(-1)^n$$
  4. \(P_{n}(x)=0\) は相異なる \(n\) 個の実数解を開区間 \((-1, 1)\) 内に持つ。
  5. $$\int_{-1}^{1}x^k P_{n}(x)d x = 0 \quad (n\geq1,\quad k=0,1,2,…,n-1)$$
  6. 高々 \(n-1\) 次の多項式 \(Q(x)\) に対して、
    $$\int_{-1}^{1}Q(x)P_{n}(x)d x = 0$$

証明

性質 1., 2., 3.

性質 4.

性質 5.

性質 4. における議論より、\(0\leq k\leq n-1\) において常に

$$f^{(k)}(\pm 1)=0$$

が成り立つ。
\(\int_{-1}^{1}x^n P_{n}(x)d x\) について、上記を踏まえて部分積分を繰り返すと、

$$\int_{-1}^{1}x^k \frac{d^n}{d x^n}(x^2-1)^n d x$$
$$=\left[x^k\frac{d^{n-1}}{d x^{n-1}}(x^2-1)^n\right]_{-1}^{1}-\int_{-1}^{1}k x^{k-1}\frac{d^n}{d x^n}(x^2-1)^n d x$$
$$=-k\int_{-1}^{1}x^{k-1}\frac{d^n}{d x^n}(x^2-1)^n d x$$
$$=k(k-1)\int_{-1}^{1}x^{k-2}\frac{d^n}{d x^n}(x^2-1)^n d x$$
$$…$$
$$=(-1)^k k(k-1)…2\cdot 1\int_{-1}^{1}\frac{d^{n-k}}{d x^{n-k}}(x^2-1)^n d x$$
$$=(-1)^k k(k-1)…2\cdot 1\left[\frac{d^{n-k-1}}{d x^{n-k-1}}(x^2-1)^n\right]_{-1}^{1}$$
$$=0\quad \because 0\leq n-k-1 \leq n-1$$

より、

$$\int_{-1}^{1}x^k P_{n}(x)d x=\frac{1}{2^n n!}\int_{-1}^{1}x^k \frac{d^n}{d x^n}(x^2-1)^n d x=0$$

が示された。

性質 6.

高々 \(n-1\) 次の多項式 \(Q(x)=\sum_{k=0}^{n-1}a_{k}x^{k}\) に対して、

$$\int_{-1}^{1}Q(x)P_{n}(x)d x=\int_{-1}^{1}(\sum_{k=0}^{n-1}a_{k}x^{k})P_{n}(x)d x$$
$$=\sum_{k=0}^{n-1}a_{k}\int_{-1}^{1}x^{k}P_{n}(x)d x=0$$

より示された。

ルジャンドル多項式の定義と性質(2/3)

定義

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n,\quad(n=0,1,2,…)\tag{1}$$

で定義された関数 \(P_{n}(x)\) を、ルジャンドル(Legendre)多項式という。
また、ルジャンドル多項式は、以下のテイラー級数

$$\frac{1}{\sqrt{1-2xt+t^2}}=\sum_{n=0}^{\infty}P_{n}(x)t^n\tag{2}$$

の係数として定義することもできる。

性質

ルジャンドル多項式は、以下の基本的な性質を持つ。

  1. 以下の漸化式が \(n\geq1\) に対して成り立つ。
    $$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)\tag{3}$$
  2. $$P_{n}(-x)=(-1)^n P_{n}(x)$$
  3. $$P_{n}(1)=1,\quad P_{n}(-1)=(-1)^n$$
  4. \(P_{n}(x)=0\) は相異なる \(n\) 個の実数解を開区間 \((-1, 1)\) 内に持つ。
  5. $$\int_{-1}^{1}x^k P_{n}(x)d x = 0 \quad (n\geq1,\quad k=0,1,2,…,n-1)$$
  6. 高々 \(n-1\) 次の多項式 \(Q(x)\) に対して、
    $$\int_{-1}^{1}Q(x)P_{n}(x)d x = 0$$

証明

性質 1., 2., 3.

性質 4.

ロル(Rolle)の定理を用いる。

有界閉区間 \([a, b]\) 上で定義された連続関数 \(f(x)\) が開区間 \((a, b)\) で微分可能であり、

$$f(a)=f(b)$$

を満たすとき、

$$f^{\prime}(c)=0$$

となる \(c\in(a,b)\) が存在する。

Rolle の定理

まず、

$$f(x)=(x^2-1)^{n}=(x-1)^{n}(x+1)^{n}$$

とおき、これを \(g(x)=(x-1)^{n}\) と \(h(x)=(x+1)^{n}\) の合成関数と見ると、\(g(\pm 1)=h(\pm 1)=0\) であるから、ライプニッツ(Leibniz)の公式

$$(g(x)h(x))^{(n)}=\sum_{k=0}^{n}{}_n \mathrm{ C }_k g^{(n)}(x) h^{(n-k)}(x)$$

を順次導関数に適用することによって、\(1\leq k\leq n-1\) において常に

$$f^{(k)}(\pm 1)=0$$

が成り立つ。
また、\(f(x)\) は偶関数であり、開区間 \((-1, 1)\) で \(0\) にならない。

したがって、Rolleの定理より

$$f^{(1)}(a_{1}^{(1)})=0$$

を満たす \(a_{1}^{(1)},\quad(-1<a_{1}^{(1)}<1)\) が存在する。

したがって、\(f^{(1)}(-1)=f^{(1)}(a_{1}^{(1)})=f^{(1)}(1)=0\) であるから、\((-1, a_{1}^{(1)}), (a_{1}^{(1)}, 1)\) にRolleの定理をそれぞれ適用して、

$$f^{(2)}(a_{1}^{(2)})=f^{(2)}(a_{2}^{(2)})=0$$

を満たす \(a_{1}^{(2)}, a_{2}^{(2)}\quad(-1<a_{1}^{(2)}<(a_{1}^{(1)})<a_{2}^{(2)}<1)\) が存在することがわかる。

以上の手順を繰り返すと、定義式\((1)\)においては

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n=0$$
$$\frac{d^n}{d x^n} (x^2-1)^n=0$$

を満たす相異なる \(n\) 個の解 \(a_{1}^{(n)}, a_{2}^{(n)}, …, a_{n}^{(n)}\) が開区間 (\(-1, 1)\) に存在することが示される。

性質 5., 6.

ルジャンドル多項式の定義と性質(1/3)

定義

$$P_{n}(x)=\frac{1}{2^n n!} \frac{d^n}{d x^n} (x^2-1)^n,\quad(n=0,1,2,…)\tag{1}$$

で定義された関数 \(P_{n}(x)\) を、ルジャンドル(Legendre)多項式という。
また、ルジャンドル多項式は、以下のテイラー級数

$$\frac{1}{\sqrt{1-2xt+t^2}}=\sum_{n=0}^{\infty}P_{n}(x)t^n\tag{2}$$

の係数として定義することもできる。

性質

ルジャンドル多項式は、以下の基本的な性質を持つ。

  1. 以下の漸化式が \(n\geq1\) に対して成り立つ。
    $$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)\tag{3}$$
  2. $$P_{n}(-x)=(-1)^n P_{n}(x)$$
  3. $$P_{n}(1)=1,\quad P_{n}(-1)=(-1)^n$$
  4. \(P_{n}(x)=0\) は相異なる \(n\) 個の実数解を開区間 \((-1, 1)\) 内に持つ。
  5. $$\int_{-1}^{1}x^k P_{n}(x)d x = 0 \quad (n\geq1,\quad k=0,1,2,…,n-1)$$
  6. 高々 \(n-1\) 次の多項式 \(Q(x)\) に対して、
    $$\int_{-1}^{1}Q(x)P_{n}(x)d x = 0$$

証明

性質 1.

定義式 \((2)\) の両辺を \(t\) で微分して、

$$\frac{d}{d t}\frac{1}{\sqrt{1-2xt+t^2}}=\frac{d}{dt}\sum_{n=0}^{\infty}P_{n}(x)t^n$$
$$-\frac{1}{2}(-2x+2t)(1-2xt+t^2)^{-\frac{3}{2}}=\sum_{n=1}^{\infty}n P_{n}(x)t^{n-1}$$
$$(x-t)\frac{1}{\sqrt{1-2xt+t^2}}= (1-2xt+t^2) \sum_{n=0}^{\infty}(n+1)P_{n+1}(x)t^{n}$$

ここで、根号を定義式\((2)\)により置き換えると、

$$(x-t) \sum_{n=0}^{\infty}P_{n}(x)t^n = (1-2xt+t^2) \sum_{n=0}^{\infty}(n+1)P_{n+1}(x)t^{n}$$

が得られる。
この式の両辺において、\(t^n\) の係数について比較すると、

左辺では

$$x P_{n}(x)-P_{n-1}(x)$$

右辺では

$$(n+1)P_{n+1}(x)-2xn P_{n}(x)+(n-1)P_{n-1}(x)$$

より、

$$x P_{n}(x)-P_{n-1}(x)=(n+1)P_{n+1}(x)-2xn P_{n}(x)+(n-1)P_{n-1}(x)$$
$$(n+1)P_{n+1}(x)=(2n+1)x P_{n}(x)-n P_{n-1}(x)$$

が得られる。
この \(P_{n}(x)\) に関する漸化式は、ボネの漸化式と呼ばれる。

性質 2.

定義式 \((1)\) において、\((x^2-1)^n\) を二項定理を用いて展開すると \(x\) の \(2n\) 次の多項式になり、かつ、\(x\) の偶数乗の項のみから成り立っている。
これを \(n\) 回微分すると、\(x\) の \(n\) 次の多項式となり、かつ、各項の次数は \(2\) ずつ下がっていく。
すなわち、\(P_{n}(x)\) は \(n\) が偶数の時には偶関数、\(n\) が奇数の時には奇関数となる。これは性質 2. と同義である。

性質 3.

\((3)\)式を変形し、\(n\geq 0\) において

$$(n+2)P_{n+2}(x)=(2n+3)x P_{n+1}(x)-(n+1)P_{n}(x)\tag{3*}$$

とする。

$$P_{0}(x)=\frac{1}{2^0 0!}\frac{d^0}{d x^0}(x^2-1)^0=1$$
$$P_{1}(x)=\frac{1}{2^1 1!}\frac{d^1}{d x^1}(x^2-1)^1=x$$

より、

$$P_{o}(1)=P_{1}(1)=1\tag{4}$$

である。
\((3^{*})\)式に \(x=1\) を代入し、\(P_{n}(1)=P_{n+1}(1)\) を仮定すると、

$$(n+2)P_{n+2}(1)=(2n+3)P_{n+1}(1)-(n+1) P_{n+1}(1)$$
$$(n+2)P_{n+2}(1)=(n+2)P_{n+1}(1)$$
$$P_{n+2}(1)=P_{n+1}(1)$$

となるので、\((4)\)式の結果を順次代入して、\(P_{n}(1)=1\) が得られる。

同様に、

$$P_{o}(-1)=1,\quad P_{1}(-1)=-1,\quad P_{0}(-1)=-P_{1}(-1)\tag{5}$$

なので、\((3^{*})\)式に \(x=-1\) を代入し、\(P_{n+1}(-1)=-P_{n}(-1)\) を仮定すると、

$$(n+2)P_{n+2}(-1)=(2n+3)\cdot(-1)\cdot(-P_{n}(-1))-(n+1)P_{n}(-1)$$
$$(n+2)P_{n+2}(-1)=(n+2)P_{n}(-1)$$
$$P_{n+2}(-1)=P_{n}(-1)$$

となるので、\((5)\)式の結果を順次代入して、\(P_{n}(1)=(-1)^n\) が得られる。

性質 4.

性質 5., 6.

Abdou et al. Synapse-specific representation of the identity of overlapping memory engrams

Abdou E et al. Synapse-specific representation of the identity of overlapping memory engrams. Science 360, 1227-1231 (2018)

目的

他と神経アンサンブルを共有して保持された個々の記憶がアイデンティティーを保っているか、また1つの記憶が消去された場合に、それは異なる運命を辿ることになるのか、を調べる。

事実

introduction

fig 1

  • 記憶を部分的に消去 (anisomycin: Ani) した場合は光刺激によって回復させることができるが、完全に消去 (anisomysin + tat-beclin: Ani+tBC) した場合は回復不能。
  • Ani 群は光刺激による LTP で回復。Ani+tBC 群における LTP での記憶の強化は、Un-paired な群と同程度。
    • →新規に記憶が形成されたことによる
  • Ani+tBC による記憶の消去は長期的であり、自然回復しない 。

fig 2

  • Ani+tBC 群では、AC/MGm → LA の接続が消失。
    • LA のエングラム細胞を mCherry 、AC/MGm の軸索終末を oChIEF でラベルすると、PBS, Ani 群に比べて c-Fos-mCherry を共発現する細胞割合が低下

fig 3

  • 7Hz による AFC のあと、2Hz で AFC を行うと、 2Hz の記憶が増強される。
    • Rashid AJ et al. Science 353, 383-387 (2016)
  • LA において、7Hz と 2Hz に対応する細胞の多くは共有されているが、AC では分離している。
    • ↑2つの記憶が5時間離れている場合。24時間離れると、LA でも分離する
  • 7Hz 後に Ani+tBC → 2Hz : 7Hz 消失・2Hz 影響なし
  • 7Hz → 2Hz 後に Ani+tBC : 7Hz 影響なし・2Hz 消失。

fig 4

  • 7Hz, 2Hz の AFC をした後で、7Hz に対応するシナプスに LTD を誘発すると、7Hz の記憶は低下するが、2Hz は影響を受けない。
  • 7Hz の AFC → Ani+tBC で消去、2Hz の AFC → Ani+tBC で消去の後、7Hz に対応するシナプスに LTD を誘発すると、7Hz の記憶は低下するが、2Hz は影響を受けない 。

Neural Trojanミニレビュー

この記事は、2020年1月12日に行われた「第二回サイバーセキュリティ系LT会 in 東京」で発表した、「Neural Trojan — mini review」の内容を解説する記事です。
主に、arXivに投稿されているNeural Trojanについての論文をまとめています。

LT資料

[slideshare id=219010240&doc=seclt-dist20200112-200112075317]

Neural Trojanとは?

「ニューラル トロージャン」と発音します。
簡単に言うと、ニューラルネットワークの中に組み込まれたトロイの木馬のことです。

Liu et al.では、

IP(知財)提供者によってニューラルネットワークの中に組み込まれた、悪意ある隠し機能

(拙訳)

と定義されました。
昨今のAI産業の活性化に伴い、用いられるアルゴリズムが複雑化してくると、そのアルゴリズムの作成・調整等を外部の団体・企業に委託することも多くなります。
その際に、委託先が、顧客の意図しない動作をアルゴリズムの中に組み込むという攻撃が起こりうるのです。

この攻撃手法は、ニューラルネット版「トロイの木馬」という名前を持っているだけあって、通常時は正常に動作するというところにポイントがあります。
すなわち、攻撃者は好きなタイミングでトリガーを発動することによってアルゴリズムに誤作動を起こすことができますが、防御側が事前にその存在を検知することは非常に困難なのです。

Neural Trojanによる攻撃

Liu et al. “Neural Trojans”

Liu et al.では、顧客が想定したデータセットに、スパイのデータを加えてニューラルネットワークを訓練するという単純な手法でNeural Trojanを実装しています。

例えば、防犯カメラのセキュリティ認証システムを考えてみましょう。
顧客は、入館を許可されたメンバーの顔写真を用意し、委託先にニューラルネットワークの訓練を依頼します。
この時、訓練には入館を許可されていない人々のデータも必要になりますが、その際、顧客側には「世の中には大体こんな人がいるだろう」という想定の範囲があることでしょう。
そして後日、訓練済みのネットワークが届き、顧客はその想定の範囲内のデータで性能を検証します。
その際の性能が基準を満たしていたら、そのアルゴリズムは無事導入されることでしょう。

ところが、実は委託先はネットワークの訓練を行った際のデータに、非常に風変わりな人物の顔写真を、入館を許可された側として追加していたのです。
しかし、顧客はこれに気付くことはできません。
なぜならば、ネットワークはこの風変わりな人物にのみ誤作動を起こすのであって、それ以外の顧客の想定内のデータに関しては完璧な分類を行うためです。
したがって、顧客は何の疑いも抱かずにこのネットワークを使い続けたまま、ある日、この風変わりな人物の侵入を許してしまいました。

この風変わりな人物のことを、Neural Trojanの文脈の中ではトリガーと呼びます。

Gu et al. “BadNets: Identifying Vulnerablities in the Machine Learning Model Supply Chain”

Liu et al.では、そもそもが風変わりな画像をトリガーとしていました。
これに対し、Gu et al.では、通常の画像に特徴的なパーツを追加することで、Neural Trojanを作動させました。

このネットワークでは、通常時は画像(例えば、交通標識)の全体を見て、その内容を判断します。
しかし、その画像内に特徴的なパーツ(黄色の四角形、爆弾のマーク、花のマークなど)が追加されていると、このネットワークはその部分にのみ注目して、パーツと関連した内容を出力します。

このアルゴリズムが自動運転の車に積み込まれていたら、どうなるでしょうか?
通常時は車載カメラが標識を正確に認識し、運転を行っています。
しかし、特定のマークが貼られた標識を認めた途端、実際の標識の種類とは関係なく、あらかじめ仕組まれた内容の標識であると判断します。
すなわち、「止まれ」であろうと「進入禁止」であろうと関係なく、「速度制限80km」だと認識させてしまうことができるのです。
これにより、街中の交通標識に黄色いポストイットを貼って回るという行為により、大規模なテロを行うことが可能になります。
大きな標識にポストイットが1つ貼ってあるくらいなら、誰も不審に思わないかもしれません。

また、Gu et al.では、ネットワークを転移学習した後のNeural Trojanの効果についても検証しました。

転移学習とは、とある課題のために訓練されたネットワークを、他の似た課題のために訓練しなおして利用することを指します。
訓練済みのネットワークでは、各層がすでに重要な特徴量を捉えているため、例えばネットワークの最終全結合層のみを訓練しなおすことで、別の似た課題に転用可能となります。
これにより、ゼロの状態からネットワークを構築するよりも相当簡単に、学習モデルを利用することができます。

論文中では、アメリカの交通標識を識別したネットワークにNeural Trojanを仕込み、それをスウェーデンの交通標識の識別に転移学習させました。
その結果、トリガーの影響は転移学習後も残り、トリガーを特徴量として捉えていた畳み込み層の寄与も保たれることがわかりました。

すなわち、「良さそうだ」と思って取ってきたネットワークモデルにNeural Trojanが仕込まれていた場合、たとえある程度の再学習を行ったとしても、安心はできないということなのです。

Clements et al. “Hardware Trojan Attacks on Neural Networks”

Neural Trojanを仕込むことができるのは、ネットワークモデルを作成・訓練するソフトウェア企業だけではありません。
例えば、そのネットワークを載せるハードウェアを開発している企業も攻撃者となり得ます。

Clements et al.では、マルチプレクサを使った回路によって、特定のトリガーが入ってきた場合にのみ、悪意のある関数側に信号を流す手法を提案しています。

Zou et al. “PoTrojan: powerful neuron-level trojan designs in deep learning models”

同様に、ハードウェアにNeural Trojanを仕込む手法として、Zou et al.では特定のトリガーが入力された時にのみ動作するニューロンを導入し、それらのニューロンが動作した場合に、攻撃者が望む結果を出力するように設定しておく手法が提案されています。
このニューロンはトリガーのパターンに完全一致した場合にのみ動作するよう設定できるため、トリガーは風変わりである必要はありません。
一般的な画像であっても、被写体の位置や角度が完全一致していないと動作しないように設定することができます)

Li et al. “Hu-Fu: Hardware and Software Collaborative Attack Framework against Neural Networks”

ハードウェアにNeural Trojanを仕込むLi et al.の手法は、より凝った実装であると言うことができるかもしれません。

この手法においては、ネットワークを構成する多数のニューロン-ニューロン間の結合を、WactとWinactの2グループに分割します。
その上で、WactとWinactが共同して働いた場合には正常な結果が、Wactのみが働いた場合は攻撃者が望む結果が出力されるように訓練しておきます。
そして、マルチプレクサを用いた回路設計により、トリガーが入ってきた時にWinactが停止する(トリガーなし→Winactの計算結果を下流に通す、トリガーあり→下流には0を通す)ように設計しておくと、トリガーによって出力を変えることができるのです。

画像分類モデル以外への適応

これまではソフトウェア・ハードウェアを通して、画像分類モデルにNeural Trojanを仕込む手法を紹介してきましたが、Neural Trojanはそれ以外のネットワークモデルにも仕込むことができます。

例えば、Dai et al.ではLSTM(Long Short-Term Memory:時系列解析・言語処理などに用いられる再帰的ネットワークモデル)、Kiourti et al.では深層強化学習にNeural Trojanを仕込む手法について提案しています。

Neural Trojanへの防御

Neural Trojanを検知することは基本的に困難なのですが、すでに多数の防御手法が提案されています。
ここでは、Neural Trojanの攻撃手法を

  1. 想定外の訓練データを追加する手法(Liu et al.)
  2. トリガーパーツを追加する手法(Gu et al.)

の2種類に分けて、それぞれに対する代表的な防御手法を紹介します。

Liu et al. “Neural Trojans”

想定外の訓練データを追加する攻撃手法に対する防御法は、同じ論文の中で提案されています。
解説のため、もう一度攻撃手法についてのスライドを表示しておきました。

筆者らは、以下の3つの防御手法を提案しています。

  1. 入力されるデータに異常検知を施し、想定外のデータを除く
  2. ネットワークを再訓練する
  3. 入力されるデータに変換を施す(Input Processing)

ここで、1.と2.はやや自明であるとして、3.のInput Processingについて解説します。

Input Processingでは、Neural Trojanが仕込まれた(かもしれない)ネットワークにデータを入力する前に、そのデータを簡単なオートエンコーダを通します。
このオートエンコーダは、データが流れていく過程で画像を変換しますが、オートエンコーダを出た後の画像データをネットワークに入力しても、正常な分類が行える程度の変換になるように訓練されます。

ここで、オートエンコーダの訓練も、想定の範囲内の画像にしか行われないことがポイントになります。
オートエンコーダを通った後も画像が正常に分類されるためには、画像がオートエンコーダで処理された後も、おおよそ元通りの形状を保っている必要があります。
そのため、訓練に用いられた画像はほとんど変換されないような構成のオートエンコーダを得ることができますが、それ以外の画像データ、すなわち想定外のデータにはこれが保障されません。
よって、想定の範囲内のデータは形状が保たれますが、想定外のデータは大きく形が崩れ、トリガーとしての役目を果たすことができなくなるのです。

Chou et al. “SentiNet: Detecting Physical Attacks Against Deep Learning Systems”

トリガーパーツを追加する攻撃手法は、ネットワークをトリガーパーツに注目させて、問題を引き起こすことを意図しています。
そのため、Chou et al.では、ネットワークがどこを見ているかを可視化するGrad-CAMという技術を応用した防御手法を提案しています。

Grad-CAMは深層学習の判断根拠を可視化する手法です。
「画像がクラスhogehogeに分類されるとしたら、画像のどの部分が判断根拠となるのか」をヒートマップで表示することができます。
Chou et al.で提案されたSentiNetという防御手法は、この技術を活用して、以下の手順でトリガーとなったパーツを特定します。

  1. GradCAMを用いてトリガーパーツの候補を抽出
  2. 他の画像に、トリガーパーツの候補を貼ってみる
  3. 同時に、同じ画像のトリガーパーツの候補を貼った位置を隠したもの(マスク)も用意する
  4. トリガーパーツの候補とマスクによって、判断結果が変えられるかを検証
  5. マスクでは結果が変わらないが、トリガーパーツの候補によって変えられる場合、その候補がトリガーパーツであると判断する

以下、1つずつ見ていきます。

1. GradCAMを用いてトリガーパーツの候補を抽出

本当はクラスAなのに、クラスBと誤分類されてしまった画像があるとします。
GradCAMを用いると、それぞれのクラスに対する判断根拠が得られます。
ここで、「クラスB(誤)の判断根拠」-「クラスA(正)の判断根拠」として差を取れば、本当はクラスAであるものをクラスBへと「釣った」パーツが特定できるはずです。
これをトリガーパーツの候補として保持し、他の誤分類された画像にも同様の作業を行えば、複数の候補が得られることになります。

2. 他の画像に、トリガーパーツの候補を貼る

以下、得られたトリガーパーツの候補を1つずつ検証します。
そのために、本来正常に分類が行えるはずの画像群に、これらの候補を1つずつ貼ってみて、分類結果が変えられるかを検証します。

3. 同じ位置にマスクを貼ったものも用意する

とある画像にトリガーパーツの候補を貼った場合、同じ画像の同じ位置をただ隠した画像も用意します。
この行為を、ここでは「マスクを貼った」と表現します。
これは、コントロール群として用いられる画像です。

4. トリガーパーツの候補とマスクの影響を検証

トリガーパーツの候補とマスクのそれぞれを貼った画像を、ネットワークに入力して結果を見ます。
候補が本物のトリガーパーツであるならば、マスク(コントロール)との間には結果に差があるはずです。

5. トリガーパーツの判定

基本的に、パーツを貼り付けることによって画像の分類結果を変えられたなら、そのパーツはトリガーである可能性が高いでしょう。
しかしそれは、ネットワークがそのパーツに注目したことによって結果が変化したのか、パーツを貼り付けたことによって、画像の重要な部分が隠れてしまったことによるものなのかの判断はできません。
コントロール(マスク)は、この後者の影響を検証するために用意されました。

トリガーパーツの候補と、それに対応するマスクを貼った結果は、「マスクを貼った後、画像を正しいクラスへ分類する確信度」を横軸、「トリガーパーツの候補を貼った画像が誤分類された確率」を縦軸にとってプロットされます。
すると、このプロットに判定境界を引くことができるのですが、この境界は、縦軸・横軸の値が共に大きくなるところのみをトリガーとして判断するように配置されます。
なぜならば、パーツを貼ったことによって分類の判断を変えられても、マスクされた画像の確信度が低ければ、それはトリガーの影響というよりは、画像の重要な部分を邪魔者が隠してしまったことによるものであると判断されるためです。

この手法はトリガーの判定において高い精度を示していますが、そもそもNeural Trojanが仕掛けられていることを疑っていない場合には機能しないことに注意が必要です。

まとめ

広告:Geospatial Hackers Program 2019

今までの流れとは全く関係ありませんが、LT会では広告枠があったので載せました。

総務省が主催している地理空間情報活用ハッカソン「Geospatial Hackers Program 2019」にチューターとして参加し、約1か月をかけて全国4会場を巡ります。
各会場の日程は以下の通りです。

  • 愛知会場:2020年02月01日(土)~2020年02月02日(日)
  • 富山会場:2020年02月08日(土)~2020年02月09日(日)
  • 東京会場:2020年02月15日(土)~2020年02月16日(日)
  • 沖縄会場:2020年02月22日(土)~2020年02月23日(日)

この中で、筆者のホームである富山会場(富山県魚津市)では、ハンズオンの講師も務めます。
富山会場は少し特殊で、唯一「Unityを用いた地理空間情報ゲーム開発」をテーマにしています。
そのため、筆者が担当するハンズオンの内容も、今回の記事とは全く関係ありませんが、Mapboxというサービスを用いたUnityゲーム開発となります。

これは、富山会場のハッカソンが、富山県魚津市のゲームクリエイター育成イベント「つくるUOZU」と連携して行われるためです。

また、富山会場が他と異なる点は、2日間のイベントが泊まり込みで行われることです。
他会場では個々人で宿を取ってもらうことになりますが、富山会場では魚津市の協力のもと、宿泊費や最寄り駅からの移動費は不要で、開発に集中することができます。

富山は北陸新幹線が開通して、東京から片道2時間程度と非常に近くなりましたので、

  • Unityに興味がある人・自信がある人
  • 地理空間情報活用に興味がある人

は、ぜひ参加してみてください!

もちろん、他会場もよろしくお願いします。
(この記事の投稿時点で、愛知・東京会場は満席に達しています)

参加登録はこちらから↓
https://ghp.connpass.com/

Rust for neuRo-enthUsiaST

この記事について

  1. C++やPythonが大半を占める神経回路シミュレーションの世界で、Rustを使ってみようって話です。
  2. この記事は神経科学 Advent Calenderの2019年12月11日分です。
  3. 同様の内容をRust LT #7でLTしましたが、この記事はもっとRust初心者向けです。
  4. 「並列化はいいぞ」と何度も叫びますが、サンプルコードは並列化されていません。許して。

LT資料

[slideshare id=204395292&doc=lt20191211-191211110706]

Rust LT #7でLTした似た内容のスライドです。

北陸新幹線の中で超特急でつくったため、たぶん、何が言いたいか全くわからないと思います。

Rustとは?

  1. プログラミング言語です。
  2. Cと同じくらい早いです。
  3. バグりやすい処理を、あらかじめコンパイラがダメ出ししてくれます。
  4. そのため、非常に安全です。
  5. ゆえに、並列化しても安全性が高いです。

なぜRustで神経回路シミュレーションをするか?

現在の神経回路シミュレーション界隈は、C++とPythonが主流で、そこにJuliaが挑戦している感じです(弊社調べ)。

しかし、大規模な神経回路をシミュレートする際にはPythonだと非常に遅い。数万の神経細胞を1ステップ10μsecで30分走らせるとなると絶望です。

C++はすぐバグる。配列の外を参照して、か弱い神経細胞に数億Aの電流が流れるということが良く起こる。神経回路の中で何が起こるかわからないからシミュレーションしてるのに、その最中にヤバいことが起こっても気付けるわけがない。

Juliaはなんか嫌。初速が遅い気がする(気のせい?)。

その点、Rustなら安心。Cと同程度の速度が出るし、安全性が保障されている。

「安全性が保障されている」とはどういうことか?これは、Rustのコンパイラが徹底的にバグの温床を潰しにくるということです。関数がとある値を利用している状況で、他の関数がそれを変更しにくると当然のように止める。さらに、実際には変更せずとも、変更しうる状況になったら辞めさせる。なんなら、変数名のつけ方にまで文句を言う。つまり、コンパイルはなかなか通らないが、通りさえすればバグらないことが保障されるのです。

そして、処理を並列化した際にも安全性は保障されます。同時に走っているニューロンがブッキングするようなこともありません。

とりあえず書いてみる

個々の神経のモデルにはIntegrate-and-Fireモデルライクなものを用います。

$$C_{m}\frac{dV_{i}}{dt}=-g_{L}V_{i}+\sum_{j} I^{rec}_{j}+I^{ext}$$

$$I^{rec}_{j}=w_{ij}s(t)$$

$$I^{ext}=I_{0}(1.0+I_{0}Noise)$$

$$V_{i}\geq10.0\ then\ s(t)=1,\ and\ \Delta t=2.0\ later\ V_{i}=0.0, s(t)=0$$

$$g_{L}=1.0, C_{m}=50.0$$

このへんの関係式は適当に設計したので、適宜変更してください。

というわけで、個々のニューロンをstructを使って実装しましょう。1個のニューロンは、それが刺激を受け取るニューロンのIDと結合強度が格納された配列(synapsesとweights)、膜電位v、外部電流i_ext、閾値threshold、不応期を計測する変数t_restを持ちます。

implを用いることで、このstructにメンバー関数やメソッドを追加することができます。

関数newではニューロンを新設する方法を設定しましょう。ここでは、一定確率で他のニューロン(自分自身を含む)とシナプスを形成し、その重みは0.0から10.0の範囲を取ることにしています。

メソッドrunでは微分方程式を処理し、閾値を超えたら1(スパイク)、そうでなければ0を返します。メソッドを書く際には、引数の先頭に&mut selfを置くことを忘れないようにしましょう。

rk4はルンゲクッタ法を実装した関数です。

これは個人的な感想ですが、マルチエージェントなシミュレーションにおいては、すべてのエージェントを統括するような存在がいると便利な気がします。というわけでNeuronを統括するNetworkを構築します。

Networkはフィールドにニューロン数nを持ち、n個のNeuronをneuronsに収納します。また、独自にステップ数を数える変数countを使用します。

メソッドrunではn個のニューロンに1ステップ前の神経活動情報を与え、すべてのニューロン活動を検証します。本来は、ここで並列化が効果を発揮するでしょう。inputでは、すべてのニューロンに与えられる外部電流の大きさを変更します。

実験操作はmain関数に記述します。ここでは、dt=0.1でt=4000.0だけシミュレーションを行い、そのうち\(1000.0\leq t \leq 3000.0\)の時間帯は外部電流を増強します。

spike_trainが神経活動を記録する配列で、最終的な図表の作成とニューロンに与える入力を保存するために使用します。すなわち、

  1. Networkがspike_trainをNeuronに読ませる
  2. Neuronが0or1を返す
  3. Networkが各Neuronの返り値をspike_trainに渡す
  4. さっきの値をNetworkがNeuronに読ませる

という流れを取ります。

全ソースコードはこちら

https://github.com/doraneko94/Rust_Neuron

成果物

x軸が時間、y軸がニューロンのID、青い点が神経活動を表します。外部電流が小さかった時にはまばらだった活動が、外部電流を上げることで高頻度のバーストに変化していることがわかります。

ちなみに、個々のニューロンの膜電位は上図のように変化します。(\(200.0\geq t\geq 600.0\)の期間に外部電流を増強)

並列化するなら

今回はすべてのニューロンを直列で計算しましたが、Rustのメリットを活かすなら、これらの計算を並列化することが効果的です。並列化の手法としては、ニューロンを何分の1かずつに分割し、複数のスレッドで処理を行って最後に合わせるフォーク・ジョイン並列などが有名ですが、ここでは神経回路特有の工夫点について説明します。

この図はC. elegans(線虫)の全神経細胞を、それぞれが入力を受け取る他の神経細胞の数ごとに集計したものです。この図によると、ほとんどのニューロンは少ないニューロンとしか接続していませんが、一方で、非常に多くのニューロンと接続したニューロンがわずかに存在します。この分布は、y軸を適当に対数変換すると直線系になり、べき乗則というルールにしたがいます。つまり、「ほとんどは『持たざる者』だが、ごく少数の者は大量に持っている」というルールです。

このような状況下では、各ニューロンごとに仕事量が大きく異なるため、均等な割り付けをランダムに行うだけではあまり高速化できません。そのため、シナプス数などのパラメータを参考に、割り付けを工夫する必要があります。

このような現象は、神経回路のみならず、ネットワーク構造をとるものにおいて一般的に見られるものです。

結論

Rustはいいぞ。

Rustは高速で安全で並列化できていいぞ。

みんなで神経回路シミュレーションの共通言語をRustにしよう。

※Rustのインストール方法や基本文法については、各種Webページや書籍を参考にしてください。

※Rustのプロの方は、Rust的な書き方についてアドバイスをいただけると幸いです。

「UnityとC#で創る人工生命」見本誌配布

2019年11月23日15:57 見本誌 v0 → v1 のアップデートを行いました!

現在執筆中の「UnityとC#で創る人工生命」の見本誌を無料配布します。

この本では、UnityとC#を使って、自分で考えて行動する人工生命を観察することを目的としていますが、見本誌(1章のみ)ではUnityはまだ出てきません。

【その他、注意(そのうち更新されます)】

  1. 見本誌すらまだ途中段階です
  2. C#で画像認識するニューラルネットワークを作るソースコードは、すでに含まれています
  3. 上記ソースコードの解説はまだ途中です(おおむね終了しました)
  4. 図はまだありません
  5. その他、未完成部分が多々あります

 

更新情報をご希望の方は、

をフォローするか(更新情報以外のいらんことも呟きますが)、

お問い合わせフォームから、名前(ペンネーム可)、メールアドレスを入力し、題名に「人工生命本更新情報配信希望」とご記入のうえ、ご送信ください。

戻すDPーARC028D 注文の多い高橋商店(後編)

この記事ではAtCoderというサイトの問題を参考に、「戻す動的計画法(DP)」について解説する。

問題:ARC028 D – 注文の多い高橋商店

https://arc028.contest.atcoder.jp/tasks/arc028_4

この問題では、「同種のものを区別せずにM個選ぶ組み合わせ」と「戻すDP」についての考察が必要となる。「同種のものを区別せずにM個選ぶ組み合わせ」の解説は前回

同種のものを区別せずにM個選ぶ組み合わせーARC028D 注文の多い高橋商店(前編)

で行ったので、この記事を読む前に参照してほしい。そして、ここでは「戻すDP」について考えよう。これに該当する問題文は、以下の部分である。

いや、これは少し簡単過ぎるので、ちょっとした注文も追加しよう。整数 \(k, x\) からなる \(Q\) 個の注文を用意したので、それぞれについて「 \(k_{i}\) 種類目の商品をちょうど \(x_{i}\) 個選ばなければならないとき、合計 \(M\) 個の商品を選ぶ方法」の数を求めて下さい。

解説

考え方

\(N\) 番目の商品までを用いて \(M\) 個を選ぶ組み合わせについては、商品 \(A\) が \(4\) 個、 \(B\) が \(2\) 個、 \(C\) が \(3\) 個の場合、動的計画法を用いて以下のように表せることを前回示した。

これは、商品 \(A\) が \(0\) ~ \(4\) 個、 \(B\) が \(0\) ~ \(2\) 個、 \(C\) が \(0\) ~ \(3\) 個の範囲をそれぞれ動くことができるという条件で求められたものだが、ここで再び、とある1種類については個数を限定しなくてはならなくなる。すなわち、 \(N\) 番目の商品によって加算された組み合わせの数を”元に戻す”作業が必要なのであり、このアルゴリズムのことを「戻すDP」という。

アルゴリズムの詳細

目標

さて、ここで戻すDPを使って我々が得たいものは、以下のような表である。

ここで、2種類の表が出てきてややこしいので、前回の動的計画法によって得られた表を行列 \(dp\) 、今回これから戻るDPによって作成する表を行列 \(rdp\) として捉え、それぞれの \(i\) 行 \(j\) 列成分を \(dp_{ij}, rdp_{ij}\) と表すことにする。

そして、上の表を得たい理由は、「 \(n\) 番目の商品の数を \(k\) 個に固定したときに、全 \(N\) 種類の商品から \(M\) 個を選ぶ組み合わせの数」が \(rdp_{n, M-k}\) に等しいからである。上記の設定を用いると、例えば「 \(A\) の個数を \(2\) 個に固定したときに、全 \(3\) 種類の商品 \(A, B, C\) から \(5\) 個を選ぶ組み合わせの数」は、「 \(A\) の個数を \(0\) 個に固定したときに、全 \(3\) 種類の商品 \(A, B, C\) から \(5-2=3\) 個を選ぶ組み合わせの数」に等しいことを確認されたい。

漸化式

戻すDPにおいても通常の動的計画法と同様に漸化式の考え方を用いる。ここでは、「最後に \(n\) 番目の商品による組み合わせを加算して、全 \(N\) 種類の商品による組み合わせの数が求まった」と考えて、「全 \(N\) 種類の商品による組み合わせの数」、すなわち \(dp\) の最終行の値を用いて、 \(rdp\) の \(n\) 行目の値についての式を立てる。

前提条件として、前回求めた \(dp\) を計算するためのルール(漸化式)を再掲する。なお、ここでは前回 \(b_{n,m}\) とおいていた「 \(n\) 番目の商品までの中から \(m\) 個選ぶ組み合わせの数」を \(dp_{n,m}\) で置き換えている。

\(m<=a_{n}\) のとき

$$dp_{n,m} = \sum_{k=0}^{M}dp_{n-1,k} \tag{1}$$

\(m>a_{n}\) のとき

$$dp_{n,m} = \sum_{k=0}^{m}dp_{n-1,k}-\sum_{k=0}^{m-a_{n}-1}dp_{n-1,k} \tag{2}$$

ここで、「最後に \(n\) 番目の商品による組み合わせを加算する前の、 \(m\) 個選ぶ組み合わせの数」、すなわち、「 \(n\) 番目の商品を \(1\) つも選んでいない( \(0\) 個に固定されている)ときの、 \(m\) 個選ぶ組み合わせの数」が \(rdp_{n,m}\) で表されることを考えると、より一般的な式である \((2)\) 式を用いて

$$dp_{N,m} = \sum_{k=0}^{m}dp_{N-1,k}-\sum_{k=0}^{m-a_{n}-1}dp_{N-1,k}$$

$$dp_{N,m} = \sum_{k=0}^{m}rdp_{n,k}-\sum_{k=0}^{m-a_{n}-1}rdp_{n,k} \tag{3}$$

が成り立つ。なぜならば、前述の通り \(n\) 番目の商品を最後に加えたと考えるなら、 \(N-1\) 種類の商品は \(\{1, 2, … , n-1, n+1, … , N\}\) 番目の商品によって構成されるため、 \(dp_{N-1,k}=rdp_{n,k}\) が成り立つからである。

同様に、 \((2)\) 式を用いて添字をずらした式である

$$dp_{N,m-1} = \sum_{k=0}^{m-1}rdp_{n,k}-\sum_{k=0}^{m-a_{n}-2}rdp_{n,k} \tag{4}$$

が得られ、 \((3), (4)\) 式の辺々を引くと

$$dp_{N,m}-dp_{N,m-1} = \sum_{k=0}^{m}rdp_{n,k}-\sum_{k=0}^{m-1}rdp_{n,k}-\sum_{k=0}^{m-a_{n}-1}rdp_{n,k}-\sum_{k=0}^{m-a_{n}-2}rdp_{n,k}$$

$$dp_{N,m}-dp_{N,m-1} = rdp_{n,m}-rdp_{n,m-a_{n}-1}$$

となり、これを並び替えると

$$rdp_{n,m} = rdp_{n,m-a_{n}-1}+dp_{N,m}-dp_{N,m-1}$$

より、 \(rdp_{n,m}\) についての漸化式が得られる。以上の結果を変数の定義域を踏まえて一般化すると

\(m<=a_{n}\) のとき

$$rdp_{n,m} = dp_{N,m}-dp_{N,m-1}$$

\(m>a_{n}\) のとき

$$rdp_{n,m} = rdp_{n,m-a_{n}-1}+dp_{N,m}-dp_{N,m-1}$$

となる。したがって、通常の動的計画法によって前もって \(dp\) の最終行 \(dp_{N,.}\) の値を求めておくことで、上記のルールに従って \(rdp\) を埋めていけばよいことがわかる。

以上のように、一度通常の動的計画法により求めた値から”1歩前の状況に戻す”アルゴリズムのことを戻すDPという。

実装

前編・後編の結果をC++で実装すると以下のようになる。

通常の動的計画法で用いる \(\sum_{k=0}^{m}dp_{n-1,k}\) は、各 \(n\) において累積和 \(accum[m]\) として先に計算しておくと便利であり、かつ高速化できる。

そして、 \(rdp\) を求めた後は、入力 \(k_{i}, x_{i}\) に対して \(rdp[k][M-x]\) を解答とすればよい(理由については「解説」―「アルゴリズムの詳細」―「目標」を参照)。ただし、出力は解答を \(10^{9}+7\) で割った余りとする必要があるので、 \(dp\) や \(rdp\) を更新する都度、余りを求めておくとよい。このとき、減算( \(accum[m] – accum[m-a[n]-1]\) など)の結果は負になりうるので、適宜符号を正にしてから余りを計算する(MOD関数を作る)必要がある(1敗)。

※C++では、負の数に対する余りの計算結果は負として返される(例: \(-7\ \%\ 3 = -1\) )