Keywords: réseaux de neurones, attaques adverses, constante de Lipschitz Stability Quantification of Neural Networks Robustness, stability, neural networks, adversarial attacks, constant A.1 Results on Combined Cycle Power Plant Data

Quantification de la stabilité des réseaux de neurones Mots clés: robustesse, stabilité,

Abstract: This thesis contains three main contributions concerning the stability of neural networks. First, we present a simple and easily implementable white box adversarial attacker for regression models operating on tabular data. We rely upon the Jacobian properties of the neural networks to generate the worst attack on the sample. We propose three metrics for quantifying the success of the attacks analogous to metrics used in classification scenarios. In our second contribution, we propose a multivariate sensitivity analysis of the neural network inputs based on its Lipschitz properties. Such a study helps to understand better the effect of individual inputs on the dynamics of the outputs. The presented method works on individual as well as any group of inputs. We also propose a "Lipschitz star" representation to display this sensitivity analysis of the inputs. In the last contribution, we design a stability control method to achieve the best trade-off between accuracy and stability during the training phase of neural network models. This method is based on a spectral normalization inserted in a training loop that allows us to achieve given performance and stability targets. We validate this approach on an Unmanned Aerial Vehicle (UAV) application. To support our contributions, we test them both on open-source, public datasets and an industrial dataset provided by Thales LAS France that is related to the domain of safety-critical applications.

Abstract

Artificial neural networks are at the core of recent advances in Artificial Intelligence. One of the main challenges faced today, especially by companies like Thales designing advanced industrial systems is to ensure the safety of new generations of products using these technologies. In 2013 in a key observation, neural networks were shown to be sensitive to adversarial perturbations, raising serious concerns about their applicability in critically safe environments. In the last years, publications studying the various aspects of this robustness of neural networks, and rising questions such as "Why adversarial attacks occur?", "How can we make the neural network more robust to adversarial noise?", "How to generate stronger attacks?" etc., have grown exponentially. The contributions of this thesis aim to tackle such problems. The adversarial machine learning community concentrates majorly on classification scenarios, whereas studies on regression tasks are scarce. Our contributions bridge this significant gap between adversarial machine learning and regression applications. Also, all the datasets used as part of the contributions are tabular, which has not been highly researched by the deep learning community and only recently received more attention.

The first contribution in Chapter 3 proposes a white-box attackers designed to attack regression models. The presented adversarial attacker is derived from the algebraic properties of the Jacobian of the network. We show that our attacker successfully fools the neural network and measure its effectiveness in reducing the estimation performance. We present our results on various open-source and real industrial tabular datasets. Our analysis relies on the quantification of the fooling error as well as different error metrics. Another noteworthy feature of our attacker is that it allows us to optimally attack a subset of inputs, which may help to analyze the sensitivity of some specific inputs. We also, show the effect of this attacker on spectrally normalised trained models which are known to be more robust in handling attacks.

The second contribution of this thesis (Chapter 4) presents a multivariate Lipschitz constant analysis of neural networks. The Lipschitz constant is widely used in the literature to study the internal properties of neural networks. But most works do a single parametric analysis, which do not allow to quantify the effect of individual inputs on the output. We propose a multivariate Lipschitz constantbased stability analysis of fully connected neural networks allowing us to capture the influence of each input or group of inputs on the neural network stability. Our approach relies on a suitable re-normalization of the input space, intending to perform a more precise analysis than the one provided by a global Lipschitz constant. We display the results of this analysis by a new representation designed for machine learning practitioners and safety engineers termed as a Lipschitz star. We perform experiments on various open-access tabular datasets and an actual Thales Air Mobility industrial application subject to certification requirements.

The use of spectral normalization in designing a stability control loop is discussed in Chapter 5. A critical part of the optimal model is to behave according to specified performance and stability targets while in operation. But imposing tight Lipschitz constant constraints while training the models usually leads to a reduction of their accuracy. Hence, we design an algorithm to train "stable-by-design" neural network models using our spectral normalization approach, which optimizes the model by taking into account both performance and stability targets. We focus on Small Unmanned Aerial Vehicles (UAVs). More specifically, we present a novel application of neural networks to detect in real-time elevon positioning faults to allow the remote pilot to take necessary actions to ensure safety.

Résumé (French)

Les réseaux de neurones artificiels sont au coeur des avancées récentes en Intelligence Artificielle. L'un des principaux défis auxquels on est aujourd'hui confronté, notamment au sein d'entreprises comme Thales concevant des systèmes industriels avancés, est d'assurer la sécurité des nouvelles générations de produits utilisant cette technologie. En 2013, une observation clé a révélé que les réseaux de neurones sont sensibles à des perturbations adverses. Ceci soulève de sérieuses inquiétudes quant à leur applicabilité dans des environnements où la sécurité est critique. Au cours des dernières années, des publications ont étudiées les différents aspects de la robustesse des réseaux de neurones, et des questions telles que "Pourquoi des attaques adverses se produisent?", "Comment pouvons-nous rendre les réseaux de neurones plus robustes à ces bruits ?", "Comment générer des attaques plus fortes ?", etc., se sont posées avec une acuité croissante. Cette thèse vise à apporter des réponses à de telles questions. La communauté s'intéressant aux attaques adverses en apprentissage automatique travaille principalement sur des scénarios de classification, alors que les études portant sur des tâches de régression sont rares. Nos contributions comblent le fossé existant entre les méthodes adverses en apprentissage et les applications de régression. De plus, tous les ensembles de données utilisés dans le cadre de nos travaux sont tabulaires. De telles données ont fait l'objet de peu de recherches approfondies dans la communauté de l'apprentissage et n'ont reçu que récemment une attention plus soutenue.

Notre première contribution, dans le chapitre 3, propose un algorithme de type "boîte blanche" pour attaquer les modèles de régression. L'attaquant adverse présenté est déduit des propriétés algébriques du Jacobien du réseau. Nous montrons que notre attaquant réussit à tromper le réseau de neurones et évaluons son efficacité à réduire les performances d'estimation. Nous présentons nos résultats sur divers ensembles de données tabulaires industriels en libre accès et réels. Notre analyse repose sur la quantification de l'erreur de tromperie ainsi que différentes métriques. Une autre caractéristique remarquable de notre algorithme est qu'il nous permet d'attaquer de manière optimale un sous-ensemble d'entrées, ce qui peut aider à identifier la sensibilité de certaines d'entre elles. Nous montrons également l'effet de cet attaquant sur des modèles entraînés avec une normalisation spectrale, connus pour être plus robustes aux attaques.

La deuxième contribution de cette thèse (Chapitre 4) présente une analyse de la constante de Lipschitz multivariée des réseaux de neurones. La constante de Lipschitz est largement utilisée dans la littérature pour étudier les propriétés intrinsèques des réseaux de neurones. Mais la plupart des travaux font une analyse mono-paramétrique, qui ne permet pas de quantifier l'effet des entrées individuelles sur la sortie. Nous proposons une analyse multivariée de la stabilité des réseaux de neurones entièrement connectés, reposant sur leur propriétés Lipschitziennes. Cette analyse nous permet de saisir l'influence de chaque entrée ou groupe d'entrées sur la stabilité du réseau de neurones. Notre approche repose sur une re-normalisation appropriée de l'espace d'entrée, visant à effectuer une analyse plus précise que celle fournie par une constante de Lipschitz globale. Nous visualisons les résultats de cette analyse par une nouvelle représentation conçue pour les praticiens de l'apprentissage automatique et les ingénieurs en sécurité appelée "étoile de Lipschitz". Nous menons des expérimentations sur différents ensembles de données tabulaires en libre accès et sur une application industrielle réelle de Thales Air Mobility soumise à des exigences de certification.

L'utilisation de la normalisation spectrale dans la conception d'une boucle de contrôle de stabilité est abordée au chapitre 5. Une caractéristique essentielle du modèle optimal consiste à satisfaire aux objectifs de performance et de stabilité spécifiés pour le fonctionnement. Cependant, contraindre la constante de Lipschitz lors de l'apprentissage des modèles conduit généralement à une réduction de leur précision. Par conséquent, nous concevons un algorithme permettant de produire des modèles de réseaux de neurones "stable dès la conception" en utilisant une nouvelle approche de normalisation spectrale, qui optimise le modèle, en tenant compte à la fois des objectifs de performance et de stabilité. Nous nous concentrons sur les petits drones aériens (UAV). Plus précisément, nous présentons une nouvelle application des réseaux de neurones pour détecter en temps réel les défauts de positionnement d'aileron, afin de permettre au télépilote de prendre les mesures nécessaires pour assurer la sécurité du vol.

Acknowledgement

First of all, I would like to thank my supervisors, Prof. Jean-Christophe Pesquet, Dr. Fateh Kaakai, and Dr. Beatrice Pesquet-Popsescu, for giving me this opportunity. I have learned a lot from each of them over the years, things I will carry forward in my research career and would make me a better researcher. I want to extend my thanks to my thesis examiners Prof. Hamid Krim and Prof. Mathieu Serrurier, for going through my thesis manuscript and giving valuable feedback, as well as my jury members Prof. Frederic Jurie and Dr. Juliette Mattioli, for their time and effort towards evaluating my thesis work. I want to thank all the CVN faculty for building a lab research atmosphere -Prof. Fragkiskos Malliaros, Prof. Marc Castella, and Prof. Hughes Talbot. I also want to thank Jana for helping me immensely with administrative tasks and making life smoother. Also, I would like to thank Prof. Angshul Majumdar for encouraging to pursue research as career and teaching me everything. I can't thank enough my friends Sanjay and Preethi for making my survival easy at the plateau and guiding me throughout for everything. I am really grateful to them. I want to thank Hyrule and Freja (cats) for helping with timely mental health support and just making life beautiful, and also to the people who live in their house Younes and Solen for being so welcoming, kind and generous. I want to thank Aymen for being my translator and Jhony for being so kind and motivating. I would like to thank my cite friends Surabhi, Shalu and Pragya for being such a great female support in Paris and also for the great food everytime and to Mihir for intently listening to my rants/stories and always providing with the best advices and solutions. I want to thank my parents and brother for believing and not giving up on me throughout the Ph.D. journey. Their blessings and support miles apart kept me going in my low phases. I want to thank Sagar for being my constant support system, chief advisor for everything, biggest critic, and biggest cheerleader. He has encouraged me to be a better researcher and made me a stronger and braver person. This thesis is the manifestation of my late grandmother, who used to address me as a Dr. from when I was two years old. I hope I made her happy and proud.

A Lipschitz Constant Analysis of Tabular Data

A.

Chapter 1 Introduction

. Context and Motivations

Neural networks (NN) are at the core of recent advances in Artificial Intelligence (AI). One of the main challenges faced today, especially by companies like Thales designing advanced industrial systems, is ensuring the safety of new products using these technologies. In 2013, neural networks were first shown to be sensitive to adversarial perturbations [START_REF] Szegedy | Intriguing properties of neural networks[END_REF], raising serious concerns about their suitability of use in mission-critical or safety-critical products and leading to the development of a new field of study concerning the robustness of neural networks.

In the literature, various classes of AI desirable properties have been defined by adding some objectives, components, or constraints to AI such as robustness, ethics, trustworthiness aiming at making neural networks more reliable for realworld applications as described in Table 1.1.

Concept Definition

Robust

AI systems with the ability to cope with errors during execution and cope with erroneous input Ethical AI systems that do what is right, fair, and just and prevent harm.

Trustworthy

AI systems that achieve their full potential if trust can be established in the development, deployment, and use [START_REF] Thiebes | Trustworthy artificial intelligence[END_REF] Fair AI systems absent from any prejudice or favoritism toward an individual or a group based on their inherent or acquired characteristics [START_REF] Mehrabi | A survey on bias and fairness in machine learning[END_REF] Safe AI systems deployed in ways that do not harm humanity Dependable AI systems that focus on reliability, verifiability, explainability, and security [START_REF] Singh | Trustworthy AI[END_REF] Human-centered AI systems that are "continuously improving because of human input while providing an effective experience between human and robot"

Socially Responsible

AI systems that follow ethical and philanthropic responsibilities alongwith legal and performance responsibilities [START_REF] Cheng | Socially responsible AI algorithms: Issues, purposes, and challenges[END_REF]. Among these different classes of AI properties, the contributions presented in this thesis come under the category of Robust AI. In particular, we deal with errors (wrong predictions) during execution due to the inability of the neural network to cope with small perturbations applied to nominal inputs within the Operational Design Domain (ODD). These perturbations can result from natural system degradation, such as aging sensors, or a deliberate and malicious attack (cyber crimes).

In Szegedy et al. [START_REF] Szegedy | Intriguing properties of neural networks[END_REF] the concept of adversarial attacks was first proposed to fool DNNs. Adding a subtle perturbation to the input of the neural network produces an incorrect output, while human eyes cannot recognize the difference in the modification of the input data. For example, after adding a slight perturbation to the original image labelled as "panda", it is classified as a "gibbon" by the same model with 99.3% confidence, while the human eyes cannot distinguish the differences between the original image and the adversarial image in Figure 1.1 [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF]. This behavior of DNNs is possible even when the model is well-trained with good accuracy. The capability of handling adversarial noise is known as the robustness of the Neural Network. Even though different models have different architectures and might use different training data, the same kind of adversarial attack strategies can attack related models (transferability of Adversarial attacks). These attacks pose a huge threat to the performance of DNNs and cannot be ignored in safetyand security-critical AI applications. Analyzing the cause of adversarial attacks can help researchers to fix the vulnerability effectively, to provide defenses, and to guide in generating more robust and safe models. . In the recent years, the number of works devoted to understanding the instability and increasing robustness of NNs have grown exponentially. Many approaches have been proposed, more or less dedicated to some specific architectures, particularly networks using only ReLU activation functions, and grounded on more or less empirical techniques [START_REF] Virmaux | Lipschitz regularity of deep neural networks: analysis and efficient estimation[END_REF][START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF]. Broadly we can classify the works in literature as follows:

i) Purely computational approaches consist of attacking a neural network and observing its response to such attacks [START_REF] Kurakin | Adversarial examples in the physical world[END_REF][START_REF] Farnia | Generalizable adversarial training via spectral normalization[END_REF].

ii) Methods based on heuristics for testing / promoting the robustness of a neural net [START_REF] Wong | Provable defenses against adversarial examples via the convex outer adversarial polytope[END_REF][START_REF] Raghunathan | Certified defenses against adversarial examples[END_REF][START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF].

iii) Studies that aim to establish mathematical certificates and provide theoretical robustness guarantees [START_REF] Leino | Globally-robust neural networks[END_REF][START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF][START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF].

The three strategies might ultimately be necessary to effectively build robust neural networks and certify them. However, we believe the only one providing strong robustness guarantees is the last one. Hence, we concentrate more on providing theoretical guarantees. In the next chapter, we give details on various state-of-the-art methods used for verifying and promoting robustness of NN.

slope = -L slope = L f(x) x x x+z f(x) f(x+z) f(x)+L||z|| p f(x)-L||z|| p Figure 1
.2: Intuition of using Lipschitz constant as a stability property of neural network [START_REF] Weng | Evaluating the robustness of neural networks: An extreme value theory approach[END_REF].

Among the possible mathematical approaches, we focus on those relying upon the analysis of the Lipschitz properties of neural networks. The Lipschitz modulus of any function is defined as the upper bound on the ratio of variations in the output to the variations in the output. In the context of neural networks, it can be seen as a slope of the changes in a high dimensional space. An intuition of using Lipschitz bounds for any neural network function f is given in Figure 1.2. Very high changes or steep slope values can be sign of highly unstable models, hence not robust to input perturbations. Although neural networks not being robust can be attributed to various factors, such properties play a fundamental role in understanding the internal mechanisms governing these complex non-linear systems. Besides, these methods make very few assumptions on the type of nonlinearities used and are thus valid for a wide range of networks. We base our work and contributions on estimating and controlling the Lipschitz bounds of neural networks and hence enable engineers to develop "Stable-by-design" models which are safe under deployment phase.

. Contributions of the Thesis

The contributions presented in this thesis aim to bridge the gaps encountered in adversarial machine learning literature especially when applied to safety-and mission-critical systems. One major limitation of the existing works is related to information on the nature of adversarial attacks in regression tasks. The literature and context of adversarial attacks and their defenses move around classification tasks and notion of clear boundaries between the classes. Technically, adversarial attacks occur for neural networks when the logit of "wrong/target" class crosses the "true" class. Similarly, the defenses are designed or models are trained so that the logit of any class will be less likely to cross the logit of "true" class. More specifically, these notions are modelled using cross-entropy losses which are primarily used for classification scenarios. Regression applications cover a major chunk of industrial applications, but the adversarial perturbations are not defined as easily in this context. Also, the first industrial application presented in the thesis is a regression problem. To bridge this major gap in the literature, we provide a simple "white" box attacker for generating adversarial attacks for regression tasks. We define some metrics similar to those used in classification scenarios for checking the effectiveness of attacking the neural networks. It is worthy to note here that adversarial attacks in the literature mostly correspond to ℓ ∞ or ℓ 2 norm while the attacker designed as part of our contribution can work for various combinations of norms on input and output space. Also, we provide a way to attack only a group of inputs at a time.

The work on estimation of Lipschitz constant in the literature aims to find a single Lipschitz constant parameter that does not highlight individual contributions of the inputs in the prediction of the output. In other words, they do not provide any information on the sensitivity of the inputs. Such analysis is important so that engineers can give importance to control the perturbations arising in certain sensors over others. To assess this information on the sensitivity of the inputs, our next contribution provides a multivariate Lipschitz analysis of the neural networks by introducing the concepts of partial/weighted Lipschitz constants. As part of our contribution, we provide a representation called Lipschitz star, a graphical and practical tool to analyze the sensitivity of a neural network model during it's development, concerning different combinations of inputs.

The Lipschitz constrained neural networks often perform poorly on the clean samples if constrained to have very tight Lipschitz bounds. Such poor performance is not acceptable for a working model. For achieving the most optimally trained models it is necessary to attain a trade-off between the clean accuracy and the Lipschitz constant value. To achieve this trade-off in our next contribution, we designed a spectral normalized stability control loop which takes care of both predefined stability and performance targets.

To support our contributions, we have validated and tested various open-source and real world applications, including a use-case provided by Thales LAS France for air traffic management and another one related to fault detection in UAVs.

. Organisation of the Thesis

This thesis is organised as follows.

i) In Chapter 2 we provide a literature review of the different strategies used to improve neural networks' stability. We have divided this chapter primarily into two parts : (a) Estimation of stability and (b) Controlling the stability of neural networks. Since the main focus of this work is based on a Lipschitz constant analysis we provide an in-depth study of works utilizing the Lipschitz properties of neural networks to improve the robustness.

ii) In Chapter 3 we focus on designing an adversarial attacker for regression tasks which are somewhat neglected in the adversarial machine learning community. We show that by utilizing the Jacobian properties of neural networks it is possible to design an attacker which generates worst attack on the sample. To achieve this, we provide some metrics for quantifying the effectiveness of the presented attacker. We show the effect of our proposed attacker on standard trained model and spectral normalized trained models which have much lower Lipschitz bounds.

iii) In Chapter 4, we present a multivariate Lipschitz analysis of neural networks. We provide the concept of partial/weighted Lipschitz constant for quantifying the sensitivity of individual inputs w.r.t to outputs. We present a new representation "Lipschitz star" to graphically represent the introduced partial Lipschitz constant. We study the effect of spectral normalization method and adversarial learning on the sensitivity of each input. Note: For this thesis we have used words stability and robustness interchangeably and they thus represent the same concept i.e. neural network's ability to handle any kind of perturbations and adversarial attacks. The stability criterion considered here highlights the fact that small perturbations in the inputs do not produce high variations of the outputs. It should be noted that the safety community has defined robustness and stability of models as two different notions. The robustness of the model is the extent to which the system can continue to operate correctly despite abnormal inputs and conditions outside the defined Operational Design Domain (ODD). In contrast, stability is the network's ability to handle small perturbations within its inputs ODD 1 . The research community working on Robust AI currently does not differentiate between the two concepts.

. List of Publications

Chapter 2 A Survey on the Stability of Neural Networks

In this chapter we provide a literature review for various works in Robust AI domain related to verification and certification of neural networks. We concentrate on the methods utilizing the Lipschitz properties of the neural network as a formal guarantee of robustness.

. Provable Defenses

Neural networks (NNs) are being deployed in a wide range of applications, including safety-and mission-critical tasks, certifying robustness of a neural network against perturbations has become an important research topic in machine learning community. There are various approaches and strategies to test the robustness of neural networks. Traditional coverage-based approaches may be irrelevant for testing neural network systems. In such cases, code certifiability can be trivially satisfied while providing only limited guarantees on the system safe behavior when facing situations that have not been strictly met during the training process.

Deep Neural Network(DNN) verification is experimentally beyond the capabilities of tools such as Linear Programming (LP) solvers or Satisfiability Modulo Theories (SMT) solvers [START_REF] Pulina | Challenging SMT solvers to verify neural networks[END_REF][START_REF] Bastani | Measuring neural net robustness with constraints[END_REF]. By far, these dedicated tools have only been able to handle very small networks (a single hidden layer with 10/20 neurons [START_REF] Pulina | Challenging SMT solvers to verify neural networks[END_REF][START_REF] Pulina | An abstraction-refinement approach to verification of artificial neural networks[END_REF]). The difficulty in using and estimating robustness properties for DNNs is the consequence of the presence of activation functions at various parts of the architecture. An attempt towards property verification of neural networks with ReLU Activations in safety critical applications was made in [START_REF] Katz | Reluplex: An efficient SMT solver for verifying deep neural networks[END_REF]. The idea is to find the exact value of minimum perturbation distance (say ∆) i.e. the trained classifier is safe against any perturbations with norm less than ∆. Finding this minimum distortion/perturbation of adversarial examples with ReLU activations is known to be an NP-complete problem [START_REF] Sinha | Certifying some distributional robustness with principled adversarial training[END_REF][START_REF] Katz | Reluplex: An efficient SMT solver for verifying deep neural networks[END_REF]. This makes formal verification technique in Reluplex [START_REF] Katz | Reluplex: An efficient SMT solver for verifying deep neural networks[END_REF] computationally demanding even for small-sized NNs and hence they suffer from scalability issues. In another line of work, the piece-wise linear nature of ReLU is exploited to efficiently compute a non-trivial certified lower bound of the minimum distortion [START_REF] Wong | Provable defenses against adversarial examples via the convex outer adversarial polytope[END_REF][START_REF] Raghunathan | Certified defenses against adversarial examples[END_REF][START_REF] Weng | Towards fast computation of certified robustness for relu networks[END_REF]. Zhang et al. [START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF] introduced an approximation technique called CROWN for computing a certified lower bound of minimum adversarial distortion given any input data point with general activation functions beyond ReLU such as tanh, sigmoid, and arctan for larger networks.

Another popular certifiable defense against ℓ ∞ norm bounded inputs via the convex outer adversarial polytope [START_REF] Wong | Provable defenses against adversarial examples via the convex outer adversarial polytope[END_REF] which is applicable to small networks with ReLU activation functions. Mirman et al. [START_REF] Mirman | Differentiable abstract interpretation for provably robust neural networks[END_REF] takes a similar approach via abstract interpretation. These methods use linear relaxations of neural networks to compute an outer bound at the final layer. However, because the convex relaxations employed are relatively expensive, these methods are typically slow to train. A simple and fast certifiable defense for ℓ ∞ norm bounded inputs is Interval Bound Propagation (IBP) [START_REF] Gowal | On the effectiveness of interval bound propagation for training verifiably robust models[END_REF][START_REF] Mirman | Differentiable abstract interpretation for provably robust neural networks[END_REF]. In [START_REF] Gowal | On the effectiveness of interval bound propagation for training verifiably robust models[END_REF] authors demonstrated that IBP bound can be quite loose for general networks. But it can be used to train large provable NNs by using appropriate loss and clever hyper-parameters to allow the network to adapt such that the IBP bound is tight. CROWN-IBP [START_REF] Zhang | Towards stable and efficient training of verifiably robust neural networks[END_REF] outperforms previous methods by combining IBP bound in a forward bounding pass and a tighter linear relaxation bound in a backward bound pass. Shi et al. [START_REF] Shi | Fast certified robust training with short warmup[END_REF] improved IBP with better initialization to accelerate training. The certification methods discussed so far provide deterministic robustness guarantees. Additionally, randomized smoothing [START_REF] Cohen | Certified adversarial robustness via randomized smoothing[END_REF][START_REF] Li | Certified adversarial robustness with additive noise[END_REF][START_REF] Lecuyer | Certified robustness to adversarial examples with differential privacy[END_REF] provides probabilistic guarantees to certify ℓ 2 norm robustness with arbitrarily high confidence. The prediction of a randomized smooth classifier is the most likely prediction returned by the base classifier that is fed by samples from a Gaussian distribution. Salman et al. [START_REF] Salman | Provably robust deep learning via adversarially trained smoothed classifiers[END_REF] further improves the performance of randomized smoothing via adversarial training. All these defenses are computationally expensive and increases the run-time.

. Estimation of Stability

. Quantification of Lipschitz constant

An m-layered feed-forward network Figure 2.1 can be modelled by the following recursive equations:

(∀i ∈ {1, . . . , m}) x i = T i (x i-1 ) = R i (W i x i-1 + b i ) (2.1)
where, at the i th layer, W i ∈ R N i ×N i-1 is the weight matrix, b i ∈ R N i is the bias vector, and R i : R N i → R N i is the activation operator. This operator may consist of the application of basic nonlinear functions, e.g. ReLU or tanh, to each component of the input. Alternatively, it may consist of a softmax operation or group sorting operations which typically arise in max pooling.

As shown in [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF], the problem is related to the choice of the weight matrices. One way of quantifying the stability of the system is to calculate a Lipschitz constant of the network. A Lipschitz constant is an upper bound on the ratio between the variations of the outputs and the variations of inputs of a function T and thus it is a measure of sensitivity of the function with respect to input perturbations. This means that, if θ ∈ [0, +∞[ is a Lipschitz constant of T, then

∥T(x + z) -T(x)∥ ⩽ θ∥z∥ (2.2)
for every input x ∈ R N 0 and perturbation z ∈ R N 0 . Note that the same notation is used here for the norms on R N 0 and R N m , but different norms can be used. If not specified, the standard Euclidean norm will be used. An important assumption is that the operators (R i ) 1⩽i⩽m are non-expansive, i.e. 1-Lipschitz. This assumption is satisfied for all the standard choices of activation operators.

Remark 1 For every (x, z) ∈ (R N 0 ) 2 , if T is continuous on the segment [x, x + z] and differentiable on ]x, x + z[, it follows from the mean value inequality that

∥T(x + z) -T(x)∥ ⩽ ∥z∥ sup α∈]0,1[ ∥T ′ (x + α(z -x))∥ S , (2.3)
where T ′ is the Jacobian of T and ∥ • ∥ S is the spectral norm. We deduce that, when T is differentiable on R N 0 , the optimal Lipschitz constant is θ = sup x∈R N 0 ∥T ′ (x)∥ S . According to the differentiation chain rule, if the activation operators (R i ) 1⩽i⩽m are differentiable, the Jacobian matrix of T at x 0 ∈ R N 0 is expressed as

T ′ (x 0 ) = R ′ m (y m )W m R ′ m-1 (y m-1 ) • • • R ′ 1 (y 1 )W 1 (2.4)
where (∀i ∈ {1, . . . , m}) y i = W i x i-1 + b i .

(2.5)

The variables (x i ) 1⩽i⩽m and (y i ) 1⩽i⩽m-1 are computed during the forward pass. It is important to note that the differentiability assumption is not satisfied by the ReLU activation function or its variants (leaky ReLU, capped ReLU,...).

Remark 2 From a practical viewpoint, it appears sometimes more relevant to determine a local Lipschitz constant of the network. Indeed, the vector of inputs of the network usually varies in a restricted subset C of R N 0 . In addition, it appears reasonable to assume that the perturbations z lie in a ball B(0, ϵ) centered at 0 and with radius ϵ > 0. We can thus define the local Lipschitz constant as

θ C,ϵ = sup x∈C sup z∈B(0,ϵ)\{0} ∥T(x + z) -T(x)∥ ∥z∥ . (2.6)
From this definition, it is clear that θ C,ϵ is lower than the optimal global Lipschitz constant, whose expression is recovered as a limit case when C = R N 0 and ϵ → +∞.

Two point must however be emphasized: i) This local Lipschitz constant may not significantly differ from the global one.

To support this claim, assume that T is differentiable on R N 0 and that there exists x ∈ C such that ∥T ′ ( x)∥ = sup x∈R N 0 ∥T ′ (x)∥. Then it follows from Eq. (2.3) and Eq. (2.6) that θ C,ϵ reduces to the optimal global Lipschitz constant whatever the choice of ϵ > 0.

ii) Computing θ C,ϵ (or an upper bound of it) often is a more difficult problem than the computation of a global Lipschitz constant of T since it corresponds to solving a twice constrained maximization problem. In particular, numerical approaches can be developed where samples within the set C are drawn and adversarial examples are generated by identifying the worst choices for z for each sampled value of the input [START_REF] Gehr | AI2: Safety and robustness certification of neural networks with abstract interpretation[END_REF]. Such approaches may appear better suited for designing networks with improved robustness than for certification purposes. Alternative measures to the local Lipschitz constant have also been considered for multi-class classification problems [START_REF] Yang | Adversarial robustness through local Lipschitzness[END_REF].

. Trivial Upper Bound

The first upper-bound on the Lipschitz constant of a neural network was derived by analyzing the effect of each layer independently and considering a product of the resulting spectral norms [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF]. This leads to the following Upper Bound:

θ m = ∥W m ∥ S ∥W m-1 ∥ S • • • ∥W 1 ∥ S .
(2.7)

Although easy to compute, this upper bound turns out be over-pessimistic and loose.

. Improved Bounds

In [START_REF] Virmaux | Lipschitz regularity of deep neural networks: analysis and efficient estimation[END_REF],the problem of computing the exact Lipschitz constant of a differentiable function is pointed out to be NP-hard. A first generic algorithm (AutoLip) for upper bounding the Lipschitz constant of any automatically differentiable function is proposed. This bound however reduces to Eq. (2.7) for standard feed-forward neural networks. Additionally, the authors proposed, an algorithm SeqLip for sequential neural networks, which shows significant improvement over AutoLip. A sequential neural network is a network for which the activation operators are separable in the sense that, for every i ∈ {1, . . . , m},

(∀x i = (ξ i,k ) 1⩽k⩽N i ∈ R N i ) R i (x) = (ρ i (ξ i,k )) 1⩽k⩽N i , (2.8)
where the activation function ρ i : R → R. 1 In [START_REF] Virmaux | Lipschitz regularity of deep neural networks: analysis and efficient estimation[END_REF], it is assumed that the functions (ρ i ) 1⩽i⩽m are differentiable, increasing, and their derivative are upper bounded by one. By using Eq. (2.4), one can conclude that a Lipschitz constant of the network is

ϑ m = sup Λ 1 ∈D N 1 ([0,1]),...,Λ m-1 ∈D N m-1 ([0,1]) ∥W m Λ m-1 • • • Λ 1 W 1 ∥ S , (2.9)
where D N (I) designate the set of diagonal matrices of dimension N × N with diagonal values in I ⊂ R. This bound simplifies as

ϑ m = sup Λ 1 ∈D N 1 ({0,1}),...,Λ m-1 ∈D N m-1 ({0,1}) ∥W m Λ m-1 • • • Λ 1 W 1 ∥ S , (2.10)
which shows that 2 N i values of the diagonal elements of matrix Λ i have to be tested at each layer i ∈ {1, . . . , m}, so that the global complexity amounts to

2 N 1 +•••+N m-1
and thus grows exponentially as a function of the number of neurons.

Estimating the Lipschitz constant using this method is intractable even for mediumsize networks and authors use a greedy algorithm to compute a bound, which may under-approximate the Lipschitz constant. This does not provide tight upper bounds.

. CPLip

In [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF] various bounds on the Lipschitz constant of a feed-forward network are derived by assuming that, for every i ∈ {1, . . . , m} the activation operator R i is α i -averaged with α i ∈]0, 1]. We recall this means that there exists a non-expansive (i.e. 1-Lipschitz) operator

Q i such that R i = (1 -α i )Id + α i Q i . The following inequality is then satisfied (∀(x, y) ∈ R N i )∥R i x -R i y∥ 2 ⩽ ∥x -y∥ 2 - 1 -α i α i ∥(Id -R i )x -(Id -R i )y∥ 2 .
(2.11)

We thus see that the smaller α i , the more "stable" R i . In the limit case when α 1 = 1, R i is non-expansive and, when α i = 1/2, R i is said to be firmly nonexpansive. An important sub-class of firmly non-expansive operators is the class of proximity operators of convex functions which are proper and lower-semicontinuous. Let Γ 0 (R N ) the be the class of such function defined from R N to ]-∞, +∞]. The proximity operator of a function f ∈ Γ 0 (R N ), at some point x ∈ R N , is the unique vector denoted by prox f (x) such that

prox f (x) = argmin p∈R N 1 2 ∥p -x∥ 2 + f (p).
(2.12)

The proximity operator is a fundamental tool in convex optimization. As shown in [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF], most of the activation functions (sigmoid, ReLU, leaky ReLU, ELU,...) currently used in neural networks are the proximity operators of proper lower-semicontinuous convex functions. This property is also satisfied by activation operators which are not separable like softmax or the squashing function used in capsule networks. The few activation operators which are not proximity operators (e.g. convex combinations of a max pooling and an average pooling) can be viewed as over-relaxations of proximity operators and correspond to a value of the averaging parameter greater than 1/2.

Based on these averaging assumptions, a first estimation of the Lipschitz constant is given by

θ m = β m; p∥W m • • • • • W 1 ∥ + m-1 ∑ k=1 ∑ (j 1 ,...,j k )∈J m,k β m;{j 1 ,...,j k } σ m;{j 1 ,...,j k } , (2.13)
where

(∀ J ⊂ {1, . . . , m -1}) β m;J = ∏ j∈J α j ∏ j∈{1,...,m-1}∖J
(1α j ), (2.14) for every k ∈ {1, . . . , m -1},

J m,k = (j 1 , . . . , j k ) ∈ N k 1 ⩽ j 1 < • • • < j k ⩽ m -1 , if k > 1; {1, . . . , m -1}, if k = 1 (2.15)
and, for every (j 1 , . . . , j k ) ∈ J m,k ,

σ m;{j 1 ,...,j k } = ∥W m • • • W j k +1 ∥ S ∥W j k • • • W j k-1 +1 ∥ S • • • ∥W j 1 • • • W 1 ∥ S . (2.16)
When, for every i ∈ {1, . . . , m -1}, R i is firmly nonexpansive, the expression simplifies as

θ m = 1 2 m-1 ∥W m • • • W 1 ∥ S + m-1 ∑ k=1 ∑ (j 1 ,...,j k )∈J m,k σ m;{j 1 ,...,j k } .
(2.17)

If, for every i ∈ {1, . . . , m -1}, R i is separable,2 a second estimation is provided which reads

ϑ m = sup Λ 1 ∈D N 1 ({2α 1 -1,1}), . . ., Λ m-1 ∈D N m-1 ({2α m-1 -1,1}) ∥W m Λ m-1 • • • Λ 1 W 1 ∥ S , (2.18)
We thus see that, when

α 1 = • • • = α m-1 = 1/2
, we recover Eq. (2.10) without making any assumption on the differentiability of the activation functions. This estimation is more accurate than the previous one in the sense that

∥W m • • • W 1 ∥ S ⩽ ϑ m ⩽ θ m . (2.19) It is proved in [9] that, if the network is with non-negative weights, that is (∀i ∈ {1, . . . , m}) W i ∈ [0, +∞[ N i ×N i , the lower bound in Eq. (2.19) is attained, i.e. ϑ m = ∥W m • • • W 1 ∥ S .
(2.20)

Another interesting result which is established in [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF] is that similar results hold if other norms than the Euclidean norm are used to quantify the perturbations on the input and the output. For example, for a given i ∈ {1, . . . , m}, for every p ∈ [1, +∞], we can define the following norm:

∀x i = (ξ i,k ) 1⩽k⩽N i ∈ R N i ∥x∥ p =          N i ∑ k=1 |ξ i,k | p 1/p , if p < +∞ sup 1⩽k⩽N i |ξ i,k |, if p = +∞.
(2.21)

If (p, q) ∈ [1, +∞] 2
, the input space R N 0 is equipped with the norm ∥ • ∥ p , and the output space R N m is equipped with the norm ∥ • ∥ q a Lipschitz constant for a network with separable activation operators is

ϑ m = sup Λ 1 ∈D N 1 ([2α 1 -1,1]), . . ., Λ m-1 ∈D N m-1 ([2α m-1 -1,1]) ∥W m Λ m-1 • • • Λ 1 W 1 ∥ p,q (2.22)
= sup

Λ 1 ∈D N 1 ({2α 1 -1,1}), . . ., Λ m-1 ∈D N m-1 ({2α m-1 -1,1}) ∥W m Λ m-1 • • • Λ 1 W 1 ∥ p,q , (2.23)
where ∥ • ∥ p,q is the subordinate L p,q matrix norm induced by the two previous norms. The ability to use norms other than the Euclidean one may be sometimes more meaningful in practice (especially for the ℓ 1 or the sup norm). However, computing such a subordinate norm is not always easy [START_REF] Lewis | A top nine list: Most popular induced matrix norms[END_REF].

. SDP based approach

The work in [START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF] focuses on neural networks using separable activation operators. It assumes that the activation function ρ i used at a layer i ∈ {1, . . . , m} is slope-bounded, i.e. there exists non-negative parameters σ min and σ max such that

(∀(ξ, ξ ′ ) ∈ R 2 ) ξ ̸ = ξ ′ ⇒ σ min ⩽ ρ i (ξ) -ρ i (ξ ′ ) ξ -ξ ′ ⩽ σ max .
As noted by the authors, most activation functions satisfy this inequality with σ min = 0 and σ max = 1. In other words, the above inequality means then that ρ i is an increasing function and it is non-expansive. But a known result [41, Proposition 2.4] states that a function ρ i satisfies these properties if and only if it is the proximity operator of some proper lower-semicontinuous convex function.

So it turns out that we recover similar assumptions to those made in [START_REF] Combettes | Deep neural network structures solving variational inequalities[END_REF]. By use of the firm non-expansivity of the operators R i Let us thus assume that σ min = 0, σ max = 1, and m ⩾ 2. It is well-known that R i is firmly nonexpansive if and only if

(∀(x, y) ∈ (R N i ) 2 (x -y) ⊤ (R i x -R i y) ⩾ ∥R i x -R i y∥ 2 .
(2.24)

The point is that, if R i is a separable operator, this inequality holds in a more general metric associated with a matrix

Q i = Diag(q i,1,1 , . . . , q i,N i ,N i ), (2.25)
where (∀k ∈ {1, . . . , N i } 2 ) q i,k,k ⩾ 0. In the following, the set of such matrices (Q i ) 1⩽i⩽m-1 will be denoted by Q. This means that

(∀(x, y) ∈ (R N i ) 2 (x -y) ⊤ Q i (R i x -R i y) ⩾ (R i x -R i y) ⊤ Q i (R i x -R i y).
(2.26)

For every (x i , y i ) ∈ (R N i ) 2 , let x i = R i (W i x i-1 + b i ) and y i = R i (W i y i-1 + b i ). It follows from Eq. (2.26) that W i (x i-1 -y i-1 ) ⊤ Q i (x i -y i ) ⩾ (x i -y i ) ⊤ Q i (x i -y i ).
(2.27)

Summing for the first m -1 layers yields

m-1 ∑ i=1 W i (x i-1 -y i-1 ) ⊤ Q i (x i -y i ) ⩾ m-1 ∑ i=1 (x i -y i ) ⊤ Q i (x i -y i ).
(2.28)

On the other hand, ϑ m > 0 is a Lipschitz constant of the neural network T if ϑ 2 m ∥x 0 -y 0 ∥ 2 ⩾ ∥W m (x m-1 -y m-1 )∥ 2 .
(2.29)

For the latter inequality to hold, it is thus sufficient to ensure that,

ϑ 2 m ∥x 0 -y 0 ∥ 2 -∥W m (x m-1 -y m-1 )∥ 2 ⩾ 2 m-1 ∑ i=1 W i (x i-1 -y i-1 ) ⊤ Q i (x i -y i ) -2 m-1 ∑ i=1 (x i -y i ) ⊤ Q i (x i -y i ). (2.

30)

This inequality can be rewritten in matrix form as

   x 0 -y 0 . . . x m-1 -y m-1    ⊤ M(ρ m , Q 1 , . . . , Q m-1 )    x 0 -y 0 . . . x m-1 -y m-1    ⩾ 0 (2.31) with ρ m = ϑ 2 m and M(ρ m , Q 1 , . . . , Q m-1 ) =          ρ m Id N 0 -W ⊤ 1 Q 1 0 -Q 1 W 1 0 . . . . . . . . . . . . . . . 0 -W ⊤ m-1 Q m-1 0 -Q m-1 W m-1 2Q m-1 -W ⊤ m W m          . (2.32)
In the case of a network having just one hidden layer, which is mainly investigated in [START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF], the above matrix reduces to

M(ρ 2 , Q 1 ) = ρ 2 Id N 0 -W ⊤ 1 Q 1 -Q 1 W 1 2Q 1 -W ⊤ 2 W 2 .
(2.33)

Condition (2.31
) is satisfied, for every (x 0 , . . . , x m-1 ) and (y 0 , . . . , y m-1 ) if and only if

M(ρ m , Q 1 , . . . , Q m-1 ) ⪰ 0.
(2.34)

It is actually sufficient that this positive semidefiniteness constraint be satisfied for any matrix Q ∈ Q N i for √ ρ m to be a Lipschitz constant. The smallest possible value of the resulting constant can be obtained by solving the following semidefinite positive programming (SDP) problem:

minimize (ρ m ,Q 1 ,...,Q m-1 )∈C ρ m , (2.35)
where C is the closed convex set

C = (ρ m , Q 1 , . . . , Q m-1 ) ∈ [0, +∞[×Q Eq. (2.34) holds . (2.36)
Although there exists efficient SDP solvers, the method remains computationally intensive. A solution to reduce its computational complexity at the expense of a lower accuracy consists of restricting the optimization of the metric matrices Q 1 , . . . , Q m-1 to a subset of Q.

One limitation of this method is that it is taylored to the use of the Euclidean norm.

Remark 3

In [START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF], it is claimed that Eq. (2.26) is valid for every metric matrix

Q i = N i ∑ k=1 q i,k,k e k e ⊤ k + ∑ 1⩽k<ℓ⩽N i q i,k,l (e k -e ℓ )(e k -e ℓ ) ⊤ , (2.37)
where (e k ) 1⩽k⩽N i is the canonical basis of R N i and (∀(k, ℓ) ∈ {1, . . . , N i } 2 ) with k ⩽ ℓ, q i,k,ℓ ⩾ 0 Unfortunately, this turns out to be incorrect as shown by the counterexample next. Assume that all the coefficients (q i,k,ℓ ) 1⩽k⩽ℓ⩽N i are zero, except one in the second summation set to 1. We have then

Q = (e k -e ℓ )(e k -e ℓ ) ⊤ (2.38) for some (k, ℓ) ∈ {1, . . . , N i } 2 with k < ℓ. Consider the simple case when ρ i is the ReLU function, that is (∀ξ ∈ R) ρ i (ξ) = max{ξ, 0}.
ρ i is an example of slope-bounded function with σ min = 0 and σ max = 1. Then, for every

x = (ξ k ) 1⩽k⩽N i ∈ R N i and y = (υ k ) 1⩽k⩽N i ∈ R N i , (e k -e ℓ ) ⊤ (x -y) = ξ k -ξ ℓ -υ k + υ ℓ (e k -e ℓ ) ⊤ (R i (x) -R i (y)) = ρ i (ξ k ) -ρ i (ξ ℓ ) -ρ i (υ k ) + ρ i (υ ℓ ).
(2.39)

We have thus

(x -y) ⊤ Q(R i x -R i y) = (ξ k -ξ ℓ -υ k + υ ℓ )(ρ i (ξ k ) -ρ i (ξ ℓ ) -ρ i (υ k ) + ρ i (υ ℓ )) (R i x -R i y) ⊤ Q(R i x -R i y) = (ρ i (ξ k ) -ρ i (ξ ℓ ) -ρ i (υ k ) + ρ i (υ ℓ )) 2 .
(2.40)

If we now assume that ξ k < 0, υ k = 0, and υ ℓ > ξ ℓ > 0, we deduce from the expression of the ReLU function that

(x -y) ⊤ Q(R i x -R i y) = (ξ k + υ ℓ -ξ ℓ )(υ ℓ -ξ ℓ ) (R i x -R i y) ⊤ Q(R i x -R i y) = (υ ℓ -ξ ℓ ) 2 .
(2.41) Inequality Eq. (2.26) then becomes

ξ k + υ ℓ -ξ ℓ ⩾ υ ℓ -ξ ℓ , (2.42)
which contradicts the fact that ξ k has been chosen negative.

The erroneous statement comes from a flaw in the deduction of Lemma 1 from Lemma 2 in [START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF]. Another counterexample was also recently provided in [START_REF] Pauli | Training robust neural networks using Lipschitz bounds[END_REF].

. Polynomial optimization based approach

The approach in [START_REF] Latorre | Lipschitz constant estimation of neural networks via sparse polynomial optimization[END_REF] applies to neural networks having a single output (i.e. N m = 1)3 . The authors state that their approach is restricted to differentiable activation functions, but it is actually valid for any separable firmly non-expansive activation operators. Indeed, when N m = 1, the Lipschitz constant in Eq. (2.22) reduces to

ϑ m = sup Λ 1 ∈D N 1 ([0,1]), . . ., Λ m-1 ∈D N m-1 ([0,1]) ∥W ⊤ 1 Λ 1 • • • Λ m-1 W ⊤ m ∥ p * , (2.43)
where p * ∈ [1, +∞] is the dual exponent of p (such that 1/p + 1/p * = 1). Recall that p ∈ [1, +∞] is the exponent of the ℓ p -norm equipping the input space. This shows that ϑ m is equal to

ϑ m = sup{Φ(x, λ 1 , . . . , λ m-1 ) | ∥x∥ p ⩽ 1, (λ i ) 1⩽i⩽m-1 ∈ [0, 1] N 1 +•••+N m-1 }, (2.44)
where, for every

x ∈ R N 0 and (λ i ) 1⩽i⩽m ∈ R N 1 +•••+N m-1 , Φ(x, λ 1 , . . . , λ m-1 ) = x ⊤ W ⊤ 1 Diag(λ 1 ) • • • Diag(λ m-1 )W ⊤ m .
(2.45)

Function Φ is a multivariate polynomial of the components of its vector arguments. Therefore if the unit ball associated with the ℓ p norm can be described via polynomial inequalities, which happens when p ∈ N \ {0} and p = +∞, then finding ϑ m turns out to be a polynomial constrained optimization problem. Solving such an optimization problem can be achieved by solving a hierarchy of convex problems. However, the size of the hierarchy tends to grow fast and if the order of the hierarchy is truncated to a too small value, the delivered result becomes inaccurate.

Leveraging the sparsity properties that might exist for the weight matrices may be helpful numerically. Note that the approach is further improved in [START_REF] Chen | Polynomial optimization for bounding Lipschitz constants of deep networks[END_REF] by using Lasserre's hierarchy.

Remark 4

i) Any polynomial optimization problem can be reduced to a quadratically constrained quadratic program. This approach is followed in [START_REF] Raghunathan | Semidefinite relaxations for certifying robustness to adversarial examples[END_REF] and the problem is then solved via an SDP relaxation. This leads to an upper bound on ϑ m .

ii) In [START_REF] Latorre | Lipschitz constant estimation of neural networks via sparse polynomial optimization[END_REF], the authors argue that their method allows a local Lipschitz constant to be computed. Indeed, in the case when the activation functions are differentiable, Eq. (2.4) holds where

(∀i ∈ {1, . . . , m -1}) R ′ i (y i ) = Diag(ρ ′ i (υ i,1 ), . . . , ρ ′ i (υ i,N i )) (2.46)
and y i = (υ i,k ) 1⩽k⩽N i is defined by Eq. (2.5). Consequently, if the input x of the network belongs to a bounded convex set C, then, for every i ∈ {1 . . . , m -1} y i varies in a bounded set too, which implies that (ρ ′ i (υ i,k )) 1⩽k⩽N i belongs to a polytope C i (say a hyperrectangle) of [0, 1] N i . This allows us to define a local Lipschitz constant as

ϑ m,C = sup{Φ(x, λ 1 , . . . , λ m-1 ) | ∥x∥ p ⩽ 1, (λ i ) 1⩽i⩽m-1 ∈ C 1 × • • • × C m-1 }.
(2.47) Computing this constant is still a a polynomial optimization problem. So, the same approach as in the global case can be followed, with an extra effort for the determination of polytopes (C i ) 1⩽i⩽m-1 .

. Local Lipschitz Bound

As discussed briefly in Remark 2 local definitions of Lipschitz bounds are also possible. Few works in literature utilize local Lipschitz bounds to improve the robustness of the NNs. The advantage of using a local bound is that we may theoretically expect tighter bounds since of course the local Lipschitz constant cannot be any larger than the global Lipschitz constant. However, using a local bound also has its cons. Firstly, a local bound is more computationally expensive, as each instance has its own bound, hence the required memory grows with the batch size. This in turn reduces the amount of parallelism that can be exploited when using a local bound, reducing the model throughput. Furthermore, because the local Lipschitz constant is different for every point, it must be computed every time the network sees a new point. By contrast, the global bound can be computed in advance, meaning that verification via the global bound is essentially free. This makes the global bound advantageous, assuming that it can be effectively leveraged for verification. Finding an exact local Lipschitz constant for a neural network is generally NP-hard [START_REF] Virmaux | Lipschitz regularity of deep neural networks: analysis and efficient estimation[END_REF], so most works focus on finding a sound upper bound. Recent work has explored methods for obtaining upper bounds on the local Lipschitz constant [START_REF] Weng | Towards fast computation of certified robustness for relu networks[END_REF][START_REF] Zhang | Efficient neural network robustness certification with general activation functions[END_REF][START_REF] Lee | Lipschitz-certifiable training with a tight outer bound[END_REF]. In [START_REF] Weng | Evaluating the robustness of neural networks: An extreme value theory approach[END_REF], the authors convert the robustness analysis problem into a local Lipschitz constant estimation problem, where they estimate this local constant by a set of independently and identically sampled local gradients. This algorithm is scalable but is not guaranteed to provide upper bounds. In a similar work, the authors of [START_REF] Weng | Towards fast computation of certified robustness for relu networks[END_REF] exploit the piece-wise linear structure of ReLU activation functions to estimate the local Lipschitz constant of neural networks. In [START_REF] Fazlyab | Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming[END_REF], the authors use quadratic constraints and semidefinite programming to analyze local (point-wise) robustness of neural networks. Hein et al. [START_REF] Hein | Formal guarantees on the robustness of a classifier against adversarial manipulation[END_REF] derived an analytical bound for two layer neural networks and found that local Lipschitz bounds could be much tighter than the global one and provide better robustness certificates. RecurJac [START_REF] Zhang | Recurjac: An efficient recursive algorithm for bounding jacobian matrix of neural networks and its applications[END_REF] is a recursive algorithm that analyzes the local Lipschitz constant in a neural network using a bound propagation [START_REF] Zhang | Towards stable and efficient training of verifiably robust neural networks[END_REF] based approach. FastLip [START_REF] Weng | Towards fast computation of certified robustness for relu networks[END_REF] is a special and weaker form of RecurJac. Jordan et al. [START_REF] Jordan | Exactly computing the local Lipschitz constant of relu networks[END_REF] formulated the computation of the local Lipschitz constant as a Mixed Integer Linear Programming (MILP) problem and they were able to solve it. Although the authors of these approaches claim that they can obtain reasonably tight and sound local Lipschitz constants, none of them have been demonstrated effective for training a certifiably robust network, which requires high efficiency and scalability.

. Dealing with scalability issues

To handle scalability issues arising when estimating a Lipschitz constant of a deep or broad network, several techniques can be employed as we show next.

. Serial splitting of the network

A feed-forward network can be decomposed as a cascade of smaller sub-networks (U j ) 1⩽j⩽k , where the j-th sub-network comprises m j weight operators (hence, ∑ k j=1 m j = m). Let θ U j ,m j be a Lipschitz constant of the j-th sub-network. An overall Lipschitz constant for T is given by

θ T,m = k ∏ j=1 θ U j ,m j .
(2.48) Various decompositions can be performed which lead to more or less tight bounds. For example, if we constrain the size of the sub-networks to be equal to m max ⩽ m, except possibly the first and the last one, the number of possible decompositions is equal to m max while the maximum number of sub-networks in a decomposition is equal to ⌊m/m max ⌋ + 1. If m = 7 and m max = 3, we can choose

     U 1 = T 1 U 2 = T 4 • T 3 • T 2 U 3 = T 7 • T 6 • T 5 , (2.49)      U 1 = T 2 • T 1 U 2 = T 5 • T 4 • T 3 U 3 = T 7 • T 6 , (2.50) or      U 1 = T 3 • T 2 • T 1 U 2 = T 6 • T 5 • T 4 U 3 = T 7 .
(2.51)

If we look at the complexity reduction based on the computation of Eq. (2.23), we see that for the j-th sub-network, the computational complexity is of the order of 2 N U j where N U j is the number of neurons in its hidden layers. For the previous example splittings, the complexity are thus reduced to 2

N 2 +N 3 + 2 N 5 +N 6 , 2 N 1 + 2 N 3 +N 4 + 2 N 6 , and 2 N 1 +N 2 + 2 N 4 +N 5
, respectively (to be compared with An additional degree of flexibility consists of factorizing the weight matrix W i at a given layer i as

2 N 1 +•••+N 6 ).
W i = V i V i (2.52)
where

V i ∈ R N i ×M i and V i ∈ R M i ×N i-1
. Such a factorization may result from a singular value decomposition of W i . In [START_REF] Virmaux | Lipschitz regularity of deep neural networks: analysis and efficient estimation[END_REF], it is proposed to perform such a decomposition and to consider k = m -1 sub-networks having one hidden layer, of the form

(∀i ∈ {1, . . . , m -2}) U i = V i+1 • R i • ( V i • +b i ) (2.53) U m-1 = R m • (V m • +b m ) • R m-1 • ( V m-1 • +b m-1 ) (2.54) with V 1 = W 1 and V m = W m .
This yields a Lipschitz constant of the form Eq. (2.48) where

(∀i ∈ {1, . . . , m -1}) θ U i ,N i = sup Λ i ∈D N i ∥V i+1 Λ i V i ∥ p,q , (2.55)
thus inducing a global computational complexity of the order of

2 N 1 + • • • + 2 N m-1 .

. Parallel splitting of the network

When the number of neurons at a given layer is large, a parallel splitting of the network may be preferable. This technique mainly applies to the computation of Eq. (2.23). The basic idea consists of block-decomposing the matrix Λ i ∈ D N i associated with the activation operator used at layer i ∈ {1, . . . , m -1} in J i diagonal matrices Λ i,j ∈ D N i,j with j ∈ {1, . . . , J i }. For every i ∈ {1, . . . , m -1}, we have thus

Λ i =       Λ i,1 0 . . . 0 0 Λ i,2 . . . . . . . . . . . . . . . 0 0 . . . 0 Λ i,J i       (2.56)
where the j-th block contains a reduced number N i,j of diagonal elements. Then

∑ J i j=1 N i,j = N i . Set J m = 1, J 0 = 1, N m,1 = N m , and N 0,1 = N 0 .
For every i ∈ {1, . . . , m}, each weight matrix can be similarly block-decomposed as

W i =    W i,1,1 . . . W i,1,J i-1 . . . . . . W i,J i ,1 . . . W i,J i ,J i-1    (2.57)
where, for every j ′ ∈ {1, . . . , J i } and j ∈ {1, . . . ,

J i-1 }, W i,j,j ′ is a matrix of dimension N i,j ′ × N i-1,j . We have then W m Λ m-1 W m-1 • • • W 2 Λ 1 W 1 = J 1 ∑ j 1 =1 • • • J m-1 ∑ j m-1 =1 W m,1,j m-1 Λ m-1,j m-1 W m-1,j m-1 ,j m-2 • • • W 2,j 2 ,j 1 Λ 1,j 1 W 1,j 1 ,1 (2.58)
By using now the triangle inequality, we deduce that

∥W m Λ m-1 W m-1 • • • W 2 Λ 1 W 1 ∥ p,q ⩽ J 1 ∑ j 1 =1 • • • J m-1 ∑ j m-1 =1 ∥W m,1,j m-1 Λ m-1,j m-1 W m-1,j m-1 ,j m-2 • • • W 2,j 2 ,j 1 Λ 1,j 1 W 1,j 1 ,1 ∥ p,q .
This leads to

ϑ m ⩽ J 1 ∑ j 1 =1 • • • J m-1 ∑ j m-1 =1 ϑ m,j 1 ,...,j m-1 , (2.59)
where

ϑ m,j 1 ,...,j m-1 = sup Λ 1,j 1 ∈D N 1,j 1 ({2α 1 -1,1}), . . ., Λ m-1,j m-1 ∈D N m-1,j m-1 ({2α m-1 -1,1}) ∥W m,1,j m-1 Λ m-1,j m-1 • • • Λ 1,j 1 W 1,j 1 ,1 ∥ p,q . (2.

60)

The complexity of the computation of the latter constant is of the order of

2 N 1,j 1 +•••+N m-1,j m-1 .
The complexity of the Lipschitz constant defined by the upperbound in Eq. (2.59) thus amounts to

J 1 ∑ j 1 =1 • • • J m-1 ∑ j m-1 =1 2 N 1,j 1 +•••+N m-1,j m-1
For example, if for every layer i ∈ {1, . . . , m -1}, a decomposition in 2 blocks of equal size is performed (i.e. J i = 2 and N i,1

= N i,2 = N i /2 ∈ N), the resulting complexity is 2 (N 1 +...+N m-1 )/2+m-1 .
As an illustration, when m = 2, we get

W 1 = W 1,1,1 W 1,2,1
(2.61)

W 2 = [W 2,1,1 W 2,1,2 ].
(2.62)

By omitting redundant indices these decompositions can be rewritten as

W 1 = W 1,1 W 1,2
(2.63)

W 2 = [W 2,1 W 2,2 ]
(2.64)

and Eq. (2.59) becomes

ϑ m ⩽ sup Λ 1,1 ∈D N 1,1 ({2α 1 -1,1}) ∥W 2,1 Λ 1,1 W 1,1 ∥ p,q + sup Λ 1,2 ∈D N 1,2 ({2α 1 -1,1}) ∥W 2,2 Λ 1,2 W 1,2 ∥ p,q .
(2.65)

. Serial and parallel splitting of the network

A combination of the two previous splitting techniques can be used for an efficient estimation of the Lipschitz constant of a feed-forward neural network. [START_REF] Szegedy | Intriguing properties of neural networks[END_REF] proposed L-BFGS (Limited-memory; Broyden, Fletcher, Goldforb, Shanno) to construct adversarial attacks and since then, there has been a plethora of works introducing various adversarial attacks and their defenses for DNNs. Goodfellow et al. [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] proposed a simpler and faster method to construct adversarial attacks (FGSM-Fast Gradient Sign Method). The generated images are misclassified by adding perturbations and linearizing the cost function in the gradient direction. This is a non-iterative attack; hence it has a lower computation cost than the previous method. The Fast Gradient Sign Method (FGSM) is an ℓ ∞ bounded attack and is often prone to label leaking.

. Control of Stability

It may be difficult for FGSM to control the perturbation level in constructing attacks. Kurakin et al. [START_REF] Kurakin | Adversarial examples in the physical world[END_REF] proposed an optimized FGSM, termed Iterative Gradient Sign Method (IGSM), which adds perturbations in multiple smaller steps and clips the results after each iteration ensuring that the perturbations are restricted to the neighborhood of the example. [START_REF] Dong | Boosting adversarial attacks with momentum[END_REF] added momentum to IGSM attacks. [START_REF] Papernot | The limitations of deep learning in adversarial settings[END_REF] proposed the Jacobian-based Saliency Map Attack (JSMA), based on the ℓ 0 sparsity measure. The basic idea is to construct a saliency map with the gradients and model the gradients based on the impact of each image pixel.

Moosavi et al. [START_REF] Moosavi-Dezfooli | Deepfool: a simple and accurate method to fool deep neural networks[END_REF] proposed a non-targeted attack method based on the ℓ 2norm, called DeepFool. It tries to find the decision boundary closest to the sample in the image space and then uses the classification boundary to fool the classifier. FGSM, JSMA, and DeepFool are designed to generate adversarial attacks corresponding to a single image to fool the trained classifier model. Moosavi et al. [START_REF] Moosavi-Dezfooli | Universal adversarial perturbations[END_REF] proposed a universal image-agnostic perturbation attack method that deceives the classifier by adding a single perturbation to all images in the dataset. Carlini et al. [START_REF] Carlini | Towards evaluating the robustness of neural networks[END_REF] proposed a powerful attack based on L-BFGS. The attack can be generated according to ℓ 1 , ℓ 2 , and ℓ ∞ norm, which can be targeted or non-targeted. Liu et al. [START_REF] Liu | Delving into transferable adversarial examples and black-box attacks[END_REF] proposed an ensemble attack method combining multiple models to construct adversarial attacks. Rony et al. [START_REF] Rony | Augmented Lagrangian adversarial attacks[END_REF] proposed a method to generate minimally perturbed adversarial examples based on Augmented Lagrangian for various distance metrics. In [START_REF] Balda | Perturbation analysis of learning algorithms: A unifying perspective on generation of adversarial examples[END_REF], authors propose a general framework for generating adversarial examples in both classification and regression tasks for applications in the image domain. Most of the methods in the literature about adversarial example generation belong to the class of white box attackers, i.e., the attacker has access to the information related to the trained neural network model, including the model architecture and its parameters. A black box attacker is introduced in [START_REF] Su | One pixel attack for fooling deep neural networks[END_REF]. Such attackers do not know the model but can interact with it. [START_REF] Su | One pixel attack for fooling deep neural networks[END_REF] proposed a one-pixel attack method that only changes one pixel for each image to construct an adversarial attack to fool DNNs. A byproduct of black-box attack is grey-box attack, where attackers might have limited information regarding the model. Weng et al. [START_REF] Weng | Evaluating the robustness of neural networks: An extreme value theory approach[END_REF] proposed a computationally feasible method called Cross Lipschitz Extreme Value for nEtwork Robustness (CLEVER), which applies extreme value theory to estimate a lower bound of the minimum adversarial perturbation required to misclassify the image. CLEVER is the first attack-independent method and can also evaluate the intrinsic robustness of neural networks.

Most of the introduced attack methods are applied digitally, where the adversary supplies input images directly to the DNN. However, this is not always the case in scenarios that use cameras, microphones, or other sensors to receive the signals as input. These systems can still be attacked by generating physical-world adversarial objects. For example, Eykholt et al. [START_REF] Eykholt | Robust physical-world attacks on deep learning visual classification[END_REF] attached stickers to road signs that can severely threaten autonomous car sign recognizer and in [START_REF] Athalye | Synthesizing robust adversarial examples[END_REF] introduced the generation of physical 3D adversarial objects. These kinds of adversarial objects are more destructive for deep learning models because they can directly challenge many practical applications of DNN, such as face recognition, autonomous vehicles, etc.

Regression: In regression tasks, there are no natural margins as in the case of classification tasks, and difficulties hinder adversarial learning in a regression setting to define the adversarial attacks, their success, and evaluation metrics. Despite the number of works in adversarial attack generation, there are few articles dealing with regression tasks. Tong et al. [START_REF] Tong | Adversarial regression with multiple learners[END_REF] looked at adversarial attacks in the setting of an ensemble of multiple learners, investigating the interactions between these linear learners and an attacker in regression tasks, modeled as a Multi-Learner Stackelberg Game (MLSG). However, the investigated linear case cannot capture the larger class of non-linear models. The focus only on specific applications of regression is common. Ghafouri et al. [START_REF] Ghafouri | Adversarial regression for detecting attacks in cyber-physical systems[END_REF] examined an important problem: selecting an optimal threshold for each sensor against an adversary for regression tasks in cyberphysical systems. Deng et al. [START_REF] Deng | An analysis of adversarial attacks and defenses on autonomous driving models[END_REF] introduced the concept of adversarial threshold, which is related to a deviation between the original prediction and the prediction of an adversarial example, i.e., an acceptable error range in driving models. In a regression context, [START_REF] Nguyen | Adversarial attacks, regression, and numerical stability regularization[END_REF] introduced a generically useful defense to reduce the effectiveness of adversarial attacks. They consider adversarial attacks as a potential symptoms of numerical instability in the learned function.

. Adversarial Defenses

Adversarial Training : All ML models are vulnerable to adversarial attacks [START_REF] Papernot | Practical black-box attacks against machine learning[END_REF]. Therefore, defending against adversarial examples is urgent for ML security. Goodfellow et al. [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] proposed the concept of adversarial training to improve the robustness of the model. The idea is to add adversarial attacks to the training data and continuously generate new examples at each training step. The number and relative weight of adversarial examples at each batch is controlled by the loss function independently. Adversarial training can be viewed as the process of minimizing classification error rates when the data is maliciously perturbed and it can be formulated as a minimax problem. Although a model can be robust to whitebox attacks after the adversarial training step, it is still vulnerable to the adversarial attacks generated from other models, i.e., the model is not robust to black-box attacks. Based on this, Tramer et al. [START_REF] Tramèr | Ensemble adversarial training: Attacks and defenses[END_REF] proposed the concept of ensemble adversarial training. It consists in augmenting the training data constructed not only from the model being trained but also from the other pre-trained models, increasing the diversity of adversarial attacks and improving the generalization ability.

Defensive Strategies to improve robustness, such as adversarial training only provide empirical robustness, without any formal guarantees on the safety of trained models. Stronger and more aggressive attacks have successfully broken many existing adversarial defenses [START_REF] Athalye | Obfuscated gradients give a false sense of security: Circumventing defenses to adversarial examples[END_REF]. Similarly, Carlini and Wagner [START_REF] Carlini | Towards evaluating the robustness of neural networks[END_REF] show that the method of Defensive Distillation is still vulnerable to their adversarial examples. In contrast, certified defenses using bounds such as the Lipschitz constant provide formal robustness guarantees that any norm-bounded adversary cannot alter the prediction of a given network.

. Control using Lipschitz constant

Utilizing the Lipschitz constant to certify robustness has been studied at several instances in the literature. Neural networks trained without any robustness constraint usually have very large global Lipschitz constant bounds [START_REF] Szegedy | Intriguing properties of neural networks[END_REF], so most existing works train the network with a criterion that promotes small Lipschitz bound.

Parseval Networks: [START_REF] Cisse | Parseval networks: Improving robustness to adversarial examples[END_REF] designed networks with orthogonal weights, whose Lipschitz constants are exactly 1. This can be too restrictive and later works mostly use the power iteration to obtain per-layer induced norms, whose product is a Lipschitz constant Spectral Normalization: [START_REF] Miyato | Spectral normalization for generative adversarial networks[END_REF] showed control on the Lipschitz constant using spectral normalization for GANs. We discuss Spectral Normalization in Section 2.3.2.1.

Hinge Regularisation: The authors of [START_REF] Serrurier | Achieving robustness in classification using optimal transport with hinge regularization[END_REF] propose to learn 1-Lipschitz networks for binary classification using a new loss which is a hinge regularized version of the Kantorovich-Rubinstein dual formulation for the Wasserstein distance estimation. They prove that the proposed loss function has a direct interpretation in terms of adversarial robustness together with certifiable robustness bound.

Prior work seeks to use global or local Lipschitz bounds during training to promote robustness. Computing and imposing loose global Lipschitz bounds is often easy but over-regularize the network during training and decrease its accuracy on clean data.

Lipschitz Margin training(LMT): In [START_REF] Tsuzuku | Lipschitz-margin training: Scalable certification of perturbation invariance for deep neural networks[END_REF] authors train models that are certifiably robust by constructing a new loss on worst logits using the global Lipschitz bounds. They add a value of √ 2ϵL glob , where ϵ is the perturbation radius, to all logits other than the logit corresponding to the ground-truth class and calculate the loss thereafter.

Box constrained propagation (BCP): Lee et al. [START_REF] Lee | Lipschitz-certifiable training with a tight outer bound[END_REF] achieves a tighter outer bound than global Lipschitz-based outer bound, by considering local information via interval bound (box) propagation. They also compute the worst-case logit based on the intersection of a (global) ball and a (local) box.

GloRo: Leino et al. [START_REF] Leino | Globally-robust neural networks[END_REF] bounds the upper bounds on the worst margins using the global Lipschitz constant. It constructs a new logit with a newly constructed class called the bottom class and determines if the sample can be certified.

Local-Lip: [START_REF] Huang | Training certifiably robust neural networks with efficient local Lipschitz bounds[END_REF] utilizes the interactions between activation functions (e.g. ReLU, MaxMin) and weight matrices. By eliminating the corresponding rows and columns where the activation function output is constant, they guarantee a lower provable local Lipschitz bound than the global Lipschitz bound for the neural network.

LipNet1: Bethune et al. [START_REF] Béthune | The many faces of 1-Lipschitz neural networks[END_REF] concentrate on classifiers constructed with 1-Lipschitz networks. They use DeelLip Library [START_REF] Serrurier | Achieving robustness in classification using optimal transport with hinge regularization[END_REF] for constructing 1-Lipschitz neural networks. They show that the classifiers trained with 1-Lipschitz constrained can be as expressive as those trained without any Lipschitz constraints and are certifiably robust. They introduce an implicit parameter τ, which controls the network expressiveness. A careful examination of the algorithm reveals that controlling this parameter is equivalent to relaxing the Lipschitz constant of the neural network, which is the approach we follow in a more direct manner in this work.

Known Lipschitz Constants: In Piat et al. [START_REF] Piat | Towards an evaluation of Lipschitz constant estimation algorithms by building models with a known Lipschitz constant[END_REF], authors learn general neural networks for which the target Lipschitz constants are known before hand. They propose a PGD-like algorithm to approximate the target fixed Lipschitz constant by maximizing a spectral norm.

. Spectral Normalization

Spectral normalization stabilizes DNN training by constraining the Lipschitz constant of the objective function. Spectrally normalized DNNs have also been shown to generalize well [START_REF] Bartlett | Spectrally-normalized margin bounds for neural networks[END_REF], which is an indication of stability in machine learning. For fully-connected layers, the spectral norm of a given weight matrix W can be computed by using the approach described by [START_REF] Miyato | Spectral normalization for generative adversarial networks[END_REF] using the power iteration method. For each W, we initialize a vector ũ and approximate both the left and right singular vectors by iterating the update rules:

ṽ ← W ũ/∥W ũ∥ 2 ũ ← W ⊤ ṽ/∥W ⊤ ṽ∥ 2
The final singular value can be approximated with σ(W) ≈ ṽT W ũ. The final weight matrix after spectral normalization W SN is given as

W SN = W/σ(W).
In the case of a CNN consisting of convolutional layers, dense layers, and ReLU activations, the spectral norm of each of the ReLU layers can be calculated using the power iteration method [START_REF] Farnia | Generalizable adversarial training via spectral normalization[END_REF][START_REF] Gouk | Regularisation of neural networks by enforcing Lipschitz continuity[END_REF]. Gouk et al. [START_REF] Gouk | Regularisation of neural networks by enforcing Lipschitz continuity[END_REF] also give a procedure for bounding the spectral norm of skip connections, and batch normalization layers enable this approach on ResNet architectures. The effectiveness of spectral normalization in training safety-critical systems was shown in [START_REF] Shi | Neural lander: Stable drone landing control using learned dynamics[END_REF][START_REF] O'connell | Neural-fly enables rapid learning for agile flight in strong winds[END_REF]. The DNNs are trained with layer-wise spectrally normalized weight matrices. In practice, we can apply spectral normalization to the weight matrices in each layer during training as follows:

W SN = W σ(W) γ 1/m (2.66)
where γ is the intended global Lipschitz constant of the DNN having m layers.

. Conclusion

Several developments related to various state-of-art approaches have been presented in this chapter. In the following chapters, we present some contributions and discuss the limitations of state-of-the-art techniques. We show how to tackle a few of these shortcomings so as to better understand and improve the stability of neural networks. We evaluate our proposed methods on diverse open-source and industrial datasets in the field of safety-critical applications.

Chapter 3 Adversarial Attacks on Regression Tasks

. Introduction

Adversarial attacks against neural networks and their defenses have been primarily investigated in classification tasks. However, adversarial attacks in a regression setting remain understudied, although regression tasks are an important part of mission-critical and safety-critical use cases in many domains like aviation, space, transportation, and defense. In this work, we present an adversarial attacker for regression tasks derived from the algebraic properties of the Jacobian of the neural network. We show that our attacker successfully fools the neural network and we measure its effectiveness in reducing the estimation performance. This white-box adversarial attacker aims to support engineers in designing safety-critical regression machine learning models. We present our results on various open-source and real industrial tabular datasets. In particular, the proposed adversarial attacker outperforms attackers based on random perturbations of the inputs. Our analysis relies on the quantification of the fooling error as well as various error metrics. A noteworthy feature of our attacker is that it allows us to optimally attack a subset of inputs, which may help to analyze the sensitivity of some specific inputs, as we will see in the next chapter.

. Limitations of Previous Works

As pointed out in Section 2.3.1, adversarial machine learning has received increased attention in the past decade. Adversarial attacks cause vulnerability in model deployment and especially need to be considered in the deployment of security-critical AI applications. Despite the newfound interest of the research community in trustworthy and explainable AI, only a few works are investigating adversaries in the case of regression tasks.

Current advances in the adversarial machine learning field evolve around the issue of designing attacks and defenses with a focus on the use of neural networks in image analysis and computer vision [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF], [START_REF] Kurakin | Adversarial examples in the physical world[END_REF]. Much fewer works concern tabular data. However, most machine learning tasks in the industry rely on tabular data, e.g., fraud detection, product failure prediction, anti-money laundering, recommendation systems, click-through rate prediction, or flight arrival time prediction. To the best of our knowledge, the only work dealing with adversarial attacks in white box settings for tabular data has been proposed in [START_REF] Ballet | Imperceptible adversarial attacks on tabular data[END_REF], and this work handles only classification tasks. In this work, we focus on generating adversarial attacks for neural networks in the specific scenario when i) a regression task is performed and ii) tabular data are employed.

. Proposed Attacker

. Objective

As already mentioned, the problem of adversarial attacks is closely related to the robustness issue of a neural network, i.e., its sensitivity to perturbations. Let T : R N 0 → R N m be the considered neural network having N 0 scalar inputs and N m scalar outputs. If x ∈ R N 0 is a given vector of inputs for some data for which y is the associated target output, the network has been trained to produce an output T(x) close to y. If the input is now perturbed by an additive vector e ∈ R N 0 , the perturbed output is T(x + e). Attacking the network then amounts to finding a perturbation e of preset magnitude, which makes the output of the network maximally deviate from a reference output. This reference output may be the model output T(x) or the ground truth output y. Since our purpose is to develop an efficient approach even if the network accuracy is not very high, we choose y as the reference output when available. In this context, measures of deviation and magnitude of the perturbation play an important role in the mathematical formulation of the problem. As a common choice, the measure of perturbation magnitude will be here an ℓ p -norm where p ∈ [1, +∞]. For measuring the output deviation, we will similarly consider an ℓ q -norm where q ∈ [1, +∞]. It must be emphasized that this choice makes sense when dealing with regression problems. In this context, the ℓ 2 or the ℓ 1 norms are frequently used as training loss functions. On the other hand, the ℓ +∞ norm is also a popular measure when dealing with reliability issues.

. Optimization formulation

In the described setting, the design of the attacker can be formulated as the problem of finding the "worst perturbation" e such that e ∈ Argmax where C p,δ is the closed and convex set defined as

C p,δ = {e ∈ R N 0 | ∥Σ -1/2 e∥ p ⩽ δ}. (3.2)
Σ ∈ R N 0 ×N 0 is symmetric positive definite matrix. δ is a parameter that controls the maximum allowed perturbation, and Σ is a weighting matrix typically corresponding to the covariance matrix of the inputs. For instance, if we assume it to be a diagonal matrix, it simply introduces a normalization of the perturbation components with respect to the standard deviations of the associated inputs.

For standard choices of activation functions, T is a continuous function. Using the Weierstrass theorem, the existence of a solution (not necessarily unique) to Problem 3.1 is then ensured. Although C p,δ is a relatively simple convex set, this problem appears as a difficult non-convex problem due to the fact that i) T is a complex nonlinear operator, ii) we maximize an ℓ q measure which, in addition, leads to a non-smooth cost function when q = 1 or q = +∞. A further difficulty is that we usually need to attack a large dataset to evaluate the robustness of a network, and the provided optimization algorithm should therefore be fast.

. Algorithm

We propose to implement a two-step approach.

• Step 1. We first perform a linearization based on the following first-order Taylor expansion:

T(x + e) ≃ T(x) + J(x)e, (3.3)
where J(x) ∈ R N m ×N 0 is the Jacobian of the network at x.1 Note that J(x) can be computed by classical back-propagation techniques. We will make a second approximation, that is y ≃ T(x). Based on these two approximations and after the variable change e ′ = δ -1 Σ -1/2 e, Problem 3.1 simplifies to

maximize e ′ ∈B p ∥J(x)Σ 1/2 e ′ ∥ q , (3.4)
where B p is the closed ℓ p ball centered at 0 and with unit radius. Note that the optimal cost value in Eq. (3.4) is the subordinate norm of matrix J(x)Σ 1/2 when the input space is equipped with the ℓ p norm and the output space with the ℓ q one. We recall that this subordinate norm is defined, for every matrix M ∈ R N m ×N 0 , as

∥M∥ p,q = sup z∈R N 0 \{0}
∥Mz∥ q ∥z∥ p .

(3.5)

Problem 3.4 is thus equivalent to find a vector e ′ for which the value of the cost function is equal to ∥J(x)Σ 1/2 ∥ p,q . For values of (p, q) listed below the expression of such vector has an explicit form.

-If p = q = 2, e ′ is any unit ℓ 2 norm eigenvector of Σ 1/2 J(x) ⊤ J(x)Σ 1/2 associated with the maximum eigenvalue of this matrix. This vector can be computed by performing a singular value decomposition of J(x)Σ 1/2 .

-If p = 2 and q = +∞, e ′ is any unit ℓ 2 norm vector colinear with a row of J(x)Σ 1/2 having maximum ℓ 2 norm.

-If p = +∞ and q = +∞, e ′ is a unit norm vector whose elements are equal to (ϵ (i) ) 1⩽i⩽N 0 where, for every i ∈ {1, . . . , N 0 }, ϵ i ∈ {-1, 0, 1} is the sign of the i-th element of a row of J(x)Σ 1/2 with maximum ℓ 1 norm.

-If p = 1 and q = 1, e ′ is a vector which has only one nonzero component equal to ±1, the index of this component corresponds to the column of J(x)Σ 1/2 with maximum ℓ 1 norm.

-If p = 1 and q = 2, e ′ is a vector with only one nonzero component equal to ±1. The index of this component corresponds to a column of J(x)Σ 1/2 with maximum ℓ 2 norm.

-If p = 1 and q = +∞, e ′ is again a vector with only one nonzero component equal to ±1. The index of this component corresponds to a column of J(x)Σ 1/2 where is located an element of maximum absolute value.

• Step 2. In the previous optimization step, the optimal solution is not unique.

Indeed if e = δΣ 1/2 e ′ is a solution to Problem 3.4, thene is also a solution. In addition, there may exist other reasons for the multiplicity of the solutions. For example, there may be several maximum norm rows for matrix J(x)Σ 1/2 . Among all the possible choices, we propose to choose the solution e leading to the maximum deviation w.r.t. the ground truth, that is such that ∥T(x + e) -y∥ q is maximum. This requires to perform a search on a small number of possible candidates. Note that no approximation error is involved in this step. If the ground truth for the output is not available, it can be replaced by the model output.

• Post-optimization. If 1 < q < +∞ and T is assumed to be differentiable, e → ∥T(x + e) -y∥ q q is a differentiable function. A further refinement consists of minimizing this function over C p,δ by using a projected gradient algorithm with Armijo search for the step-size. The previous estimates of e can then be used to initialize the algorithm. According to our numerical tests, implementing this strategy when q = 2 only brings a marginal improvement. Moreover, this approach cannot be used when q = 1 or q = +∞.

. Attacking a group of inputs

It can also be interesting to attack only a selected subset of inputs. It may help in identifying the more sensitive inputs of the network. Also, for some inputs like unsorted categorical ones, attacks are often meaningless since they introduce a main change in the informative contents of the dataset, which can be easily detected. Our proposed approach can be adapted to generate such partial attacks. In Problem 3.4, it is indeed sufficient to replace matrix Σ 1/2 by DΣ 1/2 D, where D a masking diagonal matrix whose diagonal elements are equal to 1 when the input is attacked and 0 otherwise. The optimal solutions e ′ and e = δDΣ 1/2 D e ′ have then their components equal to 0 for the non-attacked inputs. Note that the naive approach which would consist in solving Eq. (3.4) and setting to zero the resulting perturbation components for non-attacked inputs would be sub-optimal.

. Error Metrics

Analogous to metrics such as fooling rates/success rates in classification tasks we propose three metrics to quantify the performance of the proposed adversarial attacker in regression tasks. These metrics formulations are given in Table 3.1. Let K be the number of samples in the dataset.

• Mean Accuracy Error (MAE) : It is the average of the distance between the ground truth of the sample and the prediction of the attacked sample over K samples.

• Fooling Error (E) : It is the average of the distance between the prediction of the original sample and the prediction of the attacked sample over K samples.

• Symmetric Mean Accuracy Error (SMAPE) : This measure takes into account relative error distances of perturbed samples and original samples from the ground truth. This value is averaged over samples where distance of ground truth is farther from perturbed samples than the original sample. For all the above metrics, a high value implies a strong attack.

Mean Accuracy Error

MAE = 1 K K ∑ k=1 ∥T(x k + e k ) -y k ∥ q Fooling Error E = 1 K K ∑ k=1 ∥T(x k + e k ) -T(x k )∥ q Symmetric Mean Accuracy Percentage Error SMAPE = 2 K + K+ ∑ k=1 ∥T(x k + e k ) -y k ∥ q -∥T(x k ) -y k ∥ q ∥T(x k + e k ) -y k ∥ q + ∥T(x k ) -y k ∥ q Table 3
.1: Error metrics used for evaluation. The mean value computed for SMAPE is limited to the K + positive values of the elements in the summation. e k is the perturbation generated by the adversarial attacker on the k-th sample in the dataset of length K.

. Experiments

. Open Source Datasets

We run our experiments on three open source regression datasets. The Combined Cycle Power Plant [START_REF] Tüfekci | Prediction of full load electrical power output of a base load operated combined cycle power plant using machine learning methods[END_REF] dataset has 4 features with 9,568 instances. The task is to predict the net hourly electrical energy output using hourly average ambient variables. The Red Wine Quality dataset [START_REF] Cortez | Modeling wine preferences by data mining from physicochemical properties[END_REF] contains 1,599 total samples and each instance has 11 features. The features are physio-chemical and sensory measurements for wine. The output variable is a quality score ranging from 0 to 10, where 10 represents for best quality and 0 for least quality. For the Abalone dataset, the task is to model an Abalone's age based purely on its physical measurements. This would allow Abalone's age estimation without cutting its shell. There are in total 4,177 instances with 8 input variables including one categorical variable. The datasets are divided with a ratio of 4:1 between training and testing data. The categorical attributes are dealt with by using one hot encoding based on the number of categories. The input attributes are normalized by removing their mean and scaling to unit variance.

We train fully connected networks for the estimation of variables from the datasets. The network architecture for the dataset are given below. The values represent the number of hidden neurons in the layers. The activation function at each layer is ReLU except for the last layer.

• Combined cycle Power Plant dataset -(10, 6, 1)

• Red Wine Quality dataset -(100, 100, 100, 10, 1)

• Abalone Data set -(256, 256, 256, 256, 1)

. Industrial Dataset -Safety Critical Application

An industrial application dataset is also considered with 2,219,097 training, 739,639 validation, and 739,891 test samples. The description of the input/output variables of the dataset is given in Table 3.2. The variable to be predicted is the Estimation of Arrival time (ETE) of a flight, given variables including the distance and speed, and also an initial estimate of ETE. The dataset is related to flight control, an activity area where safety is critical. The input attributes are normalized by removing their mean and scaling to unit variance. For models, we build fully connected networks with ReLU activation function on all the hidden layers except the last one. The network architecture is shown in the Figure 3.1.

. Adversarial attacks on Standard Training

For checking the efficacy of the proposed adversarial attacker, we first train the neural networks without any constraints using the network architecture presented in the previous section with the aim of reducing the prediction/performance loss on the train dataset. This is refered to as a standard training procedure. We compare the proposed adversarial attacker with random noise attackers generated by i.i.d. perturbations. We use three additive noise distributions: Gaussian, uniform, and binary for comparisons. The output of these attackers have been normalized so as to meet the desired bound on the norm of the perturbation. The metrics (MAE, E and SMAPE) are computed on the test samples where K is the total number of samples in the test set. The results on the 4 datasets for varying noise levels for ℓ 2 perturbations are shown in Table 3.3. Our proposed attacker successfully attacks the neural network models, and is significantly better than random noise attackers in terms of distance metrics.

We also show the histograms of (∥T( [START_REF] Mehrabi | A survey on bias and fairness in machine learning[END_REF], where (e k ) 1⩽k⩽K have been generated from various noise distributions and the proposed adversarial attacker. From the histograms, we observe the proposed attacker is guaranteed to attack the sample as it generates the worst positive deviation, while the random noise attackers are symmetric around zero. This implies that random attackers sometimes make the prediction better rather than making it deviate from the ground truth predictions. We observe similar trends for all four test datasets in the histograms.

x k + e k ) -y k ∥ q -∥T(x k ) -y k ∥ q ) 1⩽k⩽K in Figures 3.2-3.
Noise

1 × 10 -1 2 × 10 -1 1 × 10 -1 2 × 10 -1 5 × 10 -2 1 × 10 -1 1 × 10 -1 2 × 10 -

. Attacks on Spectral Normalization training

In these experiments, we train our networks while using a spectral normalization technique as shown in section 2.3.2.1 which has been proven to be very effective of Fooling Error (E) and Symmetric Mean Accuracy Percentage error (SMAPE) for various datasets and varying perturbation magnitude under ℓ 2 perturbations is given in Table 3.4. We observe that the spectral normalized models are more resilient to attacks than the standard trained models. Thus spectral normalization constraints help to generate more robust and safer models.

Noise

. Attacking only continuous variables

All the previous results have been obtained with attack and noise addition on all the input features present in the datasets. As pointed in Section 3.3.4, the introduced adversarial attacker is capable of attacking a group of inputs. While generating an adversarial attack, we avoid attacking the categorical input variables [START_REF] Ballet | Imperceptible adversarial attacks on tabular data[END_REF]. Attacking categorical inputs may introduce a major change in the informative contents of the dataset, which is often easily identifiable and detectable. This is clearly in contradiction with the purpose of generating adversarial attacks which are hardly noticeable by the user. Hence in Abalone and industrial datasets, we attack only the continuous variables. For the Combined Power plant dataset, we attack 3 out of 4 continuous variables since it does not contain any categorical variables. Similarly, for the Red-wine dataset we attack 8 continuous variables out of 11. The performance of the adversarial attacker, when attacking only few inputs, is shown in Table 3.5. We observe that, from the attacker viewpoint, attacking only few inputs has an advantage since some inputs can be more sensitive than others and can contribute more in generating perturbed outputs. Attacking all the inputs might add redundancy to the generated perturbation and hence render the attack more easily detectable.

. Comparison on different measures of deviations

As emphasized in Section 3.3.3, our adversarial attacker is applicable for various measures of perturbation on input and output deviations. The previous results have been obtained for the value p = q = 2 termed as ℓ 2 attacks here. We further show results for p = q = 1 termed as ℓ 1 attacks and for p = q = +∞ termed as ℓ ∞ attacks in Table 3.6. We observe that our proposed attacker us successful in attacking under different measures of perturbations and remains better than random noise attackers in all cases and for all metrics.

Noise

. Summary

The contributions of this work are summarised as follows:

• We propose a simple, novel, and easily implementable Jacobian-based adversarial attacker for regression tasks in the context of white box attacks.

• We define three error metrics: Mean Accuracy Error, Fooling Error, and Symmetric Mean Accuracy Percentage Error that help to analyze the attacker's effectiveness.

Noise • Our attacker is generic in the sense that it can handle any measure (ℓ 1 , ℓ 2 , ℓ ∞ ) on input or output perturbations according to the target application.

1 × 10 -1 2 × 10 -1 1 × 10 -1 2 × 10 -1 1 × 10 -1 2 × 10 -1 Attacks ℓ 2 ℓ 1 ℓ ∞ MAE gauss (×10 -
• We show that the proposed attacker allows us to attack any given subset of input features optimally. This feature may be helpful when handling specific tabular datasets and also be insightful when information regarding sensitivity or the ability to control some inputs is available.

• We evaluate our results on open-source regression datasets and an industrial dataset (output and input features described in Table 3.2) which lies in the domain of safety-critical applications.

• Our tests concentrated on fully connected networks, but it is worth pointing out that the proposed approach can be applied to any network architecture.

Chapter 4 Multivariate Estimation of Stability of Neural Networks

. Introduction

In this chapter, we introduce a multivariate Lipschitz constant-based stability analysis of fully connected neural networks allowing us to capture the influence of each input or group of inputs on the neural network stability. Our approach relies on a suitable re-normalization of the input space, with the objective to perform a more precise analysis than the one provided by a global Lipschitz constant. We investigate the mathematical properties of the proposed multivariate Lipschitz analysis and show its usefulness in better understanding the sensitivity of the neural network with regard to groups of inputs. We display the results of this analysis by a new representation designed for machine learning practitioners and safety engineers termed as a Lipschitz star. The Lipschitz star is a graphical and practical tool to analyze the sensitivity of a neural network model during its development, with regard to different combinations of inputs. By leveraging this tool, we show that it is possible to build stable-by-design models using spectral normalization techniques for controlling the stability of a neural network.Thanks to our multivariate Lipschitz analysis, we can also measure the efficiency of adversarial training in inference tasks. We used the adversarial attacker proposed in Chapter 3 to generate adversarial attack and perform adversarial training. We perform experiments on various open access tabular datasets, and also on a real Thales Air Mobility industrial application subject to certification requirements, also considered in the previous chapter.

. Limitations of Previous Works

The main limitations of the Lipschitz constant, defined in either global or local context, is that it only provides a single parameter to quantify the stability of a neural network. Such a single-parameter analysis does not facilitate the understanding of potential sources of instability. In particular, it may be insightful to identify the inputs which have the highest impact in terms of sensitivity on the output. In the context of tabular data mining, the inputs often have quite heterogeneous characteristics. Some of them are categorical data, often encoded in a specific way (e.g., one-hot encoder [START_REF] Hancock | Survey on categorical data for neural networks[END_REF]) and among them, one can usually distinguish those which are unsorted (like labels identifying countries) or those which are sorted (like severity scores in a disease). So, it may appear useful to analyze in a specific manner each type of inputs of a NN and even sometimes to exclude some of these inputs (e.g., unsorted categorical data for which the notion of small perturbation may be meaningless) from the performed sensitivity analysis. A comparison of the state-of-the-art and proposed approach is presented in Table 4

. Partial/Weighted Lipschitz Constant

We present in this section a new approach based on a suitable weighting operation performed in the computation of Lipschitz constants. This enables a multivariate sensitivity analysis of the neural network stability for individual inputs or groups of inputs. We will start by motivating this weighting from a statistical standpoint. Then we will define it in a more precise manner, before discussing its resulting mathematical properties.

. Statistical motivations

For tractability, assume that the perturbation at the network input is a realization of a zero-mean Gaussian distributed random vector z with N 0 × N 0 covariance matrix Σ ≻ 0. Then, its density upper level sets are defined as

C η = {z ∈ R N 0 | z ⊤ Σ -1 z ⩽ η}, (4.1)
for every η ∈]0, +∞[. The set C η defines an ellipsoid where the probability density takes its highest values. More precisely, the probability for z to belong to this set is independent of Σ and is equal to

P(z ∈ C η ) = γ(N 0 /2, η/2) Γ(N 0 /2) , (4.2)
where Γ is the gamma function and γ the lower (unnormalized) incomplete gamma function.

Proof of Eq. (4.2): Let B(0, √ η) be the closed ball of R N 0 with center 0 and radius √ η and let S N 0 -1 be the surface of the unit-radius sphere in R N 0 -1 . We have

P(z ∈ C η ) = 1 (2π) N 0 /2 (det(Ω)) 1/2 C η exp - 1 2 z ⊤ Ω -1 z dz = 1 (2π) N 0 /2 B(0, √ η) exp - 1 2 ∥u∥ 2 du = 1 (2π) N 0 /2 S N 0 -1 √ η 0 ρ N 0 -1 exp - ρ 2 2 dρ = 1 (2π) N 0 /2 2(π) N 0 /2 Γ(N 0 /2) √ η 0 ρ N 0 -1 exp - ρ 2 2 dρ = 1 Γ(N 0 /2) η 2 0 ζ N 0 /2-1 exp(-ζ) dζ = γ(N 0 /2, η/2) Γ(N 0 /2) . (4.3)
On the other hand, let us assume that the maximum standard deviation σ max of the components of z (i.e., square root of the maximum diagonal element of matrix Σ) is small enough. If we suppose that the network T is differentiable in the neighborhood of a given input x ∈ R N 0 , as the input perturbation is small enough, we can approximate the network output by the following first-order expansion:

T(x + z) ≃ T(x) + T ′ (x)z. (4.4)
Let us focus our attention on perturbations in C η . By doing so, we impose some norm-bounded condition, which may appear realistic for adversarial perturbations. Then, we will be interested in calculating

sup z∈C η ∥T(x + z) -T(x)∥ ≃ sup z∈C η ∥T ′ (x)z∥. (4.5)
By making the variable change z ′ = z/ √ η and using Eq. (4.1),

sup z∈C η ∥T(x + z) -T(x)∥ ≃ √ η sup z ′ ∈C 1 ∥T ′ (x)z ′ ∥ = √ η sup z∈R N 0 z ′ ̸ =0 ∥T ′ (x)z ′ ∥ ∥z ′ ∥ Σ -1 = √ ησ max sup z ′ ∈R N 0 z ′ ̸ =0 ∥T ′ (x)z ′ ∥ ∥z ′ ∥ Ω -1 , (4.6)
where Ω = Σ/σ 2 max and ∥ • ∥ Ω -1 = (•) ⊤ Ω -1 (•). This suggests that, in this context, the suitable subordinate matrix norm for computing the Lipschitz constant is obtained by weighting the Euclidean norm in the input space with Ω -1 . We can also deduce from Eq. (4.6), by setting z ′′ = Ω -1/2 z ′ , that

sup z∈C η ∥T(x + z) -T(x)∥ ≃ √ ησ max sup z ′′ ∈R N 0 z ′′ ̸ =0 ∥T ′ (x)Ω 1/2 z ′′ ∥ ∥z ′′ ∥ = √ η∥T ′ (x)Σ 1/2 ∥ S . (4.7)
On the other hand, based on the first-order approximation in Eq. (4.4), T(x + z) is approximately Gaussian with mean T(x) and covariance matrix

T ′ (x)ΣT ′ (x) ⊤ . As ∥T ′ (x)ΣT ′ (x) ⊤ ∥ S = ∥T ′ (x)Σ 1/2 ∥ 2 S
, we see that another insightful interpretation of Eq. (4.7) is that, up to the scaling factor √ η, it approximately delivers the square root of the spectral norm of the covariance matrix of the output perturbations.

. New definition of a weighted Lipschitz constant

Based on the previous motivations, we propose to employ a weighted norm to define a Lipschitz constant of the network as follows:

Definition 5 Let Ω be an N 0 × N 0 symmetric positive definite real-valued matrix. We say that θ Ω m is an Ω-weighted norm Lipschitz constant of T as described in Fig. 2.1 if

(∀(x, z) ∈ (R N 0 ) 2 ) ∥T(x + z) -T(z)∥ ⩽ θ Ω m ∥z∥ Ω -1 . (4.8)
The above definition can be extended to non Euclidean norms by making use of exponents (p, q) ∈ [1, +∞] 2 and by replacing inequality Eq. (4.8) with

(∀(x, z) ∈ (R N 0 ) 2 ) ∥T(x + z) -T(z)∥ q ⩽ θ Ω m ∥Ω -1/2 z∥ p .
(4.9)

By changes of variable, this inequality can also be rewritten as

(∀(x ′ , z ′ ) ∈ (R N 0 ) 2 ) ∥T(Ω 1/2 (x ′ + z ′ )) -T(Ω 1/2 z ′ )∥ q ⩽ θ Ω m ∥z ′ ∥ p . (4.10)
Therefore, we see that calculating θ Ω m is equivalent to deriving a Lipschitz constant of the network T where an additional first linear layer Ω 1/2 has been added. Throughout the rest of this section, it will be assumed that, for every i ∈ {1, . . . , m -1} the activation operator R i is separable and α i -averaged. It then follows from Eq. (2.23) that an Ω-weighted norm Lipschitz constant of T is

ϑ Ω m = sup Λ 1 ∈D N 1 ({2α 1 -1,1}), . . ., Λ m-1 ∈D N m-1 ({2α m-1 -1,1}) ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 ∥ p,q . (4.11)
Although all our derivations were based on the fact that Ω is positive definite, from the latter expression we see that, by continuous extension, ϑ Ω m can be defined when Ω is a singular matrix.

. Sensitivity with respect to a group of inputs

In this section, we will be interested in a specific family of weighted norms associated with the set of matrices

{Ω ϵ,K | ∅ ̸ = K ⊂ {1, . . . , N 0 }, ϵ ∈]0, 1]},
defined, for every nonempty subset K of {1, . . . , N 0 } and for every ϵ ∈]0, 1], as

Ω ϵ,K = Diag(σ 2 ϵ,K,1 , . . . , σ 2 ϵ,K,N 0 ), (4.12)
where

(∀ℓ ∈ {1, . . . , N 0 }) σ ϵ,K,ℓ = 1 if ℓ ∈ K ϵ otherwise. (4.13)
If we come back to the statistical interpretation in Section 4.3.1, Ω ϵ,K is then (up to a positive scale factor) the covariance matrix of a Gaussian random vector z with independent components. 1 The components with indices in K have a given variance σ2 max while the others have variance ϵ 2 σ 2 max . Such a matrix thus provides a natural way of putting emphasis on the group of inputs with indices in K. Thus, variables ϑ Ω ϵ,K m will be termed partial Lipschitz constants in the following. The next proposition lists the main properties related to the use of such weighted norms for calculating Lipschitz constants.

Proposition 6 Let (p, q) ∈ [1, +∞] 2 . For every nonempty subset K of {1, . . . , N 0 } and for every ϵ ∈]0, 1], let Ω ϵ,K be defined as above and let ϑ Ω ϵ,K m be defined by Eq. (4.11). Let K 0 and K 1 be nonempty subsets of {1, . . . , N 0 }. Then the following hold:

i) As ϵ → 0, ϑ Ω ϵ,K 0 m
converges to the Lipschitz constant of a network where all the inputs with indices out of K 0 are kept constant.

ii) ϑ Ω 1,K 0 m is equal to the global Lipschitz constant ϑ m defined by Eq. (2.23).

iii) Let (ϵ, ϵ ′ ) ∈]0, 1] 2 . If Ω ϵ,K 0 ⪯ Ω ϵ ′ ,K 1 , 2 then ϑ Ω ϵ,K 0 m ⩽ ϑ Ω ϵ ′ ,K 1 m . iv) Function ϑ Ω •,K 0 m : ]0, 1] → [0, +∞[ : ϵ → ϑ Ω ϵ,K 0 m is monotonically increasing. v) Let ϵ ∈]0, 1]. If K 0 ⊂ K 1 , then ϑ Ω ϵ,K 0 m ⩽ ϑ Ω ϵ,K 1 m . vi) Let ϵ ∈]0, 1], let K ∈ N \ {0}, and let ω K,ϵ = N 0 -1 K -1 1 + N 0 K -1 ϵ . (4.14)
We have

max K⊂{1,...,N 0 } card K=K ϑ Ω ϵ,K m ⩽ ϑ m ⩽ 1 ω K,ϵ ∑ K⊂{1,...,N 0 } card K=K ϑ Ω ϵ,K m . (4.15)
vii) Let ϵ ∈]0, 1], let P be a partition of {1, . . . , N 0 }, and let ω P,ϵ = 1 + (card P -1)ϵ.

We have

max K∈P ϑ Ω ϵ,K m ⩽ ϑ m ⩽ 1 ω P,ϵ ∑ K∈P ϑ Ω ϵ,K m . (4.16) viii) Let K 2 be such that K 1 ∩ K 2 ̸ = p and K 1 ∪ K 2 = K 0 . Let p * ∈ [1, +∞] be such that 1/p + 1/p * = 1 Then ϑ Ω ϵ,K 0 m ⩽ (ϑ Ω ϵ,K 1 m ) p * + (ϑ Ω ϵ,K 2 m ) p * 1/p * + o(ϵ).
(4.17)

Proofs: i): As ϵ → 0, the diagonal elements of Ω ϵ,K 0 with index ℓ ̸ ∈ K tend to zero, which in Eq. (4.11) amounts to assuming that the corresponding components of the input perturbation are zero. The existence of the limit Ω 0,K 0 is secured based on the remark at the end of Section4.3.2.

ii): If ϵ = 1, then Ω ϵ,K 0 = Id N 0 and Eq. (4.11) reduces to Eq. (2.23).

iii): For every i ∈ {1, . . . , m -1}, let Λ i ∈ D N i ({2α i -1, 1}). Then, by using the triangle inequality,

∥W m Λ m-1 • • • Λ 1 W 1 (Ω ϵ,K 0 ) 1/2 ∥ p,q ⩽ ∥W m Λ m-1 • • • Λ 1 W 1 (Ω ϵ ′ ,K 1 ) 1/2 ∥ p,q ∥(Ω ϵ ′ ,K 1 ) -1 Ω ϵ,K 0 ∥ 1/2 p,p . (4.18)
By taking the supremum of both sides with respect to the (Λ i ) 1⩽i⩽m-1 matrices, we deduce that

ϑ Ω ϵ,K 0 m ⩽ ∥(Ω ϵ ′ ,K 1 ) -1 Ω ϵ,K 0 ∥ 1/2 p,p ϑ Ω ϵ ′ ,K 1 m . ( 4

.19)

On the other hand,

(Ω ϵ ′ ,K 1 ) -1 Ω ϵ,K 0 = Diag σ 2 ϵ,K 0 ,1 σ 2 ϵ ′ ,K 1 ,1 , . . . , σ 2 ϵ,K 0 ,N 0 σ 2 ϵ ′ ,K 1 ,N 0 , (4.20)
where, by using the fact that

Ω ϵ,K 0 ⪯ Ω ϵ ′ ,K 1 , (∀ℓ ∈ {1, . . . , N 0 }) σ ϵ,K 0 ,ℓ σ ϵ ′ ,K 1 ,ℓ ⩽ 1. (4.21)
Since (Ω ϵ ′ ,K 1 ) -1 Ω ϵ,K 0 is a diagonal matrix with elements lower than or equal to 1, ∥(Ω ϵ ′ ,K 1 ) -1 Ω ϵ,K 0 ∥ p,p ⩽ 1 and it follows from Eq. (4.19) that

ϑ Ω ϵ,K 0 m ⩽ ϑ Ω ϵ ′ ,K 1 m . (4.22) iv): Let (ϵ, ϵ ′ ) ∈]0, 1] 2 with ϵ < ϵ ′ . We have Ω ϵ,K 0 ⪯ Ω ϵ ′ ,K 0 and, according to iii), ϑ Ω ϵ,K 0 m ⩽ ϑ Ω ϵ ′ ,K 0 m . (4.23) v): If K 0 ⊂ K 1 , then Ω ϵ,K 0 ⪯ Ω ϵ,K 1 and the result follows from iii). vi): We have ∑ K⊂{1,...,N 0 } card K=K Ω 1/2 ϵ,K = N 0 -1 K -1 + N 0 K - N 0 -1 K -1 ϵ Id N 0 .
(4.24)

By using the relation

N 0 K = N 0 K N 0 -1 K -1 , (4.25)
we deduce that

∑ K⊂{1,...,N 0 } card K=K Ω 1/2 ϵ,K = ω K,ϵ Id N 0 . (4.26)
For every i ∈ {1, . . . , m -1}, let Λ i ∈ D N i ({2α i -1, 1}). Then, ∥W m Λ m-1 • • • Λ 1 W 1 ∥ p,q = 1 ω K,ϵ W m Λ m-1 • • • Λ 1 W 1 ∑ K⊂{1,...,N 0 } card K=K Ω 1/2 ϵ,K p,q ⩽ 1 ω K,ϵ ∑ K⊂{1,...,N 0 } card K=K ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 ϵ,K ∥ p,q . (4.27)
We deduce that

ϑ m = sup Λ 1 ∈D N 1 ({2α 1 -1,1}), . . ., Λ m-1 ∈D N m-1 ({2α m-1 -1,1}) ∥W m Λ m-1 • • • Λ 1 W 1 ∥ p,q ⩽ 1 ω K,ϵ ∑ K⊂{1,...,N 0 } card K=K sup Λ 1 ∈D N 1 , . . ., Λ m-1 ∈D N m-1 ∥W m • • • Λ 1 W 1 Ω 1/2 ϵ,K ∥ p,q = 1 ω K,ϵ ∑ K⊂{1,...,N 0 } card K=K ϑ Ω ϵ,K m .
(4.28)

Furthermore, according to ii) and iv), for every K ⊂ {1, . . . , N 0 },

ϑ Ω ϵ,K m ⩽ ϑ Ω 1,K m = ϑ m . (4.29)
This yields max

K⊂{1,...,N 0 } card K=K ϑ Ω ϵ,K m ⩽ ϑ m . (4.30)
vii): The proof is similar to that of vi) by noticing that, if P is a partition of {1, . . . , N 0 }, then

∑ K∈P Ω 1/2 ϵ,K = ω P,ϵ Id N 0 . (4.31) viii): For every ∅ ̸ = K ⊂ {1, . . . , N 0 }, ϑ Ω ϵ,K m = ϑ Ω 0,K m + o(ϵ).
It is thus sufficient to prove this inequality in the limit case when ϵ → 0. For every i ∈ {1, . . . , m -1}, let Λ i ∈ D N i ({2α i -1, 1}). Let x ∈ R N 0 and let x K be the projection of x onto the space of vectors whose components indexed by {1, . . . , N 0 } \ K are zero. We have thus

x K 0 = x K 1 + x K 2 . Then, ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 0,K 0 x∥ q = ∥W m Λ m-1 • • • Λ 1 W 1 x K 0 ∥ q ⩽ ∥W m Λ m-1 • • • Λ 1 W 1 x K 1 ∥ q + ∥W m Λ m-1 • • • Λ 1 W 1 x K 2 ∥ q = ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 0,K 1 x K 1 ∥ q + ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 0,K 2 x K 2 ∥ q ⩽ ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 0,K 1 ∥ p,q ∥x K 1 ∥ p + ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 0,K 2 ∥ p,q ∥x K 2 ∥ p .
(4.32)

By using Hölder's inequality, we deduce that

∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 0,K 0 x∥ q ⩽ ∥W m • • • Λ 1 W 1 Ω 1/2 0,K 1 ∥ p * p,q + ∥W m • • • Λ 1 W 1 Ω 1/2 0,K 2 ∥ p * p,q 1/p * (∥x K 1 ∥ p p + ∥x K 2 ∥ p p ) 1/p . (4.33) Since ∥x K 1 ∥ p p + ∥x K 2 ∥ p p = ∥x K 0 ∥ p p , it follows that ∥W m Λ m-1 • • • Λ 1 W 1 Ω 1/2 0,K 0 ∥ p,q ⩽ ∥W m • • • Λ 1 W 1 Ω 1/2 0,K 1 ∥ p * p,q + ∥W m • • • Λ 1 W 1 Ω 1/2 0,K 2 ∥ p * p,q 1/p * (4.34)
Taking the supremum with respect to (Λ i ) 1⩽i⩽m-1 and majorizing the supremum of the sum in the right-hand side by the sum of the suprema yield Eq. (4.17).

Let us comment on these results. According to Property i) in the limit case when ϵ → 0, only the inputs with indices in K 0 are used in the computation of the associated Lipschitz constant. In turn, Property ii) states that, when ϵ = 1, we recover the classical expression of a Lipschitz constant where the perturbations on all the inputs are taken into account. In addition, based on Property iv), the evolution of ϑ Ω ϵ,K 0 m when ϵ varies from 1 to 0 provides a way of assessing how the group of inputs indexed by K 0 contributes progressively to the overall Lipschitz behaviour of the network. Although one would expect that summing the Lipschitz constants obtained for each group of inputs would yield the global Lipschitz constant, Properties vi) and vii) show that this does not hold in general whatever the way the entries are split (possibly overlapping groups of given size K or disjoint groups of arbitrary size). Instead, after suitable normalization, such sums provide upper bounds on ϑ m . Furthermore, it follows from ii), Eq. (4.15), and Eq. (4.16) that the difference between these normalized sums and ϑ m tends to vanish when ϵ increases. Note that, when looking at the sensitivity with respect to individual inputs, i.e., when the considered set of indices are singletons, both vi) (with K = 1) and vii) (with P = {{k} | k ∈ {1, . . . , N 0 }} lead to the same inequality

max k∈{1,...,N 0 } ϑ Ω ϵ,{k} m ⩽ ϑ m ⩽ 1 1 + (N 0 -1)ϵ N 0 ∑ k=1 ϑ Ω ϵ,{k} m . (4.35)

. Lipschitz Star

We propose a new representation for displaying the results of the Lipschitz analysis of a neural network utilizing radial plots. More precisely, we plot the values of (ϑ

Ω ϵ,{k} m
) 1⩽k⩽N 0 on a star or radar chart where each branch of the star corresponds to the index k of an input. For each value of ϵ, a new plot is obtained which is displayed in a specific color. Note that, according to Proposition 6iv), the plots generated for different ϵ values cannot cross. When ϵ = 1, we obtain an "isotropic" representation whose "radius" corresponds to the global Lipschitz constant ϑ m of the network. This representation is called a Lipschitz star. All the results of our analysis will be displayed using this representation. To verify the correctness of the partial Lipschitz computation, we first study simple synthetic examples of polynomial systems for which we can calculate explicitly the partial Lipschitz constants. We generate input-output data for the defined systems, and train a fully connected model using a standard training, i.e., without any constraints. We compare this approach with a training subject to a spectral norm constraint on the layers.

. Experiments

Spectral Normalization: For safety critical tasks, Lipschitz constant and performance targets can be specified as engineering requirements, prior to network training. A Lipschitz target can be defined by a safety analysis of the acceptable perturbations for each output knowing the input range and it constitutes a current practice in many industries. As mentioned in the previous chapter, 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. In our experiments, we train networks while using a spectral normalization technique. We recall from section 2.3.2.1 for practical application of spectral normalization, given an m layer fully connected architecture and a Lipschitz target θ, we can constrain the spectral norm of each layer to be less than m √ θ. According to Eq. (2.7), this ensures that the upper bound on the global Lipschitz constant is less than θ.

For each training, we study the effect of input variables on the stability of the networks. As proposed in Section 4.3.3, for a given group of inputs with indices in K, we will quantify the partial Lipschitz constant ϑ

Ω ϵ,K m . The obtained value of ϑ Ω ϵ,K m
allows us to evaluate how the corresponding group of variables may potentially affect the stability of the network. For simplicity, in the next section, we will focus on the limit case when ϵ = 0 (see the last remark in Section 4.3.2).

Partial Lipschitz constant values ϑ Ω 0,K m , for all possible choices for K, are computed using the numerical method described in Section 2.2.1.4 and compared with the theoretical values derived in the following subsection. More details on the models are also provided in these sections.

. Polynomial systems

We consider regression problems where the data is synthesized by a secondorder multivariate polynomial. The system to be modelled is thus described by the following function:

(∀(ξ 1 , . . . , ξ N 0 ) ∈ R N 0 ) f (ξ 1 , . . . , ξ N 0 ) = N 0 ∑ k=1 a k ξ k + N 0 ∑ k=1 N 0 ∑ l=1 b k,l ξ k ξ l , (4.36)
where (a k ) k∈N 0 and (b k,l ) 1⩽k,l⩽N 0 are the real-valued polynomial coefficients. Note that, such a polynomial system is generally not Lipschitz-continuous. The Lipschitzcontinuity property only holds on every compact set. Subsequently, we will thus study this system on the hypercube [-M, M] N 0 with M > 0.

The explicit values of the partial Lipschitz constant on this domain can be derived as follows. We first calculate the gradient of f

∇ f (ξ 1 , . . . , ξ N 0 ) = ∂ k f (ξ 1 , . . . , ξ N 0 ) 1⩽k⩽N 0 , (4.37)
where, for every k ∈ {1, . . . , N 0 }, ∂ k f denotes the partial derivative w.r.t. the k-th variable given by

∂ k f (ξ 1 , . . . , ξ N 0 ) = a k + N 0 ∑ l=1 (b k,l + b l,k )ξ l .
(4.38)

For every K ⊂ {1, . . . , N 0 }, the partial Lipschitz constant θΩ 0,K of the polynomial system (restricted to [-M, M] N 0 ) w.r.t. the group of variables with indices in K is then equal to

θΩ 0,K = sup (ξ 1 ,...,ξ N 0 )∈[-M,M] N 0 λ Ω 0,K (ξ 1 , . . . , ξ N 0 ), (4.39)
where, for every diagonal matrix

Λ = Diag(ε 2 1 , . . . , ε 2 N 0 ) with (ε 1 , . . . , ε N 0 ) ∈ [0, +∞[ N 0 , λ Λ (ξ 1 , . . . , ξ N 0 ) = ∥(∇ f (ξ 1 , . . . , ξ N 0 )) ⊤ Λ 1/2 ∥ 2 = N 0 ∑ k=1 ε k ∂ k f (ξ 1 , . . . , ξ N 0 ) 2 .
(4.40)

Since the partial derivatives in Eq. (4.38) are affine functions of the variables (ξ 1 , . . . , ξ N 0 ), λ Λ is a convex function. We deduce that the supremum in Eq. (4.39) is attained when ξ 1 = ±M, . . . , ξ N 0 = ±M, so that θΩ 0,K can be computed by looking for the maximum of a finite number (2 N 0 ) of values.

. Numerical results

In our numerical experiments, we consider a toy example corresponding to

N 0 = 3 and (∀(ξ 1 , ξ 2 , ξ 3 ) ∈ R 3 ) f (ξ 1 , ξ 2 , ξ 3 ) = ξ 1 + 100ξ 3 -ξ 2 2 + γξ 1 ξ 3 , (4.41)
where γ ∈ [0, +∞[. We deduce from Eq. (4.40) that

λ Λ (ξ 1 , ξ 2 , ξ 3 ) = ε 1 (1 + γξ 3 ) 2 + 4ε 2 ξ 2 2 + ε 3 (100 + γξ 1 ) 2 (4.42)
and, consequently,

sup (ξ 1 ,ξ 2 ,ξ 3 )∈[-M,M] 3 λ Λ (ξ 1 , ξ 2 , ξ 3 ) = ε 1 (1 + γM) 2 + 4ε 2 M 2 + ε 3 (100 + γM) 2 .
(4.43)

By looking at the seven possible binary values of (ε 1 , ε 2 , ε 3 ) ̸ = (0, 0, 0), we thus calculate the Lipschitz constant of f with respect to each group of inputs. For example,

• if ε 1 = 1, ε = 0, ε 3 = 0, we calculate θΩ 0,K with K = {1}, i.e., evaluate the sensitivity w.r.t. the first variable;

• if ε 1 = ε 2 = 1, ε 3 = 0, we calculate θΩ 0,K with K = {1, 2}, i.e.
, evaluate the joint sensitivity w.r.t. the first and second variables;

• if ε 1 = ε 2 = ε 3 = 1,
we calculate θΩ 0,K with K = {1, 2, 3}, i.e., evaluate the sensitivity w.r.t. all the variables (global Lipschitz constant).

These Lipschitz constants allow us to evaluate the intrinsic dynamics of the system, that is how it responds when its inputs vary. Our interest will be now to evaluate how this dynamics is modified when the system is modelled by a neural network. To do so, three systems are studied by choosing γ ∈ {0, 1 10 , 1} and M = 50. We generate 5,000 data samples from each system, the input values being drawn independently from a random uniform distribution. While training the neural networks, the dataset is divided with a ratio of 4:1 into training and testing samples. The input is normalized using its mean and standard deviation, while the output is max-normalized. We build neural networks for approximating the systems using two hidden layers (m = 3) with a number of hidden neurons equal to 30 in each layer and ReLU activation functions. The training loss is the mean square error.

For different values of γ, we report the values of the partial Lipschitz constants in Tables 4.2-4.4. The variable θ K corresponds to θΩ 0,K for the analytical value we derived from previous formulas, whereas it corresponds to the Lipschitz constant ϑ Ω 0,K 3

, when computed for the neural network trained either in a standard manner or with a spectral normalization constraint. The value of θ used in the spectral normalization was adjusted to obtain a similar global Lipschitz constant to the polynomial system. In the caption, we also indicate the accuracy in terms of normalized mean square error (NMSE) and normalized mean absolute error (NMAE). These values are slightly higher for constrained training, but remain quite small.

Partial

Comments on the results:

• In general, ξ 3 impacts the output of this system the most, and (ξ 2 , ξ 3 ) mainly account for the global dynamics of the system.

• With standard training, we see that there exists a significant increase of the sensitivity with respect to the input variations, so making the neural network vulnerable to adversarial perturbations.

• By using spectral normalization, it is possible to constrain the global Lipschitz constant of the system to be close to the analytical global value while keeping a good accuracy. One may however notice an increase of the sensitivity to ξ 1 and ξ 3 , and a decrease of the sensitivity to ξ 2 with respect to the original system.

• For all the three models, the values obtained with neural networks follow the same trend, for different groups of inputs, as those observed with the analytical values.

• Although the Lipschitz constant of the neural networks is computed on the whole space and the one of the system on [-50, 50] 3 , our Lipschitz estimates appear to be consistent without resorting to a local analysis.

These observations emphasize the importance of controlling the Lipschitz constant of neural network models through specific training strategies. In addition, we see that evaluating the Lipschitz constant with respect to groups of inputs allow us to have a better understanding of the behaviour of the models.

In this section, we have discussed the proposed method for synthetic datasets. In the next section, the sensitivity analysis will be made on widely used open source datasets and an industrial dataset.

. Open-source use-cases

We study regression problems involving tabular datasets to showcase our proposed multivariate analysis of the stability of neural networks. Tabular data take advantage of heterogeneous sources of information coming from different sensors or data collection processes. We apply our methods on widely used tabular datasets: i) Combined Cycle Power Plant dataset3 which has 4 attributes with 9,568 instances; ii) Auto MPG dataset4 consists of 398 instances with 7 attributes; iii) Boston Housing dataset5 consists of 506 instances with 13 attributes. For Combined Power Plant and Auto MPG datasets, we solve a regression problem with a single output, whereas for Boston Housing dataset we consider a two-output regression problem with "price" and "ptratio :pupil-teacher ratio by town" as the output variables. The attributes in the dataset are a combination of continuous and categorical ones. The datasets are divided with a ratio of 4:1 between training and test data. Information on the input and output attributes are provided in Tables 4.5

. Industrial Dataset

Thales Air Mobility industrial application represents the prediction of the Estimated Time En-route (ETE), meaning the time spent by an aircraft between the take-off and landing, considering a number of variables as described in Table 3 The application is important in air traffic flow management, which is an activity area where safety is critical. The purpose of the proposed sensitivity analysis is thus to help engineers in building safe by design models complying with given safety stability targets. The dataset consists of 2,219,097 training, 739,639 validation, and 739,891 test samples.

For all the models, we build fully connected networks with ReLU6 activation function on all the hidden layers, except the last one. The models are trained on Keras with Tensorflow backend. The initializers are set to Glorot uniform. The network architecture of the different models, number of layers, and neurons are tabulated in Table 4.8. Combined Cycle Power Plant dataset with (10, 6) network architecture is trained with two hidden layers having 10 and 6 hidden neurons, respectively. For Thales Air Mobility industrial application (10 × (30)) implies that the neural network has 10 hidden layers with 30 neurons each.

. Sensitivity analysis with respect to each input

In this section we study the effect of input variables on the stability of the networks. More specifically, we study the effect of input variations on the stability of the networks by quantifying ϑ Ω ϵ,K m with ϵ ∈]0, 1], for various choices of K, instead of a global Lipschitz constant accounting for the influence of the whole set of inputs. The obtained value of ϑ Ω ϵ,K m allows us to evaluate how the corresponding group of variables may potentially affect the stability of the network. By performing this analysis for several choices of K, we thus generate a multivariate analysis of the Lipschitz regularity of the network.

As shown by Proposition 6, varying the ϵ parameter is also insightful since it allows us to measure how the network behaves when input perturbations are gradually more concentrated on a given subset of inputs.

Although our approach can be applied to groups of inputs, for simplicity in this section, we will focus on the case when the sets K reduce to singletons.

For each dataset, we first perform a standard training when designing the network. To facilitate comparisons, the Lipschitz star of the network trained in such standard manner is presented as the first subplot of all the figures. Next, we show the variation in terms of input sensitivity, when i) a Lipschitz target is imposed, and ii) when an adversarial training of the networks is performed. The network architecture remains unchanged, for all our experiments and each dataset. All the Lipschitz constants for each value of ϵ are calculated using LipSDP-Neuron [START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF]. Since an increased stability may come at the price of a loss of accuracy [START_REF] Tsipras | Robustness may be at odds with accuracy[END_REF], we also report the performance of the networks on test datasets in terms of MAE (Mean Absolute Error) for each of the Lipschitz star plot. Comments on the Thales Air Mobility industrial application: From the star plots, it is clear that the various variables have a quite different effect on the Lipschitz behavior of the network. This is an expected outcome since these variables carry a different amount of information captured by learning. From Figure 4.5 we observe that variables 1 -Flight Distance and 3 -Initial ETE play a prominent role, while variables 5 -Longitude Origin, and 8 -Longitude Destination are also sensitive. Some plausible explanations for these facts are mentioned below.

• Flight distance: The impact of a change of this input can be significant since, as a result of air traffic management separation rules, commercial aircrafts cannot freely increase their speed to minimize the impact of a longer flight distance.

• Initial ETE: Modifying this input is equivalent to changing the initial conditions, which will have a significant impact. It is possible, in the worst case scenario, to accumulate other perturbations coming from other cou- pled inputs and parameters (e.g., weather conditions) and this is probably the reason why the partial Lipschitz constant is very high, and close to the global Lipschitz constant.

• Longitude origin and destination parameters: These parameters are related to different continents and even countries of the origin and destination airports, probably with different qualities of air traffic equipment.

. Effect of adversarial training

Generating adversarial attacks and performing adversarial training constitute popular methods in designing robust neural networks. However, these techniques have received less attention for regression tasks, since most of the works deal with classification tasks [START_REF] Eykholt | Robust physical-world attacks on deep learning visual classification[END_REF][START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF][START_REF] Kurakin | Adversarial examples in the physical world[END_REF]. Also, most of the existing works in the deep learning literature are for standard signal/image processing problems, whereas there are only few works handling tabular data [START_REF] Ke | TabNN: A universal neural network solution for tabular data[END_REF][START_REF] Zhang | Deep learning over multi-field categorical data[END_REF]. One noticeable exception is [START_REF] Ballet | Imperceptible adversarial attacks on tabular data[END_REF] which investigates problems related to adversarial attacks for classification tasks involving tabular data.

Since our applications are related to regression problems for which few existing works are directly applicable, we designed a specific adversarial training method. More specifically, for a given amplitude of the adversarial noise and for each sample in the training set, we generate the worst attack based on the spectral properties of the Jacobian of the network, computed by backpropagation at this point. At each epoch of the adversarial training procedure, we solve the underlying minmax problem [START_REF] Tu | Theoretical analysis of adversarial learning: A minimax approach[END_REF]. More details on the generation of adversarial attacks for regression attacks can be found in [START_REF] Gupta | An adversarial attacker for neural networks in regression problems[END_REF].

The generated adversarial attacks from the trained model at the previous epoch are successively concatenated to the training set for the next training epoch, much like in standard adversarial training practices using FGSM [START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] and Deepfool [START_REF] Moosavi-Dezfooli | Deepfool: a simple and accurate method to fool deep neural networks[END_REF] attacks. While generating adversarial attacks on tabular data, some of the variables may be more susceptible to attacks than others. The authors of [START_REF] Ballet | Imperceptible adversarial attacks on tabular data[END_REF] take care of this aspect by using a feature importance vector. They also only attack the continuous variables, disregarding categorical ones while generating attacks. For the Power plant and Boston Housing datasets, we attack all the four input variables, while on the MPG dataset, we attack only the continuous variables. For the industrial dataset, we generate attacks for the five most sensitive input variables. We also tried attacking all the variables of the dataset but this was not observed to be more efficient. The results in form of Lipschitz star are given in Figures 4.6-4.9.

As expected, adversarial training leads to a shrinkage of the star plots, which indicates a better control on the stability of the trained models, while also improving slightly the MAE. In the test we did, we observe however that our adversarial training procedure is globally less efficient than the spectral normalization technique.

. Sensitivity w.r.t. pair of variables

We now consider the case when the set K contains pairs of elements. We first show the corresponding Partial Lipschitz constants using a Lipschitz star, for the different datasets we have discussed in the article. Vertices in the Lipschitz star represent the obtained Lipschitz constant value ϑ Ω ϵ,K m for all possible combinations of pair of variables with varying values of ϵ, i.e., it represents the sensitivity w.r.t. to that particular pair.

As shown by Figure 4.10-4.13, this Lipschitz star representation can be useful for displaying the influence of groups of variables instead of single ones. This may be of high interest when the number of inputs is large, especially if they can be grouped into variables belonging to a given class having a specific physical meaning (e.g., electrical variables versus mechanical ones). Such Lipschitz star representation might however not be very insightful for identifying the coupling that may exist between the variables within a given group. For example, it may happen that, considered together, two variables yield an increased sensitivity than the sensitivity of each of them individually. The reason why we need to find a better way for highlighting these coupling effects is related to Proposition 6v) which states that, for every ϵ ∈]0, 1] and (k, ℓ) ∈ {1, . . . ,

N 0 } 2 , max{ϑ Ω ϵ,{k} m , ϑ Ω ϵ,{ℓ} m } ⩽ ϑ Ω ϵ,{k,ℓ} m . (4.44)
This property means that, when considering a pair of inputs, the one with the highest partial Lipschitz constant will "dominate" the other. To circumvent this difficulty and make our analysis more interpretable, we can think of normalizing the Lipschitz constant in a suitable manner. Such a strategy is a common practice in statistics when, for example, the covariance of a pair of variables is normalized by the product of their standard deviations to define their correlation factor. Once again, we can take advantage of the properties established in Proposition 6 to provide us a guideline to perform this normalization. In addition to Eq. (4.44), according to Property viii),

ϑ Ω ϵ,{k,ℓ} m ⩽ (ϑ Ω ϵ,{k} m ) p * + (ϑ Ω ϵ,{ℓ} m ) p * 1/p * + o(ϵ) ⩽ 2 1/p * max{ϑ Ω ϵ,{k} m , ϑ Ω ϵ,{ℓ} m } + o(ϵ). (4.45)
The two previous inequalities suggest to normalize the Lipschitz constant for pairs of inputs by defining

ϑ Ω ϵ,{k,ℓ} m = 1 2 1/p * -1 ϑ Ω ϵ,{k,ℓ} m max{ϑ Ω ϵ,{k} m , ϑ Ω ϵ,{ℓ} m } -1 . (4.46)
Indeed, when ϵ is close to zero, Eq. (4.44)-Eq. (4.46) show that ϑ

Ω ϵ,{k,ℓ} m ∈ [0, 1].
Note that, for the diagonal terms, ϑ Ω ϵ,{k,k} m = 0. The higher ϑ Ω ϵ,{k,ℓ} m , the higher the gain in sensitivity due to the coupling between k and ℓ. The normalized values for the different datasets are reported in Table 4.9.

. Interpretation of the results

We summarize some observations/properties concerning the stability of the NNs which can be drawn from training on different datasets and leveraging the quantitative tools we have proposed in this article.

(a) Combined Power Plant Dataset

• '3 -Exhaust Vacuum' is the most sensitive variable out of the four variables.

• We observe for any variable coupled with '3' gives a higher partial Lipschitz constant.

• From Table 4.9(a), we see that the effect is mostly caused by the sensitivity of '3' and there is no gain when coupled with other variables. Hence, '3' dominates the overall sensitivity of the NN. • On the other hand, we observe that, '0' when coupled with '1' and '2' becomes more sensitive as evidenced by the gain in Table 4.9(a).

(b) Auto MPG Dataset

• Variable '6 -Origin' and '3 -Weight' are the most sensitive variables.

• The values of partial Lipschitz constant peak when the other variables are coupled with '3' or '6'.

• From Table 4.9(b), we see that most of the values coupled with either '3' or '6' are close to zero, except when '3' and '6' are coupled together. Also, we see an exception when '5' is coupled with either '3' or '6'. This suggests that altogether '3', '5', and '6' have a higher impact on the stability of the network.

(c) Boston Housing Dataset

• Variable '7 -DIS' and '11 -LSTAT' are the most sensitive variables.

• We observe a high partial Lipschitz constant when coupling any variables with '7' or '11'. • From Table 4.9(c), we see that all the values for both '7' and '11' coupled with other variables are close to zero, except when '7' and '11' are jointly considered. Hence, '7' and '11' dominate the sensitivity of the NN.

• We observe from the table of normalized values, that '2-9' have a higher impact on the sensitivity of the NN when coupled. Similar observation can be made for pairs '2-8','1-4','3-4'.

(d) Thales Air Mobility industrial application

• Variable '1 -Flight distance','3 -Initial ETE', and '8 -Longitude Destination' are the most sensitive variables.

• We see peaks in the partial Lipschitz constant values when these highly sensitive variables are coupled with other variables.

• But when analyzing the normalized tables, it becomes clear that the gain is mostly due to these sensitive variables. • We also observe from Table 4.9(d), an increased sensitivity of '0' when coupled with other variables '5', '7', '10', '11', and '13'.

0 -1 0 -2 0 -3 0 -4 0 -5 0 -6 0 -7 0 -8 0 -9 0 -1 0 0 -1 1 1 -2 1 -3 1 -4 1 -5 1 -6 1 -7 1 -8 1 -9 1 -1 0 1 -1 1 2 -3 2 -4 2 -5 2 -6 2 -7 2 -8 2 -9 2 -1 0 2 -1 1 3 -4 3 -5 3 -6 3 -7 3 -8 3 -9 3 -1 0 3 -1 1 4 -5 4 -6 4 -7 4 -8 4 -9 4 -1 0 4 -1 1 5 -6 5 -7 5 -8 5 -9 5 -1 0 5 -1 1 6 -7 6 -8 6 -9 6 -1 0 6 -1 1 7 -8 7 -9 7 -1 0 7 -1 1 8 -9 8 -1 0 8 -1 1 9 -1 0 9 -1 1 1 0 -1 1

. Summary

The contributions of this work are summarised as follows:

• This chapter proposes a new multivariate analysis of the Lipschitz regularity of a neural network. The theoretical foundations of our approach are given in Section 4.3.

• A multivariate analysis of the Lipschitz properties of NNs is performed by generating a set of partial Lipschitz constants. This opens a new dimension to studying the stability of NNs.

• Our sensitivity analysis allows us to capture the behaviour of an individual input or group of inputs on the output of the neural network. • The results of this analysis are displayed by a new graphical representation termed as a Lipschitz star. This representation is helpful in displaying how each input or group of inputs contributes to the global Lipschitz behaviour of a network.

• Using the proposed analysis, we also study quantitatively the effect of spectral normalization constraint and adversarial training on the stability of NNs.

• We showcase our results on various open-source datasets along with a real industrial application in the domain of Air Traffic Management. Chapter 5 Spectral Normalization Loop for controlling stability

. Introduction

As discussed in the previous chapter, a form 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. With this chapter we move to controlling and tightening the Lipschitz bounds of the neural networks to promote robustness. We show the efficiency of a formal method based on the Lipschitz constant for quantifying the stability of neural network models. We present how this formal 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 (accuracystability trade-off). We first, show the efficacy of the proposed control loop on a public tabular dataset related to German Credit risk and then move to safety critical application. We evaluate our approach on an application related to Small Unmanned Aerial Vehicles (sUAVs) 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 a novel application of neural networks to detect in real-time elevon positioning faults to allow the remote pilot to take necessary actions to ensure the safety (e.g.remote triggering of the drone emergency parachute). In both scenarios, we show the effectiveness of our proposed algorithm in reaching a optimal model following both accuracy and stability targets that are user specified.

. Limitations of Previous Works

As introduced in the Section 2.3.2, the Lipschitz constant can be used to train more robust models. The existing methods are hindered with several limitations. One being the lack of accuracy-vs-stability trade-off. Randomly constraining neural networks to have very tight Lipschitz bounds or maybe making them 1-Lipschitz as in the case of Spectral Normalization will lead the accuracy to drop drastically on clean data. It is possible that, under very tight constraints, the networks fails to learn anything. Such a reduction in the performance is not acceptable in applications such as UAVs. Engineers need to maintain respectable levels of accuracy and stability to deploy any solution in that case. To mitigate this effect, we propose a clever use of Spectral Normalization to design a stability control loop allowing us to lower the Lipschitz bounds while maintaining a good performance.

. Proposed Design method

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. 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 Eq. (2.18), 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 makes the minimization of the loss function hard. 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.

Consider an m-layer fully connected NN architecture (FCN) denoted by T with a desired Lipschitz constant θ target and an accuracy target acc target . We can constrain the spectral norm of each layer to be less than m θ target . According to Eq. (2.7), this ensures that the upper bound on the global Lipschitz constant is less than θ target . This bound is however over-pessimistic, 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 θ target .

The approach is summarized in Algorithm 1. Module Train SN (θ, . . .) is a spectral normalization 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 √ θ), 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 (classification accuracy) of the trained model on the validation dataset. If θ is small, the training will be over-constrained and the accuracy of the model will suffer. To relax this constraint, we iteratively increase the value of θ through a factor α close to 1 to find the best value θ opt which would be less than θ 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. the accuracy of the trained model is greater than acc target and the Lipschitz constant is less than θ target .

Algorithm 1 Spectral Normalization-based Stability Control Loop

. Experiments

. Open source dataset

We first show the efficacy of the proposed stability control loop on an Open source dataset. We use Statlog (German Credit Data) Data Set. It has a classification setting of 1,000 instances with 23 continuous features per instance. The goal is to determine credit-worthiness of customers (Good or bad) . We build fully connected networks with ReLU activation function on all the hidden layers, except the last one. The architecture has two hidden layers with 50 and 20 neurons. The input attributes are scaled between 0 and 1.

If we train the neural network model without any stability constraint, i.e., perform standard training, the network performance is 75.2% with an estimate for the Lipschitz constant of 17.58. Table 5.1 tabulates the results for basic spectral normalization (with a Lipschitz target) in columns 2 and 3. The value of θ is determined using LipSDP-neuron. The results with our proposed method, adding an accuracy target, are given in columns 5 and 6. For the first set of results, we train this model without any accuracy target and perform a simple Spectral normalization to constrain the Lipschitz norm of each layer to 1; we see a reduction in the accuracy by more than 5% with respect to the standard trained model. The accuracy improves slowly as the Lipschitz constraint is relaxed, but the values are much lower than the target value. For the second set of results, we provide both accuracy and Lipschitz targets so as to obtain a better solution by using Algorithm 1. We see an increase in the overall performance of the models. We see a similar increase of the accuracy as the stability constraints are relaxed and it can even become better than the baseline in some cases. We also observe that the upper bounds for the Lipschitz constant provided by LipSDP are much closer to the targeted values. For an accuracy target of 75% with a Lipschitz target less than 5, the algorithm fails to find a trade-off between stability and accuracy. We report the results for θ target less than 5 with acc target of 70% in Table 5.1. The best solution in terms of accuracy is obtained with θ target = 12 and acc target = 75%.

. UAV use case

Unmanned Aerial Vehicles (UAVs) are currently being used in various applications such as surveillance, home delivery of products by e-commerce websites, disaster management, emergency services, efficient farming, insurance underwriting, etc. UAVs are capable of replacing humans in potentially dangerous situations, but also reduce the cost of operations in some businesses, these being the main advantages of using drones. But systems of this kind need to be continuously monitored to detect in real time any problem (called a "fault") that could lead to a failure of the system, to the potential loss of the UAV, and thus an interruption of the mission. Broadly, a fault analysis algorithm comprises of a continuously monitoring a system during its operation for detecting faults (Fault Detection), locate faults (Fault Isolation) and predict their temporal evolution (Fault Identification). This projected surge of use in the future by all industries makes real-time fault diagnosis of these systems a priority to avoid further critical losses i.e there can be many mission or safety critical implication to the faults if left undetected. Fault diagnosis techniques are becoming increasingly important in order to ensure high levels of safety and reliability in aeronautical systems. This gives rise for the need for more intelligent and reliable systems in the future. Hence, the advances in UAV or aerial vehicle technology is accompanied by the concern of safety. Designing a system that is fault tolerant is a recognized approach in the field of dependability engineering. Formally, fault tolerance is the ability of a system to perform its intended function despite the presence or occurrence of faults, whether physical damage to the hardware, software defects [START_REF] Laprie | Definition and analysis of hardware-and software-fault-tolerant architectures[END_REF], or malicious attacks. In this work, we will only deal with the fault detection capability that is depicted in Figure 1. It is designed to detect mechanical faults like an elevon stucked or locked in a given position. Other dimensions of fault tolerance involving error recovery mechanisms or other safety measures are not addressed here.

. Limitation of AI in UAV use cases

Machine Learning 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 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). 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] X.-H. 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 data-driven 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.

. 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 = Fu com + E (5.1)
where u com is the commanded control deflection from the autopilot to the actuators, 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

u app r u app l = f r 0 0 f l u com r u com l + e r e l , (5.2)
where () r and and () l stand 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 com r , 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. Each fault category is labeled according to its fault code in order to generate the .

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, for . LipSDP, we do not use the tightest bound LipSDP-Network, since it is erroneous as mentioned in Section 2.2.1.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. The 12 th and 13 th July flights contain approximately 15 minutes of similar fault pattern of reduced efficiency of right elevon ( f r = 0.3) making it suitable for binary classification. 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 we employed sigmoid for binary classification. We use the cross entropy function as the loss function for training the model. The feature set is standardized by removing the mean and scaling to unit variance. in terms of Lipschitz constant for a varying number of hidden layers and hidden neurons is shown in Table 5.4. Notation [START_REF] Kurakin | Adversarial examples in the physical world[END_REF][START_REF] Singh | Trustworthy AI[END_REF] in Table 5.4 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. The 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 5.3 with Algorithm 1. We chose the same configuration for the neural networks as the baseline and introduced constraints by setting θ target = 4 and acc target = 75%. Table 5.5a 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. 1 It is useful to test the 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 5.6. 2 We remark that spectral normalized 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 elevons. For injecting the faults, efficiency of the right elevon is reduced ( f r ) and similarly for the left control elevon ( f l )Fault codes for different classes are : '0' -nominal , '1'f r = 0.3 , '2' -f l = 0.9 and '3'f l = 0.3, defining 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. For training the neural network, flight data are selected between 400 and 2600 seconds. 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. It is to be noted here again that the flights made on different days will be affected by the environmental factors, which will impact the real-time prediction of the faults. 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 training time, and also on the flight data of 23 rd July. We show the prediction and generalization capabilities of the neural networks over stateof-the-art SVM and logistic regression. We tabulate the classification metrics Chapter 6 Conclusion and Future Work

. Conclusion

Although they may appear at the top of their advancements, neural networks are hindered by security, privacy and safety issues due to their sensitivity to various attacks and perturbations they might encounter during their operation. It thus becomes essential to understand the causes of neural network instability, identify the concern areas, and provide solutions aiming to improve their stability, and ensure that AI-based systems work as desired during their whole operational phase.

In this thesis, we have discussed various limitations and challenging problems in the emerging field of robust AI, especially related to safety critical applications. We have emphasized that the Lipschitz regularity of the network is a useful concept to address stability issues from a mathematically rigorous perspective. In a nutshell, we have studied various concepts and latest developments towards formalizing the robustness and Lipschitz stability of neural network models. As shown in Chapter 2, the concepts of robustness is quite broad and varied. For example, according to experts in safety at Thales, we have to distinguish as:

• ML Model Robustness: The capacity of an ML model to preserve its expected/intended performance under well-characterized abnormalities (e.g. lack of generalization) or deviations to its inputs and operating conditions outside its Operational Design Domain (ODD).

• ML Model Stability: The capacity of an ML model to preserve its expected/intended output(s) under well-characterized and bounded perturbations to its inputs and operating conditions within its operational design domain .

The perturbations encountered can be natural (e.g. sensor noise, bias,...), due to failures (e.g. invalid data from degraded sensors) or maliciously inserted (e.g. pixels modified in an image) to fool the model predictions. Perturbations can also be simply defined as true data locally different from the original data used for the model training and that might lead to a wrong prediction and an incorrect behaviour of the system.

Below are the summarised contributions of this thesis before discussing potential future work to promote the ideas presented in this thesis.

• We discussed a major gap in the field of adversarial machine learning and regression tasks. Regression tasks cover a wide range of safety-and missioncritical applications, and significantly fewer efforts have been put into understanding the vulnerabilities in regression models. To overcome this limitation, we have proposed a generic white-box attacker for regression neural networks. The Jacobian of the trained neural network is used to drive the attacks. Analogous to classification scenarios, useful metrics are proposed to quantify the effectiveness of these attacks.

• Our proposed attacker has several unique features such as attacking few or group of inputs and can address different norm perturbations on input and output. Our proposed attacker is shown to be more powerful than random attackers.

• Moving towards the notion of Lipschitz certificates which is the backbone of the stability analysis carried out in this thesis, we propose a new multivariate analysis of the Lipschitz regularity of a neural network by defining partial Lipschitz constants. The theoretical foundations of this approach are presented in section 4.3. Most current works utilizing Lipschitz certificates offer a single parameter analysis and limits the information that can be extracted from existing theoretical guarantees of stability.

• Our sensitivity analysis allows us to capture the influence of an individual input or group of inputs on the output of the neural network. The results of this analysis are displayed by a novel graphical representation termed as a Lipschitz star. This representation is helpful in displaying how each input or group of inputs contributes to the global Lipschitz behaviour of a network.

• For controlling the Lipschitz bound of the neural network, we present a stability control loop based on spectral normalization to train optimal neural models achieving performance and stability targets. Model trained with tight Lipschitz certificates are more robust to adversarial attacks and also maintain good accuracy.

• All of the algorithms and methods presented in this thesis are tested on various open-source tabular datasets, an air traffic management case provided by Thales LAS France and a real dataset of UAV safety application which is publicly available.

• Following the insights presented in the work, we have shown that it is possible to build neural networks which are "Stable-by-Design" and with guaranteed stability through Lipschitz certificates. These results regarding robustness are strong developments towards understanding the vulnerabilities in ML models and provide defenses against attacks.

. Future work

The contributions of the thesis being summarised, a few possible future research directions based on the contributions and new concepts related to the work are mentioned in this section.

. Stronger Adversarial Attacks for Regression

The adversarial attacker presented in the thesis is effective, however stronger attacks can be achieved aligning to the classification scenarios. For generating adversarial attacks in chapter 3, we use fooling error to make the perturbed sample deviate from the ground truth sample. We can use stronger or multiple objectives to create stronger attacks. More sophisticated optimization algorithms such as augmented Lagrangian principles [START_REF] Rony | Augmented Lagrangian adversarial attacks[END_REF] and proximal splitting techniques [START_REF] Rony | Proximal splitting adversarial attacks for semantic segmentation[END_REF] could be explored for this purpose in the context of regression tasks. Our tests concentrated on fully connected networks; it is worth pointing out that the proposed adversarial attacker can be applied to any network architecture, but it needs to be validated on other architectures. Also, there is a growing body of work on intermediate-level attacks (ILAP) [START_REF] Huang | Enhancing adversarial example transferability with an intermediate level attack[END_REF][START_REF] Li | Yet another intermediate-level attack[END_REF] which try to eliminate or reduce the need for soft or hard labels. An initial direction -in the latent space of the hidden, intermediate layers -is determined via some suitable attack. Gradient steps are then used to maximally perturb the intermediate hidden representation to find a suitable adversarial attack. Interestingly, such a method can often outperform the method used to set the initial direction itself, while being more transferable across different models. These work eliminating the use of labels can help to design adequate adversarial attackers for regression tasks. Overall, the literature studying adversarial attacks and defenses for regression tasks is limited and needs more attention in future works.

. MIMO and Volterra Networks

One obvious extension of our analysis presented in chapter 4, where we have concentrated mostly on MISO (Multi-Input Single-Output) systems, is the sensitivity analysis of MIMO ((Multi-Input Multi-Output) systems. In chapter 4, we provide partial Lipschitz constants w.r.t to input when all the outputs are active-ignoring how the individual or group inputs can affect individual or a group of outputs. To get further analysis, we can easily calculate these partial Lipschitz constants corresponding to only specific outputs. Such analysis is important for understanding the sensitivity of inputs, but it can also help to control unintended couplings between variables that can be introduced by the neural network model. We also observe a cross-link between the weighting matrices allowing us to obtain partial Lipschitz constants and explainability masks used in images, highlighting the regions of the image which contribute most to the decision of CNNs to predict a certain label. Similarly, the partial Lipschitz constant highlights the importance of certain inputs w.r.t to output.

Another possible direction is to design new neural network architectures for tabular data learning using the concept of Volterra neural networks [START_REF] Ghanem | Latent code-based fusion: A volterra neural network approach[END_REF] with reduced or no use of activation functions. A neural network with a lower number of activation functions should be intuitively more robust with less regions where the approximation function exhibits steep slopes (i.e. higher local Lipschitz bounds), hence it is more stable to perturbations. We can design a neural network using second-order Volterra kernels to learn all relationships as in Eq. (4.36).

. Spectral Normalization and its variants

In chapter 5, we use the same level of Lipschitz normalization at each layer in the network. This might not be best way to normalize the Lipschitz constant of the layer because different layers in the neural network architecture learn more or less informative features, which implies that Lipschitz normalization should depend on the importance of the layer or its width. To mitigate this, we can add a learnable parameter, associated with each individual layer, to the spectral normalization algorithm which can adjust itself while training to attain an optimal robust model. Also, we have used fully connected layers in our experiments in chapter 5 but Spectral normalization is well suited and tested for convolutional neural networks as well. We are limited by the use of LipSDP module or any other module to estimate Lipschitz bounds that work well only for dense networks with non-expansive operators. Obtaining tight and numerically tractable Lipschitz bounds for convolutional networks thus remains a main challenge. In [START_REF] Gouk | Regularisation of neural networks by enforcing Lipschitz continuity[END_REF], author proposes a relaxed approach to normalize weights of different layers which is not as restrictive as Spectral normalization. They do so for improving the performance of the neural networks. Such a relaxation of normalization can also be explored to have better control on the stability. Also, another means of attaining an accuracy-stability trade-off can be achieved by utilizing two objective functions in order to maintain accuracy (i.e. first achieve the best possible accuracy and keep the accuracy in the tolerance levels) and gradually changing the weights to achieve tighter Lipschitz bounds.

. Other architectures and application domains

Our work concentrated on using simple fully connected layers, understanding them, and controlling their stability. More prevalent architectures such as ResNet, EfficinetNet, graph networks, transformers, and generative models such as GANs need to be explored. Even in case of tabular data, though researchers may rely on methods such as XgBoost [START_REF] Chen | Xgboost: A scalable tree boosting system[END_REF], deeper learning models, especially attention-based mechanisms have emerged, such as TabTransformer [START_REF] Huang | Tabtransformer: Tabular data modeling using contextual embeddings[END_REF], FT-Transformer [START_REF] Gorishniy | Revisiting deep learning models for tabular data[END_REF], TabNet transformer [START_REF] Arik | Tabnet: Attentive interpretable tabular learning[END_REF], and most recently Hopfield networks [START_REF] Schäfl | Hopular: Modern Hopfield networks for tabular data[END_REF], successful in learning such data. Also, a point to be noticed for these architectures is the use of activation functions such as GELU in architectures, which are not firmly non-expansive, while the architecture in itself consists of simple fully connected layers. In addition, transformers using standard DP product attention mechanism were shown to be non Lipschitz continuous [START_REF] Kim | The Lipschitz constant of selfattention[END_REF], again questioning the use of such networks in real-world applications. A Lipschitz continuous attention mechanism was proposed, but the robustness of this block is still unknown. Another area which needs to be studied is the use of batch and layer normalization layers and how they contribute to the overall robustness of the neural networks. Similarly, [START_REF] Dasoulas | Lipschitz normalization for selfattention layers with application to graph neural networks[END_REF] proposes Lipschitz normalization for improving the performance of attention based graph neural networks without discussion of their adversarial and certified robustness. These architectures must be studied for robustness and vulnerabilities to deliver more trustworthy models. Overall, there have been plethora of works designing and defending such networks against adversarial attacks, but handful that certify theoretical guarantees against perturbations.

. Other notions of robustness

A very close notion to Lipschitz bounds is related to curvature, which is a mathematical quantity that encodes the deviation of a curve from 'flatness'. The high degree of non-linearity present in deep neural networks, or equivalently their large curvature, is critical to achieving good performance in difficult tasks. Networks learnt by penalizing the curvature keep the trained model as linear as possible. This implies that networks are learnt to obtain certain upper bounds on the curvature similar to Lipschitz bounds. When a function is twice differentiable, the maximum of the norm of the Hessian is the local Lipschitz constant of the Jacobian of the network. In deep learning, the curvature of a function at a point is quantified as the norm of the Hessian at that point. So, it is equivalent to look at the Lipschitz properties of the Jacobian. Low curvature bounds is yet another way to verify that the linearity is indeed the source of increased robustness. In literature, a low curvature is imposed to induce robustness in the networks, and this is done by penalizing the norms of the Hessian. Moosavi et al. [START_REF] Moosavi-Dezfooli | Robustness via curvature regularization, and vice versa[END_REF] penalizes the Frobenius norm of the Hessian, and [START_REF] Qin | Adversarial robustness through local linearization[END_REF] introduces a local linearity regularizer, which also implicitly penalizes the Hessian. Relationships and analysis can be explored between Lipschitz bounds and curvature since the two concepts are complementary and closely related. It would thus be interesting to explore training strategies utilizing curvature as the objective.

A.1 . Effect of Regularization Techniques

We present the results for three regularisation techniques: ℓ 1 , ℓ 2 , and Dropout [START_REF] Srivastava | Dropout: a simple way to prevent neural networks from overfitting[END_REF]. We study the relationship between the parameters associated with each regularisation and the NN stability quantified by its Lipschitz constant. Attention also to be paid to the resulting accuracy. We apply these regularization techniques while training our NNs, then compute a Lipschitz constant associated with the obtained weights. We use the three state-of-the-art estimation methods which have been described in Section 2.2.1.2. The first one is LipSDP [START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF] which uses Euclidean norms for both the input and output spaces (L 2,2 spectral norm). The second one is the polynomial based approach LipOpt [START_REF] Chen | Polynomial optimization for bounding Lipschitz constants of deep networks[END_REF] where the input and output spaces are equipped with the sup norm while estimating the Lipschitz constant. This estimation approach is thus linked to the L +∞,+∞ subordinate matrix norm. The third estimation method is CPLip [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF] which can work for any norm on the input and output spaces.

Each experiment was run 10 times and we chose the model with best performance (least MAE or highest classification accuracy) and computed a Lipschitz constant for this model. The results for the datasets using ℓ 1 and ℓ 2 regularization are reported in Tables A.1, A.2, and A.3 for varying values of the regularization parameter α which controls the strength of the ℓ 1 and ℓ 2 penalty on the weights. Similar results with varying drop-rates are presented in Tables A. [START_REF] Thiebes | Trustworthy artificial intelligence[END_REF], A.5, and A.6. Droprate corresponds to the proportion of the neurons which will be shut-off while training a neural network.

A.2 . Positive Weighted Networks

Next we analyse the stability of NNs when the weights are constrained to be non-negative. The comparison between arbitrary signed network and positively signed network for all the three datasets is shown in Table A.7.

A.3 . Addition of Noise to the Dataset

We also perform a stability analysis when the original dataset is corrupted with random noise which is a standard practice while training NNs having continuous input variables. From the original dataset, we generated a dataset 20 times larger by including noisy samples.

More precisely, to all the normalised input features, we added a noise value drawn from a random i.i.d. zero-mean Gaussian distribution with a small standard deviation. We combined the original training set with the generated noisy samples and trained our model on the augmented dataset. The results on Combined Cycle Power plant with variation of standard deviation are given in Table A

A.4 . Comments on the results

A first observation is that CPLip provides slightly tighter bounds than LipSDP and LipOpt. The two latter approaches may however be more scalable when applied to deeper networks. Another remark is that similar behaviours can be seen when using different norms. While designing a network for deployment, it also appears that there is a trade-off between stability and performance.

• ℓ 1 and ℓ 2 regularization increase the stability of the network consistently, but the performance is maintained up to a certain value of α, from where the accuracy drops. ℓ 1 usually yields better results, confirming the robustness of this norm as training measure.

• We globally observe an increasing stability of the neural networks as we increase the value of droprate. The results are however less consistent than with regularization.

• The positive constraint leads to a significant loss of performance in terms of MAE and classification accuracy, but the stability of the networks is improved by a significant margin.

• As expected, adding noise leads to a drop in the value of the Lipschitz constant as the noise level increases, without having a negative impact on the accuracy. b is the graphical representation of the presented concept. As it can be observed, the weights are split to emulate a MIMO system. Each kernel is associated with a specific shape of the matrix, which is depicted in Figure 1.c. The lower and upper triangular null parts are displayed in dark gray, while the light gray central part contains overlap and may use different weight values.

B.1.1 . Learning algorithm

For training the proposed ACNN (Adaptive Convolutional Neural Network), we use a stochastic gradient-like optimization based on Adam. Consider the vector of parameters of the network, η = (η i ) 1⩽i⩽m , such that, for each layer i ∈ {1, . . . , m}, η i represents a vector of dimension N i (N i-1 + 1), composed of the elements of weight matrix W i and the components of the bias vector b i .

To secure the conditions of robustness while imposing the desired structure for our network, the parameter vector η is projected onto a closed set S that expresses all these constraints. The parameter update at epoch n > 0 is performed for mini- batches (M q,n ) 1⩽q⩽Q . If the training data are denoted by (z k ) 1⩽k⩽K , where z k is the k-th pair of inputs and their associated outputs, the operations performed during the n-th epoch reads: Partition {1, . . . , K} into mini-batches (M q,n ) 1⩽q⩽Q (B.1)

For every q ∈ {1, . . . , Q} (B.2) where the square, the square root, and the division are performed component-wise, and

               t = (n -1)Q + q For every i ∈ {1, . . . , m}           g i,t = ∑ k∈M q,n ∇ i ℓ z k , (η i,t 1⩽i⩽m ) µ i,t = β 1 µ i,t-1 + (1 -β 1 )g i,t ν i,t = β 2 ν i,t-1 + (1 -β 2 )g 2 i,t γ t = γ 1 -β t 2 /
S i,t = η i | [(η ⊤ j,t+1 ) j<i η ⊤ i (η ⊤ j,t ) j>i ] ⊤ ∈ S . (B.4)
Here-above, ℓ denotes the loss function, ∇ i represents the gradients with respect to η i ; µ i,t and ν i,t represent the first and second momentum estimates at the iteration t, initialized with µ i,0 = ν i,0 = 0. P S i,t designates the projection onto the constraint set S i,t . Although the set S is non-convex, the sets S i,t will be defined as the intersection of three closed and convex constraint sets, as detailed next. The parameters used for learning are set to β 1 = 0.9, β 2 = 0.999, γ = 0.001, and ϵ = 10 -12 .

To ensure the computation of a tight robustness bound we impose non-negative weights for every i ∈ {1, . . . , m} by considering the constraint set:

D i = {W i ∈ R N i ×N i-1 | W i ⩾ 0} (B.5)
Let R i (resp. Q i ) be the number of output (resp. input) channels used in layer i ∈ {1, . . . , m}. The proposed algebraic structure of each weight operator is guaranteed by splitting the corresponding matrix

W i into R i × Q i sub-matrices of dimension N ′ i × M ′ i (with N ′ i = N i /R i and M ′ i = N i-1 /Q i )
, denoted by (W (r,q) i ) 1⩽r⩽R i ,1⩽q⩽Q i . The desired band structure of the sub-matrix W (r,q) i is ensured by imposing that it belongs to the following vector space:

E i = {(V u,v ) 1⩽u⩽N ′ i ,1⩽v⩽M ′ i ∈ R N ′ i ×M ′ i | ∀(u, v) s.t. |(M ′ i -1)(u -1) -(N ′ i -1)(v -1)| ⩾ d i , V u,v = 0}.
Herein, d i is an integer and, if

N ′ i = M ′ i > 1, 2⌊d i /(N ′
i -1)⌋ -1 plays a role similar to a kernel length in standard CNNs. Finally, to control the robustness, we need to limit the Lipschitz constant of

B.2.2 . Experimental setup

The noisy data for training, validating, and testing is generated by adding white Gaussian noise to the original samples. The noise has zero mean and its standard deviation is randomly chosen so that the resulting signal-to-noise ratio (SNR) varies between 5 and 30dB. The dataset samples are normalized between 0 and 1. We extract the frequency features from the audio signal using a STFT. The network estimates the STFT coefficients of the samples and an ISTFT is performed as the post-processing step. We consider a Hanning sliding analysis window of length T = 23 ms, with an overlap between two consecutive windows of 50%. The STFT is performed on 1024 points. In total, from each audio segment, a vector of length L = 513 frequency coefficients is obtained, constituting the input of our ACNN.

The denoising is performed using a 6-layer ACNN architecture, as presented in Figure B.1. The network has an encoder-decoder structure. Each layer employs ReLU as activation function followed by a batch normalization step that acts as a regularizer and prevents the model from overfitting.

B.2.3 . Simulations and results

In order to measure the performance of our proposed ACNN architecture, we perform two sets of experiments. In the first set, we control the Lipschitz constant

2. 1

 1 m-layered feed-forward neural network architecture. For the i th -layer, W i is the linear weight operator, b i the bias vector, and R i the activation operator. . . . . . . . . . . . . . . . . 2.2 A network T is split into three sub-networks U 1 , U 2 , and U 3 . . . . . . . . . . . . . . . . . 3.1 Network Architecture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Error distribution of random attacks and proposed adversarial attack on Combined cycle Power Plant dataset for perturbation level of 2 × 10 -1 for ℓ 2 attack. . . . . . . . . . . . . 3.3 Error distribution of random attacks and proposed adversarial attack on Red-wine dataset for perturbation level of 2 × 10 -1 for ℓ 2 attack. . . . . . . . . . . . . . . . . . . . . . . . 3.4 Error distribution of random attacks and proposed adversarial attack on Abalone dataset for perturbation level of 1 × 10 -1 for ℓ 2 attack. . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Error distribution of random attacks and proposed adversarial attack on industrial dataset for perturbation level of 2 × 10 -1 for ℓ 2 attack. . . . . . . . . . . . . . . . . . . . . . . . 4.1 Lipschitz Star: Representing Partial Lipschitz constants w.r.t. to the inputs as the vertices of the radial plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Sensitivity w.r.t. to each input on Combined Cycle Power Plant dataset. Influence of a spectral normalization constraint. a) Standard training: Lipschitz constant = 0.66, MAE = 0.007 , b) With spectral normalization: Lipschitz constant = 0.25, MAE = 0.0066. . . . . 4.3 Sensitivity w.r.t. to each input on Auto MPG dataset. Influence of a spectral normalization constraint. a) Standard training : Lipschitz constant = 2.75, MAE = 0.05 , b) With spectral normalization: Lipschitz constant = 0.76, MAE = 0.04. . . . . . . . . . . . . . . . . . . 4.4 Sensitivity w.r.t. to each input on Boston Housing dataset. Influence of a spectral normalization constraint. a) Standard training : Lipschitz constant = 18.56, MAE(y 1 ) = 2.45, MAE(y 2 ) = 1.41, b) With spectral normalization: Lipschitz constant = 8.06 , MAE(y 1 ) = 2.96, MAE(y 2 ) = 1.35. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Sensitivity w.r.t. to each input on Thales Air Mobility industrial application. Influence of a spectral normalization constraint. a) Standard training: Lipschitz constant = 45.46, MAE = 496.37 (s), b) With spectral normalization constraint: Lipschitz constant = 16.62, MAE = 478.88 (s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Sensitivity w.r.t. to each input on Combined Cycle Power Plant dataset. Effect of adversarial training. a) Standard training: Lipschitz constant = 0.657, MAE = 0.007, b) Adversarial training: Lipschitz constant = 0.37, MAE = 0.0068. . . . . . . . . . . . . . . . . . . . . . 4.7 Sensitivity w.r.t. to each input on Auto MPG dataset. Effect of adversarial training. a) Standard training: Lipschitz constant = 2.75, MAE = 0.05, b) Adversarial training: Lipschitz constant = 1.84, MAE = 0.042. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Sensitivity w.r.t. to each input on Boston Housing dataset. Effect of adversarial training. a) Standard training : Lipschitz = 18.56, MAE(y 1 ) = 2.45, MAE(y 2 ) = 1.41, b) Adversarial training: Lipschitz constant = 16.50, MAE(y 1 ) = 2.35 MAE(y 2 ) = 1.32. . . . . . . . . . . 4.9 Sensitivity w.r.t. to each input on Thales Air Mobility industrial application. Effect of adversarial training. a) Standard training: Lipschitz = 45.47, MAE = 496.37 (s), Adversarial training. b) Lipschitz = 34.26, MAE = 494.7 (s). . . . . . . . . . . . . . . . . . . . . . . 4.10 Sensitivity w.r.t to pair of variables on Combined Power Plant dataset . . . . . . . . . . . 4.11 Sensitivity w.r.t to pair of variables on Auto MPG dataset. . . . . . . . . . . . . . . . . . 4.12 Sensitivity w.r.t to pair of variables on Boston Housing Dataset. . . . . . . . . . . . . . . 4.13 Sensitivity w.r.t to pair of variables on Thales Air Mobility industrial application. . . . . . . 5.1 Fault Detection in drones/UAVs : Inertial sensors of the drone: accelerometer measures translational acceleration, and the gyrometer measures angular velocities (rotational motion). By combining these measurements the flight controller is able to calculate the drone current attitude (angle of flying ) and perform necessary corrections to the measurements. 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Time variations in the features of the tabular dataset for small fixed-wing UAV on 12 th and 13 th July. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Time variations in the features of the tabular dataset for mall fixed-wing UAV on 21 st and 23 rd July. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.1 Proposed architecture of Adaptive Convolutional Neural Network (ACNN). a) An encoderdecoder architecture composed of a 6-layer FCN followed by ReLU activation function. b) Relation between proposed FCNs and CNNs; the weights are split into sub-matrices simulating convolutive filters in CNNs c) Each of the sub-matrices is constrained to have a band structure as shown in this example. The dark gray area marks the zero-entries, while the light-gray colour corresponds to the ones that are allowed to be non-zero. . . . . . . . . . B.2 Convergence profile of the proposed method. . . . . . . . . . . . . . . . . . . . . . . . . . List of Tables 1.1 Different classes of AI systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Error metrics used for evaluation. The mean value computed for SMAPE is limited to the K + positive values of the elements in the summation. e k is the perturbation generated by the adversarial attacker on the k-th sample in the dataset of length K. . . . . . . . . . . . 3.2 Input and output variables description for the Thales Air Mobility industrial application dataset. 3.3 Comparison on evaluation metrics random attacker vs. proposed adversarial attacker with variation in perturbation level (ℓ 2 attack). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Standard training vs Spectral Normalization training on ℓ 2 attacks. . . . . . . . . . . . . . 3.5 Standard training attacking all inputs vs standard training attacking few inputs on ℓ 2 attacks. 3.6 Comparison on industrial dataset for ℓ 2 , ℓ 1 and ℓ ∞ attacks with variation in perturbation levels. MAE std (MAE for Standard training) = 9.2 × 10 -3 . . . . . . . . . . . . . . . . . . 4.1 Comparison of state-of-the-art Lipschitz estimation approaches vs the proposed one. . . . . 4.2 Comparison of Lipschitz constant values when γ = 0. Test performance for standard training: NMSE = 0.007, NMAE = 0.005, for spectral normalization: NMSE = 0.011, NMAE = 0.009. 4.3 Comparison of Lipschitz constant values when γ = 1/10. Test performance for standard training: NMSE = 0.006, NMAE = 0.005, for spectral normalization: NMSE = 0.009, NMAE = 0.007. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Comparison of Lipschitz constant values when γ = 1. Test performance for standard training: NMSE = 0.006, MAE =0.005, for spectral normalization: NMSE = 0.014, NMAE = 0.009. 4.5 Input and output attributes of Combined Cycle Power Plant dataset. . . . . . . . . . . . . 4.6 Input and output attributes of Auto MPG dataset. . . . . . . . . . . . . . . . . . . . . . . 4.7 Input and output attributes of Boston Housing dataset. . . . . . . . . . . . . . . . . . . . 4.8 Network Architecture and training setup for different datasets. . . . . . . . . . . . . . . . 4.9 2 nd order normalized coupling matrix with ϵ = 0.001 on a) Combined Power Plant Dataset b) Auto MPG Dataset c) Boston Housing Dataset and d) Thales Air Mobility industrial application. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Results on German Credit dataset using the proposed algorithm. . . . . . . . . . . . . . . 5.2 Input and output attributes of UAV dataset for fault detection. a * -Linear acceleration , ω * -angular velocity , u * -autopilot command controls. . . . . . . . . . . . . . . . . . . . . . 5.3 Binary Classification : Comparison between SVM, logistic regression, and best NN trained with baseline training on 12 th July. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 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. . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Comparison between models trained with baseline and spectral normalization constraint (θ target = 4). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Binary Classification: Results on 12 th July clean data, FGSM and PGD attacks using model trained on 12 th July. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 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 . . . . . . . . . . . . . 5.8 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. . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: By adding an unnoticeable perturbation, an image labelled as "panda" is classified as "gibbon" [1].

Figure 2 . 1 :

 21 Figure 2.1: m-layered feed-forward neural network architecture. For the i th -layer, W i is the linear weight operator, b i the bias vector, and R i the activation operator.

Figure 2 . 2 :

 22 Figure 2.2: A network T is split into three sub-networks U 1 , U 2 , and U 3 .

  Distillation: Adversarial training needs adversarial examples to train the model; thus, the defense is related to the process of adversarial examples construction. In Papernot et al. [52] proposed a universal defensive method for neural networks called defensive distillation. The distillation method uses a small model to simulate a large and computationally intensive model without affecting the accuracy and can solve the problem of missing information. Different from the traditional distillation technique, defensive distillation aims to smooth the model during the training process by generalizing examples outside the training data. The basic idea of defensive distillation is to generate smooth classifiers that are more resilient to adversarial examples, reducing the sensitivity of the DNN to the input perturbation.

  e∈C p,δ ∥T(x + e) -y∥ q , (3.1)

Figure 3 . 1 :

 31 Figure 3.1: Network Architecture.

Figure 3 . 2 :

 32 Figure 3.2: Error distribution of random attacks and proposed adversarial attack on Combined cycle Power Plant dataset for perturbation level of 2 × 10 -1 for ℓ 2 attack.

4. 4 . 1 .

 41 Validation on Synthetic Datasets 4.4.1.1 . Context

Figure 4 . 1 :

 41 Figure 4.1: Lipschitz Star: Representing Partial Lipschitz constants w.r.t. to the inputs as the vertices of the radial plot.

4. 4 . 5 .

 45 Effect of training with specified Lipschitz target Spectral norm constrained training is performed as explained in Section 4.4.1.1. The results are shown for our three datasets in Figures 4.2-4.5. On these plots, we can observe a shrinkage of the Lipschitz stars following the reduction of the target Lipschitz value. Interestingly, improving stability does not affect significantly the performance of the networks. Let us comment on the last use case in light of the obtained results.

Figure 4 . 2 :Figure 4 . 3 :

 4243 Figure 4.2: Sensitivity w.r.t. to each input on Combined Cycle Power Plant dataset. Influence of a spectral normalization constraint. a) Standard training: Lipschitz constant = 0.66, MAE = 0.007 , b) With spectral normalization: Lipschitz constant = 0.25, MAE = 0.0066.

Figure 4 . 4 :

 44 Figure 4.4: Sensitivity w.r.t. to each input on Boston Housing dataset. Influence of a spectral normalization constraint. a) Standard training : Lipschitz constant = 18.56, MAE(y 1 ) = 2.45, MAE(y 2 ) = 1.41, b) With spectral normalization: Lipschitz constant = 8.06 , MAE(y 1 ) = 2.96, MAE(y 2 ) = 1.35.

Figure 4 . 5 :

 45 Figure 4.5: Sensitivity w.r.t. to each input on Thales Air Mobility industrial application. Influence of a spectral normalization constraint. a) Standard training: Lipschitz constant = 45.46, MAE = 496.37 (s), b) With spectral normalization constraint: Lipschitz constant = 16.62, MAE = 478.88 (s).

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: Sensitivity w.r.t. to each input on Combined Cycle Power Plant dataset. Effect of adversarial training. a) Standard training: Lipschitz constant = 0.657, MAE = 0.007, b) Adversarial training: Lipschitz constant = 0.37, MAE = 0.0068.

Figure 4 . 8 :

 48 Figure 4.8: Sensitivity w.r.t. to each input on Boston Housing dataset. Effect of adversarial training. a) Standard training : Lipschitz = 18.56, MAE(y 1 ) = 2.45, MAE(y 2 ) = 1.41, b) Adversarial training: Lipschitz constant = 16.50, MAE(y 1 ) = 2.35 MAE(y 2 ) = 1.32.

Figure 4 . 9 :

 49 Figure 4.9: Sensitivity w.r.t. to each input on Thales Air Mobility industrial application. Effect of adversarial training. a) Standard training: Lipschitz = 45.47, MAE = 496.37 (s), Adversarial training. b) Lipschitz = 34.26, MAE = 494.7 (s).

Figure 4 . 10 :

 410 Figure 4.10: Sensitivity w.r.t to pair of variables on Combined Power Plant dataset

Figure 4 . 11 :

 411 Figure 4.11: Sensitivity w.r.t to pair of variables on Auto MPG dataset.

Figure 4 . 12 :

 412 Figure 4.12: Sensitivity w.r.t to pair of variables on Boston Housing Dataset.

Figure 4 . 13 :

 413 Figure 4.13: Sensitivity w.r.t to pair of variables on Thales Air Mobility industrial application.

Figure 5 . 1 :

 51 Figure 5.1: Fault Detection in drones/UAVs : Inertial sensors of the drone: accelerometer measures translational acceleration, and the gyrometer measures angular velocities (rotational motion). By combining these measurements the flight controller is able to calculate the drone current attitude (angle of flying ) and perform necessary corrections to the measurements. 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 .

Figure 5 . 2 :

 52 Figure 5.2: Time variations in the features of the tabular dataset for small fixed-wing UAV on 12 th and 13 th July.

Figure 5 . 3 :

 53 Figure 5.3: Time variations in the features of the tabular dataset for mall fixed-wing UAV on 21 st and 23 rd July.

Figure B. 1 :

 1 Figure B.1: Proposed architecture of Adaptive Convolutional Neural Network (ACNN). a) An encoder-decoder architecture composed of a 6-layer FCN followed by ReLU activation function. b) Relation between proposed FCNs and CNNs; the weights are split into sub-matrices simulating convolutive filters in CNNs c) Each of the sub-matrices is constrained to have a band structure as shown in this example. The dark gray area marks the zero-entries, while the light-gray colour corresponds to the ones that are allowed to be non-zero.

Figure B. 2 :

 2 Figure B.2: Convergence profile of the proposed method. The dataset is divided into training, validation, and test sets. The training set contains 90 clips with variable length, ranging from 6 s up to 150 s, having in total over an hour (67 min) of audio recordings. Some songs are repeated on different octaves to obtain a minimum number of occurrences for all notes. The validation dataset contains 9 songs with length between 12 and 120 s, forming around 10 minutes of audio signals. The test set has one 2-minute long clip.

  1 Effect of Regularization Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Positive Weighted Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Addition of Noise to the Dataset . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Comments on the results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B Robustness of a New Type of Neural Networks using Positive Weights B.1 Making the bridge between CNNs and Fully Connected Networks (FCN) . . . . . . . . . . B.1.1 Learning algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Experimental Evaluation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2.1 Dataset Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2.2 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2.3 Simulations and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . By adding an unnoticeable perturbation, an image labelled as "panda" is classified as "gibbon" [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Intuition of using Lipschitz constant as a stability property of neural network [2]. . . . . . .

	List of Figures
	1.1
	Bibliography

Table 1 .

 1

1: Different classes of AI systems

Table 3 .

 3 2: Input and output variables description for the Thales Air Mobility industrial application dataset.

		0	Speed
		1	Flight Distance
		2	Departure Delay
		3	Initial ETE
		4	Latitude Origin
		5	Longitude Origin	continuous
	Input	6 7	Altitude Origin Latitude Destination
		8 Longitude Destination
		9	Altitude Destination
		10	Arrival Time Slot	7 slots (categorical)
		11	Departure Time Slot	7 slots (categorical)
		12	Aircraft Category	6 classes (categorical)
		13	Airline Company	19 classes (categorical)
	Output 3	Refinement ETE	continuous

  1

	Dataset	Power Plant	Wine		Abalone	Industrial
	MAE std	6.4 × 10 -3 6.4 × 10 -3	0.47	0.47	1.68	1.68	9.2 × 10 -3 9.2 × 10 -3
	MAE gauss MAE uni MAE bin MAE adv	6.5 × 10 -3 6.8 × 10 -3 6.5 × 10 -3 6.8 × 10 -3 6.5 × 10 -3 6.9 × 10 -3 10.3 × 10 -3 14.2 × 10 -3	0.46 0.46 0.47 0.58	0.47 0.47 0.48 0.66	1.68 1.68 1.68 2.04	1.68 1.68 1.68 2.40	9.6 × 10 -3 10.7 × 10 -3 9.6 × 10 -3 10.7 × 10 -3 9.6 × 10 -3 10.7 × 10 -3 20.9 × 10 -3 32.5 × 10 -3
	E guass E uni E bin E adv	1.3 × 10 -3 2.5 × 10 -3 1.3 × 10 -3 2.5 × 10 -3 1.4 × 10 -3 2.7 × 10 -3 4.0 × 10 -3 8.0 × 10 -3	0.04 0.05 0.04 0.12	0.09 0.09 0.09 0.21	0.02 0.02 0.03 0.36	0.05 0.05 0.05 0.72	2.6 × 10 -3 5.1 × 10 -3 2.6 × 10 -3 5.2 × 10 -3 2.7 × 10 -3 5.4 × 10 -3 11.8 × 10 -3 24.0 × 10 -3
	SMAPE gauss	0.33	0.34	0.34	0.49	0.05	0.09	0.45	0.65
	SMAPE uni	0.34	0.52	0.33	0.51	0.04	0.08	0.46	0.66
	SMAPE bin	0.36	0.56	0.29	0.48	0.05	0.09	0.47	0.67
	SMAPE adv	0.62	0.87	0.41	0.56	0.38	0.58	0.96	1.24

Table 3 . 3

 33

: Comparison on evaluation metrics random attacker vs. proposed adversarial attacker with variation in perturbation level (ℓ 2 attack).

  × 10 -1 11.8 × 10 -3 9.1 × 10 -3

		E adv	E spec	SMAPE adv SMAPE spec
			Power Plant
	1 × 10 -1 4.0 × 10 -3 3.9 × 10 -3	0.62	0.60
	2 × 10 -1 8.0 × 10 -3 7.7 × 10 -3	0.87	0.84
			Wine
	1 × 10 -1	0.12	0.02	0.41	0.12
	2 × 10 -1	0.21	0.03	0.56	0.17
			Abalone
	5 × 10 -2	0.36	0.11	0.38	0.12
	1 × 10 -1	0.72	0.23	0.58	0.21
			Industrial
	1 0.91	0.49
	2 × 10 -1 24.0 × 10 -3 17.5 × 10 -3	1.24	0.72

Table 3 . 4

 34

: Standard training vs Spectral Normalization training on ℓ 2 attacks.

Table 3 .

 3

5: Standard training attacking all inputs vs standard training attacking few inputs on ℓ 2 attacks.

Table 3 .

 3 3 ) 6: Comparison on industrial dataset for ℓ 2 , ℓ 1 and ℓ ∞ attacks with variation in perturbation levels. MAE std (MAE for Standard training) = 9.2 × 10 -3

		9.6	10.7	9.2	9.4	10.5	13.5
	MAE uni (×10 -3 ) MAE bin (×10 -3 ) MAE adv (×10 -3 )	9.6 9.6 20.9	10.7 10.7 32.5	9.2 9.2 18.5	9.3 9.3 28.1	11.1 13.0 31.1	15.3 22.0 52.5
	E gauss (×10 -3 ) E uni (×10 -3 ) E bin (×10 -3 ) E adv (×10 -3 )	2.6 2.6 2.7 11.8	5.1 5.2 5.4 24.0	0.84 0.81 0.72 9.5	1.7 1.6 1.4 20.0	4.7 6.0 9.9 22.0	9.4 12.0 19.0 45.0
	SMAPE gauss	0.45	0.65	0.22	0.35	0.63	0.84
	SMAPE uni	0.46	0.66	0.21	0.34	0.71	0.92
	SMAPE bin	0.47	0.67	0.20	0.32	0.87	1.08
	SMAPE adv	0.96	1.24	0.87	1.15	1.22	1.47

Table 4 .

 4 .1.

	Method	Properties	Sensitivity of Inputs
	Naive upper Bound [1]	spectral bound, loose bound, univariate	No
	SDPLip [16]	ℓ 2 norm, more scalable to broad networks, univariate	No
	CPLip [9]	ℓ p ∈ [1, +∞], not scalable to broad networks, univariate	No
	LipOpt-k [43]	ℓ p ∈ [1, +∞] , univariate	No
	Proposed	scalable to broad networks,multivariate	Yes

1: Comparison of state-of-the-art Lipschitz estimation approaches vs the proposed one.

LC Analytical Standard Spectral Normalized

	θ {1}	1	133.9	6.75
	θ {2}	100	211.7	76.3
	θ {3}	100	299.7	136.0
	θ {1,2}	100.0	229.0	102.2
	θ {1,3}	100.0	303.1	136.0
	θ {2,3}	141.4	314.2	141.2
	θ {1,2,3}	141.4	315.3	141.2

Table 4 .

 4

	Partial LC Analytical Standard Spectral Normalized
	θ {1}	6	138.7	10.1
	θ {2}	100	219.2	90.0
	θ {3}	105	302.7	138.9
	θ {1,2}	100.2	231.6	108.1
	θ {1,3}	105.2	306.3	139.0
	θ {2,3}	145	316.4	147.2
	θ {1,2,3}	145.1	316.5	147.2

2: Comparison of Lipschitz constant values when γ = 0. Test performance for standard training: NMSE = 0.007, NMAE = 0.005, for spectral normalization: NMSE = 0.011, NMAE = 0.009.

Table 4 .

 4 3: Comparison of Lipschitz constant values when γ = 1/10. Test perfor-

	mance for standard training: NMSE = 0.006, NMAE = 0.005, for spectral normalization:
	NMSE = 0.009, NMAE = 0.007.

Table 4 .

 4 4: Comparison of Lipschitz constant values when γ = 1.

	Partial LC Analytical Standard Spectral Normalized
	θ {1}	51	274.7	59.5
	θ {2}	100	298.9	80.3
	θ {3}	150	388.7	183.7
	θ {1,2}	112.6	337.0	119.4
	θ {1,3}	158.4	392.2	183.7
	θ {2,3}	180.3	400.1	188.9
	θ {1,2,3}	187.4	400.5	189.0
				Test performance
	for standard training: NMSE = 0.006, MAE =0.005, for spectral normalization: NMSE
	= 0.014, NMAE = 0.009.

Table 4 .

 4 -4.7, respectively. 5: Input and output attributes of Combined Cycle Power Plant dataset.

		0		Temperature
	Input	1 2		Ambient Pressure Relative Humidity	continuous
		3		Exhaust Vacuum
	Output 4 Net hourly electrical energy output continuous
			0	Cylinders
			1 Displacement
		Input	2 3	Horsepower Weight	continuous
			4	Acceleration
			5	Model Year
			6	Origin	unsorted categorical
	Output 7	MPG	continuous

Table 4 .

 4 6: Input and output attributes of Auto MPG dataset.

Table 4 .

 4 .2. 7: Input and output attributes of Boston Housing dataset.

		0	CRIM
		1	ZN
		2	INDUS
		3	CHAS
		4	NOX
		5	RM
	Input	6	AGE	continuous
		7	DIS
		8	RAD
		9	TAX
		10	B
		11	LSTAT
	Output	12 13 PTRATIO continuous MEDV continuous

Table 4 .

 4 8: Network Architecture and training setup for different datasets.

	Dataset	layers & neurons EPOCHS Optimizer	lr
	Combined Cycle Power Plant	(10, 6)	100	Adam	0.01
	Auto MPG	(16, 8)	1000	RMSprop	0.001
	Boston Housing	(10, 5)	500	RMSprop	0.001
	Thales Air Mobility App.	(10 × (30))	100	Adam	0.01

Table 4 .

 4 9: 2 nd order normalized coupling matrix with ϵ = 0.001 on a) Combined Power

	Plant Dataset b) Auto MPG Dataset c) Boston Housing Dataset and d) Thales Air Mo-
	bility industrial application.

1 :

 1 Inputs : training dataset (X train , Y train ), validation dataset (X valid , Y valid ), target Lipschitz constant θ target , target accuracy acc target , neural network T, multiplicative factor α > 1, maximum iteration number n max 2: Output : optimal neural network T opt , Lipschitz constant θ opt , and accuracy acc opt : L ← θ target , n ← 0 5: T ← Train SN (θ, X train , Y train ) 6: θ est ← LipEst(T) 7: acc ← Predict(T, X valid , Y valid ) 8: while (n < n max ) and (θ est < θ target ) and (acc < acc target ) do Train SN (L, X train , Y train ) Predict(T, X valid , Y valid ) T opt ← T 16: θ opt ← θ est 17: acc opt ← acc 18: if n = n max then

		θ target	θ	acc(%) acc target (%) θ prop acc prop (%)
		1	1.02	69.4	70	1.00	70.0
		2	1.56	70.4	70	1.95	71.2
		3	1.90	70.9	70	2.73	72.8
		4	2.49	71.2	70	3.50	74.2
		5	2.78	72.3	70	4.73	74.4
		8	3.42	72.6	75	7.59	75.2
		12	5.2	73.1	75	10.11	76.4
		14	5.8	73.7	75	13.29	75.6
	9:	θ ← α θ
	10:
	11:	θ est ← LipEst(T)
	12:
	13:	n ← n + 1
	14: end while
	15: 19:	Notification: No trade-off found between Stability and Accuracy
	20: else
	21:	Notification: Trade-off found between Stability and Accuracy
	22: end if

3: Algo : 4T ← acc ←

Table 5 .

 5 1: Results on German Credit dataset using the proposed algorithm.

Table 5 .

 5 2: Input and output attributes of UAV dataset for fault detection. a

* -Linear acceleration , ω * -angular velocity , u * -autopilot command controls.

multi-labeled data. Data have been split into two parts containing 80% training data and 20% test data.

  While training, the dataset Results on 13 th July trained on 12 th July data.

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

(b) Multi-class Classification: Results on 23 rd July trained on 21 st July data.

Table 5 .

 5 5: Comparison between models trained with baseline and spectral normalization constraint (θ 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 5 .

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

Table 5 .

 5 7: 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 .

Table A .

 A 1: Results on Combined Cycle Power Plant Data Set for ℓ 1 and ℓ 2 regulariza-

									ℓ 1			ℓ 2
		α		0	0.00001 0.0001 0.001	0.01	0.1	0.00001 0.0001 0.001	0.01	0.1
		MAE	0.0069 0.0071 0.0069 0.0094 0.0300	0.0300	0.0073	0.007 0.0077 0.0299	0.030
	L 2,2	LipSDP [16] CPLip[9]	0.657 0.638	0.705 0.681	0.26 0.25	0.024 0.024 5.33e-11 3.99e-11 ≈ 0 ≈ 0	1.01 0.96	0.741 0.73	0.027 0.027 1.26e-09 1.09e-17 ≈ 0 ≈ 0
	L+∞,+∞	Lipopt-k [43] CPLip[9]	1.41 1.23			1.39 1.17	0.47 0.46	0.028 0.028 9.16e-11 6.75e-11 ≈0 ≈ 0	1.89 1.486	1.26 1.26	0.040 1.07e-09 0.040 2.35e-09 1.95e-17 ≈ 0
	tion
									ℓ1			ℓ2
		α	0	0.00001 0.0001 0.001	0.01	0.1	0.2	0.00001 0.0001 0.001	0.01	0.1	0.2
		MAE	0.0418 0.0444 0.0405 0.0443 0.0505 0.1490	0.1490	0.0515 0.0404 0.0424 0.0454 0.1489	0.1489
	L2,2	LipSDP [16] CPLip[9]	2.75 2.747	1.74 1.705	0.38 0.373	0.16 0.16	0.11 0.110 8.75e-09 7.44e-09 ≈ 0 ≈ 0	2.13 2.11	0.78 0.721	0.201 0.201	0.089 0.089 9.96e-09 9.72e-09 ≈ 0 ≈ 0
	L+∞,+∞	Lipopt-k [43] CPLip[9]	9.47 6.98	5.66 4.36		1.04 1.03	0.286 0.1642 1.66e-08 1.83e-08 0.29 0.16 2.30e-08 2.06e-08	7.048 4.97	3.088 0.4365 0.2099 1.41e-08 1.45e-08 1.93 0.44 0.21 2.9e-08 2.8e-08

Table A .

 A 2: Results on Auto MPG Data Set for ℓ 1 and ℓ 2 regularization

							ℓ1					ℓ2
		α	0	0.00001 0.0001 0.001	0.01	0.1	0.2	0.00001 0.0001 0.001 0.01	0.1	0.2
		Acc	84.94	85.16	85.57 85.54	76.30	76.30	76.30	85.32	85.26 85.46 84.73	76.32	76.32
	L2,2	LipSDP[16] CPLip[9]	8.21 8.208	6.29 6.29	5.15 5.15	3.45 3.45 5.81e-10 1.84e-10 3.02e-10 ≈ 0 ≈ 0 ≈ 0	9.15 9.15	5.32 5.32	4.19 4.19	1.91 1.91 1.839e-10 6.04e-11 ≈ 0 ≈ 0
	L+∞,+∞	Lipopt-k [43] 56.22 CPLip[9] 56.22	31.77 31.77	20.53 12.36 20.53 12.36 4.33e-09 1.34e-09 2.21e-09 ≈ 0 ≈ 0 ≈ 0	42.74 42.74	22.48 17.90 9.06 22.48 17.90 9.06	2.43e-11 9.25e-10 3.85e-10 ≈ 0

Table A .

 A 3: Results on Adult Data Set for ℓ 1 and ℓ 2 regularization

  Lipopt-k[START_REF] Latorre | Lipschitz constant estimation of neural networks via sparse polynomial optimization[END_REF] 56.22 46.07 49.32 43.58 40.21 39.87 41.26 CPLip[9] 56.22 46.07 49.32 43.58 40.21 39.87 41.26 Table A.6: Results on Adult Data Set with Dropout

		Drop-rate	0	0.05	0.1	0.2	0.3	0.4	0.5
		MAE	0.0069 0.0069 0.0074 0.0086 0.0115 0.0111 0.016
	L 2,2	LipSDP [16] CPLip[9]	0.66 0.64	0.16 0.16	0.23 0.23	0.56 0.54	0.45 0.42	0.27 0.27	0.40 0.39
	L +∞,+∞	Lipopt-k [43] CPLip[9]	1.41 1.23	0.26 0.26	0.32 0.32	1.01 1.01	0.79 0.69	0.42 0.42	0.69 0.69
	Table A.4: Results on Combined Cycle Power Plant Data set with Dropout
		Drop-rate	0	0.05	0.1	0.2	0.3	0.4	0.5
		MAE	0.042 0.039 0.0364 0.043 0.044 0.045 0.050
	L 2,2	LipSDP [16] CPLip[9]	2.75 2.75	1.86 1.73	2.1 2.1	1.89 1.87	2.28 2.28	1.88 1.88	1.41 1.42
	L +∞,+∞	Lipopt-k [43] 9.47 CPLip [9] 6.98	6.41 4.45	5.49 4.631	4.89 4.89	5.2 5.19	3.98 3.98	3.07 3.07
		Table A.5: Results on Auto MPG Data set with Dropout
		Drop-rate	0	0.05	0.1	0.2	0.3	0.4	0.5
		Acc	84.94 85.13 85.14 85.08 85.03 85.09 84.81
	L 2,2	LipSDP [16] CPLip[9]	8.21 8.21	7.17 7.17	7.39 7.39	6.82 6.82	6.53 6.53	6.68 6.68	6.63 6.63
	L +∞,+∞
				106

Table A .

 A .8. 8: Results on Combined Cycle Power Plant Data set with added noise

			Dataset 1	Dataset 2	Dataset 3
			Arbitrary Positive Arbitrary Positive Arbitrary Positive
		MAE/ACC	0.0069	0.021	0.042	0.08	84.94	83.54
	L 2,2	LipSDP [16] CPLip [9]	0.66 0.64	0.03 0.03	2.75 2.75	0.57 0.57	8.21 8.21	3.65 3.65
	L +∞,+∞	Lipopt-k [43] CPLip[9]	1.41 1.23	0.06 0.06	9.47 6.98	3.29 1.24	56.22 56.22	18.81 18.81
		Table A.7: Results with positive constraint on the weights
		std	No Noise	0.01	0.05	0.1	0.2
		MAE	0.0069	0.0064 0.0061 0.0065 0.0071
	L 2,2	LipSDP [16] CPLip[9]	0.66 0.64	0.61 0.61	0.34 0.33	0.37 0.36	0.10 0.09
	L +∞,+∞	Lipopt-k [43] CPLip[9]	1.411 1.234	1.12 1.12	0.58 0.57	0.61 0.59	0.24 0.16

  (1β t 1 ) η i,t+1 = P S i,t η i,tγ t µ i,t /( √ ν i,t + ϵ) ,

	(B.3)

This definition of model stability and robustness comes from the standardization working group EUROCAE WG-114 and SAE G34.

More generally, a function ρ i,k can be applied to each component ξ i,k but this situation rarely happens in standard neural networks.

The result remains valid if different scalar activation functions are used in a given layer.

This can be extended to multiple output network, if the output space is equipped with the ℓ +∞ norm.

We assume that J(x) is defined at x, see[START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient methods and deep learning[END_REF] for a justification of this assumption in the non-smooth case.

Recall that this interpretation is valid when p =

in Eq. (4.11).[START_REF] Weng | Evaluating the robustness of neural networks: An extreme value theory approach[END_REF] ⪯ designates the Loewner order.

https://archive.ics.uci.edu/ml/datasets/Combined+Cycle+Power+Plant

https://archive.ics.uci.edu/ml/datasets/auto+mpg

https://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html

We present the results only for ReLU, but we tested our approach with other activation functions such as tanh as well and found the trends in sensitivity of inputs to be similar.

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.

• We also demonstrate the effectiveness of spectral normalized models in handling adversarial tasks.

https://speed.pub.ro/downloads/

in controlling Lipschitz properties in GANs. Recall that, given an m layer fully connected architecture and a Lipschitz target θ, we constrain the spectral norm of each layer to be less than m √ θ. This ensures that the upper bound on the global Lipschitz constant is less than θ. We keep the network architectures exactly the same as the standard training procedures in the previous section. The performance of adversarial attacker on standard and spectrally normalized trained model in terms 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 com r , 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 5.3. 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 that it is possible to train neural networks with much better generalization performance on different days with varying weather conditions.

In the 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

. Summary

The contributions of this work are summarised as follows:

• We 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 normalization constraints to achieve an accuracystability trade-off.

• We first show the efficacy of our proposed algorithm on an open source dataset. Our algorithm achieved an optimal solution with a trade-off between performance and stability.

• Next, we analyse real flight data application for fault detection function embedded in an UAV. 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.

Appendix A Lipschitz Constant Analysis of Tabular Data

This appendix describes one of our preliminary studies related to Lipschitz analysis on tabular data published in [START_REF] Gupta | A quantitative analysis of the robustness of neural networks for tabular data[END_REF].

Deep learning on tabular data has received much less attention than deep learning for standard signal/image processing problems particularly seen in the area of computer vision and natural language processing. Tabular data allows to take advantage of heterogeneous sources of information coming from different sensors or registered variables (like altitude of an aircraft, departure and destination airport, duration of the flight, company type). Tabular data analysis covers a wide variety of applications, e.g. fraud detection,product failure prediction, anti-money laundering, recommendation systems [START_REF] Covington | Deep neural networks for youtube recommendations[END_REF], click-through rate prediction [START_REF] Guo | Deepfm: a factorization-machine based neural network for ctr prediction[END_REF] etc. A generalised NN framework for tabular data is presented in [START_REF] Ke | TabNN: A universal neural network solution for tabular data[END_REF]. Most of the mentioned tasks may be hampered with safety concerns and require reliability in the performance of the NN used for predicting or classifying data. In [START_REF] Ballet | Imperceptible adversarial attacks on tabular data[END_REF], authors presented a method for generation of imperceptible adversarial attacks for tabular data.

In our stability analysis, we study three widely used tabular datasets from UCI Repository: 1) Combined Cycle Power Plant Data Set has 4 attributes with 9568 instances , 2) Auto MPG Data Set consists of 398 instances with 7 attributes, 3) Adult Data Set consists of 48842 instances in total with 14 attributes. Dataset 2 and 3 include both continuous and categorical attributes, whereas dataset 1 contains only continuous attributes. The datasets are divided with a ratio of 4:1 between training and testing data. The categorical attributes are dealt with by using one hot encoding based on the number of categories. The input attributes are normalised by removing their mean and scaling to unit variance. All the architectures are made up of two hidden fully connected layers with the following characteristics:

• Combined Cycle Power Plant Data set -(4, 10, 6, 1) • Auto MPG Data set - [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF][START_REF] Fazlyab | Efficient and accurate estimation of Lipschitz constants for deep neural networks[END_REF][START_REF] Virmaux | Lipschitz regularity of deep neural networks: analysis and efficient estimation[END_REF][START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] • Adult Data set - [START_REF] Zhang | Deep learning over multi-field categorical data[END_REF][START_REF] Singh | Trustworthy AI[END_REF][START_REF] Singh | Trustworthy AI[END_REF][START_REF] Goodfellow | Explaining and harnessing adversarial examples[END_REF] Appendix B

Robustness of a New Type of Neural Networks using Positive Weights

This appendix presents the work done in collaboration with another PhD student at CVN. The work is published in [START_REF] Neacsu | Signal denoising using a new class of robust neural networks[END_REF].

B.1 . Making the bridge between CNNs and Fully Connected Networks (FCN)

We use the fact that, if all the weights of a feed-forward neural network are non-negative, an optimal Lipschitz constant reduces to ∥W m • • • W 1 ∥ S [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF]. In other words, it has the same Lipschitz constant as a linear network where the identity function is substituted for all the activation operators. For a network with weights having arbitrary signs, ϑ m = ∥W m • • • W 1 ∥ S only constitutes a lower bound on the Lipschitz constant [START_REF] Combettes | Lipschitz certificates for neural network structures driven by averaged activation operators[END_REF].

In this work, we aim at filling the gap between FCNs and CNNs. In signal processing context, a convolutive layer is a Multiple-Input Multiple-Output (MIMO) filter. For one-dimensional signals, each of these filters can be viewed as a Toeplitz matrix generated by the impulse response of the filter, which is applied to the vector of signal samples. If the filter length is short, large upper and lower triangular parts of this matrix are null. In our proposed approach, we will keep this band structure for the weight matrix, which is equivalent to performing local processing at each time within a sliding window. However, in order to add more flexibility in this architecture, we will allow all the nonzero coefficients of this matrix to be fully optimized. the network to a given value ϑ > 0. The related constraint can be expressed as:

with the convention that for the first layer (i = 1), B i,t is the N i × N i identity matrix Id N i and for the last layer (i = m), A i,m = Id N i . Above, (W j,t ) 1⩽j⩽m designates the estimates of the weight matrices at iteration t of the projected ADAM optimizer.

The projection onto the intersection of the above three closed convex sets has no closed-form expression. The intersection D i ∩ E R i ×Q i i is however quite simple to handle since the projection onto this set reduces to

) of a weight matrix W i , we use the dual forward-backward algorithm, as presented in [START_REF] Abboud | Dual block-coordinate forward-backward algorithm with application to deconvolution and deinterlacing of video sequences[END_REF].

B.2 . Experimental Evaluation

The proposed network has been evaluated for denoising music signals.

B.2.1 . Dataset Description

We train our proposed ACNN on a dataset consisting of musical exercises and songs performed on a Ronald organ. The organ covers 5 octaves (range C2-C7), each octave having 12 semitones, generating a total of 61 different possible notes. For the recordings, the whole range of notes is used. The songs are recorded in MIDI format using MidiEditor, in the following manner: the recording mode from MidiEditor is activated before each song being played and is stopped after the song is finished (so there is silence at the beginning and end of each recording). In total, the dataset contains 100 MIDI recordings, with a sampling frequency F s = 44100 Hz, constituting 1 h and 17 min of audio. The data set is available online 1 . of the architecture for three values ϑ equal to 1, 5, and 10. In the second experiment, we test our architecture by varying the number of channels, i.e. the way we split each weight matrix. Note that for both scenarios we consider a network with m = 6 layers and ∀i ∈ {1, . . . , 6}, [START_REF] Tsipras | Robustness may be at odds with accuracy[END_REF][START_REF] Nguyen | Adversarial attacks, regression, and numerical stability regularization[END_REF][START_REF] Nguyen | Adversarial attacks, regression, and numerical stability regularization[END_REF][START_REF] Tsipras | Robustness may be at odds with accuracy[END_REF]237). We evaluate the performance on 3 standard metrics: Peak to Signal Noise Ratio (PSNR), Mean squared error (MSE), and Cross-correlation (CC), as shown in the Table B.1. We compare our method with a standard denoising technique, based on a Wavelet decomposition. We employ a 5-level decomposition using Symlet8 filters combined with SureShrink thresholding. We also compare ACNN with a FCN implementation, for which we ensured the Lipschitz bound ϑ = 1. Here, the FCN has no structural constraints but, in addition to the imposed Lipschitz property, it has positive weights.

Table B.1 reports the quantitative results obtained by all three approaches, evaluated over the test set. Our method outperforms the baseline on all measures by a significant margin. ACNN also outperforms classical FCN, inferring that structure imposed on the weight matrix for ACNN leads to a model with better generalization power, whereas FCN implementation may be prone to overfitting. We observed that, without any Lipschitz constraint, FCNs tend to have a very high Lipschitz constant of the order 10 5 -10 6 , which emphasizes the importance of controlling the Lipschitz behaviour of the system. The Lipschitz constant of the network is also closely related to the expressiveness of the trained model. Architectures with tight robustness constraint have fewer degrees of freedom and are thus expected to be less accurate, and lead to a slower convergence. The convergence profile of our proposed method is shown in Figure B.2.