HSIC and MMD

---

1. HSIC

博客链接：HSIC简介：一个有意思的判断相关性的思路

前几天，在机器之心看到这样的一个推送《彻底解决梯度爆炸问题，新方法不用反向传播也能训练ResNet》，当然，媒体的标题党作风我们暂且无视，主要看内容即可。机器之心的这篇文章，介绍的是论文《The HSIC Bottleneck: Deep Learning without Back-Propagation》的成果，里边提出了一种通过HSIC Bottleneck来训练神经网络的算法。

坦白说，这篇论文笔者还没有看明白，因为对笔者来说里边的新概念有点多了。不过论文中的“HSIC”这个概念引起了笔者的兴趣。经过学习，终于基本地理解了这个HSIC的含义和来龙去脉，于是就有了本文，试图给出HSIC的一个尽可能通俗（但可能不严谨）的理解。

背景

HSIC全称“ Hilbert-Schmidt independence criterion”，中文可以叫做“希尔伯特-施密特独立性指标”吧，跟互信息一样，它也可以用来衡量两个变量之间的独立性。

度量相关

我们知道，互信息的基本形式是：

$$I(X,Y)=\iint p(x,y)\log \frac{p(x,y)}{p(x)p(y)}dxdy$$如果$I(X,Y)=0$那么就说明$p(x,y)\equiv p(x)p(y)$，也就是两个变量是相互独立的，否则是相关的。但$\log \frac{p(x,y)}{p(x)p(y)}$这一项意味着我们要用某种方式对概率密度进行估计。

HSIC的作用跟互信息类似，但是跟互信息不一样的是，它不需要估计两个变量的概率密度，而是直接转化为采样的形式。

长期关注本博客的读者都知道，“互信息”是本博客经常出现的概念，我们可以用互信息做新词发现（比如《基于切分的新词发现》），也可以用互信息做无监督学习（比如《深度学习的互信息：无监督提取特征》），互信息的重要性可见一斑。如果说有一个指标可以取代互信息、比互信息还方便，那肯定是笔者必须去学习的对象了。

问题定义

一般来说，我们将问题定义为：

有数据$(x_1,y_1),(x_2,y_2),\dots ,(x_n,y_n)\sim p(x,y)$，判断$p(x,y)$是否恒等于$p(x),p(y)$，即 $x,y$ 是否独立。

严格来讲，如果是对于连续变量，这里的“恒等于”指的是“几乎处处等于”，但我们这里不严格声明这一点。
为了描述的规范，这里设 $x\in X,y\in Y$，而 $f(x),g(y)\in \mathbb{R}$。注意 $x,y$ 可能是两个含义完全不一样的变量，比如 $x$ 可能是“星期一”，$y$可能是“上班”，$p(x,y)$就是“今天是星期一，且今天要上班”的概率。鉴于此，$X,Y$可能是两个完全不一样的域。
基本的思路是去计算互信息，但很多问题中我们都无法很好地估计概率或概率密度。一种可能的方案是转化为对偶问题，用类似对抗的思路去学习互信息（infomax的思路），但这种方法可能会不稳定，而且受到采样方案的影响。最好的方案就是能有一个类似“相关系数”的指标，让我们可以显式地计算和优化这个指标。
HSIC就是冲着这个目标而来的～

HSIC

这里我们尽可能清晰地引入HSIC的概念。然而，“尽可能清晰”不等价于篇幅尽可能短，事实上，下面的篇幅依然会比较长，而且有不少数学公式，但是相对于标准教程里边一上来就引入希尔伯特空间、再生核、各种算子等做法，这里的介绍应该算是对很多不了解相关概念的读者来说都是友好的了。

基本思想

HSIC留意到：

$p(x,y)\equiv p(x)p(y)$ 当且仅当对于任意的 $f,g$，式
$$\begin{array}{rl}C[f,g]=& \iint p(x,y)f(x)g(y)dxdy-\iint p(x)p(y)f(x)g(y)dxdy\\ =& {\mathbb{E}}_{(x,y)\sim p(x,y)}[f(x)g(y)]-{\mathbb{E}}_{x\sim p(x)}[f(x)]{\mathbb{E}}_{y\sim p(y)}[g(y)]\end{array}$$都等于0。

这个结论显然不难理解。有意思的是，等号右边是采样的形式，也就是说我们将这个指标转化为了采样的形式，避免了直接估算概率密度。
所以，我们就有一个判断独立性的方法：选取“足够多”的 $f,g$，然后计算：
$$L_H=\sum _{f,g}(C[f,g]{)}^2$$看 $L_H$ 与 0 的接近程度；反过来，如果在优化问题中，我们希望特征 $x,y$ 尽可能相互独立，那么我们就可以将 $L_H$ 加入到损失函数中。

抽丝剥茧

其实 $L_H$ 的形式已经很好地体现了HSIC的判别思想。下面我们就沿着这个思路，继续抽丝剥茧，逐步地走向HSIC最终的形式。
首先我们把$(C[f,g]{)}^2$算一算：

$$\begin{array}{rl}HSIC(X,Y)=(C[f,g]{)}^2=& ({\mathbb{E}}_{(x,y)\sim p(x,y)}[f(x)g(y)]{)}^2+({\mathbb{E}}_{x\sim p(x)}[f(x)]{)}^2({\mathbb{E}}_{y\sim p(y)}[g(y)]{)}^2\\ & -2({\mathbb{E}}_{(x,y)\sim p(x,y)}[f(x)g(y)])({\mathbb{E}}_{x\sim p(x)}[f(x)])({\mathbb{E}}_{y\sim p(y)}[g(y)])\end{array}$$

然后我们用一个技巧：我们知道 ${\mathbb{E}}_{x\sim p(x)}[f(x)]={\mathbb{E}}_{x^{\prime }\sim p(x^{\prime })}[f(x^{\prime })]$，说明了这个期望值的结果跟随机变量的记号没啥关系。所以我们有
$$\begin{array}{rl}({\mathbb{E}}_{x\sim p(x)}[f(x)]{)}^2=& ({\mathbb{E}}_{x_1\sim p(x)}[f(x_1)])({\mathbb{E}}_{x_2\sim p(x)}[f(x_2)])\\ =& {\mathbb{E}}_{x_1\sim p(x),x_2\sim p(x)}[f(x_1)f(x_2)]\end{array}$$

把其余的项都这样变换，最终我们就可以得到
$$\begin{array}{rl}(C[f,g]{)}^2=& {\mathbb{E}}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}[f(x_1)f(x_2)g(y_1)g(y_2)]\\ & +{\mathbb{E}}_{x_1\sim p(x),x_2\sim p(x),y_1\sim p(y),y_2\sim p(y)}[f(x_1)f(x_2)g(y_1)g(y_2)]\\ & -2{\mathbb{E}}_{(x_1,y_1)\sim p(x,y),x_2\sim p(x),y_2\sim p(y)}[f(x_1)f(x_2)g(y_1)g(y_2)]\end{array}$$这样一来，每一项都是$f(x_1)f(x_2)g(x_1)g(x_2)$的期望，只不过变量的采样分布不一样。

特征函数

现在的问题是：要选择哪些 $f,g$ 呢？怎样才算“足够多”呢？
类比向量空间的知识，所有可能的 $f(x)$ 能组成一个向量空间 $\mathcal{F}$，所有的 $g(y)$ 也一样组成一个向量空间 $\mathcal{G}$。如果能把这两个空间的所有“基底”都遍历一遍，那肯定就够了。那问题就是：如何找到所有的基底呢？
这时候“核函数”就登场了。所谓核函数，那就是——呃，其实说起来很复杂，我也不大懂。简单来说，核函数是类似于线性代数中“正定矩阵”的存在，就是一个定义在 $X\times X$ 的二元对称函数$K(x_1,x_2)=K(x_2,x_1)$，然后我们把一元函数 $f(x)$ 类比为一个向量，那么
$$\int K(x_1,x_2)f(x_2)d{x}_2$$就相当于一个矩阵乘以向量的矩阵运算。跟矩阵的特征值和特征向量一样，核函数也能定义特征值和特征函数，满足下述恒等式的一元函数 $\psi$ 就称为这个核函数的特征函数：
$$\int K(x_1,x_2)\psi (x_2)d{x}_2=\alpha \psi (x_1)$$上面的内容都是铺垫的，其严格定义则是属于“再生核希尔伯特空间“范畴。后面我们用到的，实际上是两点性质：

1、核函数的所有特征函数 ${\psi }_1,{\psi }_2,\dots$ 构成该空间的一组正交基；
2、核函数的所有特征值 ${\alpha }_1,{\alpha }_2,\dots$ 都是正的，且满足
$$K(x_1,x_2)=\sum _i{\alpha }_i{\psi }_i(x_1){\psi }_i(x_2)$$

HSIC登场

经过上述铺垫，HSIC基本上就可以登场了～
首先，假如我们已经有定义在 $X\times X$ 的核函数 $K_X(x_1,x_2)$，那么我们就可以算出$K_X(x_1,x_2)$ 对应的特征值 ${\alpha }_1,{\alpha }_2,\dots$ 和特征函数 ${\psi }_1,{\psi }_2,\dots$；同样地，有了定义在 $Y\times Y$ 的核函数 $K_Y(y_1,y_2)$ 后，也可以算出 $K_Y(y_1,y_2)$ 对应的特征值 ${\beta }_1,{\beta }_2,\dots$ 和特征函数 ${\varphi }_1,{\varphi }_2,\dots$。
然后，因为特征函数构成了基底，所以我们可以把 $f,g$ 换成对应特征函数 ${\psi }_i,{\varphi }_j$
$$L_H=\sum _{i,j}(C[{\psi }_i,{\varphi }_j]{)}^2$$
因为所有的特征值都是正的，所以我们还可以用特征值为权重进行加权求和，而不改变$L_H$的作用：
$$L_H=\sum _{i,j}{\alpha }_i{\beta }_j\cdot (C[{\psi }_i,{\varphi }_j]{)}^2$$
现在我们代入到上面去，就得到
$$\begin{array}{rl}{L}_H=& {\mathbb{E}}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}\left[\sum _{i,j}{\alpha }_i{\beta }_j{\psi }_i(x_1){\psi }_i(x_2){\varphi }_j(y_1){\varphi }_j(y_2)\right]\\ & +{\mathbb{E}}_{x_1\sim p(x),x_2\sim p(x),y_1\sim p(y),y_2\sim p(y)}\left[\sum _{i,j}{\alpha }_i{\beta }_j{\psi }_i(x_1){\psi }_i(x_2){\varphi }_j(y_1){\varphi }_j(y_2)\right]\\ & -2{\mathbb{E}}_{(x_1,y_1)\sim p(x,y),x_2\sim p(x),y_2\sim p(y)}\left[\sum _{i,j}{\alpha }_i{\beta }_j{\psi }_i(x_1){\psi }_i(x_2){\varphi }_j(y_1){\varphi }_j(y_2)\right]\end{array}$$
最后，方括号里边的实际上就是$K_X(x_1,x_2)K_Y(y_1,y_2)$，于是，HSIC就登场了：
$$\begin{array}{rl}HSIC(X,Y)=& {\mathbb{E}}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right]\\ & +{\mathbb{E}}_{x_1\sim p(x),x_2\sim p(x),y_1\sim p(y),y_2\sim p(y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right]\\ & -2{\mathbb{E}}_{(x_1,y_1)\sim p(x,y),x_2\sim p(x),y_2\sim p(y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right]\end{array}$$
这就是我们最重要寻找的度量相关性的指标，它纯粹是采样的形式，而且 $K_X,K_Y$ 都是事先给定的、通常是可微的，因此这就是一个可以明确采样计算、可以直接优化的指标！

在实际计算中，我们可选的核函数有很多，比较常用的是
$$K(x_1,x_2)=\exp \left(-\frac{\Vert x_1-x_2{\Vert }_2^2}{{\sigma }^2}\right)$$其中 $\sigma >0$ 是一个常数，本文开头提到的论文《The HSIC Bottleneck: Deep Learning without Back-Propagation》也是用这个核函数。不同的核函数效果有点不一样，但是都能保证$HSIC(X,Y)=0\iff p(x,y)\equiv p(x)p(y)$。

矩阵形式

最后，我们来推导一下在有限样本下的矩阵形式。

按照采样求期望的思想，${\mathbb{E}}_{(x_1,y_1)\sim p(x,y)}$ 实际上就是对所有的样本对$(x_i,y_i)$的结果求平均，而 ${\mathbb{E}}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}$ 其实就是将这个平均操作做两次，所以：
$${\mathbb{E}}_{(x_1,y_1)\sim p(x,y),(x_2,y_2)\sim p(x,y)}\left[K_X(x_1,x_2)K_Y(y_1,y_2)\right]=\frac{1}{n^2}\sum _{i=1}^{\beta }\sum _{j=1}^{\beta }\left[K_X(x_i,x_j)K_Y(y_i,y_j)\right]$$其中 $K_X(x_i,x_j),K_Y(y_i,y_j)$ 实际上都是 $\beta \times \beta$ 的对称矩阵中 $i$ 行 $j$ 列的元素，分别记为 $K_X{,}_y$，那么上述运算可以写成矩阵乘法 $\frac{1}{{\beta }^2}\text{Tr}(K_X{K}_Y)$，其中 $\text{Tr}$ 表示矩阵的迹。基于同样的思想，第二项实际上就是 “$K_X$ 所有元素的平均乘以 $K_Y$ 所有元素的平均”，如果非要写成矩阵形式的话，那就是 $\frac{1}{{\beta }^4}\text{Tr}(K_{X}1K_{Y}1)$，其中加粗的 $1$ 表示大小为 $\beta \times \beta$ 的全1矩阵。相应地，最后一项是“$K_X{K}_Y$ 所有元素平均值的 $1/\beta$ 的两倍”，即 $\frac{2}{{\beta }^3}\text{Tr}(K_X{K}_{Y}1)$。

所以，如果用矩阵形式表示HSIC，那就是
$$\begin{array}{rl}HSIC(X,Y)& =\frac{1}{{\beta }^2}\text{Tr}(K_X{K}_Y)+\frac{1}{{\beta }^4}\text{Tr}(K_{X}1K_{Y}1)-\frac{2}{{\beta }^3}\text{Tr}(K_X{K}_{Y}1)\\ & =\frac{1}{{\beta }^2}\text{Tr}(K_X{K}_Y+K_X\frac{1}{\mathbf{\beta }}{K}_Y\frac{1}{\mathbf{\beta }}-2K_X{K}_Y\frac{1}{\mathbf{\beta }})\\ & =\frac{1}{{\beta }^2}\text{Tr}(K_{X}J{K}_{Y}J)\end{array}$$其中 $J=\mathbf{I}-1/\mathbf{\beta }$
则无偏估计为：
$$\begin{array}{rl}HSIC(X,Y)& =\frac{1}{(\beta -1{)}^2}\text{Tr}(K_{X}J{K}_{Y}J)\end{array}$$

---

2. MMD (Maximum Mean Discrepancy)

MMD 的含义为 $p$、$q$ 两个概率分布的差距。令 $\mathcal{F}$ 为由 $f:X\to \mathbb{R}$ 组成的一类函数。$p,q$ 是两个概率分布，$X,Y$ 是 从 $p,q$ 采样得到的样本的集合。定义：
$$MMD(\mathcal{F},p,q)=\underset{f\in \mathcal{F}}{\sup }|E_{x\sim p}[f(x)]-E_{y\sim q}[f(y)]|$$

设 $x,y$ 投射到 $\mathcal{F}$ 上的值的均值元素分别为：
$$\begin{gathered} {⟨{\mu }_x,f⟩}_{\mathcal{F}}:={\mathbf{E}}_x\left[⟨\varphi (x),f{⟩}_{\mathcal{F}}\right]={\mathbf{E}}_x[f(x)] \\ {⟨{\mu }_y,f⟩}_{\mathcal{F}}:={\mathbf{E}}_y\left[⟨\varphi (y),f{⟩}_{\mathcal{F}}\right]={\mathbf{E}}_y[f(y)] \end{gathered}$$此外：
$${\Vert {\mu }_x\Vert }_{\mathcal{F}}^2={\mathbf{E}}_{x,x^{\prime }}\left[{⟨\varphi (x),\varphi \left(x^{\prime }\right)⟩}_{\mathcal{F}}\right]={\mathbf{E}}_{x,x^{\prime }}\left[k\left(x,x^{\prime }\right)\right]$$其中，$x$ 和 $x^{\prime }$ 是取自分布 $p_x$ 的独立副本。

则：
M M D 2 ( F , p , q ) = { sup ⁡ f ∈ F ∣ E x ∼ p [ f ( x ) ] − E y ∼ q [ f ( y ) ] ∣ } 2 = { sup ⁡ f ∈ F ∣ < μ x − μ y , f > F ∣ } 2 = ∣ ∣ μ x − μ y ∣ ∣ 2 = < μ p , μ p > F + < μ q , μ q > F − 2 < μ p , μ q > F = E x , x ′ < ϕ ( x ) , ϕ ( x ′ ) > F + E y , y ′ < ϕ ( y ) , ϕ ( y ′ ) > F      − 2 E x , y < ϕ ( x ) , ϕ ( y ) > F = E x , x ′ k ( x , x ′ ) + E y , y ′ k ( y , y ′ ) − 2 E x , y k ( x , y ) \begin{aligned} MMD^2(\mathscr{F}, p, q) &=\{\sup_{f\in\mathscr{F}}| E_{x\sim p}[f(x)] - E_{y\sim q}[f(y)]|\}^2 \\ &=\{\sup_{f\in\mathscr{F}}| <\mu_x-\mu_y, f>_\mathscr{F}|\}^2\\ &=||\mu_x-\mu_y||^2 \\ &=<\mu_{p}, \mu_{p}>_{\mathscr{F}}+<\mu_{q}, \mu_{q}>_{\mathscr{F}}-2<\mu_{p}, \mu_{q}>_{\mathscr{F}} \\ &=E_{x, x^{\prime}}<\phi(x), \phi\left(x^{\prime}\right)>_{\mathscr{F}}+E_{y, y^{\prime}}<\phi(y), \phi\left(y^{\prime}\right)>_{\mathscr{F}} \\ &\ \ \ \ -2 E_{x, y}<\phi(x), \phi(y)>_{\mathscr{F}} \\ &=E_{x, x^{\prime}} k\left(x, x^{\prime}\right)+E_{y, y^{\prime}} k\left(y, y^{\prime}\right)-2 E_{x, y} k(x, y) \end{aligned} MMD2(F,p,q)​={f∈Fsup​∣Ex∼p​[f(x)]−Ey∼q​[f(y)]∣}2={f∈Fsup​∣<μx​−μy​,f>F​∣}2=∣∣μx​−μy​∣∣2=<μp​,μp​>F​+<μq​,μq​>F​−2<μp​,μq​>F​=Ex,x′​<ϕ(x),ϕ(x′)>F​+Ey,y′​<ϕ(y),ϕ(y′)>F​    −2Ex,y​<ϕ(x),ϕ(y)>F​=Ex,x′​k(x,x′)+Ey,y′​k(y,y′)−2Ex,y​k(x,y)​

MMD 的无偏估计计算公式为:

$$\begin{array}{rl}MM{D}^2[\mathcal{F},p,q]& =\frac{1}{m(m-1)}\sum _{i=1}^m\sum _{j\ne i}^{m}k\left(x_i,x_j\right)+\frac{1}{n(n-1)}\sum _{i=1}^n\sum _{j\ne i}^{n}k\left(y_i,y_j\right)\\ & -\frac{2}{mn}\sum _{i=1}^m\sum _{j=1}^{n}k\left(x_i,y_j\right)\end{array}$$