Learning Dynamic Causal Graphs Under Parametric Uncertainty via Polynomial Chaos Expansions
Authors: Liang Cao
Venue: ICLR 2026 Poster (2026)
Keywords: Causal Discovery
Polynomial Chaos Expansion
Parametric Uncertainty
Functional Causal Models
Uncertainty Quantification
Primary area: causal reasoning

Abstract: Existing causal discovery methods are fundamentally limited by the assumption of a static causal graph, a constraint that fails in real-world systems where causal relationships dynamically vary with underlying system parameters. This discrepancy prevents the application of causal discovery in critical domains such as industrial process control, where understanding how causal effects change is essential. We address this gap by proposing a new paradigm that moves beyond static graphs to learn functional causal representations. We introduce a framework that models each causal link not as a static weight but as a function of measurable system parameters. By representing these functions using Polynomial Chaos Expansions (PCE), we develop a tractable method to learn the complete parametric causal structure from observational data. We provide theoretical proofs for the identifiability of these functional models and introduce a novel, provably convergent learning algorithm. On a large-scale chemical reactor dataset, our method learns the dynamic causal structure with a 90.9% F1-score, nearly doubling the performance of state-of-the-art baselines and providing an interpretable model of how causal mechanisms evolve.

Published as a conference paper at ICLR 2026
LEARNING DYNAMIC CAUSAL GRAPHS UNDER PARA- METRIC UNCERTAINTY VIA POLYNOMIAL CHAOS Ex- PANSIONS
Liang Cao
Department of Chemical and Biological Engincering The University of British Columbia Vancouver, BC V6T 173, Canada
clube1 9 @ma 1. .bc.ca
ABSTRACT
Existing causal discovery methods are fundamentally limited by the assumption of a static causal graph, a constraint that fails in real-world systems where causal rela- tionships dynamically vary with underlying system parameters. This discrepancy prevents the application of causal discovery critical domains such as industrial process control, where understanding how causal effects change is essential. We address this gap by proposing a new paradigm that moves beyond static graphs learn functional causal representations. We introduce framework that mod- els each causal link not as static weight but function of measurable system parameters. By representing these functions using Polynomial Chaos Expansions (PCE), we develop a tractable method to learn the complete parametric causal structure from observational data. We provide theoretical proofs for the identi- fiability of these functional models and introduce novel, provably convergent learning algorithm. On a large-scale chemical reactor dataset, our method learns the dynamic causal structure with 90.9% Fi-score, nearly doubling the perfor- mance of state-of-the tht-hhe-rtt art baselines and providing an interpretable model of how causal mechanisms evolve.
INTRODUCTION
Industrial process control systems generate massive volumes sensor data requiring automated analysis for optimization and predictive maintenance (Fang et al., 2024; Zhou et al., 2016). Under- standing causal relationships between process variables is essential for root cause analysis, anomaly detection, and adaptive control strategies (Zhang al., 2016). However, industrial processes hibit unique challenges that violate assumptions of existing causal discovery methods: causal rela- tionships may vary systematically with operating conditions, sensors exhibit complex multi-modal and heavy-tailed noise distributions, and safety-critical applications demand rigorous uncertainty quantification (Cao et al., 2025; Wang et al., 2025).
Many widely used causal discovery methods for observational data are formulated in terms of a single, static causal graph whose edge strengths do depend on observed context or operating pa- rameters, even though there is a growing body of work on time-varying and context-specific causal structures (Song et al.. 2009; Huang et al., 2019). In reality, industrial causal effects are functions of measurable parameters. For instance, in chemical reactors, the influence of feed temperature on product quality depends strongly on catalyst activity, which degrades over time. Heat exchanger effectiveness varies with fouling levels, fundamentally altering thermal control loops. These para- metric dependencies are not mere nuisances but contain critical information for process optimization and predictive maintenance. Our goal in this paper is therefore not replace existing approaches to epistemic or aleatoric uncertainty, but to complement them with a representation which each causal edge is explicit function of a -iw-dinnnsional vector of operating parameters.
The field of causal discovery has evolved through three major directions, each with distinct limita- tions for industrial applications. Constraint-basse methods such as the Peter-Clark (PC) algorithm

---

Published as a conference paper at ICLR 2026
Discovery
Uncertainty Quantification
Industrial Data
Discovered Dynamic Causal Graph
Imail
MMluy York
PABS
Edge Rafinament
Figure 1: Overview of the polynomial chaos theory for causal discovery framework in dynamic uncertainty systems.
(Spirtes et al., 2000), Fast Causal Inference (FCI) (Spirtes et al., 1999), and Really Fast Causal Inference (RFCI) (Colombo et al., 2012) rely on conditional independence tests infer causal structures. While theoretically sound, these methods struggle with finite sample sizes and become computationally intractable for high-dimensional industrial data with hundreds sensors. Recent advances like PC-stable (Colombo & Maathuis, 2014) improve stability but fail when faced with complex noise distributions common in industrial sensors.
Score-based methods including Greedy Equivalence Search (GES) (Chickering. 2002) and Fast Greedy Equivalence Search (FGES) (Ramsey et al., 2017) optimize scoring functions over possible Directed Acyelic Graph (DAG) structures. The breakthrough NOTEARS algorithm (Zheng et al., 2018) reformulated structure learning as continuous optimization with differentiable acyclicity con- straints, spawning variants like Directed Acyelic Graph Graph Neural Network (DAG-GNN) (Yu et al., 2019) for nonlinear relationships, Reinforcement Learning- Bayesian Information Criterion (RL-BIC) (Zhu et al., 2020) using reinforcement learning. In their most common instantiations, these approaches return a single best-scoring DAG and point estimates of edge weights; uncertainty over graphs and parameters typically handled by separate Bayesian or bootstrap procedures (e.g., Lorch et al., 2021; Cundy et al., 2021) rather than being integrated with an explicit model of how edge strengths vary with operating parameters. This limitation is particularly critical for safety- critical industrial applications where confidence in causal recommendations directly impacts opera- tional decisions and safety outcomes.
Functional causal models exploit asymmetries in data distributions for identifiability. The Linear Non-Gaussian Acyclic Model (LiNGAM) (Shimizu et al., 2006) proved that linear models with non-Gaussian noise yield unique causal structures, later extended to DirectLiNGAM (Shimizu et al., 2011) and Vector Autoregressive LiNGAM (VAR-LiNGAM) (Hyvärinen et al., 2010) for time series. Nonlinear extensions include additive noise models (ANM) (Hoyer et al., 2008), post- nonlinear models (PNL) (Zhang& & Hyvärinen, 2009), and the general identifiable functional causal model framework (Peters et al., 2014). However, in their standard form these models do not repre- sent how causal effects change as an explicit function of observed operating parameters.
Recent industrial applications (Sui al., 2025; Runge et al., 2019; Yang et al., 2025) have high- lighted these limitations, often requiring extensive preprocessing or domain-specific modifications lacking theoretical justification. Bayesian approaches such as Differentiable Bayes for Structure Learning (DiBS) (Lorch et al., 2021) and Bayesian Causal Discovery with Neural Networks (BCD Nets) (Cundy et al., 2021) explicitly quantify posterior uncertainty over graphs and parameters, ad- dressing epistemic uncertainty due to finite data, but they still treat each edge as static and do not model how its strength varies with operating parameters.
Polynomial Chaos Expansion (PCE), introduced by Wiener (Wiener, 1938) and generalized by Xiu (Xiu & Karniadakis, 2002), provides mathematically rigorous framework for representing and propagating uncertainty through complex systems. PCE has been successfully applied in uncertainty quantification for engineering systems (Sudret, 2008), sensitivity analysis (Crestaux et al., 2009), and process monitoring (Cao et al., 2026). Recent algorithmic advances including sparse PCE (Jiang et al., 2025), adaptive basis selection (Dai et al., 2025), and multi-fidelity approaches (Liu al., 2021) have made PCE computationally tractable for high-dimensional problems. Despite this
2

---

Published as a conference paper at ICLR 2026
success in forward uncertainty propagation, PCE has, to the best of our knowledge, not yet been systematically exploited for causal discovery.
This paper introduces PCT-CD, bridging uncertainty quantification and causal discovery. Figure provides an overview of the proposed framework. Our contributions are summarized as follows:
1. From Static Graphs Dynamic Functions:
We formalize an industrial structural causal model in which a single underlying DAG is equipped with edge weights that are explicit functions of operating conditions, and we prove identifiability these parametric mecha- nisms under mild assumptions.
2. An End-to-End Solution for Dynamic Systems: We propose PCT-CD, an integrated algo- rithm specifically designed for parametric causal discovery. It translates complex process data into an interpretable model of how causal links evolve, providing actionable intelli- gence beyond simple correlation.
3. Demonstrated Real-World Viability: We empirically validate PCT-CD on controlled benchmark with parameter-varye varriig mechanisms. PCT-CD achieves high F1-scores and equips engineers with uncertainty-awee tools that are essential for making robust deci- sions in high-stakes environments.
The remainder of this paper is organized as follows. Section 2 details the PCT-CD methodology, in- cluding parametric SEM formulation, PCE representation, and our novel conditional independence test. Section 3 establishes theoretical foundations with identifiability and convergence proofs. Sec- tion 4 validates our approach on synthetic benchmark with controlled parametric variation and an industrial process dataset, comparing against a broad set of baseline methods. Section 5 concludes with implications for industrial deployment and future research directions.
2
METHODOLOGY
Our proposed framework, PCT-CD, introduces a novel paradigm for causal discovery by explic- itly modeling how causal relationships vary functions of measurable system parameters. This achieved by integrating the theory of PCE into a hybrid structure learning algorithm. While Bayesian and bootstrap-basse methods typically quantify epistemic uncertainty arising from finite data (for example via posterior distributions over static graphs and parameters), PCT-CD is designed address parametric uncertainty by representing causal edges explicit functions of measurable system parameters, and is conceptually complementary to these existing approaches.
The methodology unfolds in four stages: first, we formulate a Structural Equation Model (SEM) where causal coefficients are functions of a parameter vector Ş. Second, we represent these functions using PCE, transforming the non-parametric problem into a tractable parametric one. Third, we develop novel conditional independence test tailored to this representation to discover an initial causal skeleton. Finally, we refine this structure and quantify edge strengths using score-based optimization with natural gradient approach, ensuring both accuracy and computational efficiency.
2.1 PROBLEM FORMULATION
We consider a complete probability space (E,F,P) where all random quantities are defined. The core innovation of our framework is the explicit modeling of parametric uncertainty through ran- dom vector & E = C Rd This vector represents known, measurable operating conditions (e.g., ambient temperature, catalyst age, feedstock quality) with a joint probability distribution jęg that has finite moments of all orders. In our theoretical analysis we assume that læ is known so that stan- dard PCE basis adapted to μlg can be chosen; in practice, when only samples of g are available, an empirical orthogonal basis can be constructed from the observed parameter values (see Appendix for details). This formulation emphasizes a dimension that is often implicit in causal discovery, where operating conditions typically treated as fixed and uncertainty is mainly modeled arising from finite data and stochastic noise.
We observe process variables, collected in vector X these variables are generated by linear SEM (equivalently. a linear structural causal model in the sense of structural causal inference) where the causal relationships are functions of the parameter
XT,)T E R We assume that
3

---

Published as a conference paper at ICLR 2026
vector
x, Σ (i)(E) 4 1=1,...,, JEPA;
where PA, C (I....,n) (i) denotes the set of causal parents of variable X,, the functions 6 are unknown, -nnereeittegrable functions capturing the parameter-dependen causal effects, and e; are mutually independent, centered, sub- Gaussiann noise terms. The underlying causal structure forms DAG G (V, E), where V = (1....,m) n) and an edge (j,t) E E exists if and only if j E PA,, We make the standard assumptions of causal sufficiency (no unmeasured common causes) and faithfulness (all conditional independencies in the data are consequences of d-separation G). Throughout, edge set E does not depend on É: only the edge weights biy (E) vary with operating conditions. We assume m i.idd samples {(X(1), {()) from this model, where & indexes samples rather than time.
2.2 POLYNOMIAL. CHAOS REPRESENTATION
The central challenge is to learn the functions (E)(E), the true causal strength varying with system parameters ç. We address this by representing each causal coefficient function using PCE. For many choices of [Ag. including the classical Wiener-Askey scheme, there exists a corresponding basis of orthogonal multivariate polynomials (V) adapted to μr (e.g.. Xiu & Karniadakis, 2002; Sudret, 2008; Crestaux et al., 2009). Common examples include Hermite polynomials for Gaussian parameters, Legendre for uniform and Laguerre for exponential.
Any square-integrabl function (5) can be expanded in this basis. By truncating the expansion at a total polynomial degree Np. we obtain a finite-dimensional approximation:
(E)(E) ~ Σ
Ozj.a Va(EE)
OEANP
where AN, == fa E N |a| = Ei=1ar a I N2} } the set of multi-indices, and the coefficients Ozj.a are the spectral projections of the function onto the basis, given by Ozjaa (Et)(E), (E)))) r... The cardinality of the basis is P AN, (No+t)). This representa- tion converts the infinite-dimensional problem of learning functions big (E) into finite-dimensional problem of estimating the spectral coefficients lij.a.
For functions that are continuously differentiable s times, the spectral error decays polynomi- ally: [bbi - 1I, k < CN,² where IIN, is the projection operator. For analytic func- tions, which are common in physical systems, convergence is exponential: l|brg - II,, biglls <
C Cexp(->NV (e.g., Sudret, 2008; Crestaux et al., 2009). For nigh-dimensional parameter spaces (d >> 1). the basis size P can become computationally prohibitive. We employ hyperbolic trun- cation schemes, which prioritize low-order interaction terms and significantly reduce the basis size while often retaining high accuracy for functions with decaying importance of higher-order interac- tions (Jiang et al., 2025; Dai et al., 2025).
2.3 PCT-CONDITIONAL INDEPENDENCE TEST
Once have established the PCE representation, our initial goal is to identify the causal skele- ton. Standard conditional independence (CI) tests that operate on the marginal distribution of XB, .z) can fail in the presence of parameter-varying mechanisms: a causal relationship may be strong for many values of but cancel out in the marginal (c.g., sign-changing effects), lead- ing marginal CI tests to falsely conclude independence. We therefore test conditional dependence as a function of the operating parameters.
Let
CABIZ(E) = Cov Xn Xz, denote the conditional covariance function over the parameter space.
Definition 1 (PCT-conditional independence). We say that XA and X8 are PCT-conditionally in- dependent given Xz if
||CABiZllLA(u))= = B [C(CAAIZ)(E) = =0,
(3)
equivalently, (E) = 0 for pe-almost every &
4

---

Published as a conference paper at ICLR 2026
Under mild moment conditions, CA) L' (pgg sion
admits (truncated) polynomial chaos expan-
CABIZ(E) =
Σ CABIZ (E)
P:= P:=|AN,
(4)
cGEAN,
By orthogonality of ( J. the null hypothesis is equivalent (under truncation) to the vector of PCE coefficients being zero:
Ho CABIZ= (CABiz,alaAAN, 0.
To estimate CABIZ(E), first remove the dependence (Xz; by regressing XA and X on the interaction featurrs {Xx V. (XxVa(E))ae2) (E))
BA. x v v +
(2) A
X
Σ Σ
BB (
kEZ AEAAN
KEZ aE,AN,
The residual product (t) (t)
A B|2,c is then an empirical estimate of the conditional covariance signal (()), Define the basis vector (E) = (Va(E))EEAN, ER", and
m
(+0) :0 (1)
AlZ.E Biz, wb (5(()) R.P. Ō:
(t)
m
Let Σ be the sample covariance of [v(t))m
It=l"
Σ =
- -C)
Under Ho and standard regularity conditions (i.i.d. samples conditional on &, finite fourth moments, and nonsingularity of Σ), multivariate CLT yields vm Ò NN(0, E). Consequently, the Wald statistic
TPCT == m CT E C S xar-
E := Σ + EregI,
(5)
where df rank(E) (typically = P) and Ereg > is a small ridge term for numerical stability. We use this test the CI oracle inside PC-style skeleton search (Algorithm I).
2.4 SCORE-BASED LEARNING WITH NATURAL GRADIENT
Although constraint- based methods are effective for skeleton discovery, they can be unstable with finite data. We therefore use the output of the constraint-based phase initialization for more robust score-basee based optimization. We formulate structure learning as the optimization of penalized likelihood score over the space DAGs and PCE coefficients and define the PCT-BIC score as:
S(E
m
X(()
Σ JEPA; AEAN,
Σ
Big V ((()))) )
+ AII(EE, O)
(6)
20t
where ④ = ((...aa) is the collection of all PCE coefficients. The group sparsity penalty (E θ) Ea, 1[1 l2 > 0} encourages sparse DAGs by penalizing the number of non-zero causal links. Here >>0 is a regularization weight and in practice we treat as a tunable hyperparameter.
Optimizing this score is challenging due to the combinatorial nature of the graph space and the high dimensionality of the parameter space e We employ a greedy search strategy combined with efficient gradient-basee optimization of the coefficients for given graph structure, accepting edge additions or deletions only when they preserve acyclicity of g.
Natural-gradient preconditioned updates for coefficient optimization. For a fixed DAG G (hence fixed parent sets (PA,)). the sparsity penalty in equation 6 is constant, and optimizing over θ reduces to n separate least-squares problems, one per node i. For each node define the design matrix &; E (Mx((PAA)))) with entries
(Dll.taa)) = x Vo E PA; E AN,:
5

---

Published as a conference paper at ICLR 2026
and let 0, I PA.|P stack (0ij,a)jEPA. Then the data-fit fit term for node is
((
1 x
0.0:112
x; x: (x.)) p.... x: X(m))TT
20²
Under a Chnaa-Gausssin (or quasi-likelihood) view, the empirical Fisher information is
F
o.m $ 4
which generally not diagonal because the regressors depend on & and can be correlated. We therefore apply Pisher-preconndtioneed update (natural gradient)
0: - 1 ( Φ; Ereg 1 (Φ,0; x,),
with small ridge Ereg 0 for numerical stability and warm starts the graph is updated during greedy search.
The complete PCT-CD algorithm, summarized in Algorithm 2 in the appendix, integrates these components into a multi-phase procedure that ensures both structural accuracy and robust parameter estimation. This provides not only the final graph and functional relationships but also confidence intervals for causal strengths and probabilities for the existence of each edge.
THEORETICAL ANALYSIS
In this section, we establish the theoretical foundations the PCT-CD framework. We prove that under reasonable conditions, the true parametric causal DAG is uniquely identifiable from observa- tional data. Furthermore, we provide finite-sample guarantees for the recovery of the causal structure and analyze the convergence properties of optimization procedure. Formal statements and proofs of the main results are deferred to the Appendix for clarity.
3.1 ASSUMPTIONS AND PRELIMINARIES
Assumption (Data-generating process). The variables X (XI,..., obey the SEM equa- tion with a DAG G. where each biz(E) E (.3())) admits the chaos expansion equation 2 with trancation bias controlled by polynomial convergence theory. The parameter vector & is indepen- dent of noises {E:} and has a known distribution μę with finite moments of all orders.
Assumption (Noise). The disturbances e; are mutually independent, centered, sub- Sab-Gaussian with proxy of and finite fourth moments. Moreover, most one E, is Gaussian (equivalently, all but at most one are non- non-Gaussian).
Assumption 3 (Faithfulness and stability). For Pg-aammst every &, the conditional distribution
P(X 2) is faithful to the DAG G. The operator norm satisfies E[][B(())lop) < 1, where X B + e is the matrix form with [B(E)];j = bir(E) for j PA; and zero otherwise.
Assumption 4 (Design regularity). Fix truncation order Np and basis size P |Ax, and let be an upper bound on in-degree. For each node and any candidate parent set S, with {S, S s, define the interaction feature vector
Φ... := (X, Va((E) jES,, dEAx, E S...P
and the population Gram matrix Gt,s, El ST Assume there exists a constant y >> 0 such that
Amin (G,s.) 2 for all i and all S; considered by the algorithm.
Moreover, let Di.9. ppxx([S,)P) be the empirical design matrix whose t-th row is (B) the empirical Gram matrices concentrate that with probability least 1 5/4,
i... Assume
1 T 1.5,
Amin
$1.5
simultaneouully for all and all S, considered by the algorithm.
m

---

Published as a conference paper at ICLR 2026
3.2 IDENTIHIABILITY OF PARAMETRIC CAUSAL STRUCTURES
Identifiability is the cornerstone of any causal discovery method, ensuring that the underlying causal structure can, in principle, be recovered from the joint distribution of the observed variables. We extend the classical results of LiNGAM to our parametric setting.
Theorem 1 (PCT Identifiability. Under Assumptions 1-3, the DAG and the collection of co- efficient functions {Er are identifiable from the joint distributiin of (X, up to the usual permutation and scaling indeterminacies of independent component models. In particular, bey(E) is identifiable as an element of T.² (1E), and therefore its PCE coefficients
Ozjaa (bir(E), L((e))
LP))) uniquely determined by orthogonal projection for every AN,-
Remark (all-Gaussian case). If all e; are Gaussian, then for each fixed & the model reduces to GaussinG SEM, for which the DAG is in general identifiable only up to Markov equivalence without additional assumptions (e.g., heteroscedasticity across environments, multi-envvronment in- terventions, or other asymmetries). In this work we do claim full DAG identifiability in the all- all-Gaussian case.
3.3 FINITE-SAMPE GUARANTEES AND CONSISTENCY
While identifiability ensures recovery from the distribution, practical algorithms operate on finite data. We state a finite-sample guarantee for recovering the truncuted parametric mechanisms under fixed polynomial order Np-
Truncation and effective edge strength. Let NN, denote the ( -othhoonall projection onto (Nx)
the span of { a JoEANN, Define the truncated edge function big (E) (1),)))))EE. Let the truncation bias be
Strunc max |bg (U))E
(N+) (1ZA(H)))
and define the minimum effective (truncated) edge strength
KN, min [1b [[(())))
(Np)
(3,))
Group edge detection. O.j (0ij.a)EEAN, R by least-squares regression interaction features (E)). Define the estimated edge strength
For each candidate edge + estimate the truncated coefficient vector
( 10
(22((e)) = Σ
tj. (V22).
oEAN,
and declare (j,z) É iff |0 (()) for threshold F > 0.
Theorem 2 (Sample complexity for recovering truncated edge functions). Assume Assumptions 1-4 and let be the maximum in in-degree considered by the algorithm. Fix the truncation order Np (hence P |AN, and suppose KNp 0. Choose the threshold =KN,/2 Then there exists constant
C. > 0 such that if
m C o2 (sP) log
2nn2P 8
7 TNN
the trou--thresoolding rule recovers the correct edge set E (for the truncated model) with prob- 2etruee where min(..))
ability least 5. Moreover,
then KN, Etrune 2 k/2, and the same bound guarantees recovery of the true (non-truncated) edge set as well.
The bound captures the correct dependencies: sample complexity scales with the number of in- teraction features per node (sP), deteriorates with ill-conditioned designs (y small), and increases quadratically the weakest effective edge HN,, decreases.

---

Published as a conference paper at ICLR 2026
3.4 CONVERGENCE ANALYSIS
Finally, we analyze the convergence of the coefficient optimization phase, which is critical for the efficiency of the score-based refinement. The use of natural gradient descent key achieving rapid convergence.
Theorem (Natural-gradien (preconditioned LS) convergence). Fix DAG and consider optimiz- ing the coefficient-only objective for fixed graph. For each node define the least- ceast-squarss loss
(.
20 x. 2.0.1.
where 4, is the interaction-feature design matrix and 0; stacks the coefficients for node i. Assume $ V; is positive definite and 0 the unique minimizer.
Consider the (ridge- -precoodiiioooee update
0(1+2) =0! 11 (8, O. Erog Erog!) 1 "0, x:),
Grog 2 0.
Then the error evolves linearly:
0 ( 0 = (1, D1 + Eveg!)
(e, Di (2:) (0, 01).
(8)
In particular, for 0 < 1j the iteration converges linearly; if Ereg 0 and rj 1, it reaches the least-squares optimum in one step.
4 EXPERIMENTAL VALIDATION
We validate PCT-CD on a comprehensive industrial process dataset comprising 10,000 samples from chemical reactor network at Parkland Refinery in Canada. The system monitors critical process variables including feed temperatures, reactor pressures, product quality indicators, and flow rates, with 11 established causal relationships verified through process engineering principles and operational expertise. The system exhibits parametric uncertainty from three primary sources: heat transfer coefficients (L1), reaction rate constants (E2). and yield factors (Ez), making it ideal for demonstrating the advantages of modeling parameter-dependen causal relationships.
4.1 EXPERIMENTAL SETUP
The dataset represents a hierarchical chemical process where variables form network structure with source nodes (feed streams), intermediate processing stages (reactors and separators), and ter- minal outputs (product quality metrics). Each sample includes simultancous measurements of all process variables along with the corresponding parameter values, captured under varying opera- tional conditions over a 6-month period. The ground truth causal structure was established through combination of process flow diagrams, material balance equations, and expert knowledge from plant operators. We compare PCT-CD against 23 state-of-the-at methods spanning six categories.
selected through cross-validation: N, 4 (polynomial degree), Ossig 0.05 (significance level for conditional independence tests), X = (regularization parameter), and B 200 (bootstrap samples).
PCT-CD parameters
4.2 PERFORMANCE RESULTS
Table presents comprehensive performance metrics across all methods. PCT-CD achieves ex- ceptional performance with 90.9% F1-score, correctly identifying out of 11 true edges (True Positives) with only false positive and 1 false negative, resulting in structural Hamming distance (SHD) of 2. This represents nearly double the performance of the best bascline methods.
Analyzing the results by method category reveals systematic patterns. Constraint-based methods (PC, FCT) achieve moderate precision (55.6%) but suffer from low recall (45.5%), indicating con- servative edge detection that misses many true relationships. Score-based approaches (GES, GIES, NOTEARS) show balanced precision and recall around 50-54%, but their static graph assumption fundamentally limits performance. Among functional causal models, traditional I. LiNGAM performs
8

---

Published as a conference paper at ICLR 2026
Table 1: Performance Comparison Across All Methods
Method
TP
FP
FN
Prec.
Recall
F1
SHD
ICA-LiNGAM
1
14
10
0.067
0.091
0.077
24
DirectLiNGAM
2
13
9
0.133
0.182
0.154
22
VAR-LINGAM
3
7
8
0.300
0.273
0.286
15
RECI
9
7
0.308
0.364
0.333
16
PCMCI
4
8
7
0.333
0.364
0.348
15
CCD
5
11
6
0,312
0.455
0.370
17
LiNGAM
5
10
6
0.333
0.455
0.385
16
ElasticNet
5
8
6
0.385
0.455
0.417
14
Entropy-Based
5
6
6
0,455
0.455
0.455
12
GP-Based
5
5
6
0,500
0.455
0,476
NOTEARS
5
5
6
0.500
0.455
0.476
CGNN
7
10
4
0,412
0.636
0.500
14
Dasso-Granger
7
5
0,462
0.545
0.500
12
FCI
5
4
6
0.556
0.455
0.500
10
PC
4
6
0.556
0.455
0.500
10
ANM
6
6
5
0,500
0.545
0.522
PNL
6
5
0,500
0.545
0.522
GIES
6
5
5
0.545
0.545
0.545
10
GES
6
5
5
0.545
0.545
0.545
10
CAM
8
6
3
0.571
0.727
0.640
9
GraNDAG
8
6
3
0.571
0.727
0.640
9
SAM
8
6
3
0.571
0.727
0.640
9
PCT-CD
10
1
1
0.909
0.909
0.909
2
Table Parametric Causal Effects with Uncertainty Quantification
Edge
Mean
95% CI
Boot Prob
Dominant
XI + X2
0.642
[0.411. 0.873]
0.95
E1 (heat)
XI + X3
0.465
[0.305, 0.669]
0.88
& (reaction)
XI X4
0.359
[0.217, 0.499]
0.84
{1 (heat)
X2 + X5
0.275
[0.151, 0.401]
0.78
&1 (heat)
X2 X7
0.213
[0. 14.142,00 0.284]
0.89
Ç (heat)
X3 + X4
0.185
[0.098, 0.272]
0.85
&2 (reaction)
X3 X6
0.214
[0.081, 0.285]
0.81
Ea (reaction)
X4 + X7
0.176
[0.092, 0.244]
0.91
Ea (yield)
X5 X8
0.133
[0.071, 0.183]
0.82
& (heat)
X6 + X9
0.109
[0.034, 0.156]
0.93
&u (reaction)
X7 X9
0.452
[0.317, 0.632]
0.95
Es (yield)
poorly (38.5% Fl-score) while ICA-LiNGAM shows the worst performance (7.7% Fl-score), sug- gesting severe model misspecificationn under parametric variation. Figure 2 visualizes the discovered causal structures across 23 methods, providing comprehensive comparison of graph recovery quality. The performance gap between PCT-CD (90.9% F1 F1-score) and the next best methods (CAM, GraNDAG, SAM at 64.0%) highlights the value explicit uncertainty modeling.
4.3 PARAMETRIC UNCERTAINTY QUANTIFICATION
Table 2 reveals PCT-CD's unique capability to quantify how causal relationships vary with system parameters. Each edge's strength represented as a continuous function of the parameter vector with confidence intervals capturing both estimation uncertainty and parametric variation. The strongest relationship X1-XX2 varies by over 100% depending on heat transfer conditions, while constrained variation. Heat transfer coefficients (51) predominantly influ-
weaker edges show ence feed and thermal control pathways, reaction rate constants (5x) govern intermediate transfor- mations, and yield factors (Ex) control product quality paths.

---

Published as a conference paper at ICLR 2026
Deseneer
Wiiti Sdge
Figure 2: Discovered causal structures for all 23 methods. PCT-CD (top-left) accurately recovers the ground truth with minimal errors, while baseline methods show varying degrees of false positives (red edges) and false negatives (missing edges)
These results highlight how PCT-CD captures parameter-dependeer variations in causal strength. Rather than assigning single static weight, each edge represented as function of &, with con- fidence bands quantifying estimation and parametric uncertainty. This enables the method to dis- tinguish edges that are consistently strong (e.g., XI +22 under heat transfer variation) from those whose influence is highly context-specific (e.g.. X2-xX5 with non-monotonic dependence on &2).
From a methodological perspective, these results demonstrate that PCE-based representations low the discovery algorithm to recover not only the existence of edges but also their functional sensitivity to operating conditions. Such functional profiles provide interpretable evidence of how causal mechanisms adapt to process variability, feature not accessible to static graph models. This property is particularly important in industrial domains where safe control requires anticipating how interventions may propagate differently under changing parameters.
5 CONCLUSION
This paper addressed the critical limitation of static assumptions in industrial causal discovery by introducing framework to model dynamic, parameter-dependeer relationships. Our proposed method, PCT-CD, successfully learns these functional causal links, demonstrating superior perfor- mance with 90.9% F1-score real-world chemical process dataset. The core contribution lies in establishing theoretical identifiability for parametric causal structures and providing a robust algo- rithmic solution. This work provides a significant step towards building more realistic and reliable causal models for smart manufacturing, enabling enhanced process control and more accurate root cause analysis under varying operating conditions. Future research could extend this framework to handle unobserved confounders, incorporate more complex nonlinear interactions, and explore its application in online, adaptive control systems.
10

---

Published as a conference paper at ICLR 2026
REFERENCES
Liang Cao, Jingyi. Wang, Jianping. Su, Yi. Luo, Yankai. Cao, Richard D. Braatz, and R. Bhushan. Gopaluni. Comprehensive analysis on machine learning approaches for interpretable and stable soft sensors. IEEE Transactions Instrumentation and Measurement, 14:1-17, 2025.
Liang Cao, Jianping Su, and Yi Luo. Causal structure discovery for industrial process monitoring: A novel approach for nonlinear systems with complex noise. IEEE Internet of Things Journal, pp. 1-11, 2026. doi: 10. 0.11090000.200266633115
David M. Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3:507-554, 2002.
Diego Colombo and Marloes H. Maathuis. Order-independent constraint-based causal structure learning. Journal of Machine Learning Research, 15(1):3921-396,, 2014.
Diego Colombo, Marloes H. Maathuis, Markus Kalisch, and Thomas S. Richardson. Learning high- dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, 40(1):294-321, 2012.
Thomas Crestaux, Olivier Le Maître, and Jean-Mare Martinez. Polynomial chaos expansion for sensitivity analysis. Reliability Engineering & System Safery, 9417):111-1112,, 2009.
Chris Cundy, Aditya Grover, and Stefano Ermon. BCD nets: Scalable variational approaches for Bayesian causal discovery. In Proceedings the International Conference on Neural Information Processing Systems, pp. 7095-7110, 2021.
Wei Dai, Dewen Li, Hui Liu, and Yuelin Liu. A cost surrogate model for TSO-DSO coordination based polynomial chaos expansion. IEEE Transactions on Power Systems, 40(4):3580-3583 2025.
Qihang Fang, Gang Xiong, Mengchu Zhou, Tariku Sinshaw Tamir, Chao-bo Yan, Huaiyu Wu, Zhen Shen, and Fei- Yue Wang. Process monitoring. diagnosis and control additive manufacturing. IEEE Transactions Automation Science and Engineering, 21(1):1041-1067 2024.
Patrik O. Hoyer, Dominik Janzing, Joris M. Mooij. Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Proceedings of the International Conference Neural Information Processing Systems, pp. 689-696, 2008
Biwei Huang, Kun Zhang, Mingming Gong, and Clark Glymour. Causal discovery and forecasting in nonstationary environments with state-space models. In Proceedings of the 36th International Conference on Machine Learning, volume 97, pp. 2901-2910. PMLR, 2019.
Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O. Hoyer. Estimation of a structural vector autoregression model using non-Gaussianity. Journal of Machine Learning Research, 11: 1709 1731, 2010.
Haolin Jiang, Francesco Ferranti, and Giulio Antonini. Enhanced sparse polynomial chaos expan- sion for electromagnetic compatibility uncertainty quantification problems. IEEE Transactions on Electromagnetic Compatibility, 67(6)) 1751, 2025.
Zhanlin. Liu, Ashis. G. Banerjee, and Youngjun. Choe. Identifying the influential inputs for network output variance using sparse polynomial chaos expansion. IEEE Transactions on Automation 2021.
Science and Engineering,
Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. Dibs: Differentiable Bayesian structure learning. In Proceedings of the International Conference on Neural Infor- mation Processing Systems, pp. 24111-24123, 2021.
Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 15(1):2009-253, 2014.
11

---

Published as a conference paper at ICLR 2026
Joseph Ramsey, Madelyn Glymour, Ruben Sanchez-Romenn and Clark Glymour. A million vari- ables and more: The fast greedy equivalence search algorithm for learning high- high-dimensional graphical causal models, with application to functional magnetic resonance images. Interna- fional Journal of Data Science and Analytics, 3(2): 2017.
Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps- Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D. Mahecha, Jordi Muñoz-Marí, Egbert H. van Nes, Jonas Peters, Rick Quax, Markus Reichstein, Marten Scheffer, Bernhard Schölkopf, Peter Spirtes, George Sugihara, Jie Sun, Kun Zhang, and Jakob Zscheischler. Inferring causation from time series Earth system sciences, Nature Communications, 10:2553, 2019.
Shohei Shimizu, Patrik O. Hoyer, Aapo Hyvärinen, and Antti Kerminen. A linear moon-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7:2003-2030, 2006.
Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O. Hoyer, and Kenneth Bollen. DirectLiNGAM: A direct method for learning linear non- Gan-aassi structural equation model. Journal of Machine Learning Research, 12(33):1225-- 1248, 2011.
Le Song, Mladen Kolar, and Eric P. Xing. Time-varying dynamic Bayesian networks. In Proceedings of the International Conference on Neural Information Processing Systems. pp. 1732-1740, 2009.
Peter Spirtes, Christopher Meek, and Thomas Richardson. An algorithm for causal inference in the presence of latent variables and selection bias. In Computation, Causation, and Discovery, pp- 211-252. MIT Press, Cambridge, MA, 1999.
Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search. MIT Press, Cambridge, MA, 2nd edition, 2000.
Bruno Sudret. Global sensitivity analysis using polynomial chaos expansions. Reliability Engineer- ing& & System Safety, 93(7):964-979, 2008
Qingkai Sui, Yalin Wang, Chenliang Liu, Kai Wang, and Bei Sun. Attribution-aide nonlinear granger causality discovery method and its industrial application. IEEE Transactions Industrial Informatics, 21(8):6147-6157, 2025.
Xiaoming Wang, Chaobo Zhang, Yiqi Liu, Xiaojun Liang, Chunhua Yang, and Weihua Gui. Ad- vancing industrial process control with deep Geaming-enhancee model predictive control for non- linear time-delay systems. IEEE Transactions on Industrial Informatics, 21(9):6823-6833 2025.
Norbert Wiener. The homogeneous chaos. American Journal of Mathematics, 60(4):897-936, 1938.
Dongbin Xiu and George Em Karniadakis. The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24(2):6196444, 2002.
Xing Yang, Tian Lan, Hao Qiu, and Chen Zhang. Nonlinear causal discovery via dynamic latent variables. IEEE Transactions on Automation Science and Engineering, 22:10381-10391, 2025.
Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In Proceedings of the 36th International Conference on Muchine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 7 7154-7163. PMLR, 2019.
Kun Zhang and Aapo Hyvärinen. On the identifiability of the post-nonlinear causal model. In Proceedings of the 25th Conference on Uncertainty in Artificial Intelligence, pp. 647-655, 2009.
Yingwei Zhang, Yunpeng Fan, and Wenyou Du. Nonlinear process monitoring using regression and reconstruction method. IEEE Transactions Automation Science and Engineering. 13(3): 1343-1354, 2016.
Xun Zheng, Bryon Aragam, Pradeep Ravikumar, and Eric Xing. DAGs with no tears: Continuous optimization for structure learning. In Proceedings of the International Conference on Neural Information Processing Systems, pp. 9472-9483, 2018.
12

---

Published as a conference paper at ICLR 2026
Cheng Zhou, Kaibo Liu, Xi Zhang, Weidong Zhang, and Jianjun Shi. An automatic process mon- itoring method using recurrence plot in progressive stamping processes. IEEE Transactions on Automation Science and Engineering, (12):1102-1111, 2016.
Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning. In International Conference on Learning Representations, 2020.
A APPENDIX
In this appendix we collect proofs, algorithmic details, and additional experimental material, in- cluding the full setup and results for the synthetic benchmark with parameter-varying mechanisms discussed in Section 4.
A.1 PROOF SKETCH OF PCT IDENTIFIABILITY (THEOREM I)
We sketch the main argument and highlight the role of conditioning on the operating parameters. For Ieealmost every fixed value & = &, the SEM equation 1 becomes standard linear acyclic model
X B(E)X + E,
equivalently
X = A()E,, (E)) = (I - -B)
Assumption 3 (stability) guarantees I - (5) is invertible for pg-a.e. &, so A is well-defined.
Step 1: Identifiability for each fixed & (ICA/LINGAM). Under Assumption 2, the components of E are mutually independent and all but at most one are non non-Gaussian Therefore, for each fixed
ç. the linear mixing model X = e is identifialle up to permutation and scaling: if X d (g) with independent components é, then there exist permutation matrix P and an invertible diagonal matrix D(E) such that
(E) = PD((E,, == DPTT
This is standard consequence of ICA identifiability (c.g., via the Darmois-Skitovich theorem).
Step 2: The permutation/scaling cannot depend on E. Because & is independent of the noise vector e (Assumption 1). the distribution of does not change with &. If D were to vary with on a of positive me-measure, then the induced distribution of D(E) P e would necessarily vary with for any non-Gaussian component (scaling a non-Gaussian random variable changes its distribution), contradicting the requirement that the noise law is Ç-invariant and independent &. Hence DE)) must be constant (a.e. in E), and similarly the permutation P must be constant.
Step Recovering the DAG and the coefficient functions. For a DAG model, B({) strictly upper (or lower) triangular under topological ordering. and A (I B(E)) inherits corresponding triangular structure with ones on the diagonal. This structural constraint fixes the remaining constant permutation//saaling ambiguity, yielding unique A and hence unique B for /lc-a.e Ç. Therefore each edge function biy(E) is identifiable as an I.a (14g) function.
Step 4: Uniqueness of PCE coefficients. Once (E) E ((( is identified, its polynomial chaos coefficients are uniquely given by orthogonal projection:
(b:;(E) V. (E) E AN,.
Oija
(V2)
This completes the proof sketch.
A.2 PROOF SKETCH SAMPLE COMPLEXITY (THEOREM
Fix node and a candidate parent set S; with |S,| < Let t; E mRxx be the design matrix with (t) V, (E()),, and O, stack the corresponding d, = [.S|P columns corresponding to features (x) truncated coefficients. Write
y: = v.0, 6i,
6 (c /=1-
(1)
13

---

Published as a conference paper at ICLR 2026
The OLS estimator satisfies
0, 0; = -0, -
By Assumption 4, the population Gram matrix well-conditioned with Amin (alos,,, S,,,,) 2 Standard concentration implies that for sufficiently large Amin ( 0 Φ.) 2 2/2 with high probability, hence
my
Therefore,
2 1 7
110,
0,112
+ E
Since of are centered sub-Gart sub-Gaussian and independent of the regressors conditional on g(,, stan- dard concentration bound for sums of sub-Gaussian random vectors yields
-
1 E,
d 1
2n2 8
Φ m
S
C1
log m
with probability least - 6/(2.). A union bound over i I,..., gives the stated rate with d; S sP.
Each directed edge j corresponds to P-dimensional block 0, inside 0,, the same bound controls |0 - 0,y ll2 uniformly over at most n² candidate edges. If satisfies equation then this uniform error is at most /2, which implies exact recovery by thresholding at HN,/22 D
A.3 PROOF OF NATURAL GRADIENT CONVERGENCE (THEOREM 3)
The convergence analysis of natural gradient descent fundamentally differs from standard gradient methods due the incorporation of the Fisher information metric, which provides more appro- priate geometry for the parameter space. We establish the convergence rate by analyzing how the algorithm behaves in the Riemannian manifold defined by the Fisher matrix.
For a fixed DAG structure, the PCT-BIC score S()) becomes a quadratic function of the PCE coefficients θ The gradient in Euclidean space is:
Ox
Ves =
1
( - -x(( i=1
(9)
de
where (1(1) represents the model prediction. The natural gradient transforms this direction using the inverse Fisher matrix:
=F
(10)
The Fisher information matrix captures the local curvature the log-likelihood surface. It decom- poses into blocks indexed by edges, where the block for edge ji, has entries:
F Fiyla,a - EX v (.(E) (E)
(11)
Note that although the PCE basis satisfies Eg V (E) (E)) Serar ()) this orthogonality alone does not make the Fisher matrix diagonal, because X, depends on & through the parametric edge functions (E), so the expectation does not factorize as E[X,3] Ee [ Wor] ] in general. This is consistent with the main text observation and is precisely why the preconditioned update requires
solving linear system with ridge regularization rather than simple diagonal rescaling.
To establish the convergence rate, we analyze the evolution of the error in the Fisher norm. Let (!) (t) 2" denote the error iteration f. The natural gradient update yields:
14

---

Published as a conference paper at ICLR 2026
A (1+1) A -
(12)
Using the Taylor expansion of the gradient around O* and the fact that VeS(e*) 0:
Vg HA (4) O(IIA(1)33)
(13)
where H is the Hessian matrix at the optimum. For the quadratic objective arising from linear models, the Hessian is constant and equals H F (A), where the perturbation term comes from the regularization.
Substituting this into the update equation:
A (++1) (I- A
(14)
The strong convexity parameter μ and smoothness constant L in the original Euclidean metric trans- late to corresponding parameters paf and LE in the Fisher metric through the eigenvalue bounds:
H
LF
(15)
Amax (F)
Amin (F)
The spectral radius of the iteration matrix (I 11-1FF- H) determines the convergence rate. With the optimal step size 17" = 2/(Fr LF), we achieve:
LF MF
1-p p +
(I - 1f F HH)
(16)
LF
+ 4t
where p p= /F / T. Amin (F) (m/x is the condition number of the Fisher matrix.
Therefore, the error contracts at each iteration according to:
| |F 11-1 LF
|A
(17)
This linear convergence rate represents a significant improvement over standard gradient descent, whose convergence rate depends the condition number of the Hessian in Euclidean space. Al- though the Fisher/Gram matrix is generally not diagonal (since X, depends on {), the precondition- ing absorbs this ill-conditioning into the linear solve step, making natural gradient descent particu- larly effective for high-dimensional PCE coefficient estimation.
A.4 ALGORITHM 1: PCT--II TEST ALGORITHM
PCT Conditional Independence Test
Input:
Samples ((X(1, { (())) index sets A, B, Z, basis index set AN,, significance level Osig, ridge Greg >
Output: p-value for PCT conditional independence // 1) Build basis vectors
for t 1 to m do (6()) + (V.((""))) E R
end for // 2) Residualize conditioning on (Xz,E) {) via interaction features
Build design matrix Φz E (Iz(() with columns {x((² (E(( (E)
Fit OLS of XA on Φz to obtain residuals (1) (AAlz,g) t==
Fit (ridge-)OLS of X on 2 to obtain residuals (t) 1t
3) Form moment vectors and estimate mean/covariance
15

---

Published as a conference paper at ICLR 2026
for 1 to do
v(0) t AlZ B|2 20
(t)
e
ERP
end for mm
C- m -1=1 vl0)
E M-T I
Erog + Erog
1/4 4) Wald statistic and p-value
TPCT C C
df < rank(Sree)
Return: p-value- =1 Fx (Trct)
p-value < asig then
Reject Ho: dependencies detected
else
Accept Ho: conditionally independent
end
A.5 ALGORITHM PCT-CD MAIN ALGORITIM
PCT-CD: Polynomial Chaos Theory for Causal Discovery
Input: Data {(X (1), { parameters N,, X
Output: Final graph G, coefficients {0tz)- uncertainty measures
// Phase 1: Initial Structure Discovery
Use PCT-CI test obtain initial DAG Go // Phase 2: Score-Based Refinement Initialize with g < Go // Forward Phase while score improves
Find edge (i,j) that maximally improves PCT-BIC score
if adding (i,j) maintains acyclicity then
G < G U {(,
Re-optimize coefficients θ using natural gradient end if end while // Backward Phase
while score improves do
Find edge (i,j) whose removal maximally improves score (i,,5)
G G end while
// Phase 3: Edge Orientation Refinement
for each edge (i.j) e 9 do
Verify orientation using non Gaussianity/residuul methods end for
for each non-ad coj-adjacent pair (i.j) & do
Test for nonlinear relationship using MI and residual analysis
if criteria met AND acyclicity preserved then
Consider adding edge (i.j)
end end for
// Phase 4: Uncertainty Quantification
Generate B bootstrap samples from original data
for b 1 B do
Rerun Phases 1-3 bootstrap sample Obtain (b) and (b)
end for
Compute edge probabilities: P(i ()
EE 1((,5) E F(0)
16

---

Published as a conference paper at ICLR 2026
Table 3: Ablation Study Results
Configuration
Precision
Recall
Fl Score
Full PCT-CD
0.909
0.909
0.909
Without PCE
0.611
0.636
0.623
Without Multi-criteria
0.647
0.818
0.722
Without Bootstrap
0.769
0.727
0.747
PCE Order N₂ 2 0.667
0.545
0.600
PCE Order N, = 3
0.786
0.727
0.755
PCE Order N₂ = 4
0.909
0.909
0.909
PCE Order Nº = 5
0.846
0.818
0.832
Penalty À 0.01
0.611
0.818
0.700
Penalty À 0.1
0.733
0.818
0.773
Penalty À
0.909
0.909
0.909
Penalty À 10
0.857
0.727
0.787
Table 4: Computational Scaling Analysis
Variables
Samples
Runtime
Memory
Fl Score
9
10,000
42.38
892MB
0.909
20
10,000
4.2min
3.1GB
0.795
50
10,000
28.5min
9.8GB
0.741
100
10,000
2.3hr
24.2GB
0.698
9
1,000
4.8s
218MB
0.636
9
5,000
23.1s
564MB
0.773
9 10,000
42.3s
892MB
0.909
9
50,000
3.7min
3,8GB
0,945
100,000 9
7.2min
7.3GB
0.964
Compute confidence intervals for each Oij.a Calculate Sobol indices from PCE coefficients Return: G, (0.g edge probabilities, confidence intervals, Sobol indices
A.6 ABLATION STUDIES AND COMPUTATIONAL SCALING
Table quantifies each component's contribution overall performance. Removing PCE causes the largest performance drop (28.6% Fl decrease), confirming polynomial chaos representation as fundamental to capturing parametric uncertainty. Multi-criteria refinement improves precision from 64.7% to 90.9% by preventing false positives. Bootstrap uncertainty quantification contributes 16.2% Fl improvement through better threshold calibration.
Parameter sensitivity analysis reveals optimal settings: PCE order Np = 4 balances expressiveness and overfitting, while regularization PCE orders lack sufficient flexibility, while higher orders overfit given finite samples.
optimally trades model complexity against fit. Lower
Table evaluates scalability across different problem sizes and sample counts. Runtime scales quadratically with variable count and linearly with samples, remaining tractable for industrial ap- plications. Performance improves monotonically with sample size, reaching 96.4% F1-score at 100,000 samples, demonstrating effective utilization of large industrial datasets. The method scales to var-abbe systems 2.3 hours, confirming practical applicability to complex industrial pro- cesses.
17

---

Published as a conference paper at ICLR 2026
Precision Recall Trade-off
Predisoo
003
00
02
0.4
26
0.6
1.0
Recall
Figure 3: Precision-recall trade-off across all methods. PCT-CD achieves both high precision and recall simultaneously.
Parametric Causal Effectss 96%
1007
1336
32-77
23:00
83428
Hed Berder
Offect Strength
Figure 4: Forest plot of parametric causal effect strengths with 95% confidence intervals for each discovered edge.
A.7 ADDITIONAL EXPERIMENTAL RESULTS
Figure 3 illustrates the precision-recall trade-off across all methods. PCT-CD occupies unique position in the high-performance region (top-right), achieving both high precision (90.9%) and high recall (90.9%) simultaneously. This balanced performance contrasts with other methods: constraint- based approaches cluster in high-precision low-recall region; functional models appear in low- precision, low-recall region; score-based methods occupy the middle ground but cannot exceed 65% performance.
Figure 4 provides forest plot of parametric causal effect strengths with 95% confidence intervals for each discovered edge. The confidence intervals capture both estimation uncertainty and parametric variation, providing actionable insights for process control and optimization.
18

---

Published as a conference paper at ICLR 2026
A.8 LARGE LANGUAGE MODEL USAGE DISCLOSURE
We acknowledge the use of large language models to assist in grammar checking and language polishing throughout this manuscript.
19


Reviews:
[{"review_id": "f7BxfScJBq", "reviewer": "Reviewer_Vtcr", "summary": "Summary. The paper proposes PCT-CD, which models edge weights as functions of operating parameters via Polynomial Chaos Expansion, claims identifiability, and reports large gains on a single industrial dataset.", "questions": "1. You state that “traditional methods provide only point estimates and ignore uncertainty.” Which classes of methods are you referring to specifically, and which uncertainty-aware baselines did you exclude and why?\n\n2. Is the distribution of the operating parameters ξ assumed known or estimated from data? If estimated, how sensitive is your method to misspecification?\n\n3. Does the causal graph change with ξ (structure-varying), or are only edge weights functional in ξ while the graph is fixed? Clarify the formal model class.", "limitations": "", "rating": 2, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "Strengths.\n- Clear statement of goal: parameter-dependent causal effects. \n\n- Simple decomposition of pipeline: CI testing in the PCE basis, then score-based refinement. \n\n- Readability: paper is generally well written and easy to follow.", "weaknesses": "1. The introduction asserts that existing score-based methods “provide only point estimates without uncertainty quantification,” and that constraint-based methods relying on independence tests are insufficient in practice. These are too narrow and partially incorrect. There is a long line of Bayesian and bootstrap approaches for DAG posteriors and uncertainty over graphs and parameters; many constraint-based methods are paired with robust CI procedures and stability devices. The paper needs careful scoping and citations when criticizing “traditional” methods.\n\n2. The text (line 185) suggests that “traditional CI tests… are insufficient,” then replaces them with a bespoke PCT-CI definition without an empirical comparison against strong CI tests matched to data characteristics. The use of  PCT-CI is not justified methodologically. The argument reads as assertion rather than evidence. Provide comparisons to well-tuned kernel-based CI, partial correlation with robust estimation, or recent conditional mutual information estimators, on synthetic and real settings.\n\n3. Key background on parameter-varying or context-specific causal discovery, time-varying DAGs, covariate-dependent SEMs, and distribution shift causal structure is largely missing. The PCE/QoI citations cluster in early classics (Wiener 1938; Xiu & Karniadakis 2002) and very recent 2024–2025 engineering pieces, leaving a gap spanning two decades. This pattern suggests coverage of a few base methods plus several very recent items rather than a thorough survey. Expand the related work to include context-specific independence, conditional DAGs, regime-switching causal models, and Bayesian DAG posteriors with uncertainty. \n\n4. All results come from one proprietary refinery dataset with 9 variables and 11 asserted ground-truth edges. This is neither diverse nor standard. The head-to-head table mixes methods with very different assumptions and tuning needs, without transparent hyperparameter search, preprocessing, or split protocol. Strong directed baselines might perform well on such process-control data if tuned with domain knowledge; a single domain does not support claims of general superiority. Use public benchmarks (synthetic with controlled parametric variation, cause-me style datasets, time-varying DAG suites), release code, and report robust model selection.\n\n5. The comparison list includes methods that assume static graphs and others that model nonlinearity, but the paper then generalizes conclusions about “traditional methods.” If the claim is about advantages under parametric variation, show controlled synthetic experiments where the ground-truth edge functions vary with ξ, sweep PCE order and noise types, and compare to alternatives explicitly designed for nonstationary or context-specific structure. Current evidence is insufficient to support broad claims.\n\nMinor issues:\n- Typographical placeholders (“Algorithm ??”) and inconsistent section cross-references. \n\n- Ambiguity about whether ξ’s joint law is known a priori or estimated, and how misspecification impacts tests and scores. \n\n- The statement that score-based methods “provide only point estimates” ignores common uncertainty add-ons; rephrase and cite carefully.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "mvLkAjPzHW", "reviewer": "Reviewer_uJS9", "summary": "The paper presents a method for learning causal relationships between a\ncollection of measured process variables, where the strength of causal links\ndepend on uncertain parameters.\nThis type of model is relevant in process control applications, for instance,\nwhere operating conditions can change the causal graph of process variables.\nThe authors propose to model the relationships between each process variable\nand its `causal parents' as linear equations, where the coefficients are\nfunctions of the parameters. They use Polynomial Chaos Expansion (PCE)\nto represent the coefficient functions.\nBased on this model, the authors develop a conditional independence test used\nto initialise the learned causal graph. This is improved using a natural\ngradient method to optimise a likelihood score.\nGradient updates are weighted with the Fisher matrix to improve convergence.\nOn the theoretical front, the authors provide an identifiability theorem, a\nsample complexity bound and an analysis of the convergence of the gradient\nmethod, using standard assumptions.\nFinally, the method is validated on an industrial dataset, where it achieves\nstate-of-the-art performance across a variety of standard metrics.", "questions": "* The wording of Assumption 2 is unclear. In \"or the collection of coefficient\n  functions is non-degenerate\", what is meant by \"or\"? Is this an equivalent\n  reformulation of the previous statement, or an alternative but orthogonal\n  assumption that will lead to the same results?\n\n* Related to the above, in the proof of Thm. 1, it seems to me that\n  non-degeneracy is explicitly used to derive eq. (25). However, the text refers\n  to Assumption 1, the relevance of which I am not able to see here. Is this a\n  typo?\n\n* Regarding eq. (15) in the proof of Thm. 1: is it clear that the elements of\n  $A$ are in $L^2(\\Xi)$, so that it can be expanded in this way?\n\n* The text of the proof of Theorem 2 could be improved, e.g., by writing an\n  expression for $\\hat{\\theta}_{ij,\\alpha}$; avoid writing that the expectation\n  is \"approximately $P\\sigma_\\epsilon / m$, being more explicit; do not redefine\n  $\\kappa$ but refer back to Section 3.3.\n\n* In Section 2.4, it is written that the Fisher matrix is block-diagonal, but in\n  the proof of Thm. 3, this is refined to diagonal. This should be corrected for\n  consistency. More generally, some more details on the derivation of the\n  expression for the Fisher matrix would be appreciated.\n\n* Please provide references in Section 2.2 for the principal results\n  mentioned there, such as the choice of basis depending on the distribution,\n  the decay of the spectral error, and the hyperbolic truncation schemes.\n\n* Note the several broken references to the Algorithm environment.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "* Original formulation of the problem in a new setting, together with a\n  computationally tractable method. The integration of PCE is original and\n  appears to be a very good fit for the addressed problem.\n\n* Comprehensive theoretical analysis of the proposed method.\n\n* Experimental validation demonstrates strong performance compared to\n  alternative methods in the literature.", "weaknesses": "* The assumption of the parameters $\\xi$ having a known distribution seems quite\n  restrictive; the paper would benefit from a discussion of this and possible\n  workarounds.\n\n* As per the paper the polynomial basis functions used in PCE are dependent on\n  the distribution of $\\xi$. Section 4 does not discuss what assumptions were in\n  place regarding the distribution of $\\xi$, or in fact which basis was used for\n  the experimental results.\n\n* It's peculiar that the main algorithm, presented in appendix A.5, is not part\n  of the main body of the paper. Additionally, there are a lot of steps there\n  which are not discussed elsewhere in the paper: \"non-Gaussianity/residual\n  methods\", \"Test for nonlinear relationship using MI and residual analysis\"\n  (the abbreviation MI seems to be undefined), \"Consider adding edge\n  (i, j)\" (what is meant by \"consider\" here?). Unless these steps are explained, at\n  least in the appendix, the paper seems to lack reproducibility.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "fB7SAymuGR", "reviewer": "Reviewer_daxv", "summary": "The authors leverage polynomial chaos expansion to identify time-varying causal graphs with parametric uncertainty. Polynomial chaos expansion appears to be an adequate tool for a relevant but relatively understudied topic in causal inference. Therefore, the paper certainly has merit. The proposed approach outperforms several standard baselines.", "questions": "1) What are the \"standard regularity conditions\" that you assume for the independence test?\n\n2) From where does it follow that under these standard regularity conditions, the estimates are asymptotically independent and follow a standard normal distribution?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1) The paper is overall well-written and easy to follow.\n\n2) The considered problem setting is clear and relevant, and nicely embedded into the bigger streams of causal inference literature.\n\n3) The usage of polynomial chaos expansion is innovative.\n\n4) The experiments show a clear advantage over existing methods.", "weaknesses": "1) While I agree that time-varying causal graphs are not studied very often, I think it is not fair to say that all existing methods are limited by the assumption of a static causal graph. There is prior work that addresses this problem setting and deserves discussion. For instance, work on dynamic Bayesian networks (Song et al., \"Time-varying dynamic Bayesian networks,\" NeurIPS 2009), or Huang et al., \"Causal discovery and forecasting in nonstationary environments with state-space models,\" ICML 2019, address similar problem settings. Thus, I think the claims should be toned down a bit as well.\n\n2) A similar comment on the claim that traditional causal discovery assumes that all uncertainty is epistemic and rises solely from finite data samples. There is a lot of work on additive noise models, where the uncertainty is clearly not purely epistemic, even if it is not about parameter uncertainty.\n\n3) Some claims when developing the independence test require a bit more justification, see questions. It would also be good to state that proofs can be found in the appendix.\n\n4) Comparing to some methods that themselves consider time-varying graphs would strengthen the evaluation.\n\n5) References to the algorithm are shown as ??", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "gooSOGRD0u", "reviewer": "Reviewer_WoqE", "summary": "This paper introduces a novel and elegant method for learning dynamic causal graphs, where the causal relationships change as a function of system parameters. The problem of discovering dynamic causal structures is both timely and challenging, and the proposed algorithm, which cleverly combines concepts from different families of causal discovery methods, represents a relevant and interesting contribution.\n\nWhile the core idea is strong and I am inclined to recommend acceptance, the paper in its current form could be further improved by addressing several smaller issues concerning ambiguity, adding more detailed discussions, and expanding the empirical evaluation. I believe that incorporating these suggestions will further increase the paper's impact and clarity.", "questions": "1.  The proposed method uses a hybrid approach: skeleton discovery via a constraint-based method, followed by a score-based method for orientation and parameterization. Is this two-stage process strictly necessary? Have you explored whether a purely constraint-based approach (e.g., a PC-style algorithm using your proposed conditional independence test) is theoretically sound or empirically viable?\n2.  A discussion on the limitations and potential failure modes of the proposed algorithm would be highly valuable for readers. Under what conditions (e.g., type of dynamics, data scarcity, violation of assumptions) would you expect the method to perform suboptimally?", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "*   **Novelty and Significance:** The paper tackles the important and under-explored problem of learning causal graphs that are not static but change dynamically with system parameters.\n*   **Elegant Methodological Contribution:** The proposed algorithm provides a novel synthesis of constraint-based and score-based causal discovery principles, creating a powerful and interesting hybrid approach.\n*   **Timeliness:** The research area is of high interest to the machine learning community, and this work is a valuable early contribution to a nascent field.\n*   **Theoretical rigor:** The authors put visible effort into providing a solid theoretical analysis of their proposed algorithm.", "weaknesses": "### **Areas for Improvement**\n\nThe paper is promising, but its clarity and empirical validation could be strengthened in several key areas.\n\n**1. Clarity and Methodological Precision:**\n\n*   **Missing Details and Citations:** Several sections are very dense and would benefit from additional detail and citations to support the claims made. This is particularly true for:\n    *   Line 161 onwards: A citation is needed for the statement made here.\n    *   Line 205: Please explicitly name the specific optimization algorithm used in your work.\n    *   Line 223: The origin or justification for this property is unclear. Please elaborate.\n    *   Section 220 onwards: This section is difficult to follow. Please expand on the concepts presented and provide citations for established results to improve comprehensibility.\n    *   The index `α'` appears to be used without a formal definition in the main text.\n\n*   **Ambiguity of Temporal Dimension:** There seems to be a discrepancy in the problem formulation. Equation (1) presents a static formulation without a time dimension. However, Equation (5) and later references to time-series causal discovery algorithms imply that the data `X` has a temporal component. Please clarify whether the proposed method is designed for time-series data, i.i.d. samples, or both. This is a crucial detail for understanding the method's scope and applicability.\n\n*   **Acyclicity Constraint:** Score-based causal discovery algorithms typically deploy an explicit acyclicity constraint in addition to a sparsity regularizer. The objective function only appears to include a sparsity penalty ($\\lambda \\Vert (E, \\Theta) \\Vert_0$). Could you clarify if an acyclicity constraint is used, and if so, how it is enforced? If not, please explain how the acyclicity of the learned graph is guaranteed.\n\n**2. Experimental Evaluation and Analysis:**\n\n*   **Baseline Performance Analysis:** In the experiments, the static baseline algorithms perform surprisingly poorly, even though the causal parameters in the dataset have a strictly positive range. In such a case, one might expect a static method to converge to average coefficients, which, however, should still represent the proper causal connections.  Could the authors comment on why the baselines fail so decisively? Is this potentially due to hyperparameter choices? A more in-depth discussion would improve the reliability of the empirical results.\n\n*   **Value of Synthetic Experiments:** While the chosen real-world benchmark is interesting, the paper would be strengthened by including synthetic experiments. This would allow for a more controlled, fine-grained analysis to:\n    *   Demonstrate precisely under which conditions (e.g., speed of parameter change, noise levels, sample size) the proposed method outperforms static baselines.\n    *   Identify the limitations and potential failure modes of the algorithm.\n\n*   **Additional Datasets:** The empirical evaluation currently relies on a single ground truth graph. To demonstrate broader applicability and robustness, the authors should consider evaluating on additional benchmarks. There are several well-established datasets that could be suitable, such as:\n    *   Time-series benchmarks: [4], [5]\n    *   I.I.D. sample benchmarks with varying contexts: [6], [7]\n\n*   **Missing Related Work:** The discussion of prior work could be expanded. While the field is new, there are some relevant works that should be cited and discussed, such as [2] and [3].\n\n### **Minor Weaknesses**\n\n*   **Terminology:** As suggested by [1], the authors might consider using \"Structural Causal Model\" (SCM) instead of \"Structural Equation Model\" (SEM) to avoid confusion with the different usage of the term in the social sciences.\n*   **Cross-Referencing:** Please add backlinks to algorithm listings (e.g., on lines 196, 234) and to proofs in the appendix to improve readability and navigation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)