Learning Theory for Distribution Regression

Introduction

We address the learning problem of distribution regression in the two-stage sampled setting, where we only have bags of samples from the probability distributions: we regress from probability measures to real-valued (Póczos et al., 2013) responses, or more generally to vector-valued outputs (belonging to an arbitrary separable Hilbert space). Many classical problems in machine learning and statistics can be analysed in this framework. On the machine learning side, multiple instance learning (Dietterich et al., 1997; Ray and Page, 2001; Dooly et al., 2002) can be thought of in this way, where each instance in a labeled bag is an i.i.d. (independent identically distributed) sample from a distribution. On the statistical side, tasks might include point estimation of statistics on a distribution without closed form analytical expressions (e.g., its entropy or a hyperparameter).

Can the distribution regression problem be solved consistently under mild conditions? What is the exact computational-statistical efficiency trade-off implied by the two-stage sampling?

In our work the estimated mapping ($\hat{f}$) is the analytical solution of a kernel ridge regression (KRR) problem.Beyond its simple analytical formula, kernel ridge regression also allows efficient distributed (Zhang et al., 2015; Richtárik and Takác̆, 2016), sketch (Alaoui and Mahoney, 2015; Yang et al., 2016) and Nyström based approximations (Rudi et al., 2015). The performance of $\hat{f}$ depends on the assumed function class ($\mathscr{H}$), the family of $\hat{f}$ candidates used in the ridge formulation. We shall focus on the analysis of two settings:

Well-specified case ($f_{*}\in\mathscr{H}$): In this case we assume that the regression function $f_{*}$ belongs to $\mathscr{H}$. We focus on bounding the goodness of $\hat{f}$ compared to $f_{*}$. In other words, if $\mathcal{R}[f_{*}]$ denotes the prediction error (expected risk) of $f_{*}$, then our goal is to derive a finite-sample bound for the excess risk, $\mathcal{E}(\hat{f},f_{*})=\mathcal{R}[\hat{f}]-\mathcal{R}[f_{*}]$ that holds with high probability. We make use of this bound to establish the consistency of the estimator (i.e., drive the excess risk to zero) and to derive the exact computational-statistical efficiency trade-off of the estimator as a function of the sample number ($l$, $N=N_{i},\,\forall\,i$) and the problem difficulty (see Theorem 5 and its corresponding remarks for more details).

Misspecified case ($f_{*}\in L^{2}\backslash\mathscr{H}$): Since in practise it might be hard to check whether $f_{*}\in\mathscr{H}$, we also study the misspecified setting of $f_{*}\in L^{2}$; the relevant case is when $f_{*}\in L^{2}\backslash\mathscr{H}$. In the misspecified setting the ’richness’ of $\mathscr{H}$ has crucial importance, in other words the size of $D_{\mathscr{H}}^{2}=\inf_{f\in\mathscr{H}}\|f_{*}-f\|_{2}^{2}$, the approximation error from $\mathscr{H}$. Our main contributions consist of proving a finite-sample excess risk bound, using which we show that the proposed estimator can attain the ideal performance, i.e., $\mathcal{E}(\hat{f},f_{*})-D_{\mathscr{H}}^{2}$ can be driven to zero. Moreover, on smooth classes of $f_{*}$-s, we give a simple and explicit description for the computational-statistical efficiency trade-off of our estimator (see Theorem 9 and its corresponding remarks for more details).

One can define kernel learning algorithms on bags based on set kernels (Gärtner et al., 2002) by computing the similarity of the sets/bags of samples representing the input distributions; set kernels are also called called multi-instance kernels or ensemble kernels, and are examples of convolution kernels (Haussler, 1999). In this case, the similarity of two sets is measured by the average pairwise point similarities between the sets. From a theoretical perspective, nothing is known about the consistency of set kernel based learning method since their introduction in 1999 (Haussler, 1999; Gärtner et al., 2002): i.e. in what sense (and with what rates) is the learning algorithm consistent, when the number of items per bag, and the number of bags, are allowed to increase?

It is possible, however, to view set kernels in a distribution setting, as they represent valid kernels between (mean) embeddings of empirical probability measures into a reproducing kernel Hilbert space (RKHS; Berlinet and Thomas-Agnan, 2004). The population limits are well-defined as being dot products between the embeddings of the generating distributions (Altun and Smola, 2006), and for characteristic kernels the distance between embeddings defines a metric on probability measures (Sriperumbudur et al., 2011; Gretton et al., 2012). When bounded kernels are used, mean embeddings exist for all probability measures (Fukumizu et al., 2004). When we consider the distribution regression setting, however, there is no reason to limit ourselves to set kernels. Embeddings of probability measures to RKHS are used by Christmann and Steinwart (2010) in defining a yet larger class of easily computable kernels on distributions, via operations performed on the embeddings and their distances. Note that the relation between set kernels and kernels on distributions was also applied by Muandet et al. (2012) for classification on distribution-valued inputs, however consistency was not studied in that work. We also note that motivated by the current paper, Lopez-Paz et al. (2015) have recently presented the first theoretical results about surrogate risk guarantees on a class (relying on uniformly bounded Lipschitz functionals) of soft distribution-classification problems.

Our contribution in this paper is to establish the learning theory of a simple, mean embedding based ridge regression (MERR) method for the distribution regression problem. This result applies both to the basic set kernels of Haussler (1999); Gärtner et al. (2002), the distribution kernels of Christmann and Steinwart (2010), and additional related kernels. We provide finite-sample excess risk bounds, prove consistency, and show how the two-stage sampled nature of the problem (bag size) governs the computational-statistical efficiency of the MERR estimator. More specifically, in the

derive finite-sample bounds on the excess risk: We construct $\mathcal{R}[\hat{f}]-\mathcal{R}[f_{*}]\leq r(l,N,\lambda)$ bounds holding with high probability, where $\lambda$ is the regularization parameter in the ridge problem ($\lambda\rightarrow 0$, $l\rightarrow\infty$, $N=N_{i}\rightarrow\infty$).

establish consistency and computational-statistical efficiency trade-off of the MERR estimator on a general prior family $\mathcal{P}(b,c)$ as defined by Caponnetto and De Vito (2007), where $b$ captures the effective input dimension, and larger $c$ means smoother $f_{*}$ ($1<b$, $c\in(1,2]$). In particular, when the number of samples per bag is chosen as $N=l^{a}\log(l)$ and $a\geq\frac{b(c+1)}{bc+1}$, then the learning rate saturates at $l^{-\frac{bc}{bc+1}}$, which is known to be one-stage sampled minimax optimal (Caponnetto and De Vito, 2007). In other words, by choosing $a=\frac{b(c+1)}{bc+1}<2$, we suffer no loss in statistical performance compared with the best possible one-stage sampled estimator.

Note: the advantage of considering the $\mathcal{P}(b,c)$ family is two-fold. It does not assume parametric distributions, yet certain complexity terms can be explicitly upper bounded in the family. This property will be exploited in our analysis. Moreover, (for special input distributions) the parameter $b$ can be related to the spectral decay of Gaussian Gram matrices, and existing analysis techniques (Steinwart and Christmann, 2008) may be used in interpreting these decay conditions.

misspecified case: We establish consistency and convergence rates even if $f_{*}\notin\mathscr{H}$. Particularly, by deriving finite-sample bounds on the excess risk we

prove that the MERR estimator can achieve the best possible approximation accuracy from $\mathscr{H}$, i.e. the $\mathcal{R}[\hat{f}]-\mathcal{R}[f_{*}]-D_{\mathscr{H}}^{2}$ quantity can be driven to zero (recall that $D_{\mathscr{H}}=\inf_{f\in\mathscr{H}}\|f_{*}-f\|_{2}$). Specifically, this result implies that if $\mathscr{H}$ is dense in $L^{2}$ ($D_{\mathscr{H}}=0$), then the excess risk $\mathcal{R}[\hat{f}]-\mathcal{R}[f_{*}]$ converges to zero.

Due to the differences in the assumptions made and the loss function used, a direct comparison of our theoretical result and that of Póczos et al. (2013) remains an open question, however we make three observations. First, our approach is more general, since we may regress from any probability measure defined on separable, topological domains endowed with kernels. Póczos et al.’s work is restricted to compact domains of finite dimensional Euclidean spaces, and requires the distributions to admit probability densities; distributions on strings, graphs, and other structured objects are disallowed. Second, in our analysis we will allow separable Hilbert space valued outputs, in contrast to the real-valued output considered by Póczos et al. (2013). Third, density estimates in high dimensional spaces suffer from slow convergence rates (Wasserman, 2006, Section 6.5). Our approach mitigates this problem, as it works directly on distribution embeddings, and does not make use of density estimation as an intermediate step.

The principal challenge in proving theoretical guarantees arises from the two-stage sampled nature of the inputs. In our analysis of the well-specified case, we make use of Caponnetto and De Vito (2007)’s results, which focus (only) on the one-stage sample setup. These results will make our analysis somewhat shorter (but still rather challenging) by giving upper bounds for some of the objective terms. Even the verification of these conditions requires care since the inputs in the ridge regression are themselves distribution embeddings (i.e., functions in a reproducing kernel Hilbert space).

In the misspecified case, RKHS methods alone are not sufficient to obtain excess risk bounds: one has to take into account the “richness” of the modelling RKHS class ($\mathscr{H}$) in the embedding $L^{2}$ space. The fundamental challenge is whether it is possible to achieve the best possible performance dictated by $\mathscr{H}$; or in the special case when further smoothness conditions hold on $f_{*}$, what convergence rates can yet be attained, and what computational-statistical efficiency trade-off realized. The second smoothness property could be modelled for example by range spaces of (fractional) powers of integral operators associated to $\mathscr{H}$. Indeed, there exist several results along these lines with KRR for the case of real-valued outputs: see for example (Sun and Wu, 2009a, Theorem 1.1), (Sun and Wu, 2009b, Corollary 3.2), (Mendelson and Neeman, 2010, Theorem 3.7 with Assumption 3.2). The question of optimal rates has also been addressed for the semi-supervised KRR setting (Caponnetto, 2006, Theorem 1), and for clipped KRR estimators (Steinwart et al., 2009) with integral operators of rapidly decaying spectrum. Our results apply more generally to the two-stage sampled setting and to vector valued outputs belonging to separable Hilbert spaces. Moreover, we obtain a general consistency result without range space assumptions, showing that the modelling power of $\mathscr{H}$ can be fully exploited, and convergence to the best approximation available from $\mathscr{H}$ can be realized.Specializing our result, we get explicit rates and an exact computational-statistical efficiency description for MERR as a function of sample numbers and problem difficulty, for smooth regression functions.

There are numerous areas in machine learning and statistics, where estimating vector-valued functions has crucial importance. Often in statistics, one is not only confronted with the estimation of a scalar parameter, but with a vector of parameters. On the machine learning side, multi-task learning (Evgeniou et al., 2005), functional response regression (Kadri et al., 2016), or structured output prediction (Brouard et al., 2011; Kadri et al., 2013) fall under the same umbrella: they can be naturally phrased as learning vector-valued functions (Micchelli and Pontil, 2005). The idea underlying all these tasks is simple and intuitive: if multiple prediction problems have to be solved simultaneously, it might be beneficial to exploit their dependencies. Imagine for example that the task is to predict the motion of a dancer: taking into account the interrelation of the actor’s body parts is likely to lead to more accurate estimation, as opposed to predicting the individual parts one by one, independently. Successful real-world applications of a multi-task approach include for example preference modelling of users with similar demographics (Evgeniou et al., 2005), prediction of the daily precipitation profiles of weather stations (Kadri et al., 2010), acoustic-to-articulatory speech inversion (Kadri et al., 2016), identifying biomarkers capable of tracking the progress of Alzheimer’s disease (Zhou et al., 2013), personalized human activity recognition based on iPod/iPhone accelerometer data (Sun et al., 2013), finger trajectory prediction in brain-computer interfaces (Kadri et al., 2012) or ecological inference (Flaxman et al., 2015); for a recent review on multi-output prediction methods see (Álvarez et al., 2011; Borchani et al., 2015). A mathematically sound way of encoding prior information about the relation of the outputs can be realized by operator-valued kernels and the associated vector-valued RKHS-s (Pedrick, 1957; Micchelli and Pontil, 2005; Carmeli et al., 2006, 2010); this is the tool we use to allow vector-valued learning tasks.

Finally, we note that the current work extends our earlier conference paper (Szabó et al., 2015) in several important respects: we now show that the MERR method can attain the one-stage sampled minimax optimal rate; we generalize the analysis in the well-specified setting to allow outputs belonging to an arbitrary separable Hilbert spaces (in contrast to the original scalar-valued output domain); and we tackle the misspecified setting, obtaining finite sample guarantees, consistency, and computational-statistical efficiency trade-offs.

The paper is structured as follows: The distribution regression problem and the MERR technique are introduced in Section 2. Our assumptions are detailed in Section 3. We present our theoretical guarantees (finite-sample bounds on the excess risk, consistency, computational-statistical efficiency trade-offs) in Section 4: the well-specified case is considered in Section 4.1, and the misspecified setting is the focus of Section 4.2. Section 5 is devoted to an overview of existing heuristics for learning on distributions. Conclusions are drawn in Section 6. Section 7 contains proof details. In Section 8 we discuss our assumptions with concrete examples.

The Distribution Regression Problem

Below we first introduce our notation (Section 2.1), then formally define the distribution regression task (Section 2.2).

We use the following notations throughout the paper:

Functional analysis: Let $(N_{1},\left\|\cdot\right\|_{N_{1}})$ and $(N_{2},\left\|\cdot\right\|_{N_{2}})$ denote two normed spaces, then $\mathscr{L}(N_{1},N_{2})$ stands for the space of $N_{1}\rightarrow N_{2}$ bounded linear operators; if $N_{1}=N_{2}$, we will use the $\mathscr{L}(N_{1})=\mathscr{L}(N_{1},N_{2})$ shorthand. For $M\in\mathscr{L}(N_{1},N_{2})$ the operator norm is defined as $\left\|M\right\|_{\mathscr{L}(N_{1},N_{2})}=\sup_{0\neq h\in N_{1}}\left\|Mh\right\|_{N_{2}}/\left\|h\right\|_{N_{1}}$, $Im(M)=\{Mn_{1}\}_{n_{1}\in N_{1}}$ denotes the range of $M$, $Ker(M)=\{n_{1}\in N_{1}:Mn_{1}=0\}$ is the null space of $M$. Let $\mathscr{K}$ be a Hilbert space. The adjoint operator $M^{*}\in\mathscr{L}(\mathscr{K})$ of an operator $M\in\mathscr{L}(\mathscr{K})$ is the operator such that $\left<Ma,b\right>_{\mathscr{K}}=\left<a,M^{*}b\right>_{\mathscr{K}}$ for all $a$ and $b$ in $\mathscr{K}$. $M\in\mathscr{L}(\mathscr{K})$ is called positive if $\left<Ma,a\right>_{\mathscr{K}}\geq 0$ ($\forall a\in\mathscr{K}$), self-adjoint if $M=M^{*}$, and trace class if $\sum_{j\in J}\left<|M|e_{j},e_{j}\right>_{\mathscr{K}}<\infty$ for an $(e_{j})_{j\in J}$ ONB (orthonormal basis) of $\mathscr{K}$ ($|M|:=(M^{*}M)^{\frac{1}{2}}$), in which case $Tr(M):=\sum_{j\in J}\left<Me_{j},e_{j}\right>_{\mathscr{K}}<\infty$; compact if $cl\left[Ma:a\in\mathscr{K},\left\|a\right\|_{\mathscr{K}}\leq 1\right]$ is a compact set. Let $\mathscr{K}_{1}$ and $\mathscr{K}_{2}$ be Hilbert spaces. $M\in\mathscr{L}(\mathscr{K}_{1},\mathscr{K}_{2})$ is called Hilbert-Schmidt if $\left\|M\right\|_{\mathscr{L}_{2}(\mathscr{K}_{1},\mathscr{K}_{2})}^{2}=Tr(M^{*}M)=\sum_{j\in J}\left<Me_{j},Me_{j}\right>_{\mathscr{K}_{2}}<\infty$ for some $(e_{j})_{j\in J}$ ONB of $\mathscr{K}_{1}$. The space of Hilbert-Schmidt operators is denoted by $\mathscr{L}_{2}(\mathscr{K}_{1},\mathscr{K}_{2})=\{M\in\mathscr{L}(\mathscr{K}_{1},\mathscr{K}_{2}):\left\|M\right\|_{\mathscr{L}_{2}(\mathscr{K}_{1},\mathscr{K}_{2})}<\infty\}$. We use the shorthand notation $\mathscr{L}_{2}(\mathscr{K})=\mathscr{L}_{2}(\mathscr{K},\mathscr{K})$ if $\mathscr{K}:=\mathscr{K}_{1}=\mathscr{K}_{2}$; $\mathscr{L}_{2}(\mathscr{K})$ is separable if and only if $\mathscr{K}$ is separable (Steinwart and Christmann, 2008, page 506). Trace class and Hilbert-Schmidt operators over a $\mathscr{K}$ Hilbert space are compact operators (Steinwart and Christmann, 2008, page 505-506); moreover,

the set of mean embeddings (Berlinet and Thomas-Agnan, 2004) of the distributions to the space $H$.The $x\mapsto\mu_{x}$ mapping is defined for all $x\in\mathscr{M}^{+}_{1}(\mathscr{X})$ if $k$ is bounded, i.e., $\sup_{u\in\mathscr{X}}k(u,u)<\infty$. Let $Y$ be a separable Hilbert space, where the inner product is denoted by $\left<\cdot,\cdot\right>_{Y}$; the associated norm is $\left\|\cdot\right\|_{Y}$. $\mathscr{H}=\mathscr{H}(K)$ is the $Y$-valued RKHS (Pedrick, 1957; Micchelli and Pontil, 2005; Carmeli et al., 2006, 2010) of $X\rightarrow Y$ functions with $K:X\times X\rightarrow\mathscr{L}(Y)$ as the reproducing kernel (we will present some concrete examples of $K$ in Section 3; see Table 1); $K_{\mu_{x}}\in\mathscr{L}(Y,\mathscr{H})$ is defined as

Further, $f(\mu_{x})=K_{\mu_{x}}^{*}f$ $(\forall\mu_{x}\in X,f\in\mathscr{H})$.

Regression function: Let $\rho$ be the $\mu$-induced probability measure on the $Z=X\times Y$ product space, and let $\rho(\mu_{x},y)=\rho(y|\mu_{x})\rho_{X}(\mu_{x})$ be the factorization of $\rho$ into conditional and marginal distributions.Our assumptions will guarantee the existence of $\rho$ (see Section 3). Since $Y$ is a Polish space (because it is separable Hilbert) the $\rho(y|\mu_{a})$ conditional distribution ($y\in Y$, $\mu_{a}\in X$) is also well-defined (Steinwart and Christmann, 2008, Lemma A.3.16, page 487). The regression function of $\rho$ with respect to the $(\mu_{x},y)$ pair is denoted by

2 Distribution Regression

We now formally define the distribution regression task. Let us assume that $\mathscr{M}^{+}_{1}(\mathscr{X})$ is endowed with $\mathscr{S}_{1}=\mathcal{B}(\tau_{w})$, the weak-topology generated $\sigma$-algebra; thus $(\mathscr{M}^{+}_{1}(\mathscr{X}),\mathscr{S}_{1})$ is a measurable space. In the distribution regression problem, we are given samples $\hat{\mathbf{z}}=\{(\{x_{i,n}\}_{n=1}^{N_{i}},y_{i})\}_{i=1}^{l}$ with $x_{i,1},\ldots,x_{i,N_{i}}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}x_{i}$ where $\mathbf{z}=\{(x_{i},y_{i})\}_{i=1}^{l}$ with $x_{i}\in\mathscr{M}^{+}_{1}\left(\mathscr{X}\right)$ and $y_{i}\in Y$ drawn i.i.d. from a joint meta distribution $\mathcal{M}$ defined on the measurable space $(\mathscr{M}^{+}_{1}(\mathscr{X})\times Y,\mathscr{S}_{1}\otimes\mathcal{B}(Y))$, the product space enriched with the product $\sigma$-algebra. Unlike in classical supervised learning problems, the problem at hand involves two levels of randomness, wherein first $\mathbf{z}$ is drawn from $\mathcal{M}$, and then $\hat{\mathbf{z}}$ is generated by sampling points from $x_{i}$ for all $i=1,\ldots,l$. The goal is to learn the relation between the random distribution $x$ and response $y$ based on the observed $\hat{\mathbf{z}}$. For notational simplicity, we will assume that $N=N_{i}$ ($\forall i$).

In other words, the distribution $x\in\mathscr{M}^{+}_{1}\left(\mathscr{X}\right)$ is first mapped to $X\subseteq H$ by the mean embedding $\mu$, and the result is composed with $f$, an element of the RKHS $\mathscr{H}$.

which is minimized by the $f_{\rho}$ regression function. The classical regularization approach is to optimize

instead of $\mathcal{R}$, based on samples $\mathbf{z}$. Since $\mathbf{z}$ is not available, we consider the objective function defined by the observable quantity $\hat{\mathbf{z}}$,

where $\hat{x}_{i}=\frac{1}{N}\sum_{n=1}^{N}\delta_{x_{i,n}}$ is the empirical distribution determined by $\left\{x_{i,n}\right\}_{i=1}^{N}$. The ridge regression objective function has an analytical solution: given training samples $\hat{\mathbf{z}}$, the prediction for a new $t$ test distribution is

where $\mathbf{k}=\left[K(\mu_{\hat{x}_{1}},\mu_{t}),\ldots,K(\mu_{\hat{x}_{l}},\mu_{t})\right]\in\mathscr{L}(Y)^{1\times l}$, $\mathbf{K}=[K(\mu_{\hat{x}_{i}},\mu_{\hat{x}_{j}})]\in\mathscr{L}(Y)^{l\times l}$, $[y_{1};\ldots;y_{l}]\in Y^{l}$.

It is important to note that the algorithm has access to the sample points only via their mean embeddings $\{\mu_{\hat{x}_{i}}\}_{i=1}^{l}$ in Eq. (9).

There is a two-stage sampling difficulty to tackle: The transition from $f_{\rho}$ to $f_{\mathbf{z}}^{\lambda}$ represents the fact that we have only $l$ distribution samples ($\mathbf{z}$); the transition from $f_{\mathbf{z}}^{\lambda}$ to $f_{\hat{\mathbf{z}}}^{\lambda}$ means that the $x_{i}$ distributions can be accessed only via samples ($\hat{\mathbf{z}}$).

While ridge regression can be performed using the kernel $K_{\mathscr{G}}$, the two-stage sampling makes it difficult to work with arbitrary $K_{\mathscr{G}}$. By contrast, our choice of $K_{\mathscr{G}}(x,x^{\prime})=K(\mu_{x},\mu_{x^{\prime}})$ enables us to handle the two-stage sampling by estimating $\mu_{x}$ with an empirical estimator, and using it in the algorithm as shown above.

One could also formulate the problem (and get guarantees) for more abstract $X\subseteq H\rightarrow Y$ regression tasks [see Eq. (7)] on a convex set $X$ with $H$ and $Y$ being general, separable Hilbert spaces. Since distribution regression is probably the most accessible example where two-stage sampling appears, and in order to keep the presentation simple, we do not consider such extended formulations in this work.

Our main goals in this paper are as follows: first, to analyse the excess risk

both when $f_{\rho}\in\mathscr{H}$ (the well-specified case) and $f_{\rho}\in L^{2}_{\rho_{X}}\backslash\mathscr{H}$ (the misspecified case); second, to establish consistency ($\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)\rightarrow 0$, or in the misspecified case $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)-D_{\mathscr{H}}^{2}\rightarrow 0$, where $D_{\mathscr{H}}^{2}:=\inf_{q\in\mathscr{H}}\left\|f_{\rho}-S_{K}^{*}q\right\|_{\rho}^{2}$ is the approximation error of $f_{\rho}$ by a function in $\mathscr{H}$); and third, to derive an exact computational-statistical efficiency trade-off as a function of the $(l,N,\lambda)$ triplet, and of the difficulty of the problem.

Assumptions

In this section, we detail our assumptions on the $(\mathscr{X},Y,k,K)$ quartet. Our analysis for the well-specified case uses existing ridge regression results (Caponnetto and De Vito, 2007) focusing on problem (8) where only a single-stage sampling is present, hence we have to verify the associated conditions. Though we make use of these results, the analysis still remains challenging; the available bounds can moderately shorten our proof. We must take particular care in verifying that Caponnetto and De Vito (2007)’s conditions are met, since they need to hold for the space of mean embeddings of the distributions ($X=\mu\left(\mathscr{M}_{1}^{+}(\mathscr{X})\right)$), whose properties as a function of $\mathscr{X}$ and $H$ must themselves be established.

$(\mathscr{X},\tau)$ is a separable, topological space.

$k$ is bounded, in other words $\exists B_{k}<\infty$ such that $\sup_{u\in\mathscr{X}}k(u,u)\leq B_{k}$, and continuous.

The $\{K_{\mu_{a}}\}_{\mu_{a}\in X}$ operator family is uniformly bounded in Hilbert-Schmidt norm and Hölder continuous in operator norm. Formally, $\exists B_{K}<\infty$ such that

and $\exists L>0$, $h\in(0,1]$ such that the mapping $K_{(\cdot)}:X\rightarrow\mathscr{L}(Y,\mathscr{H})$ is Hölder continuous:

$y$ is bounded: $\exists C<\infty$ such that $\left\|y\right\|_{Y}\leq C$ almost surely.

These requirements hold under mild conditions: in Section 8, we provide insight into the consequences of our assumptions, with several concrete illustrations (e.g. regression with set- and RBF-type kernels).

Error Bounds, Consistency & Computational-Statistical Efficiency Trade-off

In this section, we present our analysis of the consistency of the mean embedding based ridge regression (MERR) method.

Given the estimator ($f^{\lambda}_{\hat{\mathbf{z}}}$) in Eq. (9), we derive finite-sample high probability upper bounds (see Theorems 2 and 7) for the excess risk $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)$, and in the misspecified setting, for the excess risk compared to the best attainable value from $\mathscr{H}$, i.e., $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)-D_{\mathscr{H}}^{2}$. We illustrate the bounds for particular classes of prior distributions, and work through special cases to obtain consistency conditions and computational-statistical efficiency trade-offs (see Theorems 4, 9 and the 3rd bullet of Remark 8). The main challenge is how to turn the convergence rates of the mean embeddings into those for an error $\mathcal{E}$ of the predictor. Although the main ideas of the proofs can be summarized relatively briefly, the full details are more demanding. High-level ideas with the sketches of the proofs and the obtained results are presented in Section 4.1 (well-specified case) and Section 4.2 (misspecified case). The derivations of some technical details of Theorems 2 and 7 are available in Section 7.

We first focus on the well-specified case ($f_{\rho}\in\mathscr{H}$) and present our first main result. We derive a high probability upper bound for the excess risk $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)$ of the MERR method (Theorem 2). The upper bound is instantiated for a general class of prior distributions (Theorem 4), which leads to a simple computational-statistical efficiency description (Theorem 5); this shows (among others) conditions when the MERR technique is able to achieve the one-stage sampled minimax optimal rate. We first give a high-level sketch of our convergence analysis and an intuitive interpretation of the results. An outline of the main proof ideas is given below, with technical details in Section 7.

Let us define $\mathbf{x}=\{x_{i}\}_{i=1}^{l}$ and $\hat{\mathbf{x}}=\{\{x_{i,n}\}_{n=1}^{N}\}_{i=1}^{l}$ as the ‘x-part’ of $\mathbf{z}$ and $\hat{\mathbf{z}}$, respectively. One can express $f_{\mathbf{z}}^{\lambda}$ [Eq. (8)] (Caponnetto and De Vito, 2007), and similarly $f_{\hat{\mathbf{z}}}^{\lambda}$ [Eq. (9)], as

where $T_{\mu_{a}}=K_{\mu_{a}}K_{\mu_{a}}^{*}\in\mathscr{L}(\mathscr{H})$ ($\mu_{a}\in X$), $T_{\mathbf{x}},T_{\hat{\mathbf{x}}}:\mathscr{H}\rightarrow\mathscr{H}$, $g_{\mathbf{z}},g_{\hat{\mathbf{z}}}\in\mathscr{H}$. By these explicit expressions, one can decompose the excess risk into 5 terms (Szabó et al., 2015, Section A.1.8):

Three of the terms ($S_{1}$, $S_{2}$, $\mathscr{A}(\lambda)$) are identical to the terms in Caponnetto and De Vito (2007), hence the earlier bounds can be applied. The two new terms ($S_{-1}$, $S_{0}$) resulting from two-stage sampling will be upper bounded by making use of the convergence of the empirical mean embeddings. These bounds will lead to the following results:

if $l\geq 2C_{\eta}B_{K}\mathscr{N}(\lambda)/\lambda$, $\lambda\leq\left\|T\right\|_{\mathscr{L}(\mathscr{H})}$ and N\geq\big{(}1+\sqrt{\log(l)+\delta}\big{)}^{2}2^{\frac{h+6}{h}}B_{k}(B_{K})^{\frac{1}{h}}L^{\frac{2}{h}}/\lambda^{\frac{2}{h}}.

Below we specialize our excess risk bound for a general prior class, which captures the difficulty of the regression problem as defined in Caponnetto and De Vito (2007). This $\mathcal{P}(b,c)$ class is described by two parameters $b$ and $c$: larger $b$ means faster decay of the eigenvalues of the covariance operator $T$ [in Eq. (17)], hence smaller effective input dimension; larger $c$ corresponds to a smoother regression function. Formally:

Definition of the $\mathcal{P}(b,c)$ class: Let us fix the positive constants $R$, $\alpha$, $\beta$. Then given $1<b$, $c\in(1,2]$, the $\mathcal{P}(b,c)$ class is the set of probability distributions $\rho$ on $Z=X\times Y$ such that

a range space assumption is satisfied: $\exists g\in\mathscr{H}$ s.t. $f_{\rho}=T^{\frac{c-1}{2}}g$ with $\left\|g\right\|_{\mathscr{H}}^{2}\leq R$,

in the spectral decomposition of $T=\sum_{n=1}^{\infty}\lambda_{n}\left<\cdot,e_{n}\right>_{\mathscr{H}}e_{n}$, where $(e_{n})_{n=1}^{\infty}$ is a basis of $Ker(T)^{\perp}$, the eigenvalues of $T$ satisfy $\alpha\leq n^{b}\lambda_{n}\leq\beta\quad(\forall n\geq 1)$.

We make few remarks about the $\mathcal{P}(b,c)$ class:

Range space assumption on $f_{\rho}$: The smoothness of $f_{\rho}$ is expressed as a range space assumption, which is slightly different from the standard smoothness conditions appearing in non-parametric function estimation. By the spectral decomposition of $T$ given above [$\lambda_{1}\geq\lambda_{2}\geq\ldots>0,\lim_{n\rightarrow\infty}\lambda_{n}=0$], $T^{r}u=\sum_{n=1}^{\infty}(\lambda_{n})^{r}\left<u,e_{n}\right>_{\mathscr{H}}e_{n}\quad(r=\frac{c-1}{2}\geq 0,u\in\mathscr{H})$ and

Specifically, in the limit as $r\rightarrow 0$, we obtain $f_{\rho}\in Im(T^{0})=Im(I)=\mathscr{H}$ (no constraint); larger values of $r$ give rise to faster decay of the $(c_{n})_{n=1}^{\infty}$ Fourier coefficients. This is the concrete meaning of $f_{\rho}\in Im(T^{r})$.

Spectral decay condition: We can provide a simple illustration of when the spectral decay conditions hold, in the event that the distributions are normal with means $m_{i}$ and identical variance ($x_{i}=N(m_{i},\sigma^{2}I$)). When Gaussian kernels ($k$) are used with linear $K$, then $K(\mu_{x_{i}},\mu_{x_{j}})=e^{-c\left\|m_{i}-m_{j}\right\|^{2}}$ (Muandet et al., 2012, Table 1, line 2) (Gaussian, with arguments equal to the difference in means). Thus, this Gram matrix will correspond to the Gram matrix using a Gaussian kernel between points $m_{i}$. The spectral decay of the Gram matrix will correspond to that of the Gaussian kernel, with points drawn from the meta-distribution over the $m_{i}$. Thus, the source conditions are analysed in the same manner as for Gaussian Gram matrices: see e.g. Steinwart and Christmann (2008) for a discussion of these spectral decay properties.

In the $\mathcal{P}(b,c)$ family, the behaviour of $\mathscr{A}(\lambda)$, $\mathscr{B}(\lambda)$ and $\mathscr{N}(\lambda)$ is known: $\mathscr{A}(\lambda)\leq R\lambda^{c}$, $\mathscr{B}(\lambda)\leq R\lambda^{c-1}$, $\mathscr{N}(\lambda)\leq\beta\frac{b}{b-1}\lambda^{-\frac{1}{b}}$. Specializing Theorem 2 and retaining its assumptions, we get:

Suppose the conditions in Theorem 2 hold. Let $\rho\in\mathcal{P}(b,c)$, where $1<b$ and $c\in(1,2]$. Then

Discarding the constants in Theorem 4, the study of convergence of the excess risk $\mathcal{E}(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho})$ to boils down to finding $N$ and $\lambda$ (as a function of $l$) where $N\rightarrow\infty$, $\lambda\rightarrow 0$ and

as $l\rightarrow\infty$. Let us choose $N=l^{\frac{a}{h}}\log(l)$; in this case Eq. (19) reduces to

One can assume that $a>0$, otherwise $r(l,\lambda)\rightarrow 0$ fails to hold; in other words, $N$ should grow faster than $\log(l)$. Matching the ‘bias’ ($\lambda^{s}$) and ‘variance’ (other) terms in $r(l,\lambda)$ to choose $\lambda$, and guaranteeing that the matched terms dominate and the constraints in Eq. (20) hold, one gets the following simple description for the computational-statistical efficiency trade-off:The derivations are available in the supplement.

(Computational-statistical efficiency trade-off; well-specified case; $\rho\in\mathcal{P}(b,c)$) Suppose the conditions in Theorem 2 hold. Let $\rho\in\mathcal{P}(b,c)$ and $N=l^{\frac{a}{h}}\log(l)$, where $0<a$, $1<b$, $c\in(1,2]$. If

$a\leq\frac{b(c+1)}{bc+1}$, then $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)=\mathcal{O}_{p}\left(l^{-\frac{ac}{c+1}}\right)$ with $\lambda=l^{-\frac{a}{c+1}}$,

$a\geq\frac{b(c+1)}{bc+1}$ then $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)=\mathcal{O}_{p}\left(l^{-\frac{bc}{bc+1}}\right)$ with $\lambda=l^{-\frac{b}{bc+1}}$.

Theorem 5 formulates an exact computational-statistical efficiency trade-off for the choice of the bag size ($N$) as a function of the number of distributions ($l$) and problem difficulty ($b$, $c$).

$a$-dependence: A smaller bag size (smaller $a$; $N=l^{\frac{a}{h}}\log(l)$) means computational savings, but reduced statistical efficiency. It is not worth increasing $a$ above $\frac{b(c+1)}{bc+1}$ since from that point the rate becomes $r(l)=l^{-\frac{bc}{bc+1}}$; remarkably, this rate is minimax in the one-stage sampled setup (Caponnetto and De Vito, 2007). The sensible choice $a=\frac{b(c+1)}{bc+1}<2$ means that the one-stage sampled minimax rate can be achieved in the two-stage sampled setting with bag size $N$ sub-quadratic in $l$.

$h$-dependence: In accord with our ‘smoothness’ assumptions it is rewarding to use smoother $K$ kernels (larger $h\in(0,1]$) since this reduces the bag size [$N=l^{\frac{a}{h}}\log(l)$].

$c$-dependence: The strictly decreasing property of $c\mapsto\frac{b(c+1)}{bc+1}$ implies that for ‘smoother’ problems (larger $c$) fewer samples ($N$) are sufficient.

Below we elaborate on the sketched high-level idea and prove Theorem 2. Proof of Theorem 2 (detailed derivations of each step can be found in Section 7.1)

Decomposition of the excess risk: We have the following upper bound for the excess risk

using which upper bounds on $S_{1}$ and $S_{2}$ that hold with probability $1-\eta$ are obtained. It is known that $\mathscr{A}(\lambda)\leq R\lambda^{c}$.

Probabilistic bounds on $\|g_{\hat{\mathbf{z}}}-g_{\mathbf{z}}\|_{\mathscr{H}}^{2}$, $\|T_{\mathbf{x}}-T_{\hat{\mathbf{x}}}\|_{\mathscr{L}(\mathscr{H})}^{2}$, $\|\sqrt{T}(T_{\hat{\mathbf{x}}}+\lambda I)^{-1}\|_{\mathscr{L}(\mathscr{H})}^{2}$, $\|f_{\mathbf{z}}^{\lambda}\|_{\mathscr{H}}^{2}$: One can bound $S_{-1}$ and $S_{0}$ as

For the terms on the r.h.s., we derive upper bounds [for the definition of $\alpha$, see Eq. (24)]

The bounds hold under the following conditions:

$\|g_{\hat{\mathbf{z}}}-g_{\mathbf{z}}\|_{\mathscr{H}}^{2}$ (see Section 7.1.1): if the empirical mean embeddings are close to their population counterparts, i.e.,

This event has probability $1-le^{-\alpha}$ over all $i=1,\ldots,l$ samples; see (Altun and Smola, 2006) and (Szabó et al., 2015, Section A.1.10).

$\|T_{\mathbf{x}}-T_{\hat{\mathbf{x}}}\|_{\mathscr{L}(\mathscr{H})}^{2}$ (see Section 7.1.2): (24) is assumed.

$\|\sqrt{T}(T_{\hat{\mathbf{x}}}+\lambda I)^{-1}\|_{\mathscr{L}(\mathscr{H})}^{2}$ (Szabó et al., 2015, Section A.1.11): (24), $\bm{\Theta}(\lambda,\mathbf{z})\leq\frac{1}{2}$, and

$\|f_{\mathbf{z}}^{\lambda}\|_{\mathscr{H}}^{2}$: The bound is guaranteed to hold under the conditions of the bounds of $S_{1}$ and $S_{2}$.8

Union bound: By applying an $\alpha=\log(l)+\delta$ reparameterization, and combining the received upper bounds with Caponnetto and De Vito (2007)’s results for $S_{1}$ and $S_{2}$, Theorem 2 follows (Section 7.1.3) with a union bound.

Finally, we note that existing results/ideas were used at two points to simplify our analysis: bounding $S_{1}$, $S_{2}$, $\Theta(\lambda,\mathbf{z})$, $\left\|f_{\mathbf{z}}^{\lambda}\right\|_{\mathscr{H}}^{2}$ (Caponnetto and De Vito, 2007) and $\left\|\mu_{x_{i}}-\mu_{\hat{x}_{i}}\right\|_{H}$ (Altun and Smola, 2006). We also corrected some constants in the previous works (Altun and Smola, 2006; Caponnetto and De Vito, 2007).

2 Results for the Misspecified Case

In this section, we focus on the misspecified case ($f_{\rho}\in L^{2}_{\rho_{X}}\backslash\mathscr{H}$) and present our second main result, which was inspired by the proof technique of Sriperumbudur et al. (2014, Theorem 12). We derive a high probability upper bound for $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)$, i.e., the excess risk of the MERR method (Theorem 7) which gives rise to consistency results (3rd bullet of Remark 8) and precise computational-statistical efficiency trade-off (Theorem 9). Theorem 7 consists of two finite-sample bounds:

The first, more general bound [Eq. (27)] will be used to show consistency in the misspecified case (see the 3rd bullet of Remark 8), in other words that $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)$ can be driven to its smallest possible value determined by the “richness” of $\mathscr{H}$:

The value of $D_{\mathscr{H}}$ equals the approximation error of $f_{\rho}$ by a function from $\mathscr{H}$. Specifically, if $\mathscr{H}$ [precisely $S_{K}^{*}(\mathscr{H})=\{S_{K}^{*}q:q\in\mathscr{H}\}\subseteq L^{2}_{\rho_{X}}$] is dense in $L^{2}_{\rho_{X}}$, then $D_{\mathscr{H}}=0$.

The second, specialized result [Eq. (28)] under additional smoothness assumptions on $f_{\rho}$ will give rise to a precise computational-statistical efficiency trade-off in terms of the problem difficulty ($s$) and sample numbers ($l$, $N$); this result can be seen as the misspecified analogue of Theorem 5.

After stating our results, the main ideas of the proof follow; further technical details are available in Section 7.2. Our main theorem for bounding the excess risk is as follows:

Then for arbitrary $q\in\mathscr{H}$ with probability at least $1-\eta-e^{-\delta}$

where $D_{a}(\lambda,q)=\|f_{\rho}-S_{K}^{*}q\|_{\rho}+\max(1,\|T\|_{\mathscr{L}(\mathscr{H})})\lambda^{\frac{1}{2}}\|q\|_{\mathscr{H}}$.

We give a short insight into the assumptions of Theorem 7, followed by consequences of the theorem.

$\sqrt{\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)}$: Notice that in the bounds [(27), (28)], instead of the excess risk, its square root appears; this has technical reasons, as it is easier to have the $D_{a}(\lambda,q)$ quantity (without multiplicative constants) appear on the r.h.s. of Eq. (27) with this form.

Consistency in the misspecified case: The consequence of Theorem 7(1) is as follows. Discarding the constants in Eq. (27), we obtain the upper bound (notice that the constant multiplier of $\left\|f_{\rho}-S_{K}^{*}q\right\|_{\rho}$ in the last term was one):

By choosing $N=l^{1/h}\log l$, $\sqrt{r(l,\lambda)}$ is bounded by

Discarding the constants in Eq. (28) we get10

Our goal is to drive $r(l,N,\lambda)$ to zero with a suitable choice of the $(l,N,\lambda)$ triplet under the stronger range space assumption. Since in Eq. (29) $\min(1,s)$ appears, one can assume without loss of generality that $s\in(0,1]$; consequently $1-\frac{s}{2}\in\left[\frac{1}{2},1\right)$ and $\frac{1}{l^{\frac{1}{2}}\lambda^{\frac{1}{2}}}\leq\frac{1}{\lambda^{1-\frac{s}{2}}l^{\frac{1}{2}}}$. Let us choose $N=l^{2a/h}\log(l)$; in this case using the previous dominance note, Eq. (29) reduces to the study of

One can assume that $a>0$, otherwise $r(l,\lambda)\rightarrow 0$ fails to hold: in other words, $N$ should grow faster than $\log(l)$. Matching the ‘bias’ ($\lambda^{s}$) and ‘variance’ (other) terms in $r(l,\lambda)$ to choose $\lambda$, guaranteeing that the matched terms dominate and the constraint in Eq. (30) hold, one can arrive at the following computational-statistical efficiency trade-off:8

$a\leq\frac{s+1}{s+2}$, then $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)=\mathcal{O}_{p}\left(l^{-\frac{2sa}{s+1}}\right)$ with $\lambda=l^{-\frac{a}{s+1}}$,

$a\geq\frac{s+1}{s+2}$, then $\mathcal{E}\left(f^{\lambda}_{\hat{\mathbf{z}}},f_{\rho}\right)=\mathcal{O}_{p}\left(l^{-\frac{2s}{s+2}}\right)$ with $\lambda=l^{-\frac{1}{s+2}}$.

Theorem 9 provides a complete computational-statistical efficiency trade-off description for the choice of the bag size $(N)$ as a number of the distributions ($l$).

$a$-dependence: A smaller value of ‘$a$’ (smaller bags $N=l^{2a/h}\log(l)$) leads to a computational advantage, but one looses in statistical efficiency. As ‘$a$’ reaches $\frac{s+1}{s+2}$, the rate becomes $r(l)=l^{-\frac{2s}{s+2}}$ and one does not gain from further increasing the value of $a$. The sensible choice of $a=\frac{s+1}{s+2}\leq\frac{2}{3}$ means that $N$ can again be sub-quadratic ($2a<\frac{4}{3}<2$) in $l$.

$h$-dependence: By using smoother $K$ kernels (larger $h\in(0,1]$) one can reduce the size of the bags: $h\mapsto 2a/h$ is decreasing in $h$. This is compatible with our smoothness requirement on $f_{\rho}$.

$s$-dependence: “Easier” tasks (larger $s$) give rise to faster convergence. Indeed, in the $r(l)=l^{-\frac{2s}{s+2}}$ rate the $s\mapsto\frac{2s}{s+2}$ exponent is strictly increasing function of the problem difficulty ($s$). For example, for extremely non-smooth regression problems ($s\approx 0$) the convergence can be arbitrary slow ($\lim_{s\rightarrow 0}\frac{2s}{s+2}=0$). In the smooth case ($s=1$) $\lim_{s\rightarrow 1}\frac{2s}{s+2}=\frac{2}{3}$ and one can achieve the $r(l)=l^{-\frac{2}{3}}$ rate.

The main steps of the proof of Theorem 7 are as follows: Proof of Theorem 7 (the details of the derivation are available in Section 7.2) Steps 1-7 will be identical in both proofs,Importantly, with a slight modification of the more general, first part of Theorem 7, one can get the specialized second setting of the theorem (see Step 8). and we present them jointly.

Decomposition of the excess risk: By the triangle inequality, we have

Bound on $\left\|S_{K}^{*}\left(f^{\lambda}_{\hat{\mathbf{z}}}-f^{\lambda}_{\mathbf{z}}\right)\right\|_{\rho}$: UsingSee for example de Vito et al. (2006) on page 88 with the $(\mathscr{H},\mathscr{G},A,T):=(\mathscr{H},L^{2}_{\rho_{X}},S_{K}^{*},T)$ choice. the fact that

and the definitions of $S_{-1}$ and $S_{0}$ [see Eqs. (15)-(16)], we obtain

through an application of triangle inequality. One can derive without a $\mathcal{P}(b,c)$ prior assumption (Section 7.2.1) the upper boundSee the remark at the end of Section 7.2.1.

for the r.h.s. of Eq. (33) under the conditions that $\bm{\Theta}(\lambda,\mathbf{z})\leq\frac{1}{2}$ (which holds with probability $1-\eta$ if $\left[12B_{K}\log(2/\eta)/\lambda\right]^{2}\leq l$), and that Eqs. (24)-(25) hold.

Decomposition of \big{\|}S_{K}^{*}f^{\lambda}_{\mathbf{z}}-f_{\rho}\big{\|}_{\rho}: By the triangle inequality and Eq. (32), we have

Decomposition of \big{\|}\sqrt{T}\left(f^{\lambda}_{\mathbf{z}}-f^{\lambda}\right)\big{\|}_{\mathscr{H}}: Making use of the analytical expressions for $f^{\lambda}_{\mathbf{z}}$ and $f^{\lambda}$ [see Eq. (13) and Eq. (17)], and the operator Woodbury formula (Ding and Zhou, 2008, Theorem 2.1, page 724) we arrive at the decomposition (see Section 7.2.2)

where $g_{\rho}=S_{K}f_{\rho}$. As it is known (Caponnetto and De Vito, 2007, page 348) $\|\sqrt{T}(T_{\mathbf{x}}+\lambda I)^{-1}\|_{\mathscr{L}(\mathscr{H})}\leq 1/\sqrt{\lambda}$ provided that $\bm{\Theta}(\lambda,\mathbf{z})\leq\frac{1}{2}$.

Bound on $\left\|g_{\mathbf{z}}-g_{\rho}\right\|_{\mathscr{H}}$, $\left\|T-T_{\mathbf{x}}\right\|_{\mathscr{L}(\mathscr{H})}$: By concentration arguments the bounds

hold with probability at least $1-\eta$, each (see Section 7.2.3, 7.2.4).

where we used at the last step that $\min(1,s+1)=1$; this follows from $s\geq 0$.

Bound on $\left\|S_{K}^{*}f^{\lambda}-f_{\rho}\right\|_{\rho}$:

No range space assumption: One can construct (Section 7.2.6) the bound

which holds for arbitrary $q\in\mathscr{H}$.

Union bound: Applying an $\alpha=\log(l)+\delta$ reparameterization, changing $\eta$ to $\frac{\eta}{3}$ and combining the derived results (in case of the first statement with $s=0$) with a union bound, Theorem 7 follows.

To contrast the derivation of the well- and the misspecified cases, we note that previous results [Section 4.1, or Caponnetto and De Vito (2007)’s bound] were used at two points:

In Step 2 by using Eq. (32) and transforming the $L^{2}_{\rho_{X}}$ error $\left\|S_{K}^{*}\left(f^{\lambda}_{\hat{\mathbf{z}}}-f^{\lambda}_{\mathbf{z}}\right)\right\|_{\rho}$ to $\mathscr{H}$, we could rely on our previous bounds for $S_{-1}$ and $S_{0}$. However, we were required to use a different concentration argument to guarantee $\bm{\Theta}(\lambda,\mathbf{z})\leq\frac{1}{2}$ since we no longer assume the $\mathcal{P}(b,c)$ prior class.

In Step 4 the first term could be bounded by Caponnetto and De Vito (2007). Its $\bm{\Theta}(\lambda,\mathbf{z})\leq\frac{1}{2}$ condition was guaranteed by Step 2; and see Section 7.2.1.

We note that our misspecified proof method was inspired by Sriperumbudur et al. (2014, Theorem 12), where the authors focused on the consistency of an infinite-dimensional exponential family estimator.

Related Work

In this section we discuss existing approaches and heuristic techniques to tackle learning problems on distributions.

Methods based on parametric assumptions: A number of methods have been proposed to compute the similarity of distributions or bags of samples. As a first approach, one could fit a parametric model to the bags, and estimate the similarity of the bags based on the obtained parameters. It is then possible to define learning algorithms on the basis of these similarities, which often take analytical form. Typical examples with explicit formulas include Gaussians, finite mixtures of Gaussians, and distributions from the exponential family (with known log-normalizer function and zero carrier measure, see Kondor and Jebara, 2003; Jebara et al., 2004; Wang et al., 2009; Nielsen and Nock, 2012). A major limitation of these methods, however, is that they apply quite simple parametric assumptions, which may not be sufficient or verifiable in practise.

Methods based on parametric assumption in a RKHS: A heuristic related to the parametric approach is to assume that the training distributions are Gaussians in a reproducing kernel Hilbert space (see for example Jebara et al., 2004; Zhou and Chellappa, 2006, and references therein). This assumption is algorithmically appealing, as many divergence measures for Gaussians can be computed in closed form using only inner products, making them straightforward to kernelize. A fundamental shortfall of kernelized Gaussian divergences is the lack of their consistency analysis in specific learning algorithms.

Kernels based techniques: A more theoretically grounded approach to learning on distributions has been to define positive definite kernels on the basis of statistical divergence measures on distributions, or by metrics on non-negative numbers; these can then be used in kernel algorithms. This category includes work on semigroup kernels (Cuturi et al., 2005), non-extensive information theoretical kernel constructions (Martins et al., 2009), and kernels based on Hilbertian metrics (Hein and Bousquet, 2005). For example, the intuition of semigroup kernels (Cuturi et al., 2005) is as follows: if two measures or sets of points overlap, then their sum is expected to be more concentrated. The value of dispersion can be measured by entropy or inverse generalized variance. In the second type of approach (Hein and Bousquet, 2005), homogeneous Hilbert metrics on the non-negative real line are used to define the similarity of probability distributions. While these techniques guarantee to provide valid kernels on certain restricted domains of measures, the performance of learning algorithms based on finite-sample estimates of these kernels remains a challenging open question. One might also plug into learning algorithms (based on similarities of distributions) consistent Rényi and Tsallis divergence estimates (Póczos et al., 2011, 2012), but these similarity indices are not kernels, and their consistency in specific learning tasks remains an open question.

Multi-instance learning: An alternative paradigm in learning when the inputs are “bags of objects” is to simply treat each input as a finite set: this leads to the multi-instance learning task (MIL, see Dietterich et al., 1997; Ray and Page, 2001; Dooly et al., 2002). In MIL one is given a set of labelled bags, and the task of the learner is to find the mapping from the bags to the labels. Many important examples fit into the MIL framework: for example, different configurations of a given molecule can be handled as a bag of shapes, images can be considered as a set of patches or regions of interest, a video can be seen as a collection of images, a document might be described as a bag of words or paragraphs, a web page can be identified by its links, a group of people on a social network can be captured by their friendship graphs, in a biological experiment a subject can be identified by his/her time series trials, or a customer might be characterized by his/her shopping records. The MIL approach has been applied in several domains; see the reviews from Babenko (2004); Zhou (2004); Foulds and Frank (2010); Amores (2013).

“Bag-of-objects” methods (MIL, not classification): Beyond classification, there exist several heuristics—without consistency guarantees—for many other multi-instance problems in the literature, including regression (Ray and Page, 2001; Dooly et al., 2002; Zhou et al., 2009; Kwok and Cheung, 2007), clustering (Zhang and Zhou, 2009; Zhang et al., 2009, 2011; Chen and Wu, 2012), ranking (Bergeron et al., 2008; Hu et al., 2008; Bergeron et al., 2012), outlier detection (Wu et al., 2010), transfer learning (Raykar et al., 2008; Zhang and Si, 2009), and feature selection, -weighting and -extraction (also called dimensionality reduction, low-dimensional embedding, manifold learning, see Raykar et al., 2008; Ping et al., 2010; Sun et al., 2010; Carter et al., 2011; Zafra et al., 2013; Chai et al., 2014a, b, and references therein).

Approaches using set metrics: Adapting the bag viewpoint of MIL, one can come up with set metric based learning algorithms.Often these “metrics” are only semi-metrics, as they do not satisfy the triangle inequality. Probably one of the most well-known set metrics is the Hausdorff metric (Edgar, 1995), which is defined for non-empty compact sets of metric spaces, specifically for sets containing finitely many points. There also exist other (semi)metric constructions on points sets (Eiter and Mannila, 1997; Ramon and Bruynooghe, 2001). Unfortunately, the classical Hausdorff metric is highly sensitive to outliers, seriously limiting its practical applicability. In order to mitigate this deficiency, several variants of the Hausdorff metric have been designed in the MIL literature, such as the maximal-, the minimal- and the ranked Hausdorff metrics, with successful applications in MIC (Wang and Zucker, 2000) and multi-instance outlier detection (Wu et al., 2010); and the average Hausdorff metric (Zhang and Zhou, 2009) and contextual Hausdorff dissimilarity (Chen and Wu, 2012), which have been found useful in multi-instance clustering. Unfortunately, these methods lack any theoretical guarantee when applied in specific learning problems.

Functional data analysis techniques: Finally, the distribution regression task might also be interpreted as a functional data analysis problem (Ramsay and Silverman, 2002, 2005; Müller, 2005), by considering the probability measures $x_{i}$ as functions. This is a highly non-standard setup, however, since these functions ($x_{i}$) are defined on $\sigma$-algebras and are non-negative, $\sigma$-additive.

Conclusion

We have established a learning theory of distribution regression, where the inputs are probability measures on separable, topological domains endowed with reproducing kernels, and the outputs are elements of a separable Hilbert space. We studied a ridge regression scheme defined on embeddings of the input distributions to a reproducing kernel Hilbert space, which has a simple analytical solution, as well as theoretically sound, efficient methods for approximation (Zhang et al., 2015; Richtárik and Takác̆, 2016; Alaoui and Mahoney, 2015; Yang et al., 2016; Rudi et al., 2015). We derived explicit bounds on the excess risk as a function of the number of samples and problem difficulty. We tackled both the well-specified case (when the regression function belongs to the assumed RKHS modelling class), and the more general misspecified setup. As a special case of our results, we proved the consistency of regression for set kernels (Haussler, 1999; Gärtner et al., 2002), which was a $17$-year-old open problem, and for a recent kernel family (Christmann and Steinwart, 2010), which we have expanded upon (Table 1). We proved an exact computational-statistical efficiency trade-off for the MERR estimator: in the well-specified setting, we showed how to choose the bag size in the two-stage sampled setup to match the one-stage sampled minimax optimal rate (Caponnetto and De Vito, 2007); and in the misspecified setting, our rates approximate closely an asymptotically optimal estimator imposing stricter eigenvalue decay conditions (Steinwart et al., 2009). Several exciting open questions remain, including whether improved/optimal rates can be derived in the misspecified case, whether we can obtain consistency guarantees for non-point estimates, and how to handle non-ridge extensions.

Finally, we note that although the primary focus of the current paper was theoretical, we have applied the MERR method (Szabó et al., 2015, Section A.2) to supervised entropy learning and aerosol prediction based on multispectral satellite images.For code, see https://bitbucket.org/szzoli/ite/. In future work, we will address applications with vector-valued outputs.

Proofs

We provide proofs for our results detailed in Section 4: Section 7.1 (resp. Section 7.2) focuses on the well-specified case (resp. misspecified setting). The used lemmas are enlisted in Section 7.3.

We give proof details concerning the excess risk in the well-specified case (Theorem 2).

By (13), (14) we get g_{\hat{\mathbf{z}}}-g_{\mathbf{z}}=\frac{1}{l}\sum_{i=1}^{l}\big{(}K_{\mu_{\hat{x}_{i}}}-K_{\mu_{x_{i}}}\big{)}y_{i}; hence by applying the Hölder property of $K_{(\cdot)}$, the boundedness of $y_{i}$ ($\left\|y_{i}\right\|_{Y}\leq C$) and (24), we obtain

with probability at least $1-le^{-\alpha}$, based on a union bound.

Using the definition of $T_{\mathbf{x}}$ and $T_{\hat{\mathbf{x}}}$, and exploiting (with $\|\cdot\|_{\mathscr{L}(\mathscr{H})}$) that in a normed spaceEq. (35) holds since $\left\|\cdot\right\|^{2}$ is convex function, thus $\left\|\frac{1}{n}\sum_{i=1}^{n}f_{i}\right\|^{2}\leq\frac{1}{n}\sum_{i=1}^{n}\left\|f_{i}\right\|^{2}$. $(N,\left\|\cdot\right\|)$, $f_{i}\in N$, ($i=1,\ldots,n$)

To upper bound $\|T_{\mu_{x_{i}}}-T_{\mu_{\hat{x}_{i}}}\|_{\mathscr{L}(\mathscr{H})}^{2}$, let us see how $T_{\mu_{u}}=K_{\mu_{a}}K^{*}_{\mu_{a}}$ acts. The existence of an $E\geq 0$ constant satisfying $\left\|(T_{\mu_{u}}-T_{\mu_{v}})(f)\right\|_{\mathscr{H}}\leq E\left\|f\right\|_{\mathscr{H}}$ implies $\left\|T_{\mu_{u}}-T_{\mu_{v}}\right\|_{\mathscr{L}(\mathscr{H})}\leq E$. We continue with the l.h.s. of this equation using Eq. (35):

By Eq. (45) and the Hölder continuity of $K_{(\cdot)}$, one arrives at

Hence $\left\|(T_{\mu_{u}}-T_{\mu_{v}})(f)\right\|_{\mathscr{H}}^{2}\leq 4B_{K}L^{2}\left\|\mu_{u}-\mu_{v}\right\|_{H}^{2h}\left\|f\right\|_{\mathscr{H}}^{2}\hskip 2.84544pt\Rightarrow\hskip 2.84544ptE^{2}=4B_{K}L^{2}\left\|\mu_{u}-\mu_{v}\right\|_{H}^{2h}$. Exploiting this property in (36) with Eq. (24) we arrive to the bound

1.3 Proof: final union bound in Theorem 2

Until now, we obtained that if (i) the sample number $N$ satisfies Eq. (25), (ii) (24) holds (which has probability at least $1-le^{-\alpha}=1-e^{-[\alpha-\log(l)]}=1-e^{-\delta}$ applying a union bound argument; $\alpha=\log(l)+\delta$), and (iii) $\bm{\Theta}(\lambda,\mathbf{z})\leq\frac{1}{2}$ is fulfilled [see Eq. (22)], then

By taking into account Caponnetto and De Vito (2007)’s bounds for $S_{1}$ and $S_{2}$, $S_{1}\leq 32\log^{2}\left(\frac{6}{\eta}\right)\left[\frac{B_{K}M^{2}}{l^{2}\lambda}+\frac{\Sigma^{2}\mathscr{N}(\lambda)}{l}\right]$, $S_{2}\leq 8\log^{2}\left(\frac{6}{\eta}\right)\left[\frac{4B_{K}^{2}\mathscr{B}(\lambda)}{l^{2}\lambda}+\frac{B_{K}\mathscr{A}(\lambda)}{l\lambda}\right]$, plugging all the expressions to (21), we obtain Theorem 2 with a union bound.

2 Proofs of the Misspecified Case

We present the proof details concerning the excess risk in the misspecified case (Theorem 7).

where we made use of (2), the $\|T_{\mu_{\mathbf{x}_{i}}}\|_{\mathscr{L}_{2}(\mathscr{H})}\leq B_{K}$ identity following from the boundedness of $K$ (Caponnetto and De Vito, 2007, page 341, Eq. (13)), and the spectral theorem. Consequently, by the Bernstein’s inequality (Lemma 7.3.1 with $\mathscr{K}=\mathscr{L}_{2}(\mathscr{H})$, $B=2B_{K}/\lambda$, $\sigma=B_{K}/\lambda$) we obtain that for $\forall\eta\in(0,1)$

Thus, for $\bm{\Theta}(\lambda,\mathbf{z})\leq\frac{1}{2}$ with probability $1-\eta$ it is sufficient to have

Under these conditions, we arrived at the upper bound

where as opposed to Section 7.1.3 and Eq. (23) we used a slightly cruder $\left\|f_{\mathbf{z}}^{\lambda}\right\|_{\mathscr{H}}^{2}\leq\frac{C^{2}}{\lambda}$ bound; it holds without the $\mathcal{P}(b,c)$ assumption by the definition of $f_{\mathbf{z}}^{\lambda}$ and the boundedness of $y$ since $\lambda\left\|f_{\mathbf{z}}^{\lambda}\right\|_{\mathscr{H}}^{2}\leq\frac{1}{l}\sum_{i=1}^{l}\left\|y_{i}\right\|_{Y}^{2}\leq C^{2}$.

Remark: Notice that the price we pay for not assuming that the prior belongs to the $\mathcal{P}(b,c)$ class ($b>1$) is a slightly tighter $\frac{1}{\lambda^{2}}\leq l$ constraint [Eq. (38)] instead of $\frac{1}{\lambda^{1+\frac{1}{b}}}\leq l$ in Eq. (19), and a somewhat looser $\left\|f_{\mathbf{z}}^{\lambda}\right\|_{\mathscr{H}}^{2}$ bound.

Using the analytical formula of $f^{\lambda}_{\mathbf{z}}$ [see Eq. (13)] and that of $f^{\lambda}$ [see Eq.(17)]

one gets $(T+\lambda I)f^{\lambda}=S_{K}f_{\rho}$ $\Rightarrow$ $\lambda f^{\lambda}=S_{K}f_{\rho}-Tf^{\lambda}$ and

Let us rewrite $(T+\lambda I)^{-1}$ by the $(A+UV)^{-1}=A^{-1}-A^{-1}U\left(I+VA^{-1}U\right)^{-1}VA^{-1}$ operator Woodbury formula (Ding and Zhou, 2008, Theorem 2.1, page 724)

where the CBS (Cauchy-Bunyakovsky-Schwarz) inequality was applied. Since

using Eq. (41) and the analytical expression for $f^{\lambda}$ [see Eq. (39)] we have

Below we give upper bounds on these two terms.

exploiting the Parseval’s identity and that \big{(}\frac{\lambda_{i}}{\lambda_{i}+\lambda}-1\big{)}^{2}\leq 1.

Making use of the two derived bounds, we get \left\|S_{K}^{*}f^{\lambda}-f_{\rho}\right\|_{\rho}\leq\left\|f_{\rho}-S_{K}^{*}q\right\|_{\rho}+\max\big{(}1,\left\|T\right\|_{\mathscr{L}(\mathscr{H})}\big{)}\lambda^{\frac{1}{2}}\left\|q\right\|_{\mathscr{H}}.

3 Supplementary Lemmas

In this section, we list two lemmas used in the proofs.

3.2 Lemma on bounded, self-adjoint compact operators; Sriperumbudur et al. (2014, Proposition A.2, page 39)

Let $M$ be a bounded, self-adjoint compact operator on a separable Hilbert space $\mathscr{K}$. Let $a\geq 0$, $\lambda>0$, and $s\geq 0$. Let $f\in\mathscr{K}$ such that $f\in Im\left(M^{s}\right)$. If $s+a>0$, then

Note: specifically for $s=0$ we have $Im\left(M^{s}\right)=Im\left(I\right)=\mathscr{K}$, in other words, there is no additional range space constraint.

Discussion of Our Assumptions

We give a short insight into the consequences of our assumptions (detailed in Section 3) and present some concrete examples.

Well-definedness of $\rho$: The boundedness and continuity of $k$ imply the measurability of $\mu:(\mathscr{M}^{+}_{1}(\mathscr{X}),\mathcal{B}(\tau_{w}))\rightarrow(H,\mathcal{B}(H))$. Let $\tau$ denote the open sets on $H=H(k)$, $\left.\tau\right|_{X}=\{A\cap X:A\in\tau\}$ the subspace topology on $X$, and $\left.\mathcal{B}(H)\right|_{X}=\{A\cap X:A\in\mathcal{B}(H)\}$ the subspace $\sigma$-algebra on $X$. By noting (Schwartz, 1998, Corollary 5.2.13) that $\mathcal{B}\left(\left.\tau\right|_{X}\right)=\left.\mathcal{B}(H)\right|_{X}=\{A\in\mathcal{B}(H):A\subseteq X\}\subseteq\mathcal{B}(H)$, the H-measurability of $\mu$ guarantees the measurability of $\mu:(\mathscr{M}^{+}_{1}(\mathscr{X}),\mathcal{B}(\tau_{w}))\rightarrow(X,\left.\mathcal{B}(H)\right|_{X})$, and hence the well-definedness of $\rho$, the measure induced by $\mathscr{M}$ on $X\times Y$; for further details see (Szabó et al., 2015, Section A.1.1).Note that the referred proof also holds for separable Hilbert $Y$, and by the simplified reasoning above the original $X\in\mathcal{B}(H)$ condition could be avoided.