Optimal Kernel Choice for Score Function-based Causal Discovery

Dear Reviewer VHXH,

Thank you very much for your efforts in reviewing this paper. We identify there exist some misunderstandings. And we have reexamined the notation in the paper fixing the typos. Here are some clarification of your concerns.

Q1. The difference from the existing Gaussian-process-based framework is unclear.

A1. I believe there might be some misunderstanding here. The contribution of our method does not lie in using Gaussian processes to optimize kernel hyperparameters. Our contribution builds upon Huang et al.'s method. Huang et al. utilize the characterization of general independence relationships with conditional covariance operators and propose using the regression model for causal discovery with the regression form $\phi(X; \sigma_x) = f(PA; \theta) + \epsilon$. In Huang et al.'s model, the parameter $\sigma_x$ in $\phi(\cdot)$ is pre-defined. It assumes that the pre-defined feature $\phi(X)$ is generated from a non-linear projection $f(PA_X)$ with an independent noise $\epsilon$ added, which is hard to satisfy with such an manually selected $\phi(\cdot)$ and may not fully capture the characteristics of the data in many cases. While in our method, we extended the regression model above with trainable parameters of $\phi(\cdot)$ rather than the fixed one.

By introducing more learnable parameters, our model becomes more flexible and can more accurately uncover relationships within the data. Additionally, due to the introduction of extra parameters, we propose using the joint distribution's marginal likelihood as the score for parameter learning and causal discovery, which is different from Huang's conditional likelihood based ones. To the best of our knowledge, there is currently no other approaches using the same model and score function for causal discovery.

And our proposed score function in Eq.13 is not the same with that of [1]. Our score has an extra term (last term in Eq.13), which reflects the difference of relationship modeling mention above. And this last term arises from transforming the joint distribution of dependent and independent variables into the joint distribution of the independent variables and the estimated mixed noise. It results from simplifying the Jacobian matrix obtained from using the change-of-variables theorem.

Q2. Relationship between $k_x$ and $\phi(x)$

A2. When we introduce an RKHS space $H_\mathcal{X}$ on $\mathcal{X}$, it comes with the definition of a continuous feature mapping $\phi_\mathcal{X}$ and a measurable positive-definite kernel function $k_\mathcal{X} : \mathcal{X} \times \mathcal{X} \rightarrow \mathbb{R}$. Since $\phi_X$ is in an infinite-dimensional space without a proper probability measure. For practical purposes, given an observation $x$ of $X$, we project $\phi_\mathcal{X}(x)$ into its empirical feature map, which is $\mathbf{k}_x$.

Q3. The connection between Eq.3 and Eq.4 is unclear.

A3. We are sorry for your confusion. The Eq.3 is derived from Theorem 7 in [1], which is the theoretical foundation that allows us to transform the independence relations test into a regression model selection problem. Eq.4 is the corresponding regression model, where a continuous feature mapping $\phi(X)$ is applied to the dependent variable $X$. Due to the properties of characteristic kernels, $\phi(X)$ can be decomposed into an infinite linear combination of basis functions $[\phi_1, \cdots, \phi_i, \cdots]$, and each $\phi_i \in {\mathcal{H_X}}$ is orthogonal to each other. We associate each $\phi_i$ with the function $g$ in Eq.3, leading to the conclusion that Eq.5 holds for any given $\phi \in {\mathcal{H_X}}$. Thus, based on Eq.5, we can further transform the conditional independence test into a parameter selection problem for the regression model defined in Eq.4.

Q4. In Eq.6, there is no definition or description about parameter. What kind of parameters do the authors consider?

A4. In the later part, we introduced Gaussian processes to mitigate model overfitting. By parameterizing $f$, we aim to reduce the number of parameters in $f$, representing the parameters $\theta$ as optimizable parameters in the Gaussian process.

Reference

[1] K. Fukumizu, F. R. Bach, and M. I. Jordan. 2004. Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces. Journal of Machine Learning Research 5 (2004), 73–79.

Unfortunately, I cannot clearly understand what the authors mean, mainly due to the grammatical errors in their response, though I will not enumerate them.

I posed two questions about notations, $\varepsilon \in H_X$? and $\phi(x) = k_x$?, but they seem not be fully answered (I understand the bold $k_x$ is a concatenated feature vector of observations; thank you for the response.). Regarding the latter, in usual setup in kernel methods, we often express feature mapping $\phi$ like as $\phi(x) = k_{X}(x, \cdot)$. However, from the authors' response, I am not sure whether the notations are introduced in the same way.

As pointed out by other Reviewer, if the related work is not sufficiently mentioned, the readers cannot understand the technical significance of the proposed method. I understand that developping GP-based formulation for Bayesian network structure learning like UAI2000 paper is not the proposal of this work. However, since both address similar problem setup and use similar GP-based formulations, the authors should clearly explain the difference. Although the authors mention the formulation difference in their response (i.e., the last term in Eq. 13), I cannot clearly understand the meaning of reflects the difference of relationship modeling mention above (mention -> mentioned?; What does relationship modeling mean?). If they claim that their model formulation is more flexible than UAI2000 paper, they should perform empirical performance comparison.

I am not sure, but I remember that UAI2000 work can also adapt the empirical Bayes framework for hyperaparameter tuning. If so, the contribution of the proposed method is to make it possible to perform gradient-based optimization for hyperparameter tuning (to my understanding). However, without any empirical comparison, we cannot see the technical significance of this contribution.

For this reason, I feel that the paper requires a major revision and is not ready for the publication at such top conferences as ICLR.

Dear Reviewer SijP,

Thank you very much for your efforts in reviewing this paper. We identify there exist some misunderstandings. And here are some clarification of your concerns.

Q1. Why joint marginal likelihood successfully avoids overfitting?

A1. What we want to highlight here is that using the likelihood of the joint distribution instead of conditional distribution is not to mitigate overfitting. Given that our model introduces additional learnable parameters (kernel function parameters imposed on the dependent variable), we choose to employ the likelihood of the joint distribution to constrain these parameters and learn a non-trivial relationship. As mentioned in the last paragraph of Section 3.2, if we were to use conditional likelihood for learning these extra parameters, the model would directly fail to capture real relationships but project all samples into a very small range of the feature space with excessively high scores, which is uninformative for causal discovery purposes. This is the reason why we must use the joint distribution—it fundamentally stems from the differences in the assumed models. Building on this, to alleviate model overfitting, we utilized Gaussian processes to parameterize the projection function $f$, aiming to enhance the effectiveness of causal discovery. The corresponding score also takes the form of joint marginal likelihood.

Q2. Using more complex kernels, such as deep kernels.

A2. Your valuable insights are appreciated. Due to the fact that our approach is based on the characterization of conditional independence with conditional covariance operators[1], it requires the use of the characteristic kernel. The use of more complex deep kernels, etc., will need further exploration.

Q3. Is Equation 11 the same as the conditional likelihood?

A3. It is not the same with the conditional likelihood. If it is a conditional distribution, in Equation (11), only the first term will be present, and there won't be the term $E[\log | det J_i|]$. That is, $p(x|pa, f, \sigma_x) = p(f(pa) + \widetilde{{\pmb \epsilon}}| pa, f, \sigma_x) = p(\widetilde{\epsilon})$.

Q4. You mentioned "maximizing the likelihood function itself may potentially lead to overfitting in structure learning", what do you mean by this? Could you elaborate about it?

A4. When using models like MLPs to represent $f$ in the model and directly employing MLE or minimizing mutual information to learn the model parameters, overfitting is prone to occur. The main drawback is that this model can fit perfectly to any data by learning constant functions. Consequently, the estimated noise $\tilde{\epsilon}$ will be constant and thus always independent from $X$ with a high score, leading to incorrect connections. Therefore, here we employ Gaussian Processes to parameterize $f$, with the expectation of reducing model overfitting. There are also some other methods to avoid overfitting, such as rank-based [2] or auto-encoding-based [3].

Q5. What do you mean by "mixture of independent noise variables"?

A5. We formulate the generative relationship among presumed variables as $\phi(x) = f(pa) + \epsilon$, where $\epsilon$ represents the independent noise variable vector, and each dimension is independent of the others. Subsequently, our objective involves learning the model parameters and conducting a causal graph search by fitting the likelihood of this set of independent noise variables. That is "mixture of independent noise variables".

Q6. Why joint distribution likelihood enforces constraints on the additional parameters?

A6. As we answered in A3, according to Eq.11, the conditional probability distribution cannot constrain the kernel function parameters imposed on $X$. In contrast, the joint distribution likelihood, with the inclusion of the term $E[\log|\mathrm{det} J_i|]$, provides constraints on the kernel function parameters imposed on $X$. This term plays a role in constraining the kernel function parameters by influencing the Jacobian determinant's expectation, providing valuable constraints during the parameter learning process.

Reference

[1] K. Fukumizu, F. R. Bach, and M. I. Jordan. 2004. Dimensionality reduction for supervised learning with reproducing kernel Hilbert spaces. Journal of Machine Learning Research 5 (2004), 73–79.

[2] Keropyan G, Strieder D, Drton M. Rank-Based Causal Discovery for Post-Nonlinear Models International Conference on Artificial Intelligence and Statistics, 2023: 7849-7870.

[3] Uemura K, Shimizu S. Estimation of post-nonlinear causal models using autoencoding structure. International Conference on Acoustics, Speech and Signal Processing, 2020: 3312-3316.

Dear Reviewer oMj1,

Thank you very much for your efforts in reviewing this paper.

Q1. This seems too incremental as a conference contribution and in my opinion is of insufficient originality.

A1. The choice of the kernel is a crucial step in kernel-based methods, as it directly impacts the accuracy of the results. Unlike previous approaches, which often heuristically select kernel parameters manually, our method introduces an optimization-based approach for adaptive kernel parameter selection. Our approach dynamically considers the characteristics of the given data and selects appropriate kernel parameters accordingly. Introducing additional parameters poses a challenge for existing conditional likelihood-based methods to supervise model parameters. Therefore, we propose a supervised approach based on mutual information. Our experimental results demonstrate the stability improvement of our proposed method compared to existing approaches.

Q2. Several key references missed.

A2. Our proposed method introduces a scoring approach that models relationships in a more generalized manner, focusing on the local associations between given potential independent and dependent variables. Based on our understanding, references [2-5] predominantly focus on searching for DAGs under the assumption of a relatively simple local relationship model. While the mentioned references [6-7] appear to address extreme scenarios involving hard intervention targets or missing variables in the field of causal discovery. Given the extensive scope and diverse content of causal discovery, we have chosen to align our work with a narrower field relevant to our method.

Q3. What is the additional computational cost of the method, when a possibly non-convex objective has to be optimized instead of the simpler objective by Huang?

A3. Due to the introduction of additional parameters (parameters in the kernel applied to the dependent variable), our method has a relatively higher computational complexity compared with score function proposed by Huang et al. However, the primary computational bottleneck in Gaussian processes lies in the Cholesky decomposition of $K_{PA}$ with a computational complexity of $O(n^3)$. Therefore, the computational complexity of both our method and Huang are primarily dominated by this term.

Dear Reviewer 2GNL,

Thank you very much for your efforts in reviewing this paper. We identify there exist some misunderstandings. And here are some clarification of your concerns.

Q1. Lack of in-depth discussions on different choices of kernels.

A1. In this paper, we do not select kernel across different kernel families. What we do in this paper is to first give a specific characteristic kernel (such as Gaussian kernel used in the paper) and apply it to the dependent variable, and then optimize the parameters of the kernel (such as the width in the Gaussian kernel) to find the optimal parameters that best fit the data.

Q2. The variance parameter of Marginal Likelihood of Huang et al. 2018 is actually trainable.

A2. We believe there might be some misunderstanding here. What we want to highlight is that our method can optimize the parameters of the kernel function imposed on the dependent variable, which is pre-defined heuristically and fixed in the existing paper. Specifically, in Eq.6 of our paper (i.e., $\phi(X; \sigma_x) = f(PA; \theta) + \epsilon$), the kernel function parameter $\sigma_x$ involved in $\phi(X)$ imposed on $X$ is learnable in our approach, while it is fixed in Huang et al.'s method. Due to the introduction of additional learnable parameters, we extended the regression model in RKHS and accordingly use joint distribution likelihood $p(X, PA |\\bf{\sigma})$ to simultaneously optimize $\sigma_x$ and other learnable parameters (i.e., $\theta$ in $f$) in the model including the mentioned noise variance.

Q3. The linear function class on the RKHS induced by RBF kernel is capable to be consistent with any underlying model.

A3. When the kernel is characteristic, it indeed possesses the property of fitting any function. However, what we aim to compare here are the differences in the models between our approach and Huang et al.'s method. In general, both methods assume a relationship model, which has the form: $g(X) = f(PA_X) + \epsilon$ and $\epsilon$ is the noise. The difference lies in the fact that in Huang's method, the parameters in $g$ are heuristically pre-defined and fixed, resulting in the fixed features $g(x)$. This assumes that , for such a pre-fixed $g(x)$, it is generated by a non-linear mapping $f$ of the independent variable $x$ with its parent variables $pa$, plus the independent noise $\epsilon$. This assumption is often hard to satisfy with most manually-selected $g$ and is not related to the functional space induced by characteristic kernel. That's why there are still the dependency between the error and the regressor. In contrast, our method optimizes parameters in both $f$ and $g$ simultaneously, which can select more proper parameters of $g$ with a suitable feature space and can better estimates the independent noise. This is what we want to express in the Motivation part.

Q4. Why it leads to a spurious edge between $X_1$ and $X_3$ in Caese 2 of Motivation?

A4. The GES method searches for all possible combinations of parent variables for a target variable and identifies the combination with the highest score. In contrast, in Huang et al.'s method, the score when $PA_{X_3} = (X_1, X_2)$ is higher than when $PA_{X_3} = (X_2)$, leading to incorrect edges. We attribute this error to the inability of Huang et al.'s method to adequately capture the relationships between $X_2$ and $X_1$ as well as between $X_3$ and $X_2$. This results in the estimated noise $\widetilde{\epsilon}_3$ not being independent of $\widetilde{\epsilon}_2$, as evident in Figure 1(d).

Q5. What is the difference in their dynamics that leads to potential performance differences?

A5. The likelihood of the conditional distribution is insufficient to constrain the parameters of the kernel function applied to the dependent variable. From the conditional distribution, we can observe that $p(x|pa, f, \sigma_x) = p(f(pa) + \widetilde{{\mathbf \epsilon}}| pa, f, \sigma_x) = p(\widetilde{\epsilon})$, lacking the term $E[\log|\mathrm{det} J_i|]$ in Eq.11 of the paper. Consequently, it will fail to capture real relationships if we use $\log p(X = x | PA, f, \sigma_x)$ as the score function with trainable parameters in $\phi(X)$. Specifically, it will project all samples into a very small range of the feature space with excessively high scores, which is uninformative for causal discovery purposes. Therefore, we must use $\log p(X = x, PA = pa | f, \sigma_x)$ as our score function, driven by the differences in the assumed models.

Dear Reviewer qtuB,

Thank you very much for your efforts in reviewing this paper.

Q1. Is there any significant technical difference compared with Huang et al. (2018) in this part?

A1. Section 4 is to demonstrate that our proposed score exhibits the property of local consistency, ensuring that our score can collaborate with the GES method for accurate causal discovery. Due to differences in the assumed models between our method and Huang et al., Eq.3 in our paper represents a joint distribution $p(PA_X, X | {\pmb \sigma})$, whereas in Huang et al., it represents a conditional distribution $p(X | PA, {\pmb \sigma})$. Additionally, the score in Huang et al. proposed does not involve optimizing kernel function parameters (i.e., $\sigma_x$ in Eq.6 of the paper) imposed on $X$, whereas our method simultaneously optimizes these parameters along with other parameters in the model, selecting optimal kernel parameters for the given data. Therefore, our approach differs in terms of both the model and the corresponding score function from Huang et al.