两个多元正态分布的KL散度、巴氏距离和W距离
By 苏剑林 | 2021-07-08 | 95973位读者 |正态分布是最常见的连续型概率分布之一。它是给定均值和协方差后的最大熵分布(参考《“熵”不起:从熵、最大熵原理到最大熵模型(二)》),也可以看作任意连续型分布的二阶近似,它的地位就相当于一般函数的线性近似。从这个角度来看,正态分布算得上是最简单的连续型分布了。也正因为简单,所以对于很多估计量来说,它都能写出解析解来。
本文主要来计算两个多元正态分布的几种度量,包括KL散度、巴氏距离和W距离,它们都有显式解析解。
正态分布 #
这里简单回顾一下正态分布的一些基础知识。注意,仅仅是回顾,这还不足以作为正态分布的入门教程。
概率密度 #
正态分布,也即高斯分布,是定义在$\mathbb{R}^n$上的连续型概率分布,其概率密度函数为
\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}
这里的$\boldsymbol{x},\boldsymbol{\mu}\in\mathbb{R}^n$,$\boldsymbol{\mu}$即均值向量(本文的向量默认情况下都为列向量),而$\boldsymbol{\Sigma}\in\mathbb{R}^{n\times n}$即为协方差矩阵,它要求是正定对称的。可以看到,正态分布由$\boldsymbol{\mu}$和$\boldsymbol{\Sigma}$唯一确定,因此不难想象它的统计量都是$\boldsymbol{\mu}$和$\boldsymbol{\Sigma}$的函数。当$\boldsymbol{\mu}=\boldsymbol{0}, \boldsymbol{\Sigma}=\boldsymbol{I}$时,对应的分布称为“标准正态分布”。
基本性质 #
通常来说,基本的统计量就是均值和方差了,它们对应着正态分布的两个参数:
\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}
由此也可以推出二阶矩的结果:
\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}
还有一个常用的统计量是它的熵:
\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}
其计算过程可以参考后面KL散度的推导。
高斯积分 #
概率密度函数意味着$\int p(\boldsymbol{x}) d\boldsymbol{x} = 1$,这就可以推出:
\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}
设$\boldsymbol{\omega} = \boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}$,那么得到高斯积分
\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}
利用它我们可以算出正态分布的特征函数
\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}
特征函数可以用来算正态分布的各阶矩。
线性代数 #
这里补充一些线性代数基础,它们在后面的推导中会频繁用到。同样地,这仅仅是“回顾”,并不能作为线性代数教程。
内积范数 #
首先,我们来定义内积和范数。对于向量$\boldsymbol{x}=(x_1,\cdots,x_n)$和$\boldsymbol{y}=(y_1,\cdots,y_n)$,内积按照
\begin{equation}\langle\boldsymbol{x},\boldsymbol{y}\rangle = \sum_{i=1}^n x_i y_i\end{equation}
来定义,而模长定义为$\Vert \boldsymbol{x}\Vert = \sqrt{\langle\boldsymbol{x},\boldsymbol{x}\rangle}$。对于$m\times n$的矩阵$\boldsymbol{A}=(a_{i,j}),\boldsymbol{B}=(b_{i,j})$,我们按照类似的方式定义:
\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}
这称为Frobenius内积,对应的$\Vert \boldsymbol{A}\Vert_F = \sqrt{\langle\boldsymbol{A},\boldsymbol{A}\rangle_F}$称为Frobenius范数。不难看到,Frobenius内积和范数,事实上就是把矩阵展平为向量后,当作常规的向量来运算。
关于Frobenius内积,最关键的性质之一是成立恒等式
\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}
也就是说,矩阵的Frobenius内积可以转化为矩阵乘法的迹,并且交换相乘顺序不改变结果(不改变迹的结果,但是矩阵乘法的整体结果会改变)。
正定对称 #
接着,来看正定对称矩阵的一些性质。$\boldsymbol{\Sigma}$是一个正定对称矩阵,对称说的是$\boldsymbol{\Sigma}^{\top}=\boldsymbol{\Sigma}$,正定说的是对于任意非零向量$\boldsymbol{\xi}\in\mathbb{R}^n$,都有$\boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} > 0$。可以证明,如果$\boldsymbol{\Sigma}_1,\boldsymbol{\Sigma}_2$都是正定对称矩阵,那么$\boldsymbol{\Sigma}_1^{-1},\boldsymbol{\Sigma}_2^{-1},\boldsymbol{\Sigma}_1+\boldsymbol{\Sigma}_2$也都是正定对称矩阵。如果$\boldsymbol{C} = \boldsymbol{B}^{\top}\boldsymbol{A}\boldsymbol{B}$,$\boldsymbol{B}$是可逆阵,那么$\boldsymbol{C}$是正定对称的当且仅当$\boldsymbol{A}$是正定对称的。
此外还有半正定的概念,指对于任意非零向量$\boldsymbol{\xi}\in\mathbb{R}^n$,都有$\boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} \geq 0$,也就是说可能存在非零向量$\boldsymbol{\xi}$使得$\boldsymbol{\xi}^{\top}\boldsymbol{\Sigma}\boldsymbol{\xi} = 0$。不过考虑到正定矩阵在半正定矩阵中稠密,所以我们不严格区分正定和半正定了,统一按照正定矩阵来处理。
正定对称矩阵有一个重要的性质,那就是它的SVD分解跟特征值分解一致,即具有下述形式的分解
\begin{equation}\boldsymbol{\Sigma} = \boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{U}^{\top}\end{equation}
其中$\boldsymbol{U}$是正交矩阵,而$\boldsymbol{\Lambda}$是对角阵,并且对角线上的元素都是正的。该结果的一个直接推论是:正定对称矩阵都可以“开平方”,其平方根为$\boldsymbol{\Sigma}^{1/2} = \boldsymbol{U}\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{\top}$,其中$\boldsymbol{\Lambda}^{1/2}$是指将对角线上的元素都开平方,可以检验平方根矩阵也是正定对称的。反过来,可以开平方的对称矩阵,一定也是正定对称矩阵。
矩阵求导 #
最后,在求Wasserstein距离的时候,还需要用到一些矩阵求导公式,如果不了解的读者,可以直接参考维基百科的“Matrix Calculus”。当然,其实也不难,主要用到了
\begin{equation}\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}\end{equation}
剩下的可以结合迹的运算公式来派生出来,比如
\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散度 #
作为第一个尝试,我们来算两个高斯分布的KL散度(Kullback-Leibler divergence)。KL散度算是最常用的分布度量之一了,因为它积分之前需要取对数,这对于指数簇分布来说通常能得到相对简单的结果。此外它还跟“熵”有着比较紧密的联系。
计算结果 #
两个概率分布的KL散度定义为
\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}
对于两个正态分布来说,计算结果是
\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}
特别地,当$q$是标准正态分布时,结果简化为
\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}
推导过程 #
从KL散度的定义知道,我们主要把$\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}
现在,关于迹的恒等式就可以派上用场了:
\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}
注意$\boldsymbol{\mu}_q=\boldsymbol{\mu}_p,\boldsymbol{\Sigma}_q=\boldsymbol{\Sigma}_p$时,上式就等于$n$,此时就对应正态分布的熵。所以最终得到
\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),它定义为
\begin{equation}BD(p(\boldsymbol{x}), q(\boldsymbol{x})) = -\log \int \sqrt{p(\boldsymbol{x}) q(\boldsymbol{x})} d\boldsymbol{x}\end{equation}
与之相关的还有一个叫做“Hellinger距离”的概念,它的平方定义为$\frac{1}{2}\int \left(\sqrt{p(\boldsymbol{x})} - \sqrt{q(\boldsymbol{x})}\right)^2 d\boldsymbol{x}$,展开后就能发现跟巴氏距离本质是等价的。
计算结果 #
对于两个正态分布来说,它们的巴氏距离为
\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}
这里$\boldsymbol{\Sigma}=\frac{1}{2}(\boldsymbol{\Sigma}_p + \boldsymbol{\Sigma}_q)$。可以看到结果是对称的,这是因为巴氏距离的定义本身就是对称的。
当两者之一为标准正态分布时,结果并没有明显简化,所以这里就不单独写出来了。
推导过程 #
按照定义,两个正态分布的巴氏距离,是下述积分的负对数:
\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}
记$\boldsymbol{y}=\boldsymbol{x}-\boldsymbol{\mu}_p, \boldsymbol{\Delta}=\boldsymbol{\mu}_p - \boldsymbol{\mu}_q$,积分部分可以换元为
\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}
这里$\boldsymbol{\Sigma}=\frac{1}{2}(\boldsymbol{\Sigma}_p + \boldsymbol{\Sigma}_q)$。按照前面介绍的高斯积分公式$\eqref{eq:g-int}$,积分结果是
\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}
所以最终
\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距离 #
如果读者还想看了解更多关于概率散度的内容,可以参考书籍《Statistical Inference Based on Divergence Measures》。现在我们转向另一类概率度量——基于最优传输的W距离(Wasserstein距离)。
沿用《从Wasserstein距离、对偶理论到WGAN》中的记号,W距离的定义如下:
\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}
目前仅当$\rho=2$时求出了解析解,因此下面都是在求$\mathcal{W}_2[p,q]$,并且简单起见,记
\begin{equation}\mathcal{W}_2^2[p,q] = \left(\mathcal{W}_2[p,q]\right)^2\end{equation}
计算结果 #
有意思的是,关于两个正态分布的W距离的结果,流传着两个不同的版本,这两个版本都有一定的认知度,但却没有看到有明确说两者等价的资料。两个版本出自不同的论文,还被冠以了不同的名字。
版本1 #
首先是流传相对较广的版本,很多时候搜索“正态分布的Wasserstein距离”都是给出这个结果:
\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}关于这个结果,有的读者可能困惑于“怎么关于$p,q$不是对称的”,事实上,它关于$p,q$是对称的,因为
\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}
然后我们可以直接验证$(\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}$,所以有$\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})$。
版本2 #
第二个看起来稍微简单一些的版本,结果是:
\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}
这个版本通常被称为“Fréchet距离”,一般情况下通过关键词“正态分布的Fréchet距离”才能搜到这个结果。GAN中经常使用的评价指标FID(Frechet Inception Distance),就是基于这个公式进行计算的。可以模仿前面证明它关于$p,q$的对称性,当然也可以从下面的等价性讨论中直接得出。
等价性 #
按道理,第2个版本简洁一些,应该以第二个版本为标准传播才对。所以为什么现在依然流传着两个不同的版本,笔者也相当困惑。从理论上来看,证明两个版本的等价性也不难,根据迹的恒等式我们有:
\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}
然后直接验证$(\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$即可。
特殊时 #
特别地,如果$\boldsymbol{\Sigma}_p,\boldsymbol{\Sigma}_q$的乘法可以交换,那么将会简化为非常直观的形式:
\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}
为什么说它非常直观呢?因为正态分布的参数为$\boldsymbol{\mu},\boldsymbol{\Sigma}$,所以比较正态分布的差异其实就是比较$\boldsymbol{\mu},\boldsymbol{\Sigma}$的差异,按照机器学习的习惯,一个很容易相当想到的指标是平方误差
\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}
但从物理角度来看,这个指标是不妥的,因为如果将$\boldsymbol{\mu}$看成是长度量纲,那么$\boldsymbol{\Sigma}$就具有长度平方的量纲,所以$\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2$和$\Vert \boldsymbol{\Sigma}_p - \boldsymbol{\Sigma}_q\Vert_F^2$是具有不同量纲的两个量,不能相加。而为了使得量纲一致,直观的想法就是把$\boldsymbol{\Sigma}$“开平方”后再算平方误差,这就得到了式$\eqref{eq:w-jiaohuan}$。
特别地,当$q$为标准正态分布时,结果简化为
\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}
推导过程1 #
现在介绍第一个证明,主要参考了论文《A class of Wasserstein metrics for probability distributions》。另外《The distance between two random vectors with given dispersion matrices》也提供了一个类似的证明,也可以参考。
下面的推导过程则是经过笔者简化的,相对原论文的证明来说简单一些,但依然不可避免地会涉及到较多的线性代数知识,我们将分几个部分介绍。
去均值 #
不失一般性,我们可以只考虑均值为0的分布$p,q$。因为如果$p,q$的均值不为0,那么设对应的均值为0的分布为$\tilde{p},\tilde{q}$,此时有
\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}
该结果意味着
\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}
所以,只需要算出均值都为零时的Wasserstein距离,然后加上$\Vert \boldsymbol{\mu}_p - \boldsymbol{\mu}_q\Vert^2$就得到了一般情况的结果。
纯代数 #
现在我们假设$p,q$的均值均为0,然后计算
\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}
其中
\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}
构成联合分布$\gamma$的协方差矩阵。我们知道协方差矩阵是正定对阵矩阵,所以从代数的角度看,问题变成了:
已知$\boldsymbol{\Sigma}_{\gamma}= \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix}$为正定对称矩阵,求$\text{Tr}(\boldsymbol{C})$的最大值。
舒尔补 #
为此,我们需要利用下述关于“舒尔补”的恒等式:
\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}
其中对称矩阵$\boldsymbol{S} = \boldsymbol{\Sigma}_q - \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$称为“舒尔补(Schur Complement)”,该分解具有$\boldsymbol{B}^{\top}\boldsymbol{A}\boldsymbol{B}$的形式,要想它是正定的,那么$\boldsymbol{A}$要是正定的,而$\boldsymbol{\Sigma}_p$已经是正定的,所以$\boldsymbol{S}$需要是正定的。
分参数 #
我们尝试分离参数,即从$\boldsymbol{S} = \boldsymbol{\Sigma}_q - \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$中把$\boldsymbol{C}$解出来。首先移项得到$\boldsymbol{\Sigma}_q - \boldsymbol{S} = \boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$,由于$\boldsymbol{\Sigma}_p$是正定对称的,所以$\boldsymbol{\Sigma}_p^{-1}$也是,从而$\boldsymbol{C}^{\top}\boldsymbol{\Sigma}_p^{-1}\boldsymbol{C}$也是正定对称的,那么它具有正定对称的平方根,即存在正定对称矩阵$\boldsymbol{R}$,使得
\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}
这说明$\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{C}\boldsymbol{R}^{-1}$是正交矩阵,记为$\boldsymbol{O}$,那么$\boldsymbol{C} = \boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R}$。
乘子法 #
此时,变量分别是$\boldsymbol{O}$和$\boldsymbol{R}$,求$\text{Tr}(\boldsymbol{C})=\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R})$的最大值。我们先固定$\boldsymbol{R}$,求取最大值时的$\boldsymbol{O}$,此时相当于在$\boldsymbol{O}^{\top}\boldsymbol{O}=\boldsymbol{I}$的约束下,求$\text{Tr}(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R})$的最大值,我们用“拉格朗日乘子法”:引入新参数矩阵$\boldsymbol{W}$,转化为下述无约束极值问题
\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}
求导:
\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}
首先留意到$\boldsymbol{O}^{\top}\boldsymbol{O} - \boldsymbol{I}$是对称的,因此对应的参数矩阵$\boldsymbol{W}$也是对称的,于是我们有:
\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}
即$\boldsymbol{O}\boldsymbol{W}\boldsymbol{O}^{\top}=(\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{R}^2\boldsymbol{\Sigma}_p^{1/2})^{1/2}$,所以此时
\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}
不等式 #
最后需要把$\boldsymbol{R}$确定下来。回顾$\boldsymbol{R}$的定义,我们有$\boldsymbol{R}^2=\boldsymbol{\Sigma}_q - \boldsymbol{S}$,其中$\boldsymbol{S}$是正定矩阵。直觉上$\boldsymbol{S}=\boldsymbol{0}$时取得最大值,事实上也确实如此,这算是“Weyl不等式”的一个直接推论。
根据Weyl不等式,如果矩阵$\boldsymbol{A},\boldsymbol{B}$都是正定对称矩阵,设$\boldsymbol{A},\boldsymbol{B},\boldsymbol{A}+\boldsymbol{B}$的特征值从小到大排列分别为$0\leq\lambda_1^{(A)} \leq \dots \leq \lambda_n^{(A)}$、$\lambda_1^{(B)} \leq \dots \leq \lambda_n^{(B)}$和$0\leq\lambda_1^{(A+B)} \leq \dots \leq \lambda_n^{(A+B)}$,那么对于任意$1\leq i \leq n$,都有$\lambda_i^{(A)}\leq \lambda_i^{(A+B)}$和$\lambda_i^{(B)}\leq \lambda_i^{(A+B)}$,也就是说:
正定对称矩阵的和的特征值,一一对应地大于它们各自的特征值。
有了这个结论,那就简单了,设$(\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2})^{1/2}$的特征值为$0 \leq \lambda_1 \leq \dots \leq \lambda_n$,那么它的迹就是$\lambda_1 + \dots + \lambda_n$,对应地,$\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2}$的特征值为$0 \leq \lambda_1^2 \leq \dots \leq \lambda_n^2$,注意$\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2}$是正定对称矩阵(对称是显然的,而因为它能开平方,所以正定),$\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{S}\boldsymbol{\Sigma}_p^{1/2}$也是正定对称的(因为$\boldsymbol{S}$是正定对称的),所以它们的特征值,都不超过它们的和——也就是$\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2}$的特征值,所以说,$(\boldsymbol{\Sigma}_p^{1/2}(\boldsymbol{\Sigma}_q - \boldsymbol{S})\boldsymbol{\Sigma}_p^{1/2})^{1/2}$每个特征值的最大值(也就是迹的最大值),在$\boldsymbol{S}=\boldsymbol{0}$处取到。
至于Weyl不等式的证明,主要利用到了Rayleigh quotient和Courant–Fischer定理,有兴趣了解的读者自行查阅这两部分资料后,再查阅Wely不等式的证明就好。事实上,熟悉这两部分内容后,Weyl不等式基本上就“水到渠成”了。
推导过程2 #
这里继续介绍另一个更为简单的证明,原始证明可以在《The Fréchet distance between multivariate normal distributions》找到。相对于第一个证明而言,该证明更简单直接,尤其是不需要太多的纯线性代数知识。下面的推导过程依旧是经过笔者进一步简化的,比原始论文更好理解一些。
在这个推导过程中,“去均值”、“纯代数”两个步骤跟“推导过程1”是一样的,不再重复。所以,此时问题已经被转化为
已知$\boldsymbol{\Sigma}_{\gamma}= \begin{pmatrix} \boldsymbol{\Sigma}_p & \boldsymbol{C}\\ \boldsymbol{C}^{\top} & \boldsymbol{\Sigma}_q\end{pmatrix}$为正定对称矩阵,求$\text{Tr}(\boldsymbol{C})$的最大值。
分块阵 #
由于$\boldsymbol{\Sigma}_{\gamma}$是正定对称矩阵,所以它必然可以表达成$\boldsymbol{D}\boldsymbol{D}^{\top}$的形式,我们将$\boldsymbol{D}$表达为分块矩阵$\begin{pmatrix}\boldsymbol{A} \\ \boldsymbol{B}\end{pmatrix}$,其中$\boldsymbol{A},\boldsymbol{B}\in \mathbb{R}^{n\times 2n}$,此时
\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}
对应地有$\boldsymbol{\Sigma}_p=\boldsymbol{A}\boldsymbol{A}^{\top}, \boldsymbol{\Sigma}_q=\boldsymbol{B}\boldsymbol{B}^{\top},\boldsymbol{C}=\boldsymbol{A}\boldsymbol{B}^{\top}$。
乘子法 #
在上述参数化之下,问题转化为:
已知$\boldsymbol{A}\boldsymbol{A}^{\top}=\boldsymbol{\Sigma}_p, \boldsymbol{B}\boldsymbol{B}^{\top}=\boldsymbol{\Sigma}_q$,求$\text{Tr}(\boldsymbol{A}\boldsymbol{B}^{\top})$的最大值。
这是一个带约束的最大值问题,我们用“拉格朗日乘子法”:引入新参数矩阵$\boldsymbol{W}_p, \boldsymbol{W}_q$,转化为下述无约束极值问题
\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}
求导:
\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}
注意到$\boldsymbol{A}\boldsymbol{A}^{\top} - \boldsymbol{\Sigma}_p$和$\boldsymbol{B}\boldsymbol{B}^{\top} - \boldsymbol{\Sigma}_q$都是对称的,所以对应的参数矩阵$\boldsymbol{W}_p, \boldsymbol{W}_q$也是对称的,此时
\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}
令$\boldsymbol{W}_p=\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{R}\boldsymbol{\Sigma}_p^{-1/2}$,代入上式得$\boldsymbol{\Sigma}_q=\boldsymbol{\Sigma}_p^{-1/2}\boldsymbol{R}^2\boldsymbol{\Sigma}_p^{-1/2}$,即
\begin{equation}\boldsymbol{R} = (\boldsymbol{\Sigma}_p^{1/2}\boldsymbol{\Sigma}_q\boldsymbol{\Sigma}_p^{1/2})^{1/2}
\end{equation}
而
\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}
文章小结 #
本文详细计算了两个多元正态分布的KL散度、巴氏距离和W距离,给出了它们的显式解析解,这些结果在某些场景下可以作为隐变量的正则项使用,来规范隐变量的分布。此外,本文还可以作为比较有挑战性的线性代数练习题,供大家参考练习。
转载到请包括本文地址:https://www.spaces.ac.cn/archives/8512
更详细的转载事宜请参考:《科学空间FAQ》
如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。
如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!
如果您需要引用本文,请参考:
苏剑林. (Jul. 08, 2021). 《两个多元正态分布的KL散度、巴氏距离和W距离 》[Blog post]. Retrieved from https://www.spaces.ac.cn/archives/8512
@online{kexuefm-8512,
title={两个多元正态分布的KL散度、巴氏距离和W距离},
author={苏剑林},
year={2021},
month={Jul},
url={\url{https://www.spaces.ac.cn/archives/8512}},
}
July 15th, 2021
公式(8)下面的模长定义那里有个笔误~
谢谢反馈,已经修正。
July 26th, 2021
解决了我多年对这几个距离的困惑,非常感谢。我想请问,我可以转载到【运筹or帷幄】上吗?我会注明本网站地址。
可以~
July 26th, 2021
苏神,公式(7)怎么由公式(6)推导出来,能详细一下吗?
公式$(6)$意味着
$$\int\exp\left\{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})+\boldsymbol{\omega}^{\top}(\boldsymbol{x}-\boldsymbol{\mu})\right\}d\boldsymbol{x} = \sqrt{(2\pi)^n \det(\boldsymbol{\Sigma})}\exp\left\{\frac{1}{2}\boldsymbol{\omega}^{\top}\boldsymbol{\Sigma}\boldsymbol{\omega}\right\}
$$
两边除以$\sqrt{(2\pi)^n \det(\boldsymbol{\Sigma})}$,得到
$$\mathbb{E}_{\boldsymbol{x}}\left[\exp\left\{\boldsymbol{\omega}^{\top}(\boldsymbol{x}-\boldsymbol{\mu})\right\}\right] = \exp\left\{\frac{1}{2}\boldsymbol{\omega}^{\top}\boldsymbol{\Sigma}\boldsymbol{\omega}\right\}$$
这跟$(7)$等价。
August 11th, 2021
公式23推导得到下面公式24的结论过程中,直接应用了公式6,那么是不是隐含了一个前提,就是sigma_q的逆是等于sigma逆/2(此处simga是指公式23中的(sigma_q+sigma_p)/2);因为公式6成立是假设omiga是sigma的逆乘以mu,从而保证概率积分为1。两个公式对照,那么这样就意味着,sigma_p和sigma_q是相同的???这个假设成立吗? 还是我理解有误,请指正,谢谢!
公式$(6)$是一个恒等式,$\boldsymbol{\omega}$跟$\boldsymbol{\Sigma}$是相互独立的,没有联系。
确实如此,公式6是一个恒等式。非常感谢!指出了一个认识误区
September 15th, 2021
写的很详细,数学功底很强,但是等式(42)的求偏导似乎存在问题,少了一个系数2,在后面的推导中可以约掉,结论是正确的。
确实是,感谢指出,已经修正。
November 21st, 2021
苏神,在计算巴氏距离时套用高斯公式(推导时$\boldsymbol{\omega} = \boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}$),$\boldsymbol{\omega}$是否需要满足一定条件?例如存在$\boldsymbol{\mu}$使得$\boldsymbol{\omega}$能够分解为$\boldsymbol{\Sigma^{-1}\mu}$。
$\boldsymbol{\Sigma}$可逆的情况下,还需要加什么条件?只是一个线性变换而已。
November 22nd, 2021
苏神,“注意到$\boldsymbol{A}\boldsymbol{A}^{\top} - \boldsymbol{\Sigma}_p$和$\boldsymbol{B}\boldsymbol{B}^{\top} - \boldsymbol{\Sigma}_q$都是对称的,所以对应的参数矩阵$\boldsymbol{W}_p, \boldsymbol{W}_q$也是对称的”,后一段的所以是怎么得到的?
比如我们希望引入约束$f_{i,j}=0$,如果$f_{i,j}$是对称的,那么意味着$f_{i,j}=0$跟$f_{j,i}=0$是同一个约束。我们用拉格朗日乘数法引入$\sum_{i,j}\lambda_{i,j} f_{i,j}$,事实上需要引入$n(n+1)/2$个$\lambda_{i,j}$就够了;如果非要引入$n^2$个$\lambda_{i,j}$,那么就要让相同约束对应的$\lambda_{i,j}$相等,也就是$\lambda_{i,j}=\lambda_{j,i}$。
谢谢苏神解惑!
还有个问题,也是因为约束条件引入了$n^2$个约束,所以约束前乘了$1/2$吗,我在求导时发现刚好两个$1/2$合到一起了。
$$\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}$$
1/2是为了形式简单而引入的而已,你改为1/3也不改变结果。相当于求$f(x)$的最大值,跟求$f(x/2)$的最大值,结果是等价的(如果$x\in\mathbb{R}$的话)
November 25th, 2021
苏神,关于标量对矩阵求导,我看了一些资料发现跟你的公式有不同。
\begin{equation}\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}\end{equation}
根据微分与导数的关系求梯度(https://zhuanlan.zhihu.com/p/24709748)
\begin{equation}
d(\text{tr}(\boldsymbol{XA}))
=\text{tr}(d\boldsymbol{XA})
=\text{tr}(\boldsymbol{A}d\boldsymbol{X})
=\text{tr}(\boldsymbol{(A^T)}^Td\boldsymbol{X})
\end{equation}
然后根据导数与微分的关系
\begin{equation}
df=\text{tr}(\frac{\partial{f}}{\partial \boldsymbol{X}}^Td\boldsymbol{X})
\end{equation}
得到结果为$\boldsymbol{A}^T$。
根据维度,我觉得这个结果合理,假如$f=\text{tr}(\boldsymbol{X}_{m\times n}A_{n\times m})$,则$\frac{\partial\,f}{\partial \boldsymbol{X}_{m\times n}}=\boldsymbol{A_{n\times m}}^T$维度相等。
请问,如果改为$\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}^{\top}$,本文的结果有任何实质变化吗?
事实上,$\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}$还是$\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}^{\top}$,纯粹是人为定义的,你所谓的“导数与微分的关系”$df=\text{tr}(\frac{\partial{f}}{\partial \boldsymbol{X}}^{\top} d\boldsymbol{X})$,我要是定义为$df=\text{tr}(\frac{\partial{f}}{\partial \boldsymbol{X}} d\boldsymbol{X})$,也无任何不妥之处。
所以,不用纠结这些形式,我一开始就给出$\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}$,就是要声明我这里采用的是$\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}$这种定义而不是$\frac{\partial\,\text{Tr}\left(\boldsymbol{X}\boldsymbol{A}\right)}{\partial \boldsymbol{X}} = \boldsymbol{A}^{\top}$这种定义,只要我提前声明了,就不至于引起混淆,也就可以合理阅读下去了。
明白了,谢谢苏神~
September 27th, 2022
写的太好了!!! 完美解答了我看这个
https://djalil.chafai.net/blog/2010/04/30/wasserstein-distance-between-two-gaussians/
疑惑的点!!!
August 22nd, 2023
苏神~~~$(40)$分参数那段泥说由于$\sum_p$是正定对称的,所以$\sum^{−1}_p$也是,从而$C^⊤\sum^{−1}_pC$也是正定对称的.这么说的前提是否是:"$C$是可逆的"捏?那咱们如何证明$C$是可逆得丫~
苏神在嘛~~请问是我的理解有问题嘛
抱歉,之前看漏了这个留言。这里似乎没用到$\boldsymbol{C}$的逆矩阵?为什么要求$\boldsymbol{C}$可逆?
你是担心“正定性”?事实上这里半正定问题也不大。
是滴苏神因为在本文正定对称 #那一段(就是包含$(11)$的那一段)那个定理,$D=B^⊤AB$, $B$是可逆的,那么$D$是正定对称的当且仅当$A$是正定对称的,这么做的前提是$B$是可逆的嘛。如果$B$不可逆那么如果$A$是正定对称的,$D$就有可能是正定,但也有可能是半正定. 所以转到$(40)$这里苏神您的意思就是$C^T\sum_p^{-1}C$是半正定也没事儿是吗~~0.0~~
反复看了一下,$(40)$的最终推导结果是希望得出$\boldsymbol{C} = \boldsymbol{\Sigma}_p^{1/2}\boldsymbol{O}\boldsymbol{R}$,其中$\boldsymbol{O}$是正交矩阵,而$\boldsymbol{O}$是(半)正定对称矩阵。如果$\boldsymbol{C}$是半正定的,那么直接让$\boldsymbol{R}$也是半正定就行,似乎不大影响整个推导过程,细节我没有细抠,不过这里的证明本身就是摘录简化自原始证明,应该不会有啥问题。
好滴苏神~~~谢谢泥~~~那我也多看两遍~~~O(∩_∩)O