Post-Nonlinear Causal Relationship with Finite Samples: A Maximal Correlation Perspective

Thank you so much for your careful reading and comments. We have carefully addressed all your questions below and revised the manuscript accordingly. We sincerely hope that our answers can fully address your concerns.

W1: Presentation of Lemma 4. Should add more assumptions?

Lemma 4 is in the context of PNL learning and follows the same assumptions as Lemma 3. As we have assumed the invertibility of underlying functions $\bar{f}$ and $\bar{g}$ in Lemma 4, one can always find $\hat{f},\hat{g}$ such that " the composition of PNL functions and un-mixing nonlinearities are linear". We have added more descriptions to make Lemma 4 complete.

W2: Identification bound of PNL has not improved. Require illustration of the signification impact.

We acknowledge that the identification result has not been improved. However, we clarify the identifiability results in the context of PNL learning in Lemma 4, which has not been discussed in the existing literature. Besides, the estimation algorithm development is also fairly important, since there is a lack of effective and efficient algorithms to recover the underlying nonlinearities. Our proposed method is easy to implement and has a benign optimization structure, which provides a practical and efficient way to learn the underlying nonlinearities. We hope this can bring new thoughts to both the machine learning and causal inference communities.

W3: In the "NONLINEAR FUNCTION FITTING" experiment, there are only two group generation mechanisms are used. More results with more types of non-linear functions should be provided if possible.

The two synthetic examples were presented for illustration, and one of them was used in Uemura & Shimizu (2020). Following your suggestions, we updated the supplement with more examples. Please find additional experiments in our anonymous code link.

Q1: Have you applied this method to real-world data? The gene data from the DREAM4 competition is designed to mimic real data. In addition, we also showed the causal discovery result with uncertainty quantification using a real dataset (Tuebingen) in Figure 10 in suppl.I.2.

[ Update: Response to W2 is updated ]

Thank you so much for your careful reading and comments. Our remarks and clarifications are as follows. We sincerely hope that our answers can fully address your concerns.

W1: The use of over-parameterized neural network (NN). We need to clarify that NN has been used in some previous PNL-based methods (e.g., MLP-PNL, AbPNL), but they are not used in our method. We observed the over-fitting and optimization issues when we tried to reproduce their results using NN, which motivates the benefits of using our proposed method based on the RFF. We did use different NN architectures (over-parameterized wide NN and a relatively smaller deep-narrow NN) to illustrate the issues of previous works, see supp. C, i.e., they both produce unsatisfying results.

The reasons that we include the over-parameterized NN for comparison are as follows: 1) it works well in many practical applications, and 2) it has a nice loss landscape for optimization (no bad basin/spurious valley under mild assumptions [Sun et al. 2020]).

W2: In analogy with the single-index model, is there any identification issue when expanding $f_1$ to a purely nonlinear one?

Thanks for pointing out the relation to the single-index model. If we understood correctly, you are referring to the noiseless case, $y =m(x)= f_2(f_1(x))$, which is not discussed in our manuscript. When $f_1(x)= \beta^T\phi(x)$, the PNL model reduces to a single-index model as defined in [Lin and Kulasekera, 2007]. And $\beta$ is identifiable when it has a unit norm and $m$ is a non-constant continuous function.

If $f_1$ is extended to a more flexible nonlinear function w.r.t the parameter $\theta$ (e.g., a NN parameterized by $\theta$), similar results MAY NOT hold, i.e., $m(x) = g(f(x))=h(q(x))$ but $\{g\neq h, f\neq q \}$. A counter-example (in scalar case) is $m(x)= \sin^2(x)$. The nonlinear transformations can be $\{g(z) = (\alpha z)^2,f(x)=\frac{1}{\alpha} \sin(x) \}$ and $\{h(z)= sin^2(\beta z), q(x) = \frac{1}{\beta} x \}$, where $\alpha,\beta \neq 0$ are some scaling factors.

However, in the PNL model in our paper, there is an independent noise term, and our Lemma 2 gives the corresponding identifiability result.

Please correct us if we misunderstood your question. We'll be happy to address further questions/comments.

Q1: Why use maximal correlation? There are obviously many "nonparametric correlations" at our disposal, e.g. d-corr, maximal information coefficients (by Reshef et al.), and many others.

Maximal Correlation is a natural choice for PNL learning, as the optimized transformations $f$ and $g$ are closely related to the underlying nonlinear functions $f_1,f_2^{-1}$ in the PNL model. However, other dependence measures do not involve such a transformation learning process. Another consideration is computational efficiency. Since we need to compute dependence during optimization, many known dependence measures (e.g., maximal information coefficient) are not suitable, see the time comparison in the RDC paper (Lopez-Paz et al., 2013). For the independence test, however, other dependence measures could be considered as potential replacements.

Q2: Any connection to CCA? Certainly, some components in our method (e.g., HGR and RDC) are strongly related to nonlinear CCA. One famous work along this research line is Kernelized CCA by Bach and Jordan (2002), where the connection to CCA was thoroughly discussed. The distinction between Kernelized CCA and RDC lies in the representations of $f$ and $g$: Kernelized CCA uses $K_1\alpha_1$ for $f(X)$ (representer theorem), while RDC uses the linear combination of RFFs, $\Phi(X)^T\alpha$.

R. Sun, D. Li, S. Liang, T. Ding and R. Srikant, "The Global Landscape of Neural Networks: An Overview," in IEEE Signal Processing Magazine, vol. 37, no. 5, pp. 95-108, Sept. 2020, doi: 10.1109/MSP.2020.3004124.

Lin, Wei, and K. B. Kulasekera. “Identifiability of Single-Index Models and Additive-Index Models.” Biometrika 94, no. 2 (2007): 496–501. http://www.jstor.org/stable/20441387.

Bach F R, Jordan M I. Kernel independent component analysis[J]. Journal of machine learning research, 2002, 3(Jul): 1-48.

Thank you so much for your careful reading and comments. We have carefully addressed your questions below and revised the manuscript accordingly. We sincerely hope that our answers have fully addressed your concerns.

W1 & Q1: On Presentation clarity.

We refer to the "meaningful solutions" as those solutions that can reflect the underlying nonlinear transformations. As one can see in the last two rows in Figure 6, the learned solutions cannot reveal the underlying functions as shown in Figure. 4, and thus are less meaningful. We have added a short description in the main paper to make this clear.

The causal discovery method, IGCI, is not designed for PNL learning and cannot give explicit function transformations of $f_1,f_2$, thus cannot provide "transparent and interpretable transformations". We have rewritten this phrase to "explicit nonlinear transformations".

For the convenience of the audience, we have added the formal definition of randomized dependence coefficient (RDC) in the updated version.

W2: How HGR was used to motivate PNL learning? As $Y = f_2 (f_1(X)+\epsilon)$ in the PNL setting, it is obvious that $X$ and $Y$ are highly correlated when the noise is small. Our goal is to learn $f_1,f_2$. Let $f(X) = f_1(X),g(Y) = f_2^{-1}(Y)$, and we hope to learn $f$ and $g$ through maximizing the HGR correlation $E(g(Y)f(X))$.

W3: The main contributions seem overstated. The reported numerical gain is compared with the SOTA PNL-based method. Our focus lies on comparing with interpretable baselines, as other baselines are incapable of providing explicit transformations for $f_1$ and $f_2$.

W4: While one of the disadvantages of the PNL algorithms is claimed to be the optimization issue, thus motivating the RFF parametrization and linear HSIC kernel, neither is the experimental performance of this variant discussed, nor is its theoretical properties. MC-PNL still uses a gradient-based algorithm for the universal kernel and banded loss regularizer.

We did include the theoretical properties in the text under eq (14). We show that the BCD algorithm can converge to a critical point and each sub-problem is a quadratic program. The experimental convergence results have been shown in Figure 3. Only for fine-tuning, the gradient-based optimization is needed to handle the universal kernel and banded loss regularizer.

Q1: see W1

Q2: Typo. Yes, it should be $<-\delta$ in Line 3 in Algorithm 2.

Q3: Is there any ablation study done to determine the parameters of the RFF?

We are unsure about the specific type of ablation study you are expecting. It would be very helpful if you could specify the question a little bit more.