Two Multivariate Normal Distributions: KL Divergence, Bhattacharyya Distance and Wasserstein Distance
The normal distribution is one of the most common continuous probability distributions. It is the maximum-entropy distribution given a fixed mean and covariance (see "Entropy" Untamed: From Entropy, the Maximum Entropy Principle, to the Maximum Entropy Model (Part 2)), and it can also be viewed as the second-order approximation of any continuous distribution, playing a role analogous to the linear approximation of a general function. From this perspective, the normal distribution is arguably the simplest continuous distribution there is. And precisely because it is so simple, many estimators associated with it admit closed-form analytic expressions.
This post mainly works out several measures between two multivariate normal distributions, including the KL divergence, the Bhattacharyya distance, and the Wasserstein distance (W distance), all of which have explicit analytic solutions.
The Normal Distribution
Let's briefly review some basics of the normal distribution. Note that this is only a review, not a proper introduction to the normal distribution.
Probability density
The normal distribution, also known as the Gaussian distribution, is a continuous probability distribution defined on $\mathbb{R}^n$, whose probability density function is
\begin{equation}p(\boldsymbol{x})=\frac{1}{\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma})}}\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right\}\end{equation}more
Here $\boldsymbol{x},\boldsymbol{\mu}\in\mathbb{R}^n$, with $\boldsymbol{\mu}$ being the mean vector (throughout this post, vectors are by default column vectors), and $\boldsymbol{\Sigma}\in\mathbb{R}^{n\times n}$ being the covariance matrix, which is required to be symmetric positive definite. As we can see, the normal distribution is uniquely determined by $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$, so it's not hard to imagine that all of its statistics are functions of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$. When $\boldsymbol{\mu}=\boldsymbol{0}, \boldsymbol{\Sigma}=\boldsymbol{I}$, the corresponding distribution is called the "standard normal distribution."
Basic properties
Generally speaking, the basic statistics are the mean and variance, which correspond to the two parameters of the normal distribution:
\begin{equation}\begin{aligned} \mathbb{E}_{\boldsymbol{x}}\left[\boldsymbol{x}\right]=&\int p(\boldsymbol{x}) \boldsymbol{x} dx=\boldsymbol{\mu}\\ \mathbb{E}_{\boldsymbol{x}}\left[(\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{\top}\right]=&\int p(\boldsymbol{x}) (\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{\top} dx=\boldsymbol{\Sigma}\\ \end{aligned}\end{equation}
From this we can also derive the second moment:
\begin{equation} \mathbb{E}_{\boldsymbol{x}}\left[\boldsymbol{x}\boldsymbol{x}^{\top}\right]=\boldsymbol{\mu}\boldsymbol{\mu}^{\top} + \mathbb{E}_{\boldsymbol{x}}\left[(\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{\top}\right]=\boldsymbol{\mu}\boldsymbol{\mu}^{\top} + \boldsymbol{\Sigma}\end{equation}
Another commonly used statistic is its entropy:
\begin{equation}\mathcal{H} = \mathbb{E}_{\boldsymbol{x}}\left[-\log p(\boldsymbol{x})\right]=\frac{n}{2}(1 + \log 2\pi) + \frac{1}{2}\log \det(\boldsymbol{\Sigma}) \end{equation}
Its derivation can be found later on, as part of the discussion of the KL divergence.
Gaussian integral
The probability density function requires $\int p(\boldsymbol{x}) d\boldsymbol{x} = 1$, from which we can derive:
\begin{equation}\begin{aligned} \sqrt{(2\pi)^n \det(\boldsymbol{\Sigma})} =& \int\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right\}d\boldsymbol{x} \\ =& \int\exp\left\{-\frac{1}{2}\boldsymbol{x}^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{x}+\boldsymbol{\mu}^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{x}-\frac{1}{2}\boldsymbol{\mu}^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}\right\}d\boldsymbol{x} \end{aligned}\end{equation}
Setting $\boldsymbol{\omega} = \boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}$, we obtain the Gaussian integral
\begin{equation} \int\exp\left\{-\frac{1}{2}\boldsymbol{x}^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{x}+\boldsymbol{\omega}^{\top}\boldsymbol{x}\right\}d\boldsymbol{x} = \sqrt{(2\pi)^n \det(\boldsymbol{\Sigma})}\exp\left\{\frac{1}{2}\boldsymbol{\omega}^{\top}\boldsymbol{\Sigma}\boldsymbol{\omega}\right\}\label{eq:g-int} \end{equation}
Using this we can compute the characteristic function of the normal distribution
\begin{equation}\mathbb{E}_{\boldsymbol{x}}\left[\exp\left(\boldsymbol{\omega}^{\top}\boldsymbol{x}\right)\right]=\exp\left(\boldsymbol{\omega}^{\top}\boldsymbol{\mu}+\frac{1}{2}\boldsymbol{\omega}^{\top}\boldsymbol{\Sigma}\boldsymbol{\omega}\right)\\ \end{equation}
The characteristic function can be used to compute the moments of the normal distribution of any order.
Linear Algebra
Here we supplement some linear algebra background that will be used frequently in the derivations that follow. Again, this is only a "review" and not intended as a linear algebra tutorial.
Inner product and norm
First, let's define the inner product and the norm. For vectors $\boldsymbol{x}=(x_1,\cdots,x_n)$ and $\boldsymbol{y}=(y_1,\cdots,y_n)$, the inner product is defined as
\begin{equation}\langle\boldsymbol{x},\boldsymbol{y}\rangle = \sum_{i=1}^n x_i y_i\end{equation}
and the magnitude is defined as $\Vert \boldsymbol{x}\Vert = \sqrt{\langle\boldsymbol{x},\boldsymbol{x}\rangle}$. For matrices $\boldsymbol{A}=(a_{i,j}),\boldsymbol{B}=(b_{i,j})$ of shape $m\times n$, we define things analogously:
\begin{equation}\langle\boldsymbol{A},\boldsymbol{B}\rangle_F = \sum_{i=1}^m\sum_{j=1}^n a_{i,j} b_{i,j}\end{equation}
This is called the Frobenius inner product, and the corresponding $\Vert \boldsymbol{A}\Vert_F = \sqrt{\langle\boldsymbol{A},\boldsymbol{A}\rangle_F}$ is called the Frobenius norm. It's easy to see that the Frobenius inner product and norm are, in effect, just the ordinary vector inner product and norm applied after flattening the matrix into a vector.
Regarding the Frobenius inner product, one of the most important properties is the identity
\begin{equation}\langle\boldsymbol{A},\boldsymbol{B}\rangle_F = \text{Tr}\left(\boldsymbol{A}^{\top}\boldsymbol{B}\right) = \text{Tr}\left(\boldsymbol{B}\boldsymbol{A}^{\top}\right) = \text{Tr}\left(\boldsymbol{A}\boldsymbol{B}^{\top}\right) = \text{Tr}\left(\boldsymbol{B}^{\top}\boldsymbol{A}\right) \end{equation}
That is, the Frobenius inner product of two matrices can be converted into the trace of a matrix product, and swapping the order of multiplication does not change the result (it doesn't change the trace, though the overall result of the matrix product itself does change).
Positive definite symmetric matrices
Next, let's look at some properties of positive definite symmetric matrices. Let $\boldsymbol{\Sigma}$ be a positive definite symmetric matrix: symmetric means $\boldsymbol{\Sigma}^{\top}=\boldsymbol{\Sigma}$, and positive definite means that for any nonzero vector $\boldsymbol{\xi}\in\mathbb{R}^n$, we have $\boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} > 0$. It can be shown that if $\boldsymbol{\Sigma}_1,\boldsymbol{\Sigma}_2$ are both positive definite symmetric matrices, then $\boldsymbol{\Sigma}_1^{-1},\boldsymbol{\Sigma}_2^{-1},\boldsymbol{\Sigma}_1+\boldsymbol{\Sigma}_2$ is also positive definite symmetric. If $\boldsymbol{C} = \boldsymbol{B}^{\top}\boldsymbol{A}\boldsymbol{B}$ and $\boldsymbol{B}$ is invertible, then $\boldsymbol{C}$ is positive definite symmetric if and only if $\boldsymbol{A}$ is positive definite symmetric.
There is also the concept of positive semi-definiteness, meaning that for any nonzero vector $\boldsymbol{\xi}\in\mathbb{R}^n$, we have $\boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} \geq 0$, i.e., there may exist a nonzero vector $\boldsymbol{\xi}$ such that $\boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} = 0$. However, since positive definite matrices are dense within the set of positive semi-definite matrices, we won't strictly distinguish between them here, and will treat everything uniformly as positive definite.
Positive definite symmetric matrices have an important property: their SVD coincides with their eigendecomposition, i.e., they admit a decomposition of the form
\begin{equation}\boldsymbol{\Sigma} = \boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{\top}\end{equation}
where $\boldsymbol{U}$ is an orthogonal matrix and $\boldsymbol{\Lambda}$ is a diagonal matrix whose diagonal entries are all positive. A direct corollary of this result is that positive definite symmetric matrices can be "square-rooted," with square root $\boldsymbol{\Sigma}^{1/2} = \boldsymbol{U}\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top}$, where $\boldsymbol{\Lambda}^{1/2}$ means taking the square root of each diagonal entry; one can verify that the square root matrix is also positive definite symmetric. Conversely, any symmetric matrix that admits a square root must itself be positive definite symmetric.
Matrix calculus
Finally, when computing the Wasserstein distance we will also need some matrix calculus formulas. Readers unfamiliar with these can refer directly to the Wikipedia page on "Matrix Calculus". It's actually not too hard — we mainly need
\begin{equation}\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}\end{equation}
and everything else can be derived by combining this with the rules for the trace, for example
\begin{equation}\frac{\partial\,\text{Tr}\left(\boldsymbol{A}\boldsymbol{X}\boldsymbol{B}\right)}{\partial \boldsymbol{X}} = \frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{B}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{B}\boldsymbol{A}\end{equation}
KL Divergence
As a first attempt, let's compute the KL divergence (Kullback–Leibler divergence) between two Gaussian distributions. The KL divergence is one of the most commonly used measures between distributions, because taking the logarithm before integrating usually gives relatively simple results for distributions in the exponential family. It is also closely related to "entropy."
Result
The KL divergence between two probability distributions is defined as
\begin{equation}KL(p(\boldsymbol{x})\Vert q(\boldsymbol{x}))=\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[\log \frac{p(\boldsymbol{x})}{q(\boldsymbol{x})}\right]=\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[\log p(\boldsymbol{x})\right]+\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[-\log q(\boldsymbol{x})\right]\end{equation}
For two normal distributions, the result is
\begin{equation} KL(p(\boldsymbol{x})\Vert q(\boldsymbol{x}))=\frac{1}{2}\left[(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)-\log \det(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Sigma}_p) + \text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Sigma}_p\right) - n\right] \end{equation}
In particular, when $q$ is the standard normal distribution, this simplifies to
\begin{equation} KL(p(\boldsymbol{x})\Vert q(\boldsymbol{x}))=\frac{1}{2}\left[\Vert\boldsymbol{\mu}_p\Vert^2-\log \det(\boldsymbol{\Sigma}_p) + \text{Tr}(\boldsymbol{\Sigma}_p) - n\right]\end{equation}
Derivation
From the definition of the KL divergence, we mainly need to compute $\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[-\log q(\boldsymbol{x})\right]$:
\begin{equation}\begin{aligned} \mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[-\log q(\boldsymbol{x})\right] =&\, \mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[\frac{n}{2}\log 2\pi + \frac{1}{2}\log \det(\Sigma_q) + \frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_q)\right]\\ =&\,\frac{n}{2}\log 2\pi + \frac{1}{2}\log \det(\boldsymbol{\Sigma}_q) + \frac{1}{2}\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[(\boldsymbol{x}-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_q)\right] \end{aligned}\end{equation}
Now the trace identity mentioned earlier comes in handy:
\begin{equation}\begin{aligned} \mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[(\boldsymbol{x}-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_q)\right]=&\,\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[\text{Tr}\left((\boldsymbol{x}-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_q)\right)\right]\\ =&\,\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[\text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_q)(\boldsymbol{x}-\boldsymbol{\mu}_q)^{\top}\right)\right]\\ =&\,\text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[(\boldsymbol{x}-\boldsymbol{\mu}_q)(\boldsymbol{x}-\boldsymbol{\mu}_q)^{\top}\right]\right)\\ =&\,\text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[\boldsymbol{x}\boldsymbol{x}^{\top}-\boldsymbol{\mu}_q\boldsymbol{x}^{\top} - \boldsymbol{x}\boldsymbol{\mu}_q^{\top} + \boldsymbol{\mu}_q\boldsymbol{\mu}_q^{\top}\right]\right)\\ =&\,\text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\left(\boldsymbol{\Sigma}_p + \boldsymbol{\mu}_p\boldsymbol{\mu}_p^{\top}-\boldsymbol{\mu}_q\boldsymbol{\mu}_p^{\top} - \boldsymbol{\mu}_p\boldsymbol{\mu}_q^{\top} + \boldsymbol{\mu}_q\boldsymbol{\mu}_q^{\top}\right)\right)\\ =&\,\text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Sigma}_p + \boldsymbol{\Sigma}_q^{-1}(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)^{\top}\right)\\ =&\,\text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Sigma}_p\right) + (\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)\\ \end{aligned}\end{equation}
Note that when $\boldsymbol{\mu}_q=\boldsymbol{\mu}_p,\boldsymbol{\Sigma}_q=\boldsymbol{\Sigma}_p$, the expression above equals $n$, which corresponds exactly to the entropy of the normal distribution. So we finally obtain
\begin{equation}\begin{aligned} KL(p(\boldsymbol{x})\Vert q(\boldsymbol{x}))=&\,\frac{1}{2}\left[n\log 2\pi + \log \det(\boldsymbol{\Sigma}_q) + \text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Sigma}_p\right) + (\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)\right] \\ &\,\qquad- \frac{1}{2}\left[n\log 2\pi + \log \det(\boldsymbol{\Sigma}_p) + n\right]\\ =&\,\frac{1}{2}\left[(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{\mu}_p-\boldsymbol{\mu}_q)-\log \det(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Sigma}_p) + \text{Tr}\left(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Sigma}_p\right) - n\right] \end{aligned}\end{equation}
Bhattacharyya Distance
Next, let's look at the Bhattacharyya distance, defined as
\begin{equation}BD(p(\boldsymbol{x}), q(\boldsymbol{x})) = -\log \int \sqrt{p(\boldsymbol{x}) q(\boldsymbol{x})} d\boldsymbol{x}\end{equation}
Related to it is a concept called the "Hellinger distance", whose square is defined as $\frac{1}{2}\int \left(\sqrt{p(\boldsymbol{x})} - \sqrt{q(\boldsymbol{x})}\right)^2 d\boldsymbol{x}$; expanding this out reveals that it is essentially equivalent to the Bhattacharyya distance.
Result
For two normal distributions, the Bhattacharyya distance is
\begin{equation} BD(p(\boldsymbol{x}), q(\boldsymbol{x})) = \frac{1}{2}\log \frac{\det(\boldsymbol{\Sigma})}{\sqrt{\det(\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q)}} + \frac{1}{8}(\boldsymbol{\mu}_p - \boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_p - \boldsymbol{\mu}_q) \end{equation}
where $\boldsymbol{\Sigma}=\frac{1}{2}(\boldsymbol{\Sigma}_p + \boldsymbol{\Sigma}_q)$. Notice that the result is symmetric, which follows from the fact that the definition of the Bhattacharyya distance is itself symmetric.
When one of the two distributions is the standard normal distribution, the result doesn't simplify noticeably, so we won't write it out separately here.
Derivation
By definition, the Bhattacharyya distance between two normal distributions is the negative logarithm of the following integral:
\begin{equation}\begin{aligned} &\qquad\int \sqrt{p(\boldsymbol{x}) q(\boldsymbol{x})} d\boldsymbol{x}=\frac{1}{\sqrt[4]{(2\pi)^{2n}\det(\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q)}}\times \\ &\int \exp\left\{-\frac{1}{4}(\boldsymbol{x}-\boldsymbol{\mu}_p)^{\top}\boldsymbol{\Sigma}_p^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_p)-\frac{1}{4}(\boldsymbol{x}-\boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_q)\right\}d\boldsymbol{x} \end{aligned}\end{equation}
Writing $\boldsymbol{y}=\boldsymbol{x}-\boldsymbol{\mu}_p, \boldsymbol{\Delta}=\boldsymbol{\mu}_p - \boldsymbol{\mu}_q$, the integral part can be changed of variables to
\begin{equation}\begin{aligned} &\int \exp\left\{-\frac{1}{4}\boldsymbol{y}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{y}-\frac{1}{4}(\boldsymbol{y}+\boldsymbol{\Delta})^{\top}\boldsymbol{\Sigma}_q^{-1}(\boldsymbol{y}+\boldsymbol{\Delta})\right\}d\boldsymbol{y}\\ =&\int \exp\left\{-\frac{1}{4}\boldsymbol{y}^{\top}\left(\boldsymbol{\Sigma}_p^{-1}+\boldsymbol{\Sigma}_q^{-1}\right)\boldsymbol{y}-\frac{1}{2}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}_q^{-1}\boldsymbol{y} - \frac{1}{4}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Delta}\right\}d\boldsymbol{y}\\ =&\int \exp\left\{-\frac{1}{2}\boldsymbol{y}^{\top}\left(\boldsymbol{\Sigma}_p^{-1}\boldsymbol{\Sigma}\boldsymbol{\Sigma}_q^{-1}\right)\boldsymbol{y}-\frac{1}{2}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}_q^{-1}\boldsymbol{y} - \frac{1}{4}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Delta}\right\}d\boldsymbol{y}\end{aligned}\end{equation}
where $\boldsymbol{\Sigma}=\frac{1}{2}(\boldsymbol{\Sigma}_p + \boldsymbol{\Sigma}_q)$. Using the Gaussian integral formula $\eqref{eq:g-int}$ introduced earlier, the integral evaluates to
\begin{equation}\begin{aligned} &\,\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma}_p^{-1}\boldsymbol{\Sigma}\boldsymbol{\Sigma}_q^{-1})^{-1}}\exp\left\{\frac{1}{8}\left(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Delta}\right)^{\top}\left(\boldsymbol{\Sigma}_p^{-1}\boldsymbol{\Sigma}\boldsymbol{\Sigma}_q^{-1}\right)^{-1}\left(\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Delta}\right)-\frac{1}{4}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}_q^{-1}\boldsymbol{\Delta}\right\}\\ =&\,\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}_p)}\exp\left\{\frac{1}{8}\boldsymbol{\Delta}^{\top}\left(\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q^{-1} - 2\boldsymbol{\Sigma}_q^{-1}\right)\boldsymbol{\Delta}\right\}\\ =&\,\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}_p)}\exp\left\{\frac{1}{8}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q^{-1} - 2\boldsymbol{\Sigma}\boldsymbol{\Sigma}_q^{-1}\right)\boldsymbol{\Delta}\right\}\\ =&\,\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}_p)}\exp\left\{-\frac{1}{8}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{\Delta}\right\} \end{aligned}\end{equation}
so finally
\begin{equation}\begin{aligned} BD(p(\boldsymbol{x}), q(\boldsymbol{x})) =&\, -\log \left[\frac{\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}_p)}}{\sqrt[4]{(2\pi)^{2n}\det(\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q)}}\exp\left\{-\frac{1}{8}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{\Delta}\right\}\right] \\ =&\, -\log \left[\frac{\sqrt[4]{\det(\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q)}}{\sqrt{\det\left(\boldsymbol{\Sigma}\right)}}\exp\left\{-\frac{1}{8}\boldsymbol{\Delta}^{\top}\boldsymbol{\Sigma}^{-1}\boldsymbol{\Delta}\right\}\right]\\ =&\,\frac{1}{2}\log \frac{\det(\boldsymbol{\Sigma})}{\sqrt{\det(\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q)}} + \frac{1}{8}(\boldsymbol{\mu}_p - \boldsymbol{\mu}_q)^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}_p - \boldsymbol{\mu}_q) \end{aligned} \end{equation}
W Distance
If readers would like to learn more about probability divergences in general, they can consult the book Statistical Inference Based on Divergence Measures. We now turn to another class of probability metrics — the W (Wasserstein) distance, based on optimal transport.
Following the notation of From the Wasserstein Distance and Duality Theory to WGAN, the W distance is defined as follows:
\begin{equation}\begin{aligned} \mathcal{W}_{\rho}[p,q]=&\,\left(\inf_{\gamma\in \Pi[p,q]} \iint \gamma(\boldsymbol{x},\boldsymbol{y}) \Vert\boldsymbol{x} - \boldsymbol{y}\Vert^{\rho} d\boldsymbol{x}d\boldsymbol{y}\right)^{1/\rho}\\ =&\,\left(\inf_{\gamma\in \Pi[p,q]} \mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\gamma(\boldsymbol{x},\boldsymbol{y})} \left[\Vert\boldsymbol{x} - \boldsymbol{y}\Vert^{\rho}\right]\right)^{1/\rho} \end{aligned}\end{equation}
According to the sources the author has found, among two normal distributions, an analytic solution has only been worked out for the case $\rho=2$, so below we will focus on computing $\mathcal{W}_2[p,q]$, and for simplicity we write
\begin{equation}\mathcal{W}_2^2[p,q] = \left(\mathcal{W}_2[p,q]\right)^2\end{equation}
Result
Interestingly, there are two different versions of the result for the W distance between two normal distributions circulating in the literature. Both versions have some degree of recognition, yet I have not seen any source that explicitly states their equivalence. The two versions come from different papers, and have even been given different names.
Version 1
First, the more widely circulated version — searching for "Wasserstein distance between normal distributions" will most often turn up this result:
\begin{equation}\mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + \text{Tr}(\boldsymbol{\Sigma}_p) + \text{Tr}(\boldsymbol{\Sigma}_q) - 2\text{Tr}((\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2})\label{eq:w-v1}\end{equation}Some readers might be puzzled as to "why this isn't symmetric with respect to $p,q$" — in fact it is symmetric in $p,q$, because
\begin{equation}\begin{aligned}\text{Tr}((\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2})=&\,\text{Tr}((\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2}\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{\Sigma}_q^{-1/2}\boldsymbol{\Sigma}_q^{1/2}\boldsymbol{\Sigma}_p^{1/2})\\ =&\,\text{Tr}(\boldsymbol{\Sigma}_q^{1/2}\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2}\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{\Sigma}_q^{-1/2}) \end{aligned}\end{equation}
We can then directly verify that $(\boldsymbol{\Sigma}_q^{1/2}\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2}\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{\Sigma}_q^{-1/2})^2=\boldsymbol{\Sigma}_q^{1/2}\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q^{1/2}$, so $\text{Tr}((\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2})=\text{Tr}((\boldsymbol{\Sigma}_q^{1/2}\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q^{1/2})^{1/2})$ holds.
Version 2
The second, slightly simpler-looking version gives the result:
\begin{equation}\mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + \text{Tr}(\boldsymbol{\Sigma}_p) + \text{Tr}(\boldsymbol{\Sigma}_q) - 2\text{Tr}((\boldsymbol{\Sigma}_p\boldsymbol{\Sigma}_q)^{1/2})\label{eq:w-v2}\end{equation}
This version is usually called the "Fréchet distance", and generally one only finds this result by searching for "Fréchet distance between normal distributions." The FID (Fréchet Inception Distance) metric, commonly used to evaluate GANs, is computed using precisely this formula. Its symmetry with respect to $p,q$ can be shown in the same way as before, or it follows directly from the equivalence discussion below.
Equivalence
By rights, since the second version is a bit simpler, it should be the one that has become the standard. So why do two different versions still circulate, is something I find rather puzzling. Theoretically, proving the equivalence of the two versions isn't hard: using the trace identity, we have
\begin{equation}\begin{aligned}\text{Tr}((\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2})=&\,\text{Tr}((\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2}\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{\Sigma}_p^{1/2})\\ =&\,\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2}\boldsymbol{\Sigma}_p^{-1/2}) \end{aligned}\end{equation}
and then we just need to verify $(\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2}\boldsymbol{\Sigma}_p^{-1/2})^2=\boldsymbol{\Sigma}_p \boldsymbol{\Sigma}_q$ directly.
A special case
In particular, if $\boldsymbol{\Sigma}_p,\boldsymbol{\Sigma}_q$ commute under multiplication, the result simplifies to a very intuitive form:
\begin{equation}\mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + \Vert \boldsymbol{\Sigma}_p^{1/2} - \boldsymbol{\Sigma}_q^{1/2}\Vert_F^2\label{eq:w-jiaohuan}\end{equation}
Why is this so intuitive? Because the parameters of the normal distribution are $\boldsymbol{\mu},\boldsymbol{\Sigma}$, so comparing the difference between two normal distributions really amounts to comparing the difference between $\boldsymbol{\mu},\boldsymbol{\Sigma}$. Following common practice in machine learning, a metric that comes to mind readily is the squared error
\begin{equation}\mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + \Vert \boldsymbol{\Sigma}_p - \boldsymbol{\Sigma}_q\Vert_F^2\end{equation}
But from a physical point of view, this metric is problematic: if we regard $\boldsymbol{\mu}$ as having the dimension of length, then $\boldsymbol{\Sigma}$ has the dimension of length squared, so $\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2$ and $\Vert \boldsymbol{\Sigma}_p - \boldsymbol{\Sigma}_q\Vert_F^2$ have different dimensions and cannot simply be added together. To make the dimensions consistent, the intuitive fix is to "take the square root" of $\boldsymbol{\Sigma}$ before computing the squared error, which gives us equation $\eqref{eq:w-jiaohuan}$.
In particular, when $q$ is the standard normal distribution, the result simplifies to
\begin{equation}\mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p\Vert^2 + \Vert \boldsymbol{\Sigma}_p^{1/2} - \boldsymbol{I}\Vert_F^2\end{equation}
Derivation 1
Let us now present the first proof, mainly following the paper A class of Wasserstein metrics for probability distributions. In addition, The distance between two random vectors with given dispersion matrices provides a similar proof, which readers may also consult.
The derivation below has been simplified by the author, and is somewhat simpler than the proof in the original paper, though it still unavoidably involves a fair amount of linear algebra. We will present it in several parts.
Removing the mean
Without loss of generality, we may consider only distributions $p,q$ with zero mean. Indeed, if $p,q$ does not have zero mean, let $\tilde{p},\tilde{q}$ denote the corresponding zero-mean distribution; then we have
\begin{equation}\begin{aligned} &\,\mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\gamma(\boldsymbol{x},\boldsymbol{y})}\left[\Vert \boldsymbol{x} - \boldsymbol{y}\Vert^2\right] \\ =&\, \mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\tilde{\gamma}(\boldsymbol{x},\boldsymbol{y})}\left[\Vert (\boldsymbol{x} + \boldsymbol{\mu}_p) - (\boldsymbol{y} + \boldsymbol{\mu}_q)\Vert^2 \right]\\ =&\,\mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\tilde{\gamma}(\boldsymbol{x},\boldsymbol{y})}\left[\Vert \boldsymbol{x} - \boldsymbol{y}\Vert^2 + \Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + 2\langle\boldsymbol{x} - \boldsymbol{y}, \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\rangle\right]\\ =&\,\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + \mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\tilde{\gamma}(\boldsymbol{x},\boldsymbol{y})}\left[\Vert \boldsymbol{x} - \boldsymbol{y}\Vert^2 \right] \end{aligned}\end{equation}
This result implies that
\begin{equation}\mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + \mathcal{W}_2^2[\tilde{p},\tilde{q}]\end{equation}
So we only need to work out the Wasserstein distance in the case of zero mean, and then add $\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2$ to obtain the result for the general case.
Pure algebra
Now assume $p,q$ both have mean zero, and consider computing
\begin{equation}\begin{aligned} \mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\gamma(\boldsymbol{x},\boldsymbol{y})}\left[\Vert \boldsymbol{x} - \boldsymbol{y}\Vert^2\right] =&\, \mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\gamma(\boldsymbol{x},\boldsymbol{y})}\left[\boldsymbol{x}^{\top} \boldsymbol{x} + \boldsymbol{y}^{\top} \boldsymbol{y} - 2\boldsymbol{y}^{\top} \boldsymbol{x}\right]\\ =&\, \mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\gamma(\boldsymbol{x},\boldsymbol{y})}\left[\text{Tr}\left(\boldsymbol{x} \boldsymbol{x}^{\top} + \boldsymbol{y} \boldsymbol{y}^{\top} - 2\boldsymbol{x}\boldsymbol{y}^{\top} \right)\right]\\ =&\, \text{Tr}\left(\mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\gamma(\boldsymbol{x},\boldsymbol{y})}\left[\boldsymbol{x} \boldsymbol{x}^{\top} + \boldsymbol{y} \boldsymbol{y}^{\top} - 2\boldsymbol{x}\boldsymbol{y}^{\top}\right]\right)\\ =&\, \text{Tr}(\boldsymbol{\Sigma}_p) + \text{Tr}(\boldsymbol{\Sigma}_q) - 2\text{Tr}(\boldsymbol{C}) \end{aligned}\end{equation}
where
\begin{equation}\boldsymbol{\Sigma}_{\gamma}= \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix}=\mathbb{E}_{(\boldsymbol{x},\boldsymbol{y})\sim\gamma(\boldsymbol{x},\boldsymbol{y})}\left[\begin{pmatrix}\boldsymbol{x} \\ \boldsymbol{y}\end{pmatrix}\begin{pmatrix}\boldsymbol{x}^{\top} & \boldsymbol{y}^{\top}\end{pmatrix}\right]\end{equation}
forms the covariance matrix of the joint distribution $\gamma$. We know that a covariance matrix is positive definite symmetric, so from an algebraic point of view, the problem becomes:
Given that $\boldsymbol{\Sigma}_{\gamma}= \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix}$ is positive definite symmetric, find the maximum value of $\text{Tr}(\boldsymbol{C})$.
Schur complement
To do this, we need to use the following identity involving the "Schur complement":
\begin{equation} \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix} = \begin{pmatrix} \boldsymbol{I} & \boldsymbol{0}\\ \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1} & \boldsymbol{I}\end{pmatrix} \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{\Sigma}_q - \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}\end{pmatrix} \begin{pmatrix} \boldsymbol{I} & \boldsymbol{\Sigma}_p^{-1}\boldsymbol{C} \\ \boldsymbol{0} & \boldsymbol{I}\end{pmatrix} \end{equation}
where the symmetric matrix $\boldsymbol{S} = \boldsymbol{\Sigma}_q - \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$ is called the "Schur complement". This decomposition has the form $\boldsymbol{B}^{\top}\boldsymbol{A}\boldsymbol{B}$, and for it to be positive definite, $\boldsymbol{A}$ must be positive definite; since $\boldsymbol{\Sigma}_p$ is already positive definite, this requires $\boldsymbol{S}$ to be positive definite.
Separating the parameters
Let's try to separate the parameters, i.e., solve for $\boldsymbol{C}$ from $\boldsymbol{S} = \boldsymbol{\Sigma}_q - \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$. First, rearranging gives $\boldsymbol{\Sigma}_q - \boldsymbol{S} = \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$. Since $\boldsymbol{\Sigma}_p$ is positive definite symmetric, so is $\boldsymbol{\Sigma}_p^{-1}$, and hence so is $\boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$, which therefore has a positive definite symmetric square root; that is, there exists a positive definite symmetric matrix $\boldsymbol{R}$ such that
\begin{equation}\boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C} = \boldsymbol{R}^2\quad\Leftrightarrow\quad \left(\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{C}\boldsymbol{R}^{-1}\right)^{\top}\left(\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{C}\boldsymbol{R}^{-1}\right)=\boldsymbol{I}\end{equation}
This shows that $\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{C}\boldsymbol{R}^{-1}$ is an orthogonal matrix, which we denote $\boldsymbol{O}$; hence $\boldsymbol{C} = \boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R}$.
Lagrange multipliers
At this point, the variables are $\boldsymbol{O}$ and $\boldsymbol{R}$, and we want to find the maximum of $\text{Tr}(\boldsymbol{C})=\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R})$. Let's first fix $\boldsymbol{R}$ and find the value of $\boldsymbol{O}$ that maximizes it — this amounts to maximizing $\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R})$ subject to the constraint $\boldsymbol{O}^{\top}\boldsymbol{O}=\boldsymbol{I}$. We use "Lagrange multipliers": introduce a new parameter matrix $\boldsymbol{W}$, converting this into the following unconstrained extremum problem
\begin{equation}F = \text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R}) - \frac{1}{2}\text{Tr}(\boldsymbol{W}(\boldsymbol{O}^{\top}\boldsymbol{O} - \boldsymbol{I}))\end{equation}
Taking the derivative:
\begin{equation}\begin{aligned} &\frac{\partial F}{\partial \boldsymbol{O}} = \boldsymbol{0} \quad \Rightarrow\quad \boldsymbol{R}\boldsymbol{\Sigma}_p^{1/2} = \boldsymbol{W}\boldsymbol{O}^{\top}\\ &\frac{\partial F}{\partial \boldsymbol{W}} = \boldsymbol{0} \quad \Rightarrow\quad \boldsymbol{O}^{\top}\boldsymbol{O} = \boldsymbol{I}\\ \end{aligned}\end{equation}
First, note that $\boldsymbol{O}^{\top}\boldsymbol{O} - \boldsymbol{I}$ is symmetric, and therefore the corresponding parameter matrix $\boldsymbol{W}$ is also symmetric, so we obtain
\begin{equation}\left(\boldsymbol{O}\boldsymbol{W}\boldsymbol{O}^{\top}\right)^2=\left(\boldsymbol{W}\boldsymbol{O}^{\top}\right)^{\top}\left(\boldsymbol{W}\boldsymbol{O}^{\top}\right)=\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{R}^2\boldsymbol{\Sigma}_p^{1/2}\end{equation}
that is, $\boldsymbol{O}\boldsymbol{W}\boldsymbol{O}^{\top}=(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{R}^2\boldsymbol{\Sigma}_p^{1/2})^{1/2}$, so at this point
\begin{equation}\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R})=\text{Tr}(\boldsymbol{O}\boldsymbol{R}\boldsymbol{\Sigma}_p^{1/2})=\text{Tr}(\boldsymbol{O}\boldsymbol{W}\boldsymbol{O}^{\top})=\text{Tr}((\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{R}^2\boldsymbol{\Sigma}_p^{1/2})^{1/2})\end{equation}
Inequality
Finally, we need to pin down $\boldsymbol{R}$. Recall the definition of $\boldsymbol{R}$: we have $\boldsymbol{R}^2=\boldsymbol{\Sigma}_q - \boldsymbol{S}$, where $\boldsymbol{S}$ is a positive definite matrix. Intuitively, the maximum should be attained when $\boldsymbol{S}=\boldsymbol{0}$, and indeed this is the case — it's a direct corollary of "Weyl's inequality".
By Weyl's inequality, if matrices $\boldsymbol{A},\boldsymbol{B}$ are both positive definite symmetric, and the eigenvalues of $\boldsymbol{A},\boldsymbol{B},\boldsymbol{A}+\boldsymbol{B}$, sorted in increasing order, are $0\leq\lambda_1^{(A)} \leq \dots \leq \lambda_n^{(A)}$, $\lambda_1^{(B)} \leq \dots \leq \lambda_n^{(B)}$, and $0\leq\lambda_1^{(A+B)} \leq \dots \leq \lambda_n^{(A+B)}$ respectively, then for any $1\leq i \leq n$ we have $\lambda_i^{(A)}\leq \lambda_i^{(A+B)}$ and $\lambda_i^{(B)}\leq \lambda_i^{(A+B)}$. That is:
The eigenvalues of the sum of positive definite symmetric matrices are each, correspondingly, larger than the respective eigenvalues of the summands.
With this result in hand, things become simple: let the eigenvalues of $(\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2})^{1/2}$ be $0 \leq \lambda_1 \leq \dots \leq \lambda_n$, so its trace is $\lambda_1 + \dots + \lambda_n$. Correspondingly, the eigenvalues of $\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2}$ are $0 \leq \lambda_1^2 \leq \dots \leq \lambda_n^2$. Note that $\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2}$ is positive definite symmetric (symmetry is obvious, and since it admits a square root, it's positive definite), and $\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{S}\boldsymbol{\Sigma}_p^{1/2}$ is also positive definite symmetric (since $\boldsymbol{S}$ is positive definite symmetric). Hence their eigenvalues cannot exceed those of their sum — i.e., the eigenvalues of $\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2}$ — so the maximum of each eigenvalue of $(\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2})^{1/2}$ (and hence the maximum of the trace) is attained at $\boldsymbol{S}=\boldsymbol{0}$.
As for the proof of Weyl's inequality, it mainly relies on the Rayleigh quotient and the Courant–Fischer theorem. Interested readers can look up these two topics first, and then consult the proof of Weyl's inequality — once you're familiar with those two ingredients, Weyl's inequality basically falls into place naturally.
Derivation 2
Here we present another, simpler proof; the original can be found in The Fréchet distance between multivariate normal distributions. Compared to the first proof, this one is more direct and, in particular, doesn't require as much pure linear algebra. The derivation below has again been further simplified by the author, and should be easier to follow than the original paper.
In this derivation, the two steps "removing the mean" and "pure algebra" are the same as in "Derivation 1," so we won't repeat them. So at this point the problem has already been reduced to:
Given that $\boldsymbol{\Sigma}_{\gamma}= \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix}$ is positive definite symmetric, find the maximum value of $\text{Tr}(\boldsymbol{C})$.
Block matrix
Since $\boldsymbol{\Sigma}_{\gamma}$ is positive definite symmetric, it can necessarily be expressed in the form $\boldsymbol{D}\boldsymbol{D}^{\top}$. Let us write $\boldsymbol{D}$ as the block matrix $\begin{pmatrix}\boldsymbol{A} \\ \boldsymbol{B}\end{pmatrix}$, where $\boldsymbol{A},\boldsymbol{B}\in \mathbb{R}^{n\times 2n}$. Then
\begin{equation}\begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix} = \begin{pmatrix}\boldsymbol{A} \\ \boldsymbol{B}\end{pmatrix} \begin{pmatrix}\boldsymbol{A}^{\top} & \boldsymbol{B}^{\top}\end{pmatrix} = \begin{pmatrix} \boldsymbol{A}\boldsymbol{A}^{\top} & \boldsymbol{A}\boldsymbol{B}^{\top}\\ \boldsymbol{B}\boldsymbol{A}^{\top} & \boldsymbol{B}\boldsymbol{B}^{\top}\end{pmatrix}\end{equation}
and correspondingly $\boldsymbol{\Sigma}_p=\boldsymbol{A}\boldsymbol{A}^{\top}, \boldsymbol{\Sigma}_q=\boldsymbol{B}\boldsymbol{B}^{\top},\boldsymbol{C}=\boldsymbol{A}\boldsymbol{B}^{\top}$.
Lagrange multipliers
Under this parametrization, the problem becomes:
Given $\boldsymbol{A}\boldsymbol{A}^{\top}=\boldsymbol{\Sigma}_p, \boldsymbol{B}\boldsymbol{B}^{\top}=\boldsymbol{\Sigma}_q$, find the maximum value of $\text{Tr}(\boldsymbol{A}\boldsymbol{B}^{\top})$.
This is a constrained maximization problem, so we again use "Lagrange multipliers": introduce a new parameter matrix $\boldsymbol{W}_p, \boldsymbol{W}_q$, converting it into the following unconstrained extremum problem
\begin{equation}F = \text{Tr}(\boldsymbol{A}\boldsymbol{B}^{\top}) - \frac{1}{2}\text{Tr}(\boldsymbol{W}_p(\boldsymbol{A}\boldsymbol{A}^{\top} - \boldsymbol{\Sigma}_p)) - \frac{1}{2}\text{Tr}(\boldsymbol{W}_q(\boldsymbol{B}\boldsymbol{B}^{\top} - \boldsymbol{\Sigma}_q))\end{equation}
Taking the derivative:
\begin{equation}\begin{aligned} &\frac{\partial F}{\partial \boldsymbol{A}} = \boldsymbol{0} \quad \Rightarrow\quad \boldsymbol{B}^{\top} = \boldsymbol{A}^{\top}\boldsymbol{W}_p\\ &\frac{\partial F}{\partial \boldsymbol{B}} = \boldsymbol{0} \quad \Rightarrow\quad \boldsymbol{A}^{\top} = \boldsymbol{B}^{\top}\boldsymbol{W}_q\\ &\frac{\partial F}{\partial \boldsymbol{W}_p} = \boldsymbol{0} \quad \Rightarrow\quad \boldsymbol{A}\boldsymbol{A}^{\top} = \boldsymbol{\Sigma}_p\\ &\frac{\partial F}{\partial \boldsymbol{W}_q} = \boldsymbol{0} \quad \Rightarrow\quad \boldsymbol{B}\boldsymbol{B}^{\top} = \boldsymbol{\Sigma}_q\\ \end{aligned}\end{equation}
Noting that both $\boldsymbol{A}\boldsymbol{A}^{\top} - \boldsymbol{\Sigma}_p$ and $\boldsymbol{B}\boldsymbol{B}^{\top} - \boldsymbol{\Sigma}_q$ are symmetric, the corresponding parameter matrix $\boldsymbol{W}_p, \boldsymbol{W}_q$ must also be symmetric, so at this point
\begin{equation}\boldsymbol{\Sigma}_q = \boldsymbol{B}\boldsymbol{B}^{\top} = \left(\boldsymbol{A}^{\top}\boldsymbol{W}_p\right)^{\top}\left(\boldsymbol{A}^{\top}\boldsymbol{W}_p\right)=\boldsymbol{W}_p\boldsymbol{A}\boldsymbol{A}^{\top}\boldsymbol{W}_p=\boldsymbol{W}_p\boldsymbol{\Sigma}_p\boldsymbol{W}_p\\ \end{equation}
Setting $\boldsymbol{W}_p=\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{R}\boldsymbol{\Sigma}_p^{-1/2}$ and substituting into the above gives $\boldsymbol{\Sigma}_q=\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{R}^2\boldsymbol{\Sigma}_p^{-1/2}$, that is,
\begin{equation}\boldsymbol{R} = (\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2} \end{equation}
while
\begin{equation}\begin{aligned} \text{Tr}(\boldsymbol{A}\boldsymbol{B}^{\top}) =&\, \text{Tr}(\boldsymbol{A}\boldsymbol{A}^{\top}\boldsymbol{W}_p)=\text{Tr}(\boldsymbol{\Sigma}_p\boldsymbol{W}_p)\\ =&\,\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{R}\boldsymbol{\Sigma}_p^{-1/2}) = \text{Tr}(\boldsymbol{R}\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{\Sigma}_p^{1/2})\\ =&\,\text{Tr}(\boldsymbol{R}) \end{aligned}\end{equation}
Summary
This post has worked out in detail the KL divergence, Bhattacharyya distance, and Wasserstein distance between two multivariate normal distributions, giving explicit analytic solutions for each. These results can serve, in certain settings, as regularization terms on latent variables, to constrain the distribution of the latent variables. Beyond that, this post can also serve as a fairly challenging set of linear algebra exercises for readers to work through on their own.
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.