多変量正規分布を適用した場合のKullback Leibler Divergence(KL Divergence; カルバックライブラーダイバージェンス)の導出を丁寧に書いてみようと思います。
つまり、

\begin{aligned}
&D_{KL}(p(x) \| q(x)) =\\
&{1 \over 2}\left[ \log {|\Sigma_1| \over |\Sigma_0|} - d + {\rm Tr} (\Sigma_1^{-1}\Sigma_0 ) + 
(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0)^T \Sigma_1^{-1} (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0)\right],
\end{aligned}

\begin{aligned}
p(x) \sim N(\boldsymbol{\mu}_0, \Sigma_0), \\
q(x) \sim N(\boldsymbol{\mu}_1, \Sigma_1)
\end{aligned}

を導出したいと思います。

1.道具の準備

1-1. 行列の積のトレース

行列の積を成分表示すると

(AB)_{ij} = \left\{ \sum_{k=1}^b a_{ik} b_{kj}  \right\} \quad \cdots {\rm Eq.} (0)

なので、行列の積のトレースは

{\rm Tr}(AB) = \sum_i\sum_j a_{ij} b_{ji} \quad \cdots {\rm Eq.} (1)

と書けます。

1-2. 確率変数の積の期待値

共分散は分散共分散行列$\Sigma$の$i,j$成分なので、

\begin{aligned}
\{\Sigma\}_{ij} &= \sigma_{ij} \\
&= \mathbf{E}[(x_i-\mu_i)(x_j-\mu_j)] \\
&= \mathbf{E}[x_ix_j -\mu_ix_j - \mu_jx_i + \mu_i\mu_j] \\
&= \mathbf{E}[x_ix_j] -\mu_i\mu_j - \mu_j\mu_i + \mu_i\mu_j \\
&= \mathbf{E}[x_ix_j] -\mu_i\mu_j
\end{aligned}

よって、確率変数の積の期待値は

\mathbf{E}[x_ix_j] = \sigma_{ij} + \mu_i\mu_j \quad \cdots {\rm Eq.} (2)

のように書けます。

1-3. 2次形式

2次形式は、

\boldsymbol{x}^TA\boldsymbol{x} = \sum_i \sum_j a_{ij}x_i x_j

と書き下せます。

1-4. 確率変数の２次形式の期待値

次に、確率変数の２次形式に関する変形。

\begin{aligned}
\mathbf{E}[\boldsymbol{x}] &= \boldsymbol{\mu}, \\
\mathbf{E}[(\boldsymbol{x} - \boldsymbol{\mu})(\boldsymbol{x} - \boldsymbol{\mu})^T ] &=  \Sigma
\end{aligned}

とする。

\begin{aligned}
&\mathbf{E}[\boldsymbol{x}^TA\boldsymbol{x}] \\
&=\mathbf{E}[\sum_i \sum_j a_{ij} x_i x_j] \\
&=\sum_i \sum_j a_{ij} \mathbf{E}[x_i x_j] \\
&=\sum_i \sum_j a_{ij} (\sigma_{ij} + \mu_i\mu_j) \quad (\leftarrow using\ Eq.(2)) \\
&=\sum_i \sum_j a_{ij}\sigma_{ji} +  \sum_i \sum_j a_{ij}\mu_i\mu_j \\
&={\rm Tr}(A\Sigma) + \boldsymbol{\mu}^TA\boldsymbol{\mu} \quad (\leftarrow using\ Eq.(1)\ \ for\ 1_{st}\ term)  \quad \cdots {\rm Eq.} (3)
\end{aligned}

(5行目の左の項は$\Sigma$が対象行列なので$\sigma_{ij}=\sigma_{ji}$を適用しています。)

1-5. 中心化した確率変数の、分散共分散行列の逆行列に関する２次形式の期待値

(タイトルが長いな・・・)　つまり、$\mathbf{E}[(\boldsymbol{x}-\boldsymbol{\mu})^T\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu})]$です。

\begin{aligned}
\mathbf{E}&[(\boldsymbol{x}-\boldsymbol{\mu})^T\Sigma^{-1}(\boldsymbol{x}-\boldsymbol{\mu})] \\
&={\rm Tr}(\Sigma^{-1}\Sigma) + \boldsymbol{0}^TA\boldsymbol{0}  \quad (\leftarrow using\ Eq.(3)) \\
&= {\rm Tr}(I_d) = d   \quad \cdots {\rm Eq.} (4)
\end{aligned}

1-6. 異なる分布の平均値と逆行列を適用した場合の1-3.の２次形式の期待値

2つの分布、

p(x) \sim N(\mu_0, \Sigma_0), \\
q(x) \sim N(\mu_1, \Sigma_1)

があるとします。
また、$\mathbf{E}_{p(x)}[\cdot]$を分布$p(x)$に関する期待値とします。
$p(x)$に関する期待値をとるのですが、異なる分布$q(x)$の期待値で中心化を行い、またこの$q(x)$の分散共分散行列の逆行列を使って2次形式をとったもののことです。（下の式を見たほうが早そうです）

また、逆行列の成分表示として

\left( \Sigma_1^{-1} \right)_{ij} = s_{1, ij}

という記号を使います。

\begin{aligned}
&\mathbf{E}_{p(x)}[(\boldsymbol{x}-\boldsymbol{\mu}_1)^T\Sigma_1^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_1)] \\
&= \mathbf{E}_{p(x)}[\sum_i \sum_j s_{1, ij}(x_i-\mu_{1,i})(x_j-\mu_{1,j})] \\
&= \mathbf{E}_{p(x)}[\sum_i \sum_j s_{1, ij}(x_ix_j - x_i\mu_{1,j} - x_j\mu_{1,i} + \mu_{1,i}\mu_{1,j}  )] \\
&= \sum_i \sum_j s_{1, ij} (\mathbf{E}_{p(x)}[x_ix_j] - \mathbf{E}_{p(x)}[x_i]\mu_{1,j} - \mathbf{E}_{p(x)}[x_j]\mu_{1,i} + \mu_{1,i}\mu_{1,j} )] \\
&= \sum_i \sum_j s_{1, ij} ((\sigma_{ij} + \mu_{0,i}\mu_{0,j}) - \mu_{0,i}\mu_{1,j} - \mu_{0,j}\mu_{1,i} + \mu_{1,i}\mu_{1,j} )   \quad (\leftarrow using\ Eq.(2)) \\
&= \sum_i \sum_j s_{1, ij}\sigma_{ij} + \sum_i \sum_j s_{1, ij}(\mu_{0,i}\mu_{0,j} - \mu_{0,i}\mu_{1,j} - \mu_{0,j}\mu_{1,i} + \mu_{1,i}\mu_{1,j} )\\
&= \sum_i \sum_j s_{1, ij}\sigma_{ij} + \sum_i \sum_j s_{1, ij}
(\mu_{0,i}-\mu_{1,i} ) (\mu_{0,j} -\mu_{1,j}) \\
&= {\rm Tr}(\Sigma_0\Sigma_1^{-1}) + 
(\boldsymbol{\mu}_0 -\boldsymbol{\mu}_1)^T\Sigma_1^{-1}(\boldsymbol{\mu}_0-\boldsymbol{\mu}_1)   \quad \cdots {\rm Eq.} (5) 
\quad (\leftarrow using\ Eq.(1)\ \ for\ 1_{st}\ term )
\end{aligned}

2.正規分布の場合のKullback Leibler Divergenceの導出

以上の道具を使って多変量正規分布を適用した場合のKullback Leibler Divergenceの導出を見ていきたいと思います。
2つの分布は下記の通りです。

p(x) \sim N(\mu_0, \Sigma_0), \\
q(x) \sim N(\mu_1, \Sigma_1)

このとき、Kullback Leibler Divergence $D_{KL}(p(x) || q(x)) $は、

\begin{aligned}
&D_{KL}(p(x) \| q(x)) \\
&= \int p(x) \log {p(x) \over q(x)} dx \\
&= \int p(x)[ \log p(x) - \log q(x)] dx \\
&= \int p(x)[ 
-{d \over 2}\log 2\pi -{ 1 \over 2} \log |\Sigma_0| -{ 1 \over 2} (\boldsymbol{x}-\boldsymbol{\mu}_0)^T\Sigma_0^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_0) \\
&\quad\quad\quad\quad +{d \over 2}\log 2\pi +{ 1 \over 2} \log |\Sigma_1| +{ 1 \over 2} (\boldsymbol{x}-\boldsymbol{\mu}_1)^T\Sigma_1^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_1) 
] dx \\
&= {1 \over 2}\log {|\Sigma_1| \over |\Sigma_0|} -{1 \over 2} 
\mathbf{E}_{p(x)}[(\boldsymbol{x}-\boldsymbol{\mu}_0)^T\Sigma_0^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_0) ]
+{1 \over 2} 
\mathbf{E}_{p(x)}[(\boldsymbol{x}-\boldsymbol{\mu}_1)^T\Sigma_!^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_1) ]
\\
&= {1 \over 2}\log {|\Sigma_1| \over |\Sigma_0|} -{1 \over 2} d
+{1 \over 2} 
\mathbf{E}_{p(x)}[(\boldsymbol{x}-\boldsymbol{\mu}_1)^T\Sigma_1^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_1) ]    \quad (\leftarrow using\ Eq.(4)\ \ for\ 2_{nd}\ term) 
\\
&= {1 \over 2}\log {|\Sigma_1| \over |\Sigma_0|} -{1 \over 2} d
+{1 \over 2} 
{\rm Tr}(\Sigma_0\Sigma_1^{-1}) + 
(\boldsymbol{\mu}_0 -\boldsymbol{\mu}_1)^T\Sigma_1^{-1}(\boldsymbol{\mu}_0-\boldsymbol{\mu}_1)
\quad (\leftarrow using\ Eq.(5)\ \ for\ 3_{rd}\ term) 
\\
&= {1 \over 2}\left[ \log {|\Sigma_1| \over |\Sigma_0|} - d + {\rm Tr} (\Sigma_0 \Sigma_1^{-1}) + 
(\boldsymbol{\mu}_0 - \boldsymbol{\mu}_1)^T \Sigma_1^{-1} (\boldsymbol{\mu}_0 - \boldsymbol{\mu}_1)\right]\\
&= {1 \over 2}\left[ \log {|\Sigma_1| \over |\Sigma_0|} - d + {\rm Tr} (\Sigma_1^{-1}\Sigma_0 ) + 
(\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0)^T \Sigma_1^{-1} (\boldsymbol{\mu}_1 - \boldsymbol{\mu}_0)\right]
\end{aligned}

できました！