Score-informed Neural Operator for Enhancing Ordering-based Causal Discovery

We appreciate your valuable comments and suggestions.

W1-(a). Accuracy of Hessian diagonal approximation

To assess the accuracy of Hessian diagonal estimation, we provide a comparison between SciNO and the second-order Stein gradient estimator. Since we are unable to include images, we provide the mean and standard deviation of the estimation errors for each graph used in Figure 2 in the following table. As shown, SciNO yields more accurate and stable estimates, with errors closer to zero and lower variance than the Stein-based estimator.

Furthermore, as presented in answer for W3-(a) from Reviewer 2fQy, we compared SciNO to the SCORE, which is based on the Stein gradient estimator. SciNO generally achieves a lower Order Divergence than SCORE, particularly in settings high-dimensional datasets. This combined evidence highlights SciNO’s superiority over the second-order Stein gradient estimator in accuracy.

W1-(b). Stability of CaPS w/ SciNO on low-dimensional datasets

To assess potential performance degradation, we compared CaPS and SciNO on five low-dimensional datasets.

These results confirm that SciNO does not degrade performance in low-dimensional settings. In all five cases, SciNO either improved the performance of CaPS or maintained it.

W2. Connection between causal discovery and LLM control, with comparison to previous methods

Our research is motivated by [Vashishtha et al., 2025], which demonstrated the benefits of using LLMs for ordering tasks (L27-28). We aim to improve ordering-based causal discovery by guiding it with LLM priors. Because existing pairwise-query approaches (L67–68) cannot be directly applied to SciNO due to structural differences, we introduce a novel integration method in which the LLM sequentially predicts leaf nodes given the inferred causal order (L180–181). This method depends on the accuracy of the Hessian-diagonal statistics in Equation (8). Since Figures 2 and A.3 show that SciNO provides a more stable and accurate estimate than DiffAN, we focused our main experiments exclusively on SciNO. In response to the reviewer’s suggestion, we conducted additional experiments applying our control method to existing approaches.

1. Comparison with DiffAN-based control

   The table is provided in our response to Reviewer 2jQc’s Q3 and will be added to Table 4 of the revised manuscript. DiffAN-based control yields a 49% average improvement in order divergence over uncontrolled Llama, which is lower than SciNO’s 64%. This clearly demonstrates that our superior performance arises from both the novel control method and the SciNO’s stable, accurate evidence.

2. Comparison with SCORE-based control

   To benchmark against kernel methods (e.g. SCORE, CaPS) that cannot leverage Deep Ensembles [Lakshminarayanan et al., 2017], we used 1,000-sample bootstrapping. As noted in L186, our focus is the minimal-variance criterion eq.(3), so we conducted SCORE-based control. Although this yields large improvements on MAGIC-NIAB, it does not outperform SCORE alone on ECOLI70 and ARTH150, indicating that LLM priors lack robust synergy with kernel-based methods.

• SCORE-based control (Comparison of order divergence between uncontrolled Llama and controlled Llama)

  Method\Dataset	MAGIC-NIAB(d44)	ECOLI70(d46)	MAGIC-IRRI(d64)	ARTH150(d107)
  Llama-3.1-8b(1-norm)	9	32	18	80
  w/ control(rank)	1 (↓88.9%)	14 (↓56.3%)	11 (↓38.9%)	45 (↓43.8%)
  w/ control(CI)	1 (↓88.9%)	14 (↓56.3%)	11 (↓38.9%)	43 (↓46.3%)
  Llama-3.1-8b(0.5-norm)	14	32	17	80
  w/ control(rank)	1 (↓92.9%)	16 (↓50.0%)	10 (↓41.2%)	45 (↓43.8%)
  w/ control(CI)	2 (↓85.7%)	16 (↓50.0%)	7 (↓58.8%)	43 (↓46.3%)

• SCORE

  Method\Dataset	MAGIC-NIAB(d44)	ECOLI70(d46)	MAGIC-IRRI(d64)	ARTH150(d107)
  SCORE	8	7	17	38

Q1. Detailed explanation of what Figure 5's subgraph shows

As noted in L266-268, the bottom of the Figure 5 shows the probability of selecting a ground-truth leaf node at each causal ordering step on ECOLI70. When multiple ground-truth leaf nodes exist at a given step, we plot the highest probability among them.

The green bars represent the initial probabilities based solely on the LLM prior, $\mathbb{P}_{\mathrm{AR}}(x_{t+1}|x_{1:t},\mathtt{context}) $ from eq.(8). In high-dimensional settings, this prior often produces low, uninformative probabilities. Our control method combines these with SciNO’s evidence via eq.(8), yielding the posterior probabilities $ \mathbb{P}(x_{t+1}|x_{1:t},\mathtt{stat},\mathtt{context}) $ (shown as pink bars).

As discussed in L268-270, the dramatic rise in the ground-truth node’s probability demonstrates that, even if the LLM prior is inaccurate, the data-driven evidence compensates for this inaccuracy and enables reliable causal ordering. Crucially, during the initial steps when errors would otherwise accumulate, these probabilities are effectively boosted. We observe the same corrective effect consistently across high-dimensional datasets (Figure A.5).

---

[Vashishtha et al., 2025] — Causal Order: The Key to Leveraging Imperfect Experts in Causal Inference. ICLR, 2025.

[Lakshminarayanan et al., 2017] — Simple and scalable predictive uncertainty estimation using deep ensembles. NeurIPS, 2017.

We are grateful for your detailed and helpful comments.

Q1. Details on how noise injection destroys functional information

The functional information refers to the underlying functional relationships between variables, such as the derivatives of the score function. Diffusion models learn data distributions by injecting noise into the input, but this process can destroy these functional relationships. As shown in Figure 2, MLP-based diffusion models poorly estimate score derivatives, failing to preserve functional information. Furthermore, Figure 3 shows a weak correlation between MMD and OD when using a standard MLP architecture, providing additional evidence that MLP-based diffusion models fail to capture the functional information effectively.

Q2. Clarification of 'robust performance' (L150)

We agree that the phrase 'robust performance' can be vague. Our use of 'robust performance' was intended to describe the trend of performance as the number of nodes increases, rather than its standard deviation as shown in Table 1. To accurately reflect this, we will revise the sentence to state that SciNO 'achieves a performance advantage over baselines that grows with the graph's dimensionality, highlighting its superior scaling properties'.

Q3. Applicability of the control method to DiffAN

The probabilistic control method proposed in Section 3.3 uses Hessian‐diagonal statistics, therefore it can be applied to DiffAN as well. Fundamentally, its performance hinges on the quality of the evidence term in eq.(8), i.e., the Hessian‐diagonal statistics. As we shown in Figure 2 and Figure A.3, SciNO provides a more stable and accurate estimation of the Hessian than DiffAN, so we exclusively use SciNO for all control experiments in the paper.

The table below reports results obtained under the same settings as Table 4, but with the evidence term extracted from DiffAN. On high-dimensional graphs, using SciNO's evidence leads to an average 64% reduction in order divergence (Table 4), versus only 49% with DiffAN. This clearly demonstrates that our superior performance arises from both the novel control method and the stable, accurate evidence provided by SciNO.

• DiffAN-based control (Comparison of order divergence between uncontrolled Llama and controlled Llama)

  Method\Dataset	Physics(d7)	Sachs(d8)	MAGIC-NIAB(d44)	ECOLI70(d46)	MAGIC-IRRI(d64)	ARTH150(d107)
  Llama-3.1-8b(1-norm)	3	3	9	32	18	80
  w/ control(rank)	4 (↑33.3%)	3	2 (↓77.8%)	22 (↓31.3%)	7 (↓61.1%)	25 (↓68.8%)
  w/ control(CI)	3	4 (↑33.3%)	6 (↓33.3%)	16 (↓50.0%)	13 (↓27.8%)	52 (↓35.0%)
  Llama-3.1-8b(0.5-norm)	1	3	14	32	17	80
  w/ control(rank)	4 (↑300%)	3	1 (↓92.9%)	20 (↓37.5%)	8 (↓52.9%)	25 (↓68.8%)
  w/ control(CI)	4 (↑300%)	5 (↑66.7%)	8 (↓42.9%)	18 (↓43.8%)	14 (↓17.6%)	52 (↓35.0%)

• SciNO-based control (Comparison of order divergence between uncontrolled Llama and controlled Llama, Table 4)

  Method\Dataset	Physics(d7)	Sachs(d8)	MAGIC-NIAB(d44)	ECOLI70(d46)	MAGIC-IRRI(d64)	ARTH150(d107)
  Llama-3.1-8b(1-norm)	3	3	9	32	18	80
  w/ control(rank)	1 (↓66.7%)	5 (↑66.7%)	4 (↓55.6%)	11 (↓65.6%)	9 (↓50.0%)	21 (↓73.8%)
  w/ control(CI)	1 (↓66.7%)	3	2 (↓77.8%)	15 (↓53.1%)	10 (↓44.4%)	18 (↓77.5%)
  Llama-3.1-8b(0.5-norm)	1	3	14	32	17	80
  w/ control(rank)	1	4 (↑33.3%)	2 (↓85.7%)	10 (↓68.8%)	9 (↓47.1%)	21 (↓73.8%)
  w/ control(CI)	1	4 (↑33.3%)	2 (↓85.7%)	14 (↓56.3%)	10 (↓41.2%)	18 (↓77.5%)

Q4. Performance gap in SHD between SciNO and CAM

We would like to clarify that our work concentrates on identifying and addressing the limitations of Ordering-based Causal Discovery methods that identify leaf nodes based on the Hessian diagonal.

To better contextualize the performance gap, we evaluated CAM [Bühlmann et al., 2014] under the same experimental setup. Note that our datasets differ from those in [Lachapelle et al., 2020], [Wang et al., 2021] as they are newly generated. SciNO achieves consistently lower SHD on the high-dimensional ER datasets and superior performance on the semi-synthetic datasets for all datasets except ARTH150.

Therefore, our method enhances upon existing score-matching based Causal Ordering methods and is also highly competitive with established Causal Discovery benchmarks like CAM. As detailed in our response to Reviewer Vgze (W3) our method also significantly outperforms other classic algorithms like PC [Spirtes et al., 2000] and GES [Chickering, 2002]. We appreciate the reviewer’s suggestion and will include these results in the paper.

Q5-(a). Placement of case study in method section

We agree with the reviewer’s point about improving the paper’s flow. Section 3.2 includes a critical ablation study and analytical findings that support the design choices and effectiveness of our proposed method. Table 1 presents an ablation study comparing our Learnable Time Encoding (LTE) module with Positional Encoding (PE), justifying its use. Figures 2 and 3 further support our approach by showing that SciNO achieves more stable Hessian approximation and learns causal relationship. Including these results directly in the Method section was intended to help readers understand why and how the method works as they encounter each design component. On other hand, we agree that the memory efficiency analysis in Figure 4 fits better in the Experiments section. We will move it to Section 4 in the final version to improve the overall organization. We appreciate the reviewer’s helpful suggestion.

Q5-(b). Table caption formatting

Thank you for pointing this out. We will move all table captions to the top to adhere to the standard format.

Q5-(c). Abbreviation formatting consistency

We appreciate the reviewer’s attention to detail in identifying this inconsistency. To correct this, we will adopt a consistent rule for the first use of any abbreviation, for example: Fourier Neural Operator (FNO) [22] and Positional Embedding (PE) [40].

---

[Bühlmann et al., 2014] — CAM: Causal Additive Models, High-Dimensional Order Search and Penalized Regression. The Annals of Statistics, 2014.

[Lachapelle et al., 2020] — Gradient-Based Neural DAG Learning. ICLR, 2020.

[Wang et al., 2021] — Ordering based causal discovery with reinforcement learning. IJCAI, 2021.

[Spirtes et al., 2000] — Causation, Prediction, and Search. MIT press, 2000.

[Chickering, 2002] — Optimal Structure Identification with Greedy Search. JMLR, 2002.

We sincerely appreciate the reviewer’s valuable feedback.

W1. Deeper discussion on when and why LLM control method works

In high-dimensional settings, the LLM prior alone can be uninformative and may even assign high probability to incorrect nodes. As shown in Table 4, both GPT-4o and uncontrolled Llama struggled significantly with causal ordering when relying solely on semantic information. This is illustrated at the bottom of Figure 5, where the LLM prior probability for the ground-truth leaf node (green bars) on the ECOLI70 dataset remains very low across most steps, highlighting its lack of confidence in large variable spaces.

Our probabilistic control method addresses this issue by integrating the LLM’s semantic prior with SciNO’s data-informed evidence. As described in L268–270, this integration dramatically boosts the posterior probability of the true leaf node (pink bars at the bottom of Figure 5). Crucially, our method delivers substantial improvements in the early ordering stages when the risk of error accumulation is highest.

To investigate the impact of variable semantics, we conducted an ablation study that masked variable names and varied the proportion of available descriptions (Supplementary Section A.5.2). As noted in L1019–1021, Supplementary Table A.1 demonstrates a robust synergy between the LLM’s semantic information and SciNO’s data-based evidence. Even with only 10% of variable descriptions provided, our controlled model outperforms using SciNO alone in the majority of cases. This result confirms that our method can improve performance even with limited semantic information. Furthermore, by tuning the temperature parameter τ, we can adjust the balance between prior information and evidence to optimize performance across different scenarios.

W2. Clarification on architectural contributions

As stated in L133-141, our architecture incorporates two key modifications to the Neural Operator from [Lim et al., 2023], which was designed for image generation tasks. First, we decompose complex-valued Fourier features into real and imaginary parts to better capture functional structures. Second, we integrate a Learnable Time Embedding (LTE) module which is similar to [Li et al., 2021] designed for multi-dimensional spatial Positional Embedding (PE) in Vision Transformer (L145).

Since our model is based on a continuous-time diffusion process, we adopt LTE as a more suitable alternative to standard PE, a choice validated by our results in Table 1 where LTE shows a significant performance advantage (L137–138). Additionally, as shown in our response to Reviewer hiE3 (Q4-(b)), our multiplicative use of LTE offers substantial gains over a standard additive approach.

To further clarify the effectiveness of our modifications, we conducted a new ablation study applying LTE to DiffAN. These results demonstrate that simply applying LTE to DiffAN fails to produce the same performance gain, especially on high-dimensional datasets (e.g., 82.0 vs. 20.8 on ARTH150).

To summarize, rather than naively combining existing components, we integrate architectural modifications, guided by Neural Operator theory (Section 3.1), to construct a functional diffusion model that effectively estimates the Hessian diagonal (Figure 2) and captures causal relations (Figure 3).

W3. Comparison with non-ANM baselines

We would like to clarify that our work aims to resolve key issues in ordering-based Causal Discovery methods that utilize Hessian diagonal estimation, rather than proposing a framework that outperforms all existing Causal Discovery methods.

In response to the reviewer’s suggestion, we additionally conducted experiments comparing our method with the constraint-based PC algorithm [Spirtes et al., 2000] and the score-based GES algorithm [Chickering, 2002]. While our primary evaluation metric is Order Divergence, we converted the outputs of all methods to CPDAGs and measured the Structural Hamming Distance(SHD-C) to ensure a fair comparison. These results demonstrate that our method consistently outperforms both PC and GES across ER datasets.

As we show in our response to Reviewer 2jQc (Q4), our method is also highly competitive with CAM [Bühlmann et al., 2014]. This highlights the advantages of our approach against representative methods beyond the ANM family: it achieves superior accuracy while also being more scalable, overcoming the high computational complexity that limits these classic methods. We will add this comparison to our paper to better position our work within the wider causal discovery literature.

---

[Lim et al., 2023] — Score-based Generative Modeling through Stochastic Evolution Equations in Hilbert Spaces. NeurIPS, 2023.

[Li et al., 2021] — Learnable Fourier Features for Multi-Dimensional Spatial Positional Encoding. NeurIPS, 2021.

[Spirtes et al., 2000] — Causation, Prediction, and Search. MIT press, 2000.

[Chickering, 2002] — Optimal Structure Identification with Greedy Search. JMLR, 2002.

[Bühlmann et al., 2014] — CAM: Causal Additive Models, High-Dimensional Order Search and Penalized Regression. The Annals of Statistics, 2014.

We would like to thank the reviewer for very helpful suggestions.

W1-(a). Strong Sobolev space assumption

We appreciate detailed questions about our theoretical results.

The condition $k > 2 + \frac{D}{2}$ for Theorem 3.1 (pointwise approximation) appears strong in high dimensions. However, our method's practical success and theoretical grounding rely on a weaker condition. The key is that for our proposed causal discovery algorithms, which compute statistics like the variance eq.(3) or mean eq.(4) of the Hessian diagonal, we only need convergence in the energy norm, not the stronger pointwise sense. We provide a specific guarantee for this in Theorem A.1, which does not require the demanding $k>D/2$ condition. Therefore, Theorem A.1 provides the sufficient theoretical result for our approach. We presented Theorem 3.1 to show that under stronger conditions (which are still comparable to the infinite-differentiability assumption in the baseline SCORE method ), our operator also achieves the pointwise approximation.

W1-(b). Inapplicability near t=0 for time-conditioned scores

Thanks for your remarkable point. Our theory is designed to approximate the target score function $\mathbf{S}\_{\ast} = \nabla\_{\mathbf{x}} \log P\_{\text{data}}$ (L830-836 and Theorem A.1) via neural operators $\widehat{\mathbf{S}}^{\theta}(\cdot)$. The connection to the marginal score model follows a fundamental procedure in diffusion models. Approximation to marginal score functions with a time variable via neural operators in Hilbert space, $\widehat{\mathbf{S}}^{\theta}(t,\cdot) \approx \mathbf{S}(t, \cdot)$ for each $t$, is proved in prior work on HDM [Lim et al., 2023]. Hence the completeness of Hilbert space implies that the distance between $\widehat{\mathbf{S}}^{\theta}(\cdot)$ and the sequence of functional diffusion models will close to zero as we train score models accurately. We will clarify this by adding a paragraph to the appendix.

W2-(a). Limited number of runs in Section 4.1-4.3

To address these concerns about statistical robustness, we performed additional experiments, and the new results will be included in our revision.

1. Experiments in Section 4.1-4.2

   We increased the number of independent runs on real datasets from 5 to 10 to align with our ER setup. As shown in Table 2, SciNO consistently outperforms DiffAN across nearly all datasets and metrics (with a minor exception for SHD on Sachs). In Table 3, SciNO continues to achieve performance comparable to CaPS and exhibits improved order divergence on higher-dimensional graphs such as MAGIC-IRRI and ARTH150.

   • Table 2

     Dataset\Method(metric)	DiffAN(OD)	DiffAN(SHD)	DiffAN(SID)	SciNO(OD)	SciNO(SHD)	SciNO(SID)
     Physics(d7)	3.3±0.78	8.6±2.06	16.8±4.73	1.9±0.54	5.8±1.60	8.2±4.62
     Sachs(d8)	5.7±2.69	22.8±4.56	26.1±10.34	5.6±2.65	23.1±3.86	23.3±9.10
     MAGIC-NIAB(d44)	8.8±6.85	89.1±25.85	227.2±157.64	4.1±0.70	72.0±4.52	125.5±28.08
     ECOLI70(d46)	21.8±4.26	111.0±11.47	684.8±80.67	12.8±3.19	90.3±7.13	500.9±58.16
     MAGIC-IRRI(d64)	12.0±3.87	146.9±13.46	321.9±64.64	10.6±1.69	144.5±4.43	254.4±39.66
     ARTH150(d107)	43.3±15.77	613.3±72.49	2456.0±953.08	21.4±2.76	515.6±15.26	1186.7±104.19

   • Table 3

     Dataset\Method(metric)	CaPS w/ SciNO(OD)	CaPS w/ SciNO(SHD)	CaPS w/ SciNO(SID)
     MAGIC-NIAB(d44)	37.5±0.50	165.4±0.5	1185.9±5.61
     ECOLI70(d46)	25.4±0.80	141.3±4.78	882.4±13.02
     MAGIC-IRRI(d64)	41.6±0.80	183.5±2.97	1280.7±36.83
     ARTH150(d107)	48.2±0.87	380.0±4.63	2799.1±73.26

2. Experiments in Section 4.3

   For Llama-3.1 and SciNO, we evaluated performance using fixed-weight models. SciNO was trained using Deep Ensembles [Lakshminarayanan et al., 2017], which inherently account for initialization variance and thus reduce variability across trials. Accordingly, we recorded a single measurement per dataset.For GPT-4o, we initially reported the best result from 3 runs to approximate the model’s potential performance ceiling. To address concerns about generalizability, we have increased the number of runs to 10 and will report both the mean and standard deviation. Our original conclusion that GPT-4o exhibits limitations on high-dimensional causal graphs remains valid.

   Method\Dataset	Physics(d7)	Sachs(d8)	MAGIC-NIAB(d44)	ECOLI70(d46)	MAGIC-IRRI(d64)	ARTH150(d107)
   GPT-4o	1.1±0.54	4.4±1.11	11.9±0.94	41.1±2.30	20.9±2.30	77.1±7.46

W2-(b). Overgeneralization of performance claims

To better reflect the results, where SciNO primarily excels in the Order Divergence (OD) metric (on all datasets but Sachs), we will revise the text to be more precise: "Table 2 shows that SciNO enhances the performance of DiffAN in most metrics, especially in terms of OD."

W2-(c). Overstatement of SciNO’s advantages in MAGIC-IRRI and ARTH150

We agree that your point regarding Table 3 is valid. Our original intention was to emphasize the performance gains in terms of Order Divergence (OD), which directly measures the accuracy of causal ordering. To avoid confusion, we will revise the sentence to explicitly state “in terms of OD.”

W3-(a). Missing empirical comparison with SCORE

We agree that a direct comparison with the SCORE baseline is essential, and we include the results below. As the results show, SciNO performs comparably to SCORE on low-dimensional graphs and demonstrates a significant advantage on high-dimensional datasets. This supports our paper's claim that SciNO enhances scalability and performance in higher-dimensional settings.

W3-(b). Terminology for 'Real-World' vs. 'Semi-Synthetic' Data

In this paper, we followed the terminology used in prior works (SCORE, DiffAN, and CaPS), which refer to datasets constructed from real-world graphs as real data. However, we agree that classifying them as 'real data' may lead to confusion. We will update the paper to consistently define these datasets as 'semi-synthetic'.

W3-(c). Missing empirical runtime comparison with CaPS

While our original Figure 4 focused on memory efficiency, we have conducted the experiments comparing the runtime of CaPS against SciNO on the MAGIC NIAB(d44) dataset, varying the number of samples. Unlike CaPS, SciNO is trained in a multi-step fashion, so its total runtime depends on the number of epochs. Although GPU parallelization can reduce sample-dependent runtime, it does not mitigate the cost associated with iterative training. Therefore, we report the runtime per epoch for clarity. As demonstrated in the 100,000-sample setting, where CaPS runs out of memory while SciNO completes successfully.

We will add a plot of these results to the experimental section of our revised paper. We thank the reviewer again for this valuable suggestion.

W3-(d). PE vs. LTE ablation on semi-synthetic datasets

To further validate our design choice, we performed the ablation study comparing the Order Divergence (OD) of Positional Embedding (PE) and our proposed Learnable Time Embedding (LTE) on semi-synthetic datasets. As the results show, the performance advantage of our proposed LTE is most pronounced on high-dimensional datasets, which is consistent with our findings on ER graphs.

---

[Lakshminarayanan et al., 2017] — Simple and scalable predictive uncertainty estimation using deep ensembles. NeurIPS, 2017.

[Lim et al., 2023] — Score-based Generative Modeling through Stochastic Evolution Equations in Hilbert Spaces. NeurIPS, 2023.

• Table 3: note that the CaPS is a kernel-based method which produces identical results when we run multiple experiments. As the reviewer suggested, we provide the complete Table 3 below with the re-evaluated results on SciNO. We will include the result in the revised manuscript.

  Dataset\Method(metric)	CaPS (OD)	CaPS (SHD)	CaPS (SID)	CaPS w/ SciNO(OD)	CaPS w/ SciNO(SHD)	CaPS w/ SciNO(SID)
  MAGIC-NIAB(d44)	36	160	1024	37.5±0.50	165.4±0.5	1185.9±5.61
  ECOLI70(d46)	25	134	786	25.4±0.80	141.3±4.78	882.4±13.02
  MAGIC-IRRI(d64)	42	179	1192	41.6±0.80	183.5±2.97	1280.7±36.83
  ARTH150(d107)	49	406	3051	48.2±0.87	380.0±4.63	2799.1±73.26

• Table 4: note that the Llama’s pre-trained weights are unchangeable, so the source of randomness is due to the SciNO’s initialization. In the control experiment, we have used $M=30$ (L262) independently trained models for computing eq.(10) and eq.(12) via Deep Ensemble, which greatly reduces the effect of random initialization. Hence the rest of Table 4 results remain unchanged even though we run multiple experiments. We provide the complete Table 4 below with the re-evaluated results on GPT-4o. The confidence intervals of SciNO's results will also be finalized in the revised manuscript when the multiple deep ensemble training (which takes about couple of weeks) is done.

  Method\Dataset	Physics(d7)	Sachs(d8)	MAGIC-NIAB(d44)	ECOLI70(d46)	MAGIC-IRRI(d64)	ARTH150(d107)
  GPT-4o	1.1±0.54	4.4±1.11	11.9±0.94	41.1±2.30	20.9±2.30	77.1±7.46
  Llama-3.1-8b(1-norm)	3	3	9	32	18	80
  w/ control(rank)	1	5	4	11	9	21
  w/ control(CI)	1	3	2	15	10	18
  Llama-3.1-8b(0.5-norm)	1	3	14	32	17	80
  w/ control(rank)	1	4	2	10	9	21
  w/ control(CI)	1	4	2	14	10	18

For clarification, as the reviewer suggested, we promise to include the mean and standard deviation over 10 independent experiments in Tables 3 and 4.

To alleviate the reviewer’s concern on multiple experiments, we have conducted supplementary multiple experiments with the reduced number of ensembling (M=3), because running the full experiment with M=30 takes a couple of weeks. We will include the complete results after the discussion period.

We provide the modified experiments below with the re-evaluated results. The results confirms our method's consistency. Under rank-based control, the performance with M = 3 closely matches the results in Table 4. For the CI-based control, we observe slight drops on some datasets, an expected result of using fewer models for CI estimation. Importantly, even with M = 3, our method still yields significant improvements over the uncontrolled Llama and outperforms standalone SciNO (reported in Table 2 of our response to W2-(a)) on high-dimensional datasets (ARTH150 being the only exception).

These findings suggest our main results with a full M = 30 ensemble will be even more robust. We will include the complete 10-run results for M = 30 in the final manuscript.

We sincerely thank the reviewer for the constructive feedback.

Q1. Concern over whether the method belongs in causal discovery literature

As the reviewer pointed out, score approximation in function space is a core component of our method. We want to clarify that our primary goal is to enhance ordering-based causal discovery (e.g., DiffAN, CaPS) motivated by [Vashishtha et al., 2025], which studies that predicting causal orderings rather than full graph structures can improve consistency and reduce structural errors when combined with LLMs (L27–28). Therefore, SciNO ultimately serves the purpose of estimating the causal order, rather than a generic score approximation. Hence, we focus on ordering-based causal discovery algorithms that require Hessian diagonal approximation under the ANM assumption, as we clarify in Introduction (L30-46) and Supplementary Section A.1.1 (L762-764). We also included non-score matching related works in A.1.1.

To address the concern on the comparison against other methods, we have conducted additional experiments, such as PC [Spirtes et al., 2000], GES [Chickering, 2002], and CAM [Bühlmann et al., 2014]. See the answers to Reviewer Vgze’s W3 and Reviewer 2jQc’s Q4 (We apologize that we are unable to include the table due to the length constraint). Our results demonstrate that SciNO achieves superior performance on ER datasets compared to PC and GES, and also outperforms CAM in high-dimensional settings.

We will revise the related work section to better contextualize our work. We will add a discussion of the ANM methods that build the graph from the roots and contrast them with our approach, clarifying our motivation. We appreciate the reviewer’s valuable comment.

Q2. Connect between causal discovery and LLM control

As we mentioned in the answer to Q1, [Vashishtha et al., 2025] highlighted the benefits of using LLMs for ordering tasks, providing a key motivation for our research (L27–28).

Recent causal discovery studies have shown that integrating LLM prior knowledge can boost performance. Existing LLM-based methods, however, rely on querying pairwise relationships between variables (L66–69), which suffers quadratic complexity (L73). This prompting strategy is structurally incompatible with ordering-based approaches (e.g., DiffAN, CaPS, SciNO) that sequentially predict leaf nodes. We resolve this problem by treating the LLM’s generation as sequential leaf node prediction given the inferred causal order (L180–181), thereby enabling us to seamlessly combine the LLM priors and SciNO-based control with linear complexity (L73).

Therefore, our contribution is summarized as follows: enhancing ordering-based causal discovery algorithms, and bridging the gap when combining LLM prior knowledges and causal ordering algorithms. We aim to improve ordering-based causal discovery by guiding it with LLM priors, motivated by [Vashishtha et al., 2025], which concludes that SciNO and LLM control are not standalone topics. We will revise the introduction to make this clearer.

Q3. Clarification of Sections 3.1–3.2

We deeply appreciate the reviewer’s comment, that help us improve the paper’s clarity. Contrary to MLP-based approach (DiffAN), SciNO is founded on neural operator theory. Since MLP-based models are numerically unstable when estimating the higher-order derivatives, DiffAN has limitation in high-dimensional settings despite of its motivation. Hence Section 3.1 explains (L112-115) our theoretical foundation why score matching in function spaces is necessary for approximating the Hessian diagonal eq.(6) and Theorem 3.1 demonstrates approximation power of a family of neural operators including SciNO. Due to the page limitation, we present architectural and operational details in Section A.3. For clarity, we will move the Remark A.4 (L863-867) to the main part, which explains the theoretical connection between two sections and supplement high-level description.

1. FFT (inverse FFT) module in each Fourier layer receives (outputs) $\mathbb{R}^{H}$-valued tensor.
2. Theorem 3.1 presents a universal approximation of neural operator family for estimating score function and its derivatives. As approximating eq.(6) requires higher-order derivatives, $k$ in Theorem 3.1 denotes the order of differentiability when we define the Sobolev space. We will add a sentence clarifying $k$ in eq.(7).
3. Baseline structure of SciNO is motivated from FNO, which uses Fourier transform to formulate the operator learning in Sobolev spaces. Thus, Fourier transform is used to ensure high-order derivative approximation. We will add this explanation for clarity.

Q4-(a). Missing baseline details

We had previously omitted a detailed description of the baseline architecture from [Lim et al., 2023] since it is detailed in the original paper. Following your feedback, we will add a description of the [Lim et al., 2023] architecture and a comparison with ours in Supplementary Section A.3.

Q4-(b). Ablation study: multiplicative vs. additive LTE

We have conducted an experiment based on Order Divergence to validate our choice of a multiplicative approach over the more common additive one. As the results demonstrate, multiplicative LTE combination consistently yields better performance, providing empirical validation for our design.

Q4-(c). Clarification on the novelty of our contribution

To address the concern that our contribution may be incremental, we want to clarify that our main novelty does not lie in the architecture itself. Our key contribution is the proposal of SciNO as a functional diffusion framework which can accurately estimate the Hessian diagonal (L76-78). This capability is non-trivial, since standard diffusion models inherently fail at this task as we demonstrate in Figure 2. Furthermore, as shown in Figure 3, SciNO can learn the underlying causal relationships in the data.

Our proposed framework is not simply the application of a known heuristic. As discussed in our response to Reviewer Vgze (W2), simply applying the LTE to a standard diffusion model does not guarantee the performance improvements we achieved. Our work shows that the beneficial combination of our Neural Operator architecture and the LTE leads to our effectiveness on this task.

We validate this framework through a theoretical proof (Section 3.1) and empirical results in our experiments, which we believe constitutes a significant contribution.

Q5. Misuse of the term 'causal representation'

We agree that our use of 'causal representation' could create confusion with the distinct field of causal representation learning. Our intention was to emphasize that SciNO learns the underlying functional information of the data distribution. To improve clarity, we will replace 'causal representation' with 'causal relationship' throughout the manuscript.

Q6. Clarification of Section 3.3: Node Set V and i.i.d. Assumption

1. Definition of Node Set $\mathcal{V}$

   The node set $\mathcal{V}$ is defined in Section 2 (L90). To enhance clarity, we will restate this definition in Section 3.3 of the revised manuscript.

2. Justification for the i.i.d. Assumption

   The Central Limit Theorem applies to a set of independent and identically distributed (i.i.d.) random variables. In our method, we consider each variance estimate $ (\sigma_i^{(m)})^2 $ from the M models as a random variable.

   As stated in L190–191 and L193, each of the M models is trained and evaluated independently, so their predictions are independent. Moreover, since all models share the same architecture and are trained on the same data with identical hyperparameters, their outputs are drawn from the same distribution.

   Therefore, the variance estimates $ (\sigma_i^{(m)})^2 $ satisfy the i.i.d. condition.

Q7. Missing nonlinear ANM details in Synthetic Experiments

The reviewer is correct that nonlinearity is crucial for identifiability in Gaussian ANMs. As stated on L226, we generate the data using a Gaussian Process with an RBF kernel, which ensures that the functional relationships are nonlinear. To improve clarity, we will revise the text to state: '...we sample node values according to a nonlinear additive noise model (ANM) with Gaussian noise, where the functions are generated by a Gaussian Process with an RBF kernel.'

---

[Spirtes et al., 2000] — Causation, Prediction, and Search. MIT press, 2000.

[Chickering, 2002] — Optimal Structure Identification with Greedy Search. JMLR, 2002.

[Bühlmann et al., 2014] — CAM: Causal Additive Models, High-Dimensional Order Search and Penalized Regression. The Annals of Statistics, 2014.

[Vashishtha et al., 2025] — Causal Order: The Key to Leveraging Imperfect Experts in Causal Inference. ICLR, 2025.

[Lim et al., 2023] — Score-based Generative Modeling through Stochastic Evolution Equations in Hilbert Spaces. NeurIPS, 2023.

For example, the authors claim that the main contribution of the paper is to propose a novel causal discovery algorithm.

For clarification, we did not claim that SciNO is a novel causal discovery algorithm. I wonder why the reviewer thinks so, since we didn't claim that in the paper or in our rebuttal. One of our paper's contributions is enhancing the performance of Hessian diagonal-based causal discovery algorithms. I wonder what the reviewer's rationale is for claiming that CaPS, which is a recent algorithm among Hessian diagonal-based causal discovery algorithms, is not suitable as a baseline.

authors did not provide the details of the outcome of the experiments

We copy the answer to Reviewers 2jQc and Vgze here. To better contextualize the performance gap, we evaluated CAM [Bühlmann et al., 2014] under the same experimental setup. Note that our datasets differ from those in [Lachapelle et al., 2020], [Wang et al., 2021] as they are newly generated. SciNO achieves consistently lower SHD on the high-dimensional ER datasets and superior performance on the semi-synthetic datasets for all datasets except ARTH150. Therefore, our method enhances upon existing score-matching based Causal Ordering methods and is also highly competitive with established Causal Discovery benchmarks like CAM.

In response to the reviewer’s suggestion, we additionally conducted experiments comparing our method with the constraint-based PC algorithm [Spirtes et al., 2000] and the score-based GES algorithm [Chickering, 2002]. While our primary evaluation metric is Order Divergence, we converted the outputs of all methods to CPDAGs and measured the Structural Hamming Distance(SHD-C) to ensure a fair comparison. These results demonstrate that our method consistently outperforms both PC and GES across ER datasets.

I find the current write-up in the paper on LLM extremely unclear. I am not sure how the authors are considering the uncertainty coming from the LLM prediction.

If you could point out which parts are "extremely unclear" in the LLM experiments presented in this paper, we would be able to clarify them. The LLM control part proposed in this paper is designed to manage the possibility that LLM predictions could be uncertain or incorrect, and this is explained in Section 3.3. While the reviewer claimed that comparison with existing techniques is necessary, the appropriate comparison target for the LLM control technique proposed in this paper is the LLM prediction results without the control applied.

We would appreciate considering the feedback if the reviewer can provide the detailed reason for claiming that our experimental results distract readers.