行列式の微分と近似

<script type="text/x-mathjax-config">
MathJax.Hub.Config({
  TeX: { 
  	equationNumbers: { 
  		autoNumber: "AMS"
  	}
  }
});
</script>
$$
\newcommand{\b}[1]{\mathbf{#1}}
$$
行列式の微分が必要になることってたまにありますよね、ってことで行列式の微分とテイラー展開を求めてみましょう。

行列式の微分に機械学習系の勉強をしている時に遭遇するのはおそらく多次元正規分布のパラメーターの推定（最尤推定）が一番初めだと思います。昔この問題にぶち当たった時 PRML の2章の部分で自分は[アンチョコ](http://prml.in/)を見ながらフムフムと計算していたのですが、この本の行列式の導出方法は洗練されていて思いつきにくく若干遠回りだと感じました。

もう少しいい方法がないかと探していたら [Terence Tao](https://terrytao.wordpress.com/2013/01/13/matrix-identities-as-derivatives-of-determinant-identities/) さんのブログで説明されてました。
なのでここで説明されているやり方に沿って行列式の微分を求めてみましょう。

## 行列式の微分

微分の定義通り計算する。

$$
\begin{eqnarray}
\frac{d}{dx} det \b{A}(x) & = & \lim_{\delta x \to 0} \frac{det(\b{A}(x+\delta x)) - det(\b{A}(x))}{\delta x} \\\
& = & \lim_{\delta x \to 0} \frac{det(\b{A}(x) + \frac{\partial \b{A}}{\partial x} \delta x + O((\delta x)^2)) - det(\b{A}(x))}{\delta x} \\\
& = & \lim_{\delta x \to 0} \frac{det(\b{A}(x) + \frac{\partial \b{A}}{\partial x} \delta x) - det(\b{A}(x)) + O((\delta x)^2)}{\delta x} \\\
& = & \lim_{\delta x \to 0} \frac{det(\b{A}(x)) (det(\b{I} + \b{A}^{-1}\frac{\partial \b{A}}{\partial x} \delta x) - 1) + O((\delta x)^2)}{\delta x} 
\end{eqnarray}
$$

と式変形が出来る。ここで問題になるのが $ det (\b{I} + \b{A}^{-1}\frac{\partial \b{A}}{\partial x} \delta x) $ の項である。
$ \b{N} = \b{A}^{-1}\frac{\partial \b{A}}{\partial x} $ とおけば求めるべきは $ det(\b{I} + \b{N}\delta x) $ である。ここで行列の中身がどのようになっているのかを具体的に書き出すと

$$
 \b{I} + \b{N}\delta x = \left(
 		\begin{array}{c}
 			1 + N_{11}\delta x & N_{12} \delta x & \ldots & N_{1n} \delta x \\\
 			N_{21} \delta x & 1 + N_{22} \delta x & \ldots & N_{2n} \delta x \\\
 			\vdots & \vdots & \ddots & \vdots \\\
 			N_{n1} \delta x & N_{n2} \delta x & \ldots & 1 + N_{nn} \delta x
 		\end{array}
  \right)
$$

である。行列式は

$$
det(A) = \sum_{\sigma \in S_n} sgn(\sigma) A_{1\sigma(1)} A_{2\sigma(2)} \cdots A_{n\sigma(n)}
$$

で与えられることから$ \delta x $ が非常に小さい時、**行列の対角成分の項以外は行列式の値にほとんど影響を与えない**事がわかる。また置換群の性質から対角成分を除いた項を全て $ (\delta x)^2 $ に比例する。
つまり具体的に書くと

$$
\begin{eqnarray}
det(\b{I} + \b{N} \delta x) & = & \prod_{i=1}^{n} (1 + N_{ii}\delta x) + O((\delta x)^2) \\\
& = & 1 + \sum_{i=1}^{n} N_{ii} \delta x + O((\delta x)^2) \\\
& = & 1 + tr(\b{N}) \delta x + O((\delta x)^2) \\\
\end{eqnarray}
$$

であることが分かる。これを先ほどの式に代入すると

$$
\begin{eqnarray}
\frac{d}{dx} det \b{A}(x) & = & \lim_{\delta x \to 0} \frac{det(\b{A}(x)) (1 + tr(\b{A}^{-1} \frac{\partial \b{A}(x)}{\partial x})\delta x + O((\delta x)^2) - 1) + O((\delta x)^2)}{\delta x} \\\
& = & \lim_{\delta x \to 0} \frac{det(\b{A}(x)) ~ tr(\b{A}^{-1} \frac{\partial \b{A}(x)}{\partial x}) \delta x +  O((\delta x)^2)}{\delta x} \\\
& = & det(\b{A}(x)) ~ tr(\b{A}^{-1}(x) \frac{\partial \b{A}(x)}{\partial x})
\end{eqnarray}
$$

として行列式の微分を求められた。


### 本当に成り立っているか検証

以下の行列で本当にこの公式が成り立っているか以下の行列で確認してみる。

$$
\begin{eqnarray}
\b{A}(x) = \left(
 		\begin{array}{c}
 			2 x & x^2 \\\
 			x^3 + x & \frac{1}{x} \\\
 		\end{array}
  \right)
\end{eqnarray}
$$

まず普通に行列式を計算してみる

$$
\begin{eqnarray}
det (\b{A}(x)) & = & 2x \cdot \frac{1}{x}  - x^2 (x^3 + x) \\\
& = & 2 - x^5 - x^3 \\\
\end{eqnarray}
$$

より

$$
\begin{eqnarray}
\frac{d}{dx} det (\b{A}(x)) = - 5 x^4 - 3 x^2
\end{eqnarray}
$$

一方で行列式の微分の公式を用いて計算すると

$$
\begin{eqnarray}
\b{A}^{-1}(x) & = & \frac{1}{2 - x^5 + x^3} \left(
 		\begin{array}{c}
 			\frac{1}{x} & -x^2 \\\
 			- x^3 - x & 2x \\\
 		\end{array}
  \right) \\\
\frac{\partial}{\partial x} \b{A}(x) & = & \left(
 		\begin{array}{c}
 			2 & 2x \\\
 			3x^2 + 1 & -\frac{1}{x^2} \\\
 		\end{array}
  \right)
\end{eqnarray}
$$

これらを掛け合わせることで以下の式を得る。

$$
\b{A}^{-1} \cdot \frac{\partial}{\partial x} \b{A}(x) = 
\frac{1}{2-x^5-x^3} \left(
 		\begin{array}{c}
 		 	\frac{2}{x} - 3 x^4 - x^2 & ... \\\
 			... & -2 x^4 -2 x^2 - \frac{2}{x} \\\
 		\end{array}
  \right)
$$

したがって

$$
\begin{eqnarray}
det(\b{A}(x)) ~ tr(\b{A}^{-1}\frac{\partial}{\partial x} \b{A}(x)) & = & -5x^4 - 3x^2
\end{eqnarray}
$$

これより具体的な行列でも公式が正しいことが確認できた。

## 行列式のテイラー展開

行列式の近似値を求める際に具体的な行列式の計算を用いて以下の1次近似を出した。

$$
\begin{eqnarray}
det(\b{I} + \b{N} \delta x) & = &  1 + tr(\b{N}) \delta x + O((\delta x)^2) \\\
\end{eqnarray}
$$

さて、上の行列のさらに高次の微小量はどのような値になっているのか気にならないだろうか？ということで3次までの微小量を求めてみました。

この計算を前回と同様に具体的な行列式を展開して求めるのはエレガントではないので、以下のように $ f(x) $ を定義し、そのテイラー展開を行うことで高次の微小量を求めることにする。

$$
\begin{eqnarray}
f(x) = det(\b{I} + \b{A} x)
\end{eqnarray}
$$

まず微分の定義に従って一階微分を求める。

$$
\begin{eqnarray}
f'(x) & = & \lim_{\delta x \to 0} \frac{det(\b{I} + (x + \delta x)\b{A}) - det(\b{I} + x\b{A})}{\delta x} \\\
& = & \lim_{\delta x \to 0} \frac{det(\b{I} + x\b{A})\\{det(\b{I} + (\b{I} + x \b{A})^{-1}\b{A} \delta x) - 1\\}}{\delta x} \\\
& = & \lim_{\delta x \to 0} \frac{det(\b{I} + x\b{A})\\{1 + tr((\b{I} + x\b{A})^{-1}\b{A}) \delta x - 1\\}}{\delta x} \\\
& = & det(\b{I} + x\b{A})tr(\b{F}^{-1}(x)) \\\
& = & f(x) ~ tr(\b{F}^{-1}(x))
\end{eqnarray}
$$

ここで $ \b{F}(x) $ を以下のように定めた。

$$
\b{F}(x) = \b{A}^{-1}(\b{I} + x\b{A})
$$

ここで具体的に $ x = 0 $ を上式に代入すると

$$
\begin{eqnarray}
f(0) = 0 \\\
\b{F}^{-1}(0) = \b{A}
\end{eqnarray}
$$

より $ f'(0) = tr(\b{A}) $  となる。
同様にして二階微分を求めると

$$
\begin{eqnarray}
f^{(2)}(x) & = & f'(x)tr(\b{F}^{-1}(x)) - f(x) tr(\b{F}^{-2}(x)) \\\
& = & f(x) (tr(\b{F}^{-1}(x))^2 - tr(\b{F}^{-2}(x)) \\\
\end{eqnarray}
$$

同じ様にして三階微分も

$$
\begin{eqnarray}
f^{(3)}(x) & = & f'(x) (tr(\b{F}^{-1})^2 - tr(\b{F}^{-2}) + \nonumber \\\
& & ~~~~ f(x)(-2 ~ tr(\b{F}^{-1})tr(\b{F}^{-2}) + 2 tr(\b{F}^{-3})) \\\
& = & f(x) (tr(\b{F}^{-1})^3 - 3tr(\b{F}^{-1})tr(\b{F}^{-2}) + 2tr(\b{F}^{-3}))
\end{eqnarray}
$$

これらの結果から

$$
\begin{eqnarray}
f^{(1)}(0) & = & tr(\b{A}) \\\
f^{(2)}(0) & = & tr(\b{A})^2 - tr(\b{A}^2) \\\
f^{(3)}(0) & = & tr(\b{A})^3 - 3tr(\b{A})tr(\b{A}^2) + 2tr(\b{A}^3) \\\
\end{eqnarray}
$$

この結果を用いいて $ f(x) $ のテイラー展開する。

$$
\begin{eqnarray}
f(x) & = & f(0) + f^{(1)}(0)x + \frac{1}{2!} f^{(2)} x^2 + \frac{1}{3!} f^{(3)} x^3 +  ... \\\
f(x) & = & 1 + tr(\b{A}) x \\\
& & ~~~~ + \frac{1}{2} (tr(\b{A})^2 - tr(\b{A}^2)) x^2 \\\
& & ~~~~~~~~~ + \frac{1}{6} (tr(\b{A})^3 - 3tr(\b{A})tr(\b{A}^2) + 2tr(\b{A}^3))x^3 + .... \\\
\end{eqnarray}
$$

$ x \to \delta x $ とすることで目的の高次の微小量が得られていることが分かる。

## 行列式のテイラー展開と固有値との関係
さて 4x4 の逆行列の係数に先ほど求めたテイラー展開の値が入っているのがwikipediaの[Inversion of 4×4 matrices](https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_4.C3.974_matrices) で確認できる。

$ det(\b{I} + x\b{A}) $  の展開と $ \b{A}^{-1} $ の値という一見別に見えるものが繋がっているのは面白いですね。実はこれはケイリーハミルトンの定理が成り立つこを考えれば自明なのですが...