こんしゅーの統計memo

統計学のmemoをしていきます

統計で使う行列の微分

統計で使いそうな行列の性質と公式を記載しておきます． 詳しくは統計のための行列代数(上,下)をご参照． 全体を通して実数の範囲で考えます．

ベクトルによる微分

ベクトル \boldsymbol{x}=(x _ 1,...,x _ N)^\top による微分は，微分の対象となる関数 f(\boldsymbol{x}) を \boldsymbol{x} の各要素 x _ i で微分したものを並べたベクトル

\displaystyle{ \dfrac{\mathrm{d}f(\boldsymbol{x})}{\mathrm{d}\boldsymbol{x}} =\left(\dfrac{\mathrm{d}f(\boldsymbol{x})}{\mathrm{d}x _ 1},...,\dfrac{\mathrm{d}f(\boldsymbol{x})}{\mathrm{d}x _ N}\right)^\top }

と定義します． ここで，関数 f(\boldsymbol{x}) はスカラー値であることに注意します． この定義に従って，以下の典型的なベクトルによる微分の公式を導出します．

内積に対する微分

係数ベクトル \boldsymbol{a}=(a _ 1,...,a _ N)^\top と \boldsymbol{x}=(x _ 1,...,x _ N)^\top の内積 \boldsymbol{a}^\top\boldsymbol{x} の \boldsymbol{x} による微分は，各 x_i の微分が

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{a}^\top\boldsymbol{x}}{\mathrm{d}x_i}=\dfrac{\mathrm{d}}{\mathrm{d}x _ i}\sum _ {j=1}^Na _ jx _ j=a _ i }

となることから，

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{a}^\top\boldsymbol{x}}{\mathrm{d}\boldsymbol{x}} =\boldsymbol{a} }

となります． 微分結果の係数ベクトル \boldsymbol{a} に転置がかかっていることに注意します．

二次形式に対する微分

N\times N の係数行列

\displaystyle{ \boldsymbol{A} =\begin{pmatrix}A _ {11} & \cdots & A _ {1N}\\\vdots & \ddots & \vdots\\A _ {N1} & \cdots & A _ {NN}\\\end{pmatrix} }

と \boldsymbol{x}=(x _ 1,...,x _ N)^\top の二次形式 \boldsymbol{x}^\top\boldsymbol{Ax} の \boldsymbol{x} による微分を求めます． 行列 \boldsymbol{A} による微分ではなくベクトル \boldsymbol{x} による微分であることに注意します． 内積の微分と同様に各 x_i の微分は

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{x}^\top\boldsymbol{Ax}}{\mathrm{d}x _ i} =\dfrac{\mathrm{d}}{\mathrm{d}x _ i}\sum _ {j,k=1}^Nx _ jA _ {jk}x _ k =\sum _ {j=1}^Nx _ jA _ {ji}+\sum _ {k=1}^NA _ {ik}x _ k }

となるので，これを列ベクトルとして並べると

\displaystyle{ \begin{pmatrix} \sum _ {j=1}^NA _ {j1}x _ j+\sum _ {k=1}^NA _ {1k}x _ k\\ \vdots\\ \sum _ {j=1}^NA _ {jN}x _ j+\sum _ {k=1}^NA _ {Nk}x _ k\\ \end{pmatrix} = \boldsymbol{Ax}+\boldsymbol{A}^\top\boldsymbol{x} =(\boldsymbol{A}+\boldsymbol{A}^\top)\boldsymbol{x} }

となります．よって，二次形式の微分は

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{x}^\top\boldsymbol{Ax}}{\mathrm{d}\boldsymbol{x}} =(\boldsymbol{A}+\boldsymbol{A}^\top)\boldsymbol{x} }

となることがわかります．

加えて，行列 \boldsymbol{A} が対称な行列 A _ {ij} = A _ {ji} ならば \boldsymbol{A}=\boldsymbol{A}^\top なので，

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{x}^\top\boldsymbol{Ax}}{\mathrm{d}\boldsymbol{x}} =2\boldsymbol{Ax} }

となります．

行列による微分

ベクトルによる微分と同様に， M\times N 行列

\displaystyle{ \boldsymbol{X} =\begin{pmatrix} X _ {11} & \cdots & X _ {1N} \\\vdots & \ddots & \vdots\\ X _ {M1} & \cdots & X _ {MN} \\ \end{pmatrix} }

による微分も，微分の対象となる関数 f(\boldsymbol{X}) を \boldsymbol{X} の各要素 X _ {ij} で微分したものを並べた行列

\displaystyle{ \dfrac{\mathrm{d}f(\boldsymbol{X})}{\mathrm{d}\boldsymbol{X}} =\begin{pmatrix} \dfrac{\mathrm{d}f(\boldsymbol{X})}{\mathrm{d}X _ {11}} & \cdots & \dfrac{\mathrm{d}f(\boldsymbol{X})}{\mathrm{d}X _ {1N}} \\\vdots & \ddots & \vdots\\ \dfrac{\mathrm{d}f(\boldsymbol{X})}{\mathrm{d}X _ {M1}} & \cdots & \dfrac{\mathrm{d}f(\boldsymbol{X})}{\mathrm{d}X _ {MN}} \\ \end{pmatrix} }

と定義します． ここでも，関数 f(\boldsymbol{X}) はスカラー値であることに注意します． この定義に従って，以下の典型的な行列による微分の公式を導出します．

行列の Trace

まず，行列積のトレースを，行列の要素の和で表しておきます． \boldsymbol{A} を N\times M 行列，\boldsymbol{X} を M \times L 行列，\boldsymbol{B} を L \times K 行列とします． 行列 \boldsymbol{A} の i,j 要素を [ \boldsymbol{A} ] _ {ij} と表すと，行列 \boldsymbol{A} のトレース \mathrm{tr}(\boldsymbol{A}) は

\displaystyle{ \mathrm{tr}(\boldsymbol{A}) = \sum _ i [\boldsymbol{A}] _ {ii} = \sum _ i A _ {ii} }

で定義されます．

行列積 \boldsymbol{AX} を和で表すと

\displaystyle{ [\boldsymbol{AX}] _ {ij} = \sum _ {k=1}A _ {ik}X _ {kj} }

と表されます．したがって，

\displaystyle{ \mathrm{tr}(\boldsymbol{AX}) =\sum _ i[\boldsymbol{AX}] _ {ii} =\sum _ i\sum _ kA _ {ik}X _ {ki} }

となります．

同様に行列積 \boldsymbol{AXB} は

\displaystyle{ [\boldsymbol{AXB}] _ {ij} = \sum _ k\sum _ lA _ {ik}X _ {kl}B _ {lj} }

と表されるので，

\displaystyle{ \mathrm{tr}(\boldsymbol{AXB}) =\sum _ i[\boldsymbol{AXB}] _ {ii} =\sum _ i\sum _ k\sum _ lA _ {ik}X _ {kl}B _ {li} }

となります．

Trace の微分

行列積 \boldsymbol{AX} のトレース \mathrm{tr}(\boldsymbol{AX}) を X _ {mh} で微分すると， X _ {mh} の係数を抜き出せばよいから，

\displaystyle{ \dfrac{\mathrm{d}}{\mathrm{d}X _ {mh}}\mathrm{tr}(\boldsymbol{AX}) =\dfrac{\mathrm{d}}{\mathrm{d}X _ {mh}}\sum _ i\sum _ kA _ {ik}X _ {ki} =A _ {hm} }

となるので，

\displaystyle{ \dfrac{\mathrm{d}\mathrm{tr}(\boldsymbol{AX})}{\mathrm{d}\boldsymbol{X}} =\boldsymbol{A}^\top }

となります． 同様に，\boldsymbol{XB} のトレースの微分から

\displaystyle{ \dfrac{\mathrm{d}}{\mathrm{d}X _ {mh}}\mathrm{tr}(\boldsymbol{XB}) =\dfrac{\mathrm{d}}{\mathrm{d}X _ {mh}}\sum _ i\sum _ kX _ {ik}B _ {ki} =B _ {hm} }

となるので，

\displaystyle{ \dfrac{\mathrm{d}\mathrm{tr}(\boldsymbol{XB})}{\mathrm{d}\boldsymbol{X}} =\boldsymbol{B}^\top }

がわかります．

また，行列積 \boldsymbol{AXB} のトレース \mathrm{tr}(\boldsymbol{AXB}) を X _ {mh} で微分すると

\displaystyle{ \dfrac{\mathrm{d}}{\mathrm{d}X _ {mh}}\mathrm{tr}(\boldsymbol{AXB}) =\dfrac{\mathrm{d}}{\mathrm{d}X _ {mh}}\sum _ i\sum _ k\sum _ lA _ {ik}X _ {kl}B _ {li} =\sum _ iA _ {im}B _ {hi} =[\boldsymbol{BA}] _ {hm} }

と計算できるので，

\displaystyle{ \dfrac{\mathrm{d}\mathrm{tr}(\boldsymbol{AXB})}{\mathrm{d}\boldsymbol{X}} =(\boldsymbol{BA})^\top=\boldsymbol{A}^\top\boldsymbol{B}^\top }

がわかります．

二次形式の微分

ベクトルによる微分と表記を変え，N\times N の係数行列

\displaystyle{ \boldsymbol{X} =\begin{pmatrix}X _ {11} & \cdots & X _ {1N}\\\vdots & \ddots & \vdots\\X _ {N1} & \cdots & X _ {NN}\\\end{pmatrix} }

と \boldsymbol{a}=(a _ 1,...,a _ N)^\top の二次形式 \boldsymbol{a}^\top\boldsymbol{Xa} の \boldsymbol{X} による微分を求めます． X _ {mh} で微分すると，

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{a}^\top\boldsymbol{Xa}}{\mathrm{d}X _ {mh}} =\dfrac{\mathrm{d}}{\mathrm{d}x _ i}\sum _ {j,k=1}^Na _ jX _ {jk}a _ k =a _ ma _ h }

となります． したがって，これらの要素を行列に表すと

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{a}^\top\boldsymbol{Xa}}{\mathrm{d}\boldsymbol{X}} =\begin{pmatrix} a _ 1a _ 1 & \cdots & a _ 1a _ N \\ \vdots & \ddots & \vdots\\ a _ Na _ 1 & \cdots & a _ Na _ N \\ \end{pmatrix} }

となるので，

\displaystyle{ \dfrac{\mathrm{d}\boldsymbol{a}^\top\boldsymbol{Xa}}{\mathrm{d}\boldsymbol{X}} =\boldsymbol{aa}^\top }

となることがわかります．

カーネル法の基礎

この投稿では，私が後で見返したときに簡単にカーネル法を思い出せるように，その基礎を数学的な詳細に余りこだわらずにmemoします．

カーネル法の考え方

まず線形予測問題を考えます．すなわち， 入力する変数を表すベクトル \boldsymbol{x} と基底関数によるベクトル

\phi(\boldsymbol{x})=\left(\phi_1(\boldsymbol{x}),...,\phi_D(\boldsymbol{x})\right)^\top\in\mathbb{R}^ D

を使った，パラメータベクトル \boldsymbol{\beta}\in\mathbb{R}^D の線形モデル

f(\boldsymbol{x})=\boldsymbol{\beta}^\top\phi(\boldsymbol{x})

で出力となる変数 y を予測します．

データ (\boldsymbol{x} _ i,y _ i),...,(\boldsymbol{x} _ n,y _ n) が得られた時に，y_i を f(\boldsymbol{x} _ i) で近似する事を考えます． 最小二乗法で推定することにすると，二乗誤差

E(\beta)=\displaystyle{\sum_{i=1}^ n\left(y _ i-\beta^\top \phi(\boldsymbol{x} _ i)\right)^ 2}

を最小にする \betaを求めることになります． 変数のデータを表す n\times D 行列を \boldsymbol{X} と，出力のデータを表すベクトルを \boldsymbol{Y} と表すと，E を最小にする \hat{\beta} は

\hat{\beta}=\left(\boldsymbol{X}^ \top\boldsymbol{X}\right)^ {-1}\boldsymbol{X}^ \top\boldsymbol{Y}

と求めることができます．

ここで，視点を変えて，\hat{\beta} の別の表現方法を考えます． そのために，最小二乗法の目的関数 E(\beta) における \beta を次のように分類します． まず， \phi(\boldsymbol{x} _ 1),...,\phi(\boldsymbol{x} _ D)^\top で張られる \mathbb{R}^ D の部分空間を S と表します．すなわち，\beta _ S が S に含まれていれば，それは n 個のスカラー値の係数 \alpha _ 1,...,\alpha _ n を使って，

\displaystyle{\beta _ S=\sum _ {i=1}^ n\alpha _ i\phi(\boldsymbol{x} _ i)}

と表されることを意味します．そして \mathbb{R}^ Dの直交補空間を S^ \bot と表すと，S^ \bot の要素 \beta _ {S^ \bot} を使って，\beta は

\beta=\beta _ S+\beta _ {S^ \bot}

と表され，\beta _ {S^ \bot} は S の直交補空間に含まれることから，任意の i=1,...,n について，\beta _ {S^ \bot}^ \top\phi(\boldsymbol{x} _ i)=0 をみたします．

これらを使って，最小二乗法の目的関数 E(\beta) を書き換えると，次のようになります．

\displaystyle{E(\beta)=\sum _ {i=1}^ n\left(y _ i-\beta _ S^ \top\phi(\boldsymbol{x} _ i)\right)^ 2+\sum _ {i=1}^ ny _ i^ 2. }

二項目は \beta _ {S^ \bot}^ \top\phi(\boldsymbol{x} _ i)=0 から得られる項ですが，最小化に寄与しない項となるので，最適化の候補となるパラメータは \phi(\boldsymbol{x} _ i) の線形和となるパラメータ \beta _ S だけを考えれば十分であるということがわかりました． さらに関数 k(\boldsymbol{x} _ i,\boldsymbol{x} _ j) を

\displaystyle{k(\boldsymbol{x} _ i,\boldsymbol{x} _ j)=\phi(\boldsymbol{x} _ i)^ \top\phi(\boldsymbol{x} _ j) }

とおくと，\beta _ S が \phi(\boldsymbol{x} _ i) の線形和で表されることと合わせて，最小二乗法の目的関数を

\displaystyle{E(\beta)=\sum _ {i=1}^ n\left(y _ i-\sum _ {j=1}^ n\alpha _ j\phi(\boldsymbol{x} _ j)^ \top\phi(\boldsymbol{x} _ i)\right)^ 2=\sum _ {i=1}^ n\left(y _ i-\sum _ {j=1}^ n\alpha _ jk(\boldsymbol{x} _ j,\boldsymbol{x} _ i)\right)^ 2}

と書き換えることができます． この係数 \alpha _ i を並べたベクトル \alpha と表したとき，E(\beta) を最小にする \hat{\alpha} は，i,j 要素に K _ {ij} = k(\boldsymbol{x} _ i,\boldsymbol{x} _ j) を持つ n\times n 行列を K をとおくと，

\hat{\alpha}=K^ {-1}Y
と表すことができ，この結果から，回帰関数の推定量を

\displaystyle{f(\boldsymbol{x})=\sum _ {i=1}^ n\hat{\alpha} _ ik(\boldsymbol{x} _ i,\boldsymbol{x})}

と表現することができます．

こうすると，いままで線形予測問題で使っていた基底関数 \phiを考える代わりに，関数 k(\boldsymbol{x},\boldsymbol{x}') を考えれば十分であることがわかりました．また，k(\boldsymbol{x},\boldsymbol{x}') を決めて，データから行列 K を求めたとき，\hat{\alpha} を求めるのに必要な計算量は \phi の次元 D ではなく，データのサンプルサイズ n に依存する事がわかります． これから逆に，「先に関数 k(\boldsymbol{x},\boldsymbol{x}') を決めたときに，いままでと同じ様に予測問題を考えることができるだろうか」という疑問が湧いてきます．この疑問に答えるのが「Representer 定理」と呼ばれるもので，その道具立てとして「再生核Hilbert空間」というものが必要になってきます．

再生核Hilbert空間とRepresenter 定理

以降ではデータ点を x _ i，x _ j などのように表します． 上記で定めた k(x _ i,x _ j) は正定値カーネルと呼ばれます． 元々，k(x _ i,x _ j)=\phi(x _ i)^ \top\phi(x _ j) と定めていたことに対応させて，正定値カーネルが次の2つの条件をみたすものと定めます．

  • (対称性):k(x _ i,x _ j)=k(x _ j,x _ i)
  • (正値性):任意の c _ 1,...,c _ n に対して，\sum _ {i,j=1}^ nc _ ic _ jk(x _ i, x _ j) \ge0

この性質から，上で定めたような n\times n 行列 K=(K _ {ij})=(k(x _ i,x _ j)) は半正定値な対称行列となります．このような行列はGram行列と呼ばれます．

この正定値カーネルはいくつか性質を持ちますが，特に重要なものは以下のつです

先程の議論で回帰関数の推定量が

\displaystyle{f(\boldsymbol{x})=\sum _ {i=1}^ n\hat{\alpha} _ ik(\boldsymbol{x} _ i,\boldsymbol{x})}

で得られていたことを踏まえて，この正定値カーネルの線形和で表される関数がどれだけ表現力があるか考えます． まずはその集合 \mathcal{H} を以下の様に定めます:

\displaystyle{\mathcal{H}=\left\{f(x)=\sum _ {i=1}^ n\alpha _ ik(z _ i,x)\right\} }

ここで，n は任意の自然数で， z _ i は任意の入力データ変数で観測された変数でなくても構いません． 異なる \mathcal{H} の要素 f,g の距離を考えるために，それらの内積を定めます．そのため，正定値カーネル同士の内積を

\displaystyle{\langle k(x _ i,\cdot),k(x _ j,\cdot)\rangle=k(x _ i,x _ j) }

と定めると，\mathcal{H} の要素

\displaystyle{f(\cdot)=\sum\alpha_ik(x _ i,\cdot)}，\displaystyle{g(\cdot)=\sum\beta_jk(x _ j,\cdot) }

の内積を

\displaystyle{\langle f,g\rangle=\left\langle\sum _ i\alpha _ ik(x _ i,\cdot),\sum _ i\beta _ jk(x _ j,\cdot)\right\rangle=\sum _ {i,j}\alpha _ i\beta _ j\langle k(x _ i, \cdot),k(x _ j, \cdot)\rangle=\sum _ {i,j}\alpha _ i\beta _ jk(x _ i, x _ j) }

と計算することができます．したがって

\displaystyle{\langle f,g\rangle=\sum _ {i,j}\alpha _ i\beta _ jk(x _ i, x _ j)=\sum _ j\beta _ jf(x _ j)=\sum _ i\alpha _ ig(x _ i) }

となることから，\langle f,g\rangle を片方の関数値から計算することができ，線形和の表し方に依存しないことがわかります． また上式で g(\cdot)=k(x,\cdot) とおくことで，

\displaystyle{\langle f,k(x,\cdot)\rangle=\sum _ i\alpha _ i\langle k(x_i,\cdot),k(x,\cdot)\rangle=\sum _ i\alpha _ ik(x_i,x)=f(x) }

と計算することができ，この性質はカーネルの再生性と呼ばれます． 実は，再生性を持つ正定値カーネルをもつ関数の空間を考えることは，正定値カーネルの線形性から \mathcal{H} を考えることと同等になります． このことから，次に定める再生核Hilbert空間の性質を調べることを目的にすればよいことがわかりました．

再生核Hilbert空間 集合 \mathcal{X} 上の関数からなるHilbert空間を \left(\mathcal{H},\langle\cdot,\cdot\rangle _ {\mathcal{H}}\right) と表します． 関数 k:\mathcal{X}^ 2\rightarrow\mathbb{R} が存在して，任意の x\in\mathcal{X} と f\in\mathcal{H} に対して

\displaystyle{k(x,\cdot)\in\mathcal{H}}，\displaystyle{\langle f,k(x,\cdot)\rangle _ {\mathcal{H}}=f(x) }

が成り立つとき，\mathcal{H} を集合 \mathcal{X} 上の再生核Hilbert空間と呼び，関数 k を再生核と呼びます．

このように定めた再生核Hilbert空間の再生核は正定値カーネルとなり，また再生核は一意になります． このことから，再生核Hilbert空間を定めることから議論を始めることは，正定値カーネルから始めることと同等になることがわかります．

再生核Hilbert空間の性質はいくつかあるのですが，次のRepresenter定理が統計モデルの推定に関連するものです．

Representer定理 データ点として (x _ 1,y _ 1),...,(x _ n,y _ n) が得られたとき，損失関数 L と正則化項 \Psi の和

\displaystyle{ L\left(\left\{f(x_i)+b \right\}\right)+\Psi(||f|| _ {\mathcal{H}}) }

を最小化する f について，係数 \alpha _ i が存在して

\displaystyle{ f(x)=\sum _ {i=1}^ n\alpha _ ik(x _ i,x) }

と表すことができる．

参考文献