Permutation-Free High-Order Interaction Tests
Authors: Zhaolu Liu
Robert Peach
Mauricio Barahona
Venue: ICML 2025 poster (2025)
Keywords: high-order interaction; kernel method; permutation-free
Primary area: general_machine_learning->scalable_algorithms

Abstract: Kernel-based hypothesis tests offer a flexible, non-parametric tool to detect high-order interactions in multivariate data, beyond pairwise relationships. Yet the scalability of such tests is limited by the computationally demanding permutation schemes used to generate null approximations. Here we introduce a family of permutation-free high-order tests for joint independence and partial factorisations of $d$ variables. Our tests eliminate the need for permutation-based approximations by leveraging V-statistics and a novel cross-centring technique to yield test statistics with a standard normal limiting distribution under the null. We present implementations of the tests and showcase their efficacy and scalability through synthetic datasets. We also show applications inspired by causal discovery and feature selection, which highlight both the importance of high-order interactions in data and the need for efficient computational methods.

Permutation-Pree High-Order Interaction Tests
Zhaolu Liu' Robert I., Peach 23 Mauricio Barahona
Abstract
servational data. A growing body of work has approached this challenge using tools from information theory (Rosas et al., 2020; 2024), persistent homology (Santoro et al., 2023; 2024), Bayesian inference (Young et al., 2021) and Hamiltonian-basee functions (Schneidman et al., 2003).
Kernel-basee hypothesis tests offer a flexible, non- parametric tool to detect high-order interactions in multivariate data, beyond pairwise relationships. Yet the scalability of such tests is limited by the computationally demanding permutation schemes used to generate null approximations. Here we introduce family of permutation- high-order tests for joint independence and partial factorisa- tions of d variables. Our tests eliminate the need for permutation-base aging V-statistics and a novel cross-centring tech- nique to yield test statistics with a standard normal limiting distribution under the null. We present implementations of the tests and showcase their efficacy and scalability through synthetic datasets. We also show applications inspired by causal dis- covery and feature selection, which highlight both the importance of high-order interactions in data and the need for efficient computational methods.
Within the statistical and machine learning literature, kernel- based tests provide a flexible, non-parametric framework for identifying relationships in data. The detection of pairwise independence between two variables using HSIC (Gretton et al., 2007) has found broad applicability in tasks such as generative adversarial networks (Binkowski et al., 2018), feature selection (Song al., 2012), and transfer learn- ing (Long et al., 2015). Beyond pairwise independence, kernel-based tests have been extended to assess joint inde- pendence among d variables via dHSIC (Pfister et al., 201 8), as well as Lancaster interactions among three variables (Se- jdinovic et al., 2013a), which were later generalised to d- order Streitberg interactions (Liu al., 2023b). However, these tests typically rely on computationally expensive per- mutation or bootstrap schemes (involving 100 s p S 1000 permutations), reducing their applicability to datasets with larger number of variables, where the subsets of variables to be tested for interactions explode combinatorially.
approximations by lever-
1. Introduction
Many complex real-world systems involve interactions be- yond pairwise relationships, necessitating an explicit under- standing of interactions among three or more entities (Bat- tiston et al., 2020; 2021; Rosas et al., 2022; Thibeault et al., 2024). Such high-order interactions amongst d > 2 vari- ables have been shown to affect processes in social (Benson et al., 2016), ecological (Grilli et al., 2017) and biological systems (Luff et al., 2024; Arnaudon et al., 2022), often in ways that cannot be reduced to combinations of pair- wise effects (Schaub et al., 2020; Bick et al., 2023). Given the relative scarcity of relational datasets that record high- order phenomena, it is important to develop unsupervised, mathematically grounded methods for their detection in ob-
Here, we tackle this limitation on scalability by introducing family of permutation-fre high-order tests that retain the flexibility of kemel-based methods while eliminating the computational burden of null-distribution approximation. To note, our approach applies to the general settings where one seeks to determine factorisation properties among sub- sets of variables (e.g., Lancaster and Streitberg interactions), rather than merely rejecting or accepting joint independence.
Our main contributions can be summarised as follows: (1) We define the permutation-fre -rre joint inde- pendence test, xdHSIC for any d. (2) To address the com- binatorial challenges of testing for high-order factorisation, we introduce the cross-centring technique to develop the permutation-free d-order Lancaster factorisation test, xLI. (3) Building XLI, we derive a permutation-fre Streitberg interaction test, xSI. (4) We show that, for the special case d 2, 2, our high-order tests have a more concise mathemat- ical formulation than the existing xHSIC (Shekhar et al., 2023). (5) Through empirical evaluations in machine leamn ing applications, we demonstrate that the permutation- ffee
Contributions
'Department of Mathematics, Imperial College London, United Kingdom DDpaatmentt of Neurology, University Hospi- tal Würzburg, Germany DDepartment of Brain Sciences, Imperial College London, United Kingdom. Correspondence to: Mauricio Barahona < m. barahhoooa @ imperial.ac.uk>.
Proceedings of the nddd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Permutation-Pri Free High-Order Interaction Tests
tests achieve substantial computational speedup, exceeding 100-fold, over their permutation- baaedd counterparts.
d = 2
Pairwise Independence A AiP- APP A3P
Paper Structure Section 2 introduces joint independence and partial factorisation as two different generalisations of pairwise relationships, and briefly reviews kernel- -ased tests. In Section 3, we develop our permutation-Fes joint independence test for any d, XdHSIC, and our permutation- free factorisation tests for any xLI and XSI. We relate our proposed methods to existing work in Section 4, and provide numerical validations in Section 5. We conclude by discussing directions for future research in Section 6.
==3
Joint Independence
47P
Factorisation
Partial
4 A3P
Joint Independence
4
Lancaster Interaction
A
Streitberg Interaction
4
more factorisations Pr Ploout are included
Figure 1. Generalisations of pairwise dependence high-order interactions. All three measures of statistical interaction (joint independence, Lancaster, Streitberg) are the same for d 2. Joint independence captures less information for d 2 3. Only Streitberg interaction considers all factorisations for
2. Preliminaries
2.1. Statistical Interactions
Throughout this paper. we consider multivariate probability distributions of a set of d variables D {X {X.....,X*) x4) Pi..d-
If 123 can be factorised into any of P P2,, P2P13 PaP or PiP2Ps, then A 0 (see Appendix B). Lancaster (1969) rewrote and generalised this measure for d-order interactions as
denoted P xx
Pairwise Interaction The simplest statistical interaction emerges for d = 2. In this case, lack of statistical in- teraction is equivalent to pairwise independence, which is determined by the difference between the joint distribution P12 and the product of the marginals P,P2:
A P
(P P.).
(4)
where the P.; are defined implicitly by P.P...PP. Pig... and IP; PP, P..
A 2 = P12 - P2
(1)
However, for d 4, A fails to capture certain fac- torisations, e.g.. if P1231 F1eP31 then A (P12 I (1Px .31P,) which 1990). Despite this limitation, Liu et al. (2023b) proved that for d 4, A 0 only if the factorisation of Fi... includes singleton variables, highlighting its partial informativeness in capturing structural dependencies.
Hence A = X 1 x (pairwise independence).
zero in general (Streitberg.
Joint Independence One way to generalise the above notion to d > 2 variables is to focus on joint independence:
A IP 1...d
1 P.
Clearly, AfF vanishes if and only if the d variables are jointly independent, property that has been employed in Independent Component Analysis (Hyvärinen & Oja, 2000), causal discovery (Pfister et al., 2018; Laumann et al., 2023), and variational autoencoders (Lopez et al., 2018). How- ever, for 3, joint independence is uninformative as criterion for the presence of high-order interactions. For example, three of the four factorisations of Pizs (IP P23.
Streitherg Interaction To account for all factorisations of Pi..d the Streitberg interaction was introduced (Streitberg, 1990). Let π = by|bel 16 denote a partition of the set D of d variables into ||| blocks bj. and let P. == 11 Phy. The Streitberg interaction for d variables is:
A Σ (-1) Tt()((
(5)
PaPia and 12) are neither full 3-way interactions nor jointly independent like P2P3.
where the sum extends over IT(D), the set of all partitions of the set D. Crucially, A/P 0 only If PP... can be factorised in any way (Streitberg, 1990).
Lancaster Interaction An alternative generalisation in- terprets the vanishing condition of TP as factorisation of the joint distribution into disjoint components. To capture this perspective, the Lancaster interaction for variables defined as:
Figure 1 summarises these three perspectives of high-order interactions, showing how they all coincide for d = 2; begin to diverge for d= 3; and Streitberg interactions are unique in capturing all possible factorisations for 2 4. For details on derivations via lattice theory, see Liu al. (2025).
A 3 TP Pi23
P. P23
P2P113 P3IP12 + 2P, P2P3.
(3)
2

---

Permutation- Free High-Order Interaction Tests
2.2. Kernel-Based Tests
where Pit) and Pi) are the empirical embeddings in dHSIC, estimated from two disjoint sample sets: the first half S {(1) "., ..,,,")))) sssisn) n} and the second half S {(0), ",,,,,,)))) in 1 1 2m) of the i.i.d. samples.
To transform the interaction measures from Section 2.1 into non-parametric test statistics, embed the under- lying distributions into reproducing kernel Hilbert spaces (RKHS). Recall that symmetric, positive-definite ker- nel k² : X' x X > R induces an associated RKHS H with the reproducing property. For any X' E X', let (x+) = k" (x denote the canonical feature map. The kernel mean embedding of a distribution P; is then given by ut; = Ex([d [db (x*)]- This embedding satisfies Mr,
Unlike &HSIC in (Shekhar et al., 2023), which was based U-statistics, here adopt V-statistics to reduce com- putational overhead, while retaining robust estimation ac- curacy (Pfister et al., 2018). Specifically, we form the 2n x 2n kernel matrices K' with entries Kab = ka af}) for 1 S a, b S 2n where {=1) are i.i.d. samples for the i-th variable. In our construction, the n x n submatrices on the diagonal of K (with entries K where both sam- ples a, b come from either Sf from Sy) do not contribute to the estimate. Instead, we only use the n n upper-right off- -if-oaaa block (lower-left block is redundant due to sym- metry), which measures distances between samples from the two disjoint subsets. We then define XdHSIC, the nor- malised permutation free dHSIC.
Ex f(X') (f, /rr) hES for all in H', where (ss is the inner product. If k is characteristic, the map μe; is injective, ensuring that the norm of the signed measure is zero if and only if the measure itself is zero (Gret- ton et al., 2007; Muandet et al., 2017).
Using such embeddings, all high-orda order interaction measures in Section 2.1 can be naturally converted into kernel-based test statistics by constructing kernel mean embeddings of the relevant distributions, and computing their Hilbert-Schmidt norms. Adopting this methodology, the measures in Eq. (2), Eq. (4) and Eq. (5) have been developed into, respectively, the dHSIC (Pfister et al., 2018), the Lancaster interaction test (LI) (Sejdinovic al., 2013a; Liu et al., 2023b), and the Streitherg interaction test (SI) (Liu et al., 2023b).
Definition 3.1. Let KC' K
[[mm,-+1:2]] be the off- off-diagonal n n submatrix of K2, Then xdHSIC is defined as:
xdHSIC- 1R-xHISIC
with
87
xdHSIC = no
T
Okx)
d
1+
1
d
n2d
A key practical issue is the fact that all of these test statistics are degenerate under their respective null hypotheses. To obtain valid rejections, one typically employs a permuta- tion procedure to approximate the null distribution, or uses heuristic strategies such as a Gamma approximation (as in HSIC (Gretton et al., 2007) and dHSIC (Pfister et al., 2018)). Permutation-basse approaches involve p runs (typi- cally p> > 100) achieve an accurate approximation of the null distribution, and are thus computationally expensive for large sample sizes or multiple variables. This computational burden underscores the need for alternative methods that are both theoretically principled and computationally efficient.
f=1
(=1
1
(
d
1
(6ITK))2
no+1
j-1
no+1
(-1
57 - xx) 1+ )
1
k (
(0
d
nd (Oxr
d -
xdHSIc1 2
nd (
where the operator O represents the Hadamard product and 1 is the n x 1 vector of ones.
xdHSIC can be used to test d-order joint independence: Hypothesis 3.2. (Joint Independence)
Recently, Shekhar et al. (2023) introduced xHSIC, a variant of HSIC for d = 2 that avoids the need for permutations by establishing that its null distribution converges to a standard normal. Inspired by this approach, we develop family of permutation-feet high-order tests in Section 3.
Ho P PI,,..d .d IT IP;
H₂ PI,
PI,,..d d F IP.
Then we have: Ho = ApP-
x = 0.
3. Permutation-Free High-Order Tests
The construction of xdHSIC ensures it follows N/ 1) under the null (Shekhar et al., 2023), thus enabling the test with a single computation in contrast with p permutations for standard dHSIC. Importantly, its time complexity matches that of single dHSIC computation.
3.1.Joint Independence
To extend the permutation-free pairwise independence mea- sure of Shekhar et al. (2023) to the multivariate setting, we first define the unnormalised permutation-free dHSIC using the data-spliting technique:
Proposition 3.3. xdHSIC can be computed in O
To prove this, it is sufficient to show that both the numer- ator, xdHSTC, and the denominator, 87. can be computed
xdHSIC= = (pbb... - p P
3

---

Vermutation-Free High-Order Interaction Tests
in (dn'), which follows from an argument in Pfister al.
factorisations involving at least one singleton variable (Liu et al., 2023b). To formalise the test, we introduce the fol- lowing null hypothesis:
(2018) (see Appendix A).
3.2. Partial Factorisation
Hypothesis 3.6. (Lancaster Factorisation)
A limitation of the joint independence test is that, if the null hypothesis is rejected, it does not distinguish between partial factorisation, where disjoint subsets of variables indicate low-order interactions, and true fattriitoii,ttion, where the joint distribution cannot be factorised in any way, implying that all variables interact. To address this limita- tion, more sophisticated criteria have been developed (Stre- itberg, 1990), and the corresponding kerne--based tests have been introduced (Sejdinovic et al., 2013a; Liu al., 2023b). However, these methods encounter degenerate cases, which necessitate permutation- braadd resampling a step that is climinated in the permutation-fri Lancaster and Streitberg interaction tests introduced below.
Ho : V HEL
V | =2
XL H1 A M1-H1 -I
v V|TT|= =
where Hx, is the event such that the joint distribution fac- torises according to partition π containing exactly one singleton, whereas =H, indicates that such a factorisation does not hold. If any of the d subhypotheses in Ho is true, then A = 0. Hence if all d subhypotheses are rejected, a d-order Lancaster interaction is detected.
Note that even though Hy contains only d subhypotheses, which correspond to the factorisations with one singleton, all other factorisations that satisfy the Lancaster vanishing condition are actually subsumed by them (Liu et al., 2023b).
Centring: Our approach relies on matrix centring, which widely used in both kernel-based (Gretton et al., 2007: Sejdinovic et al., 2013a) and distance-based (Chakraborty & Zhang, 2019; Böttcher, 2020) high-order interaction tests to improve computational efficiency and simplify notation. For and x n kernel matrix K, its centred version is K = CKC, where C = I, 1 11",, where I is the n x n identity matrix.
We can now define the normalised permutation- statistic associated with a subhypothesis HT, of Ho.
Definition 3.7. Let H : PPi..d Pmll Pmm.i,....m- "....m-1,m-- 1m=1, P for some S m S d. Under the null of H Ha the normalised permutation-fice Lancaster interaction, xLI, is:
Equivalence: Two kernels are said to be equivalent if they induce the same semimetric on the domain (Sejdinovic et al., 2013b). In particular, the centred kernel is translation- invariant and thus is equivalent to the original kernel.
XLI = Vn.xLI xLI
with
8L.
Cross- Cross-centring: To extend the permutation-free scheme to a broader class of measures, we introduce an analogous cross-centring, which arises naturally under data splitting.
xL1 = 1
K
O
K 1
n² u
ie D/m
1 diag(1.)
Definition 3.4. (Cross-centring) Let 1,, ()) @ and 1, (9) 8 1, where @ denotes the Kronecker product, and define the matrices C, 12n 11,1, T, and Ce- 12m- 41,1,", where I2n is the 2n x 2n identity matrix. For a 2nx 2n kernel matrix K with indices 1 n associated with S and indices n + 2n associated with Sy its cross-centred version is given by K = CC,
$7 = n13
K
K
1
xL12
seD/m
where D/m denotes all d variables except XT the opera- tors o, O represent Hadamard products, and the overline (.) indicates cross-centring.
Remark: In our permutation-fis procedure we have a dif- ferent test statistic for each subhypothesis, all with fixed null distribution (standardised (0, 1)).
Lemma 3.5. The cross-centred kernel
k(I,²') = ((x)) ((00))=(00) 1 th 1=1
[dz)) 0(x') Σ E
2n
To reject the null Hy for Lancaster factorisation, we carry out d subtests (one for each of the singleton factorisations), each with complexity Odn?).
j=n+1
where 66-) is the canonical feature map, is equivalent to k(x,z') as n + 0o. See Appendix A.
Proposition 3.8. The time complexity of xL1 for all Haz across any d is (dn²).
Having established this asymptotic equivalence, we express our test statistics in terms of cross-centred kernel matrices.
The composite hypothesis Ho is rejected only if all sub- hypotheses are individually rejected, and, as this has been proven to be less conservative than the Bonferroni correc- tion (Rubenstein et al., 2016), we use this correction for our tests. When the subtests are executed sequentially rather
3.2. 1. LANCASTER INTERACTION TEST
We now derive the d-order permutation-Fis Lancaster in- teraction test. The Lancaster criterion can be used to detect
4

---

Permutation- Free High-Order Interaction Tests
than in parallel, the procedure can terminate early if any subtest fails to reject its corresponding null, thereby saving computational time. We provide a detailed description of the test procedure of XLI for d = 3 in Appendix B.
In our numerics, we test the complete factorisation hypoth- esis (3.9) via two-step process denote XIII+XII: we first use XLI to test the d subhypotheses involving TT. (bi- partitions with singleton), followed by applying XSI for the subhypotheses involving the non-singleton bipartitions, Tg. Together, these two sets cover all the (2 1) bi- partitions needed to test for the complete factorisation, i.e.. T|=-2 TL wg. Although XSI could, in principle, be used to perform all subtests in Hypothesis (3.9) directly, we found this approach to be computationally less advanta- geous. Indeed, given that Liu et al. (2023b) showed that the permutation-ban
3.2.2. STREITBERG INTERACTION TEST
As discussed, the Lancaster interaction is unable to detect factorisations that do not involve singletons. The measure that does vanish for all factorisations of the joint distribution is the Streitberg interaction in Eq. (5). To test if Pi--d can be factorised in any way, we have the following hypothesis:
SI test involves only non-singleton parti- tions, if we were to use XSI, the terms in the permutation- free xST exceed this scope, further contributing to its higher computational cost if it were applied for all subtests. As a result, using xSI alone does not always guarantee improved computational efficiency. In contrast, XLI+XSI ensures effi- ciency, by maintaining quadratic complexity with the num- ber of samples n regardless of order d. This advantage becomes increasingly significant as d grows, making our method more suitable for practical applications involving high-order data.
Hypothesis 3.9. (Complete Factorisation)
Ho: V Z
V ||| 2
H H
Y |
where H denotes the event such that the joint distribution Pi..d is factorised as P== Pb, Pbg.
Note that only partitions with exactly two blocks need to be considered in the composite tests, as all further factorisations Pi.-.- = br ...... are subsumed in bipartition subtests. The partitions with two blocks correspond the second level of the d-order partition lattice (Liu et al., 2023b).
3.3. Computational Complexity
Table provides a summary of the time complexity of the permutation-baxxe and permutation tests, highlighting the computational advantage of the tests proposed this paper. For the joint independence (dHSIC) and Lancaster interaction (I.I) tests, the permutation-fiet tests offer p- fold reduction in complexity by eliminating the need for resampling. For the complete factorisation test (SI), our strategy of combining XSI and XLI is far more efficient than the permutation- based test in two ways: (i) p-fold reduction for each each subtest; (ii) n² scaling decoupled from d. Empirical evaluation of these scalings is shown in Section 5 using synthetic data with ground truths.
The Streitberg (complete) factorisation hypothesis includes the Lancaster factorisation hypothesis plus, additionally, the Hag where Ts are partitions into two blocks with no singletons. For example, for d = 6, the Streitberg subtests must check specifically for P P123P456. whereas, in contrast, F P1P23156 is already included in the Lancaster factorisation tests (due to the singleton), and IP12TP34PP36 does not need to be checked, since it has cardinality 3 and is subsumed in 2-block factorisations P F1234 56. 11255 34 and F12
Definition 3.10. Under the null subhypothesis of Hag. where ns by (b2, the permutation- ffree Streitberg inter- action is defined as
Table 1. Time complexity of high-order tests. Here p is the num- ber of permutations, n is the number of samples, d is the number of variables and Fa is the number of partitions without singletons for d variables, a combinatorial number. The permutation-free tests developed in this paper are highlighted in blue.
Vi- xSI
XSI=
with
$3
1
XSI
1T
O
K
K"
1t,
n
a
Permutation- based
Permutation-fice free
péb,
by
dHSIC: (pdn") LI: (pd?n")
xdHSIC: Odd")
1 (1))
xII: (dP
85-
O geb2
K"
1
xSI
3
S1: Opt
XSI + +LII: O(2
pebs
Unlike the permutation- based Streitberg statistic, whose time complexity grows combinatorially with d (Liu al., 2023b), xSI grows only quadratically with
4. Related Work
Our work is inspired by Shekhar et al. (2023) who in- troduced xHSIC to make HSIC (Gretton et al., 2007) permutation-fre by approximating HSIC using the data-
Proposition 3.11. The time complexity xSI for all H across any d (dn?).
5

---

Permutation-Fiet FFee High-Order Interaction Tests
splitting technique and U-statistic (Kim & Ramdas, 2024; Shekhar et al., 2022), and normalising by the empirical stan- dard deviation, thereby achieving standard normal null distribution. However, to the best of our knowledge, these methods do not extend readily to more than two variables. In contrast, our kernel-based formulation leverages V-statistics and cross-centring to define permutation-fiee tests for Lan- caster and Streitberg interactions for d 2 variables.
4-way MVG: power
4-way MVG: time
1.0
0.8-
0.8
0.6
0.6
0.4
0.4
dHSIC RdHSIC
dHSIC
RdHSII
500
1000
0.0
10-2 time
10
sample size,
Remark: For d 2, our formulation recovers pairwise in-
time vares
time d' varies
dependence (Fig. 1). Yet by using cross-centring and V- statistics we obtain a permutation-fre statistic XHSICy := xSI = xL_,--, which is simpler, more compact, and with higher power than the original XHSIC in Shekhar et al. (2023). See derivation and numerics in Appendix C.
0.012
SdNSIc [sacena]
0.010
1.0
diiic (ecoond) Relisic (ecand)
0.03
CHSIC (secomd)
0.008
-0.8
0.006
-066
0.02
0.004
0.4
0.02
0.002
0.2
5. Experiments
0.000
250
500
750
1000
0.00
2.5
10.0
sample
order,
We evaluate our proposed permutation-free tests using range of experiments. For all cases, under the null hy- pothesis of each test, we confirm standard normality and controlled type-I errors of the corresponding xdHSIC, xL1 XSI statistics (see Appendix D, Figs. 10-11 for examples). Below focus on the statistical power and computational efficiency of these methods compared to their permutation- based counterparts. Unless noted otherwise, the significance level is set to a 0.05, and we use Gaussian kernels with bandwidth given by the median heuristic. Additional details of each experiment are available in Appendix E.
Figure 2. Joint independence test. (a)-(c) Sampling from d variable MVG (covariance matrix, inset of (a)), we compare dHSIC (with p = 100 permutations) and xdHSIC: (a) statisti- cal power as function of the number of samples, (b) trade-off between power and computational time (circle diameter propor- tional to sample size n);; (c) CPU time as increases. (d) CPU time for MVGs of the same form as in (a) but with increasing number of variables, d (with
500 samples). Note that the right scale of the CPU time in (c) and (d) is 100 times larger, thus reflecting the p-fold computational reduction over dHSIC.
a
b
1.0-
1.0
5.1. Validation on Synthetic Examples
type-1l error
type-1 error
Joint Independence Test We first demonstrate the RdHSIC statistic for permutation- joint independence testing. Figure 2 shows the application to d = 4 multi- variate Gaussian (MVG) distribution (covariance matrix in the inset of (a), with 3 0.5). The power of the xdHSIC test has similar sensitivity compared to the permutation-based dHSIC test with p = 100 permutations, but xdHSIC reaches the maximum power roughly 100 times faster (Fig. 2b-c), as expected from the p-fold reduction in complexity. The same reduction in computation time obtained when applied to d-variable MVGs of varying d (Fig, 2d).
0.6-
0.6
XdHSIC
XU
0.2
XSI+XUI
0.2
0.0
0.5
1.0
0.0
0.5
1.0
interaction strength, 8
interaction strength, p
Figure 3. Partial factorisations of MVG with d 5 variables. (a) xdHSIC fails to detect the partial factorisation with singleton, P. P2P4s regardless of interaction strength. (b) Only xS1+XLI identifies the factorisation with no singletons, P12Ip4ss MVG covariance matrices are insets in each case.
Factorisation Tests Next, we demonstrate that our permutation-fres high-order tests successfully preserve the vanishing conditions for different factorisation tests. In Figure 3a, we consider data from a MVG with d 5 vari- ables and covariance matrix such that P1 P254s (see inset). As expected, xdHSIC fails to detect this partial (singleton) factorisation, which is correctly detected by both ÆLI and XSI+XLI Figure 3b shows that, when analyse data from a MVG that is factorisable as F12Pa15 (no singletons), only XLI+XLI succeeds in detecting this interaction.
We then simulate data from a XOR gate with d 5 variables, in which the interaction exclusively of order without 2- 3-, or 4-order interactions. Fig. 4a-c shows that the power of the permutation- xSI+xLI is superior to the permutation-
SI with substantial reduction in com- putation time. For example, for n 500 samples, xSI+*LI can be computed within 0.5 seconds whereas SI requires roughly 6 minutes. In Figure 4d, we consider a series of XOR gates with increasing number of variables, d. Whereas XSI+XLI can be computed in less than 1 second for d = 10,
6

---

Permutation- High-Order Interaction Tests
SI becomes already computational infeasible for d = 6, due to the rapid combinatorial increase in the number of non- singleton partitions (Fa) (Sloane et al., 2018), which makes the computation quickly intractable as d increases (Table 1).
making xLII an attractive choice for detecting V-structures when conventional conditional independence tests become unreliable or too slow.
Furthermore, in Appendix B we thoroughly evaluate the sen- sitivity of the permutation-free method by applying it to V- structures encoded by three nonlinear forms (sinc, log, poly- nomial) using three different kernel functions (RQ. Laplace, RBF). We show that regardless of the non-linear form or kernel, our proposed method achieves similar accuracy to the permutation-based counterpart (Fig. 8).
b
5-way XOR: n effect
5-way XOR: 8 effect
1.0
1.0
0.8
Bower
0.6
Bower
0.6
0.4
0.4
0.2
0.2
0.0
0.0
a
dataset A
b
dataset B
50
100
0.00
0.25
0.50 0.75
1.0-
1.0
sample size,
interaction strength. $
time varies
d
time as d varies
type-il entor
0.8
type-1l eror
0.8
500
0.6
0.6
VS|+XUI (sacond)
0.8
RSI+RLI
400
XSI+ XSI+XLI (secand)
10
Si (second))
0.6-
SI (second)
0.4
xlp:))))")
0.6
300
Ol:O)p"))
200
0.4
10c
0.2-
0.2
XL: ((o?)
0.0
0.0
0.2
100
0.24
10
noise dimension
10
15
noise dimension
10 1.5
0.0
0.0
200
400
2
sample
order,
Figure 5, V-structure detection. Type-II errors of XLI, LI and KCI for detecting the conditional dependence as the noise dimen- sions vary in dataset (a) A and (b) B from Sejdinovic et al. (2013a).
Figure 4. 5-way XOR factorisation. Power of xSI+xLI and ST tests for an XOR gate with d = 5 variables as a function of: (a) sam- ple size, n;; (b) interaction strength, 3. Computation time of both tests for: (c) increasing sample size, n (for d = 5); (d) increasing number of variables, d (for 100). The ST test is computation- ally infcasible above ==3 5 (dashed grey line). Note the different loft/rght scales for CPU times in (c) and (d), highlighting the large reduction in computational cost achieved by xSI+xLI.
DAG Causal Structural Learning An alternative route for causal learning is provided by cassed methods (Pe- ters et al., 2014; Laumann et al., 2023), whereby each vari- able is modelled as a function of its direct causes (i.e., parent nodes in directed acyclic graph (DAG)) plus an additive noise term. If the noise terms for all regressions are jointly independent, then the proposed DAG cannot be rejected.
5.2. Applications to Causal Discovery
Following closely Simulation from Pfister et al. (2018), we simulate n samples from one DAG with d = 4 variables (see Appendix D). We then compare the accuracy of our permutation free xdHSIC against the permutation based dHSIC (using p = 200) to recover the ground truth DAG. Figure 6a shows that xdHSIC identifies the correct DAG more frequently than dHSIC for n > 500, and it reaches perfect accuracy using roughly two-thirds of the samples needed by dHSIC, all of it with a 200-fold computational speed-up. Similarly, the structural Hamming distance (SHD) of xdHSIC reaches zero faster than dHSIC (Fig. 6b).
V-Structure Detection A V-structure (or collider) is fundamental component in bonstraint-based causal discov- ery, where two variables X and Y both causally influence a third variable Z. If X 1 Y, rejecting the Lancaster fac- torisation hypothesis for = 3 is equivalent to rejecting the conditional independence X u Y|Z, thereby indicating the presence of V-structure (Sejdinovic et al., 2013a).
Here, we use two datasets (A and B) from Sejdinovic et al. (2013a), in which X a Y. and we increase the noise dimen- sion to make the causal relationship more difficult to detect. Figures 5a-b shows the comparison of our permutation-fis XLI to the permutation- based LI and kernel-based condi- tional independence KCI (Zhang et al., 2011) tests. Each experiment uses n = 500 samples and the time complexity is On²), (pn²) and Opn"", respectively. KCI rapidly fails as the noise dimension grows, whereas Lancaster-basee based approaches remain more robust. Although XLI is marginally less powerful than LI, it runs approximately 100 times faster (0.082s vs. 7.6s for p = 100 permutations), providing a trade- off between statistical power and computational efficiency,
In real-world causal discovery, one often needs to search over all possible DAGs (see Appendix D). For 3, 4, 5, and 6 nodes, there are 25, 543, 29281, and 3781503 DAGs, respec- tively. Since the ratio of DAG counts between consecutive node sets (3, 4, 5, 6} is less than 200, the computation time required for HHSTC (with p 200 permutations) to test joint independence in DAGs of 3 (or 4, 5) nodes exceeds that of xdHSIC applied to DAGS with an additional node, i.e., 4 (or 5, 6). This efficiency gain can facilitate experiments to uncover causal relationships more effectively.
7

---

Permutation-Free High-Order Interaction Tests
b
DAG detection SHD
DAG detection accuracy
Table 2. Feature selection. Only xSI+XLI identifies the ground truth variables that generate Y. Other methods either include confounding variables or miss generating variables (O).
12
dHSIC XdHISIC
accuracy
0.6-
SHD
Method
x² x²
X
X²
X
X
X
X
X
MSF.
Variance
5.01
0.4
Univariate
5.01
0.2-
dHSIC RdHSIC
RFE
Sequential
...
....
4.47
4.69
SHAP &HSIC ROHSIC RLI
4.67
500 750
1000 sample sice, n
300 400 500 600 700 800 sample size, n
10
4.67
5.01
Figure 6. DAG causal discovery. (a) Accuracies of the correct DAG being detected among all fully connected DAGS of four nodes and (b) the distributions of the SHD of the detected DAG using xdHSIC and dHSIC out of 100 experiments.
xSI+xL1
4.35
Ground truch
4.35
Figure 7 shows the percentage of 2-, 3-, 4-, and 5-way inter- actions for all sectors, where we compute interactions for 500 sets of stocks within the same sector and, as a compari- son, for 5000 sets of stocks randomly drawn from different sectors. As expected, 2- 3-, 4-, and 5-way interactions are more common within a sector than across sectors, in con- cordance with the GICS industrial taxonomy that underpins the defimition of the sectors Interestingly, both Utilities and Energy exhibit very high 3-, 4-, and 5-way within-sector interactions. This likely stems from the highly regulated na- ture of these sectors, resulting in stocks with similar returns and, consequently, high (within-sector) redundancy. Con- versely, Information Technology, Consumer Discretionary and Health Care display lower high-order interactions. This suggests that companies within these sectors may be more closely connected to firms in other sectors as they are likely influenced by external factors, indicating (cross-sector) syn- ergistic relationships. A summary heatmap of the high-order interaction profiles of all sectors shown in Figure 13 in Appendix E. These high-order interactions could be of help when considering portfolio diversification.
5.3. Application to Feature Selection
Feature selection is a crucial step in many machine learn- ing pipelines. The task is to identify the most informative features to enhance model performance and interpretability. Most techniques, however, focus solely on pairwise rela- tionships between a feature and the target variable (Póczos et al., 2012; Song al., 2012), overlooking high-order in- teractions (Fumagalli et al., 2024; Muschalik et al., 2024).
To illustrate the application to feature selection, we con- struct a dataset (n = 500) with nine features {x),... XX),
of which {Xr ,..., } are i.i.d. uniform variables interact- ing through an XOR gate to generate a target variable Y. and {X*, X²) are MVG's with covariance [[1, 0.9], [0.9, I]]- The ground truth interaction between the variables is thus (Y.X'... ""....,",?; a. {X*, ;x hence no single variable or pair of variables is predictive of Y. In Table 2, we compare the performance of feature selection methods based on high- order tests with popular methods based on univariate and pairwise measures, including the sequential method, recur- rent feature selection (RFE) (Pedregosa et al., 2011) and SHAP (Lundberg & Lee, 2017). The permutation-fret test XSI+XLI is the only method able to detect the ground truth feature set with lowest mean squared error (MSE). All other methods, including pairwise independence, joint indepen- dence and Lancaster tests, suffer from errors in selecting the correct features due to their inability to identify the complex high-order interactions between the input features and the target variable in this example.
6. Discussion
The use of fine-tuned permutation schemes has been shown to be an essential ingredient for interaction tests to produce reliable results, but this process can be both time-consuming and challenging (Rindt et al., 2021). The permutation- free methods introduced here, xIHSIC, XLI and XSI, offer a more efficient approach to handling complex data structures circumventing the need for permutations and the tuning of associated hyperparameters, thus simplifying the workflow and ensuring consistent performance.
5.4. Application to Dataset of Stock Daily Returns
To demonstrate the scalability and applicability of our permutation-free tests, consider large real-world fi- nancial dataset: daily returns of stocks in the S&P 500 from 2020 to 2024. This dataset comprises 504 variables (stocks) and 1004 samples (returns), spanning 11 sectors (Fig. 7). As is standard in the finance literature, we assume the daily returns to be i.i.d. (Ali& & Giaccotto, 1982).
Regarding when to use permutation-fiee vs. permutation- based tests, our experiments show that for small sample size, permutation-based methods perform comparably, and occasionally slightly better. This is likely due to the need for data splitting in permutation-frr methods, which limits their ability to capture the full data structure in low-lample regimes as test statistics rely on inner products between
8

---

Permutation-Fres Free High-Order Interaction Tests
1.0
d-Order
HOI Percentage
2way
3way
4way
5way
0.2
0.0
Information Tectnology
Steples
Indurtials
Financials
Health
Utiitiea
Deccmitonary
Real Escate
Materials
Communication Services
trandomisd)
Figure 7. Stock interactions in S&P500. Utilities and Energy show high 3-, 4- and 5-way interactions, likely due regulation redundancy. contrast, Information Technology, Health Care, and Consumer Discretionary exhibit lower high-order interactions, suggesting stronger cross-sector connectivity.
embeddings estimated from the separate halves. Hence, rule of thumb, for low number of samples and low order d the trade-off between computation time and sta- tistical power would not be unfavourable for permutation- based methods, e.g., for n S 50 and d S permutation- based methods perform well with just = 100 permutations. Outside this regime, however, we strongly recommend the permutation-free methods introduced here.
Acknowledgements
MB acknowledges support by EPSRC through grants EP/W024020/1, funding the project "Statistical physics of cognition", and EP/NO14529/1, funding the EPSRC Cen- tre for Mathematics of Precision Hcalthcare, and by the Nuffield Foundation, under the project "The Future of Work and Well-being: The Pissarides Review". RP acknowledges funding from the Deutsche Forschungsgemeinschatt (DFG) Project-ID 424778381- 295.
Nevertheless, several challenges and limitations remain. The data-splitting technique that underpins our permutation- free strategy relies on the assumption of i.i.d. samples, which can be restrictive when dependencies exist across sam- ples (e.g., time series, network data, spatial data). Although permutation based techniques have been developed ad- dress stationary (Chwialkowski & Gretton, 2014; Ruben- stein et al., 2016) and non-stationary time-series data (Liu et al., 2023a), the development of permutation-free meth- ods for such data remains an open challenge. Additionally, while we focus here on joint independence and factorisa- tion hypotheses, other forms of higher-order structure, such as conditional independence or more intricate composite hypotheses, might benefit from similar permutation- free schemes. Finally, although the permutation-fret free approach substantially improves the efficiency of high-order tests, there remains potential to further reduce computational over- head by incorporating scalable kernel approximation tech- niques such as random Fourier features and the Nyström method (Zhang et al., 2018). Investigating these directions could further broaden the applicability of kernel-based high- order interaction tests in large-scale machine learning and statistical settings.
The authors thank Felix Laumman, Jianxiong Sun, Shub- hanshu Shekhar and Zhuoyu Li for valuable discussions.
Impact Statement
The societal impact of our work stems from its ability to address the critical need for efficient and scalable meth- ods to detect high-order interactions in multivariate data, By eliminating the reliance on computationally intensive permutation schemes, our permutation free tests enable re- searchers to uncover complex relationships beyond pairwise interactions, which are crucial in domains like causal dis- covery and feature selection. These advancements facilitate more effective data analysis in fields such as healthcare, eco- nomics, and environmental science, where understanding intricate interactions can lead to breakthroughs in diagnosis, decision-maken making, and policy formulation.
References
Ali, M. M. and Giaccotto, C. The identical distribution hypothesis for stock market prices-location-andd scale- shift alternatives. Journal of the American Statistical Association, 77(377) 1982.
Code Availability Code to implement the permutation- free tests in this paper available at https: //github. cow/barahona-aaeeee tahhon-research-groupp ree-HOr -o..-HHI.glll. .git.
Arnaudon, A., Peach, R. L., Petri, G., and Expert, P. Con- necting hodge and sakaguchi kuramoto through a mathe- matical framework for coupled oscillators on simplicial complexes. Communications Physics, 5(1):211, 2022.
-akaguchi-kuramotto

---

Permutation-Fre High-Order Interaction Tests
Battiston, F., Cencetti, G., lacopini, I,, Latora, V., Lucas, M., Patania, A., Young. J.-G., and Petri, G. Networks beyond pairwise interactions: structure and dynamics, Physics Reports, 874:1-92, 2020.
Laumann, F,, Von Kügelgen, J., Park, J., Schölkopf, B., and Barahona, M. Kernel-based independence tests for causal structure learning on functional data. Entropy, 25(12): 1597, 2023.
Battiston, F., Amico, E., Barrat, A., Bianconi, G., Ferraz Arruda, G., Franceschiillo,, B., lacopini, I., Kéfi, S., La- tora, V., Moreno, Y., et al. The physics of higher-order interactions in complex systems. Nature Physics, 17(10): 1093-1098, 2021.
Liu, Z., Peach, R. L.,, Laumann, F., Vallejo Mengod, S., and Barahona, M. Kernel-basee joint independence tests for multivariate stationary and non-stationary time series, Royal Society Open Science, 10(11))23087, 2023a.
Liu, Z., Peach, R. L., Mediano, P. A., and Barahona, M. Interaction measures, partition lattices and kernel tests for high-order interactions. Advances in Neural Processing Systems, 36, 2023b.
Benson, A. R., Gleich, D. F., and Leskovee, Higher-order organization of complex networks. Science, 353(6)
2016.
Bick, C., Gross, E., Harrington, H. A., and Schaub, M. T. What are higher-order networks? SIAM Review, 65(3): 106-731,220 2023.
Liu, Z., Barahona, M., and Peach, R. L. Information- theoretic measures on lattices for high-order interactions. AISTATS - Proceedings of The 28th International Con- ference on Artificial Intelligence and Statistics, PMLR, 2025.
Binkkosski, M., Sutherland. D. J., Arbel, M., and Gret- ton, A. Demystifying mmd gans. arXiv preprint arXiv: (ath.1801.1111,,11 2018.
Long, M., Cao. Y., Wang. J., and Jordan, M. Learning transferable features with deep adaptation networks. In International conference on machine learning, pp. 97- 105. PMLR, 2015.
Böttcher, B. Dependence and dependence structures: esti- mation and visualization using the unifying concept of distance multivariance. Open Statistics, 1(1):1-48, 2020.
Lopez, R., Regier, J., Jordan, M. I., and Yosef, N. Infor- mation constraints on auto-encoding variational bayes. Advances neural information processing systems, 31, 2018.
Chakraborty, S. and Zhang, X. Distance metrics for measur- ing joint dependence with application to causal inference. Journal of the American Statistical Association, 2019.
Chwialkowski, K. and Gretton, A. kernel independence test for random processes, In International Conference Machine Learning, pp. 1422-1430 1430. PMLR, 2014.
Luff, C. E.., Peach, R., Mallas, E.-J., Rhodes, E., Laumann, F., Boyden, E. S., Sharp, D. J., Barahona, M., and Gross- man, N. The neuron mixer and its impact on human brain dynamics. Cell Reports, 43(6), 2024.
Fumagalli, F., Muschalik, M., Hüllermeier, E., Hammer, B., and Herbinger, J. Unifying feature-based explanations with functional anova and cooperative game theory. arXiv preprint arXiv:2412.17155 2024.
Lundberg, S. M. and Lee, S.-L unified approach to inter- preting model predictions. Advances in Neural Informa- tion Processing Systems, 30, 2017.
Muandet, K., Fukumizu, K., Sriperumbudur, B., and Schölkopf, B. Kernel mean embedding of distributions: A review and beyond. Found. Trends Mach. Learn.. 10(1 (011-))11-44, June 2017. ISSN 1935-8237. doi: 10. 10.1561/22000000660 URL httpt:t////// /ttp:///oooooo..//////// 1561/2200000060.
Gretton, A., Fukumizu, K., Teo, C., Song, L., Schölkopf, B., and Smola, A. A kernel statistical test of independence. Advances in neural information processing systems, 20, 2007.
Grilli, J., Barabás, G., Michalska-Smiin M. J., and Allesina, S. Higher-order interactions stabilize dynamics in com- petitive network models. Nature, 548(76666:200-13, 2017.
Muschalik, M., Baniecki, H., Fumagalli, F., Kolpaczki, P., Hammer, B., and Hüllermeier, E. shapiq: Shapley inter- actions for machine learning. Advances in Neural Infor- mation Processing Systems, 37:130324-130357 2024.
Hyvärinen, A. and Oja, E. Independent component analysis: algorithms and applications. Neural networks, 13(4-5): 411-430, 2000.
Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cour- napeau, D., Brucher, M., Perrot, M., and Duchesnay, E. Scikit-learn Machine learning in Python. Journal of Machine Learning Research, 12:2825-2830, 2011.
Kim, I. and Ramdas, A. Dimension Dimension-agnostic inference using cross u-statistics. Bernoulli, 30(1):683-711 2024.
Lancaster, H. The Chi- Chi-Squared Distribution. Wiley, 1969.
10

---

Permutation-Fret High-Order Interaction Tests
Peters, J., Mooij, J. M., Janzing, D., and Schölkopf, B. Causal discovery with continuous additive noise models. J. Mach. Learn. Res., 15(1):2009-2553, January 2014. ISSN 1532-4435.
Sejdinovic, D., Gretton, A., and Bergsma, W. A kernel test for three-variable interactions. Advances in Neural Information Processing Systems, 26, 2013a.
Sejdinovic, D., Sriperumbudur, B., Gretton, A., and Fuku- mizu, K. Equivalence of distance-based and rkhs-based
Pfister, N., Bühlmann, P., Schölkopf, B., and Peters, J. Kernel-based tests for joint independence. Journal of the Royal Statistical Society: Series B (Statistical Methodol- ogy), 80(1):5-31, 2018.
statistics in hypothesis testing. The annals of stafistics, pp. 2263-2291, 2013b.
Shekhar, S., Kim, I., and Ramdas, A. permutation-free kernel two-sample test. Advances in Neural Information Processing Systems, 35:18168-18118 18180, 2022.
Póczos, B., Ghahramani, Z., and Schneider, J. Copula- based kernel dependency measures. In Proceedings of the 29th International Coference on International Conference on Machine Learning, ICML' 12, pp. 1635-1642, Madison,
Shekhar, S., Kim, I., and Ramdas, A. A permutation-fre
kernel independence test. Journal of Machine Learning Research, 24(369):1-63 2023.
WI, USA, 2012. Omnipress. ISBN 97814 1781450312851.
Rindt, D., Sejdinovic, D., and Steinsaltz, Consistency of permutation tests of independence using distance covari- ance, hsic and dhsic. Stat, 10(1):e364, 2021.
Sloane, N. J. et al. The on-line encyclopedia of integer sequences, sequence a000296. Published electronically at https://veis. org, 2018.
Rosas, F. E., Mediano, P. A., Rassouli, B., and Barrett, A. B. An operationna information decomposition via synergis- tic disclosure. Journal of Physics A: Mathematical and Theoretical, 53(48):485001, 2020.
Song, L., Smola, A., Gretton, A., Bedo, J., and Borgwardt, K. Feature selection via dependence maximization. Jour- nal of Machine Learning Research, 13(5), 2012.
Streitberg, B. Lancaster interactions revisited. The Annals of Statistics, pp. 1878-1885, 1990.
Rosas, F. E., Mediano, A., Luppi, A. Varley, T. F., Lizier, J. T., Stramaglia, S., Jensen, H. J., and Marinazzo, D. Disentangling high-order mechanisms and high-order behaviours in complex systems. Nature Physics, 18(5): 476-477, 2022.
Thibeault, V., Allard, A., and Desrosiers, The low-rank hypothesis of complex systems. Nature Physics, 20(2): 294-302, 2024.
Rosas, F.,, Gutknecht, A., Mediano, P. A., and Gastpar, M. Characterising high-order interdependence via entropic conjugation. arXiv preprint arXiv:2410. 10485, 2024.
Young, J.-G., Petri, G., and Peixoto, T. Hypergraph recon- struction from network data. Communications Physics, (1):1-11, 2021.
Zhang, K., Peters, J., Janzing, D., and Schölkopf, B. Kernel-based conditional independence test and appli- cation in causal discovery. In Proceedings of the Twenly- Seventh Conference on Uncertainty in Artificial Intelli- gence, UAI'1, pp. 804-813, Arlington, Virginia, USA, 2011 1. AUAI Press. ISBN 9780974903972.
Rubenstein, P. K., Chwialkowski, K. P., and Gretton, A. A kernel test for three-variable interactions with random processes. arXiv preprint 1rXiv:1603.00900 2016.
Santoro, A., Battiston, F., Petri, G., and Amico, E. Higher- order organization of multivariate time series. Nature Physies, 19(2):221-229 2023.
Zhang, Q., Filippi, S., Gretton, A., and Sejdinovic, D. Large- scale kernel methods for independence testing. Statistics and Computing, 28:1 130, 2018. 8.
Santoro, A., Battiston, F., Lucas, M.. Petri, G., and Amico, E. Higher-order connectomics of human brain function reveals local topological signatures of task decoding, in- dividual identification, and behavior. Nature Communi- cations, 15214, 2024.
Schaub, M. T., Benson, A. R., Horn, P., Lippner, G., and Jadbabaie, A. Random walks simplicial complexes and the normalized hodge I-laplacian. SIAM Review, 62 (2):353-39)) 3913 2020.
Schneidman, E., Still, S., Berry, M. J., Bialek, W., et al. Net- work information and connected correlations. Physical review letters, 91(23):238701, 2003.

---

Permutation-Fret High-Order Interaction Tests
Appendices
A. Proofs for Section 3
A.1. Proof of Proposition 3.3
First observe that each term in the numerator can be computed in (dn²) since the expression only contains Hadamard product and vector matrix multiplication Pfister (2018). Similarly for the denominator, the variance estimation introduces significant computational overhead beyond simple vector compared xdHSIC.
A.2. Proof of Lemma 3.5
x The conventional centring of matrix K denoted by K CKC, where C' -----4 111 and 1 is an vector of ones, has been shown that the respective centred kernel function & is an equivalent kernel to k. Indeed, as shown in Sejdinovic et al. (2013b), a kernel y equivalent to ls it can expressed as:
by (2,2))= (k(,,)) f,k(,2)) f)
for some function f. Now for conventional centring k. ppx which can be approximated by Tix 1/ Ei di (xi) and thus is an equivalent kernel.
Here, the cross-centred Te is also equivalent as Ei ø (x;) and E & (x;) can be seen as approximation of pex-
A.3. Proof of Proposition 3.8
Here, consider the Hadamard product of all the kernels in D except K* individual kernel, then the XLI can be seen special xHSIC and therefore they share the same time complexity.
A.4. Proof of Proposition 3.11
Here, consider the two Hadamard products two kernels themselves, then the xSI can be seen as special xHSIC and therefore they share the same time complexity.
B. Extended analysis of the d = Lancaster interaction test
B.1. Lancaster interaction hypothesis and test
TP in Equation (3) is the extension as d increases from to 3 in the direction of the factorisation route as illustrated in Figure 1. Now instead of testing for joint independence through the full factorisation I2123 P222,33, we would like understand if the joint distribution Pizs is partially factorisable summarised in the composite hypothesis below:
Hypothesis B.1. (Order-3 Factorisation)
Ho (P123 P12Pss v (P123 = PPsPP) v (P123 PP23P1)
H1 : (2123 P32P3s) 123 # P2) A (P123 123 P13P1)
If any of the subhypotheses in the null hypothesis Ho is true, then A 0. For example, if P123 PiglPs, this implies F13 = PyPy and P223 = F2Ps, and, as a result, we have:
P123 Pi2P3
A Pizz - PI2IPs - P13I22 P23P1 2P
And, similarly, Pizz P1PP2 P 0 and Piz3 PzsP, APP O. Hence Hy A 0. Note that
though Ho only considers the three partial factorisations, the factorisation P. P2Ip' is actually subsumed by them, i.e. P123
P, P2I33 implies that P123 P12P3 and 1P IPa and P123 P23IP1. In other words, P123 P2P3 = A P = 0 also.
By the contrapositive, this means: P = Ho = H1. Hence, when a d 3 interaction is detected via the alternative hypothesis Hi, all three subhypotheses in Hp, as well as the complete factorisation Pi23 P22,, need be rejected.
12

---

Permutation-Fre High-Order Interaction Tests
B.2. Construction of the permutation-free normalised test statistic
By computing the RKHS embedding and applying the permutation-free scheme, the unnormalised test statistics of three- variable Lancaster interaction be immediately written out in terms of the cross-centred kernels:
2x Σ
xLI
K o K
K
1ij
Note that key distinction between the permutation--asee methods and the permutation-free approaches discussed here lies in how multiple tests are conducted. In permutation-based methods, the test statistic remains constant, while the permutations adapt to changes in the subhypotheses. contrast, with the permutation-Fi approach, where the null distributions are all standardised as N (O, 1), the test statistics must instead vary.
Without loss of generality, under the subhypothesis P123 P121P3, xLI can be simplified into:
2n
xL [123]]
K
R
K
j=n+1
1j
The normalised test statistic is
Vn-xLI
[
[LI[223]=
with
8 (1q111)
2n
$ [1219]]
Σ
Σ
K
K
K
jan+1
.ij
Σ
Σ
2n
K o K
R
xIIII
2xLI3
Σ
2n
K o K
K
1-1
j=m+|
ni g=nl
Σ
Σ
K
K
K
-XLI3
j=n+1
Similarly, for the other two subhypotheses, we just move the double cross-centring to the pairwise variables. The whole testing procedure summarised in Algorithm
Algorithm 1 Lancaster Interaction Test at 3
Input: Centred kernel matrices K R2. K computed from observational data for each subhypothesis (total three subhypotheses) do Compute the test statistic XLI
xLL > (1 a) -qaantile of 10,, then
Reject H; and proceed to the next subhypothesis
else
Terminate and do not reject composite hypothesis end if end for if all three subhypotheses rejected then
Reject the composite hypothesis
else
Do not reject the composite hypothesis end
The composite hypothesis Ho is rejected only if all subhypotheses are individually rejected. As this has been proven to be less conservative than the Bonferroni correction (Rubenstein et al., 2016), we use this correction for the factorisation tests for any d. When the subtests are executed sequentially rather than in parallel, the procedure can terminate early if any subtest fails to reject its corresponding subhypothesis, thereby saving computational time.
13

---

Permutation-Fres Free High-Order Interaction Tests
B.3. Application to additional examples: nonlinear form and different kernels
In addition the examples from Sejdinovic et al. (2013a), we also apply our method to three more examples with diverse nonlinear forms. In all cases, our proposed permutation-fret method achieves similar accuracy compared with the permutation-basee counterpart. Moreover the results are robust with different kernel functions, Rational-Quadern (RQ) kernel, Laplace Kernel and Radial Basis Function (RBF) kernel.
C
sinc
log
polynomial
1.0
0.8
0.8
0.8
type-ll errer
0.6
0.6
0.6
30 (R8F))
S1 (AB)
XJ (MBF)
0.2
231 O.apaac) KL (12)
0.2
21 (lapick)
0.2
X3 Das
XU 010)
XU (1Q)
10
5.0
10.0
noise dimension
noise dimension
noise diimension
Figure Structural causal models with diverse nonlinear forms and kernel functions. We constructed three V-structures using (a) sinc (b) log, and (c) polynomial nonlinear forms (sce inset).
C. Definition of the xHSIC pairwise test using V-statistics
We formulate the corresponding V-statistic version of the permutation-free pairwise independence test statistic, XHSICv using the cross-centring technique:
VI XHSICV
RHSICy
SV
1
2n
xHSICv
Σ
K
K
n12 =1 1
1ij
1 22
2n
Σ
K
o
K
- xHSICy
s=1
f=++1
ij
This is simpler formulation compared with the original XHSIC in Shekhar et al. (2023).
Figure 9 shows that not only is xHSICy theoretically grounded, like xHSIC, but it also has higher power compared with both XHSIC and HSIC.
XHSIC,
2-way power
1.0
0.4
Indep. L-disin
Indepp. I-distn y=0.05
0.3
N(D.1)
0.8
0.8
density
0.6
power
0.6
0.5
0.4
0.4
HSIC
0.3
0.2
RHSIC
XHSIC,
0.1
0.0
10
-5
5 10
50
100
150
200 sample
100
Figure 9. Experiments of pairwise independence tests: xHSICy follows a standard normal distribution with two independent (- distributions; controls the type-T error with the data from t-distrihutions, and outperforms both HSTC (Gretton et al., 2007) and XHSIC (Shekhar et al., 2023) with dependent MVG from the example in Equation (13) from Shekhar et al. (2023) with parameters b = 2, ndim 10, and 0.5.
14

---

Permutation- Free High-Order Interaction Tests
D. Further numerical checks of the test statistics
RdHSIC
XLI
X5I
Indep. I-distn
0.4
Indep. MVG
0.4
Indeo.
MVG
0.3
N(0.1)
0.3
N(0.1)
0.3-
N(0.1)
density
denstty
0.2
0.2
0.2
0.1
0.1
0.1
0.0
0.0
0.0
Figure 10. Null distributions of permutation-free test statistics. All three test statistics introduced in this paper follow a standard normal distribution, (0, 1) under their respective nulls.
XdHSIC
XSI
1.0
1.0
1.0
Indep. t-distn
Indep. MVG
Indep MVG
0.8
y=0.05
0.8
0=0.05
0.8
y=0.05
0.6
0.6
0.6
0.5
0.44
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.1
0.0
0.0
50
100
150
50
100
150
100
150
sample sice. n
semple size,
sample size,
Figure 11. Controlled type-I error of permutation-Fee test statistics. All three test statistics introduced in this paper have a valid level of type-I error for 0.05, indicated by the dashed line.
E. Further details of datasets used in Section 5
The code to implement the permutation free tests introduced in this paper available at https: com/ barahona- rethhonnrreeaarr-grrooooooo/pppp PwwpccuupppeeePePeeeeeee ree- HOT git.
E.1. XOR example from Liu et al. (2023b)
We generate n samples of V,W,X, X, Y, Z 4(0,4), Z = +V mod 4 and Z ~ ((0,4)) (where the samples [i n] act as noise). We then gradually increase interaction proportion, 0 s i/n
E.2. V-structure dataset A and B from Sejdinovic et al. (2013a)
Dataset A: For (X,Y, Z) R*XR*xR* X,Y N (0, W Exp ( Zi sign (X₁Y)) W, and Zzp N (0, Ip-1).
/2
Dataset B: For {X, Z) IR R x R.
X,Y (0, Zzn N (0, 1)11,, and
XP+C,
w.p. 1/3,
71
Y2+
w.p. 1/3
X,Y++,
w.p. 1/3
where ~~ 00,00.0 and W.P. stands for with probability.
In both datasets, the 3-way interaction becomes increasingly difficult to detect as the noise dimension p increases.
15

---

Permutation- Free High-Order Interaction Tests
E.3. Causal discovery example from Pfister al. (2018)
For an additive noise model over random variables X ,...,,
X
X
E PaPrret
falk (x*) + No, j (I.....dd
(6)
with corresponding DAG G. The noise variables N normally distributed and are jointly independent with standard deviation sampled uniformly between /2 and 2. Nodes without parents follow Gaussian distribution with standard deviation sampled uniformly between 5 and 10. The functions f5t are sampled from a Gaussian process with Gaussian kernel and bandwidth one.
DAG
DAG
DAG 7
DAG
DAG 9
DAG 10
DAG 11
DAG 12
DAG
DAG
17
DAG 19
DAG 20
DAG 21
DAG
DAG 23
DAG
0
Figure 12, Fully connected DAG of nodes. We use DAG as the ground truth causal structure to simulate the data.
In this example we simulate the data from DAG and compare all the possible fully connected DAGs shown in Figure 12. The procedure the causal discovery is outlined as the following: (i) Run the generalised additive noise model for each DAG and get the residuals. (ii) Check the joint independence between the d residuals. (iii) Report the DAG with the largest p-value.
16

---

Permutation-Pre Free High-Order Interaction Tests
E.4. Real-World Financial Data
The CPU time for computing high-order interaction percentages took hours without any parallelisation on a 2015 iMac with 4 GHz Quad-Core Intel Core i7 processor and 32 GB 1867 MHz DDR3 memory. We estimate the time taken for the permutation-basee methods would have been least 2 weeks. This highlights the scalability of our proposed method.
Information Technology
Consumer Staples
Industrials
0.8
Financials
Health
Utilities
Consumer Discretionary
Real Estate
Materials
Energy
Communication Services
Inter-sector (randomised)
2way 3way 4way Sway order
Figure 13. Percentages of high-order interactions between stocks within sectors in S&P500. Heatmap of the percentages shown in Figure summarising the differences across sectors in their high-order interaction structure, also relative to randomised sets (last row).
17


Reviews:
[{"review_id": "jOsEW73XRh", "reviewer": "Reviewer_s8B8", "summary": "This paper introduces permutation-free high-order interaction tests for joint independence and partial factorization of d variables. Traditional kernel-based hypothesis tests, such as HSIC and its extensions (e.g., dHSIC), rely on computationally expensive permutation-based null approximations. The authors propose a new family of tests that eliminate permutation-based computations by leveraging\nV-statistics and a novel cross-centring technique to yield test statistics with a standard normal limiting distribution under the null.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their reviews. Below we refer to the Figure S1, S2 & S3 in the link: https://imgur.com/a/uju1pHq.\n\n**Re: SCMs with diverse nonlinear forms and alternative kernels:** For the example in Fig. 5, we have added three new analyses where the V-structures are encoded by sinc, log and polynomial nonlinear forms. Furthermore, we have also included three kernels for each new nonlinear example: Rational-Quadratic kernel, Radial Basis Function kernel and Laplace kernel. Regardless of the non-linear form or kernel, our proposed method achieves similar accuracy compared with the permutation-based counterpart. See Figure S1 (a)(b)(c).\n\nSimilarly, for the DAG causal discovery application (Fig. 6), we have added a new nonlinear example that combines all three nonlinear forms simultaneously (sinc, log and polynomial) and compares three kernels. We find that our method is robust to these choices and achieves similar accuracy to the permutation-based counterpart. See Figure S1 (d).\n\n**Re: real dataset and large sample size:** We have added a large real financial dataset to demonstrate the scalability and broad applicability of our proposed method. Specifically, we compute high-order interactions between stocks in the S&P 500 from 2020 to 2024 using their iid. daily returns (standard assumption in finance [Ali & Giaccatto (1982)]). This dataset comprises 504 variables (stocks) and 1004 samples (returns), spanning 11 sectors. We compute 2-, 3-, 4-, and 5-way interactions for 500 sets of stocks within the same sector and 5000 sets of stocks randomly drawn from different sectors. The percentages of detected high-order interactions are reported in the accompanying bar plot (Figure S3). As expected, all 2-, 3-, 4-, and 5-way interactions are significantly more common within sectors than across sectors. For nearly all sectors, 2-way interactions dominate, which aligns with the GICS industrial taxonomy upon which these sectors are based upon. Interestingly, we observe that Utilities and Energy exhibit exceptionally high 3-, 4-, and 5-way interactions. This likely stems from the highly regulated nature of these sectors, resulting in stocks with similar returns and, consequently, high (within-sector) redundancy. Conversely, Information Technology, Consumer Discretionary and Health Care display lower high-order interactions. This suggests that companies within these sectors may be more closely connected to firms in other sectors as they are likely to be influenced by diverse external factors, indicating (across-sector) synergistic relationships rather than redundancy. These high-order interactions may help investors build a more diverse portfolio. These findings suggest that our method is capable of providing valuable insights into complex, real-world datasets. See Figure S3 (d) using the link.\n\nImportantly, these experiments took ~8 hours without any parallelisation on a 2015 iMac with 4 GHz Quad-Core Intel Core i7 processor and 32 GB 1867 MHz DDR3 memory. We estimate the time taken for the permutation-based methods would have been at least 2 weeks. This highlights the scalability of our proposed method.\n\nAli and Giaccotto. Journal of the American Statistical Association 77.377 (1982): 19-28.\n\n**Re: ablation study of cross-centring:** Thank you for this suggestion. We constructed an order-4 jointly independent MVG and show that only cross-centring makes the test statistic follow a standard normal distribution. Using traditional centring, or not using any centring techniques at all, results in undesired null distributions. See Figure S2.", "claims_and_evidence": "Yes.\n\nHowever, to better demonstrate the effectiveness of the proposed method, it would be beneficial if the authors could:\n\t1.\tEvaluate datasets generated by nonlinear SCMs with diverse nonlinear forms, such as polynomial, sinc, and log functions.\n\t2.\tExplore alternative kernel choices beyond the Gaussian kernel.\n\t3.\tPresent performance results for larger sample sizes.", "methods_and_evaluation_criteria": "Yes. The used datasets include (1) Synthetic datasets (Multivariate Gaussian, XOR gates), (2) Feature selection datasets (high-order interactions in multivariate data). The evaluation metrics includes Type-I error control (Figure 9), Statistical power of permutation-free vs. permutation-based tests (Figure 2, Figure 4), and Computational efficiency comparison (Table 1).\n\nIt misses  (1) SCMs with different nonlinear forms, (2) real-world datasets (e.g., finance, genomics) used for validation, and (3) ablation studies on cross-centering’s impact on variance estimation.", "theoretical_claims": "I did not check the proofs.", "experimental_designs_or_analyses": "Strengths:\nComputational efficiency evaluations provide practical insights into scaling.\nComparisons with existing methods (HSIC, dHSIC, SHAP) demonstrate the framework’s advantages.\n\nWeaknesses:\nLack of sensitivity analysis on kernel choice.\nSCMs with different nonlinear forms are not considered.", "relation_to_broader_scientific_literature": "This paper extends kernel-based independence test and builds on high-order interaction testing.", "essential_references_not_discussed": "It seems most references have been discussed.", "other_strengths_and_weaknesses": "Strengths:\n\t•\tThis paper proposes a novel permutation-free high-order dependency test.\n\t•\tSignificant computational speedup (100x faster than permutation-based methods).\n\t•\tStrong theoretical grounding and empirical validation.\n\nWeaknesses:\n\t•\tIn the experiments, the authors did not consider nonlinear SCMs with different nonlinear forms.", "other_comments_or_suggestions": "No further comments", "questions_for_authors": "1. Consider datasets generated by nonlinear SCMs with different nonlinear forms, such as polynomial, sinc, and log.\n\n2. Consider other kernels besides the Gaussian kernel\n\n3. Show the performance when the sample size is large", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "nGsMlxbk4G", "reviewer": "Reviewer_vR6R", "summary": "This paper proposes a set of permutation-free high-order interaction tests, addressing the computational inefficiency of existing permutation-based kernel hypothesis tests. It uses cross-centering techniques to derive test statistics that follow a standard normal distribution under the null hypothesis. The authors provide theoretical justification, empirical validation on synthetic data, and comparisons with existing methods. The results demonstrate significant computational advantages over traditional approaches, making high-order interaction testing feasible for larger datasets.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their positive review. Below we refer to Figure S1, S2 & S3 in the link: https://imgur.com/a/uju1pHq.\n\n**Re: iid assumption:** As the reviewers correctly pointed out, our methods depend on the data being iid and indeed in many application areas, the nature of the data is often, for example, temporal. If many realisations of each time series are present (and the realisations are iid) one can still use  permutation-free tests acting on the repeated realisations. Concretely, this means that one computes $k(x_i, x_j)$ where $x_i$ and $x_j$ are the time series vectors. Indeed, such an approach would even enable the detection of interactions between non-stationary time-series when there are multiple iid realisations. However, in situations where only one realisation of the time series is present, then one would have to rely on transformations to make the data points iid (as in stock returns in the example below), or one must check the mixing conditions in Chwialkowski et al. (2016) to develop the permutation-free tests for time series data.\n\n**Re: real and large dataset:** We have added a large real financial dataset (504 variables and 1004 samples) to demonstrate the scalability and applicability of our proposed method. In summary, we show that in highly-regulated sectors such as Energy and Utilities, there are consistently high 3-, 4- and 5-way interactions due to redundancy. Moreover we observe that the companies in Information Technology, Consumer Discretionary and Health Care are prone to external factors therefore have low intra-sector interactions. Importantly, these experiments took ~8 hours without any parallelisation on a 2015 iMac. We estimate the time taken for the permutation-based methods to be at least 2 weeks. See details in response to Reviewer s8B8.\n\n\n**Re: discuss when to use permutation-based methods and why $n=500$ for causal discovery example:** In our experiments, we observe that, when the sample size is low, the permutation-based methods achieve very similar performance, and in some cases slightly more powerful, than permutation-free methods (as noted by the reviewers for the causal discovery example in Fig. 6). We believe this is because in the low sample regime permutation-free methods are not able to capture the full structure of the data due to the data splitting—specifically, the test statistics are computed as the inner product between empirical embeddings estimated from two separate sample halves. Therefore, as a rule of thumb, if both the data available and the order of interest are relatively small, i.e., both $n$ and $d$ are small, then we could consider using the permutation-based methods. In such a case, the trade-off between computation time and statistical power would not be unfavourable. More concretely, when $n\\leq 50$ and $d\\leq 5$ the permutation-based method could be approximated well enough with just 100 permutations making it both computationally feasible and statistically powerful. However, outside of this regime we would strongly recommend the permutation-free method introduced in this paper.  \n\n \n**Re: Figure 3:** The reason that both xdHSIC and xLI fail follows from their vanishing conditions, i.e., xdHSIC becomes zero in the presence of joint independence and xLI becomes zero when the factorisation contains at least one singleton. The partial factorisation in Figure 3b is constructed precisely such that the ground truth is $P_{12345}=P_{12}P_{345}$ so it does not fall into these two categories: (i) it is neither jointly independent nor (ii) contains singletons. Therefore, both measures are unable to vanish in the presence of this partial factorisation hence their type-II errors increase as the interaction strength increases. In contrast, the Streitberg test stays fixed near alpha=0.05 with a controlled type-I error. We apologise for this lack of clarity and will make this clearer in the revised manuscript.", "claims_and_evidence": "I did not find any problematic claims.", "methods_and_evaluation_criteria": "The evaluation criteria are appropriate, with comparisons to state-of-the-art permutation-based methods (e.g., dHSIC, LI, SI) and applications to real-world tasks (e.g., causal discovery, feature selection).", "theoretical_claims": "The proofs are sound, and the theoretical results are consistent with the empirical observations.", "experimental_designs_or_analyses": "The experiments are comprehensive and support the claims made in the paper. \n\nIn Section 5.4, why the proposed method cannot work well when $n$ is less than 500? The authors should elaborate on it.", "relation_to_broader_scientific_literature": "The proposed kernel-based test can be applied in causal discovery and feature selection, bridging the gap between theoretical developments and practical applications.", "essential_references_not_discussed": "None.", "other_strengths_and_weaknesses": "Strengths:\n\n- The paper is well-written, with clear explanations of the theoretical and empirical contributions.\n- The authors conduct comprehensive synthetic experiments demonstrating the effectiveness of the approach, including comparisons with existing kernel-based tests.\n- The proposed methods are demonstrated in real-world applications, such as causal discovery and feature selection.\n\n\nWeaknesses:\n\n- The proposed methods assume i.i.d. samples, which may not hold in some real-world scenarios (e.g., time series or network data). The authors also acknowledge this limitation.\n\n- While the experiments are comprehensive, additional real-world datasets and larger simulated datasets could further validate the practical utility of the proposed methods.\n\n- The proposed permutation-free method is claimed to be an improvement over permutation-based tests, but there is limited discussion on potential drawbacks or cases where permutation-based methods might still be preferable. It might be better to discuss cases when we opt to use permutation-free or permutation-based methods.", "other_comments_or_suggestions": "Please see Other Strengths And Weaknesses.", "questions_for_authors": "In Figure 3, why $xdHSIC$ and $\\bar{x}SI$ fail to detect the partial factorization, with different interaction strengths? Please explain and give a more intuitive analysis of it.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lAU9ONsGjv", "reviewer": "Reviewer_eBnj", "summary": "The paper studies testing for independence across many variables. The Streitberg interaction is a way to test for any factorization of a joint distribution. While it had been kernelized before, now it is kernelized and centered with sample splitting, in a way that avoids permutations.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewers for their positive review and recognising our effort of making the technical literature accessible. We have corrected the typo in our revision and address your concerns below.\n\n**Re: formal statements of normality:**  The normality results follow directly  from our construction and previous results, hence we had left them out for space considerations and to focus on the new results on computational efficiency of the permutation-free tests used for high-order interactions. But we agree it will be better to have these statements formalised briefly. \n\nProposition 1 (Null normality of $\\overline{\\mathrm{x}}\\mathrm{dHSIC}$)\nLet $h_I$ be the core function of the V-statistic for $\\overline{\\mathrm{x}}\\mathrm{dHSIC}$. We assume that $E(h_I^2)<\\infty$, i.e. the second moment of the core function is finite. Then under the null hypothesis in Hypothesis 3.2 (joint independence), $\\overline{\\mathrm{x}}\\mathrm{dHSIC}\\sim \\mathcal{N}(0,1)$ \n\nProposition 2 (Null normality of $\\overline{\\mathrm{x}}\\mathrm{LI}$)\nLet $h_L$ be the core function of the V-statistic for $\\overline{\\mathrm{x}}\\mathrm{LI}$. We assume that $E(h_L^2)<\\infty$, i.e. the second moment of the core function is finite. Then under the null hypothesis in Hypothesis 3.6 (Lancaster Factorisation), $\\overline{\\mathrm{x}}\\mathrm{LI}\\sim \\mathcal{N}(0,1)$ \n\nProposition 3 (Null normality of $\\overline{\\mathrm{x}}\\mathrm{SI}$)\nLet $h_S$ be the core function of the V-statistic for $\\overline{\\mathrm{x}}\\mathrm{SI}$. We assume that $E(h_S^2)<\\infty$, i.e. the second moment of the core function is finite. Then under the null hypothesis in Hypothesis 3.9 (Complete Factorisation), $\\overline{\\mathrm{x}}\\mathrm{SI}\\sim \\mathcal{N}(0,1)$ \n\nSketch of the proofs: After sample splitting, the V-statistics are no longer degenerate, and hence follow a standard normal distribution under the null hypotheses. The existence of second moment assumptions are sufficient here and have been used as a standard assumption in kernel-based methods.\n\nTyping long and complex equations in OpenReview is cumbersome, so we have omitted the more technical details here, but we will include full proofs in the revised manuscript.\n\n**Re: causal paper linked:** We thank the reviewer for suggesting this interesting connection to testing the kernel mean embeddings of the counterfactual distributions. The counterfactual distribution and its kernel mean embeddings formalised in Singh et al. indeed presents us an opportunity to explore the high-order interactions associated with counterfactual outcome in a multivariate regression setting. For instance, it would be interesting to compare the interactions between the covariates and the outcomes before and after the intervention. Additionally, if there are multiple outcomes variables of interest, one can assess the interactions between the outcomes before and after interventions. We will discuss this interesting connection in the revised manuscript.", "claims_and_evidence": "I find the claim to be clearly presented and well supported.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "It would be nice to have a formal statement of the limiting distribution under the null. If I understand correctly, the authors describe it verbally and suggest it is immediate from their constructions and previous results, but it would be good to formalize it.\n\nI do not see any technical issues in the results that are currently presented, which seem like straightforward summaries of computational complexity.", "experimental_designs_or_analyses": "The experimental designs seem reasonable, and inherited from earlier works in the literature.", "relation_to_broader_scientific_literature": "The key techniques are from Shekhar et al. (2023) and Liu et al. (2023b)", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "While the result is not groundbreaking, it is clearly presented. The paper also makes a technical literature accessible in a nice way.\n\nThe main weakness is that the paper lacks a formal statement of the limiting distribution under the null. I would like the authors to address this point.", "other_comments_or_suggestions": "I think a natural connection could be to testing kernel mean embeddings of counterfactual distributions as formalized in \"Kernel methods for causal functions: dose, heterogeneous and incremental response curves\".\n\n“kerne” on page 4.", "questions_for_authors": "N/A", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "lFSA3AJUCK", "reviewer": "Reviewer_wudW", "summary": "This paper introduces permutation-free kernel-based tests for detecting high-order interactions in multivariate data. The proposed methods, xdHSIC, xLI, and xSI, leverage V-statistics and cross-centring to achieve a standard normal null distribution, eliminating computationally intensive permutations. Empirical evaluations demonstrate the efficiency and scalability of proposed methods while maintaining comparable statistical power. Applications to causal discovery and feature selection highlight the methods' scalability and utility in real-world scenarios involving complex dependencies.", "questions": "", "limitations": "", "rating": 2, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their review. Below we discuss the concerns raised by the reviewer. \n\n**Re: iid assumption:** We have added a discussion of how our method might be applicable to time series in the response to Reviewer vR6R. In summary, our method is still applicable if multiple iid realisations of a times series (stationary or non-stationary) are available. If only one time series is observed, then one might have to carry out some transformations, (e.g. stock returns in the real-world financial example below shown in Figure S3 in link: https://imgur.com/a/uju1pHq) or assume certain mixing conditions.\n\n**Re: real and large dataset:** We have added a large real financial dataset (504 variables and 1004 samples) to demonstrate the scalability and broad applicability of our proposed method. In summary, we show that in highly-regulated sectors such as Energy and Utilities, there are consistently high 3-, 4-and 5-way interactions due to redundancy. Moreover we observe that companies in Information Technology and Health Care are prone to external factors therefore have low intra-sector interactions. Importantly, these experiments took ~8 hours without any parallelisation on a 2015 iMac. We estimate the time taken for the permutation-based methods to be at least 2 weeks. See further details in response to Reviewer s8B8.\n\n**Re: only small $d$ considered, example of thousands of variables:** In a system consisting of $M$ (thousands of) variables, high-order behaviours are captured up to order $d << M$. In general, detecting $d$-order interactions among a group of $M$ variables is inherently combinatorial and quickly becomes infeasible as $M$ and/or $d$ increase.  However, recent studies across various application areas have focused on identifying interactions beyond pairwise ($d>2$), and have shown that even low-order interactions ($d=3,4$) significantly impact network structure and dynamics, see  Battiston et al. (2020) and Santoro et al. (2023). These findings underscore the practical benefits of high-order interaction tests, even to moderate $d$, in revealing new structural and functional properties that would be overlooked by pairwise analyses alone, as our financial example also highlights.\n\n**Re: the accuracy and computational complexity:** In our experiments, we observe that, when the sample size is low, the permutation-based methods achieve very similar performance, and in some cases slightly more powerful, than permutation-free methods (as noted by the reviewers for the causal discovery example in Fig. 6). We believe this is because in the low sample regime permutation-free methods are not able to capture the full structure of the data due to the data splitting—specifically, the test statistics are computed as the inner product between empirical embeddings estimated from two separate sample halves. Therefore, as a rule of thumb, if both the data available and the order of interest are relatively small, i.e., both $n$ and $d$ are small, then we could consider using the permutation-based methods. In such a case, the trade-off between computation time and statistical power would not be unfavourable. More concretely, when $n\\leq 50$ and $d\\leq 5$ the permutation-based method could be approximated well enough with just 100 permutations making it both computationally feasible and statistically powerful. However, outside of this regime we would strongly recommend the permutation-free method introduced in this paper.  \nWe note that we do not make additional assumptions compared with Liu et al., (2023b). Both methods rely on iid data.\n\n**Re: ablation studies:** Thank you for the suggestion. Constructing the test statistics using V-statistics not only enhances computational efficiency but is also essential before applying cross-centring to achieve normality under the null. Gretton et al., (2007) and Liu et al.,(2023b) employ V-statistics to enable traditional centring for these purposes. For the importance of cross-centring, please see the ablation study of the cross-centring technique in response to Reviewer s8B8.\n\n**Re: random fourier features:** We thank the reviewer for bringing up Random Fourier Features (RFF). RFF is a technique that constructs a low-dimensional feature map such that  its inner product is close to the true kernel matrix. As a result, it reduces the computational complexity of computing the kernel matrix from $O(n^2)$ to $O(nm)$ where $m$ is the number of RFFs. This could be advantageous, as it is linear in $n$, but the hyperparameter $m$ would have to be tuned to make tradeoffs between efficiency and accuracy. For each individual subtest, RFF could reduce the  $O(dn^2)$ to $O(dnm)$.  This means that in future work, we can also apply RFF on top of the permutation-free strategy (which has already eliminated the number of permutations $p$) in order to achieve better performance. \n\nWe will include this discussion in the revised manuscript and pursue it in our future research, with thanks.", "claims_and_evidence": "1. The methods assume i.i.d. samples, which is acknowledged as a limitation but not validated against dependent data (e.g., time series). \n2. The paper claims that the theoretical computational complexity of its permutation-free tests (e.g., xdHSIC, xLI, xSI) scales as O(dn^2), where d is the number of variables and n is the sample size. However, the experiments primarily focus on smaller values of d rather than demonstrating scalability to extremely large d or real-world datasets with thousands of variables.", "methods_and_evaluation_criteria": "1. The paper claims that the theoretical computational complexity of its permutation-free tests (e.g., xdHSIC, xLI, xSI) scales as O(dn^2), where d is the number of variables and n is the sample size. However, the experiments primarily focus on smaller values of d rather than demonstrating scalability to extremely large d or real-world datasets with thousands of variables. \n2. From Figure 6, it seems that the permutation-free methods require a larger sample size. The paper doesn't seem to discuss what factors the accuracy of the test methods is related to. It's also unclear whether the permutation-free methods can achieve the same accuracy as the permutation-based methods under reasonable assumptions or conditions. Additionally, since the computational complexity of the permutation-free methods is related to the sample size n, if more samples are needed to meet the accuracy requirements, the computational cost will also increase.", "theoretical_claims": "The paper introduces novel definitions (e.g., Definition 3.1 for xdHSIC) and hypotheses (e.g., Hypothesis 3.2 for joint independence) to enable permutation-free high-order interaction tests. While the theoretical framework is logically consistent and supported by mathematical derivations (e.g., V-statistics, cross-centring), the justification for these definitions and hypotheses could benefit from additional evidence to enhance their solidity. For example: rigorous validation of assumptions and ablation studies.", "experimental_designs_or_analyses": "The paper claims that the theoretical computational complexity of its permutation-free tests (e.g., xdHSIC, xLI, xSI) scales as O(dn^2), where d is the number of variables and n is the sample size. However, the experiments primarily focus on smaller values of d rather than demonstrating scalability to extremely large d or real-world datasets with thousands of variables.", "relation_to_broader_scientific_literature": "The paper positions these definitions as extensions of Shekhar et al. (2023).", "essential_references_not_discussed": "It would be valuable for the authors to elaborate on the connection between permutation-free tests and random Fourier feature-based methods in the context of high-order interaction detection, given their shared focus on computational efficiency and kernel approximation.", "other_strengths_and_weaknesses": "N/A", "other_comments_or_suggestions": "N/A", "questions_for_authors": "N/A", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper introduces permutation-free high-order interaction tests for joint independence and partial factorization of d variables. The null hypothesis states that the joint distribution factorizes into a product of marginals, and the alternative hypothesis states that it does not. Unlike state-of-the-art approaches such as HSIC which requires expensive permutations to sample from the null distribution, the new approach does not, and significantly saves computation.

The core proposal is built on previous works of Shekhar et al. (2023) and Liu et al. (2023b). The paper is well received by the reviewers, noting the novelty of the proposed permutation-free test and its well-grounded theory (s8B8, vR6R, eBnj). The claim on computation saving is clearly presented and supported by the empirical results. Initial concerns raised in the official reviews that **were addressed** include 

1. lack of sensitivity analysis on kernel choice [s8B8], 
2. lack of experiments on real datasets [vR6R], 
3. limited discussion on when permutation-based approaches can be better [vR6R, wudW]
4. lack of experimental results on non-i.i.d. data [wudW]

The AC sees the last concern as a minor one, considering other contributed results. The authors **submitted a strong rebuttal** that addresses all the above concerns. In the next revision, please include the new results on the financial dataset, structural causal models, and other results from the rebuttal.

AC recommendation: accept.
Decision: Accept (poster)