# density_estimation_in_infinite_dimensional_exponential_families__189f7129.pdf

Journal of Machine Learning Research 18 (2017) 1-59 Submitted 1/16; Revised 4/17; Published 7/17

Density Estimation in Infinite Dimensional Exponential Families

Bharath Sriperumbudur bks18@psu.edu Department of Statistics, Pennsylvania State University University Park, PA 16802, USA.

Kenji Fukumizu fukumizu@ism.ac.jp The Institute of Statistical Mathematics 10-3 Midoricho, Tachikawa, Tokyo 190-8562 Japan.

Arthur Gretton arthur.gretton@gmail.com ORCID 0000-0003-3169-7624 Gatsby Computational Neuroscience Unit, University College London Sainsbury Wellcome Centre, 25 Howland Street, London W1T 4JG, UK

Aapo Hyv arinen a.hyvarinen@ucl.ac.uk Gatsby Computational Neuroscience Unit, University College London Sainsbury Wellcome Centre, 25 Howland Street, London W1T 4JG, UK

Revant Kumar rkumar74@gatech.edu College of Computing, Georgia Institute of Technology 801 Atlantic Drive, Atlanta, GA 30332, USA.

Editor: Ingo Steinwart

In this paper, we consider an infinite dimensional exponential family P of probability densities, which are parametrized by functions in a reproducing kernel Hilbert space H, and show it to be quite rich in the sense that a broad class of densities on Rd can be approximated arbitrarily well in Kullback-Leibler (KL) divergence by elements in P. Motivated by this approximation property, the paper addresses the question of estimating an unknown density p0 through an element in P. Standard techniques like maximum likelihood estimation (MLE) or pseudo MLE (based on the method of sieves), which are based on minimizing the KL divergence between p0 and P, do not yield practically useful estimators because of their inability to efficiently handle the log-partition function. We propose an estimator ˆpn based on minimizing the Fisher divergence, J(p0 p) between p0 and p P, which involves solving a simple finite-dimensional linear system. When p0 P, we show that the proposed estimator is consistent, and provide a convergence rate of n min{ 2

2β+2} in Fisher divergence under the smoothness assumption that log p0 R(Cβ) for some β 0, where C is a certain Hilbert-Schmidt operator on H and R(Cβ) denotes the image of Cβ. We also investigate the misspecified case of p0 / P and show that J(p0 ˆpn) infp P J(p0 p) as n , and provide a rate for this convergence under a similar smoothness condition as above. Through numerical simulations we demonstrate that the proposed estimator outperforms the non-parametric kernel density estimator, and that the advantage of the proposed estimator grows as d increases.

c 2017 Bharath Sriperumbudur, Kenji Fukumizu, Arthur Gretton, Aapo Hyv arinen and Revant Kumar.

License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/. Attribution requirements are provided at http://jmlr.org/papers/v18/16-011.html.

Sriperumbudur et al.

Keywords: density estimation, exponential family, Fisher divergence, kernel density estimator, maximum likelihood, interpolation space, inverse problem, reproducing kernel Hilbert space, Tikhonov regularization, score matching

1. Introduction

Exponential families are among the most important classes of parametric models studied in statistics, and include many common distributions such as the normal, exponential, gamma, and Poisson. In its natural form , the family generated by a probability density q0 (defined over Ω Rd) and sufficient statistic, T : Ω Rm is defined as

Pfin := n pθ(x) = q0(x)eθT T(x) A(θ), x Ω: θ Θ Rmo (1)

where A(θ) := log R

ΩeθT T(x)q0(x) dx is the cumulant generating function (also called the log-partition function), Θ {θ Rm : A(θ) < } is the natural parameter space and θ is a finite-dimensional vector called the natural parameter. Exponential families have a number of properties that make them extremely useful for statistical analysis (see Brown, 1986 for more details).

In this paper, we consider an infinite dimensional generalization (Canu and Smola, 2005; Fukumizu, 2009) of (1),

P = n pf(x) = ef(x) A(f)q0(x), x Ω: f F o ,

where the function space F is defined as

F = n f H : e A(f) < o , with A(f) := log Z

Ω ef(x)q0(x) dx

being the cumulant generating function, and (H, , H) a reproducing kernel Hilbert space (RKHS) (Aronszajn, 1950) with k as its reproducing kernel. While various generalizations are possible for different choices of F (e.g., an Orlicz space as in Pistone and Sempi, 1995), the connection of P to the natural exponential family in (1) is particularly enlightening when H is an RKHS. This is due to the reproducing property of the kernel, f(x) = f, k(x, ) H, through which k(x, ) takes the role of the sufficient statistic. In fact, it can be shown (see Section 3 and Example 1 for more details) that every Pfin is generated by P induced by a finite dimensional RKHS H, and therefore the family P with H being an infinite dimensional RKHS is a natural infinite dimensional generalization of Pfin. Furthermore, this generalization is particularly interesting as in contrast to Pfin, it can be shown that P is a rich class of densities (depending on the choice of k and therefore H) that can approximate a broad class of probability densities arbitrarily well (see Propositions 1, 13 and Corollary 2). This generalization is not only of theoretical interest, but also has implications for statistical and machine learning applications. For example, in Bayesian non-parametric density estimation, the densities in P are chosen as prior distributions on a collection of probability densities (e.g., see van der Vaart and van Zanten, 2008). P has also found applications in nonparametric hypothesis testing (Gretton et al., 2012; Fukumizu et al., 2008) and dimensionality reduction (Fukumizu et al., 2004, 2009) through the mean and

Density Estimation in Infinite Dimensional Exponential Families

covariance operators, which are obtained as the first and second Fr echet derivatives of A(f) (see Fukumizu, 2009, Section 1.2.3). Recently, the infinite dimensional exponential family, P has been used to develop a gradient-free adaptive MCMC algorithm based on Hamiltonian Monte Carlo (Strathmann et al., 2015) and also has been used in the context of learning the structure of graphical models (Sun et al., 2015).

Motivated by the richness of the infinite dimensional generalization and its statistical applications, it is of interest to model densities by P, and therefore the goal of this paper is to estimate unknown densities by elements in P when H is an infinite dimensional RKHS. Formally, given i.i.d. random samples (Xa)n a=1 drawn from an unknown density p0, the goal is to estimate p0 through P. Throughout the paper, we refer to case of p0 P as wellspecified, in contrast to the misspecified case where p0 / P. The setting is useful because P is a rich class of densities that can approximate a broad class of probability densities arbitrarily well, hence it may be widely used in place of non-parametric density estimation methods (e.g., kernel density estimation (KDE)). In fact, through numerical simulations, we show in Section 6 that estimating p0 through P performs better than KDE, and that the advantage of the proposed estimator grows with increasing dimensionality.

In the finite-dimensional case where θ Θ Rm, estimating pθ through maximum likelihood (ML) leads to solving elegant likelihood equations (Brown, 1986, Chapter 5). However, in the infinite dimensional case (assuming p0 P), as in many non-parametric estimation methods, a straightforward extension of maximum likelihood estimation (MLE) suffers from the problem of ill-posedness (Fukumizu, 2009, Section 1.3.1). To address this problem, Fukumizu (2009) proposed a method of sieves involving pseudo-MLE by restricting the infinite dimensional manifold P to a series of finite-dimensional submanifolds, which enlarge as the sample size increases, i.e., p ˆf(l) is the density estimator with

ˆf(l) = arg max f F(l) 1 n

a=1 f(Xa) A(f), (2)

where F(l) = {f H(l) : e A(f) < } and (H(l)) l=1 is a sequence of finite-dimensional subspaces of H such that H(l) H(l+1) for all l N. While the consistency of p ˆf(l) is proved in Kullback-Leibler (KL) divergence (Fukumizu, 2009, Theorem 6), the method suffers from many drawbacks that are both theoretical and computational in nature. On the theoretical front, the consistency in Fukumizu (2009, Theorem 6) is established by assuming a decay rate on the eigenvalues of the covariance operator (see (A-2) and the discussion in Section 1.4 of Fukumizu (2009) for details), which is usually difficult to check in practice. Moreover, it is not clear which classes of RKHS should be used to obtain a consistent estimator (Fukumizu, 2009, (A-1)) and the paper does not provide any discussion about the convergence rates. On the practical side, the estimator is not attractive as it can be quite difficult to construct the sequence (H(l)) l=1 that satisfies the assumptions in Fukumizu (2009, Theorem 6). In fact, the impracticality of the estimator, ˆf(l) is accentuated by the difficulty in efficiently handling A(f) (though it can be approximated by numerical integration).

A related work was carried out by Barron and Sheu (1991) also see references therein where the goal is to estimate a density, p0 by approximating its logarithm as an expansion in terms of basis functions, such as polynomials, splines or trigonometric series. Similar to

Sriperumbudur et al.

Fukumizu (2009), Barron and Sheu proposed the ML estimator p ˆfm, where

ˆfm = arg max f Fm 1 n

a=1 f(Xa) A(f)

and Fm is the linear space of dimension m spanned by the chosen basis functions. Under the assumption that log p0 has square-integrable derivatives up to order r, they showed that KL(p0 p ˆfm) = Op0(n 2r/(2r+1)) with m = n1/(2r+1) for each of the approximating families, where KL(p q) = R p(x) log(p(x)/q(x)) dx is the KL divergence between p and q. Similar work was carried out by Gu and Qiu (1993), who assumed that log p0 lies in an RKHS, and proposed an estimator based on penalized MLE, with consistency and rates established in Jensen-Shannon divergence. Though these results are theoretically interesting, these estimators are obtained via a procedure similar to that in Fukumizu (2009), and therefore suffers from the practical drawbacks discussed above.

The discussion so far shows that the MLE approach to learning p0 P results in estimators that are of limited practical interest. To alleviate this, one can treat the problem of estimating p0 P in a completely non-parametric fashion by using KDE, which is well-studied (Tsybakov, 2009, Chapter 1) and easy to implement. This approach ignores the structure of P, however, and is known to perform poorly for moderate to large d (Wasserman, 2006, Section 6.5) (see also Section 6 of this paper).

1.1 Score Matching and Fisher Divergence

To counter the disadvantages of KDE and pseudo/penalized-MLE, in this paper, we propose to use the score matching method introduced by Hyv arinen (2005, 2007). While MLE is based on minimizing the KL divergence, the score matching method involves minimizing the Fisher divergence (also called the Fisher information distance; see Definition 1.13 in Johnson (2004)) between two continuously differentiable densities, p and q on an open set Ω Rd, given as

Ω p(x) log p(x) log q(x) 2 2 dx, (3)

where log p(x) = ( 1 log p(x), . . . , d log p(x)) with i log p(x) := xi log p(x). Fisher divergence is closely related to the KL divergence through de Bruijn s identity (Johnson, 2004, Appendix C) and it can be shown that KL(p q) = R 0 J(pt qt) dt, where pt = p N(0, t Id), qt = q N(0, t Id), denotes the convolution, and N(0, t Id) denotes a normal distribution on Rd with mean zero and diagonal covariance with t > 0 (see Proposition B.1 for a precise statement; also see Theorem 1 in Lyu, 2009). Moreover, convergence in Fisher divergence is a stronger form of convergence than that in KL, total variation and Hellinger distances (see Lemmas E.2 & E.3 in Johnson, 2004 and Corollary 5.1 in Ley and Swan, 2013).

To understand the advantages associated with the score matching method, let us consider the problem of density estimation where the data generating distribution (say p0) belongs to Pfin in (1). In other words, given random samples (Xa)n a=1 drawn i.i.d. from p0 := pθ0, the goal is to estimate θ0 as ˆθn, and use pˆθn as an estimator of p0. While the

Density Estimation in Infinite Dimensional Exponential Families

MLE approach is well-studied and enjoys nice statistical properties in asymptopia (i.e., asymptotically unbiased, efficient, and normally distributed), the computation of ˆθn can be intractable in many situations as discussed above. In particular, this is the case for pθ(x) = rθ(x)

A(θ) where rθ 0 for all θ Θ, A(θ) = R

Ωrθ(x) dx, and the functional form of r is known (as a function of θ and x); yet we do not know how to easily compute A, which is often analytically intractable. In this setting (which is exactly the setting of this paper), assuming pθ to be differentiable (w.r.t. x), and R

Ωp0(x) log pθ(x) 2 2 dx < , θ Θ, J(p0 pθ) =: J(θ) in (3) reduces to

2 ( i log pθ(x))2 + 2 i log pθ(x) dx + 1

Ω p0(x) log p0(x) 2 2 dx,

(4) through integration by parts (see Hyv arinen, 2005, Theorem 1), under appropriate regularity conditions on p0 and pθ for all θ Θ. Here 2 i log pθ(x) := 2 x2 i log pθ(x). The main advantage of the objective in (3) (and also (4)) is that when it is applied to the situation discussed above where pθ(x) = rθ(x)

A(θ) , J(θ) is independent of A(θ), and an estimate of θ0 can be obtained by simply minimizing the empirical counterpart of J(θ), given by

2 ( i log pθ(Xa))2 + 2 i log pθ(Xa) + 1

Ω p0(x) log p0(x) 2 2 dx.

Since Jn(θ) is also independent of A(θ), ˆθn = arg minθ Θ Jn(θ) may be easily computable, unlike the MLE. We would like to highlight that while the score matching approach may have computational advantages over MLE, it only estimates pθ up to the scaling factor A(θ), and therefore requires the approximation or computation of A(θ) through numerical integration to estimate pθ. Note that this issue (of computing A(θ) through numerical integration) exists even with MLE, but not with KDE. In score matching, however, numerical integration is needed only once, while MLE would typically require a functional form of the log-partition function which is approximated through numerical integration at every step of an iterative optimization algorithm (for example, see (2)), thus leading to major computational savings. An important application that does not require the computation of A(θ) is in finding modes of the distribution, which has recently become very popular in image processing (Comaniciu and Meer, 2002), and has already been investigated in the score matching framework (Sasaki et al., 2014). Similarly, in sampling methods such as sequential Monte Carlo (Doucet et al., 2001), it is often the case that the evaluation of unnormalized densities is sufficient to calculate required importance weights.

1.2 Contributions

(i) We present an estimate of p0 P in the well-specified case through the minimization of Fisher divergence, in Section 4. First, we show that estimating p0 := pf0 using the score matching method reduces to estimating f0 by solving a simple finite-dimensional linear system (Theorems 4 and 5). Hyv arinen (2007) obtained a similar result for Pfin where the estimator is obtained by solving a linear system, which in the case of Gaussian family

Sriperumbudur et al.

matches the MLE (Hyv arinen, 2005). The estimator obtained in the infinite dimensional case is not a simple extension of its finite-dimensional counterpart, however, as the former requires an appropriate regularizer (we use 2 H) to make the problem well-posed. We would like to highlight that to the best of our knowledge, the proposed estimator is the first practically computable estimator of p0 with consistency guarantees (see below).

(ii) In contrast to Hyv arinen (2007) where no guarantees on consistency or convergence rates are provided for the density estimator in Pfin, we establish in Theorem 6 the consistency and rates of convergence for the proposed estimator of f0, and use these to prove consistency and rates of convergence for the corresponding plug-in estimator of p0 (Theorems 7 and B.2), even when H is infinite dimensional. Furthermore, while the estimator of f0 (and therefore p0) is obtained by minimizing the Fisher divergence, the resultant density estimator is also shown to be consistent in KL divergence (and therefore in Hellinger and total-variation distances) and we provide convergence rates in all these distances.

Formally, we show that the proposed estimator ˆfn is converges as

f0 ˆfn H = Op0(n α), KL(p0 p ˆfn) = Op0(n 2α) and J(p0 p ˆfn) = Op0

n min n 2 3 , 2β+1

if f0 R(Cβ) for some β > 0, where R(A) denotes the range or image of an operator A, α = min{1

4, β 2β+2}, and C := Pd i=1 R

Ω ik(x, ) ik(x, ) p0(x) dx is a Hilbert-Schmidt operator on H (see Theorem 4) with k being the reproducing kernel and denoting the tensor product. When H is a finite-dimensional RKHS, we show that the estimator enjoys parametric rates of convergence, i.e.,

f0 ˆfn H = Op0(n 1/2), KL(p0 p ˆfn) = Op0(n 1) and J(p0 p ˆfn) = Op0(n 1).

Note that the convergence rates are obtained under a non-classical smoothness assumption on f0, namely that it lies in the image of certain fractional power of C, which reduces to a more classical assumption if we choose k to be a Mat ern kernel (see Section 2 for its definition), as it induces a Sobolev space. In Section 4.2, we discuss in detail the smoothness assumption on f0 for the Gaussian (Example 2) and Mat ern (Example 3) kernels. Another interesting point to observe is that unlike in the classical function estimation methods (e.g., kernel density estimation and regression), the rates presented above for the proposed estimator tend to saturate for β > 1 (β > 1

2 w.r.t. J), with the best rate attained at β = 1 (β = 1

2 w.r.t. J), which means the smoothness of f0 is not fully captured by the estimator. Such a saturation behavior is well-studied in the inverse problem literature (Engl et al., 1996) where it has been attributed to the choice of regularizer. In Section 4.3, we discuss alternative regularization strategies using ideas from Bauer et al. (2007), which covers nonparametric least squares regression: we show that for appropriately chosen regularizers, the above mentioned rates hold for any β > 0, and do not saturate for the aforementioned ranges of β (see Theorem 9).

(iii) In Section 5, we study the problem of density estimation in the misspecified setting, i.e., p0 / P, which is not addressed in Hyv arinen (2007) and Fukumizu (2009). Using a more sophisticated analysis than in the well-specified case, we show in Theorem 12 that J(p0 p ˆfn) infp P J(p0 p) as n . Under an appropriate smoothness assumption

Density Estimation in Infinite Dimensional Exponential Families

q0 (see the statement of Theorem 12 for details), we show that J(p0 p ˆfn) 0 as n along with a rate for this convergence, even though p0 / P. However, unlike in the well-specified case, where the consistency is obtained not only in J but also in other distances, we obtain convergence only in J for the misspecified case. Note that while Barron and Sheu (1991) considered the estimation of p0 in the misspecified setting, the results are restricted to the approximating families consisting of polynomials, splines, or trigonometric series. Our results are more general, as they hold for abstract RKHSs.

(iv) In Section 6, we present preliminary numerical results comparing the proposed estimator with KDE in estimating a Gaussian and mixture of Gaussians, with the goal of empirically evaluating performance as d gets large for a fixed sample size. In these two estimation problems, we show that the proposed estimator outperforms KDE, and the advantage grows as d increases. Inspired by this preliminary empirical investigation, our proposed estimator (or computationally efficient approximations) has been used by Strathmann et al. (2015) in a gradient-free adaptive MCMC sampler, and by Sun et al. (2015) for graphical model structure learning. These applications demonstrate the practicality and performance of the proposed estimator.

Finally, we would like to make clear that our principal goal is not to construct density estimators that improve uniformly upon KDE, but to provide a novel flexible modeling technique for approximating an unknown density by a rich parametric family of densities, with the parameter being infinite dimensional, in contrast to the classical approach of finite dimensional approximation.

Various notations and definitions that are used throughout the paper are collected in Section 2. The proofs of the results are provided in Section 8, along with some supplementary results in an appendix.

2. Definitions & Notation

We introduce the notation used throughout the paper. Define [d] := {1, . . . , d}. For a :=

(a1, . . . , ad) Rd and b := (b1, . . . , bd) Rd, a 2 := q Pd i=1 a2 i and a, b := Pd i=1 aibi. For a, b > 0, we write a b if a γb for some positive universal constant γ. For a topological space X, C(X) (resp. Cb(X)) denotes the space of all continuous (resp. bounded continuous) functions on X. For a locally compact Hausdorffspace X, f C(X) is said to vanish at infinity if for every ϵ > 0 the set {x : |f(x)| ϵ} is compact. The class of all continuous f on X which vanish at infinity is denoted as C0(X). For open X Rd, C1(X) denotes the space of continuously differentiable functions on X. For f Cb(X), f := supx X |f(x)| denotes the supremum norm of f. Mb(X) denotes the set of all finite Borel measures on X. For µ Mb(X), Lr(X, µ) denotes the Banach space of r-power (r 1) µ-integrable functions. For X Rd, we will use Lr(X) for Lr(X, µ) if µ is a Lebesgue measure on X. For f Lp(X, µ), f Lr(X,µ) := R

X |f|r dµ 1/r denotes the Lr-norm of f for 1 r < and we denote it as Lr(X) if X Rd and µ is the Lebesgue measure. The convolution f g of two measurable functions f and g on Rd is defined as

(f g)(x) := Z

Rd f(y)g(x y) dy,

Sriperumbudur et al.

provided the integral exists for all x Rd. The Fourier transform of f L1(Rd) is defined as f (y) = (2π) d/2 Z

Rd f(x) e i y,x dx

where i denotes the imaginary unit 1.

In the following, for the sake of completeness and simplicity, we present definitions restricted to Hilbert spaces. Let H1 and H2 be abstract Hilbert spaces. A map S : H1 H2 is called a linear operator if it satisfies S(αx) = αSx and S(x+x ) = Sx+Sx for all α R and x, x H1, where Sx := S(x). A linear operator S is said to be bounded, i.e., the image SBH1 of BH1 under S is bounded if and only if there exists a constant c [0, ) such that for all x H1 we have Sx H2 c x H1, where BH1 := {x H1 : x H1 1}. In this case, the operator norm of S is defined as S := sup{ Sx H2 : x BH1}. Define L(H1, H2) be the space of bounded linear operators from H1 to H2. S L(H1, H2) is said to be compact if SBH1 is a compact subset in H2. The adjoint operator S : H2 H1 of S L(H1, H2) is defined by x, S y H1 = Sx, y H2, x H1, y H2. S L(H) := L(H, H) is called self-adjoint if S = S and is called positive if Sx, x H 0 for all x H. α R is called an eigenvalue of S L(H) if there exists an x = 0 such that Sx = αx and such an x is called the eigenvector of S and α. For compact, positive, self-adjoint S L(H), Sr : H H, r 0 is called a fractional power of S and S1/2 is the square root of S, which we write as

S := S1/2. An operator S L(H1, H2) is Hilbert-Schmidt if S HS := (P j J Sej 2 H2)1/2 < where (ej)j J is an arbitrary orthonormal basis of separable Hilbert space H1. S L(H1, H2) is said to be of trace class if P j J (S S)1/2ej, ej H1 < . For x H1 and y H2, x y is an element of the tensor product space H1 H2 which can also be seen as an operator from H2 to H1 as (x y)z = x y, z H2 for any z H2. R(S) denotes the range space (or image) of S.

A real-valued symmetric function k : X X R is called a positive definite (pd) kernel if, for all n N, α1, . . . , αn R and all x1, . . . , xn X, we have Pn i,j=1 αiαjk(xi, xj) 0. A function k : X X R, (x, y) 7 k(x, y) is a reproducing kernel of the Hilbert space (Hk, , Hk) of functions if and only if (i) y X, k(y, ) Hk and (ii) y X, f Hk, f, k(y, ) Hk = f(y) hold. If such a k exists, then Hk is called a reproducing kernel Hilbert space. Since k(x, ), k(y, ) Hk = k(x, y), x, y X, it is easy to show that every reproducing kernel (r.k.) k is symmetric and positive definite. Some examples of kernels that appear throughout the paper are: Gaussian kernel, k(x, y) = exp( σ x y 2 2), x, y Rd, σ > 0 that induces the following Gaussian RKHS,

Hk = Hσ := n f L2(Rd) C(Rd) : Z |f (ω)|2e ω 2 2/4σ dω < o ,

the inverse multiquadric kernel, k(x, y) = (1 + x y

c 2 2) β, x, y Rd, β > 0, c (0, ) and

the Mat ern kernel, k(x, y) = 21 β

Γ(β) x y β d/2 2 Kd/2 β( x y 2), x, y Rd, β > d/2 that

induces the Sobolev space, Hβ 2 ,

Hk = Hβ 2 := n f L2(Rd) C(Rd) : Z (1 + ω 2 2)β|f (ω)|2 dω < o ,

where Γ is the Gamma function, and Kv is the modified Bessel function of the third kind of order v (v controls the smoothness of k).

Density Estimation in Infinite Dimensional Exponential Families

For any real-valued function f defined on open X Rd, f is said to be m-times continuously differentiable if for α Nd 0 with |α| := Pd i=1 αi m, αf(x) = α1 1 . . . αd d f(x) = |α|

xα1 1 ... x αd d f(x) exists. A kernel k is said to be m-times continuously differentiable if

α,αk : X X R exists and is continuous for all α Nd 0 with |α| m where α,α := α1 1 . . . αd d α1 1+d . . . αd 2d . Corollary 4.36 in Steinwart and Christmann (2008) and Theorem 1 in Zhou (2008) state that if α,αk exists and is continuous, then αk(x, ) = α1 1 . . . αd d k(x, ) = |α|

xα1 1 ... x αd d k((x1, . . . , xd), ) Hk with x = (x1, . . . , xd) and for every f Hk, we have

αf(x) = αk(x, ), f Hk and α,αk(x, x ) = αk(x, ), αk(x , ) Hk.

Given two probability densities, p and q on Ω Rd, the Kullback-Leibler divergence (KL) and Hellinger distance (h) are defined as KL(p q) = R p(x) log p(x)

q(x) dx and h(p, q) = p q L2(Ω) respectively. We refer to p q L1(Ω) as the total variation (TV) distance between p and q.

3. Approximation of Densities by P

In this section, we first show that every finite dimensional exponential family, Pfin is generated by the family P induced by a finite dimensional RKHS, which naturally leads to the infinite dimensional generalization of Pfin when H is an infinite dimensional RKHS. Next, we investigate the approximation properties of P in Proposition 1 and Corollary 2 when H is an infinite dimensional RKHS.

Let us consider a r-parameter exponential family, Pfin with sufficient statistic T(x) := (T1(x), . . . , Tr(x)) and construct a Hilbert space, H = span{T1(x), . . . , Tr(x)}. It is easy to verify that P induced by H is exactly the same as Pfin since any f H can be written as f(x) = Pr i=1 θi Ti(x) for some (θi)r i=1 R. In fact, by defining the inner product between f = Pr i=1 θi Ti and g = Pr i=1 γi Ti as f, g H := Pr i=1 θiγi, it follows that H is an RKHS with the r.k. k(x, y) = T(x), T(y) Rr since f, k(x, ) H = Pr i=1 θi Ti(x) = f(x). Based on this equivalence between Pfin and P induced by a finite dimensional RKHS, it is therefore clear that P induced by a infinite dimensional RKHS is a strict generalization to Pfin with k( , x) playing the role of a sufficient statistic.

Example 1 The following are some popular examples of probability distributions that belong to Pfin. Here we show the corresponding RKHSs (H, k) that generate these distributions. In some of these examples, we choose q0(x) = 1 and ignore the fact that q0 is a probability distribution as assumed in the definition of P.

Exponential: Ω= R++ := R+\{0}, k(x, y) = xy.

Normal: Ω= R, k(x, y) = xy + x2y2.

Beta: Ω= (0, 1), k(x, y) = log x log y + log(1 x) log(1 y).

Gamma: Ω= R++, k(x, y) = log x log y + xy.

Inverse Gaussian: Ω= R++, k(x, y) = xy + 1

Poisson: Ω= N {0}, k(x, y) = xy, q0(x) = (x! e) 1.

Sriperumbudur et al.

Binomial: Ω= {0, . . . , m}, k(x, y) = xy, q0(x) = 2 m m c .

While Example 1 shows that all popular probability distributions are contained in P for an appropriate choice of finite-dimensional H, it is of interest to understand the richness of P (i.e., what class of distributions can be approximated arbitrarily well by P?) when H is an infinite dimensional RKHS. This is addressed by the following result, which is proved in Section 8.1.

Proposition 1 Define

P0 := n πf(x) = ef(x) A(f)q0(x), x Ω: f C0(Ω) o

where Ω Rd is locally compact Hausdorff. Suppose k(x, ) C0(Ω), x Ωand Z Z k(x, y) dµ(x) dµ(y) > 0, µ Mb(Ω)\{0}. (5)

Then P is dense in P0 w.r.t. Kullback-Leibler divergence, total variation (L1 norm) and Hellinger distances. In addition, if q0 L1(Ω) Lr(Ω) for some 1 < r , then P is also dense in P0 w.r.t. Lr norm.

A sufficient condition for Ω Rd to be locally compact Hausdorffis that it is either open or closed. Condition (5) is equivalent to k being c0-universal (Sriperumbudur et al., 2011, p. 2396). If k(x, y) = ψ(x y), x, y Ω= Rd where ψ Cb(Rd) L1(Rd), then (5) can be shown to be equivalent to supp(ψ ) = Rd (Sriperumbudur et al., 2011, Proposition 5). Examples of kernels that satisfy the conditions in Proposition 1 include the Gaussian, Mat ern and inverse multiquadrics. In fact, any compactly supported non-zero ψ Cb(Rd) satisfies the assumptions in Proposition 1 as supp(ψ ) = Rd (Sriperumbudur et al., 2010, Corollary 10). Though P0 is still a parametric family of densities indexed by a Banach space (here C0(Ω)), the following corollary (proved in Section 8.2) to Proposition 1 shows that a broad class of continuous densities are contained in P0 and therefore can be approximated arbitrarily well in Lr norm (1 r ), Hellinger distance, and KL divergence by P.

Corollary 2 Let q0 C(Ω) be a probability density such that q0(x) > 0 for all x Ω, where Ω Rd is locally compact Hausdorff. Suppose there exists a constant ℓsuch that for any ϵ > 0, R > 0 that satisfies | p(x)

q0(x) ℓ| ϵ for any x with x 2 > R. Define

Pc := p C(Ω) : Z

Ω p(x) dx = 1, p(x) 0, x Ωand p

q0 ℓ C0(Ω) .

Suppose k(x, ) C0(Ω), x Ωand (5) holds. Then P is dense in Pc w.r.t. KL divergence, TV and Hellinger distances. Moreover, if q0 L1(Ω) Lr(Ω) for some 1 < r , then P is also dense in Pc w.r.t. Lr norm.

By choosing Ωto be compact and q0 to be a uniform distribution on Ω, Corollary 2 reduces to an easily interpretable result that any continuous density p0 on Ωcan be approximated arbitrarily well by densities in P in KL, Hellinger and Lr (1 r ) distances.

Density Estimation in Infinite Dimensional Exponential Families

Similar to the results so far, an approximation result for P can also be obtained w.r.t. Fisher divergence (see Proposition 13). Since this result is heavily based on the notions and results developed in Section 5, we defer its presentation until that section. Briefly, this result states that if H is sufficiently rich (i.e., dense in an appropriate class of functions), then any p C1(Ω) with J(p q0) < can be approximated arbitrarily well by elements in P w.r.t. Fisher divergence, where q0 C1(Ω).

4. Density Estimation in P: Well-specified Case

In this section, we present our score matching estimator for an unknown density p0 := pf0 P (well-specified case) from i.i.d. random samples (Xa)n a=1 drawn from it. This involves choosing the minimizer of the (empirical) Fisher divergence between p0 and pf P as the estimator, ˆf which we show in Theorem 5 to be obtained by solving a simple finitedimensional linear system. In contrast, we would like to remind the reader that the MLE is infeasible in practice due to the difficulty in handling A(f). The consistency and convergence rates of ˆf F and the plug-in estimator p ˆf are provided in Section 4.1 (see Theorems 6 and 7). Before we proceed, we list the assumptions on p0, q0 and H that we need in our analysis.

(A) Ωis a non-empty open subset of Rd with a piecewise smooth boundary Ω:= Ω\Ω, where Ωdenotes the closure of Ω.

(B) p0 is continuously extendible to Ω. k is twice continuously differentiable on Ω Ω with continuous extension of α,αk to Ω Ωfor |α| 2.

(C) i i+dk(x, x)p0(x) = 0 for x Ωand p

i i+dk(x, x)p0(x) = o( x 1 d 2 ) as x Ω, x 2 for all i [d].

(D) (ε-Integrability) For some ε 1 and i [d], i i+dk(x, x), q

2 i 2 i+dk(x, x) and p

i i+dk(x, x) i log q0(x) Lε(Ω, p0), where q0 C1(Ω).

Remark 3 (i) Ωbeing a subset of Rd along with k being continuous ensures that H is separable (Steinwart and Christmann, 2008, Lemma 4.33). The twice differentiability of k ensures that every f H is twice continuously differentiable (Steinwart and Christmann, 2008, Corollary 4.36). (C) ensures that J in (3) is equivalent to the one in (4) through integration by parts on Ω(see Corollary 7.6.2 in Duistermaat and Kolk, 2004 for integration by parts on bounded subsets of Rd which can be extended to unbounded Ωthrough a truncation and limiting argument) for densities in P. In particular, (C) ensures that R

Ω if(x) ip0(x) dx = R

Ω 2 i f(x)p0(x) dx for all f H and i [d], which will be critical to prove the representation in Theorem 4(ii), upon which rest of the results depend. The decay condition in (C) can be weakened to p

i i+dk(x, x)p0(x) = o( x 1 d 2 ) as x Ω, x 2 for all i [d] if Ωis a (possibly unbounded) box where d = #{i [d]|(ai, bi) is unbounded}.

(ii) When ε = 1, the first condition in (D) ensures that J(p0 pf) < for any pf P. The other two conditions ensure the validity of the alternate representation for J(p0 pf) in (4) which will be useful in constructing estimators of p0 (see Theorem 4). Examples of

Sriperumbudur et al.

kernels that satisfy (D) are the Gaussian, Mat ern (with β > max{2, d/2}), and inverse multiquadric kernels, for which it is easy to show that there exists q0 that satisfies (D).

(iii) (Identifiability) The above list of assumptions do not include the identifiability condition that ensures pf1 = pf2 if and only if f1 = f2. It is clear that if constant functions are included in H, i.e., 1 H, then pf = pf+c for any c R. On the other hand, it can be shown that if 1 / H and supp(q0) = Ω, then pf1 = pf2 f1 = f2. A sufficient condition for 1 / H is k C0(Ω Ω). We do not explicitly impose the identifiability condition as a part of our blanket assumptions because the assumptions under which consistency and rates are obtained in Theorem 7 automatically ensure identifiability.

Under these assumptions, the following result proved in Section 8.3 shows that the problem of estimating p0 through the minimization of Fisher divergence reduces to the problem of estimating f0 through a weighted least squares minimization in H (see parts (i) and (ii)). This motivates the minimization of the regularized empirical weighted least squares (see part (iv)) to obtain an estimator fλ,n of f0, which is then used to construct the plug-in estimate pfλ,n of p0.

Theorem 4 Suppose (A) (D) hold with ε = 1. Then J(p0 pf) < for all f F. In addition, the following hold.

(i) For all f F,

J(f) := J(p0 pf) = 1

2 f f0, C(f f0) H , (6)

where C : H H, C := R

Ωp0(x) Pd i=1 ik(x, ) ik(x, ) dx is a trace-class positive operator with

i=1 ik(x, ) if(x) dx.

(ii) Alternatively,

2 f, Cf H + f, ξ H + J(p0 q0)

ik(x, ) i log q0(x) + 2 i k(x, ) dx H

and f0 satisfies Cf0 = ξ.

(iii) For any λ > 0, a unique minimizer fλ of Jλ(f) := J(f) + λ

2 f 2 H over H exists and is given by fλ = (C + λI) 1ξ = (C + λI) 1Cf0.

(iv) (Estimator of f0) Given samples (Xa)n a=1 drawn i.i.d. from p0, for any λ > 0, the unique minimizer fλ,n of ˆJλ(f) := ˆJ(f) + λ

2 f 2 H over H exists and is given by

fλ,n = ( ˆC + λI) 1ˆξ,

Density Estimation in Infinite Dimensional Exponential Families

where ˆJ(f) := 1

2 f, ˆCf H + f, ˆξ H + J(p0 q0), ˆC := 1

n Pn a=1 Pd i=1 ik(Xa, ) ik(Xa, ) and

ik(Xa, ) i log q0(Xa) + 2 i k(Xa, ) .

An advantage of the alternate formulation of J(f) in Theorem 4(ii) over (6) is that it provides a simple way to obtain an empirical estimate of J(f) by replacing C and ξ by their empirical estimators, ˆC and ˆξ respectively from finite samples drawn i.i.d. from p0, which is then used to obtain an estimator of f0. Note that the empirical estimate of J(f), i.e., ˆJ(f) depends only on ˆC and ˆξ which in turn depend on the known quantities, k and q0, and therefore fλ,n in Theorem 4(iv) should in principle be computable. In practice, however, it is not easy to compute the expression for fλ,n = ( ˆC + λI) 1ˆξ as it involves solving an infinite dimensional linear system. In Theorem 5 (proved in Section 8.4), we provide an alternative expression for fλ,n as a solution of a simple finite-dimensional linear system (see (7) and (8)), using the general representer theorem (see Theorem A.2). It is interesting to note that while the solution to J(f) in Theorem 4(ii) is obtained by solving a non-linear system, Cf0 = ξ (the system is non-linear as C depends on p0 which in turn depends on f0), its estimator fλ,n proposed in Theorem 4, is obtained by solving a simple linear system. In addition, we would like to highlight the fact that the proposed estimator, fλ,n is precisely the Tikhonov regularized solution (which is well-studied in the theory of linear inverse problems) to the ill-posed linear system ˆCf = ˆξ. We further discuss the choice of regularizer in Section 4.3 using ideas from the inverse problem literature.

An important remark we would like to make about Theorem 4 is that though J(f) in (6) is valid only for f F, as it is obtained from J(p0 pf) where p0, pf P, the expression f f0, C(f f0) H is valid for any f H, as it is finite under the assumption that (D) holds with ε = 1. Therefore, in Theorem 4(iii, iv), fλ and fλ,n are obtained by minimizing Jλ and ˆJλ over H instead of over F, as the latter does not yield a nice expression (unlike fλ and fλ,n, respectively). However, there is no guarantee that fλ,n F, and so the density estimator pfλ,n may not be valid. While this is not an issue when studying the convergence of fλ,n f0 H (see Theorem 6), the convergence of pfλ,n to p0 (in various distances) needs to be handled slightly differently depending on whether the kernel is bounded or not (see Theorems 7 and B.2). Note that when the kernel is bounded, we obtain F = H, which implies pfλ,n is valid.

Theorem 5 (Computation of fλ,n) Let fλ,n = arg inff H ˆJλ(f), where ˆJλ(f) is defined in Theorem 4(iv) and λ > 0. Then

fλ,n = ˆξ λ +

i=1 β(a 1)d+i ik(Xa, ), (7)

where ˆξ is defined in Theorem 4(iv) and β = (β(a 1)d+i)a,i is obtained by solving

(G + nλI) β = 1

Sriperumbudur et al.

with (G)(a 1)d+i,(b 1)d+j = i j+dk(Xa, Xb) and

(h)(a 1)d+i = ˆξ, ik(Xa, ) H = 1

j=1 i 2 j+dk(Xa, Xb) + i j+dk(Xa, Xb) j log q0(Xb).

We would like to highlight that though fλ,n requires solving a simple linear system in (8), it can still be computationally intensive when d and n are large as G is a nd nd matrix. This is still a better scenario than that of MLE, however, since computationally efficient methods exist to solve large linear systems such as (8), whereas MLE can be intractable due to the difficulty in handling the log-partition function (though it can be approximated). On the other hand, MLE is statistically well-understood, with consistency and convergence rates established in general for the problem of density estimation (van de Geer, 2000) and in particular for the problem at hand (Fukumizu, 2009). In order to ensure that fλ,n and pfλ,n are statistically useful, in the following section, we investigate their consistency and convergence rates under some smoothness conditions on f0.

4.1 Consistency and Rate of Convergence

In this section, we prove the consistency of fλ,n (see Theorem 6(i)) and pfλ,n (see Theorems 7 and B.2). Under the smoothness assumption that f0 R(Cβ) for some β > 0, we present convergence rates for fλ,n and pfλ,n in Theorems 6(ii), 7 and B.2. In reference to the following results, for simplicity we suppress the dependence of λ on n by defining λ := λn where (λn)n N (0, ).

Theorem 6 (Consistency and convergence rates for fλ,n) Suppose (A) (D) with ε = 2 hold.

(i) If f0 R(C), then fλ,n f0 H p0 0 as λ 0, λ n and n .

(ii) If f0 R(Cβ) for some β > 0, then for λ = n max n 1 4 , 1 2(β+1) o ,

fλ,n f0 H = Op0 n min n 1 4 , β 2(β+1) o as n .

(iii) If C 1 < , then for λ = n 1

2 , fλ,n f0 H = Op0(n 1/2) as n .

Remark (i) While Theorem 6 (proved in Section 8.5) provides an asymptotic behavior for fλ,n f0 H under conditions that depend on p0 (and are therefore not easy to check in practice), a non-asymptotic bound on fλ,n f0 H that holds for all n 1 can be obtained under stronger assumptions through an application of Bernstein s inequality in separable Hilbert spaces. For the sake of simplicity, we provided asymptotic results which are obtained through an application of Chebyshev s inequality.

(ii) The proof of Theorem 6(i) involves decomposing fλ,n f0 H into an estimation error part, E(λ, n) := fλ,n fλ H, and an approximation error part, A0(λ) := fλ f0 H, where fλ = (C + λI) 1Cf0. While E(λ, n) 0 as λ 0, λ n and n without any assumptions on f0 (see the proof in Section 8.5 for details), it is not reasonable to expect

Density Estimation in Infinite Dimensional Exponential Families

A0(λ) 0 as λ 0 without assuming f0 R(C). This is because, if f0 lies in the null space of C, then fλ is zero irrespective of λ and therefore cannot approximate f0.

(iii) The condition f0 R(C) is difficult to check in practice as it depends on p0 (which in turn depends on f0). However, since the null space of C is just constant functions if the kernel is bounded and supp(q0) = Ω(see Lemma 14 in Section 8.6 for details), assuming 1 / H yields that R(C) = H and therefore consistency can be attained under conditions that are easy to impose in practice. As mentioned in Remark 3(iii), the condition 1 / H ensures identifiability and a sufficient condition for it to hold is k C0(Ω Ω), which is satisfied by Gaussian, Mat ern and inverse multiquadric kernels.

(iv) It is well known that convergence rates are possible only if the quantity of interest (here f0) satisfies some additional conditions. In function estimation, this additional condition is classically imposed by assuming f0 to be sufficiently smooth, e.g., f0 lies in a Sobolev space of certain smoothness. By contrast, the smoothness condition in Theorem 6(ii) is imposed in an indirect manner by assuming f0 R(Cβ) for some β > 0 so that the results hold for abstract RKHSs and not just Sobolev spaces which then provides a rate, with the best rate being n 1/4 that is attained when β 1. While such a condition has already been used in various works (Caponnetto and Vito, 2007; Smale and Zhou, 2007; Fukumizu et al., 2013) in the context of non-parametric least squares regression, we explore it in more detail in Proposition 8, and Examples 2 and 3. Note that this condition is common in the inverse problem theory (see Engl, Hanke, and Neubauer, 1996), and it naturally arises here through the connection of fλ,n being a Tikhonov regularized solution to the ill-posed linear system ˆCf = ˆξ. An interesting observation about the rate is that it does not improve with increasing β (for β > 1), in contrast to the classical results in function estimation (e.g., kernel density estimation and kernel regression) where the rate improves with increasing smoothness. This issue is discussed in detail in Section 4.3.

(v) Since C 1 < only if H is finite-dimensional, we recover the parametric rate of n 1/2

in a finite-dimensional situation with an automatic choice for λ as n 1/2.While Theorem 6 provides statistical guarantees for parameter convergence, the question of primary interest is the convergence of pfλ,n to p0. This is guaranteed by the following result, which is proved in Section 8.6.

Theorem 7 (Consistency and rates for pfλ,n) Suppose (A) (D) with ε = 2 hold and k := supx Ωk(x, x) < . Assume supp(q0) = Ω. Then the following hold:

(i) For any 1 < r with q0 L1(Ω) Lr(Ω),

pfλ,n p0 Lr(Ω) 0, h(pfλ,n, p0) 0, KL(p0 pfλ,n) 0 as λ n , λ 0 and n .

In addition, if f0 R(Cβ) for some β > 0, then for λ = n max n 1 4 , 1 2(β+1) o ,

pfλ,n p0 Lr(Ω) = Op0(θn), h(p0, pfλ,n) = Op0(θn), KL(p0 pfλ,n) = Op0(θ2 n)

as n where θn := n min n 1 4 , β 2(β+1) o .

(ii) J(p0 pfλ,n) 0 as λn , λ 0 and n . In addition, if f0 R(Cβ) for some

Sriperumbudur et al.

β 0, then for λ = n max n 1 3 , 1 2(β+1) o ,

J(p0 pfλ,n) = Op0

n min n 2 3 , 2β+1

2(β+1) o as n .

(iii) If C 1 < , then θn = n 1

2 and J(p0 pfλ,n) = Op0(n 1) with λ = n 1

Remark (i) Comparing the results of Theorem 6(i) and Theorem 7(i) (for Lr, Hellinger and KL divergence), we would like to highlight that while the conditions on λ and n match in both the cases, the latter does not require f0 R(C) to ensure consistency. While f0 R(C) can be imposed in Theorem 7 to attain consistency, we replaced this condition with supp(q0) = Ω a simple and easy condition to work with which along with the boundedness of the kernel ensures that for any f0 H, there exists f0 R(C) such that p f0 = p0 (see Lemma 14).

(ii) In contrast to the results in Lr, Hellinger and KL divergence, consistency in J can be obtained with λ converging to zero at a rate faster than in these results. In addition, one can obtain rates in J with β = 0, i.e., no smoothness assumption on f0, while no rates are possible in other distances (the latter might also be an artifact of the proof technique, as these results are obtained through an application of Theorem 6(ii) in Lemma A.1) which is due to the fact that the convergence in these other distances is based on the convergence of fλ,n f0 H, which in turn involves convergence of A0(λ) := fλ f0 H to zero while the convergence in J is controlled by A 1

C(fλ f0) H which can be shown to behave

λ) as λ 0, without requiring any assumptions on f0 (see Proposition A.3). Indeed, as a further consequence, the rate of convergence in J is faster than in other distances.

(iii) An interesting aspect in Theorem 7 is that pfλ,n is consistent in various distances such as Lr, Hellinger and KL, despite being obtained by minimizing a different loss function, i.e., J. However, we will see in Section 5 that such nice results are difficult to obtain in the misspecified case, where consistency and rates are provided only in J.

While Theorem 7 addresses the case of bounded kernels, the case of unbounded kernels requires a technical modification. The reason for this modification, as alluded to in the discussion following Theorem 4, is due to the fact that fλ,n may not be in F when k is unbounded, and therefore the corresponding density estimator, pfλ,n may not be welldefined. In order to keep the main ideas intact, we discuss the unbounded case in detail in Section B.2 in Appendix B.

4.2 Range Space Assumption

While Theorems 6 and 7 are satisfactory from the point of view of consistency, we believe the presented rates are possibly not minimax optimal since these rates are valid for any RKHS that satisfies the conditions (A) (D) and does not capture the smoothness of k (and therefore the corresponding H). In other words, the rates presented in Theorems 6 and 7 should depend on the decay rate of the eigenvalues of C which in turn effectively captures the smoothness of H. However, we are not able to obtain such a result see the remark following the proof of Theorem 6 for a discussion. While these rates do not reflect

Density Estimation in Infinite Dimensional Exponential Families

the intrinsic smoothness of H, they are obtained under the smoothness assumption, i.e., range space condition that f0 R(Cβ) for some β > 0. This condition is quite different from the classical smoothness conditions that appear in non-parametric function estimation. While the range space assumption has been made in various earlier works (e.g., Caponnetto and Vito (2007); Smale and Zhou (2007); Fukumizu et al. (2013) in the context of nonparametric least square regression), in the following, we investigate the implicit smoothness assumptions that it makes on f0 in our context. To this end, first it is easy to show (see the proof of Proposition B.3 in Section B.3) that

i I ciφi : X

i I c2 i α 2β i <

where (αi)i I are the positive eigenvalues of C, (φi)i I are the corresponding eigenvectors that form an orthonormal basis for R(C), and I is an index set which is either finite (if H is finite-dimensional) or I = N with limi αi = 0 (if H is infinite dimensional). From (9) it is clear that larger the value of β, the faster is the decay of the Fourier coefficients (ci)i I, which in turn implies that the functions in R(Cβ) are smoother. Using (9), an interpretation can be provided for R(Cβ) (β > 0 and β / N) as interpolation spaces (see Section A.5 for the definition of interpolation spaces) between R(C β ) and R(C β ) where R(C0) := H (see Proposition B.3 for details). While it is not completely straightforward to obtain a sufficient condition for f0 R(Cβ), β N, the following result provides a necessary condition for f0 R(C) (and therefore a necessary condition for f0 R(Cβ), β > 1) for translation invariant kernels on Ω= Rd, whose proof is presented in Section 8.7.

Proposition 8 (Necessary condition) Suppose ψ, φ Cb(Rd) L1(Rd) are positive definite functions on Rd with Fourier transforms ψ and φ respectively. Let H and G be the RKHSs associated with k(x, y) = ψ(x y) and l(x, y) = φ(x y), x, y Rd respectively. For 1 r 2, suppose the following hold:

Rd ω 2 2ψ (ω) dω < ; (ii) φ

ψ < ; (iii) 2 2(ψ )2

φ L r 2 r (Rd); (iv) q0

Then f0 R(C) implies f0 G H.

In the following, we apply the above result in two examples involving Gaussian and Mat ern kernels to get insights into the range space assumption.

Example 2 (Gaussian kernel) Let ψ(x) = e σ x 2 with Hσ as its corresponding RKHS (see Section 2 for its definition). By Proposition 8, it is easy to verify that f0 R(C) implies f0 Hα Hσ for σ

2 < α σ. Since Hβ Hγ for β < γ (i.e., Gaussian RKHSs are nested), f0 R(C) ensures that f0 lies in H σ

2 +ϵ for arbitrary small ϵ > 0.

Example 3 (Mat ern kernel) Let ψ(x) = 21 s

2 2 Kd/2 s( x 2), x Rd with Hs 2(Rd)

as its corresponding RKHS (see Section 2 for its definition) where s > d

2. By Proposition 8, we have that for q0 L1(Rd), if f0 R(C), then f0 Hα 2 (Rd) Hs 2(Rd) for 1 + d

Sriperumbudur et al.

2. Since Hδ 2(Rd) Hγ 2 (Rd) for γ < δ (i.e., Sobolev spaces are nested), this

means f0 lies in H 2s 1 d

2 ϵ 2 (Rd) for arbitrarily small ϵ > 0, i.e., f0 has at least 2s 1 d

2 weak-derivatives. By the minimax theory (Tsybakov, 2009, Chapter 2), it is well known that for any α > δ 0,

inf ˆfn sup f0 Hα 2 (Rd) ˆfn f0 Hδ 2(Rd) n α δ 2(α δ)+d , (10)

where the infimum is taken over all possible estimators. Here an bn means that for any two sequences an, bn > 0, an/bn is bounded away from zero and infinity as n . Suppose

f0 / Hα 2 (Rd) for α 2s 1 d

2, which means f0 H 2s 1 d

2 ϵ 2 (Rd) for arbitrarily small ϵ > 0. This implies that the rate of n 1/4 obtained in Theorem 6 is minimax optimal if H is chosen to be H1+d+ϵ 2 (Rd) (i.e., choose α = 2s 1 d

2 ϵ and δ = s in (10) and solve for s by equating the exponent in the r.h.s. of (10) to 1

4). Similarly, it can be shown that if q0 L2(Rd), then the rate of n 1/4 in Theorem 6 is minimax optimal if H is chosen

to be H 1+ d

2 +ϵ 2 (Rd). This example also explains away the dimension independence of the rate provided by Theorem 6 by showing that the dimension effect is captured in the relative smoothness of f0 w.r.t. H.

While Example 3 provides some understanding about the minimax optimality of fλ,n under additional assumptions on f0, the problem is not completely resolved. In the following section, however, we show that the rate in Theorem 6 is not optimal for β > 1, and that improved rates can be obtained by choosing the regularizer appropriately.

4.3 Choice of Regularizer

We understand from the characterization of R(Cβ) in (9) that larger β values yield smoother functions in H. However, the smoothness of f0 R(Cβ) for β > 1 is not captured in the rates in Theorem 6(ii), where the rate saturates at β = 1 providing the best possible rate of n 1/4 (irrespective of the size of β). This is unsatisfactory on the part of the estimator, as it does not effectively capture the smoothness of f0, i.e., the estimator is not adaptive to the smoothness of f0. We remind the reader that the estimator fλ,n is obtained by minimizing the regularized empirical Fisher divergence (see Theorem 4(iv)) yielding fλ,n = ( ˆC + λI) 1ˆξ, which can be seen as a heuristic to solve the (non-linear) inverse problem Cf0 = ξ (see Theorem 4(ii)) from finite samples, by replacing C and ξ with their empirical counterparts. This heuristic, which ensures that the finite sample inverse problem is well-posed, is popular in inverse problem literature under the name of Tikhonov regularization (Engl et al., 1996, Chapter 5). Note that Tikhonov regularization helps to make the ill-posed inverse problem a well-posed one by approximating α 1 by (α + λ) 1, λ > 0, where α 1 appears as the inverse of the eigenvalues of C while computing C 1. In other words, if ˆC is invertible, then an estimate of f0 can be obtained as ˆfn = ˆC 1ˆξ, i.e., ˆfn = P i I ˆξ,ˆφi H

ˆαi ˆφi, where (ˆαi)i I and (ˆφi)i I are the eigenvalues and eigenvectors of ˆC respectively. However, ˆC being a rank n operator defined on H (which can be infinite dimensional) is not invertible and therefore the regularized estimator is constructed as fλ,n = gλ( ˆC)ˆξ where gλ( ˆC) is defined through functional calculus (see Engl, Hanke, and

Density Estimation in Infinite Dimensional Exponential Families

Neubauer, 1996, Section 2.3) as

gλ( ˆC) = X

i I gλ(ˆαi) , ˆφi H ˆφi

with gλ : R+ R and gλ(α) := (α + λ) 1. Since the Tikhonov regularization is wellknown to saturate (as explained above) see Engl et al. (1996, Sections 4.2 and 5.1) for details , better approximations to α 1 have been used in the inverse problems literature to improve the rates by using gλ other than ( + λ) 1 where gλ(α) α 1 as λ 0. In the statistical context, Rosasco et al. (2005) and Bauer et al. (2007) have used the ideas from Engl et al. (1996) in non-parametric regression for learning a square integrable function from finite samples through regularization in RKHS. In the following, we use these ideas to construct an alternate estimator for f0 (and therefore for p0) that appropriately captures the smoothness of f0 by providing a better convergence rate when β > 1. To this end, we need the following assumption quoted from Engl et al. (1996, Theorems 4.1 4.3 and Corollary 4.4) and Bauer et al. (2007, Definition 1) that is standard in the theory of inverse problems.

(E) There exists finite positive constants Ag, Bg, Cg, η0 and (γη)η (0,η0] (all independent of λ > 0) such that gλ : [0, χ] R satisfies:

(a) supα D |αgλ(α)| Ag, (b) supα D |gλ(α)| Bg

λ , (c) supα D |1 αgλ(α)| Cg and (d) supα D |1 αgλ(α)|αη γηλη, η (0, η0] where D := [0, χ] and χ := d supx Ω,i [d] i i+dk(x, x) < .

The constant η0 is called the qualification of gλ which is what determines the point of saturation of gλ. We show in Theorem 9 that if gλ has a finite qualification, then the resultant estimator cannot fully exploit the smoothness of f0 and therefore the rate of convergence will suffer for β > η0. Given gλ that satisfies (E), we construct our estimator of f0 as fg,λ,n = gλ( ˆC)ˆξ.

Note that the above estimator can be obtained by using the data dependent regularizer, 1 2 f, ((gλ( ˆC)) 1 ˆC)f H in the minimization of ˆJ(f) defined in Theorem 4(iv), i.e.,

fg,λ,n = arg inf f H ˆJ(f) + 1

2 f, ((gλ( ˆC)) 1 ˆC)f H.

However, unlike fλ,n for which a simple form is available in Theorem 5 by solving a linear system, we are not able to obtain such a nice expression for fg,λ,n. The following result (proved in Section 8.8) presents an analog of Theorems 6 and 7 for the new estimators, fg,λ,n and pfg,λ,n.

Theorem 9 (Consistency and convergence rates for fg,λ,n and pfg,λ,n) Suppose (A) (E) hold with ε = 2.

(i) If f0 R(Cβ) for some β > 0, then for any λ n 1/2,

fg,λ,n f0 H = Op0 (θn) ,

Sriperumbudur et al.

where θn := n min n β 2(β+1), η0 2(η0+1) o with λ = n max n 1 2(β+1), 1 2(η0+1) o . In addition, if k < , then for any 1 < r with q0 L1(Ω) Lr(Ω),

pfg,λ,n p0 Lr(Ω) = Op0(θn), h(p0, pfg,λ,n) = Op0(θn) and KL(p0 pfg,λ,n) = Op0(θ2 n).

(ii) If f0 R(Cβ) for some β 0, then for any λ n 1/2,

J(p0 pfg,λ,n) = Op0

n min{2β+1,2η0} min{2β+2,2η0+1}

with λ = n 1 min{2β+2,2η0+1} .

(iii) If C 1 < , then for any λ n 1/2,

fg,λ,n f0 H = Op0(θn) and J(p0 pfg,λ,n) = Op0(θ2 n)

with θn = n 1

2 and λ = n 1 min{2,2η0} .

Theorem 9 shows that if gλ has infinite qualification, then smoothness of f0 is fully captured in the rates and as β , we attain Op0(n 1/2) rate for fg,λ,n f0 H in contrast to n 1/4

(similar improved rates are also obtained for pfg,λ,n in various distances) in Theorem 6. In the following example, we present two choices of gλ that improve on Tikhonov regularization. We refer the reader to Rosasco et al. (2005, Section 3.1) for more examples of gλ.

Example 4 (Choices of gλ) (i) Tikhonov regularization involves gλ(α) = (α + λ) 1 for which it is easy to verify that η0 = 1 and therefore the rates saturate at β = 1, leading to the results in Theorems 6 and 7.

(ii) Showalter s method and spectral cut-offuse

gλ(α) = 1 e α/λ

α and gλ(α) =

( 1 α, α λ 0, α < λ

respectively for which it is easy to verify that η0 = + (see Engl, Hanke, and Neubauer, 1996, Examples 4.7 & 4.8 for details) and therefore improved rates are obtained for β > 1 in Theorem 9 compared to that of Tikhonov regularization.

5. Density Estimation in P: Misspecified Case

In this section, we analyze the misspecified case where p0 / P, which is a more reasonable case than the well-specified one, as in practice it is not easy to check whether p0 P. To this end, we consider the same estimator pfλ,n as considered in the well-specified case where fλ,n is obtained from Theorem 5. The following result shows that J(p0 pfλ,n) infp P J(p0 p) as λ 0, λn and n under the assumption that there exists f F such that J(p0 pf ) = infp P J(p0 p). We present the result for bounded kernels although it can be easily extended to unbounded kernels as in Theorem B.2. Also, the presented result for Tikhonov regularization extends easily to pfg,λ,n using the ideas in the proof of Theorem 9. Note that unlike in the well-specified case where convergence in other distances can be shown even though the estimator is constructed from J, it is difficult to show such a result in the misspecified case.

Density Estimation in Infinite Dimensional Exponential Families

Theorem 10 Let p0, q0 C1(Ω) be probability densities such that J(p0 q0) < where Ωsatisfies (A). Assume that (B), (C) and (D) with ε = 2 hold. Suppose k < , supp(q0) = Ωand there exists f F such that

J(p0 pf ) = inf p P J(p0 p).

Then for an estimator pfλ,n constructed from random samples (Xa)n a=1 drawn i.i.d. from p0, where fλ,n is defined in (7) also see Theorem 4(iv) with λ > 0, we have

J(p0 pfλ,n) inf p P J(p0 p) as λ 0, λn and n .

In addition, if f R(Cβ) for some β 0, then

J(p0 pfλ,n) r

inf p P J(p0 p) + Op0

n min n 1 3 , 2β+1

with λ = n max n 1 3 , 1 2(β+1) o . If C 1 < , then for λ = n 1

J(p0 pfλ,n) r

inf p P J(p0 p) + Op0(n 1/2).

with λ = n 1

While the above result is useful and interesting, the assumption about the existence of f

is quite restrictive. This is because if p0 (which is not in P) belongs to a family Q where P is dense in Q w.r.t. J, then there is no f H that attains the infimum, i.e., f does not exist and therefore the proof technique employed in Theorem 10 will fail. In the following, we present a result (Theorem 12) that does not require the existence of f but attains the same result as in Theorem 10, but requiring a more complicated proof. Before we present Theorem 12, we need to introduce some notation.

To this end, let us return to the objective function under consideration,

J(p0 pf) = 1

Ω p0(x) log p0

i=1 ( if if)2 dx,

where f = log p0

q0 and p0 / P. Define

W2(Ω, p0) := f C1(Ω) : αf L2(Ω, p0), |α| = 1 .

This is a reasonable class of functions to consider as under the condition J(p0 q0) < , it is clear that f W2(Ω, p0). Endowed with a semi-norm,

f 2 W2 := X

|α|=1 αf 2 L2(Ω,p0),

W2(Ω, p0) is a vector space of functions, from which a normed space can be constructed as follows. Let us define f, f W2(Ω, p0) to be equivalent, i.e., f f , if f f W2 = 0.

Sriperumbudur et al.

In other words, f f if and only if f and f differ by a constant p0-almost everywhere. Now define the quotient space W 2 (Ω, p0) := {[f] : f W2(Ω, p0)} where [f] := {f W2(Ω, p0) : f f } denotes the equivalence class of f. Defining [f] W 2 := f W2, it is easy to verify that W 2 defines a norm on W 2 (p0). In addition, endowing the following bilinear form on W 2 (Ω, p0)

[f] , [g] W 2 := Z

|α|=1 ( αf)(x)( αg)(x) dx

makes it a pre-Hilbert space. Let W2(Ω, p0) be the Hilbert space obtained by completion of W 2 (Ω, p0). As shown in Proposition 11 below, under some assumptions, a continuous mapping Ik : H W2(Ω, p0), f 7 [f] can be defined, which is injective modulo constant functions. Since addition of a constant does not contribute to pf, the space W2(Ω, p0) can be regarded as a parameter space extended from H. In addition to Ik, Proposition 11 (proved in Section 8.10) describes the adjoint of Ik and relevant self-adjoint operators, which will be useful in analyzing pfλ,n in Theorem 12.

Proposition 11 Let supp(q0) = Ωwhere Ω Rd is non-empty and open. Suppose k satisfies (B) and i i+dk(x, x) L1(Ω, p0) for all i [d]. Then Ik : H W2(Ω, p0), f 7 [f] defines a continuous mapping with the null space H R. The adjoint of Ik is Sk : W2(Ω, p0) H whose restriction to W 2 (Ω, p0) is given by

Sk[h] (y) = Z

i=1 ik(x, y) ih(x) p0(x) dx, [h] W 2 (Ω, p0), y Ω.

In addition, Ik and Sk are Hilbert-Schmidt and therefore compact. Also, Ek := Sk Ik and Tk := Ik Sk are compact, positive and self-adjoint operators on H and W2(Ω, p0) respectively where

i=1 ik(x, y) ig(x)p0(x) dx, g H, y Ω

and the restriction of Tk to W 2 (Ω, p0) is given by

i=1 ik(x, ) ih(x) p0(x) dx

, [h] W 2 (Ω, p0).

Note that for [h] W 2 (Ω, p0), the derivatives ih do not depend on the choice of a representative element almost surely w.r.t. p0, and thus the above integrals are well defined. Having constructed W2(Ω, p0), it is clear that J(p0 pf) = 1

2 [f ] Ikf 2 W2, which means estimating p0 is equivalent to estimating f W2(Ω, p0) by f F. With all these preparations, we are now ready to present a result (proved in Section 8.11) on consistency and convergence rate for pfλ,n without assuming the existence of f .

Theorem 12 Let p0, q0 C1(Ω) be probability densities such that J(p0 q0) < . Assume that (A) (D) hold with ε = 2 and χ := d supx Ω,i [d] i i+dk(x, x) < . Then the following

Density Estimation in Infinite Dimensional Exponential Families

(i) As λ 0, λn and n , J(p0 pfλ,n) infp P J(p0 p).

(ii) Define f := log p0

q0 . If [f ] R(Tk), then

J(p0 pfλ,n) 0 as λ 0, λn and n .

In addition, if [f ] R(T β k ) for some β > 0, then for λ = n max n 1 3 , 1 2β+1 o

J(p0 pfλ,n) = Op0 n min n 2 3 , 2β 2β+1 o

. (iii) If E 1 k < and T 1 k < , then J(p0 pfλ,n) = Op0 n 1 with λ = n 1

Remark (i) The result in Theorem 12(ii) is particularly interesting as it shows that [f ] W2(Ω, p0)\Ik(H) can be consistently estimated by fλ,n H, which in turn implies that certain p0 / P can be consistently estimated by pfλ,n P. In particular, if Sk is injective, then Ik(H) is dense in W2(Ω, p0) w.r.t. W2, which implies infp P J(p0||p) = 0 though there does not exist f H for which J(p0||pf ) = 0. While Theorem 10 cannot handle this situation, (i) and (ii) in Theorem 12 coincide showing that p0 / P can be consistently estimated by pfλ,n P. While the question of when Ik(H) is dense in W2(Ω, p0) is open, we refer the reader to Section B.4 for a related discussion.

(ii) Replicating the proof of Theorem 4.6 in Steinwart and Scovel (2012), it is easy to show that for all 0 < γ < 1, R(T γ/2 k ) = [W2(Ω, p0), Ik(H)]γ,2, where the r.h.s. is an interpolation space obtained through the real interpolation of W2(Ω, p0) and Ik(H) (see Section A.5 for the notation and definition). Here Ik(H) is endowed with the Hilbert space structure by Ik(H) = H/H R. This interpolation space interpretation means that, for β 1 2, R(T β k ) H modulo constant functions. It is nice to note that the rates in Theorem 12(ii) for β 1

2 match with the rates in Theorem 7 (i.e., the well-specified case) w.r.t. J for 0 β 1

2. We highlight the fact that β = 0 corresponds to H in Theorem 7 whereas β = 1

2 corresponds to H in Theorem 12(ii) and therefore the range of comparison is for β 1

2 in Theorem 12(ii) versus 0 β 1

2 in Theorem 7. In contrast, Theorem 10 is very limited as it only provides a rate for the convergence of J(p0||pfλ,n) to infp P J(p0||p) assuming that f is sufficiently smooth.

Based on the observation (i) in the above remark that infp P J(p0 p) = 0 if Ik(H) is dense in W2(Ω, p0) w.r.t. W2, it is possible to obtain an approximation result for P (similar to those discussed in Section 3) w.r.t. Fisher divergence as shown below, whose proof is provided in Section 8.12.

Proposition 13 Suppose Ω Rd is non-empty and open. Let q0 C1(Ω) be a probability density and

PFD := n p C1(Ω) : Z

Ω p(x) dx = 1, p(x) 0, x Ωand J(p q0) < o .

For any p PFD, if Ik(H) is dense in W2(Ω, p) w.r.t. W2, then for every ϵ > 0, there exists p P such that J(p p) ϵ.

Sriperumbudur et al.

6. Numerical Simulations

We have proposed an estimator of p0 that is obtained by minimizing the regularized empirical Fisher divergence and presented its consistency along with convergence rates. As discussed in Section 1, however one can simply ignore the structure of P and estimate p0 in a completely non-parametric fashion, for example using the kernel density estimator (KDE). In fact, consistency and convergence rates of KDE are also well-studied (Tsybakov, 2009, Chapter 1) and the kernel density estimator is very simple to compute requiring only O(n) computations compared to the proposed estimator, which is obtained by solving a linear system of size nd nd. This raises questions about the applicability of the proposed estimator in practice, though it is very well known that KDE performs poorly for moderate to large d (Wasserman, 2006, Section 6.5). In this section, we numerically demonstrate that the proposed score matching estimator performs significantly better than the KDE, and in particular, that the advantage with the proposed estimator grows as d gets large. Note further that the maximum likelihood approach of Barron and Sheu (1991) and Fukumizu (2009) does not yield estimators that are practically feasible, and therefore to the best of our knowledge, the proposed estimator is the only viable estimator for estimating densities through P.

In the following, we consider two simple scenarios of estimating a multivariate normal and mixture of normals using the proposed estimator and demonstrate the superior performance of the proposed estimator over KDE. Inspired by this preliminary empirical investigation, recently, the proposed estimator has been explored in two concrete applications of gradient-free adaptive MCMC sampler (Strathmann et al., 2015) and graphical model structure learning (Sun et al., 2015) where the superiority of working with the infinite dimensional family is demonstrated. We would like to again highlight that the goal of this work is not to construct density estimators that improve upon KDE but to provide a novel modeling technique of approximating an unknown density by a rich parametric family of densities with the parameter being infinite dimensional in contrast to the classical approach of finite dimensional approximation.

We consider the problems of estimating a standard normal distribution on Rd, N(0, Id) and mixture of Gaussians,

2φd(x; α1n, Id) + 1

2φd(x; β1n, Id)

through the score matching approach and KDE, and compare their estimation accuracies. Here φd(x; µ, Σ) is the p.d.f. of N(µ, ΣId). By choosing the kernel, k(x, y) = exp( x y 2 2 2σ2 )+ r(x T y+c)2, which is a Gaussian plus polynomial of degree 2, it is easy to verify that Gaussian distributions lie in P, and therefore the first problem considers the well-specified case while the second problem deals with the misspecified case. In our simulations, we chose r = 0.1, c = 0.5, α = 4 and β = 4. The base measure of the exponential family is N(0, 102Id). The bandwidth parameter σ is chosen by cross-validation (CV) of the objective function ˆJλ (see Theorem 4(iv)) within the parameter set {0.1, 0.2, 0.4, 0.6, 0.8, 1, 1.2, 1.4, 1.6} σ , where σ is the median of pairwise distances of data, and the regularization parameter λ is set as λ = 0.1 n 1/3 with sample size n. For KDE, the Gaussian kernel is used for the smoothing kernel, and the bandwidth parameter is chosen by CV from

Density Estimation in Infinite Dimensional Exponential Families

5 10 15 20 0

Socre objective function

Gaussian distribution: n = 500

Score match KDE

100 200 300 400 500 0

Sample Size

Socre objective function

Gaussian distribution: d = 7

Score match KDE

0 5 10 15 20 0

Score objective function

Gaussian mixture: n = 300

Score match KDE

100 200 300 400 500 0

Sample Size

Score objective function

Gaussian mixture: d = 7

Score match KDE

Figure 1: Experimental comparisons with the score objective function: proposed method and kernel density estimator

{0.02, 0.04, 0.06, 0.08, 0.1, 0.2, 0.4, 0.6, 0.8, 1.0} σ ; where for both the methods, 5-fold CV is applied. Since it is difficult to accurately estimate the normalization constant in the proposed method, we use two methods to evaluate the accuracy of estimation. One is the objective function for the score matching method,

2 | i log p(x)|2 + 2 i log p(x) p0(x)dx,

and the other is correlation of the estimator with the true density function,

Cor(p, p0) := ER[p(X)p0(X)] p

ER[p(X)2]ER[p0(X)2] ,

where R is a probability distribution. For R, we use the empirical distribution based on 10000 random samples drawn i.i.d. from p0(x).

Figures 1 and 2 show the score objective function ( J(p)) and the correlation (Cor(p, p0)) (along with their standard deviation as error bars) of the proposed estimator and KDE for the tasks of estimating a Gaussian and a mixture of Gaussians, for different sample sizes

Sriperumbudur et al.

5 10 15 20 0.65

Correlation

Gaussian distribution: n = 500

Score match KDE

100 200 300 400 500 0.8

Sample Size

Correlation

Gaussian distribution: d = 7

Score match KDE

0 5 10 15 20 0.5

Correlation

Gaussian mixture: n = 300

Score match KDE

100 200 300 400 500 0.7

Sample Size

Correlation

Gaussian mixture: d = 7

Score match KDE

Figure 2: Experimental comparisons with the correlation: proposed method and kernel density estimator

(n) and dimensions (d). From the figures, we see that the proposed estimator outperforms (i.e., lower function values) KDE in all the cases except the low dimensional cases ((n, d) = (500, 2) for the Gaussian, and (n, d) = (300, 2), (300, 4) for the Gaussian mixture). In the case of the correlation measure, the score matching method yields better results (i.e., higher correlation) besides in the Gaussian mixture cases of d = 2, 4, 6 (Fig.2, lower-left) and some cases of d = 7 (lower-right). The proposed method shows an increased advantage over KDE as the dimensionality increases, thereby demonstrating the advantage of the proposed estimator for high dimensional data.

7. Summary & Discussion

We have considered an infinite dimensional generalization, P, of the finite-dimensional exponential family, where the densities are indexed by functions in a reproducing kernel Hilbert space (RKHS), H. We showed that P is a rich object that can approximate a large class of probability densities arbitrarily well in Kullback-Leibler divergence, and addressed the main question of estimating an unknown density, p0 from finite samples drawn i.i.d. from it, in well-specified (p0 P) and misspecified (p0 / P) settings. We proposed a density

Density Estimation in Infinite Dimensional Exponential Families

estimator based on minimizing the regularized version of the empirical Fisher divergence, which results in solving a simple finite-dimensional linear system. Our estimator provides a computationally efficient alternative to maximum likelihood based estimators, which suffer from the computational intractability of the log-partition function. The proposed estimator is also shown to empirically outperform the classical kernel density estimator, with advantage increasing as the dimension of the space increases. In addition to these computational and empirical results, we have established the consistency and convergence rates under certain smoothness assumptions (e.g., log p0 R(Cβ)) for both well-specified and misspecified scenarios.

Three important questions still remain open in this work which we intend to address in our future work. First, the assumption log p0 R(Cβ) is not well understood. Though we presented a necessary condition for this assumption (with β = 1) to hold for bounded continuous translation invariant kernels on Rd, obtaining a sufficient condition can throw light on the minimax optimality of the proposed estimator. Another alternative is to directly study the minimax optimality of the rates for 0 < β 1 (for β > 1, we showed that the above mentioned rates can be improved by an appropriate choice of the regularizer) by obtaining minimax lower bounds under the source condition log p0 R(Cβ) and the eigenvalue decay rate of C, using the ideas in De Vore et al. (2004). Second, the proposed estimator depends on the regularization parameter, which in turn depends on the smoothness scale β. Since β is not known in practice, it is therefore of interest to construct estimators that are adaptive to unknown β. Third, since the proposed estimator is computationally expensive as it involves solving a linear system of size nd nd, it is important to study either alternate estimators or efficient implementations of the proposed estimator to improve the applicability of the method.

We provide proofs of the results presented in Sections 3 5.

8.1 Proof of Proposition 1

Sriperumbudur et al. (2011, Proposition 5) showed that H is dense in C0(Ω) w.r.t. uniform norm if and only if k satisfies (5). Therefore, the denseness in L1, KL and Hellinger distances follow trivially from Lemma A.1. For Lr norm (r > 1), the denseness follows by using the bound pf pg Lr(Ω) 2e2 f g e2 f f g q0 Lr(Ω) obtained from Lemma A.1(i) with f C0(Ω) and g H.

8.2 Proof of Corollary 2

For any p Pc, define pδ := p+δq0

1+δ . Note that pδ(x) > 0 for all x Ωand p pδ Lr(Ω) = δ p q0 Lr(Ω)

1+δ , implying that limδ 0 p pδ Lr(Ω) = 0 for any 1 r . This means, for any ϵ > 0, δϵ > 0 such that for any 0 < θ < δϵ, we have p pθ Lr(Ω) ϵ, where pθ(x) > 0 for all x Ω.

Sriperumbudur et al.

Define f := log pθ

q0 cθ where cθ := log ℓ+θ

1+θ. It is clear that f C(Ω) since p, q C(Ω). Fix any η > 0 and define

A := {x : f(x) η} = x : p(x)

q0(x) ℓ (ℓ+ θ) (eη 1) .

Since p q0 ℓ C0(Ω), it is clear that A is compact and so f C0(Ω). Also, it is easy to verify that pθ = ef A(f)q0 which implies pθ P0, where P0 is defined in Proposition 1. This means, for any ϵ > 0, there exists pg P such that pθ pg Lr(Ω) ϵ under the assumption that q0 L1(Ω) Lr(Ω). Therefore p pg Lr(Ω) 2ϵ for any 1 r , which proves

the denseness of P in Pc w.r.t. Lr norm for any 1 r . Since h(p, q) q

p q L1(Ω) for any probability densities p, q, the denseness in Hellinger distance follows.

We now prove the denseness in KL divergence by noting that

KL(p pδ) = Z

{p>0} p log p + pδ

p + q0δ dx Z

{p>0} p p + pδ

p + q0δ 1 dx

p>0 (p q0) p p + q0δ dx δ p q0 L1(Ω) 2δ,

which implies limδ 0 KL(p pδ) = 0. This implies, for any ϵ > 0, δϵ > 0 such that for any 0 < θ < δϵ, KL(p pθ) ϵ. Arguing as above, we have pθ P0, i.e., there exists f C0(Ω) such that pθ = efq0 R efq0 dx. Since H is dense in C0(Ω), for any f C0(Ω) and any ϵ > 0,

there exists g H such that f g ϵ. For pg P, since R p log pθ

pg dx log pθ

2 f g 2ϵ, we have

KL(p pg) = Z

and the result follows.

8.3 Proof of Theorem 4

(i) By the reproducing property of H, since if(x) = f, ik(x, ) H for all i [d], it is easy to verify that

i=1 f f0, ik(x, ) 2 H dx

i=1 f f0, ( ik(x, ) ik(x, )) (f f0) H dx

Ω p0(x) f f0, Cx(f f0) H dx, (11)

where in the second line, we used a, b 2 H = a, b H a, b H = a, (b b)a H for a, b H with H being a Hilbert space and

i=1 ik(x, ) ik(x, ). (12)

Density Estimation in Infinite Dimensional Exponential Families

Observe that for all x Ω, Cx is a Hilbert-Schmidt operator as Cx HS Pd i=1 ik(x, ) 2 H = Pd i=1 i i+dk(x, x) < and (f f0) (f f0) is also Hilbert-Schmidt as (f f0) (f f0) HS = f f0 2 H < . Therefore, (11) is equivalent to

Ω p0(x) (f f0) (f f0), Cx HS dx.

Since the first condition in (D) implies R

Ω Cx HSp0(x) dx < , Cx is p0-integrable in the Bochner sense (see Diestel and Uhl, 1977, Definition 1 and Theorem 2), and therefore it follows from Diestel and Uhl (1977, Theorem 6) that

(f f0) (f f0), Z

Ω Cx p0(x) dx

where C := R

ΩCx p0(x) dx is the Bochner integral of Cx, thereby yielding (6).

We now show that C is trace-class. Let (el)l N be an orthonormal basis in H (a countable ONB exists as H is separable see Remark 3(i)). Define B := P l Cel, el H so that

Ω el, Cxel Hp0(x) dx = X

i=1 el, ik(x, ) 2 H p0(x) dx

i [d],l el, ik(x, ) 2 H p0(x) dx ( ) = Z

i=1 ik(x, ) 2 H p0(x) dx < ,

which means C is trace-class and therefore compact. Here, we used monotone convergence theorem in ( ) and Parseval s identity in ( ). Note that C is positive since f, Cf H = R

Ωp0(x) f 2 2 dx 0, f H.

(ii) From (6), we have J(f) = 1 2 f, Cf H f, Cf0 H + 1

2 f0, Cf0 H. Using if0(x) = i log p0(x) i log q0(x) for all i [d], we obtain that for any f H,

f, Cf0 H = Z

i=1 if(x) if0(x) dx

i=1 if(x) ip0(x) dx Z

i=1 if(x) i log q0(x) dx

i=1 2 i f(x) dx Z

i=1 if(x) i log q0(x) dx

ξx z }| { d X

i=1 2 i k(x, ) + ik(x, ) i log q0(x)

H dx (c) = f, ξ H, (13)

where (b) follows from integration by parts under (C) and the equality in (c) is valid as ξx is Bochner p0-integrable under (D) with ε = 1. Therefore Cf0 = ξ. For the third term, f0, Cf0 H = R

Ωp0(x) Pd i=1 ( if0(x))2 dx and the result follows.

Sriperumbudur et al.

(iii) Define c0 := J(p0 q0). For any λ > 0, it is easy to verify that

2 (C + λI)1/2f + (C + λI) 1/2ξ 2 H 1

2 ξ, (C + λI) 1ξ H + c0.

Clearly, Jλ(f) is minimized if and only if (C + λI)1/2f = (C + λI) 1/2ξ and therefore fλ = (C + λI) 1ξ is the unique minimizer of Jλ(f).

(iv) Since (iv) is similar to (iii) with C replaced by ˆC and ξ replaced by ˆξ, we obtain fλ,n = ( ˆC + λI) 1ˆξ.

8.4 Proof of Theorem 5

We prove the result based on the general representer theorem (Theorem A.2). From Theorem 4(iv), we have

fλ,n = arg inf f H 1 2 f, ˆCf H + f, ˆξ H + λ

= arg inf f H 1 2n

i=1 f, ik(Xa, ) 2 H + f, ˆξ H + λ

= arg inf f H V ( f, φ1 H, . . . , f, φnd H, f, φnd+1 H) + λ

where V (θ1, . . . , θnd, θnd+1) := 1 2n Pn a=1 Pd i=1 θ2 (a 1)d+i + θnd+1, φ(a 1)d+i := ik(Xa, ), a

[n], i [d] and φnd+1 := ˆξ. Therefore, it follows from Theorem A.2 that

fλ,n = δˆξ +

i=1 β(a 1)d+iφ(a 1)d+i (14)

where δ and β satisfy

with K = G h h T ˆξ 2 H

. Since V z t

, (15) reduces to λδ + 1 = 0 and λβ +

nh = 0 yielding δ = 1

n G + λI)β = 1 nλh.

Remark Instead of using the general representer theorem (Theorem A.2), it is possible to see that the standard representer theorem (Kimeldorf and Wahba, 1971; Sch olkopf et al., 2001) gives a similar, but slightly different linear system, and the solutions are the same if K is non-singular. The general representer theorem yields that β and δ are solution to

, where F = 1

. On the other hand, by using the standard

representer theorem, it is easy to show that fλ,n has the form in (14) with δ and β being

solution to KF β δ

. Clearly, both the solutions match if K is invertible while

the latter has many solutions if K is not invertible.

Density Estimation in Infinite Dimensional Exponential Families

8.5 Proof of Theorem 6

fλ,n fλ = ( ˆC + λI) 1 ˆξ + ( ˆC + λI)fλ ( ) = ( ˆC + λI) 1 ˆξ + ˆCfλ + C(f0 fλ)

= ( ˆC + λI) 1(C ˆC)(fλ f0) ( ˆC + λI) 1(ˆξ + ˆCf0)

= ( ˆC + λI) 1(C ˆC)(fλ f0) ( ˆC + λI) 1(ˆξ ξ) + ( ˆC + λI) 1(C ˆC)f0,

where we used λfλ = C(f0 fλ) in ( ). Define S1 := ( ˆC + λI) 1(C ˆC)(fλ f0) H, S2 := ( ˆC + λI) 1(ˆξ ξ) H and S3 := ( ˆC + λI) 1(C ˆC)f0 H so that

fλ,n f0 H fλ,n fλ H + fλ f0 H S1 + S2 + S3 + A0(λ), (16)

where A0(λ) := fλ f0 H. We now bound S1, S2 and S3 using Proposition A.4. Note that C = R

ΩCx p0(x) dx where Cx is defined in (12) is a positive, self-adjoint, trace-class operator and (D) (with ε = 2) implies that

Ω Cx 2 HSp0(x) dx Z

i=1 ik(x, ) 2 H

Ω ik(x, ) 4 H p0(x) dx < .

Therefore, by Proposition A.4(i,iii),

S1 ( ˆC + λI) 1 (C ˆC)(fλ f0) H = Op0

S2 ( ˆC + λI) 1 ˆξ ξ H = Op0

where by using the technique in the proof of Proposition A.4(i), we show below that ˆξ

ξ H = Op0(n 1/2). Note that Ep0 ˆξ ξ 2 H =

Ω ξx 2 Hp0(x) dx ξ 2 H n

Ω ξx 2 Hp0(x) dx

n , where ξx H is defined in (13) and (D) (with ε = 2) implies that R

Ω ξx 2 Hp0(x) dx < . Therefore ˆξ ξ H = Op0(n 1/2) follows from an application of Chebyshev s inequality. Again using Proposition A.4(i,iii), we obtain that

S3 ( ˆC + λI) 1 (C ˆC)f0 H = Op0

Using the bounds in S1, S2 and S3 in (16), we obtain

fλ,n f0 H = Op0

1 λ n + A0(λ)

+ A0(λ). (20)

(i) By Proposition A.3(i), we have that A0(λ) 0 as λ 0 if f0 R(C). Therefore, it follows from (20) that fλ,n f0 H 0 as λ 0, λ n and n .

(ii) If f0 R(Cβ) for β > 0, it follows from Proposition A.3(ii) that

A0(λ) max{1, C β 1} C βf0 Hλmin{1,β}

Sriperumbudur et al.

and therefore the result follows by choosing λ = n max n 1 4 , 1 2(β+1) o .

(iii) Note that

S1 = ( ˆC + λI) 1(C ˆC)(fλ f0) H C( ˆC + λI) 1 C 1 (C ˆC)(fλ f0) H,

S2 = ( ˆC + λI) 1(ˆξ ξ) H C( ˆC + λI) 1 C 1 ˆξ ξ H,

S3 = ( ˆC + λI) 1(C ˆC)f0 H C( ˆC + λI) 1 C 1 (C ˆC)f0 H

and A0(λ) = fλ f0 H C 1 C(fλ f0) H.

It follows from Proposition A.4(v) that C( ˆC+λI) 1 1 for n c λ2 where c is a sufficiently large constant that depends on Pd i=1 R

Ω( i i+dk(x, x))2p0(x) dx but not on n and λ. Using the bounds on (C ˆC)(fλ f0) H, ˆξ ξ H and (C ˆC)f0 H from part (i) and the bound on C(fλ f0) H from Proposition A.3(ii), we therefore obtain

fλ,n f0 H Op0

as n and the result follows.

Remark Under slightly strong assumptions on the kernel, the bound on S1 in (17) can be improved to obtain S1 = Op0(n 1/2) while the one on S3 in (19) can be refined to obtain

where N(λ) := Tr((C + λI) 1C) is the intrinsic dimension of H. Using

the fact that N(λ) 1

λ, it is easy to verify that the latter is an improved bound than the one in (19). In addition S3 dominates S1. However, if S2 in (18) is not improved, then S2 dominates S3, thereby resulting in a bound that does not capture the smoothness of k (or the corresponding H). Unfortunately, even with a refined analysis (not reported here), we are not able to improve the bound on S2 wherein the difficulty lies with handling ξ.

8.6 Proof of Theorem 7

Before we prove the result, we present a lemma.

Lemma 14 Suppose supx Ωk(x, x) < and supp(q0) = Ω. Then F = H and for any f0 H there exists f0 R(C) such that p f0 = p0.

Proof Since supx Ωk(x, x) < , it implies that, for every f H, R

Ωef(x)q0(x) dx < and hence F = H. Also, under the assumptions on k and q0, it is easy to verify that supp(p0) = Ω, which implies

N(C) = f H : Z

Ω f 2 2 p0(x) dx = 0

is either R or {0}, where N(C) denotes the null space of C. Let f0 be the orthogonal projection of f0 onto R(C) = N(C) . Then f0 f0 R and therefore p f0 = pf0.

Density Estimation in Infinite Dimensional Exponential Families

Proof of Theorem 7. From Theorem 4(iii), fλ = (C +λI) 1Cf0 = (C +λI) 1C f0 where the second equality follows from the proof of Lemma 14. Now, carrying out the decomposition as in the proof of Theorem 6(i), we obtain fλ,n fλ = ( ˆC + λI) 1(C ˆC)(fλ f0) ( ˆC + λI) 1(ˆξ ξ) + ( ˆC + λI) 1(C ˆC) f0 and therefore,

fλ,n f0 H ( ˆC + λI) 1 (C ˆC)(fλ f0) H + ξ ˆξ H + (C ˆC) f0 H + A0(λ),

where A0(λ) = fλ f0 H. The bounds on these quantities follow those in the proof of Theorem 6(i) verbatim and so the consistency result in Theorem 6(i) holds for fλ,n f0 H. By Lemma 14, since pf0 = p f0, it is sufficient to consider the convergence of pfλ,n to p f0. Therefore, the convergence (along with rates) in Lr (for any 1 r ), Hellinger and KL distances follow from using the bound fλ,n f0 p

k fλ,n f0 H (obtained through the reproducing property of k) in Lemma A.1 and invoking Theorem 6.

In the following, we obtain a bound on J(p0 pfλ,n) = 1

C(fλ,n f0) 2 H. While one can trivially use the bound

C(fλ,n f0) 2 H

C 2 fλ,n f0 2 H to obtain a rate on J(p0 pfλ,n) through the result in Theorem 6(ii), a better rate can be obtained by carefully bounding

C(fλ,n f0) 2 H as shown below. Consider

C(fλ,n f0) H

C(fλ,n fλ) H + A 1

2 (λ) + A 1 2 (λ)

C( ˆC + λI) 1 (C ˆC)(fλ f0) H + ξ ˆξ H + (C ˆC) f0 H

2 (λ) + A 1 2 (λ),

C(fλ f0) H and A 1 2 (λ) :=

C( f0 f0) H. It follows from Theo-

rem 4(i) and Lemma 14 that A 1 2 (λ) = J(p0 p f0) = 0. Also it follows from Proposition A.4(v)

C( ˆC + λI) 1 1

λ for n c λ2 where c is a large enough constant that does not

depend on n and λ and depends only on Pd i=1 R ik(x, ) 4 H p0(x) dx. Using the bounds from the proof of Theorem 6(i) for the rest of the terms within parenthesis, we obtain

C(fλ,n f0) H Op0

2 (λ). (22)

The consistency result therefore follows from Proposition A.3(i) by noting that A 1

as λ 0. If f0 R(Cβ) for some β 0, then Proposition A.3(ii) yields A 1

max{1, C β 1

2 }λmin{1,β+ 1

2 } C βf0 H which when used in (22) provides the desired rate

with λ = n max{ 1

3 , 1 2(β+1)}. If C 1 < , then the result follows by noting

C(fλ,n f0) H

C fλ,n f0 H and invoking the bound in (21).

8.7 Proof of Proposition 8

Observation 1: By (Wendland, 2005, Theorem 10.12), we have

H = f L2(Rd) Cb(Rd) : f ψ L2(Rd) ,

Sriperumbudur et al.

where f is defined in L2 sense. Since Z

Rd |f (ω)| dω Z

Rd |f (ω)|2

Rd ψ (ω) dω 1

where we used ψ L1(Rd) (see Wendland, 2005, Corollary 6.12), we have f L1(Rd). Hence Plancherel s theorem and continuity of f along with the inverse Fourier transform of f allow to recover any f H pointwise from its Fourier transform as

f(x) = 1 (2π)d/2

Rd eix T ωf (ω) dω, x Rd. (23)

Observation 2: Since ψ L1(Rd) and ψ 0, we have for all j = 1, . . . , d, Z

Rd |ωj|ψ (ω) dω 2 = Z

Rd ψ (ω) dω 2 Z

Rd |ωj| ψ (ω) R

Rd ψ (ω) dω dω 2

Rd ψ (ω) dω Z

Rd |ωj|2ψ (ω) dω

Rd ψ (ω) dω Z

Rd ω 2ψ (ω) dω (i) < ,

where we used Jensen s inequality in ( ). This means ωjψ (ω) L1(Rd), j [d] which ensures the existence of its Fourier transform and so

jψ(x) = 1 (2π)d/2

Rd(iωj)ψ (ω)eix T ω dω, x Rd, j [d]. (24)

Observation 3: For g H, we have for all j [d],

Rd |ωj||g (ω)| dω Z

Rd |g (ω)|2

Rd |ωj|2ψ (ω) dω 1

Rd |g (ω)|2

Rd ω 2ψ (ω) dω 1

which implies ωjg (ω) L1(Rd), j = 1, . . . , d. Therefore,

jg(x) = 1 (2π)d/2

Rd(iωj)g (ω)eix T ω dω, x Rd, j [d]. (25)

Observation 4: For any g G, we have Z

Rd |g (ω)|2

ψ (ω) dω = Z

Rd |g (ω)|2

φ (ω) φ (ω) ψ (ω) dω g 2 G

which implies g H, i.e., G H.

We now use these observations to prove the result. Since f0 R(C), there exists g H such that f0 = Cg, which means

j=1 jk(x, y) j p0(x) dx

Density Estimation in Infinite Dimensional Exponential Families

Rd ei(x y)T ω(iωj)ψ (ω) dω jg(x) p0(x) dx

Rd eix T ω jg(x) p0(x) dx (iωj)ψ (ω)e iy T ω dω

(25) = 1 (2π)d/2

(i( )jg p 0 ) (ω)(iωj)ψ (ω)e iy T ω dω

which from (23) means f 0 (ω) = Pd j=1 (i( )jg p 0 ) ( ω)( iωj)ψ (ω) where we have invoked Fubini s theorem in ( ) and represents the convolution. Define Lr(Rd) := r and θ := r r 1. Consider

Rd |f 0 (ω)|2

φ (ω) dω = Z

(i( )jg p 0 ) ( ω)(iωj)

(ψ )2(ω)(φ (ω)) 1 dω

i( )jg p 0 ( ω)|ωj|

(ψ )2(ω)(φ (ω)) 1 dω

i( )jg p 0 2 ( ω) ω 2(ψ )2(ω)(φ (ω)) 1 dω

iωjg (ω) p 0 (ω) 2 ( )

2(ψ )2( )(φ ( )) 1 r 2 r

where in the following we show that Pd j=1 |iωjg (ω) p 0 (ω)|2 ( ) L θ 2 (Rd), i.e.,

iωjg (ω) p 0 (ω) 2 ( )

iωjg (ω) p 0 (ω) 2 ( ) θ

iωjg (ω) p 0 (ω) 2 θ

( ) Pd j=1 iωjg (ω) 2 1 p 0 2 θ

( ) p0 2 r Pd j=1 iωjg (ω) 2 1

where we have invoked generalized Young s inequality (Folland, 1999, Proposition 8.9) in ( ), Hausdorff-Young inequality (Folland, 1999, p. 253) in ( ), and observation 3 combined with (iv) in ( ). This shows that f0 R(C) f0 G, i.e., R(C) G.

8.8 Proof of Theorem 9

To prove Theorem 9, we need the following lemma (De Vito et al., 2012, Lemma 5), which is due to Andreas Maurer.

Lemma 15 Suppose A and B are self-adjoint Hilbert-Schmidt operators on a separable Hilbert space H with spectrum contained in the interval [a, b], and let (σi)i I and (τj)j J be

Sriperumbudur et al.

the eigenvalues of A and B, respectively. Given a function r : [a, b] R, if there exists a finite constant L such that |r(σi) r(τj)| L|σi τj|, i I, j J, then r(A) r(B) HS L A B HS.

Proof of Theorem 9. (i) The proof follows the ideas in the proof of Theorem 10 in Bauer et al. (2007), which is a more general result dealing with the smoothness condition, f0 R(Θ(C)) where Θ is operator monotone. Recall that Θ is operator monotone on [0, b] if for any pair of self-adjoint operators U, V with spectra in [0, b] such that U V , we have Θ(U) Θ(V ), where is the partial ordering for self-adjoint operators on some Hilbert space H, which means for any f H, f, Uf H f, V f H. In our case, we adapt the proof for Θ(C) = Cβ. Define rλ(α) := gλ(α)α 1. Since f0 R(Cβ), there exists h H such that f0 = Cβh, which yields

fg,λ,n f0 = gλ( ˆC)ˆξ f0 = gλ( ˆC)(ˆξ + ˆCf0) + rλ( ˆC)Cβh

= gλ( ˆC)(ˆξ ξ) + gλ( ˆC)(C ˆC)f0 + rλ( ˆC) ˆCβh + rλ( ˆC)(Cβ ˆCβ)h. (26)

fg,λ,n f0 H gλ( ˆC)(ˆξ ξ) H | {z } (A)

+ gλ( ˆC)( ˆC C)f0 H | {z } (B)

+ rλ( ˆC)Cβh H | {z } (C)

+ rλ( ˆC)(Cβ ˆCβ)h H | {z } (D)

We now bound (A) (D). Since (A) gλ( ˆC) ˆξ ξ H, we have (A) = Op0 1 λ n

used (b) in (E) and the bound on ˆξ ξ H from the proof of Theorem 6(i). Similarly, (B) gλ( ˆC) ( ˆC C)f0 H implies (B) = Op0 1 λ n

where (b) in (E) and Proposition A.4(i) are invoked. Also, (d) in (E) implies that

(C) rλ( ˆC) ˆCβ h H max{γβ, γη0}λmin{β,η0} C βf0 H.

(D) can be bounded as

(D) rλ( ˆC) Cβ ˆCβ C βf0 H.

We now consider two cases:

β 1: Since α 7 αθ is operator monotone on [0, χ] for 0 θ 1, by Theorem 1 in Bauer et al. (2007), there exists a constant cθ such that ˆCθ Cθ cθ ˆC C θ cθ ˆC C θ HS. We now obtain a bound on ˆC C HS. To this end, consider

E ˆC C 2 HS = E ˆC 2 HS C 2 HS

i=1 ik(x, ) ik(x, )

HS p0(x) dx d

Z ik(x, ) 4 p0(x) dx,

which by Chebyshev s inequality implies that

ˆC C HS = Op0(n 1/2)

Density Estimation in Infinite Dimensional Exponential Families

and therefore (D) = Op0(n β/2). Since λ n 1/2, we have (D) = Op0(λβ).

β > 1: Since α 7 αθ is Lipschitz on [0, χ] for θ 1, by Lemma 15, Cβ ˆCβ Cβ ˆCβ HS βχβ 1 C ˆC HS and therefore (C) = Op0(n 1/2).

Collecting all the above bounds, we obtain

fg,λ,n f0 H Op0

+ Op0 λmin{β,η0}

and the result follows. The proofs of the claims involving Lr, h and KL follow exactly the same ideas as in the proof of Theorem 7 by using the above bound on fg,λ,n f0 H in Lemma A.1.

(ii) We now bound J(p0 pfg,λ,n) =

C(fg,λ,n f0) 2 H as follows. Note that

C(fg,λ,n f0) = (

ˆC)(fg,λ,n f0) | {z } (I )

ˆC(fg,λ,n f0) | {z } (II )

We bound (I ) H as

ˆC fg,λ,n f0 H c 1

C ˆC HS fg,λ,n f0 H

+ Op0 λmin{β,η0}+ 1

where we used the fact that α 7 α is operator monotone along with λ n 1/2. Using (26), (II ) H can be bounded as

ˆCgλ( ˆC) ˆξ + ˆCf0 H + p

ˆCrλ( ˆC) ˆCβ C βf0 H

ˆCrλ( ˆC) Cβ ˆCβ C βf0 H

ˆCrλ( ˆC) ˆCβ (γβ+ 1

2 γη0)λmin{β+ 1

ˆCrλ( ˆC) (γ 1

2 γη0)λmin{ 1

with ˆCf0 + ˆξ and Cβ ˆCβ bounded as in part (i) above. Here (a b) := max{a, b}. Combining (I ) H and (II ) H, we obtain the required result.

(iii) The proof follows the ideas in the proof of Theorems 6 and 7. Consider fg,λ,n f0 = gλ( ˆC)( ˆCf0 + ˆξ) + rλ( ˆC)f0 so that

fg,λ,n f0 H C 1 Cgλ( ˆC)( ˆCf0 + ˆξ) H + C 1 Crλ( ˆC)f0 H

C 1 ˆCf0 + ˆξ H ˆCgλ( ˆC) + ˆC C gλ( ˆC)

+ C 1 f0 H ˆCrλ( ˆC) + ˆC C rλ( ˆC) .

Therefore fg,λ,n f0 H = Op0(n 1/2)+O λmin{1,η0} where we used the fact that λ n 1/2

and the result follows.

Sriperumbudur et al.

8.9 Proof of Theorem 10

Before we analyze J(p0 pfλ,n), we need a small calculation for notational convenience. For any probability densities p, q C1, it is clear that p

2J(p q) = log p log q 2 L2(p). We generalize this by defining p

2J(p q µ) := log p log q 2 L2(µ) .

Clearly, if µ = p, then J(p q µ) matches with J(p q). Therefore, for probability densities p, q, r C1, p

J(p q p) + p

J(q r p). (27)

Based on (27), we have r

inf p P J(p0 p) q

J(p0 pfλ,n p0) q

J(p0 pf p0) + q

J(pf pfλ,n p0)

inf p P J(p0 p p0) + q

J(pf pfλ,n p0)

inf p P J(p0 p) + 1

fλ,n f , C(fλ,n f ) H

inf p P J(p0 p) + 1

C(fλ,n f ) H (28)

inf p P J(p0 p) + 1

C(fλ,n fλ) H + 1

where A (λ) =

C(fλ f ) H. The result simply follows from the proof of Theorem 7, where we showed that

C(fλ,n fλ) H = Op0 1

and A (λ) = O(λmin{1,β+ 1

f R(Cβ) for β 0 as λ 0, n . When C 1 < , we bound

C(fλ,n f ) H in (28) as

C fλ,n f H where fλ,n f H is in turn bounded as in (21).

8.10 Proof of Proposition 11

For f H, we have

i=1 ( if)2 p0(x) dx f 2 H

i=1 ik(x, ) 2 H p0(x) dx < ,

which means f W2(Ω, p0) and therefore [f] W2(Ω, p0). Since Ikf W2 = [f] W 2 = f W2 c f H < where c is some constant, it is clear that Ik is a continuous map from H to W2(Ω, p0). The adjoint Sk : W2(Ω, p0) H of Ik : H W2(Ω, p0) is defined by the relation Skf, g H = f, Ikg W2, f W2(Ω, p0), g H. If f := [h] W 2 (Ω, p0), then

[h] , Ikg W2 = [h] , [g] W 2 = X

Ω ( αh)(x)( αg)(x) p0(x) dx.

Density Estimation in Infinite Dimensional Exponential Families

For y Ωand g = k( , y), this yields

Sk[h] (y) = Sk[h] , k( , y) H = [h] , Ikk( , y) W2 = Z

i=1 ik(x, y) ih(x)p0(x) dx.

We now show that Ik is Hilbert-Schmidt. Since H is separable, let (el)l 1 be an ONB of H. Then we have

l Ikel 2 W2 = X

i=1 ( iel(x))2 p0(x) dx = Z

l el, ik(x, ) 2 H p0(x) dx

i=1 ik(x, ) 2 H p0(x) dx < ,

which proves that Ik is Hilbert-Schmidt (hence compact) and therefore Sk is also Hilbert Schmidt and compact. The other assertions about Sk Ik and Ik Sk are straightforward.

8.11 Proof of Theorem 12

By slight abuse of notation, f is used to denote [f ] in the proof for simplicity. For f F, we have

J(p0 pf) = 1

2 Ikf f 2 W2 = 1

2 Ekf, f H Skf , f H + 1

Since k satisfies (C) it is easy to verify that Skf , f H = f, ξ H, f H (see proof of Theorem 4(ii)). This implies Skf = ξ and

J(p0 pf) = 1

2 Ekf, f H + f, ξ H + 1

2 f 2 W2, (29)

where ξ is defined in Theorem 4(ii), and Ek is precisely the operator C defined in Theorem 4(ii). Following the proof of Theorem 4(ii), for λ > 0, it is easy to show that the unique minimizer of the regularized objective, J(p0 pf) + λ

2 f 2 H exists and is given by

fλ = (Ek + λI) 1ξ = (Ek + λI) 1Skf . (30)

We would like to reiterate that (29) and (30) also match with their counterparts in Theorem 4 and therefore as in Theorem 4(iv), an estimator of f is given by fλ,n = ( ˆEk+λI) 1ˆξ. In other words, this is the same as in Theorem 4(iv) since ˆEk = ˆC, and can be solved by a simple linear system provided in Theorem 5. Here ˆEk is the empirical estimator of Ek. Now consider q

2 J(p0 pfλ,n) = Ikfλ,n f W2 Ik(fλ,n fλ) W2 + Ikfλ f W2

Ek(fλ,n fλ) H + B(λ), (31)

where B(λ) := Ikfλ f W2. The proof now proceeds using the following decomposition, equivalent to the one used in the proof of Theorem 6(i), i.e.,

fλ,n fλ = ( ˆEk + λI) 1ˆξ fλ

Sriperumbudur et al.

= ( ˆEk + λI) 1(ˆξ + ˆEkfλ + λfλ)

( ) = ( ˆEk + λI) 1(ˆξ + ˆEkfλ + Skf Ekfλ ˆSkf + ˆSkf ),

where we used (30) in ( ). ˆSkf is well-defined as it is the empirical version of the restriction of Sk to W 2 (p0). Since Skf Ekfλ = Sk(f Ikfλ) and ˆSkf ˆEkfλ = ˆSk(f Ikfλ), we have fλ,n fλ = ( ˆEk + λI) 1(ˆξ + ˆSkf ) + ( ˆEk + λI) 1( ˆSk Sk)(f Ikfλ)

Ek(fλ,n fλ) H p

Ek( ˆEk+λI) 1 ˆξ + ˆSkf H + ( ˆSk Sk)(f Ikfλ) H . (32)

It follows from Proposition A.4(v) that

Ek( ˆEk + λI) 1 1

for n c λ2 where c is a sufficiently large constant that does not depend on n and λ. Following the proof of Proposition A.4(i), we have

E ˆξ + ˆSkf 2 H = n 1

n ξ + Skf 2 H + 1

i=1 ik(x, ) if + ξx

wherein the first term is zero as Skf + ξ = 0 and since

i=1 ik(x, ) if + ξx

H 2 ξx 2 H + 2χ f 2 2,

the integral in the second term is finite because of (D) and f W2(Ω, p0). Therefore, an application of Chebyshev s inequality yields

ˆξ + ˆSkf H = Op0(n 1/2). (34)

We now show that ( ˆSk Sk)(f Ikfλ) H = Op0(B(λ)n 1/2). To this end, define g := f Ikfλ and consider

Ep0 ˆSkg Skg 2 H =

Ω Pd i=1 ik(x, ) ig(x) 2 Hp0(x) dx Skg 2 H n χ

which therefore yields the claim through an application of Chebyshev s inequality. Using this along with (33) and (34) in (32), and using the resulting bound in (31) yields

2 J(p0 pfλ,n) Op0

+ B(λ). (35)

(i) We bound B(λ) as follows. First note that

B(λ) = Ik(Sk Ik + λI) 1Skf f W2 = (Tk + λI) 1Tkf f W2

Density Estimation in Infinite Dimensional Exponential Families

and so for any h H, we have

B(λ) = (Tk + λI) 1Tkf f W2 ((Tk + λI) 1Tk I)(f Ikh) W2 | {z } (I)

+ (Tk + λI) 1Tk Ikh Ikh W2 | {z } (II)

Since Tk is a self-adjoint compact operator, there exists (αl)l N and ONB (φl)l N of R(Tk) so that Tk = P l αl φl, W2φl. Let (ψj)j N be the orthonormal basis of N(Tk). Then we have

αl αl + λ 1 2 f Ikh, φl 2 W2 + X

j f Ikh, ψj 2 W2

l f Ikh, φl 2 W2 + X

j f Ikh, ψj 2 W2 = f Ikh 2 W2. (37)

From (Tk + λI) 1Tk = Ik(Ek + λI) 1Sk and Sk Ikh = Ekh, we have

(II) = Ik(Ek + λI) 1Ekh Ikh W2 = p

Ek(Ek + λI) 1Ekh p

where the inequality follows from Proposition A.3(ii). Using (37) and (38) in (36), we obtain B(λ) f Ikh W2 + h H

λ, using which in (35) yields

2 J(p0 pfλ,n) f Ikh W2 + Op0

Since the above inequality holds for any h H, we therefore have

2 J(p0 pfλ,n) inf h H

λ h H + Op0

λ, W2(p0), Ik(H)) + Op0

where the K-functional is defined in (A.6). Note that Ik(H) = H/H R is continuously embedded in W(p0). From (A.6), it is clear that the K-functional as a function of t is an infimum over a family of affine linear and increasing functions and therefore is concave, continuous and increasing w.r.t. t. This means, in (39), as λ 0,

λ, W2(p0), Ik(H)) inf h H f Ikh W2 = r

2 inf p P J(p0 p).

Since J(p0 pfλ,n) infp P J(p0 p), we have that J(p0 pfλ,n) infp P J(p0 p) as λ 0, λn and n .

(ii) Recall B(λ) from (i). From Proposition A.3(i) it follows that B(λ) 0 as λ 0

if f R(Tk). Therefore, (35) reduces to q

2 J(p0 pfλ,n) Op0 1

+ B(λ) and the

consistency result follows. If f R(T β k ) for some β > 0, then the rates follow from

Sriperumbudur et al.

Proposition A.3 by noting that B(λ) max{1, Tk β 1}λmin{1,β} T β k f W2 and choosing

λ = n max n 1 3 , 1 2β+1 o .

(iii) This simply follows from an analysis similar to the one used in the proof of Theorem 6(iii).

8.12 Proof of Proposition 13

For any p PFD, define f := log p

q0 , which implies that [f] W2(p). Since Ik(H) is dense in W2(p), we have for any ϵ > 0, there exists g H such that [f] Ikg W2

2ϵ. For a given g H, pick pg P. Therefore,

J(p pg) = 1

Ω p(x) log p log pg 2 2 dx = 1

2 [f] Ikg 2 W2 ϵ

and the result follows.

Acknowledgments

We thank the action editor and two anonymous reviewers for their detailed comments that significantly improved the paper. Theorem A.2 and the proof are due to an anonymous reviewer. We also thank Dougal Sutherland for a careful reading of the paper which helped to fix minor errors. A part of the work was carried out while BKS was a Research Fellow in the Statistical Laboratory, Department of Pure Mathematics and Mathematical Statistics, University of Cambridge. BKS thanks Richard Nickl for many valuable comments and insightful discussions. KF is supported in part by MEXT Grant-in-Aid for Scientific Research on Innovative Areas 25120012.

A. Appendix: Technical Results

In this appendix, we present some technical results that are used in the proofs.

A.1 Bounds on Various Distances Between pf and pg

In the following result, claims (iii) and (iv) are quoted from Lemma 3.1 of van der Vaart and van Zanten (2008).

Lemma A.1 Define P := pf = ef A(f)q0 : f ℓ (Ω) , where q0 is a probability density on Ω Rd and ℓ (Ω) is the space of bounded measurable functions on Ω. Then for any pf, pg P , we have

(i) pf pg Lr(Ω) 2e2 f g e2 min{ f , g } f g q0 Lr(Ω) for any 1 r ;

(ii) pf pg L1(Ω) 2e f g f g ;

Density Estimation in Infinite Dimensional Exponential Families

(iii) KL(pf pg) c f g 2 e f g (1 + f g ) where c is a universal constant;

(iv) h(pf, pg) e f g /2 f g .

Proof (i) Define B(f) := R efq0 dx. Consider

pf pg Lr(Ω) = efq0 B(f) egq0

efq0B(g) egq0B(f) Lr(Ω) B(f)B(g)

efq0 (B(g) B(f)) + ef eg q0B(f) Lr(Ω) B(f)B(g)

efq0 (B(g) B(f)) Lr(Ω) B(f)B(g) +

ef eg q0B(f) Lr(Ω) B(f)B(g)

|B(g) B(f)| efq0 Lr(Ω)

(ef eg)q0 Lr(Ω) B(g) . (A.1)

Observe that

|B(f) B(g)| Z

Ω |ef eg|q0 dx = Z

Ω eg|ef g 1|q0 dx e f g f g B(g)

since |eu v 1| |u v|e|u v| for any u, v R. Similarly,

(ef eg)q0 Lr(Ω) e f g f g egq0 Lr(Ω).

Using these above, we obtain

pf pg Lr(Ω) e f g f g

B(f) + egq0 Lr(Ω)

Since efq0 Lr(Ω) e f q0 Lr(Ω) and B(f) e f , from (A.2) we obtain

pf pg Lr(Ω) e f g f g q0 Lr(Ω) e2 f + e2 g .

2e f g f g q0 Lr(Ω)e2 max{ f , g }

2e2 f g f g q0 Lr(Ω)e2 min{ f , g }

where we used max{a, b} min{a, b} + |a b| for a, b 0 in the last line above.

(ii) This simply follows from (A.2) by using r = 1.

A.2 General Representer Theorem

The following is the general representer theorem for abstract Hilbert spaces.

Sriperumbudur et al.

Theorem A.2 (General representer theorem) Let H be a real Hilbert space and let (φi)m i=1 Hm. Suppose J : H R be such that J(f) = V ( f, φ1 H, . . . , f, φm H) , f H where V : Rn R is a convex differentiable function. Define

fλ = arg inf f H J(f) + λ

where λ > 0. Then there exists (αi)m i=1 Rm such that fλ = Pm i=1 αiφi where α := (α1, . . . , αm) satisfies the following (possibly nonlinear) equation

λα + V (Kα) = 0,

with K being a linear map on Rm and (K)i,j = φi, φj H, i [m], j [m].

Proof Define A : H Rm, f 7 ( f, φi H)m i=1. Then fλ = arg inff H V (Af) + λ

2 f 2 H. Therefore, Fermat s rule yields

0 = A V (Afλ) + λfλ fλ = A 1

( α Rm) fλ = A α, α = 1

( α Rm) fλ = A α, α = 1

λ V (AA α),

where A : Rm H is the adjoint of A which can be obtained as follows. Note that

i=1 αi f, φi H =

H ( f H) ( α Rm)

and thus A α = Pm i=1 αiφi. Therefore AA α = Pm j=1 αj Aφj = Pm j=1 αj( φj, φi H)m i=1, and so for every i [m], (AA α)i = Pm j=1 φj, φi Hαj and hence AA = K.

A.3 Bounds on Approximation Errors, A0(λ) and A 1

The following result is quite well-known in the linear inverse problem theory (Engl et al., 1996).

Proposition A.3 Let C be a bounded, self-adjoint compact operator on a separable Hilbert space H. For λ > 0 and f H, define fλ := (C + λI) 1Cf and Aθ(λ) := Cθ(fλ f) H for θ 0. Then the following hold.

(i) For any θ > 0, Aθ(λ) 0 as λ 0 and if f R(C), then A0(λ) 0 as λ 0.

(ii) If f R(Cβ) for β 0 and β + θ > 0, then

Aθ(λ) max{1, C β+θ 1}λmin{1,β+θ} C βf H.

Density Estimation in Infinite Dimensional Exponential Families

Proof (i) Since C is bounded, compact, and self-adjoint, the Hilbert-Schmidt theorem (Reed and Simon, 1972, Theorems VI.16, VI.17) ensures that C = P l αlφl φl, H, where (αl)l N are the positive eigenvalues and (φl)l N are the corresponding unit eigenvectors that form an ONB for R(C). Let θ = 0. Since f R(C),

A2 0(λ) = (C + λI) 1Cf f 2 H =

αi αi + λ f, φi Hφi X

i f, φi Hφi

λ αi + λ f, φi Hφi

2 f, φi 2 H 0 as λ 0

by the dominated convergence theorem. For any θ > 0, we have

A2 θ(λ) = Cθ(C + λI) 1Cf Cθf 2

Let f = f R+f N where f R R(Cθ), f N R(Cθ) if 0 < θ 1 and f R R(C), f N R(C)

if θ 1. Then

A2 θ(λ) = Cθ(C + λI) 1Cf Cθf 2

H = Cθ(C + λI) 1Cf R Cθf R 2

α1+θ i αi + λ f R, φi Hφi X

i αθ i f R, φi Hφi

λαθ i αi + λ f R, φi Hφi

λαθ i αi + λ

2 f R, φi 2 H 0

(ii) If f R(Cβ), then there exists g H such that f = Cβg. This yields

A2 θ(λ) = Cθ(C + λI) 1Cf Cθf 2

H = Cθ(C + λI) 1Cβ+1g Cθ+βg 2

λαθ+β i αi + λ g, φi Hφi

λαθ+β i αi + λ

!2 g, φi 2 H. (A.3)

Suppose 0 < β + θ < 1. Then

αβ+θ i λ αi + λ = αi αi + λ

β+θ λ αi + λ

1 θ β λβ+θ λβ+θ.

On the other hand, for β + θ 1, we have

αβ+θ i λ αi + λ = αi αi + λ

αβ+θ 1 i λ C β+θ 1λ.

Using the above in (A.3) yields the result.

Sriperumbudur et al.

A.4 Bound on the Norm of Certain Operators and Functions

The following result is used in many places throughout the paper. We would like to highlight that special cases of this result are known, e.g., see the proof of Theorem 4 in Caponnetto and Vito (2007) where concentration inequalities are obtained for the quantities in Proposition A.4 using Bernstein s inequality. Here, we provide asymptotic statements using Chebyshev s inequality.

Proposition A.4 Let X be a topological space, H be a separable Hilbert space and L+ 2 (H) be the space of positive, self-adjoint Hilbert-Schmidt operators on H. Define R := R

X r(x) d P(x)

and ˆR := 1 n Pm a=1 r(Xa) where P M1 +(X), (Xa)m a=1 i.i.d. P and r is a L+ 2 (H)-valued measurable function on X satisfying R

X r(x) 2 HS d P(x) < . Define gλ := (R + λI) 1Rg for g H, λ > 0, and A0(λ) := gλ g H. Let α 0 and θ > 0. Then the following hold:

(i) ( ˆR R)(gλ g) H = OP A0(λ) m

(ii) Rα(R + λI) θ λα θ.

(iii) ˆRα( ˆR + λI) θ λα θ.

(iv) (R + λI) θ( ˆR R) = OP q

(v) Rα( ˆR + λI) 1 λα 1 for m c λ2 where is c is a sufficiently large constant that depends on R r(x) 2 HS d P(x) but not on m and λ.

Proof (i) Note that for any f H,

EP ( ˆR R)f 2 H = EP ˆRf 2 H + Rf 2 H 2EP ˆRf, Rf H,

where EP ˆRf, Rf H = 1 n Pn a=1 EP r(Xa)f, Rf H = 1 n Pn a=1 EP r(Xa), f Rf HS. Since R

X r(x) 2 HS d P(x) < , r(x) is P-integrable in the Bochner sense (see Diestel and Uhl, 1977, Definition 1 and Theorem 2), and therefore it follows from Diestel and Uhl (1977, Theorem 6) that EP r(Xa), f Rf HS = R

X r(x) d P(x), f Rf HS = Rf 2 H. Therefore,

EP ( ˆR R)f 2 H = EP ˆRf 2 H Rf 2 H,

EP ˆRf 2 H = EP

a,b=1 EP r(Xa)f, r(Xb)f H.

Splitting the sum into two parts (one with a = b and the other with a = b), it is easy to verify that EP ˆRf 2 H = 1

X r(x)f 2 H d P(x) + m 1

m Rf 2 H, thereby yielding

EP ( ˆR R)f 2 H = 1

X r(x)f 2 H d P(x) Rf 2 H

X r(x)f 2 H d P(x)

X r(x) 2 H d P(x).

Density Estimation in Infinite Dimensional Exponential Families

Using f = gλ g, an application of Chebyshev s inequality yields the result.

(ii, iii) Rα(R + λI) θ = supi γα i (γi+λ)θ = supi h γi γi+λ α 1 (γi+λ)θ α i supi 1 (γi+λ)θ α λα θ, where (γi)i N are the eigenvalues of R. (iii) follows by replacing (γi)i N with the eigenvalues of ˆR.

(iv) Since (R+λI) θ( ˆR R) (R+λI) θ( ˆR R) HS, consider EP (R+λI) θ( ˆR R) 2 HS, which using the technique in the proof of (i), can be shown to be bounded as

EP (R + λI) θ( ˆR R) 2 HS 1

X (R + λI) θr(x) 2 HS d P(x). (A.4)

(R + λI) θr(x) 2 HS = (R + λI) θr(x), (R + λI) θr(x) HS = (R + λI) 2θ Tr (r(x)r(x)) = (R + λI) 2θ r(x) 2 HS λ 2θ r(x) 2 HS, (A.5)

where the last inequality follows from (iii). Using (A.5) in (A.4), we obtain

EP (R + λI) θ( ˆR R) 2 HS 1 mλ2θ

X r(x) 2 HS d P(x).

The result therefore follows by an application of Chebyshev s inequality.

(v) We use the idea in Step 2.1 of the proof of Theorem 4 in Caponnetto and Vito (2007), where Rα( ˆR+λI) 1 is written equivalently as follows: Note that ˆR+λI = ( ˆR R)+(R+λI), which implies

( ˆR + λI) 1 = ( ˆR R) + (R + λI) 1 = (R + λI) 1 I (R ˆR)(R + λI) 1 1

and so Rα( ˆR+λI) 1 = Rα(R+λI) 1 I (R ˆR)(R + λI) 1 1 . Using the Von Neumann series representation, we have

Rα( ˆR + λI) 1 = Rα(R + λI) 1 X

(R ˆR)(R + λI) 1 j

Rα( ˆR + λI) 1 Rα(R + λI) 1

j=0 (R ˆR)(R + λI) 1 j HS

j=0 (R ˆR)(R + λI) 1 j HS.

From the proof of (iv), we have that for any δ > 0, with probability at least 1 δ,

(R ˆR)(R + λI) 1 HS q R

X r(x) 2 HS d P(x) mλ2δ . Suppose m

X r(x) 2 HS d P(x) s2λ2δ where s < 1.

Then P j=0 (R ˆR)(R + λI) 1 j HS P j=0 sj = 1 1 s. This means for m c λ2 where c is sufficiently large, we obtain Rα( ˆR + λI) 1 λα 1.

Sriperumbudur et al.

A.5 Interpolation Space

In this section, we briefly recall the definition of interpolation spaces of the real method. To this end, let E0 and E1 be two arbitrary Banach spaces that are continuously embedded in some topological (Hausdorff) vector space E. Then, for x E0 + E1 := {x0 + x1 : x0 E0, x1 E1} and t > 0, the K-functional of the real interpolation method (see Bennett and Sharpley, 1988, Definition 1.1, p. 293) is defined by

K(x, t, E0, E1) := inf{ x0 E0 + t x1 E1 : x0 E0, x1 E1, x = x0 + x1}.

Suppose E and F are two Banach spaces that satisfy F , E (i.e., F E and the inclusion operator id : F E is continuous), then the K-functional reduces to

K(x, t, E, F) = inf y F x y E + t y F . (A.6)

The K-functional can be used to define interpolation norms, for 0 < θ < 1, 1 s and x E0 + E1, as

( R t θK(x, t) s t 1 dt 1/s , 1 s < supt>0 t θK(x, t), s = .

Moreover, the corresponding interpolation spaces (Bennett and Sharpley, 1988, Definition 1.7, p. 299) are defined as

[E0, E1]θ,s := {x E0 + E1 : x θ,s < } .

B. Appendix: Miscellaneous Results

In this appendix, we present the proofs of some claims that we made in Sections 1, 4 and 5.

B.1 Relation between Fisher and Kullback-Leibler Divergences

The following result provides a relationship between Fisher and Kullback-Leibler divergences.

Proposition B.1 Let p and q be probability densities defined on Rd. Define pt := p N(0, t Id) and qt := q N(0, t Id) where N(0, t Id) denotes a normal distribution on Rd with mean zero and diagonal covariance with t > 0. Suppose pt and qt satisfy

ipt(x) log pt(x) = o ( x α 2 ) , ipt(x) log qt(x) = o ( x α 2 ) and i log qt(x)pt(x) = o ( x α 2 )

as x 2 for all i [d] where α = 1 d. Then

KL(p q) = Z

0 J(pt qt) dt, (B.1)

where J is defined in (3).

Density Estimation in Infinite Dimensional Exponential Families

Proof Under the conditions mentioned on pt and qt, it can be shown that

d dt KL(pt qt) = J(pt qt). (B.2)

See Theorem 1 in Lyu (2009) for a proof. The above identity is a simple generalization of de Bruijn s identity that relates the Fisher information to the derivative of the Shannon entropy (see Cover and Thomas, 1991, Theorem 16.6.2). Integrating w.r.t. t on both sides

of (B.2), we obtain KL(pt qt)

t=0 = R 0 J(pt qt) dt which yields the equality in (B.1) as

KL(pt qt) 0 as t and KL(pt qt) KL(p q) as t 0.

B.2 Estimation of p0: Unbounded k

To handle the case of unbounded k, in the following, we assume that there exists a positive constant M such that f0 H M, so that an estimator of f0 can be constructed as

fλ,n = arg inf f H ˆJλ(f) subject to f H M, (B.3)

where ˆJλ is defined in Theorem 4(iv). This modification yields a valid estimator p fλ,n as

long as k satisfies R

k(x,x)q0(x) dx < , since this implies fλ,n F. The construction of fλ,n requires the knowledge of M, however, which we assume is known a priori. Using the representer theorem in RKHS, it can be shown (see Section B.2.1) that

fλ,n = δˆξ +

j=1 β(b 1)d+j jk(Xb, )

where δ and β are obtained by solving the following quadratically constrained quadratic program (QCQP),

( β, δ) =: Θ = arg min Θ Rnd+1 1 2ΘT HΘ + ΘT subject to ΘT KΘ M2,

with := (h, ˆξ 2 H), Θ := (β, δ) and K, H being defined in the proof of Theorem 5 and the remark following it. The following result investigates the consistency and convergence rates for p fλ,n.

Theorem B.2 (Consistency and rates for p fλ,n) Let M f0 H be a fixed constant,

and fn,λ be a clipped estimator given by (B.3). Suppose (A) (D) with ε = 2 hold. Let

supp(q0) = Ωand R

k(x,x)q0(x) dx < . Define η(x) = p

k(x,x). Then, as λ n , λ 0 and n ,

(i) p fλ,n p0 L1(Ω) 0, KL(p0 p fλ,n) 0 if η L1(Ω, q0);

(ii) for 1 < r , p fλ,n p0 Lr(Ω) 0 if ηq0 L1(Ω) Lr(Ω) and e M

k( , )q0 Lr(Ω);

Sriperumbudur et al.

(iii) h(p fλ,n, p0) 0 if p

k( , )η L1(Ω, q0);

(iv) J(p0 p fλ,n) 0.

In addition, if f0 R(Cβ) for some β > 0, then p fλ,n p0 Lr(Ω) = Op0(θn), h(p0, p fλ,n) =

Op0(θn), KL(p0 p fλ,n) = Op0(θn) and J(p0 p fλ,n) = Op0(θ2 n) where θn := n min n 1 4 , β 2(β+1) o

with λ = n max n 1 4 , 1 2(β+1) o assuming the respective conditions in (i)-(iii) above hold.

Proof For any x Ω, since |f0(x)| f0 H p

k(x, x) M p

k(x, x) and | fλ,n(x)| M p

k(x, x), we have

e fλ,n(x) ef0(x) e M

k(x,x) fλ,n(x) f0(x) η(x) fλ,n f0 H, (B.4)

where we used the fact that |ex ey| ea|x y| for x, y [ a, a] and η(x) := p

k(x,x). In the following, we obtain bounds for p fλ,n p0 Lr(Ω) for any 1 r , h(p fλ,n, p0)

and KL(p0 p fλ,n) in terms of fλ,n f0 H. Define B(f) := R

Ωefq0 dx. Since k satisfies R

k(x,x)q0(x) dx < , then it is clear that fλ,n F as B( fλ,n) < since Z

Ω e fλ,n(x)q0(x) dx Z

k(x,x)q0(x) dx Z

k(x,x)q0(x) dx < .

Similarly, it is easy to verify that B(f0) < .

(i) Recalling (A.1), we have

p fλ,n p0 Lr(Ω) |B( fλ,n) B(f0)| ef0q0 Lr(Ω)

B( fλ,n)B(f0) + (ef0 e fλ,n)q0 Lr(Ω)

If r = 1, we obtain

p fλ,n p0 L1(Ω) |B( fλ,n) B(f0)|

B( fλ,n) + (ef0 e fλ,n)q0 L1(Ω)

Using (B.4), we bound |B( fλ,n) B(f0)| as

|B( fλ,n) B(f0)| Z

e fλ,n(x) ef0(x) q0(x) dx η L1(Ω,q0) fλ,n f0 H.

Also for any f H with f H M, we have B(f) R

k(x,x)q0(x) dx =: θ, where θ > 0. Again using (B.4), we have

(ef0 e fλ,n)q0 Lr(Ω) ηq0 Lr(Ω) fλ,n f0 H

and ef0q0 Lr(Ω) e M

k(x,x)q0 Lr(Ω). Therefore,

p fλ,n p0 Lr(Ω) η L1(Ω,q0) e M

k(x,x)q0 Lr(Ω) fλ,n f0 H θ2

Density Estimation in Infinite Dimensional Exponential Families

+ ηq0 Lr(Ω) fλ,n f0 H

and for r = 1, p fλ,n p0 L1(Ω) 2 η L1(Ω,q0) fλ,n f0 H θ .

KL(p0 p fλ,n) = Z

Ω p0 log p0 p fλ,n dx = Z

ef0 fλ,n B( fλ,n)

f0 fλ,n + log B( fλ,n)

|B( fλ,n) B(f0)|

B(f0) + fλ,n f0 L1(Ω,p0) 2 ηq0 L1(Ω)

θ fλ,n f0 H.

(iii) It is easy to verify that

h(p fλ,n, p0) =

e fλ,n/2 L2(Ω,q0) ef0/2 ef0/2 L2(Ω,q0)

2 e fλ,n/2 ef0/2 L2(Ω,q0) ef0/2 L2(Ω,q0)

where the above inequality is obtained by carrying out and simplifying the decomposition as in (A.1). Using (B.4), we therefore have

h(p fλ,n, p0)

Ωk(x, x)e M

k(x,x)q0 dx θ fλ,n f0 H.

(iv) As f0, fλ,n F, by Theorem 4, we obtain J(p0 p fλ,n) = 1 2

C( fλ,n f0) 2 H

C 2 fλ,n f0 2 H.

Note that we have bounded the various distances between p fλ,n and p0 in terms of fλ,n

f0 H. Since fλ,n = fλ,n with probability converging to 1, the assertions on consistency are proved by Theorem 6(i) in combination with Lemma 14 as we did not explicitly assume f0 R(C) and the rates follow from Theorem 6(iii).

Remark The following observations can be made while comparing the scenarios of using bounded vs. unbounded kernels in the problem of estimating p0 through Theorems 7 and B.2. First, the consistency results in Lr, Hellinger and KL distances are the same but for additional integrability conditions on k and q0. The additional integrability conditions are not too difficult to hold in practice as they involve k and q0 which can be chosen appropriately. However, the unbounded situation in Theorem B.2 requires the knowledge of M which is usually not known. On the other hand, the consistency result in J in Theorem B.2 is slightly weaker than in Theorem 7. This may be an artifact of our analysis as

Sriperumbudur et al.

we are not able to adapt the bounding technique used in the proof of Theorem 7 to bound J(p0 p fλ,n) = 1

C( fλ,n f0) 2 H as it critically depends on the boundedness of k. There-

fore, we used a trivial bound of J(p0 p fλ,n) = 1

C( fλ,n f0) 2 H 1

C 2 fλ,n f0 2 H, which yields the result through Theorem 6(i). Due to the same reason, we also obtain a slower rate of convergence in J. Second, the rate of convergence in KL is slower than in Theorem B.2, which again may be an artifact of our analysis. The convergence rate for KL in Theorem 7 is based on the application of Theorem 6(ii) in Lemma A.1, where the bound on KL in Lemma A.1 critically uses the boundedness to upper bound KL in terms of squared Hellinger distance.

B.2.1 Derivation of fλ,n

Any f H can be decomposed as f = f + f where

f span n ˆξ, ( jk(Xb, ))b,j o =: H ,

which is a closed subset of H and f H := g H : g, h H = 0, h H so that H = H H . Since the objective function in (B.3) matches with the one in Theorem 5,

using the above decomposition in (B.3), it is easy to verify that ˆJ depends only on f H so that (B.3) reduces to

( f λ,n, f λ,n) = arg inf f H ,f H

f 2 H+ f 2 H M2

ˆJλ(f ) + λ

2 f 2 H + λ

2 f 2 H (B.5)

and fλ,n = f λ,n + f λ,n. Since f is of the form in (14), using it in (B.5), it is easy to show

that ˆJλ(f ) + λ

2 f 2 H = 1

2ΘT HΘ + ΘT . Similarly, it can be shown that f 2 H = ΘT KΘ. Since f appears in (B.5) only through f 2 H, (B.5) reduces to

(Θ , c ) = arg inf Θ Rnd+1,c 0 ΘT KΘ+c M2

1 2ΘT HΘ + ΘT + λ

where f λ,n is constructed as in (14) using Θ and f λ,n is such that f λ.n 2 H = c . The necessary and sufficient conditions for the optimality of (Θ , c ) is given by the following Karush-Kuhn-Tucker conditions,

(H + 2τK)Θ + = 0, λ

2 + η τ = 0 (Stationarity)

ΘT KΘ + c M2, c 0 (Primal feasibility)

η 0, τ 0 (Dual feasibility)

τc = 0, η(ΘT KΘ + c M2) = 0 (Complementary slackness)

Combining the dual feasibility and stationary conditions, we have η = τ λ

2 0, i.e., τ λ

2. Using this in the complementary slackness involving τ and c , it follows that c = 0. Since f λ,n 2 = c , we have f λ,n = 0, i.e., fλ,n is completely determined by f λ,n. Therefore f λ,n is of the form in (14) and (B.6) reduces to a quadratically constrained quadratic program.

Density Estimation in Infinite Dimensional Exponential Families

B.3 R(Cβ) and Interpolation Spaces

Proposition B.3 presents an interpretation for R(Cβ) (β > 0 and β / N) as interpolation spaces between R(C β ) and R(C β ) where R(C0) := H. An inspection of its proof shows that Proposition B.3 holds for any self-adjoint, bounded, compact operator defined on a separable Hilbert space.

Proposition B.3 Suppose (B) and (D) hold with ε = 1. Then for all β > 0 and β / N

R(Cβ) = h R(C β ), R(C β ) i

where R(C0) := H, and the spaces R(Cβ) and R(C β ), R(C β )

β β ,2 have equivalent norms.

To prove Proposition B.3, we need the following result which we quote from Steinwart and Scovel (2012, Lemma 6.3) (also see Tartar, 2007, Lemma 23.1) that interpolates L2-spaces whose underlying measures are absolutely continuous with respect to a measure ν.

Lemma B.4 Let ν be a measure on a measurable space Θ and w0 : Θ [0, ) and w1 : Θ [0, ) be measurable functions. For 0 < β < 1, define wβ := w1 β 0 wβ 1 . Then we have [L2(w0 dν), L2(w1 dν)]β,2 = L2(wβ dν)

and the norms on these two spaces are equivalent. Moreover, this result still holds for weights w0 : Θ (0, ) and w1 : Θ [0, ], if one uses the convention 0 := 0 in the definition of the weighted spaces.

Proof of Proposition B.3. The proof is based on the ideas used in the proof of Theorem 4.6 in Steinwart and Scovel (2012). Recall that by the Hilbert-Schmidt theorem, C has the following representation, C = X

i I αiφi φi, H

where (αi)i I are the positive eigenvalues of C, (φi)i I are the corresponding unit eigenvectors that form an ONB for R(C) and I is an index set which is either finite (if H is finite-dimensional) or I = N with limi αi = 0 (if H is infinite dimensional). Let (ψi)i J be an ONB for N(C) where J is some index set so that any f H can be written as

i I f, φi Hφi + X

i J f, ψi Hψi =: X

where θi := φi if i I and θi := ψi if i J with ai := f, θi H. Let β > 0. By definition, g R(Cβ) is equivalent to h H such that g = Cβh, i.e.,

i I αβ i h, φi Hφi =: X

i I biαβ i φi

Sriperumbudur et al.

where bi := h, φi H. Clearly P i I b2 i = P i I h, φi 2 H h 2 H < , i.e., (bi) ℓ2(I). Therefore

i I biαβ i φi : (bi) ℓ2(I)

i I ciφi : (ci) ℓ2(I, α 2β)

where α := (αi)i I. Let us equip this space with the bilinear form *X

i I ciφi, X

R(Cβ) := (ci), (di) ℓ2(I,α 2β)

so that it induces the norm

R(Cβ) := (ci) ℓ2(I,α 2β) .

It is easy to verify that (αβ i φi)i I is an ONB of R(Cβ). Also since R(Cβ1) R(Cβ2) for 0 < β2 < β1 < and id : R(Cβ1) R(Cβ2) is continuous, i.e., for any g R(Cβ1),

g R(Cβ2) = (ci) ℓ2(I,α 2β2) =

c2 i α2β2 i sup i I |αi|β1 β2 (ci) ℓ2(I,α 2β1)

= C β1 β2 g R(Cβ1) <

and so R(Cβ1) , R(Cβ2). Similarly, we can show that R(C) , H. In the following, we first prove the result for 0 < β < 1 and then for β > 1.

(a) 0 < β < 1: For any f H and g R(C), we have

i I J aiθi X

i I J (ai ci)θi

H = (ai ci) 2 ℓ2(I J)

where we define ci := 0 for i J. For t > 0, we find

K(f, t, H, R(C)) = inf g R(C) f g H + t g R(C)

= inf (ci) ℓ2(I,α 2) (ai ci) ℓ2(I J) + t (ci) ℓ2(I,α 2)

= K(a, t, ℓ2(I J), ℓ2(I, α 2)).

From this we immediately obtain the equivalence

f [H, R(C)]β,2 (ai) [ℓ2(I J), ℓ2(I, α 2)]β,2

where 0 < β < 1. Applying the second part of Lemma B.4 to the counting measure on I J yields [ℓ2(I J), ℓ2(I, α 2)]β,2 = ℓ2(I, α 2β).

Density Estimation in Infinite Dimensional Exponential Families

Since R(Cβ) and ℓ2(I, α 2β) are isometrically isomorphic, we obtain R(Cβ) = [H, R(C)]β,2.

(b) β > 1 and β / N: Define γ := β . Let f R(Cγ) and g R(Cγ+1), i.e., (ci) ℓ2(I, α 2γ) and (di) ℓ2(I, α 2γ 2) such that f = P i I ciφi and g = P i I diφi. Since

f g 2 R(Cγ) = (ci di) 2 ℓ2(I,α 2γ),

for t > 0, we have

K(f, t, R(Cγ), R(Cγ+1)) = inf g R(Cγ+1) f g R(Cγ) + t g R(Cγ+1)

= inf (di) ℓ2(I,α 2γ 2) (ci di) ℓ2(I,α 2γ) + t (di) ℓ2(I,α 2γ 2)

= K(c, t, ℓ2(I, α 2γ), ℓ2(I, α 2γ 2)),

from which we obtain the following equivalence

f [R(Cγ), R(Cγ+1)]β γ,2 (ci) [ℓ2(I, α 2γ), ℓ2(I, α 2γ 2)]β γ,2 ( ) = ℓ2(I, α 2β),

where ( ) follows from Lemma B.4 and the result is obtained by noting that ℓ2(I, α 2β) and R(Cβ) are isometrically isomorphic.

B.4 Denseness of Ik H in W2(Rd, p)

In this section, we discuss the denseness of Ik H in W2(Rd, p) for a given p PFD, where PFD is defined in Theorem 13, which is equivalent to the injectivity of Sk (see Rudin, 1991, Theorem 4.12). To this end, in the following result we show that under certain conditions on a bounded continuous translation invariant kernel on Rd, the restriction of Sk to W 2 (Rd, p) is injective when d = 1, while the result for any general d > 1 is open. However, even for d = 1, this does not guarantee the injectivity of Sk (which is defined on W2(Rd, p)). Therefore, the question of characterizing the injectivity of Sk (or equivalently the denseness of Ik H in W2(Rd, p)) is open.

Proposition B.5 Suppose k(x, y) = ψ(x y), x, y Rd where ψ Cb(Rd) L1(Rd), R ω 2ψ (ω) dω < and supp(ψ ) = Rd. If d = 1, then the restriction of Sk to W 2 (Rd, p) is injective for any p PFD.

Proof Fix any p PFD. We need to show that for [f] W 2 (Rd, p), Sk[f] = 0 implies [f] = 0. From Proposition 11, we have

j=1 jk(x, ) jf(x) p(x) dx = Z

j=1 jψ(x ) jf(x) p(x) dx

Rd iωjψ (ω)ei ω, x dω jf(x) p(x) dx

j=1 φj(ω)ψ (ω)ei ω, dω

Sriperumbudur et al.

where φj(ω) := 1 (2π)d/2

Rd(iωj)e i ω,x jf(x) p(x) dx.

Skf = 0 implies Pd j=1 φj(ω)ψ (ω) = 0 for all ω Rd. Since supp(ψ ) = Rd, we have Pd j=1 φj(ω) = 0 a.e., i.e., for ω-a.e.,

Rd(iωj) jf(x) p(x) e i ω,x dx =

j=1 (iωj)(p jf) (ω).

For d = 1, this implies ( jf)p = 0 a.e. and so f W2 = 0.

Examples of kernels that satisfy the conditions in Proposition B.5 include the Gaussian, Mat ern (with β > 1) and inverse multiquadrics on R.

N. Aronszajn. Theory of reproducing kernels. Trans. Amer. Math. Soc., 68:337 404, 1950.

A. Barron and C-H. Sheu. Approximation of density functions by sequences of exponential families. Annals of Statistics, 19(3):1347 1369, 1991.

F. Bauer, S. Pereverzev, and L. Rosasco. On regularization algorithms in learning theory. Journal of Complexity, 23(1):52 72, 2007.

C. Bennett and R. Sharpley. Interpolation of Operators. Academic Press, Boston, 1988.

L. D. Brown. Fundamentals of Statistical Exponential Families with Applications in Statistical Decision Theory. IMS, Hayward, CA, 1986.

S. Canu and A. J. Smola. Kernel methods and the exponential family. Neurocomputing, 69: 714 720, 2005.

A. Caponnetto and E. De Vito. Optimal rates for regularized least-squares algorithm. Foundations of Computational Mathematics, 7:331 368, 2007.

D. Comaniciu and P. Meer. Mean shift: A robust approach toward feature space analysis. IEEE Transactions on Pattern Analysis and Machine Intelligence, 24(5):603 619, 2002.

T. M. Cover and J. A. Thomas. Elements of Information Theory. John Wiley & Sons, Inc., New York, 1991.

E. De Vito, L. Rosasco, and A. Toigo. Learning sets with separating kernels. http://arxiv.org/abs/1204.3573, April 2012.

R. De Vore, G. Kerkyacharian, D. Picard, and V. Temlyakov. Mathematical methods for supervised learning. IMI Preprints, University of South Carolina, 2004.

J. Diestel and J. J. Uhl. Vector Measures. American Mathematical Society, Providence, 1977.

Density Estimation in Infinite Dimensional Exponential Families

A. Doucet, N. de Freitas, and N. J. Gordon, editors. Sequential Monte Carlo Methods in Practice. Springer-Verlag, New York, 2001.

J. J. Duistermaat and J. A. C. Kolk. Multidimensional Real Analysis II: Integration. Cambridge University Press, Cambridge, UK, 2004.

H. W. Engl, M. Hanke, and A. Neubauer. Regularization of Inverse Problems. Kluwer Academic Publishers, Dordrecht, The Netherlands, 1996.

G. B. Folland. Real Analysis: Modern Techniques and Their Applications. Wiley Interscience, New York, 1999.

K. Fukumizu. Exponential manifold by reproducing kernel Hilbert spaces. In P. Gibilisco, E. Riccomagno, M.-P. Rogantin, and H. Winn, editors, Algebraic and Geometric mothods in Statistics, pages 291 306. Cambridge University Press, 2009.

K. Fukumizu, F. Bach, and M. Jordan. Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces. Journal of Machine Learning Research, 5:73 99, 2004.

K. Fukumizu, A. Gretton, X. Sun, and B. Sch olkopf. Kernel measures of conditional dependence. In J.C. Platt, D. Koller, Y. Singer, and S. Roweis, editors, Advances in Neural Information Processing Systems 20, pages 489 496, Cambridge, MA, 2008. MIT Press.

K. Fukumizu, F. Bach, and M. Jordan. Kernel dimension reduction in regression. Annals of Statistics, 37(4):1871 1905, 2009.

K. Fukumizu, L. Song, and A. Gretton. Kernel Bayes rule: Bayesian inference with positive definite kernels. Journal of Machine Learning Research, 14:3753 3783, 2013.

A. Gretton, K. Borgwardt, M. Rasch, B. Sch olkopf, and A. Smola. A kernel two-sample test. Journal of Machine Learning Research, 13:723 773, 2012.

C. Gu and C. Qiu. Smoothing spline density estimation: Theory. Annals of Statistics, 21 (1):217 234, 1993.

A. Hyv arinen. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6:695 709, 2005.

A. Hyv arinen. Some extensions of score matching. Computational Statistics & Data Analysis, 51:2499 2512, 2007.

O. Johnson. Information Theory and The Central Limit Theorem. Imperial College Press, London, UK, 2004.

G. S. Kimeldorf and G. Wahba. Some results on Tchebycheffian spline functions. Journal of Mathematical Analysis and Applications, 33:82 95, 1971.

C. Ley and Y. Swan. Stein s density approach and information inequalities. Electronic Communications in Probability, 18:1 14, 2013.

Sriperumbudur et al.

S. Lyu. Interpretation and generalization of score matching. In Proc. 25th conference on Uncertainty in Artificial Intelligence, Montreal, Canada, 2009.

G. Pistone and C. Sempi. An infinite-dimensional geometric structure on the space of all the probability measures equivalent to a given one. Annals of Statistics, 23(5):1543 1561, 1995.

M. Reed and B. Simon. Functional Analysis. Academic Press, New York, 1972.

L. Rosasco, E. De Vito, and A. Verri. Spectral methods for regularization in learning theory. Technical Report DISI-TR-05-18, DISI, Universit a degli Studi di Genova, 2005.

W. Rudin. Functional Analysis. Mc Graw-Hill, USA, 1991.

H. Sasaki, A. Hyv arinen, and M. Sugiyama. Clustering via mode seeking by direct estimation of the gradient of a log-density. In Proc. of European Conference on Machine Learning, pages 19 34, 2014.

B. Sch olkopf, R. Herbrich, and A. Smola. A generalized representer theorem. In Proceedings of the 14th Annual Conference on Computational Learning Theory and 5th European Conference on Computational Learning Theory, pages 416 426, London, UK, 2001. Springer-Verlag.

S. Smale and D.-X. Zhou. Learning theory estimates via integral operators and their approximations. Constructive Approximation, 26:153 172, 2007.

B. K. Sriperumbudur, A. Gretton, K. Fukumizu, B. Sch olkopf, and G. R. G. Lanckriet. Hilbert space embeddings and metrics on probability measures. Journal of Machine Learning Research, 11:1517 1561, 2010.

B. K. Sriperumbudur, K. Fukumizu, and G. R. G. Lanckriet. Universality, characteristic kernels and RKHS embedding of measures. Journal of Machine Learning Research, 12: 2389 2410, 2011.

I. Steinwart and A. Christmann. Support Vector Machines. Springer, New York, 2008.

I. Steinwart and C. Scovel. Mercer s theorem on general domains: On the interaction between measures, kernels and RKHSs. Constructive Approximation, 35:363 417, 2012.

H. Strathmann, D. Sejdinovic, S. Livingstone, Z. Szab o, and A. Gretton. Gradient-free Hamiltonian Monte Carlo with efficient kernel exponential families. In C. Cortes, N.D. Lawrence, D.D. Lee, M. Sugiyama, R. Garnett, and R. Garnett, editors, Advances in Neural Information Processing Systems 28, pages 955 963. Curran Associates, Inc., 2015.

S. Sun, M. Kolar, and J. Xu. Learning structured densities via infinite dimensional exponential families. In C. Cortes, N.D. Lawrence, D.D. Lee, M. Sugiyama, R. Garnett, and R. Garnett, editors, Advances in Neural Information Processing Systems 28, pages 2278 2286. Curran Associates, Inc., 2015.

L. Tartar. An Introduction to Sobolev Spaces and Interpolation Spaces. Springer-Verlag, Berlin, 2007.

Density Estimation in Infinite Dimensional Exponential Families

A. B. Tsybakov. Introduction to Nonparametric Estimation. Springer, New York, 2009.

S. van de Geer. Empirical Processes in M-Estimation. Cambridge University Press, Cambridge, UK, 2000.

A. W. van der Vaart and J. H. van Zanten. Rates of contraction of posterior distributions based on Gaussian process priors. Annals of Statistics, 36(3):1435 1463, 2008.

L. Wasserman. All of Nonparametric Statistics. Springer, 2006.

H. Wendland. Scattered Data Approximation. Cambridge University Press, Cambridge, UK, 2005.

D.-X. Zhou. Derivative reproducing properties for kernel methods in learning theory. Journal of Computational and Applied Mathematics, 220:456 463, 2008.