【强基固本】MMD：最大均值差异

“强基固本，行稳致远”，科学研究离不开理论基础，人工智能学科更是需要数学、物理、神经科学等基础学科提供有力支撑，为了紧扣时代脉搏，我们推出“强基固本”专栏，讲解AI领域的基础知识，为你的科研学习提供助力，夯实理论基础，提升原始创新能力，敬请关注。

来源：知乎—卡兵

地址：https://zhuanlan.zhihu.com/p/163839117

01

怎么描述一个随机变量

去描述一个随机变量，最直接的方法就是给出它的概率分布函数  。一些简单的分布可以这么干，比如正太分布给出均值和方差就可以确定，但是对于一些复杂的、高维的随机变量，我们无法给出它们的分布函数。

这时候我们可以用随机变量的矩来描述一个随机变量，比如一阶中心矩是均值，二阶中心矩是方差等等。如果两个分布的均值和方差都相同的话，它们应该很相似，比如同样均值和方差的高斯分布和拉普拉斯分布。但是很明显，均值和方差并不能完全代表一个分布，这时候我们就需要更高阶的矩来描述一个分布[1]。

举个例子，就好比描述人一样，如果两个人身高、轮廓都一样，我们会说这两个人很像。但是如果要说这两个人是一个人的话，我们如要更多的信息，比如血型、DNA等更加复杂的信息

而MMD的基本思想就是，如果两个随机变量的任意阶都相同的话，那么两个分布就是一致的。而当两个分布不相同的话，那么使得两个分布之间差距最大的那个矩应该被用来作为度量两个分布的标准。

02

如何衡量两个随机变量的差异

MMD常被用来度量两个分布之间的距离，是迁移学习中常用的损失函数。定义如下，  的分布为    ，    的分布为    [2]：

这个公式里面有四个符号，第一个是    求上界（简单理解就是求最大值），第二个是    表示求期望，第三个是  表示映射函数，第四个表示函数    在再生希尔伯特空间中的范数应该小于等于1。

举个简单的例子，    就相当于将    映射到高阶上去，比如    ，那么对应的期望值    就相当于在求    阶矩（    可以是1,2,3,...），然后将他们的上确界（最大值）作为MMD的值。注意这里举的例子只是便于理解。

03

如何表示一个随机变量的任意阶矩

刚才讲到，两个分布应该是由高阶矩来描述的，那么如何获得一个随机变量的高阶矩呢？核函数。

在支持向量机中我们都知道有一个高斯核函数，它对应的映射函数恰好可以映射到无穷维上，映射到无穷维上再求期望，正好可以得到随机变量的高阶矩，这个方法有一个更高大上的名字，叫做kernel embedding of distributions[2]，这个简单理解就是将一个分布映射到再生希尔伯特空间（每个核函数都对应一个RKHS）上的一个点，

这样两个分布之间的距离就可以用两个点的内积进行表示!

这样两个分布之间的距离就可以用两个点的内积进行表示!

这样两个分布之间的距离就可以用两个点的内积进行表示! 如下图[3]

那么一个随机变量的经过映射    之后，期望就可以用表示为，如果对下面的表述方式不熟悉的话可以参考这位大佬的博客[4]：

第一个等号就是利用期望展开；第二个等号利用的是RKHS的再生性，用人话来讲就是RKHS中的    都可以写成一个无穷维的向量  与基底向量  的内积；第三个等号利用的是内积的性质；最后一个等号 表示的就是kernel mean embedding。

这个等价的概念就是将    利用  映射到无穷维上，然后在每一个维度上都求期望[5]。

04

MMD的数学化简

MMD中需要求一个sup，我们要将它用更直观的形式表达出来。

这里的第一个式子就是MMD的直接定义，第二个等号利用的是kernel embedding那一节中最后的结果。注意到，

这里利用的是内积的性质，，而对限制的目的也是为了上界避免随意取到无穷大。因此，可以得到，

由于  无法直接得到，但是可以用均值替代计算，因为均值是期望的无偏估计。假设，并且    有    个样本，    有  个样本。

到这里我们还算是没有办法求，因为  是无穷维的。但是核技巧的关键就在于不需要显式地表示映射函数  来求两个向量的内积。因此我们对MMD进行平方，化简得到内积并用核函数表达[6]。

注意这里的  之类矩阵的就是周志华老师《机器学习》里提到的核矩阵（也叫做Gram矩阵）

05

代码实现

在实际的应用中一般使用的是多核MMD，在这里实现的时候利用的就是多个不同RBF-bandwith的高斯核。

pytorch版本 easezyc/deep-transfer-learning github [7]

import numpy as np
import torch

def guassian_kernel(source, target, kernel_mul=2.0, kernel_num=5, fix_sigma=None):
    """计算Gram核矩阵
    source: sample_size_1 * feature_size 的数据
    target: sample_size_2 * feature_size 的数据
    kernel_mul: 这个概念不太清楚，感觉也是为了计算每个核的bandwith
    kernel_num: 表示的是多核的数量
    fix_sigma: 表示是否使用固定的标准差
        return: (sample_size_1 + sample_size_2) * (sample_size_1 + sample_size_2)的
                        矩阵，表达形式:
                        [   K_ss K_st
                            K_ts K_tt ]
    """
    n_samples = int(source.size()[0])+int(target.size()[0])
    total = torch.cat([source, target], dim=0) # 合并在一起

    total0 = total.unsqueeze(0).expand(int(total.size(0)), \
                                       int(total.size(0)), \
                                       int(total.size(1)))
    total1 = total.unsqueeze(1).expand(int(total.size(0)), \
                                       int(total.size(0)), \
                                       int(total.size(1)))
    L2_distance = ((total0-total1)**2).sum(2) # 计算高斯核中的|x-y|

    # 计算多核中每个核的bandwidth
    if fix_sigma:
        bandwidth = fix_sigma
    else:
        bandwidth = torch.sum(L2_distance.data) / (n_samples**2-n_samples)
    bandwidth /= kernel_mul ** (kernel_num // 2)
    bandwidth_list = [bandwidth * (kernel_mul**i) for i in range(kernel_num)]

    # 高斯核的公式，exp(-|x-y|/bandwith)
    kernel_val = [torch.exp(-L2_distance / bandwidth_temp) for \
                  bandwidth_temp in bandwidth_list]

    return sum(kernel_val) # 将多个核合并在一起

def mmd(source, target, kernel_mul=2.0, kernel_num=5, fix_sigma=None):
    n = int(source.size()[0])
    m = int(target.size()[0])

    kernels = guassian_kernel(source, target,
                              kernel_mul=kernel_mul, kernel_num=kernel_num, fix_sigma=fix_sigma)
    XX = kernels[:n, :n] 
    YY = kernels[n:, n:]
    XY = kernels[:n, n:]
    YX = kernels[n:, :n]

    XX = torch.div(XX, n * n).sum(dim=1).view(1,-1)  # K_ss矩阵，Source<->Source
    XY = torch.div(XY, -n * m).sum(dim=1).view(1,-1) # K_st矩阵，Source<->Target

    YX = torch.div(YX, -m * n).sum(dim=1).view(1,-1) # K_ts矩阵,Target<->Source
    YY = torch.div(YY, m * m).sum(dim=1).view(1,-1)  # K_tt矩阵,Target<->Target

    loss = (XX + XY).sum() + (YX + YY).sum()
    return loss

if __name__ == "__main__":
    # 样本数量可以不同，特征数目必须相同

    # 100和90是样本数量，50是特征数目
    data_1 = torch.tensor(np.random.normal(loc=0,scale=10,size=(100,50)))
    data_2 = torch.tensor(np.random.normal(loc=10,scale=10,size=(90,50)))
    print("MMD Loss:",mmd(data_1,data_2))

    data_1 = torch.tensor(np.random.normal(loc=0,scale=10,size=(100,50)))
    data_2 = torch.tensor(np.random.normal(loc=0,scale=9,size=(80,50)))

    print("MMD Loss:",mmd(data_1,data_2))

# MMD Loss: tensor(1.0866, dtype=torch.float64)
# MMD Loss: tensor(0.0852, dtype=torch.float64)

参考

1. 随机变量的矩和高阶矩有什么实在的含义？ - 姚岑卓的回答 https://www.zhihu.com/question/25344430/answer/64509141

2. abA Hilbert Space Embedding for Distributions https://www.cc.gatech.edu/~lsong/papers/SmoGreSonSch07.pdf

3. Kernel Mean Embedding of Distributions: A Review and Beyond https://arxiv.org/abs/1605.09522

4. A Story of Basis and Kernel http://songcy.net/posts/story-of-basis-and-kernel-part-2/

5. Kernel Distribution Embedding - 李新春的文章 - 知乎 https://zhuanlan.zhihu.com/p/114264831

6. MMD计算的核技巧公式推导 - 王晋东不在家的文章 https://zhuanlan.zhihu.com/p/63026435

7. https://github.com/easezyc/deep-transfer-learning https://github.com/easezyc/deep-transfer-learning

本文目的在于学术交流，并不代表本公众号赞同其观点或对其内容真实性负责，版权归原作者所有，如有侵权请告知删除。

“强基固本”历史文章

更多强基固本专栏文章，

请点击文章底部“阅读原文”查看

分享、点赞、在看，给个三连击呗！