Spectral Ranking of Causal Influence in Complex Systems.

Similar to natural complex systems, such as the Earth's climate or a living cell, semiconductor lithography systems are characterized by nonlinear dynamics across more than a dozen orders of magnitude in space and time. Thousands of sensors measure relevant process variables at appropriate sampling rates, to provide time series as primary sources for system diagnostics. However, high-dimensionality, non-linearity and non-stationarity of the data are major challenges to efficiently, yet accurately, diagnose rare or new system issues by merely using model-based approaches. To reliably narrow down the causal search space, we validate a ranking algorithm that applies transfer entropy for bivariate interaction analysis of a system's multivariate time series to obtain a weighted directed graph, and graph eigenvector centrality to identify the system's most important sources of original information or causal influence. The results suggest that this approach robustly identifies the true drivers or causes of a complex system's deviant behavior, even when its reconstructed information transfer network includes redundant edges.

1. Introduction

Semiconductor lithography systems are extremely complicated electromechanical systems, capable of sub-nanometer positioning and sub-milliKelvin temperature control, while generating extreme ultraviolet light from laser-pulsed tin plasma. It is notoriously difficult to fully understand the complex interactions among thousands of observed variables affecting the output of such systems. Model-based approaches alone are inadequate to effectively diagnose rare or new system issues, as they inherently do not model abnormal behavior. To efficiently, yet reliably, reduce the search space of potential causes, we consider both local and global causal influences of a complex system’s components, as measured by Schreiber’s transfer entropy [1] and Bonacich’s eigenvector centrality [2], respectively. Transfer entropy quantifies predictive information transfer as a potential signature of causality between two stationary time series, providing direction, strength and delay of linear or non-linear interactions. For time series with a Gaussian distribution, transfer entropy reduces to Granger causality [3]. However, being a bivariate measure, transfer entropy naturally disregards the multivariate nature of interactions in complex systems. Therefore, multivariate approaches use conditional transfer entropy [4] to separate true cause–effect relations from mere correlations, i.e., direct from so-called transitive indirect or semi-metric, and thus redundant relations, by iteratively conditioning out (subsets of) all other time series. Unfortunately, such an information decomposition is infeasible in (real time) diagnosis or prognosis of high-dimensional technological complex systems, due to its exponential scaling of computational costs with time-series dimension. Instead, we use a ranking algorithm that relies on standard transfer entropy in exhaustive, but computationally feasible bivariate interaction analysis of a complex system’s multivariate time series, resulting in an information transfer network that is likely to contain redundant edges. However, Benzi et al. [5] show that eigenvector centrality represents the upper limit of Katz centrality [6], where the global influence of nodes is essentially determined by the longest possible paths in the network, which implies its insensitivity to shortest paths, including transitive indirect or semi-metric, and thus redundant edges. In addition, empirical studies by Kalavri et al. [7] revealed that PageRank (or eigenvector centrality) yields similar rankings when computed with, or without, a graph’s (first-order) semi-metric edges. For diagnostic purposes as envisioned, we suffice to apply standard transfer entropy for approximative network inference of a complex system’s true causal structure and then use graph eigenvector centrality for computationally efficient, yet consistently accurate ranking of probable causes or key drivers of the system under disturbance.

This study relates most closely to the work of Streicher et al. [8], who implemented this approach into a two-step ranking algorithm for (chemical) plant-wide fault detection and diagnosis.

In this paper, we introduce the measures of transfer entropy and eigenvector centrality and describe the ranking algorithm as used. In rolling window analysis of a simulated time series representing two bidirectionally coupled Lorenz systems, we compare its causal network inference results to those of a constraint-based algorithm for multivariate causal analysis. The rolling analysis approach also allows the ranking algorithm to estimate the global influence of each state variable as exerted on the coupled Lorenz systems over time. Finally, we assess the ranking algorithm in diagnosing a real-world industrial system issue, using a higher-dimensional time series.

2. Applied Algorithms and Methods

2.1. Transfer Entropy

Transfer entropy () is an information-theoretic implementation of Wiener’s notion of causality applied to time series [9], whereby the cause precedes—and contains unique information about—the effect. Consider two stationary ergodic Markov processes X and Y (or and as below) and their corresponding time series and of M samples. Transfer entropy quantifies reduction in uncertainty about future states of a source process X, when passed states of a target process Y are observed in addition to passed states of X itself. As an asymmetric measure based on transition probabilities, transfer entropy naturally incorporates directional and dynamic information, which may imply causation between X and Y: where and represent states of X and Y at time t, while indicates maximized information transfer from X to Y, computed across a range of embedding delays , such that . The embedding dimensions k and l denote the number of passed states in X and Y used to condition the probabilities of transition to the next state of Y (or X) represented by and . The logarithm base defines the informational unit of transfer entropy in bits, i.e., reduction in average code-length required to optimally encode the target variable (effect), given passed states of the source variable (cause) and target variable. Herein, we keep the embedding dimensions at the commonly used , mostly for computational reasons. For every pair from a multivariate time series , we apply Equation (1) to estimate information transfer and thereby determine the candidate source times series and target time series . We then generate surrogates , which share their amplitude distribution and power spectrum with original time series , using the iterative amplitude adjusted Fourier transform (iAAFT) proposed by Schreiber et al. [10]. Following this study, we estimate information transfer from each surrogate to (original) target time series and obtain . If , information transfer is considered to be significant or non-significant otherwise. The resultant information transfer network, given by a directed weighted graph , comprises a set of N nodes and of edges. Each edge connects a source node to target node with strength , i.e., information transfer .

2.2. Eigenvector Centrality

To diagnose performance issues or even failures within technological complex systems, we wish to locate where perturbations enter the system and propagate, causing downstream effects throughout the system. Therefore, we consider an information transfer network wherein the (main) sources of original information can be identified by measurement and ranking of each node’s global network influence termed centrality or importance. Centrality, the basic principle of Google’s search engine [11], has proven useful to measure and rank a page’s relevance based on its inbound links. Here, we use out-degree eigenvector centrality, which defines the centrality of a node as proportional to the summed centralities of its outbound neighbors: where is a proportionality factor and c is the eigenvector of centralities associated with eigenvalue of adjacency matrix W, whose entries denote information transfer from node to node in graph G. If W is non-negative and irreducible, the Perron-Frobenius theorem [12,13] ensures there is a unique vector of N centralities associated with the largest positive eigenvalue or spectral radius of W, satisfying Equation (2). Usually is normalized, such that each entry indicates the centrality or importance of a node in graph G on a relative scale from 0 to 1. Alternatively, matrix W is normalized to a transition matrix P, whose entries denote probabilities of transition from node to node in a random walk on graph G or Markov chain while , such that:

Following Google’s PageRank approach [14], we modify matrix P by adding a teleportation probability and an all-ones matrix J to obtain an , irreducible, positive matrix : where a random walker at node follows an edge with probability or jumps to any other node in the N-nodes network with probability . Herein, we use PageRank’s typical value of . Considering the Markov chain associated with matrix , the Perron–Frobenius theorem ensures an unique stationary probability distribution that matches the eigenvector of associated with eigenvalue , such that . Eigenvector is usually computed via power iteration [15] in k steps: , where denotes an initial distribution.

We use the FaultMap algorithm (https://github.com/SimonStreicher/FaultMap, accessed on 20 May 2015) from Streicher [8], which, to our knowledge, is the only open-source implementation of the method as described above and summarized in the pseudo-code of Algorithm 1, as shown below.

Input:N-dimensional time series , of M samples from system X represented by: , statistical significance threshold in, e.g., a rank-order test using iAAFT surrogates, embedding delay range and embedding history lengths k and l for time series i and time series

Output: adjacency matrix , where entry () represents information transfer from node i to node j

for to N do

for to do

for to do

compute by Equation (1)

if then

else

end if

end for

end for

end for

Input: matrix where , ranking distance

Output: node centrality score vector

initialize with probabilities

whiledo

compute eigenvector of matrix associated

with eigenvalue , such that

end while

2.3. Validation

In what follows, we assess the method’s accuracy in causal inference and node ranking using simulated time series, followed by root cause analysis from real-world diagnostic data. Firstly, we consider a system of two bidirectionally coupled Lorenz systems and , investigated by Wibral et al. [16] and given by: The bidirectional coupling is governed by time-delayed quadratic terms. We used Pydelay [17] to generate a multivariate time series of 150K samples by numerically integrating Equations (5a)–(5f) with step size and initial conditions , and . To assess FaultMap’s accuracy in causal network reconstruction against other model-free approaches, we use PCMCI (v4.0) from Runge [18] as a distinct, constraint-based, multivariate alternative. FaultMap estimates edge-weight as using the Kraskov–Stögbauer–Grassberger estimator for transfer entropy in the Java Information Dynamics Toolkit from Lizier [19]. PCMCI performs condition selection at an optimized significance level using Akaike’s information criterion, followed by a conditional independence test where we use the linear ParCorr test (at ) instead of the computationally expensive nonlinear CMI test. Following the recommended sample size of Bauer [20], we use a rolling analysis window of 2K samples to define 50 adjacent time series slices for causal network reconstruction of the coupled Lorenz systems by both algorithms and node importance ranking by FaultMap only.

3. Results and Discussion

3.1. Coupled Lorenz Systems

Figure 1 depicts the rolling window analysis approach, the butterfly-shaped attractors of the coupled Lorenz systems in phase space and cause–effect detection count heatmaps for both algorithms. Figure 2 exemplifies FaultMap’s rolling window analysis results of which Figure 3 specifically reports information transfer (delay) via coupling and global influence of X-, Y- and Z-variables on the coupled Lorenz systems. The heatmaps in Figure 1b show the total count of each cause–effect relation detected by FaultMap or PCMCI over 50 adjacent time series slices. Firstly, PCMCI’s linear ParCorr test detected nonlinear coupling remarkably well, with a rate of 85% vs. 89% by FaultMap and achieved a notably higher detection rate (95%) of (direct) causal links throughout the coupled Lorenz systems than FaultMap (83%). ParCorr’s sensitivity to nonlinear interactions has been previously discussed by Krich et al. [21]. PCMCI’s higher detection rate of direct links compared to FaultMap is mainly due to its 100% detection rate of self-influence at each system state variable, whereas FaultMap only reached a 100% detection rate of self-influence at the coupled system state variables and . We could not find an obvious explanation for the algorithms’ differing detection rates of self-influence (see the heatmaps’ diagonal), but may argue that self-influence as subset of the previously defined edge-set E, or formally , is irrelevant in the context of Equation (2). Lastly, note that unlike what we would expect from its multivariate causal network reconstruction designed to distinguish direct from indirect effects, PCMCI detected 23% more transitive indirect links than FaultMap.

Figure 1

(a) Time series of system composed of bidirectionally delay-coupled Lorenz systems and , generated from Equations (5a)–(5f). (b) Heatmaps of detection count per cause–effect relation, for FaultMap (left) and PCMCI (right). Direct cause–effect relations are denoted by (*).

Figure 2

Information transfer network of two bidirectionally delay-coupled Lorenz systems. Edges indicate direct (⟶) or transitive indirect (⤏) information transfer. Edge annotations denote information transfer delay (sec). Node importance indicates a node’s global network influence. Edge-weight represents level of information transfer ().

Figure 3

Dynamic information transfer via bidirectional delay-coupling of Lorenz systems and , and dynamic importance of the Lorenz system state variables. (a) Dynamic information transfer (delay) via bidirectional delay-coupling between Lorenz systems and . (b) Distribution of information transfer (delay) in 3a. (c) Dynamic importance of state variables in delay-coupled Lorenz systems and . (d) Distributions of Lorenz system state variable importance in 3c.

Due to the transitivity of bivariate information transfer, most (if not all) networks inferred by FaultMap include direct and indirect connections, all of which passed strict significance tests, as and shown in Figure 2. The coupled Lorenz systems are easily discernible by two subnetworks and that exhibit distinct levels of information transfer with similar time delays in the range of milliseconds. At a rate of at least 98%, FaultMap detected the same edges , , and indicating (normalized) net influence within both Lorenz subsystems, as Gencaga et al. [22] found for a single Lorenz system. PCMCI shows similar detection results for these edges (see heatmaps). It is important to note that the rolling window approach took PCMCI 495 min of runtime on a 24-core HPC node with 192 Gb for causal analysis vs. 58 min by FaultMap for a two-step causal analysis and node ranking. PCMCI’s considerably higher computational costs compared to FaultMap obviously limit its applicability for the diagnostic purposes considered above.

In Figure 3a,b, we focus on reconstruction of time delays in coupling , given the modeled delays in Equations (5b) and (5e). Figure 3a reveals the dynamic nature of bidirectional information transfer in terms of strength and delay. The median difference (≈0.05) of information transfer distributions in Figure 3b suggests that Lorenz system predominantly drives , which is reasonable to expect since the coupling strength () in Equation (5b) is twice the coupling strength () in Equation (5e). Given coupling delays of 3 and 5 s, the distribution of reconstructed interaction delays seems realistic but may be impacted by lag synchronisation of the Lorenz systems, as suggested by Coufal et al. [23]. Figure 3c shows highly dynamic global influence of particularly X- and Y-variables, while the Y-variables remain the driving force within their respective Lorenz subsystem at all times. Grey-colored bars highlight time windows in which is identified as driving, and as driven subsystem. Since the coupling strengths are constant, is likely to dominate in strength when the state reaches vanishingly low values relative to the state. The median difference across all node importance distributions in Figure 3d reflects the aforementioned drive–response relations in, and between, the interacting Lorenz systems. It might explain FaultMap’s 100% detection rate of Y-variable self-loops vs. lower rates of all other self-loops. The outcome of both Y-variables as Lorenz system key driver complies with the Lorenz model of Rayleigh–Bénard convection [24], where temperature difference drives convective heat transfer in addition to conduction at Rayleigh number (see Equations (5b) and (5e)). To our knowledge, this is the first rolling window analysis to date, capturing the dynamics of time-varying information transfer or global influence (importance) of state variables within interacting Lorenz systems. Our findings may enable automated identification of monitoring observables for performance (anomaly) diagnostics or predictive maintenance within technological complex systems. The ability to capture time-varying importance of a complex system’s state variables is also relevant in time-series analysis of natural complex systems, including the Earth’s climate.

3.2. Technological Complex Systems

Nanolithography systems are among the most complex technological systems today, capable of sub-nanometer positioning and sub-milliKelvin temperature control, even as system modules accelerate at up to 15 Gs. Such systems are particularly challenging for model-based diagnosis of rare or new issues, due to nonlinear interactions across multiple time and spatial scales. To assess FaultMap’s potential in diagnosis of issues within such systems, we investigate temperature, flow and pressure instability within an ASML subsystem. Therefore, we use a multivariate time series of 315 binary samples from 366 parameters related to the problem. As shown in Figure 4, FaultMap identified parameter as a primary source of original information, i.e., most probable cause leading to event , through a network of collateral effects . The indicated root cause is confirmed to be correct by a series of automatically logged system events as well as service actions. Interestingly, the event log messages and follow the network’s direction of time i.e., from cause to effect (causal inference), while the logged service actions related to and follow the reversed time order, i.e., from effect to cause (diagnostic inference). Hence, the last logged service action appears as a direct effect of root cause in the network. This observation is promising, with regard to the automation of reliable data-driven diagnostics for technological complex systems.

Figure 4

Top-ranked node (root cause) transfers original information towards event via a network of collateral effects within an ASML subsystem. (Figure 2 legend applies here.)

4. Conclusions

To fully understand a complex system’s dynamical behavior, it is essential to identify its main sources of causal influence affecting downstream elements throughout the system. We empirically show that spectral centrality analysis of its causal network as approximated by standard transfer entropy allows one to accurately and consistently identify the most important node(s) of original information representing the most probable cause(s) or driver(s) of disturbance in the system. The ranking algorithm we use compares favorably against the alternative algorithm for multivariate information transfer estimation in causal analysis of two nonlinearly coupled Lorenz systems. In addition, it shows to be accurate, consistent and efficient, identifying the alternately driving and driven Lorenz subsystems, as well as the driving force within either subsystem, over time. Finally, the ranking algorithm correctly traces back the original disturbance within a high-dimensional technological complex system from sampled time series of several hundreds of parameters. Considering the high-dimensionality of observations across multiple time and spatial scales from such systems, we conclude that the inherent robustness of spectral centrality to semi-metricity of directed networks makes it a viable option for reliable and scalable diagnostics. Additionally, spectral centrality ranking allows for feature selection and is particularly useful in identifying long-term effects.

Regardless of the computational costs, state-of-the-art multivariate causal inference methods may be the better choice to account and control for unobserved variables or capture synergistic interactions. However, comprehensive comparison of our method with these approaches is beyond the scope of this study and therefore recommended for future research.

We thank David Sigtermans for our many inspiring discussions and Leonardo Barbini for his valuable feedback on the causal analysis results of the Lorenz systems. Finally, we are grateful to Simon Streicher for his implementation and helpful suggestions using it.