def MMD2u_estimator(K, m, n):
    """
    Compute the MMD^2_u unbiased statistic.
    This is an implementation of an unbiased MMD^2_u estimator:
    Equation (3) in Gretton et al. Journal of Machine Learning Research 13 (2012) 723-773

    :param K: numpy-array
        the pair-wise kernel matrix
    :param m: int
        dimension of the first data set
    :param n: int
        dimension of the second data set

    :return: float
        an unbiased estimate of MMD^2_u
    """
    K_x = K[:m, :m]
    K_y = K[m:, m:]
    K_xy = K[:m, m:]
    return 1.0 / (m * (m - 1.0)) * (K_x.sum() - K_x.diagonal().sum()) + \
            1.0 / (n * (n - 1.0)) * (K_y.sum() - K_y.diagonal().sum()) - \
            2.0 / (m * n) * K_xy.sum()