Mariia Vladimirova
email: mariia.vladimirova@inria.fr

Julyan Arbel
email: julyan.arbel@inria.fr

Pablo Mesejo
email: pmesejo@decsai.ugr.es

Bayesian neural networks become heavier-tailed with depth

.

Introduction

Bayesian approaches investigate models by assuming a prior distribution on their parameters. Bayesian machine learning refers to extending standard machine learning approaches with posterior inference, a line of research pioneered by the works [START_REF] Neal | Bayesian training of backpropagation networks by the hybrid Monte Carlo method[END_REF]; [START_REF] Mackay | A practical Bayesian framework for backpropagation networks[END_REF] on Bayesian neural networks. The interest of the Bayesian approach to NNs is at least twofold. First, it offers a principled approach for modeling uncertainty of the training procedure, which is a limitation of standard NNs that only provide point estimates. A second main asset of Bayesian models is that they represent regularized versions of their classical counterparts. For instance, maximum a posteriori (MAP) estimation of a Bayesian regression model with double exponential (Laplace) prior is equivalent to Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF], while a Gaussian prior leads to ridge regression. When it comes to neural networks, the regularization mechanism is also appreciated in the literature, since neural networks traditionally suffer from overparameterization, resulting in overfitting.

Central in the field of regularization techniques is the weight decay penalty [START_REF] Krogh | A simple weight decay can improve generalization[END_REF], which is equivalent to MAP estimation of a Bayesian neural network with independent Gaussian priors on the weights. [START_REF] Srivastava | Dropout: A simple way to prevent neural networks from overfitting[END_REF] have suggested dropout as a regularization method in which neurons are randomly turned off. [START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning[END_REF] proved that the neural network trained with dropout is equivalent to a probabilistic model, i.e. a deep Gaussian process [START_REF] Damianou | Deep Gaussian processes[END_REF]. It leads to the consideration of such neural networks as Bayesian models.

This study is devoted to the investigation of hidden units prior distributions in Bayesian neural networks under assumption of independent Gaussian weights. We first describe a fully connected neural network architecture as illustrated in Figure 2, Appendix D. Given an input x ∈ R N , the -th hidden layer unit activations are defined as

g ( ) (x) = W ( ) h ( -1) (x), h ( ) (x) = φ(g ( ) (x)), (1)
where W ( ) is a weight matrix including the bias vector. A nonlinear activation function φ : R → R is applied element-wise, which is called nonlinearity, g ( ) = g ( ) (x) is a vector of prenonlinearities, and h ( ) = h ( ) (x) is a vector of post-nonlinearities. When we refer to either pre-or post-nonlinearities, we will use the notation U ( ) .

2 Heavy-tailed result Definition 2.1 (Sub-Weibull random variable). A random variable X, that satisfies

P(|X| ≥ x) ≤ exp -x 1/θ /K for all x ≥ 0.
for K > 0, is called a sub-Weibull random variable with the tail parameter θ > 0, which is denoted by X ∼ subW(θ).

Informally, the tails of a subW(θ) distribution are dominated by (i.e. decay at least as fast as) the tails of a Weibull variable with the shape parameter equal to 1/θ [START_REF] Rinne | The Weibull distribution: a handbook[END_REF]. The larger tail parameter θ, the heavier the tails of the sub-Weibull distribution. We refer to [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF] for more details about sub-Weibull distribution.

Assumptions on neural network

Let all weights (including biases) be independent and have zero-mean normal distribution

W ( ) i,j ∼ N (0, σ 2 w ), (2)
for all 1 ≤ ≤ L, 1 ≤ i ≤ H -1 and 1 ≤ j ≤ H . Given some input x, such prior distribution induces by forward propagation (1) a prior distribution on the pre-nonlinearities and post-nonlinearities, whose tail properties are the focus of this section. To this aim, the nonlinearity φ is required to span at least half of the real line as follows. We introduce an extended version of the nonlinearity assumption from [START_REF] Matthews | Gaussian process behaviour in wide deep neural networks[END_REF]: Definition 2.2 (Extended envelope property for nonlinearities). A nonlinearity φ : R → R is said to obey the extended envelope property if there exist c 1 , c 2 , d 1 , d 2 ≥ 0 such that the following inequalities hold

|φ(u)| ≥ c 1 + d 1 |u| for all u ∈ R + or u ∈ R -, |φ(u)| ≤ c 2 + d 2 |u| for all u ∈ R. (3)
The interpretation of this property is that φ must shoot to infinity at least in one direction (R + or R -) at least linearly (first line of (3)), and also at most linearly (second line of (3)). Compactly supported nonlinearities such as sigmoid and tanh do not satisfy the extended envelope property but the majority of other nonlinearities do, including ReLU, ELU, SELU (see [START_REF] Klambauer | Self-normalizing neural networks[END_REF] for details), and others. Lemma 2.1. Let a nonlinearity φ : R → R obey the extended envelope property. Then for any symmetric random variable X the following asymptotic equivalence3 holds

E[φ(X) k ] E[X k ], for k → ∞. (4)

Main theorem

This section postulates the rigorous result with a proof sketch. In Supplementary material one can find proofs of intermediate lemmas and a covariance theorem which states the non-negative covariance between post-nonlinearities. Theorem 2.1 (Sub-Weibull units). Consider a feed-forward Bayesian neural network with Gaussian priors (2) with nonlinearity φ satisfying the extended envelope condition of Definition 2.2. Then conditional on the input x, the marginal prior distribution induced by forward propagation (1) on any unit (pre-or post-nonlinearity) of the -th hidden layer is sub-Weibull with optimal tail parameter θ = /2. That is for any 1 ≤ ≤ L, and for any 1 ≤ m ≤ H , U ( ) m ∼ subW( /2), where a subW distribution is defined in Definition 2.1, and U ( )

m is either a pre-nonlinearity g ( ) m or a post-nonlinearity h ( ) m .
Before proving the result, we illustrate it in Figure 1 which represents the first three hidden layers pre-nonlinearity marginal distributions (left panel). These densities are obtained as kernel density estimators from a sample of size 10 5 from the prior on the pre-nonlinearities, which is itself obtained by sampling 10 5 sets of weights W from the Gaussian prior (2) and forward propagation via (1). The three hidden layers of neural network have H 1 = 25, H 2 = 24 and H 3 = 4 hidden units, respectively. Being a linear combination involving symmetric weights W , pre-nonlinearities g are also symmetric, thus we visualize only their positive part. The input vector x ∈ R 50 is sampled from a standard normal distribution once for all at the start. The nonlinearity φ is the ReLU function. The prior distribution of post-nonlinearities has a Dirac mass at zero with a coefficient of 1/2 and they are no more symmetric. But the post-nonlinearity prior distribution tails remain the same as of pre-nonlinearities on R + , and is represented on the right panel of Figure 1. Proof of Theorem 2.1. The idea is to prove by induction with respect to hidden layer depth that preand post-nonlinearities satisfy the asymptotic moment equivalence g ( )

k
k /2 and h ( ) k k /2 . The statement of the theorem then follows by the moment characterization of optimal sub-Weibull tail coefficient in Proposition A.1.

Base step: the first hidden layer pre-nonlinearity g follows normal distribution, since g = W T i x ∼ N (0, σ 2 w x 2 ). Then, for normal zero-mean variable g from Lemma B.1: g k √ k. As activation function φ obeys extended envelope property (Definition 2.2), according to Lemma B.2, nonlinearity moments are asymptotic equivalent to symmetric variable moments φ(g) k g k ∼ √ k. It implies that first hidden layer post-nonlinearities h is sub-Gaussian or sub-Weibull with tail parameter θ = 1/2 (Definition 2.1).

Inductive step: let us suppose the post-nonlinearity of ( -1)-th hidden layer satisfies the moment condition. Hidden units satisfy the non-negative covariance theorem (Theorem C.1):

Cov h ( -1) s , h( -1) t ≥ 0, for any s, t ∈ N.

Let the number of hidden units in ( -1)-th layer equals to H. Then according to Lemma B.3, under assumption of zero-mean Gaussian weights, pre-nonlinearities of -th hidden layer g ( ) = H j=1 W

( -1) i,j h ( -1) j also satisfy the moment condition, but with θ = /2 g ( )

k k /2 .
Using Lemma B.2 from Supplementary material, one can show that post-nonlinearities h ( ) satisfy the same moment condition as pre-nonlinearities g ( ) . This finishes the proof.

A Equivalent sub-Weibull distribution properties

Proposition A.1 (Sub-Weibull distribution). Let X be a random variable. Then the following properties are equivalent; the parameters K i > 0 appearing in these properties differ from each other by at most an absolute constant factor.

1. The tails of X satisfy

P(|X| ≥ x) ≤ 2 exp -x 1/θ /K 1 for all x ≥ 0.
2. The moments of X satisfy

X k = E[|X| k ] 1/k ≤ K 2 k θ for all k ≥ 1. 3. The MGF of X 1/θ satisfies E exp λ 1/θ X 1/θ ≤ K 2 exp(K 1/θ 3 λ 1/θ )
for all λ such that |λ| ≤ 1 K3 .

4. The MGF of X 1/θ is bounded at some point, namely

E exp X 1/θ /K 4 ≤ 2.
Proof. 1 ⇒ 2. Assume property 1 holds. Applying the integral identity for |X| k , we obtain

E |X| k = ∞ 0 P |X| k > x dx = ∞ 0 P |X| > x 1/k dx ≤ ∞ 0 2 exp -x 1/(kθ) /K 1 dx = 2K kθ 1 kθ ∞ 0 e -u u kθ-1 du = 2K kθ 1 kθ Γ (kθ) ∼ K kθ 1 kθ (kθ -1) kθ-1 ∼ (K 1 kθ) kθ .
Taking the k-th root of the expression above yields property 2

X k (K 1 θ) θ k θ ≤ K 2 k θ , with K 2 = (K 1 θ) θ .
2 ⇒ 3. Assume property 2 holds. Recalling the Taylor series expansion of the exponential function, we obtain

E exp λ 1/θ X 1/θ = E 1 + ∞ k=1 (λ 1/θ |X| 1/θ ) k k! = 1 + ∞ k=1 λ k/θ E[|X| k/θ ] k! . Property 2 guarantees that E[|X| k ] ≤ K 2 k k/θ and E[|X| k/θ ] ≤ K 2 (k/θ) k for some K 2 . Stirling s approximation yields k! ≥ (k/e) k .
Substituting these two bounds, we get

E exp λ 1/θ X 1/θ ≤ ∞ k=1 λ k/θ K 2 (k/θ) k (k/e) k = ∞ k=0 K 2 (eλ 1/θ /θ) k = K 2 1 -eλ 1/θ /θ ,
provided that eλ 1/θ /θ < 1, in which case the geometric series above converges. To bound this quantity further, we can use the numeric inequality 1 1-x ≤ e 2x , which is valid for

x ∈ [0, 1/2]. It follows that E exp λ 1/θ X 1/θ ≤ K 2 exp 2eλ 1/θ /θ
for all λ satisfying |λ| ≤ θ 2e θ . This yields property 3 with K 3 = (2e/θ) θ .

3 ⇒ 4. Assume property 3 holds. Take λ = 1/K 4 , where K 4 ≥ K 3 /(ln 2ln K 2 ) θ . This yields property 4. 4 ⇒ 1. Assume property 4 holds. We may assume that K 4 = 1. Then, by Markov's inequality and property 3, we obtain

P |X| > x = P e |X| 1/θ > e x 1/θ ≤ E[e |X| 1/θ ] e x 1/θ ≤ 2e -x 1/θ /K1 .
This proves property 1 with K 1 = 1.

Remark A.1. The constant 2 that appears in some properties in Proposition A.1 does not have any special meaning. It is chosen for simplicity and can be replaced by other absolute constants.

B Intermediate lemmas

Introduce the definition of asymptotic equivalence between numeric sequences: (5)

Lemma B.1 (Gaussian moments). Let X be a normal random variable such that X ∼ N (0, σ 2 ), then the following asymptotic equivalence holds

X k √ k.
Proof. The moments of central normal absolute random variable |X| are equal to

E[|X| k ] = R |x| k p(x) dx = 2 ∞ 0 x k p(x) dx = 1 √ π σ k 2 k/2 Γ k + 1 2 . ( 6
)
We have the expression for the Gamma function

Γ(z) = 2π z z e z 1 + 1 12z + o 1 z . (7)
Substituting ( 7) into the central normal absolute moment (6), we obtain

E[|X| k ] = 1 √ π σ k 2 k/2 4π k + 1 k + 1 2e (k+1)/2 1 + 1 6(k + 1) + o 1 k = 2σ k √ 2e k + 1 e k/2 1 + 1 6(k + 1) + o 1 k .
Then the roots of absolute moments can be written in the form of

X k = σ e 1/(2k) k + 1 e 1 + 1 6(k + 1) + o 1 k 1/k = σ e √ k + 1 e 1/(2k) 1 + 1 6(k + 1)k + o 1 k 2 = σ e c k √ k + 1.
Here the coefficient c k denotes

c k = 1 e 1/(2k) 1 + 1 6(k + 1)k + o 1 k 2 → 1,
with k → ∞. Thus, asymptotic equivalence holds

X k √ k + 1 √ k.
Lemma B.2 (Nonlinearity moments). Let X be symmetric random variable, φ(x) be nonlinear function satisfying extended envelope property, then the following asymptotic equivalence holds φ(X) k X k .

Proof. According to extended envelope property,

E[φ(X) k ] E[X k ].
That means there exist constants d and D such that for all k ∈ N it holds

d ≤ E[φ(X) k ] E[|X| k ] ≤ D.
Observing that

d ≤ d 1/k ≤ φ(X) k X k ≤ D 1/k ≤ D ,
the bounding constants are d = min{1, d}, D = max{1, D}. It yields asymptotic equivalence φ(X) k X k . The lemma is proved.

Lemma B.3 (Multiplication moments).

Let W and X be independent random variables such that W ∼ N (0, σ 2 ) and for some p > 0 it holds

X k k p . (8)
Let W i be independent copies of W , and X i be copies of X, i = 1, . . . , H with non-negative covariance between moments of copies Cov X s i , X t j ≥ 0, for i = j, s, t ∈ N.

(9)

Then we have the following asymptotic equivalence

H i=1 W i X i k k p+1/2 . ( 10
)
The statement proof is based on mathematical induction, see [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF].

C Covariance theorem

Theorem C.1 (Non-negative covariance between hidden units). Consider the deep neural network described in Section 3 with assumptions from Section 3.1. The covariance between hidden units of the same layer is non-negative. Moreover, for given -th hidden layer units h ( ) and h( ) , it holds Cov h ( ) s h( ) t ≥ 0, where s, t ∈ N.

For first hidden layer = 1 there is equality for all s and t.

We refer to [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF]

Figure 1 :

 1 Figure 1: Illustration of the first three layers hidden units marginal prior distributions.

  Definition B.1 (Asymptotic equivalence). Two sequences a k and b k are called asymptotic equivalent and denoted as a k b k if there exist constants d > 0 and D > 0 such that d ≤ a k b k ≤ D, for all k ∈ N.

  for the proof. Figure 2: Neural network architecture and characterization of the -layer units prior distribution as sub-Weibull distribution with tail parameter /2 (see Definition 2.1).

	D Illustration
	input	1 st hid.	2 nd hid.	3 rd hid.	th hid.
	layer	layer	layer	layer	layer
	x	h (1)	h (2)	h (3)	. . . h ( )
					. . .
	. . .	. . .	. . .	. . .	. . . . . .
		subW( 1 2 ) subW(1) subW( 3 2 ) subW( 2 )

* Andalusian Research Institute in Data Science and Computational Intelligence (DaSCI), University of Granada, Granada, Spain Third workshop on Bayesian Deep Learning (NIPS 2018), Montréal, Canada.

See Definition B.1 for the asymptotic equivalence definition in Supplementary material.


Mariia Vladimirova
email: mariia.vladimirova@inria.fr

Julyan Arbel
email: julyan.arbel@inria.fr

Pablo Mesejo
email: pmesejo@decsai.ugr.es

Bayesian neural network priors at the level of units

Keywords: Bayesian neural network, heavy-tailed prior, sparsity

come

Introduction

Neural networks (NNs), and their deep extensions [START_REF] Goodfellow | Deep learning[END_REF], have largely been used in many research areas such as image analysis [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], signal processing [START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF], or reinforcement learning [START_REF] Silver | Mastering the game of go with deep neural networks and tree search[END_REF], just to name a few. The performances of such neural networks have greatly strengthened the line of research that aims at better understanding the driving mechanisms behind the effectiveness of deep neural networks. One important aspect of this line of research that has recently gained much attention is the study of distributional properties of the NNs through Bayesian inference.

Bayesian approaches investigate models by assuming a prior distribution on their parameters. Bayesian machine learning refers to extending standard machine learning approaches with posterior inference, a line of research pioneered by the works Neal (1992); MacKay (1992) on Bayesian neural networks which now extends to a broad class of models, including Bayesian Generative adversarial network [START_REF] Saatci | Bayesian GAN[END_REF]. See [START_REF] Nicholas | Deep learning: A Bayesian perspective[END_REF] for a review. The interest of the Bayesian approach to NNs is at least twofold. First, it offers a principled approach for modeling uncertainty of the training procedure, which is a limitation of standard NNs which only provide point estimates. A second main asset of Bayesian models is that they represent regularized versions of their classical counterparts. For instance, maximum a posteriori (MAP) estimation of a Bayesian regression model with double exponential (Laplace) prior is equivalent to Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF], while a Gaussian prior leads to ridge regression. When it comes to neural networks, the regularization mechanism is also well appreciated in the literature, since neural networks traditionally suffer from overparameterization, resulting in overfitting.

Central in the field of regularization techniques is the weight decay penalty [START_REF] Krogh | A simple weight decay can improve generalization[END_REF], which is equivalent to MAP estimation of a Bayesian neural network with independent Gaussian priors on the weights. [START_REF] Srivastava | Dropout: A simple way to prevent neural networks from overfitting[END_REF] have suggested dropout as a regularization method in which neurons are randomly turned off. [START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning[END_REF] proved that the neural network trained with dropout is equivalent to a probabilistic model, i.e. a deep Gaussian process [START_REF] Damianou | Deep Gaussian processes[END_REF]. It leads to the consideration of such neural networks as Bayesian models.

Recent papers study various distributional properties of Bayesian neural networks. Matthews et al. (2018b), or its extended version Matthews et al. (2018a), and [START_REF] Lee | Deep neural networks as Gaussian processes[END_REF] showed that deep neural networks tend in distribution to the Gaussian process when the number of hidden units grows to infinity, under the assumption of Gaussian weights for properly rescaled prior variances. The work by [START_REF] Bibi | Analytic expressions for probabilistic moments of pl-dnn with gaussian input[END_REF] provides the expression of the first two moments of the output units of a one layer neural network. Obtaining the moments is a first step to characterizing a whole distribution, however the methodology of [START_REF] Bibi | Analytic expressions for probabilistic moments of pl-dnn with gaussian input[END_REF] is limited to the first two moments and to one layer neural networks. [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF] investigate the marginal prior distribution of the network units, stating that as the depth increases, the distribution becomes more heavy-tailed. For the sake of completeness, we reproduce the result here:

Theorem 1 (Sub-Weibull units, [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF] Consider a feed-forward Bayesian neural network with Gaussian priors on weights and activation function φ satisfying the extended envelope condition of Definition A.1, Appendix A. Consider the marginal prior distribution, conditional on the input x, induced by forward propagation on any unit before or after activation of the -th hidden layer. It is sub-Weibull with optimal tail parameter θ = /2. That is for any layer with H hidden units, a unit U ( )

m of the layer (1 ≤ m ≤ H ) U ( ) m ∼ subW( /2
), where a subW distribution is defined in Definition A.2, Appendix A.

We provide in this note an interpretation of Theorem 1 in terms of sparsity-inducing mechanism at the level of the units. We stress that the term unit and notation U ( ) refer indistinctly to units before or after activation. The notation W ( ) stands for the weight matrix including the bias vector.

Sparsity-inducing prior on the units

Shrinkage is performed by imposing a penalty on the size of the parameters. Denote the parameters by W , the loss function by R(W ), and the penalty by λL (W ), where L is some norm on the weight space and λ some positive tuning parameter. Then, the regularized problem is

min W R(W ) + λL (W ).
The choice of the L norm has considerable effects on the problem, as can be sensed geometrically. Consider for instance L q norms, with q ≥ 0. For any q > 1, the associated L q norm is differentiable and contours have a round shape without sharp angles. In that case, the penalty effect is to shrink the W coefficients towards 0. The most well-known estimator falling in this class is the ridge regression obtained with q = 2, see Figure 1, Layer 1 and 2. In contrast, for any q ∈ (0, 1], the L q norm has some non differentiable points along the axis coordinates, see Figure 1, Layer 2, 3 and 10. Such critical points are more likely to be hit by the level curves of the loss function R(W ), thus setting exactly to zero some of the parameters. A very successful approach in this class is the Lasso obtained with q = 1. Note that the problem is computationally much easier in the convex situation which occurs only for q ≥ 1.

MAP on weights is weight decay

These penalized methods have a simple Bayesian counterpart in the form of the maximum a posteriori (MAP) estimator. In this context, the objective function R is the negative loglikelihood, while the penalty L is the negative log-prior. The objective function takes on the form of sum-of-squared errors for regression under Gaussian errors, and of cross-entropy for classification.

For neural networks, it is well-known that an independent Gaussian prior on the weights

π(W ) ∝ L =1 i,j e -1 2 (W ( ) i,j ) 2
, is equivalent to the weight decay penalty-aka ridge penalty in regression problem-with negative log-prior:

L (W ) ∝ L =1 i,j (W ( ) i,j ) 2 = W 2 2 ,
where products and sums involving i and j above are over 1 ≤ i ≤ H -1 and 1 ≤ j ≤ H , L is a number of neural network layers, H 0 and H L representing respectively the input and output dimensions.

MAP on units induces sparsity

Now moving the point of view from weights to units leads to a radically different shrinkage effect. Let U ( ) m denote the m-th unit of the -th layer (either before or after activation). As stated in Theorem 1, [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF] show that conditional on the input x, a Gaussian prior on the weights translates into some prior on the units U ( ) m that is marginally sub-Weibull with optimal tail index θ = /2. This means that the tails of U ( )

m satisfy P(|U ( ) m | ≥ u) ≤ exp -u 2 / /K for all u ≥ 0, (1)
for some positive constant K. The exponent of u in the exponential term above is optimal in the sense that Equation ( 1) is not satisfied with some parameter θ smaller than /2. Thus, the marginal density of U ( )

m is approximately proportional to π ( ) m (u) ≈ e -|u| 2/ /K , u ∈ R,
where K can be chosen the same for all -th layer units as the smallest constant satisfying equations ( 1) for all m. The joint prior distribution for all the units U = (U ( ) m ) 1≤ ≤L,1≤m≤H can be expressed from all the marginal distributions by Sklar's representation theorem as

π(U ) = L =1 H m=1 π ( ) m (U ( ) m ) C(F (U )), (2)
where C represents the copula of U (which characterizes all the dependence between the units) while F denotes its cumulative distribution function. The penalty incurred by such a prior distribution is obtained as the negative log-prior L (U ),

- L =1 H m=1 log π ( ) m (U ( ) m ) -log C(F (U )) ≈ U (1) 2 2 + • • • + U (L) 2 /L 2 /L -log C(F (U )). ( 3
)
The first L terms in the right-hand-side of (3) indicate that some shrinkage operates at the units level: at layer , the penalty term U (L) 2 / 2 / takes the form of the L 2 / norm. Thus, the deeper the layer, the stronger the sparsity at the level of the units.

Conclusion

We offer an interpretation in terms of sparsity-inducing mechanism of the heavy-tailed units distribution result derived by [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF] and reproduced here as Theorem 1. Heavy-tailed priors are known to induce a sparse model representation, such as the Lasso in a regression problem. We extrapolate this finding to the setting of a deep Bayesian neural network, showing that the units tend to be more sparsely represented as layers become deeper.

Figure 1 :

 1 Figure 1: L 2/ -norm unit balls (in dimension 2) for layers = 1, 2, 3 and 10.

© M. Vladimirova, J. Arbel & P. Mesejo.

Appendix A. Definitions

Definition A.1 (Extended envelope property for nonlinearities) A nonlinearity φ : R → R is said to obey the extended envelope property if there exist c 1 , c 2 , d 2 ≥ 0, d 1 > 0 such that the following inequalities hold

The interpretation of this property is that φ must shoot to infinity at least in one direction (R + or R -, at least linearly (first line of ( 4)), and also at most linearly (second line of ( 4)).

Of course, compactly supported nonlinearities such as sigmoid and tanh do not satisfy the extended envelope property but the majority of other nonlinearities do, including ReLU, ELU, SELU (see [START_REF] Klambauer | Selfnormalizing neural networks[END_REF] for details), and others.

Definition A.2 (Sub-Weibull random variable) A random variable X, that satisfies

for K > 0, is called a sub-Weibull random variable with the tail parameter θ, which is denoted by X ∼ subW(θ).

Informally, the tails of a subW(θ) distribution are dominated by (i.e. decay at least as fast as) the tails of a Weibull variable [START_REF] Rinne | The Weibull distribution: a handbook[END_REF]; in the same way as sub-Gaussian or sub-Exponential distributions correspond to distributions with tails lighter than Gaussian and Exponential distributions, respectively. Sub-Weibull distributions are parameterized by a positive tail index θ and equivalent to sub-Gaussian for θ = 1/2 and sub-Exponential for θ = 1. The larger tail parameter θ, the heavier the tails of the sub-Weibull distribution.


Mariia Vladimirova

Julyan Arbel

Pablo Mesejo

Bayesian neural network priors at the level of units

come

Introduction

We investigate deep Bayesian neural networks with Gaussian priors on the weights and ReLU-like nonlinearities. See [START_REF] Vladimirova | Bayesian neural networks increasingly sparsify their units with depth[END_REF]

Notations

Given an input x ∈ R N , the -th hidden layer unit activations are defined as

g ( ) (x) = W ( ) h ( -1) (x), h ( ) (x) = φ(g ( ) (x)).

Assumptions

• Gaussian prior on weights:

W ( ) i,j ∼ N (0, σ 2 w ),
• A nonlinearity φ : R → R is said to obey the extended envelope property if there exist

c 1 , c 2 , d 2 ≥ 0, d 1 > 0 such that |φ(u)| ≥ c 1 + d 1 |u| for u ∈ R + or u ∈ R -, |φ(u)| ≤ c 2 + d 2 |u| for u ∈ R.

Sub-Weibull

A random variable X, such that

P(|X| ≥ x) ≤ exp -x 1 /θ /K
for all x ≥ 0 and for some K > 0, is called a sub-Weibull random variable with tail parameter θ > 0:

X ∼ subW(θ).

Moment property:

X ∼ subW(θ) implies

X k = E|X| k 1 /k k θ ,
Meaning for all k ∈ N and for some constants d, D > 0,

d < X k /k θ < D.

Covariance theorem

The covariance between hidden units of the same layer is nonnegative. Moreover, for any -th hidden layer units h ( ) and h( ) , for

s, t ∈ N it holds Cov h ( ) s , h( ) t ≥ 0.

Penalized estimation

Regularized problem:

min W R(W ) + λL (W ), (1)
where R(W ) is a loss function, λL (W ) is a penalty, λ > 0.

For Bayesian models with prior distribution π(W ), the maximum a posteriori (MAP) solves (1) with: The -th hidden layer units U ( ) (pre-activation g ( ) or post-activation h ( ) ) of a feed-forward Bayesian neural network with:

L (W ) ∝ -log π(W )
• Gaussian priors on weights and

• extended envelope condition activation function φ have sub-Weibull marginal prior distribution with optimal tail parameter θ = /2, conditional on the input x:

U ( ) ∼ subW( /2),
Prior distributions of layers = 1, 2, 3

Illustration of units marginal prior distributions from the first three hidden layers. Neural network parameters: (N, H 1 , H 2 , H 3 ) = (50, 25, 24, 4).

Proof sketch

Induction w.r.t. layer depth :

h ( ) k k /2 ,
which is the moment characterization of sub-Weibull variable.

• Extended envelope property implies h ( )

k g ( ) k • Base step: g ∼ N (0, σ 2 ), g k √ k. Thus, h k = φ(g) k g k √ k. • Inductive step: suppose h ( -1) k k ( -1) /2 .
• Lower bound: non-negative covariance theorem:

Cov h ( -1) s , h( -1) t ≥ 0.

• Upper bound: Holder's inequality

• g ( ) = H j=1 W ( -1) i,j h ( -1) j implies h ( ) k g ( ) k k /2 .

Sparsity interpretation

MAP on weights is L2-reg. Independent Gaussian prior

π(W ) ∝ L =1 i,j e -1 2 (W ( ) i,j ) 2 ,
is equivalent to the weight decay penalty with negative log-prior:

L (W ) ∝ L =1 i,j (W ( ) i,j ) 2 = W 2 2 ,

MAP on units induces sparsity

The joint prior distribution for all the units can be expressed by Sklar's representation theorem as

π(U ) = L =1 H m=1 π ( ) m (U ( ) m ) C(F (U )),
where C is the copula of U (characterizes all the dependence between the units), F is its cumulative distribution function. The penalty is the negative log-prior:

L (U ) ≈ U (1) 2 2 + • • • + U (L) 2 /L 2 /L -log C(F (U )).
Layer W -penalty U -penalty

1 W (1) 2 2 , L 2 U (1) 2 2 , L 2 2 W (2) 2 2 , L 2 U (2) , L 1 W ( ) 2 2 , L 2 U ( ) 2/ 2/ , L 2/

Conclusion

We prove that the marginal prior unit distributions are heavier-tailed as depth increases. We further interpret this finding, showing that the units tend to be more sparsely represented as layers become deeper. This result provides new theoretical insight on deep Bayesian neural networks, underpinning their natural shrinkage properties and practical potential.

  .

	input	1 st hid.	2 nd hid.	th hid.
	layer	layer	layer	layer
	x	h (1)	h (2)	. . . h ( )
				. . .
	. . .	. . .	. . .	. . . . . .
		subW( 1 2 ) subW(1) subW( 2 )


Caroline Lawless

Julyan Arbel

Chinese restaurant process from stick-breaking for Pitman-Yor

Chinese restaurant process from stick-breaking for

Pitman-Yor

Introduction

• The Chinese restaurant process and the stick-breaking process are the two most commonly used representations of the Pitman-Yor process.

• However, the usual proof of the connection between them is indirect.

• [START_REF] Miller | An elementary derivation of the Chinese restaurant process from Sethuraman's stick-breaking process[END_REF] proved directly that the stick-breaking process gives rise to the Chinese restaurant process representation of the Dirichlet process.

• The Dirichlet process is a special case of the Pitman-Yor process.

• We extend Miller's proof to Pitman-Yor process random measures.

Pitman-Yor & Dirichlet processes

• The Dirichlet Process (DP) and the Pitman-Yor process (PY, Pitman and Yor, 1997) are discrete random probability measures.

• The PY is parametrized by d ∈ (0, 1), α > -d, and a base probability measure P 0 . The DP is recovered by letting d = 0.

• The stick-breaking representation [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF] is given by

v i ∼ Beta(1, α) for DP Beta(1 + d, α + id) for PY π k =v k k-1 i=1 (1 -v i ), φ k iid ∼ P 0 .
We define the random process P by

P = ∞ i=1 π k δ φ k .
• The Chinese restaurant process [START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF] is the distribution induced on random partitions C given by

P (C = C) =    α |C| Γ(α) Γ(n+α) c∈C Γ(| c |) for DP d t ( α d ) (t) (α) (n) t j=1 (1 -d) (|c j |-1) for PY.

Theorem

Suppose π follows the PY stick-breaking, and

z 1 , . . . , z n |π = π iid ∼ π, that is, P(z i = k | π) = π k ,
and C is the partition of [n] induced by z 1 , . . . , z n . Then C follows the PY Chineses restaurant process.

Technical Lemmas

Our proof relies on the following lemmas, which here we will state without proof. Let us abbreviate z = (z 1 , . . . , z n ). Given z ∈ N n , let C z denote the partition [n] induced by z. We define m(z) = max {z 1 , • • • , z n }, and g k (z) = #{i: z i ≥ k}. Lemma 1 For any z ∈ N n ,

P(z = z) = 1 (α) (n) c∈C z Γ(|c| + 1 -d) Γ(1 -d) m(z) k=1 α + (k -1)d g k (z) + α + (k -1)d .
Lemma 2 For any partition C of [n],

z∈N n 1(C z = C) m(z) k=1 α + (k -1)d g k (z) + α + (k -1)d = d t ( α d ) (t) c∈C (|c| -d)
.

Proof of Theorem

P(C = C) = z∈N n P(C = C|z = z)P(z = z) (a) = z∈N n 1(C z = C) 1 (α) (n) c∈C z Γ(|c| + 1 -d) Γ(1 -d) m(z) k=1 α + (k -1)d g k (z) + α + (k -1)d = 1 (α) (n) c∈C z Γ(|c| + 1 -d) Γ(1 -d) z∈N n 1(C z = C) m(z) k=1 α + (k -1)d g k (z) + α + (k -1)d (b) = 1 (α) (n) c∈C z Γ(|c| + 1 -d) Γ(1 -d) d t ( α d ) (t) c∈C (|c| -d) (c) = 1 (α) (n) c∈C (1 -d) (|c|-1) c∈C (|c| -d) d t ( α d ) (t) c∈C (|c| -d) = d t ( α d ) (t) (α) (n) t j=1 (1 -d) (|c j |-1)
where (a) is by Lemma 1, (b) is by Lemma 2, and (c) is since Γ

(|c|+1-d) = (|c| -d)Γ(|c -d|).

Further research

• The Dirichlet process and the Pitman-Yor process are only special cases of a broad class of random measures called Gibbs-type random measures.

• An interesting further study would be to investigate the possibility of extending this proof to Gibbs-type random measures.


Aleksandra Malkova

Julyan Arbel

Maria Laura

Delle Monache

DATASAFE: understanding Data Accidents for TrAffic SAFEty Acknowledgments

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Aim and data

Our aim is to understand if there exists a particular link between speed and density of cars at which collisions are more likely to occur. We intend to make use of raw sensor data collected on the Rocade Sud in the framework of the Grenoble Traffic Lab (https://gtl.inrialpes.fr/ data_download) using data for different days of the week and different day time. The sensor data consists of these quantities collected every 15 seconds: flow, occupancy and speed.

Future work

We will derive traffic density and flow rate from the sensor data by first reducing the data dimensionality (PCA). Then using state of the art Bayesian statistical learning, building upon [START_REF] Polson | Bayesian analysis of traffic flow on interstate I-55: The LWR model[END_REF], and clustering techniques we will analyze traffic behavior during periods of time that precede the accidents from which we will derive traffic patterns.

Fundamental diagram

The fundamental diagram is the main tool used in traffic flow to understand traffic characteristics. It is the graph that links the flow and the density. In this project, we study how the fundamental diagram depends on time and day of the week, weather conditions, characteristics of the road and local traffic rules. In particular, we are interested in understanding if certain particular road conditions can be linked to the generation of accidents.

Note that the two fundamental diagrams above are not radically different. In contrast, the accident is quite visible on the morning of the speed VS time plot.

Days with and without accident

Figure 1 .

 1 Figure 1. Raw speed data for one sensor.

Figure 2 .

 2 Figure 2. Distribution of the occupancy, speed and flow for one sensor.

Left:

  Fundamental diagrams. Right: Speed VS time. First row: Normal day without accident. Second row: Day with accident.They are different days, but same week day.

Figure 1 .

 1 Figure 1. Speed raw data represented according to day of the month VS time of the day, with color code from black-red for low speed to green for high speed.

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, M2SIAM, 38000 Grenoble, France

Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK,

Grenoble, France 3 Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, GIPSA-LAB, 38000 Grenoble, France {julyan.arbel, ml.dellemonache, aleksandra.malkova}@inria.fr

Acknowledgments

This project was developed in the framework of the Grenoble Alpes Data Institute, supported by the French National Research Agency under the "Investissements d'avenir" program (ANR-15-IDEX-02).

DATASAFE: understanding Data Accidents for TrAffic SAFEty


Julyan Arbel

Bayesian Nonparametric Mixtures Why and How?

published or not. The documents may come

Bayesian Nonparametric Mixtures Why and How?

Introduction

• Discrete random probability measure prior [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF] Predictive: Chinese Restaurant Process

• Random weights (p i ) i and locations (θ i ) i P = ∞ i=1 p i δ θ i → Dirichlet process DP(α, G 0 )
P(X n+1 ∈ • | X n ) = α α + n G 0 + 1 α + n k n j=1 n j δ X * j → Or for varying P(X n+1 new ) . . .

Density Estimation

Ecological risk assessment [START_REF] Arbel | Bayesian nonparametric modelling of species sampling distributions[END_REF] • Data are species critical effect concentrations (CEC), possibly censored

• Estimation of species sensitivity distribution (SSD), the density of CEC

• Safe concentration which protects most of the species: 5th percentile of the SSD (HC 5 )

• Very moderate sample sizes, ∼ 10 -50

→

Species Modeling

Data can be species, microbes, words, genes. . . (Arbel et al., 2016a) • Estimation of -discovery D = P(X n+1 is a species seen times)

Discovery probabilities

• Comparison with Good-Turing estimator → Closed form posterior and estimators → Uncertainty quantif., unavailable for GT → 2nd order (fast) approximations Diversity in ecology [START_REF] Arbel | Application of a Bayesian nonparametric model to derive toxicity estimates based on the response of Antarctic microbial communities to fuel contaminated soil[END_REF][START_REF] Arbel | Bayesian nonparametric dependent model for partially replicated data: the influence of fuel spills on species diversity[END_REF] • Assess impact of pollution on microbial community via study of diversity Div =i p i log p i → Model detects an hormetic effect → Uncertainty quantification → Prediction across full range of covariates q q q q q q q q q q q q q q q q q q q q q q Fuel concentration (g TPH/kg soil)

•

  Inference on curves: pdf, cdf, hazard, link. . .• Mixtures, exchangeable data X n = (X 1 , . . . , X n ) X 1 , . . . X n | P ∼ P?• Learn about data through posterior dist.

→

  BNP describes well variability of the data, without being prone to over-fitting → Species clustering as an outcome[START_REF] Arbel | Full Bayesian inference with hazard mixture models[END_REF] • Data are (remission) times possibly censored• Prior on hazard rate h(t) for every time t• Induces prior on survival function S (t)→ Availability of post. mean, median, mode → Smooth estimator VS Kaplan-Meyer Proper

Funding

This work is supported by the European Research Council (ERC) through StG "N-BNP" 306406.

Open Questions

• How to best use underlying clustering? [START_REF] Wade | Bayesian cluster analysis: Point estimation and credible balls[END_REF] • Find consistent estimator of number of clusters: posterior inconsistent [START_REF] Miller | Inconsistency of Pitman-Yor process mixtures for the number of components[END_REF], what about posterior mode?

• Devise efficient posterior sampling, truncation error [START_REF] Arbel | A moment-matching Ferguson & Klass algorithm[END_REF]


Olivier Marchal

Julyan Arbel

Beta and Dirichlet sub-Gaussianity

Optimal proxy variance σ 2 opt for the sub-Gaussianity of Beta distribution, improves recent conjecture σ 2 0 by Elder (2016): σ 2 opt ≤ σ 2 0 = 1 4(α+β+1) .

• Provide different proof techniques for (i) symmetrical case (α = β): direct coef. comparison of entire series (ii) non-symmetrical case (α β): ordinary differential equation satisfied by moment-generating function, aka confluent hypergeometric function

• Derive optimal proxy variance for Dirichlet.

Sub-Gaussianity & optimal proxy variance

• A random variable X with finite mean µ = E[X] is sub-Gaussian if there is a positive number σ 2 (called a proxy variance) such that:

E[exp(x(X -µ))] ≤ exp x 2 σ 2 2 for all x ∈ R. (1)
• If X is sub-Gaussian, the optimal proxy variance is:

σ 2 opt (X) = min{σ 2 ≥ 0 such that X is σ 2 -sub-Gaussian}.
Variance lower bounds optimal proxy variance: Var

[X] ≤ σ 2 opt (X). When σ 2 opt (X) = Var[X],
X is said to be strictly sub-Gaussian. • Top: varying α on x-axis, α + β = 1 (left) and α + β = 10 (right).

Optimal proxy variance illustrated & compared

• Left: α, β ∈ [0.2, 4]. • Magenta: σ 2 opt (α, β) • Green: Var[Beta(α, β)] • Black (dots): σ 2 0 = 1 4(α+β+1)
Looking for connections with literature

Optimal proxy variance for the Beta

Theorem Beta(α, β) is σ 2 opt (α, β)-sub-Gaussian with:

                   σ 2 opt (α, β) = α (α + β)x 0 1 F 1 (α + 1; α + β + 1; x 0 ) 1 F 1 (α; α + β; x 0 ) -1
where x 0 is the unique solution of the equation

log( 1 F 1 (α; α + β; x 0 )) = αx 0 2(α + β) 1 + 1 F 1 (α + 1; α + β + 1; x 0 ) 1 F 1 (α; α + β; x 0 ) . Simple explicit upper bound to σ 2 opt (α, β) is σ 2 0 (α, β) = 1 4(α+β+1) : Var[Beta(α, β)] ≤ σ 2 opt (α, β) ≤ σ 2 0 (strict when α β)

Sketch of proof

• α = β: direct coef. comparison of entire series representations of (1)

• α β: study the MGF aka confluent hypergeometric function or Kummer's function

y(x) def = E[exp(xX)] = 1 F 1 (α; α + β; λ) = ∞ j=0 Γ(α + j)Γ(α + β) ( j!)Γ(α)Γ(α + β + j) x j .
using the ordinary differential equation it satisfies

xy (x) + (α + β -x)y (x) -αy(x) = 0 via the difference u t (x) def = exp µx + σ 2 t x 2 /2 -E[exp(xX)], σ 2 t = t Var[X] + (1 -t)σ 2 0 0 1 2 3 -1 0 1 2 x10 -3 q x 0 u 0 u t opt u t non opt u 1
• For t = 0 (σ 2 0 ), dotted black curve remains > 0 • For t = t opt (σ 2 opt ), magenta curve has minimum = 0 at x 0 • For t = 1 (Var[X]), dashed green curve has negative second derivative at x = 0, directly negative around 0

• For t non opt ∈ (t opt , 1), orange, dash and dots curve is first positive, then negative, and positive again


Hongliang Lü

Julyan Arbel

Florence Forbes

Bayesian Nonparametric Priors for Hidden Markov Random Fields: Application to Image Segmentation

I. BACKGROUND AND MOTIVATION

Figure 1: Illustration of challenges for unsupervised image segmentation: blur, noise, color/contrast imperfection, partial volume effect (large slice thickness), anatomic variability and complexity, number of segments... Image segmentation in real-world applications is typically performed on noisy images. To achieve better segmentation performance, several extensions of the usual Bayesian nonparametric (BNP) mixture model with spatial regularization are therefore necessary.

II. BNP PRIORS

The Dirichlet process (DP) is one of the most commonly used BNP priors. It is a random process G defined over a probability space Y and characterized by a concentration parameter α and a base distribution G 0 , such that for any finite partition {A 1 , . . . , A p } of Y, the random vector (P (A 1 ), . . . , P (A p )) is Dirichlet distributed:

(P (A 1 ), . . . , P (A p )) ∼ Dir(αG 0 (A 1 ), . . . , αG 0 (A p )) which is often denoted by G ∼ DP(α, G 0 ).

III. STICK-BREAKING CONSTRUCTION

The DP has almost surely discrete realizations. It can be built by the stick-breaking construction:

G = ∞ k=1 π k (τ ) δ θ * k = ∞ k=1 τ k l<k (1 -τ l ) δ θ * k where θ * k iid ∼ G 0 and τ k iid ∼ B(1, α).

IV. DP-POTTS MIXTURE MODEL

The usual DP mixture model assumes that a set of data points y = {y 1 , . . . , y N } with y i ∈ R D (e.g., pixels) can be generated through the following hierarchical representation:

• G ∼ DP(α, G 0 ) • θ i |G ∼ G, i = 1, . . . , N • y i |θ i ∼ F (y i |θ i ), i = 1, . . . , N
where θ = {θ 1 , ..., θ N } denotes a set of model parameters. To take into account spatial constraints, we introduce a Potts model component using a set of assignment variables z = {z 1 , . . . , z N } with z i = z(θ i ) so as to favor spatial agregation [1]:

M (θ) ∝ exp   β i∼j δ z(θ i )=z(θ j )  
with β being the regularization parameter. The DP mixture model is thus extended to become the DP-Potts mixture model:

• G ∼ DP(α, G 0 ) • θ|M , G ∼ M (θ) × i G(θ i ) • y i |θ i ∼ F (y i |θ i ), i = 1, . . . , N
Accordingly, the stick-breaking construction of the DP-Potts mixture model can be summarized as follows:

• θ * k |G 0 ∼ G 0 and τ k |α ∼ B(1, α), k = 1, 2, . . . • π k (τ ) = τ k k-1 l=1 (1 -τ l ), k = 1, 2, . . . • p(z|τ ; β) ∝ i π z i (τ ) exp(β i∼j δ z i =z j ), z i = 1, 2, . . . • y i |z i , θ * ∼ F (y i |θ * z i ), i = 1, . . . , N

V. VARIATIONAL BAYES

In a Bayesian setting, we need to evaluate the intractable posterior distribution p(z, τ , α, θ * |y; φ) (φ denotes a set of hyperparameters) which can be estimated by means of the mean-field approximation:

q(z, τ , α, θ * ) q z (z)q τ (τ )q α (α)q θ * (θ * )

Variational Bayes (VB) consists of alternating maximization of free energy F(q z , q τ , q α , q θ * ; φ) = E q z q τ q α q θ * log p(z, τ , α, θ * , y; φ) q z q τ q α q θ * which implies [2] • E-steps: VE-z, VE-α, VE-τ and VE-θ * .

• M-steps: φ updating is straightforward except for β.

Here, M-β step leads to the estimation of β:

β = arg max β E q z q τ [log p(z|τ ; β)]
which involves p(z|τ ; β) = K(β, τ ) -1 exp (V (z; τ , β)) with the normalization constant K(β, τ ) and the potential function

V (z; τ , β) = i log π z i (τ ) + β i∼j δ z i =z j
To find the optimal value of β, further approximations, such as the mean-field-like approximation [3] of q z and replacing τ with a fixed τ = E q τ [τ ], are required.

VI SOME EXPERIMENTS AND RESULTS

Original image DP mixture model (about 1000 superpixels) DP-Potts mixture model (about 1000 superpixels)

Experiments were performed using superpixels on a subset (154 images) of the Berkeley segmentation data set (BSDS) [4]. Regarding the performance evaluation, the probabilistic rand index (PRI) was computed under different conditions: We also compared our best results (K = 50 and about 500 superpixels) with those given in the literature:

Proposed models

Results given in [

VII CONCLUSIONS AND FUTURE WORK

A general DP-Potts mixture model and the associated VB algorithm were proposed. The model was successfully applied to image segmentation on different types of datasets. We also investigated the impact of β on the segmentation results and presented an estimation procedure for β.

In the sequel, we plan to survey how β affects the inferred number of components.

Other types of priors (Pitman-Yor process, Gibbs-type priors, etc.) and other variational approximations (truncation-free) will also be considered. On the other hand, it is crucial to study theoretical properties of BNP priors under structural constraints (temporal or spatial). Other applications may also be possible, such as discovery probability and community detection in graphs.

  and DP-Potts mixture models on BSDS DP-Potts (500 superpixels) DP (500 superpixels) DP-Potts (1000 superpixels) DP (1000 superpixels)


Sara Ifqir

Vicenc ¸puig

Naima Ait-Oufroukh

Dalil Ichalal

Saïd Mammar

Robust Fault Detection for Vehicle Lateral Dynamics: A Zonotope-based Set-membership Approach

Keywords: Robust fault detection, interval models, zonotopes, set-membership, switched uncertain systems, LMIs, input-to-state stability, arbitrary switching

In this work, a model-based fault detection layout for vehicle lateral dynamics system is presented. The major focus in this study is on the handling of model uncertainties and unknown inputs. In fact, the vehicle lateral model is affected by several parameter variations such as longitudinal velocity, cornering stiffnesses coefficients and unknown inputs like wind gust disturbances. Cornering stiffness parameters variation is considered to be unknown but bounded with known compact set. Their effect is addressed by generating intervals for the residuals based on the zonotope representation of all possible values. The developed fault detection procedure has been tested using real driving data acquired from a prototype vehicle.

I. INTRODUCTION

Detection of malfunctions and faults in vehicle lateral dynamics is one of the most important tasks allowing to improve the performance and reliability of active safety systems. An early warning indicating the occurrence of faults can help to take appropriate actions and avoid damage or even fatal accidents. Therefore, it is necessary to develop a Fault Detection (FD) procedure to be able to detect potential faults as early as possible. Then, the system can maintain stability and keep acceptable performances by means of this early warning including some fault tolerant mechanism. Among numerous FD techniques, model-based design has attracted an increasing amount of attention over the last decades, and some approaches have been proposed (see, e.g. [START_REF] Chen | Design of unknown input observers and robust fault detection filters[END_REF], [START_REF] Ingimundarson | Robust fault detection using zonotope-based setmembership consistency test[END_REF], [START_REF] Nejjari | Robust fault detection for lpv systems using interval observers and zonotopes[END_REF], [START_REF] Varrier | Robust fault detection for uncertain unknown inputs lpv system[END_REF], [START_REF] Zhong | An lmi approach to design robust fault detection filter for uncertain lti systems[END_REF] and references therein). They perform fault detection by generating a residual signal based on the difference between the real and mathematical model outputs, which in the fault-free case should be zero, otherwise, a fault is detected and an alarm is generated. However, since there exist some mismatches between mathematical model and real plant due to parameter uncertainties, external disturbances and unknown inputs, the residual signal becomes nonzero, which generates false alarms, corrupting the results and making the FD system useless. Hence, robustness against these uncertainties and disturbances is an important and challenging problem in model-based FD design. In the literature, two main distinct approaches have been addressed. The first robust approach, known as active, is {sara.ifqir, naima.aitoufroukh, dalil.ichalal, said.mammar}@univ-evry.fr 2 V. Puig is with Automatic Control Department, Universitat Politècnica de Catalunya, Barcelona, Spain. vicenc.puig@upc.edu based on generating residuals that are sensitive to faults, while being insensitive to uncertainties. This approach has been extensively developed using Unknown Input Observer (UIO) [START_REF] Varrier | Robust fault detection for uncertain unknown inputs lpv system[END_REF], parity-space and identification based methods [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF]. Whereas, the second approach, known as passive, is based on enhancing residual robustness at the decisionmaking stage using an adaptive threshold [START_REF] Nejjari | Robust fault detection for lpv systems using interval observers and zonotopes[END_REF]. In other words, the passive approach consists on propagating the uncertainties to the residuals limits. When the residual is outside these limits, the alarm is not due to uncertainties but to a fault occurrence. In the Passive Fault Detection (PFD) technique, uncertain parameters are bounded and the resulting model known as interval model is used for design. In this case, the measured output is bounded by an interval, providing, therefore, an interval for the residual vector. Then, fault detection test consists in checking if the residual value is inside or not this interval after uncertainty propagation. In this paper, zonotopes are used to provide the interval enclosing the unknown true residual vector. The advantage of using zonotopes is that monotonicity property is no longer needed as in [START_REF] Lamouchi | Interval observer design for linear parameter-varying systems subject to component faults[END_REF], [START_REF] Lamouchi | Actuator fault compensation in a set-membership framework for linear parameter-varying systems[END_REF], [START_REF] Raïssi | Interval observer design for consistency checks of nonlinear continuous-time systems[END_REF] allowing to consider the whole set of possible trajectories. The zonotope-based approaches in fault detection has already been investigated in the literature, see for instance [START_REF] Guerra | Robust fault detection with unknown-input interval observers using zonotopes[END_REF], [START_REF] Ingimundarson | Robust fault detection using zonotope-based setmembership consistency test[END_REF]. The main contribution of this paper is that it provides a robust PFD schema based on the set-membership approach using zonotopes to bound the parameter uncertainty affecting the vehicle lateral dynamics. For this purpose, it is important to use a detailed and accurate vehicle model. In this study and in order to fully make use of the range variations of all vehicle uncertain parameters, such as, vehicle longitudinal velocity, cornering stiffness coefficients, and to take into account external disturbances such as lateral wind gusts, a switched uncertain discrete-time system subject to unknown inputs is proposed. In this model, the longitudinal velocity is considered piecewise constant, different local sub-model are then derived and a switching strategy dependent on the measured velocity is implemented. Furthermore, an adaptive dynamic tire model that is physical and flexible enough to permit timevarying tire performances is used. Indeed, a readjustment variable assumed to be unknown but bounded is added to correct the value of the cornering stiffness coefficients in case that a change occur in road adhesion or vehicle loading condition. Using Switched Unknown Input Observer (SUIO) principle and Input-to-State Stability concept, the residual vector is made robust against unknown lateral wind gusts disturbances and aforementioned uncertainties. This work shows the successful application of PFD method in order to detect faults in vehicle lateral dynamics. It should be pointed out that the use of adaptive threshold in vehicle lateral diagnosis systems is an area that is still fairly undeveloped. Moreover, for the best of the author knowledge, the PFD method has not been intensively investigated for switched systems. This paper, presents perhaps the first attempt to deal with such systems. The remainder of this paper is organized as follows: Section II presents the vehicle lateral model under consideration. The design of the PFD schema is proposed in Section III. Application to vehicle lateral dynamics with validation through real data is given in Section IV. The main conclusions and future works are drawn in Section V.

II. VEHICLE LATERAL DYNAMICS MODELING

In the literature, there are several mathematical models describing the vehicle lateral dynamics. In the present work, the well known two degree of freedom bicycle model representing the lateral and yaw motions is used for implementing FD procedure. The vehicle mass is equal to m. The moment of inertia is noted I z . The front and rear axles are located respectively at distances l f and l r from the center of gravity. The front and rear lateral tire forces F yf and F yr are functions of wheel slip angles α f and α r , respectively. The force F w is due to the effect of the wind gusts that are acting at a distance l w . v y and v x are lateral and longitudinal velocity respectively, and δ f stands for the front steering angle. The vehicle lateral dynamics is described by the following differential equations with lateral velocity v y and yaw rate ψ as states [START_REF] Rajamani | Vehicle dynamics and control[END_REF]:

m vy = F yf + Fyr -mvxr + Fw Iz ψ = l f F yf -lrFyr + lwFw (1)
For normal driving situations, F yf and F yr are linear with respect to wheel slip angles

F yf = c f α f , Fyr = crαr (2)
where c f and c r represent the front and rear wheel cornering stiffness parameters respectively. The front and rear wheel slip angles are defined as

α f = δ f - vy vx - l f vx ψ, αr = - vy vx + lr vx ψ (3)
Substituting ( 2) and ( 3) into (1), equations (1) are expressed in the following state-space form

ẋ(t) = Ax(t) + Bu(t) + Ed(t) y(t) = Cx(t) (4)
where

x(t) = vy ψ T , u(t) = δ f , d(t) = F w , y(t) = ψ A = - c f +cr mvx cr lr -c f l f mvx -vx cr lr -c f l f Iz vx - cr l 2 r +c f l 2 f Iz vx , B = c f m c f l f Iz E = 1 m lw Iz
Note that, the tire model ( 2) is valid only in the case of low lateral acceleration. However, when lateral acceleration becomes high, the tire forces are no longer linear with respect to the wheel slip angles due to the tire saturation property. Accordingly, the cornering stiffness parameters c f and c r used in the linear tire model (2) will vary when the road friction changes or when the nonlinear tire domain is reached. Taking these variations into account, the following linear adaptive tire model is used to correct the cornering stiffness values via the uncertain variables ∆c f and ∆c r as

F yf = (c f 0 + ∆c f )α f , Fyr = (cr0 + ∆cr)αr (5)
where c i0 , i ∈ {r, f }, represents the known nominal value and ∆c i , i ∈ {r, f }, is assumed to be unknown but bounded with a priori known bounds. Furthermore, it is clear that the uncertain bicycle model ( 4) is non-linear because of the term1 vx . In this paper, the longitudinal velocity is considered piecewise constant and a switching rule depending on the measured velocity is deduced to select the active subsystem. Then, the uncertain non-linear bicycle model ( 4) is transformed into a switched uncertain discrete-time model 1 :

   x k+1 = (A 0,σ(k) + ∆A σ(k) (θ k ))x k + (B 0,σ(k) + ∆B σ(k) (θ k ))u k + E σ(k) d k y k = Cx k (6)
where

x k ∈ R n , u k ∈ R m , y k ∈ R p , d k ∈ R q
are the state, the control input, the output vector and the unknown input respectively. σ : R + → I = {1, 2, ..., N } is the switching signal assumed to be available in real time. The matrices A 0,σ(k) ∈ {A 0,1 , A 0,2 , ..., A 0,N } and B 0,σ(k) ∈ {B 0,1 , B 0,2 , ..., B 0,N } are the state space matrices assumed to be constant and known a priori. E σ(k) and C stand for the unknown input distribution matrix and the output equation matrix, respectively. The matrices

∆A σ(k) (θ k ) ∈ {∆A 1 (θ k ), ∆A 2 (θ k ), ..., ∆A N (θ k )} and ∆B σ(k) ∈ {∆B 1 (θ k ), ∆B 2 (θ k ), .
.., ∆B N (θ k )} are assumed to be unknown but bounded representing the uncertainty of the system caused by the time-varying parameter θ k = ∆c f ∆c r T ∈ Θ, where Θ is an interval box given by:

Θ = {θ k ∈ R r |θ - k ≤ θ k ≤ θ + k } (7)
For easiness of further developments, the state space representation ( 6) is rewritten in the following equivalent form:

x k+1 = A 0,σ(k) x k + B 0,σ(k) u k + E σ(k) d k + δ σ(k) y(t) = Cx(t) (8)
where

δ σ(k) = ∆A σ(k) (θ k )x k + ∆B σ(k) (θ k )u k .

III. VEHICLE FAULT DETECTION

A. UI decoupling

For the system (8), a switched observer that reconstructs the state x k without any information on the UI d k can be constructed using the input u k and the measured output y k as

xk+1 = N σ(k) xk + K σ(k) y k + G σ(k) u k -H σ(k) y k+1 ŷk = C xk (9)
with the initial state x0 ∈ R n . N σ(k) , K σ(k) , G σ(k) and H σ(k) are gain matrices with appropriate dimensions to be determined. Theorem 1. The switched observer ( 9) is a switched UIO for system [START_REF] Lamouchi | Interval observer design for linear parameter-varying systems subject to component faults[END_REF] if and only if the following conditions are satisfied:

1)

e k+1 = N σ(k) e k + ∆ σ(k) is ISS-stable under arbitrary switching; 2) N σ(k) = P σ(k) A 0,σ(k) -K σ(k) C; 3) P σ(k) B 0,σ(k) -G σ(k) = 0; 4) P σ(k) E σ(k) = 0. where P σ(k) = I n + H σ(k) C and ∆ σ(k) = P σ(k) δ σ(k) .
Proof. The dynamics of the state estimation error e k = x k -xk is given by

e k+1 = x k+1 -xk+1 = (In + H σ(k) C)x k+1 -N σ(k) xk -K σ(k) y k -G σ(k) u k By denoting P σ(k) = I n + H σ(k)
C and using ( 6), the dynamics of the state estimation error is

e k+1 = (P σ(k) A 0,σ(k) -K σ(k) C)x k -N σ(k) xk + (P σ(k) B 0,σ(k) -G σ(k) )u k + P σ(k) E σ(k) d k + P σ(k) δ σ(k) (10) If the following relationships hold N σ(k) = P σ(k) A 0,σ(k) -K σ(k) C (11a) P σ(k) B 0,σ(k) -G σ(k) = 0 (11b) P σ(k) E σ(k) = 0 (11c)
then, ( 10) is reduced to

e k+1 = N σ(k) e k + ∆ σ(k) (12)
which completes the proof. From [START_REF] Rajamani | Vehicle dynamics and control[END_REF], one can see that d k has been decoupled under the conditions 2-4 of Theorem 1, but the effect of uncertainties ∆ σ(k) still persists. Therefore, the switched observer design consists on finding the observer gain K σ(k) such that the observer error dynamics ( 12) is asymptotically stable when ∆ σ(k) = 0 and is ISS-stable with respect to uncertainties when ∆ σ(k) = 0.

B. LMI formulation

The Lemma 1 is useful for the proof of ISS-stability of the SUIO (9). Lemma 1. The necessary and sufficient conditions for the existence of the SUIO (9) for the system (8) are:

1) rank(CE σ(k) ) = rank(E σ(k) ) = q , ∀σ(k);
2) The pair (P σ(k) A 0,σ(k) , C) is detectable ∀σ(k). Remark 1. The first condition in Lemma 1 ensures that equation (11c) is solvable, and a particular solution can be given as follows:

H σ(k) = -E σ(k) (CE σ(k) ) T (CE σ(k) ) -1 (CE σ(k) ) T (13)
Moreover, the second condition, sometimes referred as "strong detectability condition" is equivalent to that the transmission zeros from the UI to the output must be stable, i.e.

rank sIn -A 0,σ(t) E σ(k) C 0 = n + q, (14)
holds ∀σ(k) and for all complex number s with Re(s) ≥ 0.

The next theorem provides sufficient conditions for the state estimation error [START_REF] Rajamani | Vehicle dynamics and control[END_REF] to be ISS-stable under arbitrary switching. Before proceeding, we begin by stating the following equivalence. Lemma 2. The following conditions are equivalent:

(E 1 ) There exists a symmetric matrix Q such that

N T QN + (α -1)Q N T Q QN Q -γI ≺ 0 (15)
(E 2 ) There exist a symmetric matrix Q and a matrix S such that

  (α -1)Q 0 N T S T ( * ) -γIn S T ( * ) ( * ) -S T -S + Q   ≺ 0 (16)
Proof. If we apply the Schur complement with respect to the block (3,3) of ( 16), we recover directly (15) by choosing S = S T = Q, hence (E 1 ) implies (E 2 ). Moreover, from the first block of (16), we have (α -1)Q ≺ 0. Then, multiplying (16) by T = I n 0 N T 0 I n I n on the left and T T on the right, we get [START_REF] Zhong | An lmi approach to design robust fault detection filter for uncertain lti systems[END_REF], which establishes that (E 2 ) implies (E 1 ) and the proof is complete. Theorem 2. System ( 12) is uniformly2 ISS-stable with respect to the switching signal σ(k) if there exist a symmetric matrix Q, matrices S and W i , ∀i ∈ I, a constant γ > 0 for given positive scalars α 1 and 0 < α < 1 such that the following condition holds

min Q,S,Wi γ   (α -1)Q 0 A0,i T P T S T -C T W T i ( * ) -γIn S T ( * ) ( * ) -S T -S + Q   ≺ 0 (17)
with W i = SK i for any switching signal σ(k). Furthermore the state estimation error [START_REF] Rajamani | Vehicle dynamics and control[END_REF] satisfies

lim k→∞ e k 2 ≤ γ α1(1 -α) ∆ σ(k) ∞ (18)
where ∆ σ(k) = P δ σ(k) and the maximum norm is given by

∆ σ(k) ∞ = max{|∆ 1 |, |∆ 2 |, . . . , |∆ N |}. Proof.
For the stability analysis, we use the following common quadratic ISS-Lyapunov function

V (e k ) = e T k Qe k , Q = Q T 0 ( 19
)
Setting ∆V (e k ) V (e k+1 ) -V (e k ), we have

∆V (e k ) = e T k+1 Qe k+1 -e T k Qe k using (12), we get

∆V (e k ) = N σ(k) e k + ∆ σ(k) T Q N σ(k) e k + ∆ σ(k) -e T k Qe k = e T k N T σ(k) QN σ(k) e k + ∆ T σ(k) QN σ(k) e k +e T k N T σ(k) Q∆ σ(k) + ∆ T σ(k) Q∆ σ(k) -e T k Qe k
By adding and subtracting the terms αe T k Qe k and -γ∆ T σ(k) ∆ σ(k) , we obtain that

∆V (e k ) = e T k ∆ T σ(k) N T σ(k) QN σ(k) + (α -1)Q N T σ(k) Q ( * ) Q -γIn × e k ∆ σ(k) -αe T k Qe k + γ∆ T σ(k) ∆ σ(k) (20)
Then, the satisfaction

N T σ(k) QN σ(k) + (α -1)Q N T σ(k) Q ( * ) Q -γIn ≺ 0 (21)
is equivalent to satisfy (17) according to Lemma 2. Then, from (20) we get

∆V (e k ) < -αe T k Qe k + γ∆ T σ(k) ∆ σ(k) (22)
that implies

V (e k+1 ) < (1 -α)V (e k ) + γ ∆ σ(k) 2 (23)
Integrating (23) over the interval [k 0 , k), we have

V (e k ) ≤ (1-α) k-k 0 V (e0)+γ k-1 m=0 (1-α) k-m-1 ∆ σ(m) 2 (24)
knowing that (19) satisfies for some α 2 > α 1 > 0 the following inequality

α1 e k 2 ≤ V (e k ) ≤ α2 e k 2 (25)
allows to deduce that

e k 2 ≤ 1 √ α 1 α 2 (1 -α) k-k0 e 0 2 2 + γ 1 -α k-1 m=0 [1 -α] k-m ∆ σ(m) 2 2 1 2 (26)
Hence, when k → ∞ and using the maximum norm of ∆ σ(k) , ( 18) is obtained which ends the proof. Remark 2. The uniform ISS stability is usually used when there is no restriction on the switching signal σ(k) and requires that all the subsystems are ISS-stable. This condition is ensured by the choice of the gain K σ(k) . Note that the subsystems stability assumption is not sufficient to guarantee the stability of the overall switched system. However, the existence of the common quadratic ISS-Lyapunov function ( 19) is necessary and sufficient for the ISS-Stability of a the switched system (12) [START_REF] Rn Shorten | On the stability and existence of common lyapunov functions for stable linear switching systems[END_REF]. Remark 3. The existence of a common Lyapunov function for the switched system ( 12) is very conservative. The vehicle lateral dynamics system under consideration possesses a common Lyapunov function as will be shown in simulation, and remains stable under arbitrary switching signals. But despite this fact and in order to derive less conservative conditions our attention will be paid in future works on a less conservative class of Lyapunov functions, called multiple Lyapunov functions.

C. Fault detection using interval model

In the sequel, we are interested in actuator fault detection in vehicle lateral dynamics. To do this, the proposed switched UIO scheme is employed for actuator FD as described in the following. The fault detection procedure is based on testing whether the vehicle measured output y k is consistent with its estimation ŷk provided by the UIO [START_REF] Lamouchi | Actuator fault compensation in a set-membership framework for linear parameter-varying systems[END_REF]. When an inconsistency is detected, the existence of a fault is indicated. This consistency check is based on generating a fault indicator, known as residual r k and given by

r k = y k -ŷk (27)
which can be also expressed using state estimation error e k as follows

r k = Ce k (28)
Note that the residual vector ( 28) is expected to be zero in fault-free case and deviate from zero in the presence of faults [START_REF] Akhenak | Design of observers for takagi-sugeno fuzzy models for fault detection and isolation[END_REF], [START_REF] Chen | Robust model-based fault diagnosis for dynamic systems[END_REF]. However, due to the presence of parameter uncertainties ∆ σ(k) , the residual vector r k will be different from zero, even in a non-faulty situation. We propose in this paper to use the so called passive fault detection method, that enhances the robustness of the fault detection schema at the decision making step. The initial state uncertainties and system parameter uncertainties are propagated to the residual vector, allowing to determine residual reachable interval. When the residual is outside the given interval, the occurrence of a fault is proved, i.e. while the residual

r k ∈ R k = [r - k , r + k ] (29)
no fault occurred, and the residual value is simply due to the parameter uncertainty. r -and r + are lower and upper bounds of the residual r k , respectively. The set R k is obtained in real time using zonotopes.

D. Implementation of the FD schema using zonotopes

In this paper, zonotopes are used to implement the fault detection algorithm. The following definitions, assumptions and Lemma are introduced in advance. Definition 1. The Minkowski sum of two sets X 1 and X 2 is given by

X 1 ⊕ X 2 = {x 1 + x 2 : x 1 ∈ X 1 , x 2 ∈ X 2 }. Definition 2.
[10] Given a center vector x c ∈ R n and a matrix H ∈ R n×m , the Minkowski sum of the segments defined by the columns of matrix H is called zonotope of order m. This set is defined by:

X = x c ⊕ Hβ m = {x c + Hz : z ∈ β m } (30)
where β m is an unitary box composed by m unitary intervals. Definition 3. The interval hull X 1 of a closed set X 1 is the smallest interval box that contains X 1 . Assumption 1. The initial state x 0 is inside an initial zonotope X 0 , implying that the initial state estimation error e 0 is inside the zonotope E 0 = CX 0 , i.e. e 0 ∈ E 0 . Assumption 2. It is assumed that the state and input vectors

x k and u k are bounded as follows

x - k ≤ x k ≤ x + k , u - k ≤ u k ≤ u + k , where, x - k , x + k , u - k and u +
k are a priori known bounds reflecting the vehicle physical limits.

E k-1 E k e k+1 = N σ(k) e k + ∆ σ(k) r k = Ce k R k
+ k , x - k ∈ R n such that x - k ≤ x k ≤ x + k 1) If M ∈ R n×n is a constant matrix, then M x - k -M x + k ≤ M x k ≤ M x + k -M x - k (31) 2) If M ∈ R n×n is a variable such that M -≤ M ≤ M + for some M -, M + ∈ R n×n , then M -x k --M + x k --M -x k + + M + x k + ≤ M x k ≤ M + x k + -M -x k + -M + x k -+ M -x k - (32)
where = max(0, ) and = -. Note that since θ k ∈ Θ, the system matrices ∆A σ(k) (θ k ) and ∆B σ(k) (θ k ) are bounded based on [START_REF] Ingimundarson | Robust fault detection using zonotope-based setmembership consistency test[END_REF] as

∆A - σ(k) ≤ ∆A σ(k) (θ k ) ≤ ∆A + σ(k) (33a)
∆B - σ(k) ≤ ∆B σ(k) (θ k ) ≤ ∆B + σ(k) (33b)
Therefore, using Assumption 2 together with (33) and applying Lemma 3, the additive term ∆ σ(k) in ( 12) can be bounded by means of an interval box D given by:

D = {∆ σ(k) ∈ R n | ∆ - σ(k) ≤ ∆ σ(k) ≤ ∆ + σ(k) , ∆ - σ(k) , ∆ + σ(k) ∈ R n } (34) Furthermore, D can be rewritten as a zonotope D = ∆ c ⊕ H ∆ + σ(k) β n , where H ∆ + σ(k)
∈ R n×n is a diagonal matrix with the main diagonal being ∆ + σ(k) and ∆ c is a known and constant vector. Then, using Definitions 2 and 3, the state estimation error uncertainty is bounded using zonotopes and ( 12) is transformed to the following equivalent form

e c k+1 = N σ(k) e c k (35a) H k+1 = N σ(k) H k H ∆ + σ(k) (35b)
where e c k+1 and H k+1 are the center and segment matrix of E k+1 , respectively. The equivalent compact description of (35) is given as follows

E k+1 = N σ(k) E k ⊕ D (36)
Thus, the residual zonotope at instant k can be written as

R k = CE k (37)
by using the iterative algorithm presented in Figure 1. The zonotope R k is obtained for time instants k > 0. Subsequently, the interval for residual vector R k is obtained by computing the interval hull of the zonotope (37). By means of testing the robust FD criterion (29), the occurrence or absence of faults is detected in real time.

IV. SIMULATION RESULTS USING REAL DATA

In this section, the robust fault detection scheme proposed in Section III is applied to the actuator fault detection problem of vehicle lateral dynamic system described in Section II. The real data used in the validation process are acquired using a prototype vehicle. The run was performed on a test track located in the city of Versailles-Satory (France). The track is 3.5km length with various curve profiles allowing vehicle dynamics excitation. Several sensors are implemented on the vehicle. The yaw rate r is measured using an inertial unit. An odometer is used to provide the vehicle longitudinal speed while the steering angle δ f is measured by an absolute optical encoder. A wind lateral force acting as unknown input is added in simulation and presented in Figure 2. The steering angle, the yaw rate and longitudinal velocity profiles are plotted in Figure 3. Measured vehicle specifications for the simulation and experimental validation can be found in [START_REF] Ifqir | Robust interval observer for switched systems with unknown inputs : application to vehicle dynamics estimation[END_REF]. For the simulation scenario, three subsystems are defined for v 1 x = 8.5m/s, v 2 x = 13.55m/s and v 3 x = 18.05m/s. In addition, the following switching law is considered:

σ(t) =    1 if vx ∈ [V 0 x , V 1 x [ 2 if vx ∈ [V 1 x , V 2 x [ 3 if vx ∈ [V 2 x , V 3 x ] (38)
with

v k x = V k x -V k-1 x 2
for k = 1, 2, 3. Furthermore, it is assumed that the cornering stiffness parameters are affected by ±10% uncertainty in their nominal values. The proposed algorithm is evaluated by introducing an actuator bias into the data provided by the encoder. Two fault zones are considered and represented with the switching signal in Figure 4. The first fault zone is in a stationary space where the second mode remains unchanged. The second zone is in a space where the vehicle dynamics switches quickly from the second to third mode and from the third to second mode, alternatively. According to the procedure described in Section III, the corresponding gain matrices of the UIO (9) for the lateral dynamics system (8) are designed. Due to space limitations, the gain matrices are omitted. The resulting attenuation level is γ = 0.0145 for an actuator bias of 0.25rad, α = 0.5 and a sampling time T s = 0.6s. Simulation results of the actuator fault detection are depicted in Figure 5. The gains which are obtained by solving the optimization problem in Theorem 2 allow to attenuate the effect of the uncertainties and the residual is almost zero. On the other hand, regarding the fault zone, where the vehicle dynamics switches alternatively between two modes, and the zone where no switch occurs, it can be noticed that the residual vector r k is outside the estimated set R. The fault is well detected in both cases until it disappears. Furthermore, there is no delay in the fault indication. The fault detection accuracy is equal in both zones, and, the switching does not influence the detection quality.

V. CONCLUSION

In this paper, a new vehicle robust fault detection scheme, which combines switched UIO and zonotopes has been proposed. A multiple model switching structure is considered to take into account the longitudinal velocity variation. Then, Multiple SUIO-based FD are constructed to decouple the effect of UI from the residual vector. The gain matrices have been designed using a common ISS-Lyapunov function and an LMI formulation is obtained. As presented on simulation, the fault detection results are encouraging. The use of input-to-state stability concept guarantees uncertainties effect attenuation. Moreover, it has been shown that fast vehicle longitudinal velocity switching does not influence the fault detection quality. Future work will be devoted to extend the fault detection method to fault isolation. Furthermore, another trivial extension of the method is to consider the variation of the longitudinal velocity and therefore adopt a Takagi-Sugeno representation of the vehicle lateral model. Reduction of the conservatism of the proposed LMIs using a multiple quadratic ISS-Lyapunov function will also be considered.

Fig. 1 .

 1 Fig. 1. Main steps for residual interval computation

Fig. 2 .

 2 Fig. 2. Lateral wind gust input Fw.

Fig. 3 .

 3 Fig. 3. Steering angle δ f , yaw rate ψ and longitudinal velocity vx.

Fig. 4 .

 4 Fig. 4. Switching law and faults zone.

Fig. 5 .

 5 Fig. 5. Residue

Sara Ifqir, Naima Ait-Oufroukh, Dalil Ichalal and Saïd Mammar are with IBISC Laboratory, Univ Evry, Paris-Saclay University, Evry,France.

The zero order hold method is used to obtain[START_REF] Ifqir | Robust interval observer for switched systems with unknown inputs : application to vehicle dynamics estimation[END_REF].

The word "uniformly" refers to uniformity with respect to switching signals.