Kernelized Cumulants: Beyond Kernel Mean Embeddings
Authors: Patric Bonnier
Harald Oberhauser
Zoltán Szabó
Venue: NeurIPS 2023 spotlight (2023)
Keywords: kernel
cumulant
mean embedding
Hilbert-Schmidt independence criterion
maximum mean discrepancy
TL;DR: We define cumulants for random variables in RKHSs and demonstrate their properties and usefulness.

Abstract: In $\mathbb{R}^d$, it is well-known that cumulants provide an alternative to moments that can achieve the same goals with numerous benefits such as lower variance estimators. In this paper we extend cumulants to reproducing kernel Hilbert spaces (RKHS) using tools from tensor algebras and show that they are computationally tractable by a kernel trick. These kernelized cumulants provide a new set of all-purpose statistics; the classical maximum mean discrepancy and Hilbert-Schmidt independence criterion arise as the degree one objects in our general construction. We argue both theoretically and empirically (on synthetic, environmental, and traffic data analysis) that going beyond degree one has several advantages and can be achieved with the same computational complexity and minimal overhead in our experiments.

Kernelized Cumulants: Beyond Kernel Mean Embeddings
Zoltán Szabó² Mattematical Institute, University of Oxford *Department of Statistics, London School of Economics bonnier, oberhauser@maths. ox. ac. uk
Patric Bonnier
Harald Oberhauser
2. szabo@lse.
Abstract
In IRd, it well-known that cumulants provide an alternative to moments that can achieve the same goals with numerous benefits such as lower variance estimators. In this paper extend cumulants to reproducing kernel Hilbert spaces (RKHS) using tools from tensor algebras and show that they are computationally tractable by a kernel trick. These kernelized cumulants provide new of all-purpose statistics; the classical maximum mean discrepancy and Hilbert-Schmidt indepen- dence criterion arise the degree one objects in general construction. We argue both theoretically and empirically (on synthetic, environmental, and traffic data analysis) that going beyond degree has several advantages and be achieved with the same computational complexity and minimal overhead in our experiments.
kernel, cumulant, mean embedding, -chbbrtt independence criterion, kernel Lancaster interaction, kernel Streitberg interaction, maximum mean discrepancy, maximum variance discrepancy
Keywords:
Introduction
The moments of random variable are arguably the most popular all-purpose statistic. However, cumulants are often more favorable statistics than moments. For example, pm == E[x"" denotes the moments of real-valued random variable X, then μl2 = μ4 + Var(X) and hence the variance that directly measures the fluctuation around the mean is much better statistic for scale than the second moment /12, see Appendix A. Cumulants provide a systematic way to only record the parts of moment sequence that are not already captured by lower-order moments. While the moment and cumulant sequences (pm) ) and (Km)m carry the same information, cumulants have several desirable properties that generalize to Rd-valued random variables (McCullagh, 2018). Among these properties of cumulants, the ones that are important for our paper is that they can characterize distributions and statistical (in)dependence.
Kernel embeddings. Mean and covariance arise naturally in the context of kernel-enrice eerrellenrichedl domains. Kernel techniques (Schölkopf and Smola, 2002; Steinwart and Christmann, 2008: Saitoh and Sawano, 2016) provide principled and powerful approach for lifting data points to so-called reproducing kernel Hilbert space (RKHS; Aronszajn 1950). Considering the mean of this feature referred as kernel mean embedding (KME; Berlinet and Thomas-Agnas 2004; Smola al. 2007) also enables one to represent probability measures, and to induce a semi-metric referred to as maximum mean discrepancy (MMD; Smola et al. 2007; Gretton et al. 2012) and an independence measure called the
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Hilbert- Schemit independence criterion (HSIC"; Gretton et al. 2005). It is known that MMD is metric when the underlying kernel is characteristic (Fukumizu al., 2008; Sriperumbudur et al.. 2010). HSIC captures independence for d = 2 components with characteristic kernels (Lyons, 2013, Theorem 3.11) 1) and for > components ( (Quadrianto et al., 2009; Sejdinovie et al., 2013a; Pfister et al., 2018) with universal ones (Szabó and Sriperumbudur, 2). MMD belongs to the family of integral probability metrics (IPM; Zolotarev 1983; Müller 1997) when in the IPM the underlying function class is chosen to be the unit ball of the RKHS. MMD and HSIC (with d = 2) are known to be equivalent (Sejdinovic et al.. 2013b) to the notions of energy distance (Baringhaus and Franz, 2004; Székely and Rizzo, 2004, 2005) also called N-distance (Zinger et al., 1992; Klebanov, 2005)- and distance covariance (Székely et al.. 2007; Székely and Rizzo, 2009; Lyons, 2013) of the statistics literature. Both MMD and HSIC can be expressed in terms of expectations of kernel values which can be leveraged to design efficient estimators for them; we will refer to this trick as the expected kernel trick. A recent survey on
embedding and their applications is given by Muandet et al. (2017). The closest previous work to ours is by Makigusa (2020) who considers the variance in the RKHS which in setting can be identified as the first kemeliied cumulant after the kernel mean embedding - for two-sample testing like we do in parts of this paper. Unfortunately, Makigusa (2020) does not provide conditions on the validity of the resulting maximum variance discrepancy, and is not formulated in the context of cumulant embeddings.
Contribution. The main contribution of our paper is to introduce cumulants of random variables in RKHSs and to show that under mild conditions the proposed kernelized cumulants character- ize distributions (Theorem 2) and independence (Theorem 3). Thanks the RKHS formulation, kernelized cumulants have computable estimators (Lemma 2 and Lemma 3) and they show strong performance in two-sample and independence testing on various benchmarks (Section 4). Although cumulants are a classic tool in multi-variate statistics, they have not received attention in the kernel literature. The primary technical challenge to cireumvent in the derivation of fundamental properties of cumulants is the rich combinatorial structure which already arises in RRd from their definition via moment- -ennerainn function which is closely linked to the partition lattice (Speed, 1983, 1984). In an RKHS, even the definition of cumulants is non-straighfforward. The key insight for our extension is that the combinatorial expressions for cumulants in Ed can be generalized by using tools from tensor algebras. This in turn allows to derive the main properties of the RKHS cumulants that underpin their statistical properties.
Broader impact & limitations. We do not any direct negative societal impact arising from the proposed new set of all-purpose kernel-based divergence and dependence measure. Choosing the underlying kernels in an optimal fashion for MMD in two-sample testing (Hagrass et al., 2022) or goodness-of-ff testing (Hagrass et al., 2023)-and showing optimal rates-even for MMD with radial kernels Rd (Tolstikhin et al., 2016) are quite challenging problems requiring dedicated analysis, and are not addressed here.
The paper is structured as follows: In Section 2 we formulate the notion of cumulants of random variables in Hilbert spaces. In Section 3 we prove a kernelized version of the classical result of valued random variables on the characterization of distributions and independence using cumulants. We show that one can leverage the expected kernel trick to derive efficient estimators for our novel statistics, with MMD and HSIC arising specifically as the "degree 1 objects. In Section we demonstrate numerically that going beyond degree is advantageous. We provide a technical background (on cumulants, tensor products and tensor algebras), proofs, further details on numerical experiments, and our V-statistic based estimators in the Appendices.
Outline.
2 Moments and cumulants
We briefly revisit classical cumulants and define cumulants of random variables in Hilbert spaces.
'HSIC is MMD with the tensor product kernel evaluated on the joint distribution and the product of the marginals, equivalently HSIC equals to the Hilbert Schmidt norm of the cross- covariance operator.

---

Moments. Let y be probability measure on R and (Xi...., Xx) m. The moments μ()) = (1"()))EeN) o are defined as
( x3] ER:
where i (i1,..., ........dd (d) E Nd denotes a d-tuple of non-negative integers (i,...., ,.. 2 0). The degree of an element E N4 is defined as (ig(1):== i .... For m E N, let μ (7) = (4 (4²())) which we refer to as the m-th moments of y with the convention that μ (7) = =1
Cumulants. Cumulants x(7) = (N*(()))EEN can defined by the moment generating function as
0¹
Σ iEN4
(1)
log Σ
(7) iEN4
O = -........ Oa) E ER
i!
i!
where denote i! = in!! and 0 ... 0
equivalent definition of cumulants is via a combinatorial expression of partitions (elaborated in Appendix C.1). Cumulants have several attractive properties, the following forms main motivation.
Theorem (Characterization of distributions with cumulants on I& from Proposition I in Jammala- madaka et al. 2006). Let 7y be a probability measure on a bounded subset of with cumulants k(T)) and let (X1,..., Xa) 7. Then
1. 7 + x(?)) injective.
2. X1,..., Xa are jointly independent if and only if vd (7) 0 for all d-tuples of positive integers
i e N
2.1 Moments in Hilbert spaces
Instead of directly considering the law of a tuple of random variables (X1, 1X1.... X2) product
space X; x ... X N,, it can be advantageous to use feature maps o; X; - H4 and instead study the distribution
the Hi x ... x Ha-valued random variable (X (Xi))..., Px(X))). Motivated by this lifting, we study here moments of Hilbert-space valued random variables and assume in this subsection (with a slight abuse of notations) that one has already applied the lifting and X; Hi where i 1,...,d d. In Section 3 we specialize the construction RKHSs, and use these moments (Def. 1) to define kernelized cumulants.
Moments. In the tions of products of the coordinates of the underlying random variable. For the fiinite-dimensionll case, it convenient to develop a coordinate-free definition which can be accomplished by using tensors. To do so we make use of the following results about Hilbert spaces: for real Hilbert spaces Hi and H2 the tensor product Hi H2 is the Hilbert space given by completion of the tensor product of Hi and H2 as vector space; we also write H, Hi Ø ... & Hi. Similarly, the direct sum
case (1) we defined the moment sequence by taking expecta-
m- m-limess
H4 @ Hg is a Hilbert space. It is natural to consider E Hi-valued random variable X1 where the integral in the expectation is meant in Bochner sense. Consequently the natural state space for all moments of Hlyvaleed random variable is the tensor algebra Ti == IIm Hhm where by convention H100 = RR. See Appendix for more details on tensor products of Hilbert spaces and tensor algebras.
[x f0m
E H as the m-th moment of
Example 2.1 (Hi = Rd, m = 2). X, = (Xt...., X() is Hi = RRd-vluued then E [X E (Rd)62 can be identified with a (d d)-sized matrix whose (i, entry is E Xixi
Since we are interested in the general case of H4 x ... x Ha-valued random variable X = (XI..... Xx) arrive at the definition below.
Definition (Moments in Hilbert spaces). Let r be probability measure on H: Hi ... x Ha and let (XI,...,XX Xa) ~ y. We define
(((r) = E @ ... X E Hi,
H0: H.g 2 ... & H
(3)
3

---

for every i E No whenever the above expectation exists. The moment sequence is defined as the element
m(7) = (T))EEE E ET T1 @ ... & Ta. with T, = II
H 8m
m20
and for m N we refer to (>) =i (z as the m-moments of 7
In case of H; = R. both definitions (1) and (3) apply for (²()). Henceforth, we always refer to (3) when we write (²(7). Even in the fnnee-dimensional case, Def. is useful, for instance when X1 Hi and X2 E H2 have different state space (H1 H2).
3 Kernelized cumulants
We lift a random variable X= (X,,..., E x= x x...X Xa via a feature map 2 : + H into Hilbert space valued random variable ((X). For the rest of the paper (i) . X' will denote a collection of Polish spaces, but the reader is invited to think of them as fnite-dimensional Euclidean spaces, (ii) 7 is an RKHS with kernel &t and canonical feature map ((x) k(z, ))² and (iii) all kernels are assumed to be bounded. Our main results (Theorem 2 and Theorem 3) are that in this the expected kernel trick applies both items in the kernelized version of Theorem 1. The key to these results is an expression for inner products of cumulants in RKHSs (Lemma I).
combinatorial expression of cumulants. Classical cumulants can be defined via the moment generating function via combinatorial sums over partitions (Appendix C.1). To generalize cumu- lants to RKHSs the combinatorial definition is the most efficient way. A partition π of m elements is family of non-empty, disjoint subsets -1.--.., TT of {1,.... m) whose union is the whole set; formally U²=1= Jj = [I,...,,,, and म; n MT = Ø for + j. We call b the number blocks of the partition ı and use the shorthand ||| to denote it. The set of all partitions of m is denoted with P(m). To formulate main results, it is convenient associate with measure y and a partition the so-called partition measure 7/x that is given by permuting the marginals of 7.
Definition 2 (Partition measure). Let y be probability measure on X1 X ... x Na and x E P(d).
Define
Yr 7=:=/<x, 60 ... 10 7xx12;
where Xx, denotes the product space IīE X, and ylx,, is the corresponding marginal distribution of We call "i the partition measure induced by
We also associate with y and multi-index i the so-called diagonal measure 7² that is given by repeating marginals according to
Definition (Diagonal measure). Let 7 a probability measure on Xi ... x Xa and (in. (i1....,id) id) E Nd. Define
-::= Law X1, (I,,..,, X1, X2....,X. X Xe....
Xa....,Xu Xa),
$1 times
ig times
La times
where (X1..... Na) ~ We call 7² the diagonal measure induced by i.
In general, the partition measure " and the diagonal measure are not probability measures on Xi x ... x X₂ but spaces that constructed by permuting repeating , Xa. Formally, Ya is probability measure Xx, .... Xr. and 7 is probability measure on X x X" thus, 7> has coordinates and 7 has dig(() coordinates. These constructions be combined, writing 7 for the measure (7²), which makes sense whenever P(deg(i)). We can now write down our generalization of cumulants.
(ax -) denotes function + k( x') with E fixed.
A kerne k: X x X R. is called bounded if there exists B R such that sup k(x, ) B.

---

Definition 4 (Kernelized cumulants). Let be probability measure on X1 X .... X Xa and let (Hi,ki), (Ha ka) be RKHSs on N....., ,., respectively. We define the kernelized cumulants
Rk Kkk.....aa = (n (a)
iENu ET
as follows
()b.....a( (7))
Σ
CxlEnng 2 k (* 1......... Xm)).),
xEP(m)
where m = deg(i), Cx == (-1) 1)1, 1/4 = (7') and
(x1.. Fm); (111...,, :m)))= ;m)) ... ka(am-isth Um-- kaa Hm)
kn (1111111 111 his (3 Mhi)
(4)
is the reproducing kernel of Hi where H H1 x Ha.
Def. 4 is the natural generalization of the combinatorial definition of cumulants in H and Appendix C.2 gives an equivalent definition via a generating function analogous (2). However, our posthoc justification that these are the "right" definitions for cumulants in an RKHS are Theorems and that show that these kernelized cumulants have the same powerful properties classic cumulants in Rd (Theorem 1).
Example 3.1 (Kernelized cumulants). Let be a probability measure on X x X2, with the RKHSs
(H1, ki)s (H2,k2) given. Denote the random variables K₁ (x K2 (xx where
(X1,Xs X2) ~ Then the degree two kernelized cumulants are given as K Kk (7) = E [K?
(2,0)
E [K]] 222 (1,1)
(0,2) Kk, (7) = E [K @ K2] E [K,]] @ [K2] Kk,, (n) = k E K2 82
Inner products of cumulants.
Computing inner products of moments is straightforward thanks to a nonlinear kernel trick, see Lemma the Appendix. For example, given two probability measures 71, with corresponding random variables (Xi,..., .a) 71, (Yi...., Ya) 722 on X x ....x.x X
and RKHSs (H1, ki),.... (Hay ka) on X1..... Ns with bounded kernels, and = H4 x x Ha.
can express:
(phi k (12)) on (X,, ..........aaaaaaaa (a Yx),
(5)
where "e,,.. _ka is defined in Def. 1, and the expectation is taken over the product measure M '2- Example 3.2. the purticular case of d (5) reduces to the well well-known formula for the inner product of mean embeddings (p) (T), (12))
k(X, (X,,)) Y).
Lemma (Inner product of cumulants). Let (H1,, kı),....... (Ha, ka) be RKHSs with bounded kernels on X1,- .......4 .....,,,, respectively, and let y and n probability measures Xi x Xz i (ii....,ia) ia) E Nd such that deg(i) = m. Then
(n (7), M (1))
Σ 7.7EP(m
(Yi,..., C+ Byg Jec.Bygon en 1. (Y,, Ym)).
Point separating kernels. In the classic MMD setting the injectivity of the mean embedding + [x>,[k((,, :)] on probability measures (known as the characteristic property of the kernel k) is equivalent to the MMD being metric; this property central in applications. We formulate our theoretical results in the next section using the much weaker property of what we term "point- separating" which satisfied for essentially all popular kernels.
Definition 5 soint--eparating kernel). We call kernel X X R point-separating if the canonical feature map & x + k(I, -) is injective.
3.1 (Semi-)metrics for probability measures
In this section we use cumulants to characterize probability measures and show how to compute the distance between kernelized cumulants with the expected kernel trick.
5

---

Theorem 2 (Characterization of distributions with cumulants). Let ry and d be two probability measures on X that for every 1 I j s d kj is bounded, continuous, point-separating kernel. Then
Xa. (Hı,kı)..... .H.....H.. (H.aka) ka) RKHSs on the Polish spaces Xi...., Nx such
y 1y if and only if (kk.....ka
KKk.....k (1)))
Moreover, the expected kernel trick applies and for E N with deg(i) m, and k and 71 in (4)
d((r,n) == ||rc | (7)
N Kk (n)
(6)
Σ cxc
E 1 81 ((11, :m); (Yi,..., Ym))
M.TEP(m)
y enn! en! k. k ((1,, Xm), Mm))((n, (V,,.. Y.n)
2 so! 1 ((X1,..., Xm),(Y (.,,.,, Ym))
We recall Example 3.1 and now give examples of distances between such expressions Example 3.3 (m = 1). Applied with m = and d = 1, (6) becomes MMD (7, 7) Ki (0) (n) = EK(X, X') + EE((Y, Y") 2Ek(X, Y),
| [x) (7)
where X, X' denotes independent copies of and Y, Y' denotes independent copies of 1I.
Example 3.4 (m 2). For 2 and d = 1, (6) reduces to
1 (7) Kik (n)l = (EX(XXX')CCC) X")"(X, Xx") + E&(Y, Y')A(Y Y"") BA(X, x1) +EEK(Y + 2X ()x(X'XXXX)))))) 2EX(X, Y) Y") 2EK(X, Y") -X, Y)² 2XA(X, (XX, 2ER(Y,
where X, X'X', x' denotes independent copies of " and Y, Y.Y",Y", Y",Y" denotes independens copies of rj. This expression compares the variances in the RKHS instead of the means. This is an example of the kernel variance embedding defined in the next subsection.
The price for the weak assumption of a point- separatin kernel is that without any stronger assump- tions one does not get a metric in general, and the pul-pee way to achieve a metric is to take an infinite
over all d''s. If we only the degree MMD formula which requires characteristicness to become a metric (see Example 3.3). There are two reasons why working under weaker assumptions is useful: firstly, if the underlying kernel is not characteristic this sum gives a structured way to incorporate finer information that discriminates the two distributions; an extreme case is the linear kernel k(x, y) = (a, v) which is point-separating, and in this case the sum reduces to the differences of classical cumulanns. Secondly, under the stronger assumption of characteristicness one already has a metric after truncation degree m (the classical MMD). However, in the finite-sample case adding higher degree terms can lead to increased power. Indeed, our experiments (Section 4) show that even just going one degree further (i.e. taking m 2), can lead to more powerful tests.
term d' reduces the well-known
3.2 A characterization of independence
Here characterize independence in terms of kernelized cumulants.
Theorem 3 (Characterization of independence with cumulants). Let be a probability measure on X1 X X Xz,, and (H,,k1), kı),.... (Ha, ka) RKHSs on Polish spaces Ni...., Na such that for every 1 15 j <d is a bounded, continuous, point-separating kernel. Then
= 7xx, x 8.....7..
87xx if and only if k K (7) = 0
for every E NN. Moreover, the expected kernel trick applies in the sense that for i E Nd.
Σ
Cec, Eog
C.
Xm), (Y Ym)),
T.TEP(m)
where = dig(i), and &l and 01 are defined as in (4).

---

Applied to i = (1, the expression (7) reduces to the classical HSIC for two components, see Example 3.5 below. But for general this construction leads to genuine new statistics in RKHSs.
Example 3.5 (Specific case: HSIC, kernel Lancaster interaction, kernel Streitberg interaction). = 2 there is only one order 2 index in N4, namely = (1, 1); in this case (7) reduces to the classical HSIC equation
(1.1) 1k (1ll
(X, + (X,
- 2EA1 (X,2))
where (X and (X',Y, Y') are independent copies of the same random variable following "/- More generally, with i 1a gets the kernel Streitberg interaction (Streitberg, 1990; Sejdinovic et al., 2013a; Liu et al, 2023), and specifically the kernel Lancaster interaction (Sejdinovic et al., 2013a) for d e {2, 3}; the latter reduces to HSIC for two random variables (d 2).
3.3 Finite-sample statistics
To apply Theorem 2 and Theorem in practice, one needs to estimate expressions such as EA: ((X1, ((X1,... x,), (Yi,-.. Y..). One could use classical estimators such U-statistic (Van der Waart, 2000) which lead to unbiased estimators. However, we follow Gretton et al. (2008) and use V-statistic which is biased but conceptually simpler, easier, and efficient compute. We note that the estimators presented here all have quadratic complexity like MMD and HSIC, see Appendix E.
two-sample test for coracatractriistic feature maps. i characteristic then MMD 17)
O exactly when 7 = but we can still increase testing power by considering the distance between the kernel variance and skewness embeddings, which leads us to use our semi-metries dd (7, 1) and (( 17) as defined in (6). An efficient estimator for ((3) is given in detail in Appendix E; we provide the full expression for d(2) here.
Lemma (d( estimation, see (6)). The V-statistic for d
(7, (7, = -11e (7) Ki (n)
2
N2 Tr (K,N)
Tr (K NM
Tr KayJ M M [yN]:
M
where Tr denotes trace, (xn) i.d.
lid. n, K, [k(I; F. )1. RNN, K, = [K(,,, w;)l E IR MxM, Kz.y = [K(I,, E ERNM Jn = In 11,11 E R"*",
7. (Um)
with 1, = (1
A kernel independence test. By Theorem 3. if 7 7xx, @ 7/x, then (2)))))) 0 and x....) (7) 0. We may compute the magnitude of either y.(2.1 (7) or x (1),2) we will refer to these quantities as cross skewness independence criterion (CSIC). Note that these criteria are asymmetric. When d 2 we have probability measure y on X; x X2 and two kernels k X² + R, L: X3 R. Assume that we have samples (I,, 34 and use the shorthand notation K = K., = Ly (similarly to Lemma 2) and = Hr = $ 1yly y E N Denote by o the Hadamard product and (-) the sum over all elements of a matrix. Then one can derive the following CSIC (Nott that matrix multiplication takes precedence over the Hadamard product.)
Lemma (CSIC estimation). The V-statistic for (1.22 (1)
1
KoKoL
4K KH o -K o K LH + 4KH K LH
N2
K
+ 2K o L
+ 2KH O HK о L + 4K o HK o LH K K
N
No
K
L
K
- 8K o LH
4K o HK
+
L
N2
No
N
Remark (computational complexity w.nt. degree m). We saw that the computational complexity of the cumulant based measures is quadratic w.r.t. the sample size. Let Bm. [P(m)] be the m--h
7

---

Bell number, in other words the number of elements in P(m). The Bell numbers follow recursion:
Bm+1 = Pm + 1) = Σm=
k BR, with the first elements of the sequence being Bg = B₁ = 1, (m)
B3 B3 = 5. = Bg 52. By (6)-(7), in the worst case the number of operations to compute d' -1) or || (1)) (m = deg(i)) is proportional to B² (it equals to 3 and to B² respectively). Though asymptotically Bn grows quickly (de Bruijn, 1981; Lovász, 1993), for reasonably small degrees the computation is still manageable. In addition, merging various terms in the estimator often be carried out, which leads to computational saving. For instance, the estimator of dd(2) (see Lemma 2, Example E.I,, CSIC (Lemma 3, Example E.2) and d (Example E.3) consists of only 3. 11 and 10 + 2 x 24 terms compared to the predicted worst-case setting of 3 12, B 25, and 3B3 75 terms, respectively. On practical side, we found that using m E {2,3) is a good compromise between gain in sample efficiency and ease of implementation.
4 Experiments
In this section, we demonstrate the efficiency of the proposed kernel cumulants two-sample and independence testing.*
Two-sample test: Given N N samples from two probability measures y and v on a space X, the goal was to test the null hypothesis Ho 7 v) against the alternative Hi : 7 * vy. The compared test statistics (S) were MMD, dd²), and d(3).
Independence test: Given N paired samples from probability measure on product space Xi X2. the aim was to the test the null hypothesis Ho : 7 = 71 & 7/ against the alternative Hi 1 + 71 @ 72. The compared test statistics (S) were HSIC and CSIC.
In our experiments H1 held, and the estimated power of the tests reported. Permutation test was applied to approximate the null distribution and its 0.95-quantile (which corresponds to the level choice a 0.05): We first computed our test statistic S using the given samples (Su S), and then permuted the samples 100 times. If So was in high percentile (2 95% in our case) of the resulting distribution of S under the permutations, we rejected the null. We repeated these experiments 100 times to estimate the power of the test. This procedure was in turn repeated 5 times and the 5 samples are plotted as a box plot along with a line plot showing the mean against the number of samples
||x-y11
(N) used. All experiments were performed using the rbf- ref-kernel rbf (x, y) parameter a is called the bandwidth. We performed all experiments for every bandwidth of the form a10" where bandwidths was chosen for each method and sample size. The experiments were carried out on a laptop with an i7 CPU and 16GBs of RAM.
20
where the
1, 2.5, 5, 7.5 and b = -5, -4 -3, 2,, -1, 0 and the optimal value across the
4.1 Synthetic data
For synthetic data we designed two experiments.
2-sample test: We compared a uniform distribution with mixture of two uniforms.
Independence test: We considered the joint measure of uniform and a correlated x² random variable. We also use this same benchmark to compare the efficiency of classical and kernelized cumulants in Appendix D.
Comparing a uniform with a mixture of uniforms. Even for simpler distributions like mixtures of uniform distributions it can be hard pick up higher-order features, and ((2) can outperform MMD even when provided with moderate number of samples. Here we compared one uniform distribution U1- 1] with an equal mixture of (0.335, 0.778) and U|30.35, 0.778). The endpoints in the mixture were chosen to match the first three moments of 1]. The number of samples
used ranged from 5 to 50, and the results are summarized in Fig. 1. One can see that with d(2) the power approaches 100% much faster than with using MMD.
'All the code replicating our experiments is available at https://githu// com http://hhhw...ce..aaaacooooeeeee Kernelized- Cumulants.
8

---

100
100
Test power (%6)
Test power (%)
80
20
20
28
N
N
Figure 2: Test power function of the sample size (N) of independence testing using HSIC (red) and CSIC statistics (blue), with the independence testing between Yo.5 and X.
Figure 3: Two-sample testing using MMD (red)
and d(2) (blue) on the Scoul bicycle data set.
Independence between uniform and x². Let X ~ 00 1] and Z ~ NO, 1) be indepen- dent of X. Denote by & the c.d.f. of standard normal distribution and define Yp to be a mixture with weight p at (X) and weight 1 - 1-par pa Z so that Y, Z and Y1 + (X). We test for independence for = 0.5 between Y Y2-which will be x² distributed with degree of freedom- and X. As the statistical dependence of Y and X is more complicated than a simple correlation we expect that higher-order features of the data will help in the independence testing. The num- 60, with the results summarized in Fig. 2. One can see that CSIC supersedes HSIC for every sample size, and the difference is more pronounced for smaller ones.
80
Test power
20
N
ber of samples used ranged from
Figure Test power as function of the sam- ple size (N) of two-sample test using the MMD (red) and d(2) statistics (blue), with U[-1,1] and an equal mixture of (00.35, 0.778] and U Ul-0.35, -00 00.778) compared.
4.2 Real-world data
We demonstrate the efficiency of the kernelized cumulants on real-world data. We designed two experiments.
2-sample test: Here the goal was test if environmental data two different seasons, and traffic data at different speeds can be distinguished.
test if two distributions describing traffic flow and other traffic independence test: The aim factors are independent.
To improve performance of both test statistics, all features are standardized to lie between 0 and 1.
Scoul bicycle data. The Seoul bicycle data set (E et al., 2020) consists of environmental data along with the number of bicycle rentals. The environmental data consists of 9 numerical values, categorical value (season), and two binary values. We compare the distribution of the environmental data in the winter and the fall, 2-sample testing on two measures 7, 77 on 11 and assume that 7 + 7) where is the distribution of the environmental data in winter, and is that of the data in the fall. Permutation testing was performed for N between and 44, with results summarized in Fig. 3. As it can be observed that d(2) outperforms MMD terms of test power, For the Type error, i.e. the probability of falsely rejecting the null hypothesis (comparing winter data with itself), it hovers between 5 10% for both statistics,
expect these distributions to be different. Concretely, do

---

100
100
Test power (%6
Test power (%)
80
20
20
50
N
N
Figure 4: Independence testing using HSIC (red) Figure 5: Two-sample testing using MMD (red)
and CSIC (blue) on the Sao Paulo traffic dataset. and d(3) (blue) on the Sao Paulo traffic dataset.
which is admittedly slightly higher than the desired 5% due to the small sample size, but very similar for both statistics; for further details, the reader is referred to Fig. 8 in Appendix D.
Brazilian traffic data. We used the Sao Paulo traffic benchmark (Ferreira, 2016) to perform independence testing. The dataset consists of 16 different integer-valued statistics about the hourly traffic in Sao Paulo such blockages, fires and other reasons that might hold up traffic. This is combined with a number that describes the slowness of traffic at the given hour; so X; = R X₂ = R. One expects strong dependence between the two sets or equivalently, for the null hypothesis to be false and for the statistics are heavily skewed towards as it naturally sparse. For independence testing performed permutation testing for N between 4 and 40. The resulting test powers are summarized in Fig. 4. As can be seen, HSIC and CSIC performs similarly for very low sample sizes, but for anything else CSIC is the favorable statistic in terms of test power. For two-sample testing, we sampled N between 5 and 50 and compared the distribution of slow moving traffic with the fast moving traffic. The results are summarized in Fig. 5. It clear that d(3) performs similarly MMD in terms of test power for very small sample sizes, but significantly better for larger ones.
Conclusion
We defined cumulants for random variables in RKHSs by extending the algebraic characterization of cumulants on Id. This construction results in a structured description of the law of random variables that goes beyond the classic kernelized mean and covariance. A kernel trick allows to compute the resulting kernelized cumulants. We applied our theoretical results to two-sample and sufficient for this task, the higher-
independence testing; although kernelized mean and covariance order kernelized cumulants have the potential to increase the test power and to relax the assumptions the kernel. Our experiments real and synthetic data show that kernelized cumulants lead to significant improvement of the test power. A disadvantage of these higher-order statistics is that their theoretical analysis requires more mathematical machinery although we emphasize that the resulting estimators are simple V-statistics.
indeed
References
Arcones, M. A. and Giné, E. (1992). On the bootstrap U and V statistics. Annals of Statisties, 20:655-674
Aronszajn, N. (1950). Theory of reproducing kernels. Transactions of the American Mathematical Society, 68:337-404.
Baringhaus, L. and Franz, C. (2004). On multivariate two-sample test. Journal of Multivariate Analysis, 88:190-206.
Berlinet, A. and Thomas-Agnan, C. (2004). Reproducing Kernel Hilbert Spaces in Probability and Statistics. Kluwer.
10

---

Bogachev, V. I. (2007). Measure theory. Vol. I, II. Springer.
Bonnier, P. and Oberhauser, H. (2020). Signature cumulants, ordered partitions, and independence of stochastic processes. Bernoulli, 26(4):2727-2757
Chevyrev, and Oberhauser, (2022). Signature moments to characterize laws of stochastic processes. Journal of Machine Learning Research, 23(17)):1-42.
Chung, E. and Romano, J. P. (2013). Exact and asymptotically robust permutation tests. Annals of Statistics, 4142):484-507
de Bruijn, G. (1981). Asymptotic Methods in Analysis. Dover.
Dudley, R. (2004). Real Analysis and Probability. Cambridge University Press.
E, V., Park, J., and Cho, Y. (2020). Using data mining techniques for bike sharing demand prediction in metropolitan city. Computer Communicutions, 153:353-366.
Ferreira, R. P. (2016). Combination of artificial intelligence techniques for prediction the behavior of urban vehicular traffic in the city São Paulo. In Anais do 10. Congresso Brasileiro de Inteligência Computacional, pages 1-5.
Fukumizu, K., Gretton, A., Sun, X., and Schölkopf, B. (2008). Kernel measures of conditional dependence. In Advances in Neural Information Processing Systems (NIPS), pages 498-496.
Gärtner, T. (2003). A survey of kernels for structured data. ACM SIGKDD Explorations Newsletter, 5(1):49-58.
Gretton, A., Borgwardt, K., Rasch, M., Schölkopf, B., and Smola, A. (2012). A kernel two-sample test. Journal of Machine Learning Research, 13(25):723-773.
Gretton, A., Bousquet, O., Smola, A., and Schölkopf, B. (2005). Measuring statistical dependence with Hilbert-Schmidt norms. In Algorithmic Learning Theory (ALT), pages 63-78.
Gretton, A., Fukumizu, K., Teo, C. H., Song, L.,, Schölkopf, B., and Smola, A. (2008). A kernel statistical test of independence. In Advances in Neural Information Processing Systems (NIPS), pages 585-592.
Hagrass, O., Sriperumbudur, B. K., and Li, B. (2022). Spectral regularized kernel two-sample tests. Technical report. fattps://arre org/ abaaa/as///// /222.. 09200)..
Hagrass, O., Sriperumbudur, B. K., and Li, B. (2023). Spectral regularized kernel goodness-of-fit tests. Technical report. (https:/ //apxii///rxvv...///// org/ abs/2300 04561).
Jammalamadaka, S. R., Rao, T. and Terdik, G. (2006). Higher order cumulants of random vectors and applications to statistical inference and time series. Indian Journal of Statistics, 68(2):326-356.
Kiraly, F. J. and Oberhauser, H. (2019). Kernels for sequentially ordered data. Journal of Machine Learning Research, 20.
Klebanov, L. (2005). N-Distances and Their Applications. Charles University, Prague.
Kriege, N. M., Johansson, F. D., and Morris, C. (2020). A survey on graph kernels. Applied Network Science, 5(1).
Lang, S. (2002). Algebra. Springer.
Liu, Z., Peach, R. L., Mediano, P. A., and Barahona, M. (2023). Interaction measures, partition lattices and kernel tests for high-order interactions. Technical report. (https://arxiv, org/ ftte/2306.00900
Lodhi, H., Saunders, C.,, Shawe-Taylor J., Cristianini, N., and Watkins, C. (2002). Text classification using string kernels. Journal of machine learning research, 2:419-444.
Lovász, (1993). Combinatorial Problems and Exercise. 2nd ed. Amsterdam, Netherlands: Forth-Holland

---

Lyons, R. (2013). Distance covariance in metric spaces. The Annals of Probability,
Makigusa, N. (2020). Two-sample test based on maximum variance discrepancy. Technical report. (attps://arxix) httpaaawww.....rraaaaaaaaaaaaae wwwp/o///g/......///// 00980).
McCullagh, (2018). Tensor Methods in Statistics. Courier Dover Publications.
Muandet, K., Fukumizu, K., Sriperumbudur, B., and Schölkopf, B. (2017). Kernel mean embedding of distributions: A review and beyond. Foundations and Trends in Machine Learning. 10(1-2):1-141
Müller, A. (1997). Integral probability metrics and their generating classes of functions. Advances in Applied Probability, 29:429-444 443.
Pfister, N., Bühlmann, P., Schölkopf, B., and Peters, J. (2018). Kernel-based tests for joint indepen- dence. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(1):5-31.
Quadrianto, N., Song. L., and Smola, A. (2009). Kernelized sorting. In Advances in Neural Information Processing Systems (NIPS), pages 1289-1296.
Rudin, W. (1953). Principles of mathematical analysis. McGraw-Hil Book Company, Inc., New York-Toronto-London
Saitoh, S. and Sawano, Y. (2016). Theory of Reproducing Kernels and Applications. Springer Singapore.
Schölkopf, B. and Smola, A. (2002). Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT Press.
Sejdinovic, D., Gretton, A., and Bergsma, W. (2013a). A kernel test for varialla interactions. In Advances in Neural Information Processing Systems (NIPS), pages 1124-1132.
Sejdinovic, D., Sriperumbudur, B., Gretton, A., and Fukumizu, K. (2013b). Equivalence of distance- based and RKHS- based statistics in hypothesis testing. Annals of Statistics, 41:2263-2299
Sertling, R. (1980), Approximation Theorems of Mathematical Statistics. Wiley Series in Probability and Statisties.
Smola, A., Gretton, A., Song, L., and Schölkopf. B. (2007). A Hilbert space embedding for distributions. In Algorithmic Learning Theory (ALT), pages 13-31.
Speed, T. P. (1983). Cumulants and partition lattices. Australian Journal of Statistics, 25(2):378 388.8 Speed, T. P. (1984). Cumulants and partition lattices II. Australian Journal of Statistics, 4(1):34-53.
Sriperumbudur, B., Gretton, A., Fukumizu, K., Schölkopf, B., and Lanckriet, G. (2010). Hilbert space embeddings and metrics on probability measures. Journal of Machine Learning Research, 11:1517-1561.
Steinwart, I. and Christmann, A. (2008). Support Vector Machines. Springer.
Streitberg, B. (1990). Lancaster interactions revisited. Annals of Statistics, (8(4):1878-1885
Szab6, 7. and Sriperumbudur, B. K. (2018). Characteristic and universal tensor product kernels. Journal of Machine Learning Research, 18(233):1-29.
Székely, G. and Rizzo, M. (2004). Testing for equal distributions in high dimension. InterStat, 5:1249-1272.
Székely, G. and Rizzo, M. (2005). A new test for multivariate normality. Journal of Multivariate Analysis, 93:58-80
Székely, G. J. and Rizzo, M. L. (2009). Brownian distance covariance. The Annals of Applied Statistics, 3:1236-1288
Székely, G. J., Rizzo, M. L., and Bakirov, K. (2007). Measuring and testing dependence by correlation of distances. The Annals of Statistics, 35:2769-2794,
12

---

Tolstikhin. I., Sriperumbudur, B., and Schölkopf, B. (2016). Minimax estimation of maximal mean discrepancy with radial kernels. Advances in Neural Information Processing Systems (NIPS), pages 1930-1938.
Van der Waart, A. A. (2000), Asymptotic Statistics. Cambridge Series in Statistical and Probabilistic Mathematics.
Willard, S. (1970). General Topology. Addison- Addison-Wesley.
Zinger, A., Kakosyan, A., and Klebanov, L.. (1992). A characterization distributions by mean values of statistics and certain probabilistic metrics. Journal of Soviet Mathematics.
Zolotarev, V. (1983). Probability metries. Theory of Probability and its Applications. 28:278-300 -322.
13

---

These appendices provide additional background and elaborate on some of the finer points in the main text. In Appendix A we illustrate that cumulants have typically lower variance estimators compared to moments. Technical background on tensor products and tensor sums of Hilbert spaces, and on tensor algebras is provided in Appendix B. We present our proofs in Appendix C. In Appendix D additional details on our numerical experiments are provided. Our V-statistic based estimators are detailed in Appendix
A Moments and cumulants
Already for real-valued random variable X, moments have well-known drawbacks that make cumu- lants often preferable as statistics. For detailed introduction to the use of cumulants in statistics we refer to McCullagh (2018). Here we just mention that
1. the moment generating function f(() = E = Mm /mm describes the law of X with sequence (lm)) of moments /mm = [ E R. However, since the function + f()) is the expectation of an exponential, one would often expect that f is also "exponential in t". hence 1m should be simpler to describe as a power series. For example, for g(t) = log f(t) = 2 Ki my!
Gaussian ft) e ( Var(X) and while jmm can be in this case explicitly calculated and uneven moments vanish, the m-moments are fairly complicated compared to the power series expansion of g(t) = kyt + 2 which just consists of Kj (mean) and K2 (variance).
2.
In the moment sequence 14m, lower moments can dominate higher moments. Hence, natural idea to compensate for these "different scales" " is to systematically subtract lower moments from higher moments. As mentioned in the introduction, this is in particular troublesome if finite samples available. Even in dimension d = 1 the second moment is dominated by the squared mean, that is for real-valued random variable X 7
12(1) (u' (T)) + Var(X),
where Var(X) [ -x (1)) It is well known that the minimum variance unbiased estimators for the variance more efficient than that for the second moment: denoting them by e and R respectively, one can show (Bonnier and Oberhauser, 2020) that given N samples from X, the following holds
(EXX) (EX)) Var()))aa) -2 2 Var(X)= N
2
Var (22 Var
+
N -
This means that when X has large mean, it is more efficient to estimate its variance than its second moment since the last term the above expression dominates. Hence, the variance Var(X) is typically much more sensible second-order statistic than 2 (7) However, emphasize that there are many other reasons why cumulants can have better properties as estimators
3. Cumulants characterize laws and the independence of two random variables manifests itself simply as vanishing of cross-cumulants In view of the above item 2, this means for example that testing independence be preferable in terms of vanishing cumulants rather than testing if moments and similarly for testing if distributions are the same.
factor E EX" = E x [
The caveat to the above points is that it is not true that cumulants are always preferable. For example, there are distributions for which (a) the moment generating function is not naturally exponential in I, (b) lower moments do not dominate higher moments, (c) consequently independence or two-sample testing become worse with cumulants. While one can write down conditions under which for example, the variance of the kernelized cumulants is lower, the use of cumulants among statisticians is to simply regard cumulants as arising from natural motivations which leads to another estimator in their toolbox.
The main idea of our paper simply that for the same reasons that cumulants can turn out to be powerful for real or vector-valued random variables, cumulants of RKHS- random variables are a natural choice of statistics. The situation is more complicated since requires formalizing moment- and cumulant-generaate the computational bottleneck of working in infinite dimensions and leads to computable estimators for independence and two-sample testing.
functions in RKHS but ultimately a kernel trick allows for circumventing
14

---

Further, we note that although cumulants are classic for vector-valued data, there seems to be not much work done about extending their properties to general structured data. Our kernelized cumulants apply to any set X' where a kernel is given. This includes many practically relevant examples such as strings (Lodhi al., 2002), graphs (Kriege et al., 2020), or general sequentially ordered data (Király and Oberhauser, 2019; Chevyrev and Oberhauser, 2022); survey of kernels for structured data is provided by Gärtner (2003).
B Technical background
In Section B.1 the tensor products (8) j=1 Hyy) and direct sums of Hilbert spaces (Dier H,) are
f=
recalled. Section B.2 is about tensor algebras over Hilbert spaces (I Hmm)
B.1 Tensor products and direct sums of Banach and Hilbert spaces
Tensor products of Hilbert spaces. For Hilbert spaces H...., Ha and (h1...., ha) E Hi x x.. X H.. the multi-linear operator hn ..... .....
ha Hi @ .... H. is defined as
d II 1=1
(h₂ 8 ... 2 (a)
Sa)
JjH
for all fai E Hi x ... x Ha. By extending the inner product
(a1 x
ad, b 2 ... @ bu) =
II
finite linear combinations of a1 @ ... & a
c 8 10 ang a E R.R. E Hl, n 2
by linearity, and taking the topological completion one arrives at Hi @ ... H.. Specifically, if
(V),,,),... (Ha; ka) are RKHSs, then so is Hi X ... 6 Ha = Hzg_rt (Berlinet and Thomas- Agnan, 2004, Theorem 13) with the tensor product kernel
(0 k;)) ((x1, (()), (2/, x()..,,)))))))) II hy (z),2)),
where (II,..., ra), ( r'd) E X1 x Xx--
Tensor products of Banach spaces. For Banach spaces Bi,... BB, the construction of 69 a little more involved (Lang, 2002) one cannot rely inner product.
Direct sums of Hilbert and Banach spaces. Let (H,)er be Hilbert or Banach spaces where is some index set. The direct sum of H written as O Ho-consists of ordered tuples / and h; for all but finite number of E
= (hę)ez such that h; E H, for all Operations (addition, scalar multiplication) are performed coordinate-wise, and the inner product of a, E Dter H, is defined as (a, b) e ke M. -ier a;b,
B.2 Tensor algebras
The tensor algebra T, Hilbert space H, is defined the topological completion of the space
m20
15

---

Note that it can equivalently be defined as the subset of (ho, h1 h2....) E Im H2n such that
Em20 mm 1 Hem oo, and such it is Hilbert space with norm
[lh,. h1, h2 ... II m2o Ammo 23
Ellmm p10mm m20
Т is also algebra, endowed with the tensor product over Hj its product. For (oo, (1, 42 02 ......)) (bo. b,, b2, b2 b2...) E Ty, their product can be written down in coordinates as
m
a b=
Σ
a4 @ bm-t
m20
For sequence H....., Ha of Hilbert spaces, we define
T T ..... Td:
where T, = 1m20 HTm I1 d) Let Hi x x Ha. and recall that given tuple of
integers i (i1,....) id) No we define Hi H a multi-grading for T as
@ H&ge This allows to write down
T
II H iENd
(8)
Note that this gives credence to us using multi-indices i N to describe elements of the tensor algebra, as the multi-indices form its multi-grading.
Furthermore, T is multi qualti algebra when endowed with the (linear extension of the) following multiplication defined on the components of T
01 708 +19
(9)
(x) 8 ... 8....822 xx)) * (y1 8 ... 8 ya) = (x1 y1) @ @ (aa ya),
so that for a = (at)
LEN4 b = (bt) IENd E T, their product can be written down as
(a*
Σ a 812
(10)
i'P==
where addition of tuples 12 E N4 is defined as + = of tuple defined as deg(i) + id, T also graded algebra, with the grading written down as
+2) With the degree
II
H
(iEN"sdog(i)=m)
so that if one multiplies two elements together, the degree of their product is the sum of their degree. Finally we note that T is unital algebra and the unit has the explicit form
(1, O.O.O....)).
i.c. the element consisting of only 1 at degree 0.
C
Proofs
This section is dedicated to proofs. The equivalence between the combinatorial expressions of cumulants and the definition via a moment generating function proved in Section C.2. The derivation of our main results (Theorem and Theorem 3) are detailed in Section C.3.
16

---

C.1 Equivalent definitions of cumulants in R
Here we introduce a classical definition of cumulants via a moment generating function and its equivalence to the combinatorial expressions. If X = (X,,..., Xa) is an R valued random variable distributed according X 7, then recall that
1 (n) = E[X ER
for (i), ( id) E No.
The following definition of the cumulants x² (7) of y are equivalent
1. EIeNs (7) 6¹ log EIENe (1)) o² i
2. x²(1) E (+H) rEP(d)
where 0 = (0),..., 0a) --D) The equivalence between these two definitions of cumulants via a generating function and via their combinatorial definition, classical (McCullagh, 2018) even if our notation here is non-standard in the classical case. This equivalence is also at the heart of many proofs about properties cumulants since some properties are casier to prove via one or the other definition.
C.2 Equivalent definitions of cumulants in RKHS
In the main text, we defined cumulants RKHS by mimicking the combinatorial definition of cumulants in Rd It is natural and useful to also have the analogous definition via "generating function" for RKHS-valued random variables. However, to generalize the definition via the logarithm of the moment generating function to random variables in RKHS, requires to define logarithm for tensor series of moments. In this part, we show that this can be done and that indeed the two definitions are equivalent.
We use the shorthand k(7) == Kkj.... (r)) (7), u((r)) == (k....... (7), and we overload the notation (X1...., X2) with (1) ka(:, Xa)). With this notation, show that given coordinates i E N one may express the generalized cumulant s²(7) as either combinatorial sum over moments indexed by partitions, or by using the cumulant generating function.
More specifically, we show that the generalized cumulant of a probability measure on H4 X x.. x Ha defined as
xt
Σ F(P(m)
Be (X
where c = (-]) - 1)! can also be expressed as coordinates in the tensorized logarithm of the moment series. Motivated by the Taylor series expansion of the classic logarithm, we define
(-1)--1
log: : T
I
(x 1),².
n
n21
where denotes the product as defined in (9) and for fort= E T, tm is defined as
coordinate- wise (tmm31 Σ +:.
t for E No. Note that unlike the classical
logarithm log R, + R, the tensorized logarithm is defined on the whole space as a formal expression. This can be summarized in the following lemma:
Lemma 4.
x(())=
Σ xEP(m)
& B (X = (log Im
(11)
where 1m (I...., n
17

---

By iterating (10) we can express (11) as
Σ
Σ
(7))
+m=1m
and our goal is to express this sum over partitions. We will use the notation [n]] We can achieve our goal in two parts:
n...
1. Show that for a fixed E Nd with (ig(1) = m we can express (11) as a sum over all surjective functions from [m] to [j].
2. Show that this sum over functions reduces to a sum over partitions.
Part 1. Note that given i' +.... + 1m we may define h: [m]] + [j] by the relation (i*(*)), that is, we take h(n) to be the index for which the multi-index 1 is 1 at n. Note that this function necessarily surjective since the sum taken over non-zero multi-indices. Equivalently, for any surjective function h [m] + [j] we may define multi-indices by setting
{
if E h (c)
(i) 10=
otherwise
Note that any such multi-index will be non-zero since the function is assumed to be surjective, With this identification we can write
(log ((>1))
wi ta x
pt" 16n (7)).
[a[m]-Bb
Part 2. Recall that given a function h [m]] - [5] we can associate it its corresponding partition (P(m) by considering the set (1)...., h (-)), functions corresponding to a given partition, which are given by re-ordering the values I.....J. This reordering of the blocks does not change the summands since the marginals of the partition measure are always copies of each other and hence self-commute, hence a product of moments like
and there are exactly j! different
1t (h- (7²) * ptt +0) (7') can always be written as p (* ), the i-th coordinate of the moment
sequence of the partition measure 7H.- With this in mind we can write
(-1)-1
(-1)
(log)(7'))* 1m
Σ
((())))
Σ
1
EEP(m)
Σ
Σ
Bng (X*#).
=
Ca
TEP(m)
TEP(m)
From this immediately follows that for two probability measures Y 17 we can write
= Σ
E (x0I),
Σ
E (Y no:
TPP(m)
(EP(m)
C. B (X y Hai.
(.TEP(m)
Lemmal then follows from the definition of the tensor products.
C.3 Proof of Theorem and Theorem 3
In this section we present the proofs of Theorem 2 and Theorem We do this in slightly more abstract setting where the feature maps take values in Banach spaces for clarity, until the end when we again restrict our attention to RKHSs. We start out by showing that polynomial functions of the feature maps characterize measures (Lemma 5). From this result we show that cumulants have the same property (Theorem 4), and lastly that this also holds when working directly with the kernels (Proposition 1).
18

---

A monomial on separable Banach spaces Bi....,
.....,, B& is any expression of the form
M ( md) = 1 1y Ea)
for some (i1...., E Nd, where f E B are elements of the dual space B2 and 8; B,,² Finite linear combinations of monomials are called the polynomials. Recall that a set of functions F on set said to separate the points of S if for every x + y E there exists E such that f(x) # f(y).
Lemma 5 (Polynomial functions of feature maps characterize probability measures). Let x.....; Xa be Polish spaces, B Ba separable Banach spaces and şi 24 + B; be continuous, bounded, and injective functtions. Then the set of functions on the Borel probability measures P (Ti= of
II X
/
II X3
+R,
7H
where p ranges over all polynomials, separates the points of (IT-1
Proof. We first show that the pushforward map
d II ps: P
d II
II X, - P
B
injective. This done two parts, first we show that every Borel measure on II X; is Radon measure, then show that the pushforward map injective on Radon measures. To the first part, note that since Xi...., X. are Polish spaces, Theorem 2.5.7; Willard 1970, Theorem 16,4c), and since Borel measures on Polish spaces are Radon measures (Bogachev, 2007, Theorem 7.1.7), any y P(II X.) must be Radon measure.
their product space I X4 (Dudley 2004,
For the second part, that
II
P: II
X
II Bi i=1
II Pi (si,...,,
II (a(a))
is a norm bounded, continuous injection. Since IT B, is Hausdorff space, I is homeo- morphism on compacts since continuous injections into Hausdorff spaces are homeomorphisms on compacts (Rudin, 1953, Theorem 4.17). Let u, (IT2 be two Radon measures such that
their pushforwards are the same I[ (1) = I[ (i((), then for any compact C C I X; we have naC) v(C) as He 9i C II (C((C)) is homeomorphism. Since Radon measures are characterized by their values on compacts, this implies that μ v. Hence the pushforward map injective.
Denote by K the image of I X; under the mapping IT P in IT B, Note that K is bounded Polish space. It is enough to show that the polynomials separate the points of P(K). To see this, note that the polynomials form an algebra of continuous functions that separate the points of I B, and when restricted K they are bounded, since K norm bounded. Since K is Polish, any Borel measure is Radon, and we can apply the Stone-Weierstrn Weierstrass theorem for Radon measures (Bogachev, 2007, Exercise 7.14.79) 79) to get the assertion.
In what follows will use the following index notation for linear functionals. Fix some tuple i= (f.....,id) E Nd with deg(i) = m. Given separable Banach spaces B...,, B3 we use the notation
& BBg
'These monomials naturally extend the classical ones.
19

---

and given an element æ = (z1...., za) I B, we write x': x1 & ... @ x so that x! E B If we have functions (pp:) such that &; - Bi on some Polish spaces ..,..., .. Na.
then we write
w II X, - B
pd
Given collection of linear functionals F [I (B)) such that F = ffa we write
Fsi:= fi & ..... Sa.
p (Bai)*.
Note following trick: the monomials on I B, are exactly functions of the form
x (F
for (fi,..., Sa). this will be used in the proofs. We can now restate and prove the our theorem. Note that the cumulants here are defined like in Definition 4 which is sensible definition even if the feature maps are not associated to kernels.
Theorem (Generalization of Theorem 2 and Theorem 3). Let X,,..., Na Polish spaces and P; : X, + B; be continuous, bounded and injective feature maps into separable Banach spaces B; for = 1....d. d. Let y and r be probability measures on X₁ x ... x Ns. Then
1. if and only if n(()) = =(1])-
2. ii 7xx; if and only if the cross cumulants vanish, that is n5( 0 for all i E N4.
Proof.
Item 2: We want to show that the cross cumulants vanish if and only if - 1/2, By Lemma
5 it enough to show that
(())(())........pe (xa() = o xla, xxx(
for any monomial function p. Let us take linear functionals F (fi....,fa) and note that
(F))
Σ B,
xEP(d)
which is the classical cumulant of the vector-valued random variable
((f ox))())), (fa (Xa))
where (X Xa) ~ 7. Hence by classical results (Speed, 1983), all cross cumulants of ((1 .., (fa (x)) vanish and only if the cross moments split, that to say
()))
Ey (p(( (x))(X)), ..... (Ja pa)(()) = EE nlx, [p(Ci (i)(()), (Sa pa)
for any monomial p dRd. Since fi.... fa were arbitrary this holds for all monomials, which shows the assertion.
Item 1: By assumption n²(n) n1( for every i N4, this implies that E,p(p Pa) (Pa) for any polynomial p. so we can apply Lemma 5.
Proposition 1 (Theorem 2 and Theorem 3). Let Xi...., Xa be Polish spaces and k; X? R be collection of bounded, continuous, point-separating kernels. Let y and n be be probability measures on Xi x x Xa. Then
1. ? 11 if and only if lik. (ek....ka.??) Eke....ks (1).
= 8 d 7xx, if and only if nd. ki.....ke (7) = 0 for all i E No.
20

---

Proof. We reduce the proof to the checking of the conditions of Theorem 4. Let yi denote the canonical feature map of the kemnel k, and let B,: = Hs, be the RKHS associated to k; ( (1. (1.....,d)) dd). For all i E (1...., d). p4 is (i) bounded by the boundedness of k; since [|se((=))) (x, x) I supzex, (x,, < co, (ii) continuous by the continuity of kę (Stein- wart and Christmann, 2008, Lemma 4.29), (iii) injective by the point-separating property of k; The separability of Hi; follows (Steinwart and Christmann, 2008, Lemma 4.33) from the separability of X and the continuity of k; (i E {1.....d) dd). Note: Details on the expected kernel trick part of Theorem and Theorem 3 are provided in Section
.........
D Additional experiments and details
Here we give additional details on the experiments that were performed, and discuss some further experiments that did not fit into the main text.
Background on permutation testing. Permutation testing works by bootstrapping the distribution of a test statistic under the null hypothesis. This allows the user to estimate confidence intervals under the null, which is powerful all-purpose way of doing when analytic expressions are unavailable. As an example, assume we have two probability measures 7, y on X' with i.i.d. samples Fi...., NN 7 y1. .... yN m7. If the null hypothesis is that y 7 then may set
(z1, 22N) = (#1, NN ,I....,,''' yN)
so that for any permutation a on 2N elements, we get two different set of of i.i.d. samples from y by using the empirical measures
Ta : (Zo(1):..... Ro(N)), ;o (Zo(N+1)..... (0(2N))
and for any statistic S P(X)² + R. may estimate S(y, under the null by sampling from S(%,, (z). If the null hypothesis were true, we might expect S(7, 77 to lie in a region with high probability of the permutation estimator, and fairly weak assumptions, this yields a test at the appropriate level (Chung and Romano, 2013).
this as criteria for rejecting the null. Under
Comparing a uniform and mixture. Any uniform random variable over symmetric interval will have O mean and skewness, so symmetric mixture only needs to match the variance. X is 50/50 mixture of Ula, b) and Ul-a, bb then
Var(X)
3 (b² + ba + d²)
so if Y is distributed according to Ul-c, c) then we only need to solve
b² + ba + a = c²
which is straightforward for a given a and c.
The V-statistic for d(2) as written Lemma 2 is bottlenecked by the matrix multiplications. We may note however that for two matrices A, B it holds that
Computational complexity of estimators.
T o
where (.) denotes the sum elements and denotes the Hadamard product. We also note that for for H,, 11,1, we have (AH,),, _ Σ. Ar..- Using both of these tricks may compute both d and CSIC without any matrix multiplications, which brings the computational complexity down ( for both. For a comparison of actual computation time, see Fig. and Fig. where the average computational times for out methods are compared to the KME and and HSIC for N between 50 and 2000.
Type I error on the Seoul Bicycle data. The results when comparing the winter data to itself is presented in Fig. 8. As we see the performance similar for both estimators and lies between 5 and 10%.
21

---

Average (s)
025
820
250
500
750
1002
1203
2500
1750
2000
Figure 6: Average computational time in seconds for KME (red) and d( (blue) for sample size N between 50 and 2000.
Average (9)
03
0.
Figure 7: Average computational time in seconds for HSIC (red) and CSIC (blue) for sample size N between 50 and 2000.
Classical vs. kernelized cumulants. Using the same distributions as in the synthetic independence testing experiment, we now compare X with Y to contrast independence testing with classical cumulants with their kernelized counterpart. The results are summarized are displayed as the median value + half the difference between the 75th and 25th percentile, We consider every combination of classical vs. kernelized, variance vs. skewness, and two different sample sizes. One can observe that the classical variance based test performs poorly compared to classical skewness test, the kernelized variance test almost as powerful as the kernelized skewness test, and in all cases the kernelized tests deliver higher power.
Table I where they
E
Kernel trick computations
Here we show how to arrive the expressions used for the V-statistics used in the experiments.
20.0
150
Type
25
00
72
Figure 8: Type errors using MMD (red) and (2) (blue) on the Seoul bicycle data set.
22

---

Table 1: Comparison of classical and kernelized cumulants for independence testing with both variance and skewness.
N-20
Variance
Skewness
N=30
Variance
Skewness
Classical
19% + 3.0%
56% = 3.5%
Classical Rbf kernel
17% + 0.5%
68% + 1.0%
Rbf kernel
39% 4.5%
59% = 3.0%
65% + 3.5%
79% 1.5%
Given a real analytic function [(x, ................................ sd)) = =ierge Jixt in variables with nonzero radius of convergence and Hilbert spaces H1, Hd we may (formally) extend f to a function
f0 II H. T, (e)(
Fd)
II J IEN
Moreover, if the Hilbert spaces are RKHSs then we have the following result.
Lemma 6 (Nonlinear kernel trick). For any collection of RKHSs Hi...., .Ha with feature maps PP: : X, - H4, assume that and g are real analytic functions with radii of convergence r(f) and r(g) such that maxigisasup.er [p((z)) min(r(f), Then
(fs ((())) (
Σ iEN
figkt "1) ... ka(Fa, ya)'
Proof. Since the image the yas lie inside the radius of convergence of fer and gse the power series converge absolutely and we can write
(Se (v0 (x))) ge0 ( 8i ((')))
- Σ (1r))
Σ
0i (s²))T
ieNs
iENd
Σ
fi 11 |
Σ
(19 y1)" ka(aa,
IEN4
IEN4
where H. = 7/1 x ... x Hd-
Using Lemma 6, we can choose kernels k; : X?2 - with associated RKHSs H, and feature maps P and some i Ne with deg(i) m. We make the observation that with X (. Xa) ~ 7 Y .Y....., Ya) ~ v and g and अ in (4), one has
(x), (1))
Σ (mP(m)
Co B
s1 (X)).
Σ TEP(m)
cy s (Y*))
Σ
C (Ent (
(x)), E2y (Y+))
Σ
(x))) 10 ((l())))) ))
x,TEP(m)
Σ
G E Eeeen 1 ((),.... Xm), (Y,,..,, Vnn),
Since
[T)( =
[]x(()) -0 (m))l (wt()), n² (T)) (xd (1)) 2(k²(7), ot (x))
one gets the expected kernel trick statements of Theorem and Theorem 3.
We are now interested in explicitly computing the expression [| (7) H 2 || (2) (7)
(2) (1))) and || (7) (3)())) nik (n)ll and their corresponding V-statistics. Recall that for
Nok
23

---

a (w.l.o.g.) symmetric, measurable function h(2 2m), the V-statistic of h with N samples Z....., Nn is defined as
V(h; L1,.................... ZN) N
Σ
MZZ
Zim) ).
Under fairly general conditions, the V-statistic converges in distribution to E|A Zl.... Zm)] 2i,...,,m.. well-developed theory describes this convergence (Van der Waart, 2000; Serfling, 1980; Arcones and Giné, 1992).
and
Example E.1 (Estimating (r) nee (m)l Let X, X, ,",,""" X denote independent copies of and Y, Y", Y Y denote independent copies of n. The full expression for 1|=² (7) pa is Pik
(2)
(2) -Kk 1 1 (7) - Ki (7)1 = XX(X,X'))(X" X')(XX,,, xX") + BR(Y, X')((" + E&(X, )² + EX(Y,
(12)
+ 2EA((X Y { (xXXxXX)))) YY + 2EX ()((X,
-EEA(X, () 2EA(X, Y)²
2EE(X, X'),(XXX")) 2E (")((Y
Given samples (x:)
(1.) from y and respectively the corresponding V statistic is
1
Σ
1
Σ
k(yi, My ((NNs m1)
(13)
N
MT
k(a, 22)
+
k(yi,
N2
M
1==1
N
M
N
M
Σ
k(ri, My (K( M1) +
k( Mo )k(xi, Mn)
NºM
NM
1,k=1 j=1
i=1
j.K=1
N
N
M
Σ
Σ
k(a,,1;; (()
2
k( n3)
N
NM
i11= 1===
2
Σ
kly, M1 M.(e
N3
M3
$jA==
f.j.==1
Let us define the Gram matrices K, [k(a;, + E NXN K,, [k(yi, ((((),;)) E MMM,
Ky,y = [k( U;)] j-j=1 centering, then (13) can be rewritten
N.M. and let H = N IN IN E RNN,, HM = las E RMX be the
1
1
Tr(H K H,K,) +
Tr(H K K w) +
N2 Tr(K3) + M2 Tr(K3)
N2
M
2
2
2 NM
Tr(K zy HN K (q))
Tr(K Hy
Tr HT
Tr(K2,) my)
NM
NM
NM
2
Hr.))
Tr H,K,,
N2
M
which simplifies to
2
N2 (N-H)) M²
Tr [K[ H
Tr [K By (I H HM)KJ by II- Hx)]
NM
This estimator can be computed in quadratic time.
Example E.2 (Estimating [| (1) 22 ( wxze)) Let denote the kernel on X1 and & denote the kernel on X2. Let (X,Y), (X', Y"), (X", Y"), (x((), (3)), (x()) (())), (x(5), (()) denote in-
24

---

dependent copies of E P(G x X2). The full expression for 1 (1.2) (7) 2001 60 is
E&(X, (X), X -EA, X")(XXX))
v(O1)
+ +X)
(e)))) (()))
+ 4X, X X
y(ai)
SIEA(XX, X')(("", )E((()))))) Y) (EK X')(X", XX) X')e(r(o) y(o)
+ LIEA(XXX))ELL .X (e(r)) y(s)).
Given samples (Ii, yt) i from y the corresponding V-statistic for this expression
N
N
M(r
Σ
20), U1)
2
N
=
iJM=1
2
Σ
k( ((1,,
Σ
k( 2, w)
N3
N4
2
vi) +
k( 2 (K) (1) 21) v))
N
N
sja,l=1
1.j,s,111
k( (2,, vn) +
k(a,, xz (x), m1)
N
NA
1,,,11=11
k(
Σ
k(1;,
(8 91)
N
N5
Vn).
N6
Using the shorthand notation K = K., Ly and H Hn and denoling by the Hadamard product [A
BJs, BB,, A and (-) the sum over all elements of a matrix (A) Eig=1 Aig. the V-statistic above can be written in the simpler form
K - Ko 4K o KH - 2K o K o LH
N2
K
+ 4KH K KOLH 2K o L
+ 2KHH o HK L
N
L
K N2
+ KK O HK o LH + K K
SK o LH
N2
L
K
2
4K HK
+ 4
L
N2
N2
Again this estimator be computed in quadratic time.
Example E.3 (Estimating [|+e (7) xx (pa)l In order to estimate d(3) can write
(3)
(7,7) 7) note that one
|| (3) (7) 1x
(3) (m)) = ||m (3) (7) 72 -11 (3) (n) 7000
2/m (7), 1ik (n)) (1))>>>>
(3)
25

---

We can estimate the first two terms like in Example E.2, and the third term can be expressed as (3) (3) (Ki (7). ME (g) (X Y) 3Ek(X, k)xk(X
3IE((X, ()), + 6BK(X,
+ 3 Y)²k(X", Y") 2EK(X, F)k(X",)) + (X (),(X, Y"),(X 6EA(X, ()K(X, ()R(X, Y") 6IEK(X, ()((X', ()((X", Y')+ 4 () (')(",, ,Y",.
For simplicity will assume that we have equal number of samples (N) from both measures The V-statistic for (K& (1), (1), Hi (3)(()) (n)) can be expressed as
(xra) ) N
E Y and (3)) Eq.
N
1
3
Mas, yj) (1)
22
N3
fJ=1
ij.mml
N
N
3
6
Σ
k( Uj (= y/4
Σ
k( y3 (a (a),
N3
N
1.A==1
N
3
k(x; yj) (k(xxx) vt)
+
k( 3j (2) U,)
N4
N
i.A.ll=1
N
N
k(r, k H.)
Σ
k(r; Vs (( Ue )a(a mm)
N4
No
1.j.8,I.m=1
N
k( U1 ( ) Mm)
k( (F) (a vn).
N
N76
k.j.d.mm11
k...J.mmmm1
Using the notation K
= [k( Uo ) this estimator simplifies to
1
K Ky Kay 3K.y Kry HK xy
N22
3K o Kay o K H + 6K my o ay
Ho HK
ay
Kay N2
+ 3Kay o Kay
+ +2K
HK
HK ay
K ay
+ 2K o Kay Ho Kzy
(K o K
H
N2
K
2
- my o HK
4
K
N2
N2
We mention also that the first two terms
|| (3) (1)
Hesn 11
(3)
(1))l
be computed little
more simply than in Example E.2 since the expressions have more symmetry, using the notation
K, [k[2,, x1 2))l we can write down the V-statistic for || [?) H as
(3)
1
K2 K, o K, 6K, K,HooK,
N2
K
44, K2 K,H 3K, K,
N2
K
66, H HK v K2 12K HK,
N22
K
+4
K
N2
with a similar expression for ||
(1)
2
The estimator can be computed in quadratic time.
26


Reviews:
[{"review_id": "C2JbdRwoe0", "reviewer": "Reviewer_Q1mu", "summary": "The paper generalises the notion of kernel mean embeddings to higher-order cumulants. It proposes kernelled cumulants in the RKHS. While kernelled cumulants reside in tensor product space of the RKHS, the paper shows that Hilbert space metric between cumulants can be exactly computed using the kernel trick. Based on this construction, the paper proposes: (1) a two-sample test statistics, that generalises MMD test statistics by considering distance between cumulants, and (2) a generalisation of HSIC statistic for independence testing, again by considering distant between cumulants of joint distribution and product of marginals. \nThe advantages of the construction and proposed tests include: (i) the test are applicable for a broader class of \"point-separating\" kernels (unlike MMD/HSIC that can only useful for characteristic kernels); (ii) the new statistics can be computed in quadratic time (same as MMD); and (iii) empirically achieve higher power that classical MMD/HSIC statistics (both on synthetic and real data).", "questions": "- see weakness\n- in addition, a precise statement on computational complexity (at least for d2, d3) would be useful\n\nMinor remarks: H in Lemma 3 is not defined in main paper (but in appendix), and E et al citation seems incorrect (there is a full name for E)", "limitations": "the paper does not have immediate negative societal impact (although conclusions from hypothesis tests can always have). Hence, the work could benefit from:\n- consistency results (similar to kernel two-sample tests)\n- characterisation of whether higher-order cumulant based statistics typically tend to be larger than MMD (even under null). The comment is about sample estimates and not expected value (hence, more tied to concentration/consistency of test statistics)", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "- The idea of considering higher order moments/cumulants in RKHS is quite natural, and yet unexplored in the literature (apart from the recent work of Makigusa (2020) that consider only 2nd order moments). Hence, the contribution is novel and quite timely\n- The main strength of the paper is that the computational cost of the proposed statistics is still quadratic in sample size (same as MMD), which implies the advantages of higher order moments does not come at significant additional cost\n- The work is technically sound, and the construction of cumulant and use of kernel trick is reasonably involved.", "weaknesses": "- While it is easy to imagine in general that higher-order cumulants can distinguish between more distributions, the advantage of kernelled cumulants is difficult to grasp. If a characteristic kernel is used wouldn't mean embedding (MMD) suffice?\n- The line of argument used in the paper to demonstrate advantage of kernel cumulants is that (empirically) they have show higher power (can detect small differences better). Is there a theoretical justification for this? It would be sufficient if the authors provide a justification / reference that standard (non-kernel) cumulants are more sample efficient in some cases, where means already show separation\n\n-  For the synthetic experiments, the null rejection rate should also be plotted to show that cumulants do not have higher tendency to reject than MMD/HSIC. While this is true for the real data, unfortunately both 1st order and higher order terms reject at a rate higher than significance level\n\n- Overall, it is not clear when tests based on higher-order cumulants are indeed needed in practice. I still feel the work is relevant, but some discussions on this would certainly increase the significance of the work for the broader community \n- The paper, although well-written, is quite dense and at times bit difficult to follow, but this can be attributed to the content of the paper", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time and effort invested, and for the kind review. Below, we answer the questions in detail.\n\n - ”If a characteristic kernel is used wouldn’t mean embedding (MMD) suffice?”\nThis is true of course, a characteristic kernel is theoretically sufficient for 2-sample testing (in MMD), and a universal one\nis sufficient for independence testing (in HSIC). The main use case of cumulants are for when either one is in a situation\nwhere a non-characteristic kernel is preferable, or one does have a characteristic kernel but the data is structured in a way\nwhere one gets additional test power from using the cumulants and hence estimators which are more sample-efficient, like\nin the experiments shown in the main text.\n\n - ”Is there a theoretical justification for the higher testing power?”\nTesting power for these kinds of tests do not normally come with satisfactory descriptions since the tests are infinite dimensional. Closed form expressions of power are usually expressed as infinite convergent series, see for example Theorem 3.1 [16]. because of this, trying to study theoretical properties of testing power in this setting is a very challenging problem, and while an interesting one, we believe that the experimental results are convincing, and that further theoretical study is beyond the scope of the article.\n\n - ”both 1st order and higher order terms reject at a rate higher than significance level”\nIt is true that the rejection rate is higher than 5%, but as it is consistent for both methods we do not feel that it gives any\nunfair advantages either way. The reason for this is most likely due to the way we choose the hyperparameter by optimising\nfor test power, which would incur a bias for the rejection rate. We did include rejection rates for the bicycle data in Figure\n8 in the appendix where it can be seen that it is around 8% for both methods, and we can include the rejection rates for\nthe synthetics experiments here too.\n\n - ”(. . . ) when tests based on higher-order cumulants are indeed needed in practice. I still feel the work is relevant, but some\ndiscussions on this would certainly increase the significance of the work for the broader community” \nOne answer is that kernelized cumulants provide guarantees even when non-characteristic kernels are used. However, even for characteristic kernels the presented kernelized cumulants are advantageous since they lead to more sample-efficient estimators and hence tests. To see this on the simplest ”unkernelized” example would be to test if two normal distributed random variables\nhave the same distribution. Clearly, here the best answer is to just compute the sample means and sample covariances\nand then reject/accept if they are close enough; in contrast computing mean and second moments would lead to much\nmore noisy tests and more samples would be needed (see the discussion in Appendix A). The same arguments carry over\nwhen we use kernels for testing. Our experiments confirm this intuition and show that the tests given by statistics derived\nfrom kernelized cumulants are usually more efficient than the classic kernelized statistics. Hence, kernelized cumulants\nare usually preferrable. However, the cost of the somewhat more complicated test statistics get harder to implement, the\nhigher m becomes. On a practical side, we believe using m = 2 is a good compromise between gain in sample efficiency\n(hence better tests) and ease of implementation. We will add a sentence to emphasize this.\n\n - ”A precise statement on computational complexity”\n\n**Sample size dependence:** The exact statistic used are expanded on in examples E.1, E.2, and E.3 in the appendix. Since\nall these statistics can be phrased as polynomials of Gram matrices and some multiplication with centering matrices, the\nmost expensive computation is computing the Gram matrices themselves. For example, when computing $d^{(2)}$, assume that\none has N samples from one measure and M from the other. One first computes the 3 relevant Gram matrices which is\n$O(N^2 +NM + M^2)$ and depends on the specific kernel used. After that one computes the centered matrices which involves\n$N^2 + 2NM + M^2$ additions and the same number of divisions. One then computes the Hadamard square of the centered\nmatrices which involves $N^2 +NM + M^2$ multiplications and finally sums over all the elements and normalizes the output\nwhich involves a further $N^2 +NM + M^2$ additions and 3 divisions, and then 3 more additions to add them up. Taking the\ntotal number of operations (in addition to computing the Gram matrices), this means $2N^2 + 2M^2 + 3NM + 3$ additions\nand the same number of multiplications/divisions.\n\n**Degree dependence:** There is no free lunch. The $m$-th Bell number $B_m$ (https://en.wikipedia.org/wiki/Bell_number)\nis defined as the number of elements in $P (m)$. The Bell numbers follow a recursion: $B_{m+1} = |P(m+1)|  = \\sum_{k=0}^m \\binom{m}{k} B_{k}$, with  the first elements of the sequence being $B_0 = B_1 = 1$, $B_2 = 2$, $B_3 = 5$, $B_4 = 15$, $B_5 = 52$, $B_6 = 203$, $B_7=877$, $B_8=4140$. By (6)-(7), in the worst case the number of operations to compute $d^{(i)}(\\gamma,\\eta)$  or $\\Vert\\kappa^{i}_{k_1,\\ldots,k_d}(\\gamma)\\Vert^2$ is proportional to $B_m^2$ (it equals to $3B_m^2$ and to $B_m^2$, respectively). It is true that asymptotically $B_m$ gets very large [4, 12], but\nfor reasonably small degrees the computation is still manageable. In addition, merging various terms in the estimator can\noften be carried out, which leads to computational saving. For instance, the estimator of $d^{(2)}$ (see Lemma 2, Example E.1),\nCSIC (Lemma 3, Example E.2) and $d^{(3)}$ (Example E.3) consists of only $2$, $11$ and  $10 + 2\\times 7 = 24$ terms compared to the\npredicted worst-case setting of $3B_2^2 =12$, $B_3^2 = 25$, and $3B_3^2=75$ terms, respectively.\n\n - Minor remarks: \nThank you for pointing these out, we will address them.\n\nWe hope that our answers resolve all the questions.\n"}, {"review_id": "p0lk2idYXs", "reviewer": "Reviewer_Lvaf", "summary": "This papers revisits advances in cumulants on real data, by extending them to provide cumulants for random variables in an RKHS.", "questions": "see question in weaknesses", "limitations": "OK", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "The idea of this work is interesting, as the paper proposed to go beyond conventional kernel mean and kernel covariance. Moreover, a proposed kernel trick allow to obtain efficiently the kernelized cumulants.\nThe paper well describes the contributions. The derivations seem to be good, and some experiments on synthetic and real datasets allow to understand the relevance of the proposed approach", "weaknesses": "While the paper includes some experimental results, there are only two settings : independence test and two-sample test (MMD-like). Is there any other  tests (or beyond test) where the proposed kernelized cumulants would be relevant ?  ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time and effort invested, and for the kind review.\n\n - ”Is there any other tests (or beyond test) where the proposed kernelized cumulants would be relevant?”\nThe proposed kernelized-cumulant based measures are general divergence and independence measures, hence can be used\npractically in any application which rely on information theoretical objective and expectedly with improved sample-\nefficiency (as our hypothesis testing illustrations in the submission suggest). Examples include feature selection [2, 18, 9, 22],\ncausal discovery [14, 13, 16, 3, 17], distribution classification [11, 23] and regression [20, 19, 8, 6, 15], or generative adversarial\nnetworks [5, 10, 1].\n\nWe hope that our answer resolves the reviewer’s question."}, {"review_id": "WrecoWvoJZ", "reviewer": "Reviewer_952h", "summary": "The authors introduce the kernelized cumulant and show that it can characterize distributions and statistical (in)dependence.", "questions": "1. Is the independence criterion kernel dependent? What happens if the kernel is misspecified?\n2. What are the convergence rates for the finite-sample estimators?", "limitations": "N/A", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The kernelized cumulant provides a natural generalization of the popular maximum mean discrepancy (MMD) as well as the Hilbert-Schmidt independence criterion (HSIC).\n2. The authors illustrate several interesting properties of the kernelized cumulants and provide a two-sample test for non-characteristic feature maps as well as an independence test.\n3. The authors demonstrate the utility of the proposed method on a variety of datasets with competitive results.", "weaknesses": "1. Estimating the kernelized cumulant may require a large sample. Convergence rate of the estimator is not discussed, but the empirical results suggest comparable performance as the HSIC.\n2. The proposed method could potentially be vulnerable to kernel misspecification.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time and effort invested, and for the kind review. Below, we answer the questions in detail.\n\n - ”Is the independence criterion kernel dependent?”\nSince we are working in a very general framework, the independence criterion works for a very wide family of kernels,\neven for simple ones such as the linear kernel on $\\mathbb{R}^d$. However, our proposed methods suffer from the same drawbacks\nas independence testing with HSIC or 2-sample testing with MMD in that the testing power can be severely affected by\npicking an inappropriate kernel.\n\n - ”Convergence rates”\nAs the estimator of our proposed kernelized cumulant is a sum of V-statistics, the rates are inherited from that of V-\nstatistics.\n\n - ”kernel misspecification”\nIn our analysis we imposed minimal assumptions on the kernel (point-separating) to allow for the wide applicability of\nthe proposed framework. Understanding the optimal choice of kernels (we assume that this is what the reviewer means\nby kernel misspecification) or even the involved hyperparameters specifically for MMD (which corresponds to a degree one\nobject in our work) in specific downstream tasks is an important and highly non-trivial problem. Probably the simplest\ntask in this context is 2-sample testing, where one can phrase the goal to achieve minimax optimality. In this very specific\ncase, the analysis for constructing (almost) minimax optimal MMD-based adaptive tests can be worked out in 63 pages\n[7]. Extension of their results could mean the first step of the grounding of kernelized cumulants in downstream tasks.\nSimilarly, getting optimal rates even for the classic MMD and for radial universal kernel is non-straightforward; to our best\nknowledge, the only available result in this domain is [21]. Our investigated assumptions are weaker (hence the analysis is\nexpected to be harder) in multiple aspects: the space is Polish (hence even the notion of radial property is not defined),\ncumulants need not be 1st degree objects (as MMD), and our results (Theorem 2-3) hold for point-separating kernels\n(which requirement is significantly weaker than characteristic property, which itself is a specific case of universality).\nWe agree (in line with the penultimate paragraph of Section 1) that these are important future research directions.\n\nWe hope that our answers resolve all the questions."}, {"review_id": "vkFoCjQqas", "reviewer": "Reviewer_h8uP", "summary": "This paper generalizes the notion of cumulants to Hilbert-space-valued random variables. \nWhen these Hilbert spaces are RKHSs, the kernel trick applies so that computations can be performed with the kernel function.\nIt leads to higher-order two-sample and independence tests, which generalizes MMD and HSIC.\nThe efficiency of these tests is demonstrated numerically on both synthetic and real data.\n\nEDIT: I have read the author's rebuttal, which partially addressed my concerns.", "questions": "- I am suspicious of the fact that the computational complexity of MMD/HSIC and the higher-order extensions are the same, and are bottlenecked by the computation of the kernel matrix. Surely the complexity must increase with $m = \\mathrm{deg}(\\mathbf{i})$, even if it remains quadratic in the sample size?\n\n- Is there a relationship between going to higher-order cumulants and changing the kernel? Can we achieve the same performance than the higher-order cumulants by adapting the kernel?\n\n- Isn't the the $V$ statistic estimator a straightforward replacement of the expectation with empirical averages over the sample? I think this could be mentioned in the text.", "limitations": "The authors mentioned the two main limitations of their approach:\n- It is not clear how to design the kernels to maximize the performance of the statistical tests.\n- There is no theoretical analysis of the introduced tests.\n\nI agree that a complete resolution of these issues should be left to future work. However, I think the second weakness above, which is related to the first limitation, should be at least acknowledged.\n\nThere is no foreseeable negative societal impact.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well written and clear, despite the fact that the notation is heavy due to the complexity of the considered objects.\nI appreciated that equivalent definitions of the kernelized cumulants are given, as well as special cases to build intuition and relate it to classical notions.\n\nBecause of this, I think the paper should be accepted.", "weaknesses": "A first weakness in my opinion is that the main contribution of the paper is a straightforward combination of two known concepts: (i) cumulants and (ii) the use of the kernel methods in statistical testing. Also, rather than using the tensor algebra framework used by the authors in Section 3 and Appendix C.2, it seems to me that the simpler route would be to use the generating function $K(\\theta_1, \\dots, \\theta_d) = \\log \\mathbb E[ \\exp(\\sum_{i=1}^d \\langle \\theta_i, X_i \\rangle ) ]$ and define the cumulants from its series expansion, or equivalently derivatives at zero (which would be the correct equivalent of item 1. in Appendix C.1, rather than eq. (11)) . The definitions are of course equivalent, but it would lighten the formalism in the main text.\n\nThe second main weakness is that the relationship between MMD/HSIC and the higher-order versions *when one is allowed to change the kernel* is not discussed. If I'm not mistaken, the second-moment embedding of a probability distribution with kernel $k$ coincides with the mean embedding with the squared kernel $k^2$. So one should expect a relationship between $d^{(2)}$ with kernel $k$ and MMD with some combination of $k$ and $k^2$ as kernel. If that is the case, then the use of higher-order cumulants can be equivalently rephrased as the use of different kernels. In particular, it is not clear whether one could be better off by sticking to MMD/HSIC but with well-designed kernels, as the estimation of higher-order moments has a higher variance. I think this point should be discussed in the text.\n\nTwo additional minor remarks:\n- Repetition of \"energy distance\" in lines 37-38, which I'm guessing is a typo.\n- Missing related work: [1] considers a second-moment kernel embedding and defines kernel information-theoretic quantities.\n\n[1] Bach, Francis. \"Information theory with kernel methods.\" IEEE Transactions on Information Theory 69.2 (2022): 752-775.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time and effort invested, and for the kind review. Below, we answer the questions in detail.\n\n - ”use the generating function”\nWe agree that this is a more intuitive way to introduce the concept, but it is not very instructive when it comes to\nactually computing the statistics or writing down estimators, and in order to do so it is typical to work with combinatorial\ndescriptions of cumulants. In order to save space we introduced only the combinatorial one as it was necessary to derive\nthe statistics in the later parts of the paper.\n\n - ”higher-order cumulants can be equivalently rephrased as the use of different kernels”\nThis is true for higher order moments, but we stress that the cumulants are not equivalent to mean embeddings of\ndifferent kernels. To illustrate this, the embedding for the second cumulant of a random variable $X$ is $ \\mathbb{E}\\big[k(X,\\cdot)\\otimes k(X,\\cdot)\\big] - \\mathbb{E}\\big[k(X,\\cdot)\\otimes k(Y,\\cdot)\\big] $ where $Y$ is an i.i.d copy of $X$. The first term here is indeed a mean embedding of the product kernel\n$k^{\\otimes 2}$, but the second term acts on a different product measure than the first. This structured cancellation of the various\nmoments of the measure is central to the cumulants ability to pick up the higher order features, and makes them distinct\nfrom ”only” using higher order moments of kernels.\n\n - ”Surely the complexity must increase with the degree”\nThis is true, we apologize that this was not made sufficiently clear in the text. The complexity in the sample size is\nquadratic for all the proposed statistics, but the number of operations to compute them does increase with the degree and\nis determined by the number of partitions corresponding to that degree. The number of partitions does grow quickly as\nthe degree gets large, but for reasonably sized degrees it is not an issue. See also our reply to R1 on this topic, particularly\nthe paragraph starting with ’Degree dependence’.\n\n - ”Isn’t the the V-statistic estimator a straightforward replacement of the expectation”\nIt is very similar but there is some subtlety in which indices are summed over when computing the V-statistics, as an\nexample consider the two generalised moments of a random variable $X$, $\\mathbb{E}(X^2)$ and $[\\mathbb{E}(X)]^2$. Given $N$ samples $x_1, \\ldots, x_N$ the V-statistic for the first one is $\\frac{1}{N}\\sum_{i=1}^N x_ix_i$ and the V-statistic for the other is $\\frac{1}{N^2}\\sum_{i=1}^N\\sum_{j=1}^N x_ix_j$.\n\n - Minor remarks: \nThank you, we will address this and add the relevant citation. Regarding the reference, (i) it falls under the umbrella of\n$k^{\\otimes 2}$-based information theoretical quantities, (ii) it relies on the uncentered covariance operator, (iii) it requires the kernel\nto be universal (in contrast to our significantly relaxed point-separating assumption), and (iv) its complexity (justified in\nthe submission for the case of kernel entropy) is high, cubic in the sample size (whereas our estimator for fixed degree is\nquadratic in the sample size). We are happy to include a citation to the mentioned paper as a related work.\n\nWe hope that our answers resolve all the questions."}, {"review_id": "AKPr9MmfJH", "reviewer": "Reviewer_6dmY", "summary": "This paper proposes kernelized cumulants to extend classical cumulants in $\\mathbb{R}^d$ and shows that the kernelized cumulants provide a new set of all-purpose statistics and are computationally tractable. \nThis paper also show advantages of kernelized cumulants both theoretically and empirically.", "questions": "1. In the experiments, what is the criteria of the 'optimal value' of $\\sigma$ (Line 265)?", "limitations": "The authors have adequately addressed the limitations.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "* The paper is well-written and easy to follow.\n\n* The method, that is, kernelized cumulant, proposed in this paper is novel and might be useful in real-world applications.\nFor example, this method can be used to provide metrics for distributions and measure independence.\nMoreover, this method includes traditional MMD and HSIC as its special case and can even outperform the traditional ones.\n", "weaknesses": "* In the experiments, more kernels, such as the neural tangent kernel, could be considered.\n* It would be better if the authors provide theoretical guarantees of the effectiveness of kernelized cumulants, such as consistency.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the time and effort invested, and for the kind review. Below, we answer the questions in detail.\n\n - ”more kernels, such as the neural tangent kernel, could be considered.”\nYes, essentially any kernel can be used and in our experiments we focused on the standard kernel choices by practitioners.\nWe are not sure what the reviewer exactly means by using the neural tangent kernel in this context; it would be certainly interesting to\nuse the presented kernelized cumulants to study how dependencies arise in the training of neural networks, although this\nis beyond the scope of the current work.\n\n - ”provide theoretical guarantees”\nThis is a good point, as the estimators are V-statistics they inherit all of their nice properties like consistency, so this boils\ndown to citing the standard references for V-statistics. We will clarify this in the main text.\n\n- ”what is the criteria of the ’optimal value’ of $\\sigma$?”\nBy ’optimal value’ we mean the value for which the test has the highest power; we will clarify this in the main text.\n\nWe hope that our answers resolve all the questions."}]

Author rebuttal:
**References for the rebuttal:**

[1] Mikolaj Binkowski, Danica Sutherland, Michael Arbel, and Arthur Gretton. Demystifying MMD GANs. In International
Conference on Learning Representations (ICLR), 2018.

[2] Gustavo Camps-Valls, Joris M. Mooij, and Bernhard Sch ̈olkopf. Remote sensing feature selection by kernel dependence
measures. IEEE Geoscience and Remote Sensing Letters, 7(3):587–591, 2010.

[3] Shubhadeep Chakraborty and Xianyang Zhang. Distance metrics for measuring joint dependence with application to causal
inference. Journal of the American Statistical Association, 114(528):1638–1650, 2019.

[4] N. G. de Bruijn. Asymptotic Methods in Analysis. Dover, 1981.

[5] Gintare Karolina Dziugaite, Daniel M. Roy, and Zoubin Ghahramani. Training generative neural networks via maximum
mean discrepancy optimization. In Conference on Uncertainty in Artificial Intelligence (UAI), pages 258–267, 2015.

[6] Zhiying Fang, Zheng-Chu Guo, and Ding-Xuan Zhou. Optimal learning rates for distribution regression. Journal of
Complexity, page 101426, 2020.

[7] Omar Hagrass, Bharath K. Sriperumbudur, and Bing Li. Spectral regularized kernel two-sample tests. Technical report,
2022. (https://arxiv.org/abs/2212.09201).

[8] Ho Chung Leon Law, Danica Sutherland, Dino Sejdinovic, and Seth Flaxman. Bayesian approaches to distribution regression.
International Conference on Artificial Intelligence and Statistics (AISTATS), 84:1167–1176, 2018.

[9] Runze Li, Wei Zhong, and Liping Zhu. Feature screening via distance correlation learning. Journal of the American
Statistical Association, 107(499):1129–1139, 2012.

[10] Yujia Li, Kevin Swersky, and Richard Zemel. Generative moment matching networks. In International Conference on
Machine Learning (ICML), pages 1718–1727, 2015.

[11] David Lopez-Paz, Krikamol Muandet, Bernhard Schölkopf, and Iliya Tolstikhin. Towards a learning theory of cause-effect
inference. International Conference on Machine Learning (ICML), 37:1452–1461, 2015.

[12] Laszlo Lovasz. Combinatorial Problems and Exercise. 2nd ed. Amsterdam, Netherlands: North-Holland, 1993.

[13] Joris Mooij, Jonas Peters, Dominik Janzing, Jakob Zscheischler, and Bernhard Schölkopf. Distinguishing cause from effect
using observational data: Methods and benchmarks. Journal of Machine Learning Research, 17:1–102, 2016.

[14] Krikamol Muandet, Kenji Fukumizu, Francesco Dinuzzo, and Bernhard Schölkopf. Learning from distributions via support
measure machines. In Advances in Neural Information Processing Systems (NIPS), pages 10–18, 2011.

[15] Nicole Mucke. Stochastic gradient descent meets distribution regression. In International Conference on Artificial Intelli-
gence and Statistics (AISTATS), pages 2143–2151, 2021.

[16] Niklas Pfister, Peter Buhlmann, Bernhard Schölkopf, and Jonas Peters. Kernel-based tests for joint independence. Journal
of the Royal Statistical Society: Series B (Statistical Methodology), 80(1):5–31, 2018.

[17] Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua
Bengio. Toward causal representation learning. Proceedings of the IEEE, 109(5):612–634, 2021.

[18] Le Song, Alex Smola, Arthur Gretton, Justin Bedo, and Karsten Borgwardt. Feature selection via dependence maximization.
Journal of Machine Learning Research, 13(1):1393–1434, 2012.

[19] Danica Sutherland, Junier Oliva, Barnabas Poczos, and Jeff Schneider. Linear-time learning on distributions with approxi-
mate kernel embeddings. In AAAI Conference on Artifical Intelligence (AAAI), pages 2073–2079, 2016.

[20] Zoltan Szabo, Bharath K. Sriperumbudur, Barnab ́as P ́oczos, and Arthur Gretton. Learning theory for distribution regression.
Journal of Machine Learning Research, 17(152):1–40, 2016.

[21] Ilya Tolstikhin, Bharath Sriperumbudur, and Bernhard Schölkopf. Minimax estimation of maximal mean discrepancy with
radial kernels. In Advances in Neural Information Processing Systems (NIPS), pages 1930–1938, 2016.

[22] Andi Wang, Juan Du, Xi Zhang, and Jianjun Shi. Ranking features to promote diversity: An approach based on sparse
distance correlation. Technometrics, 64(3):384–395, 2022.

[23] Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnab ́as P ́oczos, Ruslan Salakhutdinov, and Alexander Smola. Deep
sets. In Advances in Neural Information Processing Systems (NIPS), pages 3394–3404, 2017.

Meta-review: The reviewers unanimously recommended to accept the paper.
Decision: Accept (spotlight)