PHSIC against Random Consistency and Its Application in Causal Inference

1 Introduction

Association analysis is the basis of data mining and machine learning, which has been applied in many fields. Examples include feature selection [Li et al., 2021; Song et al., 2007b], gene analysis [Yamanishi et al., 2004], causal inference [Hoyer et al., 2008], etc. There are various forms of correlation measurement, which can be roughly classified into two categories: the construction method based on the reduction error ratio [Puth et al., 2015; Hotelling, 1992] and the construction method based on the independence test [Jaworski et al., 2010; Schweizer and Wolff, 1981]. The typical representative of the first construction method is the person correlation coefficient, but it can only recognize simple linear relationships. The typical representatives of the second construction method are mutual information and HSIC. Although both mutual information and HSIC can identify different forms of dependency relationships, the empirical estimation of mutual information is expensive, while the empirical estimate of HSIC is much simpler-just the trace of the product of gram matrices.

∗The corresponding author

The good nature of HSIC drives the development of many fields. For example, in terms of neural network optimization, since the maximum information of continuous variables is difficult to calculate, Ma et al. [Ma et al., 2020] used the Hilbert-Schmidt information bottleneck (HSIC Information Bottleneck) as the regular term in the neural network optimization objective to improve the adversarial robustness of the model. In clustering problems, HSIC is introduced many times to measure the independence of cluster centers [Niu et al., 2010; Song et al., 2007a]. In the field of feature selection [Liaghat and Mansoori, 2019; Song et al., 2007b], HSIC is used to calculate the degree of association between features and labels. In the field of causal inference [Hoyer et al., 2008], HSIC is used to measure the independence between independent variables and residuals in additive noise models, so as to judge the causal direction. These are attributed to its three good characteristics: (1) it has smaller bias than other indicators; (2) it can capture many complex nonlinear dependencies without explicitly considering random variables ; (3) its empirical estimate is the trace of the product of grammatrices, which makes the objective function easy to solve.

With sufficient sample size, HSIC undoubtedly shows good performance. However, in some application scenarios, data collection is difficult due to some complex factors, so the collected samples are relatively sparse. At the same time, this random error will be amplified in high-dimensional and noisy data. HSIC generates random errors in limited samples, which causes a series of problems. (1) Firstly, as shown in Figure 1(a), in an independent situation, there is no 0 baseline, which lacks interpretability. And we can see that the smaller the sample, the higher the score of HSIC, and there is sample bias. (2) Secondly, as shown in Figure 1(b), due to the bias of limited samples, random consistency can also occur when comparing different correlations, which can seriously affect the performance of the application task. (3) Finally, as shown in Figure 1(c), in independent situations, as the dimension gets higher, the value of HSIC gets larger, which leads to the dimension deviation problem.

Not only does HSIC exhibit random consistency, but random consistency is a common problem in statistics and machine learning. In classification tasks, due to the limited number of samples and the presence of noise, the results generate random consistency, which makes machine learning algorithms lack generalization ability. In order to alleviate the random consistency in classification tasks, Wang et al. [Wang et al., 2022] proposed a pure consistency index, which verified the low bias, learning substitutability, and high recognizability of pure accuracy frameworks compared to accuracy. And Romano et al. [Romano et al., 2016] improved the classification performance of decision trees by eliminating random factors in the Gini index. In clustering tasks, the results may have biases due to factors such as the number of clusters and sample size. Vinh et al.[Vinh et al., 2009] emphasized the importance of eliminating mutual information random consistency in clustering communities, especially when the data size is relatively small relative to the number of clusters.

Based on the above observations, we propose a pure HSIC framework that can alleviate the drawbacks of HSIC and ensure that the original HSIC properties remain unchanged. PHSIC is obtained by separating the mean of random parts from the original HSIC values. To emphasize the advantages of PHSIC, we have demonstrated theoretically and experimentally that PHSIC has smaller deviations compared to HSIC, especially in independent cases where PHSIC has an error order of m-2 and HSIC has an error order of m-1 , where m represents the number of samples. At the same time, it has also been proven that the random part is the reason for the deviation of HSIC in dimensions, samples, and noise.

Association is an important issue in causal discovery. Current causal inference methods can generally be summarized as 'association + hypothesis'. For example, in constraintbased causal inference methods [Spirtes et al., 2000; Ogarrio et al., 2016], these methods mine the causal relationship between variables through the conditional independence test under the Markov faithful assumption. In function-based causal inference methods [Shimizu et al., 2006; Zhang and Hyvarinen, 2012], under the assumption of independent causal mechanism, the causal direction is further inferred by measuring the correlation between residuals and causes through correlation indicators. It can be seen that a robust correlation indicator is the premise to promote the performance of causal inference methods. Currently, identifying causal directions is difficult in small sample and noisy scenarios. Therefore, we applied PHSIC to the classical causal models to further improve the performance of causal inference model. The main code and supplementary material have been made available online 1 . The contributions of this paper are as follows:

• We analyze the random consistency phenomenon of HSIC from both theoretical and experimental perspectives.
• We propose a PHSIC framework to eliminate random consistency, proving theoretically and experimentally that PHSIC has smaller deviations than HSIC, and discussing the properties of PHSIC.
• We apply PHSIC to multiple classic causal inference tasks and experimentally verify the robustness of PHSIC based causal models to small sample and noisy data.

2 Hilbert-Schmidt Independence Criterion

In this section, we briefly introduce HSIC and demonstrate the bias of HSICb(Z).

Suppose X and Y are two random variables. Their values belong to χ and γ respectively. F presents reproducing Hilbert space for each x ∈ X. We define that φX(x) ∈ F as an element in a reproducing Hilbert space corresponding to each point x ∈ χ . Assuming the existence of a continuously bounded positive definite kernel kX : χ2 → R, where

For example, the Gaussian kernel k(x1, x2) = exp(-||x1 - x2||2 /σ2 ) is the primary choice in many literature. It is a continuously bounded kernel function, and we mainly default it to the selected kernel function in this article. Similarly, for a random variable Y , let G be a reproducing kernel Hilbert space with kernel kY : γ 2 → R and feature mapping φY : γ → G.

In order to observe the correlation between X and Y , Gretton et al. [Gretton et al., 2005] defined cross covariance operator,

CXY means a linear operation from G → F. And ⊗ represents the tensor product, µX = EX(φX), µY = EY (φY ) represents the average value of the corresponding elements of random variables X and Y in a reproducing Hilbert space, respectively. This equation cannot quantify correlation, so the square of Hilbert normal form is used instead of the cross covariance operator and it can be expressed in terms of kernels:

1 https://github.com/lijue688/main.git

which is a measure of the statistical dependence of X and Y . When random variables X and Y are independent then HSIC(F, G, PX,Y ) = 0. If X and Y are correlated, then HSIC(F, G, PX,Y ) > 0. And as the values of X and Y become more correlated, the value of HSIC becomes larger.

Gretton et al. [Gretton et al., 2007] define the following empirical HSIC estimator for an i.i.d. sample Z = {(xi , yi )}i=1,...,m:

where m is the number of data points, K is the m×m kernel matrix for X and L is that for Y :

and H is the m × m matrix defined by:

where I is an identity matrix of m×m and I is a vector of m × 1.

Gretton et al. show that for m → ∞, HSICb(Z) → 0 if and only if X is independent of Y [Gretton et al., 2007]. Next, we explored the inherent bias of HSICb(Z).

Theorem 1 ( Bias of HSICb(Z)). Let E(Z) denote the expectation over m independent copies(xi, yi) drawn from Pxy, then we have

We provide the proof of this theorem in the appendix. It can be observed that although the empirical estimates of HSIC proposed by Gretton et al. [Gretton et al., 2005; Gretton et al., 2007] in 2005 and 2007 have different coefficients, one is 1/(m- 1)2 , and the other is 1/(m)2 , their deviation orders are both O(m-1 ). From its deviation order, it can be seen that when the sample size is relatively small, it will have a certain impact on the performance of HSIC.

3 Pure HSIC to Alleviate Random Consistency

In order to alleviate the random factors generated by HSICb, we first propose a framework for eliminating random consistency, and then explore the properties of PHSIC.

Definition 1 (Eliminating Random Consistency Framework). Due to the degree of consistency between some random variables being caused by their random distribution, we introduce a framework PCM(Z1,Z2) for eliminating random consistency:

where Z1, Z2 are random variables. CM(Z1, Z2) represents the degree of consistency between Z1 and Z2, which usually represents the degree of correlation or similarity between two variables. And RCM(Z1, Z2) denotes the degree of random consistency between Z1 and Z2, which causes by factors such as limited number of samples, high dimensionas, and uniform noise. Assuming a set Z ′1 is a uniform random variable with the same value and distribution as Z1, the expected values of random consistency in random variables Z ′1 and Z2 can be used to calculate RCM(Z1,Z2). Therefore, PCM will characterize pure consistency by the difference between the degree of consistency and the degree of random consistency. Definition 2 (Empirical PHSIC). With the framework of equation (8), we define PHSIC:

where E(HSIC0(Z, F, G)) represents its expected value under the null distribution. The analysis formula is shown in equation (10), which can obtain from [Gretton et al., 2007].

Theorem 2 (Bias of PHSIC(Z)). The bias of empirical PHSIC(Z) is as follows:

Proof. Gretton et al. proposes that the unbiased estimator of HSIC can be replaced by three u-statistics [Gretton et al., 2007]:

set imn indicate that there are no duplicate combinations, the empirical averages can be used instead of expected values:

Same idea as Gretton et al. [Gretton et al., 2007], we can obtain deviation between biased estimation of HSICb(Z) and unbiased estimation of HSIC(Z) through the three difference terms, as shown in equation (14), (15), and (16):

Combining these three items, then seeking expectations on both sides, the deviation of HSICb(Z) is:

In the independent case,-3Exx′ yy′kl + 10Exx′ yy"kl - 6Exx′kEyy′ l can be merged into one item. Ignoring the bias of O(m-2 ), [Gretton et al., 2007] derived the mean of HSIC under the null hypothesis, as shown in equation (18):

Therefore, based on equations(13), (17) and (18), we can conclude:

Corollary 1. The deviation of PHSIC(Z, F, G) in independent cases is O(m-2 ).

Proof. Theorem 2 has already mentioned that in the independent case, -3Exx′ yy′ kl + 10Exx′ yy"kl - 6Exx′ kEyy′ l can be merged into one item. Therefore, the bias of HSICb in Equation (17) can be written as:

Therefore, the deviation of PHSIC (Z) in independent cases as:

Corollary 2. The deviation of PHSIC(Z,F, G) in dependent cases is O(m-1 ).

Proof. In dependence cases, the deviation of PHSIC is shown in equation (19), we can see that

3.1 Advantages of PHSICCompared to HSIC

In this section, we mainly analyze the advantages of PHSIC compared to HSICb in independent situations.

PHSIC has zero baseline. In the independent case, we can see that the mean of PHSIC is as follows:

And the average value of HSICb under independent conditions is shown in equation (7), its mean is not 0. Therefore, PHSIC is more interpretable.

PHSIC is closer to unbiased in sample size. From Theorem 1, we deduce that HSICb has a deviation order of O(m-1 ) in the independent case, while PHSIC has a deviation order of O(m-2 ) in the independent case (as shown in equation (21)). Therefore, the latter has a smaller error order than the former and is closer to unbiased.

PHSIC may closer to unbiased in dimensions. From equation (21), it can be seen that the deviation term of O(m-1 ) in HSICb comes from E(HSIC0). We observe the dimensional variation of this term from Monte Carlo perspective. Assuming that Xm×d , Y m×d is independent and uniform multidimensional random variables, d represents dimension. X = [x(1), x(2)..., x(d)]. The values of each dimension of X are the same, but the positions are different. Y = X(π), π indicates disrupting the order of column X. In this way, X and Y can remain independent, and the values of the number of columns are the same, but the positions are different. for ||uxd||2 = (m)-1 2 ∑(i,j)∈im2 ki,j, when m is fixed, the individual term kij decreases as the dimension d increases, since the value of ||uxd||2 reduce. Therefore, for E(HSIC0(Zn|X, Y )) = 1m (1 + ||ux||2 ||uy||2 - ||ux||2 -||uy||2 ), (m)-1 2 ∑(i,j)∈im2 ki,j ∈ (0, 1]. We think of it as a function of ||ux||2 (because the values of ||ux||2 and ||uy||2 are the same). Let ||ux||2 = a, E(HSIC0(Z)) = b is a function of a, b = 1 + a2 - 2a, derive from it, b′ = 2(a-1), a ∈ (0, 1], b is a monotonically decreasing function, When a decreases, b is monotonically increasing. Therefore, the E(HSIC0(Z) is monotonically increasing as the variable dimension increases. PHSIC subtracts E(HSIC0) from HSICb , making the dimension closer to unbiased.

According to the above analysis, from the perspective of random consistency, PHSIC has a smaller deviation compared to HSIC and has stronger properties.

3.2 Toy Example

We establish toy examples to further understand the advantages of PHSIC. We compare PHSIC with pHSIC (p-value of HSICb) and a popular depolarization framework SHSIC. Equation (24) is the calculation formula for SHSIC, a is constant with the value of 0.5.

Figure 2 shows the scores of pHSIC, HSICb, SHSIC and PHSIC under different noise levels in the nonlinear function y = 4(x - 5)2 in the case of binary univariate and random kernel parameters. It can be found from the four figures that compared with the other three indicators, PHSIC not only has zero baseline, but also converges quickly with a small sample. We find that p-values are only effective under global distribution, and in some smaller samples, p-values will be higher under correlated conditions, which leads to a type II error and judges some correlated values as independent. From Figure 3, we observe the sample selection fairness of PHSIC in a more intuitive way. We generated 5 groups of samples with different quantities to simulate the correlation between variables. The estimated values of Sn for different samples are calculated to calculate their selection probabilities. We select samples according to HSICb, PHSIC and SHSIC, and observe the fairness of these estimators in selecting samples of different quantities. For Figure 3 (a) and (b), we set Sn=[120 140 160 180 200] and simulated the correlations of x and y under correlated and independent conditions. Observing the selection probabilities of SHSIC, PHSIC, and HSICb for these five groups of samples, we can clearly see that PHSIC has more selection fairness.

In order to verify the unbiased nature of dimensions, we conduct experiments on independent uniform distribution simulation data with a sample size of 300 and a value range of 11-12. We increase the dimensions of x and y from 1 to 5 and compare the scores generated by different dimensions. From Figure 3 (c), it can be seen that PHSIC is more unbiased compared to HSICb and SHSIC, and its selection probability in different dimensions is closer to 0.2, HSICb and SHSIC tend to be more high-dimensional, which proves that PHSIC has more fair selectivity in different dimensions.

To further verify the fairness of PHSIC, we evaluate it under different samples and functional relationships. As shown by the red line in Figure 4 (a), we use HSICb to measure the correlation metric scores of variables x and y simulated by the function y = x2 under different noise levels for 60 samples. The green line represents the correlation metric scores of variables x and y simulated by the function y = 4(x - 5)2 under different noise levels when HSICb measures 40 samples. Figure 4 (b) using PHSIC measurement as a control. From the functional relationship, it can be seen that the scores of the variables x and y simulated by y = x2 are always higher than those simulated by y = 4(x - 5)2 in non independent cases, and the values should be the same in independent cases. From Figure 4(a), it can be seen that the independent values of HSICb measures are higher than those in some non independent cases, leading to false associations. From Figure 4 (b), it can be seen that PHSIC effectively solves this situation.

In the above toy example, we mainly demonstrate the fairness of PHSIC in different and identical associations under different sample sizes. Next, we further demonstrate its robustness in the same sample and different association relationships. We use HSICb and PHSIC to measure the scores of variables x and y in independent and non independent cases under the same sample size. We simulate 10000 times and observe the probability that the scores of two indicators were higher in independent cases than in non independent cases. From Figure 4 (c), we can see that HSICb are more prone to misselection in the case of limited samples, further demonstrating the good performance of PHSIC.

4 PHSIC-Based Causal Inference Models

In this section, we first introduce several typical causal inference models, then describe the data and experimental details we used for evaluation, and finally present the experimental results.

4.1 Introduce Typical Causal Inference Models

In this paper, we mainly focus on four classic causal inference models: ANM [Hoyer et al., 2008], ANM-MM[Hu et al., 2018], KIKO [Assaad et al., 2019], HANM [Zhao et al., 2023].

• ANM-MMmodel. The additive noise model is extended to a mixed model consisting of a finite number of ANMs, and its causal identifiability conditions are given.
• KIKO model. In the context of ANM model, KIKO model uses one regression variable instead of two to accelerate causal inference.
• HANM model. Identify many to one causal relationships and use asymmetric forward and backward models HANM to identify causal direction.
• ANM model. Inferring causal directions through the asymmetry of forward and backward models using nonlinear functions.

All these models are based on HSIC to measure the independence of causes and residuals, which is a key part of these methods. Therefore, next we will use PHSIC to further improve the performance of these models.

4.2 Experimental Setup

In this section, we introduce the specific experimental details of the four models in the appeal:

ANM-MM model. we use eight real datasets: (Rings, Viscera), (mpg, Weight) and (age, weight), etc. In order to highlight the advantages of our PHSIC, we conducted experiments on 50 samples with the same experimental parameters as the original ANM-MM.

KIKO model. we use eight real datasets:(Temperature, Co2 flux),(GNI, life expectancy) and (NEP, PPFDdif), etc. And we conducted experiments on 30 samples with the same experimental parameters as the original KIKO model.

HANM model. we use four simulated data and four real data. For the first four data (y = ex , y = x, y = sin(x), concrete), we conducted experiments on 6, 7, and 8 causes respectively, and for the last four data (ozon, y = x2 , mpgacc, mpg-mpg), we use 2 and 3 causes respectively. Exponential data(y = ex ) and concrete data were tested on 50 samples, while others were tested on 30 samples.

ANM model. We use one real dataset(CEP) and three simulated datasets(SIM, SIM-G, SIM-IN). We compare the methods based on HSIC score (ANM-HSIC), Entropy score (ANM-Entropy), Gaussian score (ANM-Gauss), and empirical Bayesian score (ANM-FN) with our proposed PHSIC score (ANM-PHSIC). The introduction of the above comparison methods can be found in [Mooij et al., 2016].

A more detailed introduction to the above dataset and experimental parameters will be provided in the appendix. For HANM and ANM models, we use accuracy indicators for evaluation, and the calculation method is as follows:

where wm is the weight of each pair of data, because there is a strong correlation between variables in the real dataset, CEP dataset decouples each pair of variables by assigning different weights. All pairs from the same dataset are weighted equally and sum to 1. dm and ˆ dm are the real direction and predicted direction of the mth pair, respectively, if the two directions are the same, it is recorded as 1, otherwise it is 0.

4.3 Experimental Result

Table 1 shows the experimental results of the HSIC based ANM-MM model and the PHSIC based ANM-MM on eight datasets. Compared with GROUND TRUTE, it can be found that the ANM-MM model based on PHSIC can be correctly recognized on all eihgt datasets, while the ANM-MM model based on HSIC cannot be correctly recognized due to its high bias in small samples. Table 2 shows the experimental results of the HSIC based KIKO model and the PHSIC based KIKO model on eight datasets. Compared with GROUND TRUTE, it can be observed that the KIKO model based on PHSIC can correctly recognize on most datasets. Table 3 shows the experimental results of the HSIC based HANM model and the PHSIC based HANMmode on eight datasets. It can be found that the HANM model based on PHSIC has significant advantages in identifying many to one causal relationships on small samples.

Figure 5 shows the experimental results on the ANM model. Figures 5 (a), (b), (c), and (d) show the experimental results of different methods on different sample sizes on the datasets CEP, SIM, SIM-G and SIM-ln, respectively. The accuracy of ANM-PHSIC in small samples is mostly higher than other indicators. As the number of samples increases, the accuracy of ANM-HSIC is consistent with that of ANMPHSIC. To verify the robustness of ANM-PHSIC under different noises, we add 10% gaussian noise and independent uniform noise to CEP and SIM data respectively. As shown in Figure 5 (e), (f), (g), (h), ANM-PHSIC also has strong robustness compared to other methods in most cases. In summary, ANM-PHSIC is more robust than other methods in cases of sparse samples and noise data.

5 Conclusion

This paper proposes a pure HSIC metric (PHSIC) to against random consistency, proving that PHSIC has tighter bounds than HSIC, and verifying its fairness and stability through theoretical analysis and simulation experiments. In addition, we apply PHSIC to multiple causal inference tasks, improving the robustness of these models in small sample and noisy scenarios. Furthermore, PHSIC can be applied to machine learning, deep learning and other fields. This is also the direction we will further apply in the future.

In conclude, it also gave us an inspiration. Unlike the traditional definition of independence when the sample tends to infinity, independence should be data-driven. From the perspective of random consistency, we can redefine the independence between random variables to construct robust and fair correlation indicators.

Acknowledgments

This work was supported by National Key Research and Development Program of China (No. 2021ZD0112400), the National Natural Science Foundation of China (Nos.62136005), the Science and Technology Major Project of Shanxi (No.202201020101006), the Key R&D Program of Shanxi Province, China (Grant no.202202020101004), the Natural Science Foundation of Shanxi Province, China (Grant nos.20210302123455).