Frédéric Alexandre

Jéremy Fix

Nicolas Rougier

Thierry Viéville

Mai

Algorithmic adjustment of neural field parameters

Keywords: Discrete neural fields, meta-parameters Champs neuronaux discrets, Méta-paramètres

Revisiting neural-field calculation maps in the discrete case, we propose algorithmic mechanisms allowing to choose a right set of parameters in order to both (i) guaranty the stability of the calculation and (ii) tune the shape of the output map.

These results do not "prove" the existence of stable bump solutions, this being already known and extensively verified numerically, but allow to calculate algorithmically the related parameters.

The results apply to scalar and vectorial neural-fields thus allowing to bypass the inherent limitations brought by mean frequency models and also take the laminar structure of the cortex or high-level representation of cortical computations into account.

We obtain an easy to implement procedure that guaranty the convergence of the map onto a fixed point, even with large sampling steps. Furthermore, we report how rectification is the minimal required non-linearity to obtain usual neural-field behaviors. We also propose a way to control and tune these behaviors (filtering, selection, tracking, remanence) and optimize the convergence rate.

This applies to both non parametric profiles, i.e. adjusting the weight values directly, or to parametric profiles and thus adjusting their parameters (e.g. Mexican-hat profiles).

Beyond these algorithmic results, the idea of studying neural computations as discrete dynamical systems and not only the discretization of a continuous system is emphasized here.

The outcome is shared as an open-source plug-in module, called EnaS (http://enas.gforge.inria.fr), to be used in existing simulation software.

Introduction

Using DNF and its extensions. The Dynamic Neural Field (DNF) theory is mainly concerned with the functional modeling of neural structures where information is considered to be encoded at the level of the population rather than at the level of single neurons. Such models first appeared in the 50s, but the theory really took off in the 70s with the works of Wilson and Cowan [START_REF] Wilson | A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue[END_REF] and Amari [2]. At the level of a single neuron, the model that is used is a mean frequency model, while the electrical activity of a neuron is approximated by a single potential. However, there also exists several spiking neuron models that represent both a finer and more accurate model of a real biological neuron. In the framework of the DNF, such models allow to bypass the inherent limitations brought by mean frequency models [START_REF] Rulkov | Modeling of spiking-bursting neural behavior using two-dimensional map[END_REF], this being taken into account via vectorial neural state [START_REF] Viéville | One step towards an abstract view of computation in spiking neural-networks[END_REF]. Furthermore, neural-fields also allow to consider high-level representation of cortical computations [START_REF] Viéville | How do high-level specifications of the brain relate to variational approaches?[END_REF].

These models most generally use excitatory recurrent collateral connections between the neurons as a function of the distance between them and global inhibition is used to ensure the uniqueness of the bubble of activity within the field [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. They exhibit so-called bump patterns, which have been observed in the prefrontal cortex and are involved in working memory tasks [START_REF] Goldman-Rakic | Topography of cognition: parallel distributed networks in primate association cortex[END_REF] or high-level cognition [START_REF] Taylor | Neural networks for consciousness[END_REF].

The dynamics of pattern formation in lateral-inhibition type neural fields with global or local inhibition has been extensively studied in a number of works where it has been demonstrated that these kinds of fields are able to maintain a localized packet of neuronal activity that can, for example, represent the current state of an agent in a continuous space or reflect some sensory input feeding the field [START_REF] Rougier | Emergence of attention within a neural population[END_REF]. A step further, the linear response of neural fields to localized inputs subject to finite propagation speeds has been studied in [START_REF] Hutt | Spontaneous and evoked activity in extended neural populations with gamma-distributed spatial interactions and transmission delay[END_REF]. The characteristic behavior of such fields is the formation of very localized packets of neural activity that tend to represent some consistent information which is present at the level of the input [START_REF] Grimbert | Mesoscopic models of cortical structures[END_REF].

However, the exact shape of these output bumps is quite difficult to predict since it is a generally a non linear consequence of both the lateral connectivity pattern and the input pattern. Thus the question remains on how to control the shape, first in a purely mechanical way (i.e. finding the proper mathematics behind), second in a more biological plausible way. This second step, which is beyond the scope of this paper, requires to complete the first step, which is also going to provide some description of the space of interesting parameters.

The main computational properties (see e.g., [START_REF] Rougier | Emergence of attention within a neural population[END_REF] for a review) of such neural field are: . -1-filtering of the output bump shape: increase versus decrease of size, volume, etc.. . -2-selection of the output bump among several input bumps, . -3-tracking of a moving input bump at the output level, . -4-remanence of the output bump after the partial or total suppression of the input bump. We are going to consider filtering and selection behaviors in this contribution and discuss the tracking and remanence mechanisms also.

The parameter adjustment problem. Here, we want to study which set of lateral weights allow to achieve a specific function, and obtain one unique controlled bump on the output. More precisely, the goal of this paper is to provide algorithmic results to Cortex this parameter adjustment problem. This contrasts with the almost universal practice to adjust empirically the parameters of neural field (even if authors are unobtrusive about it), although constructive results exist and should be used.

In [START_REF] Grimbert | Mesoscopic models of cortical structures[END_REF], for instance, homogeneous stationary solutions (i.e independent of the spatial variable) and bump stationary solutions (i.e. localized areas of high activity) in two-dimensional neural field models composed of excitatory and inhibitory neurons have been studied, in the continuous case, after [START_REF] Coombes | Waves, bumps, and patterns in neural fields theories[END_REF][START_REF] Folias | Breathing pulses in an excitatory neural network[END_REF], providing an explicit solution in terms of Bessel functions and Hankel transform. Basic properties of bounded neural fields (well-posedness, stability of their solutions in the homogeneous or locally homogeneous case) have also been studied, accounting for partial or global synchrony among the cortical columns composing the field [START_REF] Grimbert | Bifurcation analysis of Jansen's neural mass model[END_REF]. These results are based on the compactness of the related operators yielding stability results related to contracting mappings.

These results however are not easy to implement in practical codes, because, in order to obtain analytical results they consider formulations (using e.g. Heaviside nonlinearity and/or unlimited precision assumptions) which do not correspond to numerical implementations. Furthermore, in the discrete case, we may obtain more interesting results than simply approximating using "small-ǫ" results obtained in the continuous case. This means not simply consider the Euler approximation of the continuous equation but also consider more efficient discrete schemes. This also means not derive parameter adjustment algorithmic mechanisms based on approximations of continuous solutions, but build on the discrete solution directly.

What is the paper about. In the next section we review the formalism usually considered for neural field implementation. We discuss the choice of the output-state non linear relation and the way bumps are defined. We also make explicit one link between the continuous and discretized approaches.

In the first result section, we derive a sufficient condition to guaranty the convergence of the neural map towards a fixed point, given the proposed non-linearity, and discuss its applications.

In the next result section, we introduce a linear framework to define the optimal neural field weights given an input-output relation. We instantiate this mechanism for parametrized weights and local bounded weights.

In the conclusion we briefly discuss the limits and applications of the methods proposed here.

Methods

Defining a continuous map.

Let us consider1 a dynamic vectorial map v : S n × R → R m (e.g., for a 2D scalar map n = 2, m = 1) defined from:

   τ v(x, t) = -L v(x, t) + y W(x -y) a(y, t) dy + G i(x), a(x, t) = ρ (v(y, t)) , v(x, 0) = G i(x), (1)
where at location x and time t, v(x, t) represents the unit state (related to the membrane potential in scalar models) and a(x, t) accounts for the unit output positive activity (usually related to the firing rates in scalar models).

Here, i(x) models the input (related to a membrane current in scalar models). It is supposed constant in the sequel, unless explicitly specified (i.e. when considering tracking and remanence).

The spatial domain of the model is a hyper-sphere S n though the domain is assumed to be large enough to neglect boundary conditions. Numerical verifications are going to validate this approximation.

Here, we do not consider only a scalar map but also a vectorial map in order to take more complex neural models (e.g., [START_REF] Rulkov | Modeling of spiking-bursting neural behavior using two-dimensional map[END_REF]) or several cortical layer (see [START_REF] Grimbert | A brief overview of intracortical circuits[END_REF] for an illustration) or more complicated functional mechanisms (see [START_REF] Viéville | How do high-level specifications of the brain relate to variational approaches?[END_REF] for a discussion) into account.

For the input gain G, we assume without loss of generality G = I in the sequel (i.e. considering the G i → i transform).

Similarly, for the diagonal matrix L ≥ 0 (i.e. a positive matrix, thus with positive diagonal elements) which accounts for the neural unit leak, without loss of generality (i.e. considering the v → √ L -1 v transform), we assume L = I in the sequel.

The matrix L is diagonal, but in the vectorial case, the interactions between layers, i.e. between the components of v, are defined via the non-diagonal elements of the W(d) matrices.

The map is parametrized by the recurrent connection weights W. One example are radial "Mexican hat" connections (i.e. with exponential excitatory/inhibitory connections) which writes in the scalar case:

W φ (y) = A + e -|y| 2 /σ + -A -e -|y| 2 /σ 2 -,
parametrized by φ = (A + , A -, σ + , σ -), accounting for the excitatory and inhibitory gains and geometries, respectively. Other parametric instances are also considered in the present framework.

In any case, the goal is to tune the weights W in order to control the map's output. This shall be done either by nonparametric profiles, i.e. thus adjusting the weight values W directly, or by parametric profiles and thus adjusting their parameters φ.

Defining the output-state relation.

In this context, the output-state relation is modeled as a point-like non-linearity, as illustrated in Fig. 1.

Piece-linear Rectification Sigmoid Heaviside

Figure 1: Schematic representation of usual output-state non-linear relations in neural field equations: The sigmoid profile corresponds to the average relation between membrane potential and firing rate. The binary Heaviside profile is considered in order to ease analytical derivations. The piece-wise linear profile allows to take thresholding and saturation into account in conjunction with a linear response. The rectification profile takes into account the fact that the neural activity is constrained to positive values. In this contribution, we consider rectification and, to some extends, piece-linear profiles.

Cortex

The sigmoid profile corresponds to the average relation between membrane potential and firing rate (assuming slow variation of the potential with respect to firing-rate sampling measure and Gaussian distributions of membrane potential and firing threshold, see e.g., [START_REF] Dayan | Theoretical Neuroscience : Computational and Mathematical Modeling of Neural Systems[END_REF] for a review) and it is very often simplified as a piece-wise linear profile (thus with thresholding and saturation bounding a linear response) or even a drastic Heaviside profile used for analytical derivation purposes (e.g., [START_REF] Grimbert | A brief overview of intracortical circuits[END_REF]). The reason for this over-simplification is the need to derive analytical results with the caveat to model the input activity down to a simple on-off process, though analytical studies considering the sigmoidal shape of the transfer function are also available [START_REF] Ermentrout | Existence and uniqueness of travelling waves for a neural network[END_REF][START_REF] Bressloff | Front bifurcations in an excitatory neural network[END_REF]. Another choice is to accept the fact that, up to our current knowledge, the level of modelization of a neural field is rather far from the biological complexity. Applying the Occam razor principle, the fact that the neural activity, viewed as a firing-rate, is thus constrained to positive values can be considered as the minimal assumption related to biological constraints (see, e.g. [START_REF] Grossberg | Nonlinear neural networks: Principles, mechanisms, and architectures[END_REF] for an early contribution and a deep discussion on the subject). This minimal choice also results from the following numerical fact: It has been observed that (see e.g. Fig. 2), with "suitable" parameters, interesting properties of neural fields (filtering, selection, tracking) [20,[START_REF] Rougier | Emergence of attention within a neural population[END_REF] can be obtained using rectification (instead of usual sigmoid-like profiles). Figure 2: An example of input (left-view) output (right-view) transformation via a non-linear discrete map with rectification. The input is noisy (more than 20% of the signal) with three bumps. The output "filters" the noise, and select the highest bump even if the lower bumps have more energy (being wider). This result is obtained, from (2) with Mexican-Hat profiles φ = (0.0015, 0.0015, 0.45 N, N ), N = 100 and δ = 0.99. Furthermore, with the same parameters, if the bump is moving, the mechanism tracks the highest bump.

We thus are going to consider rectification as output-state relation in the sequel. We are also going to explain how the following results still apply to piece-wise linear sigmoid-like profiles.

Using a numeric non-linear discrete map.

Since we want to derive algorithmic results, we must make explicit the numeric implementation of the previous theoretic equations.

INRIA

Let us thus consider the following non-linear synchroneous discrete scheme in time and space, thus, now2 for x ∈ {0, N { n and t ∈ {0, T {, T < +∞:

u + (x, t) = ρ (u(x, t)) , u(x, t + 1) = u + (x, t) + δ -u + (x, t) + y W(x -y) u + (y, t) + i(x) , u + (x, 0) = i(x), (2)
providing 0 < δ < 1, while3 the non-linearity ρ(.) is a rectification, thus with u + (x, t) ≥ 0.

This writes in compact matrix form:

ū+ (t + 1) = F (ū + (t)) = ρ ū+ (t) + δ -ū + (t) + W • ū+ (t) + ī = ρ K ū+ (t) + δ ī , K = 1 -δ 1 -W ,
where ū+ is a R n×m vector, while W and K are R (n×m)×(n×m) Toeplitz matrices.

This corresponds to what is implemented in simulation codes (e.g., [20,[START_REF] Rougier | Emergence of attention within a neural population[END_REF]). However, the link between (1) (with K = L = I) and ( 2) is not straightforward and has to be carefully derived, as performed now.

The derivation is not straightforward because (2) is not the Euler discretization scheme of (1). This latter scheme writes:

u(x, t + 1) = u(x, t) + ∆T u(x, t) = u(x, t) + δ " -u(x, t) + P y W(x -y) u + (y, t) + i(x) " (3)
for a sampling period ∆T , with δ = ∆T /τ and providing u def = v and u + def = a. If v(x, t) > 0, v being differentiable (thus continuous), then in a neighborhood s ∈]t -ǫ, t + ǫ[ v(x, s) > 0, thus u(x, s) = u + (x, s) and u(x, s) = v(x, s) so that both schemes correspond, because (1) and (2) behave as linear equations.

Otherwise, if -v(x, t) ≥ 0 then for identical u(x, t) and u + (x, t) = 0 schemes (2) and (3) verify:

δ (τ v(x, t) + v(x, t)) = u(x, t + 1)| (3) + (1 -δ) (-v(x, t)) = u(x, t + 1)| (2) ≥ u(x, t + 1)| (3) so that: If v(x, t) ≤ 0 and τ v(x, t) ≤ -v(x, t) both schemes still correspond, since 0 ≥ u(x, t + 1)| (2) ≥ u(x, t + 1)| (3) ⇒ u(x, t + 1) + ˛(2) = u(x, t + 1) + ˛(3) = 0.
This conditions occurs if v(x, t) is small enough, in particular in a neighborhood of fixed points. As a consequence, if both schemes converge towards a fixed point, the fixed point is the same (same kind and same value).

However if v(x, t) ≤ 0 and τ v(x, t) > v(x, t) both schemes differ, in the same way the Gauss-Seidel scheme differs from (improving it) the Jacobi scheme for iterative solutions of linear systems. Here the scheme (2) make use of the u + rectified values as soon as available, instead of still using the non rectified values as for the scheme (3). See [START_REF] Porshing | Jacobi and gauss-seidel methods for nonlinear network problems[END_REF] for a technical development and elements showing the interest 4 of the "Gauss-Seidel-like" scheme (2) with respect to the "Jacobi-like" scheme (3).

As a conclusion on this aspect, the series (2) converges to a discrete approximation of (1) corresponding to the Euler discretization scheme, the standard numerical scheme being improved here.

Anyway, it is easy to verify that the derivations proposed in the sequel are obvious to adapt to the scheme (3), if used instead of the scheme (2).

Cortex

Bump as required output.

In this context, the desired output profile is a bump, i.e. a bell-like profile, more precisely a positive, decreasing, radial symmetric profile, as observed for bubble of activities in the cortex [2,[START_REF] Goldman-Rakic | Topography of cognition: parallel distributed networks in primate association cortex[END_REF]. A Gaussian profile is one example of a bump. In order to capture the notion of bump at a more general level, we have proposed to consider Gaussian enumerable linear combinations [START_REF] Frédéric Alexandre | On practical neural field parameters adjustment[END_REF], which is suitable in the linear case, but does not appear to be easily translatable to the non-linear case.

Following another track, we write a radial symmetric profile (without loss of generality around x 0 = 0) as: u(x) = u(r 2 ), r = |x| with u(r) > 0 (positivity) and u(r) ≤ 0 (decrease).

We further notice that the convolution:

(W * u)(|x| 2 ) = y W(|x -y| 2 ) u(|y| 2
) dy of two radial symmetric functions is radial symmetric (easily shown, for instance, considering the invariance by a rotation R of (W * u)(|R x| 2 ) = (W * u)(|x| 2 )) while in the general case, W(x) must be radial symmetric in order the convolution to preserve radial symmetry for any radial function u. We thus are going to consider radial symmetric weights in the present study, thus only function on the magnitude r.

There is however a tricky point here: convolution of radial symmetric functions are not magnitude convolutions. More precisely, in the 2D scalar case 5 :

(W * u)(r = |x|) = s W (s, r) u(s = |y|) ds with W (s, r) = s θ W |r -s| 2 + 2 r s (1 -cos(θ)
) dθ, whereas, for a convolution, an expression of the form W (s -r) would have been expected. This misleading fact is overcome in the second next section where a constructive method is proposed to control the bump radial symmetric profiles.

Results about the convergence

Let us first discuss how to control the discrete implementation convergence.

Implementation as a contracting mapping.

In order to obtain the convergence of the iterative equation ( 2) towards a fixed point we simply require the related mapping to be contracting, i.e.

∂F ∂u = Σ ′ K < 1 where Σ ′ = diag(• • • , ρ ′ i , • • • ) is the diagonal matrix with ρ ′ i = ρ ′ K ū(t) + δ ī ∈ {0, 1}
whether the value is negative or positive, respectively. This allows to write pragmatic conditions of convergence, as discussed now. This also restrains our study to convergence to a fixed point, without taking other asymptotic dynamics into account [START_REF] Coombes | Waves, bumps, and patterns in neural fields theories[END_REF]. 5 The same results indeed occurs in the n-dimensional vectorial case:

v(r) = R s W(s, r) u(s) with: W(s, r) = s n-1 R n dµn W " q |r -s| 2 + 2 r s (1 -e T 1 n) « , e 1 = (1, 0 • • • 0) T , writing x = r m, |m| = 1 and y = s n, |n| = 1 with dy = s n-1 ds dµn, and dµn = dn {•••j••• },j =i n i
, while ds = n T dy, thus not of the form W(sr), because of the s r term in the previous expression.

INRIA

Convergence in the linear case.

Assume that all output values are positive, i.e. ρ i = 1, thus no rectification. This occurs typically with excitatory weights only. This corresponds to a pure linear transformation.

In this case, since ∂F ∂u = K , the convergence condition now writes: | K| < 1 ⇔ | W| < 1 as soon as 0 < δ < 1, i.e. in words: the system is contracting if the magnitude of the global weight matrix is lower than one. If higher than one, the system diverges, in the general case.

Since the matrix is symmetric, because we consider radial symmetric weight profiles, its magnitude (defined as | W| = max ū,|ū=1| | W ū|) corresponds to its spectral radius i.e., its eigen-value of highest magnitude.

The previous result is then a consequence of the fact, easy to verify, that K and W have the same eigen-vectors, with corresponding eigen-values λ K and λ W verifying:

λ K = 1 -δ (1 -λ W ).
This condition is simple to check6 at the numerical level, using the power method [START_REF] Gantmatcher | Matrix Theory[END_REF], which efficiently calculates the highest matrix eigen-value, i.e., the magnitude of the matrix in our case, via the converging series (ū s , λ s ):

ūs+1 = W ūs /λ s λ s = |ū s | ū0 = (1, • • • , 1) T λ 0 = |ū 0 | yielding | W| = λ ∞ .
We numerically observed that convergence of the previous power-iteration method to estimate | W| up to a precision of 10 -3 is obtained in about 2 to 5 iterations for the different results reported here, very often 2 to 3, about 10 in the worst cases. This is due to the choice of ū0 (never orthogonal to the final eigen-vector, as required by the method) which is close to the final value ū∞ . This value has been chosen considering that W is not only a Toeplitz matrix 7 , but close to a circulant matrix (considering the previous boundary conditions), thus with eigen-elements in closed form: This value for ū0 corresponds to a eigen-value close to the weight's average value. It is easy to verify, as a rule of thumb, that the magnitude is always higher than the weight's average value, while we have numerically observed that both values are closed for high magnitude | W|.

Tuning the parameters in the linear case.

The key fact at this point, is that convergence does not depend on δ (as soon as lower than 1). It means that we can use high δ values and speed up the convergence. It also means that it is not always an optimal idea to make all the analysis in the continuous case and then blindly "discretize with a small ǫ".

More than that, the calculation not only tells if the system is contracting or not, but since the equation is linear, allows to calculate the maximal weight magnitude which guarantees the convergence. In other words, the previous algorithm tells how to rescale the gain parameters (e.g., |A -| and |A + | for Mexican-hat profiles) in order to guaranty the convergence.

Although this is not going to be true in the non-linear case, in the linear case, the final result does not depend on δ and the closed-form solution writes, reintroducing the Cortex input gain matrix G defined in [START_REF] Frédéric Alexandre | On practical neural field parameters adjustment[END_REF] for the purpose of the discussion:

ū(∞) = [1 -W] -1 Ḡī , (4)
since ū(t + 1) = K ū(t) + δ Ḡī in this case, thus:

ū(∞) = t Kt δ Ḡī = 1 -K -1 δ Ḡī = 1 -W -1 Ḡī .
This formula is instructive because it explains why using a linear map without nonlinear output-state relation does not allow the recurrent connections parameterized by W to produce a selection, unless Ḡ is singular. The simple reason is that selection means that a whole set of non similar inputs (with the desired bump added to different kind of distractors) map onto the same output (without the distractors). This linear filtering has a kernel reduced to {0} thus is not able to project several inputs onto a unique output. If Ḡ is singular, selection can be produced, but without any need of the recurrent connections. In fact, in the linear case, the DNF is equivalent to a simple feed-forward filter.

Furthermore, e.g. for scalar Mexican-hat profiles, we can describe how convergence varies with |A -|/|A + | and with σ -/σ + , as detailed in Fig. 3.

Using bounded excitation in the non-linear case.

The previous result is quite informative, but does not correspond to our need, since rectification has still to be taken into account. However, we can make use of the previous result and derive the following sufficient boundary condition, in the non-linear case:

| W+ | < 1, W+ = max(0, W) (5)
i.e. in words: the non-linear solution remains bounded and is convergent if the magnitude of the positive (excitatory) weights is lower than one.

INRIA

A few algebra yields, as soon as ū > 0, i.e., for t > 0 :

K ū + δ ī = [1 -δ [1 -W+ ]] | {z } K+ ū + δ ī + δ min `W, 0 ´ū | {z } ū<
with ū< ≤ 0 and ū> = K + ū + δ ī ≥ 0, thus:

˛F(ū) = ρ `K ū + δ ī´˛≤ ˛ρ `K + ū + δ ī´= K+ ū + δ ī˛,
since ˛ρ(ū > + ū< ) ˛≤ ˛ρ(ū > ) ˛. The upper-bounding series defined from:

0 ≤ ū(t + n) ≤ ū> n = " K+ ū> n-1 + δ ī" , ū> 0 = ū(t) is linear and convergent if and only if | K+ | < 1 ⇔ | W+ | < 1, reusing the previous result.
Since ū is positive and bounded by a convergent series, it remains bounded (i.e., convergent or oscillatory inside the bounds). However, the previous derivation yields a stronger result because: This captures the intuitive fact that the negative (inhibitory) weights action is always bounded by the rectification, thus only positive weights matter. An example of practical adjustment curve is given in Fig. 4. Here, not a Mexican hat but a step-wise constant profile has been used, in order to illustrate the fact that the derivation can be performed with any kind of weight profile.

∀i, 0 ≤ [Σ ′ Ku]i ≤ [ K+ u]i since if [ Ku]i ≤ 0 ⇒ [Σ ′ Ku]i = 0 while [ Ku]i ≤ [ K+ u]i,
Without this condition, as corroborated by numerical experiments, the system is in usual conditions divergent. This is easy to figure out: as soon as an input is closed to the W+ eigen-vector related to the eigen-value higher than one, it is going to induce a divergence, unless the inhibitory weights balance this effect. Since positive excitatory weights are concentrated in a short-range neighborhood, a localized bump essentially activates excitatory connections without inhibition, leading to divergence if | W+ | >> 1. This corresponds to what has been observed numerically, as detailed in Fig. 5. This corresponds to the fact that the proposed condition is sufficient but not necessary. Furthermore, this also does not mean that, given a restrained set of inputs, higher parameter's values are not going to yield convergence.

The system is also convergent if | W| < 1 but this is a stronger condition which bounds both excitatory / inhibitory connections whereas the proposed condition does not constraint inhibition since the rectification non-linearity allows to stabilize inhibition.

This condition still does not depend on δ, as soon as δ < 1. It does not mean that the way it converges does not depend on δ: Depending on δ the convergence might occur with dumped oscillations or not, be either faster or slower, etc.. This has been numerically explored in Fig. 6 leading to a simple rule of thumb to adjust δ. Figure 6: Convergence rate measured as the number of iterations to obtain a fixed point up to a mean value variation below 10 -3 , for the same experimental conditions as reported in Fig. 2 (corresponding to the 1st line of results), except that the excitatory weights are reinforced in order to increase the magnitude (subsequent lines of result). We clearly observed the interest of of using non-negligible values of δ, and the influence of the weights magnitude. For small values of the magnitude, a value δ ≃ 1/2 seems optimal with a switch towards δ → 1 for highest values of the magnitude. For convergence rate higher than 10 3 the system is mainly oscillating (dumped oscillations) while always bounded as predicted.

Application to tracking.

The fact convergence does not depend on δ, as soon as δ < 1, allows to introduce high convergence rates as illustrated in Fig. 6. This does not only "save computation time" but allows to track a moving input at a high rate since convergence to the modified state is improved. As far as biological neural maps are understood the convergence of the neural field towards a new state is not the result of the convergence of a progressive iterative process, but a tip over from one state to another when the input changes.

Here, our study shows how to adjust δ in order to obtain optimal convergence rate, thus optimally track a moving input.

What we have numerically experimented with the conditions reported in Fig. 2 (thus with | W+ | = 0.0), is that the tracking of a moving bump is optimally tracked for δ ≃ 1/2.

INRIA

Application to remanence.

Short-term memory behavior is related to the capability to maintain the output for a short period of time, even if the input vanishes. This is a behavior inverse from the previous one, since we do not require the system to quickly converge to the new state (corresponding to a null input), but remains in the old state (related to the previous input).

In our framework this may correspond to | W+ | ≥ 1, as shown, in the linear case, by the formula (4) and experimented numerically in Fig. 5, for 1 ≤ | W+ | ≤≃ 1.5. If | W+ | ≥ 1, the linear transfer function has "poles" (corresponding to resonances in oscillatory regimes), i.e. unbounded output for infinitesimal inputs. Numerically, this corresponds to high-gains (corresponding to dumped resonances in oscillatory regimes), and slow time decays, as required for such behavior. Tuning the scheme parameters in this range allows to approximate remanence, although this yields solutions at the edge of instability.

Another strategy is to consider not a rectification but a piece-wise linear sigmoid profile as shown in Fig. 1 and consider | W+ | ≥ 1. In that case the scheme tends to increase the output or yield divergence, but because of the saturation, the output remains bounded. This allows to easily generate remanence, as it is observed in numerical simulations [START_REF] Rougier | Emergence of attention within a neural population[END_REF], and the behavior is more robust. The add-on of the previous derivations allows to know the scheme parameters range to obtain this behavior.

Results about the bumps

Analysis of the fixed point equation.

Let us now discuss how to obtain, at the implementation level, bumps with controlled shapes. Since there is a fixed point converging dynamics, we simply have to consider the fixed point equation in this section. The scheme (2) or (3) fixed point equation writes, given a radial input i(|x|) and output u(|x|, +∞):

z = Ū+ W + ( ī -ū+ -ū-/δ) = 0 (6)
with ū+ = ρ(ū) = max(0, ū) and ū-= ρ(-ū) = min(0, ū), thus ū = ū+ + ū-, while: U + k (|x|) = y,|y|=r k u + (|x -y|) defines the convolution in the discrete case. In the 2D case for instance, for weights defined in a bounded circular neighborhood of radius R, this convolution writes:

ˆŪ + W ˜(x, y) = P K k=1 W [k] P (i,j),i 2 +j 2 =r 2 k ≤R 2 u(x -i, y -j) = W [1] u(x, y) + W [2] (u(x + 1, y) + u(x -1, y) + u(x, y + 1) + u(x, y -1)) + • • •
as represented in Fig. 7 for R = 5, while it is obvious, from a piece of symbolic computation, to derive for any R the K related terms, with: showing a supra-linear and sub-quadratic (numerically K = O(R ≃1.7 )) increase of the number of radial coefficients. This fixed-point equation splits in two parts: Figure 7: Representation of the radial symmetric discrete convolution indexing in the 2D case for a bounded circular neighborhood of radius R = 5. The 2D locations of the indexes are drawn in the top array and the correspondence between index and radius in the bottom table. This indexing has been derived using a piece of symbolic computation maple code, also allowing to automatically generates the optimal numerical C/C++ code for any convolution of radius R.

R
u(|x|) > 0 z|x| = Ū+ W + ( ī -ū+ ) |x| ⇒ [1 -W] u + (|x|) = i(|x|) u(|x|) ≤ 0 z|x| = Ū+ W + ( ī -ū-/δ) |x| ⇒ [1 -1/δ -W] u + (|x|) ≤ i(
• A linear response for points inside the bump (with u(|x|) > 0), which does not depend on δ, and is related to the output response adjustment.

• An adaptive threshold cut-off mechanism outside the bump (with u(|x|) ≤ 0), which allows to perform selection, and is related to input inhibition control.

-The key point is that, given a desired output u + , the equation is linear with respect to the weights W.

-Furthermore, δ simply scales the inhibited part of the output u -and is thus not an independent parameter.

-Here u -/δ is a hidden parameter of the problem: It allows to specify the level of inhibition against input perturbations.

-There is a scale invariance, in the sense that multiplying input and output by the same scale factor, does not change the equation with respect to W. We thus are going to consider, without any loss of generality, normalized input in the sequel, i.e. considering the volume:

V (i) = x |i(x)| (7)
input with V (i) = 1.

-Finally, the output response adjustment is entirely defined by the input-output difference [ ī-ū+ ] |x|,u(|x|)>0 , which is yet another simplification of the problem structure.

Given these properties, let us now formalize the output control.

Position of the problem.

The fact we want to "tune" the weights W, as illustrated in Fig. 8, can be formulated as follows: Given a set of input īl , indexed by l and a set of corresponding desired output ū+ l , the goal is to find the desired optimal weights W that maps the input onto the output.

INRIA

Cortex

gain" for this point and sample). A zero value allows to eliminate the related measure in the estimation process, a high value to make it important. For scalar maps (m = 1) Λ |x|,l ≥ 0 is a number, otherwise it is a m × m positive symmetric matrix. As far as Λ+ l is concerned (i.e., for point with u(|x|) > 0), this corresponds to pondering the output value, which is an interesting feature. As far as Λl is concerned (i.e., for point with u(|x|) < 0), this corresponds to pondering the "absence" of the output value, i.e. the inhibition level. This inhibition level is thus not tuned here though a value but though a kind of "estimation gain". An alternative, at the specification level, would be to simply specify u -/δ, i.e. the inhibition value itself: This is very easy to introduce in [START_REF] Folias | Breathing pulses in an excitatory neural network[END_REF] which is now the linear equation defining W.

This specification choice has also a computational advantage, since using ( 9) is numerically more stable than the previous equation: If u -/δ is not specified as an income of the problem, the previous derivation corresponds to eliminating u -/δ (solving [START_REF] Ermentrout | Existence and uniqueness of travelling waves for a neural network[END_REF] where u(|x|) < 0, thus considering input inhibition) and reintroduced the result in order to find w (solving [START_REF] Ermentrout | Existence and uniqueness of travelling waves for a neural network[END_REF] where u(|x|) > 0, thus considering output adjustment). If u -/δ is directly specified, then (9) allows us to linearly combine both output control and adjustment inhibition in the estimation.

A step further, equation ( 9) is a "Hebbian" equation in the following sense: it is local and only depends on the output value and linearly on the input value. In the present study, for the sake of simplicity, we implicitly assume that the weights W do not depend on |x|. However, this is not a limitation of the method: The criterion ( 8) is still well-defined in the more general case of non-constant weights. Then, equation ( 9) is "local" in the sense that Ū+ l only depends on the output values "around" x (more precisely in a neighborhood corresponding to the size of the output bump) while it only linearly depends on the input at point x. This is not the case if we do not introduce u -/δ as an income of the problem.

As a consequence, we propose to use [START_REF] Folias | Breathing pulses in an excitatory neural network[END_REF] in the sequel. Furthermore, for the sake of simplicity, we are going to consider Λ = I (no pondering), while we simply consider a singleton of input-output couple. The idea is to validate and study the method in a relatively simple case, before making profit of it at its whole level of generality.

Parametric 2D scalar piece-wise constant weights adjustment.

In order to apply the previous framework we start considering the simplest possible radial symmetric weights profile: W (|y|) = if |y| < S then + W 0 else -W 1 in words an excitatory flat constant profile of magnitude W 0 in a central area of radius S and a spatially unbounded flat constant inhibitory profile of magnitude W 1 outside the disk of radius S, as illustrated in Fig. 9. This is thus a parametric weight profile, parametrized with φ = (W 0 , S, W 1 ).

We further assume that S ≪ R, R being the bump radius, i.e., that the excitatory weights radius is small with respect to the bump size, which allows to derive an approximate closed-form for the input-output relation.

More precisely, this allows to write after a few algebra:

u(r = |x|) = |x-y|<S u(|y|) = W [0,R-S] (r) π S 2 u(r) -π S 4 2 + O(S 6 + S 4 r 2 ) + W [R-S,R+S] (r) υ(r 2 ) , with W [a,b] (r) = if a ≤ r ≤ b then 1 else 0 while υ(r 2
) is a term neglected in the sequel, since S ≪ R. This leads to the following closed-form formulas:

u(r) = α -1 (i(r) -β ǐ) +O(S 4 ) r ∈ [0, R -S] linear response i(r) ≤ β ǐ +O(S 4 ) r ∈ [R + S, +∞] threshold inhibition with α def = 1-π (W 0 +W 1 ) S 2 and β def = (α+π W 1 R 2 ) -1 W 1 while ǐ = |y|<R i(|y|).
Let us derive the linear response equation. If

u(r = |x|) > 0: u(x) = P y W(x -y) u(y) + i(x) = W0 P y,|x-y|<S u(y) -W1 P y,|x-y|≥S u(y) + i(x) = (W0 + W1) P y,|x-y|<S u(y) -W1 P y u(y) + i(x) = (W0 + W1) û(x) -W1 ǔ + i(x) writing ǔ = P y u + (y) = P y,|y|<R u(y). Considering the approximation û(r = |x|) ≃ W [0,R-S] (r) π S 2 u(r) + O(S 4 ) allows to further write: u(x) = π (W0 + W1) S 2 u(x) -W1 ǔ + i(x) + O(S 4 ) thus: ǔ = P x,|x|<R u(x) = π (W0 + W1) S 2 ǔ -W1 π R 2 ǔ + ǐ + O(S 4 ) considering P x,|x|<R 1 ≃ R x,|x|<R 1 = π R 2
, leads to the linear response equality, written above in the scalar case.

Otherwise, if u(r = |x|) ≤ 0 and |x| > R + S: P y W(xy) u(y) + i(x) ≤ 0 W1 ǔ + i(x) ≤ 0 leads to the threshold inhibition inequality.

In order this equation to be well defined, we must have: π (W 0 + W 1 ) S 2 < 1, which is in direct relation with the convergence condition (but does not correspond to it, since we use an approximation here).

This derivation is informative and shows that, up to the proposed approximation, the response is mainly a thresholding mechanism: the effect of the computation is to cut the input under the given threshold, thus eliminating spurious distractors of the main input bump and reducing also the output volume. The calculation is non trivial because the threshold β ǐ is adaptive and related to the input volume around the output bump. The radius R of the output bump is an implicit parameter depending on the weight parameters, thus not easy to read on the formula. These aspects are further detailed and verified numerically in Fig. 10 and Fig. 11, providing also some abacus to choose the desired output bump.

Non-parametric 2D scalar spatial-bounded weights adjustment.

In this second numerical application, we consider non-parametric 2D scalar spatialbounded symmetric radial weights, as made explicit in Fig. 7. . When W1 = 0 there is no filtering (output equals input). If W1 is too small the excitatory weights diffuse the input all over the map. If W1 is too high (as for S = 2 and W1 > 0.073) inhibition is high enough to cancel the whole input. Larger S increase diffusion of the input thus enlarge the output radius. Higher inhibition magnitudes W1 increase the thresholding thus decrease the output radius. The output volume is also obviously controlled by the input gain (taken as 1 in the derivation to ease the notations).

Input

Output Radius Output Volume Now the profile is parametrized by the radial values W [i], i ∈ {1, K}. In order to instantiate the general method proposed in ( 8) and ( 9), we fix a form for the input and output, presently:

INRIA u(r) = A H(R -r) 2/π [R 2 -r 2 ]/R 4 + H(r -R) T e -(4/π A/(R 3 T )) (r-R) -1 i(r) = H(C -r) 2/π [C 2 -r 2 ]/C 4
thus quadratic input and output profiles with predefined parameters, as detailed in Fig. 12, while the inhibited profile is exponential and chosen in order u to remains C 1 . The numerical approximation of the convolution writes:

U s (r) ≃ 2 π s θ u + (r 2 -2 r s cos(θ) + s 2 ) / θ 1 as illustrated in Fig. 12.

The general criterion [START_REF] Ermentrout | Existence and uniqueness of travelling waves for a neural network[END_REF] simply writes in this case: min W |z| + L s dW[s] 2 for a regularization parameter L.

Numerical results are reported and discussed in Fig. 13.

Application to inhibition control of the output radius.

As a final numerical application, let us consider parametric 2D scalar Mexican-hat symmetric radial weights, as made explicit in the section 2.

We reuse the previous form of input and output (bounded quadratic response and exponential inhibition) with A = 1 since the output volume control is easily tuned by the input gain. We fix T = 1 since it has been observed that this parameter value is not critical. We further do not consider regularization (L = 0) since the profile is parametric thus de-facto regular.

Since we are interested here in controlling the output radius, we are going to explore the relation between the input radius C and the output radius R.

Regarding the Mexican hat profile, we consider that the weight width parameters σ + and σ -are structural parameters, thus fixed by the connectivity. A step further, since we are interested in selection and remanence we choose |W + | = 0.9 so that the excitatory weights magnitude |A + | is fixed from σ + as detailed previously.

We are thus left with controlling the inhibitory weights magnitude |A -| obtaining results given in Fig. 14. The interesting conclusion is that Mexican-hat profiles (compared to the previous profiles) are far from being optimal to control the output radius. Sharp profiles (i.e. with excitatory/inhibitory weights clearly separated) seem a better choice. Furthermore, not the magnitude but the width of the inhibition is the best control variable to adjust the output radius. The top-left curve corresponds to a standard input-output transfer, and the optimal weights correspond to pure inhibitory profiles, which corresponds to the usual practice of experimenters. However, the obtained shape might be counter-intuitive in the sense that inhibition only increases with eccentricity up to the weights spatial bounds, without any decrease for high eccentricities. The top-right curve corresponds to a standard input-output transfer, but estimated with a much higher regularization parameter L. The nice result is that the shape is a bit smoother as expected, but in fact only marginally depends on this extra parameter, which is non critical as soon as set to a reasonable value. The middle-left curve corresponds to an output radius reduction with respect to the input radius, while the middle-right curve corresponds to an output radius expansion, which appears to be obtained by a non monotonic shape, according to the proposed criterion. The bottom-left curve corresponds to the choice of a lower inhibition level T , again showing the non-critical impact of this parameter, except an expected reduction of the inhibition weights magnitude. The bottom-right curve correspond to bump magnitude amplification, which is coherent with the apparition of excitatory weights in this case. All profiles must verify |W| < 1 since they map bounded input onto bounded output, thus in a case where the scheme is converging.

Discussion and conclusion

The previous results do not "prove" the existence of stable bump solutions, this being already known [START_REF] Wilson | A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue[END_REF]2] and extensively verified numerically. Here the goal is to calculate algorithmically the related parameters. This modest objective is attained since we now are able to tell which parameter set yields convergence and how to control the output and the shape of the bump.

INRIA R |A -| R |A -| R σ- σ+

Further applications of bump control.

One track is to consider how tuning distributed propagation delays and also feedback delays [START_REF] Atay | Neural fields with distributed transmission speeds and long-range feedback delays[END_REF] can, at the implementation level, control bumps but also traveling fronts [START_REF] Hutt | Effects of nonlocal feedback on traveling fronts in neural fields subject to transmission delay[END_REF] .

Application to other paradigms.

In a recent work [20] a model has been designed that performs global competition, only using local connections, with diffusion of the inhibition throughout the network. This is far quicker to have a few local interactions when computing activity within the network and makes the model a real candidate for distributed computations. This paradigm is taken into account here, with non-parameteric weight adjustments. A slight modification of the previous derivation, allows to not only consider synchronous but asynchronous sampling (see [START_REF] Rougier | Emergence of attention within a neural population[END_REF] for details) since the map is still contracting when using asynchronous sampling. This has been numerically experimented, as shown in Fig. 15 and has to be further investigated. Interesting enough is the fact that this simulation has been implemented on an event-based neural network simulator, in order to consider neural maps not as simple "analog" calculations, but both analog and event-based [START_REF] Rochel | An event-driven framework for the simulation of networks of spiking neurons[END_REF].

We have proposed the numerical results only for 2D scalar maps, for the simple reason we are publishing these results . . on a 2D scalar support. All derivations have been developed for the vectorial case and are thus available for higher-dimensional paradigms.

Here we have restrained the discrete scheme to (2) with 0 < δ < 1. It is an interesting perspective to consider convergence and output control for 1 ≤ δ < 2 or higher values.

Beyond parameter adjustment.

Beyond such technical variants, the idea is to get rid of continuous schemes of the form (1), whereas directly study to which extends a discrete dynamical system of the form ( 2) is of interest to analyze and simulate the neural maps behaviors. This is, for instance, the case in [START_REF] Rulkov | Modeling of spiking-bursting neural behavior using two-dimensional map[END_REF] where the dynamics of spiking and spiking-bursting activity of real biological has been studied via a discrete two-dimensional map defined with one fast and one slow variable. As a consequence, discretized neural fields power of description is still to be explored for both distributed computational paradigms and biologically plausible computational map models A Tuning bumps in the linear case Bump as required output In this context, the desired output profile is a bump, i.e. a positive, decreasing, radial symmetric profile,as observed for bubble of activities in the cortex [2,[START_REF] Goldman-Rakic | Topography of cognition: parallel distributed networks in primate association cortex[END_REF]. A Gaussian profile is one example of bump. In order to capture the notion of bump at a more general level, we propose to consider Gaussian enumerable linear combinations, i.e. Gaussian series: b(p) = s∈{s1,••• },s>0 b * (s) e -s |p| 2 , writing s = 1 2 σ 2 for an enumerable set of "width" s, which is of common use (e.g. [START_REF] Wennekers | Separation of spatio-temporal receptive fields into sums of gaussian components[END_REF]). Small s for flat Gaussian (s = 0 is the "constant part"), Large s for tight Gaussian.

The goal is thus to find the parameters allowing to generate bumps, controlling their width, amplitude and shape.

Let us now discuss how to obtain, at the implementation level, bumps with controlled shapes, in the linear case. This observation means that bumps can be approximated using Gaussian series for the simple reason that Gaussian series coefficients are nothing but the discrete approximation of the bump β profile's Laplace transform. This representation thus benefits from all the Laplace transform properties: Here, not a Mexican hat but a step-wise constant profile has been used, since the derivation can be performed with any W profile.

Gaussian series yield "nice" bumps. The previous result is useful because it is very easy to shape the bump using suitable coefficients, for instance using rational profiles of the form β 1 (r) = 1/(1 + |r| n/2 ), n even since: b 1 (s) * = 1 n n i=1 e -s sin(α i n ) sin s sin β i n + α i n α i n = π n 2 i + n 2 -1 β i n = π n (1 -2 i) or of the form β 2 (r) = 1/r n , n > 0 since: b 2 (s) * = s n-1 /n -1! as illustrated in Fig. 17. Such profiles correspond to desired bumps in neural field applications, among many other suitable alternatives. They come with closed-form formula and show how malleable are such shapes.

Note that a simple change of scale allows to adjust the width and the amplitude of the bump, the challenge here being to adjust the shape.

Gaussian series are linear map fix points. Linear map fixed points u • are defined by -u • + W * u • + i = 0, where * stands for the convolution operator. Since convolution and linear combination of Gaussian are Gaussian, providing that the input i and the weights profile W are Gaussian series, the output u • is a Gaussian series. This result is made explicit by the following formula for the Gaussian series coefficients:

-u • (ν) + π r>ν,s>ν, r s r+s =ν W (r) u • (s) r+s + i(ν) = 0,

Figure 3 :

 3 Figure 3: Representation of | W| magnitude variations for Mexican-hat profiles. Here the 2Dscalar map size is 100×100 and σ+ = 10, while A+ = 1/20, for magnitudes between 0 (black) and 1 (white). Results are obtained in less than 5 iterations of the power-method. We observe the non trivial increase and decrease of the magnitude depending on the inhibitory/excitatory gain and width ratios: It increases with σ-/σ+ and decreases with |A-|/|A+| when σ-< σ+ and |A-| < |A+| which corresponds to the intuition that inhibition "balances" excitation. A dual effect is observed when σ-> σ+ and |A-| > |A+|: Inhibition, when not thresholded by the rectification mechanism, also induces high magnitudes, thus instabilities. In other conditions the magnitude is not monotonic.

  as derived in the previous lines. Thus ˛∂F ∂u ˛= ˛Σ′ K˛≤ ˛K ˛< 1 as soon as | W+ | < 1 and the series of scheme (2) is contracting as desired.

Figure 4 :

 4 Figure 4: An example of maximal gain curve in the non-linear case. Given σ+ in abscissa (in pixel unit for a 100×100 2D scalar map), the normalized maximal value of log 10 (A+) is drawn.Here, not a Mexican hat but a step-wise constant profile has been used, in order to illustrate the fact that the derivation can be performed with any kind of weight profile.

Figure 5 :

 5 Figure 5: Behavior of the non-linear discrete scheme, when | W+ | > 1, for the same experimental conditions as reported in Fig. 2 except that the excitatory weights are reinforced in order to increase | W+ |. Given | W+ | in abscissa, the input/output gain G is drawn. We numerically observe the divergence for | W+ | >≃ 1.5. In the intermediate range | W+ | ∈ [1.. ≃ 1.5] the scheme is still converging to fixed point, but with an exponential gain increase and small residual oscillations observed during the numerical experiment.

Figure 9 :

 9 Figure 9: Left view: Radial symmetric piece-wise constant weight parameters profile, W0 being the central excitatory magnitude, W1 the peripheral inhibitory magnitude, while S is the radius of the circle separating both weight's regions. Right view: Schematic representation of the convolution of a bump u of radius R by a radial symmetric sum of radius S ≪ R. The result û is decomposed in four zones, as made explicit in the right part of the diagram, and obtained when deriving the closed form approximation of this convolution.

1 Figure 10 :

 110 Figure 10: Numerical verification of the capability to control the output radius and volume, for a normalized input profile of radius 20 and a 2D 100 × 100 scalar map. The weight radius values are S = {2 (red), 5 (green), 10 (blue), 20 (magenta)} while the excitatory magnitude W0 is adjusted from S in order to have |W| = 0.99. The inhibitory magnitude varies in W1 ∈ [0, 0.1]. When W1 = 0 there is no filtering (output equals input). If W1 is too small the excitatory weights diffuse the input all over the map. If W1 is too high (as for S = 2 and W1 > 0.073) inhibition is high enough to cancel the whole input. Larger S increase diffusion of the input thus enlarge the output radius. Higher inhibition magnitudes W1 increase the thresholding thus decrease the output radius. The output volume is also obviously controlled by the input gain (taken as 1 in the derivation to ease the notations).

1 Figure 11 :

 111 Figure 11: Numerical verification of the capability to perform selection between the main input and distractors, for a normalized input profile of radius 20 and a 2D 100 × 100 scalar map, plus a distractor of width 10, height B and eccentricity D. Here the radius is fixed S = 5 and the excitatory magnitude W0 to have |W| = 0.99. The inhibitory magnitude varies in W1 ∈ [0, 0.1]. The green curves correspond to D = 20, B = 0.9, the blue curve to D = 15, B = 0.9, and the magenta curve to D = 10, B = 0.9. The selection mechanisms (comparing the results with those of Fig.10) clearly works and we obtain similar results, as predicted by the previous derivation, since the input dependence is mainly related to the input volume in the output radius disk. Small variations qualitatively correspond to this radius variability. The red curves correspond to D = 20, B = 10 thus the distractor is higher than the central bump and the selection fails, as expected. In this latter case, when inhibition increases, the overall amount of input yields an output extinction.

Figure 12 :

 12 Figure12: Schematic representation of the two main elements of the numerical set-up. Left view: An "generic" input of radius C and volume 1 in correspondence with a generic output of height H corresponding to a volume A and radius R is chosen. The level of inhibition is fixed by an additional parameter T corresponding to the asymptotic value of the inhibited activity. Right view: A schematic representation of the numerical integration of the convolution at a given radial location.

RR = 10 ,

 10 = 10, T = 1, A = 1, C = 10, L = 1 R = 10, T = 1, A = 1, C = 10, L = 10 basic input-output transfer input-output transfer with higher regularization W k W k R = 10, T = 1, A = 1, C = 20, L = 1 R = 20, T = 1, A = 1, C = 10, L = 1 T = 0.1, A = 1, C = 10, L = 1 R = 10, T = 1, A = 5, C = 10, L = 1 input threshold reduction bump magnitude increase

Figure 13 :

 13 Figure 13: Numerical results for non-parametric weights adjustment, here for W [k], k ∈ {1, 50}.The top-left curve corresponds to a standard input-output transfer, and the optimal weights correspond to pure inhibitory profiles, which corresponds to the usual practice of experimenters. However, the obtained shape might be counter-intuitive in the sense that inhibition only increases with eccentricity up to the weights spatial bounds, without any decrease for high eccentricities. The top-right curve corresponds to a standard input-output transfer, but estimated with a much higher regularization parameter L. The nice result is that the shape is a bit smoother as expected, but in fact only marginally depends on this extra parameter, which is non critical as soon as set to a reasonable value. The middle-left curve corresponds to an output radius reduction with respect to the input radius, while the middle-right curve corresponds to an output radius expansion, which appears to be obtained by a non monotonic shape, according to the proposed criterion. The bottom-left curve corresponds to the choice of a lower inhibition level T , again showing the non-critical impact of this parameter, except an expected reduction of the inhibition weights magnitude. The bottom-right curve correspond to bump magnitude amplification, which is coherent with the apparition of excitatory weights in this case. All profiles must verify |W| < 1 since they map bounded input onto bounded output, thus in a case where the scheme is converging.

Figure 14 :

 14 Figure 14: Numerical verification of the capability to control the output radius with inhibition adjustment. The input radius values are C = {10 (red), 20 (green), 35 (blue), 50 (magenta)}. The inhibitory magnitude varies in |A -| ∈ [0, 0.02] while other parameters are adjusted as discussed in the text. Left view: σ-/σ+ = 5. The result is somehow negative: the output radius is always smaller than the input radius and has only tight variations (in the range where the inhibition is high enough to perform filtering and small enough to avoid the output extinction). Middle view: σ-/σ+ = 1.5. The result is now less negative: the output radius can be controlled by inhibition adjustment. However the output radius is always higher than the input radius. Right view: Here |A -| = 0.01 and the σ-/σ+ ∈ [1, 5] is varied. It is appears clearly that this controlled variable is much more efficient to adjust the output radius.

Figure 15 :

 15 Figure 15: An example of asynchronous sampling of such maps (event-based implementation), applying convergence criteria derived here. We have numerically verified the conjecture that the present results apply when using asynchronous sampling. Left view: intermediate result, the fact asynchronous sampling yields randomization is visible. Right view: final result, after convergence.

  Gaussian series approximate bumps. Radial symmetric profile means that b(p) = β(|p| 2 ), normalized iff p b i (p) = 1. They are positive iff β i (r) ≥ 0 and decreasing iff β ′ i (r) ≤ 0. This puts our specifications in equation. For Gaussian series: β(r) = s b * (s) e -s r = L(b * )(r) where L stands for the Laplace transform, while b * has been extended to the continuum INRIA using a Dirac comb, while the Post-Bryan inversion formula 8 [5] allows to consider the Laplace transform of such positive functions.

  Bump are normalized iff s>0 b * i (s)(π/s) n/2 = 1. Several additional properties characterize such bumps: they are flat at 0 unless singular (i.e. |b(0)| < +∞ ⇒ |∇b(0)| = 0), vanish at infinity (lim r→+∞ β(r) = 0) thus without constant component (lim s→0 b * (s) = 0), etc..

Figure 16 :

 16 Figure 16: An example of maximal gain curve in the non-linear case. Given σ+ in abscissa (in pixel unit for a 100 × 100 2D map), the normalized maximal value of log 10 (|A+|) is drawn.Here, not a Mexican hat but a step-wise constant profile has been used, since the derivation can be performed with any W profile.

RR n°6923

Notations: Scalar are written in italic, vectors are written with lower-case bold letters and matrices with upper-case bold letters. The identity matrix is written I. For a scalar function f () applied to a vector u, v = f (u) means ∀i, v i = f (u i ).INRIA

We write {a, b{= {n ∈ N , a ≤ n < b}.

Here we write ∀i ∈ {1, m}, u + i ≥ 0 as u + ≥ 0 and ∀i ∈ {1, m},u + i = ρ(u i ) as u + = ρ(u). Furthermore: ρ(u) = max(0, u) = (if u ≥ 0 then u else 0) = lim ǫ→0 ρǫ(u)with ρǫ(u) def = H(u) u e -ǫ/u , H() being the Heaviside function, are equivalent definitions of a rectification. The ρǫ() function is the C ∞ mollified form of the rectification.

[START_REF] Bressloff | Front bifurcations in an excitatory neural network[END_REF] Interesting enough is the fact that some co-authors of the paper use this scheme for years thanks to their "computer-scientist" intuition, before the following discussion has been made explicit thanks to the "physicist" of the group.RR n°6923

Numerical implementations used in this paper are available at http://enas.gforge.inria.fr.

The previous condition, although very efficient in practice, does not directly make profit of the fact that W is a Toeplitz matrix. It might however be interesting to explore this algebraic property and analyze to which extends this yields a stronger result about convergence, especially in the non-linear case.RR n°6923

Post-Bryan result: A continuous function on R + which is of exponential order for some c (i.e. ||b(r)|/e c r | ∞ < ∞) with Laplace transform b * (s) is non-negative iff ∀k ≥ 0, ∀s ≥ c, (-1) k b * (k) (s) ≥ 0.RR n°6923

Acknowledgment: Partially supported by the ANR MAPS project.

the related desired output u + , the goal is to find the optimal weights W. The inhibited output u -is a hidden parameter of the problem.

This can be turned into a variational specification: Given an input-output set {• • • ( īl , ū+ l ) • • • }, the optimal weights minimize:

where:

-The output error z |x|,l magnitude at point |x| for the input-output sample of index l (defined in [START_REF] Coombes | Waves, bumps, and patterns in neural fields theories[END_REF]) is minimized, up to some pondering matrix Λ |x|,l allowing to take into account the relative importance of a sample of index l and a map location x.

-The weights are regularized in the sense that the weight's gradient ∇W is also minimized, up to some pondering matrix L, allowing to avoid ill-posed specifications.

Here, for the sake of simplicity, we consider quadratic norms, yielding linear equations for the extrema, although in this finite discrete case, other norms are rather easy to introduce. This specification allows to immediately derive:

writing Λl the block-diagonal matrix defined by Λ |x|,l , while ∆ L W is the discrete Laplace operator accounting for the regularization term. Here we deliberately use oversimplified notations for this very standard regularization mechanism ∆ L W, see e.g., [START_REF] Viéville | How do high-level specifications of the brain relate to variational approaches?[END_REF] for further technical details.

Specification of the estimation parameters.

The hidden variable ūl /δ is defined at points where u(|x|) ≤ 0 and can be derived from ( 8):

0), i.e., taking only points with u(|x|) < 0, thus where ūis defined, into account. This yields after a few algebra:

writing Λ+ l the block-diagonal matrix defined by Λ + |x|,l = (if u(|x|) > 0 then Λ |x|,l else 0), while M † is the pseudo-inverse of the matrix M. This last equation, entirely defined from [START_REF] Ermentrout | Existence and uniqueness of travelling waves for a neural network[END_REF], specifies the weights W as a function of i and u + .

Let us discuss this last point. The pondering matrix Λ |x|,l ≥ 0 defines, for each point and sample, the relative weight of the estimation (in other words the "estimation Cortex This is true for rather general weight profiles, not necessarily "Mexican-hat" profile only.

In this specific case, it is possible to derive a closed-form expansion, convergent as soon as |A ± | < 1, which is the case in practice. We obtain up to a given order r:

with η r (p, q) = (p+q)! p! q! easily derived by symbolic calculation. Note that the previous two formulas are written in the scalar case for the sake of simplicity, whereas it is straightforward to derive them in the vectorial case. It is also straightforward to derive (using symbolic calculation) for more general kernels.

Practical implementation.

From the previous results it becomes possible to design numerical routines in order to adjust at will the bump's shape. The previous result development means that we do not have to adjust the profile in the 2D space but simply in the 1D radial profile, thanks to the Gaussian series representation. This fact very likely generalizes to non-linear profiles.

For instance, given an input i(ν) and a desired output b • 1 (ν) the following program:

Λ allows to find the optimal weight parametric value θ, using e.g. the explicit formulas given previously. The choice of the measure ν and of the metric |.| Λ is application dependent.

Another application is the input control of such calculation maps: Given parameters θ, find the input i(ν) yielding a given output b(ν). The solution is straightforward in the present framework, as illustrated in Fig. 18.


Kavya Gupta
email: kavya.gupta100@gmail.com

Fateh Kaakai

Béatrice Pesquet-Popescu

Jean-Christophe Pesquet

Safe Design of Stable Neural Networks for Fault Detection in Small UAVs

Keywords: safety, stability, Lipschitz constant, verification, machine learning, neural networks, UAV, tabular data, adversarial attacks

Stability of a machine learning model is the extent to which a model can continue to operate correctly despite small perturbations in its inputs. A formal method to measure stability is the Lipschitz constant of the model which allows to evaluate how small perturbations in the inputs impact the output variations. Variations in the outputs may lead to high errors for regression tasks or unintended changes in the classes for classification tasks. Verification of the stability of ML models is crucial in many industrial domains such as aeronautics, space, automotive etc. It has been recognized that data-driven models are intrinsically extremely sensitive to small perturbation of the inputs. Therefore, the need to design methods for verifying the stability of ML models is of importance for manufacturers developing safety critical products. In this work, we focus on Small Unmanned Aerial Vehicles (UAVs) which are in the frontage of new technology solutions for intelligent systems. However, real-time fault detection/diagnosis in such UAVs remains a challenge from data collection to prediction tasks. This work presents application of neural networks to detect in real-time elevon positioning faults. We show the efficiency of a formal method based on the Lipschitz constant for quantifying the stability of neural network models. We also present how this method can be coupled with spectral normalization constraints at the design phase to control the internal parameters of the model and make it more stable while keeping a high level of performance (accuracy-stability trade-off).

Introduction

Machine Learning (ML) methods and in particular neural networks are rapidly paving their way in various mission-critical and safety-critical domains such as aeronautics, aerospace, automotive, railways, and nuclear plants. In particular, UAVs are currently being used in various applications such as surveillance of critical infrastructures, delivery of goods in densely populated urban areas, disaster property of neural network [START_REF] Weng | Evaluating the robustness of neural networks: An extreme value theory approach[END_REF].

management, emergency services, etc. UAVs are capable of replacing humans in potentially dangerous situations, but also reduce the cost of operations in some businesses. But these kind of systems need to be continuously monitored to detect in real-time any problem ("fault") that could lead to a failure of the UAV and with a potential interruption of the mission, a potential collision with a manned aircraft (air risk), a potential collision with people on ground or with a critical infrastructure like a nuclear plant (ground risk). The need is to develop fault-tolerant systems. Formally, fault tolerance is the ability of a system to perform the intended function despite the presence or occurrence of faults, whether it is a physical damage to the hardware, software defects, or malicious attacks. In this paper, we deal with a Machine Learning based UAV fault detection (depicted in the red component of Figure 1a) that implements a Deep Neural Network (DNN) to detect mechanical faults like an elevon stucked or locked in a given position. This real-time fault detection mechanism allows the remote pilot to take necessary actions (e.g. remote triggering of the drone emergency parachute) to avoid the loss of control of the UAV which could induce an air risk or a ground risk. This Machine Learning based fault detection function should be safe to claim compliance with the applicable regulations. In particular, it should be stable with respect to noise perturbations in the input data coming from the inertial sensors and the autopilot i.e. the model continues to operate correctly despite small perturbations in its inputs within the Operational Design Domain (ODD). In other words, small perturbations that are considered as part of the real operating conditions should be considered in the designing phase of the model making the model 'stable-by-design'. However, the concept of adversarial attacks intro-duced in [START_REF] Szegedy | Intriguing properties of neural networks[END_REF] shows that adding a designed imperceptible noise to the sample is able to fool a neural network even if the model is trained to achieve high accuracy and performance. There have been studies that design different attacks [START_REF] Moosavi-Dezfooli | Deepfool: a simple and accurate method to fool deep neural networks[END_REF][START_REF] Gupta | An adversarial attacker for neural networks in regression problems[END_REF] on different model architectures and applications. Such design weakness in DNNs questions the desired stability property of neural networks and hence its deployment in mission-critical or safety-critical products like UAVs. Many works providing counter measures to adversarial attacks such as adversarial training [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] and defensive distillation have appeared in literature. But it is hard to demonstrate the completeness of these counter measures, i.e. their ability to address all types of foreseeable perturbations in the given operational context.

An example of formal measure of stability of neural network is the Lipschitz constant of the model which allows to evaluate how small perturbations in the inputs impact the output variations as illustrated in Figure 1b. A sensitivity analysis of inputs using Lipschitz constant for neural networks in presented in [START_REF] Gupta | Multivariate lipschitz analysis of the stability of neural networks[END_REF]. The contributions of this work is to demonstrate the efficiency of formal guarantees based on the Lipschitz constant for checking the stability of DNNs. We present a stability control loop grounded on spectral normalisation constraints to achieve an accuracy-stability trade-off. We show the results in various experimental designs on a small UAV dataset for classification tasks. We start with presenting generalisation capabilities of neural networks over other machine learning methods such as SVM and logistic regression, and then give Lipschitz constant estimates of different neural network configurations. We show that it is possible to train a neural network satisfying both Lipschitz and accuracy targets using our stability control loop. We also demonstrate the effectiveness of spectral normalised models in handling adversarial tasks.

The next section discusses related work in the field of fault detection using neural networks and stability analysis of these architectures. Following this, a spectral normalisation based stability control loop is described. Then the real flight data used to evaluate the fault detection model is presented. In the following section, results and observations for various scenarios related to binary and multi-class classification tasks are presented. The last section summarizes the main contributions of this work.

Related Works

Fault Detection and Diagnosis(FDD)

FDD strategies can be primarily divided into two categories model-based and data-driven. Model-based approaches require an accurate model of the aircraft for successful diagnosis of the faults. But a small UAV system is susceptible to various perturbations and lacks an accurate model. Hence, model based approaches fail in such scenarios and studies are then concentrated on using datadriven methods.

AI based data-driven techniques are gaining more and more interest in the field of fault detection. Fault detection in UAVs has been studied by using different sensors and/or for detecting faults in different fault prone parts of the vehicle. Freeman et al. [START_REF] Freeman | Model-based and data-driven fault detection performance for a small uav[END_REF] compares model-based and data-driven fault detection methods by using data from small, low-cost UAVs. Both approaches were shown to be capable of detecting different elevon faults in real data during maneuvers and in the presence of environmental variations. Guo et al. [START_REF] Guo | A hybrid feature model and deep learning based fault diagnosis for unmanned aerial vehicle sensors[END_REF] proposed a hybrid feature model and deep learning-based fault diagnosis for UAV sensors. The model detects different sensor faults, such as GPS, IMU and ADS. The authors of [START_REF] Xie | Grnn model for fault diagnosis of unmanned helicopter rotor's unbalance[END_REF] modelled fault detection using Generalized Regression Neural Network with a particle swarm optimization algorithm. Eroglu et al. [START_REF] Eroglu | Autolanding control system design with deep learning based fault estimation[END_REF] proposed a datadriven fault estimation method for estimating actuator faults from aircraft state trajectories. The authors propose to use an architecture with 1-D CNN followed by LSTM layers for learning state trajectories on simulated data of fault injected aircraft. Sadhu et al. [START_REF] Sadhu | On-board deep-learning-based unmanned aerial vehicle fault cause detection and identification[END_REF] proposed a deep learning based real-time fault detection and identification on UAV IMU sensor data using CNN and Bi-LSTM. Iannace et al. [START_REF] Iannace | Fault diagnosis for uav blades using artificial neural network[END_REF] presented a single layered neural network for detecting unbalanced blades in a UAV propeller. Bronz et al. [START_REF] Bronz | Real-time fault detection on small fixed-wing uavs using machine learning[END_REF] uses an SVM classifier with RBF kernel to classify the different faults for small fixed-wing UAVs. We use the same real flight dataset for our analysis and make comparisons with the reported results for SVM classifier.

Stability Quantification through Lipschitz Constant

The number of works studying the stability issue of neural networks has increased in manifolds with those defining and handling the stability of these networks in various manners. A way of quantifying the stability of a trained neural network model consists of studying the Lipschitz properties of the neural network. A Lipschitz constant of a multi-input multi-output function f is an upper bound on the ratio between the variations in the outputs and the variations in inputs of f , and thus it constitutes a measure of sensitivity of the function with respect to input perturbations.

For a feed-forward neural network represented by f where m is the number of hidden layers, N i the number of neurons in layer i ∈ {1, . . . , m}, let x ∈ R N0 be the input and f (x) ∈ R Nm be the associated output. If L ∈ [0, +∞[ is such that, for every input x ∈ R N0 and perturbation z ∈ R N0 , we have

∥f (x + z) -f (x)∥ ≤ L∥z∥, (1)
then L is a Lipschitz constant of the neural network f . The smaller this constant value, the more stable the neural network.

A trivial upper bound on the Lipschitz constant of a neural network was defined in [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] as the product of spectral norms (S) of each layer independently. Mathematically it can be written as

L m = ∥W m ∥ S ∥W m-1 ∥ S • • • ∥W 1 ∥ S . (2)
where W i are the weights/kernels at each layer i ∈ {1, . . . , m}. This bound was found to be too loose for some practical applications.

In [START_REF] Virmaux | Lipschitz regularity of deep neural networks: analysis and efficient estimation[END_REF], it was proved that the problem of computing the exact Lipschitz constant of a neural network with differentiable activation functions is NP-hard. In [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF], the authors proposed various bounds on the Lipschitz constant of a feedforward network under the assumption that the activation operator at each layer is α i -averaged with α i ∈]0, 1]. This assumption is satisfied by most of the widely used activation functions, which are 1/2-averaged [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF]. If the activation function at each layer is separable, a general Lipschitz constant estimate, CPLip, is defined as

L m = sup Λ1∈D1,...,Λm-1∈Dm-1 ∥W m Λ m-1 • • • Λ 1 W 1 ∥ S , (3)
where, for every i ∈ {1, . . . , m -1}, D i is the set of diagonal matrices of size N i × N i with diagonal elements equal to 2α i -1 or 1. This estimate remains true even with other norms than the Euclidean, applied to the input and output perturbation spaces. In [START_REF] Latorre | Lipschitz constant estimation of neural networks via sparse polynomial optimization[END_REF] authors proposed a polynomial constrained optimization (LipOpt) based technique for estimating Lipschitz constant. This estimate can be used for any ℓ p norm for input and output perturbations, but it is valid for neural network with single output. Solving such an optimization problem can be achieved by solving a hierarchy of convex problems such as Lasserre's hierarchy. SDP based Lipschitz Constant estimates were explored in [START_REF] Fazlyab | Efficient and accurate estimation of lipschitz constants for deep neural networks[END_REF]. This work focuses on neural networks using separable activation operators. It assumes that the activation functions at all layers are slope-bounded. By using the nonexpansiveness property of the activation function, the optimization problem of estimating a Lipschitz constant is recast as solving a semi-definite positive programming (SDP) problem. This estimation can only be used when the Euclidean norm is employed at both input and output spaces. Also, the most optimized technique for estimating the Lipschitz constant in [START_REF] Fazlyab | Efficient and accurate estimation of lipschitz constants for deep neural networks[END_REF] is erroneous, as mentioned and proved in [START_REF] Pauli | Training robust neural networks using lipschitz bounds[END_REF]. Hence, for all our results we use the valid version LipSDPneuron for all our Lipschitz constant estimates.

Spectral Normalization Control of the Model Stability

For safety critical tasks, Lipschitz constant and performance targets can be specified as engineering requirements, prior to network training. A Lipschitz constant target can be defined by a safety analysis of the acceptable perturbations for each output knowing the input ranges and it constitutes a current practice in many industrial processes. The Lipschitz constant of a neural network can be arbitrarily high when unconstrained, making it sensitive to adversarial attacks. Imposing this Lipschitz target can be done either by controlling the Lipschitz constant for each layer or for the whole network, depending on the application at hand. Such a work for controlling the Lipschitz constant has been presented in [START_REF] Serrurier | Achieving robustness in classification using optimal transport with hinge regularization[END_REF] using Hinge regularization. Spectral normalization techniques [START_REF] Miyato | Spectral normalization for generative adversarial networks[END_REF][START_REF] Pan | Fast approximate spectral normalization for robust deep neural networks[END_REF] have been proved to be very effective in controlling stability properties of DNNs and counter adversarial attacks. But they require an exact computation of the spectral norm in (3), which is computationally expensive, and therefore it becomes intractable for deep neural networks. On the other hand, it has been shown that constraining the Lipschitz constant of DNNs enhances their stability, but it makes the minimization of the loss function harder. Hence, a trade-off between stability and prediction performance needs to be reached [START_REF] Zhang | Theoretically principled trade-off between robustness and accuracy[END_REF]. We subsequently propose an iterative procedure allowing us to find the best trade-off.

Given an m-layer fully connected NN architecture (FCN) denoted by T with a Lipschitz constant L target and an accuracy target acc target . We can constrain the spectral norm of each layer to be less than m L target . According to [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF], this ensures that the upper bound on the global Lipschitz constant is less than L target . This bound is however overpessimistic, which means that, by monitoring accurately the value of the actual Lipschitz constant of the network, we can relax the bound on individual linear layers so as to get an improved accuracy. In our proposed approach, this process is done in a closed loop with two targets : good performance (acc target ) and desired stability L target .

The approach is summarized in Algorithm 1. Module Train SN (L, . . .) is spectral normalisation training module which trains the defined neural network model with the Lipschitz constant constraint, i.e. each weight matrices W i with i ∈ {1, . . . , m} of the neural network is divided by (σ i / m √ L), where σ i is the largest singular value of W i . Module LipEst estimates accurately the Lipschitz constant of the trained neural network using LipSDP-Neuron [START_REF] Fazlyab | Efficient and accurate estimation of lipschitz constants for deep neural networks[END_REF]. Module Predict checks the performance (here classification performance) of the trained model on the validation dataset. If L is small, training will be over-constrained and the accuracy of the model will suffer. To relax this constraint we iteratively increase the value of L through small α factor to find the best value of L opt which would be less than L target and ensure an increase in the performance on the validation dataset to achieve acc opt greater than acc target . We stop the iterations when both conditions are met i.e. accuracy of the trained model is greater than acc target and Lipschitz constant is less than L target .

Validation on Drone Flights Data

Dataset Description

For our fault detection and stability analysis, we use the real-flight dataset provided in [START_REF] Bronz | Real-time fault detection on small fixed-wing uavs using machine learning[END_REF]. The authors use Paparazzi Autopilot system [START_REF] Hattenberger | Using the Paparazzi UAV System for Scientific Research[END_REF] for possible flight trajectories. The flight test data is captured in outdoor setting on different days with varying environments especially wind speeds and also imperfections in the geometry due to manufacturing of the UAV.

The actuator fault model (F ) as proposed in [START_REF] Tandale | Fault-tolerant structured adaptive model inversion control[END_REF] is given by

u app = F u com + E (4)
where u com is the commanded control deflection from the autopilot to the actuators and u app is the applied control deflection (i.e. final movement of the elevon) and E is the offset. In vector form, it can be expressed as where () r and and () l stands for the right and left elevon respectively. The feature set is a vector of length 8 consisting of linear accelerations at three coordinates (a x , a y , a z ), angular rates (ω x , ω y , ω z ), and auto-pilot commanded controls (u comr , u com l ) for the two aerodynamic actuators. Depending on the values of f r and f l , the fault detection problem can be translated into two kinds of classification problems: 1) Binary Class Classification and 2) Multi-class Classification.

u appr u app l = f r 0 0 f l u comr u com l + e r e l , (5)
We use three performance measures: classification accuracy (%), F1 score, and Matthews Correlation coefficient (MCC). For stability measure we use Lipschitz constant estimate as calculated using LipSDP-Neuron [START_REF] Fazlyab | Efficient and accurate estimation of lipschitz constants for deep neural networks[END_REF]. We recall that we do not use LipSDP most tight bound LipSDP-Network since it is erroneous as mentioned in the section 2.2.

Binary Classification Results

We train the data on 12 th July and test the trained model on 12 th july unseen and 13 th July data. It is interesting to use 12 th July flight data for training and predict fault for the next day, i.e 13 th July since 12 th and 13 th July had different atmospheric conditions at the time of data collection. For further checking the efficacy and generalization capabilities of the methods, we checked the trained models on flight data for 21 st July. The faults were injected in the right elevon f r = 0.3. Fault code 0 implies 'Nominal fault' and fault code 1 implies fault in the right elevon. Implementation details: We propose to use fully connected feed-forward networks (FCN) for solving the problem. We train such neural networks of varying depth and width. We use ReLU as the activation function except for the last layer where it is employed with sigmoid for binary and softmax for multi-class classification. We use the cross entropy losses as the loss function for training the model. The feature set is standardized by removing the mean and scaling to unit variance. While training, the dataset is split into 4:1 ratio for train and validation sets. The models are trained using Keras with Tensorflow backend. The initializers are set to Glorot uniform and we use Adam as the optimizer. To capture the temporal dynamics of the system, feature set of length 8 (a x , a y , a z , ω x , ω y , ω z , u comr , u com l ) is concatenated T times. Therefore, the input feature vector to the model is of length 8T . We observe a general increasing trend in the testing performance as we increase the value of T , which suggests that concatenating more feature lists from the time history of the flight data improves the performance of the real-time prediction. Hence, we found the value T = 20 as optimal for our experiments, adding time history beyond this level did not contribute to the learning and the neural network starts to overfit.

In the first set of experiments for binary classification, we train on 12 th July flight data and test on of 12 th July, 13 th , and 21 st July and directly compare our NN best models to SVM and logistic regression. SVM classifier is trained with a RBF kernel and is optimized using grid-search. The results are given in Table 1. The results obtained using SVM classifier are satisfactory, better than for logistic regression, but they fail to generalize on the unseen data from 13 th and 21 st July. We remark it is possible to train neural networks with much better generalization performance on different days with varying weather conditions. In second set of experiments we discuss the stability performance of neural networks which is our primary concern pertaining to stable model design. The comparison of the classification task accuracy as well as the stability performance in terms of Lipschitz constant for a varying number of hidden layers and hidden neurons is shown in Table 2. Notation [START_REF] Gupta | An adversarial attacker for neural networks in regression problems[END_REF][START_REF] Freeman | Model-based and data-driven fault detection performance for a small uav[END_REF] in Table 2 implies that the FCN has two hidden layers with 10 and 6 hidden neurons and so on. In the baseline training, no constraints are employed while training the NN. Lipschitz constant (L) is computed using LipSDP-neuron [START_REF] Fazlyab | Efficient and accurate estimation of lipschitz constants for deep neural networks[END_REF] on the trained network.

Third, we show the performance and stability effect of constraining the training of the neural network by spectral normalization, as introduced in section 3 with Algorithm 1. We chose the same configuration for the neural networks as baseline and introduced constraints, by setting L target = 4 and acc target = 75%. Table 3a shows a comparison of the classification accuracy and Lipschitz constant value for baseline training and the new training procedure described in Algorithm 1. We remark that the proposed approach improves both accuracy and stability of neural network models.

We also test and compare the results of the trained neural networks when they are attacked with widely known adversarial attacks such as FGSM [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] and PGD [START_REF] Madry | Towards deep learning models resistant to adversarial attacks[END_REF]. The attacks are generated using ART toolbox 3 . An adversarial attack is a special case of addition of noise to the input samples. It is useful to test robustness of the neural network since the adversarial noise is crafted in a manner to generate misclassification errors. We show the results on 12 th July test dataset in Table 4 4 . We remark that spectral normalised training is better in handling the adversarial attacks, leading to more stable models.

Multi-class Classification Results

We use the flight data captured on the days of 21 st and 23 rd July for detecting the faults in right and left elevon. For injecting the faults, efficiency of the right elevon is reduced (f r = 0.3) and similarly for the left control elevon (f l = 0.9/0.3). Fault codes: '0' -nominal , '1'f r = 0.3 , '2'f l = 0.9 and '3'f l = 0.3, making a 4-class classification problem. The emphasis has been put on detecting faults on the left elevon which has the effect of a geometric twist. Identical faults are injected in the flight data for both days. For capturing the temporal dynamics in the prediction performance, the optimal value of T is taken to be 20 as in the previous section. For the experiments, we train a FCN with ReLU activation function for all layers, except the last one for which softmax is used. We train the neural network on the flight data of 21 st July and test on data of 21 st July not seen by the model at the time of training, and also on the flight data of 23 rd July. We show the prediction and generalization capabilities of the neural networks over state-of-the-art SVM and logistic regression. We tabulate the classification metrics in Table 5.Varying configurations of model layers and neurons classification metrics along with the stability metric (L) have been considered in Table 6. We consider acc target = 70% and L target = 4 as the model performance and stability requirements. The results are shown in Table 3b.

Observations of the performance and stability analysis:

• Neural networks perform better than SVM classifier and logistic regression in terms of prediction performance indicating better generalization capabilities of neural network models, when trained and tested on different days. • We observe the trend of increased Lipschitz constant as the number of layers and neurons in the neural network increases. • There is a steep decrease in the Lipschitz constant value when the models are trained with a Lipschitz target which satisfies the performance target as well Table 3. We also observe an increase of classification performance with spectral normalization constraint. • Spectral normalised trained models are better against adversarial attacks implying since they are designed to be more stable.

Conclusion

Our work provides a formal method by means of estimating Lipschitz constant to evaluate and control the stability of a neural network, which is necessary for many mission-critical and safety-critical systems. This work analyses a real flight data collected on a small fixed-wing UAV in varying flight conditions to design a fully connected feed-forward neural networks that implements a fault detection function embedded in the UAV. The task is to detect mechanical faults like an elevon stucked or locked in a given position even in the presence of noise in the input data coming from the inertial sensors and the autopilot and to allow the remote pilot to take necessary actions to ensure the safety. The main contributions of this work are the following: Through the evaluation of classification performance we have shown that our approach based on DNN compares favorably with SVM and logistic regression techniques. Unlike many works related to the use of neural networks in industrial applications, we address the stability of our trained neural networks using a formal verification method based on the quantification of the Lipschitz constant of the model. We show the effectiveness of spectral normalization techniques in controlling the Lipschitz constant of the network to reach a Lipschitz safety target which can be defined at system level. Following the insights presented in the paper, we show that it is possible to build neural networks which are "Stable-by-Design" and verified through Lipschitz certificates. Such designed models are more robust to adversarial attacks and present a good accuracy. These results regarding stability are strong developments toward learning assurance artefacts in ML models.

Fig. 1 :

 1 Fig. 1: (a) Fault Detection in drones/UAVs :The component called "Fault Detection Component" computes a flight condition status (nominal or faulty) using the outputs of the inertial sensors and the elevon deflection control variable produced by the autopilot. (b) Intuition of using Lipschitz constant as a stability property of neural network [24].

  Algorithm 1 Spectral Normalisation-based Stability Control Loop 1: Inputs : training dataset (Xtrain, Ytrain), validation dataset (X valid , Y valid ), target Lipschitz constant Ltarget, target accurcay acctarget, neural network T , multiplicative factor α > 1, maximum iteration number nmax 2: Output : optimal neural network Topt, Lipschitz constant Lopt, and accuracy accopt

	3: Algo :
	4: L ← Ltarget, n ← 0
	5: T ← TrainSN(L, Xtrain, Ytrain)
	6: Lest ← LipEst(T )
	7: acc ← Predict(T, X valid , Y valid )
	8: while (n < nmax) and (Lest < Ltarget) and (acc < acctarget) do
	9:	L ← α L
	10:	T ← TrainSN(L, Xtrain, Ytrain)
	11:	Lest ← LipEst(T )
	12:	acc ← Predict(T, X valid , Y valid )
	13:	n ← n + 1
	14: end while
	15: Topt ← T
	16: Lopt ← Lest
	17: accopt ← acc
	18: if n = n,max then
	19:	Notification: No trade-off found between Stability and Accuracy
	20: else
	21:	Notification: Trade-off found between Stability and Accuracy
	22: end if

Table 1 :

 1 Binary Classification : Comparison between SVM, logistic regression, and best NN trained with baseline training on 12 th July.

	Model	12 th July Acc. (%) F1 MCC Acc. (%) F1 MCC Acc. (%) F1 MCC 13 th July 21 st July
	SVM	98.7	0.99 0.98	60.5	0.59 0.38	85.3	0.85 0.72
	LogReg	89.8	0.90 0.79	58.0	0.59 0.31	56.6	0.45 0.18
	Baseline	98.6	0.99 0.97	73.3	0.74 0.52	93.1	0.93 0.86
	Model	12 th July Acc. (%) F1 MCC Acc. (%) F1 MCC 13 th July	L
		(10)	98.3	0.983 0.967	66.4	0.67 0.41 12.26
	(10,6)	98.4	0.985 0.969	71.0	0.71 0.47 27.56
	(50,10,6)	98.7	0.987 0.974	73.3	0.74 0.52 37.76

Table 2 :

 2 Binary Classification : Results on 12 th and 13 th July using NN trained with baseline training on 12 th July for various configurations of FCN.

  Results on 13 th July trained on 12 th July data. July trained on 21 st July data.

	Model (10) (10,6) (50,10,6) (a) Binary Classification: Model Baseline Spectral Acc. (%) L Acc. (%) L 66.4 12.26 75.1 3.7 71.0 27.56 75.3 3.8 73.3 37.76 76.7 3.7 (50,10) (100,10) (100,75,10) (b) Multi-class Classification: Results Baseline Spectral Acc. (%) L Acc. (%) L 69.7 31.44 72.65 3.8 69.54 34.17 71.07 3.6 68.9 69.06 71.93 3.7
	on 23 rd

Table 3 :

 3 Comparison between models trained with baseline and spectral normalization constraint (L target = 4).

	Model Training	No attack Acc. (%) MCC Acc. (%) MCC Acc. (%) MCC FGSM (ϵ = 0.2) PGD(ϵ = 0.2)
	(50,10,6)	Baseline Spectral	98.6 98.7	0.97 0.97	77.6 81.0	0.62 0.66	76.0 79.2	0.59 0.62

Table 4 :

 4 Binary Classification: Results on 12 th July clean data, FGSM and PGD attacks using model trained on 12 th July.

Table 5 :

 5 Multi-class Classification: Comparison between SVM, logistic regression and best NN trained with baseline training on 21 st July and tested on 21 st and 23 rd July .

	Model	21 st July Acc. (%) F1 MCC Acc. (%) F1 MCC 23 rd July	L
	(50,10)	95.2	0.95 0.93	69.7	0.68 0.58 31.44
	(100,10)	95.3	0.95 0.93	69.5	0.68 0.57 34.17
	(100,75,10)	95.4	0.95 0.93	68.9	0.67 0.57 69.06

Table 6 :

 6 Multi-class Classification : Results on 21 st July and 23 rd July using NN trained with baseline training on 21 st July for various configurations of FCN.

https://github.com/Trusted-AI/adversarial-robustness-toolbox

In case of SVM model, classification metric MCC had negative values for FGSM and PGD attacks implying model performs very poorly against adversarial attacks, hence we discard them from comparison.