情報幾何学を嗜む ~Bregmanダイバージェンスとその双対~

最近、情報幾何学の勉強をしている。情報幾何学は日本の甘利先生という方が切り開いてきた分野で、主には確率分布のパラメータが成す空間をリーマン多様体と捉えることで、確率分布族に対して幾何学的な解釈を与えるものである。

情報幾何学情報科学の一分野でありながら、微分幾何学の理解を要する難解なものである。はっきり言って、情報系の人間で可微分多様体やらリーマン計量やら接続やらを理解している人は一握りであろう。私も情報系、それも工学部の出身であるから、甘利先生の本を初めて手に取った修士2年のときは、あまりの難しさに一瞬で心が折れたのを覚えている。

しかし時は流れ、私も今ではわずかばかり数学の心が分かるようになってきた。そこで、いよいよこの難攻不落の要塞に攻めいってみようというわけである。

というわけで、本稿から始まるいくつかの記事の中で、情報幾何学における主要なトピックについて私が理解したところを書き連ねてみようと思う。本稿ではその第一歩として、Bregmanダイバージェンスとその双対ダイバージェンスについて考えてみる。

ダイバージェンス

情報幾何学を語る上でダイバージェンスの存在は外せない。ダイバージェンスの定義を[1]から引用する*1*2

次の3条件を満たす2点関数 D[P : Q]ダイバージェンスと呼ぶ.
1)  D[P : Q] \ge 0.
2)  P = Qのとき, このときに限り,  D[P : Q] = 0.
3)  P点と Q点が近いとし, それぞれの座標を,  {\boldsymbol \xi},\ {\boldsymbol \xi} + d{\boldsymbol \xi}とする. このとき,  D[{\boldsymbol \xi} : {\boldsymbol \xi} + d{\boldsymbol \xi}]テイラー展開すると,
 \displaystyle{
D[{\boldsymbol \xi} : {\boldsymbol \xi} + d{\boldsymbol \xi}] = \frac{1}{2} \sum g_{ij}({\boldsymbol \xi}) d\xi_i d\xi_j
}
と2次の項が最初に出るが, 行列
 \displaystyle{
G({\boldsymbol \xi}) = (g_{ij}({\boldsymbol \xi}))
}
は正定値対称である.

ただし、 {\boldsymbol \xi}は有限次元ベクトルである。

ダイバージェンスは距離の公理を満たしていない。すなわち、一般には D[P : Q] = D[Q : P]とならない。この醜い非対称性が後に華麗な蝶へと変貌を遂げるのであるが、それは双対ダイバージェンスこところで説明する。

Bregmanダイバージェンス

ダイバージェンスの中でも特に重要なものの1つにBregmanダイバージェンスがある。これは滑らかな狭義凸関数 \psi({\boldsymbol \xi})*3を用いて以下のように定義される[1]。

まず、点 {\boldsymbol \xi}'における接超平面の方程式は以下のようになる。

 \displaystyle{
z = \psi({\boldsymbol \xi}') + \nabla \psi({\boldsymbol \xi}') ({\boldsymbol \xi} - {\boldsymbol \xi}')
}

 {\boldsymbol \xi}において、この接超平面と元の関数 \psi({\boldsymbol \xi})の差は以下のようになる。

 \displaystyle{
D[{\boldsymbol \xi} : {\boldsymbol \xi}'] = \psi({\boldsymbol \xi}) - \psi({\boldsymbol \xi}') - \nabla \psi({\boldsymbol \xi}') ({\boldsymbol \xi} - {\boldsymbol \xi}')
}

これを凸関数 \psiから導かれる {\boldsymbol \xi}から {\boldsymbol \xi}'へのBregmanダイバージェンスと呼ぶ。

Bregmanダイバージェンスダイバージェンスの定義を満たすことは自明ではないため、本来であれば証明すべきである。実際、ノートで計算して確かめることは出来たのだが、それをブログに書き起こす気力と時間がなくなってしまったため、ここでは割愛する。

参考までに計算の指針だけ述べておく。まず、 D[ {\boldsymbol a} : {\boldsymbol \xi} ]について {\boldsymbol \xi} = {\boldsymbol a}における多変数のTaylor展開を計算する。0次の項は D[{\boldsymbol a} : {\boldsymbol a}] = 0である。1次の項も計算すると0になる。2次の項は計算すると \psiのHesse行列に等しくなる。 \psiは滑らかな凸関数と仮定しているため、これは正定値対称となる。

Legendre変換による双対空間

Bregmanダイバージェンスを定義するために凸関数が登場した。凸関数といえば皆さん何を思い浮かべるだろうか?いろいろあると思うが、凸関数にまつわる重要な概念としてLegendre変換が挙げられる。Legendre変換を行うことで、凸関数が定義された空間の双対空間、及び双対凸関数を得ることができる。

最終的には多変数の場合を考える必要があるが、まずは1変数の場合から考えてみよう。

1変数凸関数のLegendre変換

1変数の滑らかな凸関数 \psi(x)を考える。滑らかな凸関数の導関数は異なる xに対して必ず異なる値を取る。逆に適当な実数 pを与えると、それを導関数の値とするような点 xが一意に決まる。

すなわち、 \psi'(x)の定義域を S、値域を Dとすると、 \psi' : S \ni x \mapsto p \in D全単射となる。 Dのことを双対空間と呼ぶ。

ここで、 \psi(x)に以下のような変換を施すことで双対空間に対して新たな関数 \psi^{*}(p)を定める事ができる。

 \displaystyle{
\psi^{*}(p) = \max_{x}(px - \psi(x))
}

これをLegendre変換と呼ぶ。

右辺の px - \psi(x)を最大にする xについて考えてみよう。最大値を与える xにおいては、この式を x微分したものが0となる必要がある。すなわち、以下が成立する。

 \displaystyle{
\begin{eqnarray}
\frac{d}{dx}(px - \psi(x)) &=& 0 \\
p - \psi'(x) &=& 0 \\
p &=& \psi'(x)
\end{eqnarray}
}

すなわち、 \psi'(x)の値が pとなるような xにおいて px - \psi(x)は最大となる。これはつまり、 Dの元と一対一に対応する Sの元を選べば良いということを意味する。

詳細は後述するが、実は \psi^{*}(p)も凸関数であり、これを双対凸関数と呼ぶ。そのため、 \psi^{*}(p)に対して再度Legendre変換を施すことができるが、その結果は元の関数 \psi(x)と一致する[2]。つまり、Legendre変換の逆変換はLegendre変換そのものである。これより、双対性というのはあくまで相対的な概念に過ぎないことが分かる。

双対凸関数が凸関数であることの証明

双対凸関数が凸関数であることは自明ではないので、以下で証明してみよう。 \psi^{*}(p)の定義式において x pの関数であると考えると、2階導関数は以下のようになる。

 \displaystyle{
\begin{eqnarray}
\frac{\partial^2 \psi^{*}}{\partial p^2} &=& \frac{\partial^2}{\partial p^2} (px - \psi(x)) \\
&=& \frac{\partial}{\partial p} (x + px' - \psi'(x) x') \\
&=& 2x' - \psi''(x){x'}^2 + x''(p - \psi'(x))
\end{eqnarray}
}

 p = \psi'(x)を代入すると以下のようになる。

 \displaystyle{
\frac{\partial^2 \psi^{*}}{\partial p^2} = 2x' - \psi''(x){x'}^2
}

 p = \psi'(x)の両辺をpで微分して整理すると x' = \frac{1}{\psi''(x)}となる。これを代入すると以下のようになる。

 \displaystyle{
\frac{\partial^2 \psi^{*}}{\partial p^2} = \frac{1}{\psi''(x)}
}

2階導関数が瞬間的にでも0になることがあるような関数は面倒なので除外して考えると、 \psi''(x)は上に凸なら常に負、下に凸なら常に正となる。よって \psi^{*}{''}(p)の符号も一定となるため、 \psi^{*}(p)は凸関数である。

多変数凸関数のLegendre変換

少々くどいかもしれないが、1変数のときと同じ議論を多変数についても行ってみよう。

2つ以上の変数を持つ滑らかな凸関数 \psi({\boldsymbol \xi})を考える。滑らかな凸関数の勾配ベクトルは異なる {\boldsymbol \xi}に対して必ず異なるベクトルとなる。逆に適当な実数値ベクトル {\boldsymbol \xi}^{*}を与えると、それを勾配ベクトルとするような点 {\boldsymbol \xi}が一意に決まる。

この対応関係により、定義域と \psi({\boldsymbol \xi})の勾配ベクトルが取り得る値の間に一対一の対応関係が得られる。 {\boldsymbol \xi}^{*}が成す空間のことを双対空間と呼ぶ。

すなわち、 \nabla \psi({\boldsymbol \xi})の定義域を S、値域を Dとすると、 \nabla \psi : S \ni {\boldsymbol \xi} \mapsto {\boldsymbol \xi}^{*} \in D全単射となる。 Dのことを双対空間と呼ぶ。

ここで、 \psi({\boldsymbol \xi})に以下のような変換を施すことで双対空間に対して新たな関数 \psi^{*}({\boldsymbol \xi}^{*})を定める事ができる。

 \displaystyle{
\psi^{*}({\boldsymbol \xi}^{*}) = \max_{{\boldsymbol \xi}}({\boldsymbol \xi} \cdot {\boldsymbol \xi}^{*} - \psi({\boldsymbol \xi}))
}

これをLegendre変換と呼ぶ。

右辺の {\boldsymbol \xi} \cdot {\boldsymbol \xi}^{*} - \psi({\boldsymbol \xi})を最大にする {\boldsymbol \xi}について考えてみよう。最大値を与える {\boldsymbol \xi}においては、勾配ベクトルが零ベクトルとなる必要がある。すなわち、以下が成立する。

 \displaystyle{
\begin{eqnarray}
\nabla ({\boldsymbol \xi} \cdot {\boldsymbol \xi}^{*} - \psi({\boldsymbol \xi})) &=& 0 \\
{\boldsymbol \xi}^{*} - \nabla \psi({\boldsymbol \xi}) &=& 0 \\
{\boldsymbol \xi}^{*} &=& \nabla \psi({\boldsymbol \xi})
\end{eqnarray}
}

すなわち、 \nabla \psi({\boldsymbol \xi})の値が {\boldsymbol \xi}^{*}となるような {\boldsymbol \xi}において {\boldsymbol \xi} \cdot {\boldsymbol \xi}^{*} - \psi({\boldsymbol \xi})は最大となる。これはつまり、 Dの元と一対一に対応する Sの元を選べば良いということを意味する。

証明は大変そうなので諦めるが、1変数の場合と同じく \psi^{*}({\boldsymbol \xi}^{*})も凸関数であり、これを双対凸関数と呼ぶ。そのため、 \psi^{*}({\boldsymbol \xi}^{*})に対して再度Legendre変換を施すことができるが、その結果が元の関数 \psi({\boldsymbol \xi})と一致するというのも1変数の場合と同様である。

双対ダイバージェンス

双対凸関数は凸関数なので、これを用いると双対空間にもBregmanダイバージェンスを定義できる。

 \displaystyle{
D^{*}[{\boldsymbol \xi}^{*} : {\boldsymbol \xi}'^{*}] = \psi^{*}({\boldsymbol \xi}^{*}) - \psi^{*}({\boldsymbol \xi}'^{*}) - \nabla \psi^{*}({\boldsymbol \xi}'^{*}) ({\boldsymbol \xi}^{*} - {\boldsymbol \xi}'^{*})
}

これを双対ダイバージェンスと呼ぶ。

元のダイバージェンスとの関係

双対ダイバージェンスと元のダイバージェンスとの間には重要な関係がある。以下でそれを導いてみよう。

 {\boldsymbol \xi}, {\boldsymbol \xi}'について、双対空間において対応する点がそれぞれ {\boldsymbol \xi}^{*}, {\boldsymbol \xi}'^{*}であるとする。このとき以下の式が成り立つ。

 \displaystyle{
\begin{eqnarray}
\psi^{*}({\boldsymbol \xi}^{*}) &=& {\boldsymbol \xi} \cdot {\boldsymbol \xi}^{*} - \psi({\boldsymbol \xi}) \\
\psi^{*}({\boldsymbol \xi}'^{*}) &=& {\boldsymbol \xi}' \cdot {\boldsymbol \xi}'^{*} - \psi({\boldsymbol \xi}') \\
\nabla \psi^{*}({\boldsymbol \xi}'^{*}) &=& {\boldsymbol \xi}'
\end{eqnarray}
}

これらを双対ダイバージェンスの式に代入すると以下のようになる。

 \displaystyle{
\begin{eqnarray}
D^{*}[{\boldsymbol \xi}^{*} : {\boldsymbol \xi}'^{*}] &=& ({\boldsymbol \xi} \cdot {\boldsymbol \xi}^{*} - \psi({\boldsymbol \xi})) - ({\boldsymbol \xi}' \cdot {\boldsymbol \xi}'^{*} - \psi({\boldsymbol \xi}')) - {\boldsymbol \xi}' ({\boldsymbol \xi}^{*} - {\boldsymbol \xi}'^{*}) \\
&=& \psi({\boldsymbol \xi}') - \psi({\boldsymbol \xi}) - {\boldsymbol \xi}^{*}({\boldsymbol \xi}' - {\boldsymbol \xi})
\end{eqnarray}
}

これに  {\boldsymbol \xi}^{*} = \nabla \psi({\boldsymbol \xi})を代入すると以下のようになる。

 \displaystyle{
\begin{eqnarray}
&& \psi({\boldsymbol \xi}') - \psi({\boldsymbol \xi}) - \nabla \psi({\boldsymbol \xi}) ({\boldsymbol \xi}' - {\boldsymbol \xi}) \\
&=& D[{\boldsymbol \xi}' : {\boldsymbol \xi}]
\end{eqnarray}
}

結局、以下の式が得られた。

 \displaystyle{
D^{*}[{\boldsymbol \xi}^{*} : {\boldsymbol \xi}'^{*}] = D[{\boldsymbol \xi}' : {\boldsymbol \xi}]
}

ダイバージェンスの定義を説明した際、ダイバージェンスは対称性を満たさないということを述べた。しかし、上式が示す通りBregmanダイバージェンスについては2つの引数を入れ替えたものは双対ダイバージェンスに一致するのである。元の空間だけでは対称性がないように見えるが、双対空間まで広げて考えるとこのように美しい対称性が現れるというのは非常に面白い。

ナブラを使わない表現方法

Bregmanダイバージェンスの定義式にはナブラ ( \nabla) が含まれており少々複雑である。実はこれはちょっとした式変形で回避できる。

これまでの議論が追えていれば簡単なので、以下に式変形だけ示す。

 \displaystyle{
\begin{eqnarray}
D[{\boldsymbol \xi} : {\boldsymbol \xi}'] &=& \psi({\boldsymbol \xi}) - \psi({\boldsymbol \xi}') - \nabla \psi({\boldsymbol \xi}') ({\boldsymbol \xi} - {\boldsymbol \xi}') \\
&=& \psi({\boldsymbol \xi}) - \psi({\boldsymbol \xi}') - {\boldsymbol \xi}'^{*} ({\boldsymbol \xi} - {\boldsymbol \xi}') \\
&=& \psi({\boldsymbol \xi}) + ({\boldsymbol \xi}' \cdot {\boldsymbol \xi}'^{*} - \psi({\boldsymbol \xi}')) - {\boldsymbol \xi}'^{*} \cdot {\boldsymbol \xi} \\
&=& \psi({\boldsymbol \xi}) + \psi^{*}({\boldsymbol \xi}'^{*}) - {\boldsymbol \xi}'^{*} \cdot {\boldsymbol \xi}
\end{eqnarray}
}

まとめ

本稿では情報幾何学のトピックのうち、Bregmanダイバージェンスとその双対ダイバージェンスに関する事柄について述べた。ダイバージェンスは対称性を持たないが、BregmanダイバージェンスについてはLegendre変換による双対空間まで考えることで美しい対称構造が得られることを確認した。

本稿ではまだ微分幾何学らしい概念は登場しなかった。つまり、ここで述べたことは情報幾何学の中ではまだまだ序の口ということである。次回以降、少しずつ幾何学的な内容に踏み込んでいきたいと思う。

*1:甘利先生の本[1]ではダイバージェンス微分可能性などに触れられないままいきなりTaylor展開しているところがもやもやする。あまり細かい数学的議論に重きを置いた本ではないので、これについては滑らかな関数であり、かつ剰余項は収束すると仮定を置いてしまうしかないのだろう。

*2:ダイバージェンスの引数に点を入れたり点の座標を入れたりと記号がぶれているが、本[1]に合わせた結果なので好意的に解釈して頂けるとありがたい。

*3:本[1]ではBregmanダイバージェンスの定義に用いる凸関数の性質について厳密な条件が記載されていない。議論を簡単にするために、ここでは滑らかな狭義凸関数であるとした。

クーポンコレクター問題の確率分布を解き明かす

クーポンコレクター問題というものをご存知だろうか?これは、例えば6種類のおもちゃが出るガシャポン*1があったとして、何回くらい引けば全種類引き当てる事ができるか?というようなことを考える問題である。

この問題に対して、平均や分散がどうなるかということは非常によく語られることである[1]。また、時として不等式による評価について議論している記事を見かけることもある[2]。

しかし、肝心の確率分布については議論される事が非常に少ない。あまりにも記事を見かけないので、私は当初クーポンコレクター問題の確率分布を求めることは不可能なのではないかとさえ思っていた。

それでもめげずに調べ続けた結果、私はついに確率分布について結論を出しているページを見つけた。本稿ではそれを紹介し、喜びを分かち合いたいと思う。

第2種スターリング数

クーポンコレクター問題の確率分布を求めるためには、第2種スターリング数について理解しておく必要がある。少し長いが、定義をWikipediaから引用する[3]。

第2種スターリング数
第2種スターリング数 (Stirling number of the second kind)  \{{\textstyle {n \atop k}}\}は、 x^{n}を下降階乗冪 x^{\underline {k}}\equiv x\,(x-1)(x-2)\cdots (x-k+1)級数:

 \displaystyle{
x^{n}=\sum _{k=0}^{n}\left\{{n \atop k}\right\}\,x^{\underline {k}}
}

で展開したときの展開係数として定義される。この定義では、 0\leq k\leq nである。便宜上、 \{{\textstyle {0 \atop 0}}\}=1と定義する。第2種スターリング数は

 \displaystyle{
\left\{{n \atop k}\right\}=\left\{{n-1 \atop k-1}\right\}+k\,\left\{{n-1 \atop k}\right\}
}

なる漸化式で計算できる。

定義も大切なのだが、第2種スターリング数は以下のように特徴付けられるということが重要である[3]。

第2種スターリング数の特徴付け
第2種スターリング数 \{{\textstyle{n \atop k}}\}は、組合せ数学において、番号づけされた n個の要素をグループ k個に分割する組み合わせの数を与える。分割する要素は番号付けされているので個別に区別できるが、グループは順序を特に区別しないものとする。

クーポンコレクター問題の確率分布

準備が整ったので、確率分布を求めてみよう。以下の議論は全て[4][5]を参考にした。

集めるクーポンの種類を k種類とし、初めて全種類のクーポンを取得できるまでの試行回数を表す確率変数を Nとする。この時、 n回目の試行で初めて全種類のクーポンを取得できる確率 P(N = n)を求める。

 n回の試行により出現し得るクーポンの出方の総数は k^nである。もし n回目の試行で初めて全種類のクーポンが出るパターンの総数が分かれば、それを k^nで割ったものが求める確率である。

 n回目の試行で初めて全種類のクーポンが出揃うということは、 n-1回目の時点で k-1種類のクーポンがすでに出ている必要がある。このパターンの総数を求めるために、 n-1回目までの試行に対して1から順に番号を割り振る。そして、これらの番号付けられた試行を k-1個のグループに分ける。同じグループに分けられた試行については同じクーポンが得られたと考える。この分け方のパターンの総数は第2種スターリング数 \{\textstyle{n-1 \atop k-1}\}となる。

第2種スターリング数は分けられたグループの間の順序は区別しないが、今はクーポンの種類は区別されるので、 (k-1)!を掛ける。さらに、 n-1回目までに出る k-1種類のクーポンは k種類のうちどれであるかは問わないので、全てのパターンをカウントする必要がある。そのため、 \left(\textstyle{k \atop k-1}\right) = kを掛ける。結局、 n-1回目の時点で k-1種類のクーポンが出現するパターンの総数は k! \{\textstyle{n-1 \atop k-1}\}となる。

 n回目の試行でまだ出ていない最後のクーポンが出るパターンの総数は1なので、 n回目の試行で初めて全種類のクーポンが出るパターンの総数は k! \{\textstyle{n-1 \atop k-1}\}となる。以上の議論により次の式を得る。

 \displaystyle{
P(N = n) = \frac{k!}{k^n} \{\textstyle{n-1 \atop k-1}\}
}

計算してみよう

 n, kにいくつか具体的な数値を入れて P(N = n)を計算してみよう。手計算ではやっていられないので、pythonを利用する。第2種スターリング数はsympyのstirling関数を用いて計算することが出来るので、それを利用する[6]。数値だけ示しても分かり辛いので、結果をmatplotlibでグラフ化した。Python自体はpydroid3を使って実行した[7]。なお、sympyなどの必要なパッケージはあらかじめpipでインストールしておく必要がある。

計算に用いたソースコードを以下に示す。

from sympy.functions.combinatorial.numbers import stirling
import matplotlib.pyplot as plt
import math

startCouponKind = 2
numCouponKind = 6
numTrial = 30
for k in range(startCouponKind, numCouponKind+1):
  x = []
  y = []
  for n in range(1, numTrial+1):
    stir = stirling(n-1, k-1, kind = 2)
    prob = math.factorial(k) * stir / (k**n)
    x.append(n)
    y.append(prob)
  plt.plot(x, y, label='k = {}'.format(k))
plt.legend()
plt.savefig('/storage/emulated/0/coupon.png')

得られたグラフを以下に示す。

f:id:peng225:20190314084353j:plain

途中でピークを持ち、右側に裾野が広がっている様子が見て取れる。 nが小さいときに確率が0となる領域があるが、これはクーポンの種類よりも試行回数が少ないときの様子が表れているものである。

また、 kの増え方に対して、ピークの位置が右に移動していく速度の方がやや速いように見える。実際、全種類コンプリートするまでの試行回数の期待値は O(k \mathrm{log} k)なので、これは理論的にも辻褄が合っている[1]。

まとめ

本稿ではクーポンコレクター問題の確率分布を明らかにし、実際に計算を行った結果を示した。クーポンコレクター問題ではパターンを数え上げる事が難しかったが、その難しさを既知の概念である第2種スターリング数に押し込めることで、理論的にスッキリとした結論を得ることができた。

これまで組み合わせ数学には興味がなかったが、本稿を書き上げるうちにその面白さを垣間見る事ができた。またいずれ体系的に勉強しよう。

*1:これの呼び方は地方によって差があるような気がするが、自分の流儀で呼ばせて頂く。

いくつかのLie群がLie群であることを定義に戻って確かめる

Lie 群は難しい。この理由の1つは、議論の前提となる領域が広いことにあると思われる。Lie群とは群であり多様体であるような数学的対象である。そのため、定義を理解するだけで群論多様体の知識が求められる。また、Lie群の教科書で最初に扱われるような基本的なLie群は行列群である。しかも、そのコンパクト性に着目した議論も多い。そのため、線形代数位相空間の基礎的な事項も理解しておくことが望ましい。

繰り返すが、Lie群は難しい。私はここ最近Lie群を勉強し始めて、この事実を痛感している。こういう時は足元を一歩ずつ踏み固めて行くしかない。その一環として、本稿ではいくつかの基本的なLie群について、それらが本当にLie群になっていることを定義に照らし合わせて確認してみる。

本稿では私の独断で以下の2つのLie群を扱う。

準備

Lie群の定義

Lie群の定義を[1]より引用する。

多様体 Gが群構造を持ち、その群演算
 \displaystyle{
G \times G \to G;\ (x, y) \mapsto xy,
\ \ G \to G;\ x \mapsto x^{-1}
}
 C^{\infty}写像になるとき、 GLie群と呼ぶ。

多様体上の写像 C^{\infty}級であるということ

Lie群の定義の中で C^{\infty}写像という言葉が出てきた。この定義を[2]より引用する。

1点においてCs
連続写像 f: M \to Nが, 1点 p \in Mにおいて C^s級であるとは,  pを含む M C^r級座標近傍 (U; x_1, \cdots , x_m) f(p)を含む N C^r級座標近傍 (V; y_1, \cdots , y_n)が存在して,
(1)  f(U) \subset Vかつ
(2)  (U; x_1, \cdots , x_m) (V; y_1, \cdots , y_n)に関する fの局所座標表示が C^s級である,
この2つの条件がなりたつことである. ただし,  1 \le s \le r.
Cs写像
 f: M \to N C^s写像 (mapping of class  C^s) であるとは,  Mの各点 pにおいて f C^s級であることである.

上記定義において s = \inftyとすれば C^{\infty}写像の定義となる。

部分多様体に関する定理

後で使う定理について説明するために、正則点・臨界点、及び正則値・臨界値の定義を述べておく[2]。

正則点・臨界点
 Mの点 pにおける f微分
 \displaystyle{
(df)_p : T_p (M) \to T_{f(p)} (N)
}
が‘上へ’の線形写像であるとき,  p fの正則点 (regular point) とよぶ.  (df)_pが‘上へ’の線形写像でないとき,  p fの臨界点 (critical point) という.
正則値・臨界値
 f : M \to Nの臨界点全部の集合 (臨界点集合) を C_fで表す ( C_f Mの部分集合である) .
 C_fの像 f(C_f)に属する Nの点 q fの臨界値 (critical value) とよぶ. 臨界値でない Nの点を fの正則値 (regular value) という.

以下の定理を後の議論で使用する[2]。

部分多様体に関する定理
 q \in N C^r写像 f: M \to Nの正則値で,  f^{-1}(q) \neq \phiであるとすると, 逆像 f^{-1}(q) M (m-n)次元 C^r級部分多様体である.

Lie群であることの確認

以上で準備が整ったので、Lie群であることの確認に移る。本稿では群であることの確認はサボり、それぞれのLie群について以下の3点を確認した。

  •  C^{\infty}多様体であること
  • 群演算が C^{\infty}写像であること
  • 逆元を取る操作が C^{\infty}写像であること

一般線形群 \mathrm{GL}(n, \mathbb{R})

 C^{\infty}多様体であること

 \mathbb{R} n次正方行列全体の集合を  \mathrm{M}(n, \mathbb{R})と書く。 \mathrm{GL}(n, \mathbb{R}) \mathrm{M}(n, \mathbb{R})の部分集合のうち、以下のように表されるものである。

 \displaystyle{
\{A \in \mathrm{M}(n, \mathbb{R}) | \mathrm{det} (A) \neq 0\}
}

実は、 \mathrm{GL}(n, \mathbb{R}) \mathrm{M}(n, \mathbb{R})の開部分集合となる。詳細は[3]の命題1.17に譲り、ここでは概要だけ説明する。まず、 \mathrm{det} :  \mathrm{M}(n, \mathbb{R}) \to \mathbb{R}連続写像である。このとき、 \mathrm{GL}(n, \mathbb{R}) = \mathrm{det}^{-1}\left(\{x \in \mathbb{R} | x \ne 0\}\right)と書ける。 \{x \in \mathbb{R} | x \ne 0\} \mathbb{R}の開集合なので、連続写像 \mathrm{det}による逆像も開集合となる。

ここで、 \mathrm{M}(n, \mathbb{R})に属する行列の各成分を座標と見なすと、これは \mathbb{R}^{n^2}と同一視できる。そのため、
 \mathrm{M}(n, \mathbb{R}) C^{\infty}多様体となる。その開部分集合も C^{\infty}多様体となるので、 \mathrm{GL}(n, \mathbb{R}) C^{\infty}多様体となる。このような多様体を開部分多様体と呼ぶ[2]。

群演算が C^{\infty}写像であること

 \mathrm{GL}(n, \mathbb{R})の群演算は通常の行列の掛け算である。この写像 \mathrm{mul} : \mathrm{GL}(n, \mathbb{R}) \times  \mathrm{GL}(n, \mathbb{R}) \to \mathrm{GL}(n, \mathbb{R})と表す。以下で \mathrm{mul} C^{\infty}級であることを示す。

 \mathrm{GL}(n, \mathbb{R})はそれ自身 \mathbb{R}^{n^2}の開集合と見なせるため、そのアトラスとして自分自身だけを座標近傍として含む、要素数1の集合族を取れる。  \mathrm{GL}(n, \mathbb{R}) \times  \mathrm{GL}(n, \mathbb{R})についても、積多様体のアトラスの定義を考えると、やはりアトラスとして自分自身だけを要素として含む集合族を取れる。よって \mathrm{GL}(n, \mathbb{R})  \mathrm{GL}(n, \mathbb{R}) \times  \mathrm{GL}(n, \mathbb{R})は共に C^{\infty}多様体である。

任意の点 (A, B) \in \mathrm{GL}(n, \mathbb{R}) \times  \mathrm{GL}(n, \mathbb{R})について、 \mathrm{mul}( (A, B) ) = ABの各成分は A, Bの成分の多項式となる。よって \mathrm{mul} (A, B)において C^{\infty}級である。 (A, B)は任意なので、 \mathrm{mul} C^{\infty}写像である。

逆元を取る操作が C^{\infty}写像であること

 \mathrm{GL}(n, \mathbb{R})の元の逆元を取る操作とは、逆行列を求める操作を意味する。この写像 \mathrm{inv} : \mathrm{GL}(n, \mathbb{R}) \to \mathrm{GL}(n, \mathbb{R})と表す。以下で \mathrm{inv} C^{\infty}級であることを示す。

任意の点 A \in \mathrm{GL}(n, \mathbb{R})について、 \mathrm{inv}(A) = A^{-1}である。逆行列は余因子行列を行列式で割ることによって得られるため、 A^{-1}の各成分は Aの成分の有理式となる。よって \mathrm{inv} Aにおいて C^{\infty}級である。 Aは任意なので、 \mathrm{inv} C^{\infty}写像である。

 \mathrm{inv}の定義域は \mathrm{GL}(n, \mathbb{R})なので、行列式は常に0ではないことに注意されたい。

直交群 \mathrm{O}(n)

 C^{\infty}多様体であること

これを示すのは思いのほか難しいため、証明のアウトラインだけ述べることにする。詳細は[4]を参照されたい。

 n次実対称行列全体の集合を \mathrm{S}(n)とする。 \mathrm{S}(n) \mathbb{R}^{n(n+1)/2}と同一視できるため、 n(n+1)/2 C^{\infty}多様体である。

この時、以下のような写像を考える。

 \displaystyle{
f : \mathrm{M}(n, \mathbb{R}) \to \mathrm{S}(n), A \mapsto AA^T
}

 f C^{\infty}写像である。 Aが直交行列のとき、その逆行列 A^Tとなる。そのため、 \mathrm{O}(n) = f^{-1}(I)となる。ここで、 I単位行列である。

もし I fの正則値であれば、先ほど提示した定理により \mathrm{O}(n) n(n-1)/2 C^{\infty}多様体となる。そのためには、 \mathrm{O}(n)の全ての点における微分全射になれば良い。具体的に微分計算を行い、それが全射であることを確かめるアプローチになるが、体力の限界なのでこれ以降は[4]にお任せする。

群演算が C^{\infty}写像であること

 \mathrm{O}(n)の群演算は  \mathrm{GL}(n, \mathbb{R})における群演算を  \mathrm{O}(n) \times  \mathrm{O}(n)に制限した写像 \mathrm{mul} |_{\mathrm{O}(n) \times \mathrm{O}(n)}によって与えられる。これは包含写像 i : \mathrm{O}(n) \times \mathrm{O}(n) \to \mathrm{GL}(n, \mathbb{R}) \times  \mathrm{GL}(n, \mathbb{R}) \mathrm{mul}によって以下のように書ける。

 \displaystyle{
\mathrm{mul} |_{\mathrm{O}(n) \times \mathrm{O}(n)} = \mathrm{mul} \circ i
}

 i \mathrm{mul}はどちらも C^{\infty}級であるため、それらの合成写像 C^{\infty}級となる。よって \mathrm{O}(n)の群演算は C^{\infty}写像である。

逆元を取る操作が C^{\infty}写像であること

群演算とほぼ同じように制限写像を考えれば良い。詳細は割愛する。

まとめ

本稿では一般線形群と直交群が共にLie群であることを確認した。どちらも非常に有名なLie群でありながら、それを示すのは意外と大変だった。Lie群を理解するまでの道のりは険しい。

今こそJordan標準形と向き合う

線形代数を勉強して、Jordan標準形という言葉を耳にしたことがない人はいないだろう。Jordan標準形とは、ざっくり言えば行列の対角化を一般化したようなものである。行列の対角化はいつでもできるとは限らず、報われない行列たちが存在する。一方、Jordan標準形は常に存在することが知られており、まさに行列界の救世主と言える。

私もかつて大学院入試の際にJordan標準形について勉強したことはある。しかし、専攻が情報系というのもあって、当時勉強したのは主にJordan標準形への変形方法だけで、その背後にある数学的な面白さについては理解していなかった。

そこで、本稿ではJordan標準形とは何なのか、その理論的な詳細について考えてみる。

行列は線形写像の映し鏡

始めに、本稿を読み進めるにあたって重要となる考え方について述べておく。

ベクトル空間上の線形写像は行列によって表現することができる。特に、線形変換は線形写像であるから、これも行列によって表現可能である。対角化やJordan標準形を考える上では、この逆を考えることが重要である。すなわち、ある行列に対して、それが何らかの線形変換の表現行列なのだと捉えるのである。

行列というのは線形変換をある基底に対して表現したものに過ぎず、基底の取り方によって姿を変え得る。唯一不変なのはその背後にある線形変換であり、線形変換を通して行列の振る舞いを考えることがJordan標準形の理解へと繋がる。

線形変換を通したベクトル空間の分解

固有空間への分解

 n次元ベクトル空間 Vに対して、ある線形変換 Tが定められており、 m個の固有値 \lambda_1, \lambda_2, \cdots ,  \lambda_mを持つとする。この時、各固有値について、それに対応する固有空間 W_1, W_2, \cdots, W_nが存在する。ここでは線形変換の固有値、及び固有空間について述べているのであって、行列については一切触れていないことに注意されたい。

固有空間の基底は固有ベクトルであるから、固有空間に属する任意のベクトルに Tを作用させると、それらは単に元のベクトルの固有値倍となる。ゆえに、固有空間は T不変な部分空間となる。

全ての固有空間の次元の和が nに等しい場合、以下の式が成立する。

 \displaystyle{
\begin{eqnarray}
&& V = W_1 \oplus W_2 \oplus \cdots \oplus W_n \\
&& T(W_i) \subseteq W_i
\end{eqnarray}
}

広義固有空間への分解

さて、いつもこのようになっていれば良いのだが、残念ながらそうはならないケースが存在する。すなわち、全ての固有空間の次元の和が nを下回るような場合である。このような場合でも、何とかうまく Vを直和分解出来ないだろうか?

実は、これは可能である。 W_iが固有空間であるという条件を緩和して、 T不変性と直和分解されるという性質だけを担保することで、 Vを先ほどと全く同じような直和分解の形に持ち込むことができる。

では、具体的にどのような部分空間に分解してやれば良いだろうか?それを考えるために、まずは固有ベクトル {\bf x} \in W_iについて、以下のような式変形を行う。

 \displaystyle{
\begin{eqnarray}
T{\bf x} &=& \lambda_i {\bf x} \\
(T - \lambda_i I){\bf x} &=& {\bf 0}
\end{eqnarray}
}

これより、固有ベクトルというのは T - \lambda_i Iを作用させると零ベクトルになるようなベクトルであると言える。または、簡潔に W_i = \mathrm{Ker} (T - \lambda_i I)である。

ここから着想を得て、 T - \lambda_i Iだけでなく、 (T - \lambda_i I)^k ( kは任意の自然数) を作用させて零ベクトルになるようなベクトルの集合を考えてみる。これが実は Vの部分空間となっており、さらに T不変性と直和分解されるという性質を満たしている。このようにして得られる空間を広義固有空間と呼ぶ。これを W_i'と表記すると、やはり簡潔に W_i' =  \mathrm{Ker} \{(T - \lambda_i I)^k\}と書ける。

少し確認してみよう。まず T不変性について、以下の式を考える。

 \displaystyle{
\begin{eqnarray}
(T - \lambda_i I)^k (T{\bf x}) &=& (T - \lambda_i I)^{k-1}(T^2 - \lambda_i T) {\bf x} \\
&=& (T - \lambda_i I)^{k-1}T(T - \lambda_i) {\bf x} \\
&=& \cdots \\
&=& T(T - \lambda_i I)^k {\bf x} \\
&=& {\bf 0}
\end{eqnarray}
}

 (T - \lambda_i I)^kを作用させて零ベクトルになったので、 T{\bf x} \in W_i'が言える。

続いて Vが各 W_i'の直和になることについてだが、これは思ったより証明が大変なので、ここではサボって[1]に譲ることにする。

直和分解について一点だけ注意事項を述べておく。ここまで広義固有空間への直和分解について説明したが、実は広義固有空間はさらに直和分解できる場合がある。具体的には、後述するJordan鎖の数だけさらに直和分解可能である。詳細はここでは述べないが、[2]などが参考になるだろう。

広義固有空間の基底

次に、 W_i'の基底について考えてみよう。 W_i'の次元を m_iとする。固有ベクトル {\bf x}_1 \in W_i'について、 (T - \lambda_i I) {\bf x}_2 = {\bf x}_1を満たすベクトル {\bf x}_2を考える。このような {\bf x}_2は存在するかもしれないし、しないかもしれない。もし存在すれば、 (T - \lambda_i I)^2 {\bf x}_2 = {\bf 0}となるため、 {\bf x}_2 \in W_i'が言える。しかも、これは {\bf x}_1スカラー倍でもない。そのため、 \{{\bf x}_1, {\bf x}_2\}は一次独立である。

帰納的に、 (T - \lambda_i I) {\bf x}_{j+1} = {\bf x}_jとなるベクトル {\bf x}_{j+1}を考える。すると、 (T - \lambda_i I)^{j+1} {\bf x}_{j+1} = {\bf 0}となるため {\bf x}_{j+1} \in W_i'が言える。しかも、 \{{\bf x}_1, {\bf x}_2\, \cdots , {\bf x}_{j+1}\}は一次独立になる。

基底は最大で m_i個しか取れないので、この操作はどこかで頭打ちとなる。このようにして得られるベクトルの列をJordan鎖という。Jordan鎖は1つ以上存在し、全てのJordan鎖を構成する全ベクトルを寄せ集めると、これは W_i'の基底となる。これを広義固有ベクトルと呼ぶ。

広義固有空間の構造を完全に明らかにするためには、Jordan鎖は長さいくつのものが何本存在するのかを知る必要がある。これは \mathrm{Ker} \{(T - \lambda_i I)\}, \mathrm{Ker} \{(T - \lambda_i I)^2\}, \cdotsの次元を順に計算し、それらがどのように増えていくかを調べれば分かる。が、込み入った話になるので詳細は[3]を参照されたい。

広義固有ベクトルを基底とした場合の表現行列

以上、線形変換の性質についていろいろと述べたが、ここからいよいよJordan標準形の話に入っていく。

いきなり天下り的だが、線形変換 Tについて、全ての固有値に対する広義固有ベクトルを全て集め、それらを基底とした場合の表現行列について考えてみよう。同じ広義固有空間から抽出した広義固有ベクトルは隣り合うように並べて、その中でさらに同じJordan鎖に属するベクトルも順に並べて添字をつけたものを \{{\bf x}_1, {\bf x}_2\, \cdots , {\bf x}_{n}\}とする。ここで、以下のような写像 \phi : V \to \mathbb{R}^nを定義する。

 \displaystyle{
\phi({\bf x}_i) = {\bf e}_i\ (i=1, 2, \cdots , n)
}

ここで、 {\bf e}_i \mathbb{R}^nの標準基底である。これによって Tの表現行列 Aが定まる。 Aによって実現される \mathbb{R}^n上の線形変換を T_Aとすると、以下のような可換図式が得られる。

 \displaystyle{
\xymatrix{
    V \ar[r]^{T} \ar[d]^{\phi} & V \ar[d]^{\phi} \\
    \mathbb{R} \ar[r]^{T_A} & \mathbb{R}
}
}

これを成立させるためには、 Aはどのような行列であれば良いだろうか?これを一般的な状況で説明するのは非常に煩雑なので、ここでは以下のような具体的な設定の元で議論を進める。

  •  n = 6
  • 固有値 \lambda_1, \lambda_2の2つであり、それぞれの代数的重複度は5, 1である。
  •  \lambda_1について、以下の2つのJordan鎖がある。
    •  {\bf x}_1, {\bf x}_2, {\bf x}_3
      •  (T - \lambda_1 I){\bf x}_1 = 0
      •  (T - \lambda_1 I){\bf x}_2 = {\bf x}_1
      •  (T - \lambda_1 I){\bf x}_3 = {\bf x}_2
    •  {\bf x}_4, {\bf x}_5
      •  (T - \lambda_1 I) {\bf x}_4 = 0
      •  (T - \lambda_1 I){\bf x}_5 = {\bf x}_4
  •  \lambda_2に対応する固有値ベクトルは {\bf x}_6とする。
    •  (T - \lambda_2 I) {\bf x}_6 = 0

まず、上で示した可換図式により T_A = \phi \circ T \circ \phi^{-1}となる。また、具体的な設定の中で示した式を変形すると以下のようになる。

 \displaystyle{
\begin{eqnarray}
&&T{\bf x}_1 = \lambda_1 {\bf x}_1 \\
&&T{\bf x}_2 = {\bf x}_1 + \lambda_1 {\bf x}_2 \\
&&T{\bf x}_3 = {\bf x}_2 + \lambda_1 {\bf x}_3 \\
&&T{\bf x}_4 = \lambda_1 {\bf x}_4 \\
&&T{\bf x}_5 = {\bf x}_4 + \lambda_1 {\bf x}_5 \\
&&T{\bf x}_6 = \lambda_2 {\bf x}_6
\end{eqnarray}
}

さらに、これらに \phiを作用させると以下のようになる。

 \displaystyle{
\begin{eqnarray}
&&(\phi \circ T){\bf x}_1 = \lambda_1 {\bf e}_1 \\
&&(\phi \circ T){\bf x}_2 = {\bf e}_1 + \lambda_1 {\bf e}_2 \\
&&(\phi \circ T){\bf x}_3 = {\bf e}_2 + \lambda_1 {\bf e}_3 \\
&&(\phi \circ T){\bf x}_4 = \lambda_1 {\bf e}_4 \\
&&(\phi \circ T){\bf x}_5 = {\bf e}_4 + \lambda_1 {\bf e}_5 \\
&&(\phi \circ T){\bf x}_6 = \lambda_2 {\bf e}_6
\end{eqnarray}
}

上式に {\bf x}_i = \phi^{-1}({\bf e}_i)を代入することで、 T_Aは以下のような変換であることが分かる。

 \displaystyle{
\begin{eqnarray}
&&T_A{\bf e}_1 = \lambda_1 {\bf e}_1 \\
&&T_A{\bf e}_2 = {\bf e}_1 + \lambda_1 {\bf e}_2 \\
&&T_A{\bf e}_3 = {\bf e}_2 + \lambda_1 {\bf e}_3 \\
&&T_A{\bf e}_4 = \lambda_1 {\bf e}_4 \\
&&T_A{\bf e}_5 = {\bf e}_4 + \lambda_1 {\bf e}_5 \\
&&T_A{\bf e}_6 = \lambda_2 {\bf e}_6
\end{eqnarray}
}

よって、 Aは以下のような行列であることが分かる。

 \displaystyle{
A = 
\left(
\begin{array}{cccccc}
\lambda_1 & 1 & 0 & 0 & 0 & 0  \\
0 & \lambda_1 & 1 & 0 & 0 & 0  \\
0 & 0 & \lambda_1 & 0 & 0 & 0  \\
0 & 0 & 0 & \lambda_1 & 1 & 0  \\
0 & 0 & 0 & 0 & \lambda_1 & 0  \\
0 & 0 & 0 & 0 & 0 & \lambda_2
\end{array}
\right)
}

これはJordan標準形そのものである。結局、線形変換の基底として広義固有ベクトルを選んだ場合の表現行列こそが、Jordan標準形の正体なのである。

広義固有ベクトルへの座標変換

最後に、適当な行列をJordan標準形に変形するとはどういうことなのか、その意味を考えてみよう。

Jordan標準形になっていない n次正方行列 Bを考える。 Bをあるベクトル空間 V上の線形変換 Tの表現行列であると考えると、これは基底として広義固有ベクトル以外のものを選んだ場合であると解釈できる。これをJordan標準形に変形する事は、基底を広義固有ベクトルに取り替えることを意味する。

 Bが表現行列となるような Vの基底を {\bf y}_1, {\bf y}_2, \cdots, {\bf y}_nとし、写像 \psi: V \to \mathbb{R}^nを以下のように定める。

 \displaystyle{
\psi({\bf y}_i) = {\bf e}_i\ (i=1, 2, \cdots , n)
}

すると、以下の可換図式が得られる。

 \displaystyle{
\xymatrix{
    \mathbb{R} \ar[r]^{T_B} & \mathbb{R} \\
    V \ar[r]^{T} \ar[d]^{\phi} \ar[u]_{\psi} & V \ar[d]^{\phi} \ar[u]_{\psi} \\
    \mathbb{R} \ar[r]^{T_A} & \mathbb{R}
}
}

これより、広義固有ベクトルを基底とした場合の変換 T_Aは以下のようにして得られる。

 \displaystyle{
\begin{eqnarray}
T_A &=& \phi \circ \psi^{-1} \circ T_B \circ \psi \circ \phi^{-1} \\
&=& (\psi \circ \phi^{-1})^{-1} \circ T_B \circ (\psi \circ \phi^{-1})
\end{eqnarray}
}

 \psi \circ \phi^{-1}は基底の取り替えを表す写像である。この写像の表現行列を Pとすると、よく見慣れたJordan標準形への変換式 A = P^{-1}BPが得られる。

ちなみに、 Pを具体的に求めると、これは広義固有ベクトルを列ベクトルとして順に並べた行列になっている。そのようになる理由は説明すると長くなるので、私の体力の都合により割愛する。

まとめ

本稿ではJordan標準形の理論的側面について述べた。結論として、Jordan標準形とは広義固有ベクトルを基底としたときの線形変換の表現行列であることが分かった。また、Jordan標準形への変換とは、基底を広義固有ベクトルに取り替える操作であることが分かった。

本稿ではJordan標準形の応用面について触れることが出来なかった。これについてはまた機会があれば調べてみたい。

固有ベクトルの観点から線形変換を掌握する

線形代数において、固有値固有ベクトルの話題は花形である。これらは理論的に美しいだけでなく、応用上様々な場面に登場し、重宝されている。中でも固有値が力を発揮するのは、やはり行列の対角化を行う時であろう。すなわち、n次正方行列 Aがある条件を満たすと、ある変換行列 Pが存在して P^{-1} A Pが対角行列になるのである。

では、対角化が可能となるための「ある条件」とは一体なんだろう?それは、固有値の (固有方程式の解としての) 重複度と、その固有値に対応する固有空間の次元が全ての固有値に対して等しいことである。これは対角化可能であるための必要十分条件となっている。

そう、固有値というのは実に厄介で、例えば固有方程式を解いた結果、 \lambdaが重解として得られたとしても、 \lambdaに対応する固有空間は必ずしも2次元になるとは限らないのである。

以上の事実から、行列は固有値の数と重複度、および固有ベクトルの次元によって、いくつかの種類に分類出来そうな気がしてくる。そのように分類された行列のクラスは、それぞれ何か特徴的な性質を持つのだろうか?

本稿ではこのぼんやりとした疑問の答えに迫るべく、前半は対角化について調べ、後半は具体例を頼りに固有値・固有空間が反映する行列の性質について調べてみようと思う。

対角化とは何か?

まずは対角化について考えてみよう。対角化と言えば P^{-1} A Pという変換の仕方が特徴的である。これの意味するところを考えることによって、対角化について理解を深めていこう。

線形変換の表現

行列は線形変換と密接な関係にある。すなわち、任意の線形変換は行列によって表現することができる。 P^{-1} A Pの意味を理解する上では、Aをある線形変換の表現行列と捉えるのがよい。

まず、線形変換の定義を[1]より引用する。

線形変換
 K上の線形空間 Vから K上の線形空間 V'への写像 Tが二条件
 \displaystyle{
\begin{eqnarray}
&& T({\bf x} + {\bf y}) = T({\bf x}) + T({\bf y}) \\
&& T(a{\bf x}) = aT({\bf x})
\end{eqnarray}
}
を充すとき,  T Vから V'への線形写像と言う.
・・・ (中略) ・・・
とくに,  Vから V自身への線形写像を,  V線形変換と言う.

 Vの次元をnとする。 Vの1つの基底を {\bf e}_1,  {\bf e}_2, \cdots ,  {\bf e}_nとすると、 Vの任意の元は以下のように基底の線型結合で表すことができる。

 \displaystyle{
{\bf v} = \sum_{i=1}^n x_i {\bf e}_i
}

ここで、 {\bf v} Vの線形変換 Tで写したものはやはり Vの元となる。よって以下のように同じ基底で表すことができる。

 \displaystyle{
T({\bf v}) = \sum_{i=1}^n x_i' {\bf e}_i
}

すると、変換前後の係数の間には以下のような関係が成立する。

 \displaystyle{
 \left( \begin{array}{c}
      x_1' \\
      x_2' \\
      \vdots \\
      x_n'
    \end{array} \right)_e
  = \left( \begin{array}{cccc}
      a_{11} & a_{12} & \ldots & a_{1n} \\
      a_{21} & a_{22} & \ldots & a_{2n} \\
      \vdots & \vdots & \ddots & \vdots \\
      a_{n1} & a_{n2} & \ldots & a_{nn}
    \end{array} \right)
 \left( \begin{array}{c}
      x_1 \\
      x_2 \\
      \vdots \\
      x_n
    \end{array} \right)_e
}

上の行列 A = (a_{ij})が線形変換 Tの表現行列である。両辺のベクトルの右下にeと書いたのは、これが基底 {\bf e}_iによるベクトルであることを示している。表現行列を簡潔に Aと表した式は以下のようになる。

 \displaystyle{
 \left( \begin{array}{c}
      x_1' \\
      x_2' \\
      \vdots \\
      x_n'
    \end{array} \right)_e
  = A
 \left( \begin{array}{c}
      x_1 \\
      x_2 \\
      \vdots \\
      x_n
    \end{array} \right)_e
}

基底の取り替え

次に、線形変換 T Vの別の基底で表現するとどうなるか考えてみよう。そのような基底を {\bf f}_1,  {\bf f}_2, \cdots ,  {\bf f}_nとする。すると、先ほどと同じように変換前後のベクトルは以下のように表される。

 \displaystyle{
{\bf v} = \sum_{i=1}^n y_i {\bf f}_i
}

 \displaystyle{
T({\bf v}) = \sum_{i=1}^n y_i' {\bf f}_i
}

基底の取り替えにおいて重要な考え方は、元の基底に対する表現行列 Aをできる限り活かすことである。そのために、 Aではなく入力として与えるベクトルの方を変換することを考える。これを実現するためには、以下のような行列 Pを利用する。

 \displaystyle{
 \left( \begin{array}{c}
      x_1 \\
      x_2 \\
      \vdots \\
      x_n
    \end{array} \right)_e
=
P \left( \begin{array}{c}
      y_1 \\
      y_2 \\
      \vdots \\
      y_n
    \end{array} \right)_f
}

このような行列 Pを基底の取り替え行列と呼ぶ。これより、 Tは以下のように表現できる。

 \displaystyle{
P \left( \begin{array}{c}
      y_1' \\
      y_2' \\
      \vdots \\
      y_n'
    \end{array} \right)_f
  = A
P \left( \begin{array}{c}
      y_1 \\
      y_2 \\
      \vdots \\
      y_n
    \end{array} \right)_f
}

異世界を旅して帰還する

基底の変換後の式の両辺に P^{-1}をかけると以下のようになる。

 \displaystyle{
 \left( \begin{array}{c}
      y_1' \\
      y_2' \\
      \vdots \\
      y_n'
    \end{array} \right)_f
  = P^{-1}A P
 \left( \begin{array}{c}
      y_1 \\
      y_2 \\
      \vdots \\
      y_n
    \end{array} \right)_f
}

 P^{-1}A Pが基底 {\bf f}_iを用いた場合の Tの表現行列となる。

普通はここで話はおしまいである。基底を変換した後の表現行列が求められたのだから、めでたしめでたしというわけだ。しかし、 P^{-1}A Pという形の変換は線形代数に限らずいろんな数学の分野でよく見かける形であり、もう一段抽象的なレベルでの意味がある。これを少し掘り下げて考えてみよう。

まず、行列 Aは基底 {\bf e}_iの世界で表現されるベクトルに対して適用可能な行列である。これを別の世界のベクトル、すなわち、基底 {\bf f}_iで表現されるベクトルに適用したいと思ったら、一度 {\bf f}_iの世界から {\bf e}_iの世界に移動する必要がある。すると、移動後のベクトルには Aを適用することができる。最後に、 Aを適用して得られたベクトルを元の世界に戻してやれば、結局 {\bf f}_iの世界のベクトルに Aを適用出来たことになる。以下に図を示す。

f:id:peng225:20180917213651p:plain

このように P^{-1}, Pで変換を挟み込むことは、その変換が適用可能な世界へと移動し、変換が終わったら元の世界に戻してやるような効果がある。これが線形変換の表現行列を別の基底で表したベクトルに適用できるからくりである。

対角化は座標変換

ここまで座標変換の話ばかりしてきた訳だが、座標変換と対角化の関係について述べておく。ずばり対角化とは、線形変換の基底を固有ベクトルから成る基底に変換する操作のことである。そして、対角化可能であるとは、基底となるような固有ベクトルの組が存在することを意味している。

逆に、もし固有値の重複度より固有ベクトルの次元が小さいものが存在すると、対象となる線形空間の次元に対して固有ベクトルの本数が足りなくなってしまい、基底を作れなくなる。この場合は対角化不可能となる。

基底の取り替え行列 Pとしては、 Aのn個の列固有ベクトル {\bf p}_1, {\bf p}_2, \cdots, {\bf p}_nを並べた行列を取れば対角化できる。計算の過程を以下に示す。

 \displaystyle{
\begin{eqnarray}
P^{-1} A 
 \left( \begin{array}{cccc}
      {\bf p}_1 & {\bf p}_2 & \cdots & {\bf p}_n
    \end{array} \right)
 &=& P^{-1} \left( \begin{array}{cccc}
      \lambda_1 {\bf p}_1 & \lambda_2 {\bf p}_2 & \cdots & \lambda_n {\bf p}_n
    \end{array} \right) \\
 &=& \left( \begin{array}{cccc}
      \lambda_1 & 0 & \ldots & 0 \\
      0 & \lambda_2 & \ldots & 0 \\
      \vdots & \vdots & \ddots & \vdots \\
      0 & 0 & \ldots & \lambda_n
    \end{array} \right)
\end{eqnarray}
}

固有値固有ベクトルは何に紐づくものか?

ここで注意点を1つ。世の中ではよく「行列の固有値」とか「行列の固有ベクトル」という言い方をする。これは別に間違っていない。実は、線形変換に対しても固有値固有ベクトルという概念が存在する。

まず固有値についてだが、行列の固有値 P^{-1}, Pのサンドイッチ変換に対して不変となる。そのため、線形変換の固有値とは「任意の基底に対する表現行列の固有値である」と定義しておけば、それは一意に定まる。

多様体微分幾何学なんかが代表的であるが、座標系に依存しない量や概念というのは、その数学的対象の本質的な性質を表していると考えられるため、一般にとても重要である。

一方、固有ベクトルは座標系の取り方によって変わる。このように、固有値固有ベクトルでは不変となる範囲が異なるので注意が必要である。

固有値・固有空間に着目した線形変換の分類

ここまでで固有値固有ベクトルとかなりお友達になれたはずなので、本題に入ろう。線形変換が固有値の数や重複度、また固有空間の次元とどのような関係にあるのかを調べてみよう。簡単のために、話を \mathbb{R}上2次の正方行列に絞る。

 \mathbb{R}上2次の正方行列は以下のように分類できる。

  1. 異なる2つの実固有値を持つ場合
  2. 1つの実固有値を持ち、固有空間の次元が2の場合
  3. 1つの実固有値を持ち、固有空間の次元が1の場合
  4. 固有値を持たない場合

以下では分類された各項目の線形変換がどのような振る舞いをするのか、例を用いながら調べてみる。

異なる2つの実固有値を持つ場合

例として以下の行列を考える。

 \displaystyle{
A = \left( \begin{array}{cc}
     1  & -1/4 \\
  -1/2 & 5/4 \\
    \end{array} \right)
}

 A固有値固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[1, -1/4], [-1/2, 5/4]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(3/2, [   # 固有値1つ目
(1, -2)   # 固有ベクトル
], 1)     # 重複度
(3/4, [   # 固有値2つ目
(1, 1)    # 固有ベクトル
], 1)     # 重複度

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002223p:plain

ここで、黒の矢印は基底、グレーの点線矢印は基底に Aを作用させて得られるベクトルである。赤の矢印は固有ベクトル、オレンジの点線矢印は固有ベクトル Aを作用させて得られるベクトルである。また、水色の矢印は、矢印の根元の位置ベクトルが、 Aによって矢印の先端に写ることを意味している。

これを見ると、Aの作用によって空間が2つの固有ベクトルの方向に伸び縮みしているのが分かる。これより、このタイプの線形変換は固有ベクトル方向の拡大・縮小を表すと考えられる。ただし、伸びるのか縮むのか、またその程度がどれくらいであるかは固有値によって決まる。

1つの実固有値を持ち、固有空間の次元が2の場合

例として以下の行列を考える。

 \displaystyle{
A = \left( \begin{array}{cc}
     2  & 0 \\
     0 & 2 \\
    \end{array} \right)
}

 A固有値固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[2, 0], [0, 2]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(2, [     # 固有値
(1, 0),   # 固有ベクトル1つ目
(0, 1)    # 固有ベクトル2つ目
], 2)     # 重複度

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002458p:plain

これを見ると、Aの作用によって空間が原点を中心に拡大されているのが分かる。これより、このタイプのの線形変換は拡大・縮小を行うものと考えられる。

一点補足だが、このケースの線形変換の表現行列は単位行列の定数倍しかあり得ない。 Aのただ1つの固有値 \lambdaとすると、 Aを対角化した行列は \lambda Eとなる。対角化のための変換行列を Pとすると、結局以下のようになる。

 \displaystyle{
\begin{eqnarray}
P^{-1} A P &=& \lambda E \\
A &=& P (\lambda E) P^{-1} \\
A &=& \lambda E
\end{eqnarray}
}

よって A単位行列固有値倍となる。

1つの実固有値を持ち、固有空間の次元が1の場合

例として以下の行列を考える。

 \displaystyle{
A = \left( \begin{array}{cc}
     11/6 & -1/3 \\
       1/3  & 7/6
    \end{array} \right)
}

 A固有値固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[11/6, -1/3], [1/3, 7/6]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(3/2, [   # 固有値
(1, 1)    # 固有ベクトル
], 2)     # 重複度

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002520p:plain

これを見ると、原点を通る固有ベクトルによって引かれる直線より下側では右向きの、上側では左向きの流れがあるように見える。これより、このタイプの線形変換は固有ベクトルと平行な方向に歪みを与える働きがあると考えられる。

固有値を持たない場合

例として以下の行列を考える。

 \displaystyle{
A = \left( \begin{array}{cc}
     1   & -1/2 \\
    1/2 &   1
    \end{array} \right)
}

 A固有値固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[1, -1/2], [1/2, 1]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(1 - 0.50000000000000000?*I, [(1, 1*I)], 1)    # 固有値が虚数
(1 + 0.50000000000000000?*I, [(1, -1*I)], 1)

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002536p:plain

これを見ると、螺旋のような渦巻きが現れている。これより、このタイプの線形変換は回転と拡大・縮小の組み合わせになっていると考えられる。

まとめ

以上、線形変換についていろいろと調べてみた。その結果、線形変換は固有値固有ベクトルによっていくつかの種類に分類できそうだということが分かった。本稿では2次正方行列で表現される線形変換しか詳細に調べなかったが、3次以上についても類似の考察が可能だろうと推測される。3次までなら頑張って可視化できないこともないかもしれないので、機会があればトライしてみたい。

参考文系

[1]

線型代数入門 (基礎数学1)

線型代数入門 (基礎数学1)

区間推定の理論と実践

統計学における推定の考え方は仮説検定と双璧を成す重要な理論であり、初学者がまず目指すべき最高到達点である。推定は強力な手法であるが故に、恐らく大抵の教科書では記載があるだろうし、ググればいくらでもWebページがヒットする。

そんな状況下ではあるが、推定の中でも特に区間推定は個人的に興味をそそられる内容であり、ぜひともブログを書き起こしてみたくなった。そこで、本稿ではできるだけ二番煎じにならないよう、私が得た区間推定に対する理解について説明してみたいと思う。

本稿では正規母集団から抽出した標本データを用いた母平均の区間推定のみをターゲットとする。また、適当なサンプルデータを用いて実際に区間推定を行うことで、区間推定の理論に対して確かな実感を得てみたいと思う。

なお、本稿を執筆するに辺り、全体的に[1][2]を参考にした。

母平均の区間推定

母分散が既知の場合

正規母集団 N(\mu,\ \sigma^2)に従う独立な確率変数 X_1,\ X_2,\ \cdots .\ X_nの標本平均 \overline{X}を考える。このとき、 \overline{X}自体も確率変数となっており、これが従う確率分布は N(\mu,\ \sigma^2/n)となる。また、これを標準化した確率変数 Z = \frac{\overline{X} - \mu}{\sigma/\sqrt{n}}は、当然だが標準正規分布 N(0,\ 1)に従う。

今、母分散 \sigma^2は分かっているが、母平均 \muの値が不明な状況を考える。もちろん、実際にはこんなことはまずあり得ないのだが、このケースを考えることは後々の議論において意味がある。

この時、 \mu = E(\overline{X})という関係が成り立つので、標本平均をそのまま母分散だと考えてもあながち間違いではない。しかし、当然だが \overline{X}の値にはぶれがあり、どんなサンプルを抽出するかによって値が変わってしまう。しかも、極端なケースでは \overline{X}はとんでもなく大きくなったり小さくなったりすることがあり得る。そのため、 \overline{X}の値だけをみても、母平均について何も確実なことは言えない。

では、何か確率的な議論はできないだろうか?例えば「95%の確率で \muはAからBの間に入ります」というようなことが言えないだろうか?

結論から言うとこれは可能である。標準正規分布において、それより大きい値を取る確率が \alpha\ (0 < \alpha < 1)となるような点 ( 100 \alphaパーセント点) を Z_\alphaと表すと、先ほど導入した確率変数Zは確率 1 - \alphaで以下の不等式を満たす。

 \displaystyle{
 -Z_{\alpha/2} < Z < Z_{\alpha/2}
}

2点注意事項がある。まず、これは正規分布が左右対称だから言えることである。さもなければ、  -Z_{\alpha/2}の部分はそれより小さい値を取る確率が \alpha/2となるような点を別途求めて、それと置き換えてやる必要がある。

次に、 \alphaを2で割っている理由だが、これは上側の \alpha/2と下側の \alpha/2を両方排除し、残りが 1 - \alphaとなるようにするためである。

この不等式を変形すれば、以下のように \muが確率 1 - \alphaで取る値の範囲が分かる。

 \displaystyle{
\overline{X} - \frac{Z_{\alpha/2} \sigma}{\sqrt{n}} <  \mu < \overline{X} + \frac{Z_{\alpha/2}\sigma}{\sqrt{n}}
}

母分散が未知の場合

先ほども述べた通り、母平均は未知なのに母分散は分かっているということは普通はあり得ない。実際に得られるのは不偏分散 s^2だけというケースがほとんどである。不偏分散の定義式を以下に示す。

 \displaystyle{
s^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \overline{X})^2
}

不偏分散に対しては \sigma^2 = E(s^2)という関係式が成り立つので、これを母分散の代わりに用いることはそれほど不自然ではない。しかし、やはりサンプル数が少ない場合には母分散とかけ離れた値を取ることがあり得る。そのため、先ほど考えた確率変数Zの式に対して単に \sigma = sを代入すると、これはもはや正規分布には従わない。

では、これはどのような分布になるのだろうか?それを調べるために、Zの式に \sigma = sを代入した式を以下のように変形してみよう。

 \displaystyle{
\begin{eqnarray}
t &=& \frac{\overline{X} - \mu}{\frac{s}{\sqrt{n}}} \\
&=& \frac{\overline{X} - \mu}{\frac{\sigma}{\sqrt{n}} \cdot \frac{s}{\sigma}} \\
&=& \frac{Z}{\frac{\sqrt{\frac{1}{n-1} \sum_{i=1}^n (X_i - \overline{X})^2}}{\sigma}} \\
&=& \frac{Z}{\sqrt{\frac{1}{n-1} \sum_{i=1}^n \left(\frac{X_i - \overline{X}}{\sigma}\right)^2}}
\end{eqnarray}
}

ここで、分母に \sum_{i=1}^n \left(\frac{X_i - \overline{X}}{\sigma}\right)^2という式が現れた。この式自体も確率変数であり、これは標準正規分布に従う確率変数の二乗和を表す。このような確率変数が従う確率分布は自由度n-1の \chi^2分布と呼ばれる。

まとめると、標準正規分布に従う確率変数Zと自由度n-1の \chi^2分布に従う確率変数Wによって、元々の確率変数tは以下のように表される。

 \displaystyle{
t = \frac{Z}{\sqrt{\frac{W}{n-1}}}
}

このような確率変数tが従う確率分布は自由度n-1のt分布と呼ばれる。結論として、母分散未知のケースにおける母平均の推定にはt分布を使用する必要があることが分かった。

自由度n-1のt分布を t(n-1)と書くと、母分散既知のケースと同様の計算を行うことで、 \muは確率 1 - \alphaで以下の不等式を満たすことが分かる。

 \displaystyle{
\overline{X} - \frac{t_{\alpha/2}(n-1) s}{\sqrt{n}} <  \mu < \overline{X} + \frac{t_{\alpha/2}(n-1) s}{\sqrt{n}}
}

ただし、 t_{\alpha/2}(n-1)は自由度n-1のt分布の 100 \alpha/2パーセント点を表す。ここでもt分布が左右対称であることを利用している。

サンプル数に関する注意

t分布は自由度が大きくなると正規分布に収束していく。そのため、t分布を用いる必要があるのはサンプル数が少ない場合に限る。実用上はサンプル数が30以上あれば正規分布を使用しても差し支えないようである[1]。

区間推定の実践

ここまでに述べたことを実践してみよう。統計処理をするにはpythonかRが簡単そうであるが、Android上で環境を構築できたのでpythonを使うことにする[3]。また、データセットとしてIrisは見飽きたので、ここではAbarone (アワビ) を使用している[4]。

Abaroneデータセットには、1列目に性別、3列目に直径の情報が入っている。今回はオスのアワビの直径データについて、母平均の区間推定を行ってみた。サンプル数は敢えて少なめに10とし、t分布を使った場合と正規分布を使った場合で比較を行った。もちろん、正しいのはt分布を使う方である。

以下に今回作成したソースコードを示す。作成の際には[5]を参考にさせて頂いた。日本語コメントは本稿のために後から付したものであり、このままでは実行できないので注意されたい。

import pandas as pd
from scipy import stats
import math

# 全部の列に名前を付けるのは諦めた...
abalone = pd.read_csv('file:///storage/emulated/0/Download/abalone.data', 
                 header=None, names=['sex', 'length', 'diameter', 'height', 'a', 'b', 'c', 'd', 'e'])

# オスの情報だけ抽出
male = abalone[abalone['sex'] == 'M']

n = 10
# 自由度n-1のt分布の95%信頼区間を取得
tleft, tright = stats.t.interval(alpha=0.95, df=n-1)
# 正規分布の95%信頼区間を取得
nleft, nright = stats.norm.interval(alpha=0.95)

mean = male.iloc[:n]['diameter'].mean()
# デフォルトで不偏分散の平方根が計算される
stddev = male.iloc[:n]['diameter'].std()

# 本文中の式と若干符号が異なるので注意
tleft = mean + tleft*stddev/math.sqrt(n)
tright = mean + tright*stddev/math.sqrt(n)

nleft = mean + nleft*stddev/math.sqrt(n)
nright = mean + nright*stddev/math.sqrt(n)

print("mean = {:.3f}, stddev = {:.3f}".format(mean, stddev))
print("t-dist: left = {:.3f}, right = {:.3f}".format(tleft, tright))
print("n-dist: left = {:.3f}, right = {:.3f}".format(nleft, nright))

実行結果を以下に示す。

mean = 0.339, stddev = 0.046
t-dist: left = 0.306, right = 0.372
n-dist: left = 0.310, right = 0.368

実行結果から分かるように、正規分布を使った場合は区間が狭めに出力されてしまう。そのため、正規分布での結果を元に考察を行うと判断を誤る恐れがある。

最後に特大の注意事項がある。今回はオスのアワビの直径データが正規分布に従うと仮定して分析を行ったが、本来であればこの仮定の妥当性は良く確かめる必要がある。誤った仮定を元に一生懸命計算したとしても、その結果には何の意味もないどころか、誤った帰結へと導かれる危険性がある。

データが正規分布に良く当てはまっているかどうかを判断する方法について私には知見がない。軽くググってみたらいろいろと情報がありそうだったので、興味のある方は調べてみると面白いかもしれない。

その他の話題

本稿では正規母集団から抽出したサンプルから母平均の区間推定を行った。この他にも、例えば母分散の区間推定を行う方法などは有名である[1]。また、実用上のタスクとしては正規母集団以外の母集団を考えることもあるだろう。しかし、どんな場合であれ、基本的な区間推定の考え方は本稿に述べたものと大きく異なることはないであろう。

まとめ

以上、区間推定に関する理論と実践について述べた。統計学は応用先が非常にたくさんある分野なので、理解を深めておいて損はない。また機会があれば、多変量解析などの進んだ内容も勉強し直したい。