两个多元正态分布的KL散度、巴氏距离和W距离 (opens new window)

最佳排版请看原链接：

两个多元正态分布的KL散度、巴氏距离和W距离 - 科学空间|Scientific Spaces (opens new window)

正态分布是最常见的连续型概率分布之一。它是给定均值和协方差后的最大熵分布（参考**《“熵”不起：从熵、最大熵原理到最大熵模型（二）》 (opens new window)**），也可以看作任意连续型分布的二阶近似，它的地位就相当于一般函数的线性近似。从这个角度来看，正态分布算得上是最简单的连续型分布了。也正因为简单，所以对于很多估计量来说，它都能写出解析解来。

本文主要来计算两个多元正态分布的几种度量，包括KL散度、巴氏距离和W距离，它们都有显式解析解。

正态分布

这里简单回顾一下正态分布的一些基础知识。注意，仅仅是回顾，这还不足以作为正态分布的入门教程。

概率密度

正态分布，也即高斯分布，是定义在${\mathbb{R}}^n$ 的连续型概率分布，其概率密度函数为

$$\text{p(boldsymbol{x})=frac{1}{sqrt{(2pi)^n det(boldsymbol{Sigma})}}expleft{-frac{1}{2}(boldsymbol{x}-boldsymbol{mu})^{top}boldsymbol{Sigma}^{-1}(boldsymbol{x}-boldsymbol{mu})right}\$\$ 这里的\$boldsymbol{x},boldsymbol{mu}inmathbb{R}^n\$，\$boldsymbol{mu}\$ 均值向量（本文的向量默认情况下都为列向量），而\$boldsymbol{Sigma}inmathbb{R}^{ntimes n}\$ 为协方差矩阵，它要求是正定对称的。可以看到，正态分布由\$boldsymbol{mu}\$ \$boldsymbol{Sigma}\$ 一确定，因此不难想象它的统计量都是\$boldsymbol{mu}\$ \$boldsymbol{Sigma}\$ 函数。当\$boldsymbol{mu}=boldsymbol{0}, boldsymbol{Sigma}=boldsymbol{I}\$ ，对应的分布称为“标准正态分布”。 \#\#\# 基本性质 通常来说，基本的统计量就是均值和方差了，它们对应着正态分布的两个参数： \$\$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}}$$

由此也可以推出二阶矩的结果：

$${\mathbb{E}}_{\mathit{x}}\left[\mathit{x}{\mathit{x}}^{\mathrm{\top }}\right]=\mathit{\mu }{\mathit{\mu }}^{\mathrm{\top }}+{\mathbb{E}}_{\mathit{x}}[(\mathit{x}-\mathit{\mu })(\mathit{x}-\mathit{\mu }{)}^{\mathrm{\top }}]=\mathit{\mu }{\mathit{\mu }}^{\mathrm{\top }}+\mathbf{\Sigma }$$

还有一个常用的统计量是它的熵：

$$\mathcal{H}={\mathbb{E}}_{\mathit{x}}[-\log p(\mathit{x})]=\frac{n}{2}(1+\log 2\pi )+\frac{1}{2}\log det(\mathbf{\Sigma })$$

其计算过程可以参考后面KL散度的推导。

高斯积分

概率密度函数意味着$\int p(\mathit{x})d\mathit{x}=1$，这就可以推出：

$$\begin{array}{rl}\sqrt{(2\pi {)}^{n}det(\mathbf{\Sigma })}=& \int \exp \{-\frac{1}{2}(\mathit{x}-\mathit{\mu }{)}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}(\mathit{x}-\mathit{\mu })\}d\mathit{x}\\ =& \int \exp \{-\frac{1}{2}{\mathit{x}}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\mathit{x}+{\mathit{\mu }}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\mathit{x}-\frac{1}{2}{\mathit{\mu }}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\mathit{\mu }\}d\mathit{x}\end{array}$$

设$\mathit{\omega }={\mathbf{\Sigma }}^{-1}\mathit{\mu }$，那么得到高斯积分

$$\int \exp \{-\frac{1}{2}{\mathit{x}}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}\mathit{x}+{\mathit{\omega }}^{\mathrm{\top }}\mathit{x}\}d\mathit{x}=\sqrt{(2\pi {)}^{n}det(\mathbf{\Sigma })}\exp \left\{\frac{1}{2}{\mathit{\omega }}^{\mathrm{\top }}\mathbf{\Sigma }\mathit{\omega }\right\}$$

利用它我们可以算出正态分布的特征函数

$${\mathbb{E}}_{\mathit{x}}[\exp \left({\mathit{\omega }}^{\mathrm{\top }}\mathit{x}\right)]=\exp ({\mathit{\omega }}^{\mathrm{\top }}\mathit{\mu }+\frac{1}{2}{\mathit{\omega }}^{\mathrm{\top }}\mathbf{\Sigma }\mathit{\omega })$$

特征函数可以用来算正态分布的各阶矩。

线性代数

这里补充一些线性代数基础，它们在后面的推导中会频繁用到。同样地，这仅仅是“回顾”，并不能作为线性代数教程。

内积范数

首先，我们来定义内积和范数。对于向量$\mathit{x}=(x_1,\cdots ,x_n)$ $\mathit{y}=(y_1,\cdots ,y_n)$，内积按照

$$⟨\mathit{x},\mathit{y}⟩=\sum _{i=1}^n{x}_i{y}_i$$

来定义，而模长定义为$\Vert \mathit{x}\Vert =\sqrt{⟨\mathit{x},\mathit{y}⟩}$。对于$m\times n$ 矩阵$\mathit{A}=(a_{i,j}),\mathit{B}=(b_{i,j})$，我们按照类似的方式定义：

$$⟨\mathit{A},\mathit{B}{⟩}_F=\sum _{i=1}^m\sum _{j=1}^n{a}_{i,j}{b}_{i,j}$$

这称为Frobenius内积，对应的$\Vert \mathit{A}{\Vert }_F=\sqrt{⟨\mathit{A},\mathit{B}{⟩}_F}$ 为Frobenius范数。不难看到，Frobenius内积和范数，事实上就是把矩阵展平为向量后，当作常规的向量来运算。

关于Frobenius内积，最关键的性质之一是成立恒等式

$$⟨\mathit{A},\mathit{B}{⟩}_F=\text{Tr}\left({\mathit{A}}^{\mathrm{\top }}\mathit{B}\right)=\text{Tr}\left(\mathit{B}{\mathit{A}}^{\mathrm{\top }}\right)=\text{Tr}\left(\mathit{A}{\mathit{B}}^{\mathrm{\top }}\right)=\text{Tr}\left({\mathit{B}}^{\mathrm{\top }}\mathit{A}\right)$$

也就是说，矩阵的Frobenius内积可以转化为矩阵乘法的迹，并且交换相乘顺序不改变结果（不改变迹的结果，但是矩阵乘法的整体结果会改变）。

正定对称

接着，来看正定对称矩阵的一些性质。$\mathbf{\Sigma }$ 一个正定对称矩阵，对称说的是${\mathbf{\Sigma }}^{\mathrm{\top }}=\mathbf{\Sigma }$，正定说的是对于任意非零向量$\mathit{\xi }\in {\mathbb{R}}^n$，都有${\mathit{\xi }}^{\mathrm{\top }}\mathbf{\Sigma }\mathit{\xi }>0$。可以证明，如果${\mathbf{\Sigma }}_1,{\mathbf{\Sigma }}_2$ 是正定对称矩阵，那么${\mathbf{\Sigma }}_1^{-1},{\mathbf{\Sigma }}_2^{-1},{\mathbf{\Sigma }}_1+{\mathbf{\Sigma }}_2$ 都是正定对称矩阵。如果$\mathit{C}={\mathit{B}}^{\mathrm{\top }}\mathit{A}\mathit{B}$，$\mathit{B}$ 可逆阵，那么$\mathit{C}$ 正定对称的当且仅当$\mathit{A}$ 正定对称的。

此外还有半正定的概念，指对于任意非零向量$\mathit{\xi }\in {\mathbb{R}}^n$，都有${\mathit{\xi }}^{\mathrm{\top }}\mathbf{\Sigma }\mathit{\xi }\ge 0$，也就是说可能存在非零向量$\mathit{\xi }$ 得${\mathit{\xi }}^{\mathrm{\top }}\mathbf{\Sigma }\mathit{\xi }=0$。不过考虑到正定矩阵在半正定矩阵中稠密，所以我们不严格区分正定和半正定了，统一按照正定矩阵来处理。

正定对称矩阵有一个重要的性质，那就是它的SVD分解跟特征值分解一致，即具有下述形式的分解

$$\mathbf{\Sigma }=\mathit{U}\mathbf{\Lambda }{\mathit{U}}^{\mathrm{\top }}$$

其中$\mathit{U}$ 正交矩阵，而$\mathbf{\Lambda }$ 对角阵，并且对角线上的元素都是正的。该结果的一个直接推论是：正定对称矩阵都可以“开平方”，其平方根为${\mathbf{\Sigma }}^{1/2}=\mathit{U}{\mathbf{\Lambda }}^{1/2}{\mathit{U}}^{\mathrm{\top }}$，其中${\mathbf{\Lambda }}^{1/2}$ 指将对角线上的元素都开平方，可以检验平方根矩阵也是正定对称的。反过来，可以开平方的对称矩阵，一定也是正定对称矩阵。

矩阵求导

最后，在求Wasserstein距离的时候，还需要用到一些矩阵求导公式，如果不了解的读者，可以直接参考维基百科的“Matrix Calculus (opens new window)”。当然，其实也不难，主要用到了

$$\frac{\partial \text{Tr}\left(\mathit{X}\mathit{A}\right)}{\partial \mathit{X}}=\mathit{A}$$

剩下的可以结合迹的运算公式来派生出来，比如

$$\frac{\partial \text{Tr}\left(\mathit{A}\mathit{X}\mathit{B}\right)}{\partial \mathit{X}}=\frac{\partial \text{Tr}\left(\mathit{X}\mathit{B}\mathit{A}\right)}{\partial \mathit{X}}=\mathit{B}\mathit{A}$$

KL散度

作为第一个尝试，我们来算两个高斯分布的**KL散度（Kullback-​Leibler divergence） (opens new window)**。KL散度算是最常用的分布度量之一了，因为它积分之前需要取对数，这对于指数簇分布来说通常能得到相对简单的结果。此外它还跟“熵”有着比较紧密的联系。

计算结果

两个概率分布的KL散度定义为

$$KL(p(\mathit{x})\Vert q(\mathit{x}))={\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[\log \frac{p(\mathit{x})}{q(\mathit{x})}]={\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[\log p(\mathit{x})]+{\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[-\log q(\mathit{x})]$$

对于两个正态分布来说，计算结果是

$$KL(p(\mathit{x})\Vert q(\mathit{x}))=\frac{1}{2}[({\mathit{\mu }}_p-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}({\mathit{\mu }}_p-{\mathit{\mu }}_q)-\log det({\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_p)+\text{Tr}\left({\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_p\right)-n]$$

特别地，当$q$ 标准正态分布时，结果简化为

$$KL(p(\mathit{x})\Vert q(\mathit{x}))=\frac{1}{2}[\Vert {\mathit{\mu }}_p{\Vert }^2-\log det({\mathbf{\Sigma }}_p)+\text{Tr}({\mathbf{\Sigma }}_p)-n]$$

推导过程

从KL散度的定义知道，我们主要把${\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[-\log q(\mathit{x})]$ 出来就行了：

$$\begin{array}{rl}{\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[-\log q(\mathit{x})]=& {\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[\frac{n}{2}\log 2\pi +\frac{1}{2}\log det({\mathrm{\Sigma }}_q)+\frac{1}{2}(\mathit{x}-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}(\mathit{x}-{\mathit{\mu }}_q)]\\ =& \frac{n}{2}\log 2\pi +\frac{1}{2}\log det({\mathbf{\Sigma }}_q)+\frac{1}{2}{\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[(\mathit{x}-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}(\mathit{x}-{\mathit{\mu }}_q)]\end{array}$$

现在，关于迹的恒等式就可以派上用场了：

$$\begin{array}{rl}{\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[(\mathit{x}-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}(\mathit{x}-{\mathit{\mu }}_q)]=& {\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}\left[\text{Tr}((\mathit{x}-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}(\mathit{x}-{\mathit{\mu }}_q)\right]\\ =& {\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[\text{Tr}({\mathbf{\Sigma }}_q^{-1}(\mathit{x}-{\mathit{\mu }}_q)(\mathit{x}-{\mathit{\mu }}_q{)}^{\mathrm{\top }}]\\ =& \text{Tr}\left({\mathbf{\Sigma }}_q^{-1}{\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[(\mathit{x}-{\mathit{\mu }}_q)(\mathit{x}-{\mathit{\mu }}_q{)}^{\mathrm{\top }}]\right)\\ =& \text{Tr}\left({\mathbf{\Sigma }}_q^{-1}{\mathbb{E}}_{\mathit{x}\sim p(\mathit{x})}[\mathit{x}{\mathit{x}}^{\mathrm{\top }}-{\mathit{\mu }}_q{\mathit{x}}^{\mathrm{\top }}-\mathit{x}{\mathit{\mu }}_q^{\mathrm{\top }}+{\mathit{\mu }}_q{\mathit{\mu }}_q^{\mathrm{\top }}]\right)\\ =& \text{Tr}\left({\mathbf{\Sigma }}_q^{-1}({\mathbf{\Sigma }}_p+{\mathit{\mu }}_p{\mathit{\mu }}_p^{\mathrm{\top }}-{\mathit{\mu }}_q{\mathit{\mu }}_p^{\mathrm{\top }}-{\mathit{\mu }}_p{\mathit{\mu }}_q^{\mathrm{\top }}+{\mathit{\mu }}_q{\mathit{\mu }}_q^{\mathrm{\top }})\right)\\ =& \text{Tr}({\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_p+{\mathbf{\Sigma }}_q^{-1}({\mathit{\mu }}_p-{\mathit{\mu }}_q)({\mathit{\mu }}_p-{\mathit{\mu }}_q{)}^{\mathrm{\top }})\\ =& \text{Tr}\left({\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_p\right)+({\mathit{\mu }}_p-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}({\mathit{\mu }}_p-{\mathit{\mu }}_q)\end{array}$$

注意${\mathit{\mu }}_q={\mathit{\mu }}_p,{\mathbf{\Sigma }}_q={\mathbf{\Sigma }}_p$ ，上式就等于$n$，此时就对应正态分布的熵。所以最终得到

$$\begin{array}{rl}KL(p(\mathit{x})\Vert q(\mathit{x}))=& \frac{1}{2}[n\log 2\pi +\log det({\mathbf{\Sigma }}_q)+\text{Tr}\left({\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_p\right)+({\mathit{\mu }}_p-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}({\mathit{\mu }}_p-{\mathit{\mu }}_q)]\\ & -\frac{1}{2}[n\log 2\pi +\log det({\mathbf{\Sigma }}_p)+n]\\ =& \frac{1}{2}[({\mathit{\mu }}_p-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}_q^{-1}({\mathit{\mu }}_p-{\mathit{\mu }}_q)-\log det({\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_p)+\text{Tr}\left({\mathbf{\Sigma }}_q^{-1}{\mathbf{\Sigma }}_p\right)-n]\end{array}$$

巴氏距离

然后，我们来看看**巴氏距离（Bhattacharyya distance） (opens new window)**，它定义为

$$BD(p(\mathit{x}),q(\mathit{x}))=-\log \int \sqrt{p(\mathit{x})q(\mathit{x})}d\mathit{x}$$

与之相关的还有一个叫做“Hellinger距离 (opens new window)”的概念，它的平方定义为$ \frac{1}{2}\int\left(\sqrt{p(\boldsymbol{x})} - \sqrt{q(\boldsymbol{x})}\right)^2 d\boldsymbol{x}$，展开后就能发现跟巴氏距离本质是等价的。

计算结果

对于两个正态分布来说，它们的巴氏距离为

$$BD(p(\mathit{x}),q(\mathit{x}))=\frac{1}{2}\log \frac{det(\mathbf{\Sigma })}{\sqrt{det({\mathbf{\Sigma }}_p{\mathbf{\Sigma }}_q)}}+\frac{1}{8}({\mathit{\mu }}_p-{\mathit{\mu }}_q{)}^{\mathrm{\top }}{\mathbf{\Sigma }}^{-1}({\mathit{\mu }}_p-{\mathit{\mu }}_q)$$

这里$\mathbf{\Sigma }=\frac{1}{2}({\mathbf{\Sigma }}_p+{\mathbf{\Sigma }}_q)$。可以看到结果是对称的，这是因为巴氏距离的定义本身就是对称的。

当两者之一为标准正态分布时，结果并没有明显简化，所以这里就不单独写出来了。

推导过程

按照定义，两个正态分布的巴氏距离，是下述积分的负对数： $$\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}\$$ 记$\mathit{y}=\mathit{x}-{\mathit{\mu }}_p,\mathbf{\Delta }={\mathit{\mu }}_p-{\mathit{\mu }}_q$，积分部分可以换元为 $$\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}\$$ 这里$\mathbf{\Sigma }=\frac{1}{2}({\mathbf{\Sigma }}_p+{\mathbf{\Sigma }}_q)$。按照前面介绍的高斯积分公式$\text{(???)}$，积分结果是 $$\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}\$$ 所以最终 $$\begin{aligned} BD(p(\boldsymbol{x}), q(\boldsymbol{x})) =&, -\log \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} \ =&, -\log \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}\ =&,\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} \$$

W距离

如果读者还想看了解更多关于概率散度的内容，可以参考书籍**《Statistical Inference Based on Divergence Measures》 (opens new window)**。现在我们转向另一类概率度量——基于最优传输的W距离（Wasserstein距离）。

沿用**《从Wasserstein距离、对偶理论到WGAN》 (opens new window)**中的记号，W距离的定义如下：

$$\mathcal{W}[p,q]=\underset{\gamma \in \mathrm{\Pi }[p,q]}{inf}\iint \gamma (\mathit{x},\mathit{y})d(\mathit{x},\mathit{y})d\mathit{x}d\mathit{y}=\underset{\gamma \in \mathrm{\Pi }[p,q]}{inf}{\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \gamma (\mathit{x},\mathit{y})}[d(\mathit{x},\mathit{y})]$$

不同的$d(\mathit{x},\mathit{y})$ 得到不同的结果，为了得到较为简单的解，这里选择

$$d(\mathit{x},\mathit{y})=\Vert \mathit{x}-\mathit{y}{\Vert }^2$$

计算结果

有意思的是，关于两个正态分布的W距离的结果，流传着两个不同的版本，这两个版本都有一定的认知度，但却没有看到有明确说两者等价的资料。两个版本出自不同的论文，还被冠以了不同的名字。

版本1

首先第一个流传相对较广的版本（很多文献包括维基百科也使用这个版本）：

$$\begin{gathered} \mathcal{W}[p,q]=\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2+\text{Tr}({\mathbf{\Sigma }}_p)+\text{Tr}({\mathbf{\Sigma }}_q)-2\text{Tr}(({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}) \\ \$\$\mathrm{于}\mathrm{这}\mathrm{个}\mathrm{结}\mathrm{果}，\mathrm{有}\mathrm{的}\mathrm{读}\mathrm{者}\mathrm{可}\mathrm{能}\mathrm{困}\mathrm{惑}\mathrm{于}“\mathrm{怎}\mathrm{么}\mathrm{关}\mathrm{于}\$p,q\$\mathrm{是}\mathrm{对}\mathrm{称}\mathrm{的}”，\mathrm{事}\mathrm{实}\mathrm{上}，\mathrm{它}\mathrm{关}\mathrm{于}\$p,q\$\mathrm{对}\mathrm{称}\mathrm{的}，\mathrm{因}\mathrm{为}\$\$\begin{array}{rl}\text{Tr}(({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2})=& \text{Tr}(({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}{\mathbf{\Sigma }}_p^{-1/2}{\mathbf{\Sigma }}_q^{-1/2}{\mathbf{\Sigma }}_q^{1/2}{\mathbf{\Sigma }}_p^{1/2})\\ =& \text{Tr}({\mathbf{\Sigma }}_q^{1/2}{\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}{\mathbf{\Sigma }}_p^{-1/2}{\mathbf{\Sigma }}_q^{-1/2})\end{array} \end{gathered}$$

然后我们可以直接验证$({\mathbf{\Sigma }}_q^{1/2}{\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}{\mathbf{\Sigma }}_p^{-1/2}{\mathbf{\Sigma }}_q^{-1/2}{)}^2={\mathbf{\Sigma }}_q^{1/2}{\mathbf{\Sigma }}_p{\mathbf{\Sigma }}_q^{1/2}$，所以有$\text{Tr}(({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2})=\text{Tr}(({\mathbf{\Sigma }}_q^{1/2}{\mathbf{\Sigma }}_p{\mathbf{\Sigma }}_q^{1/2}{)}^{1/2})$。

版本2

第二个看起来稍微简单一些的版本，结果是：

$$\mathcal{W}[p,q]=\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2+\text{Tr}({\mathbf{\Sigma }}_p)+\text{Tr}({\mathbf{\Sigma }}_q)-2\text{Tr}(({\mathbf{\Sigma }}_p{\mathbf{\Sigma }}_q{)}^{1/2})$$

这个版本通常被称为“Fréchet距离 (opens new window)”。GAN中经常使用的评价指标FID（Frechet Inception Distance (opens new window)），就是基于这个公式进行计算的。可以模仿前面证明它关于$p,q$ 对称性，当然也可以从下面的等价性讨论中直接得出。

等价性

事实上，证明两者的等价性并不难：

$$\begin{array}{rl}\text{Tr}(({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2})=& \text{Tr}(({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}{\mathbf{\Sigma }}_p^{-1/2}{\mathbf{\Sigma }}_p^{1/2})\\ =& \text{Tr}({\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}{\mathbf{\Sigma }}_p^{-1/2})\end{array}$$

然后直接验证$({\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}{\mathbf{\Sigma }}_p^{-1/2}{)}^2={\mathbf{\Sigma }}_p{\mathbf{\Sigma }}_q$ 可。

特殊时

特别地，如果${\mathbf{\Sigma }}_p,{\mathbf{\Sigma }}_q$ 乘法可以交换，那么将会简化为非常直观的形式：

$$\mathcal{W}[p,q]=\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2+\Vert {\mathbf{\Sigma }}_p^{1/2}-{\mathbf{\Sigma }}_q^{1/2}{\Vert }_F^2$$

为什么说它非常直观呢？因为正态分布的参数为$\mathit{\mu },\mathbf{\Sigma }$，所以比较正态分布的差异其实就是比较$\mathit{\mu },\mathbf{\Sigma }$ 差异，按照机器学习的习惯，一个很容易相当想到的指标是平方误差

$$\mathcal{W}[p,q]=\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2+\Vert {\mathbf{\Sigma }}_p-{\mathbf{\Sigma }}_q{\Vert }_F^2$$

但从物理角度来看，这个指标是不妥的，因为如果将$\mathit{\mu }$ 成是长度量纲，那么$\mathbf{\Sigma }$ 具有长度平方的量纲，所以$\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2$ $\Vert {\mathbf{\Sigma }}_p-{\mathbf{\Sigma }}_q{\Vert }_F^2$ 具有不同量纲的两个量，不能相加。而为了使得量纲一致，直观的想法就是把$\mathbf{\Sigma }$“开平方”后再算平方误差，这就得到了式$\text{(???)}$。

特别地，当$q$ 标准正态分布时，结果简化为

$$\mathcal{W}[p,q]=\Vert {\mathit{\mu }}_p{\Vert }^2+\Vert {\mathbf{\Sigma }}_p^{1/2}-\mathit{I}{\Vert }_F^2$$

推导过程1

现在介绍第一个证明，主要参考了论文**《A class of Wasserstein metrics for probability distributions》 (opens new window)。另外《The distance between two random vectors with given dispersion matrices》 (opens new window)**也提供了一个类似的证明，也可以参考。

下面的推导过程则是经过笔者简化的，相对原论文的证明来说简单一些，但依然不可避免地会涉及到较多的线性代数知识，我们将分几个部分介绍。

去均值

不失一般性，我们可以只考虑均值为0的分布$p,q$。因为如果$p,q$ 均值不为0，那么设对应的均值为0的分布为$\tilde{p},\tilde{q}$，此时有

$$\begin{array}{rl}& {\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \gamma (\mathit{x},\mathit{y})}[\Vert \mathit{x}-\mathit{y}{\Vert }^2]\\ =& {\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \tilde{\gamma }(\mathit{x},\mathit{y})}[\Vert (\mathit{x}+{\mathit{\mu }}_p)-(\mathit{y}+{\mathit{\mu }}_q){\Vert }^2]\\ =& {\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \tilde{\gamma }(\mathit{x},\mathit{y})}[\Vert \mathit{x}-\mathit{y}{\Vert }^2+\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2+2⟨\mathit{x}-\mathit{y},{\mathit{\mu }}_p-{\mathit{\mu }}_q⟩]\\ =& \Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2+{\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \tilde{\gamma }(\mathit{x},\mathit{y})}[\Vert \mathit{x}-\mathit{y}{\Vert }^2]\end{array}$$

该结果意味着

$$\mathcal{W}[p,q]=\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2+\mathcal{W}[\tilde{p},\tilde{q}]$$

所以，只需要算出均值都为零时的Wasserstein距离，然后加上$\Vert {\mathit{\mu }}_p-{\mathit{\mu }}_q{\Vert }^2$ 得到了一般情况的结果。

纯代数

现在我们假设$p,q$ 均值均为0，然后计算

$$\begin{array}{rl}{\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \gamma (\mathit{x},\mathit{y})}[\Vert \mathit{x}-\mathit{y}{\Vert }^2]=& {\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \gamma (\mathit{x},\mathit{y})}[{\mathit{x}}^{\mathrm{\top }}\mathit{x}+{\mathit{y}}^{\mathrm{\top }}\mathit{y}-2{\mathit{y}}^{\mathrm{\top }}\mathit{x}]\\ =& {\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \gamma (\mathit{x},\mathit{y})}\left[\text{Tr}(\mathit{x}{\mathit{x}}^{\mathrm{\top }}+\mathit{y}{\mathit{y}}^{\mathrm{\top }}-2\mathit{x}{\mathit{y}}^{\mathrm{\top }})\right]\\ =& \text{Tr}\left({\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \gamma (\mathit{x},\mathit{y})}[\mathit{x}{\mathit{x}}^{\mathrm{\top }}+\mathit{y}{\mathit{y}}^{\mathrm{\top }}-2\mathit{x}{\mathit{y}}^{\mathrm{\top }}]\right)\\ =& \text{Tr}({\mathbf{\Sigma }}_p)+\text{Tr}({\mathbf{\Sigma }}_q)-2\text{Tr}(\mathit{C})\end{array}$$

其中

$${\mathbf{\Sigma }}_{\gamma }=\left(\begin{array}{cc}{\mathbf{\Sigma }}_p& \mathit{C}\\ {\mathit{C}}^{\mathrm{\top }}& {\mathbf{\Sigma }}_q\end{array}\right)={\mathbb{E}}_{(\mathit{x},\mathit{y})\sim \gamma (\mathit{x},\mathit{y})}\left[\left(\begin{array}{c}\mathit{x}\\ \mathit{y}\end{array}\right)\left(\begin{array}{cc}{\mathit{x}}^{\mathrm{\top }}& {\mathit{y}}^{\mathrm{\top }}\end{array}\right)\right]$$

构成联合分布$\gamma$ 协方差矩阵。我们知道协方差矩阵是正定对阵矩阵，所以从代数的角度看，问题变成了： 已知${\mathbf{\Sigma }}_{\gamma }=\left(\begin{array}{cc}{\mathbf{\Sigma }}_p& \mathit{C}\\ {\mathit{C}}^{\mathrm{\top }}& {\mathbf{\Sigma }}_q\end{array}\right)$ 正定对称矩阵，求$\text{Tr}(\mathit{C})$ 最大值。

舒尔补

为此，我们需要利用下述关于“舒尔补”的恒等式：

$$\left(\begin{array}{cc}{\mathbf{\Sigma }}_p& \mathit{C}\\ {\mathit{C}}^{\mathrm{\top }}& {\mathbf{\Sigma }}_q\end{array}\right)=\left(\begin{array}{cc}\mathit{I}& \mathbf{0}\\ {\mathit{C}}^{\mathrm{\top }}{\mathbf{\Sigma }}_p^{-1}& \mathit{I}\end{array}\right)\left(\begin{array}{cc}{\mathbf{\Sigma }}_p& \mathbf{0}\\ \mathbf{0}& {\mathbf{\Sigma }}_q-{\mathit{C}}^{\mathrm{\top }}{\mathbf{\Sigma }}_p^{-1}\mathit{C}\end{array}\right)\left(\begin{array}{cc}\mathit{I}& {\mathbf{\Sigma }}_p^{-1}\mathit{C}\\ \mathbf{0}& \mathit{I}\end{array}\right)$$

其中对称矩阵$\mathit{S}={\mathbf{\Sigma }}_q-{\mathit{C}}^{\mathrm{\top }}{\mathbf{\Sigma }}_p^{-1}\mathit{C}$ 为“舒尔补（Schur Complement） (opens new window)”，该分解具有${\mathit{B}}^{\mathrm{\top }}\mathit{A}\mathit{B}$ 形式，要想它是正定的，那么$\mathit{A}$ 是正定的，而${\mathbf{\Sigma }}_p$ 经是正定的，所以$\mathit{S}$ 要是正定的。

分参数

我们尝试分离参数，即从$\mathit{S}={\mathbf{\Sigma }}_q-{\mathit{C}}^{\mathrm{\top }}{\mathbf{\Sigma }}_p^{-1}\mathit{C}$ 把$\mathit{C}$ 出来。首先移项得到${\mathbf{\Sigma }}_q-\mathit{S}={\mathit{C}}^{\mathrm{\top }}{\mathbf{\Sigma }}_p^{-1}\mathit{C}$，由于${\mathbf{\Sigma }}_p$ 正定对称的，所以${\mathbf{\Sigma }}_p^{-1}$ 是，从而${\mathit{C}}^{\mathrm{\top }}{\mathbf{\Sigma }}_p^{-1}\mathit{C}$ 是正定对称的，那么它具有正定对称的平方根，即存在正定对称矩阵$\mathit{R}$，使得

$${\mathit{C}}^{\mathrm{\top }}{\mathbf{\Sigma }}_p^{-1}\mathit{C}={\mathit{R}}^2\iff {\left({\mathbf{\Sigma }}_p^{-1/2}\mathit{C}{\mathit{R}}^{-1}\right)}^{\mathrm{\top }}\left({\mathbf{\Sigma }}_p^{-1/2}\mathit{C}{\mathit{R}}^{-1}\right)=\mathit{I}$$

这说明${\mathbf{\Sigma }}_p^{-1/2}\mathit{C}{\mathit{R}}^{-1}$ 正交矩阵，记为$\mathit{O}$，那么$\mathit{C}={\mathbf{\Sigma }}_p^{1/2}\mathit{O}\mathit{R}$。

乘子法

此时，变量分别是$\mathit{O}$ $\mathit{R}$，求$\text{Tr}(\mathit{C})=\text{Tr}({\mathbf{\Sigma }}_p^{1/2}\mathit{O}\mathit{R})$ 最大值。我们先固定$\mathit{R}$，求取最大值时的$\mathit{O}$，此时相当于在${\mathit{O}}^{\mathrm{\top }}\mathit{O}=\mathit{I}$ 约束下，求$\text{Tr}({\mathbf{\Sigma }}_p^{1/2}\mathit{O}\mathit{R})$ 最大值，我们用“拉格朗日乘子法 (opens new window)”：引入新参数矩阵$\mathit{W}$，转化为下述无约束极值问题

$$F=\text{Tr}({\mathbf{\Sigma }}_p^{1/2}\mathit{O}\mathit{R})-\text{Tr}(\mathit{W}({\mathit{O}}^{\mathrm{\top }}\mathit{O}-\mathit{I}))$$

求导：

$$\begin{array}{rl}& \frac{\partial F}{\partial \mathit{O}}=\mathbf{0}\Rightarrow \mathit{R}{\mathbf{\Sigma }}_p^{1/2}=\mathit{W}{\mathit{O}}^{\mathrm{\top }}\\ & \frac{\partial F}{\partial \mathit{W}}=\mathbf{0}\Rightarrow {\mathit{O}}^{\mathrm{\top }}\mathit{O}=\mathit{I}\end{array}$$

首先留意到${\mathit{O}}^{\mathrm{\top }}\mathit{O}-\mathit{I}$ 对称的，因此对应的参数矩阵$\mathit{W}$ 是对称的，于是我们有：

$${\left(\mathit{O}\mathit{W}{\mathit{O}}^{\mathrm{\top }}\right)}^2={\left(\mathit{W}{\mathit{O}}^{\mathrm{\top }}\right)}^{\mathrm{\top }}\left(\mathit{W}{\mathit{O}}^{\mathrm{\top }}\right)={\mathbf{\Sigma }}_p^{1/2}{\mathit{R}}^2{\mathbf{\Sigma }}_p^{1/2}$$

即$\mathit{O}\mathit{W}{\mathit{O}}^{\mathrm{\top }}=({\mathbf{\Sigma }}_p^{1/2}{\mathit{R}}^2{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}$，所以此时

$$\text{Tr}({\mathbf{\Sigma }}_p^{1/2}\mathit{O}\mathit{R})=\text{Tr}(\mathit{O}\mathit{R}{\mathbf{\Sigma }}_p^{1/2})=\text{Tr}(\mathit{O}\mathit{W}{\mathit{O}}^{\mathrm{\top }})=\text{Tr}(({\mathbf{\Sigma }}_p^{1/2}{\mathit{R}}^2{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2})$$

不等式

最后需要把$\mathit{R}$ 定下来。回顾$\mathit{R}$ 定义，我们有${\mathit{R}}^2={\mathbf{\Sigma }}_q-\mathit{S}$，其中$\mathit{S}$ 正定矩阵。直觉上$\mathit{S}=\mathbf{0}$ 取得最大值，事实上也确实如此，这算是“Weyl不等式 (opens new window)”的一个直接推论。

根据Weyl不等式，如果矩阵$\mathit{A},\mathit{B},\mathit{A}+\mathit{B}$ 是正定对称矩阵，它们的特征值从小到大排列分别为$0\le {\lambda }_1^{(A)}\le \cdots \le {\lambda }_n^{(A)}$、${\lambda }_1^{(B)}\le \cdots \le {\lambda }_n^{(B)}$ $0\le {\lambda }_1^{(A+B)}\le \cdots \le {\lambda }_n^{(A+B)}$，那么对于任意$1\le i\le n$，都有${\lambda }_i^{(A)}\le {\lambda }_i^{(A+B)}$ ${\lambda }_i^{(B)}\le {\lambda }_i^{(A+B)}$，也就是说： 正定对称矩阵的和的特征值，一一对应地大于它们各自的特征值。

有了这个结论，那就简单了，设${({\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_q-\mathit{S}){\mathbf{\Sigma }}_p^{1/2})}^{1/2}$ 特征值为$0\le {\lambda }_1\le \cdots \le {\lambda }_n$，那么它的迹就是${\lambda }_1+\cdots +{\lambda }_n$，对应地，${\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_q-\mathit{S}){\mathbf{\Sigma }}_p^{1/2}$ 特征值为$0\le {\lambda }_1^2\le \cdots \le {\lambda }_n^2$，注意${\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_q-\mathit{S}){\mathbf{\Sigma }}_p^{1/2}$ 正定对称矩阵（对称是显然的，而因为它能开平方，所以正定），${\mathbf{\Sigma }}_p^{1/2}\mathit{S}{\mathbf{\Sigma }}_p^{1/2}$ 是正定对称的（因为$\mathit{S}$ 正定对称的），所以它们的特征值，都不超过它们的和——也就是${\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}$ 特征值，所以说，${({\mathbf{\Sigma }}_p^{1/2}({\mathbf{\Sigma }}_q-\mathit{S}){\mathbf{\Sigma }}_p^{1/2})}^{1/2}$ 个特征值的最大值（也就是迹的最大值），在$\mathit{S}=\mathbf{0}$ 取到。

至于Weyl不等式的证明，主要利用到了**Rayleigh quotient (opens new window)和Courant–Fischer定理 (opens new window)**，有兴趣了解的读者自行查阅这两部分资料后，再查阅Wely不等式的证明就好。事实上，熟悉这两部分内容后，Weyl不等式基本上就“水到渠成”了。

推导过程2

这里继续介绍另一个更为简单的证明，原始证明可以在**《The Fréchet distance between multivariate normal distributions》 (opens new window)**找到。相对而言该证明简单不少，尤其是不需要太多的纯线性代数知识。下面的推导过程依旧是经过笔者进一步简化的，比原始论文更好理解一些。

在这个推导过程中，“去均值”、“纯代数”两个步骤跟“推导过程1”是一样的，不再重复。所以，此时问题已经被转化为 已知${\mathbf{\Sigma }}_{\gamma }=\left(\begin{array}{cc}{\mathbf{\Sigma }}_p& \mathit{C}\\ {\mathit{C}}^{\mathrm{\top }}& {\mathbf{\Sigma }}_q\end{array}\right)$ 正定对称矩阵，求$\text{Tr}(\mathit{C})$ 最大值。

分块阵

由于${\mathbf{\Sigma }}_{\gamma }$ 正定对称矩阵，所以它必然可以表达成$\mathit{D}{\mathit{D}}^{\mathrm{\top }}$ 形式，我们将$\mathit{D}$ 达为分块矩阵$\left(\begin{array}{c}\mathit{A}\\ \mathit{B}\end{array}\right)$，其中$\mathit{A},\mathit{B}\in {\mathbb{R}}^{n\times 2n}$，此时

$$\left(\begin{array}{cc}{\mathbf{\Sigma }}_p& \mathit{C}\\ {\mathit{C}}^{\mathrm{\top }}& {\mathbf{\Sigma }}_q\end{array}\right)=\left(\begin{array}{c}\mathit{A}\\ \mathit{B}\end{array}\right)\left(\begin{array}{cc}{\mathit{A}}^{\mathrm{\top }}& {\mathit{B}}^{\mathrm{\top }}\end{array}\right)=\left(\begin{array}{cc}\mathit{A}{\mathit{A}}^{\mathrm{\top }}& \mathit{A}{\mathit{B}}^{\mathrm{\top }}\\ \mathit{B}{\mathit{A}}^{\mathrm{\top }}& \mathit{B}{\mathit{B}}^{\mathrm{\top }}\end{array}\right)$$

对应地有${\mathbf{\Sigma }}_p=\mathit{A}{\mathit{A}}^{\mathrm{\top }},{\mathbf{\Sigma }}_q=\mathit{B}{\mathit{B}}^{\mathrm{\top }},\mathit{C}=\mathit{A}{\mathit{B}}^{\mathrm{\top }}$。

乘子法

在上述参数化之下，问题转化为： 已知$\mathit{A}{\mathit{A}}^{\mathrm{\top }}={\mathbf{\Sigma }}_p,\mathit{B}{\mathit{B}}^{\mathrm{\top }}={\mathbf{\Sigma }}_q$，求$\text{Tr}(\mathit{A}{\mathit{B}}^{\mathrm{\top }})$ 最大值。

这是一个带约束的最大值问题，我们用“拉格朗日乘子法 (opens new window)”：引入新参数矩阵${\mathit{W}}_p,{\mathit{W}}_q$，转化为下述无约束极值问题

$$F=\text{Tr}(\mathit{A}{\mathit{B}}^{\mathrm{\top }})-\text{Tr}({\mathit{W}}_p(\mathit{A}{\mathit{A}}^{\mathrm{\top }}-{\mathbf{\Sigma }}_p))-\text{Tr}({\mathit{W}}_q(\mathit{B}{\mathit{B}}^{\mathrm{\top }}-{\mathbf{\Sigma }}_q))$$

求导：

$$\begin{array}{rl}& \frac{\partial F}{\partial \mathit{A}}=\mathbf{0}\Rightarrow {\mathit{B}}^{\mathrm{\top }}={\mathit{A}}^{\mathrm{\top }}{\mathit{W}}_p\\ & \frac{\partial F}{\partial \mathit{B}}=\mathbf{0}\Rightarrow {\mathit{A}}^{\mathrm{\top }}={\mathit{B}}^{\mathrm{\top }}{\mathit{W}}_q\\ & \frac{\partial F}{\partial {\mathit{W}}_p}=\mathbf{0}\Rightarrow \mathit{A}{\mathit{A}}^{\mathrm{\top }}={\mathbf{\Sigma }}_p\\ & \frac{\partial F}{\partial {\mathit{W}}_q}=\mathbf{0}\Rightarrow \mathit{B}{\mathit{B}}^{\mathrm{\top }}={\mathbf{\Sigma }}_q\end{array}$$

注意到$\mathit{A}{\mathit{A}}^{\mathrm{\top }}-{\mathbf{\Sigma }}_p$ $\mathit{B}{\mathit{B}}^{\mathrm{\top }}-{\mathbf{\Sigma }}_q$ 是对称的，所以对应的参数矩阵${\mathit{W}}_p,{\mathit{W}}_q$ 是对称的，此时

$${\mathbf{\Sigma }}_q=\mathit{B}{\mathit{B}}^{\mathrm{\top }}={\left({\mathit{A}}^{\mathrm{\top }}{\mathit{W}}_p\right)}^{\mathrm{\top }}\left({\mathit{A}}^{\mathrm{\top }}{\mathit{W}}_p\right)={\mathit{W}}_p\mathit{A}{\mathit{A}}^{\mathrm{\top }}{\mathit{W}}_p={\mathit{W}}_p{\mathbf{\Sigma }}_p{\mathit{W}}_p$$

令${\mathit{W}}_p={\mathbf{\Sigma }}_p^{-1/2}\mathit{R}{\mathbf{\Sigma }}_p^{-1/2}$，代入上式得${\mathbf{\Sigma }}_q={\mathbf{\Sigma }}_p^{-1/2}{\mathit{R}}^2{\mathbf{\Sigma }}_p^{-1/2}$，即

$$\mathit{R}=({\mathbf{\Sigma }}_p^{1/2}{\mathbf{\Sigma }}_q{\mathbf{\Sigma }}_p^{1/2}{)}^{1/2}$$

而

$$\begin{array}{rl}\text{Tr}(\mathit{A}{\mathit{B}}^{\mathrm{\top }})=& \text{Tr}(\mathit{A}{\mathit{A}}^{\mathrm{\top }}{\mathit{W}}_p)=\text{Tr}({\mathbf{\Sigma }}_p{\mathit{W}}_p)\\ =& \text{Tr}({\mathbf{\Sigma }}_p^{1/2}\mathit{R}{\mathbf{\Sigma }}_p^{-1/2})=\text{Tr}(\mathit{R}{\mathbf{\Sigma }}_p^{-1/2}{\mathbf{\Sigma }}_p^{1/2})\\ =& \text{Tr}(\mathit{R})\end{array}$$

文章小结

本文详细计算了两个多元正态分布的KL散度、巴氏距离和W距离，给出了它们的显式解析解，这些结果在某些场景下可以作为隐变量的正则项使用，来规范隐变量的分布。此外，本文还可以作为比较有挑战性的线性代数练习题，供大家参考练习。