{'fraction_non_alphanumeric': 0.05678995657536793, 'fraction_numerical': 0.011692546846222732, 'mean_word_length': 4.034793704806465, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 26, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Information Theory with Kernel Methods

We consider the analysis of probability distributions through their associated covariance operators from reproducing kernel Hilbert spaces. We show that the von Neumann entropy and relative entropy of these operators are intimately related to the usual notions of Shannon entropy and relative entropy, and share many of their properties. They come together with efficient estimation algorithms from various oracles on the probability distributions. We also consider product spaces and show that for tensor product kernels, we can define notions of mutual information and joint entropies, which can then characterize independence perfectly, but only partially conditional independence. We finally show how these new notions of relative entropy lead to new upper-bounds on log partition functions, that can be used together with convex optimization within variational inference methods, providing a new family of probabilistic inference methods.


Introduction
Characterizing and studying probability distributions through moments has a long history. For distributions supported in a vector space, only an infinite number of moments can characterize the distribution. Beyond traditional polynomial moments, kernel methods based on reproducing kernel Hilbert spaces (RKHS) [48,10] have emerged as a natural tool for studying the interactions between moments and other properties of the underlying distributions (such as independence or conditional independence), as well as providing algorithmic tools to estimate these moments from data.
A natural way of using kernel methods is to consider a feature map ϕ : X → H from the underlying space X and a specific Hilbert space H, and to consider the mean element [53], also referred to as the mean embedding, for a probability distribution p, defined as which is an element of the Hilbert space H. The key idea behind kernel methods is to study properties of H and its interaction with X only through the kernel function k : X × X → R defined as k(x, y) = ϕ(x), ϕ(y) .
Under classical universality conditions on the kernel function k [53,34], the mapping p → µ p is injective (and thus characterizes the distribution p), and can be estimated from independent and identically distributed (i.i.d.) samples efficiently, with convergence rates proportional to 1/ √ n where n is the number of observations.
Representing a probability distribution p through µ p ∈ H defines explicitly a metric between distributions as (p, q) → µ p − µ q (for the H-norm), which has been extensively used in machine learning and data science [35], in particular for measuring the dependence between two random variables, leading to algorithms for independent component analysis [17,6] and for feature selection [51], or as a model fitting criterion when estimating parameters of probability distributions [13].
A classical drawback of the underlying geometry is the lack of straightforward connection with classical information theory tools. For example for discrete data, this leads to Euclidean norms between probability mass functions, which is rarely seen as an appropriate geometry for the simplex.
In this paper, we consider the second-order moment, which we call the covariance operator 1 which is an operator from H to H, defined through f, Σ p g = f, ϕ(x) g, ϕ(x) dp(x), self-adjoint and positive semi-definite. 2 As we show in this paper, it shares similar nice properties with the mean element in terms of universality, estimation from finite samples, and general applicability to all sets X where positive definite kernels can be defined (e.g., structured discrete objects). Furthermore, owing to the many tools from quantum information theory, and in particular the von Neumann entropy tr Σ p log Σ p and its associated divergence (relative entropy), we define these notions and demonstrate their many nice properties and applications, as well as explicit quantitative links with regular Shannon entropies (in fact, we can see the regular notion of entropy as the limit of our kernel entropies when the bandwidth of the kernel goes to zero, and for discrete data, they are equivalent).
The paper is organized as follows: • We review in Section 2 properties of covariance operators and reproducing kernel Hilbert spaces.
• In Section 3, we review some of the key mathematical results on quantum relative entropies, that we will leverage in further sections.
• In Section 4, we define our notions of kernel entropy and Kullback-Leibler divergence (relative entropy), and give their main properties, which are shared with the classical notions, such as convexity, or additivity for independent variables. In particular our notions of relative entropies will always be lower bounds on the Shannon relative entropy (in this paper, we use interchangeably the terms "relative entropy" and "Kullback-Leibler (KL) divergence").
• In Section 5, we show how we can leverage existing work in kernel methods to estimate these quantities from finite samples, with a convergence rate proportional to 1/ √ n from n observations, and without the need for regularization. We also propose statistically more efficient estimators if integrals of kernel functions under p can be computed.
• In Section 6, we consider product spaces X = X 1 × · · · × X d , and show that for tensor product kernels, we can define notions of mutual information and joint entropies. We can then characterize independence perfectly, but only partially conditional independence, and extend the submodularity property of the regular entropy.
• In Section 7, we show how these new notions of relative entropy lead to new upper-bounds on log partition functions, that can be used together with convex optimization within variational inference methods [58], providing a new family of probabilistic inference methods. Illustrative applications to [0, 1] and {−1, 1} d are explicitly developed.

Covariance operators
In this section, we review the notions of covariance operators, starting with reproducing kernel Hilbert spaces. For more details, see [10,48,49].

Kernels and RKHSs
In this paper we consider a compact set X and probability distributions on X (which we assume to be Borel probability measures), as well as a kernel function k : X × X → R, such that (A1) k is a continuous positive definite kernel on the compact set X, with k(x, x) = 1 for all x ∈ X.
Positive definite kernel functions are functions for which all matrices of pairwise kernel evaluations are positive semi-definite. They are known to define a Hilbert space H (with dot product ·, · ) of functions f : X → R, called a reproducing kernel Hilbert space (RKHS), as well as a map ϕ : X → H such that, for all x, y ∈ X, (a) f (x) = f, ϕ(x) , (b) k(·, x) = ϕ(x), and (c) k(x, y) = ϕ(x), ϕ(y) .
A useful representation of functions in H is through linear combinations of kernel functions, that is, for any α ∈ R n and arbitrary elements x 1 , . . . , x n of X, the function f = n i=1 α i k(·, x i ) is in H and its norm is equal to n i,j=1 α i α j k(x i , x j ). It turns out that all functions in H can be generated as appropriate limits of such finite linear combinations.
Note that from k(x, x) = 1, we get ϕ(x) = 1 for all x ∈ X. Note that most of the developments in this paper do not require further structure on the set X, except when characterizing regularity of density functions is needed in Section 4.2, where will assume the existence of a distance.
Universal and characteristic kernels. We note that our covariance operator Σ p happens to be a mean embedding for the feature map x → ϕ(x)ϕ(x) * , from X to the space of operators from H to H. Equipped with the Hilbert-Schmidt norm, we obtain a Hilbert space which happens to be isomorphic 3 to the RKHS obtained from the kernel (x, y) → k(x, y) 2 .
Therefore, injectivity of the map p → Σ p can be obtained from sufficient conditions for injectivity of the mean element map, a property for the kernel referred to as characteristic [52]. As shown in [26,Prop. 5], this is equivalent for the associated RKHS to be dense in L 2 (q) for all probability measures q on X, a sufficient condition being here (because of our compactness assumption) that the RKHS is dense in the set of continuous functions equipped with the uniform norm, a property referred to as universality [54]. Note that if (A1) is satisfied, then the universality of k implies the universality of k 2 .
We will see in Section 2.3 below that many kernels are indeed universal, and thus the covariance operator will be characteristic of the corresponding probability distribution.

Covariance operators
Given a probability distribution p on X, we consider the operator Σ p : H → H, defined as For any f, g ∈ H, we have, by definition: g, Proposition 1 (Properties of covariance operators) Assume (A1).
(a) For any probability distribution p on X, the operator Σ p : H → H defined in Eq. (2), is self-adjoint, positive-semidefinite and has unit trace.
(b) If the kernel k 2 is universal, then the mapping p → Σ p is injective from probability distributions to self-adjoint positive-semidefinite, unit trace operators.
(c) If p has full support in X, then the operator Σ p has a trivial null space.
(d) Given the mean element defined in Eq. (1), we have Σ p µ p ⊗ µ p .
(b) If Σ p = Σ q , then for all y ∈ X, X k(x, y) 2 [dp(x) − dq(x)] = 0. The representation of functions in the RKHS associated with k 2 as linear functions of x → k(x, y) 2 , and the universality of k 2 implies that Since the operators are trace-class, they have a discrete spectrum [41,16], composed of a summable sequence of strictly positive eigenvalues tending to zero, and potentially zero (if p does not have full support in X).
The spectral decay is a key quantity for the study of kernel methods (see [21] and references therein), as used in Section 5 to estimate entropies from data.
Throughout this paper, we will consider spectral functions of trace-class self-adjoint operators A, such as log A, exp(A), or more generally f (A), for a function f : R → R. These are defined through the eigenvalue decomposition by keeping eigenvectors and replacing eigenvalues λ by f (λ) (for polynomials, we then recover the usual matrix polynomials).
Uniform distribution and symmetric sets.
In order to define the notion of kernel entropy from kernel relative entropy, we will need to use a specific "base" distribution τ , that will play the role of the uniform distribution when this notion is classically defined (e.g., X is a finite set or a subset of R d , like in all examples in Section 2.3).
For thi base distribution τ on X, we use the notation (instead of Σ τ ): In this paper, some concepts will be easier for symmetric sets, that is, sets where there exists a set T of bijective transformations t : X → X, such that for all x, y ∈ X, ∃t ∈ T such that t(x) = y. We will then assume that the kernel k and the base distribution τ are compatible with this notion of symmetry, that is, τ (A) = τ (t(A)) for any Borel set A and t ∈ T, and k(x, y) = k(t(x), t(y)) for all transformations t ∈ T.
Then the covariance operator for the base distribution will also exhibit symmetries. More precisely, all eigensubspaces will be invariant by symmetry, that is, if e 1 , . . . , e m ∈ H is an orthonormal basis of an eigensubspace corresponding to an eigenvalue λ, then 4 m i=1 e i (t(x)) 2 = m i=1 e i (x) 2 for all transformations t ∈ T and x ∈ X. This implies that the function x → m i=1 e i (t(x)) 2 is constant, equal to its expectation under τ (which is itself equal to mλ).
As a consequence, for a spectral function Σ → f (Σ) based on a real-valued function f : R → R, the dotproduct ϕ(x), f (Σ)ϕ(x) can be expanded as the sum of the contributions coming from all eigensubspaces, with the one associated to the eigenvalue λ leading to f (λ) m i=1 e i (x) 2 being independent of x (and equal to mλf (λ)). This leads notably to ϕ(x), f (Σ)ϕ(x) = tr[Σf (Σ)] for all real-valued functions f defined on R, and to tr[Σf (Σ)] = tr[Σ p f (Σ)] for any probability distributions p.
Surjectivity. By definition, every covariance operator is in the closure M of the convex hull of all ϕ(x)ϕ(x) * , for x ∈ X. We will also need the closure A of the span of these vectors. By definition, the mapping p → Σ p is surjective from probability distributions to M. Moreover, if the kernel k is universal, For orthonormal embeddings, where k(x, y) = 0 if x = y, M is (in the appropriate basis) the set of non-negative diagonal operators, while A is the set of diagonal operators (in the same basis).
Note that in most cases A is strictly included in the set of trace class operators from H to H, which is equivalent to the existence of bounded self-adjoint operators M : H → H such that ∀x ∈ X, ϕ(x), M ϕ(x) = 0. Since ϕ(x) 2 = 1 for all x ∈ X, this is equivalent to the existence of a positive bounded operator M so that ∀x ∈ X, ϕ(x), M ϕ(x) = 1, which, with an eigenvalue decomposition leads to a finite or countable family of functions (f i ) i∈I in H such that i∈I f i (x) 2 = 1, that is, a "partition of unity", which are well known to exist in many spaces (see, e.g., [9,Chapter 3] or [29, Section 1.4]).

Examples
Finite set. If X = {x 1 , . . . , x m } is a finite set of cardinality m, then the Gram matrix K ∈ R m×m defined as K ij = ϕ(x i ), ϕ(x j ) completely characterizes the kernel notions up to rotation. If K = I (orthonormal embedding), then there exist m orthonormal vectors u 1 , . . . , u m such that Σ p = m i=1 P(X = i)u i u * i , and thus, all covariance matrices commute and therefore share the same eigenbasis, and we will recover exactly the regular notions of discrete entropy.
We get a symmetric set if K is invariant by permutations, that is, K = αI + (1 − α)11 ⊤ , and as soon as K is invertible, a universal kernel. 6 Polynomial kernels. When X is a subset of R d , then we can consider kernels of the form k(x, y) = (α + β x ⊤ y) r , which correspond to ϕ(x) composed of all monomials of order up to r, and the covariance operator is thus composed of traditional moments. Here, we can have k(x, x) = 1 for all x ∈ X only if X is included in a centered sphere. Alternatively, we may only impose k(x, x) 1 for all x ∈ X, which still preserves many of our results in later sections.
Translation-invariant kernels on a compact subset of R d . These are the usual kernels k(x, y) of the form k(x, y) = q(x − y) where q is a function with strictly positive Fourier transform, which are all universal [54]. Classical examples are the Gaussian kernels k(x, y) = exp(− x − y 2 2 /σ 2 ), where H is a space of infinitely differentiable functions, or the exponential kernel k(x, y) = exp(− x − y 2 /σ), where H is the Sobolev space of functions with square-integrable derivatives up to order (d + 1)/2. However, they typically do not respect the symmetry of the underlying set (for example for spheres, like below).
Torus. If X = [0, 1], we can consider kernels of the form k(x, y) = k(x − y) where k : R → R is 1-periodic and has a strictly positive Fourier series (note that we use the same notation for the kernel and the underlying function). We then get a universal kernel and a symmetric set (invariant by periodic 5 For the non trivial inclusion, we consider A = X ϕ(x)ϕ(x) * dp(x) with X dp(x) = 1, such that A 0. We have for all f ∈ H, X f (x) 2 dp(x) 0. If k is universal, then all positive continuous functions g are limits of functions f 2 , so that X g(x)dp(x) 0, which implies that p is a probability measure. 6 Note that K invertible implies that the Hadamard (i.e., pointwise) product K • K is invertible, but not vice-versa. translation); moreover, the base distribution will be the usual uniform distribution. For example, we will often consider: for which the Fourier series isk(ω) = This extends naturally to [0, 1] d , by considering tensor products, that is, k(x, y) = d i=1 k(x i − y i ), which for example leads tok(ω) = tanh d ( σ 2 )e −σ ω 1 for ω ∈ Z d in the example above. We will also need an explicit formula for ω∈Zk (ω) logk(ω) = log tanh σ 2 − σ sinh σ .
Spheres. We can also consider X as the unit ℓ 2 -sphere in dimension d, with kernels that are invariant by rotation, leading also to symmetric sets. Kernels are then defined through their expansion in spherical harmonics [50], with links with neural networks [4].
We could also imagine extensions to Stiefel or Grassman manifolds [23], and of course to all sets X where positive definite kernels can be defined (such as graphs, trees, permutations, sets, etc.) [49].

Review of quantum information theoretic quantities
We consider here complex-valued Hilbert spaces, but all results will apply to operators on real Hilbert spaces. In quantum mechanics (see, e.g., [30]), physical states are represented as elements of a Hilbert space H, and a physical system is characterized by a convex combination of orthogonal projections on these states, which is often called a density matrix. Our covariance operators defined as Σ p = X ϕ(x)ϕ(x) * dp(x) can thus be seen as "density operators". It turns out that in quantum information theory (see, e.g., [59]), information-theoretic quantities are defined based on these density operators. In this section, we review these notions from a mathematical point of view.
Entropy. For a compact positive self-adjoint operator A on the complex Hilbert space H with finite trace, we can define the negative entropy as [2]: where Λ(A) is the set of eigenvalues of A (with the convention that 0 log 0 = 0). It may not be finite if the sequence of eigenvalues is not decreasing sufficiently fast (note that this is unlikely since they are summable). As defined, tr[A log A] is always less than (tr A) log(tr A) and can be equal to −∞.
Kullback-Leibler (KL) divergence (relative entropy). It is defined as The last property is classical in quantum information theory as the mapping A → n i=1 C i AC * i is usually referred to as a "quantum operation". Note that quantum operations are usually defined as "completely positive" trace-preserving maps, and that owing to Stinespring representation theorem (see [12, Theorem 3.1.2]), they are essentially all of that form.
One particular application of monotonicity is to consider n positive self-adjoint operators D i : H → H, i = 1, . . . , n, such that n i=1 D i = I, and consider µ, ν ∈ R n defined through µ i = tr(D i A) and ν i = tr(D i B) for A, B a positive self-adjoint operator with unit trace, then µ and ν are on the simplex, and the monotonicity leads to (see detailed proof in Appendix A.2): Further results. In Appendix A.1, we provide finer results beyond expectation with respect to a probability measure with finite support, as well as equality cases.
Integral representation. For the KL divergence, we have, by direct integration [1]: which will be used in Section 5.2 for the estimation from a finite sample. As shown in Appendix A, we also have another integral representation formula, that helps in proving additional convexity results.

Kernel entropy and Kullback-Leibler divergence
We define the kernel Kullback-Leibler (KL) divergence, or kernel relative entropy, between two probability distribution p and q as In this section, we list some of its properties, which are similar to the properties of (p, q) → µ p − µ q , but with the added benefit that is has a direct link with Shannon entropy, which is useful in itself when a precise information measure is needed (e.g., for differential privacy [22]), or when used within variational inference (see Section 7).
First, the kernel relative entropy is finite under general conditions.
Proposition 3 (Finiteness of kernel KL divergence) Assume (A1). If p is absolutely continuous with respect to q, with dp dq ∞ α, then both the regular and kernel KL divergences are between zero and log α.
Proof We have then Σ p αΣ q , leading to (note that we must have α 1): The positivity follows from Prop. 2, property (b).
This implies that if p has a bounded density with respect to the base measure τ on X, then the relative entropy D(Σ p Σ) is well-defined, and always non-negative (recall that Σ is the covariance operator for the base probability measure τ on X).
We can then define the kernel entropy as Note that we add terms on top of − tr Σ p log Σ p , one linear term in Σ p , and one constant term, to extend properties of the regular Shannon discrete entropy in all cases. However, everything simplifies in some cases. Indeed, when there are some invariances that are shared by the kernel and the base distribution τ (that is, a symmetric set as defined in Section 2.2), x → ϕ(x), (log Σ)ϕ(x) is constant on X, and thus equal to its expectation under p, which implies min x∈X ϕ(x), (log Σ)ϕ(x) = tr[Σ p log Σ], and then the definition above simplifies to H(Σ p ) = − tr Σ p log Σ p .
The kernel information quantities have natural properties mimicking the traditional quantities, as well as a bound by the regular KL divergence.
(a) The kernel KL divergence defined in Eq. (6) is non-negative, and equal to 0 for p = q. If the kernel k 2 is universal, it is zero if and only if p = q.
(b) The kernel entropy defined in Eq. (7) is always between 0 and − min x∈X ϕ(x), (log Σ)ϕ(x) (which is equal to − tr Σ log Σ for symmetric sets). If p = τ (base measure), it is equal to the upper-bound, while for Dirac measures with a point mass at any minimizers of x → ϕ(x), (log Σ)ϕ(x) , it is equal to 0 (for symmetric sets, this is all of X). If the kernel k 2 is universal, these conditions are also necessary.
(c) The mapping (p, q) → D(Σ p Σ q ) is convex, and thus the mapping p → H(Σ p ) is concave.
(d) For all probability measures p and q, D(Σ p Σ q ) D(p q). Moreover, the differential entropy with respect to the base measure τ , that is, − X log dp dτ (x)dp(x), is less than H(Σ p )+min x∈X ϕ(x), (log Σ)ϕ(x) .
(e) For all probability measures p and q, D( where · * denotes the nuclear norm and · HS the Hilbert-Schmidt norm. Proof Property (a) is the consequence of the non-negativity of the relative entropy when the two arguments have the same trace (here equal to 1), and injectivity of the map p → Σ p when k 2 is universal.
For property (b), the upper-bound results follow similarly. For the non-negativity, we notice that H(Σ p ) defined in Eq. (7) is the sum of two non-negative terms − tr[Σ p log Σ p ] (negative because all eigenvalues of Σ p are less than one) and tr Σ p log Σ − min x∈X ϕ(x), (log Σ)ϕ(x) . It is equal to zero if and only the two terms are equal to zero. For the first term, this is equivalent to Σ p being rank one, while for the second one, this is equivalent to p being supported in the set of minimizers. The sufficient condition for being zero thus follows; for the necessary condition, if k 2 is universal, k(x, x ′ ) 2 < 1 if x = x ′ , and thus Σ p being rank one implies that p is a Dirac measure. Property (c) is a consequence of the convexity of the map (A, B) → tr A(log A−log B) and of the linearity of the mapping p → Σ p . Property (d) is obtained from using joint convexity of the quantum relative entropy, and the fact that ϕ(x) = 1 for all x ∈ X. Note that we can only have equality if for all x = y, k(x, y) = 0, that is, we have an orthonormal embedding, which is only possible for discrete X. Note also that the upper-bound remains valid as soon as k(x, x) 1 for all x ∈ X.
Property (e) is simply a consequence of Pinsker inequality for the quantum relative entropy (see [62]).
Note that while the relative entropy results do not depend on the choice of the base measure, the notion of kernel entropy depends on this choice (however, for the symmetric sets presented in Section 2.3, the natural choice of the base measure is the uniform measure).
As detailed in Appendix A.4, the properties outlined in Prop. 4, in particular the lower-bound on the relative entropy, is valid for other "f -divergences" [19,55], such as the squared Hellinger distance, or the Pearson χ 2 -divergence.

Lower bound on kernel relative entropy
Beyond the upper-bound by the Shannon relative entropy (property (d) in Prop. 4), we can also get a lower bound on the kernel KL divergence, by defining for all y ∈ X, D(y) = Σ −1/2 ϕ(y)ϕ(y) * Σ −1/2 , which is a positive self-adjoint operator, so that X D(y)dτ (y) = I. We can then apply the monotonicity of the relative entropy with operators D(y) 1/2 , y ∈ X, that is, Eq. (4) extended to general expectations. We have which is non negative and such that The function h can thus be seen as a smoothing kernel 7 . Thus, by monotonicity of the relative entropy (Appendix A.2), we get: h(x, y)dp(x) is a smoothed version of p (with the same definition forq). Note that we have directly D(p q) D(p q) from standard Markov chain arguments in information theory 8 [18, Section 2.9]. Overall, we have the following sequence of inequalities: These can lead to quantitative bounds between D(Σ p Σ q ) and D(p q), in particular when the smoothing function h is putting most of its mass on pairs (x, y) where x is close to y. This is made explicit below for distributions on a metric space and distributions with Lipschitz-continuous densities.

Small-width asymptotics for metric spaces
We now provide a bound between D(p q) and D(p q) in Eq. (8), and thus between D(Σ p Σ q ) and D(p q), when the set X is equipped with a distance d, used to characterize the regularity of density functions. 9 This corresponds to bounding the difference in the classical data processing inequality for the Shannon entropy, using regularity properties of the densities of p and q. The following proposition is shown in Appendix B.
Proposition 5 (Joint bound on relative entropy) We assume that p, q have strictly positive Lipschitzcontinuous densities (denoted p and q) with respect to the base measure τ . Using definitions from Section 4.1, we have: Thus, for Lipschitz-continuous densities, the approximation of D(p q) by D(Σ p Σ q ) is controlled by a kernel dependent quantity, that we now study in special cases.
Discrete orthonormal embeddings. For these kernels where k(x, y) = 1 x=y , then we have Σ = I and thus h(x, y) = 1 x=y ; therefore the quantity sup x∈X X h(x, y)d(x, y) 2 dy is always zero, which is natural here since in this situation, kernel and Shannon entropies are equal.
For a kernel on the torus [0, 1] of the form k(x, y) = k(x − y) where k has Fourier seriesk, then we have: and for the particular choice ofk(ω) ∝ e −σ|ω| , we have: h(x, y) = tanh σ 2 tanh 2 σ 4 1 + sin 2 π(x − y) sinh 2 (σ/4) −2 , and for d(x, y) = | sin π(x − y)|, then, using exact integration, we get: For the d-dimensional torus, and, d(x, y) 2 = d j=1 sin 2 π(x j − y j ), we get the bound dσ 2 16 , and hence a nice scaling in dimension (non-exponential). Overall, in this situation, we thus get an approximation of entropy from Eq. (9) up to O(σ 2 ), where σ is the width of the kernel. Thus, if our goal is to estimate the true relative entropy, and not only have a lower bound, we need to have σ going to zero, which corresponds to RKHSs which are bigger and bigger. See discussion in Section 5.3 when estimating quantities from a finite sample.
Note that this extends to all translation invariant kernels on [0, 1] d (and not only fork(ω) ∝ e −σ ω 1 ), and more generally to translation-invariant kernels on compact spaces.
Alternative proof through leverage scores. While we provided a proof of the O(σ 2 ) approximation of the entropy through quantum data processing inequalities, the relationship between the regular entropy and the kernel entropy can also be studied through the integral representation in Eq. (5). We indeed have, using Eq. (5): We now consider that both p and q have densities with respect to the base measure τ (which we also denote p and q). The main idea, using and extending results from [38], is that for all x and λ, and p smooth enough, the quantity ϕ(x), (Σ p + λI) −1 ϕ(x) , often referred to as a "leverage score" can be approximated We thus consider the potential approximation: which is exactly the traditional KL divergence D(p q). In order to show a bound similar to Eq. (9), a non-asymptotic extension of results from [38] could be carried out.

Estimation from finite sample
Given that the covariance operator is an infinite-dimensional operator, naively we would need eigenvalue decompositions of such operators to compute the kernel information quantities, which is computationally hard. Thus finite-dimensional estimation algorithms are needed. In this section, we provide algorithms for estimating the kernel entropy H(Σ p ) from independent and identically distributed (i.i.d.) samples from p (as required for model fitting [13] or for independent component analysis [35]), as well as from certain kernel integrals. Estimators from other oracles on p or q, as well as estimators for D(Σ p Σ q ) could derived and analyzed similarly. 10

Estimators
Given Σ p = X ϕ(x)ϕ(x) * dp(x) and x 1 , . . . , x n sampled i.i.d. from p, we consider the natural estimator The natural estimate for tr Σ p log Σ p is tr Σ p logΣ p , that can be computed from the kernel matrix Proposition 6 (entropy for empirical covariance estimators) WithΣ p defined in Eq. (10), and the kernel matrix defined above, we have Proof The non-zero eigenvectors ofΣ p belong to the image space ofΣ p and are thus linear combinations , which implies K 2 = nλKα and then Kα = nλα since Kα = 0. Thus, the non-zero eigenvalues ofΣ p are exactly the ones of 1 n K, and we thus get tr Σ p logΣ p = tr 1 n K log 1 n K .
More generally, for any approximation of Σ p as n i=1 η i ϕ(x i )ϕ(x i ) * , with positive weights η i that sum to one (for example obtained from a finite grid and not sampling), we would get Running-time complexity. In order to compute the eigenvalue decomposition of K above (needed to compute the matrix logarithm), we need a running time of O(n 3 ). This can be reduced by using column sampling techniques [14], also known as "Nyström method" in this context [60]: by projecting the kernel matrix to the span of m of its columns, it leads to a running time in O(m 2 n), with some approximation that could be controlled explicitly [43].

Analysis
In order to quantify the difference between tr Σ p logΣ p and tr Σ p log Σ p , we could use regular perturbation arguments for eigenvalues of self-adjoint operators [32]. However, the function t → t log t has diverging derivatives around zero, which prevents their direct use. Instead, given the integral representation in Eq. (5), we see that it will be sufficient to estimate quantities of the form tr (Σ p + λI) −1 Σ p , which are usually referred so as "degrees of freedom" in kernel methods [21]. The natural estimator is tr (Σ p + λI) −1Σ p and its performance has already been thoroughly studied (see [45]), and an important quantity there is the maximal leverage score for the base distribution sup x∈X ϕ(x), (Σ + λI) −1 ϕ(x) , which could be replaced by the the average one (which is then equal to the degrees of freedom) if "leverage score sampling" is used (see [21] for more details). For simplicity, we do not consider leverage score sampling, noting that for symmetric sets, it does not change anything.
The next proposition shows that these estimation results for fixed λ can be integrated over λ to get an estimation bound for the entropy (see proof in Appendix C.4), without the need for extra regularization, leading to an estimation in O(1/ √ n).
Proposition 7 Assume (A1) and that p has a density with respect to the base measure which is greater from p, and the estimator defined by Eq. (11), we have: Note that as a consequence of Jensen's inequality, we always have E tr Σ p logΣ p tr Σ p log Σ p , as can be seen in simulations (see Figure 1).

Examples
We consider the d-dimensional As shown in Appendix C.6, we have the bound tr Σ(Σ + λI) Thus, the constant c in Prop. 7 above can be upper bounded by a constant times σ −d (see Appendix C.6 for details). We thus get an overall estimation rate for the kernel entropy proportional to σ −d/2 √ n.
When balancing with the O(σ 2 ) estimation of the entropy for Lipschitz-continuous strictly positive densities in Section 4.2, we can take σ ∝ n −1/(d+4) (lower bandwidth with more observations), leading to a rate of n −2/(d+4) , which is to be compared to the optimal rate equal to n −4/(d+4) for entropy estimation for densities with the same regularity [28]. Although our primary goal was not to estimate entropies, we get a rate which is the square root of the optimal rate.
Illustrative experiments. We consider X = [0, 1] and p(x) = 4 x− 1 2 , with negative differential entropy See left plot in Figure 1, where we compare two values of σ. Note that the estimation is empirically here in expectation from above, as expected because of Jensen's inequality (as opposed to the right plot).  from kernel integrals evaluated as uniform samples (method from Section 5.4). As expected we estimate the kernel entropy (limit when n tends to infinity) from above on the left plot, and from below on the right plot, with a faster convergence rate on the right plot.

Estimation from other oracles
If we have access to other oracles on p and q, we can consider other estimators. 11 For example, if we can estimate X k(x, y)k(x, z)dp(x) and X k(x, y)k(x, z)dq(x) for any y, z ∈ X, we can simply sample n points x 1 , . . . , x n from the base measure on X, and, given the projection operatorΠ on the span of and the estimator D(Σ p Σ q ) which is now a pointwise lower bound 12 on D(Σ p Σ q ) (and thus on D(p q)), which can be beneficial in variational inference, as presented in Section 7.
Noting that ϕ(x i ), Σ p ϕ(x j ) = X k(x i , x)k(x, x j )dp(x) (and similarly for Σ q ), and that K is the matrix of dot-products of ϕ(x 1 ), . . . , ϕ(x n ), the orthogonal projection may be done as where k(X, x) ∈ R n is the vector with components k(x i , x), i = 1, . . . , n, and K ∈ R n×n the kernel matrix, typically with some added regularization for numetical stability, and incomplete Cholesky decomposition so as to minimize the number of kernel evaluations (see, e.g., [7] and references therein), with an overall running time complexity of O(n 3 ). A key property is that it always underestimates the kernel relative entropy (and hence the Shannon relative entropy). The following proposition shows that it provides a better estimate (better dependence in n) than with i.i.d. samples (see proof in Appendix C.5).
Proposition 8 Assume (A1) and that p has a density with respect to the base measure which is bounded. Given uniform i.i.d. samples x 1 , . . . , x n from the base distribution, and the estimator defined by Eq. (12), we have, if for all λ > 0, sup x∈X ϕ(x), (Σ + λI) −1 ϕ(x) B log(C/λ) d , for some constants B, C: The proposition above applies to the torus situation in [0, 1] d and the kernel with exponential Fourier series, with now an improved dependence in n and in d, Other decays for the maximal degrees of freedom are considered in Appendix C.5.
Illustrative experiments. We consider as above p(x) = 4 x − 1 2 , and consider estimation from various numbers of i.i.d. observations for the projection. The required integrals are computed by numerical integration. See right plot in Figure 1, where we compare two values of σ. Note that the estimation is as expected here from below, and with a faster rate of convergence than for the estimation method from Section 5.1.

Multivariate probabilistic modelling
Information-theoretic quantities are particularly useful at characterizing statistical dependence between random variables. This is for example the basis of independent component analysis [17], and can be used for feature selection [51]. We explore how the new notions of entropy can be used in this goal.

Joint entropy, marginal entropy and mutual information
We assume that we have a product set X = X 1 × X 2 , and a feature space H = H 1 ⊗ H 2 , corresponding to the kernel k(x, y) = k 1 (x 1 , y 1 )k 2 (x 2 , y 2 ), and feature map ϕ(x) = ϕ 1 (x 1 ) ⊗ ϕ 2 (x 2 ). We also consider a base distribution on X 1 × X 2 which is the product of the base distributions. The tensor product is here defined as the span of all functions (x 1 , Given a joint probability p X 1 X 2 and its marginal probability distributions p X 1 and p X 2 , we can define the covariance operators Σ p X 1 X 2 on H 1 ⊗ H 2 , Σ p X 1 on H 1 , and Σ p X 2 on H 2 , as well as Σ 12 , Σ 1 and Σ 2 the covariance operators of the base distributions on X 1 × X 2 , X 1 and X 2 : We then define the marginal and joint entropies as: The product distribution p X 1 p X 2 (with independent components with same marginal distributions as X), which is different in general from the joint distribution p X 1 X 2 , has covariance operator Σ p X 1 p X 2 = Σ p X 1 ⊗ Σ p X 2 , and its entropy is: We have log Σ p X 1 ⊗Σ p X 2 = log Σ p X 1 ⊗I +log I ⊗Σ p X 2 . Moreover, because the base distribution is separable by assumption, we have ϕ(x), (log Σ 12 )ϕ(x) = ϕ 1 (x 1 ), (log Σ 1 )ϕ 1 (x 1 ) + ϕ 2 (x 2 ), (log Σ 2 )ϕ 2 (x 2 ) , thus, we exactly have: We can therefore define the kernel mutual information (not to be confused with the regular Shannon mutual information) as From the developments above, it is equal to mimicking the traditional equality for regular entropies. Moreover, it is equal to zero if and only if the variables are independent.
Proposition 9 (Characterization of independence) Assume (A1). Then I(X 1 , X 2 ) is a lower bound on the Shannon mutual information, and I(X 1 , X 2 ) = 0 if and only if Σ p X 1 X 2 = Σ p X 1 ⊗ Σ p X 2 . If the kernel k 2 is universal, this is equivalent to the random variables X 1 and X 1 being independent.
Proof The first statement is straightforward because of the interpretation as a relative entropy. The second statement is simply the consequence of the fact that Σ p X 1 ⊗ Σ p X 2 is the covariance operator of the product distribution p X 1 p X 2 .
We can make the following observations: • The extension to mutual information between more than two random variables is straightforward by considering tensor products between all marginals, with the same characterization property of independence.
• Our use of covariance operators is related to the "kernel generalized variance" from [6], which considers a joint covariance operator , and compute the mutual information for a Gaussian vector with this covariance operator. The kernel generalized variance can be generalized to more than two random variables, but then the independence characterization property is not satisfied. Moreover, the relationship with the usual mutual information is not as straightforward.
• The extension to conditional entropies is not straightforward. We could define Σ p X 2 |X 1 =x 1 in the natural way, but then we do not have the joint entropy being the sum of the marginal and conditional entropy (otherwise, because of axiomatic characterizations of entropy [20], we would get back exactly the Shannon entropy). There is also a potential for exponential families and conditional models akin to generalized models, which remains to be explored.

Data processing inequality
We have a form of data processing inequality, which is classical in quantum information theory (see, e.g., [40]), which we extend here.
Proposition 10 (Data processing inequality) Assume (A1). We have with equality implying (when the kernel k is universal) that the conditional distribution of X 2 given X 1 is the same for p and q, but the reverse implication is not true in general.
Proof The marginal operators Σ p X 1 and Σ q X 1 are obtained by taking "partial traces" which are quantum operations as defined in Section 3, hence the inequality. Following [46], we have inequality if and only which is unfortunately not always satisfied when conditional distributions are equal. However, in case of equality, the extra condition (c) from Prop. 2, becomes Σ p X 1 X 2 (Σ q X 1 ⊗ I) = Σ q X 1 X 2 (Σ p X 1 ⊗ I), and can be leveraged. Indeed, for all test functions f 1 , g 1 , f 2 , g 2 in H, it leads to: Considering f 1 and g 1 tending to Dirac functions atx 1 (which is possible since the kernel is univeral) we get X 2 f 2 (x 2 )g 2 (y 2 )dp X 1 X 2 (x 1 , x 2 )dq X 1 (x 1 ) = X 2 f 2 (x 2 )g 2 (y 2 )dq X 1 X 2 (x 1 , x 2 )dp 1 (x 1 ), which implies that the conditional distributions of x 2 given x 1 is the same for p and q.

Submodularity and conditional independence
We can now consider three variables X 1 , X 2 , X 3 , and use strong sub-additivity of the quantum relative entropy [46,40].
Proposition 11 (Submodularity) Assume (A1). We have: with equality (for universal kernels) implying that X 1 and X 3 are independent given X 2 , but in general the reverse implication is not true.
Note that we do not obtain a necessary and sufficient characterization of conditional independence, which can be obtained with other tools based on covariance operators [27].

Convex duality, log-partition functions and variational inference
Given that our kernel notions of relative entropies are lower bounds on the regular notions, by the traditional convex duality results between maximum entropy and log-partition functions (see, e.g., [58]), we should obtain upper-bounds on log-partition functions. In this section, we show how such bounds can be obtained and computed. We then show that such bounds can also be optimized with respect to the positive definite kernel.
This leads to a new family of variational inference methods, which can then be used in various probabilistic inference tasks [58], and in particular in Bayesian inference [42].
In this section, we will sometimes only assume that for our kernels satisfy (on top of being positive definite) k(x, x) 1 for all x ∈ X, noting that the upper-bound result on the regular Shannon relative entropy is still valid in this case.

Convex duality between operators
We have for self-adjoint operators (regardless of the link with input spaces and probabilities): with the constraint that A 0 is automatically satisfied (but not tr A = 1), with the optimal M equal to log A − log B.

Bounds on the log-partition function
We can now apply the duality relationships above to covariance operators. Given a bounded function f : X → R and a distribution q on X, our goal is to obtain upper-bounds on the log-partition function log X e f (x) dq(x) which is equal to, by convex duality: log X e f (x) dq(x) = sup p probability X f (x)dp(x) − D(p q).
In our situation, our operators are subject to belong to the subspace A, equal to the span of all ϕ(x)ϕ(x) * , for x ∈ X, which leads to slight adaptations. We do not necessarily assume that the kernels are universal, so that the hull of all ϕ(x)ϕ(x) * may not be equal to the intersection of A and positive operators with unit trace.
Isotropic kernels. If k(x, x) = 1, for all x ∈ X, then tr Σ p = 1 is equivalent to X dp(x) = 1. We define for a bounded function f : X → R, Using convex duality tools from Section 7.1, we have: Note that we can get rid of the constraint X dp(x) = 1 because it is equivalent to tr Σ p = 1.
A relaxation is to relax p not to be a non-negative measure (and use the positivity of Σ p as the only constraint), leading to By construction, we have a(f ) b(f ), and by convex duality, = inf M log tr exp(M + log Σ q ) such that ∀x ∈ X, f (x) = ϕ(x), M ϕ(x) .
Note that that if the problem above is non-feasible because ϕ is too small (so that f cannot be represented as a quadratic form in ϕ), then we get a value +∞.
If we can write (usually non-uniquely), f (x) = ϕ(x), F ϕ(x) for some self-adjoint operator F , then we have another representation as: As expected, we thus obtain upper bounds on the usual log-partition function.
Note that we could also get lower bounds using tools from Section 4.1, in particular for functions f that can be written as f (x) = X h(x, y)g(y)dτ (y).
Non isotropic kernels. If we only assume that k(x, x) 1 for all x ∈ X, then b(f ) defined in Eq. (13) and Eq. (14) is not any more an upper-bound but we can define instead with log(a) = log a if a > 1, and a − 1 otherwise.
If f can be represented through the operator F , and the constant function through the operator U , then we get We can check that if U = I, then the constraints in the definition of b (c) (f ) are equivalent to tr Σ p , and we exactly recover the expression of b(f ) in Eq. (14).

Properties of upper-bounds. The definition of upper-bounds on the log-partition function naturally leads to several questions and applications:
• Can these upper-bounds be tight? In situations like in Section 4.2 where the entropy lower-bound can be made tight (for example when the kernel bandwidth tends to zero), we should also get tight upper-bounds on the log-partition functions, that is, upper-bounds that tend to the classical logpartition function when the bandwidth goes to zero. We leave for future work a detailed study of these approximations.
• Can these upper-bounds be estimated efficiently from simple oracles on the function f ? We provide efficient algorithms in Section 7.4.
• Can these upper-bounds be minimized? We show in Section 7.5 that these upper-bounds happen to be convex in the kernel, so minimizing with respect to the kernel is possible.

Relationship with optimization
When adding a temperature parameter ε > 0, we can extend the traditional link between between optimization and log-partition functions. This corresponds essentially to considering a potential 1 ε f , or defining By convex duality, we get: When ε tends to zero, then b ε (f ) converges to Given that ϕ(x) = 1 for all x ∈ X, by writing M = cI − A for A positive, this is equal to which is exactly the optimization formulation of [44] based on "kernel sums-of-squares". Following the traditional relationship between log-partition functions and maxima, we can thus consider our bound on log-partition functions as a smoothed version on the maximum. Therefore, some of the approximation techniques developed for the optimization formulation [44] can be extended to our set-up as well. Moreover, as common in convex optimization, the entropy could be used for smoothing and accelerated optimization algorithms [37].

Computable bounds
In order to compute or approximate b(f ), which depends on the function f : X → R and the feature map ϕ : X → H, we will consider a finite-dimensional approximationφ : X →H whereH will be a finite-dimensional subspace of H. All of our approximations will always be upper-bounds on the true log-partition function.
We will always need an efficient representation ofÃ the span of allφ(x)φ(x) * , for x ∈ X, and sometimes ofM the hull of allφ(x)φ(x) * for x ∈ X. See examples in Section 7.6 and Section 7.7 for the torus and the hypercube.
In order to estimate the relative entropy D(Σ p Σ q ), we will use D(Σ p Σ q ) instead. If the kernel (x, y) → k(x, y)−k(x, y) is positive definite, then this is a lower bound, and thus the bounds on the regular quantities are preserved.
In order to compute D(Σ p Σ q ), we assume we can compute explicitly X ϕ(x)ϕ(x) * dp(x) and Ifφ is obtained from projection on a span of some ϕ(x 1 ), . . . , ϕ(x n ) like in Section 5.4, then this can be done if we can compute expectations of k(x, y)k(x, z) for any y, z ∈ X. We can also build set-specific approximate feature maps (as done in Section 7.6 and Section 7.7).
With an explicit representation of f . If we have an explicit (non-unique) representation of f as f (x) = φ(x),Fφ(x) . We then have if φ(x),φ(x) = 1 for all x ∈ X (a similar approximation can be obtained for b (c) (f ) for non-isotropic approximations): which is a finite-dimensional convex optimization problem. Since the functionM → log tr exp(M + logΣ q ) is smooth, the optimization problem above can be approximately solved by projected gradient descent (with unit step-size), with iteration: .
We then obtain the (approximated) optimalC asC = exp(M +logΣq) tr exp(M +logΣq) . It can be accelerated using extrapolation [8]. Note that for any feasibleM we obtain a lower bound onb(f ).
Without an explicit representation of f . If an explicit representation of f as f (x) = φ(x),Fφ(x) is not available, we can build an approximate representation of f such that f − φ(·),Fφ(·) ∞ is as small as possible, and then use the previous algorithm, with this extra approximation factor. On the torus, we could also use the same technique as [61] based on approximating the L ∞ -norm with Fourier transforms.
, and the set A is exactly the set of Toeplitz matrices. See Figure 2 for an illustration with several values of σ.

Kernel learning
The key property we will leverage is the concavity of D(Σ p Σ q ) with respect to the kernel. This allows to maximize efficiently with respect to the kernel to get as close to the Shannon relative entropy D(p q) as possible.
Concavity with respect to the kernel. Given two distributions p and q, then D(Σ p Σ q ) depends on the kernel function k : X × X → R. This happens to be a concave function of k (which is only assumed to be positive definite with no bounds on k(x, x)). Indeed, if k 1 and k 2 are two kernels with feature maps ϕ 1 : X → H 1 and ϕ 2 : X → H 2 , with the associated covariance operators Σ (1) p , Σ and η 1 , η 2 0 such that η 1 + η 2 = 1, then with H = H 1 ⊗ H 2 and ϕ(x) = ( √ η 1 ϕ 1 (x), √ η 2 ϕ 2 (x)), using E 1 = I 0 0 0 and E 2 = 0 0 0 I , using the monotonicity of relative entropy for quantum operations, we have: which shows the concavity.
Moreover, the relative entropy is monotonic, that is, if . Therefore, when optimizing over the kernel, we only need to look at maximal elements with respect to the positive definite order.
Linear parametrization of kernels. Given a finite-dimensional feature map ψ : X → R d , we look at kernels of the form where Λ ∈ R d×d is a symmetric positive semidefinite matrix. In order to have the upper-bound on entropy, we consider the set L of matrices Λ such that for all x ∈ X, ψ(x) ⊤ Λψ(x) 1. This is a finite-dimensional convex set, for which we need to know inner approximations (to preserve directions of bounds). We then have that is convex in Λ (see a direct proof in Appendix A.3). As mentioned earlier, by monotonicity, we only need to consider maximal elements of L.
We can decompose it as the difference of two convex functions as: This naturally leads to an optimization algorithm where we lower-bound the first convex function by its tangent and then maximize this affine function minus the other convex functions. This is a parameter-free globally and monotonically convergent algorithm.
Log-partition minimization. Given q and f : X → R, we can try to minimize the lower bound on the log-partition function instead. We assume for simplicity that we know a representation of f as f (x) = ψ(x) ⊤ Gψ(x), that the span B of all ψ(x)ψ(x) ⊤ for x ∈ X is manageable, and that (for simplicity), ψ(x) ⊤ Λψ(x) = 1 for all Λ ∈ L.
We thus need to minimize with respect to Λ ∈ L the following function (which is convex in Λ as the supremum of convex functions): The optimal C and M can be obtained by accelerated projected gradient descent. We can then minimize with respect to Λ by updating Λ like for the entropy maximization, once M (and thus C) have been estimated.
We now provide two illustrating examples, on the torus and the hypercube.
Log-partition minimization. In order to compute (and then minimize) the upper-bound on the logpartition function, we can write Note here that by the classical Caratheodory interpolation theorem (see, e.g., [24,Chapter 2]), this is exactly a(f ) which is computed.
It can be minimized by (accelerated) projected gradient descent in M , since the function M → log tr exp(M ) is 1-smooth. We simply need to be able to project on the affine subspace, which can be done by solving a linear system. 13 In order to update the kernel η, we use the exact same update as for maximizing the entropy. This is illustrated in Figure 4, where we show benefit of estimating the kernel rather than using a uniform one.
This is to be compared to the upper-bound of [31], which is 1 2 log det C(p) + 1 3 In Figure 5, we compare these bounds on matrices C corresponding to independent components with mean vectors uniformly distributed, showing the potential benefits of the new upperbound on entropy.

Conclusion
In this paper, we presented properties of the von Neumann entropy and relative entropy of covariance operators obtained from reproducing kernel Hilbert spaces. These notions are intimately related to the usual notions of Shannon entropy and relative entropy, share many of their properties, and can be applied to any type of data where positive definite kernels can be defined. They also come with additional computational benefits, with several estimation algorithms with theoretical guarantees. We have also highlighted several properties in terms of multivariate probabilistic modeling or variational inference, but several interesting avenues for future research are worth exploring. First, exploring further approximation guarantees for logpartition function estimation can probably lead to new improved algorithms for variational inference [58] in continuous and discrete domains. Second, we could extend to Rényi entropies [36,25] with similar developments (along the lines of Appendix A.4). Finally, the optimization with respect to the probability measure, which is at the core of log-partition estimation, can be extended to other optimization problems on probability measures, such as already done with similar tools for optimal transport [57] or optimal control [11].

A Proof of convexity results on operators
In this section, we provide simple direct proofs of convexity results that we use in this paper. All of these proofs are obtained from the literature, and in particular [33,47,46].

A.1 Joint convexity of the relative entropy
We consider a function f : R * + → R, which is operator-convex, that is, for which the function A → f (A) is convex as an operator-valued function, i.e., f These functions can exactly be represented as: where γ, δ > 0 and ν is a positive measure with finite mass on R + [33]. We have f (1) = α. The interesting case for this paper 2 dλ, but there are others, such as, for p ∈ (0, 1]: Given A and B two Hermitian operators, we denote by R A the right multiplication by A and L A the left multiplication by A. We have, since L B and R A commute, and log L B = L log B and log R A = R log A : where M, N = tr(M * N ) denotes the usual dot-products between operators. Following [47], we thus have: Note we can ensure the finiteness of D(A B) as soon as B −1/2 AB −1/2 is bounded. Given that the mapping (u, M ) → u, M −1 u is convex, we obtain the joint convexity of the relative entropy. It can also be shown more directly as shown below.

Proposition 12 Given positive invertible self-adjoint operators
with equality if and only if for all j such that λ j > 0, Proof Following [46, Theorem 7], we have: which is greater than zero since the function C → tr exp(K + log C) is concave. There is equality if and only for all i such that λ i > 0, A i = exp(log A − log B + log B i ), because of the Bregman divergences being equal to zero.
More generally, we get for random operators A and B, E D(A B) D(EA EB), with equality if and only if almost surely, log A − log B = log EA − log EB. Note that the concavity of the function B → tr exp(M + log B), which is crucial in deriving concentration inequalities for matrices and operators [56], is simply obtained by partial Fenchel dualization of the relative entropy.
We can also provide more refined relationships at equality, following [39].

Proposition 13
Given positive invertible self-adjoint operators A i , B i , and strictly positive real numbers λ i that sum to 1, then Proof Given the representation above, and the fact that (u, M ) → u, M −1 u is convex, we get the joint convexity of (A, B) → D(A B) is convex. In order to study the equality sign, we can use the lemma below, and for each λ > 0, . By letting λ go to zero, we get the result.
Note that following [39], we can also consider that Proof We consider the function f (u, M ) = u ⊤ M −1 u, with Taylor expansion We now use the following result on the remainder in Jensen's inequality (which itself is a direct consequence of the Taylor expansion with integral remainder taken at X, while expanded at EX): In our particular situation, by the invertibility assumption, the matrix M is invertible almost surely, and having equality in Jensen's inequality imposes to have (X − EX) ⊤ f ′′ (tX + (1 − t)EX)(X − EX) = 0 almost surely in t and X. Thus, at t = 0 and X = (u i , M i ), we get: A.2 Monotonicity of the relative entropy Proof We have, with ω = exp(2iπ/n), using the joint convexity of the relative entropy and the expression 1 n n k=1 ω k(i−j) = 1 i=j : Consider the matrices M and N defined by blocks M ij = C i AC * j and N ij = C i BC * j , as well as, the matrix U k defined by blocks (U k ) ij = 1 i=j ω ki I, we get from the equation above: because each U k is unitary. Moreover, we can write M = V AV * and N = V AV * , with a an operator V defined by a column of blocks C i , i = 1, . . . , n, so that , which shows the desired bound.
Using propositions above, there is is equality if and only: log We can also obtain necessary conditions for equality, that is, This extends to a random operator C, such that E C * C = I, for which we get: A particular example is obtained from n Hermitian positive operators D 1 , . . . , D n such that n i=1 D i = I, and we get, with µ i = tr(D i A) and ν i = tr(D i B): This corresponds to, if

A.3 Concavity of relative entropy
In this section, we consider two positive self-adjoint operators, M and N such that M αN for some α.
We consider the function: It happens to be convex in (M, N ). It is also concave in Λ in the set of positive bounded operators. This is due to the representation, valid for invertible operators:

A.4 Extension to other divergences
As mentioned in [47], both the joint convexity and the monotonicity of the relative entropy D(A B) in Prop. 12 and Prop. 14 can be extended to most operator-convex functions on R * + , and with the relative entropy being D − log . Among other interesting cases, we have: q(x) , we can then bound as We assume Lipschitz-continuity for p and q in the following form: for all x, y, q(x) q(y) − 1 d(x, y)C q , d(x, y)C p/q for some distance d on X. Then the first term in the bound above is less than 2C 2 p/q sup h(x, y)d(x, y) 2 dτ (y). For the second term, we write: q(y)q(y)

C Convergence of estimation algorithms
In this Appendix, we give convergence proofs for the estimation algorithms from Section 5, starting with concentration inequalities for operators and bounds on degrees of freedom.

C.1 Concentration of sums of random self-adjoint operators
We recall in this Appendix classical results on the concentation of covariance operators.

Thus,
We have: This leads to With u = 3/4, we get: Using df(λ) df max (λ) = d, we have the bound: .
If d n, we further get:
Proof We follow a more direct strategy as the previous proof, following estimation results for random feature expansions and column sampling [43,5].