\(\newcommand{\ci}{\perp\!\!\!\perp}\) \(\newcommand{\cA}{\mathcal{A}}\) \(\newcommand{\cB}{\mathcal{B}}\) \(\newcommand{\cC}{\mathcal{C}}\) \(\newcommand{\cD}{\mathcal{D}}\) \(\newcommand{\cE}{\mathcal{E}}\) \(\newcommand{\cF}{\mathcal{F}}\) \(\newcommand{\cG}{\mathcal{G}}\) \(\newcommand{\cH}{\mathcal{H}}\) \(\newcommand{\cI}{\mathcal{I}}\) \(\newcommand{\cJ}{\mathcal{J}}\) \(\newcommand{\cK}{\mathcal{K}}\) \(\newcommand{\cL}{\mathcal{L}}\) \(\newcommand{\cM}{\mathcal{M}}\) \(\newcommand{\cN}{\mathcal{N}}\) \(\newcommand{\cO}{\mathcal{O}}\) \(\newcommand{\cP}{\mathcal{P}}\) \(\newcommand{\cQ}{\mathcal{Q}}\) \(\newcommand{\cR}{\mathcal{R}}\) \(\newcommand{\cS}{\mathcal{S}}\) \(\newcommand{\cT}{\mathcal{T}}\) \(\newcommand{\cU}{\mathcal{U}}\) \(\newcommand{\cV}{\mathcal{V}}\) \(\newcommand{\cW}{\mathcal{W}}\) \(\newcommand{\cX}{\mathcal{X}}\) \(\newcommand{\cY}{\mathcal{Y}}\) \(\newcommand{\cZ}{\mathcal{Z}}\) \(\newcommand{\bA}{\mathbf{A}}\) \(\newcommand{\bB}{\mathbf{B}}\) \(\newcommand{\bC}{\mathbf{C}}\) \(\newcommand{\bD}{\mathbf{D}}\) \(\newcommand{\bE}{\mathbf{E}}\) \(\newcommand{\bF}{\mathbf{F}}\) \(\newcommand{\bG}{\mathbf{G}}\) \(\newcommand{\bH}{\mathbf{H}}\) \(\newcommand{\bI}{\mathbf{I}}\) \(\newcommand{\bJ}{\mathbf{J}}\) \(\newcommand{\bK}{\mathbf{K}}\) \(\newcommand{\bL}{\mathbf{L}}\) \(\newcommand{\bM}{\mathbf{M}}\) \(\newcommand{\bN}{\mathbf{N}}\) \(\newcommand{\bO}{\mathbf{O}}\) \(\newcommand{\bP}{\mathbf{P}}\) \(\newcommand{\bQ}{\mathbf{Q}}\) \(\newcommand{\bR}{\mathbf{R}}\) \(\newcommand{\bS}{\mathbf{S}}\) \(\newcommand{\bT}{\mathbf{T}}\) \(\newcommand{\bU}{\mathbf{U}}\) \(\newcommand{\bV}{\mathbf{V}}\) \(\newcommand{\bW}{\mathbf{W}}\) \(\newcommand{\bX}{\mathbf{X}}\) \(\newcommand{\bY}{\mathbf{Y}}\) \(\newcommand{\bZ}{\mathbf{Z}}\) \(\newcommand{\ba}{\mathbf{a}}\) \(\newcommand{\bb}{\mathbf{b}}\) \(\newcommand{\bc}{\mathbf{c}}\) \(\newcommand{\bd}{\mathbf{d}}\) \(\newcommand{\be}{\mathbf{e}}\) \(\newcommand{\bg}{\mathbf{g}}\) \(\newcommand{\bh}{\mathbf{h}}\) \(\newcommand{\bi}{\mathbf{i}}\) \(\newcommand{\bj}{\mathbf{j}}\) \(\newcommand{\bk}{\mathbf{k}}\) \(\newcommand{\bl}{\mathbf{l}}\) \(\newcommand{\bm}{\mathbf{m}}\) \(\newcommand{\bn}{\mathbf{n}}\) \(\newcommand{\bo}{\mathbf{o}}\) \(\newcommand{\bp}{\mathbf{p}}\) \(\newcommand{\bq}{\mathbf{q}}\) \(\newcommand{\br}{\mathbf{r}}\) \(\newcommand{\bs}{\mathbf{s}}\) \(\newcommand{\bt}{\mathbf{t}}\) \(\newcommand{\bu}{\mathbf{u}}\) \(\newcommand{\bv}{\mathbf{v}}\) \(\newcommand{\bw}{\mathbf{w}}\) \(\newcommand{\bx}{\mathbf{x}}\) \(\newcommand{\by}{\mathbf{y}}\) \(\newcommand{\bz}{\mathbf{z}}\) \(\newcommand{\RR}{\mathbb{R}}\) \(\newcommand{\NN}{\mathbb{N}}\) \(\newcommand{\balpha}{\boldsymbol{\alpha}}\) \(\newcommand{\bbeta}{\boldsymbol{\beta}}\) \(\newcommand{\btheta}{\boldsymbol{\theta}}\) \(\newcommand{\hpi}{\widehat{\pi}}\) \(\newcommand{\bpi}{\boldsymbol{\pi}}\) \(\newcommand{\hbpi}{\widehat{\boldsymbol{\pi}}}\) \(\newcommand{\bxi}{\boldsymbol{\xi}}\) \(\newcommand{\bmu}{\boldsymbol{\mu}}\) \(\newcommand{\bepsilon}{\boldsymbol{\epsilon}}\) \(\newcommand{\bzero}{\mathbf{0}}\) \(\newcommand{\T}{\text{T}}\) \(\newcommand{\Trace}{\text{Trace}}\) \(\newcommand{\Cov}{\text{Cov}}\) \(\newcommand{\Var}{\text{Var}}\) \(\newcommand{\E}{\mathbb{E}}\) \(\newcommand{\Pr}{\text{Pr}}\) \(\newcommand{\pr}{\text{pr}}\) \(\newcommand{\pdf}{\text{pdf}}\) \(\newcommand{\P}{\text{P}}\) \(\newcommand{\p}{\text{p}}\) \(\newcommand{\One}{\mathbf{1}}\) \(\newcommand{\argmin}{\operatorname*{arg\,min}}\) \(\newcommand{\argmax}{\operatorname*{arg\,max}}\) \(\newcommand{\dtheta}{\frac{\partial}{\partial\theta} }\) \(\newcommand{\ptheta}{\nabla_\theta}\) \(\newcommand{\alert}[1]{\color{darkorange}{#1}}\) \(\newcommand{\alertr}[1]{\color{red}{#1}}\) \(\newcommand{\alertb}[1]{\color{blue}{#1}}\)

1 Motivation: Representing a Distribution

Suppose we are given a random variable \(X\) with distribution \(P\) on a domain \(\cX\). The distribution is a complicated object: if it has a density \(p(x)\), then that density is a function on \(\cX\) and also an infinite-dimensional vector. A natural question is: how can we represent this distribution in a useful and manageable way?

There are some straightforward ideas for representing a distribution. For example, we can consider the mean or other moments of the distribution, with the definition

\[ m_P = \E_P[X] = \int x p(x) dx. \]

For some restricted families, such as Gaussian distributions, the mean and variance already characterize the distribution. However, in general, the mean and variance are not enough. Another straightforward idea is to divide the domain \(\cX\) into a finite number of bins or regions and represent the distribution by the probability masses in those bins. This is essentially what a histogram does. More precisely, if we partition \(\cX\) into \(m\) disjoint regions \(A_1,\ldots,A_m\), we can record the probability mass within each region:

\[ \begin{aligned} P(A_i) &= \int 1_{A_i}(x) p(x) dx \\ &= \int_{x \in A_i} p(x) dx \\ &= P(X \in A_i) \end{aligned} \]

for each \(i=1,\ldots,m\). The vector \(\big(P(A_1),\ldots,P(A_m)\big)^\T\) is a finite-dimensional representation of \(P\). This representation is restricted, and it can be difficult to compare two different definitions of the bins because their dimensions may not even match. So how about we define something more flexible that also lives in a common space, say an RKHS?

2 Partition Kernel

In the previous example, let’s view the partition from a different angle. What if we represent these probability masses as a function on the domain \(\cX\)? More specifically, we can redefine the previous example as

\[ \begin{aligned} \mu_{P,i}(\cdot) &= \int 1_{A_i}(\cdot)1_{A_i}(x) p(x) dx \\ &= \int_{x \in A_i} 1_{A_i}(\cdot) p(x) dx \\ &= \cases{ P(X \in A_i) & if $\,\,$ $\cdot \in A_i$ \cr 0 & otherwise } \end{aligned} \]

for each \(i=1,\ldots,m\). The function \(\mu_{P,i}(\cdot)\) is constant on \(A_i\) and takes the value of the probability within \(A_i\). We can consider all regions together by defining the feature map

\[ \Phi(\cdot) = \big(1_{A_1}(\cdot), \ldots, 1_{A_m}(\cdot) \big)^\T, \]

The induced partition kernel is

\[ K_{\mathcal A}(x,y)=\langle\Phi(x),\Phi(y)\rangle =\sum_{\ell=1}^m1_{A_\ell}(x)1_{A_\ell}(y) =1\{x\text{ and }y\text{ are in the same region}\}. \]

Then, the representation of the distribution can be defined as

\[ \begin{aligned} \mu_P(\cdot) &= \int \langle \Phi(\cdot) , \Phi(x) \rangle p(x) \, dx \\ &= \langle \Phi(\cdot), \int \Phi(x) p(x) \, dx \rangle \\ &= \begin{bmatrix} 1_{A_1}(\cdot) \\ \vdots \\ 1_{A_m}(\cdot) \end{bmatrix}^\top \begin{bmatrix} P(X \in A_1) \\ \vdots \\ P(X \in A_m) \end{bmatrix} \end{aligned} \]

What does this function effectively do? It takes the value \(P(X\in A_i)\) when the input is in \(A_i\). Thus, it is a piecewise-constant function on \(\cX\) that summarizes the probability masses in each region. It embeds the distribution into a space of piecewise-constant functions induced by the partition. If the regions are the terminal nodes of a decision tree, this construction gives a tree-based kernel, which we will discuss later.

3 Kernel Mean Embedding

In general, we know that each feature map \(\Phi\) induces a kernel \(K(x,y)=\langle\Phi(x),\Phi(y)\rangle\). We can therefore define the representation of a distribution \(P\) as

\[ \begin{aligned} \mu_P(\cdot) &= \int K(\cdot,x)p(x)\,dx \\ &= \int K(\cdot,x)\,dP(x) \\ &= \E_{X\sim P}[K(\cdot,X)] \end{aligned} \]

For this RKHS-valued expectation to be well-defined, it is enough to assume \(\E_{X\sim P}[\sqrt{K(X,X)}]<\infty\). This condition holds automatically for bounded kernels, including the Gaussian kernel.

This is called the kernel mean embedding of \(P\) into the RKHS \(\cH\) associated with the positive-definite kernel \(K\). There are some nice properties of this embedding. First, it is linear. For a mixture \(R=aP+bQ\), where \(a,b\geq0\) and \(a+b=1\), we have

\[ \begin{aligned} \mu_{aP+bQ}(\cdot) &= \int K(\cdot,x)\,d(aP+bQ)(x) \\ &= a\int K(\cdot,x)\,dP(x)+b\int K(\cdot,x)\,dQ(x) \\ &= a \mu_P(\cdot) + b \mu_Q(\cdot) \end{aligned} \]

Second, it satisfies the mean reproducing property. For any \(f\in\cH\), we have

\[ \begin{aligned} \langle f,\mu_P\rangle_\cH &= \left\langle f,\int K(\cdot,x)\,dP(x)\right\rangle_\cH \\ &= \int\langle f,K(\cdot,x)\rangle_\cH\,dP(x) \\ &= \int f(x) \,dP(x) \\ &= \E_{X \sim P} [f(X)] \end{aligned} \]

Here the term reproducing refers to the fact that the inner product with the mean embedding \(\mu_P\) “reproduces” the expectation of the function \(f\) under the distribution \(P\).

4 Example: Linear Kernel

Let’s take the linear kernel \(K(x,y)=x^\top y\) as an example, assuming \(\E\|X\|<\infty\). The associated RKHS is the space of linear functions, i.e., \(\cH=\{f(x)=w^\top x:w\in\RR^d\}\). The kernel mean embedding of a distribution \(P\) is

\[ \begin{aligned} \mu_P(\cdot) &= \int K(\cdot,x)\,dP(x) \\ &= \int \cdot^\top x \, dP(x) \\ &= \cdot^\top \int x \, dP(x) \\ &= \cdot^\top \E_{X \sim P}[X] \end{aligned} \]

Thus, the kernel mean embedding is a linear function whose coefficient is the mean of \(P\). Evaluating it at any input \(t\) gives the expectation of the linear combination \(t^\top X\).

5 Example: Quadratic Kernel

We know that a quadratic kernel can be written as \(K(x,y)=(x^\top y)^2\). Assuming \(\E\|X\|^2<\infty\), one associated feature map is

\[ \Phi(x) = (x_i x_j)_{1 \leq i,j \leq d} \in \RR^{d^2} \]

Hence the kernel mean embedding of a distribution \(P\) is

\[ \begin{aligned} \mu_P(\cdot) &= \int K(\cdot,x)\,dP(x) \\ &= \int (\cdot^\top x)^2 \, dP(x) \\ &= \int (\cdot^\top x) (x^\top \cdot) \, dP(x) \\ &= \int \cdot^\top (x x^\top) \cdot \, dP(x) \\ &= \cdot^\top \Big( \int x x^\top \, dP(x) \Big) \cdot \\ &= \cdot^\top \E_{X \sim P}[X X^\top] \cdot \end{aligned} \]

Remember that the term in the middle is the \(d\times d\) second-moment matrix of \(P\). Evaluating the mean embedding at any input \(t\) gives the expectation of the quadratic form \(t^\top XX^\top t=(t^\top X)^2\) when \(X\sim P\).

6 Example: Gaussian Kernel

For this example, let’s visualize the kernel mean embedding of a Uniform distribution \(P=\operatorname{Unif}(-1,1)\) on \(\RR\) with the Gaussian kernel \(K(x,y)=\exp\{-\|x-y\|^2/(2\sigma^2)\}\). The kernel mean embedding is

\[ \begin{aligned} \mu_P(\cdot) &= \int K(\cdot,x)\,dP(x) \\ &= \int \exp\Big(-\frac{\|\cdot - x\|^2}{2\sigma^2}\Big) \, dP(x) \\ &= \int_{-1}^1 \frac{1}{2} \exp\Big(-\frac{\|\cdot - x\|^2}{2\sigma^2}\Big) \, dx \end{aligned} \]

Hence, the embedding is the convolution of the Uniform density and the unnormalized Gaussian kernel. For an input \(t\), it simplifies to

\[ \mu_P(t)=\frac{\sigma\sqrt{2\pi}}{2} \left\{\Phi_0\left(\frac{t+1}{\sigma}\right) -\Phi_0\left(\frac{t-1}{\sigma}\right)\right\}, \]

where \(\Phi_0\) is the standard Gaussian CDF. Below is the plot of the embedding for \(\sigma=1,0.4,\) and \(0.1\), together with the original Uniform density.

  library(ggplot2)

  x_seq <- seq(-3, 3, length.out = 200)
  mu_P <- function(t, sigma) {
    sigma * sqrt(2 * pi) * (pnorm(t + 1, sd = sigma) - pnorm(t - 1, sd = sigma)) / 2
  }

  df <- data.frame(x = x_seq, 
                   y1 = mu_P(x_seq, sigma = 1),
                   y2 = mu_P(x_seq, sigma = 0.4),
                   y3 = mu_P(x_seq, sigma = 0.1))
  
  p <- ggplot(df) + 
    geom_line(aes(x = x, y = y1), color = 'blue', linewidth = 1) +
    geom_line(aes(x = x, y = y2), color = 'darkorange', linewidth = 1) +
    geom_line(aes(x = x, y = y3), color = 'darkgreen', linewidth = 1) +
    geom_area(data = subset(df, x >= -1 & x <= 1),
              aes(x = x, y = 0.5), fill = "gray", alpha = 0.5) +
    labs(x = 'x', y = expression(mu[P](x)), 
         title = 'KME of Uniform(-1,1) with Gaussian Kernel') +
    theme_minimal() +
    theme(text = element_text(size=16)) +
    ylim(0, 0.9) +
    annotate("text", x = 2, y = 0.83, label = "sigma == 1",
             parse = TRUE, color = 'blue', size = 5) +
    annotate("text", x = 2, y = 0.75, label = "sigma == 0.4",
             parse = TRUE, color = 'darkorange', size = 5) + 
    annotate("text", x = 2, y = 0.67, label = "sigma == 0.1",
             parse = TRUE, color = 'darkgreen', size = 5)
  
  print(p)

Let’s then look at the sample version. If we observe \(\{x_i\}_{i=1}^n\) from the \(\operatorname{Unif}(-1,1)\) distribution, then the empirical kernel mean embedding is

\[ \widehat{\mu}_P(\cdot)=\frac{1}{n}\sum_{i=1}^n K(\cdot,x_i). \]

Below is the empirical embedding based on \(n=50\) samples.

  set.seed(546)
  n <- 50
  
  x_sample <- runif(n, min = -1, max = 1)
  hat_mu <- function(t, x_sample, sigma) {
    sapply(t, function(x) mean(exp(-(x - x_sample)^2 / (2 * sigma^2))))
  }
  
  df2 <- data.frame(x = x_seq, 
                   y1 = hat_mu(x_seq, x_sample, sigma = 1),
                   y2 = hat_mu(x_seq, x_sample, sigma = 0.4),
                   y3 = hat_mu(x_seq, x_sample, sigma = 0.1))
  
  p2 <- ggplot(df2) +
    geom_line(aes(x = x, y = y1), color = 'blue', linewidth = 1) +
    geom_line(aes(x = x, y = y2), color = 'darkorange', linewidth = 1) +
    geom_line(aes(x = x, y = y3), color = 'darkgreen', linewidth = 1) +
    geom_area(data = subset(df2, x >= -1 & x <= 1),
              aes(x = x, y = 0.5), fill = "gray", alpha = 0.5) +
    labs(x = 'x', y = expression(hat(mu)[P](x)), 
         title = paste('Empirical KME of Uniform(-1,1)')) +
    theme_minimal() +
    theme(text = element_text(size=16)) +
    ylim(0, 0.9) +
    annotate("text", x = 2, y = 0.83, label = "sigma == 1",
             parse = TRUE, color = 'blue', size = 5) +
    annotate("text", x = 2, y = 0.75, label = "sigma == 0.4",
             parse = TRUE, color = 'darkorange', size = 5) + 
    annotate("text", x = 2, y = 0.67, label = "sigma == 0.1",
             parse = TRUE, color = 'darkgreen', size = 5)
  
  print(p2)

But isn’t this just a kernel density estimate (KDE)?

  • With the Gaussian kernel used here, the empirical KME is proportional to a Gaussian KDE for a fixed \(\sigma\). A KDE includes the normalizing factor \(1/(\sqrt{2\pi}\sigma)\), whereas our Gaussian RKHS kernel does not.
  • The main distinction is the target. A KME represents \(P\) as an element of an RKHS and can use any positive-definite kernel. A KDE estimates a density and uses a kernel that integrates to one; it need not be positive-definite.
  • If the RKHS kernel is characteristic, the KME determines the entire distribution. It can then support tasks such as two-sample testing without reconstructing a density.

7 Characteristic Kernels

A natural question to ask is, can we recover the original distribution \(P\) from its kernel mean embedding \(\mu_P\)? In general, the answer is no. This is easy to understand because a linear or quadratic kernel captures only the first or second moments, which do not determine a general distribution. However, if the kernel \(K\) is characteristic, then the answer is yes. A positive-definite kernel \(K\) is characteristic if the map \(P\mapsto\mu_P\) is injective, i.e.,

\[ \mu_P = \mu_Q \Longleftrightarrow P = Q. \]

This means that two distributions have the same embedding only when they are the same distribution. Many commonly used kernels are characteristic, including the Gaussian and Laplace kernels. The following optional Fourier argument explains why the Gaussian kernel is characteristic. Consider a bounded, continuous, real-valued, translation-invariant kernel on \(\RR^d\) of the form \[ K(x,y)=\psi(x-y), \] Bochner’s theorem1 states that \(K\) is positive-definite if and only if \(\psi\) is the Fourier transform of a finite nonnegative measure \(\Lambda\): \[ K(x,y)=\int_{\RR^d}e^{i\omega^\top(x-y)}\,d\Lambda(\omega). \]

Define the characteristic function of \(P\) by

\[ \varphi_P(\omega)=\E_{X\sim P}\left[e^{i\omega^\T X}\right]. \]

Using Bochner’s representation, the kernel mean embedding can be expressed in Fourier space as

\[ \begin{aligned} \mu_P(t) &=\int_{\RR^d}\!\!\left\{\int_{\RR^d}e^{i\omega^\T(t-z)}\,d\Lambda(\omega)\right\}dP(z)\\ &=\int_{\RR^d}e^{i\omega^\T t}\varphi_P(-\omega)\,d\Lambda(\omega). \end{aligned} \]

The squared distance between two embeddings has the corresponding frequency-domain form

\[ \|\mu_P-\mu_Q\|_\cH^2 =\int_{\RR^d}|\varphi_P(\omega)-\varphi_Q(\omega)|^2\,d\Lambda(\omega). \]

Suppose two distributions have the same embedding. The distance identity then gives

\[ \mu_P=\mu_Q \quad\Longrightarrow\quad \varphi_P(\omega)=\varphi_Q(\omega) \quad\text{for $\Lambda$-almost every }\omega. \]

For a bounded, continuous, translation-invariant kernel on \(\RR^d\), the kernel is characteristic if and only if \(\operatorname{supp}(\Lambda)=\RR^d\) (Sriperumbudur et al. 2010). Indeed, full support and continuity of characteristic functions extend the almost-everywhere equality to all \(\omega\), and uniqueness of characteristic functions then gives \(P=Q\).

For the Gaussian kernel, the spectral measure has density proportional to

\[ \exp\left\{-\frac{\sigma^2\|\omega\|^2}{2}\right\}, \]

which is positive everywhere. Its support is therefore the entire frequency space, which shows why the Gaussian kernel is characteristic.

This establishes injectivity, but numerically recovering a distribution from its embedding is still a difficult inverse problem. The main purpose of KME is to represent a distribution in the RKHS and use it for other tasks. The following is one of its most important applications.

8 Application: Maximum Mean Discrepancy

A key application of KME is to compare two distributions. If we have distributions \(P\) and \(Q\), we can first calculate their kernel mean embeddings \(\mu_P\) and \(\mu_Q\). Then we compare these two embeddings in the RKHS in the sense that

\[ \mathrm{MMD}(P,Q;\cH) \;=\; \|\mu_P - \mu_Q\|_{\cH}. \]

This is called the maximum mean discrepancy (MMD) between \(P\) and \(Q\) with respect to the RKHS \(\cH\). If the kernel \(K\) is characteristic, then \(\mathrm{MMD}(P,Q;\cH)=0\) if and only if \(P=Q\), so MMD is a metric on probability distributions. Even if the kernel is not characteristic, MMD remains a useful discrepancy measure. An equivalent variational form is

\[ \mathrm{MMD}(P,Q;\cH)=\sup_{\|f\|_{\cH}\leq1}\left\{\E_{X\sim P}[f(X)]-\E_{Y\sim Q}[f(Y)]\right\}. \]

This means that, if we are interested in functions of the variable, say \(f(X)\), then MMD is the largest difference in expectation over the RKHS unit ball. Because this ball is symmetric, including an absolute value gives the same supremum. If the two distributions are very different, then some function \(f\) can amplify this difference. We constrain \(f\) to have RKHS norm at most 1 to avoid the trivial solution of making \(f\) arbitrarily large. Let’s use the mean reproducing property:

\[ \begin{aligned} \E_{X \sim P}[f(X)] - \E_{Y \sim Q}[f(Y)] &= \langle f, \mu_P \rangle_{\cH} - \langle f, \mu_Q \rangle_{\cH} \\ &= \langle f, \mu_P - \mu_Q \rangle_{\cH} \\ &\leq \|f\|_{\cH} \|\mu_P - \mu_Q\|_{\cH} \end{aligned} \]

When \(\mu_P\neq\mu_Q\), the supremum is attained when \(f\) is the unit-norm function in the direction of \(\mu_P-\mu_Q\):

\[ \begin{aligned} &\sup_{\|f\|_{\cH}\leq1}\left\{\E_{X\sim P}[f(X)]-\E_{Y\sim Q}[f(Y)]\right\} \\ =& \langle \frac{\mu_P - \mu_Q}{\|\mu_P - \mu_Q\|_{\cH}}, \mu_P - \mu_Q \rangle_{\cH} \\ =& \|\mu_P - \mu_Q\|_{\cH}. \end{aligned} \]

If \(\mu_P=\mu_Q\), both sides are zero. To compute MMD in practice, we use the squared RKHS norm and the reproducing property:

\[ \begin{aligned} \mathrm{MMD}^2(P,Q;\cH) &= \|\mu_P-\mu_Q\|_{\cH}^2 \\ &= \langle \mu_P - \mu_Q, \mu_P - \mu_Q \rangle_{\cH} \\ &=\langle \mu_P, \mu_P \rangle_{\cH} + \langle \mu_Q, \mu_Q \rangle_{\cH} - 2 \langle \mu_P, \mu_Q \rangle_{\cH}\\ &= \left\langle \int K(\cdot,x)\,dP(x),\int K(\cdot,x')\,dP(x')\right\rangle_{\cH} + \left\langle \int K(\cdot,y)\,dQ(y),\int K(\cdot,y')\,dQ(y')\right\rangle_{\cH} \\ &\quad -2\left\langle \int K(\cdot,x)\,dP(x),\int K(\cdot,y)\,dQ(y)\right\rangle_{\cH} \\ &= \iint K(x,x')\,dP(x)dP(x')+\iint K(y,y')\,dQ(y)dQ(y') \\ &\quad -2\iint K(x,y)\,dP(x)dQ(y) \\ &= \E_{X,X'\sim P}[K(X,X')]+\E_{Y,Y'\sim Q}[K(Y,Y')]-2\E_{X\sim P,\,Y\sim Q}[K(X,Y)] \end{aligned} \]

Here, \(X\) and \(X'\) are independent draws from \(P\), \(Y\) and \(Y'\) are independent draws from \(Q\), and the two pairs are mutually independent.

Given i.i.d. samples \(\{x_i\}_{i=1}^n\sim P\) and \(\{y_j\}_{j=1}^m\sim Q\), the empirical embeddings give the biased estimator of the squared MMD

\[ \widehat{\mathrm{MMD}^2}_{b} =\frac{1}{n^2}\sum_{i,i'}K(x_i,x_{i'}) +\frac{1}{m^2}\sum_{j,j'}K(y_j,y_{j'}) -\frac{2}{mn}\sum_{i,j}K(x_i,y_j). \]

This is exactly \(\|\widehat{\mu}_P-\widehat{\mu}_Q\|_\cH^2\), so it is always nonnegative. Removing the diagonal terms and renormalizing gives an unbiased estimator of the population squared MMD when \(n,m\geq2\):

\[ \widehat{\mathrm{MMD}^2}_{u} =\frac{1}{n(n-1)}\sum_{i\neq i'}K(x_i,x_{i'})+\frac{1}{m(m-1)}\sum_{j\neq j'}K(y_j,y_{j'})-\frac{2}{mn}\sum_{i,j}K(x_i,y_j). \]

Although unbiased for \(\mathrm{MMD}^2(P,Q;\cH)\), \(\widehat{\mathrm{MMD}^2}_{u}\) can be negative in a finite sample.

Let’s look at an empirical example of using MMD to compare the same \(\operatorname{Unif}(-1,1)\) distribution with a Gaussian \(\cN(0,1)\) alternative. We use the Gaussian kernel with bandwidth \(\sigma=0.3\). Below is a plot of the two samples.

  set.seed(123)
  
  # Gaussian kernel
  kfun <- function(x, y, sigma = 1) {
    exp(-(outer(x, y, "-")^2) / (2 * sigma^2))
  }
  
  mmd2_unbiased <- function(x, y, sigma = 1) {
    n <- length(x)
    m <- length(y)

    if (n < 2 || m < 2) {
      stop("Both samples must contain at least two observations.")
    }

    Kxx <- kfun(x, x, sigma)
    Kyy <- kfun(y, y, sigma)
    Kxy <- kfun(x, y, sigma)

    term_xx <- (sum(Kxx) - sum(diag(Kxx))) / (n * (n - 1))
    term_yy <- (sum(Kyy) - sum(diag(Kyy))) / (m * (m - 1))
    term_xy <- mean(Kxy)

    term_xx + term_yy - 2 * term_xy
  }
  
  # Example: Gaussian vs Uniform
  n_x <- n_y <- 200
  x <- rnorm(n_x, mean = 0, sd = 1)      # Gaussian N(0,1)
  y <- runif(n_y, min = -1, max = 1)     # Uniform(-1,1)
  
  # plot the two sets of samples
  df_samples <- data.frame(
    value = c(x, y),
    group = rep(c('Gaussian', 'Uniform'), times = c(length(x), length(y))))
  p_samples <- ggplot(df_samples, aes(x = value, fill = group)) +
    geom_histogram(aes(y = after_stat(density)), position = 'identity', alpha = 0.5, bins = 30) +
    labs(x = 'Value', y = 'Density', title = 'Samples from Two Distributions') +
    theme_minimal() +
    theme(text = element_text(size=16)) +
    scale_fill_manual(values = c('Gaussian' = 'darkorange', 'Uniform' = 'lightblue'))
  print(p_samples)

  
  
  mmd2_val <- mmd2_unbiased(x, y, sigma = 0.3)
  print(mmd2_val)
## [1] 0.02625334

But how do we know if this value is large or small? We test

\[ H_0:P=Q \qquad\text{against}\qquad H_1:P\neq Q. \]

Under \(H_0\), the pooled observations are exchangeable, so permuting the sample labels gives the null distribution of the statistic.

  set.seed(1234)
  
  nperm <- 500
  mmd2_perm <- numeric(nperm)
  xy <- c(x, y)
  for (b in 1:nperm) {
    perm_idx <- sample(seq_len(n_x + n_y))
    x_perm <- xy[perm_idx[seq_len(n_x)]]
    y_perm <- xy[perm_idx[n_x + seq_len(n_y)]]
    mmd2_perm[b] <- mmd2_unbiased(x_perm, y_perm, sigma = 0.3)
  }
  
  p_value <- (1 + sum(mmd2_perm >= mmd2_val)) / (nperm + 1)
  print(p_value)
## [1] 0.003992016
  
  df_mmd <- data.frame(mmd2 = mmd2_perm)
  
  p_mmd <- ggplot(df_mmd, aes(x = mmd2)) +
    geom_histogram(aes(y = after_stat(density)), bins = 30, fill = 'lightblue', color = 'black') +
    geom_vline(xintercept = mmd2_val, color = 'red', linewidth = 1) +
    labs(x = expression(widehat(MMD^2)[u]),
         y = 'Density', 
         title = paste('Permutation Test for Squared MMD (p-value =', round(p_value, 4), ')')) +
    theme_minimal() +
    theme(text = element_text(size=16)) +
    annotate("text", x = mmd2_val - 0.005, y = max(density(mmd2_perm)$y) * 0.8,
             label = "Observed unbiased estimator of MMD^2", color = 'red', size = 5)
  
  print(p_mmd)

Negative permutation values are possible because \(\widehat{\mathrm{MMD}^2}_{u}\) is unbiased for \(\mathrm{MMD}^2(P,Q;\cH)\) but is not itself a squared RKHS norm.

MMD is widely used for two-sample testing. Related kernel embedding methods are also used for independence and conditional independence testing. MMD is also used in generative models such as MMD-GAN (Gretton et al. 2012; Li et al. 2017; Muandet et al. 2017).


Gretton, Arthur, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. 2012. “A Kernel Two-Sample Test.” Journal of Machine Learning Research 13 (25): 723–73. https://www.jmlr.org/papers/v13/gretton12a.html.
Li, Chun-Liang, Wei-Cheng Chang, Yu Cheng, Yiming Yang, and Barnabás Póczos. 2017. “Mmd Gan: Towards Deeper Understanding of Moment Matching Network.” Advances in Neural Information Processing Systems 30.
Muandet, Krikamol, Kenji Fukumizu, Bharath Sriperumbudur, and Bernhard Schölkopf. 2017. “Kernel Mean Embedding of Distributions: A Review and Beyond.” Foundations and Trends® in Machine Learning 10 (1-2): 1–141. https://doi.org/10.1561/2200000060.
Sriperumbudur, Bharath K., Arthur Gretton, Kenji Fukumizu, Bernhard Schölkopf, and Gert R. G. Lanckriet. 2010. “Hilbert Space Embeddings and Metrics on Probability Measures.” Journal of Machine Learning Research 11 (50): 1517–61. https://www.jmlr.org/papers/v11/sriperumbudur10a.html.

  1. Bochner’s theorem gives a spectral representation for bounded, continuous, translation-invariant kernels on \(\RR^d\). A Mercer expansion instead depends on a reference measure and additional regularity conditions.↩︎