Efficiently Learning Significant Fourier Feature Pairs for Statistical Independence Testing

Authors: yixin  ren
yewei  xia
hao  zhang
jihong  guan
shuigeng  zhou
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-12 19:40:43

Abstract: We propose a novel method to efficiently learn significant Fourier feature pairs for maximizing the power of Hilbert-Schmidt Independence Criterion~(HSIC) based independence tests. We first reinterpret HSIC in the frequency domain, which reveals its limited discriminative power due to the inability to adapt to specific frequency-domain features under the current inflexible configuration. To remedy this shortcoming, we introduce a module of learnable Fourier features, thereby developing a new criterion. We then derive a finite sample estimate of the test power by modeling the behavior of the criterion, thus formulating an optimization objective for significant Fourier feature pairs learning. We show that this optimization objective can be computed in linear time (with respect to the sample size $n$), which ensures fast independence tests. We also prove the convergence property of the optimization objective and establish the consistency of the independence tests. Extensive empirical evaluation on both synthetic and real datasets validates our method's superiority in effectiveness and efficiency, particularly in handling high-dimensional data and dealing with large-scale scenarios.

Main paper: 
 introduction: 
 Testing for independence is a crucial and challenging task in machine learning and statistics, with wide-range applications in causal inference [16, 31] , feature selection [6] and deep learning [23, 42] . Its primary objective is to determine whether two random variables, X and Y are independent, based on the observations of the underlying joint distribution P XY . While traditional independence tests, such as Pearson's correlation coefficient [9] and Kendall's τ , can only detect monotonic relationships between low-dimensional variables, more modern tests [26, 43, 7, 25, 27, 35, 19, 20] aim to deal with complex non-linear interactions in much more challenging higher-dimensional space [45, 29] .One class of nonlinear dependence measures [3, 15] aims to capture distributional characteristics using kernel embeddings [13] , primarily derived from the cross-covariance operators in the reproducing kernel Hilbert space (RKHS). Among them, Hilbert-Schmidt Independence Criterion (HSIC) [14] is the most popular one. It utilizes the squared Hilbert-Schmidt norm to detect dependence and exhibits outstanding performance across various data contexts by choosing suitable kernels. On the other hand, some other fundamental nonlinear dependence measures employ characteristic functions to detect the smoothed discrepancy between the joint distribution and the product of marginals. By employing appropriate characteristic functions, the statistic [39, 40] computes the covariance between distances of variable pairs. It has been demonstrated that these distance-based methods are equivalent to HSIC with specific kernels [33] . However, all these measures suffer from the drawback of requiring quadratic time (w.r.t. the sample size n) to compute the feature covariance and necessitating fixed kernel or distance functions, rendering them impractical on large-scale datasets due to the unaffordable time cost and lacking flexibility in handling complex scenarios.To address these challenges, a multitude of works grounded on these measures have emerged. Upon HSIC, [44] proposes some linear-time tests including a block-averaged statistic, a statistic with Nyström approximation, and one with finite-dimensional feature mappings using random Fourier features (RFF) [28] . For convenience, these tests are referred to as BHSIC, NyHSIC, and FHSIC, respectively. FHSIC and NyHSIC are observed to have a considerable advantage over BHSIC. However, a remaining drawback of these methods is that the features are not learnable. Therefore, these methods lack enough adaptability to complex settings, thus leading to performance degradation.In addition to time efficiency, another research direction [1, 30] aims to make independence tests adaptive to better capturing distributional distinctions. These methods either select/combine appropriate kernels from a predefined set or learn parameterized kernels. Nonetheless, their criteria still inherit the quadratic time complexity of HSIC, thus cannot be readily applied to large-scale data. independent and identically distributed (i.i.d) samples Z := (X, Y ) = {(x i , y i )} n i=1 with distribution P XY , we aim to test whether X, Y are independent (i.e., X ⊥ ⊥ Y ). This corresponds to a hypothesis testing problem formulated as H 0 : P XY = P X P Y versus H 1 : P XY ̸ = P X P Y .The testing procedure is as follows: First, define the statistic ρ and calculate its estimated value using the samples. Then, choose a significance level α (typically set to 0.05), which represents the probability that the sampling of ρ under H 0 is at least as extreme as the observed value. Finally, the null hypothesis H 0 is rejected if the p-value is not greater than α.Two types of errors may occur in this procedure. Type I error occurs when H 0 is falsely rejected, while Type II error happens when H 0 is incorrect but not rejected. A good test [43] needs to control Type I error within α while maximizing the testing power (1-Type II error rate).For independence tests, a commonly used statistic is HSIC, defined as follows: Definition 1. [14] . Let F be an RKHS with kernel k : X × X → R and G be a second RKHS on Y with kernel l : Y × Y → R, the HSIC between X and Y , denoted as HSIC(X, Y ) is defined asE k(X, X ′ )l(Y, Y ′ ) +E k(X, X ′ ) E l(Y, Y ′ ) -2E X ′ Y ′ E X k(X, X ′ )E Y l(Y, Y ′ ) , (1)where (X ′ , Y ′ ) is a independent copy of (X, Y ). An estimator of HSIC(X, Y ) is given byEQUATIONwhere k ij := k(x i , x j ), l ij := l(y i , y j ) are the entries of the n × n kernel matrices K, L respectively, H = I -1 n 11 T is the centering matrix and 1 is a vector of ones. 
revisiting hsic from frequency domain perspective: 
 We denote F as the Fourier transform, and F -1 as its inverse. When the kernels k, l are translationinvariant, i.e., there exist functions ψ, ψ k , ψ l such that for all (x, x ′ ) ∈ X × X and (y, y ′ ) ∈ Y × Y, ψ(x -x ′ , y -y ′ ) = ψ k (x -x ′ )ψ l (y -y ′ ) = k(x, x ′ )l(y, y ′ ).Then, according to the results of [36, Corollary 4] , the HSIC with function ψ can be formulated as HSIC(X, Y ) = R dx ×R dy ϕ P XY (ω) -ϕ P X P Y (ω)EQUATIONwhere ω = (ω x , ω y ) ∈ R dx × R dy , ω x , ω y are the frequencies of X and Y respectively, and ϕ P XY (ω) := e -i(ω T x x+ω T y y) dP XY , ϕ P X P Y (ω) := e -iω T x x dP X e -iω T y y dP Y (5) are the characteristic functions of P XY and P X P Y , respectively. Intuitively, Eq. ( 4) means that HSIC can be understood as the difference between the joint distribution and the product of the marginal distributions in the frequency domain, with different weights (F -1 ψ)(ω) being attached to different frequencies, which are determined by the kernel function. When F -1 ψ is almost everywhere nonzero, it can be shown that the kernel is characteristic [36, 10] . The characteristic condition ensures that the criterion is discriminative for discrepancies at almost all frequencies. However, with inappropriate choices of F -1 ψ, the differences may not be significant enough. We explain this with an example:Example. Consider the Sinusoid model that X × Y := [-π, π] 2 and (X, Y ) ∼ p xy (x, y) ∝ 1 + sin(ω 0 x) sin(ω 0 y), where p xy is the probability density function and ω 0 is a positive integer. Combining Eq. ( 5), we can calculate that ϕ P X P Y (ω) = δ(ω x )δ(ω y ) and ϕ P XY (ω) = δ(ω x )δ(ω y ) + [δ(ω x + ω 0 ) + δ(ω x -ω 0 )][δ(ω y + ω 0 ) + δ(ω y -ω 0 )], where δ is the Dirac delta function, thus the difference between them only relies on the frequency ω 0 . When the Gaussian kernels with width √ 2λ x and √ 2λ y are used, i.e., k(x, x ′ ) = exp(-∥x -x ′ ∥ 2 2 /(4λ 2 x )), l(y, y ′ ) = exp(-∥y -y ′ ∥ 2 2 /(4λ 2 y )), then the inverse Fourier transform of ψ is (F -1 ψ)(ω x , ω y ) = π -1 λ x λ y exp(-(λ 2x ω 2 x + λ 2 y ω 2 y )). Hence HSIC(X, Y ) = 4π -1 λ x λ y exp(-(λ 2x + λ 2 y )ω 2 0 ) whose maximum is taken at λ * x = λ * y = 1/( √ 2ω 0 ), indicating that the widths need to be adjusted to focus on some specific frequencies. If the common setting [14, 44] is adapted, which uses mid-widths (i.e., the median distance does not change with ω 0 since the marginal distributions do not change with ω 0 ), then the criterion will exponentially decline to 0 as ω 0 increases. In contrast, the criterion using the adaptive optimization width (1/ω 0 , 1/ω 0 ) decreases at a rate of O(ω -2 0 ), which is a considerable improvement. This example illustrates the loss of the discriminatory power of the criterion when an inappropriate F -1 ψ is chosen. The discriminatory power of the criterion heavily impacts the sample size required for the test to obtain significant results in practice, and existing inflexible configurations may lead to inadequate test power in the presence of reasonably large sample sizes. Consequently, it is important to design learnable F -1 ψ. To this end, we subsequently design a learnable objective and let it be optimized in a data-driven manner. Before this, we provide an approach to make the criterion be computed efficiently. This can be achieved by sampling in the frequency domain. Formally, a finite-dimensional approximation in the frequency domain of the integral in Eq. ( 4) is given as follows:HSIC ω (X, Y ) := 1 D x D y Dx i=1 Dy j=1 ϕ P XY (ω x;i , ω y;j ) -ϕ P X P Y (ω x;i , ω y;j )2 ,where {ω x;i } Dx i=1 , {ω y;j } Dy j=1 are sampled independently with the measure F -1 ψ k , F -1 ψ l , respectively. Note that F -1 ψ is a product measure, i.e., F -1 ψ = (F -1 ψ k ) ⊗ (F -1 ψ l ). This type of approximation is also called random Fourier features (RFF) [28] that had been applied to various kernel algorithms. We will incorporate this technique to efficiently perform computation later.4 Learning Significant Fourier Feature Pairs 
hsic with learnable fourier feature pairs: 
 To design F -1 ψ, we need to make sure that supp(F -1 ψ) = R dx × R dy to meet the characteristic condition and that its integral over the full space is 1 to ensure it is a probability measure. Also, for practical utility, F -1 ψ should embody a familiar probability density function, facilitating sampling procedures. Fortunately, a versatile array of options emerges through the judicious selection of kernels 2 with adjustable parameters. Take kernel k as an example, some commonly used kernels are listed in Tab. 1, and their inverse Fourier transforms are listed simultaneously. Additionally, to Table 1 : Some popular kernels (parameterized by σ, Σ) with corresponding density functions.Kernel ψ k (∆) F -1 ψ k (ω) T θ k (x) p k (ω) Gaussian e - ∥∆∥ 2 2 2σ 2 (2π) -dx/2 σe -σ 2 ∥ω∥ 2 2 /2 x/σ (2π) -dx/2 e -∥ω∥ 2 2 /2 Laplace e -∥∆∥ 1 σ 2 π d σ σ 2 +ω 2 d x/σ 2 π d 1 1+ω 2 d Mahalanobis e -1 2 ∆ T Σ -1 ∆ (2π) -dx/2 |Σ| -1/2 e -ω T Σ -1 ω/2 Σ 1/2 x (2π) -dx/2 e -∥ω∥ 2 2 /2be able to apply gradient-based optimization techniques, we invoke a method that disentangle the sampled objects and the learnable parameters. Specifically, we leverage a variable transform T θ k (a bijection function parameterized with θ k ) to convert the probability measure F -1 ψ k into a simple distribution (e.g. a standard Gaussian distribution) p k (ω). Simultaneously, we relocate the learnable component onto X. Consequently, we can focus on learning parameterized transformations T θ k and simplifying the computation by enabling sampling directly from p k (ω).The above scheme provides a broader form for designing. The mapping T θ k can be viewed as a feature extractor, which makes it possible to flexibly combine models (e.g., neural network) thus incorporating deep kernel [24] into the framework. Also, it should be noted that the single kernel example can also be extended to multi-kernel setting [11] by executing the procedure for each kernel.Next, we obtain the learnable independence criterion and utilize the sampling technique as in Eq. ( 6) to compute efficiently. Note that for simplicity, we take the same value for both D x and D y in Eq. ( 6) by default. By the definition, the kernel function can be expressed asEQUATIONBy applying the frequency sampling technique, we obtain the approximation asEQUATIONwhere {ω k;j } D/2 j=1 are sampled independently with distribution p k (ω) and the last equation is because the kernel function is real. To get a more computationally tractable form, we defineEQUATIONcalled learnable RFF of k then Eq. ( 8) becomes ψ(ω) k (T θ k x -T θ k x ′ ) = Λ k (x)Λ k (x ′ ) T .The expression with a similar form is also given in [44] , with the difference that we have added learnable parts. For Y , we define the corresponding symbols by substituting k for l and x for y. Also, for convenience, we default to keeping Y and X the same number of samples D from here on. Then the HSIC with learnable RFF pairs can be obtained by replacing k, l in Eq. ( 1) to ψ(ω) k , ψ (ω)l . Also, the corresponding estimator with sample Z can be obtained by replacing K, L in Eq. ( 2) to the matricesΛ X Λ T X , Λ Y Λ T Y , where Λ X := [Λ k (x 1 ); ...; Λ k (x n )]n×D and so as define for Λ Y . As a result,EQUATIONwhere Λ Xc := HΛ X , Λ Y c := HΛ Y . The time complexity is analyzed as follows. Since the computation of the mapping T θ k x depends on the specific design, here we default to analyzing the kernel case shown in Tab. 1. In this case, computingΛ X , Λ Y requires O nD(d x + d y ) time. Then calculate Λ Xc , Λ Y c cost O(nD). After that, calculate HSIC ω (Z) cost O(nD 2 ). Hence, the overall time complexity is O nD(d x + d y + D) , i.e. the running time is linear with n.Next, we model the behavior of HSIC ω (Z) to obtain an optimization objective for maximizing the power of the test. By utilizing the property that HSIC ω (Z) is a V-statistic, we can extend the results [14, Theorem 1, 2] for HSIC ω (Z), as shown in the following proposition with the proof given in the Appendix. To simplify, we denote (x i , y i ) as z i to represent the i-th sample and denote ψ(ω) k (T θ k x t -T θ k x u ) as k (ω) tu and ψ (ω) l (T θ l y t -T θ l y u ) as l (ω) tu . Proposition 1 (Asymptotics). Let h (ω) ijqr := 1 4! (i,j,q,r) (t,u,v,w) k (ω) tu l (ω) tu + k (ω) tu l (ω) vw -2k (ω) uv l (ω)tv , where the sum represents all ordered quadruples (t, u, v, w) drawn without replacement from (i, j, q, r). Then, Under the null hypothesis H 0 , HSIC ω (Z) coverages in distribution toEQUATIONwhere χ 2 11 , χ 2 12 , ... are independent χ 2 1 variates and λ l is the solution to the eigenvalue problem as in the right of Eq. (11) . Also, under the alternative H 1 , HSIC ω (Z) converges in distribution asEQUATION)with the simplified notation E j,q,r := E zj ,zq,zr and E Z := E zi,zj ,zq,zr .According to Proposition 1, the power of the test with HSIC ω can be formulated byEQUATION)where Φ is the standard normal CDF and r ω is the threshold, i.e. (1 -α)-quantile of distribution given in Eq. ( 11) that exactly controls Type I error rate to the nominal level α. Hence, to maximize the power of the test, a natural criterion is[nE Z HSIC ω (Z) -r ω ]/( √ nσ ω ).Next, we provide its estimation which can be computed in linear time.We first consider obtaining the estimator of the numerator part. For the term E Z HSIC ω (Z), we can estimate it with HSIC ω (Z) as in Eq. (10) . The estimation of the threshold r ω poses a challenge, primarily stemming from the lack of an explicit expression for the distribution of the infinite sum of chi-square variables. One avenue to address this challenge involves employing the permutation method [2, 38] to simulate the distribution under H 0 . However, this method necessitates a significant number of shuffles to accurately approximate the distribution. Furthermore, even with the implementation of parallel schemes, it incurs memory costs proportional to the number of permutations, rendering it impractical for resource-constrained scenarios. Here, we adopt a lightweight approach in practice, leveraging the gamma approximation as proposed by [14] . A gamma distribution is uniquely determined by its first and second-order moments. For these two moments, we present their corresponding linear-time estimators in Theorem 1. As a result, we can obtain the (1 -α)-quantile of the gamma distribution, denoted as c α , with estimated parameters γ := E 2 0 /V 0 , β := V 0 /E 0 in linear time. Formally, with the term E 0 and V 0 defined in Theorem 1, c α is calculated byEQUATIONwhere Γ(•) is the gamma function. By combining the way to estimate the gradients of c α [30] , we enable it for gradient-based optimization with automatic differentiation framework. As a result, we obtain a linear-time differentiable estimator of the numerator part. Theorem 1 (Linear-Time Estimators). Under H 0 , the estimators of mean and variance with bias ofO(n -1 ) to E Z [nHSIC ω (Z)] and Var Z [nHSIC ω (Z)],denote as E 0 and V 0 , respectively, are given byEQUATIONwhere () •2 is the entry-wise matrix power. Both E 0 and V 0 can be calculated in O(nD 2 ) time.For the remain term σ ω , we estimate it with σ ω that σ 2 ω := 161 n i ( 1 n 3 j,q,r h (ω) ijqr ) 2 -HSIC 2 ω (Z) . To calculate j,q,r h (ω)ijqr , the straightforward way is to compute each item h (ω) ijqr , which requires total O(n 4 ) of computation. Here we provide a way to enable it to be calculated in linear time by obtaining a matrix expression. The main result is given by j,q,r hEQUATIONwhere the definition of variables A to F with the calculation cost are given in the Fig. 1 and the derivation of Eq. ( 16) is given in the Appendix. By checking the complexity of the remaining matrix operations in Eq. ( 16), all the elements with index i can be calculated in O(nD 2 ). Combining the results obtained before that HSIC ω (Z) can also be calculated in O(nD 2 ), thus calculating the term σ ω cost O(nD 2 ) time. As a result, we obtain the overall linear-time optimization objective J := [HSIC ω (Z) -c α /n]/ σ ω , which is a clear contrast to the existing quadratic-time schemes [30] .Figure 1: The diagram shows the definition of the quantities in Eq. ( 16), with styles representing the time complexity of the computational process in the current box. ⊙: the element-wise product.After obtaining the differentiable optimization objective J, we can perform the training process end-to-end. In this process, the overfitting issues may happen especially with insufficient samples, which could influence both Type I and II errors. If we use the same sample for testing, the Type I errors may be uncontrollable [22] when the overfitting issues happen. To address this, we adopt the split scheme as in [24, 18] to allow our tests to maintain validity (controllable Type I errors). The split ratio is set to 0.5 to facilitate the balance between the two. Apart from controlling Type I errors, we still want to mitigate overfitting issues as much as possible in order to generalize the optimized Fourier feature pairs on test data thus improving the power of our tests. To this end, we select smooth function classes to control the model complexity (e.g., as measured by the VC dimension), specifically in this paper, we consider the two classes in Tab. 1 and implement them for experiments later. One is the choice of Gaussian classes that optimize the global scale, and the other we consider Mahalanobis classes and set Σ to be diag(σ 1 , ..., σ d ) for optimization, which corresponds to optimizing the scale in each dimension (which allows to capture high-frequency signals such as the edge in the image). These smooth choices also bring the advantage of interpretability [30] and it is experimentally proven that this simple choice is already able to handle most of the cases in different settings.More discussion about split strategy. Currently, there are two major classes of approaches for adaptive independence tests. One involves selecting kernels from a finite/countable set (discrete scenario) and the other involves performing kernel parameter searches in a continuous space (continuous scenario). For the former case, some methods [22, 32] control Type I errors by applying techniques from the selective inference literature without data splitting. However, these methods cannot be directly applied to a continuous scenario due to the uncountable set of kernels involved. To the best of our knowledge, both our scheme and existing methods [24, 18] rely on data splitting for the continuous case. Designing methods to control Type I errors in the continuous case without sample splitting remains a challenging and significant problem for future research.Algorithm. Our algorithm is outlined in Alg. 1. As a pre-processing step, we split the data into the training data Z tr and the testing data Z te (Line 1). The test contains two phases: 1) We learn the Fourier feature pairs with Adam [21] optimizer using full batches on Z tr (Lines 2-7). 2) With the learned Fourier feature pairs, we calculate the test statistic and threshold (Lines Obtain learnable Fourier feature pairs Λ X , Λ Y with parameters θ k , θ l and {ω j } D/2 j=1 .Calculate criterion J with Λ X , Λ Y then optimize J with (θ k , θ l ) ← (θ k , θ l ) + ϵ∇ (θ k ,θ l ) J. X ̸⊥ ⊥ Y if c α (Z te ) ≤ n te HSIC ω (Z te ) holds, otherwise X ⊥ ⊥ Y . 
theoretical results: 
 We first give the uniform bound results over a ball in parameter space which guarantees the convergence of our optimizing objective thus ensuring its effectiveness in modeling test power.Theorem 2 (Uniform Bound). Let θ k , θ l parameterize T θ k , T θ l in Banach spaces of dimension d k , d l .And T θ k , T θ l are Lipschitz to the parameters θ k , θ l with the non-negative constant L k , L l , respectively. Let Θ c be a set of (θ k , θ l ) for which σ ω ≥ c > 0 with a positive constant c and ∥θ k ∥ ≤ R θ k , ∥θ l ∥ ≤ R θ l . Let r denote the threshold, i.e., (1 -α)-quantile for the distribution in Eq. (11) and r (n) be the threshold with sample size n. Let {(ω k;j , ω l;j )} D/2 j=1 be the samplings of frequency with the sampling number D. Also, we define R ω k := sup j ∥ω k;j ∥, R ω l := sup j ∥ω l;j ∥, d s := max{d k , d l } and ξ ω := HSIC ω (Z). Then with probability at least 1 -δ, we havesup (θ k ,θ l )∈Θc ξ ω -r (n) ω /n σ ω - E Z ξ ω -r ω /n σ ω ∼ O 1 n log 1 δ + d s log n n + R ω k L k + R ω l L l √ n .Next, we show the consistency of the tests, i.e. the power of the test tends to 1 as the sample size increases. Let the U-statistic of HSIC ω (Z te ) be HSIC (u) ω (Z te ), then we have the following results. Theorem 3 (Consistency). Let θ * k , θ * l be the parameters after learning, Z te be the test samples of size m, when E Z HSIC (u) ω (Z te ) > 0, then the probability of the Type II errorEQUATION)Let the mapping functions with learned parameters θ * k , θ * l be T θ * k , T θ * l , and the corresponding range space be compact subsets of R d Tx , R d Ty , respectively. Also, the diameters of two range spaces are denoted by diam(T θ * k ), diam(T θ * l ), respectively. Let {(ω k;j , ω l;j )} D/2 j=1 be the frequency samplings with their second moment denoted by σ 2ω k := E p k (ω) [ω T k;j ω k;j ], σ 2 ω l := E p l (ω) [ω T l;j ω l;j ]. Additionally, we denote ξ u := HSIC(X, Y ), then under H 1 , we have E Z HSIC (u) ω (Z te ) > 0 with any constant probability when D = Ωd Tx +d Ty ξ 2 u log σω k diam(T θ * k )+σω l diam(T θ * l ) ξu .This result can be understood in two parts. The first one is about consistency, i.e., the Type II error rate tends to 0 at the rate of m -1/2 when condition E Z HSIC (u) ω (Z te ) > 0 holds. The second part provides the condition when E Z HSIC (u) ω (Z te ) > 0 holds, which requires sufficiently many frequency samplings. The theorem shows that the large value of the criterion HSIC(X, Y ) helps to reduce the required D. According to the results discussed in Sec. 3, there is an improvement in the criterion by finding the more significant features and thus helps to reduce the required D. To summarize, the significant features further help to guarantee the consistency of the test under the efficient requirements (smaller D). All proofs as well as additional results are given in the Appendix.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: We provide the experimental data/code in the supplemental material. Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
performance evaluation: 
 We compare the following tests: distance-based statistic dCor [39] , the original HSIC QHSIC [14] , the copula-based method RDC [26] , the three variants of HSIC NyHSIC [44] , FHSIC [44] , BH-SIC [44] and HSICAgg [32] , NFSIC [18] as introduced in Sec. 1. Among them, dCor and QHSIC are O(n 2 ) tests. RDC is calculated in O(n log n) time and the rest are O(n) tests. A detailed description of the comparing methods is given in the Appendix. For our methods, We provide two variants as mentioned in Sec. 4.3. We name the Gaussian class case as LFHSIC-G, and name the Mahalanobis class case (and set Σ as a diagonal matrix) LFHSIC-M. Additionally, for the comparative methods [30] that are relevant to us, due to their high time overhead and therefore inability to handle some settings of evaluation, we separately provide a comparison with our method under certain feasible experimental settings, the results are given in the Appendix.Experimental setup. The significance level α is set to 0.05. We use Gaussian kernels for both X and Y in all kernel-based methods. And QHSIC, RDC, NyHSIC, FHSIC, BHSIC are all with the kernel width being set to the Euclidean distance median of the samples. The number of random features D for FHSIC, LFHSIC-G/M, the number of induced variables for NyHSIC, the block size for BHSIC as well as the number of sub-diagonals R for HSICAgg are all kept consistent as recommended in [44, 32] for fair evaluation. Parameter settings for the rest of the methods follow the defaults in the code. More details of the setups are given in the Appendix.Evaluation protocol. We evaluate on four synthetic datasets [18, 30] and two real datasets [44, 30] . Synthetic datasets consist of Sine Dependency (SD), Sinusoid (Sin), Gaussian Sign (GSign), and independent subspace analysis (ISA) dataset [14] . On real data, we introduce high-dimensional image data and another music dataset to evaluate the capability of all methods in different data scenarios. Unless otherwise specified, we perform 100 repeated randomized experiments and report the average result of test power as default. More details of the generating process of each dataset and the details of the evaluation (including running time) are provided in the Appendix.Settings of SD, Sin, and GSign Dataset. The Sin data corresponds to the example in Sec. 3 that requires the method to focus on differences in specific frequencies. In SD, Y is dependent solely on the first two dimensions of X. In contrast, in GSign, Y is independent of any proper subset of X but dependent on X as a whole. Therefore, it requires the method to learn important local/global features based on the characteristics of the data to improve the test power. For SD and GSign, we set the dimension of X as 4 and 5, respectively, and the dimension of Y is 1 for both. For Sin, we set the frequency parameter ω = 5. For calculating the Type I error rate, we evaluate using samples (n = 2000) obtained by permutation for all three datasets.Performance. The results for D = 100 and D = 500 are shown in Fig. 2 . Except for NFSIC and BHSIC, all the other methods succeed in controlling the Type I error rate ≤ 0.05. LFHSIC-M/G, NFHSIC, and HSICAgg perform much better than other methods due to their ability to obtain more appropriate kernels/features for testing. LFHSIC-G/M performs on both settings of D and has a more significant advantage over the others when D is small, implying the optimization objective can still be successfully optimized and the criterion is still powerful under high-speed requirements. In addition, as the sample size increases the test power of LFHSIC-G/M is gradually converging to 1 in both settings, which corroborates the results of Theorem 3. results under different settings of D, our method LFHSIC-M performs consistently well and exhibits progressively better performance as D increases.Settings of Two Real Data. The first real dataset used is a high-dimensional image dataset 3Dshapes as in [30] . In our experiments, we vectorize image X to a vector with dimension 64×64×3 = 12, 288. The sample size is set as 128. We add standard Gaussian noise N (0, 1) to the angle Y to make the setting more challenging. The Type I error rate is evaluated by the samples obtained by permutation. Besides, we consider the Million Song Data (MSD) as the second real dataset. The first dimension represents the year of release of each song and is referred to as variable Y . The remaining 90dimensional features (e.g., mean timbre and timbre covariance) constitute variable X. We follow the recommended setting [44] , i.e., disturbing each entry of the X with an independent Gaussian noise N (0, 1000). For this dataset MSD, in order to fully utilize the data, we randomly select n ∈ {500, 1000} samples as the training set and other n samples from the remaining data 100 times for the evaluation and obtain the average result. The above training and testing processes are repeated 10 times to evaluate the robustness of the optimization scheme.Performance. The results of two real data with D = 10 are presented in Fig. 4 . For the results on 3Dshapes (shown in the left of Fig. 4 ), all methods except BHSIC and NFSIC control the Type I error well. The linear-time test has relatively lower power compared to the quadratic-time test except for LFHSIC-M, proving that its more significant features obtained in high-dimensional scenarios enable it to achieve outstanding performance even in scenarios with high approximation requirements (D = 10). Similar conclusions can drawn from the MSD dataset (shown in the right of Fig. 4 ). Additionally, the results for NFSIC and LFHSIC-G/M with different sample sizes indicate increased robustness of the optimization as the sample size increases (reflected in the reduction of variance), and the more flexible design also contributes to this (comparing LFHSIC-M and LFHSIC-G), thus can be more effectively applied to real-world scenarios.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes]Justification: we provide the details. See Sec. J in Appendix.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusion: 
 In this paper, we propose a novel method to efficiently learn significant Fourier feature pairs for maximizing the power of HSIC-based independence tests. By integrating a learnable Fourier feature module, we improve the flexibility of existing configurations and design a new criterion. The proposed linear-time optimization objective accurately models the power of the test and can be trained end-toend in a data-driven manner, ensuring both effectiveness and efficiency. Both theoretical results and experimental results show the effectiveness of our proposed method. Future work includes further improving the sampling method in the frequency domain.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: The experimental results are accompanied by statistical significance tests. See Sec. 6.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes]Justification: We provide sufficient information on the computer resources. See Sec. K in Appendix.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 9. Code Of Ethics Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: I read it. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: See the Sec. L in the Appendix. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: [TODO]Answer: [NA] Justification: [TODO] Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: • Section A: List of Symbols and Notations.• Section B: Assumptions.• Section C: Some Auxiliary Lemma.• Section D: Proof of Proposition 1.• Section E: Proof of Theorem 1.• Section F: Calculation of Eq. ( 16).• Section G: Proof of Theorem 2.• Section H: Proof of Theorem 3.• Section I: Smoothness of Optimization Objective.• Section J: Details of Experiment Setup.• Section K: Additional Experiment Results.• Section L: Limitations and Broader Impacts.O big O notion o small O notion i.i.d.independent and identically distributed R the set of real numbersB(R) Borel σ-algebra on R P X marginal distribution of X P XY joint distribution of X, Y F X distribution function of X E[X] expectation of X Var(X) variance of X X ⊥ ⊥ Y random variables X, Y are independent X ̸⊥ ⊥ Y random variables X, Y are not independent i n rthe set of all r-tuples drawn without replacement from the set {1, ..., n}n k number of k-combinations of n elements (n) k number of permutations, define as n! (n-k)!Tr(•) the trace of a square matrix 1 an vector of all ones H centering matrix define asH = I -1 n 11 T ⊙ element-wise product () •2 element-wise power d - → convergence in distribution ⊗ the product symbol of measure × the product symbol of topological space N (Θ, r)covering number with radii r for space Θ F, F -1 the Fourier transform, Fourier inverse transformThe following are the assumptions required. We denote the parameter spaces of θ k , θ l as Θ k , Θ l .(a) The mapping functions T θ k , T θ l are Lipschitz to the parameters θ k , θ l , i.e. for all x ∈ X , y ∈ Y and for allθ k , θ ′ k ∈ Θ k , θ l , θ ′ l ∈ Θ l , T θ k (x) -T θ ′ k (x) ≤ L k • ∥θ k -θ ′ k ∥, T θ l (y) -T θ ′ l (y) ≤ L l • ∥θ l -θ ′ l ∥ (18) with the nonnegative Lipschitz constant L k , L l . (b)The range of the mapping functions T θ k , T θ l are bounded. (c) The parameters θ 0 , θ 1 lie in Banach spaces of dimension d k , d l respectively. Also, the parametersθ k , θ l are bounded by R θ k , R θ l respectively, i.e., ∥θ k ∥ ≤ R θ k , ∥θ l ∥ ≤ R θ l .In this part, we give a useful expression for h(ω)ijqr for the subsequent proof. By the definition, we have h(ω) ijqr = 1 4! (i,j,q,r) (t,u,v,w) k (ω) tu l (ω) tu + k (ω) tu l (ω) vw -2k (ω) tu l (ω)tv . We simplify it by setting t = i, u = i, v = i and w = i in turn. Then we can show that h (ω) ijqr is equal toEQUATIONBy the definition, k(ω)tu := ψ (ω) k (T θ k x t -T θ k x u ) = Λ X (x t )Λ X (x u ) T is symmetric, i.e. k (ω) tu = k (ω)ut . And so as l (ω) tu . Hence we can merge the identical items (marked with the same color). As a result,EQUATIONtu + 2kEQUATION(EQUATIONWe will use Eq. ( 20) many times in subsequent proofs.Under the assumption (a), RFFs ψ(ω) k , ψ (ω) lare Lipschitz to the parameters θ k , θ l . Formally, Lemma 1. (Lipschitz Property of Fourier Feature). Let T θ k , T θ l be the mapping functions of X, Y that are Lipschitz to the parameters θ k , θ l with the non-negative constant L k , L l , respectively. Let {(ω k;j , ω l;j )} D/2 j=1 be the samplings of frequency with the sampling number D. Also, we define R ω k := sup j ∥ω k;j ∥, R ω l := sup j ∥ω l;j ∥, then for the RFFs ψEQUATIONwhere∆ x,x ′ := T θ k x -T θ k x ′ , ∆ ′ x,x ′ := T θ ′ k x -T θ ′ k x ′ and ∆ y,y ′ , ∆ ′ y,y ′ are defined by analogy.Proof. We prove the result for ψ(ω) k only since the proof for ψ(ω) l can be obtained in the same way. We start by recall the definition ψEQUATIONSince the cosine function is bounded by 1, by the mean value theorem, for fixed j, we haveEQUATION)Then according to the Cauchy-Schwarz inequality,EQUATIONBy the definition of ∆ x,x ′ and the Lipschitz property of the mapping functions T θ k , T θ l , we haveEQUATIONCombining the above results, we complete the proof.Under the assumption (b), we can obtain the uniform convergence property as follows. Lemma 2. (Uniform Convergence of Fourier Features). Let the mapping function of X, Y with parameters θ k , θ l be T θ k , T θ l , and the corresponding range space be a compact subset of R d Tx , R d Ty , respectively. Also, the diameter of two range spaces is denoted by diam(T θ k ), diam(T θ l ), respectively. Let {(ω k;j , ω l;j )} D/2j=1 be the samplings of frequency with the sampling number D, then for the RFFs with mapping functions T θ k , T θ l and frequency samplings {(ω k;j , ω l;j )} D/2 j=1 , we haveEQUATIONwhere the second moment of frequency samplings σ 2ω k := E p k (ω) [ω T k;j ω k;j ], σ 2 ω l := E p l (ω) [ω T l;j ω l;j ].Proof. Based on the derivation of RFFs in Sec. 4.1, We can view it as if the frequency sampling process is performed after the range space is obtained. Since the convergence bounds of the sampling process can be obtained directly through the results of [28, Claim 1], by replacing the input space in [28, Claim 1] to the range space here, then this part of the proof can be completed.Remark. Combining the technique in [37] , the constants in bounds can be further improved.Let HSIC (u) ω (Z), also denoted as ξ(u)ω , be the U-statistic that corresponding to HSIC ω (Z), i.e., HSIC (u) ω (Z) := 1(n)4 (i,j,q,r)∈i n 4 h (ω)ijqr . The population value of HSIC (u) ω (Z) is given by E Z ξ (u) ω which can be viewed as the result obtained after a frequency sampling approximation on HSIC(X, Y ) is performed. The bound of approximation error is given by the following Lemma. Lemma 3. (Approximation Error Bound). For simplify, we denoteEQUATION(y, y ′ ), respectively. Then we haveEQUATIONProof. We first representE Z ξ (u)ω in the form corresponding to Eq. ( 1), i.e.,EQUATION)Taking one of the items as an example and comparing it to the corresponding item in Eq. ( 1),EQUATIONThe results can be obtained for the other terms in a similar way, which completes the proof.In this section, we give a proof of the Proposition 1. We first restate the Proposition 1 here.(ω)ijqr := 1 4! (i,j,q,r) (t,u,v,w) k (ω) tu l (ω) tu + k (ω) tu l (ω) vw -2k (ω) uv l (ω)tv , where the sum represents all ordered quadruples (t, u, v, w) drawn without replacement from (i, j, q, r). Then, Under the null hypothesis H 0 , HSIC ω (Z) coverages in distribution toEQUATIONwhere χ 2 11 , χ 2 12 , ... are independent χ 2 1 variates and λ l is the solution to the eigenvalue problem as in the right of Eq. (30) . Also, under the alternative H 1 , HSIC ω (Z) converges in distribution asEQUATIONwith the simplified notation E j,q,r := E zj ,zq,zr and E Z := E zi,zj ,zq,zr .Proof. The proof is mainly based on [34, Chapter 5] . A proof for a similar result has been given in [14, Theorem 1, 2] . The difference is that we consider the asymptotic distributions of HSIC ω (Z) that used learnable RFF thus is a function of frequency samplings while they consider HSIC b (Z).Thus some steps need to be modified.Step 1: we show that HSIC ω (Z) is a V-statistic, this can be done since it can be expressed asHSIC ω (Z) = 1 n 4 i,j,q,r h (ω)ijqr . To facilitate the conclusions in [34] for the U-statistic, we define the U-statistic HSIC (u) ω (Z) := 1(n)4 (i,j,q,r)∈i n 4 h (ω)ijqr that corresponds to HSIC ω (Z).Step 2: Then we prove the result under H 0 . To begin with, we first show that E j,q,r h(ω) ijqr = 0 with (i, j, q, r) ∈ i n 4 under H 0 . According to Eq. ( 20), we can calculate E j,q,r hijqr asEQUATIONwhere we define simplified notions E i̸ =u̸ =v̸ =w u,v,w whose superscript indicates the restriction. For readability, we will require additional notations:E x k (ω) i := E i̸ =u u k (ω) iu and E x k (ω) := E t̸ =u t,u k (ω) tu(the notations for Y is defined by analogy). Under H 0 , X and Y are independence. HenceEQUATION)Then combining the results [34, Section 5.5.2], we can prove Eq. (30) .Step 3: Next we prove the asymptotic distribution under H 1 . We only need to show that |HSIC ω (Z)-EQUATIONtu , we can check that |kEQUATIONtu | ≤ 1 for all t, u, thus |h (ω) ijqr | ≤ 4 for all i, j, q, r. Hence we haveEQUATIONCombining the results [34, Section 5.5.1], we can prove Eq. ( 31).In this section, we give a proof of the Theorem 1. We first restate the Theorem 1 here.Theorem 1 (Linear-Time Estimators). Under H 0 , the estimators of mean and variance with bias of O(n -1 ) to E Z [nHSIC ω (Z)] and Var Z [nHSIC ω (Z)], denote as E 0 and V 0 , respectively, are given byEQUATIONwhere () •2 is the entrywise matrix power. Both E 0 and V 0 can be calculated in O(nD 2 ) time.Proof. We first prove the part of the mean. Recall the definition of HSIC ω (Z), we have HSIC ω (Z) =1 n 4 i,j,q,r h (ω) ijqr . Hence E Z [HSIC ω (Z)] = 1 n 4 i,j,q,r E Z h (ω)ijqr . When (i, j, q, r) ∈ i n 4 , then we can show that under H 0 , E i,j,q,r h (ω) ijqr = 0 by performing the same analysis as in Eqs. ( 32) and (33) . Then we consider the case where exactly two elements of i, j, q, r are the same, for a total of 6n(n -1)(n -2) terms. By the symmetry of h(ω)ijqr , the expectation of these terms all take the same value, and here we take h (ω)iiqr as an example. According to Eq. ( 20), we can represent hEQUATION(ω)EQUATION(EQUATIONUnder H 0 , X and Y are independence. Take the expectation on both sides, we haveEQUATIONwhere we define additional notationE x k (ω) := E t̸ =u t,u k (ω)tu (the notation for Y is defined by analogy) and use k ω tt = l ω tt = 1. Hence in this case, the sum of the contributions of all terms toE Z [nHSIC ω (Z)] is (1 -E x k (ω) -E y l (ω) + E x k (ω) E y l (ω) ) + O(1/n).For the remaining terms, i.e., the case where at least three of i, j, q, r are equal, combined with the boundedness of h (ω) ijqr , we can conclude that the sum of their contributions is O(1/n). As a result, we have shown thatEQUATION)The unbiased estimators of E x k (ω) , E y l (ω) are given by 1T (Λ X Λ T X -I n )1, 1 T (Λ Y Λ T Y -I n )1, respectively.Hence under H 0 , we obtain the estimator of mean with bias of O(n -1 ) asEQUATION)Next, we prove the part of the variance. We start by calculating Var Z [nHSIC (u) ω (Z)], where the Ustatistic HSIC (u) ω (Z) := 1 (n)4 (i,j,q,r)∈i n EQUATIONwhereEQUATION2 . Under H 0 , when (i, j, q, r) ∈ i n 4 , we can show that E j,q,r h (ω) ijqr = 0 by performing the same analysis as in Eqs. ( 32) and (33) , thus ζ 1 = 0. For calculating ζ 2 , we mainly focus on the term E q,r h (ω) ijqr . We use Eq. ( 20) again,EQUATIONUnder H 0 , X and Y are independence. Take the expectation E q,r on both sides, we haveEQUATIONwhere we define additional notationE x k (ω) := E t̸ =u t,u k (ω) tu , E x k (ω) j := E j̸ =u u k (ω)ju (the notations for Y are defined by analogy). For simplify, we denote k(ω) c;ij := k (ω) ij -E x k (ω) i -E x k (ω) j + E x k (ω) and l (ω) c;ij := l (ω) ij -E y l (ω) i -E y l (ω) j + E y l (ω). Hence combining Eq. ( 40), we haveEQUATIONSince under H 0 , the bias lead by the difference terms between HSIC (u) ω (Z) and HSIC ω (Z) vanish faster than Eq. ( 43), hence the variance of HSIC ω (Z) is identical. In the following part, we consider the empirical estimate of the leading term in Eq. (43) . We estimate kEQUATIONwhere we define notionsΛ k := 1 n n u=1 Λ k (x u ). Since computing the value of 1 T (Λ Xc Λ T Xc ) •2 1 requires O(n 2) time complexity, we transform it into a more computationally tractable form. We perform the following calculating asEQUATION) Recall the definition of Λ Xc := [HΛ X ] n×D that can be calculated in O(nD) time, thus the term [Λ T Xc Λ Xc ] D×D can be calculated in O(nD + nD 2 ) time. As a result, we obtain the estimatorEQUATIONc;ij ) 2 that can be calculated in O(nD 2 ) time. The only thing left to do is to determine the bias of the estimator. For readable, we define k(ω) ij := k (ω) ij , k (ω) i := 1 n u k (ω) iu and k (ω) := 1 n 2 u,v k (ω)uv , then by removing the terms with i = j, a estimate with difference O(n -1 ) to Eq. ( 44) is given byEQUATION)By comparing the difference between the expectation of Eq. ( 46) and E i,j (k (ω) c;ij ) 2 , we can show that this error is bound by O(1/n). We illustrate this by taking one of the cross terms as an example and the other terms by analogy, as shown in the following,EQUATIONSimilarly, we can obtain the results for E i,j (l(ω) c;ij ) 2 .As a result, we have shown that V 0 is a estimator of Var Z [nHSIC ω (Z)] with bias O(n -1 ) thus complete the whole proof.F Calculation of Eq. ( 16)Here, we give the computational details of Eq. ( 16). We mark colors to indicate correspondences.According to Eq. ( 20), we can calculate j,q,r hEQUATIONWe can further represent Eq. ( 48) in matrices form asEQUATIONNext, by variable substitution, we obtain the result asj,q,r h (ω) ijqr = 1 2 n1 T A1 + n 2 (A1) i + (1 T C)B i + (1 T B)C i -nE i -nF i -nD i -1 T D .(50) where the definition of variables A to F with the calculation cost are given in the Fig. 1 . For convenience, we re-show the diagram here for reference. The computational complexity of each step is illustrated in Fig. 5 . We explain some steps here. As a start, recall that the size of Λ X , Λ Y are both n × D. Therefore a time complexity O(nD 2 ) is required to compute Λ T X Λ Y by matrix multiplication operation. Further multiplying the obtained[Λ T X Λ Y ] D×D with Λ T Y requires O(nD 2 ) time complexity. Next, since both Λ X and (Λ T X Λ Y Λ T Y )T are of size n × D, the elemental product operation requires a time complexity of O(nD). In a similar way, we can check the time complexity for each remained step. After getting the variables A to F, since 1 T A1, A1 all can be calculated in O(nD) and 1 T B, 1 T C, 1 T D all can be calculated in O(n), we conclude that the results with index i in Eq. ( 50) can be obtained in O(nD 2 ) time.In this section, we give a proof of the Theorem 2. We first restate the Theorem 2 here. Theorem 2 (Uniform Bound). Let θ k , θ l parameterize T θ k , T θ l in Banach spaces of dimension d k , d l . And T θ k , T θ l are Lipschitz to the parameters θ k , θ l with the non-negative constant L k , L l , respectively. Let Θ c be a set of (θ k , θ l ) for which σ ω ≥ c > 0 with a positive constant c and ∥θ k ∥ ≤ R θ k , ∥θ l ∥ ≤ R θ l . Let r denote the threshold, i.e., (1 -α)-quantile for the distribution in Eq. ( 11) and r (n) be the threshold with sample size n. Let {(ω k;j , ω l;j )} D/2 j=1 be the samplings of frequency with the sampling number D. Also, we define R ω k := sup j ∥ω k;j ∥, R ω l := sup j ∥ω l;j ∥, d s := max{d k , d l } and ξ ω := HSIC ω (Z). Then with probability at least 1 -δ, we havesup (θ k ,θ l )∈Θc ξ ω -r (n) ω /n σ ω - E Z ξ ω -r ω /n σ ω ∼ O 1 n log 1 δ + d s log n n + R ω k L k + R ω l L l √ n .Proof. We take a similar roadmap of proof as [30] and extend it to our optimization objective. The roadmap of the proof is as follows: we first obtain the convergence results (with sample size n) for each estimator with fixed parameters θ k , θ l , and then extend the results to the entire parameter space via ϵ-net arguments. We begin the proof of the first part, which is based on bounded differences inequality (McDiarmid's inequality) [41, Theorem 2.9.1].ω -E Z ξ ω |. Recall the definition of ξ ω := HSIC ω (Z) = 1 n 4 i,j,q,r h (ω) ijqr . By the definition of k (ω) tu , l (ω) tu , we have |k (ω) tu | ≤ 1, |l (ω) tu | ≤ 1 for all t, u, thus |h (ω)ijqr | ≤ 4 for all i, j, q, r. Now we begin by showing the bounded differences property of h (ω) ijqr . Concretely, we replace the first sample z 1 = (x 1 , y 1 ) with z ′ 1 = (x ′ 1 , y ′ 1 ) and keep the remaining samples the same. The obtained samples are named as Z ′ . Then the difference terms between h (ω) ijqr and the new substitution h(ω) ijqr can only happen in the case that at least one of i, j, q, r is equal to 1. For the case that only one subscript is 1 (here take i = 1 for example), combining Eq. ( 20), we haveEQUATION)The whole contributes of remaining terms that at least two i, j, q, r are less than O(n -2 ), thusEQUATION)Hence HSIC ω (Z) satisfy the bounded differences property with O(n -1 ). Using McDiarmid's inequality, for fixed θ 0 , θ 1 , with probability at least 1 -δ, there exist a universal constant C 1 such thatEQUATIONBound of |r(n) ω -r ω |. As r (n)ω is the (1 -α) of the distribution of nξ ω with sample size n under H 0 , according to the Eq. ( 53), when n is large enough, there exist a universal constant C 2 such thatEQUATIONwhere the last inequation is because E Z ξ ω ∼ O(n -1 ) under H 0 (see Theorem 1 for a detailed explanation). Hence |r(n) ω | ∼ O n log(1/α) . Also, by definition r ω is a constant related to α. Bound of | σ 2 ω -σ 2 ω |. In this part, We first obtain the bound of | σ 2 ω -E Z σ 2 ω |, then obtain the bound of |E Z σ 2 ω -σ 2 ω |.As before we start by showing the bounded variance property of σ 2 ω . We replacez 1 = (x 1 , y 1 ) with z ′ 1 = (x ′ 1 , y ′ 1 )and keep the remaining samples the same. The obtained samples are named as Z ′ . For readable, we denote σ 2 ω with sample Z, Z ′ as σ 2 ω (Z), σ 2 ω (Z ′ ) respectively.Recall the definition σ 2 ω := 16EQUATIONEQUATIONAgain, the difference terms between h(ω)ijqr and the new substitution h(ω) ijqr can only happen in the case that at least one of i, j, q, r is equal to 1. Hence, Eqs. ( 55) and ( 56) are both O(n -1 ). As a result, σ 2 ω satisfy the bounded differences property with bound O(n -1 ). Using McDiarmid's inequality, with probability at least 1 -δ, there exist a universal constantEQUATION)By adding further restrictions that i, i ′ , j, q, r, j ′ , q ′ , r ′ are all different, we can obtain the corresponding expression for σ 2 ω . Hence the difference between them can only happen when at least one subscript in i ′ , j, q, r, j ′ , q ′ , r ′ is equal to i. Combining |h(ω) ijqr | ≤ 4, we have |E Z σ 2 ω -σ 2 ω | ∼ O(n -1). ϵ-net arguments. Next, we prove the second part with ϵ-net arguments. Take the parameter space Θ k of θ k as an example. We choose a cover withN (Θ k , r k ) points {p i } N (Θ k ,r k ) i=1such that for any point p ∈ Θ k , we have min i ∥p -p i ∥ ≤ r k . According to [41, Proposition 4.2.12] , by comparing the volumes, we haveN (Θ k , r k ) ≤ (4R Θ k /r k ) d k .As for the parameter space Θ l of θ l , we can also obtain a cover withN (Θ l , r l ) points that N (Θ k , r k ) ≤ (4R Θ l /r l ) d l . Here, we set r k = 4R Θ k / √ n, r l = 4R Θ l / √ n, thus N (Θ k , r k ) ≤ ( √ n) d k , N (Θ l , r l ) ≤ ( √ n) d l .Then combining the Lipschitz property as shown in Lemma 4, we have with probability at least 1 -δ,sup (θ k ,θ l )∈Θc ξ ω -E Z ξ ω ≤ C 1 1 n log 2N (Θ k , r k )N (Θ l , r l ) δ + 8R ω k L k • r k + 8R ω l L l r l ≤ C 1 1 n log 2 δ + (d k + d l ) log n 2n + 32R ω k L k R Θ k √ n + 32R ω l L l R Θ l √ n .Hence when n is large enough, there exists a positive constant C 3 , with probability at least 1 -δ,EQUATIONSimilar, when n is large enough, there exists a positive constant C 4 , with probability at least 1 -δ,EQUATIONOverall Bound. Now we combine the previously obtained results. We haveξ ω -r (n) ω /n σ ω - E Z ξ ω -r ω /n σ ω ≤ ξ ω -E Z ξ ω σ ω + r (n) ω -r ω n σ ω + |E Z ξ ω -r ω /n| • 1 σ ω - 1 σ ω .Since on Θ c , σ ω ≥ c, according to Eq. ( 59), we can make σ ω ≥ c/2 happen by assigning probability budget δ/2 when n is large enough. Also, combining Eq. ( 54), |E Z ξ ω | ≤ 1 and r ω is a constant related to α, then with n large enough, there exist positive constants C 5 , C 6 ,EQUATION)Note that we need to pay for probability budget δ/2 for the above conclusion to hold. Then by taking the supremum on both sides in Eq. ( 60) and assigning the remained probability budget δ/2 to Eqs. ( 58) and (59), we can show that when n is large enough, with probability at least 1 -δ,sup (θ k ,θ l )∈Θc ξ ω -r (n) ω /n σ ω - E Z ξ ω -r ω /n σ ω ∼O 1 c 3 1 n log 1 δ + (d k + d l ) log n n + R ω k L k R Θ k + R ω l L l R Θ l √ n (61)thus complete the proof.In this section, we give a proof of the Theorem 3. We first restate the Theorem 3 here. Theorem 3 (Consistency). Let θ * k , θ * l be the parameters after learning, Z te be the testing samples of size m, when E Z HSIC (u) ω (Z te ) > 0, then the probability of Type II errorEQUATION)Let the mapping functions with learned parameters θ * k , θ * l be T θ * k , T θ * l , and the corresponding range space be a compact subset of R d Tx , R d Ty , respectively. Also, the diameter of two range spaces is denoted by diam(T θ * k ), diam(T θ * l ), respectively. Let {(ω k;j , ω l;j )} D/2 j=1 be the frequency samplings with their second moment denoted by σ 2ω k := E p k (ω) [ω T k;j ω k;j ], σ 2 ω l := E p l (ω) [ω T l;j ω l;j ]. Additionally, we denote ξ u := HSIC(X, Y ), then under H 1 , we have E Z HSIC (u) ω (Z te ) > 0 with any constant probability when D = Ωd Tx +d Ty ξ 2 u log σω k diam(T θ * k )+σω l diam(T θ * l ) ξu .Proof. The proof consists of two parts, we first give the rate of convergence of Type II error under condition E Z HSIC (u) ω (Z te ) > 0, and next for condition E Z HSIC (u) ω (Z te ) > 0 we give a lower bound on the number of frequency samplings required for it to hold. To simplify, we denote the U-statistic HSIC (u) ω (Z te ) as ξ(u)ω in this proof. We begin to prove the first part. With the learned parameters θ * k , θ * l , the probability of the Type II error is given byEQUATIONCombing the result of the difference between ξ ω and ξ (u) ω as shown in Eq. ( 34), we haveEQUATIONwhere C 0 is a positive constant. To apply the rate of convergence of the Central Limit Theorem, we rewrite the right equation in Eq. ( 64) asEQUATION)where the standard deviation (defined in Proposition 1) σ ω > 0 under H 1 . Then according to the results in [34, Section 5.5.1 Theorem B], there exist nonnegative constant C 1 such thatEQUATION)whereν h := E i̸ =j̸ =q̸ =r Z |h (ω) ijqr | 3 < ∞.When m is large enough, we further haveEQUATION)where C 2 , C 3 , C 4 are positive constants and using r (m) ∼ O(m 1/2 ) we prove in Eq. ( 54). Hence whenE Z ξ (u) ω > 0, the leading term √ mE Z ξ (u)ω decrease as m increase. Further, to obtain the decrease rate when m is close to infinity, we consider the asymptotic expansion (when x is close to negative infinity) for the function Φ(x) as given byEQUATION) thus Φ C 2 -C 3 √ mE Z ξ (u) ω + C 4 / √ m ∼ O(m -1/2). As a result, the decreasing rate is at leastO(m -1/2). We have so far completed the first part of the proof, and we next begin the second part of the proof, i.e. obtain the number of frequency samplings required for the condition E Z ξ (u) ω > 0 to hold. For simplify, we denote Λ k (x) T Λ k (x ′ ), Λ l (y) T Λ l (y ′ ) as k (ω) (x, x ′ ), l (ω) (y, y ′ ), respectively. Then according to Lemma 2, we haveP sup x,x ′ ∈X |k (ω) (x, x ′ ) -k(x, x ′ )| ≥ ϵ ≤ 2 8 σ ω k diam(T θ * k ) ϵ 2 exp - Dϵ 2 4(d Tx + 2) , P sup y,y ′ ∈Y |l (ω) (y, y ′ ) -l(y, y ′ )| ≥ ϵ ≤ 2 8 σ ω l diam(T θ * l ) ϵ 2 exp - Dϵ 2 4(d Ty + 2).(69) Also, we denote the bounds in Eq. ( 69) as δ x (ϵ, D), δ y (ϵ, D), respectively. Next we get the bound between E Z ξ (u)ω and HSIC(X, Y ). According to Lemma 3, the bound is given byEQUATION)Since by the definition, for allEQUATION)Combining the results of Eqs. ( 70) and ( 71), we obtainEQUATION)Combining the results as shown in Eq. ( 69) and allocating the probability budget ϵ, we haveEQUATION)By setting ϵ = ξ u /2, and sinceξ u > 0 under H 1 , we conclude that E Z ξ (u) ω > 0 holds with any constant probability when D = Ω d Tx +d Ty ξ 2 u log σω k diam(T θ * k )+σω l diam(T θ * l ) ξu .In the proof of Theorem 3, we obtain the convergence bound of the statistic E Z ξ (u)ω . Actually, the convergence result for its estimation can also be obtained, as shown in the following Corollary. Corollary 1 (Approximation Error Bound of HSIC ω (Z te )). Maintaining the same conditions and notions as in Theorem 3, we have the uniform convergence bound of HSIC ω (Z te ) asEQUATIONProof. By the definition, we have for all i, j, q, r, |kEQUATIONThen according to the results as shown in Eq. (69), we have for all i, j, q, rEQUATIONRecall the definition of hEQUATIONtv and we further define the corresponding h ijqr = 1 4! (i,j,q,r) (t,u,v,w) k tu l tu + k tu l vw -2k uv l tv , then for all i, j, q, r,EQUATIONAfter that, using the expressions that we obtained before, i.e., HSIC ω (Z te ) := 1n 4 i,j,q,r h (ω)ijqr and HSIC b (Z te ) := 1 n 4 i,j,q,r h ijqr , we obtain the final bound that P supEQUATIONand thus complete the proof.We first prove the Lipschitz property for some functions. For ease of reference, we re-list here the definitions of the terms that related to the optimization objective:ξ ω := HSIC ω (Z) = 1 n 4 i,j,q,r h (ω) ijqr , σ 2 ω := 16 1 n i ( 1 n 3 j,q,r h (ω) ijqr ) 2 -HSIC 2 ω (Z) and σ 2 ω := 16 E i (E j,q,r h (ω) ijqr ) 2 -E Z h (ω) ijqr 2 .The Lipschitz property of these terms are shown as follows.Lemma 4 (Lipschitz Property of of ξ ω , E Z ξ ω , σ 2 ω , σ 2 ω ). Maintaining the same conditions and notions as in Lemma 1, we have the following Lipschitz propertyEQUATIONwhere we use the symbol ξ ω (θ k , θ l ) to denote ξ ω with the parameter θ k , θ l and the others by analogy.Proof. We start by obtaining the result of h(ω)ijqr for all i, j, q, r. Since for all i, j, q, r,EQUATIONwhere the property |k(ω) ij | ≤ 1 and |l (ω)qr | ≤ 1 are used. According the Lemma 1, we haveEQUATIONCombing the definition that hEQUATIONAlso, combining the definition of ξ ω := HSIC ω (Z) = 1EQUATIONBy usingEQUATIONFor the results of σ 2 ω , σ 2 ω , we first proof the following results. For all i, j, q, r, i ′ , j ′ , q ′ , r ′ , we haveEQUATIONwhere the first inequality holds due to property |h(ω) ijqr | ≤ 4.Then we use the expressionEQUATIONand combine the results in Eq. (85). As a result, we obtain thatEQUATIONIn a similar way, we can obtain the corresponding expression of σ 2 ω asEQUATIONThen we obtain a similar result as before, i.e.,EQUATIONwhich completes the proof.Also for the term associated with the estimated threshold (recall that it is computed from the first two moments), we obtain the following properties of E 0 , V 0 as defined in Theorem 1. Lemma 5 (Lipschitz Property of of E 0 , V 0 ). Maintaining the same conditions and notions as in Lemma 4, we have the following Lipschitz propertyEQUATIONwhere the constant C 0 (n) := n 2 (n-1) 2 and C 1 (n) := n(n-4)(n-5) (n-1)(n-2)(n-3) .Proof. The expression in Theorem 1, while easy to compute, is not suitable for this part of our proof. We begin by obtaining equivalent expressions for E 0 , V 0 . According to Eq. ( 39), we haveEQUATION)And for V 0 , we use Eq. ( 45) and obtainEQUATIONAlso, we define hEQUATIONtv that corresponding to h (ω) ijqr and also define hEQUATIONtv . Then we can further rewrite the term for X asEQUATIONijqr and for Y by analogy. Then the properties of h (ω) ijqr can also be obtained for h(k) ijqr and h (l)ijqr , e.g., |hEQUATIONijqr | ≤ 4 for all i, j, q, r. Next we start to prove the Lipschitz property of E 0 , V 0 . For E 0 , according to Eq. ( 81) and combining the results 0 ≤ 1 T (Λ X Λ T X )1 ≤ n 2 and that for Y , we can show thatEQUATIONwhere1 T (Λ X Λ T X )1 = i,j k (ω)ij by definition and so as for Y . And for V 0 , we can prove thatEQUATIONwhere the last inequation is obtained similar to Eq. ( 85), thus completes the proof.The following results extend the results in [30] to the more general case (we only restrict the mapping functions to satisfy the Lipschitz property, and thus include the Gaussian kernel case of their proof). Theorem 4 (Smoothness of Optimization Objective). Let the sample of size n be Z, and with a small positive constant c, let the set of the parameters beΘ c := {(θ k , θ l )| σ ω ≥ c, V 0 ≥ c, E 0 ≥ c} , then there exist a nonnegative constant L such that ∥∇ (θ k ,θ l ) J∥ ≤ L on Θ c, where the optimization objective J := [HSIC ω (Z) -c α /n]/ σ ω is that we used in practice.Proof. According to Lemma 4, we have shown that HSIC ω (Z), σ ω both fits the Lipschitz condition. Also, according to Lemma 5, we have shown that E 0 , V 0 are also Lipschitz with respect to θ k , θ l . Since the threshold c α is completely determined by these two moments, combining the smoothness property of the mapping from two moments to thresholds as analyzed in [30, Theorem 2], we obtain that c α is also Lipschitz with respect to θ k , θ l on Θ c . As a result, we complete the entire proof based on the Lipschitz property of composite mappings.Remark. E 0 and V 0 are positive is almost satisfied in practice since according to the definition only[1 T Λ •2 Xc 1][1 T Λ •2 Y c 1] and [1 T (Λ T Xc Λ Xc ) •2 1][1 T (Λ T Y c Λ Y c ) •2 1]need to be greater than 0.In this section, we give an introduction to the comparison methods in our experiments and provide the implementation details of each method.The methods of comparison used in the experiment are described below.• dCor [39] : An independence test that is based on the distance covariance.• QHSIC [14] : The original quadratic-time HSIC independence test.• RDC [26] : The randomized dependence coefficient that measures the independence using the canonical correlation between a finite set of random features of the copula. • NyHSIC [44] : A variant of HSIC that uses the Nyström method to approximate kernels.• FHSIC [44] : A variant of HSIC that uses the random Fourier feature to approximate kernels.• BHSIC [44] : A variant of HSIC with the block-based statistic.• HSICAgg [32] : An aggregated kernel test with the incomplete statistic of HSIC.• NFSIC [18] : A test uses the normalized version of the finite set independence criterion and chooses features on a hold-out validation set to optimize a lower bound on the test power.Below are the GitHub URLs for each comparison method.• dCor: https://pypi.org/project/dcor. Threshold. For dCor, QHSIC, RDC, NyHSIC, FHSIC, and BHSIC, we permute the samples 100 times to simulate the null distribution and compute the threshold. The thresholds for the remaining methods are obtained by asymptotic null distribution, i.e., we set the test threshold to the (1 -α)quantile of χ 2 (J) for NFSIC and obtain the test threshold of LFHSIC-G/M by gamma approximation.Details of Setup. The number of random features for FHSIC, LFHSIC-G/M, the number of induced variables for NyHSIC, the block size for BHSIC and the number of sub-diagonals R for HSICAgg are all kept consistent as recommended in [44] for fair evaluation. Specifically, we set the number of random mappings in RDC to 20 to ensure compatibility with large-scale datasets. The test location parameter J of NFHSIC is set as default as 10, since it differs from other approximation methods that blindly increasing J may lead to a loss of power as shown in [18] and can significantly escalate time costs due to its cubic time complexity O(J 3 ). In the optimization step, for stabilizing the training, in the implementation of NFSIC we determine the initial bandwidth by searching the best from 25 bandwidth combinations (including the median bandwidth combination). For LFHSIC-G/M, to be fair, we perform the same grid search on SD, Sin, and GSign datasets. In other experiments, we still use the median bandwidth as initialization for LFHSIC-G/M. Also, the maximum number of iterations for the optimization is set to 100 for NFSIC and LFHSIC-G/M. The default learning rate of the optimization of LFHSIC-G/M is set as 0.05 in all the experiments. As for HSICAgg, the default implementation of the predefined 25 pairs of bandwidths in its code is used. For synthetic data, we set the split ratio to 0.5 for NFSIC and LFHSIC-G/M, i.e., we randomly sample half of the data for training and use the remaining for independence testing, while the other methods use all data for testing. For real MSD data, we divide a small portion of the data for training and then extract 100 random subsets of the remaining data (disjoint from the training set) for evaluation.The details of the four synthetic datasets and two real datasets are described below.• Sine Dependency (SD): In this model, X follows a d-dimension multivariate normal distribution N d (0, I d ), and Y is defined as 20 sin 4π(X 2 1 + X 2 2 ) +Z, where X i is the i-th dimension of X, and Z ∼ N (0, 1) represents independent noise. Notably, when d > 2, Y exhibits a nonlinear relationship solely with the first two dimensions of X. , where sgn(•) represents the sign function, X i denotes the i-th dimension of X, and Z ∼ N (0, 1) is independent of X. The challenge lies in Y being independent of any proper subset of X but dependent on X as a whole, underscoring the importance of considering all dimensions of X simultaneously in independence testing.• ISA Dataset. We construct the data through the following steps: First, we generate n i.i.d samples of two univariate random variables with a mixture of Gaussian model, i.e. 1 2 N (-1, 0.01) + 1 2 N (1, 0.01). Second, we mix these random variables using a rotation matrix parameterized by an angle θ, which varies from 0 to π/4. A zero angle implies independence between the data, while a larger angle signifies stronger dependency. Third, we append noise with a distribution of N d-1 (0, I d-1 ) to each of the mixtures. Finally, we multiply an independent random ddimensional orthogonal matrix to obtain vectors dependent across all observed dimensions. The resulting random variables X and Y are dependent but uncorrelated. When d is greater than 1, the problem is associated with the independent subspace analysis (ISA) problem [14] . For the case d = 1, θ = π/10, an example visualization is shown in the middle of Fig. 6 . • 3DShapes Dataset. This dataset [5] comprises images depicting 3D scenes, complete with additional features like shadows and backgrounds. It encompasses six fundamental latent factors: floor hue, wall hue, object hue, object scale, object shape, and orientation, all adjustable to generate corresponding images. Orientation is treated as a dependency factor for independence testing, where we test the dependency between the image X and its orientation Y . To heighten the challenge, we maintain the object shape as a ball, minimizing the apparent orientation feature compared to other shapes like squares, while randomizing the remaining factors. An example visualization is given on the right side of Fig. 6 . • Million Song Dataset. The dataset, a subset of the Million Song Data3 [4] , comprises 515, 345 songs with 91-dimensional features. The first dimension represents the release year of each song, designated as variable Y , while the remaining features (e.g., timbre average and timbre covariance) form variable X. Our objective is to identify the dependency between X and Y . In this section, we provide additional experimental results, mainly including the visualization results on the Sin synthetic dataset, the results with more comparing methods (as explained in the main paper, due to the high time overhead therefore do participate in the evaluation of the main paper) as well as the running time of each method.We provide the visualization results on the Sin synthetic dataset to illustrate the performance of our optimization objective for the example mentioned in the main paper. The results are shown in Fig. 7 , where the setup follows our experiments (n = 2000, ω = 5, D = 100). For visualization, the negative of our optimization objective J is shown. As can be seen, our optimization objective guides to letting the bandwidth adapt to improve the test power, and here we can see that regions with bandwidths around 0.2 (corresponding to the theoretical optimal solution in our main paper) indicate better power, thus corroborating the validity of our optimization objective. Also, notice that the landscape is smooth over a wide range as demonstrated in Theorem 4, and thus contributes to non-convex optimization. As mentioned in the main paper, the method4 [30] is not involved in the comparison because it takes too much time to run in some settings. Here, we compare it with ours under some feasible settings to illustrate the improvement of our method on power-runtime trade-offs. The methods using Gaussian kernel with global width and Gaussian kernel with widths of each dimension (corresponding with ours) are employed, referred to QHSIC-O and QHSIC-W, respectively. For fairness, the same grid search procedure is employed as the initialization of the optimization. We perform the evaluation on the SD data and plot the results of the test power over time as shown in Fig. 8 . Also, for our methods, we provide the results under the setting D = 100 and D = 500 as in the main paper. The experiments are conducted with the same equipment, specifically a 6-core CPU with a 3080 GPU.Results. Our test consistently results in a better power-runtime tradeoff at different D settings. At D = 100, one test can be completed in less than a second when the sample size n = 6000. As D increases (D = 500), the number of samples required to achieve the same power decreases, but the increase in D leads to an overhead in runtime, and overall our test is still completed in a few seconds. In contrast, even though QHSIC-O/W requires fewer samples to reach the same power than our tests, the runtime rises rapidly as the sample size increases. When n = 1000, it already takes more than 10s to perform a test. When n = 3000, it needs nearly a minute to perform a test, which may greatly limit the practical application. Results. The experimental results on the left show that our method is faster than other methods with optimizable options (HSICAgg and NFSIC), and can complete a test within 10 seconds even with 100, 000 samples. Even though HSICAgg uses parallelism to optimize the computational efficiency of the scheme, the actual implementation of the parallelism is time-consuming, and hence leads to a high time overhead for a single practical test. For the results on the right, it can be seen that the quadratic complexity methods face a dramatic increase in time overhead as the sample size rises, and this is especially severe for the methods (QHSIC-O/W) that need to be optimized since the optimizing objective needs to perform multiple squared complexity operations. In contrast, our linear-time learning objective allows us to handle huge data samples very efficiently.Limitations. According to the experimental results in the main paper as well as in the Appendix, no one method is better than the others in all settings, so it is important to choose several appropriate tests for real scenarios and summarize their results in order to obtain a more reliable conclusion.Broader Impacts. This work proposes a novel framework for independence testing. The proposed linear-time optimization objective can be trained end-to-end in a data-driven manner, ensuring both effectiveness and efficiency in high-dimensional and large-scale scenarios. This could be beneficial for developing more reliable downstream algorithms in a variety of areas, including causal discovery, feature selection, and deep learning.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.The bounded, continuous, translation-invariant kernel satisfies the characteristic condition[12].Million Song Data subset: https://archive.ics.uci.edu/dataset/203/yearpredictionmsdThe code is downloaded from https://github.com/renyixin666/HSIC-LK 


Decision: This paper takes the Fourier perspective on the popular HSIC criteria, leading to a linear time computable novel criteria. All the reviewers agree that the work is quite novel. It will be great if the authors handle the apprehensions raised by the reviewers appropriately in their final version, if the paper is eventually accepted. I feel the linear computation aspect is perhaps what distinguishes the work from the others and needs to be highlighted appropriately.

Reviews: review: paper_summary: paper_summary: The paper presents a novel consistent estimator for HSIC, which is computationally efficient, and tries to maximize the testing power under controlled type-1 error. The idea is to begin with a known relation between HSIC and Fourier-based distance between the relevant characteristic functions. The estimator can be understood as a sample based approximation of this. For standard kernels, it is noted that the distribution from which the samples need to be generated is from a Gaussian. In general, it is assumed that a transformation of variables converts this distribution as a standard Gaussian. It is shown that this estimator has linear computational complexity, for given transformation parameters.

Further, it is proposed to learn the transformation's parameters (which in turn is same as learning the kernel parameters). And the objective for learning is proposed to be maximizing test power under controlled t1error. To this end, using asymptotic convergence arguments, an expression for this objective is analytically arrived at (prop1). It is discussed how to estimate this objective again in linear time (theorem1). Consistency of this objective and the overall estimator are are proved in theorems 2,3.

Simulations on synthetic and benchmarks datasets show the improved test power and computational complexity when compared with existing baselines.
strength_weakness: strength_weakness: 1. It would have been nice if the introduction section had more description of the methodology at an intuitive level. Currently, most of the introduction is motivation, related work and the actual methodology is restricted to one small para. This would give the reader a intuitive understanding and help to understand the later sections easily. Without this it took me some effort to understand what is happening in the overall methodology. (I may be still missing some important points?)

2. The derivations and the main algorithm in section 4.2,4.3 have a very close resemblance with [30] https://proceedings.mlr.press/v238/ren24a/ren24a.pdf . This significantly weakens the contribution and novelty. More importantly, main details in this section could have been cited from [30] or postponed to Appendix. The current presentation seems to give an impression that the sections are entirely new.

3. Since [30] is a close methodology, empirical comparison with it seems to be crucial. However this seems missing.
questions: questions: few minor comments:

1. notations at places is a bit confusing: e.g., line 121, 146 w_x means different etc.


--- after reading the author responses and other reviews I increase my score.
limitations: limitations: yes
paper_summary: paper_summary: The paper proposes a novel method of estimating HSIC to determine the independence of two random variables. The original formulation HSIC_b is reformulated and approximated using Monte Carlo integration, HSIC_w, with frequency samples. This reformulation was not derived by the authors but was borrowed from previous literature. The authors’ contribution is the sampling of frequencies proportionally to the inverse Fourier transform of kernel functions.
strength_weakness: strength_weakness: The proposed estimation is novel, but there are several downsides of the paper. 

Many explanations are simply provided without proper self-contained information. First, the purpose of constructing criterion J is unclear. The introduction of $c_\alpha$ is unclear. Second, the definition of Type I error for evaluation is not provided. Type I error is for the false positives. What is the definition of positive decisions in the experiment? What is the test power in the experiment? Those concepts are used without definition.

The authors used neural networks for T_\theta x. Conventional HSIC can also use similar transformation of x. For example, k(T_\theta x, T_\theta x’) can be calculated directly with a learned transformation T_\theta x. The results do not necessarily support the superiority of the calculation in the frequency domain.

In the experiments, I expected to see the consistency of the proposed estimator with large D, which is not provided.

In theory, the convergence is derived for E[HSIC_w], not for HSIC_b. 

Minor comments:
In the last line of Algorithm 1, HSIC_b should be HSIC_w.
questions: questions: Could you provide a detailed description of the objective J, Type I error, and test power used in the paper?

Please explain if we can use k(Tx_1, Tx_2) instead of using the proposed method.
limitations: limitations: The authors presented an interesting algorithm, but there are several notions that are unclear.
paper_summary: paper_summary: The authors propose to extend RFF with parametric transforms of the Fourier features of shift invariant kernels. A path similar to the seminal work of Gretton et al 2007 is followed, where the proposed $\text{HSIC}_\omega$ statistic is shown to converge in distribution to a Gaussian with the true mean (under the null). Then, linear time estimates of the first two moments are provided, which are needed in vanilla HSIC to analytically compute an approximation to the critical region.
strength_weakness: strength_weakness: The main weakness of this method is its limitation to translation invariant kernels, which is, however, quite a usual case.

Another weakness is that the dataset must be split in two parts to avoid the selective inference problem.

Finally, as an extension of HSIC, this method most likely also suffers from some hindrances in the latter, for instance that the Null distribution (or, in specific, the necessary quantile) has to be empirically approximated, say using bootstrapping, or a possibly looser bound needs to be used.

Some comments:
* Table 1 needs some more space around it to avoid confusion.
* ln 51: statistics
* $\xi_\omega$ needs to be defined in  Theorem 2, and for completeness also $\sigma_\omega$, $\hat\sigma_\omega$.
Perhaps these could be added in the Appendix as a per-lemma symbol table?
questions: questions: I had difficulty following the implications of Theorem 2, given its (reasonably) over-packed formulation and I believe that a couple more sentences explaining the involved quantities would make it easier to read.
limitations: limitations: This work is limited to translation invariant kernels which can be expressed as an integral using the spectral measure, as well as the dataset splitting issue.

A brief discussion could improve the proposed work by making these issues explicit.

Comments: title: Would you please check our rebuttal! Thanks a lot!
comment: Dear reviewer,

We submitted our rebuttal to your comments a few days ago. Would you please check our feedback and retrun your further comments? 

We are looking forward to your futher feedback on any issues of our paper! 

Thanks  again!
Dear Reviewer,

Thank you very much for your prompt feedback on our rebuttal. You asked us whether you missed something. Frankly, but without any intention to offense you, we think that you did miss something --- the most important things about our submission, that are, 1) contribution: we propose a new method for Statistical Independence Testing with $O(n)$ time and space complexities, while the method in [30] has time and space complexities of $O(n^2)$. This is absolutely a significant contribution to the area. 2) Novelty: to the end, we solve the problem from the frequency domain perspective, which is a novel solution, completely different from that of [30]. 

Our paper and [30] are all about independence testing by learning kernels, that is, they all try to solve the same problem but with different methods, specifically, we tried to develop a more efficient method with linear complexity. It is very normal that our paper and [30] use similar notation systems, have a similar background introduction or preliminaries, and adopt a look-alike derivation process. In our point of view, the value of a paper depends on its contribution and novelty, not its writing style and derivation process. 

If possible, we would appreciate it if you could spend a little more time to have a more careful comparison between our paper and [30], and we are sure you will see the significant difference between the two works.

We are looking forward to further discussing with you any issues about our paper. Thanks again!
Dear authors, I agree that the motivation is different from [30]. Also, I   appreciate the linear time complexity and the Fourier based idea. However my main complaint is that the derivations in sections 4.2,4.3 more or less carry forward from [30]. Infact, most of the material , it seems, can be cited from [30] and details postponed to appendix. Am I missing something ?