ブログ

上記でまとめたメルセンヌ・ツイスタ(Mersenne Twister)法ではフロベニウスの同伴行列(companion matrix)が用いられます。当記事では同伴行列の定義と応用に関して取りまとめを行いました。

・数学まとめ
https://www.hello-statisticians.com/math_basic

同伴行列の概要

同伴行列(companion matrix)の定義

フロベニウスの同伴行列(companion matrix)は下記のように定められる。
$$
\large
\begin{align}
C_{n} = \left(\begin{array}{cccccc} 0 & 1 & 0 & \cdots & \cdots & 0 \\ 0 & 0 & 1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 & 1 & 0 \\ 0 & 0 & \cdots & \cdots & 0 & 1 \\ -c_0 & -c_1 & \cdots & \cdots & -c_{n-2} & -c_{n-1} \end{array} \right)
\end{align}
$$

同伴行列と固有多項式

下記のように$n$次多項式$p(\lambda)$を定める。
$$
\large
\begin{align}
p(\lambda) = c_0 + c_1 \lambda + c_2 \lambda^2 + \cdots + c_{n-1} \lambda^{n-1} + \lambda^{n}
\end{align}
$$

以下、$p(\lambda)=0$と$\det{(\lambda I_n – C_n)}=0$が同値であることを示す。$\det{(\lambda I_n – C_n)}$は余因子展開を用いて下記のように変形できる。
$$
\begin{align}
& \det{(\lambda I_n – C_n)} = \left| \begin{array}{cccccc} \lambda & -1 & 0 & \cdots & \cdots & 0 \\ 0 & \lambda & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda & -1 & 0 \\ 0 & 0 & \cdots & \cdots & \lambda & -1 \\ c_0 & c_1 & \cdots & \cdots & c_{n-2} & \lambda+c_{n-1} \end{array} \right| \\
&= \lambda (-1)^{1+1} \left| \begin{array}{cccccc} \lambda & -1 & 0 & \cdots & \cdots & 0 \\ 0 & \lambda & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda & -1 & 0 \\ 0 & 0 & \cdots & \cdots & \lambda & -1 \\ c_1 & c_2 & \cdots & \cdots & c_{n-2} & \lambda + c_{n-1} \end{array} \right| + c_0 (-1)^{n+1} \left| \begin{array}{ccccc} -1 & 0 & \cdots & \cdots & 0 \\ \lambda & -1 & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & \lambda & -1 & 0 \\ 0 & \cdots & \cdots & \lambda & -1 \end{array} \right| \\
&= \lambda (-1)^{1+1} \left| \begin{array}{cccccc} \lambda & -1 & 0 & \cdots & \cdots & 0 \\ 0 & \lambda & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda & -1 & 0 \\ 0 & 0 & \cdots & \cdots & \lambda & -1 \\ c_1 & c_2 & \cdots & \cdots & c_{n-2} & \lambda + c_{n-1} \end{array} \right| + c_0 (-1)^{n+1+(n-1)} \\
&= \lambda \left| \begin{array}{cccccc} \lambda & -1 & 0 & \cdots & \cdots & 0 \\ 0 & \lambda & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda & -1 & 0 \\ 0 & 0 & \cdots & \cdots & \lambda & -1 \\ c_1 & c_2 & \cdots & \cdots & c_{n-2} & \lambda + c_{n-1} \end{array} \right| + c_0 \\
&= \lambda^2 (-1)^{1+1} \left| \begin{array}{cccccc} \lambda & -1 & 0 & \cdots & \cdots & 0 \\ 0 & \lambda & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda & -1 & 0 \\ 0 & 0 & \cdots & \cdots & \lambda & -1 \\ c_2 & c_3 & \cdots & \cdots & c_{n-2} & \lambda + c_{n-1} \end{array} \right| + c_1 \lambda + c_0 \\
&= \lambda^3 (-1)^{1+1} \left| \begin{array}{cccccc} \lambda & -1 & 0 & \cdots & \cdots & 0 \\ 0 & \lambda & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda & -1 & 0 \\ 0 & 0 & \cdots & \cdots & \lambda & -1 \\ c_2 & c_3 & \cdots & \cdots & c_{n-2} & \lambda + c_{n-1} \end{array} \right| + c_2 \lambda^2 + c_1 \lambda + c_0 \\
&= \; \cdots \\
&= c_0 + c_1 \lambda + c_2 \lambda^2 + \cdots + c_{n-1} \lambda^{n-1} + \lambda^{n}
\end{align}
$$

上記より$p(\lambda)=0 \iff \det{(\lambda I_n – C_n)}=0$であることが確認できる。

参考

・Wikipedia：同伴行列
・線形空間と固有値の応用

整数の合同式に基づく線形合同法はシンプルでわかりやすい乱数生成方法である一方で、多次元分布を考える際に多次元疎結晶構造を生じるという課題があることに注意が必要です。当記事では線形合同法における多次元疎結晶構造をPythonなどを用いて図示を行いました。

・乱数生成の基本アルゴリズムまとめ
https://www.hello-statisticians.com/explain-terms-cat/random_sampling1.html

「自然科学の統計学」の$11$章「乱数の性質」の解説などを参考に作成を行いました。

自然科学の統計学 (基礎統計学)

3,190円(11/05 17:08時点)

Amazon

線形合同法における多次元疎結晶構造

線形合同法に基づく多次元の乱数の作成

$$
\large
\begin{align}
x_{n+1} = ax_{n} + c \quad (\mathrm{mod} \; M), \quad n \geq 1 \quad (1)
\end{align}
$$

上記の式で表される線形合同法を用いて多次元の乱数を作成することを考える。たとえば$M=16, a=5, c=1$の設定に対し、$x_0=0$を考えると下記のように乱数が生成される。
$$
\large
\begin{align}
x_1 &= 1, \, x_2 = 6, \, x_3 = 15, \, x_4 = 12, \, x_5 = 13, \\
x_6 &= 2, \, x_7 = 11, \, x_8 = 8, \, x_9 = 9, \, x_{10} = 14, \\
x_{11} &= 7, \, x_{12} = 4, \, x_{13} = 5, \, x_{14} = 10, \, x_{15} = 3, \\
x_{16} &= 0, \, x_{17} = 1, \, x_{18} = 6, \, x_{19} = 15, \, x_{20} = 12, \\
x_{21} &= 13, \, \cdots
\end{align}
$$

上記の周期は$M=16$に一致する。このように$M=2^n, a=4l+1, c=2l+1$に基づいて乱数を生成する場合、乱数の周期は$M=2^n$に一致するが、詳しくは下記で示した。

ここで$x_0=0, y_0=1$を元に$2$次元の一様乱数を生成することを考える。このとき、$y_0=x_1$であるので、乱数列${ x_n }$から隣接する$2$つの項を取り出し、$2$次元の乱数と見なしたと考えても同義である。

上記のように$2$つの乱数$(x_i,y_i=x_{i+1})$を生成し、$2$次元の一様分布と考えるというのは乱数生成における一般的な考え方である。一方で線形合同法でこのような方法で$2$次元の乱数を生成すると次節で取り扱う多次元疎結晶構造が生じるので注意が必要である。

多次元疎結晶構造

前節のように$M=16, a=5, c=1$の設定に基づいて初期値$x_0=0$から$(1)$式に基づいて乱数を生成した場合を考える。このとき、下記を実行することで$(x_i,x_{i+1})$を$2$次元平面状にプロットすることができる。

import numpy as np
import matplotlib.pyplot as plt

x = np.array([0., 1., 6., 15., 12., 13., 2., 11., 8., 9., 14., 7., 4., 5., 10., 3., 0.])

plt.scatter(x[:-1], x[1:])
plt.show()

・実行結果

上図のように線形合同法を用いて$2$次元乱数を生成すると生成される乱数はいくつかの直線上に規則的に並ぶ。同様の事象が$k$次元の乱数を考えた際にも生じ、このことを多次元疎結晶構造という。多次元疎結晶構造であるかどうかを調べるにあたってはスペクトル検定などが用いられる。

多次元疎結晶構造やスペクトル検定に関しては「自然科学の統計学」の$11$章付節の「多次元疎結晶構造とスペクトル検定」を参照すると詳しい。

自然科学の統計学 (基礎統計学)

3,190円(11/05 17:08時点)

Amazon

Pythonを用いた確認

一様乱数を生成する手法は様々ですが、乱数の生成の際には乱数の周期に注意が必要です。当記事では一様乱数の生成にあたっての原始的かつシンプルな手法である線形合同法に基づいて乱数を生成する際に、最大周期を実現するにあたってのパラメータの設定とその導出について取りまとめを行いました。

・乱数生成の基本アルゴリズムまとめ
https://www.hello-statisticians.com/explain-terms-cat/random_sampling1.html

前提の確認

線形合同法の概要

$$
\large
\begin{align}
x_{n+1} \equiv a x_{n} + c \quad (\mathrm{mod} \; M), \quad n \geq 1
\end{align}
$$

上記の式に基づいて一様乱数を生成する手法を線形合同法(LCG)という。詳しくは下記で取り扱った。

線形合同法による乱数列の最大周期

線形合同法では$\mathrm{mod} \; M$を考えるが$M$で割った余りは$0$から$M-1$の$M$通りである。よって線形合同法による乱数列の上限は$M$である。$M$はあくまで上限であり、たとえば$M=8,a=2,c=2$の場合は偶数しか生成されないので周期は$M$に一致しないことに注意が必要である。

逆に生成される乱数列が$M$に近い周期を持つ場合はすぐれた設定であると解釈できる。一方で$a=1,c=1$の場合は周期は$M$に一致するが、$1,2,3,4,5,\cdots$のように連番であり、良い乱数列であると言えない。このように$M,a,c$の設定は乱数生成の質に関わる重要な作業であると理解しておくと良い。

最大周期をもたらす設定とその証明

$M=2^n, \, a=4l+1, \, c=2m+1$

周期が$2^{n}$であることの証明

整数$l,m,n$を用いて$M=2^n, \, a=4l+1, \, c=2m+1$のように設定される場合に、生成される乱数列の周期は最大周期$M=2^n$に一致する。以下このことを数学的帰納法を用いて証明する。

i) $n=1,M=2^1=2$の場合
線形合同法における漸化式は下記のように表される。
$$
\large
\begin{align}
a x_{i} + c &= (4l+1) x_i + 2m + 1 \\
&= 2(2l x_i + m) + x_n + 1 \\
& \equiv x_i + 1 \equiv x_{i+1} \quad (\mathrm{mod} \; M)
\end{align}
$$

上記より$\mathrm{mod} \; M$で$x_{i+1} \equiv x_i + 1$が成立することより、$2$で割った余りに$0,1$が交互に出現すると考えることができる。よって$n=1$のとき$M=2^n, \, a=4l+1, \, c=2m+1$の設定における最大周期は$M$である。

ⅱ) $n=k,M=2^k$で成立する場合に$n=k+1,M=2^{k+1}$で成立するか
$2$で割った余りが$0,1$で循環する場合、$4$で割った余りが$0,0,0,\cdots$や$1,1,1,\cdots$にはなり得ない。同様に、「①$n=k,M=2^k$で成立する場合は、$\mathrm{mod} \; 2^{k+1}$を考えても$M=2^k$の周期よりも小さくなることはない」。

次に、「②$n=k,M=2^k$で成立する場合は、$x_i=j, \, 0 \leq j \leq 2^{k}-1$と$x_i=2^{k}+j$に対する$\mathrm{mod} \; 2^{k+1}$での$x_{j+1}$が異なる」ことを示す。
・$x_i=j$のとき
$$
\large
\begin{align}
a x_{i} + c = (4l+1)j + 2m + 1
\end{align}
$$

・$x_i=2^{k}+j$のとき
$$
\large
\begin{align}
a x_{i} + c &= (4l+1)(2^{k}+j) + 2m + 1 \\
&= 2^{k+2} + 2^{k} + (4l+1)j + 2m + 1 \\
& \equiv 2^{k} + (4l+1)j + 2m + 1 \quad (\mathrm{mod} \; 2^{k+1})
\end{align}
$$

$2^{k} + (4l+1)j + 2m + 1-((4l+1)j + 2m + 1)=2^{k}$であることから、$\mathrm{mod} \; 2^{k+1}$では「$x_i=j$と$x_i=2^{k}+j$に対する$x_{j+1}$が異なる」ことが示される。

①と②が成立するので「$n=k,M=2^k$で成立するとき$M=2^{k+1}$を考えると周期$2^k$の乱数列が並ぶか周期$2^{k+1}$の乱数列が並ぶかのどちらかである」と考えることができ、下記のように図示できる。

上図が成立することから、「③$x_i=0$のとき$x_{i+2^{2^{k}}} \neq x_i$」を示せば乱数列の周期が$2^{k+1}$であることを示すことができる。$x_i=0$のとき$x_{i+2^{2^{k}}}$は下記のように表せる。
$$
\large
\begin{align}
x_{i+2^{2^{k}}} – \frac{2m+1}{4l} & \equiv (4l+1)^{2^{k}} \left( x_0 – \frac{2m+1}{4l} \right) \\
&= – \frac{2m+1}{4l}(4l+1)^{2^{k}} \\
&= – \frac{2m+1}{4l} \sum_{j=0}^{2^{k}} {}_{2^{k}} C_{j} (4l)^{j} \\
x_{i+2^{2^{k}}} &= \frac{2m+1}{4l} – \frac{2m+1}{4l} \sum_{j=0}^{2^{k}} {}_{2^{k}} C_{j} (4l)^{j} \\
&= – \frac{2m+1}{4l} \sum_{j=1}^{2^{k}} {}_{2^{k}} C_{j} (4l)^{j} \\
& \equiv – \frac{2m+1}{4l} \times {}_{2^{k}} C_{1} (4l)^{1} \\
& \equiv -2^{k}(2m+1) \\
& \equiv 2^{k}
\end{align}
$$

よって③が成立する。①、②、③より「$n=k,M=2^k$周期が$2^{k}$が成立する場合に$n=k+1,M=2^{k+1}$で周期が$2^{k+1}$」が成立する。i)、ⅱ)より任意の自然数$n$に関して「$M=2^n, \, a=4l+1, \, c=2m+1$」の設定の周期が$2^n$であることが示される。

$M=2^n, \, a=4l+1, \, c=2m+1$の設定の注意点

前項の導出では数学的帰納法を用いたが、i)で確認したように$x_i$を$2$で割った余りが$0,1,0,1,0,\cdots$のように繰り返す。よってこの設定に基づいて生成される乱数は偶数と奇数が交互に入れ替わる乱数となる。

乱数の生成にあたっては可能な限り規則性のない数の並びを生成する必要があることから、「偶奇が必ず入れ替わる」のはそれほど良い乱数列ではないことに注意が必要である。

参考

・自然科学の統計学 Ch$.11$：乱数の性質

自然科学の統計学 (基礎統計学)

3,190円(11/05 17:08時点)

Amazon

乱数生成によく用いられるアルゴリズムのメルセンヌ・ツイスタ法ではメルセンヌ数(Mersenne Number)という概念が用いられます。当記事ではWikipediaなどを参考にメルセンヌ数、メルセンヌ素数、完全数に関して簡単に取りまとめを行いました。

・メルセンヌツイスタ法
https://www.hello-statisticians.com/explain-terms-cat/mersenne_twister1.html

メルセンヌ数・メルセンヌ素数

定義

$$
\large
\begin{align}
M_n = 2^{n}-1, \quad n \in \mathbb{N}
\end{align}
$$

上記のように自然数$n$を用いて$M_n=2^{n}-1$のように表される数をメルセンヌ数という。メルセンヌ数は下記のように具体的に考えることができる。
$$
\large
\begin{align}
2^{1}-1 &= 1 \\
2^{2}-1 &= 3 \\
2^{3}-1 &= 7 \\
2^{4}-1 &= 15 \\
2^{5}-1 &= 31 \\
2^{6}-1 &= 63 \\
2^{7}-1 &= 127 \\
2^{8}-1 &= 255 \\
2^{9}-1 &= 511 \\
2^{10}-1 &= 1023 \\
2^{11}-1 &= 2047
& \vdots
\end{align}
$$

上記に基づいて$M_n$に素数が出てくる場合をメルセンヌ素数(Mersenne prime)という。$M_n$が素数である場合、$n$は素数であるが、逆の「$n$が素数 $\implies$ $M_n$が素数」は成立しない。具体的には下記のように考えられる。

判定法

桁が小さい数に関しては計算した$M_n$を割り切る数があるかを探せば良いが、「リュカ-レーマー・テスト」を用いることなどによって効率化できる。

メルセンヌ素数の一覧

「Wikipedia：メルセンヌ素数の一覧」が詳しいが、$51$個発見されている。

$n \leq 10,000$では$n=2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, 521, 607, 1279, 2203, 2281, 3217, 4253, 4423, 9689, 9941$が挙げられる。

完全数

メルセンヌ数$M_n = 2^{n}-1$が素数であるとき、$2^{n-1}(2^{n}-1)$は完全数である。完全数は「正の約数の和に等しい自然数」を定めた数である。前節で確認したメルセンヌ素数$M_n = 2^{n}-1$に関して以下、$2^{n-1}(2^{n}-1)$が完全数であることを確かめる。

$n=2, M_n=3$

$$
\large
\begin{align}
2^{n-1}(2^{n}-1) &= 2^{2-1}(2^{2}-1) \\
&= 2 \times 3 = 6 \\
1 + 2 + 3 &= 6
\end{align}
$$

$n=3, M_n=7$

$$
\large
\begin{align}
2^{n-1}(2^{n}-1) &= 2^{3-1}(2^{3}-1) \\
&= 2^2 \times 7 = 28 \\
1 + 2 + 4 + 7 + 14 &= 28
\end{align}
$$

参考

・Wikipedia：メルセンヌ数
・Wikipedia：完全数
・Wikipedia：メルセンヌ・ツイスタ