Score-informed Neural Operator for Enhancing Ordering-based Causal Discovery
Authors: Jiyeon Kang
Songseong Kim
Chanhui Lee
Doyeong Hwang
Joanie Hayoun Chung
Yunkyung Ko
Sumin Lee
Sungwoong Kim
Sungbin Lim
Venue: NeurIPS 2025 poster ()
Keywords: Causal Discovery
Score Matching
Diffusion Models
Neural Operator
Primary area: probabilistic_methods
TL;DR: We enhance score matching based causal discovery algorithms through a specially designed neural operator trained via diffusion modeling in smooth function spaces.

Abstract: Ordering-based approaches to causal discovery identify topological orders of causal graphs, providing scalable alternatives to combinatorial search methods. Under the Additive Noise Models (ANMs) assumption, recent causal ordering methods based on score matching require an accurate estimation of the Hessian diagonal of the log-densities. However, previous approaches mainly use Stein gradient estimators, which are computationally expensive and memory-intensive. Although DiffAN addresses these limitations by substituting kernel-based estimates with diffusion models, it remains numerically unstable due to the second-order derivatives of score models. To alleviate these problems, we propose Score-informed Neural Operator (SciNO), a probabilistic generative model in smooth function spaces designed to stably approximate the Hessian diagonal and to preserve structural information during the score modeling. 
Empirical results show that SciNO reduces order divergence by 42.7% on synthetic graphs and by 31.5% in real-world datasets on average compared to DiffAN, while maintaining memory efficiency and scalability. Furthermore, we propose a probabilistic control algorithm for causal reasoning with autoregressive models that integrates SciNO's probability estimates with autoregressive model priors, enabling reliable data-driven causal ordering informed by semantic information. Consequently, the proposed method enhances causal reasoning abilities of LLMs without additional fine-tuning or prompt engineering.

Score-informed Neural Operator for Enhancing Ordering-based Causal Discovery
Jiyeon Kang
Songseong Kim 1.2r, Chanhui Le Doyeong Hwang
Joanie Hayoun Chung', Yunkyung Ko², Sumin Lee", "ee",,ungooon Kim Sungbin Lim 1,2
LLG Al Research
Department of Statistics, Korea University
3 Department of Artificial Intelligence, Korea University
Abstract
Ordering-bass -based approaches causal discovery identify topological orders of causal graphs. providing scalable alternatives combinatorial search methods. Under the Additive Noise Model (ANM) assumption, recent causal ordering methods based score matching require accurate estimation of the Hessian diagonal of the log-densities. In this paper, we aim to improve the approximation diagonal of the log-densities, thereby enhancing the performance of ordering- based causal discovery algorithms. Existing approaches that rely on Stein gradient estimators are computationally expensive and weemory-intensive while diffusion- model-based methods remain unstable due to the second-order derivatives of score models. To alleviate these problems, propose Score-infor informeed Neural Operator (SciNO), a probabilistic generative model smooth function spaces designed to stably approximate the Hessian diagonal and preserve structural information during the score modeling. Empirical results show that SciNO reduces order di- vergence by 42.7% on synthetic graphs and by 31.5% on real-world datasets on average compared to DiffAN, while maintaining memory efficiency and scalability. Furthermore, we propose a probabilistic control algorithm for causal reasoning with autoregressive models that integrates SciNO's probability estimates with au- toregressive model priors, enabling reliable data-driven causal ordering informed by semantic information. Consequently, the proposed method enhances causal reasoning abilities of L.LMs without additional fine-tuning or prompt engineering.
the Hessian
Introduction
Ordering-basse baarin causal discovery aims to identify topological ordering of nodes in causal graph, typically represented as a Directed Acyclic Graph (DAG), such that the order reflects the underlying cause-effect relationships. The combinatorial search over DAG structures is known to be NP-hard, with complexity increasing sharply as the number of variables grows [6]. In contrast, ordering- approaches first infer a topological order of the causal graph, and then determine the direction of edges based on this order, which substantially reduces the number of candidate graph structures. By leveraging Large Language Models (LLMs), predicting causal orderings rather than full graph structures, improves consistency and reduces structural errors [42]. Moreover, such causal order suffices to identify valid backdoor adjustment sets, without full graph discovery [42, Proposition 3.2].
Recently, causal ordering methods based on score matching have been proposed [33, 35, 45] under the ANM (Additive Noise Model) assumption. All of these approaches determine the causal order by
'Equal contribution.
CCoresponding Authors. [swkim01, sungbin)@korea.ac.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

iteratively identifying and removing leaf nodes modes-ie., variables that do not influence any others in the causal graph. To accomplish this, they require the estimation of the Jacobian of the score function S(x) = Vx log P(x), in particular its diagonal, denoted as the Hessian diagonal function 2 (x):
D (x) == diag (Hing r) (x)
:== D).
(1)
Score matching based causal ordering methods sequentially identify Icaf nodes from a causal graph by estimating variance or mean of the Jacobian of the score function after each leaf removal. SCORE [33] and CaPS [45] use the second-order Stein gradient estimator [23]. However, kernel-based estimators exhibit cubic computational complexity with respect to the number of samples N, and suffer from numerical instability due to the kernel matrix inversion. Since the score function must be re-estimated iteratively every step, the procedure becomes computationally problematic in large-scale causal graphs as causal discovery algorithms require more data samples as the number of nodes D increases. Contrary to kernel-based approaches, DiffAN [35] estimates the Jacobian of the score through the Denoising Diffusion Model [13]. Instead of retraining the score model, DiffAN approximates residue terms by computing the second-order derivatives of initially trained score model at each ordering step. These residual approximations reduce computational burden and allow scalable causal ordering when the number of data points increases, but their performance degrade in high-dimensional settings.
Fundamentally, to identify leaf nodes via the score information, more reliable derivative estimation of score models is necessary; otherwise errors may accumulate throughout the ordering process, ultimately degrading the causal discovery performance. In particular, we focus on improving the approximation of the Hessian diagonal of the log-densities, which is critical for reliable ordering- based causal discovery. While modeling the score function via diffusion models is suitable from a generative perspective, injecting noise into the data can destroy the functional information inherent in the original score function. To address this challenge, inspired by recent advances in score-based generative modeling in Hilbert spaces [25, 26]. we propose SciNO a functional diffusion model framework which aims to stably approximate the Hessian diagonal of the data distribution by modeling the score function with neural operators [16, 24, 26]. To stably learn the derivatives of the score function, we introduce two key modifications to the conditondt
Neural Operator),
Fourier neural operator architecture used in [26]. First, instead of using positional embeddings [43], we incorporate a Learnable Time Encoding (LTE) module that enables the model to jointly learn spatiotemporal derivatives. Second, signals in Fourier layers are decomposed into their real and imaginary parts in the spectral domain, allowing more expressive representation of functional information. Consequently, SciNO enables stabilizing the estimation of both the score function and its second-order derivatives. thus it significantly improves performances of causal ordering methods, including DiffAN and CaPS, especially in nigh-dimensional settings. Compared to the original DiffAN, SciNO achieves 42.7% reduction in order divergence on synthetic datasets and 31.5% reduction on real-world datasets.
Building on the stable estimation of SciNO and motivated by [42], we introduce a method to control autoregressive generative models, e.g. LLMs and Mamba [12], for causal reasoning tasks. While recent studies leverage prior knowledge of L.LSs to improve causal reasoning accuracy [2, 8, 21, 27, 42], these approaches often treat I.LM responses as binary decisions, failing to reflect uncertainty in prior knowledge [5] and remaining structurally incompatible with ordering-basee -ased approaches that sequentially predict leaf nodes. To address these limitations, we propose a probabilistic control method which leverages score-informed statistics inferred by SciNO, guiding autoregressive models to generate more reliable causal reasoning. Applied to real-world datasets, our method reduces order divergence by 64% average, achieving up to 79% reduction over uncontrolled models. Furthermore, our approach can lower the computational complexity of LLM queries from ((V))) to O(V)) contrary to pairwise prompting methods [1. 8, 15. 20, 27]. Notably, these improvements achieved without fine-tuning or intensive prompt engineering. The main contributions of our work summarized as follows:
We propose functional diffusion model framework with neural operators, SciNO,, that enables stable estimation of both the score function and its derivatives, which are crucial statistics for score matching based causal ordering algorithms assuming ANMs.
SciNO enables accurate and scalable causal discovery in high-dimensional settings. We empirically show improved performances over existing score matching based causal ordering methods in memory-effictt way, enhancing its scalability to real-world problems.
2

---

We also propose probabilistic control method that enhances causal reasoning ability of autoregressive generative models. Without requiring fine tuning, instruction tuning, or prompt engineering, it is particularly effective on high-dimensional causal graphs.
2 Preliminaries
In this section, we briefly introduce preliminaries for understanding matching based causal ordering methods-SCCote [33], DiffAN [35], and CaPS [45] which differ in terms of identifiability assumptions, Hessian approximation, and leaf node selection criteria.
Additive Noise Models and Identifiability Given a causal Directed Acyelic Graph (DAG) G. causal order и V V-> - V is non-unique permutation of the set of nodes V/ = {1...., D} such that for any directed edge from node to j in G, holds (i) TTj) [30]. Let T C #K-1 C denote a sequence of causal order corresponding to selected nodes at each step l: = 1,...,,171
1, and write -Mt == V \ Tk. the indices of remaining nodes. identify a causal order from observational data D,
we impose functional assumption known as the Additive Noise Models (ANMs) [31]:
Z4 ((a(z())) + Ein
6 ~ Pi-
SCORE [33] Assuming (2) with Gaussian noise p; (0, of), where each fi is supposed to be nonlinear and twice-continuously differentiable, SCORE proposes an identifiability criterion for inferring a topological order of the causal graph G by investigating the variance of the diagonal Hessian iteratively. Due to [33, 35, Lemma 1]. Var [ (x-xx)] 0 holds if and only if E V Ty is
a leaf node step + SCORE proposes an ordering algorithm as follows:
leaf argmin Var [9, (x-xx))
JEV(
CaPS [45] Under conditions on variance of noises, CaPS proposes an identifiability criterion, where each fi can be either linear or nonlinear, extending its applicability to a wider range of structural functions. Under the ANM assumption [45, Theorem 1], CaPS utilizes the expectation of the Hessian diagonal to identify leaf nodes according to the following criterion:
leaf argmax Ex~D (x-x]]
JEV(EL
DiffAN [35] SCORE and CaPS require re-estimating the diagonal Hessian in (3) and (4) for each step, resulting in computational overhead when the number of nodes [ = D increases. DiffAN replaces the re-estimation procedure by approximating the deciduous score S(x referring the score function over the remaining variables after removing leaf nodes. Due to [35, Theorem 1],
(x)mm) S, + Σ 0.,, (,,,5(())
2,,5((x))
Sx(x)
je VEVTT Tg.
(5)
IERK
Based on (5), DiffAN leverages diffusion models to learn the score function S(x) via score model (t x) near t x 0 for approximating the Hessian diagonal 9 (xx at each step:
S (1,x)
(x,) Rz 0z, (1,, + Σ Oz, Ox Sp (1,
JEV VV Th-
(6)
0, (1,
LET
3 Score-informed Neural Operator
3.1 Score Matching in Function Spaces
While approximation (6) reduces computational cost, numerical instability is implicit in the computing second-order derivatives of score models. Conventional MLPs frequently struggle to accurately estimate derivatives [23], that the curvature information of score models can be deviated from that of the true score function. To address this limitation, we employ Hilbert Diffusion Model (HDM, [26]), which enables functional diffusion modeling Hilbert spaces H rather than applying diffusion models in Euclidean space. Due to Sobolev embedding theorem [17], we consider the Sobolev space H W where denotes the order of differentiabiliiyy and adopt class of neural operators [16] to leamn higher-order derivatives of the score function. See Definition A. and Remark A.2 for details.
3

---

Theorem 3.1 (Pointwise Approximation). Let l E N such that k > 2+ D. K C C X and ε > O, there exists a neural operator So such that
For any compact subset
2.
D
Σ
(g(x))00250 (x) SE, m = k 2
sup XEK Ja|Sm
Details of Theorem 3.1 are discussed in Appendix A.2.2. Since we assume > 2 + D.
2. we have
m that S is least twice continuously differentiable which can approximate the target score function S, with respect to the Hölder norm || 1-1(a==3 3. Theorems A. and A.2 are crucial for the approximation (6) the Hessian diagonal 2 (x) since they supports the numerical stability (35) and guarantees the approximation power (37) for computing the second-order derivatives of score models.
Note that our theory is designed to approximate the target score fumction S. Vx Paata via neural operators. The connection to the marginal score model follows a fundamental procedure in diffusion models. Approximation to marginal score functions with time variable via neural operators in Hilbert space, (t,-)) S(L,) -) for each Hilbert space implies that the distance between (').) and the sequence of functional diffusion models train score models accurately. Based HDM and the above approximation power of neural operators, we introduce the architecture SciNO in the following section, specially designed neural operator which aims to stably approximate the Hessian diagonal (x) to enhance score matching based causal discovery algorithms.
proved in prior work on HDM [26]. Hence the completeness of
will close to
Neural Operator)
3.2 Architecture of SciNO
Lecresble Toue Feroding
Irshioall
Uher
Layer
Fouria Lays
Pipec
ML
Cumplic
Misconwecter
54 Sppcrenction
(al Time-condiiooeed FNO
(b) SciO with LTE
Figure I: Architectures of time-conditioned FNO [26] and SciNO with Learnable Time Encoding.
Table I: Comparison of order divergence in Erdös-Rényi synthetic datasets for D E (10,30, 50, 100) between DiffAN w/ MLP [35], DiffAN w/ SciNO (PE), and DiffAN w/ SciNO (LTE).
Method Dataset
ER(d10)
ER(d30)
ER(d50)
ER(d100)
DiffAN w/ MLP [35]
3.2 + 1.99
26.7 + 8.1
45.5 = 7.51
117.0 + 12.78
DiffAN w/ SciNO (PE)
2.5 + 1.36
19.0 + 6.91
41.4 + 11.33
100.7 + 15.62
DiffAN w/ SciNO (LTE)
2.7 + 1.19
16.9 + 6.12
32.8 + 6.55
86.6 = 12.99
Score-based generative modeling necessitates tine-conditioning to model score functions S x) over continuous time. Following the approach in DDPM [13], HDM implements time-conditioned Fourier Neural Operator (FNO, [24]) by adding projected tensors from Positional Embedding (PE, [43]) module into each Fourier layer, which transforms the input into the spectral domain and performs spectral convolution to capture global patterns while supporting stable estimation of function derivatives. However, while the time-conditioned FNO performs well in functional generative modeling, it shows slight improvements over the based DiffAN in causal orderings tasks. This is due to the instability in estimating the derivatives of the score function ((x) near t = O.
To stably learn the derivatives of the score function, we propose two major modifications to the architecture used in [26]. First, signals in Fourier layers are decomposed into their real and imaginary parts in the spectral domain, allowing more expressive representation of functional information. With the modified Fourier layer, SciNO achieves improved causal ordering performance compared to the original DifTAN, However, this performance gain relatively decreases as the number of variables

---

increases (see SciNO (PE) in Table 1). To scale-up SciNO for nigh-dimensional graphs, rather than assigning fixed embedding vectors discrete positions as PE. we propose learnable encoding function defined over a continuum, Learnable Time Encoding (LTEE) module, that enables the model to jointly learn derivatives in both spatial and temporal directions. The detailed architectural specifications for time-conditioned FNO [26] and SciNO are provided in Appendix A.3.
Learnable Time Encoding in SciNO We first project the time variable f using F F-dimensional learnable weights w Wpoj IR. get a projected vector peoj- Then, we obtain a normalized signal (t) = [cos(t win(t poo))// /2F, which is passed through an MLPUIE : IR - R layer. This embedding similar to [22], which is designed for multi dineessaineni spatial PE. Then, inject the time encoding via elementwise multiplication into both the real and imaginary parts before the input of ML.Pspee spec in Fourier layers (see Figure 1). This continuous encoding aims to learn functional information, including spatiotemporal derivatives, and aligns with the objective of fitting higher-order derivatives by enabling more flexible representation. Consequently, the LTTE module contributes to more stable Jacobian estimation and allows SciNO to outperform baselines more clearly as graph dimensionality increases, highlighting its scalability (see SciNO (LTE) in Table I; additional ablation results are provided in Appendix A.5.1).
0.4
DE W/ MLP DPAA Sc
0.2
Residual
0.0
-0.4
Graph Index
Figure 2: Boxplots of residuals in estimation error of Hessian diagonals corresponding to leaf nodes
across 9 synthetic 2D datasets generated from Structural Equation Models (SEMs).
DIMMAN w/ MLP
DJTSAN w/ ScINO
ERC53001
Order
[Saerraa 1048
(SamEE y 1300
1.2
PMO
1.5
Bext MPP
L4t
0234
5,00
Deal
A002
Best MMD
525
Figure 3: Comparion of correlation between the best MMD and order divergence in Erdös-Rényi synthetic datasets for D E (30, 50} between DiffAN [35] and DiffAN w/ SciNO.
SciNO elaborates more explicit Causal Relationship Figure shows that DiffAN with SciNO consistently achieves residuals closer to zero when approximating Hessian diagonal (x) than the original DiffAN. This result indicates that functional diffusion modeling upon SciNO provides more stable approximation of the score function and its derivatives. We also investigate whether SciNO learns causal relationship explicitly by evaluating the goodness-of-fit between generated samples and data distributions, and analyzing its alignment with causal ordering. Leveraging diffusion models on Erdös-Réény synthetic datasets, we first generate samples by using trained score models implemented by DiffAN with MLP and SciNO. We then measure the Maximum Mean Discrepancy (MMD, [11]) between the generated samples and the data, and estimate the Spearman rank correlation between the best MMD and order divergence for each model. Figure 3 shows significant positive correlation (ER(d30): 0.453, ER(d50): 0.561) between the goodness-offfi fft and causal ordering by SciNO contrary
5

---

to MLP. Overall, these findings indicate that SciNO's stable approximation of the score function and its derivatives leads the model to learn an explicit causal relationship.
SciNO enables memory efficient Hessian Approximation For scalable causal ordering across both linear and nonlinear ANM settings, one can directly replace the score matching in DiffAN and CaPS by SciNO to approximate the Hessian diagonal (xx through (6). Alternatively, one apply a probing strategy built upon pretrained SciNO during the causal ordering, by freezing the Fourier layers and post-optimize MLP rimal to match the Stein gradient estimator [23]. This probing strategy proves useful when applying SciNO to CaPS, which allows fast and memory-efficient adaptation to the Hessian of marginalized density by avoiding the computation of kernel matrix inverse with whole dataset for each iteration, and does not require retraining of the entire network (Figure 4). While CaPS requires ( computational complexity depending on the number of samples N, the probing requires O(TB²N) N) with training epochs T (< N) and mini-batches of size B (< N). We also provide the runtime comparison in Figure A.5 for empirical reference.
3.3 Probabilistic Control of Autoregressive Causal Ordering
While SciNO is stable and scalable method for ordering- based causal discovery, its performance can be further enhanced by incorporating semantic information [2. 8, 21, 27, 42]. Motivated by [42] demonstrating the benefits of using LLMs for ordering tasks, we introduce probabilistic control method that bridges ordering-basee based causal discovery and LLM-based reasoning.
Our approach reinterprets the sequence generation of autoregressive models, e.g. Large Language Models (LLMss) as a causal ordering process. These models generate a sequence due to the conditional distribution PAR (I++1 context) over the next position 2{+11 given sequence Ti:t (II,.... . ). Here, context refers contextual information such as documental guidance (e.g. variable names or descriptions) or domain knowledge that helps to infer a causal structure priori among the variables. From a causal ordering perspective, Fia represents variables whose causal orders are predicted, while I++1 corresponds to candidate variable for the next leaf node. We construct an evidence (likelihood) based on statistics estimated from observed data using SciNO, and propose a method to integrate the prior and the evidence to infer the posterior probability as follows:
P( [I1:,, stat, context) x I PAR |xi: context) x |ssaao (stat//r:tt))))
(8)
Here, stat in the evidence term Pssaoo (1:1+11 denotes a statistic based on (xx), which can be used to determine the leaf node due to the identifiability criterion (3) or (4), proposed in SCORE [33] and CaPS [45], respectively. In this paper, we focus on the criterion (3) without loss of generality,
We propose estimators for computing (saaoo (11:++) based on statistics, average rank and confidence interval of Varx~D [x, of each node i E V. The node set V represents all variables in the causal graph, where |
D is the number of variables. To compute these statistics, multiple number of variance samples is essential and can be obtained by applying the MCMC, bootstrapping, or ensemble technique. We apply Deep Ensembles [19] train multiple independent models with different initializations. Note that due to the probing strategy, one can train multiple models from single pretrained SciNO, enabling memory memory-efficient ensemble strategy. Let !!")) denote / Varga>D (x)) for a given independently trained model m E fl...., M} and node E V.
Rank-based Control Rank-based control uses the statistic of the rank of [arx=
(.rank) SciNO (Varg>D
[x)]
Varx~D [2)(x)] [21., = i), JEV\ m.
(9)
Let ii denote the partial causal order fixed by Z1g-- Motivated from Plackett- Luce model [32], we approximate (9) by computing the average rank F(i) for each node e V across models:
(-F(()
M
pleasik) M (ofmm)
(m) < <a,-m)
( T4+1 = 1) Σ, Ejerr exp (-F())
mm)
(10)
M
m=1
Here r (m) denotes the rank of a: (m)) sorted in ascending order for each model. A lower average rank implies that the node consistently exhibits minimal variance across models. Rank-based estimator (10) provides a scale-invariant and outlier-robust way to compare variances across models.

---

Confidence Interval-based Control Due to the Central Limit Theorem (CLT), the sample mean of (m) )a approximates Gaussian distribution, from which we derive the confidence interval CI(i). Let Cltower (i) and Clupper() (i) denote the lower and the upper bound of the confidence interval. Confidence
interval-based control utilizes the following statistic of the uncertainty in (a. (m())) estimation:
P 1owee (i) Clupper (x) Te+1 = j. argmin Varx~D [x}(x)]
(11)
jEV
We approximate the probability (11) using the following empirical estimator:
m (Cltwer (i) & Clupper (Jum))) Fiies Te41
Σ 1
Clower (1)
Clupper
(mmm)) (U
(12)
M
Here, Jymn denotes the node j of minimal og (m²) and II [:) is the indicator function. Instead of directly comparing variance estimates, estimator (12) can account for the uncertainty across different models.
Algorithm 1: Probabilistic Control of Autoregressive Causal Ordering
Input: Variable list V = Veontext J V -onneex
Output: Causal order æ Initialize f-
2 while V~ * do
Let t (I+(1) (0
T)) and set V-+ V
foreach E V-x do Compute prior PAR(E x context)
Compute evidence Ps))oo());
// by (10) or (12)
i E Veumext then
Compute P(v) < PAR(UiX=, context). PScINo());
// soft supervision
else
10
Compute P(v) f |-== | Pa))
1
// hard supervision
Select leaf node V* argmax P(v) and update क v UU (t*} vEV_-
12 return
Algorithm presents our probabilistic control of autoregressive causal ordering, by using the esti- mators proposed (10) or (12). Each variable belongs to either Veontext, which contains variables with contextual information or domain knowledge, or V-comtext otherwise. The posterior computation differs depending on which set the variable belongs to. For Veontext the algorithm applies soft supervision which updates the prior prediction of autoregressive model by multiplying the evidence term (10) or (12). If variables in v E V V-comeett then the algorithm applies hard supervision by impos- ing a uninformative prior, which selects a leaf node based on the evidence term. These supervisions enable the autoregressive model to operate reliably with mixed levels of contextual information.
4 Empirical Results
We evaluate whether the proposed method improves existing ordering-based methods, DiffAN [35] and CaPS [45]. Additional baseline comparisons are provided in Appendix A.5.2. Given the inferred order, we apply pruning based on feature selection [4] to compare the structural accuracies achieved by the resulting causal graphs across different ordering methods. We also validate causal ordering performances by applying a probabilistic control to LLMs and analyze the resulting changes in probability distributions, Evaluation metries include Order Divergence (OD, [33]) for causal ordering, Structural Hamming Distance (SHD, [29]) and Structural Intervention Distance (SID, [41]) for causal discovery. See Appendix A.4.2 for details on each metric.
4.1 Synthetic Data: Erdös-Rent Random Graphs
Dataset We generate causal DAGs using the Frdös-Rényi (ER) model [9]. For a graph with D nodes, set the expected number of edges to 4D. We sample node values according to nonlinear
7

---

ANM with Gaussian noise, where the functions are generated by a Gaussian Process with Radial Basis Function (RBF) kernel with fixed bandwidth 1. To evaluate model's scalability, we vary the number of nodes across 2, 3, 5, 10, 30, 50, and 100. For each case, we generate 10 random graphs and produce 1,000 samples per graph.
Results Table 2 shows that SciNO improves the performance of DiffAN in most metrics, particularly in terms of OD, and remains scalable in high-ding -ohh-iimensioal settings, notably reducing OD from 117.0 to 86.6 when D 100. The top line graphs in Figure A.4 show the cumulative order divergence, which measures the inaccurate prediction over the leaf nodes. Original DiffAN accumulates OD more rapidly the ordering progresses, while SciNO maintains more accurate prediction, particularly in high-dimensional settings (D& E [30, 50}). The bottom of heatmap in Figure A.4 visualizes the ranked variance of the estimated Hessian diagonal for each variable at ordering steps (10, 20, 30) ER graphs with 30 nodes. SciNO provides more accurate leaf node predictions and assigns consistently lower variances ground-truth leaves, due to the stable estimation of the Hessian diagonal.
Table 2: Comparison of causal discovery metrics((D//HD//SI between DiffAN [35] and DiffAN with SciNO in (a) synthetic datasets and (b) real and semi-synthetic datasets: Physies, Sachs, and BNLearn with nonlinear ANM. Each score is recorded over 10 random graphs for synthetic datasets, and over 10 independent runs for real and semi-synthetic datasets.
DifTAN [35]
DifFAN w/ SciNO (Ours)
Dataset' Metric
OD (L)
SHD (1)
SID
OD (1)
SHD (3)
SID (L)
ER(d2)
0.2 0.4
0.61
0.3
0.0 0.0
0.0 0.0
0.0 0.0
ER(d3)
0.4
0.9
1.1
0.1 0.3 0.3
0.2:0.0 0.6
0.4
ER(d5)
1.1 1.22
1.78
4.32
0.9
1.7 1.95
4.8 4.94
ER(d10)
3.2 1.99
21.9
47.8 10.20
2.7 1.19
20.5 3.5
41.6 - 9.31
ER(d30)
26.7
546.3 82.80
16.9 6.12
88.8 1::+:100 16.5
492.0 91.64
BR(d50)
45.5 + 7.51
184.2 +: 14.63
1690.3 +: 133.55
32.8 6.55
180.0 +: 15.45
1622.0 110.92
ER(d100)
117.0 +: 12.78
463.3 +: 27.64
7562.2 + 435.94
86.6 + 12.99
445.9
32.84
7259.2
653.27
Physics(d7)
0.78
8.6 4.6+-2006
16.8 + 4.73
1.9 0.54
5.8 1.60
8.2 4.62
Sachs(d8)
5.79
22.8
26.1 10.34
5.6 2.65
23.1 3.86
23.3 9.10
MAGIC NIAB(d44)
6.85
89.1 1111--25.0 25.85
227.2 + 157.64
4.1 +: 0.70
72.0 0: 4.52
125.5 + 28.08
BCOL1700d46)
21.8+:4.220
111.0 11.47
684.8 + 80.67
12.8 3.19
90.3 7.13
500.9 58.16
BAGIC-IRRI(d6)
22.0+:3.80
146.9 +: 13.46
321.9 1: 64.64
10.6 1.69
144.5 +: 4.43
254.4 + 39.66
ARTHI500d107)
43.3 15.77
613.3 +: 72.49
2456.0 + 953.08
21.4 = 2.76
515.6 + 15.26
1186.7 = 101.19
Table 3: Comparison of causal discovery (Orriis(OD/SHD/SIDD between CaPS [45] and CaPS with SciNO in BNLear datasets with linear ANM. Each score is recorded over 10 independent runs.
CaPS [45]
CaPS w/ SciNO (Ours)
Dataset Metrie Metric
OD)
SHD (4)
SID
OD (L)
SHD (
SID (J)
MAGIC- (G44)))
36
160
1021
37.5 0.50
165.4 = 0.50
1185.9 5.61
ECOL (706446)
25
134
786
25.4 0.80
141.3 + 4.78
882.4 + 13:02
179
1192
41.6 0.80
183.5 = 2.97
1280.7 - 36.83
ARTHI 50Ld107)
49
406
3051
48.2 0.87
380.0 + 4.63
2799.1 + 73.26
4.2 Real and Semi-Synthetic Data
Dataset We use three datasets to evaluate causal structure learning in real-world graphs: (i) A
Physics commonsense-bant based synthetic dataset comprising 5,000 nonlinear-SEa -EM samples generated over 7-variable water-evaporation DAG [20]. (ii) The Sachs flow-cytometry dataset [34] comprises
7, 7,466 single-cell measurements of 8 proteins after removing cyclic terminal nodes. (iii) The BNLearn collection of four real-world graphs MAGIC-NIAB [38] (44 nodes/ 66 edges), ECOLI70 [36] (46 nodes/ 70 edges), MAGIC-IRRI [37] (64 nodes/ 102 edges), ARTH150 [28] (107 nodes/ 150
edges) each simulated with 10,000 samples under both linear and nonlinear Structural Equation
Models (SEMs). Detailed sampling procedures are provided in the Appendix 1.
Results Consistent with the results on ER synthetic graphs, DiffAN with SciNO outperforms the MI.P-based model across all real datasets (see Table 2). As discussed in Section 3.2, SciNO can be applied to both methods that assume nonlinear relationships, such as DiffAN [35], and methods based on linear structures, such as CaPS [45]. Table 3 compares the performances of the original
8

---

CaPS and CaPS with probed SciNO on BNLearn datasets under linear ANM. SciNO achieves comparable performances across all datasets terms of OD, with improved results observed on higher-dimensione graphs such as MAGIC-IRRI and ARTH150. As shown in Figure 4, while CaPS fails to scale to large datasets due to memory constraints, probed SciNO exhibits significantly lower memory usage.
ARTh 1500d182)
OOMY
OOM
OOM
OOMY
Cans DUSN
CMPS
DETE AUTURD
Somple Size
Stmp a Sime
Sampls Stra
Oumple Stn
Figure 4: Comparison of GPU memory usage of DiffAN [35]. CaPS [45], SciNO during the causal ordering process while scaling up the sample size. SciNO demonstrates memory efficiency while CaPS suffers out-of- -e-oof--mmrry (OOM) errors when the number of samples is larger than 100,000.
Drompt
Currer Capell Order S/
Ulorriend Variahles
by the wint Tikely leat rode in cortal uncy wantatiee in the yives "Drondered coriibllnn en
Select lesf node
Mot
ISI
FAR
Treecs
itoedared Uaciabiens pra2", saace,
True cxued gragh
colieooiinaa vaciables
with dsceat recaciose
Pizee stat, context) x (rees (inz,, context) x PseaNn (stat//ttt
proeaty "huied sulas Iasadiinnn
(haett tal all speed"),
Variable
P(rcealinn comoes)
Fuac( conturtt)
Fsunofstat|/zes>
TSI
0.1125
0.3350
0.1047
Incanteolied SIM Outges
Controlled EIM Outper
by SiN0
SAT
0.6237
0,2872
0.6772
WB X
SAT
V
WS
0.2641
0.3778
0,2180
1.0
LLM Prior vs. Posterior Prob: (ECOLI70)
Posterir
Probabitty
6.6
00
Figure 5: (Top) Overview of the LLM control. The LLM is prompted to output the most likely leaf node among unordered variables, Prior probabilities are multiplied by SciNO's 's data-informed evidence terms to compute posteriors. The variable with the highest posterior probability is selected as the leaf node. (Bottom) Updated posterior probabilities obtained by the SciNO- based control versus LLM prior probabilities of ground-trute truth leaf nodes each causal reasoning step (ECOLI70).
4.3 Controlling LLM for Causal Ordering
Setting We prompt LLMs (Llama-3.1-8B-mn --I--Insttuct [10], GPT-4o [14]) to predict the leaf node among unordered variables, then LLMs output the prior probability FLLM context) in an autogressive manner (see Figure 5). Compared to pairwise prompt methods [1, 8, 15, 20, 27], this approach reduces the number of LLM calls from ((V)2) to O(V)). To mitigate bias due to the variable name length, we apply a length normalization for a E (0.5, 1.0). Full experimental details provided in Appendix 1.
We evaluate the impact of evidence quality by computing the evidence term P((stat with three models: SCORE [33]. DiffAN [35], and SciNO. For this computation, we employ the estimators from (10) and (12). While DiffAN and SciNO are implemented with Deep Ensembles of M = 30 models, SCORE uses bootstrapping with 1000 replicates due to its incompatibility with Deep Ensembles.

---

Table 4: Comparison of the OD in BNLearn datasets between GPT-4o, uncontrolled Llama, and controlled Llama via SciNO. Differences are expressed as percent change relative to the uncontrolled result. Each score is recorded over 10 independent runs.
Method Dataset
MAGIC-NIAB(44)
ECOL170(d46)
MAGIC-IRRIíd64
ARTH110Kd107)
GPT-4o
119
41.1 2.30
20.9 2.30
77.1 + 7.46
Llama-3. L.ama-3.1-8b(L (.oom) w/ crat))l(rankk)
9
32
80
47) 3
10.6
0.49 (766.9%)
9.9 + 0.54 (745.0%)
20.4 + 0.92 74.5%)
w/ coItrol(CI)
2.0 0.00 (777.8%)
9.9
= 1.97 (969.1%)
11.1 = 0.83 (V38.3%)
18.9 = 1.14 (776.9%
Lluma-3, 3.tee,,,1--b..... 1-880(0.5-nomm
14
32
17
80
w/ control(rank)
4.1 0.83 (1 70.75)
10.5
0.92 (767.2%)
9.1
+ 0.30 (V (6.5%)
20.4 0.92 ( 74.5%)
conTrol(CT)
2.9 0.30 (979.3%)
(V69.4%)
10.1
0.83 (V (v406.%)
18.9 1.14 (7764%))
Results Llama-3.1 and GPT-4o, recognized as high- perffoerformacc reasoning LLMs, exhibit signif- icant limitations when confronted with figh-dimensional causal graphs. Figure 5 shows the LLM prior and updated posterior probabilities of ground-truth leaf nodes at each reasoning step on the ECOLI70 dataset. Prior probabilities remain predominantly below 0.2, indicating poor predictive per- formance when faced with a large number of variables. In contrast, updated posterior probabilities in overall steps show that LLM control can reduce inaccurate leaf prediction, thereby enhancing causal reasoning. Consequently. Table 4 shows that the SciNO-based control improves performance across all datasets with more than 40 nodes, achieving up to 75% improvement on ARTH150. Note that the control method also outperforms the results reported in Table 2, which indicates that incorporating semantic information contributes in causal ordering.
We benchmarked our SeiNO-based control against the same algorithm using evidence from DiffAN and SCORE. The full results presented in Table A.7. On high-dimensional datasets, the SeiNO- based control yields a 64% average OD reduction, substantially higher than the 46% reduction from DiffAN. This indicates that the performance gain is attributable to both the novel control method and reliable evidence provided by SciNO. The results in Table 4 are obtained under the setting where all variables are well-described, i.c., Veonext in Algorithm 1. See Appendix A.6.2 for other settings when V Vontext + V.
5 Conclusion and Limitation
We propose SciNO, a functional diffusion model framework that enables accurate causal ordering through stable approximation of the Hessian diagonal. SciNO consistently enhances existing score matching based causal ordering methods in both accuracy and scalability across diverse synthetic and real-world causal graph benchmarks. We further propose a probabilistic control method for supervising autoregressive generative models to generate more reliable causal reasoning by using score-infor score-informed statistics estimated by SciNO. In real- world decision making, the inaccuracy of causal discovery algorithms may cause detrimental consequences. Our method improves the reliability of causal ordering by contributing more trustworthy reasoning and safer downstream applications. While SciNO demonstrates its effectiveness in various causal ordering tasks, several aspects of its applicability remain to be explored. First, the score matching based causal discovery frameworks mainly assume ANMs, whose identifiability is restricted within additive functional relations. Second, SciNO is currently limited to continuous random variables. Future research will explore extending the framework to accommodate discrete or mixed-type datasets, as well as multimodal or structured inputs such as images, text, or temporal sequences.
Acknowledgments and Disclosure of Funding
This work was supported by LG Al Research and also supported by the NRF grant funded by the Korea government (MSIT) (No. RS--2224-004 1008, 50%) and (No. 2RF-2022M33JA 1063595), and also supported by Institute of Information & communications Technology Planning & Fyaluation(IITPP) grant funded by the Korea government
30%) (No. 2S-2002-004366857
ITRC, 20%). (No. 200.-222225-033482. Artificial Intelligence Star Fellowship Support Program, 15%) and (No. 200-S-20019111900779 Artificial Intelligence Graduate School Program (Korea University)).
10

---

References
[1] Alessandro Antonucci, Gregorio Piqué, and Marco Zaffalon. Zero-shot causal graph extrapola- tion from text via LLMs. arXiv preprint arXiv:3312.16670, 2023.
[2] Taiyu Ban, Lyuzhou Chen, Xiangyu Wang, and Huanhuan Chen. From query tools to causal architects: Harnessing large language models for advanced causal discovery from data. CoRR, 2023.
[3] Daniel Bernstein, Basil Saeed, Chandler Squires, and Caroline Uhler. Ordering- based causal structure learning in the presence of latent variables. In International conference on artificial intelligence and statistics, pages 4098-4108. PMLR, 2020.
[4] Peter Bühlmann, Jonas Peters, and Jan Ernest. CAM: Causal additive models, high-dimensional order search and penalized regression. The Annals Statistics, pages 2526-2556, 2014.
[5] Haoang Chi, He Li, Wenjing Yang. Feng Liu, Long Lan, Xiaoguang Ren, Tongliang Liu, and Bo Han. Unveiling causal reasoning in large language models: Reality or mirage? Advances in Neural Information Processing Systems, 11196664-966700 2024.
[6] David Maxwell Chickering. Learning Bayesian networks is NP-complete. Learning from data: Artificial intelligence and statistics pages 121-130, 1996.
[7] David Maxwell Chickering, Optimal structure identification with greedy search. Journal of machine learning research, 3(N0v):507-554, 2002.
[8] Victor-Alexandru Darvariu, Stephen Hailes, and Mirco Musolesi. Large Language Models are Effective Priors for Causal Graph Discovery. arXiv e-prints, 2024.
[9] Paul Erdos and Alfred Renyi. On the evolution random graphs. Publ. Math. Inst. Hung. Acad. Sei., 5(1):17-60, 1960.
[10] Aaron Grattafiori, Abhimanyu Dubey, Abhinav Jauhri, Abhinav Pandey, Abhishek Kadian, Ahmad Al-Dahle, Aiesha Letman, Akhil Mathur, Alan Schelten, Alex Vaughan, et al. The llama
3 herd of models. arXiv preprint arXiv:2447.21173, 2024.
[11] Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. Kernel Two-Sample Test. Journal of Machine Learning Research, 13(25):723-773, 2012.
[12] Albert Gu and Tri Dao. Mamba: Linear-Time Sequence Modeling with Selective State Spaces. In First Conference on Language Modeling, 2024.
[13] Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-6851, 2020.
[14] Aaron Hurst, Adam Lerer, Adam P Goucher, Adam Perelman, Aditya Ramesh, Aidan Clark, AJ Ostrow, Akila Welihinda, Alan Hayes, Alec Radford, et al. Gpt-4o system card. arXiv preprint arXiv:2410.21276 2024.
[15) Emre Kiciman, Robert Ness, Amit Sharma, and Chenhao Tan. Causal Reasoning and Large Language Models: Opening a New Frontier for Causality. Transactions on Machine Learning Research, 2024. Featured Certification.
[16] Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1-97, 2023.
[17] Nicolai Vladimirovich Krylov. Lectures on Elliptic and Parabolic Equations Sobolev Spaces, volume 96. American Mathematical Soc., 2008
[18] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien. Gradient- Based Neural DAG Learning. In International Conference Learning Representations, 2020.

---

[19] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive uncertainty estimation using deep ensembles. Advances in neural information processing systems, 30, 2017.
[20] Chanhui Lee, Juhyeon Kim, YongJun Jeong, Yoonseok Yeom, Juhyun Lyu, Jung- Hee Kim, Sangmin Lee, Sangjun Han, Hyeokjun Choe, Soyeon Park, Woohyung Lim, Kyunghoon Bae, Sunghin Lim, and Sanghack Lee. On Incorporating Prior Knowledge Extracted from Pre-trained language Models into Causal Discovery. In Causality and Large Models NeuulPP 2024, 2024.
[21] Peiwen Li, Xin Wang, Zeyang Zhang. Yuan Meng, Fang Shen, Yue Li, Jialong Wang, Yang Li, and Wenwu Zhu. RealTCD:: temporal causal discovery from interventional data with large language model. In Proceedings of the 33rd ACM International Conference on Information and Knowledge Management, pages 4669-4677, 2024.
[22] Yang Li, Si Si, Gang Li, Cho-Jui Hsieh, and Samy Bengio. Learnable Fourier Features for Multi-Dimensional Spatial Positional Encoding. Advances in Neural Information Processing Systems,
2021.
[23] Yingzhen Li and Richard E Turner. Gradient Estimators for Implicit Models. In International Conference Learning Representations, 2018.
[24] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadeneshelii Kaushik Bhattacharya, Andrew Stuart, Anima Anandkumar, et al. Fourier Neural Operator for Parametric Partial Differential Equations. In International Conference on Learning Representations, 2021.
[25] Jae Hyun Lim, Nikola Kovachki, Ricardo Baptista, Christopher Beckham, Kamyar Azizzade- nesheli, Jean Kossaifi, Vikram Voleti, Jiaming Song, Karsten Kreis, Jan Kautz, et al. Score-based diffusion models in function space. arXiv preprint arXiv:2302.07400 2023.
[26] Sungbin Lim, Eunbi Yoon, Tachyun Byun, Taewon Kang, Seungwoo Kim, Kyungjae Lee, and Sungjoon Choi. Score-based Generative Modeling through Stochastic Evolution Equations in Hilbert Spaces. In Shirth-seventh Conference on Neural Information Processing Systems, 2023.
[27] Stephanie Long, Alexandre Piché, Valentina Zantedeschi, Tibor Schuster, and Alexandre Drouin. Causal Discovery with Language Models as Imperfect Experts. In ICML 2023 Workshop on Structured Probabilistic Inference & Generative Modeling, 2023.
[28] Rainer Opgen-Rhein and Korbinian Strimmer. From correlation to causation networks: simple approximate learning algorithm and its application to high- dimennioini plant gene expression data. BMC systems biology, 1:1-10, 2007.
[29] Jonas Peters and Peter Bühlmann. Structural Intervention Distance for Evaluating Causal Graphs. Neural Computation, 27(3):771-799, 03 2015.
[30] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference: Foun- dations and Learning Algorithms. MIT Press, 2017.
[31] Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. The Journal of Machine Learning Research, 15(1):2009- 2053, 2014.
[32] Robin I. Plackett. The analysis of permutations. Journal of the Royal Statistical Society Series C: Applied Statistics, 2432):193-202, 1975.
[33] Paul Rolland, Volkan Cevher, Matthäus Kleindessner, Chris Russell, Dominik Janzing, Bernhard Schölkopf, and Francesco Locatello. Score Matching Enables Causal Discovery of Nonlinear Additive Noise Models. In International Conference on Machine Learning, pages 18741-18753. PMLR, 2022.
[34] Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal Protein-SS Signaling Networks Derived from Multiparameter Single-Cell Data. Science, 308(5721):523-529 2005.
12

---

[35] Pedro Sanchez. Xiao Liu, Alison Q O'Neil, and Sotirios A Tsaftaris. Diffusion Models for Causal Discovery via Topological Ordering. In The Eleventh International Conference on Learning Representations, 2023.
[36] Juliane Schäfer and Korbinian Strimmer. A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics. Statistical applications in genetics and molecular biology, 4(1), 2005.
[37] M Scutari. Bayesian networks, magic populations and multiple trait prediction. In Invited Talk at the 5th International Conference on Quantitative Genetics (ICQG 2016), 2016.
[38] Marco Scutari, Phil Howell, David J Balding, and Ian Mackay. Multiple quantitative trait analysis using Bayesian networks. Genetics, 198(1):12--337, 2014.
[39] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score- Generative Modeling through Stochastic Differential Equations. In International Conference on Learning Representations, 2021.
[40] Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2000.
[41] Ioannis Tsamardinos, Laura Brown, and Constantin Aliferis. The Max-Min Hill-Climbing Bayesian Network Structure Learning Algorithm. Machine Learning, 65:31-78, 10 2006.
[42] Aniket Vashishtha, Abbavaram Gowtham Reddy, Abhinav Kumar, Saketh Bachu, Vineeth N. Balasubramanian, and Amit Sharma. Causal Order: The Key to Leveraging Imperfect Experts in Causal Inference. In The Thirteenth International Conference on Learning Representations, 2025.
[43] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017.
[44] Yonghui Wu, Mike Schuster, Zhifeng Chen, Quoc V Le, Mohammad Norouzi, Wolfgang Macherey, Maxim Krikun, Yuan Cao, Qin Gao, Klaus Macherey, et al. Google's neural machine translation system: Bridging the gap between human and machine translation. arXiv preprint arXiv:1609.08144 2016.
[45] Zhuopeng Xu, Yujie Li, Cheng Liu, and Ning Gui. Ordering-Basse Based Causal Discovery for Linear and Nonlinear Relations. In The Thirty-eighth Annual Conference on Neural Information Processing Systems, 2024.
[46] Yue Yu, Jie Chen, Tian Gao, and Yu. DAG-GNN: DAG structure learning with graph neural networks. In International conference on machine learning, pages 7154-7163. PMLR, 2019.
13

---

A Supplementary Material
In this supplementary material, we provide the information for understanding SciNO with detailed ex- planations. The supplementary material is organized as follows: Section A. introduces related works in ordering based causal discovery, score matching methods, and LLM-assisted causal discovery. Section A.2 provides theoretical derivations and proofs for preliminaries and the proposed method. Section A.3 explains detailed description of the architecture of SciNO with layer-by-layer formu- lations on how to implement it. Section A.4 describes the experimental setup, datasets, evaluation metrics and implementation details of empirical results in Sections 4.1 and 4.2. Section A.5 presents additional ablation studies and comparisons with causal discovery baselines. Section A.6 presents experimental details in Section 4.3 and additional results on controlling LLMs for causal ordering.
A.1 Related Works
A.1.1 Ordering-basn Causal Discovery
Ordering-basse causal discovery methods are typically categorized into score matching based and non-score matching based approaches. SCORE [33], DiffAN [35], and CaPS [45] are score matching based methods that estimate the Hessian of the log density to identify Ieaf nodes and construct causal order. Non-score matching methods such CAM [4] and GSPo[3] recover the causal order using functional assumptions and conditional independence tests, respectively. In this work, we focus on the score matching based approaches [33, 35, 45], which estimate gradients of the log-density for causal ordering.
A.1.2 Score Matching and Score-based Generative Models
A common approach to estimating the score function of a data distribution is to use a gradient estimator based on the Stein identity [23]. Stein kernel-based methods allow score estimation without an explicit density and have been used in SCORE [33] and CaPS [45]. However, these methods require computing the inverse of kernel matrix, leading to high computational and memory costs. recent work, Diffusion models [13, 39] have emerged as a key method for score estimation in generative modeling. However, the commonly used Denoising Diffusion Model [13, 39] is defined over finite- dimensional spaces and thus struggle to handle functional data such derivatives. Functional diffusion models [25, 26] address this issue by defining the diffusion process over a Hilbert space and learning it directly via score- based modeling. [25, 26] utilize the Neural Operators [16, 24] to learn mappings between functions, allowing the model to capture both the functional form and its derivatives.
A.1.3 Causal Discovery with Large Language Models
Recent studies leverage domain knowledge inferred from LLM as priors or constraints enhance data-driven causal discovery [2, 8, 20, 21, 27, 42]. Most methods adopt pairwise prompts to construct causal graphs via repeated queries on variable pairs [1, 8, 15, 20, 27]. To address inability of LLMs to distinguish between direct and indirect effects, [42] introduces a triplet prompt that queries relationships among three variables simultancously, deriving causal order rather than graph structures. Existing methods treat LLMs a black-box and directly convert generated tokens into binary decisions. This obscures uncertainty and limits the scientific explanatory power of causal inference. [8] partially addresses this by extracting pairwise edge-direction probabilities from the LLM and combining them with mutual information, but this approach remains limited, as it does not directly reflect causal relationships, limiting its explanatory power.
A.2 Proof of Theorems
A.2.1
Derivation of the Hessian diagonal approximation
In this section, we derive approximations (5) and (6), which are suggested by DifTAN [35], for the sake of completeness. Under the ANM assumption (2), the score function S(x) can be derived as:
D
D
Sx(x) = Öz, log P(x) = 0z, log II P(x, Pa(zz)) = 0z, Σ log P(x, Pa(z:)). j V, (13)
14

---

where D denotes the number of nodes [V]. Let 6 2, ft. By the change of variable, it holds that
2z1 log P(a: Pa(x:)) 0z,
Σ log po(;: - fi)
(14)
a log (2 fz)
0f d log p. fi)
(15)
daj
0xy
Oe;
iECh(ay)
When leaf node xę is marginalized out, the resulting score S, (x_1)) for æz Pa(x)) becomes:
0f 0 Olog pı( fi)
S,(x = Sx) + fzz
(16)
del
Let oz, denote the residue term added to the score of xy upon marginalizing out æz;
65 == Oft i lo
(17)
dx;
del
Note that the residue term 5111 vanishes if xy Pa(a1). For a leaf node ny, it holds that:
θ log pi(m fi).
2 loggi( fi). 0(x: fi)
2 log pi(e - fi)
S,(x)
(18)
dri
da
dm
Oe
The partial derivatives of the score (i(x) with respect to Hy for any E V are given by:
82 log (-
0f
02 Sx))
(19)
de? Oz,
In particular, when j = I, this reduces to:
82 log pi (2) - fi)
(,))
(20)
de7
Thus, 5/11 can be rewritten as:
Of log (1 - fi)
8, dx;
(21)
dEL
8² log pm fi).
0f
Ô log prgpr((1 - -i)
0 log pi (I) fi)
(22)
del2 0z,
Oe,
Oe,²
Sx(x)
0., Six)
(23)
0,, Sx
Let Ty C. x denote the of selected leaf nodes up step and TTy V
TT denote the remaining
ones. The deciduous score on the remaining variables -(x- E D- are approximated by:
Sx(x)
mx) = S2(x) + Σ Oz, (ix)
E V
TT.
(24)
(i
IETA
Hence we obtain (5). Finally, the Hessian diagonal (x = (x mx) : V TK) is:
S1(x)
(x(x = 02, Sxx))+ + Σ Ozy
(0,, (x
(25)
0, (x)
IETTU
(1, x)
2 0,, (9, + Σ Oz
Or, S]
jE LEV\TE mE
(26)
0, (1,
IEnx
Thus, we obtain the approximation (6).
Remark A.I. Note that (25) requires the second-order derivatives of S(x) and so does approximation (6) replacing the score function by the score model 8 (t,
15

---

A.2.2 Approximation Power of Neural Operators
Score-based Generative Modeling in Function Spaces For readers who are not familiar with functional diffusion modeling, we briefly introduce time-reversal theory of score-based generative modeling in Hilbert spaces [26].
Consider the following stochastic evolution equations the Hilbert space H::
dx, dx, + G,dWe, dX, B, G,dW, E [0,
(27)
where (X,, X,) is pair of forward and reverse stochastic processes with coefficients (B,, Ge)ee(0.71
and (Br, G.) and W, is an H-valued Wiener process. [26, Theorem 2. I] provides the time- reversal formula which connects pairs (X,, Xx) to have the same marginal distribution:
Bz(u) = (B + Sr-( G, G, GT- te [0,
(28)
Here S, denotes the score operator for each marginal probability measures (14))e(o.r) of the solution the forward equation. Therefore, we can generate samples from Xy Patta and approximate the score operator dara in function spaces by training the score model S (t, and running the reverse equation with replacing the score operator S(t, -) in (28) by the trained score model (t,t,) -) for each t e [0, T]. Therefore, the learning objective of score-based generative modeling function spaces equivalent to the score matching between Stt, -) and S"(t, :) in the Hilbert space H. [25, 26].
Proof of Theorem 3.1 To explain the statement in Theorem 3.1 precisely, let introduce necessary definitions and theorems. We first introduce Sobolev spaces and Sobolev embedding theorem [17], which are main elements for understanding fundamental theory in neural operators [16].
Definition A.1 ((17, Chapter 1.3)). Let W (x') be the Sobolev space with the norm:
Ilflww
Σ lalsk
1lgfll
(29)
with the inner-product:
(f. (f.g)3 =
Σ 1025.000 122
(30)
lalsk
Remark A.2. Note that the inner-product (30) induces the norm which is equivalent to the Sobolev norm (29). Due the Sobolev embedding theorem [17, Chapter 10], it holds that
Wh c C
(31)
where C denotes the Hölder space with the following Hölder-norm:
D
| ( = max sup |a|sm X
[x))]
+ max
||a5|| 9-m
m =
(32)
=
2
Sobolev embedding (31) implies that for any f E W holds
If lffl |Ifll
(33)
where the notation & means the inequality < holds up to constant factors. Therefore, we can approximate a function including derivatives by (33) if we can control the Sobolev norm (29).
(Br, Ge)ee(o.r] as smooth coefficients defined on [0, T] and consider their function class as A TJ), which determines the probability measure μr by simulating the forward A. let us define an operator (t(a)
Let
stochastic evolution equation X, in (27). For given S. which maps smooth coefficients a to the target score function S. = Vx log Piut the data distribution Paataa on the data space X. In this paper, we assume the target score function satisfies S. E W/ for any Lipschitz. domain C X' where l: > 2 + 2. This is natural assumption for many probability densities which are continuously differentiable on Lipschitz subsets of X.
16

---

Theorem A.1. Suppose the map G T[) - (?)) is continuous for any Lipschitz subset C X' such that (a) = S.. Then for any & > O and compact set C A, there exists a neural operator g C(O - W (St) such that
sup a
(g g(a) w3((0) S.
(34)
Furthermore, the neural operator
9 (0) T) ( is discretizatiooninvariin such that
sup acA
|((a]]|
S
sup a'
[g
(35)
where the notation S means the inequality S holds up to constant factors..
Proof. Due to [24, Theorem 8] and Remark A.2, the neural operator G is discreeization-invarian since W (I2) continuously embedded in C(Q)). We apply [24, Theorem 11] to show (34) and (35). Note that Assumption holds since we set A = C(0, T and so does Assumption since we set
W (22. Therefore, for a given exists neural operator G such that
and compact set À A, by [24, Theorem 11], there
sup (a) ((a) ((a)ll S. oε
(36)
Thus, (34) is proved. Also, (35) holds automatically because U W (S2) is a Hilbert space with the inner product (30). The theorem proved.
Remark A.3. Theorem A. implics that there exists a set L-layer neural operators which can approximate the operator G which acts set of smooth coeffcients (B,, Gt)ee j,, Since (B,, G.) Go)eto.2] determines marginal densities (Mt inequality (34) in Theorem A. enhances the approximation to weak derivatives the target score function S, in the Sobolev space. To obtain the desired approximation result in Hölder space, we utilize the Sobolev embedding theorem (see Remark A.2), upon Theorem A.1.
[0, TJ) of stochastic equations (27), the
Now we prove the main result. Let us restate Theorem 3.1 formally.
Theorem A.2. Let ie E such that & > 2 2. For any compact subset C C X and > there exists neural operator 8 such that
sup XE K |m|sm
Σ
02S(x)
S Se, m = kk
D 2
(37)
Proof. Let E >> 0 be given and fix compact set K C Ω. Due to Theorem A.1, there exists neural operator g which maps smooth coeffcients (B,, Gt)rEjo,T] E A W (1)-valuud function,
((a) = S Ω + R,,, with the approximation power (34). Therefore, by (34), it holds that
|s. S9 se/C,
(38)
wh (0)
where the constant is determined by (31). Having C mind, by (31) and (32),
sup K
[g((()) ( S. - -50 c²- f S S. -0
Ss.
(39)
W
hence obtain (37). The theorem is proved.
17

---

A.3
Architectural Details
A.3.1
Architecture SciNO
Learnable Time Encoding
Combrc' complex
Elemane-ooee muliiiication
Exmeetowiee aditiun
Uiner Projection
LTE())
Fourier
Comples
SAkp conmentttio
X
Layor
Lagar
Layer
Real
X
DJ
Bath
Dr. Dats Dimensto
Hi Hidkee Dimeriion
Fofurreer Feature Dimmession
1: Number of Tourier Layees
Figure A.I: Architecture of Score-informed Neural Operator
D-dimensional input data X RDD is passed through the initial MLP iit RRD.
R layer:
x (1) MLPinid (X).
(40)
To explain the intermediate layers in SciNO, let é (1...., L. denote the index of Fourier layers.
Then we transform (() into the spectral domain to get the following C -dimensional tensor:
(0) = FFF(X(0)) O LTE(t), t [0,
(41)
where FFT means the Fast Fourier Transform, @ denotes the element-wise multiplication, and LTE: [0, 1] + R is the Learnable Time Encoding module. The real RR({()) and imaginary (((())
parts are then split to construct an 222 -vaueed tensor x ())( ())( Then we plug the tensor into an M
spec 2R + 2 layer in the spectral domain:
z MLP(1)
spec (x
(42)
Then we combine ( Z
real: inag. to get C -valud tensor ((9) = (2) + i Next, we transform the signal ((9) from the spectral domain the spatial domain with skip connection, which helps alleviate gradient vanishing during training:
( (() + (F))((("))
(43)
where iFFT(-) denotes the inverse Fourier transformation. By repeating the above procedure for E f1...., L}, we get the output of the Fourier layer (l--), Then X (I++1) is passed through the
final MLP rinal I + I layer for final processing:
(t, = MLP final (X ((+1)))
(44)
A.3.2 MLP Modules in SciNO
The MLP modules used in SciNO are organized as follows.
Initial MLP Layer.
The input data X E R is first projected into R through MLP
ILP Dropout,.2 (LayerNorm (LeakyReLU (W)) e H
(45)
where Wintt E HDD
18

---

Spectral MLP Layer. MLP spee is designed to effectively capture differential information in the frequency domain. The input (() E @ is processed follows:
MLP (x (0) BatchNorm [LeakyReLU (wieed spec
+ blpd ))
R22
(46)
The layer is parameterized by W spec and bipecs with shapes 22 and I
Final MLP Layer. The output ((++1) R from the final Fourier layer is mapped to RP via MLP fina:
h1 Lea (1) final X (1++1
(47)
h2 = LeakyReLU
final h1 + 1
(48)
XX (((+))))+ = wx(3)
final h2 + ((3)
MLPanaa(2
(49)
Here, w E RHxH (1) E TR. wx(2) denote the weights and biases of the MLPPinal The dimension S denotes an intermediate hidden dimension between and H.
final RRSH ( E RR.. Wr Dx5 and (
R
fimal
Time Encoding MLP Layer. Given the Fourier feature (() E 2F the time encoding LTE(I) is computed via MLPPTTE:
(2) MLPPITE (+())) GeL.U (Will (1) + LITE I E TR
2
(50)
where W E R22 LLE E RM W dimension of the intermediate layer in the Time Encoding MLP.
RRHM, and bLiE
(2)
RR Here, M denotes the hidden
A.3.3 Architecture Time-Conditioned FNO
We describe the architecture of the time-conditioned Fourier Neural Operator (FNO) [26] for tabular data, following the same setting SciNO for consistent comparison.
The time-conditioned FNO takes an input X E IR_ and first maps H-dimensional hidden representation via a lifting layer:
X Linear
(51)
time encoding TE E R obtained from the scalar diffusion time t using sinusoidal positional encoding [43] defined as mapping PE: [0, 1] - RH and an MI.PPE R. - RHH layer:
TE = MLPPP (PE(t)), t [0. 1].
(52)
This -idmmeenioonal encoding is added to the lifted input:
( < X (1) TE.
(53)
For each Fourier layer E {I,.... LJ. L}. transform the hidden representation X into the spectral domain via FFT. The time encoding, processed by MLP, is then injected into the Fourier coefficients through element- addition to get the C "e--mienniional tensor:
( = (1) + MLP TE TEE).
(54)
Then we apply spectral weights wiO E H
( W w(0) ((0)
(55)
The transformed " -dimensional signal brought back to the TR: -dmeensional spatial domain via inverse FFT:
Z ( = ()).
(56)
19

---

A residual connection is applied:
(0) ( + 0
(57)
Then, an MLP
IR#
- R layer is applied with a residual connection followed by GeLU
spatial
activation:
X GeL + MLP
spatial (X
(58)
After processing through I. layers, the final hidden representation X (L+1) is projected back to the original D-dimension via MI.P final RR.
RR. layer:
(t. X)
ML. Kimll (X (1+)))).
(59)
A.4 Experimental Setting
A.4.1 Datasets
Physics We generate a Physics commonsense- based synthetic dataset using the physics-based causal DAG introduced in [20], which models causal relations among 7 variables related to water evaporation. Each edge weight is sampled from the following uniform distributions:
e~ (--1,-1,--- 0.1) (2((01,1)
(60)
Using the sampled adjacency matrix A, we generate 5,000 samples via the following nonlinear Structural Equation Model (SEM):
x sin A (x 0.5 1)) + A (x 0.5 = ~(011.
(61)
SAT
Figure A.2:
Physics commonsense-baae based synthetic graph with 7 nodes: Rainfall (RNFL), Total Solar Inradiance (TSI), Surface Air Temperature (SAT), Wind Speed (WS), Evaporation Rate (ER), Moisture Content (MC), and Object Weight (Wgt).
Sachs The Sachs dataset [34] contains 7,466 single-cell measurements of proteins and phospho- proteins involved in human immune signaling, collected under various stimulation conditions using flow cytometry. To ensure a DAG structure for evaluating causal order, we exclude terminal nodes (praf, plcg, and PIP2) involved in cycles and use the remaining 8-node graph for experiments.
BNLearn We use four real-world DAGs from the BNLearn repository- MAGIC-NIAB [38] (44 nodes/ 66 edges), ECOL170 [36] (46 nodes/ 70 edges), MAGIC-IRRI [37] (64 nodes/ 102 edges), and ARTH150 [28] (107 nodes/ 1,50 edges)-and generate data based on their causal structures. For the linear case, samples are generated using the graph structures and linear functional relationships provided by BNLearn. For the nonlinear case, we generate data using nonlinear SEMs defined over the same graphs, where each functional relationship is parameterized by Multilayer Perceptron (MLP). generate 10,000 samples per graph in both linear and nonlinear settings for evaluation.
A.4.2 Metrics
Order Divergence. Given a causal order æ and a true adjacency matrix G, the Order Divergence (OD, [33]) is defined as:
Drop (TT, 9)= = Σ Σ
Gij.
(62)
JiRE
20

---

Order divergence counts the number of edges in & that are inconsistent with the order т. A smaller value indicates a closer alignment between и and G.
SHD.
Structural Hamming Distance (SHD, [41]) measures the structural discrepancy between predicted graph G and the ground-truth graph G. It is defined as the total number of edge insertions, deletions, and reversals required to convert 9 into G:
SHD(G, # (1.55) E V2 G and & do not have the same type of edge between i and 5)
(63)
A lower SHD indicates that the predicted graph & is closer to the true causal graph G in terms of structural similarity.
SID. Structural Intervention Distance (SID, (29]) quantifies how similarly predicted graph 9 and the ground-truth graph G infer causal relationships under intervention:
Pa (X.) not a valid adjustment set for } the intervention XX)))(X;) in graph g f
SID # (i.5),
(64)
lower SID indicates that the predicted graph 9 is more causally consistent with the true causal graph 9 under interventions, making SID suitable metric for evaluating the reliability of causal inference.
A.4.3 Implementation Details
For synthetic datasets, set the number of Fourier layers we use L = 1. The Fourier feature dimension in the LTE module (50) is fixed at E = 32. We set the hidden dimensions the model as H max(1024, 5D), and maxx118, 3D) in (45), (46), and (47), where D denotes the number of nodes [V]. All experiments are conducted on single NVIDIA A6000 GPU. Data, checkpoints, and all training configurations are publicly available'.
10, and for real-world datasets,
21

---

A.5
Additional Experiments
A.5.1 Ablation Studies
We conduct a series of ablation studies evaluate the impact key architectural choices.
PE vs. LTE In Table A.I, we provide additional evidence comparing standard PE and LTE on datasets beyond the ER graphs. Consistent with Table 1, LTE yields lower OD, especially as the number of variables increases, confirming its effectiveness in diverse graph settings.
Table A.I: Comparison of OD real and semi-synthetic datasets between PE and LTE in SciNO. Each score is recorded over 10 independent runs.
Method Dataset
Physics(d7)
Sachs(d8)
MAGIC.NIARGAA
FCO1.70044))
MAGIC-IRRlidH
ARTHI
DifFAN SaINO (PE)
2.7 1: 0.46
4.4 à: 1.50
4.0 :: 1.55
30.6 4: 2.42
2.42
16.6 4 1.12
33.9 3.94
DifFAN SeiNO (LTE)
0.54
5.6 + 2.65
4.1 0.70
12.8 128+:3110
10.6 1.09
21.4 + 2.76
To validate our design choice of using multiplicative form of LTE, we compare against the more common additive variant Table A.2. Across ER graph settings, the multiplicative variant consistently achieves lower OD, with the performance gap becoming more pronounced as the number of variables increases.
Additive vs. Multiplicative LTE
Table A.2: Comparison of OD in ER datasets between additive and multiplicative LTE in SciNO. Each score is recorded over random graphs.
Method Dataset
ER(d10)
ER(d30)
ER(d50)
ER(d100)
DiffAN w/ SeiNO (Additive ITTE)
2.8 + 1.74
17.67 + 5.64
38.86 + 11.42
99.2 + 14.27
DifTAN w/ SciNO (Multiplicative LTE)
2.7 + 1.19
16.9 + 6.12
32.8 + 6.55
86.6 + 12.99
DiffAN with LTE We applied LTE to DiffAN to assess whether LTE alone explains the perfor- mance gain. As shown in Table A.3, simply injecting LTE into DiffAN does not yield performance improvement. This suggests that LTE alone is insufficient, and that the performance gain results from combining LTE with our architectural design.
Table A.3: Comparison of OD on synthetic and real datasets between DiffAN with LTE and SciNO. Each score is recorded over 10 random graphs for synthetic datasets, and over 10 independent runs for real and synteett datasets.
Dataset Method
DifTAN w/ MLP (LTE)
DiffAN w/ SciNO (LTE)
(R(22)
0.2 0.4
0,0 0.0 0.1 4: 0.3
ER(d3)
0.1 0.3
ER(d5)
0.8 = 0.75
0.9
ER(d10)
3.9 1.51
2.19
ER(d30)
37.1 7.33
6.12
ER(d50)
62.0 + 11.96
32.8 - 6.55
ER(d100)
120.3 17.12
86.6 12.99
Physics(d7)
5.0
2.0
Sachs(d8)
7.2
5.8 2.93
MAGIC-NIAB(444)
10.6 + 10.6 12.72
3.8 - 0.75
ECOL170d446)
20.4 1 1.85
14.4 1 3.61
MAGIC-IRRI(164)
25.0 1 23.06
10.0 1 0.89
ARTH110((117)
82.0 1 23.28
20.8 2.32
22

---

A.5.2 Comparison with Causal Discovery Baselines
ANM-based Baselines We compare SciNO against ANM-based baselines: CAM [4] and SCORE [33]. As shown in Table A.4, SciNO achieves consistently lower SHD than CAM on the high- dimensional ER datasets and superior performance on the semi-synthetic benchmarks, except for ARTH150. Compared to SCORE, SciNO yields similar performance on low-dimensional graphs and demonstrates a significant advantage on nigh-dimensional settings, achieving higher accuracy than second-order Stein gradient methods. Overall, SciNO demonstrates reliable performance across both synthetic and real datasets, particularly high-dimensional regimes.
Table A.4: Comparison of causal discovery metrics(OD/SHD/SID)) between CAM [4], SCORE [33], and SciNO in (a) synthetic datasets and (b) real and semi-synthetic datasets: Physies, Sachs, and BNLearn with nonlinear ANM. Each score is recorded over 10 random graphs for synthetic datasets, and over 10 independent runs for real and synihetyttttt datasets.
Dataset Metric
CAM
SCORE (33]
DIEAN o/ SAINO (Chrr))
SHD
SID
GD (1)
SHD (()
SID (3)
GD
SHD
SID
(R(dZ)
0.0:1 0.0
0.0
0.3
0.2
0.5
0.0000
0.
0.0
ER(d3)
0.0
0.0
0.:000
0.0
0.3
0.4
0.3
4.04
1
412
1.19
FR((30)
450.2
ER/050)
238.5 17:29
1503.6 201.35
49.3 14.37
180.8 15.99
1527.1 179.6]
32.8
180.0 15:45
1622.0 110.92
ER(/100)
557.2 45.57
6289.9
533.87
120.0 18.49
455.3 20.34
7175.0 :: 435.0
86.6
12.99
445.9
32.84
7259.2
653.27
Physicsdd))
12
11
12
1.9
S. 1.60
8.2 4.66
Sachaid8)
37
41
16
30
2.62:2.05
23.1
23.3 9,10
MAGIC-NIAR(44)
167
184
189
4,1
7:.0+
125.5 78.05
ECOL(70(at))
184
771
3.19
90.3
500.9 à 58.16
MAGIC-IRRliGSA
10.6 1.69
144.5
254.4
ARTHI
51566
1186.7 104.19
Non-ANM Baselines We further compare SciNO with classical Causal Discovery algorithms that do not rely on ANM assumptions, the constraint- based PC [40] and the score-based GES [7]. To ensure fair comparison, convert the predicted graphs of all methods into CPDAGs and report results using Structural Hamming Distance for CPDAGs (SHD-C). As shown in Table A.5, SciNO consistently outperforms both PC and GES across ER datasets of various sizes. This highlights the advantages of our approach against these classical methods: achieves superior accuracy while also being more scalable, overcoming the high computational complexity. Beyond these statistical algorithms, we also consider gradient-basee approaches such as DAG-GNN [46] and GraN-DAG [18], which parameterize causal mechanisms via deep generative models and gradient-based optimization. The results are summarized in Table A.6. While DAG-GNN and GraN-DAG perform competitively on small graphs, their accuracy deteriorates as dimensionality increases. In contrast, SciNO maintains consistently low error even for D 50, clearly demonstrating its scalability.
Table A.5: Comparison of SHD-C between PC [40], GES [7] and SciNO in ER datasets. Each score is recorded over 10 random graphs.
Dataset' Method
PC [40]
GES [7]
DiffAN w/ SciNO (Ours)
(R))
0.1 4 0.3
0.2
0.0
ER(d3)
0.6 - 0.92
1.3
0.0
ER(d5)
4.08:0.64
8.43
1.3
ER(d10)
34.6 2.93
35.2
2..49
(R(330)
108.5 + 15.69
105.7 +
88.3 + 18.11
(R(5)
195.9 + 12.08
180.2 + 12.45
174.3 + 14.48
23

---

Table A.6: Comparison of causal discovery metrics (SHD/SID) between DAG-GNN [46], GraNDAG [18], and SciNO across ER synthetic datasets with varying node sizes. Results are averaged over 10 random graphs.
Dataset' Metric
DAG-GNN [46]
GraN-DAG [18]
DiffAN w/ SciNO (Ours)
SHD (L)
SID (J)
SHD ())
SID (U)
SHD (J)
BR(d2)
HR(d3)
2.1 0.94
2.9
0.3
0.4
0.2
0.4 1.2
ER(d5)
9.3
14.0
0.0
0.0
1.75
4,8 4.99
BR(d10)
34.4 21.62
7.26
20.7
41.6
ER(d30)
108.0 16.37
651.9 60.35
140.9 14.54
614.1
88.8 16.5
4
KR(d50)
192.8 10.81
1900.9
206.39
232.5 16.76
1802. + 40.28
180.0
15.45
1622.0 110.92
ER(d100)
396.0 22.12
7586.4
585.98
232.5 31.81
7623.4 455.73
445.9 32.84
7259.2 = 653.27
DiffAN w/ MLP
Sampled Data
Squared de, log P(x)
Squared d² log P(x)
4:
2
DiffAN w/ SciNO
Squared of o², log P(x)
Squared d log P(x)
2
150
100
22
Vet (00) 0,0458
Mar(S-g00) 0.0037
Figure A.3: Comparison of Hessian diagonal approximation between DiffAN with MLP (top row) and DiffAN with SciNO (bottom row). The left column shows sampled data and its true density. The center and right columns illustrate the squared approximation error of the second derivatives of log P(x) with respect to 2 and x2, computed over a meshgrid.
24

---

Ooder Divergencc
Dresk &f 90
DiNNi N/S0
Driw Al 5010
23
01
12
23
Step
Step
Step
DiffAN w/ MLP
15
ng 000e0 23
DiffAN w/ SciNO
10
10
Veriables
Figure A.4: (Top) Order divergence in ER datasets per each step. The solid lines represent the average order divergence, and the shaded regions indicate the 95% confidence intervals across 10 random graphs with D E {3, 30, 50}. (Bottom) Heatmaps of variance of the Hessian diagonal estimated by DiffAN with MLP and SciNO, respectively. Darker colors indicate lower estimated variance. Red boxes D denote ground-truth causal leaves, and white stars & indicate the variables selected by models. Green boxes denote cases where the model predicts lea node correctly.
MAGIC NAB1044)
BCOL796645)
VAGICIRNINGG
ARTH 15016107)
OOM
OOM
OOM
calk w/1450
12
10 17
12'
Figure A.5: Runtime in seconds of CaPS [45] and CaPS with SciNO across four real datasets (MAGIC-NIAB, ECOLI70, MAGIC-IRRI, and ARTH150) for different sample sizes. Since SciNO involves iterative training. report its runtime per epoch. CaPS fails with OOM 100,000 samples.
A.6 Controlling LLM for Causal Ordering
A.6.1 Experimental Details in Section 4.3
Implementation Details in LLM Control The LLM prompt is designed to select leaf node, variable among the unordered variables that does causally influence any other variable. Each prompt consists of an Unordered Variables es list, a Data Descriptin summarizing the overall data context, and Variable Descriptions providing natural language descriptions of each variable (as illustrated in Figure A.7). Simple input-output examples are provided within the prompt to clarify the task objective and promote consistent output behavior. To avoid formatting inconsistencies such stray whitespace mismatched capitalization, each variable name is prefixed with node_ and presented in the format node_variabe variable..e..eeee __mae.
For each variable name, the LLM generates token-level conditional probabilities for the sequence of tokens &1, f2..... tn that constitute u, where (f1....,t,, The probability of generating the full variable name v is computed using the chain rule:
PL context) = II
PLuM(te tt...,, 4411, context).
(65)
i=1
25

---

Here, context refers to contextual information, including variable descriptions or domain knowledge. In practice, we construct the full imput sequence by appending each variable name to the prompt the form of node_ aadd__vrriee variable able_name, and then feed the entire sequence to the model in a single forward pass. To compute the probability of each variable name, we consider only the tokens following the common prefix node, since tokenization often splits strings like node_v v into node and .v. This allows us to isolate the informative portion of the name and compute joint probabilities efficiently, without repeatedly modifying the prompt.
Length Normalization Since variable names typically consist of multiple tokens with varying length, using raw probabilities introduces a bias against longer names. To address this, we apply length normalization by exponentiating the average of the -iglliilshoods of the tokens composing the variable name [44]. To further mitigate werr-eenlizzttoon of longer names, we normalize the value by dividing it by length" x (0 < a S 1). In this paper, we use a E (0.5,
Pa norm context)
exp
Σ 7
n
log PLIM(4 | (1...,, 1-1, context)) )
(66)
Leaf node selection through LLM control. Llama-3 provides token-level probabilities directly, enabling us to compute the Icaf node probabilities using (66). The leaf node is selected as V* i..., the variable with the highest Length-normalize Tk denotes the set of selected leaf nodes up to step k. For GPT-4o, we treat the variables mentioned in the final response as leaf nodes, since the model's s probability distribution cannot be externally controlled during generation. We set the temperature to its default value of 1.
arg V\,. Pa Pa--nor(())
probability, where
Results Figure A.6 presents the results under the setting where all variables are well-described, i.e., Veontext = V. The top plot compares the LLM prior probabilities with the posterior probabili- uies-updated via control-of the ground-truth leaf nodes at each step of the causal reasoning process. A substantial increase in the posterior probability is consistently observed across nigh-dimensional datasets such as ECOLI70, MAGIC-NIAB, MAGIC-IRRI, and ARTH150. The effect is especially pronounced in the early stages of the causal ordering process, where the true causal structure must be inferred from large pool of candidate nodes. Furthermore, the bottom plot illustrates several cases in which the posterior successfully identifies the true leaf node, even when SciNO alone fails to do so. These results indicate that the LLM prior can be complemented through integration with SciNO's data-informed evidence, while the semantic knowledge embedded in the LLM can also compensate for incomplete uncertain statistical inference.
Table A.7 presents the complete benchmark results of the LLM control method using evidence from SCORE [33], DiffAN [35], and SciNO. SciNO shows clear superiority high-dimensional datasets, achieving an average 64% OD reduction, compared to 46% from DiffAN. While SCORE-based control yields improvements
some datasets, it fails surpass its data-only baseline on ECOLI70 and ARTH150, as shown in Table A.4. Overall, these findings indicate that the synergy between LLM priors and data-driven evidence is most robust with stable, diffusion-bart based methods like SciNO.
A.6.2 Additional Experiments on LLM Control
In many domains such as biomedicine, healthcare, and finance, variable names are deliberately abstracted and metadata is often scarce and uncertain. To account for this practical challenge, we conduct an additional experiment in which all variable names are masked and the proportion of variables accompanied by descriptions is systematically varied.
Experimental Setting We consider an experiment to simulate practical conditions involving a mixture of Voontexe and V-gontext # 0. where Veomext refers to variables that contain contextual or domain-specific information, whereas V cuntext includes variables without such information. Variable names are replaced with randomly generated four-letter alphabetic strings (e.g., 'zntr', 'lgvp', 'aqsm" "), and descriptions are provided for 10%, 40%, or 70% of the variables. Each experiment is repeated three times, with randomization applied to both the selection of context-give -gitext-givennn variables and their input order. We use semi-synthetic datasets from BNLearn: MAGIC-NIAB (7/44 variables have descriptions), ECOLI70 (41/46), MAGIC-IRRI (10/64), and ARTH150 (105/107). Due to the limited availability of variable descriptions in MAGIC-NIAB and MAGIC-IRRI, we restrict experiments on these two datasets to the 10% description setting.
26

---

Table A.7: Comparison of OD in real and semi synthsytthetic datasets between GPT-4o, uncontrolled Llama, and controlled Llama via SCORE [33], DiffAN [35] and SciNO. Differences are expressed as percent change relative to the uncontrolled result. Each score is recorded over 10 independent runs.
Mothed
Dataset
Physiesdd))
Sacha(d8)
MAGICNIARRH4)
ECOL17044))
MAGIC-IRRId64))
ARTHI500117))
GPI-40
0.54
4.4 1.11
11.9
2,30
7.46
Liama 3.1 sbb norm)
9
32
18
SO
controlkannk
00 0.00
9.0 0.00
1. 0.00000
14.7 + 0.46
11.6 0.440
42.9 + 1/1
SCORE [53)
control(CT)
=0:0000
600 0.00
0.46
0.22
Llama-3, --b(()
17
control''s anil)
control(C)
0.00
7.0
0.00
18.1 0.30
6.22
Llama-3. 1-3b((-normm
controli anik)
21.1
DIRAN [S5]
controXC1)
0:40
5.7
5.3-1: 1:10
18.3 1.42
14.3 1:27
67.9 9.50
Llama-311 3.1 86(0 norrm) w/ controkrand w/ controlf CI)
4.0 0.00
4.4 1.20
1,7 0:64
20.5 1.80
7.5 + 0.81
21.5 + 2.54
4. 0.00
5.0 0,00
8.2 1.72
18.2 1: 1,33
13.3 4: 1.85
67.9 1: 9.50
Lham- 3.1 SbC norm)
18
w/ controkrania
1.3 0.64
5.5 1.28
3.2 0.:0
10.6 1: 0.40
9.9
1: 0.54
20-4 0.02
SeINO (Curs)
controlf CI)
1320000
3.7 0.04
2.0 0.00
9.9 1,97
11.1 à 0.88
18,9 1: 1.14
Liama 3.1 Slb((lbs-aaaeeee..e
32 10.5 à 0.92
17
controkrannd contro(CT)
0.04
3.8 0.58
20.4
0.330
1.83
Posterior
LLM Prior vS. Postorier Prob. OMAAREN (NNU)
LLM Prior VoL Pasterier Preb. (ARTHINS)
SHP
3008
Probabilities Recovery
MAG
ECOL
0.75
MAGICIRRI
prob
Figure A.6: (Top) Comparison L.MM prior and posterior probabilities for the true leaf nodes across four datasets (ECOLI70, MAGIC-NIAB, MAGIC-IRRI, and ARTH150). Posterior probabilities computed by multiplying the LLM prior with SciNO computed probabilities. (Bottom) Heatmap of SciNO probabilities for the SciNO-sclecte
nodes each ordering step across datasets. Red boxes C denote steps where SciNO alone fails to select the correct leaf node, but the correct node is successfully recovered via integration with LLM priors.
Control with Hard and Soft Supervision Under the above experimental conditions, we follow the causal ordering procedure described in Algorithm The control mechanism is formalized as:
P( [aa, stat, context) x P F1.+ context) Pssano(stat [11:+1).
(67)
Let n. denote the partial causal order determined by 2iss sequence of nodes selected in the previous steps. The evidence term PsaNo (stat ( is estimated using two approaches introduced Section 3.3. The first estimator relies on the average rank F(i) for each node i E V Fe across models:
(-F(i)(i))
M
S(rank) mo-M (o!m) < <0,m0
( F1a. J1+1 = =,==,, EjEVio exp (-F(1))² F(I) = M
Σ
(m)
(68)
m=1
Here, r (m) is the rank of node in model m, based on the variance of its Hessian diagonal. The second estimator is based on confidence intervals:
M
m (Cloower (i) Clupper (mm)) ) 1 Ez+1 = =
Σ
1 Chower (1) S Clupper [mm) (69)
M
m
27

---

where {(m] is the node with the smallest variance of the Hessian diagonal in model m.
Depending on whether variable descriptions are available, we apply two supervision strategies. For variables without descriptions V-context. we use hard supervision, relying exclusively on the evidence term estimated by (10) or (12):
1
P(u)
IV (ss =))
(70)
For variables with descriptions Veomtent, we apply soft supervision by multiplying the LLM-derived prior (65) to the evidence term (10) or (12):
P(e) x Po (T+41 = v ,e, context) (stat( my 2++1 =)-
(71)
To adapt the soft supervision signal, one can use a temperature parameter 2
P. (,, = v) = (stato NT, J++1 = u)
(72)
P. Po- (statt/ 2) softmax(P, (stat) IT, J++1 u)).
(73)
We choose T E 2, 4, 5} for experiments. When 7 > 5, the evidence term grows exponentially and overwhelms the prior, leading to predictions that are nearly indistinguishable from using the data statistic alone. Applying the above 7-temperature to the soft supervision, the posterior predictive term computed by:
[e, P(v) x Pa Pa-nomm( (e, context). (stat/mmt e), V VE NEVI Tt.
(74)
Then we select the variable with the highest posterior probability as the next leaf node:
argmax P(v).
(75)
vEVi
Results Table A.8 reports the order divergence under the setting where Veontett A V, across varying levels of description proportions and values. Compared the results in Table performance without control reveals that the LLM's s ability to infer causality deteriorates as descriptive content decreases and variable names are masked. Across all datasets, increasing 7 from 0 consistently improves performance, even under low description conditions (c.g.. 10%), indicating robustness in the combination of LLM prior and SciNO evidence. When the description is limited, performance tends to plateau after т = 1, suggesting that the contribution from LLM priors has been fully exploited, leaving SciNO as the dominant signal. Conversely, with a higher description ratio, the effect of soft supervision becomes more pronounced as varies, allowing the model to better leverage both semantic and structural signals. Interestingly, a higher T does not always lead to better results. This highlights the importance of finding an optimal balance between the LLM'' s semantic knowledge and the evidence from SeiNO. We compare various combinations of description proportions and values. achieves the highest performance
Among them, the configuration using 70% descriptions with on the ARTH150 dataset. outperforms all other tested combinations, including those with lower description ratios, different values, and even the use of SciNO alone. This suggests that when richer semantic information is available, the integration with SciNO becomes especially effective. In conclusion, these results demonstrate that under limited description conditions, combining LLM semantic knowledge with data-driven causal inference can lead to improved performance.
28

---

Table A.8: Comparison of OD across semi-synthetic datasets under the control with hard and soft supervision setting with varying description proportions (70%, 40%, 10%) and temperature parameters 7 {0. 1, 2, 3, 4, 5). The lowest OD for each dataset and description proportion is highlighted in bold. Each score is recorded over 3 independent runs.
Proportion
r-2
MAGIC-NIBB 10%
1ane--3
1-bh005--omm)) contel(rankk
37
36.0.5 control(CE)
223
5.
4.0
4.3
4.3
1luu-33 1-5c:l----- L-noe)) cantnol cask)
303:221
35.7 = 2.9
87 0.06
8.3 0:0.5
8.0 6.0
8.
8.3 83
70%
Uhma-33 3.1 1--.0.0000mm.... vi contreitrant)
1hma33 (noem) conmmi())
303
24.0 2.0
0.0
0.8
2.3
8.3.8. 0.6.
12.3
12.0 1.0
833
11.7
1mm-3 1-5000 (-mmm)))))))) wi (CCtrol(C)))
2.9
3.8
Llama-3.1-80 3.1ma-3.-9b:1-oomma control
2t0005550mmm... controlizank)
367..997
0.6
0.5
8.3
8.3
11.7
contoo)
34
0.0
5.0
0:00
10%
Lluma 30/0.3 contriltrank)
Liuma
conmon C2)
2433..3.
0.0
0.6
8.0
10.7
0.0
1lam-33 1-30/00 S-anmn) wi cistrol(C()
3
12.7 3.8
11.0
10.7 0.6
10.7 1: 1.5
10.7 1.5
11.0 1.0
Lima 3.1 3b11 noem) vanko
297.1.1
11.3 122
0.0.8
9.3
9:
9.2 0.6
11am--311- 5-amm)) contrelteank
9.7
Llama-3.1 (1-8b:1-0oom) control
297
9.7
9.7
costeoliC))
Llama 3.1 8b:1 noem) contm
2.3
80.7 8: 122
20.3 1.2
20.2 =: 0.6
21.3 1.2
21.0 1.0
21.0 1.0
70%
Liama 3.1 86/0.5 somn) wr commoltranki
Liama 3.1 9611 noem) con
750
6:33.:33
77.3 1 55
73.3 = K.6
2000
19.7 = 0.6
1.7
19.7
19.7
1 0.6
1.2
20.0 2: 0.0
19.0 1.7
20.0
20.5 2 0.6
2.0
20.3 00.3-10 0.6
19.3
1.mma33 3.1 1-3600---aommm w/ (CI)
64.3 3.5
71.3 m K.1
17.7
0.6
12.0
3: 1.0
16.7
1: 2.1
16.3
1: 1.2
16.3
1:33:11 1.2
Litma 3.0 8b1l nomm) conmon caak)
16.0
$1.3 12.9
20.0 1.0
20.0 0.0
20,0 &: 0.0
2000 00
ARTH150
40%
Limma 36/0. 3-nom)) contreltrank)
Lima-3 (--o--) c1)
57.0
77.2
20:3 0: 1.2
20.0 11 0.0
20:0
20:0
0.0
1.7
200
200
Jama- 1-5600 5-tmm))) wi contrel(CT))
8777
59.3
21.3
20.7
20.7
1.2
1977
moem) controkzask)
3..1-500..-oom controltrank)
1933-5
1ano-3
5-amm)) (CC)
19.3-1
29

---

B
Miscellaneous
You Al assistant tasked with identifying the most likely leaf node in causal structure. A leaf node is variable that does not cause any other variables in the unordered variables set. Your goal is to determine the best leaf node among Unordered Variables using the given information.
Selection Criteria:
A leaf node does not act as cause for any other variable in Unordered Variables. If multiple candidates exist, select the one that influenced by others but does not influence any other variable Unordered Variables.
Important Formatting Rules:
Respond only with the variable name of the selected leaf node.
not include any punctuation, reasoning, quotes, formatting.
Leaf Node must be exactly variable
plain text, matching from the Unordered Variables
list.
Do not include any additional text before after the variable
Example 1:
Input: :
"node_Temperatuue Temperature") Data Description: The dataset contains weather data recorded hourly with multiple atmospheric variables. Variable Descriptions:
Unordered Variables: |"node__oee CloudCover", "node Hurnidity,, "node_ Pressure", "node,
sky covered by clouds.". "node Humidity": "The amount water vapor the air.". mnode_Pressuuue "The atmospheric pressure at given location.", ". "node_ Temppeaturrrr, "The measure of how hot cold the air is."
"node. CloudCover": "The fraction
Output: node_ Tend___uueerturre
Input: :
Unordered Variables: (unordered variables)
Data Description: (data description)
Variable Descriptions: [variable_description description)
Output:
Figure A.7: Prompting for Leaf Node Selection
30

---

Physics (u(=alpha)=0.5, Control: rank-based) )
== Causal Ordering Step t=1 ===
Prior: ('Weight of object': 0.303, 'Moisture Content object': 0.205, 'Rainfall': 0.144, .... 'Total Solar
Irradiance:: 0.071, 'Surface Air Temperature'; 0.068)
SciNO Likellhood: ("Weight of objject:: 0.635, 'Rainfall': 0.233, 'Molsture Content of object': 0.086, 'Wind
Speed': 0.003, 'Total Solar Irradiance': 0.002)
Updated Posterlor: ('Weight of object': 0.776, 'Rainfall:: 0.13,, 'Moisture Content of object': 0.071, 'Wind
Speed': 0.001,
Solar Irradiance'. 0.001)
Current Causal Order z; (Weight object)
Causal Ordering Step t=2
LLM Prior:
(MMosstrre Content of object': 0.288, 'Rainfall:: 0.220, 'Rate of Evaporation's 0.219. 'Total Solar
Irradiance:: 0.086, 'Surface Tempersture': 0.086)
SciNO Likelihood: 'Rainfall': 0.644, 'Moisture Content of object'; 0.214, 'Surface Air Temperature', 0.084,
'Wind Speed': 0.010, 'Total Solar Irradiance:: 0.005)
Updated Posterior: ('Rainfall: 0.841, 'Moisture Content object': 0.278, 'Rate of Evaporation': 0.040,
'Wind Speed': 0.005, 'Total Solar Irradiance:: 0.002)
Current Causal Order nz [Weight object, Rainfall]
(.....
Causal Ordering Step t=5
LLM Prior: ('Wind Speed': 0.464, Totta Solar Irradiance', : 0,300, 'Surface Air Temperature": 0.235)
Likelihood: ('Surface Air Temperature':
Speed': 'Total Solar Irradiance's 0.104)
Updated Posterior: ('Surface Air Temperature' 0.545, 'Wind Speed': 0.346, 'Total Solar Irradiance': 0.107)
Current Causal Order æzi (Weight of object, Rainfall, Surface Air Temperature]
== Causal Ordering Step t=6 ===
('Wind Speed': 0.584, 'Total Solar Irradiance 0.415)
SciNO Likelihood: ('Wind Speed': 0,704, Tootal Solar Irradiance:: 0.295)
Updated Posterlor: ('Wind Speed': 0.769, 'Total Solar Irradiance': 0.230)
Current Order nı [Weight object, Rainfall, ..... Speed]
Complete Causal Order: [Total Solar Irradiance, Wind Speed, Rainfall, Weight of object] Order Divergence:
Figure A.8: LLM control log on the Physics dataset.
31

---

Sachs (((=alpha)=100, Control: Cl-based
Causal Ordering Step
LLM Prior: ('pakts473': 0.195, p44/42:: O. 166, 'pjnk': 0.163, PPIP3:: 0.077, 'P38: 0.071)
ScINO Likellhood: (pmmek:: 0.125, 'PIP3" 0.125, 'p44/42", 0.12,, ... 'P38': 0.125, pjnkk: 0.12)
Updated Posterior: (pakts473': 0.195, pp44/22: 0.166, pink": 0.163, .... 'PIP3': 0.077, 'P38:: 0.071)
Current Causal Order zj: [pakts473]
Causal Ordering Step t=2
LLM Prior: ("p44/22:: 0.213, PKC' 0.210, 'ojnkk: 0.204, 'PIP3: 0.089, P38' 0.074)
SciNO Likelihood: ('pmek': 0. 142, 'PIP3" 042, 0p44/42": 0.143, 'P38:: 0.143, 'pink: 0.12
Updated Posterior: (p44/42: 0.213, 'PKC:: 0.210. 'pjnk': 0.204, ... 'PIP3:: 0.089, 'P88: 0.074)
Current Causal Order œzi [pakts473, p44/42)
(.....
Causal Ordering Step t=6
LLM Prior: ('P8': 0.496, 'PKA': 0.286, PIP3': 0.217)
SciNO Likelihood: (PPIP3: 0.416, P38' 0.416, 'PKA': 0.166)
Updated Posterior: (P38: 0.599, PPIP3 0.262, 'KAA: 0.138)
Current Causal Order æg: [pakts473, p44/42, P38]
Causal Ordering Step t=7
LLM Prior: ("PI33: 0.605, PKA' 0.394)
ScINO Likellhood: ("PI33: 0.500, 'PKA:: 0.500)
Updated Posterior: (PIP3:: 0.605, 'PKAA: 0.394)
Current Causal Order πt; [pakts473, p44/42, P38, PIP3]
== Final Results ==
Complete Causal Order: [PKA, PIP3, .... p44/42, pakts473] Order Divergence: 3
Figure A.9: LLM control log the Sachs dataset.
32

---

MAGIC-NIAB (<(=alpha)=0.5, Control: rank-based
Causal Ordering Step t=1
LLM Prior: GLASS 0.176, 'YLD:: 0.142, HT: 0.118. 'G1750:: 0.019, '11766; 0.016)
SciNO Likellhood: (YYR,IELL) FIELD': 0.328, 'FT: 0.296, 'YLD: 0.297, 'G1294": 0.000, 'G1276": 0.000) Updated Posterlor: ('YLD: 0.424, 'FT', 0.311, YR..FIED": 0.216, G1276: 0.000, 'G1294': 0.000)
Current Causal Order xy: [YLD]
Ordering Step
LLM Prior: 0.161, 'HT:: 0.149, 'G43;; 0.108, 'G1750': 0.022, 'G800': 0.015)
SciNO Likelihood: FIELD': 0.443, FT: 0.414, 'FUS' 0.102.
'G1276': 0.000, 1G1263: 0.000)
Updated Posterior: (FT': 0.712, YR.FIELD:: 0.171. FUS: 0.111.
G1294: 0.000, G61263:0.000)
æy: IYLD, FT]
(.....)
Causal Ordering Step t=42
LLM Prior: ('G418': 0.359, 'G1750: 0.347, G11217 0.292)
SciNO Likelihood: ('G418: o 685, 'G1750' 0.183, 1G1217 0211) Updated Posterior: (G418) 0.706, 111500 1922 'G1217:00101
Current Causal Order x42: [YLD, FT, G418)
Causal Ordering Step t=43
LLM Prior: ('G1750': 0.527, G1217: 0.472)
SciNO Likellhood: (/12217: 0.860, 'G1750' 0.339)
Updated Posterlor: (G1217' 0.635, 'G17500 0.364)
Current Causal Order 4 [YLD, FT, .... G418, G1217]
Final Results
Complete Causal Order: [G1750, G1217, .... G38, FT, YLD] Order Divergence:
Figure A.10: LLM control log on the MAGIC-NIAB dataset.
33

---

ECOLI70 d(=alpha)=0.5, Control: rank-based
=a Causal Ordering Step t=1 mm
LLM ('atpG': 'asnA': 0.054, 'atpD': 0.054, 'b1191": 'yedE": 0.027;
SciNO Likelihood: ('gbbD: 0,742, 'tnaA': 0.193, 'nmpC': 0.030, 'sucD': 0,000, 'yecO': 0.000) Updated Posterior: ('ygbD: 0,476, 'tnaA': 0.440, 'atpG': 0.058,
'yeOO;; 0.000, 'sucD': 0.000)
Current Causal Order x [ygbD]
= Causal Ordering Step t=2 ===
Prior: ('atpG': 0.060, 'yhel': 0.048, 'asnA': 0.046, ... 'b1583': 0.033, 'b1 11111 0.032)
SciNO Likelihood: ('cspA:: 0.720, 'tnaA': 0.274, "pspA': 0,003, 'dnaG': 0.000, 'yecO': 0.000)
Updated Posterior: ( ('cspA:: 0.579, 'tnaA': 0.418, 'pspA': 0.001, .... 'dnaG': 0.000, 'icdA': 0.000)
Current Causal Order xy; lygbD, cspA]
(.....)
Causal Ordering Step t=44
LLM Prior: ('b1191:: 0.359, 'csp@': 0.352, 'eutG:: 0.287)
SciNO Likelihood: (bb 191': 0,669, 'eutG': 0.214, 'ospG:: 0.115)
Updated Posterior: ('b1 1911 0.701, 'eutG': 0.979, 'espG': 0.118)
Current Causal Order x44 lygbD, cspA, b1191]
Causal Ordering Step t=45
LLM Prior: ('eutG': 0.540, 'cspG': 0.459)
SciNO Likelihood: ('eutG': 0.673, 'cspG': 0.326)
Updated Posterior: ('eutG': 0.708, 'cspG'; 0.291)
Current Causal Order Ags! [ygbD, cspA, .... b1191, eutG]
Final Results
Complete Causal Order: (cspG, eutG, cspA, ygbD]
Order Divergence: 10
Figure A.11: LLM control log the ECOLI70 dataset.
34

---

MAGIC-IRRI ( (d(=alpha)=0055, Control: rank-based )
Causal Ordering Step t=1
LLM Prior: I'YLD': 0.238, 'HT'; 0.070, 'BROWN': 0.089, G3092': 0.023, 'AMY': 0.021)
SciNO Likelihood: ('CHALK:; 0.664, 'GW; 0.158, 'YLD': 0.134, 'FT': 0.000, 'AMY:: 0.000)
Updated Posterior: ('YLD': 0.518, 'CHALK": 0.445, 'BROWN': 0.026, .... 'AMY': 0.000, 'FT': 0.000)
Current Causal Order x;: [YLD]
Causal Ordering Step t=2
LLM Prior: ('HT': 0.119, 'CHALK: 0.099, 'BROWN': 0.093, 'G3927": 0.030, 'G3098': 0.024). Likelihood: ('CHALK':
0,232, 'BROWN': 0.068, ... 'AMY': 0.000, 'GTEMP": 0.000)
Updated Posterior: ('CHALK's 0.866, 'BROWN:: 0.084, (GW: 0.045,
'FT: 0.000, 'HT': 0.000)
Current Causal Order xy: (YLD, CHALK]
(.....
Causal Ordering Step t=62
LLM Prior: ('G3925': 0.379, 'G3823': 0.360, 'G4158': 0.260)
SciNO Likellhood: ('G4156': 0.439, 'G3925': 0.284, 'G3823': 0.275)
Updated Posterior: ('G4158': 0.355,
0.335, 0.308)
Current Causal Order Mizz (YLD, CHALK, G4156]
Causal Ordering Step t=63
LLM Prior: ('G3823:: 0.690, 'G3925': 0.309)
SciNO Likelihood: ('G3925'; 0.689, 'G3823': 0.310)
Updated Posterior: (G3829': 0.500, 'G3925': 0.499)
Current Causal Order (YLD, CHALK, G4156, G3823]
=as Final Results waa
Complete Causal Order: (G3925, G3823, Order Divergence:
CHALK, YLD]
Figure A.12: LLM control log on the MAGIC-IRRI dataset.
35

---

ARTH150 (<(=alpha)=0.5, Control: Cl-based
=== Causal Ordering Step t=1 ==
Prior: ('519': 0.066,
'93': .... 0.025, 0.024)
SciNO Likelihood: ('496'; 0,500, 0.500, 4': 0.000,
'47' 0.000, '61'; 0.000]
Updated Posterior: ('496: 0.662, '677': 0.337, 0.000, '47: 0.000, '61': 0.000)
Current Causal Order mi: (496)
Causal Ordering Step t=2
LLM Prior: ('519': 0.073, '78': 0.047, '93': 0.040, '783' 0.026, '539': 0.025)
SciNO Likelihood: 1'677' 1.000, '4': 0.000, '8:: 0.000, '61:: '63': 0.000)
Updated Posterior: ('677: 1.000, 0.000, '8:: 0.000 ... 0.0000000, 0.000
Current Causal Order n [496, 677]
(....)
Causal Ordering Step t=105
LLM Prior: ('738': 0.378, '539': 0.311. '783': 0.309)
SciNO Likelihood: :339 0.333, '788: 0.333, 783': 0.333)
Updated Posterior: ('738:: 0.378, '539': 0.311, '783': 0.309)
Current Causal Order iuus [496, 677, 738]
Causal Ordering Step t=106
LLM Prior: ('783': 0.533, '539": 0.466)
SciNO Likelihood: ('539': 0.987, '73' 0.032)
Updated Posterlor: ('539: 0.963, '73': 0.036)
Current Causal Order [496, 677, 738, 539]
=as Final Results mam
Complete Causal Order: [783, 539, ..., 677,496) 496)
Order Divergence: 18
Figure 13: LLM control log on the ARTH150 datasct.
36

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The claims presented in the abstract and introduction clearly and accurately reflect the contributions and scope of the paper.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Limitations are discussed in Section 5.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model sell-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
37

---

Section A.2 provides the full set of assumptions and complete proofs for all theoretical results.
Justification:
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proof's provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided not)?
Question:
Answer: [Yes]
Justification: Section A.3 describes the full model architecture. Section A.4 and A.6.1 provide all experimental settings and dataset details required for reproduction. In addition. the code repository referenced in Section A.4.3 includes all source code and pretrained checkpoints.
Guidelines:
The answer NA means that the paper does not include experiments. If the paper includes experiments, a No answer to this question will not 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 make their results reproducible 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 a specific model and empirical evaluation, 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 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 large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is a new model (e.g., a large language model), then there should either be way access this model for reproducing the results or way to reproduce the model (e.g., with an open- open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in authors are welcome to describe the particular way they provide for reproducibility. the case of closed-source models, 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 have some path reproducing or verifying the results.
cases, in which
38

---

5. Open access to data and code
Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
The code repository linked in Section A.4.3 includes the full source code and pretrained checkpoints, along with data access scripts, preprocessing routines, and environment specifications necessary to reproduce all main experimental results.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines
cc/
public/guides aubllc//gullessDdddetttmmsssssss
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
//nips. wwwppsccppuulliiip cttppblic/eacccccccccccccccceccc ceriic/guiiee /ewedeeeeeoeemiieieece ioriissionn weeiissssonPolii onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Answer: [Yes]
Justification:Al experimental settings, including data splits, hyperparameters, and training procedures, clearly specified in the code repository linked in Section A.4.3, enabling full understanding and reproduction of the reported results.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Ys]]
Justification: The top of Figure A.4 presents 95% confidence intervals computed over multiple randomly generated graphs, providing a clear measure of statistical significance for the reported performance.
Guidelines:
The answer NA means that the paper does not include experiments.
39

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the enror bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
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.
is OK to report I-sigma error bars, one should state it. The authors should preferably report a 2-sigma error bar than state that they have 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 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.
8. 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 reproduce the experiments?
Answer: [Yes]
Justification: Section A.4.3 provides detailed information on compute resources. This allows for faithful reproduction of the experiments under similar conditions.
Guidelines:
The answer NA means that the paper does 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 failed experiments that didn't make it into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: http://eerrpp.eeeee httppuiic/ccissuiieeeeeeeoeeeeee
Answer: [Yes]
Justification:
The research presented in this paper fully conforms to the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. if there is special consid- eration due to laws or regulations 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: We discuss societal impacts of our work in Section 5.
40

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA No, they should explain why their work has societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (c.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. dep.lymell 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 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 being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there 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 system learns from feedback over time, improving the efficiency and accessibility ML).
11. Safeguards
Question:
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Answer: [NA]
Justification: The paper does not release any high risk models or datasets.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: The datasets and models used in this work properly credited and described in Sections 4.2 and 4.3. The code repository referenced in Section A.4.3 includes citations and links to the original sources, along with license information for all third-party assets used.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
41

---

The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if it has changed) should provided.
If this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Detailed documentation for the newly introduced assets, including usage instructions, model architecture, and training details, provided README file code repository linked in Section A.4.3.
Justification:
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/code/mmoel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing research with human subjects.
42

---

Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with
human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent)
may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and 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.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods this research? Note that the LLM used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [Yes]
Justification: LLMs used as integral part of the proposed method, as described in Section 3.3. The detailed experimental setup, including prompting strategy and configuration, is provided in Section 4.3 and A.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non-standard components.
Please refer to our LLM policy (https / /rtee:///nuureee //nurpppp..... (Conferenneee/2222 2L0M/LLMM for what should or should not be described.
43


Reviews:
[{"review_id": "QhdCkdT5BA", "reviewer": "Reviewer_hiE3", "summary": "This paper studies the approximation of score functions using a newly proposed method, named “score-informed neural operator.” This method was applied to causal discovery score-based methods, and applied to an LLM reasoning task.", "questions": "1. The paper starts by introducing causal discovery. However, key related work in ANM is missing. There are ANM methods that start building the causal graph from roots from example. To claim that this method works for causal discovery, it is necessary for the authors to compare against methods beyond score-based approaches. I would actually recommend the authors introducing the method as a generic score approximation method instead,\n2. The causal discovery part and the LLM part feel like two standalone topics, and it is strange to put them into the same paper.\n3. Section 3.1 is disconnected from Section 3.2. The authors need to explicitly spell out 1) the input to the Fourier transform and the output , and 2) how theorem 3.1 is leveraged to improved the previous methods. In addition, it is unclear what ‘k’ means in equation (7). Is the Fourier transform used to ensure higher-order derivatives exist? How is the model operating in the function space? I found the high-level description at the beginning of Section 3.2 unclear. \n4. The baseline architecture in [24] is not described in words in the paper, and the authors did not perform ablation to justify the added modifications. LTE seems like a popular heuristic in modern diffusion model trainings. If separating the Fourier layer into the real and imaginary part is the only novel part, I feel the contribution of the paper is incremental.\n5. The authors use the phrase causal representation to mean causal discovery at the end of page 4. I think the authors should correct that as causal representation is often referred to as a different task.\n6. I found Section 3.3 unclear: what is node set V? For the Confidence interval, the authors needs to justify why the data is I.i.d. In order to apply the central limit theorem. \n7. The setup of the synthetic experiments missing key information in the main body, and the authors should consider adding these key information back. For example, the authors said “we sample node values according to an additive noise model with Gaussian noise.” However, it is important to include the information that the ANM is nonlinear, since linear Gaussian mechanisms are not identifiable.", "limitations": "see above", "rating": 3, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "The paper introduces many topics, and I often find the details missing, hindering the understanding of the paper. Though, the proposed method seems to work better in approximating the score function than prior methods, the framing of the paper focuses on causal discovery, which unnecessarily complicates the problem setting (with key benchmarks in the causal discovery literature missing). The LLM reasoning is disconnected with the rest of the paper, and the setup description lacks clarity.", "was_revised": true, "final_justification": "While the author's rebuttal was useful in clarifying their contributions, they provide a high-level overview of what they are proposing to change in terms of writing rather than concrete changes. I think the amount of rewriting would be significant. However, I also acknowledge that the developed method in this paper improves upon prior work. Thus, I am decreasing my confidence score to 1.", "rebuttal": "We sincerely thank the reviewer for the constructive feedback.\n\n**Q1. Concern over whether the method belongs in causal discovery literature**\n\nAs 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.\n\nTo 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. \n\nWe 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.\n\n**Q2. Connect between causal discovery and LLM control** \n\nAs 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). \n\nRecent 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).\n\nTherefore, 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.\n\n**Q3. Clarification of Sections 3.1–3.2**\n\nWe deeply appreciate the reviewer’s comment, that help us improve the paper’s clarity.\nContrary 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.\n\n1. FFT (inverse FFT) module in each Fourier layer receives (outputs) $\\mathbb{R}^{H}$-valued tensor. \n2. 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).\n3. 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.\n\n**Q4-(a). Missing baseline details**\n\nWe 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.\n\n**Q4-(b). Ablation study: multiplicative vs. additive LTE**\n\nWe 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.\n\n| Dataset\\Method(metric) | SciNO (additive LTE) (OD) | SciNO (multiplicative LTE) (OD) |\n| --- | --- | --- |\n| ER(d10) | 2.8±1.74 | **2.7**±1.19 |\n| ER(d30) | 17.67±5.64 | **16.9**±6.12 |\n| ER(d50) | 38.86±11.42 | **32.8**±6.55 |\n| ER(d100) | 99.2±14.27 | **86.6**±12.99 |\n\n**Q4-(c). Clarification on the novelty of our contribution**\n\nTo 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.\n\nOur 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.\n\nWe validate this framework through a theoretical proof (Section 3.1) and empirical results in our experiments, which we believe constitutes a significant contribution.\n\n**Q5. Misuse of the term 'causal representation'**\n\nWe 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.\n\n**Q6. Clarification of Section 3.3: Node Set V and i.i.d. Assumption**\n\n1. Definition of Node Set \\\\(\\mathcal{V}\\\\)\n    \n    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.\n    \n2. Justification for the i.i.d. Assumption\n    \n    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.\n    \n    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.\n    \n    Therefore, the variance estimates \\\\( (\\sigma_i^{(m)})^2 \\\\)  satisfy the i.i.d. condition.\n\n**Q7. Missing nonlinear ANM details in Synthetic Experiments**\n\nThe 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.'\n\n---\n[Spirtes et al., 2000] — Causation, Prediction, and Search. MIT press, 2000.\n\n[Chickering, 2002] — Optimal Structure Identification with Greedy Search. JMLR, 2002.\n\n[Bühlmann et al., 2014] — CAM: Causal Additive Models, High-Dimensional Order Search and Penalized Regression. The Annals of Statistics, 2014.\n\n[Vashishtha et al., 2025] — Causal Order: The Key to Leveraging Imperfect Experts in Causal Inference. ICLR, 2025.\n\n[Lim et al., 2023] — Score-based Generative Modeling through Stochastic Evolution Equations in Hilbert Spaces. NeurIPS, 2023."}, {"review_id": "96CHgpT9Tz", "reviewer": "Reviewer_Vgze", "summary": "This paper addresses the ordering based causal discovery methods using score matching in the ANM models.\n\nThe causal ordering is identified usually by iteratively identifying the leaf nodes, by estimating Jacobian of the score function, and in particular, its diagonal, denoted as the Hessian diagonal function. However, such estimation is not easy: current methods may suffer from numerical instability, computational complexity, high dimensional issues, etc.\n\nHence, the authors introduce SciNO (Score-informed Neural Operator), a functional diffusion model framework designed to provide stable and scalable approximation of the Hessian diagonal. Score functions are modeled in Sobolev spaces using neural operators.\n\nFurthermore, the authors propose a probabilistic control framework for autoregressive models (e.g., LLMs), integrating SciNO’s score-derived statistics with LLM outputs to guide more accurate causal ordering.", "questions": "/", "limitations": "/", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths:**\n\n1. The motivation is clear and the problem is relevant. The identified challenge of unstable second-order derivative estimation is indeed important in score matching based causal ordering identification.\n\n2. The techniques introduced -- the use of Sobolev space approximations and neural operators -- are new to the causal discovery field (at least to my knowledge). While I haven’t verified the theoretical details in full, the framework appears sound. It would be valuable if the authors could comment on whether these theoretical insights could benefit other problems in causal discovery.\n\n3. The section on Probabilistic Control of Autoregressive Causal Ordering is interesting and seems to bring a new perspective to causal discovery, both for evaluation and for debiasing in inference results.\n\n\n**Weaknesses:**\n\n1. While the proposed LLM control framework is interesting and new, it would help to include a deeper discussion on when and why it works. For example, how much do variable semantics and naming conventions influence the predictions? Are the LLM priors generally informative or misleading in high-dimensional settings? \n\n2. For the techniques, could the authors clarify which parts of the proposed architecture are proposed in this paper and which are adapted from existing methods on other problems? If certain components are adopted elsewhere it would be helpful to discuss what problems they were originally designed to address.\n\n3. Since this is an empirical causal discovery method, the experimental comparisons could be more comprehensive. It would be informative to include baselines beyond the ANM category (e.g., constraint-based or score-based methods) to better contextualize.", "was_revised": true, "final_justification": "Most of my concerns have been addressed. I keep my positive score.", "rebuttal": "We sincerely appreciate the reviewer’s valuable feedback.\n\n**W1. Deeper discussion on when and why LLM control method works**\n\nIn 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.\n\nOur 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.\n\nTo 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.\n\n**W2. Clarification on architectural contributions**\n\nAs 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). \n\nSince 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.\n\nTo 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).\n\n| Dataset\\Method(metric) | DiffAN w/ MLP (LTE) (OD) | DiffAN w/ SciNO (LTE) (OD) |\n| --- | --- | --- |\n| Physics(d7) | 5.0±1.1 | **2.0**±0.0 |\n| Sachs(d8) | 7.2±2.71 | **5.8**±2.93 |\n| MAGIC-NIAB(d44) | 10.6±12.72 | **3.8**±0.75 |\n| ECOLI70(d46) | 20.4±1.85 | **14.4**±3.61 |\n| MAGIC-IRRI(d64) | 25.0±23.06 | **10.0**±0.89 |\n| ARTH150(d107) | 82.0±23.28 | **20.8**±2.32 |\n\nTo 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).\n\n**W3. Comparison with non-ANM baselines**\n\nWe 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.\n\nIn 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. \n\n| Dataset\\Method(metric) | PC(SHD-C) | GES(SHD-C) | DiffAN w/ SciNO(SHD-C) |\n| --- | --- | --- | --- |\n| ER(d2) | 0.1±0.3 | 0.2±0.4 | **0.0**±0.0 |\n| ER(d3) | 0.6±0.92 | 1.3±1.01 | **0.0**±0.0 |\n| ER(d5) | 9.03±0.64 | 8.4±1.43 | **1.3**±2.41 |\n| ER(d10) | 34.6±2.93 | 35.2±2.96 | **23.8**±5.49 |\n| ER(d30) | 108.5±15.69 | 105.7±15.3 | **88.3**±18.11 |\n| ER(d50) | 195.9±12.08 | 180.2±12.45 | **174.3**±14.48 |\n\nAs 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.\n\n---\n[Lim et al., 2023] — Score-based Generative Modeling through Stochastic Evolution Equations in Hilbert Spaces. NeurIPS, 2023.\n\n[Li et al., 2021] — Learnable Fourier Features for Multi-Dimensional Spatial Positional Encoding. NeurIPS, 2021.\n\n[Spirtes et al., 2000] — Causation, Prediction, and Search. MIT press, 2000.\n\n[Chickering, 2002] — Optimal Structure Identification with Greedy Search. JMLR, 2002.\n\n[Bühlmann et al., 2014] — CAM: Causal Additive Models, High-Dimensional Order Search and Penalized Regression. The Annals of Statistics, 2014."}, {"review_id": "VTdCgjhc2u", "reviewer": "Reviewer_2fQy", "summary": "Recent methods for causal discovery rely on iteratively finding leaf nodes using the second-order derivative of the log-densities. However, current methods show a lack of stability and accuracy in the estimation of such a derivative. This method leverages and improves upon recent advances in diffusion models for more accurate estimation of derivative. These scores can be used to iteratively find leaf nodes, or to improve LLM-driven causal discovery. Experiments show the improved performance induced by the new estimation method.", "questions": "Questions are in the Weaknesses above; I am ready to significantly raise my score if they are mostly addressed.", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths: the method is original and is a welcome technical contribution, the paper is generally clear.\n\nWeaknesses: the Theorem is a bit narrow and claims in the experiments are dubious.\n\nTheory (Theorem 3.1) :\n\n1a) The assumptions on the Sobolev space are strong; this means k>27 for D=50 and k>52 for D=100!\n\n1b) Further, it is not applied to a score with a time variable, at least when the time is close to 0.\n\nMain concerns on experiments:\n\n2a) Only 5 independent runs are run in Sections 4.1-4.2, and even more concerning, I understand that the best of 3 runs is Section 4.3.\n\n2b) l.230: \"Table 2 shows that SciNO consistently outperforms the original DiffAN across all metrics\" : I see ties in a majority of comparisons in Table 2, both for ER and \"real-world\" (see 3b) graphs.\n\n2c) \"SciNO achieves comparable performances across all datasets, with improved results observed on higher-dimensional graphs such as MAGIC-IRRI and ARTH150.\" The baseline is better with statistical signifiance in a majority of comparisons in Table 3, even most of MAGIC-IRRI.\n\n3) More minor points in experiments\n\n3a) There is no comparison with the original SCORE method while it is extensively discussed before.\n\n3b) BNLearn and Physics are not \"real-world\" data as I understand it, I would downgrade the \"real data\" to \"semi-synthetic data\"\n\n3c) A time complexity comparison between CaPS and the method is done in l.173-174, but Figure 4 does not show runtime.\n\n3d) The analysis of SciNO PE vs LTE is not done on other datasets than ER.", "was_revised": true, "final_justification": "On my side, all issues were resolved. We'll debate the paper's acceptance or not with other reviewers", "rebuttal": "We would like to thank the reviewer for very helpful suggestions.\n\n**W1-(a). Strong Sobolev space assumption**\n\nWe appreciate detailed questions about our theoretical results.\n\nThe 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.\n\n**W1-(b). Inapplicability near t=0 for time-conditioned scores**\n\nThanks 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.\n\n**W2-(a). Limited number of runs in Section 4.1-4.3**\n\nTo address these concerns about statistical robustness, we performed additional experiments, and the new results will be included in our revision.\n\n1. Experiments in Section 4.1-4.2\n    \n    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.\n    \n\n    - **Table 2**\n  \n        | Dataset\\Method(metric) | DiffAN(OD) | DiffAN(SHD) | DiffAN(SID) | SciNO(OD) | SciNO(SHD) | SciNO(SID) |\n        | --- | --- | --- | --- | --- | --- | --- |\n        | 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 |\n        | 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 |\n        | 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 |\n        | 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 |\n        | 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 |\n        | 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 |\n    \n    - **Table 3**\n    \n        | Dataset\\Method(metric) | CaPS w/ SciNO(OD) | CaPS w/ SciNO(SHD) | CaPS w/ SciNO(SID) |\n        | --- | --- | --- | --- |\n        | MAGIC-NIAB(d44) | 37.5±0.50 | 165.4±0.5 | 1185.9±5.61 |\n        | ECOLI70(d46) | 25.4±0.80 | 141.3±4.78 | 882.4±13.02 |\n        | MAGIC-IRRI(d64) | **41.6**±0.80 | 183.5±2.97 | 1280.7±36.83 |\n        | ARTH150(d107) | **48.2**±0.87 | **380.0**±4.63 | **2799.1**±73.26 |\n\n2. Experiments in Section 4.3\n    \n    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.\n    \n    | Method\\Dataset | Physics(d7) | Sachs(d8) | MAGIC-NIAB(d44) | ECOLI70(d46) | MAGIC-IRRI(d64) | ARTH150(d107) |\n    | --- | --- | --- | --- | --- | --- | --- |\n    | 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 |\n    \n**W2-(b). Overgeneralization of performance claims**\n\nTo 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.\"\n\n**W2-(c). Overstatement of SciNO’s advantages in MAGIC-IRRI and ARTH150**\n\nWe 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.”\n\n**W3-(a). Missing empirical comparison with SCORE**\n\nWe 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.\n\n| Dataset\\Method(metric) | SCORE(OD) | SCORE(SHD) | SCORE(SID) | DiffAN w/ SciNO (OD) | DiffAN w/ SciNO (SHD) | DiffAN w/ SciNO (SID) |\n| --- | --- | --- | --- | --- | --- | --- |\n| ER(d2) | 0.1±0.3 | 0.2±0.60 | 0.2±0.50 | **0.0**±0.0 | **0.0**±0.0 | **0.0**±0.0 |\n| ER(d3) | **0.0**±0.0 | **0.0**±0.00 | **0.0**±0.00 | 0.1±0.3 | 0.2±0.6 | 0.4±1.2 |\n| ER(d5) | **0.0**±0.0 | **0.1**±0.3 | **0.1**±0.3 | 0.9±1.14 | 1.7±1.95 | 4.8±4.94 |\n| ER(d10) | 4.2±2.60 | 21.5±5.24 | **41.2**±13.24 | **2.7**±1.19 | **20.5**±3.5 | 41.6±9.31 |\n| ER(d30) | 24.8±11.97 | 92.8±16.27 | 510.8±90.77 | **16.9**±6.12 | **88.8**±16.5 | **492.0**±91.64 |\n| ER(d50) | 49.3±14.37 | 180.8±15.99 | **1527.1**±179.61 | **32.8**±6.55 | **180.0**±15.45 | 1622.0±110.92 |\n| ER(d100) | 120.0±18.49 | 455.3±26.34 | **7175.0**±435.00 | **86.6**±12.99 | **445.9**±32.84 | 7259.2±653.27 |\n| Physics(d7) | **1** | **2** | **12** | 2.0±0.0 | 5.8±0.40 | **7.8**±1.94 |\n| Sachs(d8) | 6 | 16 | 30 | **5.8**±2.93 | **13.8**±1.72 | **19.4**±5.89 |\n| MAGIC-NIAB(d44) | 8 | **61** | 189 | **3.8**±0.75 | 67.4±2.87 | **156.6**±18.04 |\n| ECOLI70(d46) | **7** | **80** | 771 | 14.4±3.61 | 91.0±7.16 | **734.0**±50.73 |\n| MAGIC-IRRI(d64) | 17 | 175 | 529 | **10.0**±0.89 | **135.4**±5.12 | **333.2**±17.34 |\n| ARTH150(d107) | 38 | 460 | 2709 | **20.8**±2.32 | **454.6**±7.66 | **1740.4**±108.63 |\n\n**W3-(b). Terminology for 'Real-World' vs. 'Semi-Synthetic' Data**\n\nIn 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'.\n\n**W3-(c). Missing empirical runtime comparison with CaPS**\n\nWhile 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. \n\n| Sample Size\\Method | CaPS(sec) | CaPS w/ SciNO(sec/epoch) |\n| --- | --- | --- |\n| 100 | 2.04 | **0.07** |\n| 1000 | 12.08 | **6.94** |\n| 10000 | 77.07 | **62.96** |\n| 100000 | OOM | **156.71** |\n\nWe will add a plot of these results to the experimental section of our revised paper. We thank the reviewer again for this valuable suggestion.\n\n**W3-(d). PE vs. LTE ablation on semi-synthetic datasets**\n\nTo 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.\n\n| Method(metric)\\Dataset | Physics(d7) | Sachs(d8) | MAGIC-NIAB(d44) | ECOLI70(d46) | MAGIC-IRRI(d64) | ARTH150(d107) |\n| --- | --- | --- | --- | --- | --- | --- |\n| DiffAN w/ SciNO (PE) (OD) | 2.7±0.46 | **4.4**±1.5 | **4.0**±1.55 | 30.6±2.42 | 16.6±1.12 | 33.9±3.94 |\n| DiffAN w/ SciNO (LTE) (OD) | **1.9**±0.54 | 5.6±2.65 | 4.1±0.70 | **12.8**±3.19 | **10.6**±1.69 | **21.4**±2.76 |\n\n---\n[Lakshminarayanan et al., 2017]  — Simple and scalable predictive uncertainty estimation using deep ensembles. NeurIPS, 2017.  \n\n[Lim et al., 2023] — Score-based Generative Modeling through Stochastic Evolution Equations in Hilbert Spaces. NeurIPS, 2023."}, {"review_id": "yDnkhhI5Pu", "reviewer": "Reviewer_2jQc", "summary": "This paper studies the problem of ordering based causal discovery. Identifying that an existing method DiffAN suffers from unstable gradient estimation,  authors proposed Score-inforrmed Neural Operator (SciNO) to stably approximate the Hessian diagonal and to preserve structural information during the score modeling.  This is achieved by introducing two key modifications to the time-conditional Fourier neural operator architecture: a learnable time encoding module and that decomposition of signals in Fourier layers into real and imaginary parts. Further, the paper proposes a probabilistic control method that leverages score-informed statistics inferred by SciNO, guiding autoregressive models to generate more reliable causal reasoning. Empirical results validate that the improved stability of gradient estimation leads to better causal discovery performance compared with DiffAN.", "questions": "1. can you give more details on : \"injecting noise into the data can destroy the functional information inherent in the original score function\" (line 50-51) ？\n\n2. line 150 :\".. allows SciNO to achieve robust performance in high-dimensional causal graphs\": 'robust performance' is vague here. Do you mean the trend of SHDs along number of nodes or the standard deviation? Please make it accurate. \n\n3. regarding Section 3.3: it seems that the probabilistic control can also be applied to the original DiffAN? if so, what would be the difference?\n\n4. **Major Concern**: The major weakness is the limited empirical comparison, lacking classic/existing benchmark methods. In particular, the paper only compares the proposed method with its base method DiffAN. It is not clear how the algorithm performs within the whole causal discovery literature. For example, the CAM method shall  perform well on the synthetic data setup of this paper. See, e.g., [1, 2] for references, where their setup is the identical or at least very similar to the one in this paper. In particular, for ER4d (d=100) case, CAM achieves  98.8$\\pm$ 20.7 SHD, while the proposed method achieves 180.0$\\pm$15.45. Of course, a new algorithm does not need to beat all existing methods, but this large difference needs further clarifying. \n\n    [1] Gradient-Based Neural DAG Learning, ICLR, 2020.\n\n   [2] Ordering based causal discovery with reinforcement learning, IJCAI, 2021.\n\n5. minor on paper writing/organization: \n\n  - line 150: almost half of section 3.2 is about ablation study or case study, which distracts the reading flow. Suggest to move it to experiment section. \n  - tables 1-4: table caption is generally placed on top.\n  - line 126: \"HDM implements time-conditioned Fourier Neural Operator (FNO, [22]) by adding projected tensors from PE (Positional Embedding, [40]) \": be consistent about the format of abbreviation.", "limitations": "Yes. \n\nAnd see above the limitation on lacking important benchmark methods. At least further clarification is needed.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n- The paper identifies the key limitation of an existing algorithm, and the modification is sound. \n- The proposed method is efficient, compared with  DiffAN.\n- Writing is generally good (though can be further improved; see below.)\n\nWeakness:\n\n- Lacking comparison with classic algorithms.\n- Novelty is a bit limited.", "was_revised": true, "final_justification": "Authors address all my concerns and validate the empirical performance by new experimental results with CAM, which is beyond my expectation. I suggest authors add more benchmark methods, e.g., those gradient-based algorithms in the same table in future versions.", "rebuttal": "We are grateful for your detailed and helpful comments.\n\n**Q1. Details on how noise injection destroys functional information**\n\nThe 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.\n\n**Q2. Clarification of 'robust performance' (L150)**\n\nWe 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'.\n\n**Q3. Applicability of the control method to DiffAN**\n\nThe 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.\n\nThe 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.\n- DiffAN-based control (Comparison of order divergence between uncontrolled Llama and controlled Llama)\n    | Method\\Dataset | Physics(d7) | Sachs(d8) | MAGIC-NIAB(d44) | ECOLI70(d46) | MAGIC-IRRI(d64) | ARTH150(d107) |\n    | --- | --- | --- | --- | --- | --- | --- |\n    | Llama-3.1-8b(1-norm) | 3 | 3 | 9 | 32 | 18 | 80 |\n    | w/ control(rank) | 4 (↑33.3%) | 3 | 2 (↓77.8%) | 22 (↓31.3%) | **7** (↓61.1%) | **25** (↓68.8%) |\n    | w/ control(CI) | 3 | 4 (↑33.3%) | 6 (↓33.3%) | **16** (↓50.0%) | 13 (↓27.8%) | 52 (↓35.0%) |\n    | Llama-3.1-8b(0.5-norm) | 1  | 3 | 14 | 32 | 17 | 80 |\n    | w/ control(rank) | 4 (↑300%) | 3 | **1** (↓92.9%) | 20 (↓37.5%) | 8 (↓52.9%) | **25** (↓68.8%) |\n    | w/ control(CI) | 4 (↑300%) | 5 (↑66.7%) | 8 (↓42.9%) | 18 (↓43.8%) | 14 (↓17.6%) | 52 (↓35.0%) |\n\n- SciNO-based control (Comparison of order divergence between uncontrolled Llama and controlled Llama, Table 4)\n    | Method\\Dataset | Physics(d7) | Sachs(d8) | MAGIC-NIAB(d44) | ECOLI70(d46) | MAGIC-IRRI(d64) | ARTH150(d107) |\n    | --- | --- | --- | --- | --- | --- | --- |\n    | Llama-3.1-8b(1-norm) | 3 | 3 | 9 | 32 | 18 | 80 |\n    | w/ control(rank) | 1 (↓66.7%) | 5 (↑66.7%) | 4 (↓55.6%) | **11** (↓65.6%) | **9** (↓50.0%) | 21 (↓73.8%) |\n    | w/ control(CI) | 1 (↓66.7%) | 3 | **2** (↓77.8%) | 15 (↓53.1%) | 10 (↓44.4%) | **18** (↓77.5%) |\n    | Llama-3.1-8b(0.5-norm) | 1  | 3 | 14 | 32 | 17 | 80 |\n    | w/ control(rank) | 1 | 4 (↑33.3%) | **2** (↓85.7%) | **10** (↓68.8%) | **9** (↓47.1%) | 21 (↓73.8%) |\n    | w/ control(CI) | 1 | 4 (↑33.3%) | **2** (↓85.7%) | 14 (↓56.3%) | 10 (↓41.2%) | **18** (↓77.5%) |\n\n**Q4. Performance gap in SHD between SciNO and CAM**\n\nWe 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.\n\nTo 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. \n\n| Dataset\\Method(metric) | CAM(SHD) | CAM(SID) | DiffAN w/ SciNO(SHD) | DiffAN w/ SciNO(SID) |\n| --- | --- | --- | --- | --- |\n| ER(d2) | **0.0**±0.0 | **0.0**±0.0 | **0.0**±0.0 | **0.0**±0.0 |\n| ER(d3) | **0.0**±0.0 | **0.0**±0.0 | 0.2±0.6 | 0.4±1.2 |\n| ER(d5) | **0.0**±0.0 | **0.0**±0.0 | 1.7±1.95 | 4.8±4.94 |\n| ER(d10) | **16.6**±4.45 | **30.9**±11.96 | 20.5±3.5 | 41.6±9.31 |\n| ER(d30) | 129±14.99 | **450.2**±99.2 | **88.8**±16.5 | 492.0±91.64 |\n| ER(d50) | 238.5±17.29 | **1503.6**±201.35 | **180.0**±15.45 | 1622.0±110.92 |\n| ER(d100) | 537.2±45.37 | **6289.9**±533.87 | **445.9**±32.84 | 7259.2±653.27 |\n| Physics(d7) | 12 | 11 | **5.8**±0.40 | **7.8**±1.94 |\n| Sachs(d8) | 37 | 41 | **13.8**±1.72 | **19.4**±5.89 |\n| MAGIC-NIAB(d44) | 167 | 384 | **67.4**±2.87 | **156.6**±18.04 |\n| ECOLI70(d46) | 184 | 828 | **91.0**±7.16 | **734.0**±50.73 |\n| MAGIC-IRRI(d64) | 261 | 791 | **135.4**±5.12 | **333.2**±17.34 |\n| ARTH150(d107) | **437** | 2133 | 454.6±7.66 | **1740.4**±108.63 |\n\nTherefore, 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.\n\n**Q5-(a). Placement of case study in method section**\n\nWe 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.\n\n**Q5-(b). Table caption formatting** \n\nThank you for pointing this out. We will move all table captions to the top to adhere to the standard format.\n\n**Q5-(c). Abbreviation formatting consistency** \n\nWe 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].\n\n---\n[Bühlmann et al., 2014] — CAM: Causal Additive Models, High-Dimensional Order Search and Penalized Regression. The Annals of Statistics, 2014.\n\n[Lachapelle et al., 2020] — Gradient-Based Neural DAG Learning. ICLR, 2020.\n\n[Wang et al., 2021] — Ordering based causal discovery with reinforcement learning. IJCAI, 2021.\n\n[Spirtes et al., 2000] — Causation, Prediction, and Search. MIT press, 2000.\n\n[Chickering, 2002] — Optimal Structure Identification with Greedy Search. JMLR, 2002."}, {"review_id": "99z5IB0zTk", "reviewer": "Reviewer_uvzA", "summary": "This paper focuses on a key challenge of ordering-based causal discovery: the computationally expensive and memory-intensive issues of score's Jacobian (Hessian) diagonal estimation. Motivated by FNO and DiffAN, this paper proposes a new time-conditioning score-based generative model, Score-informed Neural Operator (SciNO), to stably approximate the Hessian diagonal by modeling the score function with neural operators. This architecture can improve the scalability of the existing algorithms, e.g., SCORE and CaPS, by reducing their cubic computational complexity. Additionally, combining LLM with the leaf node discrimination criteria in previous works, this paper proposes a probabilistic control of autoregressive causal ordering algorithm. With the approximation power of neural operators in Theorem 3.1, the idea is theoretically sounds. Some synthetic and real-world experiments are presented.", "questions": "1. Could you further explain what the subgraph below in Figure 5 intends to express? What specific distributions do the prior and posterior here refer to? What does this change indicate?\n\n2. See strengths and weaknesses. Most of my concerns are about the experimental part. I will adjust the score if the above experiments can reach a consistent conclusion.", "limitations": "yes", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 2, "strengths_and_weaknesses": "Pros:\n\n1. This paper is well-motivated. Improving the estimation efficiency of Hessian diagonal can improve the scalability of many existing ordering-based causal discovery methods.\n\n2. Theoretical motivations are clearly explained and the proof are adequately provided. Theorem 3.1 shows the approximation power of neural operators. The theorems and mathematical details appear to be correct, though I haven't checked all the details.\n\n3. Despite some formatting issues, this paper is well-written and easy to understand.\n\nCons:\n\n1. Lack of some important experiment. Although the main purpose of this paper is to reduce the computational complexity of Hessian diagonal estimation, the accuracy of the approximation is also important. Here I have two main concerns.\n\n* About the accuracy of Hessian diagonal estimation, SciNO can have a better or comparable performance compared to the second-order Stein gradient estimator? The author can illustrate this point by adding the second-order Stein gradient estimator in figure 2 or by presenting a comparison of order divergence.\n\n* Although SciNO can enhance the scalability of CaPS, it seems that it cannot stably improve the performance of CaPS. Does SciNO perform poorly on low-dimensional data? Providing more datasets (like Table 2) can be helpful for this analysis.\n\n2. This paper seems to be concatenated together by two independent parts, one is SciNO and the other is the autoregressive LLM control algorithm. I suggest that at least a comparison between the LLM-based method and the previous methods, such as DiffAN and CaPS, should be provided.", "was_revised": true, "final_justification": "See the feedback in my official comment.", "rebuttal": "We appreciate your valuable comments and suggestions.\n\n**W1-(a). Accuracy of Hessian diagonal approximation** \n\nTo 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.\n\n| Graph\\Method | Stein gradient estimator | SciNO  |\n| --- | --- | --- |\n| 1 | -**0.03**±1.15 | -0.05±0.72 |\n| 2 | -0.1±1.21 | -**0.04**±0.5 |\n| 3 | -0.22±1.63 | -**0.03**±0.21 |\n| 4 | -0.07±1.13 | -**0.01**±0.07 |\n| 5 | -**0.01**±1.01 | -0.02±0.08 |\n| 6 | -0.22±1.56 | -**0.04**±0.06 |\n| 7 | -0.1±0.92 | -**0.05**±1.28 |\n| 8 | -**0.03**±1.01 | -0.05±0.22 |\n| 9 | -0.18±1.56 | -**0.06**±0.11 |\n\nFurthermore, 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.\n\n**W1-(b). Stability of CaPS w/ SciNO on low-dimensional datasets**\n\nTo assess potential performance degradation, we compared CaPS and SciNO on five low-dimensional datasets. \n\n| Dataset\\Method(metric) | CaPS(OD) | CaPS w/ SciNO(OD) |\n| --- | --- | --- |\n| ER(d2) | 0.3±0.46 | **0.0**±0.0 |\n| ER(d3) | 0.4±0.66 | **0.2**±0.4 |\n| ER(d5) | **0.6**±0.92 | **0.6**±0.92 |\n| Physics(d7) | **0** | **0.0**±0.0 |\n| Sachs(d8) | 5 | **4.4**±2.29 |\n\nThese 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.\n\n**W2. Connection between causal discovery and LLM control, with comparison to previous methods**\n\nOur 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. \n\n1. **Comparison with DiffAN-based control**\n    \n    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.\n    \n2. **Comparison with SCORE-based control**\n    \n    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.\n\n- SCORE-based control (Comparison of order divergence between uncontrolled Llama and controlled Llama)\n\n    | Method\\Dataset | MAGIC-NIAB(d44) | ECOLI70(d46) | MAGIC-IRRI(d64) | ARTH150(d107) |\n    | --- | --- | --- | --- | --- |\n    | Llama-3.1-8b(1-norm) | 9 | 32 | 18 | 80 |\n    | w/ control(rank) | **1** (↓88.9%) | **14** (↓56.3%) | 11 (↓38.9%) | 45 (↓43.8%) |\n    | w/ control(CI) | **1** (↓88.9%) | **14** (↓56.3%) | 11 (↓38.9%) | **43** (↓46.3%) |\n    | Llama-3.1-8b(0.5-norm) | 14 | 32 | 17 | 80 |\n    | w/ control(rank) | **1** (↓92.9%) | 16 (↓50.0%) | 10 (↓41.2%) | 45 (↓43.8%) |\n    | w/ control(CI) | 2 (↓85.7%) | 16 (↓50.0%) | **7** (↓58.8%) | **43** (↓46.3%) |\n\n- SCORE\n\n    | Method\\Dataset | MAGIC-NIAB(d44) | ECOLI70(d46) | MAGIC-IRRI(d64) | ARTH150(d107) |\n    | --- | --- | --- | --- | --- |\n    | SCORE  | 8 | 7 | 17 | 38 |\n    \n\n**Q1. Detailed explanation of what Figure 5's subgraph shows**\n\nAs 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.\n\nThe 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).\n\nAs 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).\n\n--- \n[Vashishtha et al., 2025] — Causal Order: The Key to Leveraging Imperfect Experts in Causal Inference. ICLR, 2025.\n\n[Lakshminarayanan et al., 2017]  — Simple and scalable predictive uncertainty estimation using deep ensembles. NeurIPS, 2017."}]

Author rebuttal:

Meta-review: The paper contribution is SciNO, a functional diffusion model framework that enhances ordering-based causal discovery algorithms, 
building upon [Vashishtha et al., 2025], who show how LLMs can contribute to better ranking performances. 

The authors aim to address the limitations of Ordering-based Causal Discovery methods that identify leaf nodes based on the Hessian diagonal (empirically and convincingly showing that Stein gradient-based approaches suffer from similar or worse limitations as Hessian based approaches). 

Most concerns of the reviewers were related with insufficient empirical studies. These have been extensively addressed by the authors, considering more baselines, more datasets and more runs, and presenting lesion studies.

Other concerns were voiced, regarding the organization of the paper, as the parts related to SciNO and LLM seemed disconnected at first, but the authors did an excellent job in their rebuttal to explain this organization, and slightly changing the focus/presentation of the paper. 

The proposed revision plan comprehensively addresses the concerns.
Decision: Accept (poster)