A method for identifying causality in the response of nonlinear dynamical systems

## Abstract

Abstract Predicting the response of nonlinear dynamical systems subject to random, broadband excitation is important across a range of scientific disciplines, such as structural dynamics and neuroscience. Building data-driven models requires experimental measurements of the system input and output, but it can be difficult to determine whether inaccuracies in the model stem from modelling errors or noise. This paper presents a novel method to identify the causal component of the input-output data from measurements of a system in the presence of output noise, as a function of frequency, without needing a high fidelity model. An output prediction, calculated using an available model, is optimally combined with noisy measurements of the output to predict the input to the system. The parameters of the algorithm balance the two output signals and are utilised to calculate a nonlinear coherence metric as a measure of causality. This method is applicable to a broad class of nonlinear dynamical systems. There are currently no solutions to this problem in the absence of a complete benchmark model.

## 1 Introduction

Nonlinear dynamical systems are observed across a wide range of scientific disciplines (Mosekilde, ), from structural dynamics (Nayfeh & Mook, ; Worden & Tomlinson, ) and electrical engineering (Kyamakya et al., ), to neuroscience (Corinto & Torcini, ) and economics (Bischi, ). Modelling their behaviour offers meaningful insights, enables predictions, and allows for control to generate desirable outcomes.
Accurate models generally involve experimental measurements, e.g. to carry out model updating for physics-based models, or to train data-driven models. However, it can be difficult to determine whether the prediction errors are due to additional noise or inaccuracies in the model. Consider a nonlinear dynamical system subject to broadband random excitation, $x(t)$, which is especially prevalent in structural dynamics. The output of the system is measured to be $y_{n}(t)$, where

$$ $y_{n}(t)=y(t)+\epsilon_{n}(t)=\mathcal{F}\{x(t)\}+\epsilon_{n}(t)\ .$ (1) $$

Here $\mathcal{F}\{\cdot\}$ represents the true system, $\epsilon_{n}(t)$ is additive noise and $y(t)$ is the component of $y_{n}(t)$ that is caused by $x(t)$. The noise is assumed to be due to additional separate signal sources rather than measurement noise, hence there is no noise on the input. A forward model can be used to make a prediction $y_{z}(t)$, where

$$ $y_{z}(t)=\mathcal{G}\{x(t)\}=\mathcal{F}\{x(t)\}+\epsilon_{z}(t)\ .$ (2) $$

Here $\epsilon_{z}(t)$ denotes the model errors which are caused by the inability of the model, $\mathcal{G}\{\cdot\}$, to fully capture the function $\mathcal{F}\{\cdot\}$. Determining the component of $y_{n}(t)$ caused by $x(t)$ is important for applications which depend on the ability to predict $y_{n}(t)$. However, this is challenging because it is difficult to distinguish between $\epsilon_{n}(t)$ and $\epsilon_{z}(t)$.
Despite its wide ranging applications, this work is specifically motivated by feedforward active noise reduction (ANR), where $x$ is used as a reference signal to reduce the target signal, $y_{n}(t)$. ANR is very popular in modern earbuds, headphones and aviation headsets, but has recently begun implementation in cars to reduce road noise (De Brett, ; Oh et al., ). In order to reduce the unwanted noise, a model must be developed to accurately predict $y_{n}(t)$ using $x(t)$. The maximum noise reduction that can be achieved is quantified by the causal relationship between $x(t)$ and $y_{n}(t)$, which determines the potential application of the technology in new products. For example, if the causal relationship between $x(t)$ and $y_{n}(t)$ is weak, then the level of noise reduction that is possible may not be noticeable, however good the model is.
The input-output causal relationship can be quantified using a nonlinear coherence metric between $x(t)$ and $y_{n}(t)$. This is also equal to the linear coherence between $y(t)$ and $y_{n}(t)$ (Silva et al., ), but cannot be calculated directly due to the unavailability of $y(t)$, and so only a lower bound on the nonlinear coherence can be estimated using a model prediction, $y_{z}(t)$. Determining the true nonlinear coherence is challenging because it is difficult to distinguish between $\epsilon_{n}(t)$ and $\epsilon_{z}(t)$, which both appear as random noise. When the input to the system can be controlled, it is possible to estimate the contribution of $y(t)$ through repeat measurements (see Appendix A). However, when the input cannot be controlled, there are currently no methods which can infer the nonlinear coherence without a perfect model of $\mathcal{F}\{\cdot\}$, which can require large models and lots of data to compute. Ideally, it could be estimated using a small quantity of data to determine the feasibility of a project before investing in data collection and computational resources.
This work aims to estimate the nonlinear coherence between $x(t)$ and $y_{n}(t)$ for a particular class of nonlinear dynamical systems, with relatively small quantities of data and an incomplete model. These systems can be described by the following second order ordinary differential equation (ODE):

$$ $\frac{1}{\omega_{n}^{2}}\ddot{y}+\frac{2\zeta}{\omega_{n}}\dot{y}+y+\mathcal{N }\left(y,\dot{y},\ddot{y}\right)=x+\dot{x}+\ddot{x}\ ,$ (3) $$

where $\mathcal{N}\left(y,\dot{y},\ddot{y}\right)$ represents an arbitrary nonlinear function. Nonlinear dynamical systems of this type are common in many scientific disciplines and there has been significant work in modelling their response using data-based machine learning methods (Wang et al., ; Liang et al., ; Worden et al., ; Massingham et al., ). However, it is extremely difficult using feedforward modelling techniques (e.g. convolutional neural networks) due to the exponentially growing number of nonlinear monomial terms which contribute to the output. This perspective is captured by the Volterra series (Volterra, ). It is especially challenging for systems with long memory and for random broadband inputs. In this paper, the structure of Eq. 3 is exploited to estimate the nonlinear coherence for systems of this type, without a perfect model of the system. This enables the calculation of the causal component of a system response, which can be used to guide the investment of time and money into future research. To the best of the author’s knowledge, there are currently no other methods that achieve this. The key contributions of the work are as follows:

- •
A completely new strategy is introduced for estimating the maximum possible performance of time series prediction algorithms in the presence of noise.
- •
A novel method is presented for calculating the nonlinear coherence between input and output data for a general class of nonlinear ODEs subject to random input in the presence of noise.
- •
The robustness of the method is demonstrated on three simulated nonlinear dynamical systems which are widely applicable across several scientific disciplines.
- •
The applicability of the method is demonstrated on an experimental dataset with a strong nonlinearity.

## 2 Related Work

The presented method for calculating the level of noise present and inferring the maximum performance of feedforward models is very different to work in the literature. However, there are three broad areas of research in the literature that overlap with, and contribute to, the goals of this paper: the Kalman filter estimates the state of a system using the predictions of a model and noisy measurements; scaling laws are used to predict the performance of a neural network as the quantity of training data and model size increases; and the Error-in-Variable approach is used to infer the level of noise in the input to a neural network.
Kalman Filter
The Kalman filter (Kalman, ), and its extensions, aim to estimate the true state of a system by optimally combining the predictions of a state space model with noisy measurements. A comparison of the model error and noise level is used to balance the two signals optimally. This relies on knowledge of at least one of the model error or noise variance because the true system state is not available for hyperparameter tuning. In this paper, the primary aim is to estimate the noise level and the model error, but a similar principle is utilised.
Scaling Laws
Estimating scaling laws for various machine learning tasks has received significant attention in recent years (Hestness et al., , ; Henighan et al., ), with the aim of forecasting how the error of various models will improve as more data is collected and larger models with more parameters are built. However, this approach is limited, as it relies on using the performance of a model trained on a small quantity of data to extrapolate and predict the performance of a bigger model trained using large amounts of data. Extrapolation is difficult, as it assumes there is no noise; the loss will plateau at the noise level. It is therefore difficult to evaluate the maximum performance unless the model is already close to the maximum, at which point increasing the model size and number of training examples will not reduce the prediction error. For the method presented in this paper, even with a poor forward model, it is possible to estimate the maximal predictive performance of a feedforward model by identifying the level of noise.
Learning Input Noise Levels
As discussed above, if there is noise added to the output, it is difficult to determine the noise level because it is difficult to differentiate between unmodelled components of the system response and noise. However, if there is noise in the input, there are methods for inferring the true input and hence the input noise level. The Error-in-Variable approach can be used to infer the input to a machine learning algorithm and hence infer the noise level (Van Gorp J, ; Seghouane & Fleury, ). Therefore, as in the method presented in this paper, by inverting the problem and trying to predict the input to the system using the output, the output noise could be determined. The issue with applying this to time series data is it requires inferring the true system response for potentially millions of time samples, and so the inference is extremely high dimensional. This paper proposes an alternative approach to this problem, which is well suited to time series applications.
Other methods have been presented in the literature that aim to identify a ‘nonlinear’ coherence. However, these methods generally rely on a perfect forward model of the system and so model error and additive noise are combined in any estimation (Worden et al., ). Our approach is distinct in that only an approximate forward model is needed.

## 3 Methodology

### 3.1 Architecture

Consider the signals $x(t)$, $y_{n}(t)$, and $y_{z}(t)$. Here $x(t)$ is the input to the system, $y_{n}(t)$ is a noisy observation and $y_{z}(t)$ is the prediction of an imperfect forward model. The selection of the forward model is discussed in Section 3.3. The noise is represented by $\epsilon_{n}(t)$, which is assumed to be zero mean but can otherwise be of an arbitrary spectrum, and the true output is denoted $y(t)$. The time series data is split into $N$ frames (input-output pairs), each of length $M$. The data is split into training and validation datasets, with $N_{T}$ and $N_{V}$ frames respectively. At each frequency, $Y^{(i)}(f),\ Y_{n}^{(i)}(f)\ \text{and}\ Y^{(i)}_{z}(f)$ represent the Fourier transforms of the true, measured, and predicted versions of $y$ at frequency $f$ for frame $i$. The terms $\mathcal{E}_{n}^{(i)}(f)$ and $\mathcal{E}_{z}^{(i)}(f)$ represent the zero mean noise and model error respectively at frequency $f$ in frame $i$:

$$ $\displaystyle Y^{(i)}_{n}(f)=Y^{(i)}(f)+\mathcal{E}_{n}^{(i)}(f)$ $\displaystyle Y^{(i)}_{z}(f)=Y^{(i)}(f)+\mathcal{E}_{z}^{(i)}(f)\ ,$ (4) $$

where $\mathcal{E}_{n}^{(i)}(f)$ is assumed to be uncorrelated with $Y^{(i)}(f)$, whereas $\mathcal{E}_{z}^{(i)}(f)$ is correlated with $Y^{(i)}(f)$. The explicit dependence on $f$ and $t$ will now be dropped for compactness: lowercase denotes the time domain and uppercase denotes the frequency domain.
A key observation of the structure of Eq. 3 is that whilst it is difficult to learn a mapping from $x$ to $y$ in the ‘forward’ direction, it is relatively simple to map from $y$ to $x$ in the ‘reverse’ direction, even for highly nonlinear systems with long memory. This is because there is no nonlinear mixing of terms in time in the reverse direction, and so there is not an infinite Volterra series (Volterra, ) as there is in the forward direction. Therefore, even for highly nonlinear systems with long memory, there is a relatively simple mapping from $y$ back to $x$. The input, $x$, can therefore be used as a reference signal to infer the true $y$ from $y_{z}$ and $y_{n}$. The signals $y_{z}$ and $y_{n}$ can be combined optimally to find the best estimation of $y$, and therefore $x$, based on their relative noise levels. Figure 1 illustrates the proposed architecture. In the frequency domain, consider a linear combination of $Y^{(i)}_{n}$ and $Y^{(i)}_{z}$ at each frequency:

$$ $\hat{Y}^{(i)}=Y^{(i)}_{n}K+Y^{(i)}_{z}\left(1-K\right)\ ,$ (5) $$

where $0\leq K(f)\leq 1$. In the time domain, this is then used to train a machine learning model to predict $x$ by minimising the cost function:

$$ $\mathcal{L}^{T}=(1-\lambda)\mathcal{L}^{T}_{x}+\lambda\mathcal{L}^{T}_{y}\ .$ (6) $$

The losses $\mathcal{L}^{T}_{x}$ and $\mathcal{L}^{T}_{y}$ are calculated as

$$ $\displaystyle\mathcal{L}^{T}_{x}=\sum_{i=1}^{N^{T}}\frac{||\mathcal{F}_{\theta }\left(\hat{y}^{(i)}\right)-x^{(i)}||^{2}}{\sigma_{x}^{2}}$ $\displaystyle\mathcal{L}^{T}_{y}=\sum_{i=1}^{N^{T}}\frac{||\hat{y}^{(i)}-y_{n} ^{(i)}||^{2}}{\sigma_{y}^{2}}\ ,$ (7) $$

where $\mathcal{F}_{\theta}$ has memory and represents the trained model with parameters $\theta$ and the superscript $T$ indicates the loss is calculated using the training dataset. Each term in the loss is normalised using $\sigma_{x}^{2}$ and $\sigma_{y}^{2}$, which are calculated as

$$ $\displaystyle\sigma_{x}^{2}=\frac{1}{N}\sum_{i=1}^{N}||x^{(i)}||^{2}$ $\displaystyle\sigma_{y}^{2}=\frac{1}{N}\sum_{i=1}^{N}||y_{n}^{(i)}||^{2}\ .$ (8) $$

The constant $\lambda$ plays an important role, as it balances the two terms in the loss function. The loss $\mathcal{L}_{x}$ drives the architecture towards predicting $x$ accurately, whereas $\mathcal{L}_{y}$ drives the architecture towards $\hat{y}=y_{n}$. This affects the optimal $\hat{y}$, and hence the optimal $K$. When $\lambda$ is set to 0 0 then the only term contributing to the loss function is $||\mathcal{F}_{\theta}\left(\hat{y}^{(i)}\right)-x^{(i)}||^{2}$ and it is hypothesised that the optimal value of $K$ will also yield the optimal estimate of the noise-free output $Y$. This optimum $K$ corresponds to minimising:

$$ $J=E\left[|\hat{Y}-Y|^{2}\right].$ (9) $$

Minimising this expectation with respect to $K$, by setting $\frac{\partial J}{\partial K}=0$, yields:

$$ $K=\frac{1}{1+\frac{E\left[\mathcal{E}_{n}\mathcal{E}_{n}^{*}\right]}{E\left[ \mathcal{E}_{z}\mathcal{E}_{z}^{*}\right]}}\ .$ (10) $$

See the Appendix B.1 for more details. This is the true optimum value of $K$ for a given frequency, but it cannot be computed directly because $\mathcal{E}_{n}$ and $\mathcal{E}_{z}$ are unknown. However, under the hypothesis that $\lambda=0$ yields the true optimum for $K$, then it is expected that $K$ will converge towards this value in order to minimise $\mathcal{L}$. Therefore if the parameters of the model, $\theta$ and $K$, are trained simultaneously, the architecture will implicitly estimate the ratio of errors in order to optimally combine the two signals. This is used to calculate the nonlinear coherence in Section 3.2.
In practice, $\lambda=0$ does not yield the optimal $K$ and the role of $\lambda$ will be discussed in detail in Section 3.4.

Figure: Figure 1: Architecture. The fast Fourier transform (FFT) of the signals $y_{z}$ and $y_{n}$ is taken, before the two signals are combined at each frequency point using the parameter $K(f)$. The inverse fast Fourier transform (iFFT) is then applied to the resultant signal, which is then mapped to $x$ using a CNN.
Refer to caption: x1.png

In this paper, $\mathcal{F}_{\theta}$ is a 1D convolutional neural network (CNN). Due to the simplicity of the reverse mapping from $y$ to $x$, a relatively small CNN can be constructed, even when the system memory from $x$ to $y$ is relatively long. For the examples below, a CNN with 5 layers, a kernel width of 7 and 5 features was used. This network has 778 parameters and was applied in all the example test cases below.
The optimum value of $K$ could be estimated independently at each frequency, which would result in $M$ parameters, where $M$ is the frame length for the Fourier transform. However, it will vary with some structure across the frequency domain, and so could be represented using a smaller number of parameters, which will improve its learnability. There are many possible ways to represent $K$, such as using a neural network, gaussian process or a kernel based method. A linear piece-wise approximation is chosen here for simplicity, where $M_{c}$ uniformly spaced points are learnt: $K_{c}(f_{j})$ for $j=1:\frac{M}{M_{c}}:M$. Linear interpolation is then used between these points to obtain $K(f_{i})$ for $i=1:1:M$. The resulting function is passed through the sigmoid function in order to constrain it to lie between 0 and 1.

### 3.2 Nonlinear Coherence

The nonlinear coherence between $x$ and $y_{n}$ is equal to the linear coherence between $y$ and $y_{n}$:

$$ $\gamma^{2}_{YY_{n}}=\frac{|S_{YY_{n}}|^{2}}{S_{YY}S_{Y_{n}Y_{n}}}\ ,$ (11) $$

where $S_{AB}=E\left[AB^{*}\right]$ represents the cross spectral density between signals $A$ and $B$ at frequency $f$. Assuming $\mathcal{E}_{n}$ and $Y$ are uncorrelated, the following can be written.

$$ $\displaystyle S_{YY_{n}}$ $\displaystyle=E\left[YY^{*}_{n}\right]=E\left[YY^{*}\right]$ (12) $\displaystyle S_{Y_{n}Y_{n}}$ $\displaystyle=E\left[Y_{n}Y_{n}^{*}\right]$ (13) $\displaystyle S_{YY}$ $\displaystyle=E\left[YY^{*}\right]\ .$ (14) $$

Therefore

$$ $\displaystyle\gamma^{2}_{YY_{n}}=\frac{E\left[YY^{*}\right]}{E\left[Y_{n}Y_{n} ^{*}\right]}\ .$ (15) $$

Note that $E\left[Y_{n}Y_{n}^{*}\right]$ is not simplified further because this quantity can be measured from data. An estimate of $E\left[YY^{*}\right]$ is required. The estimation of $K$, using the architecture described in Section 1, can be used to reverse engineer an estimation of $E\left[YY^{*}\right]$. Rearranging Eq. 10, the ratio between the error in each signal is given as

$$ $C\equiv\frac{E\left[\mathcal{E}_{z}\mathcal{E}_{z}^{*}\right]}{E\left[\mathcal {E}_{n}\mathcal{E}_{n}^{*}\right]}=\frac{K}{1-K}\ .$ (16) $$

Therefore

$$ $\displaystyle CE\left[\mathcal{E}_{n}\mathcal{E}_{n}^{*}\right]$ $\displaystyle=E\left[\mathcal{E}_{z}\mathcal{E}_{z}^{*}\right]$ (17) $\displaystyle=E\left[(Y-Y_{z})(Y-Y_{z})^{*}\right]$ (18) $\displaystyle=E\left[YY^{*}\right]-2\Re\{E\left[YY_{z}^{*}\right]\}+E\left[Y_{ z}Y_{z}^{*}\right]$ (19) $$

In addition to this equation, three quantities can be calculated directly from the data: $E\left[Y_{n}Y_{n}^{*}\right]$, $E\left[Y_{z}Y_{z}^{*}\right]$ and $E\left[Y_{z}Y_{n}^{*}\right]$. Noting that $\mathcal{E}_{z}$ is correlated with $Y$ because it includes model error, two additional equations can be written:

$$ $\displaystyle E\left[Y_{n}Y_{n}^{*}\right]$ $\displaystyle=E\left[YY^{*}\right]+E\left[\mathcal{E}_{n}\mathcal{E}_{n}^{*}\right]$ (20) $\displaystyle E\left[Y_{z}Y_{n}^{*}\right]$ $\displaystyle=E\left[Y_{z}Y^{*}\right]\ .$ (21) $$

Combining Eq. 19 and Eq. 20 to eliminate $E\left[\mathcal{E}_{n}\mathcal{E}_{n}^{*}\right]$ gives the following.

$$ $\displaystyle C\left(E\left[Y_{n}Y_{n}^{*}\right]-E\left[YY^{*}\right]\right)$ $\displaystyle=E\left[YY^{*}\right]-2\Re\{E\left[YY_{z}^{*}\right]\}+E\left[Y_{ z}Y_{z}^{*}\right]\ .$ (22) $$

Solving for $E\left[YY^{*}\right]$ using Eq. 21 and Eq. 22 gives

$$ $E\left[YY^{*}\right]\approx\frac{CE\left[Y_{n}Y_{n}^{*}\right]-E\left[Y_{z}Y_{ z}^{*}\right]+2\Re\{E\left[Y_{n}Y_{z}^{*}\right]\}}{C+1}\ ,$ (23) $$

which can be calculated directly from the data using the inferred value of $K$ together with Eq. 16. This is then substituted into Eq. 15 to estimate the nonlinear coherence.

### 3.3 Generating the forward prediction

The architecture described in Section 1 requires a forward prediction of $y$. This may be the best forward prediction currently available but it can be a simpler model too. In this paper, the prediction of a temporal CNN trained on the available data is used. This is an arbitrary choice and any forward model could be used because the presented method uses the forward prediction to estimate an improved estimation of the nonlinear coherence. The better the available forward model, the better the estimation of the nonlinear coherence; this is demonstrated in the simulations below by applying the CNN to systems with varying strengths of nonlinearity. It has been found that it is preferable to use a more complex model, such as a CNN, so that the forward prediction cannot be easily mapped back to the input and the architecture drives towards $\hat{y}=y$, as discussed in Section 3.4.

### 3.4 Role of λ 𝜆 \lambda italic_λ

The solution for $K$ found by the architecture described in Section 1 with $\lambda$ set to 0 is usually smaller than the true value of $K$ across the entire frequency range. It is hypothesised that this is because $y_{z}$ is generated directly from $x$, and so there is often a mapping directly from $y_{z}$ back to $x$, rather than by inferring $y$ first. The term $\lambda$ was therefore included to add some bias to drive the architecture towards inferring $y$. Its value is determined using the validation dataset. As $\lambda$ increases, initially $\mathcal{L}^{V}_{x}$ remains constant, before beginning to increase rapidly. The best estimation of the nonlinear coherence is found at the point when $\mathcal{L}^{V}_{x}$ begins to increase rapidly. In order to identify this point, $\lambda$ was increased incrementally and the model was trained for 100 epochs at each step. When the normalised loss $\mathcal{L}^{V}_{x}$, averaged over the past 100 epochs, increased by 0.01 above the minimum value, the previous value of $\lambda$ was chosen. If $\mathcal{L}^{V}_{x}$ does not increase, as is the case for low levels of noise, $\lambda=0.99$ is used. Example plots of $\lambda$ against $\mathcal{L}^{V}_{x}$ from the tests in Sections 4 and 5 can be found in Appendix D.

### 3.5 Training

Both the neural network weights and $K$ are trained using the Pytorch infrastructure (Paszke et al., ). A learning rate of 0.01 was used with the Adam optimisation algorithm. For $\lambda=0.001$, the architecture was trained for 2000 epochs. Then, for each subsequent value of $\lambda$, the algorithm was trained for 100 epochs. The full details of the implementation and training procedure are available on github (see the reproducibility statement).

## 4 Simulations

The method was evaluated on three simulated systems, each with different functional forms of nonlinearities in order to demonstrate the applicability and usefulness of the method. For each system, the robustness of the method was tested for different levels of noise. The spectrum of the noise was bandlimited and effectively flat around the resonant peak. The noise level can be quantified by the nonlinear coherence: values close to 1.0 indicate low noise and values close to 0 indicate high noise. In each case, the nonlinearity of the test case is quantified by the normalised mean squared error between the true response and the linearised response. The predictions of the presented algorithm were made using just 10 frames of data.
To benchmark the performance of the algorithm, the predictions of the nonlinear coherence are compared against those of a 1D temporal convolutional neural network with dilations, as described by van den Oord et al. (). These benchmark predictions were used as $y_{z}$ to demonstrate the improved estimation of the nonlinear coherence that can be calculated using the method presented in this paper. It is likely that better forward models could be trained both now and in the future for the examples shown; the aim is to demonstrate how the presented method utilises the forward model to improve the estimation of the nonlinear coherence, and so it is not benchmarked against the forward prediction in the traditional sense and the selection of a forward model is somewhat arbitrary.
In the presence of noise, even with a perfect model it is difficult to determine whether the entire response has been captured. The lower bound on the nonlinear coherence is $\text{Co}\left(Y_{z},Y_{n}\right)$ and the upper bound is 1. The true nonlinear coherence could lie anywhere within this range, and the results in this section demonstrate that the method presented in this paper provides a relatively accurate estimation that could be useful in practice.

### 4.1 Polynomial stiffness

Systems with quadratic and cubic nonlinearities appear in many mechanical (Nayfeh & Mook, ; Worden & Tomlinson, ) and biological systems (Miller, ), and $y^{2}$ and $y^{3}$ are the next terms in the Taylor expansion of the nonlinearity of any system. This example therefore has wide ranging real world applications. The following dimensionless nonlinear oscillator was simulated to generate 1000 frames of length 6000:

$$ $\ddot{y}+2\zeta\dot{y}+y+\alpha_{2}y^{2}+\alpha_{3}y^{3}=x\ .$ (24) $$

The input, $x$, was bandlimited random noise with a fixed rms $\tau$. The spectrum of $x$ was effectively flat in the range of interest. The hyperparameters were set as defined in Appendix C.1.
The linearised response of the system captured 88 % of the total system response and a CNN (kernel width=20, hidden layers = 5, features=10, trained for 1000 epochs with a learning rate of 0.01) captured $94\%$ of the response using 10 frames of data. Without the method presented in this paper, with just 10 frames of data it would be impossible to determine how much of the unmodelled response is due to noise and how much is due to nonlinearities. Figure 2 shows the true nonlinear coherence (solid black line), the prediction of the CNN in the forward direction (blue dotted line) and the linear coherence between $x$ and $y_{n}$ (solid green line). The more noise that was added, the lower the true nonlinear coherence because the causal relationship between $x$ and $y_{n}$ is weaker. The predicted nonlinear coherence was calculated using 10 frames of data, which is represented by the red dashed line. The method presented in this paper provides an excellent prediction of the nonlinear coherence. For this example system, the forward prediction captures a large component of the response; the method confirms this as in each case the predicted coherence is only slightly closer to 1 than that of the forward model. This is especially important in the high noise example, as if the forward model was trained on this noisy data, the loss would appear large and researchers may be lead to believe that they should collect more data and build larger models to improve the performance. However, the estimation provided by the presented method correctly indicates that most of the error is due to additive noise and only relatively small improvements to the model are possible.

Figure: Figure 2: Polynomial stiffness case study.
Refer to caption: x2.png

### 4.2 Saturating Stiffness

Saturating nonlinearities are common in hardware applications, such as Opamps, transistors, springs, and motors (Mellodge, ; Liu & Michel, ). The following dimensionless nonlinear oscillator containing a saturating nonlinear stiffness term was used to generate 1000 frames of length 6000:

$$ $\ddot{y}+2\zeta\dot{y}+\alpha_{1}y+\alpha_{2}\tanh{(10^{4}y)}=x\ .$ (25) $$

The input, $x$, was bandlimited noise with an rms of $\tau$. The spectrum of $x$ was effectively flat in the range of interest. The hyperparameters are set as in Appendix C.2.
The linearised response of the system captured 58 % of the total system response and a CNN captured 65 % of the response using 10 frames of data (kernel width=10, hidden layers = 5, features=3, trained for 100 epochs at a learning rate of 0.01), hence this is a strongly nonlinear system.

Figure: Figure 3: Saturating stiffness case study.
Refer to caption: x3.png

Figures 3 shows the nonlinear coherence estimation calculated using 10 frames of data and the forward prediction. In this case the forward model does not capture a significant component of the response, which is indicated by the method. For all three levels of noise, the presented method provides a significant improvement over the lower bound $\text{Co}\left(Y_{z},Y_{n}\right)$ and the upper bound of 1. The nonlinear coherence is underestimated for the high noise example because, in this case, the estimation of $\hat{y}$ that can be calculated using $y_{n}$ and $y_{z}$ is relatively poor. This biases the mapping learnt from $y$ to $x$. The estimation of the nonlinear coherence would be improved by building a more accurate forward model, but this would require more resources. However, the estimation of the nonlinear coherence, even in the high noise case, is still relatively good and the general shape of the curve is captured.

### 4.3 Nonlinear Friction

Friction is ubiquitous across engineering, but can be strongly nonlinear. The following dimensionless nonlinear oscillator containing a Coulomb friction term (Ellis, ) was used to generate 1000 frames of length 6000:

$$ $\ddot{y}+\alpha_{1}y+\alpha_{2}\tanh{(10^{4}\dot{y})}=x\ .$ (26) $$

The input, $x$, was bandlimited noise with an rms of $\tau$. The spectrum of $x$ was effectively flat in the range of interest. The hyperparameters are set as in the Appendix C.3.
The linearised response of the system captured 91 % of the total system response and a CNN captured 88 % of the response using 10 frames of data (kernel width=10, hidden layers = 5, features=3, trained for 1000 epochs at a learning rate of 0.01), hence this is a moderately nonlinear system.

Figure: Figure 4: Nonlinear friction case study.
Refer to caption: x4.png

Figures 4 shows the nonlinear coherence estimation calculated using 10 frames of data and the forward prediction. In this case, the forward model captures most of the system response and so only small improvements to the model are possible. The presented method demonstrates this and provides an excellent estimation of the nonlinear coherence for each noise case.

## 5 Experiment

Figure: Figure 5: Diagram of the experimental set up.
Refer to caption: x5.png

The method was then evaluated on an experimental dataset consisting of 360 frames of length 6000. A cantilever was excited through a nonlinear connection using a bandlimited random input. The nonlinear connection consisted of magnets (the cantilever is steel) and a rubber tip which caused rattling. This nonlinear connection is complex and has not been characterised. A single mode of the cantilever was excited, so it effectively had one degree-of-freedom. An accelerometer was placed on the shaker as a reference signal, $x$, and a second accelerometer was placed along the cantilever as a target signal, $y_{n}$. The set up was considered to be effectively noiseless and then three levels of additional noise were introduced using post-processing to represent uncorrelated noise in the system. Noise was introduced artificially because for this experimental rig it is difficult to introduce mechanical noise that is uncorrelated with the input due to strong coupling between inputs. Two minutes of data (10 frames) was used to generate the forward prediction and then predict the nonlinear coherence. See Figure 5 for the full experimental set up.
Given the nonlinearity in the system and small quantity of data captured, in the noise-free case only 45 % of the response was captured by the linearised component of the response, and a CNN (kernel width=10, hidden layers = 5, features=5) captured just 74 % due to the complexity of the system response. When noise is added, the challenge is to determine the component of $y_{n}$ that is caused by $x$, as this represents the component of $y_{n}$ that could theoretically be reduced using an ANC system.

Figure: Figure 6: Experimental case study.
Refer to caption: x6.png

Figure 6 shows the predicted and true nonlinear coherence for three levels of noise. For all three levels of noise, the predicted nonlinear coherence was an excellent estimation of the true nonlinear coherence. In the low noise case, it demonstrates that significant improvements in the prediction could be achieved by collecting more resources and building larger models. Conversely, in the high noise case, the potential improvements are relatively small compared with the overall noise level. This demonstrates the insight provided by the presented method; in an ANR application, it may be worth investing resources in the low noise case to develop a new product, but not in the high noise case.

## 6 Limitations

The main challenge in the work is determining the optimal value for $\lambda$. The intuition is that $\mathcal{L}^{V}_{x}$ increases slowly as $\lambda$ increases until the nonlinear coherence prediction is close to the true nonlinear coherence, and then it begins to increase more rapidly. This can be seen in Appendix D. However, it is difficult to determine this transition point and accurately identify the optimal value for $\lambda$. Across the range of examples tested, thresholding $\mathcal{L}^{V}_{x}$ worked well empirically and the method was robust to small changes in the threshold value, but an interpretable method for setting $\lambda$ would improve this work.
Furthermore, another limitation relates to the forward model prediction, $y_{z}$. When learning the mapping from $y$ to $x$, noise causes bias in the learnt model and so leads to imperfect estimations of the nonlinear coherence. This can only be improved by reducing the level of noise in the estimation of $y$, which is achieved by either reducing the model error or reducing the additive noise. As a result, the better the forward model, the better the estimation of the nonlinear coherence. If the forward model is poor and the noise level is high, then the nonlinear coherence is likely to be under estimated, as in Section 4.2. However, it will still indicate that an improved forward model can be obtained.
In addition, the method relies on the availability of a relatively simple mapping of $y$ back to $x$. The type of systems investigated in this paper are widely applicable, but the limitations of applying the method to a larger range of systems, such as those with multiple degrees of freedom, are yet to be investigated.

## 7 Conclusions

This paper presents a novel method for calculating the nonlinear coherence for a nonlinear dynamical system in the presence of noise. A set of parameters are learnt to optimally combine an output prediction, calculated using an available model, with noisy measurements of the output to predict the input to the system. The nonlinear coherence is calculated using these parameters and is used as a metric of causality. Previously, the only way to estimate this was to model the system and assume the unmodelled component of $y_{n}$ to be noise. For highly nonlinear systems, this can be a poor estimation of the nonlinear coherence. The presented method improves upon this and is able to estimate the nonlinear coherence with excellent accuracy using a relatively small quantity of data. This would allow for users to collect experimental input and output data for a particular nonlinear dynamical system, and identify how much of the output is caused by the input across the frequency range, without developing expensive models. The method works for a broad class of nonlinear dynamical systems which are widely applicable across many scientific disciplines, from structural dynamics to neuroscience. The main limitation of the method is the interpretability and selection of the hyperparameter $\lambda$. Future work should focus on this and investigating the applicability of the method to a wider class of nonlinear dynamical systems.

#### Acknowledgements

This work was supported by the UK Engineering and Physical Sciences Research Council (EPSRC) Doctoral Training Partnership (DTP) PhD studentship. The authors would also like to thank Bose Corporation for their additional support and Professor Robin Langley for the many helpful discussions.

#### Reproduciblity Statement

The code and datasets required to reproduce the results in this paper are shared at https://github.com/obj22/repo-1. The hyperparameter settings for the simulations are provided in the appendix.

## Appendix A Controlled Inputs

If the inputs to a system can be controlled, the power spectrum density of the deterministic component of the response, $y$, can be calculated. Consider two noisy observations of $y(t)$, caused by the same input $x(t)$:

$$ $\displaystyle y_{n}(t)$ $\displaystyle=y(t)+\epsilon_{n}$ $\displaystyle y_{m}(t)$ $\displaystyle=y(t)+\epsilon_{m}\ .$ (27) $$

In the frequency domain, the power spectrum density of the deterministic component is defined to be

$$ $S_{yy}(f)=E\left[YY^{*}\right]$ (28) $$

Consider the cross-spectrum between $y_{n}$ and $y_{m}$:

$$ $\displaystyle S_{y_{n}y_{m}}(f)$ $\displaystyle=E\left[Y_{n}Y_{m}^{*}\right]=E\left[YY^{*}\right]+E\left[ \mathcal{E}_{n}Y_{m}^{*}\right]+E\left[Y_{n}\mathcal{E}_{m}^{*}\right]+E\left[ \mathcal{E}_{n}\mathcal{E}_{m}^{*}\right]$ (29) $\displaystyle=E\left[YY^{*}\right]\ .$ (30) $$

Therefore by taking the cross-spectrum between the two measured outputs, the deterministic component of the response can be identified. This is not applicable in this paper because repeat measurements are not available and $\mathcal{E}_{z}$, as defined in Section 3.2, is correlated with Y.

## Appendix B Architecture

### B.1 Derivation

The architecture is driven towards de-noising $\hat{Y}$ by optimally combining $Y_{n}$ and $Y_{z}$.

$$ $\displaystyle J$ $\displaystyle=E\left[(\hat{Y}-Y)(\hat{Y}-Y)^{*}\right]$ (31) $\displaystyle=E\left[(Y_{n}K+Y_{z}\left(1-K\right)-Y)(Y_{n}K+Y_{z}\left(1-K \right)-Y)^{*}\right]$ (32) $\displaystyle=K^{2}E\left[\mathcal{E}_{n}\mathcal{E}_{n}^{*}\right]+(1-K)^{2}E \left[\mathcal{E}_{z}\mathcal{E}_{z}^{*}\right].$ (33) $$

Then, differentiating with respect to $K$ gives

$$ $\frac{dJ}{dK}=2KE\left[\mathcal{E}_{n}\mathcal{E}_{n}^{*}\right]-2(1-K)E\left[ \mathcal{E}_{z}\mathcal{E}_{z}^{*}\right]=0\ .$ (34) $$

Therefore rearranging for $K$ gives

$$ $K=\frac{1}{1+\frac{E\left[\mathcal{E}_{n}\mathcal{E}_{n}^{*}\right]}{E\left[ \mathcal{E}_{z}\mathcal{E}_{z}^{*}\right]}}\ .$ (35) $$

## Appendix C Dataset

The dataset is split into training, validation and test sets. 10 frames were used for training, 10 frames were used for validation, and remaining frames were used for the test set. A batch size of 1 frame was used for training using the Adam optimiser. For each simulation, fourth order Runge-Kutta numerical integration was used in order to calculate the benchmark true $y$ using $x$ (Press et al., ).

### C.1 Polynomial Stiffness Case Study 1

The dataset contains 1000 frames of length 6000 samples. The parameters were set as: $\zeta=4.5$, $\alpha_{1}=5\times 10^{3}$, $\alpha_{2}=-10$ and $\alpha_{3}=3\times 10^{3}$. The input, $x$, was bandlimited noise with an rms of $\tau=1.0\times 10^{3}$. A step of 0.0025 was used in the numerical integration. Figure 7 shows the power spectrum density (PSD) for the system response and the linearised response.

Figure: Figure 7: Plot of the power spectrum density of $Y_{n}$ and $Y_{z}$ for the polynomial stiffness case study with no noise added.
Refer to caption: x7.png

### C.2 Saturating Stiffness Case Study

The dataset contains 1000 frames of length 6000 samples. The parameters were set as: $\zeta=3.0$, $\alpha_{1}=10^{4}$ and $\alpha_{2}=5\times 10^{3}$. The input, $x$, was bandlimited noise with an rms of $\tau=6.0\times 10^{3}$. A step of 0.005 was used in the numerical integration. Figure 8 shows the power spectrum density (PSD) for the system response and the linearised response.

Figure: Figure 8: Plot of the power spectrum density of $Y_{n}$ and $Y_{linear}$ for the nonlinear damping data with no noise added.
Refer to caption: x8.png

### C.3 Nonlinear Friction Case Study

The dataset contains 1000 frames of length 6000 samples. The parameters were set as: $\alpha_{1}=10^{5}$ and $\alpha_{2}=0.5$. The input, $x$, was bandlimited noise with an rms of $\tau=5$. A step of 0.002 was used in the numerical integration. Figure 9 shows the power spectrum density (PSD) for the system response and the linearised response.

Figure: Figure 9: Plot of the power spectrum density of $Y_{n}$ and $Y_{linear}$ for the nonlinear friction data with no noise added.
Refer to caption: x9.png

### C.4 Experiment Case Study

The dataset contains 180 frames of data of length 6000. A sampling frequency of 500 Hz was used to record the data. Figure 10 shows the power spectrum density (PSD) for the system response and the linearised response.

Figure: Figure 10: Plot of the power spectrum density of $Y_{n}$ and $Y_{linear}$ for the experimental data with no noise added.
Refer to caption: x10.png

## Appendix D λ 𝜆 \lambda italic_λ plots

Figures 11, 12, 13 and 14 show plots of how $\mathcal{L}^{V}_{x}$ varies as $\lambda$ increases. Plotted against the ratio of the two coefficients in the loss function, $\frac{\lambda}{1-\lambda}$, because it is the relative weighting between the two terms that is important. The red circle in each plot indicates the point at which $\mathcal{L}^{V}_{x}$ crosses the threshold of 0.01 above the minimum value. This is the value of $\lambda$ used to calculate the nonlinear coherence.

Figure: (a) Low Noise
Refer to caption: x11.png

Figure: (a) Low Noise
Refer to caption: x14.png

Figure: (a) Low Noise
Refer to caption: x17.png

Figure: (a) Low Noise
Refer to caption: x20.png

## Appendix E Computational Resources

The device used to run the code in this paper had the following specifications:
Processor: 12th Gen Intel(R) Core(TM) i7-12700KF 3.61 GHz
RAM: 64.0 GB
GPU: NVIDIA GeForce RTX 3090

## Peer Reviews

confidence: 2
contribution: 2
initial_rating: 3
presentation: 1
questions: 1. Why is Wavenet a good baseline? If one wants to distinguish what combining the forward model with a data inferred state brings you, can one not simply take the forward model from section 3.3?
2. Are the $\dot{x}$ and $\ddot{x}$ terms in Eq. 3 needed (they seem to be absent in all benchmarks)?
3. Eq. 7, is there an implicit Gaussian assumption?
4. Line 197 Why would $\lambda=0$ yield the optimal $K$, this is not obvious to me?
soundness: 3
strengths: - The method uses a sensible way to combine model/input based data with observed data.

- The resulting method outperforms the chosen baselines in terms of the coherence metric.
summary: The paper introduces a method for estimating the coherence $Co(y, y_n)$ between the observed noisy data $y_n$, and the “true” underlying system $y$. This estimate, termed $\hat{y}$ is obtained by combining, in frequency domain, a model/input based estimate $y_z$, with the noisy data $y_n$ (similar to a Kalman filter). The estimated $Co(y, y_n)$ is closer to $Co(y, y_n)$, than an only model/input-based estimate $Co(y_z, y_n)$.
weaknesses: 1. I found the paper hard to read, I think this can be improved by the authors more clearly explaining or motivating their motivation for the different steps (which might be obvious to them, but might not be to the reader), e.g.,  
    - 1.1 “Causality” is in the title and a key point of the paper, however I am totally missing some explicit motivation for why measuring coherence is the best way to estimate causality. The one citation for coherency (Silva et al., 2016), seems an odd choice.
    - 1.2 ANR introduced, but it, and the cited papers seem only tangentially related. This could be scrapped to make space for motivating the actual method, i.e., the use of coherency.
     - 1.3 No motivation is given for combining a model-based estimate with data in the frequency domain.
     - 1.4 Altough $\mathcal{F}$ of Eq. 1 is called a function, it is not a map from $x(t)$ to $y(t)$ (then there would be no dynamical system). Relatedly, it might also be better to use a different characters for $\mathcal{F}$ and $\mathcal{F}_\theta$, which are two very different things.

2. I find it very confusing that the paper introduces, and keeps referring to non-linear coherence, when the paper really seems to be after estimating the “standard”, linear coherence between $y$ and $y_n$. 

3. Page 9, has half its space dedicated to a Figure of an experimental setup, but this is just another benchmark, I would suggest to make this Figure smaller (or move it to appendix), and use the remaining space to improve readability.

4. The only evaluation here is $Co(\hat{y},y)$, against other coherences, which is not immediately intuitive. Could the authors not simply also plot and/or quantify how close $\hat{y}$ is to $y$?

5. Fig. 1 seems to lack a $\hat{Y}$ before, and $\hat{y}$ after the IFFT box. $X$ and a box with $\mathcal{F}$ (with some parameters $\theta$?) feeding into $y_z$ would also be helpful.

6. I find the motivation for needing $\mathcal{L}y$ and $\lambda$ in general a bit lacking in theoretical foundations. A possible perspective could potentially be that Eq. 5 gives one the variational estimate of the true distribution of $y$ as $q(\hat{y}|y_n, x)$,  $F_\theta$ gives one $p(x|\hat{y})$ as $\mathcal{N}(x; F_\theta(\hat{y}), \sigma_x^2)$. If one adds a prior on $\hat{y}$ as p$(\hat{y})=\mathcal{N}(y_n, \sigma_y^2)$ and writes out the typical ELBO loss on $p(x)$, both terms in Eq. 7 drop out (just the entropy of $q$ is missing) - although I am sure other perspectives are possible!

---

confidence: 2
contribution: 2
initial_rating: 3
presentation: 2
questions: - What is the relationship between K and \lambda? Is K a hyper parameter during network training? How did you find the optimal K. 
- Why the modeling error / measurement error affects the signal in an additive way (Eq.4)
soundness: 2
strengths: - The presented method learns the a nonlinear functional mapping through a neural network. It works under a wide range of nonlinearities form as long as y can be mapped back to x (Eq.3).  
- It presented a novel way to balance modeling / observation error and contributed to maximize signal reconstruction performance in highly nonlinear systems. 
- Even in high noise cases, the recovered coherence was still high.
summary: The manuscript presents a novel technique for separate measurement noise and modeling errors without accessing the actual input-output relationship. The calculation was performed in frequency domain and there's one weighting parameter per frequency controlling the contribution of modeling error and measurement error to achieve minimum reconstruction error. The performance, measured by the coherence between actual signal and recovered signal is high.
weaknesses: - Only applies to the specific ODE described in Eq.3. Maybe the authors could extend it further as long as there's a one-to-one mapping from y to x? 
- the method only works from 1D signal to 1D signal? Is there any possibility to extend to high dimensional nonlinear systems.

---

confidence: 2
contribution: 2
initial_rating: 3
presentation: 1
questions: Comments:

- The description of the problem and approach in the abstract is not clear enough.
- Could you provide a more thorough statistical analysis of the performance of the approach for the simulated datasets over different parameter regimes of the underlying model, to show its robustness?
- The figures are missing captions.
- It would be informative, for the simulated and experimental results, to explicitly show how a better approximation of the response (as presented in the figures) helps in the interpretation of the original systems.

Minor comments:
- ‘but has recently begun’ - the ‘but’ part of the sentence should probably be rephrased
- it should be mentioned in the abstract that this work is restricted to a specific class of ODEs.
- ‘A completely new strategy’ - rephrase?
- The related work section sounds too defensive.. e.g. ‘The presented method .. is very different to work in the literature.’
- ‘In this paper, the primary aim is to estimate the noise level and the model error, but a similar principle is utilised.’ - unclear.
- ‘The presented method.. provides an excellent estimation of the nonlinear coherence for each noise case’ - rephrase in a more quantitative way?
soundness: 2
strengths: - The approach is well characterized mathematically in the manuscript.
- The demonstration of the approach over different types of simulated and experimental data is interesting.
- There are comparisons to alternative approaches over the same task.
summary: The authors of the manuscript ‘A METHOD FOR IDENTIFYING CAUSALITY IN THE RESPONSE OF NONLINEAR DYNAMICAL SYSTEMS’ develop a method for the calculation of the nonlinear coherence between input and output data for a noisy ODE class given random input. The approach allows to learn an optimal combination of an output prediction with noisy measurements. The nonlinear coherence that is calculated is then used as a measure of causality. The method is shown to perform well on both simulated dynamical datasets as well as experimental data.
weaknesses: - The manuscript could be improved in terms of readability. It is in many cases very difficult to follow.
- While the method was evaluated over real experimental data, the noise was added synthetically.

Meta-review: The authors aim in this work to provide a measure by which to judge the best possible dynamical system model that can fit experimentally measured inputs and outputs of a system. The measure focuses on a frequency based approximation of coherence,  and applies the metric to a number of systems. 

There were a number of issues raised, primarily in the clarity and impact of the work. Multiple reviewers found the work difficult to parse, which clearly led to a number of confusions including in notation and in properly understanding the underlying assumptions the govern the work. Despite extended conversations with the reviewers, it seems that as it stands the assumptions tend to be limiting and the manuscript can be much more clearly honed to describe how the metric presented connects back to model goodness and interpretability. These seem to be significant drawbacks and a revised manuscript might fare better.

Decision: Reject