Frequency-Domain Better than Time-Domain for Causal Structure Recovery in Dynamical Systems on Networks
Authors: Mohammed Tuhin Rana
Mishfad Shaikh Veedu
James Melbourne
Murti Salapaka
Venue: ICLR 2026 Poster (2026)
Keywords: Causal Inference
Wiener Filter
Fast Fourier Transform
Graphical Model
Primary area: causal reasoning

Abstract: Learning causal effects from data is a fundamental and well-studied problem across science, especially when the cause-effect relationship is static in nature. However, causal effect is less explored when there are dynamical dependencies, i.e., when dependencies exist between entities across time. In general, it is not possible to reconstruct the causal graph from data alone. The conventional static causal structure recovery algorithms employ tests such as the Fischer-z test and the chi-square test to assess the conditional independence (CI) of data which forms the basis for recovering Markov Equivalent Graphs (MEGs) wherein causal structure can be recovered partially. For data that are dynamically related, multivariate least square estimation, based on Wiener Filters (WFs) relying on second order statistics for estimating a data stream from other streams, provides a means of recovering influence structures of the directed network underlying the data. Here, WF based projections can be determined in time-domain or in frequency-domain; the question this article sets out to answer is which is better? Here, we obtain concentration bounds on the accuracy of the WF estimation in both time and frequency-based approaches. Exploiting the computation speed of Fast Fourier Transform (FFT), we establish that the frequency domain provides distinct advantages. Moreover, frequency domain projections involve complex numbers; we establish that the phase properties of the resulting estimates can be effectively leveraged for better recovery of the MEG in a large class of networks; the time-domain has no analogue of phase. Thus we report the "Wiener-Phase" algorithm provides the best accuracy as well as computational advantages. We validate the theoretical analysis with numerical results. Performance comparison with state of the art algorithms are also provided. Further, the proposed algorithms are validated on a real field dataset known as the "river-runoff" dataset collected from the online repository of CauseMe, and on measurement data from transistor based circuits.

Published as a conference paper at ICLR 2026
FREQUENCY-DOMAIN BETTER THAN TIME-DOMAIN FOR CAUSAL STRUCTURE RECOVERY IN DYNAMICAL SYSTEMS ON NETWORKS
Mohammed Tuhin Rana University of Minnesota rana0082 umn. edu
Mishfad Shaikh Veedu
Google LLC mishfadd@google. com
James Melbourne
Murti V. Salapaka
Centro de Investigación en Matemáticas james. me eecournneeciaatee mx
University of Minnesota murt isttunn edu
ABSTRACT
Learning causal effects from data is a fundamental and well-studied problem across science, especially when the cause-effect relationship static in nature. However, causal effect is less explored when there dynamical dependencies, i.e., when dependencies exist between entities across time. In general, reconstruct the causal graph from data alone. The conventional static causal structure recovery algorithms employ tests such as the Fischer-z test and the chi- square test assess the conditional independence (CI) of data which forms the basis for recovering Markov Equivalent Graphs (MEGs) wherein causal structure can be recovered partially. For data that are dynamically related, multivariate least square estimation, based on Wiener Filters (WFs) relying on second order statistics for estimating a data stream from other streams, provides means of recovering influence structures of the directed network underlying the data. Here, WF based projections can be determined time-domain or in frequency- domain; the question this article sets out to answer is which is better? Here, we obtain concentration bounds on the accuracy of the WF estimation in both time and rrequency-base approaches. Exploiting the computation speed of Fast Fourier Transform (FFT), we establish that the frequency domain provides distinct advantages. Moreover, frequency domain projections involve complex numbers; establish that the phase properties of the resulting estimates can be effectively leveraged for better recovery of the MEG in large class of networks; the time-domain has analogue of phase.
not possible
Thus we report the "Wiener-Phase" algorithm provides the best accuracy as well computational advantages. We validate the theoretical analysis with numerical results. Performance comparison with state of the art algorithms are also provided. Further, the proposed algorithms are validated on real field dataset known as the "river-runoff" dataset collected from the online repository of CauseMe, and on measurement data from transistor based circuits.
INTRODUCTION AND LITERATURE SURVEY
Causal identification from data is an active and important research area relevant to multiple domains including climate science Pérez-Suay Suay & Camps- Valls 2018), economics Carfi & Caristi (2008). neuroscience Ramirez- Villegas al. (022)) and biology Hu et al. (2018) Lu al. (2023). There is considerable prior art on causation especially when the interactions are static in nature (see Pearl et al. (2016), Peters et al. 2017), ,, Spirtes et al. (2000), Wren et al. (2022) and the references therein). the cases where the entities are modeled as random variables it is not possible in general to reconstruct the causal graph from observational data alone. Indeed, Markov equivalent graphs (MEGs) that capture the same set of conditional independence (CI) relations can be determined from data (Spirtes et al.) (2000). Some recent works on recovering a unique causal graph when the underlying data generative system is assumed to have more structure presented in Peters et al. (2013). Shimizu et al. (200)).

---

Published as a conference paper at ICLR 2026
Causal inference is more challenging in the presence of dynamical (across time) dependencies Costanzo & Yagan (020): Krishnan et al. (2023); Peters et al. (2013). One of the established techniques for identifying causation in systems with dynamics Granger causality Granger (969), Bennett& & Yu (020). which leverages temporal structure assuming delays are present in interactions. Such an assumption is rendered problematic when the data is collected at slower rates than the time- constants at which the dynamics evolve, thus precluding a number of practical scenarios. Another approach in handling dependencies across times to consider the series at every time instant as a random variable thereby mapping the problem to static version. The difficulties with such an approach stem from the lack of information on the size of the horizon in the past and future to be considered and the combinatorial explosion of the number of variables that entail Ghahramani (2006), Lohmann et al. (012)) In a class of approaches for unveiling dynamic dependencies, models of how the data is being generated is assumed and the causal graph structure is recovered via an estimation of model parameters. Examples of such models include vector auto regressive ( (VAR) models Krishnan et al. 2023), Peters et al. (2013), additive noise models (ANM) Costanzo & Yagan (022)). Peters et al. 2013) and neural network-based models (NNM) Morafaah et al. (2021). In another recent line of work, the main gist of the approach is to exploit the asymmetry in the relationships of a dependency and the inverse of the dependency. These asymmetries can result from nonlinear maps or from linear filters with inverses that have discernible differences that, for example, are characterizable via power spectral densities Besserve et al. (2022). Shajaaiisales et al. (2015). Here, it is to be noted that works such as Besserve et al. (2022) and Shajarisales et al. (2015) try to infer cause and effect relationship between two variables; however, the causal graph structure recovery entails determining how the influence flows from variable to another with the possibility of intermediate variables between the
cause and the effected variables.
Compared to the case when interactions are static, the literature on causal structure recovery for systems with dynamic interactions is sparse. Several recent works have bridged the causality literature (Spirtes et al. (2000)) with multivariate projections of dynamically related time-series data (oostanoo & Yagan (2020), Materassi & Salapaka| (2016), Materassi & Salapaka (2019). It was shown in Materassi & Salapaka (2012) that the Wiener Filter (WF)-based projections can be employed to recover the moral graph of entities interacting via linear time-invariant dynamics. Materassi & Salapaka (2013) proposed that the WF-based tests can in principle be employed for performing CI tests in linear dynamical systems. We emphasize that effective and implementable methodologies for causal inference in linear dynamical systems remain much less developed than their static dependency counterparts. Moreover, WF based projections are often determined in the time-domain; as we demonstrate TD based approaches can be computationally burdensome.
This article develops a framework to reconstruct the MEG for system with dynamical interactions using a frequency domain (FD) approach. The main contributions are: 1) It is shown with examples that the WF-based reconstruction techniques provide superior performance in systems with linear dynamical interaction when compared to the conventional static counterpart. For the implementation of causation algorithms such as PC and FCI [Sirrtes et al. (2000). WF-based CI tests are provided. 2) To reduce the computational complexity of the WF implementation, an FFT-based approach is proposed which is proved both theoretically and in simulations to have faster run times than TD based approaches. 3) The pairwise CI test-based algorithms such PC is combinatorial in nature, which is intractable in bigger networks. To tackle this issue a computationally attractive algorithm called Wiener-Phase is proposed to reconstruct the MEG, albeit under some assumptions that is applicable in systems such as power-grid, consensus dynamics, and chemical reactions. 4) For estimating Wiener coefficients from time-series data, a non-asymptotic concentration bound and sample complexity analysis are provided for both the time and FFT-basel approaches. The results show that the FFT-based approach is computationally superior to time-based approach by a factor of OLL2 log N), where L is the number of future and past samples used for WF TD estimation and N the FFT window length. The simulation results corroborate the conclusion that the FFT approach gives better reconstruction error as well.
Notations: IR. C, Z, N denote the set of real numbers, complex numbers, integers and natural numbers respectively: [N] (O..... N 1}; For any time signal, == (x(t) i E Z). X denotes the Fourier transform of x, X ((w) = b oc [n)) -jwm; for a stochastic processes, I, a (w) denotes the power spectral density (PSD) of i;; for complex number, x = IR + ix,, La denotes arctan (#r RR)); Ω denotes the interval (0, 2Tr] and Ωy denotes (0, 2 /N,N..... ......,,((() N - 1): for transfer function, X, X # 0 means X is not identically zero. In a given directed graph, Pa(i), Ch(i), Sp(1) denotes
2

---

Published as a conference paper at ICLR 2026
the set of parents, children, and spouses of node i respectively. el denotes the space of absolutely summable sequences.
2
PRELIMINARIES AND RESULTS ON WIENER FILTER COMPUTATIONS
2.1 LINEAR DYNAMIC INFLUENCE MODELS (LDIM)
Consider a network with n nodes, each node E {1.... n) having the time-series measurements,
T; == {ai(k) governed by the convolutional model,
x (k) = Σ h(1)x(k- - ekk),
(1)
where x(k) := (k) 2., (k)] E R", h(e) E Ruxn. and ( := [eı (k) (k)] E R" is
the vector of exogenous noise sources with e;, ej jointly wide sense stationary. That is, for every
i E {1, (1.....,)), n.). there exists an a E R such that [k]= = a, for every k; and (k)) (0) Bley(k- - (O) =: Ray(k- - f). Taking Discrete Time Fourier Transform (DTFT), equation can be represented using the following linear dynamical influence model (LDIM),
X(w) = H(w) (w) + B(w), Vi E S2
(2)
where X(w) = [w) (w) E C², and E(w) = (w)
(N E C". H(w) E
Cax is well posed transfer function, (w)) is invertible as well as every entry of (I-H(w))
is analytic, and (w( = 0 for every w E f2. Moreover, it is assumed that the power spectral density
matrix (PSDM) of e, +E(W)= = ExEZ Ree(k)e -iwk, is positive definite and diagonal for every BE Z
2.2 GRAPH DEFINITIONS
The structure induced by LDIM equation 2 can be represented using a graph G -V, where V = [I.....m] and E = {(u, v): Hes * 0}. that is, there exists an edge u + 0 in E if He + for some w E Q. If u - e exists E then u is called parent of v (u E Pa(v)) and v is called child of u (v E Ch(u)). If fuuii + exists in EL then u and are spouses (u E Sp(u)). If u Sp(v) but u & Pa(v) Ch(v) then u is a strict spouse of v. u is said to be a kin of in G, denoted, u E (u), if at least one of the following exist in E: u - v, u - u, u - i + v for some i E V. king(i) denotes the set of kins of j in G. A chain from node to node j is an ordered sequence of edges in E. ((lo, (1), (2), ... (cn-1, fn)). where lo (topology). (G), of a directed graph, G = (V, is an undirected graph; Skel(G) == (V,E) where an edge (u, v) E E if either (u, v) E E or (v, u) E E. A path from node i to node j in the directed graph G is an ordered set of edges ((lg, (1), (f,, f,.....
ln = j. and (4,, (411) E E. Skeleton
(fn-1, .o).....--11 in its skeleton, EV, where {le,, lx+1)} E E. A path of the form ((ev, l)...., (fn-1. collider at Eg E - la + Ex+1 exists in G. Consider disjoint sets X, X,Y,Z C V in a directed
en) with to = i, l,
has
graph G = EV. Then, X and Y are d-separated given 7 in G. denoted d-sepa (X, 2, Y) if and only if every path between E X and y Y satisfies at least one of the following:
The path contains a non-collider node: = E Z. 2) The path contains a collider node w such that neither w nor the descendants of w are present in Z. In-nodes of i are given by the set {j (i. E E).
2.3 WIENER FILTER- BASEDD GRAPH/SKELETON LEARNING
Here we first present a method for estimating the 7th time-series from another set of time-series.
Let x = (=(1)= : E Z) and y = (v(t) E Z} be p and q dimensional random processes with xz{(= = [z) (1) (t) ... (t)] and y(() = [yn (1) y22 (2) ... (t)) Consider the problem of
finding the minimum mean square error (MMSE) estimate of y given x, defined by:
T org
min
E [ll((n) (n)) n
ã- T (m) ... [T(]] B, j=1....9
where M : spant (2z(k) :
I,...,,, -0 < k < oo The solution, 7n) = T7 (n) T2 [n(m)] is given by the Wiener filter (WF) Kailath al. (2000), described by: (i(n) = Σx 0o (k)x(((n - k),
1,...,4, where hT(k) E R8 are the WF co- (w)e
efficients, with (x(k) = + 2daH & (w) (w)) dw; y == -o--oo Ry (k)e--ww (k) = (k)r [0)]- Alternately, the optimal solution is described in the FD by P(w)
3

---

Published as a conference paper at ICLR 2026
yz (w) oz (w) (w), V w Ω, where Y(w) is the Fourier transform of 7. The solution is unique Paz (w) is positive definite for every w E Q. In the above description if we let C' S VV {i} and instantiate y = I2 and z z0, then the MMSE estimator, E, of the process 24. given the time-series Te satisfies (()
Wi.c(w)XX((4 where Wcc(w)
a (=) ( Here (w is
the multivariate WF obtained while projecting the time-series 2; onto xe Kailath et al. (200)) The entry of (w.)(w) corresponding to jth time-serie in the set C' is denoted by W (.ccl)(w).
2.4 TD APPROACH TO DETERMINING WF
Next, we briefly outline the conventional TD-based approach of estimating the multivariate WF
coefficients from finite data. Given the time-series data {x(t)] consider the projection of (t)
to the present, the past L values, and the future values of xz (4), j E
C C V (i) with m = Cl.
That is, for every 1 = I...., T - I., project (1(()) to span (2rg (k) 1. - I. S S + I,j E C). (T.L) A TD estimate of the WF, w 5.C E (R(21+1))x)) can be computed using the following least square formulation,
8 = arg
min BERm)LL+1) T+1 ||x; ycB12
where x; (T: L) [T L 1) ... 2L(l and ye is given by
Fay (T)
Se (T
(T- - 21)
(T)
Fum (2.) Zom (T 122--- -)
Ke (T 1)
seg (T 2)
mzz (T - 2L 1)
Fem (T -1
Yc
(5)
Fey (2L)
Tay (21- 1)
ecy(0) (0)
Zom (2L)
Femm (0)
The WF estimate, w (T.L)
C.C obtained by reshaping 8 into a matrix of the dimension (2L + 1) x m,
where the i-th (2L + 1) entries correspond to the coefficient of c₁ Without loss of generality, the
rows of whe are indexed from -L to L and the columns from to m. The Fourier transform of
(T.L)
Wic W (T.L))
(T.L)
(T.L) Fc (w) E C" is w (w) Σt - W.C. (kee- where w.c (k) is k-th row
(T.L)
of W.c The sample complexity of estimating W (T..L as provided in Theorem 5.1
(T.L),
&.C' (w) using equation scales as T = O (Im)
2.5 FFT-BASED COMPUTATION OF WF
Another approach for determining the WF coefficients is to first transform every time-series to its Fast Fourier Transform (FFT) representation, followed by projections in the FD. This approach involves sampling Ω [0, 2m] at N, equally spaced points, l2 = fwo. ... ..........Nl.. wN- where W/: == 2nk with N = 20, a E N, as shown in Appendix D This process of transforming the time-series to FFT representation is as follows. Consider the 7th time-series x; which is partitioned to R segments, each of length N, with the ,t segment denoted by 2. == {x:((r N),,, 2;((r- 1)........ (;((r- 1)N + N N--1) 1)}. -1)]. Thus, the time-series v, is given by (t) N-1 Ir-i, where (t) = ((((--1)++t) +t). Using the 7hth segment of the x, trajectory, given by .00),....,7""
N
(N (
the FFT, ((K) = 1 IN-I x (n)e -wen is computed. Let C' = C2..,,
S V where V is the set of nodes. Let X(((wR) C be the vector X (wx) X (Wk) ... X (wx) ]_ obtained by stacking the Fourier coefficients from the 4hh segment of the time series in the set Let
vN
(x(w)) X! (wk)= ... [[)() E CR and let Xe (wx) == (wx)) ... (((w))] E Rxx.
where C V. Then, for any We E ΩN the least square formulation described by,
N [[L(()))) - (c(
w (wne == arg min BECiCI
(6)
computes an estimate of the WF coefficients in the FD (see Doddi et al. (2022) for details). Notice
that Doddi et al. (2022) provided the FD computation of the Lasso-based representation of the WF for a special case of our general model. In this article we emphasize on 1) the computational advantages of employing FFT, and 2) the non- non-penalized MMSE estimate of the WF. In Appendix D.I an explicit convergence result of FFT to DTFT is provided, which shows a convergence rate of 1/VN. The sample complexity of estimating w-(f) (wx) scales as O(Nn) as shown in Theorem 5.2
2.6 COMPUTATIONAL COMPLEXITY OF ESTIMATING WF
The first major contribution of the article is the characterization of the computational efficiency in computing the multivariate WFs in FD when compared to the TD projections. In the TD, for m |C),
4

---

Published as a conference paper at ICLR 2026
the dimension of y. in equation 4 is (T 2L 1 (2L 1)m. Here the complexity of computing the least square estimate for node using equation 4 is since L < <. On the other hand, computing the FFT for window size of N samples takes N log N computations (Oppenheim & Verghese (2017). Notice that we compute the FFT with N samples per segment, and so the number of effective samples in the regression computation of equation @ is R = T/N. Thus the dimension of Xe (wx) is (T/N) x m. Computing the Wiener co-efficient for a given wi using the regression equation 6 ttakes O T'mm
(2 Re O(Tm² L²),
DNN computations, thus total of O(Tm² log N) computations are required to compute the Wiener coefficient for a single frequency. It follows that the FFT-based computation provides O /log N) improvement compared TD computation in the asymptotic worst case computation complexity. For the settings with large number of samples (large 7) or with longer delays (large L) the improvement is significant.
IDENTIFICATION OF MARKOV EQUIVALENCE GRAPH OF SPS
In the previous preliminaries section we described two methods of determining the WF coefficients from finite data; the first via projections in TD and the second via projections in the FD. We would like to emphasize that the focus of the article is on detection rather than estimation wherein we are interested in determining whether WF coefficient is zero not. The WF coefficient being zero or not will be used to retrieve the causal structure corresponding to the underlying data generative process. In general, it is not possible to reconstruct the exact and the complete structure of the directed graph G, (the generative graph), associated with an LDIM from data alone without actively intervening. Instead, in many scenarios the best can do is retrieve the Markov Equivalence Graph (MEG) which is the set of graphs that satisfy the same conditional dependence property as that of the true generative DAG Ghoshal Honorio (2018).
3.1 REALIZING THE ESSENTIAL GRAPH OF SPs USING WIENER FILTERING
In the static DAGs, reconstruction of MEG is well explored Spirtes et al. 2000). One of the popular approaches in the static setting is the Peter-Clarke (PC) algoritm,, which is performed using pairwise conditional independence (CI) tests Kalisch & mühlman (2007). In static settings, tests such as the Fischer-z test and the chi-square test are performed to assess Ci. However, the static CI tests fail in the dynamic setting because of the temporal dependency in the time-series shown in Table in the experimental results section. However, WF-based techniques can be applied on time-series with dynamical dependencies to obtain better results, as suggested by the following result from Materassi & Salapaka (2013).
Lemma 3.1 Materassi & Salapaka (013)). Consider a well-posed LDIM given by equation Let i, E V and let Z V (.,3)- Let We.1.2) 2) be the WF coefficients of estimating the time- series I; from the time-series Tj and zz. i and j are d-separated given Z in G = (V,E), then We15,2] 0, for every w Q. The converse holds almost always.
Inspired by the above result, we implement modified PC algorithm, where the WF is employed to test the CI on data with dynamical dependencies, to reconstruct the MEG. This proposed algorithm, called Wiener-PC (W-PC) is provided in Algorithm [ where the CI is performed by applying Lemma
3.1 Recall that the WF can be computed using either equation or equation
The following are the notations used in Algorithm D q denotes the maximum d-separating set size. combinations(( D. k)) lists all the k-length tuples in D, in a sorted order, without any repeated elements present, and DS(i,) denotes the d-separating set identified for i and j.
3.2 COMPUTATIONAL COMPLEXITY OF W-PC
As shown in Kalisch & Bühlman (207)) PC algorithm requires O(n9) tests. From Section 2.6 we
have that computing W (w) using FFT takes a complexity of ( log N) Rs O log N), since CI < n. Thus the overall worst case complexity of W-PC when implemented in FD using equa- tion 6 is OTTn+33 log N). in contrast to the TD approach, which has the complexity O L²).
4 WIENER- PHASE ALGORITHM
In this section, we propose a computationally attractive algorithm to reconstruct the MEG by exploiting certain properties of the WF, which are not present in any prior static approaches. Towards
5

---

Published as a conference paper at ICLR 2026
this, provide some useful lemmas and preliminaries below. The following Lemma shows that the set of kins in G can be reconstructed from the support of WF.
Lemma 4.1 (Materassi & Salapaka (012)). Consider a well-posed LDIM given by equation 2 Let W.. where V {t}. be the WF coefficients of estimating the time-series Is from the time-series 2y If W. ==D((w)) + then E king(j). The converse holds almost always.
In many applications (see Remark 17 in Veedu et al. (2021)) p), is possible to retrieve (k)l(G) exactly from the imaginary part of the WF. The following assumptions are required to obtain this result. Assumption 1. If node & has multiple incoming edges in G, then for every pair of in-nodes i, of and for every w E S2 ( = 2Hzz (w), where i. k E V. Assumption 2. If H,y # 0 then (()))W)) o O, for every i,j E V and w E Ω.
Assumption T says that the phase angle of all the incoming edges must be the same. This assumption seems restrictive; surprisingly a large class of problems satisfy this assumption. Applications satisfying Assumption are provided the supplementary material. Assumption is required for
the consistency of the reconstruction of the skeleton.
Lemma 4.2. Veedu et al. (2021) Consider a well-posed LDIM given by equation and satisfying
Assumptions 1 and 2 For any E Ω if ([W((w))) + 0 then either (i.5) E E or (j,t) E E. The
converse holds almost always.
Thus, by analyzing the real and imag- inary part of Wx, separately, one can extract the strict spouse edges. This information can in turn be em- G
Algorithm Wiener-PCC algorithm
Input: Data X(w)), wk E f2N, 7
Output: &
ployed to identify the colliders in an efficient way. We call this al- gorithm (Algorithm 2 Wiener-Phase here. The worst case asymptotic computational complexity of Wiener- Phase algorithm is (n (T log N (²)) (see Section 4.1.. which is ad- vantageous, especially in highly con- nected graphs (with large 4). We demonstrate the advantages of the Wiener-phase algorithm in 6
1. Initialize the sets, 8 ffuly connected undirected edge Col (), DS + 0, Fest ()
2. For i 1,...,
(a) For j = I,..., .,. n, and j i. Initialize D V {i,j}
ii For k = 0......
For in combinations(D) (D,
If |w(/ 115.a| "1b,a/Lil [il(we)
DS(i. + z
delete (i,j) from S
break
4.1 ANALYSIS OF THE WIENER-PHASE ALGORITHM
If |2| >= q
DS(i,)) + -1
3. For (i,5) E V V and < <.j:
In here, we analyze the Wiener-Phase algorithm (Algorithm 2 As shown in Lemma 4.2 in many applications Veedu et al. (202T), support of the imaginary part of the frequency depen- dent WF retrieves (k) exactly. Further, as shown in Lemma [. the support of W retrieves the Markov blanket structure in G. Combin- ing both, the Wiener-Phass algorithm identifies the colliders and thus the MEGs efficiently. Here, the of kin edges obtained using Lemma 4.1 and S the skeleton obtained from Lemma 4.2 Consider any edge (i, j) that belongs to but does not belong to S. Then i and j share common child without a direct link between them (thus are strict spouses) in G.
If E
(i, +- Franii 1
4. For (i.j) E x V and i <
For k =j I....,n
i. If (i, E & e (j.4}
1f II > -1& DS(j,))
Col C
E + Fest -
5. For every undirected edges, i.c., if (i,j) E Fast and (j.t) E Eest. fix the orientation of as many possible such that 1) there directed cycles 2) there no collider.
Return G (V, Frea.)
Following this procedure, can construct all the striet

---

Published as a conference paper at ICLR 2026
spouses. Let this set be SP. Recall from Section 2.6 that the complexity of computing W
1-j (w) is
Onn N). Repeating this process n times gives In N). Thus, the total complexity in finding K and S using FFT approach is On log N).
Now consider the skeleton S. In Step 4. we identify the common child be- tween the strict spouses i as follows. For any {i,j] SP, let C.y
Algorithm 2 Wiener-Phase algorithm
Input: Data X(w), wy f2
(k
Output: G
ki,k} E kj E S}, which is a set that contains all nodes k which has a link to both i and j in the skeleton. Note that the set Cip can contain non-colliders, which will be eliminated using the fol- lowing steps. For any (i.j} E SP. if [Cyy = 1, then c E Ciz is col- lider (because the link between and j in K. is formed by at least one col- lider). If [Cayl > 1, then for every E Caj. we compute w(()) (.)..]lo).
Initialize the ordering, +- ()
2. For 1,...,,.
(a) Compute w (w) using equation 6
(b) For 1,...,m
[ [I( > v, then K (1.})
ii If (3 o7 ([(m)) > T, then Su
3. Compute SP := KC \
4. For {i,}} E SP
If w +Gie)] * O, then e is a collider since [i,j] & S. The complexity of computing IC and S Onn N). Step 4 is repeated (nzq) times and Step 4(a) which checks for potential collid- among the common neighbors of and j takes ((n2g) operations. Steps 4(b) and 4(c) take (n)). Thus the com- plexity in computing Step 4 is O(n2q²) giving the overall complexity of O( (T N + 9²)). is sufficient compute Algorithm for O(1) number of w. Thus, the final complexity remains same.
(a) For k
1,...,7 i If (ii, ki] E and kj e
Cis t Cis [k}
(b) If Cf, = 1: Col + Col U {k}
(c) Else: for Cij compute wi() [-(A.c]
i. If W +0:4 > T: Col + Col
5. Return G
SAMPLE COMPLEXITY ANALYSIS
So far, have discussed about estimation of WFs using TD and FD approaches and its in causal inference algorithms. However, it is to be noted that in practice, due to finite sample effects, the estimated WFs using either equation or equation 6 will deviate from the true WF shown in equation 3 A quantification of these errors IS necessary for better understanding and implementation of WFbbased algorithms. In this section, we characterize the error in estimating the WF coefficients from data, and consequently, the number of samples required to estimate the WF reliably. The first step in bounding the error in the WF estimation bounding the error in estimating the PSDM.
The error in estimating the PSDM in TD is special case of the Blackman-Tukee estimator (Veddu jet al. a (202T) Lamperski (2023)) They compute the Fourier transform of the autocorrelation estimates to obtain the PSDM estimates. The following theorem provides sample complexity for estimating Wio from equation 4 by using the PSDM concentration bound from Lamperski (003). The proof provided in Appendix
Theorem 5.1. Consider a linear dynamical system governed by equation
Suppose that the auto- correlation function Rx (k) satisfies exponential decay, R| (k)] < Cp- and that there exists M
2
such that + S Amin (2x) S Amaz (2x)= & M. Then for any < e, 8 and L 2 log, 20
(O-plc and
M
Ω if
T T 0 (6n log(5)) ((2L 1)M*)
-.). -1.). then Yo E Q, P (IWie W > <8.
The following theorem, based on the PSDM estimator in Veedu al. 2024 provides the sample complexity for estimating Wic from e The proof is provided in the Appendix B Theorem 5.2. Consider a linear dynamical system governed by equation 2 Suppose that the auto- correlation function R (k) satisfies exponential decay, [x)]])]] Cp and that there exists
7

---

Published as a conference paper at ICLR 2026
M such that mmn (xx) I Amar ($x) S M. Let 0 E1,62 be such that E2 E1-
Suppose that N [1-p-1 "Yes"
20p
Then Vw E Ω and 0 < 8, if
T TEO 0 (6n (6n-loge log(E)) NM
then Ve E Q, P (Weolw) w > >e) <
g²
Comparing Theorem 5.1 with Theorem 5.2 it can be observed that number of samples required for given error performance is comparable when N L. That is, T Q(Ta) for the TD whereas = (NN) for the FD computation. In addition, simulation results in Section 66 show that the FI) approach shows better reconstruction accuracy than TD approach smaller range of samples.
EXPERIMENTAL RESULTS
In this section, we demonstrate that FFT-based approach provides superior computational advantages while having comparable reconstruction accuracy, if not better, in networks with dynamics using synthetic data as well as real field data.
The data generation and algorithm execution were performed in computer with intel core i9-1 1490KK 3200 Mhz, 24 cores, 32 logical processors. The generative model for the synthetic datasets were vector auto-regressive (VAR) model as shown below
ax((() + (1())xz( 1) + as (2) - 2) + (3( -0) 3) Σ biyzy(1 ec((), (7)
where i E {1, 2,... 6}. and e+(+) are zero mean i.i.d Gaussian noise with diagonal PSDM and ex(() ~ NO, 1). The coefficients
3D WFC
BD WPC WisFii
WisPPi
the VAR model were generated in accordance with the graphical structure of the corre- sponding DAGs i.e. bi was set to zero if link from to i does not exist in the corresponding DAG. The non-zero bi coefficients in the VAR model were randomly generated with bij (0.22 0.4), and a-(-) were fixed. A more detailed description of the data generation process is given in Appendix C For the FD based algorithms we averaged the magnitude and phase of the Wiener filters, obtained from data, over set of frequency y/ ). where N 2 E N. The average magnitudes were then compared to the chosen thresholds in the algorithms.
1.0-
Cocpperation
0.5
10
0.0
0.4
02-
107
102
101
10
Number of sample
Number of somplex
(a)
(b)
Figure Comparison of (a) average error rate (b) average computation time for 25 random networks generated accord- ing to equation
6.1 RANDOM NETWORKS
rate for 25 randomly generated DAGs with dynamics described by equation 7 for the three algorithms, namely, W-PC with TD-based WF estimation, W-PC with FD-based WF estimation, and Wiener-phase algorithm. Here, error rate is the ratio of the total number of wrong edges (both false positive and false negative) to the number of edges in the actual DAG. It can be observed that FD-based W-PC outperforms the TD-based W-PC in terms of average eтror. Moreover, the FD-based W-PC algorithm requires much smaller sample size compared to the TD-based W-PC algorithm for similar reconstruction accuracy, We note that the Wiener-phase algorithm also outperforms the TD-based W-PC algorithm. It can be observed that the average computation time for the TD-based W-PC algorithm significantly larger compared to the FD-based W-PC and Wiener-phas samples where there is around an order of magnitude difference between the computation time required for the TD-based W-PC and the FD based W-PC.
Fig [(a) shows the average
algorithms (see Fig. [(b).. This difference is stark at higher number of
To choose the thresholds for the randomly generated network we use a sparsity metric based approach. Bounds sparsity metrics such as average, maximum, and minimum degree are typically available apriori. As described in Appendix C.2, given a dataset and an algorithm, we use an iterative search method obtain the smallest value of thresholds that produce networks which satisfies the given bounds on average degree, and minimum degree.
8

---

Published as a conference paper at ICLR 2026
Table 1: Comparison of algorithms on 20 nodes network and river-runoff dataset
20 Nodes synthetic dataset
Real-worid river-runoff dataset
Algorithm
CS
TPR
FPR
Time
Fi
CS
TPR
FPR
Run Time
Fisher-Z.PCC
24.44
21.73
36.48
86.25
0.625
52.63
64.32
76.92
87.4
0.526
CD-NOD
23.80
10.93
32.25
87.68
0.421
54.05
65.16
76.92
88.24
0.1794
18.65
50.9
58.06
92.84
0,939
28.37
36.46
89.92
0.22
TPC
78.57
70.12
70.97
99.15
4016,57
43.64
67.1
92.31
74.79
1548,29
FFT-WPC
100
100
100
100
1011.62
75
67
09.2
98.
9.011
Wiener Phase
93.55
92.97
93.54
90.43
1.37
48.3
46.29
53.85
92.44
0.02693
6.2 COMPARISON WITH STATE OF THE ART METHOSS AND SCALABILITY
To further demonstrate the efficacy we present comparative results for network with 20 nodes. The details of the network under consideration is given in Appendix C.3 We compare the performance of our FD method based algorithms with state of the art methods, namely, Fisher-Z test based PC Kalisch & Bühlman (2007), CD-NOD Huang et al. (2020), Granger Causality (GC) Granger (169)). and Time Aware PC (TPC) Biiwas & Shlizerman (2022), on synthetic data. The Fisher-Z PC, and CD-NOD algorithms were obtained from causal-lear library [Zheng et al. (2024), and
the TPC and GC algorithms were obtained from the TimeAwarePC library Biswas & Shlizerman (2022). To compare the algorithms, we define the following metrics: True positive rate (TPR) FP (1- -FFFFN x 100%, Combined Score (CS)=IPR-FPP) -FR, and F-1 F-1 score (Fi) = ZTP-EP+FN X 100%, where TP, FP, TN, and FN stand for True Positive, False Positive, True Negative, and False Negative respectively. We also report the run times in seconds. A dataset containing 128000 samples was generated using the model of equation 7 for 20 nodes network with 22 edges. The accuracy metrics and the runtime of each algorithm using the dataset of 20 nodes is shown in Table The thresholds for the FFT-WPC and Wiener Phase algorithms were chosen to be 0.034 in both algorithms which were determined using the sparsity based threshold tuning algorithm in Appendix C.2 For the GC and Fisher-Z. PC significance levels of 0.106 and 0.01 were used respectively. The TPC algorithm was implemented
TP TP+FN X 100%. 1-Faes Positive Rate (FPR) 2TP
with maximum tima-delay of interaction be 5, and a significance level a = 0.1 for kernel-based conditional independence tests. For the bootstrap procedure in TPC, 50 bootstrap iterations with bootstrap window length of 50 recordings and bootstrap stability threshold y 0.01 were chosen. can be observed that FFT- WPC and Wiener Phase algorithms outperform TPC, Fisher-Z based PC, and GC terms of accuracy. Also observe that FFT based WPC beats TPC in terms of speed also. Moreover Wiener-phase beats both TPC and GC in speed of execution.
To further examine the scalability of the proposed algorithms, we applied the FFT based Wiener- PC and the Wiener phase algorithm to a network with 50 nodes and 54 edges. The details of the generative graph is given in Appendix [c.4 The FFT based Wiener PC algorithm with a threshold of 0.035 resulted in TPR = 07.62% FP) of 142652 seconds (= 40 hours). The Wiener Phase algorithm with a threshold of 0.051 resulted in TPR = 100%, (1 FPR) 99.83%, and CS = 99. 83%, with a total runtime of 19.49 seconds.
99. 99.. and CS 97.5%. with a total runtime
6.3 RESULTS ON REAL FIELD DATA In this section we demonstrate the performance of the proposed method on real- world benchmark dataset, the river-runoff dataset, obtained from online repository of CauseMe Runge et al. (2019).
The river runoff data has 4600 samples of average daily river runoff of the upper Danube basin Ji jtt al. (2024). The graphical interpretation of the river runoff data was taken from Ji et al. (2024)
as the ground truth. The dataset was used in the FFT-WPC, Wiener Phase, Fisher- Z test based PC.
CD-NOD, GC, and TPC algorithms. For the FFT-WPC and Wiener Phase the thresholds were chosen to be 0.404 and 0. 0.108. The significance level parameter for GC and Fisher-Z-PC were taken to be 0.112 and 0.01 respectively. The parameters for the TPC algorithm were taken from Appendix-D of Biswas & Shlizerman (2022). The comparative results are shown in Table Observe that, FFT-WPC was able to obtain comparable error performance with the best algorithm reported in Biswas & Shlizerman (2022). TPC, while run time was orders of magnitude faster. Wiener phase was the fastest with runtime of 26 milliseconds, but the CS was lower than FFT-WPC and TPC. However, Wiener phase was able to beat GC in CS and F-1 score.
6.4 APPLICATION ON CIRCUITS
In this section, we demonstrate the application of the proposed methods in circuits that employ transistors, capacitors, and resistors. In Rana et al. (2025) it has been shown that electronic circuits can be represented as LDIM. To demonstrate the effectiveness of the proposed algorithms in this
9

---

Published as a conference paper at ICLR 2026
Table 2: Comparison of algorithms on analog electronic circuit dataset
MOSFET Hardware Circuit Dataset
BJT Hardware Circuit Dataset
Algorithm
TPR
FPR
Time
CS
TPR
FPR
Run Time
Fisher-Z.PC
40
11.11
66.67
44.44
0.021
40
0
100
0
0,019
CD.NOD
66.67
66.67
100
66.67
0.0323
40
0
100
0.0309
100
88.89
0.0389
77.78
100
77.78
0.0391
TPC
77.78
100
77.78
20.816
60
55.56
100
55,56
23.641
FFT-WPC
100
100
100
100
0.08412
100
100
100
100
0.0798
Phase
100
100
0.0464
100
100
0.0419
article over state-of-the-art methods, we collected two datasets, one containing measured data from bipolar junction transistor (BJT) based physical hardware set-up and another containing measured data from metal-oxide-semieee seaide-semiconductor field-effect transistor (MOSFET) based hardware set-up.
First, we present results on the dataset collected from BJT based physical hardware circuit im-
plemented on printed circuit board (PCB). Fig. 2 shows the circuit schematic along with the
PCB and generative graph. To de- velop the circuit, MMBT2222ATTIG BJT units were used along with resis- tors and capacitors. The node voltages {vI,...4s were measured at IMS/s rate using a Ni-cRio data acquisition system and AD8421 instrumentation amplifiers. A detailed description of the hardware set-up is given in Ap- pendix C.9 A total 500000 sam- ples were collected for each node volt- age. The results of reconstruction are shown in Table 2 Observe that the proposed method outperforms state- of-the-art methods in terms of both accuracy and speed.
ADOAI Besd
Gravied Clivics Bounds
Arcasse Ramel Powe Secppie UT3D1 Doeed Powerr Sarplly
Cliratt Unake Tewt
Goople
Figure 2: Schematic of BJT based hardware implemented circuit showing PCB implementation, circuit under test, and generative graph
Next, we present results on the dataset collected from a MOSFET based hardware circuit implemented on a printed circuit board (PCB). Fig.
B shows the circuit schematic along
with the PCB and generative graph.
To develop the circuit, 2N7002E channel MOSFETS units from On- semi were used along with resistors and capacitors. The node voltages (v),...04) were measured 1MS/s rate using the set up described in the previous section. A total 500000 samples were collected for each node voltage. The results of reconstruction shown in Table [ Observe that similar to the BJT circuit results, the proposed method outperforms state- of-the-art methods in terms of both accuracy and speed. This demonstrates the dominance of our proposed method for the tested real-world dynamical systems.
ADMU2I Roud Chrrats
LIC1166 Busal
ETT006 Bood Power Supia
Chrrch Unier Tes
Gremneoo Graph
Figure 3: Schematic of MOSFET based hardware imple-
mented circuit showing PCB implementation, circuit under test, and generative graph
CONCLUSIONS
This article answers the question of whether time-domain (TD) or frequency-domain (FD) approach to reconstruct causal structure of dynamical networks is better. It is established through theory and simulation results that FD approach is better than the TD approach in terms of computational complexity. Further, simulation results show that the FD approach outperforms TD. In addition, the superiority of the WF-based causal inference over traditional static CI test-such as Fischer-Z test, and state of the art methods such Time aware PC, and Granger Causality is demonstrated through simulation. As a third contribution, an algorithm for causal structure recovery using magnitude and phase properties of the WFs is developed. This algorithm provides significant computational advantages over traditional causal inference algorithms by ruling out the combinatorial explosion of required CI tests.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work was supported by National Science Foundation via grant number ECCS 2412435 and by United States Department of Energy via grant number DE-CR0000040.
REFERENCES
Stefanos Bennett and Rose Yu. Rethinking neural relational inference for granger causal discovery. In NeurlPS 2022 Workshop on Causality for Real-world Impact, 2020.
Michel Besserve, Naji Shajarisales, Dominik Janzing, and Bernhard Schölkopf. Cause-effect in- ference through spectral independence in linear dynamical systems: theoretical foundations. In Conference on Causal Learning and Reasoning, pp. 10-1143. PMLR, 2022.
Rahul Biswas and Eli Shlizerman. Statistical perspective on functional and causal neural connec- tomics: The time-aware pe algorithm. PLOS Computational Biology, 18(1 10166605 2022.
David Carfi and Giovanni Caristi. Financial dynamical systems. Differential Systems, 2008
John AWB Costanzo and Osman Yajann. Data-driven i/o structure learning with contemporaneous causality. IEEE Transactions Control of Network Systems, 1(4):1929-1939, 2020.
David A. Cox, John Little, and Donal O' Shea. Ideals, Varieties, and Algorithms-A An Introduction to Computational Algebraic Geometry and Commutative Algebra. Springer, 2007.
Harish Doddi, Deepjyoti Deka, Saurav Talukdar, and Murti Salapaka. Efficient and passive learning of networked dynamical systems driven by non-white exogenous inputs. In International Conference on Artificial Intelligence and Statistics, pp. 9982-9997. PMLR, 2022.
Zoubin Ghahramani. Learning dynamic bayesian networks. Adaptive Processing of Sequences and Data Structures: International Summer School on Neural Networks "ER Caianiello" Vietri sul Mare, Salerno, Italy September 6-13, 1997 Tutorial Lectures, pp. 168-197, 2006.
Asish Ghoshal and Jean Honorio. Learning linear structural equation models in polynomial time and sample complexity. In Amos Storkey and Fernando Perez-Cruz (eds.), Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pp. 1466- 1475. PMLR, 09-11 Apr 2018. URL mlr press v// ghoshallhsaalll//a/ 11a html
Loukas Grafakos. Classical fourier analysis, volume 2. Springer, 2008.
Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral methods.
Econometrica: journal of the Econometric Society, pp. 424-438, 1969.
Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, USA, 2nd edition, 2012. ISBN 0521548233.
Pengfei Hu, Rong Jiao, Li Jin, and Momiao Xiong. Application of causal inference to genomic analysis: advances in methodology. Frontiers in Geneties. 9:238, 2018
Biwei Huang, Kun Zhang, Jiji Zhang, Joseph Ramsey, Ruben Samchez-Romero, Clark Glymour, and Bernhard Schölkopf. Causal discovery from heterogeneoussnonstattonaay data. Journal of Machine Learning Research, 21(89):1-53, 2020.
Junzhong Ji, Ting Wang. Jinduo Liu, Muhua Wang, and Wei Tang. River runoff causal discovery with deep reinforcement learning. Applied Intelligence, 54(4):3547-3565, 2024.
Thomas Kailath, Ali H Sayed, and Babak Hassibi. Linear estimation. Prentice Hall, 2000.
Markus Kalisch and Peter Bühlman. Estimating figh-dimensional directed acyclic graphs with the pe-algorithm. Journal of Machine Learning Research, 8(3), 2007.
11

---

Published as a conference paper at ICLR 2026
Joshin Krishnan, Rohan Money, Baltasar Beferull-Lozano and Elvin Isufi. Simplicial vector autore- gressive model for streaming edge flows. In ICASSP 2023-2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 1-5. IEEE, 2023.
Andrew Lamperski. Non-asympiotic pointwise and worst-case bounds for classical spectrum estima- tors. arXiv preprint arXiv:2303. 1908, 2023.
Gabriele Lohmann, Kerstin Erfurth, Karsten Müller, and Robert Turner. Critical comments on dynamic causal modelling. Neuroimage, 59(3):2322-2333 2012.
Ziyu Lu, Anika Tabassum, Shruti R. Kulkarn, Lu Mi, J. Nathan Kutz, Eric Todd SheaBrown, and Seung- Hwan Lim. Attention for causal relationship discovery from biological neural dy- namics. In Causal Representation Learning Workshop at NeurlPS 2023, 2023. URL https //openreview net (twprrrrieeeeeeteoooaaaaaaieoaoe ttte://foorrmmmiixxxxxx
Donatello Materassi and Murti V. Salapaka. On the problem of reconstructing an unknown topology via locality properties of the wiener filter. IEEE Transactions on Automatic Control, 57(7))1765 1777, July 2012. ISSN 0018-9286. doi: 10. 10.1009/TAC.2011222131777
Donatello Materassi and Murti V Salapaka. Reconstruction of directed acyclic networks of dynamical systems. In 2013 American Control Conference, pp. 4687-4692. IEEE, 2013.
Donatello Materassi and Murti V Salapaka. Graphoid-based methodologies in modeling, analysis, identification and control of networks of dynamic systems. In 2016 American Control Conference (ACC), pp. 4661-4675. IEEE, 2016.
Donatello Materassi and Murti V. Salapaka. Signal selection for estimation and identification in networks of dynamic systems: A graphical model approach. IEEE Transactions Automatic Control, pp. 1-1, 2019. ISSN 1558-2523. doi: 10. 109//7AC.2019922900000
Raha Moraffah, Paras Sheth, Mansooreh Karami, Anchit Bhattacharya, Qianru Wang, Anique Tahir, Adrienne Raglin, and Huan Liu. Causal inference for time series analysis: Problems, methods and evaluation. Knowledge and Information Systems, 63:3041-3085, 2021.
Alan Oppenheim and George C Verghese. Signals, systems & inference. Pearson London, 2017.
Judea Pearl, Madelyn Glymour, and Nicholas Jewell. Causal inference in statistics: A primer, John Wiley & Sons, 2016.
Adrián Pérez-Suay and Gustau Camps-Vall Valls. Causal inference in geoscience and remote sensing from observational data. IEEE Transactions on Geoscience and Remote Sensing. 57(3):1502-155) 2018.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Causal inference on time series using restricted structural equation models. Advances in neural information processing systems, 26, 2013.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
Juan F RRaairez-Villeeeeee Villegas, Michel Besserve, Yusuke Murayama, Henry C Evrard, Axel Oeltermann, and Nikos K Logothetis. Coupling of hippocampal theta and ripples with pontogeniculloocciptaal waves. Nature, 589(78400:96-100 102, 2021.
Mohammed Tuhin Rana, Mishfad Shaikh Veedu, and Murti V Salapaka. Uncovering the influence flow model of transistor amplifiers, its reconstruction and application. arXiv preprint arXiv:2508.040 .04977, 2025.
Jakob Runge, Sebastian Bathiany, Erik Bollt, Gustau Camps-Valls, Dim Coumou, Ethan Deyle, Clark Glymour, Marlene Kretschmer, Miguel D Mahecha, Jordi Muñoz-Marí, et al. Inferring causation from time series in earth system sciences. Nature communications, 10(1) 2:253, 2019.
Naji Shajarisales, Dominik Janzing, Bernhard Schölkopf, and Michel Besserve. Telling cause from effect deterministic linear dynamical systems. In International Conference on Machine Learning, pp. 285-294. PMLR, 2015.
12

---

Published as a conference paper at ICLR 2026
Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non-gaussin non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
Milad Siami, Sadegh Bolouki, Bassam Bamieh, and Nader Motee. Centrality measures in linear consensus networks with structured network uncertainties. IEEE Transactions on Control of Network Systems, 5(3):924-95) 2017.
Peter Spirtes, Clark N Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2000.
Saurav Talukdar, Deepjyoti Deka, Harish Doddi, Donatello Materassi, Michael Chertkov, and Murti V. Salapaka. Physics informed topology learning in networks of linear dynamical systems. Automatica, 112:108705, 2020. ISSN 10005-1095 doi: http:/o..orr.110010.o.atttttttto 108705.
Mishfad Shaikh Veedu, Harish Doddi, and Murti V. Salapaka. Topology learning of linear dynamical systems with latent nodes using matrix decomposition. IEEE Transactions on Automatic Control, Early Access:1-1, 2021. doi: 10. (0.1109/TAC.2221.11229999
Mishfad Shaikh Veedu, Deepjyoti Deka, and Murti Salapaka. Information theoretically optimal sample complexity of learning dynamical directed acyclic graphs. International Conference on Artificial Intelligence and Statistics, pp.
2024.
Andrew J. Wren, Pasquale Minervini, Luca Franceschi, and Valentina Zantedeschi. Learning discrete directed acyclic graphs via backpropagation. In NeurlPS 2022 Workshop on Causality for Real-world Impact, 2022. URL U http://peeneeeeeeeeteroorommoooo Sttpr:uuiiiisSSI1-----
wwt.
Yujia Zheng, Biwei Huang, Wei Chen, Joseph Ramsey, Mingming Gong, Ruichu Cai, Shohei Shimizu, Peter Spirtes, and Kun Zhang. Causal-learn: Causal discovery in python. Journal of Machine Learning Research, 25(60)11-8, 2024.
13

---

Published as a conference paper at ICLR 2026
APPENDIX
A
WHICH FREQUENCY TO CHOOSE?
Most practical systems are finite dimensional with finite dimensional realizations. These admit rational transfer functions. The following Lemma proves that, in identifying the CI, graph/moral graph, and skeleton, it is sufficient to work with one frequency selected arbitrarily, in order to get meaningful result, when we have access to large enough data.
Lemma A.I. Consider any rational polynomial transfer fumction
(2) C (= a2) ... (a ap)
(G-bb)(e b2) (z - by) p psq. PSq.CE CeC.
For any distinct w1, u2 E Ω Hey(e * if and only if (elwe) # almost surely w.t. continuous probability measure.
Proof:
By the fundamental theorem of algebra Cox al. (007). there exists most p complex zeros for (z), if Hz)(z) is not identically zero, which form set of Lebesgue measure zero Suppose (e)() + O, which implies that H:j is not identically zero. Thus, (e)() j 0 almost everywhere.
B PROOF OF THE RESULTS IN SECTION 5
| W.cll = 1.0 c0 2 folle
- ( 2c)ll2 + |0 2.0) 2
< @ 3cll c d.cllel
@ 3cll + M a. +.c122
S (103c +icll + [[
कco + M 0I - s2cl22
c. ccll2 कh 1 ccl [bcll &dl
+ M1 - Dicll2
3 cc $cclll| & Ml
< ff + Meine I (M + 1)e.
where we have used triangle inequality, properties of spectral norm, || and Euclidean norm || 112-
From Section 5.8 in Horn Johnson 2012).
M
|| ccoll < | I
coll
4 ct T5ll
1 1 M2
1400
M |-ق - $ 352
<E = | - c S
Me S
1 MI 4 - 1 clll
M
Therefore, guarantee that 100 2 < is sufficient to guarantee that | &ooll < since 2 I. Hence,
P W.oll > e) < (1000 Tooll
(8)
M M+1
For Blackman- Bhackman-Tuckey.
Lemma B.1 (Lamperski 2023 Consider a linear dynamical system governed by equation 2
Suppose that the auto-correlation function Ra(k) == E{X(n)X (n k)} satisfies exponential
decay, [k)(l) S Cp-lk and that there exists M such that S Amin (+x) )= S Amaz (Φx) S M.
Then for any < e and L Z -logp
(1-ple)
20
e(T L)
CC 00
(TT
#ool >() S exp
+ 6n 33(22++ +
14

---

Published as a conference paper at ICLR 2026
Table 3: Coefficients used in data generation for random networks experiment of Section 6.1
e: (R)
1.73938 26 13946666889 x 10-
1.0362877669668 550 x 10
5:21814377847440007777 x 10
(k)k)
1.41520817200220955 10
1.9600S12688000855525 x 10
6.93041163440004197715 10"
(k)
77891029271012888199 10
5.902372011897300229. x 10
3.3156411708
(k)
8.1157990006112205000 x 10
1:13229952818880255 x 10 2
2.653973971707788600 x
(k)
1.302817022311100443388
5.1376905777 10861684 x
2.58884528881--338422
kolk)
L6W229666733666211994
10
L.141275502255119577 10
4.8808906026668845551
Applying this,
IP - >6) SP SP (14000 <P(1400
-1400- -
- (TT.1)
M +1
-T)
S exp
32(2T. + 1)M2 +12 + 6n
(T L))
I exp
32(2L + 1)(M 1)4 6n
Thus, an upper bound on the sample complexity is
(T - L)
exp
32 + 1)(M 1)4 + 6n =
(T L)
32(21. 1)(M 1) (log(0) 6n)
(T - L) = (log(0) 6n) ((2L + 1)M*)
22
For Bartlett estimator,
Lemma B.2 (Veedu et al. (2024).. Consider a linear dynamical system governed by equation 2
Suppose that the auto-correlation function Rg(k) (B(((n)))) (( (n++k [h++k) k)} satisfies exponential
decay, (k)]l) C8-Ik) and that there exists M such that Amin (xx) S Amaz (xx) M.
Let O E1,E2 be such that E2 - Suppose that N ww Q
11-0-19261 where E1. Then
2Cp
2T
P pl0Cc - &ooll > S cxpp
( 128N + 6n
Then.
P
SP | - &cel
M +1
2²
s exp
+ +6n
128N + 12
2T
$ exp
6n
(10)
12 + +14
C EXPERIMENTS
In this section, we provide the details on data generation, threshold tuning, and simulation results on statistical significance.
C.I DATA GENERATION
For MEG reconstruction of random networks, we performed the experiments over random networks with n. 6 nodes. For each network iteration, random DAG is first generated with most
15

---

Published as a conference paper at ICLR 2026
2 parents for each node. In this graph, the data is generated using the AR model, with the non-zero entries that respect the graph structure of the random DAG generated, using
(1) (1) 1) (2) 2) + (3),. 3)
Σ j
bajaz - +e;))
(1I)
where (t) is zero mean i..d. Gaussian noise with diagonal covariance matrix and 1,........... ..... ................................ n. The
data is generated continuously according to the AR model in equation I1 That is, for f = 0 to T,
(t) computed based on equation 11 which satisfies Assumptions and Hence Lemma 4.2 and
Algorithm 2 can be applied to reconstruct the essential graph.
The coefficients a-(-) used for the data generation are shown in Table It is be noted that that for the random networks the bij coefficients were generated randomly from uni/ [0.2, 0.4). The coefficients are chosen this way to stop the time-series from becoming unbounded. IID exogenous noise setting, (t) N (0, 1) is used. Given the coefficients, the data is generated according to equation 11 for every E
C.2 RANDOM NETWORK RECONSTRUCTION FROM DATA
The datasets for the 25 random DAGs used in FD based W-PC, TD based W-PC, and Wiener- phase algorithm. The average error rate and computation plots are shown in Fig. [(a) and [(b] In this section, we quantify the performance of the algorithms statistically. The box and whisker plot of eтгоr rates for the 25 random networks are shown in Fig 4 As expected the error reduces as the number samples increases.
As mentioned in Section 6 the performance of the proposed algorithms in this article depend on the choice of the thresholds. To tune the thresholds we use sparsity based tuning approach. This approach uses bounds on the average, maximum, and minimum degree of the reconstructed topology choose a threshold. The algorithm is shown in Algorithm Such sparsity related information are typically be available in practice for DAG models.
Algorithm 3 Threshold tuning- based on sparsity metric
Input: Data, upper limit of average degree du, lower limit of average degree d,, maximum degree dmaz. minimum degree dmmnn lower limit threshold 7;, upper limit of thresholdr,, resolution for tuning Ç.
Output: Optimal choice threshold
1. Initialize threshold f + 7,
7
2. p + floor [(Tu 7)/5]
3. For i I,...,P
(a) Estimate & using reconstruction algorithm and
(b) Compute average degree (d), maximum degree (dmaz), minimum degree (dmin) for skeleton of G.
(c) d d < du and dmaz < dmax and dmin > dmin
7 T
break.
(d) - (T +
4. Return
C.3 DETAILS OF 20 NODES NETWORK
The 20 nodes network used in this article is directed acyclic graph shown in Fig. 5 The network consists of 20 vertices, and 22 edges, with maximum and minimum degree of the nodes being 5 and 1. The maximum number of incoming edges at any specific node 2 and maximum number of outgoing edges from any specific node is 3. As discussed earlier the data generative model for the 20 nodes network was taken to be of the form in equation The a.(.) coefficients for the 20 nodes network were chosen shown in Table
16

---

Published as a conference paper at ICLR 2026
WAPC
Meserghase
2.30
3000
6:00
23000
128000
60
1000
6300
23000
300000
623
1206
0900
20403
100000
(a)
(b)
(c)
Figure 4: Box and whisker plots showing statistical performance of the (a) FD based W-PC (b) TD based W-PC (c) Wiener-Phase over 25 networks for the results presented in Section
3
19
17
Figure 5: Generative graph of 20 nodes network
C.4 DETAILS OF 50 NODES NETWORK
The 50 nodes network used in this article is directed acyclic graph shown in Fig. 6 The network consists of 50 vertices, and 54 edges, with maximum and minimum degree of the nodes being and 1. The maximum number of incoming edges any specific node is 2 and maximum number of outgoing edges from any specific node 3. The (,).) coefficients for the 50 nodes network were chosen shown in T
17

---

Published as a conference paper at ICLR 2026
Table Coefficients used in data generation for 20 nodes network
(1(A)
1.73938 16 1596696899 02
1.0436285665681 18150e
03
(2(
3. 211220N12002USRS60
01
1. 3660052688000335500
01
13301166110011991111
02
(k(N)
8.289 02927 042288 1996
0.80237204897300 295e
3.315641 7084
(g(k)
8.846579905691 620560c
02
1.4252899524458102222
1.65997397770748000
(k(
8.56281 76295 1122436e
01
5. 1117690677106606666
01
2.56884628885543088
01
(g(A)
1.630296867365622 1944
01
11412078072280855997 01
-8908900022988888811
02
(13998126159100 161591696899
10156287669588111000
6.2181137817710001
(k)
3.415208 172002098808e
1.88003368800883077
-1
03301134400441000
02
ag(R)
8. 11289102922 10122888888888888 1996 01
0.802372068973000 01
3.315641 7160122
01
(k)(k)
2.65397397
(A)
5. 1137690577718666668
1.568845288855433388
(k)
16332298675650222 1916 01
111120780729593909777 01
189089060/29655111111
02
G1a
1.739381200 129 16968399
02
L.0466287669561811
02
(1)
3.415208
172002098808
13860033668008330222
(k)
8.2891 029271 04288 199c
01
1.8023720489733012255
01
3.315641 17844171 5502e
01
(1)
8.846579905699 620360e
02
.4252819525458112880 02
2.65397397 11S3973977701 17777777748000 48769699
02
(1)
762951810224666
-117600577106666644
2.568845288855
(16
1.630296867
077799666766666218 5656281916
01
1111307807225930555
8898896612966888888
02
(19
1.73938 (26 26139 10968599
02
L0436287669668 18 150e 02
6.21814378477
03
0go
72002098808
1.8600520850003055333
6.880416634400419776655
10
DB
Figure 6: Generative graph 50 nodes network
C.5 ROBUSTNESS OF WIENER PHASE ALGORITHM TO PHASE CONDITION VIOLATION
In addition to the results in the main article. our algorithm remains robust to assumption violation as demonstrated through example here. To study the robustness the Wiener Phase algorithm we employed the Wiener phase algorithm on synthetic data that violates the phase alignment assumption. The model was taken be of the following form: 0;,02;(()) + (t, 1) + 0;,2=;;( 2) 1) b, - 2) + bt - 3)) + (t()), where the coefficients bij,A were chosen randomly from [0.1, 0.9], and a,,A were chosen such that the time-series remain bounded. The network was taken to be a 20 nodes network shown in Appendix (C.3) For the 20-node network with 22 edges, and N=32 point FFT, and threshold of 0.15 5 the results shown in Table
Ej
18

---

Published as a conference paper at ICLR 2026
Table 5: Coefficients used in data generation for 50 nodes network
k
A = 2
(x(k)
1.73938126 159 1696859
20
L.1166287695588 18 150c
00
5.2181 10907511
03
(2(k)
3.415208 172002098808
1336003268800
6.8304 000044633440041 1976
15
02
o3(k)
8,2891029271
104288 1996
10
298023720189730000
10
3.3156111 138166111700111 708 WOC111116600 16502e
01
(4(h)
8.84657990569 620500e
02
4.4232800528455 10280e
02
2.65397397 1707486099
02
(k(k)
166281162951 181002466
1.156605571106666566
2.568845288855
(x(A)
01
1. 1411007807 295839997
01
4. 8890906120965588811
02
1.739381 119998126159166866
02
1.0066277569661111888 02
5.21814378447 10907511
03
(2(k)
3.415208 172022098808
01
1.36606326880033006
01
6.8304113 10001634000-1971116
02
(g(k)
8.289 02927) 04288 199c
5.8023720489730012229
3.315641 178441 71 5502e
(k)
(4465799056916205000
02
4. 1.42289965245581128800 02
166519799777774866999
02
(11(k) (R)
11228176295188 11024366
01
1157600577706666644 01
01888452888554338888
01
1.630296667568 7668666899666
1412078072958395977
L89089061269688
(k)
1.73938 126 594066839
02
0A462876695508 18100
02
122181377847840000
03
(k)
3.4152081 1 120020088088
01
138606326680083305226
01
6.83041 63440041499000 197615e
02
(1)
8.259 (2927 T04288 199c
1.80237204897300 1295e
3.315641 778447 5002c
(k)
8.846657990569 620560c
02
1.4252599526455102288
02
2.653993971700 170745699
02
017 (k)
1.562817629518102
01
5. 6.13799077 77769057771061166 71086168 le
01
11.668111288555 15288855-1388122
01
(k)
L.6302006 1946
01
141220807 2960399976
01
1.8008000120988880888
02
(1)
1.73938 16966699
10136887369508
1.218113781781110
(k)
3.415208 122022088088
01
1386063208800830022
01
6.83041 4344004199761
02
0ag (1)
1.73938122 126 1394666399
02
1.046627766566181155
02
2318144384784000733
03
3.1152U81 3.11520812202299906 20002098808
02
(lk
8,289 (02027 I 104288 1996
01
6.80237204897300 01
3.315641 1708417 706417
5002e
01
0:24 (E)
1.816579900665 6205600
02
1.1232595521158
02
2.65397397 777733977777770 1707/89999
02
0go
8.56281 76295 0102436
5. 5.13769057771
(Ac
1.63029686756566 1944
01
L.14120750729999359777 02
02
0gz (1))
1.73938 11593812615999 1696859e
02
1.016628766965111119
5.21811378177 1090751
03
028 (k)
3.415208177 3.4152081720020055
20022998808
01
1.8608326880053355222
01
02
(1)
(0211002221012000 199c
0J0033720299330022996
3.3156-11 170841710002
0:g (k)
8.84657900669 6205500
02
4.42328005284558 102220 02
2.65397397 01728977777479999
02
4gi (k)
8.562817622 7629518102488 4302436e
01
01
2.5688452888555 43845
01
0.sa (J:)
1.630296867566628 194c
-141207007295889597
1.89089061126968855581
02
(A:
1.739381 2818811619999 16968599
02
102362873695681811000 02
5.2151
03
0.sa (k)
3.415208
172002098808
01
1. 1.88005326880088000
01
G.8304166 133011634000-11999 900-4197616
02
0gg (A:)
8.2891 (29921 04288 199c
5.802372048973001222
3.315641 70841 5502e
(3)
9.844657990669 620560e
02
4.42326095284588 10280e
02
2.65397397 1707496099
02
(k)
1.562817629511 18 8802443e
01
5. 1.3769057710661668
01
2.5688452888554 430842e
01
U.a (A)
1.63029686 7000867766550 7565628169 1911
01
L141207807295038555
01
02
0:so
1.73938 126 15946663339
L04332877695608 18 8150e
6.21811378448-4990777
03
3.415208 10100081720000 22022098808
01
1 366083268800833662
01
02
(41(k) (k)
8.5628 76295 18 102436
01
5. L.1766065771006 1684e
01
2.16884412885553308.22
01
(Ap)
16630296667666226
.141207807295939599
1.890889661266988853311
02
(Ac
1.739381 1215566116668
02
1,0366283669691811100
02
5.2151
03
0444 (I))
1115220112202099600
01
1.98006326880083552222
01
6. 68804166440044197615
02
Osea (A)
8.289102922 211228111228111 1996
01
5.8023720199730022222
01
3.315611 11166117701111160000
01
046
8.84657990569 6205660
02
4.4232699528455 102800
10280e
02
2.65397397 1777480999
02
(k:)
8.5628176288 772777611292900000
01
5. 111769957770866666
01
2.56884528885543088222
01
0:s (k)
1.63026866765588
194e 01
1.14120880299500055
01
L.89089060029988558811
02
1.799981261590 16559199655
104562875695681111555
5.21814974474400075551
(A)
21115208120002988000
01
336606326880033302226
01
02
Table 6: Performance of Wiener Phase Algorithm on Phase Violation
Sample Count
05%
TPRC%
FPRY%
6-100
86.52
87.09
99.43
0.117
32000
88 89
90.32
98.57
0 0.2298
125000
92.69
93.55
99.14
1.279
19

---

Published as a conference paper at ICLR 2026
It can be observed from the above results that even though the phase alignment assumption is violated the Wiener phase algorithm produces highly accurate results. This shows that the algorithm robust to the assumption being violated.
C.6 APPLICATIONS SATISFYING ASSUMPTION I
Any network whose states evolve according to the following dynamics,
d' + d0
Σ Cij (zy - 14) pc()),
Omn d
satisfies Assumption | Some examples include Talukdar et al. (2020)
Consensus dynamics: In multi-agent systems, distributed decision making often follows the following first order consensus dynamics, where each agent i updates its states according
dr; dt
Σ Cjlly - 2;) + pc(t), j=1
where Pj denotes the receiver noise for agent j Siami al. (017)) Identifying the commu- nication topology of network of agents in multi- agent system is relevant objective of cyber attacker, and appropriate hardware/software tools need be designed to avert such attacks.
2. Thermal Resistance Capacitance (RC) networks: One of the popular approaches in modeling the real-time control of building is the gray box modeling with lumped parameters. Here, the discretized physical spaces are modeled as RC models, where each discretized physical space is represented by common temperature. An example lumped parameter model is
di di Σ
+ (1), -ー
T-T Ríj
C
where T; denotes the temperature, C. > 0 denotes the capacitance and p, denotes the total internal heat generated of zone Rij 0 is the thermal resistance between zones and j-
3. Power Grid Dynamics: When the disturbances in the power grid are small, the deviation in voltage phase angle dynamics at bus j. denoted by 8y. can be modeled by a swing equation, for example,
drr;
Di dr;
M,
+
Σ
cıyl2) xzi) + Pj.
dt2
dt
j=1.j
where Pj denotes the receiver noise for agent j.
4. Linearized Chemical Reaction networks: Chemical engineering processes often operate in transient states and are generally governed by ordinary differential equations (ODEs). ODEs are used in applications such as chemical stirring tank design, heat exchangers, or biological cell growth. Example, see here
C.7 EXAMPLE WHERE GRANGER CAUSALITY FAILS
An essential assumption for applying Granger causality is that the underlying influence dynamics should be temporally strict causal. That is, if x influences y, then "only" the "past" values of x in time should influence y. Granger causality will fail the present value of y depends on the present value of well. However, our approach present value of y can depend on present, past, or future values of r. We generated a 6-node synthetic network with zero lag dependency between the nodes that follow the dynamics described as below
accommodate any linear dynamical system where the
3
(+) + a. (1) - 1) + a, (2)x, - 2) a, - 3) = Σ>
(t + et)). (12)
7=0 3
20

---

Published as a conference paper at ICLR 2026
Table Algorithm performance on non-linear system
Algorithm
TPR%
FFR)%
C5%
HPT-WPC
100
100
100
Wiener Phase
100
100
100
The following results show that the GC fails to reconstruct the full network even at 10 samples. Our algorithms achieve zero reconstruction error with fewer samples.
Wiener Phase: perfect reconstruction at 6400 samples with a N = 32 point FFT and a threshold of 0.1.
FFT-WPC: perfect reconstruction at 15000 samples with a N 32 point FFT and threshold 0.074.
Granger Causality: CS 80% 3.2 10t samples, with significance level a 0.086.
C.8 APPLICATION OF ALGORITHMS TO NON-LINEAR SYSTEMS
Note that the proposed algorithms are purely data driven and thus are applicable to any set of time-series data; the underlying data generative system does not determine the applicability of our algorithms. This is in contrast to other methods, for example, where the algorithm depends parametrized model, and the parameters are determined via data. To demonstrate that our algorithms may be effective in many non-linear system present example below.
We used a generative model the form
(i.025((() + az, (t 1) a,,20; 2)) + a1,3zi(1 - Σ (1)((,(()) 1)) ex()):
(13)
where (t) is standard Gaussian noise. byg uniform random numbers, and a; chosen such that the time series remains bounded (poles of the transfer functions remain within the unit circle). The above model
used to generate data for the 20 nodes network of Appendix [C.3] The generated
data was then used in the FFT based Wiener PC and the Wiener Phase algorithm with thresholds 0.04 and N 32 point FFT in both algorithms. The results are shown in Table 7
It can be observed that the algorithm produce almost perfect reconstruction results for the nonlinear system described above.
C.9 DESCRIPTION OF SETUP FOR EXPERIMENTS ON CIRCUITS
As described in Rana et al. (2025), the experimental setup for data acquisition from transistor based circuits that was described Section 6.4 consists of the following hardware sub components:
computer station with LabVIEW 2021,
NI-eRIO 9054 equipped with Xilinx Artix-7 A100T FPGA,
an NI-9223 module with 4 + +10V isolated analog input channels, 16-bit simultaneous sampling of all channels maximum 1MS/s rate,
600W Multi-Range 150V/10A bench-top DC Power Supply (Model 9206) from BK Precision,
12V 50W AC/DC Converter (LRS-50-12) from Mean Well USA Inc.
low noise instrumentation amplifier (AD8421) based custom made measurement circuits,
hardware under test (HUT) which was either BJT or MOSFET amplifier circuit. The PCB implementations of the hardware under tests are shown in Fig and Fig.
The circuit under test energized using the bench-top DC supply from BK Precision and the NI-cRIO DAQ system was power from the 12V AC/DC converter from Mean Well USA. The node voltages, (UI, U2, V3, 14). shown in Fig.
and Fig. were then first amplified using the AD8421
21

---

Published as a conference paper at ICLR 2026
Table 8: Comparison of algorithms on simulated analog electronic circuit dataset
Algorithm
cs
TPR
FPR
Run Time
FC
37.81
37.74
70
67.71
0.1897
CD-NOD
37.84
37.74
70
67.74
0.298
GC
91.94
100
91.94
0.0928
TPC
79.03
100
79.03
97.81
FFTT-WCC
100
100
100
100
10.97
Wicner Phase
90.91
96.8
100
96.8
0.153
instrumentation amplifier. The amplified signals were recorded using the NI-cRio and NI-9223 based data acquisition system. The LabVIEW software was used to deploy the data acquisition algorithms required for the NI-cRio. A screen capture of a segment of the recorded voltage measurements are shown in Fig.
GUINSTEK
18k pts
28HSa/5
Stop
1,00 M
/4 A
kpppo A
M Ma
588nU
58BmU
28BmU
588mU
58us
B. 888888
Figure 7: Waveforms of recorded voltage measurements for the experiments on physical hardware circuits displayed oscilloscope.
C.10 RESULTS ON LARGER CIRCUIT SIMULATED IN CADENCE
In this section, we demonstrate the application of the proposed methods in larger circuits that employ transistors, capacitors, and resistors. To demonstrate the effectiveness of the proposed algorithms in this article over state-of-the-re art methods on larger circuits, we collected dataset containing high fidelity simulation data. A schematic of the circuit simulated in Cadence virtuoso, which is a high fidelity circuit simulator, is shown in Fig. g spice model of a commercial bipolar junction transistor (BJTT) (MMBT2222ATT1G) was used for the simulation. The circuit contains multiple transistor amplifier stages connected in common emitter configuration with bypass capacitors at the emitter. Also, the different stages are connected through capacitors leading to dynamic nature of the circuit. A total of 799999 samples were collected for each of the node voltages {01, shown in Fig. X The comparative results using the data collected from the simulated dynamic circuit is shown in Table methods. Note that
88 It can be observed that our proposed approach outperforms state-of- the FFT-WPC and the Wiener-Phase produce almost perfect reconstruction, and the Wiener-Phase algorithm is the fastest. The FFT-WPC achieves the highest Fl score followed by Wiener-Phase, GC, TPC, CD-NOD, and F PC in descending order.
1.11 EFFECT OF FREQUENCY DISCRETIZATION ON ALGORITHM PERFORMANCE AND CHOICE OF N
In this section, we present results showing the effects of the frequency discretization and sampling on the proposed algorithm. As described earlier, compute the FDD based Wiener filters sample
22

---

Published as a conference paper at ICLR 2026
RR
R4
R.
R,
Ra
R
SRz
11
R.
11.
RR.
31.
30
R2
R6
RR
50
Figure 8: Schematic and generative graph of simulated transistor amplifier
Rana et al. 2025)
the frequency interval Ω = [0, 2xr at N evenly spaced points. To quantify the effects of frequency discretization and sampling, we evaluated the performance metrics the algorithms at different values of N. The data generative model of equation IT used for generating data for 25 random DAGs each containing 10 nodes and randomly generated edges. The datasets, each containing 128000 samples, were used in Algorithms| [ and 2
with the sparsity based threshold tuning Algorithm 3 for different values of N. The average error rates and Fj scores over the 25 random DAGs are shown in Fig. It can be observed that the algorithms perform better the value of N increases. However, note that for N > 32 the improvement in algorithm performance with increasing N is incremental. Although increasing the value of N may increase the performance beyond 32 for the presented example, it would also increase the computational burden as the complexity of the algorithms is Ollog N). Therefore, to maintain an acceptable tradeoff between complexity and performance we choose N = 32.
Average Error Rate Comparison, n
-10
Average Score Comperison, n
1.0
FFT-WPC
0.9
Wiener Plase
0.8
0.8
Fi Soore
0.7
0.6
Errot
0.4
0.5
FFT-WPC
0.2
0.4
Wiener Phase
10'
10a
10a
10a
N
N
(a)
(b)
Figure 9: (a) Average error rate and (b) F, score for different values of N for FD Wiener PC and Wiener-Phass over 25 random networks
23

---

Published as a conference paper at ICLR 2026
D FREQUENCY DISCRETIZATION
The LDIM in equation 2 can be sampled at frequencies wt in fly = {wo, .u0.,..,,,y ,..,,,, where wi == 2mk thus converting the continuous frequency model to a practical model involving a finite set of N
frequencies. In TD this is equivalent to the finite impulse response model
N-1 Σ h(1)x(n -
x(n)
(14)
where n = O, O,...,N 1, x(n) and e(n) are periodic with period N, and h(() is non- zero for at most N lags. For the circular convolution model equation [4 ((wx) == ((n)e-2skn/N, 2xkn/N
VN EN= x(n)e
O.....
provides the Discrete Fourier Transform (DFT)
of the sequence x(n). Similar DFT expressions hold for H(w)) and E ((wk) as well, assuming h(n) and e(n) are periodic. The directed graph associated with equation 14 G = (V, E). is defined such that the edge - exists in E if (w))
for some k. It is well known that the Fourier coeffi- cients can be approximated via discretization Grafakos (2008). A uniform and explicit convergence rate for signals with bounded variation is provided below in Section D.I that shows convergence
rate of 1/ S/VN N for approximating equation with equation |4
D.1 CONVERGENCE OF DFT TO DTFT
Notation: For continuous function f. denotes the Fourier transform of f.
Proof:
Definition 1. Grafakos (2008) The total variation of continuous complex-vale valued function f. defined 79
on (a, b) C IR is the quantity, () = sup Σ | 1 (x;)), where P {rr: Fn -
pEP
partition and P is the set of all partitions of [a, 6].
A continuous complex-valued function f is said to be of bounded variation on chosen interval [a, C IR if its total variation finite, i.e. if (f) < +oc.
Lemma D.1. For a function f of bounded variation V, on 2x] and N 2 1, the estimation of f(m) = 1ea ((w))ewwad du given by (n) 1 N-1 @ibnnk/N f(2xk/N) for |n| S vN and zero otherwise, satisfies [|ʃ 1( S SC/VN.
N k
Proposition 1. For ff with bounded variation, VoL (f), on (a, bJ and t e [)
1. ((b)de-(6-a) s (b- (6
()-a)f(t)
Proof:
Since if has bounded variation, M == supzela,b) f(x) and m == inf f(x) are both
finite. Thus,
1 fxx)dr -bb (b (())dr-(- a)((t( s [0) ((t)) dr (M m)(b a) s (b- af.
{
Corollary 1. For J a function of bounded variation on [0, 1). given by Vf) the error in Riemann sum approximation,
N = 1.' M(r)d
N-1 f(a/N)
S V//N
N n=0
Proof: Note that
EN
<
Σ
++1/N }x)ddr
N
Fm/NN)
1
<
(n+1)/N
a/N
(N/N SV:
n=0
Jn/N
n=0
where the second inequality follows from Proposition 1
24

---

Published as a conference paper at ICLR 2026
Proposition 2 (Proposition 3.2.14, Grafakos (2008).
If is of bounded variation, then
[(n) - (f)
72 o2 O. 2x n
Corollary 2. For n. E N, and & >> 0 fixed. the estimation of the Fourier series f(n) given by N-1 (() k j f(k/N)
|n| S N otherwise
converges uniformly to (n)).
Proof: For |n| < -2-=,, it follows by Corollary 2 that
S 22 |1
[1 2n
|m| + VV
(((m) (m)
N
2 f
2T
v2Vg
+
Ne
N
Similarly, for n 2 N1-s
Vol (F) 2 |12
Vi Vi)(F)
f(n) - (m)
-
S
2T Nl-e'
The proof of Theorem D.1 follows by letting 1/2.
25


Reviews:
[{"review_id": "xB2kZ2TSBA", "reviewer": "Reviewer_oKCG", "summary": "The authors tackle the problem of causal structure learning for systems where causal dependencies are non-static and dynamical in nature. They  provide a theoretical framework based on Wiener Filters (WFs) and Fast Fourier Transform (FFT) to introduce a novel method for deciphering non-static causal dependencies of dynamical systems from time-series data. They empirically validate their proposed approach on synthetic toy systems of 20 variables (nodes) and a real data (river-runoff) setting. This work provides a new perspective on causal structure learning of dynamical systems using frequency domain WFs.", "questions": "- In practical applications, how large does $N$ need to be to ensure reasonable recovery of the MEG? \n- On lines 224-225: \"For the settings with large number of samples (large $T$) or with longer delays (large $L$) improvements are significant\". What are the practical settings where we would apply the proposed method and where there is sufficiently large $L, T$ such that it would yield competitive performance?\n- Is the proposed approach practical for high dimensional systems (large $n$)?\n- For the synthetic experiments in section 6, are the dynamical / non-stationary dependencies of the ground truth process/DAG random? Could the authors provide some additional clarification here?\n- Lines 422-424: \"Bounds on sparsity metrics such as average, maximum, and minimum degree are typically available ...\". In practical applications and scientific applications, this is not necessarily know apriori, and thus, is not a reasonable assumption to make. Could the authors provide some examples supporting this claim? \n- Could the authors provide further details for how the Granger baseline is implemented? Further, have the authors also considered a very simple correlation baselines? i.e. compute the correlation between samples across adjacent time-points and use this as the estimated DAG?\n- The CS metrics seems like a weird choice. Why not consider something like average prevision / area under the average precision curve? You can compute this without thresholding, and it would provide a metric that captures some sort of notion of accuracy across \"all thresholds\".\n\n\nMinor comments:\n\n- I would define the acronym FFT (Fast Fourier Transform) earlier in the text (in fact I don't think its defined anywhere currently).\n- Notation for $T$ slightly over-loaded. It is used to indicate transpose and in the number of time-series samples. \n\nReferences:\n\n[1] Huang et al. \"Causal discovery from heterogeneous/nonstationary data.\" JMLR. 2020.\n\n[2] Fujiwara et al. \"Causal discovery for non-stationary non-linear time series data using just-in-time modeling.\" Conference on Causal Learning and Reasoning. PMLR, 2023.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The authors introduce a novel, computationally efficient, method for causal structure recovery of non-static / non-stationary systems (systems with dynamic causal dependencies).\n- The paper provides a novel way of thinking about causal structure learning of dynamical causal dependencies from time-series data; i.e. causal structure learning in the frequency domain. This insight leads to the aforementioned improvements in computational efficiency.", "weaknesses": "- Although in general this work is reasonably well written, and the authors did a good job at presenting the many nuanced details present in this work, this paper can still be hard to follow at times. \n- The baselines in the empirical experiments seem to primarily be methods which operate in the static (stationary) setting. It would be useful to compare to some existing works that approach the problem of dynamic and non-stationary dependencies. To give some examples: [1, 2].\n- Similarly, empirical experiments are quite limited. Although I believe this work presents comprehensive theoretical backing for their approach, which provides a reasonable degree of contribution already, the authors could consider some of the datasets/systems from [1, 2] to add to their  experiments section. Or perhaps, simulating data in some of the example systems the authors provide in Appendix C.6 could also be considered for addition experiments. I believe additional empirical experiment would help further strengthen this work.", "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": "kR94jW1181", "reviewer": "Reviewer_TpiQ", "summary": "This paper investigates causal structure recovery in linear dynamical systems using Wiener Filter (WF) based approaches. The authors compare time-domain (TD) and frequency-domain (FD) methods for computing WF coefficients and propose two algorithms: Wiener-PC (extending PC algorithm with WF-based conditional independence tests) and Wiener-Phase (exploiting phase properties of WF for efficient MEG reconstruction). The work provides theoretical analysis including concentration bounds and sample complexity for both approaches, demonstrating computational advantages of FD methods.", "questions": "1: How can practitioners verify whether Assumptions 1 and 2 hold for their specific applications? What happens to algorithm performance when these assumptions are violated beyond the limited robustness study provided?\n\n2: Given the O(n^(q+3)) complexity of Wiener-PC, how does the method perform on networks with hundreds or thousands of nodes? What are the practical limits of the approach?\n\n3: While Lemma A.1 suggests any frequency can be used, how does the choice of frequency affect performance in practice, especially in noisy or finite-sample settings?", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1: The Wiener-Phase algorithm cleverly exploits the phase properties of frequency-domain WF coefficients to efficiently identify colliders and reconstruct MEGs, which has no analogue in time-domain approaches.\n\n2: The evaluation includes synthetic data with controlled network structures, scalability analysis up to 50 nodes, and validation on real-world river-runoff dataset, showing consistent advantages of FD approaches.", "weaknesses": "1: The Wiener-Phase algorithm relies on two strong assumptions (phase alignment of incoming edges and non-zero imaginary components) that significantly limit its applicability. While the authors claim these apply to \"a large class of networks,\" the conditions are quite specific and may not hold in many real-world scenarios.\n\n2: The comparison primarily focuses on traditional static methods (Fisher-Z, chi-square tests) and Granger causality, missing more recent state-of-the-art causal discovery methods for time series data.\n\n3: I'm unclear on where the distinction between static and dynamic approaches is specifically presented in your algorithm, and your specific contribution is unclear.\n\n4: While the paper shows results for up to 50 nodes, the computational complexity analysis suggests the methods may not scale well to truly large networks due to the combinatorial nature of PC-based approaches. Please increase the credibility and workload.", "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": "ELwAGiQ902", "reviewer": "Reviewer_eWah", "summary": "The work compares time- versus frequency-domain estimation of multivariate Wiener filters for causal structure recovery under a linear dynamical influence model, proves computational gains via FFT, provides finite-sample concentration bounds for both domains, and proposes Wiener-PC (CI testing via Wiener coefficients) and Wiener-Phase (phase/imaginary-part–based skeleton and collider discovery), validated on synthetic and real river-runoff data to support the claims.", "questions": "See Weaknesses", "limitations": "", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "- S1: Uses the zero/nonzero support of multivariate Wiener coefficients as a principled CI surrogate, with FFT yielding explicit complexity gains over time-domain least squares, and finite-sample bounds provided in both domains.\n\n- S2: Introduces Wiener-Phase, which leverages the imaginary part/phase of frequency-domain Wiener filters to recover the skeleton and identify colliders under phase assumptions, with strong empirical scalability.\n\n- S3: Provides broad empirical evidence on synthetic node networks and a real-world dataset, showing superior accuracy and runtime trade-offs versus multiple baselines.", "weaknesses": "- W1: Key guarantees for Wiener‑Phase rely on structural assumptions (e.g., equal phase of all incoming edges and a nonzero imaginary part of transfer functions), which are somewhat strict.  \n\n- W2: The approach hinges on heuristic thresholding of Wiener coefficients, tuned via sparsity/degree constraints rather than statistically calibrated procedures, which may make results sensitive to sampling and frequency discretization.  \n\n- W3: Baseline coverage omits widely used frequency-domain causal discovery families such as spectral Granger variants, making it harder to position the proposed methods relative to established spectral measures beyond time-domain Granger.", "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)