Tensor Data Imputation by PARAFAC with Updated Chaotic Biases by Adam Optimizer

The big data pattern analysis suffers from incorrect responses due to missing data entries in the real world. Data collected for digital movie platforms like Netflix and intelligent transportation systems is Spatio-temporal data. Extracting the latent and explicit features from this data is a challenge. We present the high dimensional data imputation problem as a higher-order tensor decomposition. The regularized and biased PARAFAC decomposition is proposed to generate the missing data entries. The biases are created and updated by a chaotic exponential factor in Adam's optimization, which reduces the imputation error. This chaotic perturbed exponentially update in the learning rate replaces the fixed learning rate in the bias update by Adam optimization. The idea has experimented with Netflix and traffic datasets from Guangzhou, China.

INTRODUCTION
The low quality and homogenous data require the verification and validation process, which is more timeconsuming. Data mining approaches are used to manage a large amount of data in a particular context. Data mining schemes increases the quality of data with a large amount of data set. The meaningful pattern and rules are discovering in a large amount of data by exploration and analyzing them. In data mining approaches, various numbers of models and methods are considered to filter large data. The data mining methods are divided into different categories based on objective evaluation. The main categories are classification, estimation, prediction, clustering, association, and profiling. As we studied in the literature [1][2][3][4][5][6][7][8][9][10][11], various data mining models are proposed for Big Data analytics. The missing data problem is a setback for big data analytics. In the modern era, an intelligent transportation system (ITS) is used to collect large traffic data. Various mobile and stationary sensors are placed in the ITS configuration. The correct information collection in the transportation system represents the state of the system for the development of correct analysis. The collected data is not directly connected to the transportation system due to the presence of traditional data (traffic volume) and a new kind of data (Bluetooth and GPS data). The missing data problem has commonly arrived in a large traffic dataset. The missing data can degrade the performance of the intelligent transportation system. To improve the reliability of the traffic data collection and record the data mining approaches are used in state of the art. Manuscript  Previously various data imputation methods are proposed to solve the missing traffic data problem. The traffic data imputation problems are categorized into two sections: traffic volume imputation and traffic speed data estimation. Both traffic volume and traffic speed data have similar characteristics, like temporal and spatial correlation. The high dimensional data is converted into the tensor as it inherits the high algebraic information than the matrix. The purpose of the selection of traffic and Netflix data is, these contain the temporal as well as spatial resolutions. The pattern is learned by tensor matrix factorization for a part of the data, and that learning is used to impute the missing entries. This paper aims to propose the advanced tensor factorization for data imputation.

Contribution
Our contribution to the work is two-fold: regularized tensor factorization and Adam optimization with dynamic learning rate. The data used in our work is a tensor matrix that is a multidimensional matrix. Several tensor factorization schemes are available like CANDECOMP/PARFAC, tucker, SVD, etc. [20,[22][23][24][25][26][27][28], but the PARAFAC method provides more stable tensor decomposition with faster convergence. The tensor is decomposed to rank one matrix. Few researchers used tensor decomposition for the high dimensional data imputations too [29] [30]. The PARAFAC decomposition still doesn't deal with the unexpected differences in the data values, which are counted as noise in the data. We propose here modified regularized PARAFAC, which decomposes the tensor into three factors [20] as well as into bias matrices for those tensor factors which avoid the fall into local optima. We used the dynamic learning rate to update the bias values during Adam's optimization. The fixed learning rate of 0.001 is changed to a logistic mapped chaotic exponential decaying learning rate with every iteration.

Paper Organization
Further, in this article, we have discussed previous related work in section 2. Section 3 discusses the problem in data imputation and PARAFC tensor factorization. The proposed extended PARAFC tensor factorization for imputation has been discussed in section 4, which is followed by simulation results in the next section. We conclude the results in section 5.

II. LITERATURE REVIEW
State of the art schemes is utilized for the data mining process to achieve processed data. A high-quality data mining approach was presented in [1] by using the K-nearest neighbor algorithm. A predictive data mining framework was presented with the help of various up to date algorithms schemes like machine learning [2].
The machine learning algorithms improves the computational efficiency in removal of unwanted noise from the massive dataset. The Ghost algorithm was used to reconstruct the missing data from a large dataset that recover off period segments of missing data [3]. The ghost algorithm searched in a sequential manner dataset to achieve the data segment. A caching approach was also used to minimize the search space and improves the computational complexity to linear. In [4], various deep generative models were compared in terms of data mining. A Graph Neural Network (GNN) model was enhanced for the missing data imputation (MDI) application. With the help of GNN, a Graph convolutional autoencoder was designed to reconstruct the complete dataset [5]. The adversarial loss and global information were included in the dataset during the reconstruction phase. A Distributed Neural Network (DNN) was designed for the imputation of missing value in terms of Big data context [6]. The outcomes of the DNN data imputation method were compared with the K-nearest neighbor approach. They provided a better response. A selfdeveloped competitive neural network Adaptive Response Theory 2 (ART2) was proposed for the data imputation [7]. The ensemble approach provided accurate results on intra cluster non-missing value elements. The ART2 approach also improved imputation accuracy. The Big Data K-mean and Big Data Fuzzy imputation cluster-based methods were presented in [8]. Similarly, in [9], two algorithms were used to remove the noisy data from Big Data set. The two methods were Heterogeneous and Homogenous ensemble filters. These algorithms were formed with a combination of Deep learning classification models. An MLClean approach was presented for the Big Data preprocessing to trained the Deep Learning model [10]. It provided data cleaning, unfair mitigation, and data sanitization. The MLClean approach was commonly used in Big Data analytics cases. The ROSEFW-RF algorithm was developed based on the Mapreduce function for imbalanced data mining [11]. The medical field data is imbalanced due to its number of parameters. To achieve balanced data, the relevant features, classes, category, and values are classified with the proposed ROSEFW-RF approach. In a study [12], the author provided the Bayesian CP factorization to handle the noisy tensor data. A non-deterministic model was developed using Bayesian CP factorization, which estimated the rank of CP automatically. A robust Bayesian generative model developed in [13] for tensor factorization. It provided the robustness to the boundary conditions and reduced the Gaussian noise, which minimizes the overfitting problem. The low-rank decomposition using a Bayesian model for multiway tensors data with missing observations was presented in [14]. The single and coupled tensor factorization models were tested by the full Bayesian Inference using VB [15]. The full Bayesian Inference method was easily implemented to any large scale model. A fully Bayesian Probabilistic Matrix Factorization (PMF) model was presented in [16]. The overall model parameters and hyperparameters controlled the model capacity. A factor-based algorithm was presented for the classification application known as tensor factorization. The fully Bayesian Probabilistic Tensor Factorization (PTF) was used to analyze the large scale dataset [17]. Some real-world problems of classification were analyzed by the Bayesian PTF model with accurate outcomes. A Bayesian Gaussian CANDECOMP/PARAFAC (BGCP) decomposition model was used to impute the multidimensional imbalance traffic data [18]. In tensor representation with BGCP, the accurate outcomes were found in the third-order case for both random missing and fiber missing traffic data scenarios. The multiple metrics and tensors combined joints were presented by the Bayesian Multi-tensor Factorization model [19]. The Bayesian MTF model analyzed the total number of data collection factors. A multidimensional EEG dataset was analyzed by the Bayesian Tensor Factorization model [20]. The noninformative EEG signals were removed from the BTF method, which improved the classification response in the medical research field. A fully Bayesian system was presented for the automatically learning parameters of missing values traffic data model with Variational Bayes (VB) [21]. The data was collected from the urban traffic speed data set collected in Guangzhou, China. The two methods were tested on the traffic data imputation. The Bayesian Augmented Tensor Factorization (BATF) provided better accuracy in terms of missing values imputations. A three procedure framework was proposed for the missing traffic data recovery with traffic patterns recognition. The proposed STD approach provided higher accuracy in terms of missing traffic data [22]. A matrix factorization technique based on a K-nearest neighbor method was proposed for the Netflix missing data recovery [23]. The matrix factorization based KNN model provided higher accuracy, and it can be easily implemented on realworld datasets. The collaborative filtering approaches were applied for the data filter in [24,25]. The unique properties of implicit feedback datasets were recognized by a factor model [25] III.

PROBLEM STATEMENT
Due to vast applications of tensor, its decomposition for multidimensional data has been the researcher's interest. Their work can be broadly classified into two categories: decomposition by CANDECOMP/PARAFAC (Canonical Decomposition/ Parallel Factor), also referred to as CP and tucker factorization. Both are higher-order principal component analysis (PCA). The convex hull analysis in [36] proves the CP decomposition is better than tucker decomposition on the criteria of good fit analysis. The PARAFAC itself takes a long time in reaching the local minima solution and also has the probability to stuck in the local minima. To avoid it with local optima, we regularized the PARAFAC Tikhonov regularization factor. Still, this form of PARAFAC only considers the latent features, but the data is Spatio-temporal. So, we need to add the biases for all the three dimensions of the data, which can extract the temporal features too. Besides the decomposition part, the rank selection is also uncertain in the tensor decomposition. The Bayesian approach has also been suggested by a few researchers for rank selection [12][13][14][15][16][17][18][19][20][21]. But these focused either on prediction problems or factorization with missing data. Whereas we couldn't find it suitable for data imputation problem as Bayesian optimization works preferably for automatic Rank detection for higher-order low-rank factorization. In the imputation application of CP, the normalized imputation error increases with the increase in rank. We have tested the scheme of regularized augmented PARAFACfor various ranks from 1-10. A plot in figure 1 demonstrates that rank 1 is suitable for data imputation scheme and imputed data for this rank by PARAFAC, and regularized augmented PARAFAC is shown in figure 1(b). It supports our convention that PARAFAC lacks in providing a global solution for missing data generation. The matrices ∈ ℝ × , ∈ ℝ × , ∈ ℝ × , and ∈ ℝ × are generated by the random entries drawn from the standard normal distribution. A new matrix ∈ ℝ × × is generated by a random setting of % of the entries of to be missing. So = * , here binary tensor is represented by the which has an equal dimension of tensor and , , = 0 for the set , , to missing. Figure 2 shows the tensor with random missing entries.

Fig. 2 Arrangement of incomplete tensor by random
indexing.

IV. PROPOSED TENSOR FACTORIZATION A. Regularized-PARAFAC
As discussed above, we factorize the tensor matrix by the PARAFAC method. For the data imputation in our work, the data is of third-order tensor ℝ × × and as for large rank ( > 10), the PARFAC decomposition can be shown as in equation 1 [29] [33].
≈ , , Equation 1 holds for the sufficiently large rank , but since in our case, we are interested in the decomposition up to = 1 [13]. The tensor matrix can be decomposed into N factors of rank 1 as: Where ∘ represents the cross product. For the third-order tensor, it is ∘ ∘ . The conviction of decomposition of equation 1 to the first-order tensor is: The || . || 2 Represents the Frobenius norm. This one is the non-convex optimization problem and can be solved iteratively as The equation 4 has an issue of local optimum solution due to non-uniform scaling of factors [33], so Paatero [28] suggested the addition of a Tikhonov regularization factor in equation 4 which resolved the issue of many local solutions in equation 3 and updated it as ℒ ; , , = Where is a regulation parameter with > 1. In our work, we are using the data set from China Traffic Data [13]. This kind of tensor data has three factors for road segments, the number of time slots, and recorded time intervals. Equation 5 maps the combined rating for users and content in the latent space, but it is not always the case. The temporal dependency on the data is not examined in that, due to which the imputation might be less correct. This biasing for user and content must be included in equation 5 [29]. The augmented equation 5 is written in equation 6 2 ) (7) Where ∅ , , are the bias factors for the three tensor dimensions, is the average of all three factors [31]. The biases are the deviation of each factor from the average . The proposed tensor is shown in figure 3. The factor matrix of the tensor can also be represented as , , = + ∅ + + + [ , , ].

Fig. 3 Proposed augmented PARAFAC factorization for 3 rd order tensor
The proposed modified-PARAFAC is an iterative process, and equation 7 is minimized by gradient descent optimization [13], Alternating least square method [31] recently. But these methods are easily stuck in local minima and less rate of convergence. The advancement in deep learning is dependent upon the fast and global optimized Adam optimization [35]. We factorize the tensor by Adam optimization. The bias is also updated iteratively in each stochastic step. These are initialized at some random value between 0 and 0.1. Though the initial feeding point for these biases can be random, however for our application, we set 0,0.1 , where represents the biases. Biases are updated in each iteration of Adam's optimization. Equation  6 is the objective function of adam optimization, and the error thus calculated is used to update the bias values. We multiply the error from equation 6 by an exponentially decreasing regularization parameter, which reduces with iteration to minimize the error. To add the chaotic perturbation in error, we multiply it with a less random factor, which is calculated by a logistic map. The biases are updated as: Here − , , is the deviation in the factors, is a constant whose values if fixed to 20 in our case. and represent the current iteration and maximum iteration in the adam optimization.
is the perturbed parameter by logistic mapping which is calculated as = 1 − , = 0,1,2,3 ….. The is a system parameter 0,4 0,1 . This way and are also updated. The pseudo-code for the update of biases is shown in algorithm1. The Netflix dataset provides a testset containing 2,817,131 movie pairs with ratings. The latest movie pairs ratings are selected from a recent subset of users in the training dataset.

C. Evaluation
We have followed the random missing entries scheme for missing values generation. After creating the missing values for 10%, 20%, and 50%, we use rank 1 from figure 1 for the tensor data completion. The proposed decomposition solution with bias update is compared with the state of the art schemes. We compared the imputation with Tucker decomposition [30], Bayesian inference [29] for the random missing data. The Netflix dataset is also checked for only random missing entries. The proposed scheme is also evaluated for gradient descent (GD) and Adam optimization for the optimal solution of equation 6. The evaluation of results is statistically done based on mean absolute error (MAE), mean absolute percentage error (MAPE), and root means square error (RMSE). If the and are estimated and true missing values, then these performance measures are defined as: Netflix dataset and traffic datasets are evaluated by the proposed scheme and compared with states of the art SVD combined tucker decomposition [30], BATF [29], BGCP [20], BCPF [21], and modified PARAFAC variants. The schemes in [20], [21], [29] have common Bayesian Gaussian tensor factorization. The CANDECOMP/PARAFAC factors are trained by Bayesian optimization. The [21] presented the hyperparameter of automatic rank determination. In contrast, the solution doesn't differ significantly from other variants in [20] and [29] as the researchers only used the higher number of hyperparameters for optimization. However, the Bayesian optimization of factorization parameters performed well over Tucker-SVD decomposition [30]. Table 1 shows the performance comparison of the proposed data imputation scheme for random missing entries for traffic datasets. It can be observed from table 1 that performance is somewhat good at lower missing ratio; however, the change is not significant with increasing missing values. The dynamic, chaotic learning rate and modified regularized tensor decomposition lead to superior performance than state of the art methods. Figure 4 represents the imputed tensor for different random missing ratios for the traffic dataset by the proposed scheme. The proposed scheme doesn't get affected by the change in the missing rate like other schemes.  Fig. 4 The imputed tensor with random missing entries percentage (a) 50% (b) 30 % and (c) 10% The bias to generate the missing values is updated by gradient descent (GD) optimization in [30]; however, as discussed in the previous section, the Adam optimizer is better than GD. We have followed the update as in algorithm 2. The SVD-Tucker decomposition proposed the constant learning rate to update the bias. Still, it has been proven recently that for the sparse dataset, the adaptive learning rate is better than a constant learning rate. The learning rate can be adaptive dependent upon time, step size, or exponentially changing with iteration. Other Deep learning toolboxes like Keras also provide the feature to use the adaptive learning rate during training. We have perturbed the exponentially decaying learning rate by logistic mapping. It is a linear mapping with variable = 1 − , = 0,1,2,3 ….. The is a system parameter 0,4 0,1 . The logistic mapping shows different behavior for different values of . Figure 5 shows the stationary, periodic, and complete bifurcation diagram for the logistic map. Figure 5(a) and 5(b) shows the behavior of logistic mapping if the system parameter is less than 1 and 3 ≤ < 1 + 6 respectively. For 0 ≤ < 1, updated biases with this logistic map value can lead to premature convergence of the training. Conclusively, data imputation won't be similar to extracted latent features of the rest of the data. The non-diminishing oscillatory behavior ( figure 5(b)) also doesn't converge at all. The chaotic behavior starts beyond 3.56994, and figure 5(c) shows the logistic mapping for a complete range of ∈ 0,4 . The area between ∈ 3.54409,4 is the stable oscillations area, and convergence in optimization can be achieved in this area.
The solid lines in figure 5(c) point to the stable solution. The more uniform is the bias, the less deviation is observed in variance for a dimension. Just to recall, the three biases are for the three dimensions of the Spatio-temporal data. The size of a bias matrix is equal to the size of the decomposed factor matrices. The first row in figure 6 shows the uniformity comparison for the bias matrices for the final selected solution for a 50% missing ratio. We have used the ternary plot since we have 3-dimensional tensor. The final updated bias for the minimum gradient in equation 7 is recorded for plotting. The vertices of the ternary plot represent the updated bias by Adam optimized modified regulated PARAFAC with a dynamic learning rate, with a fixed learning rate and proposed PARAFAC segmentation by GD optimization. These three plots are for each dimension of the data. The dots inside the plots are the bias values that arrange themselves in an elliptical shape with alignment towards the proposed scheme. The elliptical shape is formed due to uniformity in the bias values; the least is the standard deviation; more is the uniformity. The alignment towards proposed scheme vertex indicates the clustering of biases is greatly affected by larger bias magnitudes generated by the proposed scheme. However, the uniformity of bias < ∅ & , as indicated in figure  6(b) for bias for the number of days. The second row in figure 6 has the plots for imputed values and original entries. We have plotted these curves for a few samples. The sparse tensor's entries are also plotted in 6(d) whose entries are all zeros for these samples. (b) for = .  For the Netflix dataset, we have compared the proposed scheme with its three other variants only due to unavailability of state of the art schemes results on this data. The figure shows the performance plot comparison for the Netflix dataset. This dataset has primarily tested for the Netflix prize competition for recommending the most suitable movies to the users. The modified PARAFAC with the proposed dynamic learning rate has minimum MAPE and RMSE. The MAPE has shown no significant change with the increase in the missing entries, but with 80% missing data, all variants have degraded. However, no such behavior is recorded for MAE and RMSE curves. The proposed scheme has the least MAPE and RMSE values at various missing ratios. This dataset also strengthens the theory that Adam optimization outperforms the PARAFAC with gradient descent optimization. The bias for each dimension of the Netflix data is shown in the second row of figure 7.

VI. CONCLUSION
The data imputation is the first step in the data denoising and to process the data with deep learning. In this research article, we have profoundly compared the performance of imputed tensor data with state of the art schemes. The Spatio-temporal real-world tensor data can be decomposed into the low-rank matrix for imputation using two kinds of decomposition primarily: CD/PARAFAC and Tucker decomposition. We compared the proposed modified PARAFAC with a dynamic learning rate with the latest work of Tucker-SVD decomposition on the Traffic dataset. It has been observed that our method is able to reduce the mean absolute error by 20% for the 50% missing entries in the data. Few other researchers have focused on the Bayesian learning training for the tensor decomposition. It has also been suggested previously that automatic rank determination leads to a better generation of missing entries. Our dynamic Adam optimization scheme has powered the modified PARAFAC with a 0.48% reduced root mean square error for 50 % missing data. We have tested the algorithm on Netflix movies dataset too. The size of the dataset has backed the performance due to our hardware constraint, but it has been noticed that dynamically learned modified PARAFAC had won the competition. In the next part of our research, we will test the performance on Hadoop and Spark as the large size of the data created the havoc in the testing. We will also cluster the imputed tensor and semantically compared the imputation.


Positional Encoder Graph Neural Networks for Geographic Data

Graph neural networks (GNNs) provide a powerful and scalable solution for modeling continuous spatial data. However, in the absence of further context on the geometric structure of the data, they often rely on Euclidean distances to construct the input graphs. This assumption can be improbable in many real-world settings, where the spatial structure is more complex and explicitly non-Euclidean (e.g., road networks). In this paper, we propose PE-GNN, a new framework that incorporates spatial context and correlation explicitly into the models. Building on recent advances in geospatial auxiliary task learning and semantic spatial embeddings, our proposed method (1) learns a context-aware vector encoding of the geographic coordinates and (2) predicts spatial autocorrelation in the data in parallel with the main task. On spatial regression tasks, we show the effectiveness of our approach, improving performance over different state-of-the-art GNN approaches. We also test our approach for spatial interpolation, i.e., spatial regression without node features, a task that GNNs are currently not competitive at. We observe that our approach not only vastly improves over the GNN baselines, but can match Gaussian processes, the most commonly utilized method for spatial interpolation problems. The code for this study can be accessed via: https://github.com/konstantinklemmer/pe-gnn

INTRODUCTION
Geographic data describes all data that can be mapped onto the earth, for example using latitude and longitude coordinates. As such, geographic data is characterized by a natural geometric structure, which often strongly affects the observed spatial pattern. While traditional neural network approaches do not have an intuition to account for spatial dynamics, graph neural networks (GNNs) can represent spatial structures graphically. The recent years have seen many applications leveraging GNNs for modeling tasks in the geographic domain, such as inferring properties of a point-ofinterest [35] or predicting the speed of traffic at a certain location [5].
Nonetheless, as we show in this study, GNNs are not necessarily sufficient for modeling complex spatial effects: spatial context can be different at each location, which may be reflected in the relationship with its spatial neighborhood. The study of spatial context and dependencies has attracted increasing attention in the machine learning community, with studies on spatial context embeddings [24,34] and spatially explicit auxiliary task learning [22].
In this study, we seek to merge these streams of research. We propose the positional encoder graph neural network (PE-GNN), a flexible approach for better encoding spatial context into GNN-based predictive models. PE-GNN is highly modular and can work with any GNN backbone. PE-GNN contains a positional encoder (PE) [24,30], which learns a contextual embedding for point coordinates throughout training. The embedding returned by PE is then concatenated with other node features to provide the training data for the GNN operator. PE-GNN further predicts the Moran's I metric of spatial autocorrelation of the output as an auxiliary task in parallel to the main objective, expanding the approach proposed by [22] to continuous spatial coordinates. Lastly, we train PE-GNN by constructing a novel training graph, based on -nearest-neighborhood, from a randomly sampled batch of points at each training step. This forces PE to learn generalizable features, as the same point coordinate might have different spatial context (neighbors) at different training steps. Similarly, this training approach leads us to compute a "shuffled" Moran's I, implicitly nudging the model to learn a general representation of spatial autocorrelation which works across varying neighbor sets. Over a range of spatial regression tasks, we show that PE-GNN improves performance of different GNN backbones in most settings, except when very little training data is available.
The contributions of our study are as follows: • We propose PE-GNN, a novel GNN architecture including a positional encoder learning spatial context embeddings for each point coordinate to improve predictions. • We propose a novel way of training the positional encoder, distinct from [24]. We use the output of PE concatenated with other node features to predict an outcome variable. PE learns through backpropagation on the main regression loss in an end-to-end fashion. Training PE thus takes into account not only the eventual variable of interest, but also further contextual information at the current location-and its relation to other points. • We expand the Moran's I auxiliary task learning framework proposed by [22] for continuous spatial data. • Our training strategy involves the creation of a new training graph at each training step from the current, random point batch. This enables learning of a more generalizable PE embedding and allows computation of a "shuffled" Moran's I, which accounts for different neighbors at different training steps, thus overcoming the well-known scale sensitivity of Moran's I.
• To the best of our knowledge, PE-GNN is the first GNN based approach that is competitive with Gaussian Processes on pure spatial interpolation tasks, i.e., predicting a (continuous) output based solely on spatial coordinates, as well as substantially improving GNN performance on both spatial regression and spatial interpolation tasks.
2 RELATED WORK 2.1 Traditional and neural-network-based spatial regression modeling Our work considers the problem of modeling geospatial data. This poses a distinct challenge, as standard regression models (such as OLS) fail to address the spatial nature of the data, which can result in spatially correlated residuals. To address this, spatial lag models [3] add a spatial lag term to the regression equation that is proportional to the dependent variable values of nearby observations, assigned by a weight matrix. Likewise, kernel regression takes a weighted average of nearby points when predicting the dependent variable. Nevertheless, the most popular off-the-shelf methods for modeling continuous spatial data are based on Gaussian processes [6].
With the emergence of machine learning, there has been a rise of research on applications of neural network models for spatial modeling tasks. More specifically, graph neural networks (GNNs) are often used for these tasks with the spatial data represented graphically. Particularly, they offer flexibility and scalability advantages over traditional spatial modeling approaches. Specific GNN architectures including Graph Convolutional Networks [20], Graph Attention Networks [31] and GraphSAGE [12] are powerful methods for inference and representation learning with spatial data.
Recently, GNN approaches tailored to the specific complexities of geospatial data have been developed. The authors of Kriging Convolutional Networks [4] propose using GNNs to perform a modified kriging task. [12] apply GNNs for a spatio-temporal Kriging task, recovering data from unsampled nodes on an input graph.
We look to extend this line of research by providing stronger, explicit capacities for GNNs to learn spatial structures. Additionally, our proposed method is highly modular and can be combined with any GNN backbone.

Spatial context embeddings for geographic data
Through many decades of research on spatial patterns, a myriad of measures, metrics, and statistics have been developed to cover a broad range of spatial interactions. All of these measures seek to transform spatial locations, with optional associated features, into some meaningful embedding, for example, a theoretical distribution of the locations or a measure of spatial association. A set of point locations in 2 -space may be distributed randomly, or may follow some underlying spatial pattern. Spatial point processes are some of the oldest methods of analyzing spatial data and are a useful tool when point locations depend on one another, e.g., in the case of epidemiological or ecological data. Point processes are particularly popular for modeling geographic data, for example Cox processes [1] and Hawkes processes [29]. Empirical observations can be embedded using these point processes, for example, to make predictions on future spatial spread. If we now assume that point locations also come with a continuous feature value, the range of potential spatial effects and associated metrics is vastly expanded.
The most common metric for continuous geographic data is the Moran's I statistic, developed by [2]. Moran's I measures local and global spatial autocorrelation and acts as a detector of spatial clusters and outliers. The metric has also motivated several methodological expansions, like local spatial heteroskedasticity [26] and local spatial dispersion [33]. Measures of spatial autocorrelation have already been shown to be useful for improving neural network models through auxiliary task learning [22], model selection [21], embedding losses [23] and localized representation learning [10].
Beyond these traditional metrics, recent years have seen the emergence of neural network based embeddings for geographic information. [32] use kernel embeddings to learn social media user locations. [10] devise an approach using local point-of-interest (POI) information to learn region embeddings and integrate similarities between neighboring regions to learn mobile check-ins. [34] develop GPS2Vec, an embedding approach for latitude-longitude coordinates, based on a grid cell encoding and spatial context (e.g., tweets and images). [24] developed Space2Vec, another latitudelongitude embedding without requiring further context like tweets or POIs. Space2Vec transforms the input coordinates using sinusoidal functions and then reprojects them into a desired output space using linear layers. All of these studies highlight the potential of spatial coordinate embeddings (with or without spatial context enrichment) for improving predictive models for geographic data.

METHOD 3.1 Graph Neural Networks with Geographic Data
We elaborate on our method using Graph Convolutional Networks (GCNs) as an example backbone for our novel PE-GNN approach. Let us first define a datapoint = { , x i , c i }, where is a continuous target variable (scalar), x i is a vector of predictive features and c i is a vector of point coordinates, mapping the datapoint into geographic space (e.g., latitude and longitude values). We use the great-circle distance = ℎ (c i , c j ) between point coordinates to create a graph of all points = { 1 , . . . , } in the set. This can be done either by using a -nearest-neighbor approach or by drawing a boundary radius around each coordinate to define each point's neighborhood and create an adjacency matrix A. The graph = ( , ) consists of a set of vertices (or nodes) = { 1 , . . . , } and a set of edges = { 1 , . . . , } as assigned by A. Each vertex ∈ has respective node features x i and target variable . While the adjacency matrix A usually comes as a binary matrix (with values of 1 indicating adjacency and values of 0 otherwise), one can account for different distances between nodes and use point distances or kernel transformations thereof [4] to weight A. Now, let D define the degree matrix and I the identity matrix of our graph . We can then define the normalized adjacency matrixĀ as: As proposed by [20], a GCN layer can now be defined as where describes an activation function (e.g., ReLU) and W ( ) is a weight matrix parametrizing GCN layer . The input for the first GCN layer H (0) is given by the feature matrix X containing all node feature vectors x 1 , . . . , x n . The assembled GCN predicts the output Y = (X, Θ ) parametrized by Θ .

Context-aware spatial coordinate embeddings
In the traditional GCN setting, the only intuition for spatial context stems from connections between nodes which allow for graph convolutions, akin to pixel convolutions with image data. This restricts the capacity of the GCN to capture spatial patterns in several ways: • While defining good neighborhood structures can be crucial for GCN performance, this often comes down to somewhat arbitrary choices like selecting the nearest neighbors of each node. Without prior knowledge on the underlying data, the process of setting the right neighborhood parameters can be difficult and require extensive testing. Furthermore, a single value of might not be best for all nodes: different locations might be more or less dependent on their neighbors. • Assuming that no underlying graph connecting point locations is known, one would typically construct a graph using the distance (Euclidean or other) between pairs of points. In many real world settings (e.g., points-of-interest along a road network) this assumption is unrealistic and may lead to poorly defined neighborhoods. • GCNs contain no intrinsic tool to transform point coordinates into a different (latent) space that might be more informative for representing the spatial structure, with respect to the particular problem the GCN is trying to solve.
As such, GCNs struggle with tasks that explicitly require learning of complex spatial dependencies. For example, we show in our experiments that simple GCNs are not able to solve simple spatial interpolation tasks, i.e., predicting a continuous outcome variable from the point coordinates only. We propose a novel approach to overcome these difficulties, by devising a new positional encoder module, and learning a flexible spatial context encoding for each geographic coordinate, motivated by recent advances in transformers and spatial representation learning.
Given a batch of datapoints, we create the spatial coordinate matrix C from individual point coordinates c 1 , . . . , c n and define a positional encoder and a fully-connected neural network (Θ ), parametrized by Θ . Following the intuition of transformers [30] for geographic coordinates [24], the sinusoidal transform is a concatenation of scale-sensitive sinusoidal functions at different frequencies, so that with being the total number of grid scales and and setting the minimum and maximum grid scale (comparable to the lengthscale parameter of a kernel). The scale-specific encoder )] processes the spatial dimensions (e.g., latitude and longitude) of C separately, so that ; sin where = . The output from is then fed through the fully connected neural network (Θ ) to transform it into the desired vector space shape, creating the coordinate embedding matrix C = (C, , , Θ ).

Auxiliary learning of spatial autocorrelation
A further particularity of geographic data is that it often exhibits spatial autocorrelation: observations are related, in some shape or form, to their geographic neighbors. Measures of spatial autocorrelation have recently been successfully integrated into neural networks for discrete spatial data (images) [22,23], however no such approach exists for continuous spatial coordinates.
Spatial autocorrelation can be measured using the Moran's I metric of local spatial autocorrelation [2]. Moran's I captures localized homogeneity and outliers and thus functions as a detector of both spatial clustering and spatial change patterns. In the context of our problem, the Moran's I measure of spatial autocorrelation for our outcome variable is defined as: where , ∈ A denotes adjacency of observations and . As proposed by [22], predicting the Moran's I metric of the output can be used as auxiliary task during training. Auxiliary task learning [28] is a special case of multi-task learning, where one learning algorithm tackles two or more tasks at once. In auxiliary task learning, we are only interested in the predictions of one task; however, adding additional, auxiliary tasks to the learner might improve performance on the primary problem: the auxiliary task can add context to the learning problem that can help solve the main problem. This approach is commonly used, for example in reinforcement learning [9] or computer vision [14,15].
Translated to our GCN setting, we seek to predict the outcome Y and its local Moran's I metric (Y) using the same network, so that [Ŷ,( Y)] = (X). As [22] note, the local Moran's I metric is scale-sensitive and, due to its restriction to local neighborhoods, can miss out on longer-distance spatial effects [7,25]. But while [22] propose to compute the Moran's I at different resolutions, the GCN setting allows for a different, novel approach to overcome this issue.
Training of GCNs is often conducted using minibatches of points at each training step. Rather than constructing the graph of training points a priori, we opt for a procedure where in each training step, ℎ points are sampled from the training data as batch . A graph with corresponding adjacency matrix A is constructed for the batch, along with the Moran's I metric of the outcome variable (Y ). This approach brings a unique advantage: When training with (randomly shuffled) batches, points may have different neighbors in different training iterations. The Moran's I for point can thus change throughout iterations, reflecting a differing set of more distant or closer neighbors. This also naturally helps to overcome Moran's I scale sensitivity without the need to compute the metric at different resolutions. Altogether, we refer to this altered Moran's I as "shuffled Moran's I".

Positional Encoder Graph Neural Network (PE-GNN)
We now assemble the different modules of our method and introduce the Positional Encoder Graph Neural Network (PE-GNN). The modular approach allows for a seamless integration of the positional encoder and the spatial autocorrelation auxiliary task with any graph neural network architecture. The whole modeling pipeline of PE-GNN compared to a naive GNN approach is pictured in Figure 1. Taking GCN as the backbone architecture, PE-GCN is constructed as follows: Assuming a batch of randomly sampled points 1 , . . . , ℎ ∈ , a spatial graph is constructed from point coordinates c 1 , . . . , c ℎ using -nearest-neighborhood, resulting in adjacency matrix A . The point coordinates are then subsequently fed through the positional encoder (Θ ), consisting of the sinusoidal transform and a single fully-connected layer with sigmoid activation, embedding the 2 coordinates in a customizable latent space, returning vector embeddings c 1 , . . . , c ℎ = C . The neural network allows for explicit learning of spatial context, reflected in the vector embedding. We then concatenate the positional encoder output with the node features, to create the input for the first layer of our GCN as: The subsequent layers follow according to Equation 2.
To integrate the Moran's I auxiliary task, we compute the metric (Y ) for our outcome variable Y at the beginning of each training step according to Equation 5, using spatial weights from A . Prediction is then facilitated by creating two prediction heads, here linear layers, while the graph operation layers (e.g. GCN layers) are shared between tasks. Finally, we obtain predicted valuesŶ and(Y ). This auxiliary task learning process is also illustrated in Figure 1.
The loss of PE-GCN can be computed with any regression criterion, for example mean squared error (MSE): where denotes the auxiliary task weight. The final model is denoted as Θ ,Θ . Algorithm 1 describes a training cycle. Sample minibatch of ℎ points with features X , coordinates C and outcome Y .

5:
Construct a spatial graph with adjacency matrix A from coordinates C using -nearest neighbors 6: Using spatial adjacency A , compute Moran's I of output as (Y ) We begin training by initializing our model , for example a PE-GCN, with random weights on all layers and potential hyperparameters such as the PE embedding dimension. We also define our optimizer, setting hyper-parameters such as the learning rate. We now start the training cycle: At each training step, we first sample a minibatch of points from our training data. These points come as features X , point coordinates C and outcome variables Y . We then construct a graph from spatial coordinates C usingnearest-neighborhood, obtaining an adjacency matrix A . Next we use A as spatial weight matrix to compute local Moran's I values (Y ) from Y . As minibatches are randomly sampled, this creates a "shuffled" version of the metric. We then run inputs X , C , A through the two-headed model Θ ,Θ obtaining predictionŝ Y , (Ŷ ). We then compute the loss L (Y , (Y ),Ŷ , (Ŷ ), ), weighing the Moran's I auxiliary task according to weight parameter . Lastly, we use the loss L to update our model parameters Θ , Θ according to stochastic gradient descent. Training is conducted for (or until convergence) and after it is finished, the final model is returned.
PE-GNN, with any GNN backbone, helps to tackle the particular challenges of geographic patterns, outlined above: • While our approach still includes the somewhat arbitrary choice of -nearest neighbors to define the spatial graph, the proposed positional encoder network is not bound by this restriction, as it does not operate on the graph. This enables a separate learning of context-aware embeddings for each coordinate, accounting for neighbors at any potential distance within the batch. • While the spatial graph used still relies on pre-defined distance measure, the positional encoder embeds latitude and longitude values in a high-dimensional latent space. These high-dimensional coordinates are able to reflect spatial complexities much more flexibly and, added as node features, can communicate these throughout the learning process. • Batched PE-GNN training is not conducted on a single graph, but a new graph consisting of randomly sampled training points at each iteration. As such, at different iterations, focus is put on the relationships between different clusters of points. This helps our method to generalize better, rather than just memorizing neighborhood structures. • The differing training batches also help us to compute a "shuffled" version of the Moran's I metric, capturing autocorrelation at the same location for different (closer or more distant), random neighborhoods.

EXPERIMENTS 4.1 Data
We evaluate PE-GNN and baseline competitors on four real-world geographic datasets of different spatial resolutions (regional, continental and global): California Housing 1 : This dataset contains the prices of over 20, 000 California houses from the 1990 U.S. census and was introduced by [18]. The regression task at hand is to predict house prices using features x (e.g., house age, number of bedrooms) and location c. California housing is a standard dataset for assessment of spatial autocorrelation.
Election 2 : This dataset contains the election results of over 3, 000 counties in the United States and was proposed by [16]. The regression task here is to predict election outcomes using sociodemographic and economic features (e.g., median income, education) x and county locations c.
Air temperature 3 : The air temperature dataset [13] contains the coordinates of 3, 000 weather stations around the globe. For this regression task we seek to predict mean temperatures from a single node feature , mean precipitation, and location .
3d Road 4 : The 3d road dataset [17] provides 3-dimensional spatial co-ordinates (latitude, longitude, and altitude) of the road network in Jutland, Denmark. The dataset comprises over 430, 000 points and can be used for a regression task where altitude is predicted using latitude and longitude coordinates c. Note that this dataset does not contain any node features x, so it is only used for the spatial interpolation task.

Baselines
We compare PE-GNN with four different graph neural network backbones: The original GCN formulation, introduced by [20] and outlined in the Methods section, graph attention mechanisms (GAT) [31] and GraphSAGE [12]. Lastly, we use Kriging Convolutional Networks (KCN) [4]. KCN differs from GCN primarily in two ways: it transforms the distance-weighted adjacency matrix A using a Gaussian kernel and adds the outcome variable and features of neighboring points to the features of each node. To avoid a data leak, test set points can only access neighbors from the training set to extract these features. We compare the naive version of all these approaches to the same four backbone architectures augmented with our PE-GNN modules. Beyond GNN-based approaches, we also compare PE-GNN to the most popular method for modeling continuous spatial data: Gaussian processes. For all approaches, we compare a range of different training settings and hyperparameters, as discussed below.

Experimental setup
To allow for a fair comparison between the different GNN approaches, we equip all models with the same architecture, consisting of two GCN / GAT / GraphSAGE layers with ReLU activation and dropout, followed by a final linear layer as regression head. The KCN model also uses GCN layers, following the author specifications. We test four different auxiliary task weights = {0, 0.25.0.5, 0.75}, where = 0 implies no auxiliary task. Spatial graphs are constructed assuming = 5 nearest neighbors, following results from previous work [4,16] and our own rigorous testing. We include a sensitivity analysis of the parameter and different batch sizes in our results section. Training for the GNN models is conducted using PyTorch [27] and PyTorch Geometric [8]. We use the Adam algorithm to optimize our models [19] and the mean squared error (MSE) as loss function. Gaussian process models (exact and   Table 1: Spatial Interpolation: Test MSE and MAE scores from four different datasets, using four different GNN backbones with and without our proposed architecture. approximate) are trained using GPyTorch [11]. Due to the size of the dataset, we only provide an approximate GP result for 3d Road. All training is conducted on single CPUs. The code for PE-GNN can be accessed here: https://github.com/konstantinklemmer/pe-gnn.

Model
Cali  Table 2: Spatial Regression: Test MSE and MAE scores from three different datasets, using four different GNN backbones with and without our proposed architecture.

Results
We test our methods on two tasks: Spatial Interpolation, predicting outcomes from spatial coordinates alone, and Spatial Regression, where additional node features are available. The results of our experiments are shown in Table 1 and 2. We evaluate both models using the mean squared error (MSE) and mean absolute error (MAE) metrics.
For the spatial interpolation task, we observe that the PE-GNN approaches consistently and vastly improve performance for all four backbone architectures across the California Housing, Air Temperature and 3d Road datasets and by a small margin for the Election dataset.
For the spatial regression task, we observe that the PE-GNN approaches consistently and substantially improve performance for all four backbone architectures on the California Housing and Air Temperature datasets. Performance remains unchanged or decreases by very small margins in the Election dataset, except for the KCN backbone which benefits tremendously from the PE-GNN approach, particularly with auxiliary tasks.
Our results are interesting in several ways: Generally, PE-GNN substantially improves over baselines in regression and interpolation settings. Most of the improvement can be attributed to the positional encoder, however the auxiliary task learning also has substantial beneficial effects in some settings, especially for the KCN models. The best setting for the task weight hyperparameter seems to heavily depend on the data, which confirms findings by [22]. To our knowledge, PE-GNN is the first GNN-based learning approach that can compete with Gaussian Processes on simple spatial interpolation baselines, though especially exact GPs still sometimes have the edge. Nonetheless, PE-GNN is of course substantially more scalable than exact GPs, which rely on expensive pair-wise distance calculations across the full training dataset. Due to this problem, we do not run an exact GP baseline for the highdimensional 3d Road dataset. For KCN models, which include neighbor information as additional node features, we observe a proneness to overfitting. As the authors of KCN mention, this effect diminishes in large enough data domains [4]. For example, KCNs are the best performing method on the 3d Road dataset-by far our largest experimental dataset. We observe that in cases when KCN learns well, PE-KCN can still improve its performance, as highlighted by the 3d road experiments. The KCN experiments also highlight the strongest effects of the Moran's I auxiliary tasks: In cases when KCN overfits (Election, Cali. Housing datasets), PE-KCN without auxiliary task ( = 0) is not sufficient to overcome the problem. However, adding the auxiliary task ( = [0.25, 0.5, 0.75]) can mitigate most of the overfitting issue. This directly confirms a theory of [22] on the beneficial effects of auxiliary learning of spatial autocorrelation. Throughout our experiments, PE-GNN provides a powerful and flexible framework for learning with geogpraphic coordinates. The positional encoder allows for the learning of spatial context features to inform the downstream task, while the Moran's I auxiliary task helps to overcome potential overfitting issues.
We also assess the robustness of PE-GNN training cycles. Figure  3 highlights the confidence intervals of PE-GNN models with GCN, GAT and GraphSAGE backbones trained on the California Housing and 3d Road datasets, obtained from 10 different training cycles. We can see that training runs exhibit only little variability. These findings thus confirm that PE-GNN can consistently outperform naive GNN baselines. Lastly, Figure 4 highlights some results from our sensitivity analyses with the and ℎ (batch size) parameters. After rigorous testing, we opt for = 5-NN approach to create the spatial graph and compute the shuffled Moran's I across all models. We chose ℎ = 2048 for Cali. Housing and 3d Road datasets and ℎ = 1024 for the election dataset.

CONCLUSION
With PE-GNN, we introduce a flexible, modular new GNN-based learning framework for geographic data. PE-GNN leverages recent findings in embedding spatial context into neural networks to improve predictive models. Our empirical findings confirm a strong performance. This study highlights how domain expertise from applied academic disciplines like geography can help improve machine learning models for the particular data that is common in these disciplines. We hope to build on the foundations of PE-GNN to develop further methods for geospatial machine learning.


Model-based Safety and Security Co-analysis: a Survey

We survey the state-of-the-art on model-based formalisms for safety and security analysis, where safety refers to the absence of unintended failures, and security absence of malicious attacks. We consider ten model-based formalisms, comparing their modeling principles, the interaction between safety and security, and analysis methods. In each formalism, we model the classical Locked Door Example where possible. Our key finding is that the exact nature of safety-security interaction is still ill-understood. Existing formalisms merge previous safety and security formalisms, without introducing specific constructs to model safetysecurity interactions, or metrics to analyze trade offs.

Introduction
New technology comes with new risks: drones may drop on to people, self-driving cars may get hacked, medical implants may leak in people's body. Such risks concer both accidental failures (safety) and malicious attacks (security). Safety and security are heavily intertwined. Measures that increase safety often decrease security and vice versa: the Internet-of-Things offers ample opportunities to monitor the safety of a power plant, but their many access points are notorious for enabling hackers to enter the system. Passwords secure patients' medical data, but are a hindrance during emergencies. It is therefore widely acknowledged, also by international risk standards [22,24], that safety and security must be analyzed in combination [5,35]. To cater for this need, various risk frameworks have been developed for safety-security co-analysis. Process-oriented frameworks consider the steps needed for performing safety-security risk analysis [23,20]. Various formalisms specifically support the risk assessment phase within this process, i.e., the identification, analysis and evaluation of safety-security risks. Text-based methods include FMVEA [42], CHASSIS [41] and SAHARA [34]. Our focus is on model-based risk assessment. These formalisms provide a detailed insight in This work was partially funded by ERC Consolidator Grant 864075 (CAESAR). different failure and attack modes, mechanisms and their root causes. Understanding these is crucial for effective decision making. Model-based formalisms allow one to compute various dependability metrics. Such metrics quantify key performance indicators, like the availability and mean time to attack/failure. Problem statement. With various formalisms being proposed for the relevant purpose of safety-security co-analysis, it becomes natural to compare these, and assess their merits. Several questions emerge: (Q1) How these formalisms compare? (Q2) Which modeling constructs exist to model the interactions between safety and security? (Q3) Which analyses do these formalisms enable? (Q4) How do these formalisms compare on industrial case studies and (Q5) What would be desirable extensions?
Our literature search followed well-established methodological principles [9], focusing on peer-reviewed publications published until 2020. We queried Scopus, Microsoft Academic and Google Scholar, using the keywords "safety", "security", "dependency" and "model-based". We further explored relevant results through backward and forward search, as well as safety-security surveys [31,36,12].

Ref. Year #Citations
Fault Tree/Attack Tree Integration (FT/AT) [17] 2009 170 Component Fault Trees (CFTs) [44] 2013 52 Attack-Fault Trees (AFTs) [32] [45,30]. This is a classical example of safety-security antagonism: if locked, a door is unsafe in case of a fire; if unlocked, it is insecure in case of a burglary. While being a toy example, similar principles appear in, e.g., fire exits in buildings and emergency stops in factories. We present conclusions and directions for future research.

Related work.
A few papers have surveyed safety and security combinations: paper [31] compares non-model-based and model-based approaches, the latter being graphical or non-graphical. Survey [12] compares seven frameworks with respect to model creation, origin, stages of the risk assessment, and use cases. Our comparison is more technical than [12]. Finally, [36] summarizes practices from safety and security modeling separately. A joint model is proposed, combining the goal structuring notation GSL with Attack-Defense Trees.  Organization. Section 2 presents our findings. Section 3 introduces Fault Trees and Attack Trees. Section 4 presents our definitions for dependencies. Section 5, 6 and 7 compare the various formalisms and Section 8 concludes the paper.

Findings
Our survey revealed several noteworthy findings, summarized in Tables 1 and 2.  [44], Attack-Fault Trees (AFTs) [32]. Others, grouped in Class 2, deploy additional constructs: State/Event Fault Treess (SEFTs) [38] use Petri nets, Boolean driven Markov processes (BDMPs) [30] use Petri nets and triggers, Attack Tree Bow Ties (ATBTs) [6] extend bowties. A third class safety-security frameworks concerns security extensions of architectural languages: STAMP [19], SysML [37] and AADL [13]. We also place Bayesian networks [29] in this class. Finding #2: No novel modeling constructs are introduced. None of the formalisms for safety and security integration that we found introduces novel modeling constructs to capture safety-security interactions. Instead, they merge existing safety and security formalisms without adding new operators. Thus, one can represent safety and security features in one model, but one may wonder to what extend the interaction can be expressed appropriately. Finding #3: Safety-security interactions are still ill-understood. While the paper [31] coins the four dependency types, definitions can be improved. In particular, we propose rigorous definitions that focus on requirements and events, to then specify these for tree-based formalisms in Section 4. Finding #4: No novel metrics were proposed. No novel metrics were introduced to quantify safety-security interactions. Again, classical metrics were studied, such as the mean time to failure and attacker success probabilities. Furthermore, trade offs between safety and security, e.g., through Pareto analysis, were not studied either, even though trade offs are natural to analyze and does not necessitate novel analysis methods.
Finding #5: No large case studies were carried out. This question on whether or not additional constructs are needed becomes more pressing, since no large case studies were carried out. In general, analyzed papers present relatively small examples for illustrative purposes. A notable exception is the medium-size, but realistic, pipeline case study in [30] and [32]. This is remarkable, because for safety and security separately, such large case studies do exist [13].

Conclusion.
In conclusion, the true interaction between safety and security is still ill-understood, and ill catered for. Thus, future work should focus on a fundamental understanding between safety and security, so that dedicated modeling constructs and metrics can be proposed. We suggest starting from a realistic case study, and synthesizing the needs for usable safety-security frameworks that fit the future needs in this important area.

Attack Trees and Fault Trees
FTs and ATs are similar in nature. Section 5 compares plain combinations of FTs and ATs, while Section 6 surveys FTs and ATs extensions with additional features. We first introduce FTs and ATs themselves. Modeling. FTs and ATs model how low level failures (resp. attacks) propagate through the system and lead to system level failures (resp. attacks). FTs have been developed in the early '60s [16]. Due to their popularity, ATs were developed in the '90 as their security counterparts [43]. FTs (resp. ATs) start with a Top Level Event (TLE), modeling a system level failure (attack), which is then refined through Boolean gates: the AND-gate indicates that all children must fail (be attacked) in order for the gate to fail (be attacked).
For the OR-gate to fail (be attacked), at least one of its children need to fail (be attacked). When refining is no longer needed, one arrives at leaves of the tree: The Basic Events (BEs) in FTs model atomic failures; the Basic Attack Step (BAS) in ATs model atomic attack steps. Non-leaves nodes, such as Locked in during fire and Attacker walks through door, are called intermediate events, denoted by rectangles. Despite their name, FTs and ATs are Directed Acyclic Graphs (DAGs), rather than trees, since subtrees can have multiple parent gates. The FT in Figure 1 models that integrity is compromised if one is either locked in during a fire or the door lock fails to close. This is modeled via the top ORgate. One is locked in during a fire if there is a fire and the door is locked. This is modeled via the AND-gate connecting the BEs Fire and Locked.
Analysis. FT and AT enable numerous analysis methods [28]: qualitative analyses include Minimal Cut Sets (MCSs), indicating which combinations of BEs or BASs lead to the TLE. The set {Fire, Locked} is a cut set in Figure 1. Quantitative analyses compute dependability metrics, such as the system reliability, attack probabilities and costs. For example, by equipping the BEs and BASs with probabilities, one can compute the likelihood of a system level failure or attack to occur. Both FTs and ATs are part of international standards [21] and have been used to analyse numerous case studies [46,18,11]. FTs and ATs also feature some differences: FTs often focus on probabilities, whereas ATs consider several other attributes, like cost, effort and required skills. Further, FTs have been extended with repairs [40], and dynamic gates [25,15]; ATs with defenses, and sequential AND (SAND) gates [27,18].

Dependencies in ATs and FTs Combinations/Extensions
This section proposes formal definitions of the four dependency types from [31] in the specific context of FTs and ATs. While [31] coins these dependencies, it does not give very precise definitions. The paper focuses on dependencies between requirements and/or measures, as do most examples (e.g., firewalls increases security, but decrease safety). We refine these definitions by focusing on events in the formalisms combining FTs and ATs. Note that these events correspond to requirements: the TLE provides the main disruption to be avoided; other events refine the TLE into subrequirements to be prevented. As such, these events and their interactions are the inverse of the logical interactions between the requirements. For each dependency type, we investigate to what extend these are expressible in the various FT and AT formalisms. These dependencies are to a large extend generic, independent of the safety-security context. When fit, we specialize to the safety-security interdependencies.
Antagonism. Two undesirable events A, B are antagonistic, or conflicting, if there is no situation in which neither event has occurred. In tree-based formalisms, let A and B have a shared security/safety event C as a child and let it be connected to either A or B exclusively through a NOT-gate. We then say that A and B are antagonistic with respect to C. The events Locked in during fire and Attacker walks through door in Figure 2b are antagonistic with respect to Door unlocked. As such, expressing antagonism requires a form of negation. . By assumption, all leaves in FTs and ATs are statistically independent. Events that are (mutually) reinforcing can also satisfy this statistical independence requirement. Thus, in tree-based formalisms the absence of (mutual) reinforcement is also required to capture Independence as intended. All tree-based formalisms from Sections 3 to 6 can express independence.
Observations. Since these formalisms combine FTs and ATs rather than developing new formalisms, it is difficult to support safety-security trade offs. No novel metrics were introduced either. It is of interest to extend trees by using two roots, for safety and security goals. These formalisms can be extended with NOT-gates. To assess their modelling capabilities, larger case studies are needed.

Formalisms combining Fault Trees and Attack Trees
This section surveys three combinations of FTs and ATs: Fault Tree/Atack Tree refines a FT's BEs with ATs; Component Fault Trees merges FTs and ATs, while Attack-Fault Trees merge dynamic FTs and dynamic ATs. Our figures use blue for events related to safety, yellow for security, and green for their combination.

Fault Tree-Attack Tree Integration
Modeling. FT/ATs [17] are based on the assumption that attackers try to force a system failure, expressed by the TLE of a FT, to occur. To do so, FT/ATs investigate how the BEs of an FT can be evoked by an attacker, refining these by ATs with the BE as goal. Figure 2a refines basic event Door lock fails in Figure 1, by making this glsbe the goal of the AT depicted in yellow. There is no clear method for modeling dependencies of components in the AT on parts of the FT, so mutual reinforcement can not be modeled due to a lack of bidirectionality. Conditional dependency and antagonism can be modeled as per Section 4.
Analysis. The paper [17] computes the probability for the TLE to occur, given probabilities for the BEs and and the BASs.
Observations. Exploiting safety faults is a common method for hackers to enter a system, e.g.: triggering a fire alarm to disengage a fire safety lock. Thus, it seems reasonable to consider the leaves of a FT as vulnerabilities, refining them as ATs. The fault tree's TLE then becomes the target for hackers. This can be reasonable e.g., shutting down a factory for ransomware. However, attackers are often exploit safety faults to achieve other goals, such as stealing digital assets.
In that case, the forcing of BEs is only a starting point. Subsequent attack steps could be modeled in an AT; this could be a natural extension for FT/ATs. Paper [17] abstractly models a chemical plant as a case study

Component Fault Trees
Modeling. CFTs equip FTs with a modular structure [26], so that a large FT can be modeled and analyzed in terms of smaller components. The paper [44] extends CFTs with security aspects by introducing a new BEs type for security breaches, which are essentially BASs. CFTs make no distinction between BEs and BASs. This is exemplified in Figure 2b, where attacks and failures are freely merged. As per Section 4, CFTs model mutual reinforcement, conditional and antagonistic dependencies. Analysis. The CFTs from [44] enable the same analysis methods as FTs. The authors remark that, if a BAS only occurs in MCSs with multiple, low probability and independent BEs, then this attack is very unlikely to ever cause a disruption, as it requires multiple other unlikely events to occur. However, one may remark that the same holds for any event.
Observations. Not distinguishing between safety and security has the advantage that all existing tools remain applicable. Moreover, one may ask if it really matters whether a disruption is due to a failure or an attack. A disadvantage of merging failures and attacks, however, is that no interactions or trade offs between safety and security can be studied. Paper [44] models an Adaptive Cruise Control system as a case study.

Attack-Fault Trees
Modeling. AFTs [32] treat attacks and failure in the same way of CFTs. However, AFTs merge dynamic ATs and dynamic FTs. These provide additional gates to model dynamic behaviour: dynamic ATs include the sequential ANDgate (SAND), modelling attacks as sequence of steps [4]. Dynamic FTs include gates for modeling spare components (SPARE), functional dependencies (FDEP) and priority ANDs (PAND). The AFTs leaves can be BEs, BASs and instant failures (IFAIL). Further, attacker profiles quantify the attacker's capabilities, such as available resources, skills and damage. These appear as attributes on the BASs. Figure 2c presents the Locked Door Example as a AFTs. The various SAND-gates indicate the order of events. E.g., for the Door lock breached event to happen, first Attack initiated must happen, and then the Access event as well.
In [32], the status of the door is modelled by two IFAIL leaves, Door locked and Door unlocked: if the door is locked, then the probability of it being unlocked is 0. Dependencies can thus be modeled as described in Section 4. Analysis. AFTs are analysed by translating the AFTs to a network of stochastic timed automata, analyzed via statistical model checking. Using the attribute values in the attacker profiles, several metrics can be computed, such as the time, cost and likelihood of the attacks. Pareto frontiers elucidate trade offs between attributes, e.g., the likelihood of an attack within a given budget.
Observations. AFTs share the two disadvantages common to all FT approaches in that there is no distinction between safety or security failures once the TLE occurs, and the methods of modeling interactions may not clearly highlight the antagonistic dependencies. The attacker profiles in AFTs support a wide range of quantifiable parameters, enabling versatile analysis and trade offs. Remarkably, AFTs are used to model the medium-sized case study of an oil pipeline, in addition to modelling the Locked Door Example.

Formalisms extending Fault Trees and/or Attack Trees
This section surveys combinations of FTs and ATs with additional features:

State/Event Fault Trees join a fault tree-like model with Petri nets, Boolean Driven Markov Processes extend FTs and ATs with Petri nets and triggers and
Attack Tree Bow-ties combine Event Trees with an FT/AT-like model.

State/Event Fault Trees
Modeling. State/Event Fault Trees (SEFTs) [38] join FTs and Petri nets † , expressing that certain failures can only happen in certain states. Whereas the leaves in ATs and FTs model atomic events, SEFTs deploy Petri nets to accommodate state changes inside basic events. In Figure 3, the Door component can move between the states Unlocked and Locked. These state changes are triggered by events, depicted as black rectangles, that can be exponentially distributed, deterministic, or triggered by other events. Both states and events can be communicated via the gates of the tree, via in and out ports. In this way, Figure 3 expresses that a fire casualty can only happen if the door is locked. Figure 3 models antagonism: the door should be unlocked to escape from Fire but locked to prevent failure of the AND-gate Enter unlocked door.  Analysis. SEFT support the same dependability metrics as attack-fault tree combinations. The tool ESSaReL translates SEFTs to extended deterministic stochastic Petri nets, which can be further analyzed by the TimeNet tool, e.g., for steady state analysis.
Observations. The authors of [38] model the small-sized example of a Tyre Pressure Monitoring System. Even though not mentioned explicitly in [38], it appears that Petri nets in the SEFT leaves must be disjoint between components. E.g., it is not possible to use the Unlocked state of the Door component as a direct input for an attack step. Allowing this could increase the expressivity of SEFTs.

Boolean Driven Markov Processes
Modeling. Boolean Driven Markov Processes [8] extend FTs by equipping each leave with a Markov process (MP), representing the different modes a component can be in. Various templates provide standard MPs to model standard failure behavior. For example, the failure in operation MP contains two modes, operational and failed. One transitions from operational to failed with an exponential failure rate λ, and back to operational with a repair rate µ. The IFAIL MP models instantaneous failures. Moreover, users can define their MPs as a stochastic Petri net [30], similarly to SEFTs. In [7], BDMPs are extended with security aspects by providing additional Markov processes for attacker steps. For example, the Attacker Action MP (AA) contains tree modes: in Idle the attacker has not yet initiated the attack. The Active mode corresponds to actual attempt, requiring an exponentially distributed time to succeed, and leading to the success mode. Further, triggers, represented by dotted red arrows, allow one MP to trigger a mode change in another MP. Figure 4 shows the Locked Door Example from [30]. The triggers pointing from Attack initiated to the OR-gate means that the OR-gate is not activated until Attack initiated happens. The use of Petri nets and the presence of intermediate events allow BDMPs to model all dependencies between safety and security, as detailed in Section 4. Analysis. BDMPs allows both quantitative and qualitative analysis. The authors of [7] build BD-MPs with the KB3 modeling software platform. The computation of the overall mean time to success (MTTS), the probability of success, and the list of possible attack success sequences is possible.
Observations. BDMPs are the only formalism, compared to the others of this work, that handles a medium size case study [10]. Furthermore, the authors discuss dependencies corresponding to those mentioned in the introduction. Because our definitions are more detailed, those do not coincide with our dependencies. Modeling. Attack Tree Bowties (ATBTs) [1] combine Bow-ties with ATs. Bowties [14] themselves combine FT and Event Trees (ETs): the left part of a Bow-tie is a FT modeling the causes of an hazardous event, which is in the middle of the Bow-tie. The ET on the right models its consequences. Barriers, i.e., a measure M preventing some failure F from happening, are modeled as AND(F, M ), so that the failure F only propagates if the barrier M fails. Now, ATBTs extend regular Bow-ties by attaching ATs to the basic events of a FT, just as in FT/ATs. Thus, in ATBTs, the left part of the Bow-tie is an FT/AT, rather than an FT. Figure 5 models the Locked Door Example via ATBTs. The leftmost part is equivalent to the FT/AT from Figure 2a. The ET on the right details the consequences in several cases: if the failure was a fire or a burglary, if alternate escape routes were available or inner doors were locked. Dependencies in ATBTs are the same as in FT/ATs, as further explained in Section 4. Analysis. Just like for regular Bow-ties, [1] identifies vulnerabilities via MCSs. By assigning likelihood level to all BEs and BASs, two likelihood levels are assigned to these cut sets: one for safety and one for security. Thus, trade offs can be made. Observations. Since the left, causal part of ATBTs is similar to FT/ATs, similar observations apply. Further it would be natural to study trade offs between safety and security by equipping ATBTs with two hazardous events, one for safety and one for security. Moreover, ATBTs are one of the few formalisms that consider safety-security trade offs. In [1] the authors create a small case study of a risk scenarioin a chemical facility.

Architectural formalisms and Bayesian Networks
Formalisms not based on trees are STAMP, SysML and AADL -which are all based on the system architecture -and Bayesian networks. Modeling. The System-Theoretic Accident Model and Processes (STAMP) is rooted in the observation that system risks do not come from component failures, but from inadequate control or enforcement of safety-related constraints. "In STAMP, systems are viewed as interrelated components that are kept in a state of dynamic equilibrium by feedback loops of information and control" [33]. Each component enforces the safety and security constraints in the processes it controls, using control actions and feedback messages. Inability to enforce these constraints results in failures in safety or security. System Theoretic Process Analysis (STPA) and its extension STPA-safesec [19] systematically identify the consequences of incorrect control actions and feedback, e.g., when these happen too early, in the wrong order, or were maliciously inserted. Figure 6 shows the Person, who can lock and unlock the door. Locking mechanism is controlled by Person, and controls if the Door can open. A safety constraint is that Person must be able to unlock a door in case of fire; a security constraint is an unauthorised person must not be able to gain access. A violation is, e.g., the scenario where the person locks the locking mechanism while the door is open, forcing the door to stay open and granting unauthorised access. STPA-safesec can discover these risks in a structured manner, however this is a manual process requiring domain knowledge and cannot be automated. STAMP is not geared towards expressing dependencies. However, STPA-safesec analysis may reveal safety-security conflicts. Analysis. STPA is used to identify potential hazards and undesired behaviours: e.g., the violation previously described. Observations. STAMP models control flows of the system and further analysis identify safety and security hazards in that control flow. Domain experts are required to properly identify issues. STAMP provides a structured way of reasoning based on a high level description of the system that should ensure the identification of safety and security requirements. STAMP can identify safety and security issues in a system, it is not for documenting those interactions. It is suited to help with the syntheses of a model describing safety and security interactions. Paper [19] uses the system connecting synchronous-islanded operating microgrids, microgrids that locally generate power separately from the main power grid, to the main grid as a case study. Modeling. The Systems Modeling Language (SysML) is a general-purpose modeling framework for systems engineering, and extends the Unified Modeling Language (UML). SysML-sec [39] extends SysML to express safety and security requirements. Its functional model describes the communication channels between processes, detailing the encryption methods and their complexity overhead cost. SysML-sec also includes a system mapping model. SysML-sec enables both safety and security properties to be expressed and verified via separate model checkers. E.g., checking if confidential communications can be intercepted by compromising the bus. The interaction between safety and security properties is not modeled. Since SysML-sec is geared to embedded software, we replaced the Locked Door Example with a digital keypad lock. Figure 7 shows the mapping model. Within the keypad, the CPU is connected to the main memory and an Input-Output Bus. Connected to this bus is the physical button pad and a Digital Analog Converter that engages and disengages the lock. Figure  8 details the process of (dis)engaging the lock. The correct codes are stored in an encrypted format, and an input code is received unencrypted. Both codes are collected, and the input code is checked against the stored key, either by decrypting the stored key or encrypting the input key. On a match, a command is sent to (un)lock the door. Analysis. The TTool [37] can encode ATs and FTs in a SysML parametric model. It then uses UP-PAAL to test for reachability of the TLE. Counter measures can be annotated to the ATs and FTs models, and individual attack vectors/BEs switched on or off. UPPAAL will then indicate if preventing the chosen subset of attack vectors/BEs is sufficient to prevent the TLE. No method for capturing the ATs and FTs is included. TTool can also verify security requirements with ProVerif [2]. Requirements can be tagged as e.g. Confidentiality, Non Repudiation, Data Origin Authenticity [3]. For example, a confidentiality requirement on all communications can be created, and comm channels from the functional model linked to these requirements. An unsecured bus does not satisfy the confidential requirement, and with no other mitigation e.g., encrypted messaging, ProVerif will find a state where the confidential requirement of the comm channel is violated, and return a trace showing how this state was reached. Safety and security are modeled and analyzed separately. Thus, violation of safety and security requirements can be observed, as well as attempts at mitigation violating other requirements, but no direct interactions.

SysML
Observations. SysML is not suitable for generally visualising the interactions between safety and security in arbitrary systems, but for designing complex embedded systems with safety and security requirements. Paper [37] uses Key Masters Keying Protocol, which aims to securely distribute a randomly generated key among a group of in-car Electronic Control Units, as a case study. Modeling. AADL models the software and hardware architecture of embedded, real-time systems. It consists of a core language and several annexes: the core describes the multithreaded, distributed software architecture, and the annexes describe real-time behaviour, and error modelling. To detect software failures (safety), the AADL standard contains the EMV2 language.Safety and security are added in one AADL model in [13]. Safety analysis computes the mean time between failures (MTBF) of the system. For security, the goal is to avoid unauthorized access to sensitive data. Data is accessed through AADL ports of subprograms. So it is important to verify that the indirect data access points are secure enough. A dedicated property set is used to associate a security level, which is an integer value, to AADL data. Figure 9 models the Locked Door Example with AADL: this is uncommon for AADL because it is often used for the interaction between software and hardware components. Here we consider a Person, a Locking mechanism, and Door access as components. They are linked through ports. AADL does not model any kind of dependencies: safety and security are considered separately, although they share the same AADL model to perform the analysis. Analysis. For safety, the MTBFs of the system is calculated via the AADL Error model statements which are in various components descriptions. They are compiled together to generate a FT. For security, the goal is to avoid unauthorized access to sensitive data via the Stood for AADL tool.

The Architectural Analysis and Design Language AADL
Observations. An AADL model is scalable and appropriate for version and configuration management, because it can be fully described by its textual representation.FTs can be analyzed after transforming AADL to another tool which makes AADL an excellent application for FTs. Interesting for future work could be to integrate also ATs with AADL. There is a small case study in [13] on the locked door created by the authors. (BN) is a probabilistic graphical model that represent probabilistic dependencies between several variables via a directed acyclic graph. Each node A represents a variable, and an edge from A to B indicates that A stochastically depends on B. A conditional probability table yields the conditional probabilities P [A|B]. If the probabilities of the leaves in the BN are known, the probabilities of the root nodes can be calculated. In [29], BNs are proposed to model safety and security dependencies. The two root nodes represent system safety and security. Figure 10 expresses that Door locked is a common factor for safety and security. BNs model conditional dependency, mutual reinforcement, and antagonism. Analysis. Qualitative analysis can be done in the form of conditional independence analysis (analyzing which nodes influence other nodes and how). BNs enable quantitative analysis to calculate reliability metrics such as mean time to failure.

Modeling. A Bayesian Network
Observations. Fault trees and attack trees can be seen as special cases of Bayesian networks, where the CPT tables encode the Boolean gates, e.g., P [A = 0|B = 1, C = 1] = 1, 0 for other B and C, for an AND-gate A with children B and C. Thus, the BNs extend ATs and FTs with flexible dependencies, and separate TLEs for safety and security. However, the gates in ATs and FTs provide a clearer visualization of the gates, since in BNs these must be read from the probability tables. Paper [29] uses the pipeline example as a case study.

Conclusion and Future Work
We surveyed 10 model-based formalisms that capture safety/security interdependencies. Out of these, only AFTs, BDMPs and BNs fully capture all four interactions: antagonism, mutual reinforcement, conditional dependency and independence. None of the formalisms present specific constructs for safety/security interdependencies: they are acquired from safety-only or security-only frameworks and joint afterwards, following different strategies. The same holds for metrics: none of those is specific for safety/security interactions. To the best of our knowledge, safety/security interactions are studied inside limited scenarios and large industrial case studies are still missing. More rigorous definitions of safety/security dependencies are needed, to account for: 1. Directionality. Are safety and security directional or bi-directional and from which direction do they flow? 2. Intensity. For a quantifiable co-analysis, intensity of these interaction has to be considered. 3. Nature of the interaction. For each of the possible interactions, from influence, to dependency or antagonism, accounting for the positive or negative impact of such an interaction is fundamental. Further extensions could be considered for every FT/AT-based formalism, i.e., using two roots to account for different safety and security TLEs, similarly to BNs. Analyzed formalisms join analysis methods and metrics from the fields of security and safety: this approach does not consider the possibility to develop metrics that are instead safety-security specific. The visualization of safety/security interdependencies should also be considered to ease readability of complex models.