Cesar Torres
email: ctorres@tamps.cinvestav.mx

Huitzil Cinvestav

Tamaulipas Mexico

Oleksandr Popovych
email: oleksandr.popovych@inria.fr

Bernard Girau
email: bernard.girau@loria.fr

Fault Tolerance of Self Organizing Maps

As the quest for performance confronts resource constraints, major breakthroughs in computing efficiency are expected to benefit from unconventional approaches and new models of computation such as brain-inspired computing. Beyond energy, the growing number of defects in physical substrates is becoming another major constraint that affects the design of computing devices and systems. Neural computing principles remain elusive, yet they are considered as the source of a promising paradigm to achieve fault-tolerant computation. Since the quest for fault tolerance can be translated into scalable and reliable computing systems, hardware design itself and the potential use of faulty circuits have motivated further the investigation on neural networks, which are potentially capable of absorbing some degrees of vulnerability based on their natural properties. In this paper, the fault tolerance properties of Self Organizing Maps (SOMs) are investigated. To asses the intrinsic fault tolerance and considering a general fully parallel digital implementations of SOM, we use the bit-flip fault model to inject faults in registers holding SOM weights. The distortion measure is used to evaluate performance on synthetic datasets and under different fault ratios. Additionally, we evaluate three passive techniques intended to enhance fault tolerance of SOM during training/learning under different scenarios.

I. INTRODUCTION

Artificial neural networks are computational models that have attracted intensive research interest and enjoyed significant renewed growth in artificial intelligence related applications. In neural networks research, one of the main problems is the architecture optimization, which aims at appropriately choosing the neural architecture and its parameters for high generalization/performance at solving a given task. However, the fact that performance maximization is of primary concern does not necessarily imply that it is the only goal that has been or should be pursued [START_REF] Alippi | Selecting accurate, robust, and minimal feedforward neural networks[END_REF]. Artificial neural networks are generally assumed to acquire some other desirable features of biological systems such as their tolerance against imprecision, uncertainty and faults, which make them even harder to study or design. Fault tolerance is the attribute of a system that allows it to preserve its expected behavior after faults have manifested themselves within the system [START_REF] Avižienis | Framework for a taxonomy of fault-tolerance attributes in computer systems[END_REF]. More precisely, a faulttolerant system might be defined as one that has provisions to avoid failure after faults have caused errors within the system. When a statement about neural networks fault tolerance is made, it should be implicitly assumed a failure condition or criterion related to its intended computational task. That is, the threshold below which it cannot longer perform its function according to the specification. As such, fault tolerance in neural networks depends on the definition of the acceptable degree of performance and its intended application [START_REF] Protzel | Performance and faulttolerance of neural networks for optimization[END_REF].

According to neurobiological studies, the human brain is able to tolerate a small amount of synapse or neuron faults or even use noise as a source of computation [START_REF] Maass | Noise as a resource for computation and learning in networks of spiking neurons[END_REF]. Nervous systems are complex, highly parallel information processing architectures made of seemingly imperfect and slow, but exceptionally adaptive and power-efficient components that carry out information processing functions [START_REF] Sejnowksi | The language of the brain[END_REF]. Moreover, brains have the capability to relearn by growth of new neurons and/or neural connections and/or retraining of the existing neural architecture. Derived from these observations, it is commonly claimed that the majority of neural network models, abstracted from biological ones, have built-in or intrinsic fault tolerance properties due to their parallel and distributed structure, and the fact that they contain more neurons or processing elements than the necessary to solve a given problem, i.e., some natural redundancy due to overprovisioning. Claiming an equivalent fault tolerance only on the basis of rough architectural similarities therefore cannot hold true in general, especially for small size neural networks [START_REF] Venkatesh | Robustness in neural computation: random graphs and sparsity[END_REF]. Studies have shown that an artificial neural network has a very limited fault tolerance capability and, as a matter of fact, most neural networks cannot be considered intrinsically fault tolerant, without a proper design. Obtaining truly fault tolerant neural networks is still a very attractive and important issue both for i) artificial intelligence based solutions, where, for instance, pervasive embedded systems will require smart objects fully merged with the environment in which they are deployed to cope with unforeseeable conditions [START_REF] Wang | Overcoming computational errors in sensing platforms through embedded machine-learning kernels[END_REF], and ii) as a source for reliable computing systems built from unreliable components, as suggested in [START_REF] Neumann | Probabilistic logics and the synthesis of reliable organisms from unreliable components[END_REF]. Rooted on neural computing, a new paradigm is needed to take advantage of new emerging devices at nanoscale dimensions and deal with both manufacturing defects and transient faults. This trend might lead to even consider faults/errors as an essential/intrinsic part of a system design. In this last direction, the robustness and the potential fault-tolerant properties of neural models call for attention as permanent and transient faults, device variation, thermal issues, and aging will force designers to abandon current assumptions that transistors, wires, and other circuit elements will function perfectly over the entire lifetime of a computing system, relying mainly on digital integrated circuits [START_REF] Rahimi | Variability mitigation in nanometer cmos integrated systems: A survey of techniques from circuits to software[END_REF].

In the literature, several experimental and less analytic works have been carried out to study issues related to neural networks fault tolerance, including the analysis of the effect of noise on the output sensitivity [START_REF] Hammadi | Improving the performance of feedforward neural networks by noise injection into hidden neurons[END_REF], the weight error sensitivity [START_REF] Zeng | Sensitivity analysis of multilayer perceptron to input and weight perturbations[END_REF], and the relationship among fault tolerance, generalization and model complexity [START_REF] Alippi | Selecting accurate, robust, and minimal feedforward neural networks[END_REF], [START_REF] Venkatesh | Robustness in neural computation: random graphs and sparsity[END_REF], [START_REF] Tchernev | Investigating the fault tolerance of neural networks[END_REF]. Such works have been carried out at different levels of abstraction, from very specific low level physical implementations, but mostly, to the high level intrinsic fault masking capacity of neural paradigms. It is worth to point out that most works have been focused on feedforward neural networks and few attempts have been made to study and improve fault tolerance of other neural models, such as self-organizing maps (SOMs). Yet, we consider that the principle of self-organization itself might have an influence onto the fault tolerance of computing systems. Therefore in this paper, we are not only interested in the way pre-learned SOMs react to faults, but also in the way fault tolerance evolves during the self-organization of the neural population. This is why we chose to study SOMs because they stand as one of the most well-recognized models of self-organizing computing systems. In this effort, the aim of this work is twofold: i) to study the intrinsic fault tolerance capabilities of SOMs under bit-flip faults on weights by analyzing their performance degradation with variable fault ratios using synthetic datasets, and 2) to apply/adapt some reported techniques for enhancing fault tolerance during training/learning in neural networks to SOMs and compare their performance under different testing scenarios. We are interested in faults in weights, since when a SOM is intended to be implemented onto a digital hardware substrate, the value of each weight vector can be disturbed by various causes such as electromagnetic field or radiation, which can be modeled by a bit-flip fault model.

The paper is organized as follows. Section II presents some relevant related work on fault tolerance in SOMs. Section III presents the fundamentals of the Kohonen SOM model, and the on-line learning mechanism for self-organization. Section IV first describes the general approach for fault tolerance assessment in neural networks, and then the specificities for SOMs when targeted to a fully parallel digital hardware implementation are presented. Section V presents three techniques that have been used to improve the intrinsic fault tolerance of SOMs, which involves a modification of the learning/training mechanism. Section VI describes the datasets used for testing and the obtained experimental results for different scenarios.

II. RELATED WORK

Despite of their importance, the studies and results on fault tolerance reported in the literature are mostly concerned with feedforward neural networks. However, such results are difficult to generalize across different neural models. Among other neural models, the Kohonen SOM is a neural network which main property is to perform data quantization while preserving some topological relations among the set of learned prototypes. This model is a recognized tool for data visualization, and it has also been used for many practical applications such as pattern classification, image processing, and robotics [START_REF] Kohonen | The self-organizing map[END_REF]. The SOM model is inspired by the organizational structure of feature maps in the brain. In a SOM, each neuron iteratively learns to extract a feature vector from the input data. Simultaneously this learning maintains neighborhood relations (close neurons learn close prototypes). This learning, property of effectively creating spatially organized internal representations, is performed autonomously in the SOM and is considered as a self-organization process. Such networks self-organize, learn through associations, and scale well to very large systems due to their regularity and modularity, potentially leading to easy hardware implementation in existing technologies: either dedicated very large scale integrated (VLSI) chips or field programmable gate array (FPGA) structures.

However, the self-organizing behavior of SOMs in the presence of faults and its fault-tolerant capability have not yet been fully explored, in contrast to the more in-depth existing studies for feedforward networks. As some previous works have shown, SOMs are expected to tolerate defective neuron weights thanks to their self-organizing mechanisms. For instance, in [START_REF] Yasunaga | Fault-tolerant self-organizing map implemented by wafer-scale integration[END_REF] authors studied the fault tolerance capability of SOM in the presence of defective neurons undergoing stuck-at faults. In this study, it was shown that defective SOMs can eventually re-organize themselves, by relearning, if the defective neuron stuck-output is larger than a critical value. Only a linear array was analyzed and discussed, where defective neurons were concentrated in one place in the array, forming what they called a defective-neuron cluster. In the experiments, 100 neurons, including six defective ones, were tested to show that the SOM can learn in spite of some faulty neurons. In [START_REF] Talumassawatdi | Fault immunization concept for self-organizing mapping neural networks[END_REF], authors addressed fault tolerance improvement in SOMs by using a technique called fault immunization of the synaptic connections. Stuck-at-a faults, a is a real value, were considered. Only one neuron was faulty at any time, but no restriction on the number of faulty links of the neuron was assumed. Weights are immunized by adding a constant value that is increased or decreased as much as possible without creating any misclassification. Fault immunization was formulated as an optimization problem on finding the corresponding constant value for each neuron. The application of this study was oriented to enhancing the reliability of the lossless image compression by a SOM. Other related works are related to modifications of the Kohonen SOM on-line learning algorithm so as to improve the fault tolerance of the learned SOM. Section V describes more precisely these works.

III. SELF ORGANIZING MAPS

A. Fundamentals

Self-organizing maps (SOMs), as proposed by Kohonen [START_REF] Kohonen | The self-organizing map[END_REF], consist of neighboring neurons commonly organized on one-or two-dimensional arrays that project patterns of arbitrary dimensionality onto a lower dimensional array of neurons. All neurons receive the same input pattern and an iterative mechanism updates the neuron's weights so as to learn to quantize the input space in an unsupervised way. This mechanism first selects the neuron whose weights best fit the given input pattern, and then brings the weights of this neuron and of its neighbors closer to the current input pattern.

More precisely, each neuron in a SOM is represented by a d-dimensional weight vector, m 2 R d , also known as prototype vector, m = [m 1 , ..., m d ], where d is the dimension of the input vectors, x. Neurons are connected to adjacent ones by a neighborhood relationship, which define the structure of the map. The mechanism for selecting the winning neuron requires a centralized entity, so that the Kohonen SOM is not a fully distributed model as the cortex organization [START_REF] Rougier | Dynamic self-organising map[END_REF]. After learning, or self-organization, two vectors that are close in the input space will be represented by prototypes of the same or of neighboring neurons on the neural map. Thus, the learned prototypes become ordered by the structure of the map, since neighboring neurons have similar weight vectors.

In the literature, two main types of SOM can be distinguished, considering the number of neurons in the map with respect to the number of expected clusters: SOMs in which the number of neurons is roughly equal to the number of expected clusters to be found in the input, and SOMs with a large number of neurons (thousands or tens of thousands), which are used to allow the emergence of intrinsic structural features of the data space; these SOMs are referred to as Emergent Self-Organizing Maps (ESOMs) [START_REF] Ultsch | Data mining and knowledge discovery with emergent selforganizing feature maps for multivariate time series[END_REF].

B. On-line mode training

The basic version of SOM learning is an on-line stochastic process which has been inspired by neurobiological learning paradigms. Several other extensions of the algorithm have been proposed [START_REF] Cottrell | Theoretical and Applied Aspects of the Self-Organizing Maps[END_REF]. In the SOM on-line mode algorithm, the learning phase starts with an appropriate (usually random) initialization of the weight vectors, m i . The input vectors are presented to the neural map in multiple iterations. For each iteration, i.e. for each input vector x, the distance from x to all the weight vectors is calculated using some distance measure. The neuron whose weight vector gives the smallest distance to the input vector x is usually called the best matching unit (BMU), denoted by c, and determined according to:

kx m c k = min i kx m i k (1)
where k•k is the distance measure, typically the Euclidean distance, x is the input vector and m i is the weight vector of neuron i. The winner c and its neighboring neurons i 2 N w update their weights according to the SOM rule:

m i (t + 1) = m i (t) + ↵(t)h ci (t)[x(t) m i (t)] (2)
where t denotes the time, x(t) is an input vector randomly drawn from the input data set at time t, ↵(t) the learning rate at time t, and h ci (t) is the neighborhood kernel around c. The learning rate ↵(t) defines the strength of the adaptation, which is application-dependent. Commonly ↵(t) < 1 is constant or a decreasing scalar function of t.

The neighboring kernel h ci (t), which is a function of the distance between the winner neuron c and neuron i, can be computed using a Gaussian function:

h ci (t) = exp( kr c r i k 2 2 2 (t) ) (3)
The term kr c r i k is the distance between neuron i and the winner neuron c. The precise value of (t) does not really matter, as long as it is fairly large in the beginning of the process, for instance in the order of 20% of the longer side of the SOM array, after which it is gradually reduced to a small fraction of it, e.g. 5% of the shorter side of the array [START_REF] Kohonen | The self-organizing map[END_REF].

A batch version of SOM learning exists. It is also iterative, but instead of using a single input vector at a time, the whole dataset (batch) is presented to the map before updating any weight [START_REF] Kohonen | The self-organizing map[END_REF]. In each training step, the dataset is partitioned according to the Voronoi regions of the map weight vectors, i.e., each data vector belongs to the dataset of the closest map unit. This algorithm is deterministic, and one of its advantages is that the limit states of the prototypes depend only on the initial choices. However, in this paper, we are interested in the way fault tolerance evolves with SOM learning. The on-line mode training provides us with a more detailed evolution to observe. Moreover, we not only consider the standard Kohonen SOM learning algorithm, but also some variants among which several ones cannot be applied in batch mode. Therefore, we will only focus on on-line mode training.

IV. FAULT TOLERANCE IN NEURAL MODELS

At a high level of abstraction, neural networks fault tolerance can be analyzed by the effects of errors in the main operators, rather independent from their intended physical implementation. In a more comprehensive approach [START_REF] Piuri | Analysis of fault tolerance in artificial neural networks[END_REF], after this initial step, physical faults affecting a specific implementation can be mapped onto such errors to estimate fault tolerance of a given neural model, and by identifying critical components, complementary and ad-hoc fault tolerance policies can be further applied to enhance the properties of the neural model implementation. In order to study fault tolerance, researchers develop models of them to examine the variety of faults that need to be tolerated during the operation of a given system.

A. Faults and fault models

The following fault models have been widely used as abstractions of physical defect mechanisms in digital electronics devices/systems [START_REF] Abraham | Fault and error models for vlsi[END_REF], and also applied to neural networks:

• Stuck-at, a data or control line appears to be held exclusively high (stuck-at-1) or low (stuck-at-0).

• Random bit flips, a data or memory element has some incorrect, but random value.

The stuck-at fault model is very popular since many defects at the transistor and interconnection structures can be modeled, as permanent faults, at the logic level with reasonable accuracy. This model is a binary model that does not capture the indeterminate states that faults may induce while occurring. It is applied to model faults on communication channels and arithmetic-logic operators. The random bit-flip model models transient faults that usually happen at registers or memory elements due to external perturbations, for instance, a single event upset. Under this model, damage/corruption is done only to the data and not to the physical circuit itself. Conceptually, it consists of a register bit that is switched randomly, resulting in a memory element that holds a wrong logic value.

B. Fault injection

A fault injection method is required for gaining insights of the behavior of a system [START_REF] Hsueh | Fault injection techniques and tools[END_REF]. Faults are probabilistically introduced into a neural model and the degree of failure, impact on the performed task, is evaluated according to some measures, which basically measure the performance distance (closeness) in performing a task by a fault-free neural network and the derived faulty networks. The measure of fault tolerance from many experiments is evaluated against the number of faults injected into the neural model. The limit of the network fault tolerance is problem-dependent and determined by operating scenarios of multiple faults leading to a violation of the performance constraints. Exhaustive testing of all possible single faults is prohibitive, not to mention that even multiple faults might occur concurrently. Hence, the strategy of randomly testing a small fraction of the total number of possible faults in a network has been adopted for tractability. It yields partial fault tolerance estimates that are statistically very close to those obtained by exhaustive testing.

C. SOM implementation choice

In this paper, we analyze SOM fault tolerance when a map is targeted to a fully parallel digital hardware implementation. We consider a parallel implementation for the SOM recall phase following a general architecture as shown in figure 1. A b-bit fixed-point representation for the internal computations and storage is used, as most SOMs hardware implementations use this number system so as to optimize hardware resources. A neuron essentially consists of registers (to store weights), adders and absolute computation modules that as a whole computes the Manhattan distance. The original SOM uses the Euclidean distance to measure the distance between input patterns and weight vectors, but it requires hardware-greedy multiplications and square roots, thus the Manhattan distance is a preferable choice so that even the distance between the inputs x and the vector prototypes m can be computed concurrently at every clock cycle. The BMU can be implemented as a comparator binary tree and a set of multiplexers, which receive the Manhattan distances and spatial coordinates from all neurons, as shown in figure 2. The BMU's outputs are the spatial coordinates (x, y) of the winner neuron. The neighborhood function plays a key role in SOM learning, but the direct implementation of the exponential function is not practical due to its complexity, instead look-up tables to store precomputed values can be used for resource saving. However, herein we only consider the on-chip implementation of the recall phase and an off-chip implementation for the learning phase.

D. SOM fault tolerance assessment

For SOM fault tolerance assessment, we consider that faults occur only in weight registers in the hardware architecture shown in figure 1. This choice considers plausible hardware faults and not only the unprecise notion of faulty neuron, yet all possible hardware faults are not considered this way (see IV-A). More precisely, we assume that each bit in defective neuron weights, which are stored in b-bit registers, is flipped with some probability and remains the same during the SOM processing. The fault rate u is defined as the percentage of flipped bits (b flips such bits) in neuron weights with respect to the total number of bits in the SOM map. It is used so as to evaluate fault tolerance in SOMs. Considering a w ⇥ w SOM map, with N inputs and b-bit wordlength for weights, yields:

u = b flips w 2 ⇥ N ⇥ b (4)
It must be pointed out that the number of faulty bits and the rate of faulty neurons increase very rapidly with the abovedefined fault rate: to give an example, if b = 16, N = 4, and u = 0.05, 5120 faulty bits appear (out of 102 400) in a 40⇥40 SOM (more precisely, for this case, an ESOM), and 96 % of all neurons contain at least one faulty weight. For space consideration, we have focused on such fault rates as they represent a more devastating effect on the SOM behavior than single, double or triple faults, which can be well masked in SOM as it has been reported in related works [START_REF] Yasunaga | Fault-tolerant self-organizing map implemented by wafer-scale integration[END_REF] [START_REF] Talumassawatdi | Fault immunization concept for self-organizing mapping neural networks[END_REF], and confirmed by our experiments. Moreover, this fault model takes into account the fact that all the introduced faults are not equivalent: faults on least significant bits do not have the same effect as faults on most significant bits.

The severity of faults in the SOM when performing a task can be characterized with respect to the fault rate u: it consists in considering the artificially introduced faults at a given fault rate as the independent variable, and a selected performance measure of the SOM as the dependent variable. As SOM performance criteria for fault tolerance, the average quantization error (AQE) and distortion measure are used, which have been proposed as a way to trust the result of the SOM convergence on a given dataset [START_REF] De Bodt | Statistical tools to assess the reliability of self-organizing maps[END_REF][24] [START_REF] Beaton | Cqoco: A measure for comparative quality of coverage and organization for self-organizing maps[END_REF]. The AQE measures the representation fidelity of the training data, and it is computed as the average distance between each data vector and the prototype of its BMU that should be very close:

AQE = 1 n n X 1 x i m c(xi) 2 ( 5
)
where n is the number of input vectors, x i is an input, m c(xi) is the weight of the BMU in the map for x i . The less the AQE, the better the fidelity of the representation of input data.

The distortion measure considers both vector quantization and structure preservation of the SOM [START_REF] Rynkiewicz | Self-organizing map algorithm and distortion measure[END_REF] and it is suitable for comparing maps of the same size [START_REF] Pölzlbauer | Survey and comparison of quality measures for selforganizing maps[END_REF]. Let's consider a SOM represented by prototypes m j , the distortion measure is defined as the average over all input patterns (a set of n input data vectors x i ) of the weighted sums of the squared distances between the pattern and all neuron prototypes, with weights that depend on the neighborhood function between each neuron and the winner neuron for the given pattern. According to [START_REF] Wu | Spherical self-organizing map using efficient indexed geodesic data structure[END_REF], the original distortion measure would be biased towards the 2D SOM. If the input's distortion is redefined to be the average distortion over its BMU's neighborhood, the results will not be affected by the different topologies of the grid and the distortion measure becomes:

⇠ = 1 n n X i=1 P m j=1 h c(xi)j x i m j 2 P m j=1 h c(xi)j (6)
where c(x i ) is the BMU of x i , and h c(xi)j is the neighborhood function of neurons c(x i ) and j.

V. ENHANCING FAULT TOLERANCE IN SOM

This section briefly introduces three techniques for SOM that have been proposed to improve its fault tolerance: weights restriction, fault insertion, and noise injection during training.

A. Restricting weights during training

The basic idea of this technique is that critical synapses have high magnitude weights [START_REF] Chin | Training techniques to obtain fault-tolerant neural networks[END_REF] [START_REF] Cavalieri | A novel learning algorithm which improves the partial fault tolerance of multilayer neural networks[END_REF]. Restricting the weight magnitudes to be low during training potentially makes the SOM more fault-tolerant. For instance, by i) penalizing high magnitude weights in the learning rule;

m j i (t + 1) 2 [m min , m max ]
, the lower and upper bound of the weight vectors components, respectively, or ii) weight balancing to distribute the weight values more uniformly, e.g., to distribute the absolute weight values around the average absolute value. Here, a weight is updated only if its absolute value, computed by the learning rule, does not exceed the following threshold:

✓(t) = 1 nw 2 X m j n X i=1 m j i (t)

B. Inserting faults during training

This technique randomly inserts faults during training [30] so as to improve fault tolerance. In each training iteration, a fixed number of weights (neurons) are randomly chosen, then faults (bit-flip) are injected according to these steps:

1) Randomly choose a predefined number of weights and insert faults in them 2) Apply the learning rule of on-line (batch) training for a pattern (all patterns) in the training set 3) Restore faulty weights to their unfaulty state and repeat the process until a stop criterion is reached

In our tests, only one faulty bit is introduced (and then restored) for each learning iteration.

C. Noise and weight perturbations

During training, in each iteration, a small number of inputs/weights are selected randomly to be perturbed by adding some type of noise [START_REF] Talumassawatdi | Fault immunization concept for self-organizing mapping neural networks[END_REF] [START_REF] Allende | Robust Self-organizing Maps[END_REF]. The basic idea is to add noise to the inputs or weights of each neuron. Two main options exist, multiplicative and additive noise; additive noise being the most used. We consider such an additive noise in our tests.

VI. EXPERIMENTAL RESULTS

This section presents the experiments that have been carried out to assess SOM fault tolerance. We have focused on studying the performance of 2D maps on artificial datasets, for different bit-flips fault rates in weight registers. We have developed a software environment for fault tolerance evaluation, which handles computations using fixed-point representations (even during learning) as in the proposed hardware implementation. Experiments are limited here to 40 ⇥ 40 2D SOMs, i.e, ESOMs, whose weights are initialized to random values, evenly distributed between [0, 1] and then scaled to a fixedpoint representation Q6.10. As SOM performance depends on many other factors, to reduce the scope of the analysis, all SOMs use the same learning rate, update radius and Gaussian neighborhood function, during 150 epochs with 100 iterations each. Different map initializations are tested, but all learning techniques apply to the same initialized maps and use the same randomly chosen input pattern for each iteration.

A. Datasets

Two simple datasets, called D1 and D2, with 2D points, were used in the experiments. D1 contains uniform observations in [0, 1] ⇥ [0, 1] and D2 observations drawn from a mixture of gaussians: three clusters with a Gaussian distribution, which have equal variance and some overlap. The standard deviation is 1 in each dimension. and the clusters centers are: (0.35, 0.35), (0.65, 0.35), and (0.5, 0.65). An example of the distributions in D1 and D2 are shown in figure 3 a), and b), respectively. It is important to point out that new points for D1 and D2 are generated at each learning iteration according to the chosen distribution, instead of pre-generating a fixed-size dataset. Thus, no cross validation technique is required.

B. Results

The AQE and distortion measure profiles during learning are shown in figure 3. The presented results are averaged over several randomly initialized maps, but all learning techniques and all maps use the same learning patterns, randomly generated at each epoch. It explains why the different curves have similar irregular profiles. Nevertheless, experiments performed with other random patterns provide the same kind of results. Both measures decrease during training for all SOMs, but for D2 the convergence is faster and the measures show less variance for different random initializations. This can be explained by the usual problems known for Kohonen SOMs that unfold themselves in a D1-like input space: some maps spread well while some others remain partially twisted.

Figures 4(a) and 4(d) shows the performance metrics profiles for the SOM learning algorithms with noise injection and thresholding. SOM trained with thresholding yields larger AQE and distortion than the others, and a slower convergence rate. These results tend to show that this technique does not well extend to the kind of fault model we study here. This is even more noticeable for the fault injection based technique: its results are not even included, for readability, in the figure since this method performs worst among the implemented algorithms. This is directly linked to the number representation system that we use herein, where faulty bits can be injected in the 6-bit integer part and induce excessively large weights even after bit restoration (e.g. 000000 becomes 100000 with a faulty bit, then 011111 with learning, and recovering the initial value of the faulty bit is useless here). Using different integer part sizes induces more or less bad results for this method.

Figures 4(b) and 4(e) shows the AQE and distortion measures for D1 and D2 versus the bit-flip percentage introduced into SOMs. Recall that a 1 % fault rate already corresponds to 512 faulty bits in the kind of SOM we consider herein. As the same results are not expected for different faulty SOMs, even on the same dataset, for each fault rate many faulty versions of the trained SOMs were generated and averaged in the presented results. Distortion rapidly increases even with a small percentage of faults, and a criterion (threshold value) should be defined below which the SOM might be considered performing properly. The AQE also increases with the fault percentage, but in a more attenuated way for small percentages. The SOM trained with noise injection slightly outperforms the threshold method, the map is less affected when the number of inserted faults in weights increases. It can also be observed for dataset D1 that the AQE increases slowly at fault rates below 10%, suggesting that faults in weights can be masked in the SOM when it performs quantization, the global internal representation being distributed enough to cover the input space even with faults. However, when an error measure that considers the local behavior of the map is used, results clearly show that the structure preservation is affected even at small fault rate percentages, as the distortion measure indicates. Nevertheless, all these results confirm that SOMs may stand as quite fault-tolerant models, since any abrupt degradation of their behavior with faults is observed. Figures 4(c) and 4(f) show in more detail how fault tolerance evolves through learning for different fault rates for stochastic online learning. The level of fault tolerance is estimated as the ratio between the average (over 100 faulty versions of each map) performance obtained for the given fault rate and the performance without faults. Thus the "optimal" or "reference" value is 1, and a ratio increase means a degradation of fault tolerance. After an immediate decrease of this ratio from the initialization to epoch 1, we see that it increases quite slowly for the AQE measure, and even stabilizes from epoch 100 at a value that depends on the fault rate for D1. For the distortion measure, there is an abrupt increase of the fault tolerance ratio up to very large values, after it stabilizes or even slightly decreases after 50 epochs. Such increase is a direct consequence of the introduction of faulty bits in the integer parts of the weights, as it was the case to explain the bad performance of the fault injection learning technique. The order of magnitude of the ratio goes up to 10 2 10 2 with such faults, since the distortion measure takes into account the weights of all neurons, whereas the AQE measure only considers the winner neuron. Again, it is tightly linked to the size of the integer part, e.g. a Q4.12 instead of Q6.10 fixed point representation would induce approximately 16 times lower values for the maximum distortions. This negative effect on distortion is observed even when introducing only a few faulty bits in our other experiments, whereas the quantizationbased fault tolerance remains optimal for such scarce faults.

To summarize, these results tend to show that:

-Pre-trained SOMs are quite tolerant to faults for vector quantization since their performance gradually degrades with the fault rate. But the specific properties of structure preservation degrade more drastically with faulty bits.

-When observing the evolution of fault tolerance during learning, self-organization induces an improvement with respect to the AQE criterion only during the first learning iterations (initial unfolding of the map). If the distortion criterion is used, even the initial self-organization phenomenon is not able to induce more fault tolerance (figures 4(c) and 4(f)).

VII. CONCLUSION

This paper presents a study of fault tolerance of SOMs. We compare the performances of SOMs obtained with three variants of the on-line training algorithm: weight restriction, fault injection and noise injection. We also study how fault tolerance evolves during learning. The obtained results show that training by inserting noise injection provides slightly better results than the other techniques but it is far from producing highly fault tolerant SOMs. Nevertheless, SOM fault tolerance is globally confirmed since the quantization quality degrades gradually with increasing fault rates. SOM fault tolerance with respect to structure preservation has been studied using the distortion measure, and results from this initial study show that this particular property of SOMs is not necessarily capable of supporting a large number of faults. To address this limit, we consider that a more collective behavior of neurons is required: instead of considering the weights of each neuron separately, a more fault tolerant behavior is expected from a population-coded approach of the prototypes, where single highly faulty weights would not directly affect the global performance. To this end, we currently study how a softwinner-takes-all approach by means of gaussian filtering of prototypes may increase the fault tolerance of SOMs. Then in order to decentralize the process and reduce its dependence to new possibly faulty registers, we intend to replace gaussian filtering by the emergence of bubbles of activity in dynamic neural fields coupled to the SOMs.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Block diagram of an abstract digital hardware implementation for a two-dimensional SOM.

Fig. 3 :

 3 Fig. 3: Quantization error and distortion for SOM stochastic on-line learning for different initialized maps on D1 and D2.

ACKNOWLEDGMENT

The authors would like to thank the partial support received through the Conacyt grant program to conduct research abroad. This work is also partially supported by Erasmus+ 2016-2017 and by the French embassy of Ukraine.

-Among the different learning modifications that have been proposed to increase SOM fault tolerance, only the noise injection technique appears to be robust to the kind of fault model considered here (figures 4(b) and 4(e)).


Benjamin Girault
email: benjamin.girault@usc.edu

Shrikanth S Narayanan

Antonio Ortega
email: antonio.ortega@sipi.usc.edu

Local stationarity of graph signals: insights and experiments

Keywords: Graph signal processing, stationary graph signals, graph learning

In this paper, we look at one of the most crucial ingredient to graph signal processing: the graph. By taking a step back on the conventional approach using Gaussian weights, we are able to obtain a better spectral representation of a stochastic graph signal. Our approach focuses on learning the weights of the graphs, thus enabling better richness in the structure by incorporating both the distance and the local structure into the weights. Our results show that the graph power spectrum we obtain is closer to what we expect, and stationarity is better preserved when going from a continuous signal to its sampled counterpart on the graph. We further validate the approach on a real weather dataset.

INTRODUCTION

Graph Signal Processing (GSP) is thriving as the tool of choice for studying data lying on arbitrary irregular structures. [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF][START_REF] Sandryhaila | Discrete Signal Processing on Graphs[END_REF] By exploiting these structures, we have been able to provide reasonable extensions to classical definitions of the signal processing toolbox, including decompositions such as the graph Fourier transform, wavelet decompositions, [START_REF] Hammond | Wavelets on graphs via spectral graph theory[END_REF] filterbanks, [START_REF] Narang | Compact Support Biorthogonal Wavelet Filterbanks for Arbitrary Undirected Graphs[END_REF] or graph empirical mode decomposition, [START_REF] Tremblay | Graph Wavelets for Multiscale Community Mining[END_REF] but also tools to study the properties of said transforms such as the uncertainty principle, [START_REF] Agaskar | A Spectral Graph Uncertainty Principle[END_REF] or more recently our proposed model for stochastic graph signals: graph stationarity. [START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF][START_REF] Girault | Signal Processing on Graphs -Contributions to an Emerging Field[END_REF] In most of the literature, the graph structure is supposed to be known. Yet, the question of choosing the right graph is of utmost importance and often overlooked. Since graph signal processing is the study of signals in the context of a given graph, then working on a graph that is not adequate with respect to the actual structure of the data leads to a study whose conclusions may not be relevant. In this paper we tackle this problem with a data driven approach, and focus on providing evidence that some conventional approaches for selecting graph structures may not be ideal, even though they are widely used in the literature.

More specifically, we refer here to graphs constructed to represent data obtained by sampling a Euclidean space. The conventional approach uses a Gaussian kernel of the Euclidean distance to assign weights to the edges. Since working with complete graphs is expensive, the next step is often to sparsify the graphs: Those weights that correspond to distant nodes are removed to have reasonable algorithm complexity. This leads to a sparse graph whose vertices have a small number of neighbors compared to the total number of vertices. This setting corresponds for example to sensor network applications that appear frequently in the graph signal processing literature, and is a well studied problem from a mathematical perspective in the context of manifold sampling. [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF] Yet, theoretical results only guarantee that graph signal smoothness is related to the underlying continuous signal smoothness: The fact that the continuous Laplacian can be approximated by the graph Laplacian when applied to a signal does not guarantee that the graph Fourier transform can approximate the continuous Fourier transform, or that graph frequencies correspond to continuous frequencies. Consequently, sampling a continuous stationary signal leads to a non-stationary graph signal hinted by spectral correlations.

The further step of approximating the graph by removing those edges with small weights to make it sparse does not help either with regards to extending the Fourier transform. Observations show that this step can actually change the properties of the graph quite drastically. Ultimately, such a method to define the graph relies on trade-offs between sparsity and accuracy.

We propose here an alternative approach involving learning a graph structure from data to overcome issues arising with these conventional methods. We study this approach on two typical cases where the common question is that of finding a graph that gives a relevant spectral model for the data. The first case is that of approximating the continuous Fourier modes: Can we model continuous frequencies on a Euclidean space using an irregular sampling of said space? The second case involves studying real data, for which distances are only one of the explanation to correlations, and we answer the question of learning a relevant graph to capture a good statistical model for the data.

Several methods exist in the literature to learn a graph structure from data. [START_REF] Egilmez | Graph Learning from Data under Laplacian and Structural Constraints[END_REF][START_REF] Segarra | Network Topology Inference from Spectral Templates[END_REF][START_REF] Pasdeloup | Characterization and Inference of Graph Diffusion Processes from Observations of Stationary Signals[END_REF] We choose here the approach of Egilmez et al. [START_REF] Egilmez | Graph Learning from Data under Laplacian and Structural Constraints[END_REF] for its ability to constraint the connectivity of the graph to any set of edges and learn their weights. Doing so, we show that the learnt graphs lead to smoother graph Fourier modes in the sense of the Euclidean domain, sparse graphs, and, more importantly, better stochastic models for the data in the sense of smoother graph power spectrum and more interpretable spectral components.

The rest of this paper is organized as follows. Section 2 introduces the necessary background on graph signal processing. Section 3 recalls the classical definition of stationarity and our definitions of (global) graph stationarity [START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF] and local graph stationarity. [START_REF] Girault | Towards a definition of local stationarity for graph signals[END_REF] Section 4 describes our approach and that of the literature, and both Sec. 5 and Sec. 6 applies our approach to an ideal case and a real case.

GRAPH SIGNAL PROCESSING

Let G = (V, E, w) be a graph with V its set of vertices, E ⊆ V × V its set of edges, and w : E → R + its edges' weight function. We assume here that the graph is undirected, i.e., whenever ij is an edge from vertex i to j, ji is also an edge and w(ij) = w(ji). Let A be the adjacency (or weight) matrix defined as A ij = w(ij) whenever ij is an edge, and A ij = 0 otherwise. Note that since the graph is symmetric, A is also symmetric.

At the core of Graph Signal Processing (GSP), lies harmonic analysis of graph signals. Graph signals are functions x : V → R or C. Assuming a numbering of the vertices from 1 to N = |V|, x(i) will be denoted x i such that x can be represented by the vector x = [x 1 , . . . , x N ] T . Everything that follows is independent of the numbering chosen. The most important tool for harmonic analysis is the so-called Graph Fourier Transform (GFT) defined as the projection onto a basis of graph signals called the graph Fourier modes. Definition of these modes varies. In this paper, the only requirement is to have a unitary transform. Let F be the matrix of the GFT. Then

F -1 = [u 0 • • • u N -1 ] = F *
is the inverse GFT, where the last equality follows from unitarity and . * is the complex conjugate operator. We obtain the following operations on a signal x:

x(l) = x, u l = u * l x = N i=1 u * l (i)x i x i = N -1 l=0 x(l)u l (i). (1)
The second requirement is to have a set of graph frequencies associated to the set of graph Fourier modes such that the lowest frequency is 0, and their magnitude is upper bounded by 1. We denote λ l those frequencies, with λ 0 = 0 associated to the lowest frequency u 0 . Let Λ = diag(λ 0 , . . . , λ N -1 ) be the diagonal matrix of graph frequencies. From now on, we denote Q = F * ΛF the fundamental matrix of the graph * .

One example of modes is the set of eigenvectors of the Laplacian matrix L = D -A, where D is the diagonal degree matrix verifying D ii = j A ij , associated to the eigenvalues of L, normalized (the exact normalization is not a sensitive parameter here and the interested reader can find one our previous work [START_REF] Girault | Translation on Graphs: An Isometric Shift Operator[END_REF] ), as graph frequencies. [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF] We obtain then Q = 1 α L, with α this normalization factor. Other definitions of graph Fourier modes exist, 2 however, the rest of this paper is not dependent on any specific GFT.

The second harmonic analysis tool we are using in this paper is a localized decomposition. Let {g i,m } i,m be a family of atoms (i.e., graph signals) such that g i,m is localized on vertex i, and m indexes the span of the atom. Additionally, we supppose these atoms to be of unit norm: g i,m 2 = 1. The decomposition is then defined as a projection onto this family similarly to the GFT: W x (i, m) = x, g i,m . One such family is the Spectral Graph Wavelet Transform 3 (SGWT) defined using the spectral properties of the atoms. More precisely, g i,m is the localization around vertex i of m th scale kernel g m defining these spectral properties:

g m (l) = exp(-κs m λ l ) g i,m = c i,m T i g m , (2)
where {s m } m is a set of scales and T i is the localization operator verifying T i x(l) = δ i (l) x(l). This operator has been shown 3 to effectively localize a smooth kernel such as g m above. In this particular case of an exponential kernel, it is shown that there is also an exponential decay of the energy of an atom g i,m with the distance (on the graph) to vertex i, validating the concept of localized atom, and decomposition.

STATIONARITY

In this section we recall a classical model for stochastic signals, the stationarity model. We do so both with a Euclidean setting, and with a graph setting.

Stationarity for Euclidean Spaces

In the rest of this paper, we denote with bold character stochastic signals. Let s be a stochastic signal. Wide Sense Stationarity (WSS) is then defined for Euclidean spaces as: [START_REF] Cressie | Statistics for spatial data[END_REF] ∀x,

E [s(x)] = µ s ∀x 1 , x 2 , E [s(x 1 )s * (x 2 )] = γ s (x 1 -x 2 ). (3)
In other words, the first moment is constant, and the second moment depends only on the coordinate difference x 2 -x 1 . An equivalent definition involves the statistical invariance through translation:

∀x, dx, E [s(x)] = E [s(x + dx)] ∀x 1 , x 2 , dx, E [s(x 1 )s * (x 2 )] = E [s(x 1 + dx)s * (x 2 + dx)] . (4)
Although equivalent for Euclidean spaces, they are not in a graph setting as we will be discussed in Sec. 3.2: Only Eq. ( 4) can be directly transposed in a graph setting while Eq. ( 3) is more complex to transpose due to the coordinate difference

x 1 -x 2 .
WSS is a simplification of the Strict Sense Stationarity (SSS) definition where a signal is SSS if and only if it is statistically the same after translation. In other words, all the statistical moments of a given signal are equal, making it a particularly difficult, if not impossible, property to verify.

Intuitively, stationarity is the property that the signal does not depend on the origin of the system of coordinates: any changes of this origin leaves the signal unaffected. Conversely, a non-stationary signal is one that exhibit local properties whose location would be translated if the origin is modified.

A Global Definition for Graphs

Global stationarity, or stationarity in short, is obtained by directly extending the concept of statistical invariance to the graph setting: In a Euclidean setting, H is a shifting operator, e.g. the time shift T {s}(t) = s(t -1) for temporal signals or a combination of the spatial shifts T x {s}(x, y) = s(x -1, y) and T y {s}(x, y) = s(x, y -1) for a 2D signal such as in image † . Unfortunately, the notion of translation on a graph is not well defined, since there is no immediate system of coordinates, or any composition law on the set of vertices (if i and j are vertices, such a law would map them to a new vertex i • j, similarly to x + y being itself a new coordinate), for which a translation of the signal is equivalent to a translation of the system of coordinates.

We proposed recently a new operator [START_REF] Girault | Signal Processing on Graphs -Contributions to an Emerging Field[END_REF][START_REF] Girault | Translation on Graphs: An Isometric Shift Operator[END_REF] on graph signals whose behavior can be compared to that of a diffusion operator, [START_REF] Girault | Localization Bounds for the Graph Translation[END_REF] the Graph Translation:

T G = exp (-ıπQ) . ( 5
)
Contrary to Q, this operator is unitary: It does not alter the energy of the signal. This property is a natural assumption for an invariance operator since the second moment cannot be invariant if the energy changes. Finally, diffusion invariance is a natural parallel to shift invariance since it is the primary use of graphs and graph signals.

We introduced the graph WSS definition using WSI through the graph translation:

µ s = µ T G s R s = R T G s ,
Note that the requirement on the correlation matrix can equivalently be replaced by the same requirement on the covariance matrix Σ s . However, due to the fact that the first moment needs to be estimated prior to estimating the covariance matrix, we will favor the correlation matrix R s when studying real data. Finally, using the GFT on the signal s, we obtain a spectral characterization of WSS graph signals: [START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF] Theorem 3.3 (Graph WSS Spectral Characterization). A graph signal s is WSS if and only if:

∀λ l = 0, µ s (l) = 0 ∀λ l = λ k , S s (l, k) = 0,
where S s = R s is the spectral correlation matrix. In other words, graph frequencies are uncorrelated, and the mean is the lowest frequency of the signal (a DC component). This also gives a formal relation between the correlation matrix R s and the covariance matrix Σ s :

Σ s = R s -µ s µ * s = R s -µ s (0)u 0 u * 0 , using the fact that Σ s = S s -µ s (0)E 1,1
, where E 1,1 is the matrix with zeros everywhere except for the top left element. From Thm. 3.3, we can define the graph Power Spectrum Density as the variance of each spectral component:

γ s (l) = E | s(l)| 2 = S s (l). (6)
In turn, we can write the elements of the correlation matrix using the PSD:

R s (i, j) = (T j γ s ) (i). (7)
Finally, two alternative extensions of stationarity for stochastic graph signals exist in the litterature following our inital work. [START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF] These are summarize on Tab. 1. [START_REF] Perraudin | Stationary Signal Processing on Graphs[END_REF] This work assumes that the graph power spectrum of a WSS graph signal can be written as a function of the graph frequency. These graph signals are also WSS according to our definition (the converse is not always true).

Perraudin and Vandergheynst approach

Paper Q 1 st moment 2 nd moment Constraints [7] Q µ s = µ T s E[ss * ] = E[(T s)(T s) * ], T = exp(-ıπQ) Q normal [17] L µ s = α1 Σ s (i, j) = (T j γ s ) (i), γ s (l) = E[| s(l)| 2 ] no multiple eigenvalues, γ s (l) = f (λ l ) [18] Q µ s = 0 E[(Q a s)((Q * ) b )s) * ] = E[(Q a+c s)((Q * ) b-c )s) * ] Q normal
Marques et al. approach [START_REF] Marques | Stationary Graph Processes and Spectral Estimation[END_REF] Here, the invariance is more involved and is with respect to the fundamental matrix Q of the graph which is not an isometric operator. As a consequence, the first moment is assumed null.

The invariance is actually equivalent to Q and R s commuting (QR s = R s Q) which is itself equivalent to our second moment assumption. These graph signals are also WSS according to our definition (the converse is not true because of the first moment requirement).

A Local Definition for Graphs

Following our work on stationarity for graph signals, and inspired by the work of Perraudin and Vandergheynst, [START_REF] Perraudin | Stationary Signal Processing on Graphs[END_REF] we introduced recently a notion of local stationarity. [START_REF] Girault | Towards a definition of local stationarity for graph signals[END_REF] Given a set of localized atoms {g i,m } i,m , the local power spectrum is given by:

S s (i, m) := E |W s (i, m)| 2 = E | s, g i,m | 2 . ( 8
)
Local stationarity is then the property that no matter where we look on the graph, the signal is the same. Using the local PSD above, this translates into:

∀i, j, m, S s (i, m) = S s (j, m)
Contrary to the global definition in Sec. 3.2, this definition involves the invariance of local descriptors of the signal, hence the name of local stationarity.

GRAPHS FOR EUCLIDEAN SPACE SAMPLINGS

In this section we explore an often overlooked aspect of graph signal processing: the selection of the right graph.

In particular, we focus here on sampling Euclidean spaces. These are for example sensor networks, or feature spaces in machine learning. We recall the usual structure used, and then describe our approach using stationary graph signals and graph learning to attain the goal of a better statistical model for the data.

State of the Art

Given a set of measurement point {x i } i∈V of a Euclidean space S, how do we connect them into a graph structure? Literature on Graph Signal Processing have been using a standard result from differential geometry 9 to justify the use of a Gaussian kernel of the distance as edge weights:

w(ij) = exp - x i -x j 2 2σ 2 , ( 9
)
for some parameter σ. The result in question is stated in Thm. 4.1.

Theorem 4.1. [9, Theorem 3.1] Let s : S → R be a continuous signal and ∆ the Laplace-Beltrami operator of S. Let G (N ) be a series of graphs constructed by uniformly sampling S, with all edges weighted by a Gaussian kernel of the Euclidean distance. Let L (N ) be the associated graph Laplacians. Then, under some assumptions on σ, and rescaling of the Laplacian, we have the following convergence result for the sampled signals s (N ) of s on the graphs:

L (N ) s (N ) ----→ N →∞ ∆ s.
Appplication of this result in the context of graph signal processing is challenging since this result only gives guarantees for an arbitrary continous signal s. Yet, the ideal convergence result would be that the graph Fourier modes converge to the continuous Fourier modes. Since those are defined as eigenvectors of the Laplacian, rather than sampling of a continuous Fourier mode, this paper [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF] does not give any guarantees with that regard.

As a consequence, the graph spectrum of a graph signal does not inform on the continuous spectrum of the underlying sampled continuous signal. The only guarantee is that given a signal s on the graph, and its spectral representation s, we have that s * Ls = l | s(l)| 2 λ l is a good representative of the smoothness s, ∆ s of the underlying continuous signal s. In particular, for a continuous Fourier mode e and e the sampled graph signal, e * Le = l | e(l)| 2 λ l is a good representative of the actual frequency, but e is not necessarily localized in the graph spectral domain.

The second question with this result is the value of σ in Eq. ( 9). Using the value used in [9, Theorem 3.1] leads to a value so small that the graph is almost disconnected: Most of the edges have weight very close to 0. The question then is how to choose this crucial parameter. The convention is to just use one value that gives interesting result, not too high to avoid all weights being close to 1, and not too small to avoid having almost all weights close to 0. One heuristic states that a good choice of σ is a third of the mean distance between vertices. [START_REF] Chapelle | Semi-Supervised Learning[END_REF] We adopt this heuristic to compare the Gaussian weights approach to ours as it gives reasonable results. Furthermore, this complete graph is rarely used without sparsifying it first. In that case, those edges between distant vertices are removed. This step is often desirable for complexity reasons of GSP algorithms: Sparsity of the graph leads to sparsity of matrices and efficient matrix-vector or matrix-matrix multiplications. We investigate two approaches here: Keep the edge between vertices in a given radius, denoted THR graph, or keep the K nearest vertices, denoted KNN for K-nearest neighbors. The complete graph is denoted CMP.

The result of Thm. 4.1 does not include the sparsification case. As a consequence, there is no guarantee that those sparse graphs are good representations of the underlying Euclidean space.

Gaussian Random Fields

In this paper, we propose a data-driven approach to define the graph. More precisely, our approach aims at following more closely the result of Thm. 4.1: we define the graph from a given continuous signal. In this section, we describe a continuous signal model that is a good candidate to infer the graph from, and its corresponding graph signal model. We consider here only zero mean signals, but this can easily be extended to arbitrary means.

Continuous model From now on, bold characters designate stochastic quantities. Let s(x) be a D-dimensional continuous stochastic signal with x ∈ R D and:

s(x) = γ s (ν) cos ı2π(x T ν + φ(ν)) dν,
for some prescribed Power Spectrum Density (PSD) γ s (ν) of frequency parameter ν ∈ R D , and φ(ν) i.i.d. phases uniformly drawn in [0, 1]. s(x, y) is called a Gaussian Random Field (GRF). Such a phase randomization is known to stationarize a non-stationary continuous signal, and is used to test stationarity using surrogates. [START_REF] Borgnat | Testing stationarity with surrogates: a time-frequency approach[END_REF] s is therefore a continuous WSS signal. We can compute its autocorrelation function:

R s (x 1 , x 2 ) = E[s(x 1 )s * (x 2 )] = γ s (ν) cos ı2π(x 1 -x 2 ) T ν dν = R s (x 1 -x 2 ), (10)
and its spectral autocorrelation function using the Fourier transform of the autocorrelation function: 8, 21

S s (ν) = F { → R s (δx)} = γ s (ν). (11)
Graph signal model The equivalent of this model in a graph setting is called a Gaussian Markov Random Field (GMRF). A signal s is GMRF if it follows the distribution N (0, Σ), a multivariate Gaussian distribution with zero mean, and covariance Σ. When the covariance matrix is jointly diagonalizable with the fundamental matrix of the graph, a condition which is met when the graph Fourier matrix diagonalize Σ, then a GMRF model can be described in the graph spectral domain as:

s(i) = N -1 l=0 s(l)u l (i),
where the variables s(l) ∼ N (0, γ s (l)) independent variables. Therefore, s is graph WSS. Its spectral autocorrelation matrix is a diagonal matrix verifying:

S s (l, l) = γ s (l). ( 12
)

A Data-Driven Approach to Choosing a Graph

Observing both Eq. ( 11) and Eq. ( 12) reveals an interesting parallel: Both spectral correlations are equal to the power spectrum. Starting from this observation, we approach the problem of choosing a graph by learning it such that the continuous PSD and the graph PSD coincide as functions of the frequencies. To do that, instead of using Gaussian weights with a somewhat arbitrary parameter σ, we learn these weights using a GRF model sampled on the set of vertices. Equation (10) gives then the correlation matrix of the stochastic graph signal obtained using the coordinates of the vertices.

To learn the weights, we use the approach of Egilmez et al. [START_REF] Egilmez | Graph Learning from Data under Laplacian and Structural Constraints[END_REF] In this work, the Laplacian matrix is optimized to make it coincide with the precision matrix Since Thm. 4.1 shows a parallel between the continuous Laplacian and the graph Laplacian, and since the eigenvalues of the continuous Laplacian are squared frequencies, the graph equivalent to continuous the PSD γ(ν) = 1/ν 2 is the graph PSD γ(λ) = 1/λ. Therefore, learning the graph from a covariance matrix is here equivalent to considering the data as a continuous GRF signal s with PSD γ s (ν) = 1/ν 2 .

Ω s = Σ † s = R † s of
We propose here to learn the graph structure from a Euclidean GRF model in 1/ν 2 . In Sec. 5, we consider the ideal case for a two dimensional Euclidean plane with GRF model in 1/ν 2 , and in Sec. 6, we look at a real case, where the exact geometry of the underlying space is unknown, and for which we assume the data to be sampled from such a GRF model on this unknown space.

Finally, to solve the question of learning a sparse graph, we use the conventional approach to select vertices. In other words, in this paper, we use K-nearest neighbors or distance thresholding to limit the number of edges of the graph, and perform the learning algorithm on those weights alone. The goal is to take into account the distance, and counterweight the particular local structure of the graph by learning the graph weights. In Sec. 5 and Sec. 6 we are using the same set of edges to better compare the Gaussian and the learnt approaches.

EXPERIMENT WITH A GRF MODEL

In this section, we look at the ideal two dimensional Euclidean space with a GRF model g whose PSD is γ g (ν x , 0) = 1/ν 2

x and γ g (0, ν y ) = 1/ν 2 y for ν x , ν y ∈ {1, . . . , N/2}, with N the number of sampled vertices in the space. We compare several properties of the graph Fourier transform and of the power spectrum of the graph signal sampled from this GRF model, using six different graphs. The first group of graphs is weighted by a Gaussian kernel, and the second is weighted by our learning approach. Each group is a set of three graphs that differ only by the sparsity: complete (CMP), K-nearest neighbor (KNN), and thresholded distance (THR). In this section, the graph has N = 200 vertices, K = 6, and the radius for THR is set to have the same mean number of neighbors for both KNN and THR.

Graph Fourier Transform

Studying the spectral properties of the graph signal is relevant only if the graph Fourier modes do correspond to meaningful modes, that is, if we were to extrapolate them on the underlying Euclidean space, then the graph frequency should be a good indication of the continuous frequency of the extrapolated signal.

Qualitatively, we show on Fig. 1 several graph Fourier modes along with their closest continuous Fourier mode (leftmost column). Two remarks are in order: (i) the sparsified graphs Gauss KNN and Gauss THR show strong clustering effect (almost piecewise constant), and (ii) the graph Learnt CMP has Fourier modes that do not look like continuous Fourier modes.

The first point is important here as the number of edges removed can have a high impact on the graph Fourier modes. Here the mean number of neighbors (around 7) can often be found in the literature, but we see that drawing any conclusion on the graph Fourier spectrum of the sampled signal is difficult. The same remark applies for small values of σ (see Sec. 5.4). The second point can be interpreted as overfitting of the learning approach since the number of parameters (edge weights) is close to the size of the data (covariance matrix). We can further study these issues quantitatively using the graph Laplacian variation on Gauss CMP. We know [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF] that given a continuous signal s, the output Ls of the graph Laplacian L applied to the sampling s of s converges exponentially as the number of vertices increases to the output ∆ s of the continuous Laplacian ∆ applied to s. Consequently, s * Ls is a good proxy for the quantity s, ∆ s . Figure 1 shows the results for the various graph Fourier modes depicted. These numbers confirm the qualitative remarks we made above. In other words, if working on Gauss CMP can be afforded in terms of complexity and a good value of σ is found, then it is a good choice, otherwise, the graphs Learnt KNN and Learnt THR are also good choices.

Original Data

Spectral Correlations

The goal being that the sampled continuous GRF signal is graph WSS, we look at the correlation matrix R g of the sampled graph signal g of the GRF signal g, and the more informative matrix of Pearson correlations:

C g (k, l) = S g (k, l) S g (k, k) S g (l, l) .
This matrix removes the influence of the variance on each mode and shows how much correlation there is between any pair of graph frequencies. Figure 2 shows this matrix for the sampled continous GRF graph signal g using different graphs. From this figure, we observe that (i) for all graphs a large amount of correlation is given by the power spectrum (there less is off-diagonal energy than for the direct domain correlation matrix), and (ii) Gauss CMP has more correlations between frequencies than the other graphs.

With respect to spectral correlations, the graph Gauss CMP should be avoided, and the learnt graphs give good results, especially the sparse ones.

Graph Power Spectrum

In this section, we look at how the graph power spectrum varies with the graph frequency, compared to the underlying continuous PSD of the GRF model g. Recall that this PSD is a 1/ν 2 model, such that we are looking for a 1/λ graph PSD. On the other hand, the Learnt graphs follow more closely the target of b = 1, with less noise on the PSD. However, the graph Learnt CMP has a different behavior in the lower spectrum, with b closer to 6. This shows that the overfitting behavior shown by the GFT is at play here, and lower graph frequency modes are actually poluted by higher continuous frequencies. Consequently energy of the higher continuous frequencies is shifted to the lower spectrum, and a higher value of b is obtained.

In terms of graph selection, the learnt sparsified graphs work best for the graph power spectrum of the graph signal g to match that of the underlying continuous GRF model.

Influence of σ

Before studying the local power spectrum of the graph signal, we focus on one crucial parameter of the Gaussian weights approach: the parameter σ. On Fig. 4(a), we study its influence on spectral correlation (i.e., the energy of the off-diagonal elements in the matrices of Fig. 2). The figure shows one important point: There is an optimal value to minimize spectral correlations, around σ = 0.032 (versus σ CMP = 0.173, σ KNN = 0.0255 and σ THR = 0.0242 when using the heuristic of one third of the mean distance used on Fig. 2).

There are however arguments against using such a method to set the value of σ: The GFT for this small value of σ is very close to that of Gauss KNN and Gauss THR on Fig. 1. More precisely, decreasing the value of σ tends to disconnect clusters as many more edges will have small weights, and as such the low frequency Fourier modes are forced to be almost constant on those clusters. This is hinted by Fig. 4(b) where the weights sharply decay with the distance when using the optimal σ, such that some edges in Gauss THR are almost removed, those between clusters. As a consequence, the graph Fourier modes are not easily interpretable in the underlying Euclidean space, from which the signal g is sampled.

In other words, Fig. 4 shows that it is possible to minimize the mean spectral correlation (and on this example it is equal to mean spectral correlation of Learnt KNN and Learnt THR), but doing so is at the cost of a meaningful GFT (not shown here, but included in the matlab code associated to this paper). The learnt sparse graphs on the other hand show the best of both worlds: meaningful graph Fourier modes, and small spectral correlations. Learning the graph allows here to have weights that depend on both the distance and the neighborhoods to overcome the impact of the irregular Euclidean space sampling. This is illustrated on Fig. 4(b) with a non-linear relation between weight and distance.

Graph Local Power Spectrum

To compute the local power spectrum, we use the spectral graph wavelet transform described in Sec. 2. The exact parameters vary from one graph to another to account for the difference of eigenvalue repartitions (see Fig. 5), and can be found in the source code.

Notably, the learnt graphs have several very high eigenvalues. These correspond to those vertices that are extremely close to each other, making them highly correlated, hence being given a high weight to the edge connecting them. Although these high weights could be attenuated using regularization in the learning algorithm, they do represent an important feature of the sampling: If a given signal have a high variation between two such vertices, then there is also a high variation in the underlying continuous signal, hence a high frequency. In other words, these eigenvalues do correspond to very high frequencies in the continuous domain.

However, for the SGWT that relies on the largest eigenvalue to define the wavelets filters, that means that the spectral responses of the filters are skewed towards higher graph frequencies where there are actually few of them. To circumvent this we parametrized differently the span of the wavelet functions in the SGWT from one graph to another. We will further investigate this point in a future communication.

In terms of local stationarity, we observe on Fig. 6(a) that not all vertices have the same local power spectrum, for any graph. However, there is a link between the mean (over vertices) of the local power spectrum and the global spectrum as shown in Fig. 6(b) and studied before. [START_REF] Girault | Towards a definition of local stationarity for graph signals[END_REF] We can therefore study the global power spectrum from the local power spectrum, but deciding global stationarity from it is not possible, even with the exact local PSD matrix.

Using the SGWT, as shown previously, [START_REF] Girault | Towards a definition of local stationarity for graph signals[END_REF] we can infer the graph PSD. However, local and global stationarity are not equivalent definitions, and one does not imply the other, as shown here. One important remark though: The graph used is here important as the SGWT depends on the GFT and the graph frequencies. We will further investigate the relationship between the SGWT and the graph in a future communication to better understand local stationarity.

GRAPHS FOR REAL SIGNALS

We now look at experimental validation of this data-driven approach of selecting a graph. In a recent communication, [START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF] we used a sensor network dataset to study stationarity and used the graph Gauss CMP to study stationarity of the sensor readings. We propose to refine the study of this dataset in two ways: (i) we perform a different, less invoved, preprocessing, and (ii) we look at graph stationarity with respect to different graph constructions to select one that best describe the statistics of the data. More precisly, this dataset is called the Molène dataset, [START_REF]Données horaires des 55 stations terrestres de la zone Large Molène sur un mois[END_REF] from the name of an island off the coast of Brittany, France. As the name suggests, the sensor network is spread in the region of Brittany, and we study the temperature readings. Those readings are performed hourly, during the month of January 2014. This makes for 744 readings. Each of the 37 weather stations is associated with a geographical position such that construction of the Gaussian kernel is possible to build the Gaussian graphs. K for Gauss KNN is chosen equal to 6, and the radius for Gauss THR is chosen to obtain the same mean number of neighbors in both sparsified graphs. , (e 0.2 , 0.4), (e 0.2 , 0.3), (e 1.5 , 1), (e 0.7 , 0.6), (e 0.7 , 0.6). The yellow dash-dotted line corresponds to the case b = 1.0. Note that the lowest frequency 0 is not shown here because of the logarithmic scale on the horizontal axis.

Using the sensor readings as many realizations of the same stochastic graph signal, we can compute sample estimators of the first two moments. The goal is to ultimately give a stochastic model for this dataset.

Data Preprocessing

As shown before, [START_REF] Girault | Stationary Graph Signals using an Isometric Graph Translation[END_REF] preprocessing of the data is necessary because of obvious non-ergodicity shown by the data. Examples of those include cyclo-stationarities (day-night cycles), and trend over the month of the temperatures. We previously used the empirical mode decomposition as preprocessing step, a method empirically known to remove these issues for a more sensible ergodicity assumption. However, this method is particularly involved, and as its name suggests, it is empirical. In, [START_REF] Perraudin | Stationary Signal Processing on Graphs[END_REF] the authors used a simple mean removal from the data, but this method is not ideal since it has to be done twice making for 37 + 744 mean estimations: one for each time instant to remove the noise, and one per vertex to cancel out the first moment.

We propose here to model this data using a different assumption: temperatures follow a Brownian motion model. In other words, we assume that temperature increments are independent Gaussian variables. This model as the advantages of not performing any estimation (as opposed to mean removal), of an easy intuition (as opposed to an EMD), and is reasonable in a graph signal processing setting: Correlations between distant vertices are much smaller with temperature increments. One explanation for this phenomenon comes from the study of those stations highly correlated yet very far apart: They are on the coast line such that the sea and its rather constant temperature is attenuating the temperature variation due to the air ‡ . Further analysis reveals a first moment close to 0 on all vertices (maximum amplitude of the first moment is 2.5 × 10 -3 , with temperature increments in [-8°C, 8°C]). The first moment condition of stationarity is then verified. In the rest of this section we study spatial correlations of temperature increments.

Learnt Graph

We use the same data-driven approach to build the learnt graphs. However, instead of having to resort to synthetic graph signals, we here use the sensor readings: the temperature increments. We obtain a covariance matrix fed to the learning algorithm. [START_REF] Egilmez | Graph Learning from Data under Laplacian and Structural Constraints[END_REF] Figure 7 shows then some low frequency graph Fourier modes of the Gaussian and learnt graphs, along with the Laplacian variations with respect to Gauss CMP of those modes. The same conclusion on Gauss KNN and Gauss THR can be drawn here: There is a noticeable clustering effect. Learnt CMP is again less smooth than one would expect, and the learnt sparse graphs achieve good smoothness, making them the best choice to have good Fourier modes and sparsity.

Graph Stationarity and Graph PSD

Figure 8 shows the spectral Pearson correlation matrices for the different graphs. Results are again very good for the learnt sparse graphs, especially for Learnt THR that shows small spectral correlations in average, but also the smaller extreme correlations. Noteworthly, all graphs are able to substantially decrease the dimension from a dense correlation matrix in the vertex domain, to a highly sparse spectral correlation matrix. ‡ On this dataset one could very well work with the original temperatures and learn a graph, but the sea would need to be taken into account. Finally, Fig. 9 shows the associated graph PSD for the different graphs. First of all, except for the graph Learnt CMP, no graph PSD follows the model in 1/λ. However, these are data real for which we do not expect the graph PSD to follow this model. Note also that overfitting is at play on Learnt CMP as the graph learning approach will tend to find the graph that best matches this graph signal to a model in 1/λ. The second remark is on the fourth graph Fourier mode for the graphs Gauss KNN and Gauss THR (third data point on Fig. 9): This component is very small. Looking at Fig. 7 we see that these modes have high energy on a vertex corresponding to an island station. The temperature of this station should however not be very different than those of the close coastline, such that this mode is absent in the temperature increment graph signal, leading to a small energy. Here the clustering effect goes against the real data, and the graph PSD is not representative of the real data. Finally, the learnt sparse graph achieve good smoothness of the graph PSD, one for which graph filtering using low order polynomial is reasonable. [START_REF] Shuman | The Emerging Field of Signal Processing on Graphs: Extending High-Dimensional Data Analysis to Networks and Other Irregular Domains[END_REF]

Local Power Spectrum

To conclude this experimental section, we look at the local power spectrum on Fig. 11. The associated frequency spectrum of the SGWT is shown on Fig. 10. First of all, for no graph, the signal is local stationary as the local power spectrum varies from one vertex to another. As pointed out in a previous communication, [START_REF] Girault | Towards a definition of local stationarity for graph signals[END_REF] local stationarity is a property that is very hard to verify due to the high number of parameters to the model (number of scales times the number of vertices) compared to the number of samples (vertices). Deciding local stationarity could however be done statistically, and we will further investigate this crucial question in a future communication.

On this other hand, the mean local power spectrum on Fig. 11(b) is smooth for all graphs and can be used to estimate the global power spectrum. Notice that the scaling function coefficient (first data point) is associated to a very high energy. This is due to the fact that the data has some DC component, although it is zero in average. By construction of the SGWT atoms, only the scaling function is affected by this energy.

Gauss CMP

CONCLUSION

For sampled Euclidean spaces, we showed that the usual choice of graph weights deriving from a Gaussian kernel of the Euclidean distance is not ideal, especially with respect to the graph Fourier transform. This is true for both a complete and sparsified graphs. We showed then that learning the weights allows for more flexibility as those weights can depend not just on the distance between vertices but also on the connectivity patterns that are generated by the Euclidean space sampling. We experimentally obtained better results using our learning approach with respect to stationarity, power spectrum, and even the graph Fourier modes smoothness in the lower spectrum.

We illustrated this approach on a real dataset involving temperature measurements. With our approach, we obtained more sensible graph Fourier modes, and a notion of stationarity on a sparse graph that make sense with respect to the data. This is further validated with a smooth graph power spectrum. Preliminary results on the local power spectrum hint at a local decomposition that is interesting to obtain an estimate of the global power spectrum using the SGWT. However, sensitivity of the SGWT parameters to the distribution of eigenvalues, leads to the question of the relation between the two. We will look further into this question, and the broader question of the localized set of atoms in a future communication.

Definition 3 . 1 (

 31 Weak Statistical Invariance (WSI)). A stochastic (graph or Euclidean) signal s is said WSI with respect to the operator H if and only if its first two statistical moments are invariant through H: µ s = µ Hs R s = R Hs , with µ s = E [s] the first moment ( mean) and R s = E [ss * ] the second moment ( correlation matrix).

  the data, with . † the pseudo-inverse: min L Tr (LΣ) -log |L| subj. to. L is a valid Laplacian, with |L| the pseudo-determinant of the matrix. This optimization function is actually a specific instance of a maximum a posteriori estimation of the Laplacian given a GMRF model N 0, L † [10, Proposition 5]. This corresponds to the graph PSD γ s (0) = 0 and γ s (l) = 1/λ l .

Figure 1 .

 1 Figure 1. Comparison of the graph Fourier modes #2, #3, and #4 (from top to bottom) for various graphs. Edges are omitted for the complete graphs for clarity. Values range from -0.3 (blue) to 0.3 (red) including 0 (white). The left columns shows a corresponding 2D Fourier mode of frequency ν, rotation θ in the 2D plane and phase shift φ ∈ [0, 1]. The parameters (ν, θ, φ) are given below the figures, with the bottom left one corresponding to a modulation of two Fourier modes. Laplacian variations are computed with respect to Gauss CMP graph Laplacian.

Figure 2 .

 2 Figure 2. Theoretical Pearson correlation matrix of the sampled continuous GRF signal model g. Values range from 0 (black) to 0.5 (white). Contrast is increased for better visibility and some entries are actually greater than 0.5. For reference, the spatial Pearson correlation matrix is given (left), and the mean Pearson correlation (off-diagonal elements) is given below the matrices.

Figure 3 .

 3 Figure 3. Graph PSD for the sampled continuous GRF signal g using different graphs. The red dashed lines corresponds to a function γ(λ) = aλ -b . (a, b) parameters are from left to right: (e 1000 , 2.0), (e 0.4 , 0.5), (e 0.5 , 0.5), (e 3 , 1.0), (e, 1.0), (e, 1.0). The yellow dash-dotted line corresponds to the case b = 1.0. The plain black line for Learnt CMP corresonds to b = 6.0. Note that the lowest frequency 0 is not shown here because of the logarithmic scale on the horizontal axis.

Figure 3

 3 Figure 3 plots the graph PSD against the graph frequency in a loglog graph. Several observations can be made out of this figure. First, all graphs show a colored noise model for the PSD: γ(λ) = aλ -b . However, the sparsified Gaussian graphs Gauss KNN and Gauss THR are closer to b = 0.5 rather than the target b = 1 of a 1/λ model. The graph Gauss CMP is also closer to b = 2, however, we observe that most of the energy is contained within the low frequencies, and the target b = 1 is not far off factoring in the noise.

Figure 4 .

 4 Figure 4. (a) Influence on the mean spectral Pearson correlation of the parameter σ in the Gaussian kernel (the curves for Learnt KNN and Learnt THR almost coincide). (b) Weight against the Euclidean distance for Gauss THR with σ set by the heuristic (1/3 of the mean distance) or selected to minimize spectral correlation on (a), and for the graph Learnt THR. The vertical axis is truncated and the points corresponding to very close vertices have very high weight.

Figure 5 .Figure 6 .

 56 Figure 5. Frequency responses of the SGWT kernels gm of various graphs. Bold crosses represent the eigenvalues of the graph.

Figure 7 .

 7 Figure 7. Comparison of the graph Fourier modes #2, #3, and #4 (from top to bottom) for various graphs with the Molène dataset. Edges are omitted for the complete graphs for clarity. Values range from -0.7 (blue) to 0.7 (red) including 0 (white). One value is above 0.7 on the last row for the Learnt KNN graph. Laplacian variations are computed with respect to Gauss CMP graph Laplacian.

Figure 8 .

 8 Figure 8. Pearson correlation matrices of the temperature increments. Values range from 0 (black) to 1.0 (white). For reference, the spatial Pearson correlation matrix is given (left), and the mean Pearson correlation (off-diagonal elements) is given below the matrices.

Figure 9 .

 9 Figure 9. Graph PSD for temperature increments using different graphs. The red dashed lines corresponds to a function γ(λ) = aλ -b . (a, b) parameters are from left to right: (e 0.8 , 0.6), (e 0.2 , 0.4), (e 0.2 , 0.3), (e 1.5 , 1), (e 0.7 , 0.6), (e 0.7 , 0.6). The yellow dash-dotted line corresponds to the case b = 1.0. Note that the lowest frequency 0 is not shown here because of the logarithmic scale on the horizontal axis.

Figure 10 .

 10 Figure 10. Frequency responses of the SGWT kernels gm of various graphs for the Molène dataset. Bold crosses represent the eigenvalues of the graph.

  Figure 11. (a) Local power spectrum of the signal g in dB (horizontally: vertices, vertically: scales, coarser on top), and (b) its mean over vertices in function of the scale.

  Definition 3.2 (Graph WSS 7 ). A stochastic graph signal s is graph WSS if and only if it is WSI through the graph translation T G :

Table 1 .

 1 Comparison of graph stationarity approaches.

† For dimensions greater than 2, stationarity is actually defined as the invariance through a set of operators, e.g. {Tx, Ty} in the 2D case.

ACKNOWLEDGMENTS

This work was supported in part by NSF under grants CCF-1410009, CCF-1527874, CCF-1029373.

APPENDIX A. GRAPH WEIGHTS

Figure 12 shows using colors the weight of the edges for all the graphs studied in this paper. As discussed in the core of the communication, weights on the sparsified Gaussian graphs are very much similar to one another, making those graphs almost unweighted. The consequence is that the graph Fourier transform depends only on the connectivity, and not on the closeness of those vertices. On the other hand, the learnt graph are able to counterweight the influence of the connectivity to truely fit to the data model at hand.