Journal of Machine Learning Research 17 (2016) 1-40

Submitted 12/14; Revised 1/16; Published 9/16

Learning Theory for Distribution Regression

arXiv:1411.2066v4 [] 21 Oct 2016

Zoltán Szabó∗

zoltan.szabo@gatsby.ucl.ac.uk

ORCID 0000-0001-6183-7603
Gatsby Unit, University College London
Sainsbury Wellcome Centre, 25 Howland Street
London - W1T 4JG, UK

Bharath K. Sriperumbudur

bks18@psu.edu

Department of Statistics
Pennsylvania State University
University Park, PA 16802, USA

Barnabás Póczos

bapoczos@cs.cmu.edu

Machine Learning Department
School of Computer Science
Carnegie Mellon University
5000 Forbes Avenue Pittsburgh PA 15213 USA

Arthur Gretton

arthur.gretton@gmail.com

ORCID 0000-0003-3169-7624
Gatsby Unit, University College London
Sainsbury Wellcome Centre, 25 Howland Street
London - W1T 4JG, UK

Editor: Ingo Steinwart

Abstract
We focus on the distribution regression problem: regressing to vector-valued outputs from
probability measures. Many important machine learning and statistical tasks fit into this
framework, including multi-instance learning and point estimation problems without analytical solution (such as hyperparameter or entropy estimation). Despite the large number
of available heuristics in the literature, the inherent two-stage sampled nature of the problem makes the theoretical analysis quite challenging, since in practice only samples from
sampled distributions are observable, and the estimates have to rely on similarities computed between sets of points. To the best of our knowledge, the only existing technique
with consistency guarantees for distribution regression requires kernel density estimation
as an intermediate step (which often performs poorly in practice), and the domain of the
distributions to be compact Euclidean. In this paper, we study a simple, analytically computable, ridge regression-based alternative to distribution regression, where we embed the
distributions to a reproducing kernel Hilbert space, and learn the regressor from the embeddings to the outputs. Our main contribution is to prove that this scheme is consistent
in the two-stage sampled setup under mild conditions (on separable topological domains
enriched with kernels): we present an exact computational-statistical efficiency trade-off
analysis showing that our estimator is able to match the one-stage sampled minimax op∗. Now at Department of Applied Mathematics, CMAP, École Polytechnique, Route de Saclay, 91128
Palaiseau, France.
c 2016 Zoltán Szabó, Bharath K. Sriperumbudur, Barnabás Póczos, and Arthur Gretton.

Szabó et al.

timal rate (Caponnetto and De Vito, 2007; Steinwart et al., 2009). This result answers a
17-year-old open question, establishing the consistency of the classical set kernel (Haussler,
1999; Gärtner et al., 2002) in regression. We also cover consistency for more recent kernels
on distributions, including those due to Christmann and Steinwart (2010).
Keywords: Two-Stage Sampled Distribution Regression, Kernel Ridge Regression, Mean
Embedding, Multi-Instance Learning, Minimax Optimality

1. 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 (Pó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).
Intuitive description of our goal: Let us start with a somewhat informal definition
of the distribution regression problem and an intuitive phrasing of our goals. Suppose that
our data consist of z = {(xi , yi )}li=1 , where xi is a probability distribution, yi is its label
(in the simplest case yi ∈ R or yi ∈ Rd ) and each (xi , yi ) pair is i.i.d. sampled from a
meta distribution M. However, we do not observe xi directly; rather, we observe a sample
i.i.d.
l
i
xi,1 , . . . , xi,Ni ∼ xi . Thus the observed data are ẑ = {({xi,n }N
n=1 , yi )}i=1 . Since xi,j is
sampled from xi , and xi is sampled from M, we call this process two-stage sampling. Our
goal is to predict a new yl+1 from a new batch of samples xl+1,1 , . . . , xl+1,Nl+1 drawn from a
new distribution xl+1 ∼ M. For example, in a medical application, the ith patient might be
identified with a probability distribution (xi ), which can be periodically assessed by blood
i
tests ({xi,n }N
n=1 ). We are also given some health indicator of the patient (yi ), which might
be inferred from his/her blood measurements. Based on the observations (ẑ), we might
wish to learn the mapping from the set of blood tests to the health indicator, and the hope
is that by observing more patients (larger l) and performing a larger number of tests (larger
Ni ) the estimated mapping (fˆ = fˆ(ẑ)) becomes more “precise”. Briefly, we consider the
following questions:
Can the distribution regression problem be solved consistently
under mild conditions? What is the exact computationalstatistical efficiency trade-off implied by the two-stage sampling?
In our work the estimated mapping (fˆ) is the analytical solution of a kernel ridge regression
(KRR) problem.1 The performance of fˆ depends on the assumed function class (H), the
1. Beyond its simple analytical formula, kernel ridge regression also allows efficient distributed (Zhang et al.,
2015; Richtárik and Takác̆, 2016), sketch (Alaoui and Mahoney, 2015; Yang et al., 2016) and Nyström
based approximations (Rudi et al., 2015).

2

Learning Theory for Distribution Regression

family of fˆ candidates used in the ridge formulation. We shall focus on the analysis of two
settings:
1. Well-specified case (f∗ ∈ H): In this case we assume that the regression function
f∗ belongs to H. We focus on bounding the goodness of fˆ compared to f∗ . In other
words, if R[f∗ ] denotes the prediction error (expected risk) of f∗ , then our goal is to
derive a finite-sample bound for the excess risk, E(fˆ, f∗ ) = R[fˆ] − 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 computationalstatistical efficiency trade-off of the estimator as a function of the sample number (l,
N = Ni , ∀ i) and the problem difficulty (see Theorem 5 and its corresponding remarks
for more details).
2. Misspecified case (f∗ ∈ L2 \H): Since in practise it might be hard to check whether
f∗ ∈ H, we also study the misspecified setting of f∗ ∈ L2 ; the relevant case is when
f∗ ∈ L2 \H. In the misspecified setting the ’richness’ of H has crucial importance,
2
2 = inf
in other words the size of DH
f ∈H kf∗ − f k2 , the approximation error from 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.,
2 can be driven to zero. Moreover, on smooth classes of f -s, we give a
E(fˆ, f∗ ) − DH
∗
simple and explicit description for the computational-statistical efficiency trade-off of
our estimator (see Theorem 9 and its corresponding remarks for more details).
There exist a vast number of heuristics to tackle learning problems on distributions; we
will review them in Section 5. However, to the best of our knowledge, the only prior work
addressing the consistency of regression on distributions requires kernel density estimation
(Póczos et al., 2013; Oliva et al., 2014; Sutherland et al., 2016), which assumes that the
response variable is scalar-valued,2 and the covariates are nonparametric continuous distributions on Rd . As in our setting, the exact forms of these distributions are unknown;
they are available only through finite sample sets. Póczos et al. estimated these distributions through a kernel density estimator (assuming these distributions have a density) and
then constructed a kernel regressor that acts on these kernel density estimates.3 Using the
classical bias-variance decomposition analysis for kernel regressors, they showed the consistency of the constructed kernel regressor, and provided a polynomial upper bound on the
rates, assuming the true regressor to be Hölder continuous, and the meta distribution that
generates the covariates xi to have finite doubling dimension (Kpotufe, 2011).4
One can define kernel learning algorithms on bags based on set kernels (Gä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
2. Oliva et al. (2013, 2015) consider the case where the responses are also distributions or functions.
3. We would like to clarify that the kernels used in their work are classical smoothing kernels—extensively
studied in non-parametric statistics (Györfi et al., 2002)—and not the reproducing kernels that appear
throughout our paper.
4. Using a random kitchen sinks approach, with orthonormal basis projection estimators Oliva et al. (2014);
Sutherland et al. (2016) propose distribution regression algorithms that can computationally handle large
scale datasets; as with Póczos et al. (2013), these approaches are based on density estimation in Rd .

3

Szabó et al.

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; Gä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); Gä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
1. well-specified case: We
(a) derive finite-sample bounds on the excess risk: We construct R[fˆ] − R[f∗ ] ≤

r(l, N, λ) bounds holding with high probability, where λ is the regularization parameter in the ridge problem (λ → 0, l → ∞, N = Ni → ∞).
(b) establish consistency and computational-statistical efficiency trade-off of the MERR
estimator on a general prior family 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 ∈ (1, 2]). In particular, when the number of samples per bag is chosen
bc
− bc+1
, which is
as N = la log(l) and a ≥ b(c+1)
bc+1 , then the learning rate saturates at l
known to be one-stage sampled minimax optimal (Caponnetto and De Vito, 2007).
In other words, by choosing a = 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 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
4

Learning Theory for Distribution Regression

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.
2. misspecified case: We establish consistency and convergence rates even if f∗ ∈
/ H.
Particularly, by deriving finite-sample bounds on the excess risk we
(a) prove that the MERR estimator can achieve the best possible approximation accu2 quantity can be driven to zero (recall that
racy from H, i.e. the R[fˆ] − R[f∗ ] − DH
DH = inf f ∈H kf∗ − f k2 ). Specifically, this result implies that if H is dense in L2
(DH = 0), then the excess risk R[fˆ] − R[f∗ ] converges to zero.

(b) analyse the computational-statistical efficiency trade-off: We show that by choosing
2sa
s+1
the bag size as N = l2a log(l) (a > 0) one can get rate l− s+1 for a ≤ s+2
, and the
2s

− s+2
rate saturates for a ≥ s+1
, where the difficulty of the regression problem
s+2 at l
is captured by s ∈ (0, 1] (a larger s means an easier problem). This means that
2
easier tasks give rise to faster convergence (for s = 1, the rate is l− 3 ), the bag size
4
N can again be sub-quadratic in l (2a ≤ 2(s+1)
s+2 ≤ 3 < 2), and the rate at saturation
2s

is close to r̃(l) = l− 2s+1 , which is the asymptotically optimal rate in the one-stage
sampled setup, with real-valued output and stricter eigenvalue decay conditions
(Steinwart et al., 2009).

Due to the differences in the assumptions made and the loss function used, a direct comparison of our theoretical result and that of Pó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. Pó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 Pó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 (H) in
the embedding L2 space. The fundamental challenge is whether it is possible to achieve the
best possible performance dictated by H; or in the special case when further smoothness
conditions hold on f∗ , what convergence rates can yet be attained, and what computationalstatistical efficiency trade-off realized. The second smoothness property could be modelled
5

Szabó et al.

for example by range spaces of (fractional) powers of integral operators associated to 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 H can be fully exploited, and convergence to the best approximation
available from H can be realized.5
There are numerous areas in machine learning and statistics, where estimating vectorvalued 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
(Á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 (Szabó 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,
5. 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.

6

Learning Theory for Distribution Regression

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.

2. The Distribution Regression Problem
Below we first introduce our notation (Section 2.1), then formally define the distribution
regression task (Section 2.2).
2.1 Notation
We use the following notations throughout the paper:
 
• Sets, topology, measure theory: Let K be a Hilbert space; cl V is the closure of a set
V ⊆ K. ×i∈I Si is the direct product of sets Si . f ◦g is the composition of function f and g.
Let (X, τ ) be a topological space and let B(X) := B(τ ) be the Borel σ-algebra induced by
the topology τ . If (X, d) is a metric space, then B = B(d) is the Borel σ-algebra generated
by the open sets induced by metric d. M+
1 (X) denotes the set of Borel probability measures
on the (X, B(X)) measurable space. Given measurable spaces (U1 , S1 ) and (U2 , S2 ), the
S1 ⊗ S2 product σ-algebra (Steinwart and Christmann, 2008, page 480) on the product
space U1 × U2 is the σ-algebra generated by the cylinder sets U1 × S2 , S1 × U2 (S1 ∈ S1 ,
+
S2 ∈ S2 ). The weak topology (τw = τw (X, τ )) on M
R 1 (X) is defined as the weakest topology
+
such that the Lh : (M1 (X), τw ) → R, Lh (x) = X h(u)dx(u) mapping is continuous for
all h ∈ Cb (X) = {(X, τ ) → R bounded, continuous functions}.
• Functional analysis: Let (N1 , k·kN1 ) and (N2 , k·kN2 ) denote two normed spaces, then
L(N1 , N2 ) stands for the space of N1 → N2 bounded linear operators; if N1 = N2 , we will
use the L(N1 ) = L(N1 , N2 ) shorthand. For M ∈ L(N1 , N2 ) the operator norm is defined
as kM kL(N1 ,N2 ) = sup06=h∈N1 kM hkN2 / khkN1 , Im(M ) = {M n1 }n1 ∈N1 denotes the range
of M , Ker(M ) = {n1 ∈ N1 : M n1 = 0} is the null space of M . Let K be a Hilbert space.
The adjoint operator M ∗ ∈ L(K) of an operator M ∈ L(K) is the operator such that
hM a, biK = ha, M ∗ biK for all a and b in K. M ∈ L(K) isP
called positive if hM a, aiK ≥ 0
∗
(∀a ∈ K), self-adjoint if M = M , and trace class if
j∈J h|M |ej , ej iK < ∞ for an
1

(ej )j∈J ONB (orthonormal basis) of K (|M | := (M ∗ M ) 2 ), in which case T r(M ) :=
P
j∈J hM ej , ej iK < ∞; compact if cl [M a : a ∈ K, kakK ≤ 1] is a compact set. Let K1
and K2 be Hilbert
spaces. M ∈ L(K1 , K2 ) is called Hilbert-Schmidt if kM k2L2 (K1 ,K2 ) =
P
T r(M ∗ M ) =
j∈J hM ej , M ej iK2 < ∞ for some (ej )j∈J ONB of K1 . The space of
Hilbert-Schmidt operators is denoted by L2 (K1 , K2 ) = {M ∈ L(K1 , K2 ) : kM kL2 (K1 ,K2) <
∞}. We use the shorthand notation L2 (K) = L2 (K, K) if K := K1 = K2 ; L2 (K)
is separable if and only if K is separable (Steinwart and Christmann, 2008, page 506).
Trace class and Hilbert-Schmidt operators over a K Hilbert space are compact operators
(Steinwart and Christmann, 2008, page 505-506); moreover,
kAkL(K) ≤ kAkL2 (K) ,

∀A ∈ L2 (K),

kABkL2 (K) ≤ kAkL2 (K) kBkL(K) ,
7

∀A, B ∈ L2 (K).

(1)
(2)

Szabó et al.

I is the identity operator; Il ∈ Rl×l is the identity matrix.
• RKHS, mean embedding: Let H = H(k) be an RKHS (Steinwart and Christmann,
2008) with k : X × X → R as the reproducing kernel. Denote by
Z

+
+
k(·, u)dx(u) = Eu∼x [k(·, u)] ∈ H
X = µ M1 (X) = {µx : x ∈ M1 (X)} ⊆ H, µx =
X

the set of mean embeddings (Berlinet and Thomas-Agnan, 2004) of the distributions to
the space H.6 Let Y be a separable Hilbert space, where the inner product is denoted
by h·, ·iY ; the associated norm is k·kY . H = H(K) is the Y -valued RKHS (Pedrick,
1957; Micchelli and Pontil, 2005; Carmeli et al., 2006, 2010) of X → Y functions with
K : X × X → L(Y ) as the reproducing kernel (we will present some concrete examples
of K in Section 3; see Table 1); Kµx ∈ L(Y, H) is defined as
K(µx , µt )(y) = (Kµt y)(µx ),

(∀µx , µt ∈ X), or K(·, µt )(y) = Kµt y.

(3)

Further, f (µx ) = Kµ∗x f (∀µx ∈ X, f ∈ H).
• Regression function: Let ρ be the µ-induced probability measure on the Z = X × Y
product space, and let ρ(µx , y) = ρ(y|µx )ρX (µx ) be the factorization of ρ into conditional
and marginal distributions.7 The regression function of ρ with respect to the (µx , y) pair
is denoted by
Z
y dρ(y|µa ) (µa ∈ X)
(4)
fρ (µa ) =
Y

and for f ∈ L2ρX let kf kρ =

q

hf, f iρ := kf kL2ρ

=
X

hR

i1
2
2
kf
(µ
)k
dρ
(µ
)
. Let us
a Y
X a
X

assume that the operator-valued kernel K : X × X → L(Y ) is a Mercer kernel (that is
H = H(K) ⊆ C(X, Y ) = {X → Y continuous functions}), is bounded (∃BK < ∞ such
that kK(x, x)kL(Y ) ≤ BK ), and is a compact operator for all x ∈ X. These requirements
∗ :
will be guaranteed by our assumptions, see Section 7.2.6. In this case, the inclusion SK
2
2
H ֒→ LρX is bounded, and its adjoint SK : LρX → H is given by
(SK g)(µu ) =

Z

K(µu , µt )g(µt )dρX (µt ).

(5)

X

We further define T̃ as
∗
T̃ = SK
SK : L2ρX → L2ρX ;

(6)

in other words, the result of operation (5) belongs to H, which is embedded in L2ρX . T̃ is
a compact, positive, self-adjoint operator (Carmeli et al., 2010, Proposition 3), thus by
the spectral theorem T̃ s exists, where s ≥ 0.
6. The x 7→ µx mapping is defined for all x ∈ M+
1 (X) if k is bounded, i.e., supu∈X k(u, u) < ∞.
7. Our assumptions will guarantee the existence of ρ (see Section 3). Since Y is a Polish space (because it is separable Hilbert) the ρ(y|µa ) conditional distribution (y ∈ Y , µa ∈ X) is also well-defined
(Steinwart and Christmann, 2008, Lemma A.3.16, page 487).

8

Learning Theory for Distribution Regression

2.2 Distribution Regression
We now formally define the distribution regression task. Let us assume that M+
1 (X) is
endowed with S1 = B(τw ), the weak-topology generated σ-algebra; thus (M+
(X),
S1 ) is
1
a measurable space. In the distribution regression problem, we are given samples ẑ =
i.i.d.
+
l
l
i
{({xi,n }N
n=1 , yi )}i=1 with xi,1 , . . . , xi,Ni ∼ xi where z = {(xi , yi )}i=1 with xi ∈ M1 (X)
and yi ∈ Y drawn i.i.d. from a joint meta distribution M defined on the measurable
space (M+
1 (X) × Y, S1 ⊗ B(Y )), the product space enriched with the product σ-algebra.
Unlike in classical supervised learning problems, the problem at hand involves two levels of
randomness, wherein first z is drawn from M, and then ẑ is generated by sampling points
from xi for all i = 1, . . . , l. The goal is to learn the relation between the random distribution
x and response y based on the observed ẑ. For notational simplicity, we will assume that
N = Ni (∀i).
As in the classical regression problem (Rd → R), distribution regression can be tackled
via kernel ridge regression (using a squared loss as the discrepancy criterion). The kernel
(say KG ) is defined on M+
1 (X), and the regressor is then modelled by an element in the
RKHS G = G(KG ) of functions mapping from M+
1 (X) to Y . In this paper, we choose
(X)
and
so
that the function (in G) to describe the
KG (x, x′ ) = K(µx , µx′ ) where x, x′ ∈ M+
1
(x, y) random relation is constructed as a composition f ◦ µx , i.e.
f ∈H=H(K)

µ

→ X(⊆ H = H(k)) −−−−−−−→ Y.
M+
1 (X) −

(7)

In other words, the distribution x ∈ M+
1 (X) is first mapped to X ⊆ H by the mean
embedding µ, and the result is composed with f , an element of the RKHS H.
Let the expected risk for a f˜ : X → Y (measurable) function be defined as
 
2
R f˜ = E(x,y)∼M f˜(µx ) − y Y ,

which is minimized by the fρ regression function. The classical regularization approach is
to optimize
l

fzλ = arg min
f ∈H

1X
kf (µxi ) − yi k2Y + λ kf k2H
l

(8)

i=1

instead of R, based on samples z. Since z is not available, we consider the objective function
defined by the observable quantity ẑ,
l

fẑλ = arg min
f ∈H

1X
kf (µx̂i ) − yi k2Y + λ kf k2H ,
l

(9)

i=1

P
N
where x̂i = N1 N
n=1 δxi,n is the empirical distribution determined by {xi,n }i=1 . The ridge
regression objective function has an analytical solution: given training samples ẑ, the prediction for a new t test distribution is
(fẑλ ◦ µ)(t) = k(K + lλI)−1 [y1 ; . . . ; yl ],

(10)

where k = [K(µx̂1 , µt ), . . . , K(µx̂l , µt )] ∈ L(Y )1×l , K = [K(µx̂i , µx̂j )] ∈ L(Y )l×l ,
[y1 ; . . . ; yl ] ∈ Y l .
9

Szabó et al.

Remark 1
• It is important to note that the algorithm has access to the sample points only via
their mean embeddings {µx̂i }li=1 in Eq. (9).
• There is a two-stage sampling difficulty to tackle: The transition from fρ to fzλ represents the fact that we have only l distribution samples (z); the transition from fzλ
to fẑλ means that the xi distributions can be accessed only via samples (ẑ).
• While ridge regression can be performed using the kernel KG , the two-stage sampling makes it difficult to work with arbitrary KG . By contrast, our choice of
KG (x, x′ ) = K(µx , µx′ ) enables us to handle the two-stage sampling by estimating
µx with an empirical estimator, and using it in the algorithm as shown above.
• In case of scalar output (Y = R), L(Y ) = L(R) = R and (10) is a standard linear
equation with K ∈ Rl×l , k ∈ R1×l . More generally, if Y = Rd , then L(Y ) = L(Rd ) =
Rd×d and (10) is still a finite-dimensional linear equation with K ∈ R(dl)×(dl) and
k ∈ Rd×(dl) .
• One could also formulate the problem (and get guarantees) for more abstract X ⊆
H → 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

E fẑλ , fρ := R[fẑλ ] − R[fρ ],

both when fρ ∈ H (the well-specified case) and fρ ∈ L2ρX \H (the misspecified case); second,


2 → 0,
to establish consistency (E fẑλ , fρ → 0, or in the misspecified case E fẑλ , fρ − DH
2
∗
2 := inf
where DH
q∈H kfρ − SK qkρ is the approximation error of fρ by a function in H); and
third, to derive an exact computational-statistical efficiency trade-off as a function of the
(l, N, λ) triplet, and of the difficulty of the problem.

3. Assumptions
In this section, we detail our assumptions on the (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 = µ M+
1 (X) ),
whose properties as a function of X and H must themselves be established.
Our assumptions are as follows:
1. (X, τ ) is a separable, topological space.

10

Learning Theory for Distribution Regression

2. Y is a separable Hilbert space.
3. k is bounded, in other words ∃Bk < ∞ such that supu∈X k(u, u) ≤ Bk , and continuous.
4. The {Kµa }µa ∈X operator family is uniformly bounded in Hilbert-Schmidt norm and
Hölder continuous in operator norm. Formally, ∃BK < ∞ such that

kKµa k2L2 (Y,H) = T r Kµ∗a Kµa ≤ BK ,

(∀µa ∈ X),

(11)

and ∃L > 0, h ∈ (0, 1] such that the mapping K(·) : X → L(Y, H) is Hölder continuous:
kKµa − Kµb kL(Y,H) ≤ L kµa − µb khH ,

∀(µa , µb ) ∈ X × X.

(12)

5. y is bounded: ∃C < ∞ such that kykY ≤ 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).

4. 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.
high probability upper
Given the estimator (fẑλ ) in Eq. (9), we derive finite-sample

λ
setting,
bounds (see Theorems 2 and 7) for the excess risk E fẑ , fρ , and in the misspecified

2 . We
for the excess risk compared to the best attainable value from H, i.e., E fẑλ , fρ − DH
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 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.
4.1 Results for the Well-specified Case
We first focus on the well-specified case (fρ ∈ H) and present our
 first main result. We
λ
derive a high probability upper bound for the excess risk E fẑ , fρ 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.
l
Let us define x = {xi }li=1 and x̂ = {{xi,n }N
n=1 }i=1 as the ‘x-part’ of z and ẑ, respectively.
λ
One can express fz [Eq. (8)] (Caponnetto and De Vito, 2007), and similarly fẑλ [Eq. (9)],
11

Szabó et al.

as
l

fzλ

−1

= (Tx + λ)

1X
Tx =
Tµxi ,
l

gz ,

fẑλ = (Tx̂ + λ)−1 gẑ ,

Tx̂ =

1
l

i=1
l
X

Tµx̂i ,

i=1

l

1X
gz =
Kµxi yi ,
l
gẑ =

1
l

i=1
l
X

Kµx̂i yi ,

(13)

(14)

i=1

where Tµa = Kµa Kµ∗a ∈ L(H) (µa ∈ X), Tx , Tx̂ : H → H, gz , gẑ ∈ H. By these explicit expressions, one can decompose the excess risk into 5 terms (Szabó et al., 2015, Section A.1.8):

 
E fẑλ , fρ = R fẑλ − R [fρ ] ≤ 5 [S−1 + S0 + A(λ) + S1 + S2 ] ,

where

√
S−1 = S−1 (λ, z, ẑ) = k T (Tx̂ + λI)−1 (gẑ − gz )k2H ,
√
S0 = S0 (λ, z, ẑ) = k T (Tx̂ + λI)−1 (Tx − Tx̂ )fzλ k2H ,
√
√
A(λ) = k T (f λ − fρ )k2H , S1 = S1 (λ, z) = k T (Tx + λI)−1 (gz − Tx fρ )k2H ,
√
S2 = S2 (λ, z) = k T (Tx + λI)−1 (T − Tx )(f λ − fρ )k2H ,
Z
∗
Tµa dρX (µa ) = SK SK
: H → H.
f λ = arg min(R[f ] + λ kf k2H ), T =
f ∈H

(15)
(16)

(17)

X

Three of the terms (S1 , S2 , A(λ)) are identical to the terms in Caponnetto and De Vito
(2007), hence the earlier bounds can be applied. The two new terms (S−1 , S0 ) 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:
Theorem 2 (Finite-sample excess risk bounds; well-specified case) Let
K(·) : X → L(Y, H) be Hölder continuous with constants L, h. Let l ∈ Z+ , N ∈ Z+ ,
0 < λ, 0 < η < 1, 0 < δ, Cη = 32 log 2 (6/η),
as
√ kykY ≤ C (a.s.) and A(λ) be the residual
λ − f k2
,
T
as
in
(17),
B(λ)
=
kf
defined above. Define M = 2(C + kfρ kH BK ), Σ = M
ρ
H
2
as the reconstruction error, and N(λ) = T r[(T + λI)−1 T ] as the effective dimension. Then
with probability at least 1 − η − e−δ , the excess risk can be upper bounded as




(

E fẑλ , fρ ≤ 5


2h

p
(2Bk )h 
4L2 1 + log(l) + δ
λN h

C 2 + 4BK ×

 2



   2
B(λ) BK A(λ)
Σ2 N(λ)
64 M BK
6
24 4BK
2
+
+
× log
+ 2
+ B(λ) + kfρ kH
η
λ
l2 λ
l
λ
l2
l
)
 2
BK B(λ) BK A(λ) BK M 2 Σ2 N(λ)
+
+ 2
+
+A(λ) + Cη
l2 λ
4lλ
l λ
l
2

if l ≥ 2Cη BK N(λ)/λ, λ ≤ kT kL(H) and N ≥ 1 +
12

p

2 h+6
1
2
2
log(l) + δ 2 h Bk (BK ) h L h /λ h .

Learning Theory for Distribution Regression

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 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 P(b, c) class: Let us fix the positive constants R, α, β. Then given
1 < b, c ∈ (1, 2], the P(b, c) class is the set of probability distributions ρ on Z = X × Y
such that
c−1
1. a range space assumption is satisfied: P
∃g ∈ H s.t. fρ = T 2 g with kgk2H ≤ R,
∞
2. in the spectral decomposition of T = ∞
n=1 λn h·, en iH en , where (en )n=1 is a basis of
⊥
b
Ker(T ) , the eigenvalues of T satisfy α ≤ n λn ≤ β (∀n ≥ 1).
Remark 3 We make few remarks about the P(b, c) class:
• Range space assumption on fρ : The smoothness of fρ is expressed as a range space
assumption, which is slightly different from the standard smoothness conditions appearing in non-parametric function estimation. By theP
spectral decomposition of T
r
given above [λ1 ≥ λ2 ≥ . . . > 0, limn→∞ λn = 0], T r u = ∞
n=1 (λn ) hu, en iH en (r =
c−1
2 ≥ 0, u ∈ H) and
Im(T r ) =

n X∞

n=1

cn en :

X∞

n=1

o
c2n λ−2r
<
∞
.
n

(18)

Specifically, in the limit as r → 0, we obtain fρ ∈ Im(T 0 ) = Im(I) = H (no constraint); larger values of r give rise to faster decay of the (cn )∞
n=1 Fourier coefficients.
This is the concrete meaning of fρ ∈ 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
mi and identical variance (xi = N (mi , σ 2 I)). When Gaussian kernels (k) are used
2
with linear K, then K(µxi , µxj ) = e−ckmi −mj k (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
mi . The spectral decay of the Gram matrix will correspond to that of the Gaussian
kernel, with points drawn from the meta-distribution over the mi . 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 P(b, c) family, the behaviour of A(λ), B(λ) and N(λ) is known: A(λ) ≤ Rλc , B(λ) ≤
1
b
λ− b . Specializing Theorem 2 and retaining its assumptions, we get:
Rλc−1 , N(λ) ≤ β b−1
Theorem 4 (Finite-sample excess risk bound for ρ ∈ P(b, c))
Suppose the conditions in Theorem 2 hold. Let ρ ∈ P(b, c), where 1 < b and c ∈ (1, 2].
13

Szabó et al.

Then




(

E fẑλ , fρ ≤ 5


2h
p
"
4L2 1 + log(l) + δ
(2Bk )h
λN h

C 2 + 4BK ×

!#
#
( "
 2
)
Rλc−1 BK Rλc
Σ2 βb
3 4BK
2 M 2 BK
+
+
+ Rλc−1 + kfρ k2H
+ 2
× Cη
1
λ
l2 λ
4λ
l2
l
(b − 1)lλ b
"
#)
2 Rλc−2
c−1
2
2 βb
B
B
Rλ
B
M
Σ
K
K
K
+ Rλc + Cη
.
+
+ 2
+
1
l2
4l
l λ
(b − 1)lλ b

Discarding the constants in Theorem 4, the study of convergence of the excess risk E(fẑλ , fρ )
to 0 boils down to finding N and λ (as a function of l) where N → ∞, λ → 0 and


b+1
1
1
1
log(l)
logh (l)
1
+1+
r(l, N, λ) =
≤N
+ λc + 2 + 1 → 0, s.t. lλ b ≥ 1,
2
1+ 1b
N hλ
λ2 l 2
l
λ
lλ
lλ b
λh
(19)
a

as l → ∞. Let us choose N = l h log(l); in this case Eq. (19) reduces to
r(l, λ) =

1

l2+a λ3

+

1

la λ

+

1

1
la+1 λ2+ b

+ λc +

1

l2 λ

+

1

lλ

1
b

→ 0, s.t. lλ

b+1
b

≥ 1, la λ2 ≥ 1. (20)

One can assume that a > 0, otherwise r(l, λ) → 0 fails to hold; in other words, N should
grow faster than log(l). Matching the ‘bias’ (λs ) and ‘variance’ (other) terms in r(l, λ)
to choose λ, 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:8
Theorem 5 (Computational-statistical efficiency trade-off; well-specified case;
a
ρ ∈ P(b, c)) Suppose the conditions in Theorem 2 hold. Let ρ ∈ P(b, c) and N = l h log(l),
where 0 < a, 1 < b, c ∈ (1, 2]. If
 ac 

−
− a
λ, f
• a ≤ b(c+1)
,
then
E
f
=
O
ρ
p l c+1 with λ = l c+1 ,
ẑ
bc+1
• a≥

b(c+1)
bc+1

 bc 

b
then E fẑλ , fρ = Op l− bc+1 with λ = l− bc+1 .

Remark 6 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

• a-dependence: A smaller bag size (smaller a; N = l h log(l)) means computational
savings, but reduced statistical efficiency. It is not worth increasing a above b(c+1)
bc+1
bc

since from that point the rate becomes r(l) = l− bc+1 ; remarkably, this rate is minimax
in the one-stage sampled setup (Caponnetto and De Vito, 2007). The sensible choice
a = 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.
8. The derivations are available in the supplement.

14

Learning Theory for Distribution Regression

• h-dependence: In accord with our ‘smoothness’ assumptions it is rewarding to use
a
smoother K kernels (larger h ∈ (0, 1]) since this reduces the bag size [N = l h log(l)].
• c-dependence: The strictly decreasing property of c 7→ 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)
1. Decomposition of the excess risk: We have the following upper bound for the excess
risk

 
E fẑλ , fρ = R fẑλ − R [fρ ] ≤ 5 [S−1 + S0 + A(λ) + S1 + S2 ] .

(21)

2. It is sufficient to upper bound S−1 and S0 : Caponnetto and De Vito (2007) have
shown that for ∀η > 0 if l ≥ 2Cη BK N(λ)/λ, λ ≤ kT kL(H) , then P(Θ(λ, z) ≤ 1/2) ≥
1 − η/3, where
Θ(λ, z) = (T − Tx )(T + λI)−1

L(H)

,

(22)

using which upper bounds on S1 and S2 that hold with probability 1 − η are obtained.
It is known that A(λ) ≤ Rλc .
√
3. Probabilistic bounds on kgẑ − gz k2H , kTx − Tx̂ k2L(H) , k T (Tx̂ + λI)−1 k2L(H) , kfzλ k2H :
One can bound S−1 and S0 as
S−1 ≤
and
S0 ≤

√

√

2
kg
L(H) ẑ

T (Tx̂ + λI)−1

T (Tx̂ + λI)−1

2
kTx
L(H)

− gz k2H

− Tx̂ k2L(H) fzλ

2
.
H

For the terms on the r.h.s., we derive upper bounds [for the definition of α, see Eq. (24)]
kgẑ −

gz k2H

2

≤L C

2 (1

kTx −

+

√

2h

√

α) (2Bk )h
,
Nh

Tx̂ k2L(H)

≤

(1 +

√

T (Tx̂ + λI)−1

2h h+2
2 (B

α)

Nh

h
2
k ) BK L

L(H)

2
≤√ ,
λ

,

and
fzλ

2
H

≤6



  2

M BK
Σ2 N(λ)
16
2 6
log
+
λ
η
l2 λ
l
  2


4BK B(λ) BK A(λ)
4
2
2 6
+
+ B(λ) + kfρ kH .
+ 2 log
λ
η
l2
l

The bounds hold under the following conditions:
15

(23)

Szabó et al.

• kgẑ − gz k2H (see Section 7.1.1): if the empirical mean embeddings are close to their
population counterparts, i.e.,
kµxi − µx̂i kH ≤

(1 +

√ √
α) 2Bk
√
,
N

(∀i = 1, . . . , l).

This event has probability 1 − le−α over all i = 1, . . . , l samples;
(Altun and Smola, 2006) and (Szabó et al., 2015, Section A.1.10).

(24)
see

• kTx − Tx̂ k2L(H) (see Section 7.1.2): (24) is assumed.
√
• k T (Tx̂ + λI)−1 k2L(H) (Szabó et al., 2015, Section A.1.11): (24), Θ(λ, z) ≤ 21 , and
(1 +

√

h+6
h

2

α) 2

1

2

Bk (BK ) h L h

2

λh

≤ N.

(25)

• kfzλ k2H : The bound is guaranteed to hold under the conditions of the bounds of S1
and S2 .8
4. Union bound: By applying an α = log(l) + δ reparameterization, and combining the
received upper bounds with Caponnetto and De Vito (2007)’s results for S1 and S2 ,
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 anal2
ysis: bounding S1 , S2 , Θ(λ, z), fzλ H (Caponnetto and De Vito, 2007) and kµxi − µx̂i kH
(Altun and Smola, 2006).9
4.2 Results for the Misspecified Case
In this section, we focus on the misspecified case (fρ ∈ L2ρX \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 E fẑ , fρ , 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:
1. 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 E fẑλ , fρ can
be driven to its smallest possible value determined by the “richness” of H:
2
∗
DH
:= inf kfρ − SK
qk2ρ .

(26)

q∈H

The value of DH equals the approximation error of fρ by a function from H. Specifi∗ (H) = {S ∗ q : q ∈ H} ⊆ L2 ] is dense in L2 , then D = 0.
cally, if H [precisely SK
H
ρX
ρX
K
9.

We also corrected some
Caponnetto and De Vito, 2007).

constants

in

16

the

previous

works

(Altun and Smola,

2006;

Learning Theory for Distribution Regression

2. The second, specialized result [Eq. (28)] under additional smoothness assumptions on
fρ 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:
Theorem 7 (Finite-sample excess risk bounds; misspecified case) Let l ∈ Z+ ,

2
 
K
C
≤ l
N ∈ Z+ , 0 < λ, 0 < η < 1, 0 < δ and Cη = log η6 . Assume that 12B
η
λ
p
2 h+6
1
2
2
and 1 + log(l) + δ 2 h Bk (BK ) h L h /λ h ≤ N .
1. Then for arbitrary q ∈ H with probability at least 1 − η − e−δ
h

p
h
√

q
(2Bk ) 2 
2LC
1
+
log(l)
+
δ

2 BK
λ
1+ √
+
(27)
E fẑ , fρ ≤
√
h
λ
λN 2
(
! 
)
√
√

2Cη
C BK
σ 1q
2BK
2C B K
√
+ √
+√
+
λ kfρ kρ Da (λ, q) + Da (λ, q),
l
l
λ
l
l λ
1

∗ qk + max(1, kT k
where Da (λ, q) = kfρ − SK
ρ
L(H) )λ 2 kqkH .

2. In addition, suppose fρ ∈ Im(T̃ s ) for some s > 0, where T̃ is defined in Eq. (6). Then
with probability at least 1 − η − e−δ , we have
h

p
h
√

q
(2Bk ) 2 
2LC
1
+
log(l)
+
δ

2 BK
λ
1+ √
E fẑ , fρ ≤
+
√
h
λ
λN 2
(
! 
√
√

2BK
2Cη
2C B K
C BK
σ 1
√
+
+ √
+√
×
l
l
λ
l
l λ
s
)


s
max 1, T̃
λ T̃ −s fρ Db (λ, s) + Db (λ, s), (28)
L(L2ρ )
ρ
X
s−1
where Db (λ, s) = max(1, kT̃ kL(L
2

ρX )

)λmin(1,s) kT̃ −s fρ kρ .

Remark 8 We give a short insight into the assumptions of Theorem 7, followed by consequences of the theorem.
• Range space assumption on fρ : The range space assumption for the compact, positive, self-adjoint operator, T̃ = T̃ (K) : L2ρX → L2ρX in the 2nd part of Theorem 7
can be interpreted similarly to that on T ; see Eq. (18). One can also prove alternative descriptions for Im(T̃ s ) in terms of interpolation spaces (Steinwart and Scovel,
2012, Theorem 4.6,
 page 387), or the decay
 of the 2-approximation error function,
2
A2 (λ) = inf f ∈H(K) λ kf kH(K) + R[f ] − R[fρ ] (Smale and Zhou, 2003; Steinwart et al.,
2009).
17

Szabó et al.

q


E fẑλ , fρ : 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 Da (λ, 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
∗ qk in the last term was one):
multiplier of kfρ − SK
ρ
q
√
∗ q
h
f
−
S
+
λ kqkH
p
ρ
√
2
K
log (l)
1
ρ
∗
√
+√ +
+ kfρ − SK
qkρ + λ kqkH .
r(l, N, λ, q) =
h
lλ
λ l
N 2λ
p
By choosing N = l1/h log l, r(l, λ) is bounded by
q


p
∗ q




f ρ − SK
kqkH √
ρ
1
∗
√
+ 3 √ + λkqkH + Op √
.
inf kfρ − SK qkρ +

q∈H 
λ l
λl
λ4 l

•

√
Our goal is to investigate the behavior
of the bound
→ 0 and λ l →
n
o
q as l → ∞, λ p
∗ qk + α
∗ q
f
−
S
+
β
∞. Define K(α, β, γ) := inf q∈H kfρ − SK
.
kqk
+
γkqk
ρ
H
H
K
ρ
ρ

K(α, β, γ) is the pointwise infimum of affine functions, therefore it is upper semicontinuous and concave on R3 (Aliprantis and Border, 2006, Lemmas 2.41 and 5.40 ); it
is continuous on ×3i=1 R>0 (Rockafellar and Wets, 2008, Theorem 2.35). Moreover, by applying (Rockafellar and Wets, 2008, Corollary 2.37) it extends continuously to ×3i=1 R≥0 ;
specifically it 
is continuous at
 (α, β, γ) = 0. In other words, as l → ∞, λ → 0 and
√
√
1
1
√ , 3 √ , λ → DH and we get consistency in the misspecified case,10
λ l → ∞, K
λ l λ4

l

p

r(N, l, λ) → DH .

Discarding the constants in Eq. (28) we get10
√
h
p
λmin(1,s)
log 2 (l)
1
1
√
r(l, N, λ) =
+√ +
+ λmin(1,s) , subject to 2 ≤ l.
h
λ
lλ
λ l
N 2λ

(29)

Our goal is to drive r(l, N, λ) to zero with a suitable choice of the (l, N, λ) triplet under
the stronger range space assumption. Since in Eq. (29) min(1, s) appears,
one can assume

1
s
without loss of generality that s ∈ (0, 1]; consequently 1 − 2 ∈ 2 , 1 and 1 1 1 ≤ 1−1s 1 .
l2 λ2

λ

2l2

Let us choose N = l2a/h log(l); in this case using the previous dominance note, Eq. (29)
reduces to the study of
p
1
1
+ λs → 0, s.t. lλ2 ≥ 1.
(30)
r(l, λ) = a +
l λ λ1− 2s l 21

One can assume that a > 0, otherwise r(l, λ) → 0 fails to hold: in other words, N should
grow faster than log(l). Matching the ‘bias’ (λs ) and ‘variance’ (other) terms in r(l, λ) to
choose λ, 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
2

10. We have discarded the log(l)/λ h ≤ N constraint implied by the convergence of the first term in

18

√

r.

Learning Theory for Distribution Regression

Theorem 9 (Computational-statistical efficiency trade-off; misspecified case,
2a
fρ ∈ Im(T̃ s )) Suppose that fρ ∈ Im(T̃ s ) and N = l h log(l), where s ∈ (0, 1], a > 0.
If
 2sa 

a
−
λ, f
with λ = l− s+1 ,
• a ≤ s+1
,
then
E
f
=
O
ρ
p l s+1
ẑ
s+2
 2s 

1
λ, f
• a ≥ s+1
,
then
E
f
l− s+2 with λ = l− s+2 .
=
O
ρ
p
ẑ
s+2

Remark 10 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 = l2a/h log(l)) leads to a compus+1
, the rate
tational advantage, but one looses in statistical efficiency. As ‘a’ reaches s+2
2s

becomes r(l) = l− s+2 and one does not gain from further increasing the value of a. The
s+1
sensible choice of a = s+2
≤ 32 means that N can again be sub-quadratic (2a < 34 < 2)
in l.
• h-dependence: By using smoother K kernels (larger h ∈ (0, 1]) one can reduce the size
of the bags: h 7→ 2a/h is decreasing in h. This is compatible with our smoothness
requirement on fρ .
• s-dependence: “Easier” tasks (larger s) give rise to faster convergence. Indeed, in the
2s
2s
exponent is strictly increasing function of the problem
r(l) = l− s+2 rate the s 7→ s+2
difficulty (s). For example, for extremely non-smooth regression problems (s ≈ 0)
2s
the convergence can be arbitrary slow (lims→0 s+2
= 0). In the smooth case (s = 1)
lims→1

2s
s+2

=

2
3

2

and one can achieve the r(l) = l− 3 rate.
2s

2s

• We may compare our r(l) = l− s+2 result with the ro (l) = l− 2s+1 (one-stage sampled)
rate (Steinwart et al., 2009, β/2 := s, q := 2, p := 1 in Corollary 6), which was
shown to be asymptotically optimal on Y = R for continuous k on compact metric
X. Steinwart et al.’s result is more general in terms of q (kf kqH based regularization)
and p (kf k∞ ≤ C kf kpH kf kρ1−p , ∀f ∈ H; in our case p = 1), although it imposes
an additional eigenvalue constraint [(Steinwart et al., 2009, Eq. (6))] as well as fρ ∈
Im(T̃ s ). Moreover, one can observe that ro (l) ≤ r(l) with a small gap, and that
for s → 0 and s = 1, ro (l) = r(l); see Fig. 1. We further remind the reader that our
MERR analysis also holds for separable Hilbert output spaces Y , separable topological
domains X enriched with a bounded, continuous kernel k, and that we handle the twostage sampled setting.
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,11 and we present them jointly.
1. Decomposition of the excess risk: By the triangle inequality, we have
q


∗ λ
∗
∗ λ
fẑ − fρ ρ ≤ SK
E fẑλ , fρ = SK
fẑλ − fzλ ρ + SK
fz − fρ ρ .

(31)

11. 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).

19

Szabó et al.

0.8

Rate

0.6

0.4
−log [r (l)]=2s/(2s+1)
l o

0.2

−log [r(l)]=2s/(s+2)
l

0
0

0.2

0.4
0.6
Smoothness (s)
2s

0.8

1

2s

Figure 1: Comparison of the ro (l) = l− 2s+1 and r(l) = l− s+2 rates as function of the problem
difficulty/smoothness (s).
∗ fλ − fλ
2. Bound on SK
z
ẑ



: Using12 the fact that
√
2
∗
T h H (∀h ∈ H),
kSK
hk2ρ =
ρ

and the definitions of S−1 and S0 [see Eqs. (15)-(16)], we obtain
p
√
p


∗
SK
fẑλ − fzλ ρ =
T fẑλ − fzλ H ≤ S−1 + S0 ,

(32)

(33)

through an application of triangle inequality. One can derive without a P(b, c) prior
assumption (Section 7.2.1) the upper bound13
√

√
h 
p
p
2LC(1 + α)h (2Bk ) 2
2 BK
S−1 + S0 ≤
1+ √
√
h
λ
λN 2

for the r.h.s. of Eq. (33) under the conditions that Θ(λ, z) ≤ 21 (which holds with
probability 1 − η if [12BK log(2/η)/λ]2 ≤ l), and that Eqs. (24)-(25) hold.
∗ fλ − f
3. Decomposition of SK
ρ ρ : By the triangle inequality and Eq. (32), we have
z


∗ λ
∗
∗ λ
∗
∗ λ
SK
f z − f ρ ρ = SK
fzλ − f λ + SK
f − f ρ ρ ≤ SK
fzλ − f λ ρ + SK
f − fρ ρ
√

∗ λ
T fzλ − f λ H + SK
f − fρ ρ .
(34)
=

√

T fzλ − f λ H : Making use of the analytical expressions for fzλ
4. Decomposition of
and f λ [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)

√
√

T fzλ − f λ H ≤
T (Tx + λI)−1 L(H) kgz − gρ kH +
 

∗
SK f ρ H ,
kT − Tx kL(H) λ−1 SK fρ − (T̃ + λI)−1 SK
√
where gρ = SK fρ√
. As it is known (Caponnetto and De Vito, 2007, page 348) k T (Tx +
λI)−1 kL(H) ≤ 1/ λ provided that Θ(λ, z) ≤ 21 .
∗
12. See for example de Vito et al. (2006) on page 88 with the (H, G, A, T ) := (H, L2ρX , SK
, T ) choice.
13. See the remark at the end of Section 7.2.1.

20

Learning Theory for Distribution Regression

5. Bound on kgz − gρ kH , kT − Tx kL(H) : By concentration arguments the bounds
!
√
√
 

 

2
2
4C B K
4BK
2C BK
4σ
√
√
kgz − gρ kH ≤
log
log
+
+
, kT − Tx kL(H) ≤
l
η
l
η
l
l
hold with probability at least 1 − η, each (see Section 7.2.3, 7.2.4).
 2

∗ S f
6. Decomposition of SK fρ − (T̃ + λI)−1 SK
K ρ H : Exploiting the analytical formula
for f λ , one can construct (Section 7.2.5) the upper bound


 2

∗
∗
∗ λ
SK fρ − (T̃ + λI)−1 SK
SK fρ H ≤ T̃ fρ − (T̃ + λI)−1 SK
SK f ρ ρ SK
f − fρ ρ .



∗ S f
s
7. Bound on T̃ fρ − (T̃ + λI)−1 SK
K ρ ρ : Using our assumptions that fρ ∈ Im(T̃ )
(s ≥ 0)14 and exploiting the separability of L2ρX , by Lemma 7.3.2 (K = L2ρX , f = fρ ,
∗ S we obtain the upper bound
M = T̃ , a = 1) and T̃ = SK
K




∗
SK fρ ρ = T̃ fρ − (T̃ + λI)−1 T̃ fρ ρ
T̃ fρ − (T̃ + λI)−1 SK




λ T̃ −s fρ ρ ,
≤ max 1, kT̃ ksL(L2 ) λmin(1,s+1) T̃ −s fρ ρ = max 1, kT̃ ksL L2
( ρX )
ρX
where we used at the last step that min(1, s + 1) = 1; this follows from s ≥ 0.

∗ fλ − f
8. Bound on SK
ρ

ρ

:

(a) No range space assumption: One can construct (Section 7.2.6) the bound
 1
∗ λ
∗
SK
f − fρ ρ ≤ kfρ − SK
qkρ + max 1, kT kL(H) λ 2 kqkH ,
which holds for arbitrary q ∈ H.

∗ f λ = (T̃ + λI)−1 T̃ f identity
(b) Range space assumption in L2ρX : Using the SK
ρ
2
[see Eq. (43)], and Lemma 7.3.2 (M = T̃ , K = LρX , a = 0), we get


s−1
∗ λ
SK
f − fρ ρ = (T̃ + λI)−1 T̃ fρ − fρ ρ ≤ max 1, T̃ L(L2 ) λmin(1,s) T̃ −s fρ ρ .
ρX

9. Union bound: Applying an α = log(l) + δ reparameterization, changing η to η3 and
combining the derived results (in case of the first statement with s = 0) with a union
bound, Theorem 7 follows.
Remark 11 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:

∗ fλ − fλ
(a) In Step 2 by using Eq. (32) and transforming the L2ρX error SK
to H,
z
ẑ
ρ
we could rely on our previous bounds for S−1 and S0 . However, we were required to
use a different concentration argument to guarantee Θ(λ, z) ≤ 21 since we no longer
assume the P(b, c) prior class.
14. Note that we choose s = 0 and s > 0 in the first and second theorem part, respectively.

21

Szabó et al.

(b) In Step 4 the first term could be bounded by Caponnetto and De Vito (2007). Its
Θ(λ, z) ≤ 21 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.

5. 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 Rényi and Tsallis divergence estimates (Póczos et al., 2011, 2012),
but these similarity indices are not kernels, and their consistency in specific learning tasks
remains an open question.
22

Learning Theory for Distribution Regression

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, classification): Despite the large number of MIL
applications and the spate of heuristic solution techniques, there exist few theoretical results in the area (Auer, 1998; Long and Tan, 1998; Blum and Kalai, 1998; Babenko et al.,
2011; Zhang et al., 2013; Sabato and Tishby, 2012) and they focus on the multi-instance
classification (MIC) task. In particular, let us first consider the standard MIC assumption
(Dietterich et al., 1997), where a bag is declared to be positive (labelled with “1”) if at least
one of its instances is positive (“1”); otherwise, the bag is negative (“0”).15 In other words,
if the instances (xi,n ) in the ith bag {xi,1 , . . . , xi,N } have hidden label L(xi,n ) ∈ {0, 1}, then
the observed label of the bag is yi = h(xi,1 , . . . , xi,N ) = max(L(xi,1 ), . . . , L(xi,N )) ∈ {0, 1}.
In case of the original APR (axis-aligned rectangles; Dietterich et al., 1997) hypothesis class,
function L is equal to the indicator of an unknown rectangle R (L = IR ). In other words, a
bag is declared to be positive if there exists at least one instance in the bag, which belongs
to R.16 The goal is to learn R with high probability given the bags ({xi,1 , . . . , xi,N }-s) and
their labels (yi -s). Long and Tan (1998) proved the PAC learnability (probably approximately correct; Valiant, 1984) of the APR hypothesis class, if all instances in each bag are
i.i.d. and follow the same product distribution over the instance coordinates. On the other
hand, for arbitrary distributions over bags, when the instances within a bag might be statistically dependent, APR learning under MIC is NP-hard (Auer, 1998); the same authors
also showed that the product property (Long and Tan, 1998) on the coordinates is not required to obtain PAC results. Blum and Kalai (1998) extended PAC learnability of APR-s
to hypothesis classes learnable from one-sided classification noise. In contrast to the previous approaches (Long and Tan, 1998; Auer, 1998; Blum and Kalai, 1998), Babenko et al.
(2011) modelled the bags as low-dimensional manifolds, and proved PAC bounds. By relaxing the standard MIC assumption, Sabato and Tishby (2012) showed PAC-learnability
for general MIC hypothesis classes with extended “max” functions. Zhang et al. (2013)
derived high-probability generalization bounds in the MIC setting, when local and global
representations are combined. Our work falls outside this setting since the label and bag
generation mechanisms we consider are different: we do not assume an exact form of the
15. The motivation of this assumption comes from drug discovery: if a molecule has at least one well-binding
configuration, then it is considered to bind well.
16. In terms of drug binding prediction, this means that a molecule binds to a target iff at least one of its
configurations falls within a fixed, but unknown rectangle.

23

Szabó et al.

labelling mechanism (function L and max in h). Rather, the labelling is presumed to be
stochastically determined by the underlying true distribution, not deterministically by the
instance realizations in the bags (these are presumed i.i.d., and may be bag-specific).
“Bag-of-objects” methods (MIL, not classification): Beyond classification,
there exist several heuristics—without consistency guarantees—for many other multiinstance 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.17 Probably one of the most wellknown set metrics is the Hausdorff metric (Edgar, 1995), which is defined for nonempty 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; Müller, 2005), by considering the probability measures xi as functions. This is a
highly non-standard setup, however, since these functions (xi ) are defined on σ-algebras
and are non-negative, σ-additive.

6. 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; Richtárik and Taká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 well17. Often these “metrics” are only semi-metrics, as they do not satisfy the triangle inequality.

24

Learning Theory for Distribution Regression

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; Gä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 (Szabó et al., 2015, Section A.2) to supervised entropy
learning and aerosol prediction based on multispectral satellite images.18 In future work,
we will address applications with vector-valued outputs.

7. 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.
7.1 Proofs of the Well-specified Case
We give proof details concerning the excess risk in the well-specified case (Theorem 2).
7.1.1 Proof of the bound on kgẑ − gz k2H

P
By (13), (14) we get gẑ − gz = 1l li=1 Kµx̂i − Kµxi yi ; hence by applying the Hölder
property of K(·) , the boundedness of yi (kyi kY ≤ C) and (24), we obtain
l

kgẑ − gz k2H ≤
≤

l
L2 X

l

i=1

1 X
l
l2
i=1

kyi k2Y kµx̂i − µxi k2h
H

l

1X
2
Kµx̂i − Kµxi L(Y,H) kyi k2Y
l
i=1
2h

√ 2h
√ √
l
h
L2 C 2 X (1 + α) 2Bk
2 2 (1 + α) (2Bk )
√
=L C
≤
l
Nh
N
i=1


Kµx̂i − Kµxi yi

2
H

≤

with probability at least 1 − le−α , based on a union bound.
18. For code, see https://bitbucket.org/szzoli/ite/.

25

Szabó et al.

7.1.2 Proof of the bound on kTx − Tx̂ k2L(H)
Using the definition of Tx and Tx̂ , and exploiting (with k · kL(H) ) that in a normed space19
(N, k·k), fi ∈ N , (i = 1, . . . , n)
Xn

i=1

we get

fi

2

≤n

Xn

i=1

kfi k2 ,

l

kTx − Tx̂ k2L(H) ≤

1 X
Tµxi − Tµx̂i
l
l2
i=1

(35)

2
L(H)

.

(36)

To upper bound kTµxi − Tµx̂i k2L(H) , let us see how Tµu = Kµa Kµ∗a acts. The existence of an
E ≥ 0 constant satisfying k(Tµu − Tµv )(f )kH ≤ E kf kH implies kTµu − Tµv kL(H) ≤ E. We
continue with the l.h.s. of this equation using Eq. (35):
2

k(Tµu − Tµv )(f )k2H = Kµu Kµ∗u (f ) − Kµv Kµ∗v (f ) H


2
= Kµu Kµ∗u (f ) − Kµ∗v (f ) + (Kµu − Kµv ) Kµ∗v (f ) H
h

 2
≤ 2 Kµu Kµ∗u (f ) − Kµ∗v (f ) H + (Kµu − Kµv ) Kµ∗v (f )

2
H

By Eq. (45) and the Hölder continuity of K(·) , one arrives at


Kµu Kµ∗u (f ) − Kµ∗v (f )

2
H

≤ kKµu k2L(Y,H) Kµ∗u (f ) − Kµ∗v (f )

i

.

2
Y

2
2
kf k2H = kKµu k2L(Y,H) k(Kµu − Kµv )∗ kL(H,Y ) kf k2H
L(H,Y )
2
= kKµu k2L(Y,H) kKµu − Kµv k2L(Y,H) kf k2H ≤ BK L2 kµu − µv k2h
H kf kH ,
2
2
(Kµu − Kµv ) Kµ∗v (f ) H ≤ kKµu − Kµv k2L(Y,H) Kµ∗v (f ) Y
2
2
≤ kKµu − Kµv k2L(Y,H) Kµ∗v L(H,Y ) kf k2H ≤ BK L2 kµu − µv k2h
H kf kH .

≤ kKµu k2L(Y,H) Kµ∗u − Kµ∗v

2
2h
2
2
Hence k(Tµu − Tµv )(f )k2H ≤ 4BK L2 kµu − µv k2h
H kf kH ⇒ E = 4BK L kµu − µv kH . Exploiting this property in (36) with Eq. (24) we arrive to the bound

kTx −

Tx̂ k2L(H)

√ 2h
l
l
4BK L2 X
4BK L2 X (1 + α) (2Bk )h
2h
≤
kµxi − µx̂i kH ≤
l
l
Nh
i=1
i=1
√ 2h
(1 + α) 2h+2 (Bk )h BK L2
=
.
Nh

(37)

7.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−α = 1 − e−[α−log(l)] = 1 − e−δ applying a union bound
19. Eq. (35) holds since k·k2 is convex function, thus

26

1
n

Pn

i=1

fi

2

≤

1
n

Pn

i=1

kfi k2 .

Learning Theory for Distribution Regression

1
2

argument; α = log(l) + δ), and (iii) Θ(λ, z) ≤

is fulfilled [see Eq. (22)], then

"
√ 2h
√ 2h
4
(1 + α) (2Bk )h (1 + α) 2h+2 (Bk )h BK L2
S−1 + S0 ≤
L2 C 2
+
×
λ
Nh
Nh

 2
   2



B(λ) BK A(λ)
64 M BK
Σ2 N(λ)
24 4BK
2
2 6
+
+
× log
+ 2
+ B(λ) + kfρ kH
η
λ
l2 λ
l
λ
l2
l

√ 2h
4L2 (1 + α) (2Bk )h
C 2 + 4BK ×
=
λN h

   2
 2



B(λ) BK A(λ)
64 M BK
Σ2 N(λ)
24 4BK
2
2 6
× log
+
+
+ 2
+ B(λ) + kfρ kH .
η
λ
l2 λ
l
λ
l2
l
By taking
account Caponnetto
and DeVito
h
i
h 2 (2007)’s bounds
i for S1 and S2 , S1 ≤
  into
4BK B(λ)
Σ2 N(λ)
BK A(λ)
BK M 2
2 6
2 6
+ lλ
, S2 ≤ 8 log η
, plugging all the expres32 log η
l
l2 λ +
l2 λ
sions to (21), we obtain Theorem 2 with a union bound.
7.2 Proofs of the Misspecified Case
We present the proof details concerning the excess risk in the misspecified case (Theorem 7).
7.2.1 Proof of the bound on

√

S−1 +

√

S0 without P(b, c)

The upper bounds on S−1 and S0 [which are defined in Eqs. (15), (16)] remain valid without modification provided that (i) Θ(λ, z) = k(T − Tx )(T + λ)−1 kL(H) ≤ k(T − Tx )(T +
λ)−1 kL2 (H) ≤ 21 , where we used Eq. (1), (ii) Eq. (24) is satisfied (which has probability
1 − le−α ) and (iii) Eq. (25) holds. Our goal below is to guarantee the Θ(λ, z) ≤ 21 condition
with high probability without assuming that the prior belongs to P(b, c).
Requirement Θ(λ, z) ≤ 21 : Let us define ξi = Tµxi (T + λ)−1 ∈ L2 (H), (i = 1, . . . , l). With
P
this choice we get E[ξi ] = T (T + λ)−1 , (T − Tx )(T + λ)−1 = E[ξi ] − 1l li=1 ξi and
kξi kL2 (H) ≤ Tµxi

L2 (H)

(T + λ)−1

L(H)



≤ BK /λ ⇒ E kξi k2L2 (H) ≤ (BK )2 /λ2 ,

where we made use of (2), the kTµxi kL2 (H) ≤ BK 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 K = L2 (H), B = 2BK /λ,
σ = BK /λ) we obtain that for ∀η ∈ (0, 1)
P



(T − Tx )(T + λ)

Thus, for Θ(λ, z) ≤
2



−1

1
2

L2 (H)

≤2



2BK
BK
+√
λl
lλ



 
2
≥ 1 − η.
log
η

with probability 1 − η it is sufficient to have

BK
2BK
+√
λl
lλ



 
 

 2
2
2
12BK
2
6BK
1
log
≤ l.
log
≤ √ log
≤ ⇔
η
η
2
λ
η
lλ
27

(38)

Szabó et al.

Under these conditions, we arrived at the upper bound
s
#
"
r
√ 2h
p
p
4L2 C 2 (1 + α) (2Bk )h √
4BK
S−1 + S0 ≤
1+
λN h
λ
√

√
h 
2LC(1 + α)h (2Bk ) 2
2 BK
√
=
1+
,
√
h
λ
λN 2
2

2

where as opposed to Section 7.1.3 and Eq. (23) we used a slightly cruder fzλ H ≤ Cλ
bound; it holds without the P(b, c) assumption by the definition of fzλ and the boundedness
P
2
of y since λ fzλ H ≤ 1l li=1 kyi k2Y ≤ C 2 .

Remark: Notice that the price we pay for not assuming that the prior belongs to the
1
P(b, c) class (b > 1) is a slightly tighter λ12 ≤ l constraint [Eq. (38)] instead of 1+
1 ≤ l in
Eq. (19), and a somewhat looser fzλ

2
H

λ

7.2.2 Proof of the decomposition of
Using the analytical formula of

fzλ

b

bound.
√

T fzλ − f λ



H

[see Eq. (13)] and that of f λ [see Eq.(17)]

∗
f λ = (SK SK
+ λI)−1 SK fρ = (T + λI)−1 SK fρ

(39)

one gets (T + λI)f λ = SK fρ ⇒ λf λ = SK fρ − T f λ and
fzλ − f λ = (Tx + λI)−1 gz − f λ = (Tx + λI)−1 gz − (Tx + λI)−1 (Tx + λI)f λ



= (Tx + λI)−1 gz − (Tx + λI)f λ = (Tx + λI)−1 gz − Tx f λ − λf λ

= (Tx + λI)−1 gz − Tx f λ − SK fρ + T f λ
= (Tx + λI)−1 (gz − SK fρ ) + (Tx + λI)−1 (T − Tx )f λ

= (Tx + λI)−1 (gz − SK fρ ) + (Tx + λI)−1 (T − Tx )(T + λI)−1 SK fρ .
Let us rewrite (T +λI)−1 by the (A+U V )−1 = A−1 −A−1 U I + V A−1 U
Woodbury formula (Ding and Zhou, 2008, Theorem 2.1, page 724)

−1

(40)

V A−1 operator


−1 ∗ −1
∗ −1
∗
) = (λ−1 I) − (λ−1 I)SK I + SK
(T + λI)−1 = (λI + SK SK
(λ−1 I)SK
SK (λ I)
∗
= (λ−1 I) − λ−1 SK (λI + T̃ )−1 SK
.

−1
−1
−1
−1
By the derived expression
get (T
 + λI) SK fρ = λ SK fρ − λ SK (λI +
 for (T + λI) −1, we
∗
−1
∗
−1
T̃ ) SK SK fρ = λ SK fρ − (T̃ + λI) SK SK fρ . Plugging this result to Eq. (40), introducing the gρ = SK fρ notation, using the triangle inequality we get
√

T fzλ − f λ H =
io
h
n
√
∗
=
SK f ρ
T (Tx + λI)−1 (gz − SK fρ ) + (T − Tx )λ−1 SK fρ − (T̃ + λI)−1 SK
H

√
 

−1 ∗
−1
−1
T (Tx + λI)
f
−
(
T̃
+
λI)
S
S
f
kg
−
g
k
+
kT
−
T
k
λ
S
≤
z
ρ H
x L(H)
K ρ
K K ρ H .
L(H)

28

Learning Theory for Distribution Regression

7.2.3 Proof of the bound on kgz − gρ kH
R
P
As is known gz = 1l li=1 Kµxi yi [see Eq. (13)] and gρ = X Kµx fρ (µx )dρX (µx )
(Caponnetto and De Vito, 2007, Eq. (23), page 344). Let ξi = Kµxi yi ∈ H (i = 1, . . . , l).
P
2
In this case E[ξi ] = gρ , gρ − gz = E[ξi ] − 1l li=1 ξi , and kξi k2H = Kµxi yi H ≤


√
2
Kµxi L(Y,H) kyi k2Y ≤ BK C 2 ⇒ kξi kH ≤ C BK ⇒ E kξi k2H ≤ C 2 BK using the boundedness of kernel K ( Kµxi

2
L(Y,H)

≤ BK ) and the boundedness of output y (k|ykY ≤ C). Ap√
√
plying the Bernstein inequality (see Lemma 7.3.1 with K = H, B = 2C BK , σ = C BK )
one gets that for any η ∈ (0, 1)
!
√
√
 !
C BK
2
2C B K
+ √
≥ 1 − η.
P kgz − gρ kH ≤ 2
log
l
η
l
7.2.4 Proof of the bound on kT − Tx kL(H)
P
Let ξi = Tµxi ∈ L2 (H) (i = 1, . . . , l), then E[ξi ] = T , T − Tx = T − 1l li=1 Tµxi , kξi kL2 (H) =
h
i
2 . Applying the kT − T k
kTµxi kL2 (H) ≤ BK , E kξi k2L2 (H) ≤ BK
x L(H) ≤ kT − Tx kL2 (H)
relation [see Eq. (1)] and the Bernstein inequality (see Lemma 7.3.1 with K = L2 (H),
B = 2BK , σ = BK ), we obtain that for any η ∈ (0, 1)

 


2
2BK
σ
P kT − Tx kL(H) ≤ 2
+ √ log
≥ 1 − η.
l
η
l
∗ S f
7.2.5 Proof of the decomposition of SK fρ − (T̃ + λI)−1 SK
K ρ

∗ S a, ai = T̃ a, a
Since kSK ak2H = hSK a, SK aiH = hSK
K
ρ

ρ



2
H

(∀a ∈ L2ρX ) by the definition of

∗ S [see Eq. (6)], we can rewrite the target term as
the adjoint operator and T̃ = SK
K
i 2
h
∗
=
SK f ρ
SK fρ − (T̃ + λI)−1 SK
H
D h
E
i
∗
∗
= T̃ fρ − (T̃ + λI)−1 SK
SK f ρ
SK fρ , fρ − (T̃ + λI)−1 SK
ρ
i
h
−1 ∗
−1 ∗
fρ − (T̃ + λI) SK SK fρ
≤ T̃ fρ − (T̃ + λI) SK SK fρ
ρ

ρ

,

where the CBS (Cauchy-Bunyakovsky-Schwarz) inequality was applied. Since
∗
∗
∗
∗
(SK
SK + λI)SK
= SK
(SK SK
+ λI)

∗
∗
∗
∗
SK
(SK SK
+ λI)−1 = (SK
SK + λI)−1 SK

∗
∗
∗
∗
SK
(SK SK
+ λI)−1 SK = (SK
SK + λI)−1 SK
SK

(41)

∗
SK
(T + λI)−1 SK = (T̃ + λI)−1 T̃

(42)

using Eq. (41) and the analytical expression for f λ [see Eq. (39)] we have
∗
∗
∗
(T̃ + λI)−1 T̃ fρ = (T̃ + λI)−1 SK
SK fρ = (SK
SK + λI)−1 SK
SK f ρ
∗
∗
∗ λ
= SK
(SK SK
+ λI)−1 SK fρ = SK
f


 2

−1 ∗
∗ S f
and SK fρ − (T̃ + λI)−1 SK
K ρ H ≤ T̃ fρ − (T̃ + λI) SK SK fρ

29

(43)
ρ

∗ fλ − f
SK
ρ

ρ

.

Szabó et al.

∗ fλ − f
7.2.6 Proof of the bound on SK
ρ

q′

ρ

Let us apply (i) the Af − f = Af − f − + = (A − I)(f − q ′ ) + Aq ′ − q ′ relation with
∗ q, where q ∈ H is an arbitrary element from H, (ii)
A = (T̃ + λI)−1 T̃ , f = fρ and q ′ = SK
Eq. (43) and (iii) the triangle inequality to arrive at
∗ λ
SK
f − fρ

ρ

q′

= (T̃ + λI)−1 T̃ fρ − fρ ρ


∗
∗
∗
= (T̃ + λI)−1 T̃ − I (fρ − SK
q) + (T̃ + λI)−1 T̃ SK
q − SK
q ρ


∗
∗
∗
≤ (T̃ + λI)−1 T̃ − I (fρ − SK
q) ρ + (T̃ + λI)−1 T̃ SK
q − SK
q

ρ

.

Below we give upper bounds on these two terms.
First, notice that µx ∈ X 7→ kK(µx , µx )kL(Y ) ≤ BK . This boundedness with the strong
continuity of K(·) imply (Carmeli et al., 2006, Proposition 12) that H ⊆ C(X, Y ), i.e., K
is a Mercer kernel. Since Kµx is a Hilbert-Schmidt operator for all µx ∈ X [see Eq. (11)],
it is also a compact operator (∀µx ∈ X). The compactness of Kµx -s with the bounded and
∗ and that T̃ is a compact, positive,
Mercer property of K guarantees the boundedness of SK
self-adjoint operator (Carmeli et al., 2010, Proposition 3).
∗ q)k : Since T̃ is a compact positive self-adjoint opBound on k[(T̃ + λI)−1 T̃ − I](fρ − SK
ρ
erator, by the spectral theorem (Steinwart and
2008, Theorem 4.27, page 127)
 Christmann,

there exist an (ui )i∈I countable ONB in cl Im(T̃ ) , and a1 ≥ a2 ≥ . . . > 0 such that
P
T̃ f = i∈I ai hf, ui iρ ui (∀f ∈ L2ρX ) and let (vj )j∈J (J is also countable by the separabil

ity20 of L2ρX ) an ONB in Ker(T̃ ∗ ) = Ker(T̃ ); L2ρX = cl Im(T̃ ) ⊕ Ker(T̃ ). Thus,



∗
(T̃ + λI)−1 T̃ − I (fρ − SK
q)

2
ρ

2
X
ai
∗
∗
− 1 hfρ − SK
q, ui i2ρ +
hfρ − SK
q, vj i2ρ
ai + λ
j∈J
i∈I
X
X
2
∗
∗
∗
≤
hfρ − SK q, ui iρ +
hfρ − SK q, vj i2ρ = kfρ − SK
qk2ρ

=

X
i∈I

exploiting the Parseval’s identity and that

j∈J

λi
λi +λ

2
− 1 ≤ 1.

∗ q − S ∗ qk : By using Eq. (42), Eq. (32), and Lemma 7.3.2
Bound on k(T̃ + λI)−1 T̃ SK
K ρ
∗ , K = H, f = q, a = 1 ), the target quantity can be bounded as
(M = T = SK SK
2
∗
∗
(T̃ + λI)−1 T̃ SK
q − SK
q

ρ

∗
∗
∗
= SK
(T + λI)−1 SK SK
q − SK
q ρ
√ 

∗
T (T + λI)−1 SK SK
=
q−q H

 1
√ 

T (T + λI)−1 T q − q H ≤ max 1, kT kL(H) λ 2 kqkH .
=

Making use of the two derived bounds, we get
 1
max 1, kT kL(H) λ 2 kqkH .

∗ fλ − f
SK
ρ

ρ

∗ qk +
≤ kfρ − SK
ρ

20. L2ρX = L2 (X, B(H)|X , ρX ; Y ) is isomorphic to L2 (X, B(H)|X , ρX ; R)⊗Y , where ⊗ is the tensor product
of Hilbert spaces. The separability follows from that of Y and L2 (X, B(H)|X , ρX ; R); the latter holds
(Cohn, 2013, Proposition 3.4.5) since B(H)|X is countably generated since X ⊆ H is separable.

30

Learning Theory for Distribution Regression

7.3 Supplementary Lemmas
In this section, we list two lemmas used in the proofs.
7.3.1 Bernstein’s inequality (Caponnetto and De Vito, 2007, Prop. 2, p. 345)
Let ξi (i = 1, . . . , l) be i.i.d. realizations of a random variable on a (Ω, A, P ) probability
space with values in a separable
i space K. If there exist B > 0, σ > 0 constants
h Hilbert
2
B
such that kξ(ω)kK ≤ 2 a.s., E kξkK ≤ σ 2 , then for all l ≥ 1 and η ∈ (0, 1) we have
l

P

1X
ξi − E[ξ1 ]
l
i=1

K

≤2



B
σ
+√
l
l



 !
2
≥ 1 − η.
log
η

7.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 K. Let
a ≥ 0, λ > 0, and s ≥ 0. Let f ∈ K such that f ∈ Im (M s ). If s + a > 0, then




λmin(1,s+a) M −s f K .
M a (M + λI)−1 M f − f K ≤ max 1, kM ks+a−1
L(K)

Note: specifically for s = 0 we have Im (M s ) = Im (I) = K, in other words, there is no
additional range space constraint.

8. 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 ρ: The boundedness and continuity of k imply the measurability of µ : (M+
1 (X), B(τw )) → (H, B(H)). Let τ denote the open sets on H = H(k),
τ |X = {A ∩ X : A ∈ τ } the subspace topology on X, and B(H)|X = {A ∩ X : A ∈
B(H)} the subspace σ-algebra on X. By noting (Schwartz, 1998, Corollary 5.2.13)
that B ( τ |X ) = B(H)|X = {A ∈ B(H) : A ⊆ X} ⊆ B(H), the H-measurability of µ
guarantees the measurability of µ : (M+
1 (X), B(τw )) → (X, B(H)|X ), and hence the
well-definedness of ρ, the measure induced by M on X × Y ; for further details see
(Szabó et al., 2015, Section A.1.1).21
• Separability of X: separability of X and the continuity of k implies the separability
of H = H(k) (Steinwart and Christmann, 2008, Lemma 4.33, page 130). Also, since
X ⊆ H, X is separable.
• Finiteness of Bk : If X is compact, then the continuity of k implies Bk < ∞.
• Finiteness of BK , compact metricness of X: Let X be a compact metric space. In
this case M+
1 (X) is also compact metric (Parthasarathy, 1967, Theorem 6.4, page 55).
21. Note that the referred proof also holds for separable Hilbert Y , and by the simplified reasoning above
the original X ∈ B(H) condition could be avoided.

31

Szabó et al.

22 (not just measurable), then X is
Hence if µ : (M+
1 (X), τw ) → H(k) is continuous
compact metric and thus by the Hölder property of K(·) , it is continuous implying
that BK < ∞.

• K properties: It is known (Caponnetto and De Vito, 2007, page 339-340) that
K(µa , µb ) = Kµ∗a Kµb ,
kKµa kL(Y,H) =

(∀µa , µb ∈ X)
p
≤ BK , (∀µa ∈ X).

Kµ∗a L(H,Y )

(44)
(45)

Remark: In terms of Eq. (44), the Eq. (11) assumption means that the
{K(µa , µa )}µa ∈X operators are trace class, specifically they are compact operators.
• Separability of H: The separability of X and the continuity of K imply the separability of H. Indeed, since µa 7→ Kµa is Hölder continuous w.r.t. the Hilbert-Schmidt
norm it is also continuous. As a result it is continuous w.r.t. the operator norm, and
thus also w.r.t. the strong topology. Using this property with the finiteness of BK the
separability of H follows (Carmeli et al., 2006, Proposition 5.1, Corollary 5.2).
• Our assumptions imply Caponnetto and De Vito (2007)’s conditions (not considering
the P(b, c) prior requirement). Indeed
1. Y is a separable Hilbert space by assumption; the same property also holds for H
as we have seen.
2. The measurability of (µx , µt ) 7→ hKµx w, Kµt viH for ∀w, v ∈ Y is guaranteed by the
continuity of K(·) w.r.t. the strong topology.
R
R
3. We have X×Y kyk2Y dρX (µx , y) ≤ X×Y C 2 dρX (µx , y) = C 2 < ∞ due to the
boundedness of y, and hence ∃Σ > 0, ∃M > 0 such that for ρX -almost µx ∈ X
Z
m!Σ2 M m−2
(∀m ≥ 2).
(46)
dρ(y|µ
)
≤
ky − fρ (µx )km
x
Y
2
Y

Indeed, by (Caponnetto and De
(33)) the Bernstein condition (46)
R Vito, 2007, Eq.
2
M
holds if ky − fρ (µx )kY ≤ 2 , Y ky − fρ (µx )kY dρ(y|µx ) ≤ Σ2 . In our case using
the boundedness of y, the regression function is also bounded and the same√
holds
for ky − fρ (µx )kY by the
√ triangle inequality: ky − fρ (µx )kY ≤ C + kfρ kH BK ;
thus, M = 2(C + kfρ kH BK ), Σ = M
2 is a suitable choice.
4. The Polishness of X × Y was used by Caponnetto and De Vito (2007) to assure
the existence of ρ(y|µa ); we guaranteed this existence under somewhat milder conditions (see footnote 7).
Real-valued outputs: We now consider the specific case of Y = R, when the following
simplifications and results hold. By noting that in this case T r(Kµ∗a Kµa ) = K(µa , µa ),
Eq. (11) simplifies to the boundedness of kernel K in the traditional sense
K(µa , µa ) ≤ BK

(∀µa ∈ X).

(47)

22. For example, if k is universal, then µ metrizes the weak topology τw (Sriperumbudur et al., 2010, Theorem 23, page 1552), hence µ is continuous.

32

Learning Theory for Distribution Regression

KG
kµa −µb k2H
−
2θ 2

e

Ke
e−

h=1

KC

kµa −µb kH
2θ 2

h=

1
2



Kt
Ki
− 1
−1 
−1 
2
kµa − µb k2H + θ 2
1 + kµa − µb kθH
1 + kµa − µb k2H /θ 2
h=1

h=

θ
2

(θ ≤ 2)

h=1

Table 1: Nonlinear kernels on mean embedded distributions: K = K(µa , µb ); θ > 0. For the
Hölder continuity of ΨK , we assume that X is a compact metric space and µ is continuous.

Eq. (12) reduces to the Hölder continuity of the canonical feature map ΨK (µc ) :=
K(·, µc ) : X → H, in other words ∃L > 0, h ∈ (0, 1] such that kΨK (µa ) − ΨK (µb )kH ≤
L kµa − µb khH , ∀(µa , µb ) ∈ X × X. In case of a linear kernel, K(µa , µb ) = hµa , µb iH ,
(µa , µb ∈ X), the Hölder continuity of ΨK holds with L = 1, h = 1, and BRK = Bk is a suitable choice. Evaluating the kernel K at the empirical embeddings µx̂i = X k(·, u)dx̂i (u) =
1 PN
n=1 k(·, xi,n ) ∈ H yields the standard set kernel
N
*
+
N
N
N
1 X
1 X
1 X
K(µx̂i , µx̂j ) = µx̂i , µx̂j H =
k(xi,n , xj,m )
k(·, xi,n ),
k(·, xj,m )
= 2
N
N
N
n=1

m=1

H

n,m=1

by the bilinearity of h·, ·iH and the reproducing property of k.
Remark: One can define many nonlinear kernels (see Table 1) on mean embedded distributions. These kernels are the natural extensions to distributions of the Gaussian
(Christmann and Steinwart, 2010), exponential, Cauchy, generalized t-student and inverse
multiquadric kernels. If X is a compact metric space and µ is continuous, then the ΨK
canonical feature maps, associated to K-s in Table 1, can be shown to satisfy our Hölder
continuity requirement [Eq. (12)]; for details, see (Szabó et al., 2015, Section A.1.5-A.1.6).

Acknowledgments
We would like to thank the anonymous reviewers for their highly valuable, constructive
suggestions to improve the manuscript. This work was supported by the Gatsby Charitable
Foundation, NSF grant 1247658, and DOE grant DE-SC001114. A part of the work was
carried out while Bharath K. Sriperumbudur was a research fellow in the Statistical Laboratory, Department of Pure Mathematics and Mathematical Statistics at the University of
Cambridge, UK.

References
Ahmed Alaoui and Michael W. Mahoney. Fast randomized kernel ridge regression with
statistical guarantees. In Advances in Neural Information Processing Systems (NIPS),
pages 775–783, 2015.
Charalambos D. Aliprantis and Kim C. Border. Infinite Dimensional Analysis: A Hitchhiker’s Guide. Springer, 2006.
33

Szabó et al.

Yasemin Altun and Alexander Smola. Unifying divergence minimization and statistical
inference via convex duality. In Conference on Learning Theory (COLT), pages 139–153,
2006.
Mauricio A. Álvarez, Lorenzo Rosasco, and Neil D. Lawrence. Kernels for vector-valued
functions: A review. Foundations and Trends in Machine Learning, 4:195–266, 2011.
Jaume Amores. Multiple instance classification: Review, taxonomy and comparative study.
Artificial Intelligence, 201:81–105, 2013.
Peter Auer. Approximating hyper-rectangles: Learning and pseudorandom sets. Journal of
Computer and System Sciences, 57:376–388, 1998.
Boris Babenko. Multiple instance learning: Algorithms and applications. Technical report, Department of Computer Science and Engineering, University of California, San
Diego, 2004. (http://cms.brookes.ac.uk/research/visiongroup/talks/rg_dec_11_
09/bbabenko_re.pdf).
Boris Babenko, Nakul Verma, Piotr Dollár, and Serge Belongie. Multiple instance learning
with manifold bags. In International Conference on Machine Learning (ICML), pages
81–88, 2011.
Charles Bergeron, Jed Zaretzki, Curt Breneman, and Kristin P. Bennett. Multiple instance
ranking. In International Conference on Machine Learning (ICML), pages 48–55, 2008.
Charles Bergeron, Gregory Moore, Jed Zaretzki, Curt M. Breneman, and Kristin P. Bennett.
Fast bundle algorithm for multiple-instance learning. IEEE Transactions on Pattern
Analysis and Machine Intelligence, 34:1068–1079, 2012.
Alain Berlinet and Christine Thomas-Agnan. Reproducing Kernel Hilbert Spaces in Probability and Statistics. Kluwer, 2004.
Avrim Blum and Adam Kalai. A note on learning from multiple-instance examples. Machine
Learning, 30:23–29, 1998.
Hanen Borchani, Gherardo Varando, Concha Bielza, and Pedro Larranaga. A survey on
multi-output regression. Wiley Interdisciplinary Reviews: Data Mining and Knowledge
Discovery, 5:216–233, 2015.
Céline Brouard, Florence d’Alché Buc, and Marie Szafranski. Semi-supervised penalized
output kernel regression for link prediction. In International Conference on Machine
Learning (ICML), pages 593–600, 2011.
Andrea Caponnetto. Optimal rates for regularization operators in learning theory. Technical
report, Massachusetts Institute of Technology, 2006. (http://www6.cityu.edu.hk/ma/
doc/people/caponnettoa/regop_TR%28TR11%29.pdf).
Andrea Caponnetto and Ernesto De Vito. Optimal rates for regularized least-squares algorithm. Foundations of Computational Mathematics, 7:331–368, 2007.
34

Learning Theory for Distribution Regression

Claudio Carmeli, Ernesto De Vito, and Alessandro Toigo. Vector valued reproducing kernel
Hilbert spaces of integrable functions and Mercer theorem. Analysis and Applications, 4:
377–408, 2006.
Claudio Carmeli, Ernesto De Vito, Alessandro Toigo, and Veronica Umanitá. Vector valued
reproducing kernel Hilbert spaces and universality. Analysis and Applications, 8:19–61,
2010.
Kevin M. Carter, Raviv Raich, William G. Finn, and Alfred O. Hero. Information-geometric
dimensionality reduction. IEEE Signal Processing Magazine, 28:89–99, 2011.
Jing Chai, Hongtao Chen, Lixia Huang, and Fanhua Shang. Maximum margin multipleinstance feature weighting. Pattern Recognition, 47:2091–2103, 2014a.
Jing Chai, Xinghao Ding, Hongtao Chen, and Tingyu Li. Multiple-instance discriminant
analysis. Pattern Recognition, 47:2517–2531, 2014b.
Ying Chen and Ou Wu. Contextual Hausdorff dissimilarity for multi-instance clustering.
In International Conference on Fuzzy Systems and Knowledge Discovery (FSKD), pages
870–873, 2012.
Andreas Christmann and Ingo Steinwart. Universal kernels on non-standard input spaces.
In Advances in Neural Information Processing Systems (NIPS), pages 406–414, 2010.
Donald L. Cohn. Measure Theory: Second Edition. Birkhäuser Basel, 2013.
Marco Cuturi, Kenji Fukumizu, and Jean-Philippe Vert. Semigroup kernels on measures.
Journal of Machine Learning Research, 6:11691198, 2005.
Ernesto de Vito, Lorenzo Rosasco, and Andrea Caponnetto. Discretization error analysis
for Tikhonov regularization. Analysis and Applications, 4:81–99, 2006.
Thomas G. Dietterich, Richard H. Lathrop, and Tomás Lozano-Pérez. Solving the multiple
instance problem with axis-parallel rectangles. Artificial Intelligence, 89:31–71, 1997.
Jiu Ding and Aihui Zhou. A spectrum theorem for perturbed bounded linear operators.
Applied Mathematics and Computation, 201:723–728, 2008.
Daniel R. Dooly, Qi Zhang, Sally A. Goldman, and Robert A. Amar. Multiple-instance
learning of real-valued data. Journal of Machine Learning Research, 3:651–678, 2002.
Gerald Edgar. Measure, Topology and Fractal Geometry. Springer-Verlag, 1995.
Thomas Eiter and Heikki Mannila. Distance measures for point sets and their computation.
Acta Informatica, 34:109–133, 1997.
Theodoros Evgeniou, Charles A. Micchelli, and Massimiliano Pontil. Learning multiple
tasks with kernel methods. Journal of Machine Learning Research, 6:615–637, 2005.
Seth Flaxman, Yu-Xiang Wang, and Alex Smola. Who supported Obama in 2012? ecological
inference through distribution regression. In ACM SIGKDD International Conference on
Knowledge Discovery and Data Mining (KDD), pages 289–298, 2015.
35

Szabó et al.

James Foulds and Eibe Frank. A review of multi-instance learning assumptions. The
Knowledge Engineering Review, 25:1–25, 2010.
Kenji Fukumizu, Francis Bach, and Michael Jordan. Dimensionality reduction for supervised
learning with reproducing kernel Hilbert spaces. Journal of Machine Learning Research,
5:73–99, 2004.
Thomas Gärtner, Peter A. Flach, Adam Kowalczyk, and Alexander Smola. Multi-instance
kernels. In International Conference on Machine Learning (ICML), pages 179–186, 2002.
Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. Journal of Machine Learning Research, 13:723–773,
2012.
László Györfi, Michael Kohler, Adam Krzyzak, and Harro Walk. A Distribution-Free Theory
of Nonparametric Regression. Springer, New-york, 2002.
David Haussler. Convolution kernels on discrete structures. Technical report, Department
of Computer Science, University of California at Santa Cruz, 1999. (http://cbse.soe.
ucsc.edu/sites/default/files/convolutions.pdf).
Matthias Hein and Olivier Bousquet. Hilbertian metrics and positive definite kernels on
probability measures. In International Conference on Artificial Intelligence and Statistics
(AISTATS), pages 136–143, 2005.
Yang Hu, Mingjing Li, and Nenghai Yu. Multiple-instance ranking: Learning to rank images
for image retrieval. In IEEE Conference on Computer Vision and Pattern Recognition
(CVPR), pages 1–8, 2008.
Tony Jebara, Risi Kondor, and Andrew Howard. Probability product kernels. Journal of
Machine Learning Research, 5:819–844, 2004.
Hachem Kadri, Emmanuel Duflos, Philippe Preux, Stéphane Canu, and Manuel Davy. Nonlinear functional regression: a functional RKHS approach. International Conference on
Artificial Intelligence and Statistics (AISTATS; JMLR W&CP), 9:374–380, 2010.
Hachem Kadri, Alain Rakotomamonjy, Francis Bach, and Philippe Preux. Multiple
operator-valued kernel learning. In Advances in Neural Information Processing Systems
(NIPS), pages 2429–2437, 2012.
Hachem Kadri, Mohammad Ghavamzadeh, and Philippe Preux. A generalized kernel approach to structured output learning. International Conference on Machine Learning
(ICML; JMLR W&CP), 28:471–479, 2013.
Hachem Kadri, Emmanuel Duflos, Philippe Preux, Stéphane Canu, Alain Rakotomamonjy,
and Julien Audiffren. Operator-valued kernels for learning from functional response data.
Journal of Machine Learning Research, 17:1–54, 2016.
Risi Kondor and Tony Jebara. A kernel between sets of vectors. In International Conference
on Machine Learning (ICML), pages 361–368, 2003.
36

Learning Theory for Distribution Regression

Samory Kpotufe. k-NN regression adapts to local intrinsic dimension. Technical report, Max
Planck Institute for Intelligent Systems, 2011. (http://arxiv.org/abs/1110.4300).
James T. Kwok and Pak-Ming Cheung. Marginalized multi-instance kernels. In International Joint Conferences on Artificial Intelligence (IJCAI), pages 901–906, 2007.
Philip M. Long and Lei Tan. PAC learning of axis-aligned rectangles with respect to product
distributions from multiple-instance examples. Machine Learning, 30:7–21, 1998.
David Lopez-Paz, Krikamol Muandet, Bernhard Schölkopf, and Iliya Tolstikhin. Towards a
learning theory of cause-effect inference. International Conference on Machine Learning
(ICML; JMLR W&CP), 37:1452–1461, 2015.
André F. T. Martins, Noah A. Smith, Eric P. Xing, Pedro M. Q. Aguiar, and Mário A. T.
Figueiredo. Nonextensive information theoretical kernels on measures. Journal of Machine Learning Research, 10:935–975, 2009.
Shahar Mendelson and Joseph Neeman. Regularization in kernel learning. The Annals of
Statistics, 38:526–565, 2010.
Charles A. Micchelli and Massimiliano Pontil. On learning vector-valued functions. Neural
Computation, 17:177–204, 2005.
Krikamol Muandet, Kenji Fukumizu, Francesco Dinuzzo, and Bernhard Schölkopf. Learning
from distributions via support measure machines. In Advances in Neural Information
Processing Systems (NIPS), pages 10–18, 2012.
Hans-Georg Müller. Functional modelling and classification of longitudinal data. Scandinavian Journal of Statistics, 32:223–240, 2005.
Frank Nielsen and Richard Nock. A closed-form expression for the Sharma-Mittal entropy
of exponential families. Journal of Physics A: Mathematical and Theoretical, 45:032003,
2012.
Junier Oliva, Barnabás Póczos, and Jeff Schneider. Distribution to distribution regression.
International Conference on Machine Learning (ICML; JMLR W&CP), 28:1049–1057,
2013.
Junier Oliva, William Neiswanger, Barnabás Póczos, Eric Xing, Hy Trac, Shirley Ho, and
Jeff Schneider. Fast function to function regression. International Conference on Artificial
Intelligence and Statistics (AISTATS; JMLR W&CP), 38:717–725, 2015.
Junier B. Oliva, Willie Neiswanger, Barnabás Póczos, Jeff Schneider, and Eric Xing. Fast
distribution to real regression. International Conference on Artificial Intelligence and
Statistics (AISTATS; JMLR W&CP), 33:706–714, 2014.
Kalyanapuram R. Parthasarathy. Probability Measures on Metric Spaces. Academic Press,
1967.
George Pedrick. Theory of reproducing kernels for Hilbert spaces of vector valued functions.
Technical report, 1957.
37

Szabó et al.

Wei Ping, Ye Xu, Kexin Ren, Chi-Hung Chi, and Shen Furao. Non-I.I.D. multi-instance
dimensionality reduction by learning a maximum bag margin subspace. In AAAI Conference on Artificial Intelligence, pages 551–556, 2010.
Barnabás Póczos, Liang Xiong, and Jeff Schneider. Nonparametric divergence estimation
with applications to machine learning on distributions. In Uncertainty in Artificial Intelligence (UAI), pages 599–608, 2011.
Barnabás Póczos, Liang Xiong, Dougal Sutherland, and Jeff Schneider. Support distribution
machines. Technical report, Carnegie Mellon University, 2012. (http://arxiv.org/abs/
1202.0302).
Barnabás Póczos, Alessandro Rinaldo, Aarti Singh, and Larry Wasserman. Distribution-free
distribution regression. International Conference on Artificial Intelligence and Statistics
(AISTATS; JMLR W&CP), 31:507–515, 2013.
Jan Ramon and Maurice Bruynooghe. A polynomial time computable metric between point
sets. Acta Informatica, 37:765–780, 2001.
James O. Ramsay and Bernard W. Silverman. Applied Functional Data Analysis. Springer
Verlag, New York, 2002.
James O. Ramsay and Bernard W. Silverman. Functional Data Analysis. Springer Verlag,
New York, 2005.
Soumya Ray and David Page. Multiple instance regression. In International Conference on
Machine Learning (ICML), pages 425–432, 2001.
Vikas C. Raykar, Balaji Krishnapuram, Jinbo Bi, Murat Dundar, and R. Bharat Rao.
Bayesian multiple instance learning: Automatic feature selection and inductive transfer.
In International Conference on Machine Learning (ICML), pages 808–815, 2008.
Peter Richtárik and Martin Takác̆. Distributed coordinate descent method for learning with
big data. Journal of Machine Learning Research, 17:1–25, 2016.
R. Tyrrell Rockafellar and Roger J-B Wets. Variational Analysis. Springer, 2008.
Alessandro Rudi, Raffaello Camoriano, and Lorenzo Rosasco. Less is more: Nyström computational regularization. In Advances in Neural Information Processing Systems (NIPS),
pages 1648–1656, 2015.
Sivan Sabato and Naftali Tishby. Multi-instance learning with any hypothesis class. Journal
of Machine Learning Research, 13:2999–3039, 2012.
Laurent Schwartz. Analyse III, Calcul Integrél. Hermann, 1998.
Steve Smale and Ding-Xuan Zhou. Estimating the approximation error in learning theory.
Analysis and Applications, 1:17–41, 2003.
38

Learning Theory for Distribution Regression

Bharath Sriperumbudur, Arthur Gretton, Kenji Fukumizu, Gert Lanckriet, and Bernhard
Schölkopf. Hilbert space embeddings and metrics on probability measures. Journal of
Machine Learning Research, 11:1517–1561, 2010.
Bharath K. Sriperumbudur, Kenji Fukumizu, and Gert R. G. Lanckriet. Universality, characteristic kernels and RKHS embedding of measures. Journal of Machine Learning Research, 12:2389–2410, 2011.
Bharath K. Sriperumbudur, Kenji Fukumizu, Revant Kumar, Arthur Gretton, and Aapo
Hyvärinen. Density estimation in infinite dimensional exponential families. Technical
report, 2014. (http://arxiv.org/pdf/1312.3516).
Ingo Steinwart and Andres Christmann. Support Vector Machines. Springer, 2008.
Ingo Steinwart and Clint Scovel. Mercer’s theorem on general domains: On the interaction
between measures, kernels, and RKHSs. Constructive Approximation, 35:363–417, 2012.
Ingo Steinwart, Don R. Hush, and Clint Scovel. Optimal rates for regularized least squares
regression. In Conference on Learning Theory (COLT), 2009.
Hongwei Sun and Qiang Wu. Application of integral operator for regularized least-square
regression. Mathematical and Computer Modelling, 49:276–285, 2009a.
Hongwei Sun and Qiang Wu. A note on application of integral operator in learning theory.
Applied and Computational Harmonic Analysis, 26:416–421, 2009b.
Xu Sun, Hisashi Kashima, and Naonori Ueda. Large-scale personalized human activity
recognition using online multitask learning. IEEE Transactions on Knowledge and Data
Engine, 25:2551–2563, 2013.
Yu-Yin Sun, Michael K. Ng, and Zhi-Hua Zhou. Multi-instance dimensionality reduction.
In AAAI Conference on Artificial Intelligence, pages 587–592, 2010.
Dougal J. Sutherland, Junier B. Oliva, Barnabás Póczos, and Jeff Schneider. Linear-time
learning on distributions with approximate kernel embeddings. In AAAI Conference on
Artifical Intelligence (AAAI), pages 2073–2079, 2016.
Zoltán Szabó, Arthur Gretton, Barnabás Póczos, and Bharath Sriperumbudur. Two-stage
sampled learning theory on distributions. In International Conference on Artificial Intelligence and Statistics (AISTATS), pages 948–957, 2015.
Leslie Valiant. A theory of the learnable. Communications of the ACM, 27:1134–1142, 1984.
Fei Wang, Tanveer Syeda-Mahmood, Baba C. Vemuri, David Beymer, and Anand Rangarajan. Closed-form Jensen-Rényi divergence for mixture of Gaussians and applications to
group-wise shape registration. Medical Image Computing and Computer-Assisted Intervention, 12:648–655, 2009.
Jun Wang and Jean-Daniel Zucker. Solving the multiple-instance problem: A lazy learning
approach. In International Conference on Machine Learning (ICML), pages 1119–1126,
2000.
39

Szabó et al.

Larry Wasserman. All of Nonparametric Statistics. Springer, 2006.
Ou Wu, Jun Gao, Weiming Hu, Bing Li, and Mingliang Zhu. Identifying multi-instance
outliers. In SIAM International Conference on Data Mining (SDM), pages 430–441, 2010.
Yun Yang, Mert Pilanci, and Martin J. Wainwright. Randomized sketches for kernels: Fast
and optimal non-parametric regression. Annals of Statistics, 2016. (to appear; arXiv:
http://arxiv.org/abs/1501.06195).
Amelia Zafra, Mykola Pechenizkiy, and Sebastián Ventura. HyDR-MI: A hybrid algorithm
to reduce dimensionality in multiple instance learning. Information Sciences, 222:282–
301, 2013.
Dan Zhang and Luo Si. Multiple instance transfer learing. In IEEE International Conference
on Data Mining Workshops (ICDMW), pages 406–411, 2009.
Dan Zhang, Fei Wang, Luo Si, and Tao Li. M3IC: Maximum margin multiple instance
clustering. In International Joint Conferences on Artificial Intelligence (IJCAI), pages
1339–1344, 2009.
Dan Zhang, Fei Wang, Luo Si, and Tao Li. Maximum margin multiple instance clustering
with applications to image and text clustering. IEEE Transactions on Neural Networks,
22:739–751, 2011.
Dan Zhang, Jingrui He, Luo Si, and Richard D. Lawrence. MILEAGE: Multiple Instance LEArning with Global Embedding. International Conference on Machine Learning
(ICML; JMLR W&CP), 28:82–90, 2013.
Min-Ling Zhang and Zhi-Hua Zhou. Multi-instance clustering with applications to multiinstance prediction. Applied Intelligence, 31:47–68, 2009.
Yuchen Zhang, John C. Duchi, and Martin J. Wainwright. Divide and conquer kernel ridge
regression: A distributed algorithm with minimax optimal rates. 16:3299–3340, 2015.
Jiayu Zhou, Jun Liu, Vaibhav A. Narayan, and Jieping Ye. Modeling disease progression
via multi-task learning. NeuroImage, 78:233–248, 2013.
Shaohua Kevin Zhou and Rama Chellappa. From sample similarity to ensemble similarity:
Probabilistic distance measures in reproducing kernel Hilbert space. IEEE Transactions
on Pattern Analysis and Machine Intelligence, 28:917–929, 2006.
Zhi-Hua Zhou. Multi-instance learning: A survey. Technical report, AI Lab, Department
of Computer Science & Technology, Nanjing University, Nanjing, China, 2004. (http://
cs.nju.edu.cn/zhouzh/zhouzh.files/publication/techrep04.pdf).
Zhi-Hua Zhou, Yu-Yin Sun, and Yu-Feng Li. Multi-instance learning by treating instances
as non-I.I.D. samples. In International Conference on Machine Learning (ICML), pages
1249–1256, 2009.

40

Learning Theory for Distribution Regression

Supplement
2
H

9. Proofs of the fzλ

Bound, Theorem 5 and Theorem 9
2

This section contains the derivations of the fzλ H bound (used in Theorem 2; see Section 9.1), Theorem 5 (Section 9.2) and Theorem 9 (Section 9.3).
9.1 Bound on fzλ

2
H
2

Below we derive the stated Eq. (23) bound for fzλ H ; it is guaranteed to hold under
the conditions of the bounds for S1 and S2 obtained in (Caponnetto and De Vito, 2007,
Eq. (48), above Eq. (46), Eq. (43)); see (†)1 , (†)2 , (†)3 below.
Applying the triangle inequality and the definition of B(λ) we get
fzλ

H

≤ fzλ − f λ

H

+ f λ − fρ

H

+ kfρ kH = fzλ − f λ

H

+

p

B(λ) + kfρ kH .

fzλ − f λ can be decomposed (Caponnetto and De Vito, 2007, page 347) as
fzλ − f λ = (Tx + λI)−1 (gz − Tx fρ ) + (Tx + λI)−1 (T − Tx )(f λ − fρ ).
Thus
fzλ − f λ

H

≤ (Tx + λI)−1 (gz − Tx fρ ) H + (Tx + λI)−1 (T − Tx )(f λ − fρ ) .
H
{z
} |
|
{z
}
=:(∗)1

=:(∗)2

1. The first term can be estimated as

1

1

(T + λI)− 2 (gz − Tx fρ )
L(H)
H
#
" r
r


(†)1
1
1 M 2 BK
Σ2 N(λ)
6
−1
≤ (Tx + λI) (T + λI) 2
.
+
2 log
η
l
λ
l
L(H)
|
{z
}

(∗)1 ≤ (Tx + λI)−1 (T + λI) 2

=:(∗)3

By (Caponnetto and De Vito, 2007, page 350)

h
i
1
1 −1
1
(∗)3 = (T + λI)− 2 I − (T + λI)− 2 (T − Tx )(T + λI)− 2
1

≤ (T + λI)− 2

L(H)

h

1

L(H)
1

I − (T + λI)− 2 (T − Tx )(T + λI)− 2

i−1

(†)2
L(H)

2
≤ √ .
λ

2. The second expression can be bounded as
(∗)2 ≤ (Tx + λI)−1

L(H)

(T − Tx )(f λ − fρ )
41

(†)3
H

2
log
≤
λ

#
r
p
 "
2BK B(λ)
6
BK A(λ)
.
+
η
l
l

Szabó et al.

To sum up, we obtained that
fzλ

H

#
r
 " r 2
6
1 M BK
Σ2 N(λ)
2
+
≤ √ 2 log
η
l
λ
l
λ
#
"
r
p
 
p
2BK B(λ)
2
6
BK A(λ)
+ log
+ B(λ) + kfρ kH ,
+
λ
η
l
l

and hence by Eq. (35)

  2

2
M BK
Σ2 N(λ)
6
16
log2
+
fzλ
≤6
λ
η
l2 λ
l
H
  2


4BK B(λ) BK A(λ)
4
2
2 6
+
+ B(λ) + kfρ kH .
+ 2 log
λ
η
l2
l
9.2 Proof of Theorem 5
In the following λ is chosen to match the ’bias’ (λc ) and a ’variance’ (other) term in r(l, λ)
[see Eq. (48)], guarantee that the matched terms dominate and the constraints in Eq. (48)
are also satisfied; according to our assumptions 1 < b and c ∈ (1, 2].
r(l, λ) =

1
l2+a λ3

+

1
la λ

+

1
la+1 λ

+ λc +

2+ 1b

1
l2 λ

+

1
lλ

1
b

→ 0, s.t. lλ

b+1
b

≥ 1, la λ2 ≥ 1. (48)

Casesin terms
i of a choice:
b(c+1)
• a ∈ 0, bc+1 : Since b(c+1)
bc+1 < 2 ⇔ b(1 − c) < 2 (⇐ b > 1, c > 1) one has
r(l, λ) =

1

l2+a λ3

+

1

la λ

+

1

2+ 1b

la+1 λ

1

+ λc +

lλ

1
b

→ 0, s.t. lλ

b+1
b

1
l2 λ

≤

1
la λ ,

≥ 1, la λ2 ≥ 1.

and
(49)

– ( 1 =) 4 : [order of the 1st and 4th terms in Eq. (49) are the matched and the first
term will be discarded]:

1
l2+a λ3

2+a

= λc ⇔ λ = l− 3+c , and
2+a

2+a

1

r(l) = l−a+ 3+c + l−(a+1)+ 3+c (2+ b ) + l−
s.t. l

b+1
1− 2+a
3+c b

≥ 1, l

a−2 2+a
3+c

✄

✄

c+1+ 2b
1− 1b

> 2, which contradicts to a ≤

(2+a)c
3+c

+l

2+a
−1+ (3+c)b

→ 0,

(50)

≥ 1.

≥ −(a + 1) +
∗ ✂3 ✁ ≥ ✂2 ✁ [the 3rd term dominates the 2nd in Eq. (50)]: − (2+a)c

3+c
1
1
1
2+a
2
+
⇔
(a
+
1)(c
+
3)
≥
(2
+
a)(c
+
2
+
)
⇔
a
1
−
≥
c
+
1 + 2b ⇔ a ≥
3+c
b
b
b
– ( 2 =) 4 :

1
la λ

= λc ⇔ λ = l

a
− c+1

b(c+1)
bc+1

< 2.

, and

3a

1

r(l) = l−(2+a)+ c+1 + l−(a+1)+ c+1 (2+ b ) + l− c+1 + l
s.t. l

a(b+1)
1− b(c+1)

a

≥ 1, l

2a
a− c+1

42

≥ 1.

ac

a
−1+ b(c+1)

→ 0,

Learning Theory for Distribution Regression

✄

✄

✄

✄

ac
3a
≥ −(2 + a) + c+1
⇔ (2 + a)(c + 1) ≥ a(c + 3) ⇔ c + 1 ≥ a; this
∗ ✂3 ✁ ≥ ✂1 ✁: − c+1
property
holds
because
c
+
1
>
2
≥
a. 
✄
✄

ac
a
∗ ✂3 ✁≥ ✂2 ✁: − c+1 ≥ −(a + 1) + c+1 2 + 1b ⇔ (a + 1)(c + 1) ≥ a c + 2 + 1b ⇔ c + 1 ≥

c+1
c+1
> b(c+1)
a 1 + 1b ⇔ 1+
1 ≥ a; this holds since
bc+1 ≥ a.
1+ 1
b

b

b(c+1)
bc+1

a
ac
≥ −1 + b(c+1)
⇔ b ≥ a. Since
∗ ✂3 ✁ ≥ ✂4 ✁: − c+1
b(c+1)
property holds by our a ≤ bc+1 assumption.

≤ b ⇔ 1 ≤ b, the required

∗ Constraint-1 in r(l): It is sufficient that 1 − a(b+1)
b(c+1) > 0 ⇔ a <
b(c+1)
bc+1

b(c+1)
b+1 ;

this holds since

b(c+1)
b+1 .

a≤
<
2a
∗ Constraint-2 in r(l): It is enough to have a − c+1
> 0 ⇔ c + 1 > 2, which holds since
c > 1.
ac
a
To sum up, in this case λ = l− c+1 and the rate is r(l) = l− c+1 → 0.
– ( 3 =) 4 :

1

1
la+1 λ2+ b

r(l) = l

= λc ⇔ λ = l

−(2+a)+ 3(a+1)
1
2+

s.t. l
✄

1−

✄

b

+c

−

+l

a+1 b+1
2+ 1 +c b
b

a+1
2+ 1 +c
b

−a+

, and

a+1
2+ 1 +c
b

≥ 1, l

+ l

a− 2(a+1)
1
2+

b

+c

− c(a+1)
1
2+

a≤

b(c+1)
bc+1

– ( 5 =) 4 :

−1+

a+1 1
2+ 1 +c b
b

→ 0,



assumption.

λc

=

1

1

lλ b

⇔λ=l

−(2+a)+

s.t. l

−

3
c+ 1
b

1
c+ 1
b

+l

1− 1 1 b+1
b
c+
b

, and

−a+

1
c+ 1
b

≥ 1, l

+l

a− 2 1
c+
b

✄

∗ ✂4 ✁≥ ✂2 ✁: − c+c 1 ≥ −a + c+1 1 ⇔ a ≥
b
b
 to our a choice.
contradicts
• a ∈ b(c+1)
bc+1 , ∞ :
r(l, λ) =

+l

a+1
⇔ a 2 + 1b + c ≥ (a + 1)(c + 1) ⇔ a(1 + 1b ) ≥
2+ 1b +c
b(c+1)
c+1
However, 1+
1 >
bc+1 (⇔ c > 1) which contradicts to our
b

c+1
.
1+ 1b

r(l) = l

✄

+c

≥ 1.

c(a+1)
≥ −a +
∗ ✂3 ✁ ≥ ✂2 ✁: − 2+
1
+c
b

c+1 ⇔ a ≥

b

c+1
.
c+ 1b

−(a+1)+

1
c+ 1
b

(2+ 1b )

+ l

−

c
c+ 1
b

→ 0,

≥ 1.
✄

✄

As we have seen ( 3 = 4 : ✂3 ✁≥ ✂2 ✁) this

b+1
1
1
1
1
1
+ λc + 2 + 1 → 0, s.t. lλ b ≥ 1, la λ2 ≥ 1.
+
+
l2+a λ3 la λ la+1 λ2+ 1b
l λ lλ b

– ( 6 =) 4 : Using the previous ’ 5 = 4 ’ case, λc =

r(l) = l

−(2+a)+

s.t. l

3
c+ 1
b

+l

1− 1 1 b+1
b
c+
b

−a+

1
c+ 1
b

≥ 1, l

+l

a− 2 1
c+
b

−(a+1)+

≥ 1.
43

1
c+ 1
b

1

1

lλ b

(2+ 1b )

⇔λ=l
+ l

−

c
c+ 1
b

−

1
c+ 1
b

+l

, and

−2+

1
c+ 1
b

→ 0,

Szabó et al.

✄


b(c+3)
c+3
1
b ≥ c+3 ⇔ a ≥ c+ 1 −2 = bc+1 −2.
b
b(c+1)
b(c+3)
bc+1 ≥ bc+1 − 2 ⇔ bc + 1 ≥ b which

✄

∗ ✂4 ✁≥ ✂1 ✁: − c+c 1 ≥ −(2+a)+ c+3 1 ⇔ (2+a) c +
b
b

This requirement holds by our a choice since
✄
✄is valid.
∗ ✂4 ✁≥ ✂2 ✁: As we have seen (previous ’ 5 = 4 ’) this means a ≥ b(c+1)
bc+1 which holds
by
our
a
choice.
✄
✄


∗ ✂4 ✁ ≥ ✂3 ✁: − c+c 1 ≥ −(a + 1) + c+1 1 2 + 1b ⇔ (a + 1) c + 1b ≥ 2 + 1b + c ⇔ a ≥
2+ 1b +c
−1
c+ 1b
b(c+1)
⇔b
✄ bc+1 ✄

b

=

b

2b+1+bc
bc+1

− 1. This condition holds by our a choice since

2b+1+bc
bc+1

−1≤

≤ bc which is valid.


∗ ✂4 ✁≥ ✂5 ✁: − c+c 1 ≥ −2 + c+1 1 ⇔ 2 c + 1b ≥ c + 1 ⇔ 2 ≥ b(c+1)
bc+1 what we have already
b
b
established.
b+1
∗ Constraint-1 in r(l): It is enough to have 1 − c+1 1 b+1
b > 0 ⇔ 1 > bc+1 ⇔ bc > b,
b

which holds.
∗ Constraint-2 in r(l): It is sufficient that a −
satisfied because a >

2
c+ 1b

b(c+1)
2b
bc+1 > bc+1 .
1
− 1
b
− bc+1
c+
b

2
c+ 1b

> 0 ⇔a >

=

2b
bc+1 ;

this is

bc

, the rate is r(l) = l− bc+1 .
To sum up, in this case λ = l
=l
✄
✄
2+a
– ( 1 =) 4 : Using the previous ’✂1 ✁= ✂4 ✁’ case, l2+a1 λ3 = λc ⇔ λ = l− 3+c , and
2+a

2+a

1

r(l) = l−a+ 3+c + l−(a+1)+ 3+c (2+ b ) + l−
s.t. l

b+1
1− 2+a
3+c b

≥ 1, l

a−2 2+a
3+c

✄
✄
∗ ✂3 ✁ ≥ ✂5 ✁: − (2+a)c
3+c ≥ −1 +

– (

(2+a)c
3+c

2+a

+ l−2+ 3+c + l

2+a
−1+ (3+c)b

→ 0,

≥ 1.

(3+c)b
2+a
(3+c)b ⇔ (3 + c)b ≥ (2 + a)(1 + bc) ⇔ 1+bc − 2 ≥ a.
b(c+1)
However, (3+c)b
1+bc − 2 ≤ bc+1 ⇔ b ≤ bc + 1 holds, thus the requirement can not be
satisfied due to the b(c+1)
bc+1 < a assumption.
a
2 =) 4 : Using the previous ’ 2 = 4 ’ case we have la1λ = λc ⇔ λ = l− c+1 , and
3a

1

r(l) = l−(2+a)+ c+1 + l−(a+1)+ c+1 (2+ b ) + l− c+1 + l−2− c+1 + l
s.t. l
✄

a(b+1)
1− b(c+1)

✄

a

ac

a

2a

≥ 1, la− c+1 ≥ 1.

ac
a
⇔ b(c + 1) ≥ a(1 + bc) ⇔
≥ −1 + b(c+1)
∗ ✂3 ✁≥ ✂5 ✁: − c+1
b(c+1)
to our bc+1 < a choice.

– ( 3 =) 4 : Using the previous ’ 3 = 4 ’,
r(l) = l

−(2+a)+

s.t. l

3(a+1)
2+ 1 +c
b

a
−1+ b(c+1)

+l

b+1
1− a+1
2+ 1 +c b
b

−a+

a+1
2+ 1 +c
b

≥ 1, l

✄
✄
c(a+1)
≥ −1 +
∗ ✂3 ✁≥ ✂5 ✁: − 2+
1
+c
b

1

la+1 λ2+ b

+ l

a− 2(a+1)
2+ 1 +c
b

a+1 1
2+ 1b +c b

1

−

c(a+1)
2+ 1 +c
b

b(c+1)
bc+1

= λc ⇔ λ = l
+l

−2+

a+1
2+ 1 +c
b

−

→ 0,

≥ a; this contradicts
a+1
2+ 1 +c
b

+l

−1+

, and

a+1 1
2+ 1 +c b
b

→ 0,

≥ 1.
2b
1+bc − 1 ≥ a.
b(c+1)
bc+1 assumption

⇔ 2b + 1 + bc ≥ (a + 1)(1 + bc) ⇔

2b
However, b(c+1)
bc+1 ≥ 1+bc − 1 ⇔ 2bc + 1 ≥ b holds; thus, by the a >
the required property can not be satisfied.

44

Learning Theory for Distribution Regression

1
l2 λ

– ( 5 =) 4 : λc =

2

⇔ λ = l− c+1 , and
2

6

2

2c

1

r(l) = l−(a+2)+ c+1 + l−a+ c+1 + l−(a+1)+ c+1 (2+ b ) + l− c+1 + l
s.t. l

2(b+1)
1− b(c+1)

✄

✄

2
−1+ b(c+1)

4

→ 0,

≥ 1, la− c+1 ≥ 1.

2c
≥ −1 +
∗ ✂4 ✁≥ ✂5 ✁: − c+1

2
b(c+1)

2+2bc
b(c+1)

⇔1≥

⇔ b ≥ bc + 2 which does not hold.

9.3 Proof of Theorem 9
p
In the sequel we choose λ by matching 2 terms in r(l, λ) [Eq. (51)], guarantee that the
matched terms dominate and the constraint in Eq. (51) holds; we proceed by matching the
’bias’ (λs ) and ’variance’ (other) terms; s ∈ (0, 1].
p

• ( 1 =) 3 :

✄

1
la λ

r(l, λ) =

1
la λ

+

1
1− 2s

λ

l

+ λs → 0, s.t. lλ2 ≥ 1.

1
2

(51)

a

= λs ⇔ λ = l− s+1 , and
p

✄

1

a

s

2a

sa

r(l) = l− 2 + s+1 (1− 2 ) + l− s+1 → 0, s.t. l1− s+1 ≥ 1.

a
sa
≥ − 12 + s+1
(1 − 2s ) ⇔ s + 1 ≥ a(2 + s) ⇔ s+1
– ✂2 ✁≥ ✂1 ✁: − s+1
s+2 ≥ a.
2a
– Condition in r(l): it is sufficient to have 1 − s+1
> 0 ⇔ s+1
2 > a.
p

a
sa
− s+1
s+1
s+1
To sum up, if a ≤ s+2 < 2 then λ = l
leads to the rate r(l) = l− s+1 ; specifip
s
r(l) = l− s+2 .
cally, if a = s+1
then
s+2

• ( 2 =) 3 :

✄

1

s 1
λ1− 2 l 2

=

λs

p

✄

⇔λ=l

−1
2
s
1+ 2

1
=− s+2

, and

1

1
s+1
s
≥ −a + s+2
⇔ a ≥ s+2
.
– ✂2 ✁≥ ✂1 ✁: − s+2
– Condition in r(l): it is sufficient to have 1 −

To sum up, if

s+1
s+2

2

s

r(l) = l−a+ s+2 + l− s+2 → 0, s.t. l1− s+2 ≥ 1.

≤ a then choosing λ = l

1
− s+2

45

2
s+2

> 0 ⇔ s > 0 which always holds.
p
s
the rate is r(l) = l− s+2 → 0.

