元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

テンソルを使いこなしたい3(行列の行列微分)

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

の続き。
本記事ではテンソル微分演算についてまとめる。行列の行列微分がややこしかったので、縮約記法で計算すればスッキリできるのではと思ったのが事の発端。ようやくここまで来た。

準備

スカラーのベクトル微分、ベクトルのスカラー微分、ベクトルのベクトル微分スカラーの行列微分を以下のように定義する。また、それぞれの最右辺は縮約記法での表現を表す。 {} $$ \begin{align} \frac{\mathrm{d}f}{\mathrm{d}\boldsymbol{x}} &= \left( \begin{array}{c} \frac{\mathrm{d}f}{\mathrm{d}x_{1}} \\ \frac{\mathrm{d}f}{\mathrm{d}x_{2}} \\ \vdots \\ \frac{\mathrm{d}f}{\mathrm{d}x_{n}} \end{array} \right) = \partial{}_{i} f \\ \frac{\mathrm{d}\boldsymbol{f}}{\mathrm{d}x} &=\left(\begin{array}{c} \frac{\mathrm{d}f_{1}}{\mathrm{d}x} \\ \frac{\mathrm{d}f_{2}}{\mathrm{d}x} \\ \vdots \\ \frac{\mathrm{d}f_{n}}{\mathrm{d}x} \end{array} \right) = \partial{} f_{i} \\ \frac{\mathrm{d}\boldsymbol{f}}{\mathrm{d}\boldsymbol{x}} &= \left(\begin{array}{cccc} \frac{\mathrm{d}f_{1}}{\mathrm{d}x_{1}} & \frac{\mathrm{d}f_{1}}{\mathrm{d}x_{2}} & \ldots & \frac{\mathrm{d}f_{1}}{\mathrm{d}x_{n}} \\ \frac{\mathrm{d}f_{2}}{\mathrm{d}x_{1}} & \frac{\mathrm{d}f_{2}}{\mathrm{d}x_{2}} & \ldots & \frac{\mathrm{d}f_{2}}{\mathrm{d}x_{n}}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\mathrm{d}f_{m}}{\mathrm{d}x_{1}} & \frac{\mathrm{d}f_{m}}{\mathrm{d}x_{2}} & \ldots & \frac{\mathrm{d}f_{m}}{\mathrm{d}x_{n}} \end{array} \right) = \partial{}_{j} f_{i} \label{vec_by_vec} \\ \frac{\mathrm{d}f}{\mathrm{d}\boldsymbol{X}} &= \left(\begin{array}{cccc} \frac{\mathrm{d}f}{\mathrm{d}x_{11}} & \frac{\mathrm{d}f}{\mathrm{d}x_{12}} & \ldots & \frac{\mathrm{d}f}{\mathrm{d}x_{1n}} \\ \frac{\mathrm{d}f}{\mathrm{d}x_{21}} & \frac{\mathrm{d}f}{\mathrm{d}x_{22}} & \ldots & \frac{\mathrm{d}f}{\mathrm{d}x_{2n}}\\ \vdots & \vdots & \ddots & \vdots \\ \frac{\mathrm{d}f}{\mathrm{d}x_{m1}} & \frac{\mathrm{d}f}{\mathrm{d}x_{m2}} & \ldots & \frac{\mathrm{d}f}{\mathrm{d}x_{mn}} \end{array} \right) = \partial{}_{ij} f \label{scalar_by_matrix} \end{align} $$

また、\eqref{vec_by_vec}より  {} $$ \begin{align} \partial{}_{j} x_{i}= \delta_{ij} \end{align} $$ である(この後使う)。

なお、\eqref{vec_by_vec}は転置した形として定義されていることもある。行列表現のまま計算を進めようとすると、行列積が可換でないことに注意する必要がある。縮約記法で計算すれば気にしなくて良くなる。

さて、ベクトルの行列微分、行列のベクトル微分、行列の行列微分はどうなるのだろう。行列の行列微分がわかれば他もノリでわかると思うので、上記と同様に縮約記法を使って行列の行列微分の定義から、ニューラルネットワークの逆誤差伝播の式を導出するところまでやってみる。

行列の行列微分

行列 \displaystyle Xの自身による微分は縮約記法を使って以下のように定義される。

 {} $$ \begin{align} \frac{\partial X_{kl}}{\partial X_{ij}} = \delta_{ik}\delta_{li} \end{align} $$

要は、 \displaystyle X_{kl} \displaystyle X_{kl}微分したときは \displaystyle 1それ以外は \displaystyle 0と言っているだけである。「じゃあ全成分 \displaystyle 1になるじゃん」とか一瞬思うかもしれない(自分は思った)がそれは大間違い。 \displaystyle k,lそれぞれについて \displaystyle i,jを走査して計算するので4階テンソルになっている(添字も4つあるし)。

逆誤差伝播

以下のような単純なニューラルネットワークを考える。 Lは損失関数でスカラーである。形式的に書いてるだけなので実態は何でも良い。

f:id:hotoke-X:20190210064732p:plain

これを式で書けば

$$ \begin{align} Y &= WX + B \\ X &= \left(\begin{array}{c} x_{1} \\ \vdots \\ x_{n} \end{array} \right), W = \left(\begin{array}{ccc} w_{11} & \ldots & w_{1n} \\ \vdots & \ddots & \vdots \\ w_{m1} & \ldots & w_{mn} \end{array} \right), B = \left(\begin{array}{c} b_{1} \\ \vdots \\ b_{m} \end{array} \right) \end{align} $$

である。計算グラフにして逆誤差伝播を赤で書き入れると以下のようになる。

f:id:hotoke-X:20190210082429p:plain

ナニコレ。特に \displaystyle \left( \frac{\partial Y}{\partial X}\right)^{\mathrm{T}} \frac{\partial L}{\partial Y}。そりゃ \displaystyle Xの次元に一致しなきゃいけないから逆算すればわかるけど気持ち悪い。一方で、 \displaystyle \frac{\partial L}{\partial Y} \frac{\partial Y}{\partial W}も連鎖律としてはスッキリしているが、ベクトルの行列微分がある。しかも \displaystyle X^{\mathrm{T}}についてもなんで転置されてるんだかよくわからない。

逆誤差伝播(縮約記法でスッキリ ver.)

では縮約記法で \displaystyle \frac{\partial L}{\partial X}, \frac{\partial L}{\partial W}を計算してみよう。

 {} $$ \begin{align} Y_{i} = W_{ij} X_{j} + B_{i} \end{align} $$ として、  {} $$ \begin{align} \frac{\partial L}{\partial X_{j}} &= \frac{\partial L}{\partial Y_{k}} \frac{\partial Y_{k}}{\partial X_{j}} = \frac{\partial L}{\partial Y_{k}} \frac{\partial \left(W_{kl} X_{l} + B_{k} \right)}{\partial X_{j}} = \frac{\partial L}{\partial Y_{k}} W_{kl} \delta_{lj} \\ &= \frac{\partial L}{\partial Y_{k}} W_{kj} = W_{kj} \frac{\partial L}{\partial Y_{k}} = W^{\mathrm{T}} \frac{\partial L}{\partial Y} \\ \frac{\partial L}{\partial W_{ij}} &= \frac{\partial L}{\partial Y_{k}} \frac{\partial Y_{k}}{\partial W_{ij}} = \frac{\partial L}{\partial Y_{k}} \frac{\partial \left(W_{kl} X_{l} + B_{k} \right)}{\partial W_{ij}} = \frac{\partial L}{\partial Y_{k}} X_{l} \delta_{ik}\delta_{lj} \\ &= \frac{\partial L}{\partial Y_{i}} X_{j} = \frac{\partial L}{\partial Y} X^{\mathrm{T}} \end{align} $$

最右辺は行列表現で書いた。今まで悩んでたのが馬鹿らしくなるくらい簡単に計算できた。

まとめ

もう全部これで良くない?

とはいえ、実装と実行速度を考えると行列表現が得られた方が実用的であることは明白。なので、縮約記法で計算してから行列表現に戻してやるのが実務的なやり方なんだと思う。

追記:
前回記事でまとめたエディントンのイプシロンは全く必要なかった(笑)。まぁでも知っていることは良い事だ。

参考

http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/3274/pdf/imm3274.pdf

テンソルを使いこなしたい2(レヴィ=チヴィタ記号)

hotoke-x.hatenablog.com

の続き。一気に書くとまとまりがないので分けた。

本記事ではエディントンのイプシロン(レヴィ=チヴィタ記号)についてまとめる。これにより、ベクトル積、行列式、三重積などがスッキリ表現できる。

エディントンのイプシロン(レヴィ=チヴィタ記号)

添字が (1, 2, 3)の偶置換なら 1、奇置換なら -1、それ以外は0とするもの。 Wikipediaではテンソルであると言い切っているが、 レヴィ・チヴィタの記号 [物理のかぎしっぽ]では、疑テンソルであると説明している(ただ、あまり実用上は気にしなくて良さそう)。以下のように定義されている。

$$ \begin{align} \varepsilon_{i j k}=\left\{\begin{array}{ll} 1 & \left( (i, j, k)\in \left\{(1,2,3), (2,3,1), (3,1,2)\right\} \right) \\ -1 & \left( (i, j, k)\in \left\{(1,3,2), (3,2,1), (2,1,3)\right\} \right) \\ 0 & \left( \text {otherwise} \right) \end{array}\right. \end{align} $$

すなわち

$$ \begin {align} \varepsilon_{i j k}=\varepsilon_{j k i}=\varepsilon_{k i j}=-\varepsilon_{i k j}=-\varepsilon_{j i k}=-\varepsilon_{k j i} \label{prop_epsilon} \end {align} $$

3つの直交した単位ベクトル( \boldsymbol{e}_{i}, \boldsymbol{e}_{j}, \boldsymbol{e}_{k} )を使えば、スカラー三重積でも表現できる。

$$ \begin{align} \left| \begin{array}{ccc} \boldsymbol{e}_{i} & \boldsymbol{e}_{j} & \boldsymbol{e}_{k} \end{array} \right| :=\boldsymbol{e}_{i} \cdot\left(\boldsymbol{e}_{j} \times \boldsymbol{e}_{k}\right)=\varepsilon_{i j k} \label{triple} \end{align} $$

スカラー三重積の性質からも、\eqref{prop_epsilon}が成り立っていることがわかる。さらに、\eqref{triple}より また、クロネッカーのデルタを用いれば以下のような性質があることが示せる。

$$ \begin{align} \varepsilon_{i j k} \varepsilon_{l m n} &=\operatorname{det} \left[ \begin{array}{ccc} {\delta_{i l}} & {\delta_{i m}} & {\delta_{i n}} \\ {\delta_{j l}} & {\delta_{j m}} & {\delta_{j n}} \\ {\delta_{k l}} & {\delta_{k m}} & {\delta_{k n}} \end{array} \right] \\ &=\delta_{i l} \left(\delta_{j m} \delta_{k n}-\delta_{j n} \delta_{k m} \right)+\delta_{i m} \left(\delta_{j n} \delta_{k l}-\delta_{j l} \delta_{k n}\right)+\delta_{i n}\left(\delta_{j l} \delta_{k m}-\delta_{j m} \delta_{k l} \right) \label{mul_epsilon} \end{align} $$

証明は具体的に i, j, k, l, m, nに数字を入れて計算してみれば良い。例えば i=l=2の時、考えられる (j,k,m,n)のパターンは

$$ \begin{align} \left(j, k, m, n \right) \in \left\{\left(1, 3, 1, 3\right), \left(1, 3, 3, 1\right), \left(3,1, 1, 3\right), \left(3,1, 3,1\right)\right\} \end{align} $$

となる。なお、エディントンのイプシロンの定義より値が0となる組み合わせは無視した。 この時、\eqref{mul_epsilon}は 0となる項を省略すると

$$ \begin{align} \delta_{il} \left(\delta_{jm}\delta_{kn} - \delta_{jn}\delta_{km} \right) &= \delta_{22} \left(\delta_{jm}\delta_{kn} - \delta_{jn}\delta_{km} \right) \end{align} $$

となるので、\eqref{mul_epsilon}が成り立つことを簡単に確認できる。他の条件でも同様である。 また、\eqref{mul_epsilon}より以下のような性質も導ける*1

$$ \begin{align} \sum_{k} \varepsilon_{i j k} \varepsilon_{l m k} &=\operatorname{det} \left[ \begin{array}{ll} {\delta_{i l}} & {\delta_{i m}} \\ {\delta_{j l}} & {\delta_{j m}}\end{array}\right] \\ &=\delta_{i l} \delta_{j m}-\delta_{i m} \delta_{j l} \\ \sum_{j, k} \varepsilon_{i j k} \varepsilon_{l j k} &=2 \delta_{i l} \\ \sum_{i, j, k} \varepsilon_{i j k} \varepsilon_{i j k}&=6 \end{align} $$

ちょっと実践

ベクトルの外積がシンプルになる。

$$ \begin{align} \boldsymbol{a} \times \boldsymbol{b} = \sum_{i, j} \varepsilon_{i j k} a_j b_k \end{align} $$

3次元の行列式スカラー三重積)もコンパクトに表現できる。

$$ \begin{align} \operatorname{det} \left[ \begin{array}{lll} {a_{11}} & {a_{12}} & {a_{13}} \\ {a_{21}} & {a_{22}} & {a_{23}} \\ {a_{31}} & {a_{32}} & {a_{33}}\end{array}\right] &=\operatorname{det} \left[ \begin{array}{ccc} {\sum_{i} \delta_{1 i} a_{i 1}} & {\sum_{j} \delta_{1 j} a_{j 2}} & {\sum_{k} \delta_{1 k} a_{k 3}} \\ {\sum_{i} \delta_{2 i} a_{i 1}} & {\sum_{j} \delta_{2 j} a_{j 2}} & {\sum_{k} \delta_{2 k} a_{k 3}} \\ {\sum_{i} \delta_{3 i} a_{i 1}} & {\sum_{j} \delta_{3 j} a_{j 2}} & {\sum_{k} \delta_{3 k} a_{k 3}}\end{array}\right] \\ &= \sum_{i, j, k} \operatorname{det} \left[ \begin{array}{lll} {\delta_{1 i}} & {\delta_{1 j}} & {\delta_{1 k}} \\ {\delta_{2 i}} & {\delta_{2 j}} & {\delta_{2 k}} \\ {\delta_{3 i}} & {\delta_{3 j}} & {\delta_{3 k}}\end{array}\right] a_{i 1} a_{j 2} a_{k 3} \\ &=\sum_{i, j, k} \varepsilon_{i j k} a_{i 1} a_{j 2} a_{k 3} =\sum_{i, j, k} \varepsilon_{i j k} a_{1 i} a_{2 j} a_{3 k} \end{align} $$

最後の式変形は \operatorname{det} A = \operatorname{det} A^{\mathrm{T}}であることを利用した。

その他

実は添字の書き方には下付きと上付きで使い分けるルールがある。
本記事で扱う範囲を遥かに超えるので説明は省略した(というか説明できない)。
代わりに参考になるページ、資料へのリンクを貼っておく。

  1. 共変ベクトルと反変ベクトル [物理のかぎしっぽ]
    添字について解説してくれている
  2. http://takeno.iee.niit.ac.jp/~shige/math/lecture/basic2/data/detprf1.pdf
    導入は天下りだが、的p.9で添字の扱いに触れている。また、定理7で $$ \begin{align} |A B|=\left\{\begin{array}{ll} {0} & (n<m) \\ {|A||B|} & (n=m) \\ {\sum_{i_{1}<\cdots<i_{m}} \left|A_{i_{1}, \ldots, i_{m}}| \left|B^{i_{1}, \ldots, i_{m}}| \right.\right.} & (n>m) \end{array}\right. \end{align} $$ の証明をしている。 \varepsilon_{i j k} \varepsilon_{l m n} の計算の参考になる。
  3. ベクトル解析公式の証明 – 準備篇 | 高校物理の備忘録
    証明がたくさん書いてある。丁寧に途中式まである。

まとめ

使いこなすまでの道のりはまだまだ長そう...。

テンソルを使いこなしたい1(アインシュタインの縮約記法)

「ゼロから作るディープラーニング1」で、Affineレイヤの逆誤差伝播のために行列積の行列微分がシレっとでてきた。一旦書き下して納得はしたものの、行列の順序を意識しながら連鎖律を操るのはシンドイ。そこで、Einsteinの縮約記法でテンソルを扱えたらスッキリするのではと思ったので扱い方のメモ。以下がわかれば簡単な微分についても同様のノリで計算できるようになる。気が向いたら書く。

テンソルとは

ベクトルとか行列の多次元への拡張版。まとめてテンソルって呼ぶことにしましょうってノリで良いのでは(知らんけど)。

スカラー:0階テンソル
ベクトル:1階テンソル
行列  :2階テンソル

のような対応がある。

Einsteinの縮約記法

 \sum_iとか行列の次元とか、自明なものはどんどん省略しちゃおうぜってことらしい。

例)

 s    :0階テンソル
 \boldsymbol{u, v}  :1階テンソル
 \boldsymbol{T, A, B}:2階テンソル

とおくと、

$$ \begin{align} {\boldsymbol u}^T {\boldsymbol v} \longrightarrow \delta_{ij} u_{i} v_{j} = u_i v_i \label{tensor0} \end{align} $$

 \delta_{ij}クロネッカーのデルタ。直接対応しているわけではないが、この指標表示が意味するところは

$$ \begin{align} \mathrm{Tr}\hspace{1pt} \left({\boldsymbol u} {\boldsymbol v}^T \right) \end{align} $$

と同じようなものである。

$$ \begin{align} {\boldsymbol T} {\boldsymbol v} \longrightarrow T_{ij} v_j \label{tensor1} \end{align} $$

$$ \begin{align} {\boldsymbol u} {\boldsymbol v}^T &\longrightarrow u_i v_j \label{tensor2.1} \end{align} $$

$$ \begin{align} {\boldsymbol A} {\boldsymbol B} \longrightarrow A_{ij} B_{jk} \label{tensor2.2} \end{align} $$

ここで重要なのは、ベクトルや行列の成分を表す i, jなどの添字には特別な意味がない事で、次元の下限、上限を文脈から読み取ることになる。添字を表すのに無駄に文字を消費せずに済む。

また、計算結果に影響を与えないダミー指標がある。\eqref{tensor2.2}で言えば jがダミー指標で、計算の過程で総和をとるので計算結果に影響を与えない。つまり、計算のためだけに一時だけに導入された使い捨ての添字がダミー指標である。

テンソル積(Einsteinの縮約記法 ver.)

では、直積(クロス積*1、2つのベクトルのテンソル積)をEinsteinの縮約記法に書き直してみる。

$$ \begin{align} {\boldsymbol T} = {\boldsymbol u} \otimes {\boldsymbol v}^T \longrightarrow T_{ij} = {\boldsymbol u}_{i} {\boldsymbol v}_{j} = {\boldsymbol v}_{j} {\boldsymbol u}_{i} \label{outer} \end{align} $$

テンソルの成分を明示的に扱っているので行列積のように積の順序を気にしなくてよくなる。これだけでも大分スッキリする。

テンソル内積

そりゃそうだよねって形だけど念のため。

$$ \begin{align} {\boldsymbol A} \cdot {\boldsymbol B} = \delta_{ik} \delta_{jl} A_{ij} B_{kl} \end{align} $$

縮約記法から元々の記法へ

各成分へ基底をかけて和をとるだけ。
見たほうが早い。

$$ \begin{align} {\boldsymbol T} &= \left( \begin{array}{ccc} T_{11} & \ldots & T_{1n} \\ \vdots & \ddots & \vdots \\ T_{m1} & \ldots & T_{mn} \end{array} \right) \\ &= \sum_{i=1}^{m} \sum_{j=1}^{n} T_{ij} \left({\boldsymbol e}_i \otimes {\boldsymbol e}_j \right) \end{align} $$

なお、 {\boldsymbol e}_i,  {\boldsymbol e}_jは基底を表す。

まとめ

アインシュタイン先生に感謝

*1:内積に対して外積とも呼ばれることがある。ただし直積と外積は3次元では一致するが同義ではないことに注意が必要。また、outer productは直積の意味で使われるので日本語と英語の違いにも注意すること

Anaconda経由でTensorFlow (GPU)+Keras環境を構築する

タイトルの通りです。

これまでだと、GPU版TensorFlowのためにCUDA環境を構築するのがとにかく面倒でした。

良い方法はないかとググっていたら、Anaconda経由ですべてが完結するメチャクチャ簡単な方法を見つけたので紹介します。


以下の記事です。

Windows10とUbuntuの両方とも同じやり方で行けるらしい(本記事での環境はWindows10)。

towardsdatascience.com

大事なところだけ抜き出すと、プロンプト画面で

conda create --name tf_gpu tensorflow-gpu 

と入力して仮想環境を構築するだけとのこと。

tf_gpuの部分は仮想環境に付ける名前なので、自分の好きにしたら良いです。

これで後は勝手に仮想環境が構築され、numpyやTensorFlow, Kerasも勝手に入ります。

matplotlibやjupyterなどは入っていないので

activate tf_gpu

で仮想環境をアクティベートして、必要なパッケージを適宜インストールしましょう。

from tensorflow.python.client import device_lib
device_lib.list_local_devices()

を実行して、

device_type: "GPU"

が出力の中にあればインストール成功です。

関数の最適化について勉強メモ

ざっくりと最適化についてまとめる(少し詳しい目次程度)。
手法間のつながりとか流れを確認する用なので、内容には踏み込まない。
証明や細かい解説は書籍を参照すること*1

間違いがあればご指摘ください。

勾配法

  • 最急降下法
    今いる地点の勾配(1階微分)が最大になる方向へ突き進む方法。

  • ニュートン法
    最適化したい関数を2次の項までテイラー展開する(二次近似)。 近似した関数の勾配(1階微分)が最大になる方向へ突き進む方法。 2次の項まで反映されるので、最急降下法より収束が早く、二次収束することが知られている。 二階微分(ヘッセ行列)の逆行列を計算しなければならないのが難点。 (数値計算において逆行列は桁落ちが激しいので避けたいということ)

  • 共役勾配法
    最適化したい関数を2次の項までテイラー展開する(二次近似)。 近似した関数の勾配(1階微分)=0の式を作ると、二次の項でヘッセ行列xベクトルの形が出てくる。このベクトルが共役勾配である。 共役勾配は最適化したい関数の1階微分から得られる勾配とは、少しずれたベクトルになる。 ずれ分のベクトルを導入し、 共役勾配=勾配+αずれベクトル と線形和で表現し、αを求めることでヘッセ行列の逆行列を計算することを回避した。

最小二乗法

あまりに一般的すぎるので省略

非線形最小二乗法

最小二乗法は連立一次方程式の解を求める方法だった。 非線形連立方程式に拡張したもの。

  • ガウス・ニュートン
    複雑な関数の2階微分を計算するのは困難なので、ヘッセ行列を近似してからニュートン法を用いる方法。 解の近傍で与式=0と置ける場合、それを利用して二階偏微分の項を近似的に0とみなすことでヘッセ行列を2階微分を計算せずに近似できる。 解の近傍でしか使えないのが難点。

  • レーベンバーグ・マーカート法
    ガウス・ニュートン法は解の近傍でしか使えない。そこで、まず勾配法などで雑に探索を行い解の近傍に到達したらガウス・ニュートン法に切り替えるのがこの方法。

統計的最適化

  • 最尤推定

    • 出力誤差モデル
      誤差が正規分布に従うとすれば、最小二乗法と一致する。

    • 入力誤差モデル
      入力にも誤差がのっているとするモデル。 データ点から推定する直線までの距離を最小化する。

  • データ分類(クラスタリング

    • 教師なし学習
      ベイズ更新でクラスへの所属確率が収束するまでループ。EMアルゴリズムの特別な場合になっている。

    • EMアルゴリズム
      未知変数について周辺化した対数尤度関数の極値計算から自然と導出される。
      未知変数の周辺化がEステップ、極値を求めるのがMステップに対応する。
      推定パラメータが収束するまでEステップとMステップを反復する。
      この反復過程で対数尤度が単調増加することはJensenの不等式から証明できる。

線形計画法

1次不等式による制約付きで、1次式を最大化/最小化する。

  • スラック変数
    制約不等式を等式にするために導入する非負の変数。

  • シンプレックス法
    線形計画問題を組織的に解く(可能領域の頂点を全探索せずに済む)。 実装にはシンプレックス表(タブロー)を用いた実装が便利? 常に関数が増大するように探索するので、初期の解が非負ならその後の解も非負である。

    • 退化(縮退)
      n次元空間内の可能領域を表す多面体の頂点が、n+1枚以上の平面が交わっている点になっていると起こる。2次元平面で考えれば可能領域を表す直線3つ以上が1つの点を共有してる状態で、2直線の交点として考えたときに解釈が複数存在することになる。解釈を変えても関数の値が変化しない。この状態を退化(縮退)という。この時、シンプレックス法アルゴリズム的には退化の分無駄なステップを刻むことになる。退化から脱出する方法としては、辞書式順序による方法、および記号摂動法などがあるとのこと。

    • 人工変数を導入
      出発時の解が見つからない(見つかりにくい)場合に人工変数を導入する。
      最適解が求まった時に人工変数が0になっていればなかったことと同じなので問題ないという考え。
      可能領域が空集合になっていると、人工変数が0にならない(そのようなケースの判定に使える)

    • 双対原理
      書くと長いので省略。「裏返しの関係」ともいわれる?難しい話ではないので調べてみると良い(ただ、ややこしくはある)。 線形計画では、双対問題を解くと、元の問題の最適解が求まる。
      ただ、線形計画問題は簡単に解けるのであまりメリットを感じない(個人の感想)。

非線形計画法

制約条件も目的関数も一般の非線形関数になったもの。

  • 凸計画
    目的関数が上に凸、かつ可能領域も凸な問題。
    線形計画法も凸計画問題の一つ。

  • 2次計画
    制約条件がすべて線形
    目的関数が上に凸な問題
    (cvxoptというPythonパッケージがあって簡単に使える。)

ラグランジュ乗数

目的関数は上に凸、制約不等式はすべて下に凸、制約等式がすべて1次式。
この時、制約不等式、制約等式にラグランジュ乗数をかけて目的関数と合わせてラグランジュ関数をつくる。
ラグランジュの未定乗数法を理解しておくとわかりやすい。
ここから、KKT(カルーシュ・キューン・タッカー)条件がでてくる。

双対原理

線形計画の時の拡張になっている。
線形計画問題では双対問題を解くことで主問題の最適解が求まった。
一方で2次計画では、双対問題を解くことで主問題の最適解が求まるわけではない。
しかしながら主問題より簡単に解け、主問題の解を限定できる点で強力。
また、主問題と双対問題を突き合わせることで最適解かどうかの判定に使える(サポートベクターマシーンなどに利用されているらしい)。

動的計画法

変数が離散値をとる関数の最適解を効率的に求める方法。
大きな問題を小さな問題に分割して再帰的に解いていく。
最適経路問題として考えた方がわかりやすいかもしれない(個人の感想)。
最適性原理が成り立つ(どの部分解もその部分問題の最適解になっている)。
半順序しか与えられていない問題にも使えるらしい(半順序ってなに?って人は集合・位相論のテキストを読むとよい)。

  • ストリングマッチング
    文字列のマッチングを行う問題。
    評価関数を導入して最適化問題に帰着させる。
    情報理論を勉強しておくと呑み込みが早いと思う。

  • 制約のある多段階決定問題
    制約付きでもうまく変数変換すれば、見かけ上制約のない動的計画法に帰着できるパターンについて紹介してくれている。

参考書籍

www.amazon.co.jp

*1:金谷健一(2005)「これなら分かる最適化数学」共立出版株式会社

2次形式の最大値と最小値

前回の記事で、固有値固有ベクトルについて書いた。

hotoke-x.hatenablog.com

この記事中にある正方行列の対角化を式変形すると、行列 \displaystyle Aのスペクトル分解(固有値分解)が得られ、これを使って2次形式の標準形が得られる。

正方行列のスペクトル分解(固有値分解)

正方行列 \displaystyle Aの対角化は以下のような式だった。

$$ \begin{align} \boldsymbol{P^{-1}AP} &= \boldsymbol{D} \\ &= \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{align} $$

 \displaystyle P固有ベクトル \displaystyle Dは対応する固有値を対角成分に並べたものである。
これを式変形して、 \displaystyle Aのスペクトル分解を得る。

$$ \begin{align} \boldsymbol{A} &= \boldsymbol{P D P^{-1}} \end{align} $$

2次形式の標準形

正方行列のスペクトル分解を使うと、2次形式の標準形が得られる(2次形式と2次関数は別物なので注意)*。

まず、2次形式

$$ \begin{align} \boldsymbol{x^{T} A x} \end{align} $$

を考える。次に、ベクトル \displaystyle \boldsymbol{x}固有ベクトルを基底とする空間 \displaystyle \boldsymbol{P}に投げる。

$$ \begin{align} \boldsymbol{t} &= \boldsymbol{P x} \\ \boldsymbol{x} &= \boldsymbol{P^{T} t} \end{align} $$

これを2次形式の式に代入すれば

$$ \begin{align} \boldsymbol{x^{T} A x} &= \boldsymbol{\left(P^{T} t \right)^{T} A P^{T} t} \\ &= \boldsymbol{t^{T} P A P^{T} t} \\ &= \boldsymbol{t^{T} D t} \end{align} $$

となる。ただし

$$ \begin{align} \boldsymbol{P^{-1}} &= \boldsymbol{P^{T}} \\ \boldsymbol{P^{-1}AP} &=\boldsymbol{PAP^{-1}} = \boldsymbol{D} \end{align} $$

を利用した。これを2次形式の標準形と呼ぶ。ただ、固有値固有ベクトルを使って2次形式を書き直しただけである。

二次形式の最大値、最小値

 \displaystyle \boldsymbol{A}が半正定値対称行列なら、2次形式の最大値、最小値はそれぞれ \displaystyle \boldsymbol{A}の最大固有値、最小固有値に対応する。標準形を得ると、このように話がシンプルになる。

少々雑だったが、ここまでがよくある書籍の解説である。これですっきり理解できるならそれでよいが、なんで最大固有値、最小固有値が最大値、最小値に対応するのかイメージできないかもしれない。

2次形式の標準形の幾何学的解釈

これなら分かる最適化数学(金谷健一)、例題1.25、1.27より \displaystyle x^{2}+y^{2}=1の時

$$ \begin{align} f=6x^{2}+4xy+3y^{2} \end{align} $$

の最大値、最小値を考える。これを標準形に直すと

$$ \begin{align} f=2x' ^2 + 7y' ^2 \end{align} $$ となる。 \displaystyle x' \displaystyle y'固有ベクトル空間 \displaystyle \boldsymbol{P}で変換後の \displaystyle x \displaystyle y

$$ \begin{align} \boldsymbol{t} &= \boldsymbol{P x} \\ \boldsymbol{x} &= \left(\begin{array}{c} x \\ y \end{array}\right), \boldsymbol{t} = \left(\begin{array}{c} x' \\ y' \end{array}\right) \\ \end{align} $$

である。係数2と7は固有値になっている。

また、 \displaystyle x^{2}+y^{2}=1 x'^{2}+y'^{2}=1は同値である。

改めて標準形の式を眺めてみると、幾何学的には放物面に対応していることがわかる。実際にプロットしてみると

f:id:hotoke-X:20190103202346p:plain

となり、等高線が楕円になっていることがわかる。 \displaystyle \boldsymbol{x'^{2}+y'^{2}}=1に制限すれば、 \displaystyle x'=1のとき最小値2(最小固有値)、 \displaystyle y'=1のとき最大値7(最大固有値)になることがわかる。

*2次関数は2次以下の項からなる関数。2次形式は2次の項のみからなる関数。

参考書籍

  • 金谷健一(2005)「これなら分かる最適化数学」共立出版株式会社

固有値と固有ベクトル

固有値固有ベクトルについてメモ。
何をしているかはわかっているけど、何が起きているかがわかっていなかった。

  • 「固有」ってなんだよと思っている
  • 行列の掛け算で何が起きてるのかイメージがわかない

そんな人たちは参考になるかも。

数弱なので間違いがあれば教えてくれると助かります。

固有値固有ベクトルについておさらい

行列 {\displaystyle A} に対して $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$

が成り立つ{\displaystyle \boldsymbol{0}}でないベクトル{\displaystyle \boldsymbol{u}}を行列{\displaystyle A}固有ベクトル{\displaystyle \lambda}をその固有値という。
{\displaystyle A}{\displaystyle n \times n}対称行列なら、{\displaystyle A}{\displaystyle n}個の固有値{\displaystyle \lambda_1,\ldots,\lambda_n}を持つ。また、{\displaystyle \boldsymbol{u}}の各成分間は互いに直交する。

ここで、初めの式をちょっと変形する $$ \begin{equation} \left(\lambda \boldsymbol{I} - \boldsymbol{A} \right) \boldsymbol{u} = \boldsymbol{0} \end{equation} $$ 定義より{\displaystyle \boldsymbol{u}} \neq  \boldsymbol{0}より $$ \begin{equation} | \lambda \boldsymbol{I} - \boldsymbol{A}| = \boldsymbol{0} \end{equation} $$

この式を固有方程式(特性方程式)と呼ぶ。

つまりなんなのか

行列{\displaystyle A}を左からかけることを、行列{\displaystyle A}を作用させるなんて物理では言うらしい。
つまり、ベクトルを伸ばしたり、縮めたり、回転したりして別のベクトルへ変換することを意味している。

ここで改めて式を眺めてみる。 $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$ この式が意味するところは、「行列{\displaystyle A}による{\displaystyle \boldsymbol{u}}の変換は、結局{\displaystyle \boldsymbol{u}}{\displaystyle \lambda}倍することと同じ」ということ。
別の言い方をすれば、「行列による変換をしても、長さが{\displaystyle \lambda}倍変化するだけで向きが変わらないベクトル{\displaystyle \boldsymbol{u}}がある」ということで、そんなベクトルを固有ベクトルと呼んでいる。

基底変換

今、デカルト座標系にベクトル $$ \begin{align} \boldsymbol{v} = \left( \begin{array}{c} x \\ y \end{array} \right) \end{align} $$ が存在して、座標変換によって新たな表現を得ることを考える。
これはデカルト座標での表現なので、基底ベクトル{\displaystyle \boldsymbol{e}_x}, {\displaystyle \boldsymbol{e}_y}を使って大げさに書けば $$ \begin{align} \boldsymbol{v} &= x\boldsymbol{e}_x + y\boldsymbol{e}_y \\ &= x\left( \begin{array}{c} 1 \\ 0 \end{array} \right) +y\left( \begin{array}{c} 0 \\ 1 \end{array} \right) \\ &= \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right) \left(\begin{array}{c} x \\ y \end{array}\right) \end{align} $$ となる。
この基底ベクトルを $$ \begin{align} \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right)\left(\begin{array}{c} x \\ y \end{array}\right)=\left(\begin{array}{cc} p_x & q_x \\ p_y & q_y \end{array}\right)\left(\begin{array}{c} x' \\ y' \end{array}\right) \end{align} $$ となるような別の基底ベクトル $$ \begin{align} \boldsymbol{p}=\left(\begin{array}{c} p_x \\ p_y \end{array}\right), \;\; \end{align} \boldsymbol{q}=\left(\begin{array}{c} q_x \\ q_y \end{array}\right) $$ で置き換えて、xy座標(デカルト座標)以外の任意のpq座標で表現できることがわかる。
ここでまた、最初の式をまた眺めてみる。 $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$ この式の正方行列{\displaystyle A}は、ベクトルの変換意味しているのだった(ベクトル空間の線形写像)。
説明のために、 $$ \begin{equation} \boldsymbol{A} = \left(\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array}\right), \;\; \boldsymbol{u} = \left(\begin{array}{c} u_1 \\ u_2 \end{array}\right) \end{equation} $$ とすると、 $$ \begin{align} \boldsymbol{A u}=\left(\begin{array}{cc} a_{11}\\ a_{21} \end{array}\right)u_1 +\left(\begin{array}{cc} a_{21} \\ a_{22} \end{array}\right)u_2 \end{align} $$ と書けることから、行列{\displaystyle A}の列成分は変換後の空間での基底ベクトルを表していたことがわかる。 つまり、{\displaystyle (u_1, u_2)}という数字をそのまま行列{\displaystyle A}による変換後の空間に渡すとどうなるか計算していただけのことだった。

正方行列の対角化

さて、行列{\displaystyle A}による変換前後で固有ベクトル{\displaystyle \boldsymbol{u}}はスケールが固有値{\displaystyle \lambda}倍変わるだけで向きが変わらないのであった。
ということは、

  1. 行列{\displaystyle A}による固有ベクトルのスケール変化だけを事前に調べておく(行列Aの固有ベクトル固有値を調べておく)
  2. 変換後の空間を、固有ベクトルを基底として表現する

という手順で行列{\displaystyle \boldsymbol{A}}による変換は簡単な式に置き換えられそうだ。

{\displaystyle n\times n}正方行列{\displaystyle A}によってn次元ベクトル{\displaystyle \boldsymbol{s}}が別のn次元ベクトル{\displaystyle \boldsymbol{t}}に変換される、 $$ \boldsymbol{t} = \boldsymbol{A s} $$ という変換を考える。ここに、行列{\displaystyle A}の右側に固有ベクトルを並べた行列 $$ \begin{align} P&=\left(\boldsymbol{p}_1, \ldots, \boldsymbol{p}_n\right) \\ \boldsymbol{p}_1 &= \left(\begin{array}{c} p_{11} \\ \vdots \\ p_{1n} \end{array}\right), \;\; \boldsymbol{p}_2 = \left(\begin{array}{c} p_{21} \\ \vdots \\ p_{2n} \end{array}\right), \;\; \cdots, \;\; \boldsymbol{p}_n = \left(\begin{array}{c} p_{n1} \\ \vdots \\ p_{nn} \end{array}\right) \end{align} $$ をねじ込んでみると $$ \boldsymbol{t} = \boldsymbol{(AP)(P^{-1}s)} $$ となる。{\displaystyle \boldsymbol{AP}}の部分は、一番最初に登場した式 $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$ の左辺に対応している。異なるのは、固有ベクトルをひとつずつ扱っているのではなく、行列として並べてしまったことくらいである。じゃあ右辺はどうなってるのと思うかもしれないが、固有値と対応する固有ベクトルの積が以下のように並んでいるだけだ。 $$ \begin{equation} \boldsymbol{A P} = A \left(\boldsymbol{p}_1, \ldots, \boldsymbol{p}_n\right) = \left(\lambda_1 \boldsymbol{p}_1, \ldots, \lambda_n\boldsymbol{p}_n\right) \end{equation} $$

話を戻す。
$$ \boldsymbol{t} = \boldsymbol{(AP)(P^{-1}s)} $$ の式をよく見てみると、{\displaystyle \boldsymbol{s}\longrightarrow\boldsymbol{t}}とする変換だったのに、{\displaystyle \boldsymbol{P^{-1}s}\longrightarrow\boldsymbol{t}}のように見えてしまう。{\displaystyle \boldsymbol{AP}}を使いたいために格好が悪くなってしまった。

いっそのこと、{\displaystyle \boldsymbol{P^{-1}s}\longrightarrow\boldsymbol{P^{-1}t}}と考えてしまえばよさそうだ。つまり、{\displaystyle \boldsymbol{s}}{\displaystyle \boldsymbol{t}} をそれぞれ{\displaystyle \boldsymbol{P^{-1}}} が基底になっている空間に投げてしまってから考えるということだ。もし、元の座標系での値が知りたかったら{\displaystyle \boldsymbol{P}}を掛ければ元通りなので問題はない。 式で書くと

$$ \begin{equation} \boldsymbol{P^{-1} t} = \boldsymbol{(P^{-1}AP)(P^{-1}s)} \end{equation} $$

となる(有名な式らしい...知らんけど)。

なんだか大事になっているように見えるが、実はそうではない。\displaystyle \boldsymbol{P^{-1}AP}の部分が実は、固有値を対角成分に並べた対角行列になるからだ(\displaystyle \boldsymbol{P^{-1}AP}を対角行列にする操作を対角化と呼ぶ)。

$$ \begin{align} \boldsymbol{(P^{-1}AP)} &= \boldsymbol{D} \\ &= \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{align} $$

つまり、行列{\displaystyle A}によるベクトルの変換は、{\displaystyle \boldsymbol{P^{-1}}}が基底となる空間へ投げてしまえば、そこでは基底ベクトルの方向へ固有値倍しただけだとみなせるということである。

固有ベクトルを基底とした空間へ投げるんじゃないの?」

と思われるかもしれないがこれは書き方の問題で、結局

$$ \begin{align} \boldsymbol{P^{-1}AP} =\boldsymbol{PAP^{-1}} = \boldsymbol{D} \\ \end{align} $$

なので、

$$ \begin{equation} \boldsymbol{P^{-1} t} = \boldsymbol{(P^{-1}AP)(P^{-1}s)} \end{equation} $$

でも

$$ \begin{equation} \boldsymbol{Pt} = \boldsymbol{(PAP^{-1})(Ps)} \end{equation} $$

でも同じことである。「固有ベクトルを基底とした空間に投げてから考える」という意味では後者の方がスッキリしているのかもしれない。

【数理統計】Fisher情報量とクラメール・ラオ下限(CR下限)

Fisher情報量って何なのか理解したかっただけ。
Fisher情報量がわかれば、クラメール・ラオ下限(不偏推定量の分散の下限)がわかる。
備忘録として簡単にまとめておく。


前提として、不偏推定量についてだけ考えるものとする。
不偏推定量に関してはこちら

hotoke-x.hatenablog.com

定量が不偏性を持てば、その推定量の最適性・適切性の判定が容易になる。 互いに独立な標本{\displaystyle X_1, \ldots, X_n}の同時密度関数が{\displaystyle f(x|\theta)}で与えられるとする。すなわち $$ \begin{align} X_i \in \boldsymbol{X} &, \quad i \in \left\{1, \ldots, n \right\} \\ f(\boldsymbol{x}|\theta) := \prod_{i=1}^n p(x_i; \theta) &,\quad \boldsymbol{x} = \left\{x_1, \ldots, x_n \right\} \end{align} $$

このとき、 $$ \begin{align} \forall \theta \in \Theta,\qquad \boldsymbol{E} \left(\hat \theta_n \right) \end{align} $$

を満足する。これを{\displaystyle \theta}について微分 $$ \begin{equation} \frac{\partial}{\partial \theta} \boldsymbol{E} \left(\hat \theta_n \right) = 1 \end{equation} $$

これを展開 $$ \begin{equation} \frac{\partial}{\partial \theta} \int f(\boldsymbol{x}|\theta) ~ \hat \theta_n ~ \mathrm{d}\boldsymbol{x} = \int \hat \theta_n \frac{\partial}{\partial \theta} f(\boldsymbol{x}|\theta) \mathrm{d}\boldsymbol{x} \label{eq:fisher1} \end{equation} $$

ここで、 $$ \begin{align} \frac{\partial}{\partial \theta} \log{L} &= \frac{1}{L} \frac{\partial L }{\partial \theta} \\ \frac{\partial L}{\partial \theta} &= L \frac{\partial}{\partial \theta} \log{L} \end{align} $$

より

$$ \begin{align} \mathrm{R.H.S} &= \int \hat \theta_n \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} f(\boldsymbol{x}|\theta) \mathrm{d} x \\ &= \boldsymbol{E} \left[\hat \theta_n \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right] \\ &= Covariance \left(\hat \theta_n, \frac{\partial \log{f(\boldsymbol{x} | \theta)} }{\partial \theta} \right) \\ &\leq \sqrt{V \left(\hat \theta_n\right)} \sqrt{I_n\left(\theta \right)} \label{eq:fisherinfo} \end{align} $$

ここで、{\displaystyle I_n\left(\theta \right)}はFisher情報量と呼ばれる量で、クラメール・ラオの下限(CR下限)と密接な関係がある(後述)。なお、

$$ \begin{align} Covariance \left(\hat \theta_n, \frac{\partial \log{f(\boldsymbol{x} | \theta)} }{\partial \theta} \right) \leq \sqrt{V \left(\hat \theta_n\right)} \sqrt{I_n\left(\theta \right)} \end{align} $$ ではコーシー・シュワルツ不等式を利用した。また、

$$ \begin{align} I_n\left(\theta \right) &= V \left(\frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right) \\ &= \boldsymbol{E} \left[\left( \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta}\right) ^2 \right] \end{align} $$

特に{\displaystyle n}個の標本が{\displaystyle i.i.d}なら $$ \begin{align} I_n \left( \theta \right) &= \boldsymbol{E} \left[ \left( \frac{\partial}{\partial \theta} \sum_{i=1}^{n} \log{f(\boldsymbol{x}|\theta)} \right) ^2 \right] \\ &= n \boldsymbol{E} \left[\left(\frac{\partial}{\partial \theta} \log{f(\boldsymbol{x_1}|\theta)} \right) ^2 \right] \\ &= nI_1 \left( \theta \right) \end{align} $$

ここで $$ \begin{equation} \frac{\partial}{\partial \theta} \boldsymbol{E} \left(\hat \theta_n \right) = 1 \end{equation} $$

を思い出せば $$ \begin{align} Covariance \left(\hat \theta \left(\boldsymbol{x} \right), \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right) &= 1 \\ Covariance \left(\hat \theta \left(\boldsymbol{x} \right), \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right) &\leq \sqrt{V \left(\hat \theta_n \right)} \sqrt{I_n \left(\theta \right)} \end{align} $$

より $$ \begin{align} 1 &\leq \sqrt{I_n \left(\theta \right)} \\ V \left(\hat \theta \right) &\geq \frac{1}{nI_1 \left(\theta \right)} \label{eq:CR} \end{align} $$

となる。以上より、フィッシャー情報量{\displaystyle I_1 (\theta)}によって推定量の分散の下限が決まることがわかる(推定量としての良さがわかる)。この下限をクラメール・ラオ下限(CR下限)と呼ぶ。

【Julia】Julia-1.0.0のProjectを試す【仮想環境?】

Julia1.0.0が遂にリリースされました。

f:id:hotoke-X:20180810011534p:plain

仮想環境っぽいものを作って環境の切り替えができるようなので、試してみます。

Projectの作成

Anaconda Pythonではconda createコマンドで仮想環境を作れました。

Julia1.0.0ではProjectと呼ぶようです。

Anacondaでは

conda create -n envname

とする必要がありましたが、Juliaではプロジェクトのディレクトリに入って

activate .

とするようです。

その後、

julia>]
(projectname)pkg> add hogehoge

で、プロジェクトローカルな環境にhogehogeパッケージが入ります。

その後バックスペースでpkg環境を抜けて同フォルダで作業すれば良いっぽい。


ただ、Julia1.0.0がリリースされたのが記事作成時点で昨日の今日なので、対応が追い付いているパッケージを探すのが結構大変です(2018/08/10現在)。

私は情弱なのでおとなしく修正を待つことにします。

(ちなみに現状遊ぶなら0.6が無難)