Julia 小技集

完全に自分用
便利な小技をメモしていく

随時更新。

現在の環境

Julia Version 1.1.0
Commit 80516ca202 (2019-01-21 21:24 UTC)
Platform Info:
  OS: Linux (x86_64-pc-linux-gnu)
  CPU: AMD Ryzen 7 1700 Eight-Core Processor
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, znver1)
Environment:
  JULIA_DEPOT_PATH = /opt/julia
  JULIA_PKGDIR = /opt/julia
  JULIA_VERSION = 1.1.0

• 型を一段階アップグレードする関数
  例：Int64 -> Int128

wilden(x)

• Array{Array{hoge, 1}, 1}を2次元配列に直したい

vcat(a...)

• 配列の要素ごとに関数を適用

map([1,2,3]) do x
    何らかの処理
end

xに配列の要素を代入してdo-end間の処理をしてくれる。

• 多次元配列の成分ごとに関数を適用

mapslices(関数, 対象配列,dims=次元の指定)

次元の指定では[1,2]のように複数の次元を指定できる。

Meta programming

動的に構造体を定義

discourse.julialang.org

例）

macro make_struct(struct_name, schema...)
    fields=[:($(entry.args[1])::$(entry.args[2])) for entry in schema]
    esc(quote mutable struct $struct_name
        $(fields...)
        end
    end)
end
println(@macroexpand @make_struct STRUCT_NAME (x,Int32) (y, Float32))
@make_struct STRUCT_NAME (x, Int32) (y, Float32)

その他

デバッグ

デバッグにはRevise.jl(https://t.co/4v48ihKEOR)が便利ですよ。REPLを起動したままファイルを更新すれば自動で変更が反映され，再コンパイルも最小限で済みます。

— (｢・ω・)｢ｶﾞｵｰ (@bicycle1885) 2018年10月12日

マルチGPU

忘れないうちに昨日やったことを書いた（令和初の記事）
JuliaでもMulti-GPUsしたい - part of something https://t.co/rHKiO5TUzi @astellon_musicさんから

— Astellon (@astellon_music) 2019年4月30日

• 関数のベクトル、多次元配列への拡張

function f(x)
    xに対する処理
end


function g(x::Vector{T})  where {T<:AbstractFloat}
    xに対する処理
end

# 行列の成分ごとに上記の関数gを適用
function g(x::Matrix{T}) where {T<:AbstractFloat}
    mapslices(g, x, dims=次元の指定)
end

構造に対しては明示的に、数値型に関してはパラメトリック型で定義すると使いやすい（個人の感想） f(x)を要素ごとに適用したいならf.(x)とすれば良い。
@viewなどで参照を渡す場合は具体型ではなく抽象型で定義したほうが良い（SubArrayがわたされるので）

【初心者】Deep Learning用語まとめ

「ゼロから作るDeep Learning」を読んでいる。 一気読みできればいいが、そういう訳ににもいかないので簡単な用語を忘れることがある。 そんな時に本を開くのは億劫なので、ここに簡単な用語をまとめておく。

---

ニューラルネットワークに入る前に...

• パーセプトロン
  多入力単出力の構造 信号は0, 1 (ステップ関数) 単純パーセプトロンでAND, NAND, ORゲートを作れる。（多層パーセプトロンにすればXORも可能）

• バイアス
  上流からの重み付き和とは別に常に加算される値

• 白色化（whitening）
  データ全体の分布の形状を均一にする方法

• one-hot-label
  ラベルの表現方式のひとつ。
  ラベルが0~9まであり、正解ラベルが7の時、[0,0,0,0,0,0,1,0,0]のように表示する。

ニューラルネットワーク

• 活性化関数
  入力信号の総和を出力信号へ変換する。単純パーセプトロンではステップ関数を使っていた。ステップ関数以外を使えばニューラルネットワークと呼ぶようになるらしい（ゼロから作るDeep Learning, p44より）。活性化関数には線形関数を用いてはならないことに注意（多層にする意味がなくなる。単純な合成写像を考えてみればよい。）。

  • ステップ関数
    閾値を超えたら1、それ以外は0を出力する関数。
    all-or-nothingな信号という点で一番ニューロンっぽい気がする（個人の感想）

  • シグモイド関数
    $$ \begin{align} \frac{1}{1+\exp(-x)} \end{align} $$ 値を[0, 1]に押し込めるのが嬉しい関数。
    滑らかな関数（微分可能）であるところもポイント。

  • ReLU (Rectified Linear Unit) 関数
    閾値を超えたら入力をそのまま出力、それ以外は0を出力する関数。

• 出力層の活性化関数
  出力層だけ活性化関数を特別扱いする。

  • 恒等関数
    入力をそのまま出力する。回帰問題に用いる。

  • ソフトマックス関数
    を入力信号、を出力信号として
    $$ \begin{align} y_k=\frac{\exp(a_k)}{\sum_{i=1}^{n} \exp(a_i)} \end{align} $$
    これも出力を[0, 1]に押し込めるのがうれしい。
    オーバーフロー対策に定数 を導入して $$ \begin{align} y_k = \frac{\exp(a_k-C)}{\sum_{i=1}^{n}{\exp(a_i-C)}} \end{align} $$
    とするのが良いらしい。にはの最大値を持ってくることが一般的。

損失関数

ニューラルネットワークの学習に用いる損失関数いろいろ。 以下ではyを出力（ラベルごとの確率）、tを正解ラベルのone-hot-labelとする。kはデータの次元数とする。

• 二乗和誤差
  説明不要。お馴染み。
  $$ \begin{align} E = \frac{1}{2} \sum_{k} (y_k-t_k)^2 \end{align} $$

• 交差エントロピー誤差
  に関するの期待値を計算していることになる。情報理論の観点で言えば、に関するの平均符号長を計算していることになる。情報理論ででてくるエントロピーと異なり自然対数を用いているが、本質的には同じものなので注意。今回は細かい説明はしないので、ちゃんと理解したい場合は符号化理論を勉強しましょう。
  $$ \begin{align} E = - \sum_{k} t_k \log y_k \end{align} $$

ニューラルネットワークの学習

• ミニバッチ学習
  ニューラルネットワークの学習において、全データについて損失関数を計算するのは現実的でない。
  そこで一部のデータをランダムサンプリングして損失関数を求め、それを全データの代表値とみなして学習を行う方法のこと。
  ランダムサンプリングしているといっても、そもそもどうやってデータが集められたのかには注意する必要があるなと思った（個人の感想）。

• 勾配法
  勾配法については以下を参照

hotoke-x.hatenablog.com

• 確率的勾配降下法（SGD; Stochastic Gradient Descent）
  ミニバッチとしてランダムサンプリングされたデータを使用した勾配降下法のこと。
  以下の式に従ってパラメータの更新を行う。 $$ \begin{align} \boldsymbol{W} \leftarrow \boldsymbol{W}- η\frac{\partial L}{\partial \boldsymbol{W}} \label{s_g_d} \end{align} $$ は損失関数、は学習率である。

• エポック（epoch）
  学習の単位のこと。
  訓練データをすべて使い切った時の回数に対応。
  1000個のデータに対して50個のミニバッチで学習する場合、20回のミニバッチ学習=1 epochとなる。
  エポックごとに訓練データの認識精度とテストデータの認識精度を比較する。この時、認識精度に差があると過学習が起きていることになる。

誤差逆伝播法

数値微分を用いたニューラルネットワークの学習は実装が容易。一方で計算に時間がかかる。誤差逆伝播法で勾配の計算を効率よく行う。

• 計算グラフで考える
  ノードごとの局所的な計算で済む。 局所的な計算結果を保持できるので、逆伝播によって微分が計算できる。
  ただし、バッチごとに計算を行っている場合行列演算になるので、行列積が可換でないことに注意。
  縮約記法を使えば簡単（以下の記事を参照）。

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

hotoke-x.hatenablog.com

• 勾配確認
  数値微分による逆誤差伝播の実装と、解析解による逆誤差伝播の実装の値を比べる事。
  解析解による実装を行った時、解が間違っていないことを確認できる。

Momentum

SGDの\eqref{s_g_d}に慣性項を導入した。前の状態を引き継ぐことでSGDの欠点の一つであるパラメータの過剰更新を抑制する。パラメータ更新の式は以下の通り。

$$ \begin{align} \boldsymbol{v}_{t+1} &\leftarrow \alpha \boldsymbol{v}_{t} - η \frac{\partial L}{\partial \boldsymbol{W}} \\ \boldsymbol{W} &\leftarrow \boldsymbol{W} + \boldsymbol{v}_{t+1} \end{align} $$

AdaGrad

Adaptive Gradientの略？
学習率を減衰させていく方法。初めは大雑把に、時間が経つにつれ小さく探索する。パラメータの更新式は以下の通り。

$$ \begin{align} \boldsymbol{h} &\leftarrow \boldsymbol{h} + \frac{\partial L}{\partial \boldsymbol{W}} \odot \frac{\partial L}{\partial \boldsymbol{W}} \\ \boldsymbol{W} &\leftarrow \boldsymbol{W} -η \frac{1}{\sqrt{\boldsymbol{h}}} \frac{\partial L}{\partial \boldsymbol{W}} \end{align} $$

ここで、はアダマール積（要素ごとの積）である。が大きい要素は、次の学習率がより小さくなるように調整される。ループによって学習率が下がる一方なので、最終的には学習が進まなくなるのが欠点。これを改善したものにRMSPropという方法があるらしい。

Adam

MomentumとAdaGradの融合（AdaptiveなMomentumということだと思う）。
詳細は原著*1を読まないとわからない。

重みの初期値問題

小さい値をランダムに与える。
過学習を防ぐために小さい値で重みを学習したいという願いが込められている。
PRMLでみた正則化に似ていると思った（制約を与えたわけでは無いが、事前分布に願いを込める意味で）。

勾配消失

例えばシグモイド関数の値が0, 1付近だと、微分値は小さくなる。アクティベーションが0, 1に偏っていると、逆誤差伝播の時に勾配が消失し、学習が進まなくなる。特定の値にアクティベーションの値が集中するのも良くない。ノード全体が同じ値を出力するということなので、ノードが複数存在している意味が無くなってしまう。

• Xavierの初期値
  前層ノードの個数を使っての標準偏差を持つ分布を用いる（なんだか中心極限定理と関係がありそう）。 活性化関数が線形であることが前提らしい。sigmoid関数、tanh関数は対称で中央付近が線形とみなせるとのこと（それでいいの？？）。詳細は原著参照*2。

• Heの初期値
  ReLUを用いる場合はHeの初期値が良いとの事。の標準偏差を持つ分布を用いる。ReLUは負の値が全部0に押し付けれられるので、Xavierの初期値よりも広がりを持たせるということらしい（ReLUが青天井の値を持つことを考えると3でも5でも良い気がする...）。 こちらも詳細は原著参照*3

• Batch Normalization
  各層ごとにアクティベーション分布に広がりを持たせたいなら、各層ごとに調整したらいいじゃん。と思ったら書いてあった。データと各層ごとの出力をに標準化する。 活性化関数の前と後どちらでBatch Normalizationするかは議論になっているらしい。
  個人的には活性化関数の後でBatch Normalizationの方が自然な気がしている。理由は、各層の入力に広がりを持たせたいのが目的なら、その直前で分布の調整をした方が自然に感じるから。活性化関数の前で標準化してしまうと、活性化関数によって再び分布が偏る可能性が排除できない。だからこそXavierの初期値では線形関数である前提が必要だったのだと思う。活性化関数の後なら活性化関数が非線形でも使えるかもしれない（という期待を込めた感想）。
  histogram equalizationとか使ったらどうなるんだろう..。アルゴリズム的に破綻するかもしれないけど。

正則化

過学習を抑制する方法。

過学習

訓練データのみに最適化されてしまって、テストデータに全然対応できない状態。

Weight decay

損失関数に、重みのL2ノルム正則化項 $$ \begin{align} \frac{1}{2}\lambda\boldsymbol{W}^2 \end{align} $$ を加算して学習を行う。

はハイパーパラメータ。

（ABICのノリでの自動決定とかできないだろうか）

Dropout

ノードをランダムに除去しながら学習する方法。ノードを除去することで毎回異なるモデルを生成するので疑似的なアンサンブル学習として見れるとのこと。

検証データ（validation data）

ハイパーパラメータを調整するためのデータ。 テストデータでハイパーパラメータを調整すると、ハイパーパラメータがテストデータに関して過学習を起こすので別に用意する必要がある。

ハイパーパラメータの最適化

書籍を読んだ感じでは、グリッドサーチっぽかった。ただ、ハイパーパラメータ探索ではグリッドサーチよりもランダムに探索したほうが良い結果になるとの事。ハイパーパラメータごとに結果に与える影響度が異なるからだとか。ハイパーパラメータを探索するスケールが適当でないことが多いってことかな？

ベイズ的に最適化する方法もあるらしい*4。MCMC（マルコフ連鎖モンテカルロ法）かABC（近似ベイズ計算）だろうか。パラメトリックに最適化するなら事前分布の設定がまた問題になりそう。何はともあれ、これも原著を読まないことにはわからない。

畳み込みニューラルネットワーク（Convolutional neural network: CNN）

全結合層だと、データの位置情報を無視して学習してしまう。 CNNでは位置情報も学習できる（ということになっている）。

Convolution レイヤ

フィルタを通すことで、出力サイズを落とす（次元という言い方は正しくないかもしれない）。 例えばMNISTなら、28x28だった画像をフィルタを通すことで26x26にするような操作。 3次元データ（例えばRGB画像）に対してフィルタを適用する場合は、フィルタも3次元にし、出力される3層のz軸方向に和をとる（なので出力は2次元）。

フィルタを複数個用意すれば出力は3次元になる。バッチの個数分同じフィルタを用意すれば一括処理が可能になる。

パディング

Convolution レイヤばかり通していたら出力サイズが小さくなる一方なので、入力行列の外側を何らかの値（例えば0）で埋めて、畳み込みによる出力サイズの現象を防ぐ方法。

ストライド

フィルターを適用する位置の間隔

Pooling レイヤ

空間を小さくする演算処理。
Poolingのウィンドウサイズとストライドは同じ値にするのが通常。
学習パラメータがない。
チャネル数も変化しない。
微小な位置変化に対してロバスト。ちょっとぐらい入力データが動いても、フィルタの特性（フィルタ適用範囲の最大値をとってくるなど）によって出力結果が変化しにくくなる。

画像認識ではMax poolingが一般的らしい。Average poolingもあるが、確かに外れ値に引っ張られるのはちょっと使いにくそう（素人的感想）。 入力をFFTした画像とかにしたらどうなるのだろう。ノイズはノイズで学習しやすくなる気がする（これも素人的感想）

代表的なCNN

LeNet

• 畳み込み層とプーリング層を重ねて、最後に全結合層を使う。
• 活性化関数はReLU
  オリジナルではプーリング層はただサブサンプリングするだけで、活性化関数もsigmoidだった。現在はMax poolingが主流らしい。

AlexNet

• 活性化関数にReLUを使用
• 局所正規化を行う（上で書いたBatch Normalizationのようなものだと思う）。
• Dropoutを使用する

Deep Learning

VGG*5

• ReLUを使用
• 重みの初期値はHeの初期値を使用（ReLUを使用するから？）
• 小さなフィルタ（3x3）による畳み込み（些細な特徴でも拾いたいから？）
• Adamでパラメータ更新
• 全結合層の後にDropout（上層でメチャクチャ特徴を学習させて、最後にアンサンブル学習することで変な学習を回避？？）

データを回転させたり動かしたりして疑似的にデータを増やす「データ拡張」も認識精度向上に貢献するとのこと。

層を増やす利点

1. 単純にパラメータが減る。
   5x5のフィルターによる出力と、3x3のフィルタを2回通すのとでは、出力サイズは同じだが、後者の方がフィルタの要素数が小さくて済む。
2. より少ないデータで学習ができる 層が深いと高度な情報が抽出できるため

他のモデルの例

• GoogLeNet*6
  VGGとは違って、層のつながり自体が2次元的（インセプション構造というらしい）。1x1フィルターの畳み込みでサイズを削減することによって計算速度を向上させているとの事。
  ただ、層のつながりを2次元的にする恩恵はなんなのだろうか...。これも要原著参照ということで。

• ResNet*7
  入力層を畳み込み層をすっ飛ばして、出力と合算させる構造が入っている。
  層を深くしすぎて抽出できる特徴がスッカスカになって来たころに、昔を思い出させようとするお気持ちなのだろうか。もしそうなら、層の深ささえ適切なら同じような結果が得られそうな気がした。また、合算処理によって新たな特徴マップが得られるということなら、データ拡張でも似たような結果が得られそうだとも思った。知らんけど。

物体検出の例

1. R-CNN*8
   高速版もあるらしい*9

セグメンテーションの例

1. FCN*10
   全結合層が畳み込み層からなる。

画像キャプション生成

代表的なものとしては、NIC*11（Recurrent Neural Netwotk (RNN)の一種）がある。 再帰的な構造があり、時系列や言語などの連続データに使われる。

画像スタイル変換

絵のテイストを変えたりする*12

画像生成

スタイル変換とは異なり、完全に新しい画像を生成する。
DCGAN*13

Deep Q-Network（DQN, 強化学習の一つ）

1. エージェントが行動する
2. 行動に対する報酬を得る これを繰り返してずーっと学習させる。
   DQNでは最適行動価値関数とかいう関数をCNNで近似するらしい。これも論文読まないとわからないかも*14。

宿題

• 2層ニューラルネットワークの実装
• パーセプトロンの収束定理
  有限回の学習によって線形分離可能な問題が解ける。
• MITのディープラーニングの講義 アーカイブもあるので2017、2018年のも観れる。
• AlphaGOの論文も読めるようにする*15。

参考書籍

斎藤 康毅（2016）「ゼロから作るDeep Learning-Pythonで学ぶディープラーニングの理論と実装」株式会社オライリー・ジャパン

*1:D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.

*2:X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” Proceedings of the thirteenth international conference on artificial intelligence and statistics, 2010.

*3:K. He, X. Zhang, S. Ren, and J. Sun, “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification,” Proceedings of the IEEE international conference on computer vision, 2015.

*4:J. Snoek, H. Larochelle, and R. Adams, “Practical bayesian optimization of machine learning algorithms,” Advances in neural information processing systems, 2012.

*5:K. Simonyan and A. Zisserman, “Very Deep Convolutional Networks for Large-Scale Image Recognition,” arXiv preprint arXiv:1409.1556, 2014.

*6:C. Szegedy, W. Liu, Y. Jia, P. Sermanet, S. Reed, D. Anguelov, D. Erhan, V. Vanhoucke, and A. Rabinovich, “Going deeper with convolutions,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2015.

*7:K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2016

*8:R. Girshick, J. Donahue, T. Darrell, and J. Malik, “Rich feature hierarchies for accurate object detection and semantic segmentation,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2014.

*9:S. Ren, K. He, R. Girshick, and J. Sun, “Faster r-cnn: Towards real-time object detection with region proposal networks,” Advances in neural information processing systems, 2015

*10:E. Shelhamer, J. Long, and T. Darrell, “Fully Convolutional Networks for Semantic Segmentation,” Ieee T Pattern Anal, vol. 39, no. 4, pp. 640–651, 2016.

*11:O. Vinyals, A. Toshev, S. Bengio, and D. Erhan, “Show and tell: A neural image caption generator,” Proceedings of the IEEE conference on computer vision and pattern recognition, 2015.

*12:L. Gatys, A. Ecker, and M. Bethge, “A neural algorithm of artistic style,” arXiv preprint arXiv:1508.06576, 2015.

*13:A. Radford, L. Metz, and S. Chintala, “Unsupervised representation learning with deep convolutional generative adversarial networks,” arXiv preprint arXiv:1511.06434, 2015.

*14:Mnih, K. Kavukcuoglu, D. Silver, A. Rusu, J. Veness, M. Bellemare, A. Graves, M. Riedmiller, A. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, S. Legg, and D. Hassabis, “Human-level control through deep reinforcement learning,” Nature, vol. 518, no. 7540, p. 529, 2015.

*15:D. Silver, A. Huang, C. Maddison, A. Guez, L. Sifre, G. van den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, S. Dieleman, D. Grewe, J. Nham, N. Kalchbrenner, I. Sutskever, T. Lillicrap, M. Leach, K. Kavukcuoglu, T. Graepel, and D. Hassabis, “Mastering the game of Go with deep neural networks and tree search,” Nature, vol. 529, no. 7587, p. 484, 2016.

Windows10+Anaconda+PyTorch

Chainerを使いたかったけどWindowsは公式サポートでない上に、CUDA、cuDNNのインストールがやたら面倒くさい...。 他のライブラリでもっと簡単に使えるものはないかと探していたら、PyTorchならAnaconda経由で環境が整うらしいのでメモ。
linuxでもmacOSでも同様の手順で行けると思います（試してない）。

ちなみにtensorflowもAnacondaだけで環境が整います。

hotoke-x.hatenablog.com

---

以下の記事によると「PyTorchはCUDAをマニュアルインストールする必要はないよ！」とのこと。

discuss.pytorch.org

インストール方法

前提としてAnacondaのインストールは終わっているものとします。

PyTorchの公式

で、以下のようにポチポチと自分の環境を選択します。

pytorch

「CUDA」のところは、自分が使いたいバージョンを選べば良いようです。すると、「Run this Command:」にcondaコマンドが表示されるのでコピペして実行しましょう。これだけです。

正しくインストールできたかは以下のpythonコードで確認できます。

import torch 
print(torch.cuda.is_available())

Trueが返ってくれば正常にインストールできています。

CUDAのバージョンとかドライバのバージョンとか、ubuntuのカーネルとか気にしなくてよいのが最高にうれしいですね。

テンソルを使いこなしたい３（行列の行列微分）

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

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

行列の行列微分

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

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

要は、はで微分したときはそれ以外はと言っているだけである。「じゃあ全成分になるじゃん」とか一瞬思うかもしれない（自分は思った）がそれは大間違い。それぞれについてを走査して計算するので4階テンソルになっている（添字も4つあるし）。

逆誤差伝播

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

これを式で書けば

$$ \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} $$

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

ナニコレ。特に。そりゃの次元に一致しなきゃいけないから逆算すればわかるけど気持ち悪い。一方で、も連鎖律としてはスッキリしているが、ベクトルの行列微分がある。しかもについてもなんで転置されてるんだかよくわからない。

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

では縮約記法でを計算してみよう。

$$ \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

テンソルを使いこなしたい２（レヴィ=チヴィタ記号）

hotoke-x.hatenablog.com

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

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

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

添字がの偶置換なら、奇置換なら、それ以外は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つの直交した単位ベクトル（）を使えば、スカラー三重積でも表現できる。

$$ \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} $$

証明は具体的にに数字を入れて計算してみれば良い。例えばの時、考えられるのパターンは

$$ \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}はとなる項を省略すると

$$ \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} $$

最後の式変形はであることを利用した。

その他

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

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} $$ の証明をしている。の計算の参考になる。
3. ベクトル解析公式の証明 – 準備篇 | 高校物理の備忘録
   証明がたくさん書いてある。丁寧に途中式まである。

まとめ

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

*1:Wikipediaより

テンソルを使いこなしたい１（アインシュタインの縮約記法）

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

テンソルとは

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

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

のような対応がある。

Einsteinの縮約記法

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

例）

　　　 ：0階テンソル
　　：1階テンソル
：2階テンソル

とおくと、

• 内積（0階テンソル）

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

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

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

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

• 行列とベクトルの積（1階テンソル）

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

• 行列積（2階テンソル）

$$ \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} $$

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

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

テンソル積（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} $$

なお、, は基底を表す。

まとめ

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

*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階微分）=０の式を作ると、二次の項でヘッセ行列xベクトルの形が出てくる。このベクトルが共役勾配である。 共役勾配は最適化したい関数の1階微分から得られる勾配とは、少しずれたベクトルになる。 ずれ分のベクトルを導入し、 共役勾配=勾配+αずれベクトル と線形和で表現し、αを求めることでヘッセ行列の逆行列を計算することを回避した。

最小二乗法

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

非線形最小二乗法

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

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

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

統計的最適化

• 最尤推定

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

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

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

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

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

線形計画法

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

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

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

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

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

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

非線形計画法

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

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

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

ラグランジュ乗数

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

双対原理

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

動的計画法

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

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

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

参考書籍

www.amazon.co.jp

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

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

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

hotoke-x.hatenablog.com

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

正方行列のスペクトル分解（固有値分解）

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

$$ \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} $$

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

$$ \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} $$

を考える。次に、ベクトルを固有ベクトルを基底とする空間に投げる。

$$ \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次形式を書き直しただけである。

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

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

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

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

これなら分かる最適化数学（金谷健一）、例題1.25、1.27よりの時

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

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

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

$$ \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は固有値になっている。

また、とは同値である。

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

となり、等高線が楕円になっていることがわかる。に制限すれば、のとき最小値2（最小固有値）、のとき最大値7（最大固有値）になることがわかる。

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

参考書籍

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

固有値と固有ベクトル

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

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

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

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

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

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

が成り立つでないベクトルを行列の固有ベクトル、をその固有値という。
が対称行列なら、は個の固有値を持つ。また、の各成分間は互いに直交する。

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

この式を固有方程式（特性方程式）と呼ぶ。

つまりなんなのか

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

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

基底変換

今、デカルト座標系にベクトル $$ \begin{align} \boldsymbol{v} = \left( \begin{array}{c} x \\ y \end{array} \right) \end{align} $$ が存在して、座標変換によって新たな表現を得ることを考える。
これはデカルト座標での表現なので、基底ベクトル, を使って大げさに書けば $$ \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} $$ この式の正方行列は、ベクトルの変換意味しているのだった（ベクトル空間の線形写像）。
説明のために、 $$ \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} $$ と書けることから、行列の列成分は変換後の空間での基底ベクトルを表していたことがわかる。 つまり、という数字をそのまま行列による変換後の空間に渡すとどうなるか計算していただけのことだった。

正方行列の対角化

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

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

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

正方行列によってn次元ベクトルが別のn次元ベクトルに変換される、 $$ \boldsymbol{t} = \boldsymbol{A s} $$ という変換を考える。ここに、行列の右側に固有ベクトルを並べた行列 $$ \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)} $$ となる。の部分は、一番最初に登場した式 $$ \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)} $$ の式をよく見てみると、とする変換だったのに、のように見えてしまう。を使いたいために格好が悪くなってしまった。

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

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

となる（有名な式らしい...知らんけど）。

なんだか大事になっているように見えるが、実はそうではない。の部分が実は、固有値を対角成分に並べた対角行列になるからだ（を対角行列にする操作を対角化と呼ぶ）。

$$ \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} $$

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

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

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

$$ \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} $$

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