The normal distribution is one of the most common continuous probability distributions. It is the maximum entropy distribution given a specified mean and covariance (refer to "Entropy-less: From Entropy and the Maximum Entropy Principle to Maximum Entropy Models (II)"). It can also be viewed as a second-order approximation of any continuous distribution, holding a status equivalent to the linear approximation of a general function. From this perspective, the normal distribution can be considered the simplest continuous distribution. Precisely because of its simplicity, analytical solutions can be derived for many estimators.
This article primarily calculates several metrics between two multivariate normal distributions, including KL divergence, Bhattacharyya distance, and Wasserstein distance, all of which have explicit analytical solutions.
The Normal Distribution
Here we briefly review some basic knowledge of the normal distribution. Note that this is merely a review and is not intended to serve as an introductory tutorial.
Probability Density
The normal distribution, also known as the Gaussian distribution, is a continuous probability distribution defined on \mathbb{R}^n. Its probability density function (PDF) is: 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\} Where \boldsymbol{x}, \boldsymbol{\mu} \in \mathbb{R}^n, \boldsymbol{\mu} is the mean vector (vectors in this article are column vectors by default), and \boldsymbol{\Sigma} \in \mathbb{R}^{n \times n} is the covariance matrix, which is required to be symmetric and positive definite. As we can see, a normal distribution is uniquely determined by \boldsymbol{\mu} and \boldsymbol{\Sigma}; thus, it is easy to imagine that its statistics are functions of \boldsymbol{\mu} and \boldsymbol{\Sigma}. When \boldsymbol{\mu}=\boldsymbol{0} and \boldsymbol{\Sigma}=\boldsymbol{I}, the corresponding distribution is called the "standard normal distribution."
Basic Properties
Generally, the basic statistics are the mean and variance, which correspond to the two parameters of the normal distribution: \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} From this, the result for the second moment can also be derived: \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} Another commonly used statistic is its entropy: \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}) The calculation process can refer to the derivation of the KL divergence later.
Gaussian Integral
The probability density function implies \int p(\boldsymbol{x}) d\boldsymbol{x} = 1, which leads to: \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} Let \boldsymbol{\omega} = \boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}, then we obtain the Gaussian integral: \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} Using this, we can calculate the characteristic function of the normal distribution: \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) The characteristic function can be used to calculate the moments of the normal distribution of various orders.
Linear Algebra
Here we supplement some linear algebra basics that will be frequently used in the subsequent derivations. Similarly, this is only a "review" and does not serve as a linear algebra tutorial.
Inner Product and Norm
First, let’s define the inner product and norm. For vectors \boldsymbol{x}=(x_1,\cdots,x_n) and \boldsymbol{y}=(y_1,\cdots,y_n), the inner product is defined as: \langle\boldsymbol{x},\boldsymbol{y}\rangle = \sum_{i=1}^n x_i y_i And the norm is defined as \Vert \boldsymbol{x}\Vert = \sqrt{\langle\boldsymbol{x},\boldsymbol{x}\rangle}. For m\times n matrices \boldsymbol{A}=(a_{i,j}) and \boldsymbol{B}=(b_{i,j}), we define them in a similar way: \langle\boldsymbol{A},\boldsymbol{B}\rangle_F = \sum_{i=1}^m\sum_{j=1}^n a_{i,j} b_{i,j} 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 is easy to see that the Frobenius inner product and norm are essentially treating the flattened matrix as a regular vector for calculation.
Regarding the Frobenius inner product, one of the most critical properties is the identity: \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) That is, the Frobenius inner product of matrices can be converted into the trace of matrix multiplication, and swapping the multiplication order does not change the result (it does not change the trace result, although the overall result of matrix multiplication changes).
Positive Definite Symmetry
Next, let’s look at some properties of positive definite symmetric matrices. \boldsymbol{\Sigma} is a positive definite symmetric matrix; "symmetric" means \boldsymbol{\Sigma}^{\top}=\boldsymbol{\Sigma}, and "positive definite" means that for any non-zero vector \boldsymbol{\xi}\in\mathbb{R}^n, \boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} > 0. It can be proven 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 are also positive definite symmetric matrices. If \boldsymbol{C} = \boldsymbol{B}^{\top}\boldsymbol{A}\boldsymbol{B} and \boldsymbol{B} is an invertible matrix, then \boldsymbol{C} is positive definite symmetric if and only if \boldsymbol{A} is positive definite symmetric.
There is also the concept of semi-positive definiteness, which means for any non-zero vector \boldsymbol{\xi}\in\mathbb{R}^n, \boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} \geq 0, meaning there may exist a non-zero vector \boldsymbol{\xi} such that \boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} = 0. However, considering that positive definite matrices are dense within semi-positive definite matrices, we will not strictly distinguish between positive definite and semi-positive definite and will treat them uniformly as positive definite matrices.
Positive definite symmetric matrices have an important property: their SVD decomposition is consistent with their eigenvalue decomposition, i.e., they have a decomposition of the following form: \boldsymbol{\Sigma} = \boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{\top} where \boldsymbol{U} is an orthogonal matrix and \boldsymbol{\Lambda} is a diagonal matrix with positive diagonal elements. A direct corollary of this result is that positive definite symmetric matrices can have a "square root," defined as \boldsymbol{\Sigma}^{1/2} = \boldsymbol{U}\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top}, where \boldsymbol{\Lambda}^{1/2} refers to taking the square root of each diagonal element. It can be verified that the square root matrix is also positive definite symmetric. Conversely, a symmetric matrix that can have a square root must also be positive definite symmetric.
Matrix Calculus
Finally, when calculating the Wasserstein distance, some matrix calculus formulas are needed. Readers who are unfamiliar can refer to the Wikipedia entry for "Matrix Calculus." Of course, it is not difficult; the main formula used is: \frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A} The rest can be derived by combining the trace operation formulas, such as: \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}
KL Divergence
As a first attempt, let’s calculate the KL divergence (Kullback-Leibler divergence) between two Gaussian distributions. KL divergence is one of the most commonly used distribution metrics because it requires taking the logarithm before integration, which usually yields relatively simple results for exponential family distributions. Furthermore, it has a close connection with "entropy."
Calculation Result
The KL divergence between two probability distributions is defined as: 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] For two normal distributions, the calculation result is: 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] In particular, when q is the standard normal distribution, the result simplifies to: 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]
Derivation Process
From the definition of KL divergence, we mainly need to calculate \mathbb{E}_{\boldsymbol{x}\sim p(\boldsymbol{x})}\left[-\log q(\boldsymbol{x})\right]: \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} Now, the trace identity can be put to use: \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} Note that when \boldsymbol{\mu}_q=\boldsymbol{\mu}_p and \boldsymbol{\Sigma}_q=\boldsymbol{\Sigma}_p, the above expression equals n, which corresponds to the entropy of the normal distribution. Thus, we finally obtain: \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}
Bhattacharyya Distance
Next, let’s look at the Bhattacharyya distance (BD), which is defined as: BD(p(\boldsymbol{x}), q(\boldsymbol{x})) = -\log \int \sqrt{p(\boldsymbol{x}) q(\boldsymbol{x})} d\boldsymbol{x} Related to this is the concept of 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}. After expansion, it is essentially equivalent to the Bhattacharyya distance.
Calculation Result
For two normal distributions, their Bhattacharyya distance is: 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) Here \boldsymbol{\Sigma}=\frac{1}{2}(\boldsymbol{\Sigma}_p + \boldsymbol{\Sigma}_q). We can see that the result is symmetric, which is because the definition of Bhattacharyya distance itself is symmetric.
When one of the two is a standard normal distribution, the result does not simplify significantly, so it is not listed separately here.
Derivation Process
According to the definition, the Bhattacharyya distance between two normal distributions is the negative logarithm of the following integral: \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} Let \boldsymbol{y}=\boldsymbol{x}-\boldsymbol{\mu}_p and \boldsymbol{\Delta}=\boldsymbol{\mu}_p - \boldsymbol{\mu}_q. The integral part can be transformed by substitution: \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} Here \boldsymbol{\Sigma}=\frac{1}{2}(\boldsymbol{\Sigma}_p + \boldsymbol{\Sigma}_q). According to the Gaussian integral formula \eqref{eq:g-int} introduced earlier, the integral result is: \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} So finally: \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}
Wasserstein Distance
If readers want to learn more about probability divergences, they can refer to the book "Statistical Inference Based on Divergence Measures." Now we turn to another class of probability metrics—the Wasserstein distance based on optimal transport.
Following the notation in "From Wasserstein Distance and Duality Theory to WGAN", the Wasserstein distance is defined as: \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} According to the materials I found, an analytical solution for the Wasserstein distance between two normal distributions has currently only been derived for \rho=2. Therefore, the following calculations are for \mathcal{W}_2[p,q], and for simplicity, we denote: \mathcal{W}_2^2[p,q] = \left(\mathcal{W}_2[p,q]\right)^2
Calculation Results
Interestingly, there are two different versions of the result for the Wasserstein distance between two normal distributions. Both versions have a certain level of recognition, but I have not seen any materials explicitly stating that the two are equivalent. The two versions come from different papers and have been given different names.
Version 1
First is the relatively widely circulated version. Often, searching for "Wasserstein distance of normal distributions" yields this result: \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} Regarding this result, some readers might be confused as to "why it is not symmetric with respect to p and q." In fact, it is symmetric because: \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} Then we can 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}).
Version 2
The second version, which looks slightly simpler, is: \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} This version is usually called the "Fréchet distance." Generally, this result can only be found by searching for the keyword "Fréchet distance of normal distributions." The FID (Frechet Inception Distance), a commonly used evaluation metric in GANs, is calculated based on this formula. Its symmetry with respect to p and q can be proven similarly, or derived directly from the equivalence discussion below.
Equivalence
Logically, since the second version is more concise, it should be the standard version propagated. Why two different versions still circulate is quite confusing to me. Theoretically, proving the equivalence of the two versions is not difficult. According to the trace identity, we have: \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} Then simply verify that (\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.
Special Case
In particular, if the multiplication of \boldsymbol{\Sigma}_p and \boldsymbol{\Sigma}_q is commutative, it simplifies to a very intuitive form: \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} Why is this very intuitive? Since the parameters of a normal distribution are \boldsymbol{\mu} and \boldsymbol{\Sigma}, comparing the difference between normal distributions is essentially comparing the differences between \boldsymbol{\mu} and \boldsymbol{\Sigma}. According to machine learning habits, an easily conceived metric is the squared error: \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 However, from a physical perspective, this metric is inappropriate because if \boldsymbol{\mu} is considered to have the dimension of length, then \boldsymbol{\Sigma} has the dimension of length squared. Therefore, \Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 and \Vert \boldsymbol{\Sigma}_p - \boldsymbol{\Sigma}_q\Vert_F^2 are two quantities with different dimensions and cannot be added. To make the dimensions consistent, an intuitive idea is to "take the square root" of \boldsymbol{\Sigma} before calculating the squared error, which leads to equation \eqref{eq:w-jiaohuan}.
In particular, when q is the standard normal distribution, the result simplifies to: \mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p\Vert^2 + \Vert \boldsymbol{\Sigma}_p^{1/2} - \boldsymbol{I}\Vert_F^2
Derivation Process 1
Now I will introduce the first proof, mainly referring to the paper "A class of Wasserstein metrics for probability distributions." Additionally, "The distance between two random vectors with given dispersion matrices" also provides a similar proof for reference.
The following derivation process is simplified by me. It is simpler than the proof in the original paper but still inevitably involves a significant amount of linear algebra. We will introduce it in several parts.
Removing the Mean
Without loss of generality, we can consider only distributions p, q with a mean of 0. If the means of p, q are not 0, let the corresponding distributions with mean 0 be \tilde{p}, \tilde{q}. Then we have: \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} This result implies: \mathcal{W}_2^2[p,q]=\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 + \mathcal{W}_2^2[\tilde{p},\tilde{q}] Therefore, we only need to calculate the Wasserstein distance when the means are both zero, and then add \Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2 to get the result for the general case.
Pure Algebra
Now we assume the means of p, q are both 0, and calculate: \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} where \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] constitutes the covariance matrix of the joint distribution \gamma. We know that the covariance matrix is a positive definite symmetric matrix, so from an algebraic perspective, the problem becomes:
Given that \boldsymbol{\Sigma}_{\gamma}= \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix} is a positive definite symmetric matrix, find the maximum value of \text{Tr}(\boldsymbol{C}).
Schur Complement
To this end, we need to use the following identity regarding the "Schur complement": \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} 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}. For it to be positive definite, \boldsymbol{A} must be positive definite. Since \boldsymbol{\Sigma}_p is already positive definite, \boldsymbol{S} must also be positive definite.
Parameter Separation
Attempt 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, rearrange to get \boldsymbol{\Sigma}_q - \boldsymbol{S} = \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}. Since \boldsymbol{\Sigma}_p is positive definite symmetric, \boldsymbol{\Sigma}_p^{-1} is also positive definite symmetric, and thus \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C} is positive definite symmetric. Therefore, it has a positive definite symmetric square root, i.e., there exists a positive definite symmetric matrix \boldsymbol{R} such that: \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} This indicates that \boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{C}\boldsymbol{R}^{-1} is an orthogonal matrix, denoted as \boldsymbol{O}, so \boldsymbol{C} = \boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R}.
Lagrange Multipliers
At this point, the variables are \boldsymbol{O} and \boldsymbol{R}. We seek the maximum value of \text{Tr}(\boldsymbol{C})=\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R}). We first fix \boldsymbol{R} and find the \boldsymbol{O} that maximizes the value. This is equivalent to maximizing \text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R}) under the constraint \boldsymbol{O}^{\top}\boldsymbol{O}=\boldsymbol{I}. We use the "Lagrange multiplier method": introduce a new parameter matrix \boldsymbol{W} and transform it into the following unconstrained extremum problem: 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})) Taking the derivative: \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} First, note that \boldsymbol{O}^{\top}\boldsymbol{O} - \boldsymbol{I} is symmetric, so the corresponding parameter matrix \boldsymbol{W} is also symmetric. Thus we have: \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} 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: \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})
Inequality
Finally, \boldsymbol{R} needs to be determined. Recalling 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 value is obtained when \boldsymbol{S}=\boldsymbol{0}, and this is indeed the case. This can be seen as a direct corollary of "Weyl’s inequality."
According to Weyl’s inequality, if matrices \boldsymbol{A}, \boldsymbol{B} are both positive definite symmetric matrices, and the eigenvalues of \boldsymbol{A}, \boldsymbol{B}, \boldsymbol{A}+\boldsymbol{B} arranged from smallest to largest 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)}, 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, one-to-one, greater than their respective eigenvalues.
With this conclusion, it’s 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. Then 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 a positive definite symmetric matrix (symmetry is obvious, and since it can have a square root, it is 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). Therefore, their eigenvalues do not exceed the eigenvalues of their sum—which is \boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2}. Thus, the maximum value of each eigenvalue of (\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2})^{1/2} (and thus the maximum value of the trace) is reached at \boldsymbol{S}=\boldsymbol{0}.
As for the proof of Weyl’s inequality, it mainly utilizes the Rayleigh quotient and the Courant–Fischer theorem. Interested readers can look up these two parts of the material and then check the proof of Weyl’s inequality. In fact, after becoming familiar with these two parts, Weyl’s inequality basically follows naturally.
Derivation Process 2
Here I will introduce another simpler proof, the original of which can be found in "The Fréchet distance between multivariate normal distributions." Compared to the first proof, this one is more direct, especially since it does not require as much pure linear algebra knowledge. The following derivation process is further simplified by me, making it easier to understand than the original paper.
In this derivation process, the two steps "Removing the Mean" and "Pure Algebra" are the same as in "Derivation Process 1" and will not be repeated. Thus, the problem has been transformed into:
Given that \boldsymbol{\Sigma}_{\gamma}= \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix} is a positive definite symmetric matrix, find the maximum value of \text{Tr}(\boldsymbol{C}).
Block Matrix
Since \boldsymbol{\Sigma}_{\gamma} is a positive definite symmetric matrix, it must be expressible in the form \boldsymbol{D}\boldsymbol{D}^{\top}. We express \boldsymbol{D} as a block matrix \begin{pmatrix}\boldsymbol{A} \\ \boldsymbol{B}\end{pmatrix}, where \boldsymbol{A}, \boldsymbol{B} \in \mathbb{R}^{n \times 2n}. At this point: \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} Correspondingly, we have \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 parameterization, the problem transforms into:
Given \boldsymbol{A}\boldsymbol{A}^{\top}=\boldsymbol{\Sigma}_p and \boldsymbol{B}\boldsymbol{B}^{\top}=\boldsymbol{\Sigma}_q, find the maximum value of \text{Tr}(\boldsymbol{A}\boldsymbol{B}^{\top}).
This is a constrained maximum value problem. We use the "Lagrange multiplier method": introduce new parameter matrices \boldsymbol{W}_p, \boldsymbol{W}_q and transform it into the following unconstrained extremum problem: 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)) Taking the derivative: \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} Note that \boldsymbol{A}\boldsymbol{A}^{\top} - \boldsymbol{\Sigma}_p and \boldsymbol{B}\boldsymbol{B}^{\top} - \boldsymbol{\Sigma}_q are both symmetric, so the corresponding parameter matrices \boldsymbol{W}_p, \boldsymbol{W}_q are also symmetric. At this point: \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 Let \boldsymbol{W}_p=\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{R}\boldsymbol{\Sigma}_p^{-1/2}. Substituting this into the above equation yields \boldsymbol{\Sigma}_q=\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{R}^2\boldsymbol{\Sigma}_p^{-1/2}, i.e., \boldsymbol{R} = (\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2} And: \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}
Summary
This article has calculated in detail the KL divergence, Bhattacharyya distance, and Wasserstein distance between two multivariate normal distributions, providing their explicit analytical solutions. These results can be used as regularization terms for latent variables in certain scenarios to regulate their distribution. Furthermore, this article can serve as a challenging set of linear algebra exercises for reference and practice.
When reposting, please include the original address of this article: https://kexue.fm/archives/8512
For more detailed reposting matters, please refer to: "Scientific Space FAQ"