Bernadette Bouchon-Meunier
email: bernadette.bouchon-meunier@lip6.fr

Christophe Marsala
email: christophe.marsala@lip6.fr

Entropy and monotonicity in artificial intelligence

Keywords: Entropy, Monotonicity, Measure of fuzziness, Intuitionistic entropy measure, Divergence. Artificial Intelligence

Introduction

The concept of information is complex and corresponds to factual elements regarding an event or a situation, such as observations or news, as well as the amount of knowledge brought by these factual elements. A piece of information can be analysed either from the point of view of its content, semantics or meaning or from the point of view of its form, encoding or label. Claude Shannon [START_REF] Shannon | The mathematical theory of communication[END_REF] and Norbert Wiener [START_REF] Wiener | Cybernetics, or control and communication in the animal and the machine[END_REF] simultaneously introduced the first measure to evaluate information in 1948 in a general theory of communication. Both of them were only considering the information communicated by the observation of the occurrence of a message or an event among a set of messages or events, respectively in communications systems and in the framework of cybernetics. The numerical measure they introduced was based on probabilities and had nothing to do with semantics. In reaction to this proposal, Rudolf Carnap and Yehoshua Bar-Hillel [START_REF] Carnap | An outline of a theory of semantic information[END_REF] proposed to formalise the concept of semantic information on the basis of a logical approach of natural language, through the so-called amount of information. Later on, Edwin Thompson Jaynes [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] established the maximum entropy principle to make inferences on the basis of partial information, by looking for the maximum entropy, given the available knowledge. Solomon Kullback then introduced with Richard Leibler and developed the concept of discrimination measure or divergence enabling to compare two probability distributions [START_REF] Kullback | On information and sufficiency[END_REF][START_REF] Kullback | Information theory and statistics[END_REF] .

Artificial intelligence was born in 1956 during the well-known "Dartmouth Summer Research Project on Artificial Intelligence" [START_REF] Mccarthy | A proposal for the Dartmouth summer research project on artificial intelligence[END_REF]. Even though it was not a direct emanation of cybernetics, it was strongly interrelated with it and it would have looked natural to consider Wiener's measure of information in Artificial Intelligence. It was not the case because Artificial Intelligence originally limited itself to symbolic methods, based on classical logic, rejecting any kind of numerical treatment of information, including probabilities. It might have looked relevant to use Carnap-Bar Hillel's measure of information, but this theory has never been applied to Artificial Intelligence, because of its numerical nature.

It is only in the 80s that some branches of Artificial Intelligence established a bridge with information theory. The Kullback discrimination measure [START_REF] Lemmer | Efficient minimum information updating for bayesian inferencing in expert systems[END_REF] and the maximum entropy principle [START_REF] Cheeseman | A method of computing generalized bayesian probability values for expert systems[END_REF] were for instance investigated in the framework of Bayesian inference and updating. One of the first consistent utilisations of measures of information in machine learning was Quinlan's decision trees [START_REF] Quinlan | Induction of decision trees[END_REF],

based on a concept of information gain used to choose an attribute providing its maximum value [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. At the same time, several other seminal works changed the view of numerical methods in Artificial Intelligence. Nilsson's probabilistic logic [START_REF] Nilsson | Probabilistic logic[END_REF] promoted Shannon's measure of information to develop a method to implement probabilistic entailment. Judea Pearl's early works and his book [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF] largely promoted the use of probabilities in Artificial Intelligence, his belief networks being considered as semantics-based systems. The optimisation of an entropy or a gain of information they recommended became popular in logic (for instance [START_REF] Goldszmidt | A maximum entropy approach to nonmonotonic reasoning[END_REF]) before being more largely applied in Artificial Intelligence.

Since this period, measuring the information provided by the observation of events has been a challenge. A number of quantities has been pointed out and studied in the literature to achieve this goal, called entropies or measures of information in the original probabilistic framework. Their extension to other frameworks such as fuzzy knowledge representation or its generalisations have given rise to the study of other aspects of information, for instance fuzziness or specificity. They have often been constructed by analogy with the probabilistic case, which may look artificial, but their properties go far beyond a simple analogy. They refer to the information contained in belief functions, for instance in [START_REF] Jiroušek | A new definition of entropy of belief functions in the Dempster-Shafer theory[END_REF][START_REF] Jiroušek | On properties of a new decomposable entropy of Dempster-Shafer belief functions[END_REF], or fuzzy systems, for instance in [START_REF] Zhang | Information structures and uncertainty measures in a fully fuzzy information system[END_REF]. Such entropies are used in Bayesian networks [START_REF] Gao | Learning bayesian network parameters via minimax algorithm[END_REF] or decision making systems [START_REF] Jiang | Strategy selection under entropy measures in movement-based three-way decision[END_REF][START_REF] Allahverdyan | Adaptive decision making via entropy minimization[END_REF]. We propose to continue the analysis initiated in [START_REF] Bouchon-Meunier | Entropy measures and views of information[END_REF] and formalised in [START_REF] Bouchon-Meunier | Entropy and monotonicity, in: Information Processing and Management of Uncertainty in Knowledge-Based Systems[END_REF], in focusing on the property of monotonicity, which appears essential for the notion of entropy and its utilisation in Artificial Intelligence. Most of the measures existing in the literature point out properties which have in common to be related to different aspects of monotonicity based on an order taking into account the context, the point of view or the knowledge representation. The entropy increases when the means to perform observations is refined, the concept of refinement taking various forms, from the utilisation of a better tool providing less vagueness to the utilisation of additional features to obtain more detailed information or, to express it in a more colourful language in the case where a photo brings information on an object, from zooming on the object to taking additional photos from another side of the object. The importance of the monotonicity of entropy, considered under the angle of the coarseness of information, has been pointed out in [START_REF] Amari | Information geometry in optimization, machine learning and statistical inference[END_REF]. Other views of the monotonicity of entropy have been proposed in random variable analysis, see for instance [START_REF] Artstein | Solution of shannon's problem on the monotonicity of entropy[END_REF].

This paper does not pretend to review all entropies or quantities of information introduced in the literature and related to artificial intelligence, in any way, as this would be far beyond the size of this manuscript. We only focus on a study of the common property of monotonicity of such measures with regard to a refinement of information, showing that the main differences between these quantities come from the diversity of orders defining such a refinement. Our aim is to propose a clarification of the concept of refinement of information and the underlying monotonicity, and to illustrate this paradigm by the utilisation of such measures in Artificial Intelligence.

The paper is organised as follows. In Section 2, we introduce entropy measures and definitions of monotonicity according to three different forms associated with different visions of the refinement of information. In Section 3, we illustrate these visions on various classic entropy measures: probabilistic entropies, measures of fuzziness, similarity relation-based entropy measures, entropies in the settings of the Atanassov intuitionistic fuzzy sets, and we develop also a study on popular divergence measures. In Section 4, we review some common uses of entropies in Artificial Intelligence and their highlighted properties of monotonicity. Finally, some conclusions are drawn and a set of perspectives is presented in Section 5.

Monotonicity of entropy measures

Shannon [START_REF] Shannon | The mathematical theory of communication[END_REF] clearly based his definition of entropy, considered as "measures of information, choice and uncertainty" on a concept of monotonicity, as he states "The uncertainty of y is never increased by knowledge of x. It will be decreased unless x and y are independent events, in which case it is not changed", which can be regarded as an interpretation of recursivity he furthermore requires from entropy. Additivity and recursivity are among the most important algebraic properties of entropy [START_REF] Aczél | On Measures of Information and their Characterizations[END_REF], and they imply an increase of the entropy resulting from the refinement of information acquired on an event through observations. Later on, Renyi [START_REF] Rényi | On measures of entropy and information[END_REF] introduced the first of a long list of generalisations of Shannon entropy, still satisfying a property of additivity. It is worth mentioning Mugur-Schächter's work [START_REF] Mugur-Schächter | The general relativity of descriptions[END_REF] on the general relativity of descriptions. She considers that any process of knowledge extraction is associated with epistemic operators called a delimiter and a view, representing the influence of the context and the observation tool on the considered event.

A refinement of information results from a change in the observation tool. In his generalised information theory, Kampé de Fériet [START_REF] Kampé De Fériet | Mesures de l'information par un ensemble d'observateurs[END_REF][START_REF] Kampé De Fériet | Mesure de l'information fournie par un événement[END_REF] takes into account observers and also requests a monotonicity of information with respect to an order on events.

Information theory do not pretend to evaluate all aspects, and it provides an evaluation of the decrease of uncertainty after an observation of events by means of entropies.

Given the amount of data available in the numerical world, which is covering all aspects of modern life, evaluating information is a major issue. All tools enabling the user to compare two pieces of information, to evaluate the information available in a given environment, to make diagnosis or prediction on the basis of information provided by observations or data, to aggregate chunks of information, are useful. Unfortunately, there are many such tools and it is difficult to see their common features. This is why we propose to analyse measures of information and to revisit classic approaches of information evaluation in order to focus on monotonicity which we consider the most natural and relevant property requested from such a tool. Let us consider a set of objects or events that represent the real world. In this paper, for the sake of simplicity, we only consider finite sets, but this work could be generalised to non-countable sets. We use the notation proposed in the seminal paper by Aczél and Daróczy on the so-called inset entropy [START_REF] Aczél | A mixed theory of information. I: Symmetric, recursive and measurable entropies of randomized systems of events[END_REF] to formalise the available information on the set of objects or events and taking into account the context, the point of view and the chosen knowledge representation.

We consider an algebra B defined on a finite universe U. For any integer n > 0, we note:

• X n = {(x 1 , . . . , x n ) | x i ∈ B, x i ∩ x j = 0 if i = j, ∀i, j = 1, . . . , n}; • P n = {(p 1 , . . . , p n ) | p i ∈ [0, 1]
}, p i being associated with x i through a function p : B → [0, 1], a particular case being a probability measure defined on (U, B);

• W n = {(w x1 , . . . , w xn ) | w xi ∈ R + , ∀i = 1, . . . , n}, a family of n-tuples of weights 1 associated with n-tuples of B through a function f : B → R +n , such that f (x 1 , . . . , x n ) = (w x1 , . . . , w xn ).
Similarly to the definition of inset entropy [START_REF] Aczél | A mixed theory of information. I: Symmetric, recursive and measurable entropies of randomized systems of events[END_REF], we introduce an entropy measure as a sequence of mappings E n : X n ×P n ×W n → R + satisfying several properties among a long list, for instance available in [START_REF] Aczél | On Measures of Information and their Characterizations[END_REF] or in [START_REF] Klir | Uncertainty-Based Information[END_REF].

In the sequel, for the sake of simplicity, we use the notation

     x 1 , . . . , x n p 1 , . . . , p n w x1 , . . . , w xn     
rather than ((x 1 , . . . , x n ), (p 1 , . . . , p n ), (w x1 , . . . , w xn )) to represent an element of X n × P n × W n , according to the notation used in [START_REF] Aczél | A mixed theory of information. I: Symmetric, recursive and measurable entropies of randomized systems of events[END_REF].

We claim that the most significant properties to characterise an entropy measure are relative to monotonicity with respect to a refinement of information which can take various forms, depending on a chosen order. Such a monotonicity corresponds to the natural idea that the more details, precision, certainty we obtain from the observation of objects or events, or equivalently the more refined information we have on them, the bigger the amount of information we have on these objects or events.

To use a metaphor, we can consider that we are facing a picture of an object providing some amount of information on it. We can first improve the light on the object before taking another picture in order to decrease the fuzziness of the details, or take another picture with a higher resolution, which gives more information on the object according to an increase of the precision, both cases corresponding to a monotonicity described in 2.1. We can also select a part of the object and make several pictures of this part, in a form of weak recursivity described in 2.2. We can finally partition the object into different parts and, for each of them, make more pictures providing a bigger amount of information on the object, in a form of weak additivity presented in 2.3.

We present these three forms of monotonicity which can be adapted to the knowledge representation we choose, as highlighted in the next sections.

O-monotonicity

The first form of monotonicity, noted O-monotonicity, is defined according to a given (partial) order on the elements: the monotonicity highlights a link between the order of the elements and the order induced by the measure.

Definition

Let ≺ be a partial order on a reliable observation of the objects or events, O-monotonicity can be written as follows: if

     x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n      ≺      x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n      then E n      x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n      ≤ E n      x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n     

Particular cases

Examples of monotonicity can be based on the following partial orders:

     x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n      ≺      x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n     
if and only if one of the following conditions is satisfied: (O1), called sharpness in [START_REF] De Luca | A definition of a nonprobabilistic entropy in the setting of fuzzy sets theory[END_REF], or (O2)

(O1) ∀i = 1, . . . , n, if w i ≥ 1 2 then w i ≤ w i ; 170 (O2) ∀i = 1, . . . , n, w i ≥ w i .
These examples are based on an order related to W . Other orders, for instance related to P , could be also used. Such examples will be studied later.

R-monotonicity

The second form of monotonicity, noted R-monotonicity, is based on a 175 decrease of the coarseness of a partition of the universe of discourse.

Definition

R-monotonicity can correspond to a property of weak recursivity defined as follows:

E n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n w x1 , w x2 , . . . , w xn      ≥ E n-1      x 1 ∪ x 2 , x 3 , . . . , x n p 1 + p 2 , p 3 , . . . , p n w x1∪x2 w x3 , . . . , w xn     
A particular case of weak recursivity is what we call the ψ-recursivity, defined for a function ψ : X 2 × P 2 × W 2 → R + as:

E n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n w x1 , w x2 . . . , w xn      = E n-1      x 1 ∪ x 2 , x 3 . . . , x n p 1 + p 2 , p 3 . . . , p n w x1∪x2 , w x3 . . . , w xn      +ψ      x 1 , x 2 p 1 , p 2 w x1 , w x2      E 2      x 1 , x 2 p1 p1+p2 , p2 p1+p2 w x1 , w x2      .

Particular case

The classic property of recursivity corresponds to:

ψ 0      x 1 , x 2 p 1 , p 2 w x1 , w x2      = p 1 + p 2 ,
where the weights are not taken into account.

A-monotonicity

The third form of monotonicity, noted A-monotonicity, is based on the consideration of a secondary finite universe U and an algebra B on U providing more details on the observed phenomenon or object, through additional observations.

Definition

Similarly to the situation on U, we consider for any integer m

• X m = {(x 1 , . . . , x m ) | x i ∈ B , ∀i}; • P m = {(p 1 , . . . , p m ) | p i ∈ [0, 1]}, p i being associated with x i through a function p : B → [0, 1]; • W m = {(w 1 , . . . , w m ) | w i ∈ R + ;
∀i}, a family of m-tuples of weights associated with m-tuples of elements of B through a function f :

B → R + , such that f (x i ) = w i .
We further suppose that there exist two combination operators and • enabling us to equip the Cartesian product of U × U with similar distributions:

• P n P m = {(p 1 p 1 , . . . , p i p j , . . .) | p i p j ∈ [0, 1]}, p i p j being associated 195 with (x i , x j ) for any i and j through a function p p ,

• W n • W m = {(w 1,1 , . . . , w i,j , . . .) | w i,j ∈ R + , ∀i, j}, is defined through a function f • f : B × B → R + , such that: f • f (x 1 , x j ) = w i,j for all i = 1, . . . , n and j = 1, . . . , m.
Such a refinement leads to a property of weak additivity stating the following:

E n×m      (x 1 , x 1 ), (x 1 , x 2 ), . . . , (x i , x j ), . . . , (x n , x m ) p 1 p 1 , p 1 p 2 , . . . , p i p j , . . . , p n p m w x1,x 1 , w x1,x 2 , . . . , w xi,x j , . . . , w xn,x m      ≥ max      E n      x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n      , E m      x 1 , . . . , x m p 1 , . . . , p m w 1 , . . . , w m          

Particular case 200

The classic additivity property stands in the case where U and U are independent universes, p and p being probability distributions on (U, B) and (U , B ), weights generally not being taken into account. It yields:

E n×m      (x 1 , x 1 ), (x 1 , x 2 ), . . . , (x i , x j ), . . . , (x n , x m ) p 1 p 1 , p 1 p 2 , . . . , p i p j , . . . , p n p m w x1,x 1 , w x1,x 2 , . . . , w xi,x j , . . . , w xn,x m      = E n      x 1 , . . . , x n p 1 , . . . , p n w x1 , . . . , w xn      + E m      x 1 , . . . , x m p 1 , . . . , p m w x 1 , . . . , w x m     

Diverse entropy measures

In this section, the previous definitions of monotonicity are studied in the 205 case of classical and well-known entropy measures.

Shannon and weighted entropies

It is well-known that the classic Shannon entropy defined as:

E S n (p) = - n i=1 p i log p i ,
only taking into account X n and P n , is additive and recursive and then Rmonotonous and A-monotonous.

Its generalisation to the case where weights are associated with events to rep-210 resent a cost or an importance is a weighted entropy defined on X n × P n × W n as follows [START_REF] Guiaşu | Weighted entropy[END_REF]:

E w n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n w 1 , w 2 , . . . , w n      = - n i=1 w i p i log p i .
The weighted entropy is O-monotonous with regard to the partial order (O2).

It is also recursive, and then R-monotonous when considering

w x1∪x2 = (p 1 w 1 + p 2 w 2 ) (p 1 + p 2 ) .
In addition, the weighted entropy is A-monotonous as soon as we consider an aggregation function at least equal to the maximum:

E n×m      (x 1 , x 1 ), . . . . . . , (x n , x m ) p 1 p 1 , . . . . . . , p n p m w 1 • w 1 , w 1 • w 2 , . . . , w n • w m      ≥ max      E n      x 1 , . . . , x n p 1 , . . . , p n w 1 , . . . , w n      , E m      x 1 , . . . , x m p 1 , . . . , p m w 1 , . . . , w m          
whenever w i • w j ≥ max(w i , w j ) for all i = 1, . . . , n and j = 1, . . . , m

An entropy of the same form as the weighted entropy has been introduced

by Zadeh [33] in the case where weights are replaced by membership degrees, and the so-called entropy of a fuzzy set is defined as:

E Z n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n µ 1 , µ 2 , . . . , µ n      = - n i=1 µ i p i log p i .
It is obviously still O-monotonous (with partial order (O2)), and A-monotonous according to [START_REF] Zadeh | Probability measures of fuzzy events[END_REF].

Parameterised entropies

A number of generalisations of Shannon's entropies flourished in the 60s and in the 70s, independent of W n and preserving some of the basic properties of Shannon's entropy. The first one was the Renyi's entropy of order α, defined as:

E Rα n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n w 1 , w 2 , . . . , w n      = 1 1 -α log n i=1 p i α ,
for a parameter α strictly positive and different from 1, the Shannon entropy corresponding to the limit case when α tends to 1. This quantity is additive, as proved in [START_REF] Rényi | On measures of entropy and information[END_REF], and then A-monotonous, but not recursive. It is easy to see that Renyi's entropy of order α is nevertheless R-monotonous.

Another parameterised form of entropy is Daróczy's entropy of type β, for a parameter β strictly positive and different from 1, defined as [START_REF] Daróczy | Generalized information functions[END_REF]:

E Dβ n (p) = 1 2 1-β -1 ( n i=1 p i β -1).
The Shannon entropy corresponds to the limit case when β tends to 1. When β = 2, we obtain a quantity proportional by a factor 2 to the Gini diversity index [START_REF] Breiman | Classification and Regression Trees[END_REF] used in the construction of decision trees by the Cart method.

It was known not to be either additive or recursive and was proved to satisfy recursivity of type β, equivalent to the weak recursivity we consider, with:

ψ      x 1 , x 2 p 1 , p 2 w x1 , w x2      = (p 1 + p 2 ) β .
It also satisfies an additivity of type β, implying the weak additivity. The Daróczy's entropy of type β is then R-monotonous and A-monotonous.

Measure of fuzziness

Shortly after the weighted entropy, another entropy measure was introduced by De Luca and Termini by analogy with the Shannon entropy, but in a nonprobabilistic framework [START_REF] De Luca | A definition of a nonprobabilistic entropy in the setting of fuzzy sets theory[END_REF], and then independently of P n . It is a measure of fuzziness, in the case where f is the membership function of a fuzzy set on U :

E DLT n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n w 1 , w 2 , . . . , w n      = - n i=1 w i log w i - n i=1 (1 -w i ) log(1 -w i ).
A major property of this quantity is its O-monotonicity with respect to the above mentioned partial order (O1) defining the sharpness.

230

It can further be observed that, in the case where the weights are possibility degrees, ie. max(w 1 , . . . , w n ) = 1, this measure of fuzziness is also weakly recursive and then R-monotonous:

E DLT n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n w 1 , w 2 , . . . , w n      ≥ E DLT n-1      x 1 ∪ x 2 , x 3 , . . . , x n p 1 + p 2 , p 3 , . . . , p n max(w 1 , w 2 ), w 3 , . . . , w n     

Entropy measures under similarity relations

We consider a similarity relation S on U = {x 1 , . . . , x n }, reflexive, symmetric and min-transitive. R.R. Yager [START_REF] Yager | Entropy measures under similarity relations[END_REF] defines the following entropy measure on

X n × P n × W n : E Sim n      x 1 ,
(O3) S S ⇔ S(x i , x j ) ≤ S (x i , x j ) ∀i, j.
This entropy measure is also A-monotonous, if we define a joint similarity relation S × S on the Cartesian product U × U as follows, for two similarity relations S defined on U and S defined on U :

S × S ((x i , y j ), (x k , y l )) = min S(x i , x k ), S (y j , y l )

for any x i and x k in U , any y j and y l in U .

Ambiguity or nonspecificity

In the framework of possibility distributions, corresponding to a fuzzy setbased knowledge representation in which a membership degree is interpreted as the possibility of the observed variable to take a given value, with a maximum equal to 1, a measure of ambiguity or non-specificity has been introduced by [START_REF] Higashi | Measures of uncertainty and information based on possibility distributions[END_REF] and called U-uncertainty. It is defined under the hypothesis that the x i are ranked according to a possibility distribution:

π 1 ≥ π 2 ≥ ... ≥ π n and π 1 > 0.
This quantity, independent of P , got more popularity when pointed out by [START_REF] Yuan | Induction of fuzzy decision trees[END_REF] in the induction of fuzzy decision trees.

E HK n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n π 1 , π 2 , . . . , π n      = n-1 i=1 (π i -π i+1 ) log i.
It was proved [START_REF] Klir | Fuzzy Sets, Uncertainty, and Information[END_REF] that this quantity is O2-monotonous. It is also additive, then A-monotonous, when the two possibility distributions are non-interactive, as follows:

E HK n×m      (x 1 , x 1 ), . . . , (x n , x m ) p 1 p 1 , . . . , p n p m min(π 1 , π 1 ), . . . , min(π n , π m )      ≥ max      E HK n      x 1 , . . . , x n p 1 , . . . , p n π 1 , . . . , π n      , E HK m      x 1 , . . . , x m p 1 , . . . , p m π 1 , . . . , π m          
It is also recursive, and then R-monotonous.

Divergence

Kulback and Leibler's divergence [START_REF] Kullback | On information and sufficiency[END_REF][START_REF] Kullback | Information theory and statistics[END_REF] is another form of entropy, when W n is identical with P n , corresponding to the following:

J KL n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n q 1 , q 2 , . . . , q n      = n i=1 p i log( p i q i ),
where p 1 , p 2 , . . . , p n and q 1 , q 2 , . . . , q n are two probability distributions on the same set of events, for instance a prior and a posterior distribution. They prove that J KL n is additive, and therefore A-monotonous. This divergence is also recursive, such that:

J KL n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n q 1 , q 2 . . . , q n      = J KL n-1         x 1 ∪ x 2 , x 3 . . . , x n p 1 + p 2 , p 3 . . . , p n q 1 + q 2 , q 3 . . . , q n         + (p 1 + p 2 )J KL 2      x 1 , x 2 p1 p1+p2 , p2 p1+p2 q1 q1+q2 , q2 q1+q2      It is therefore R-monotonous.

255

Parameterised forms of divergence were introduced, in the same spirit as the parameterised entropies. The first one is divergence of order α introduced by Renyi [START_REF] Rényi | On measures of entropy and information[END_REF] as the gain of information resulting from the replacement of q 1 , q 2 , . . . , q n by p 1 , p 2 , . . . , p n , defined as follows:

J Rα n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n q 1 , q 2 , . . . , q n      = 1 α -1 log n i=1 p α i q α-1 i ,
for α ≥ 0 and α = 1. It is additive, and therefore A-monotonous.

The second parameterised divergence was introduced by Rathie and Kannappan [START_REF] Rathie | A directed-divergence function of type β[END_REF] as the directed divergence of type β:

J RKβ n      x 1 , x 2 , . . . , x n p 1 , p 2 , . . . , p n q 1 , q 2 , . . . , q n      = 1 2 β-1 -1 ( n i=1 p β i q β-1 i -1).
The authors proved that it is R-monotonous. Furthermore, they prove that it has a form of strong non-additivity entailing the following:

J RKβ n×m      (x 1 , x 1
), (x 1 , x 2 ), . . . , (x i , x j ), . . . , (x n , x m ) p 1 p 1 , p 1 p 2 , . . . , p i p j , . . . , p n p m q 1 q 1 , q 1 q 2 , . . . , q i q j , . . . ,

q n q m      = J RKβ n      x 1 , . . . , x n p 1 , . . . , p n q 1 , . . . , q n      + n i=1 p β i q 1-β i J RKβ m      x 1 , . . . , x m p 1 , . . . , p m q 1 , . . . , q m     
which implies a property of A-monotonicity

Intuitionistic entropy measures

In this section, we consider the setting of the Atanassov intuitionistic fuzzy sets (AIFS) where several entropy measures have been introduced [START_REF] Guo | On the entropy for Atanassov's intuitionistic fuzzy sets: an interpretation from the perspective of amount of knowledge[END_REF][START_REF] Szmidt | New measures of entropy for intuitionistic fuzzy sets[END_REF]. First of all, some basics of AIFS are recalled.

Let X be a universe, an Atanassov intuitionistic fuzzy set (AIFS) A of X is defined [START_REF] Atanassov | Intuitionistic fuzzy sets[END_REF] by:

A = {(x, µ A (x), ν A (x))|x ∈ X} with µ : X → [0, 1], ν : X → [0, 1] and 0 ≤ µ A (x) + ν A (x) ≤ 1, ∀x ∈ X.
Here, µ A (x) and ν A (x) represent respectively the membership degree and the non-membership degree of x in A. Given an intuitionistic fuzzy set A of X, the hesitancy lying on the membership of x to A is the intuitionistic index of x to

A defined for all x ∈ X as π A (x) = 1 -(µ A (x) + ν A (x)). It is easy to see that we always have π A (x) ∈ [0, 1].
The inclusion of AIFS is defined as:

A ⊆ B if and only if µ A (x) ≤ µ B (x) and ν A (x) ≥ ν B (x), ∀x ∈ X.
The union of two AIFS A and B is defined as the AIFS A ∪ B such that

µ A∪B (x) = max(µ A (x), µ B (x)) and ν A∪B (x) = min(ν A (x), ν B (x)), ∀x ∈ X.
The intersection of two AIFS A and B is defined as the AIFS A ∩ B such that

µ A∩B (x) = min(µ A (x), µ B (x)) and ν A∩B (x) = max(ν A (x), ν B (x)), ∀x ∈ X. It can be easily seen that µ A∪B (x) ∈ [0, 1] and µ A∩B (x) ∈ [0, 1].

Definitions of entropy measures for AIFS.

Several works in AIFS theory have proposed the definition for an entropy of an intuitionistic fuzzy set A. With our notations, we represent these quantities as:

280 E IF S n      x 1 , . . . , x n p 1 , . . . , p n (µ 1 , ν 1 ), . . . , (µ n , ν n )      , with µ A (x i ) = µ i and ν A (x i ) = ν i .
Here, the weights from W n are thus defined

as tuples from [0, 1] × [0, 1].
There exist various definitions of entropy measures in the AIFS setting [START_REF] Guo | On the entropy for Atanassov's intuitionistic fuzzy sets: an interpretation from the perspective of amount of knowledge[END_REF],

independent of P n . For instance, the entropy measure given in [START_REF] Szmidt | New measures of entropy for intuitionistic fuzzy sets[END_REF] is defined as:

E IF S,S n      x 1 , . . . , x n p 1 , . . . , p n (µ 1 , ν 1 ), . . . , (µ n , ν n )      = 1 - 1 2n n i=1 |µ i -ν i |
In [START_REF] Guo | On the entropy for Atanassov's intuitionistic fuzzy sets: an interpretation from the perspective of amount of knowledge[END_REF], the following entropy measure is also introduced:

E IF S,G n      x 1 , . . . , x n p 1 , . . . , p n (µ 1 , ν 1 ), . . . , (µ n , ν n )      = 1 2n n i=1 1 -|µ i -ν i | (1 + π i ), with π i = 1 -(µ i + ν i ).
Another way to define an entropy measure is presented in [START_REF] Burillo | Entropy on intuitionistic fuzzy sets and on intervalvalued fuzzy sets[END_REF] where the definition is based on extensions of the Hamming distance and the Euclidean distance to AIFS. For instance, the following entropy measure is proposed:

E IF S,B n      x 1 , . . . , x n p 1 , . . . , p n (µ 1 , ν 1 ), . . . , (µ n , ν n )      = n i=1 π i 3.7.2.
Entropy measures for AIFS and monotonicity.

In [START_REF] Guo | On the entropy for Atanassov's intuitionistic fuzzy sets: an interpretation from the perspective of amount of knowledge[END_REF], it is recalled that, in the AIFS setting, a monotonicity property for 285 an entropy measure could be ensured by definition. The authors present several definitions that lie on the definition of a partial order on W n and the concept of less fuzzy than. The following definitions of partial order could be used.

(O4)

A is less fuzzy than B if for all x ∈ X µ A (x) ≤ µ B (x) and

ν A (x) ≥ ν B (x) if µ B (x) ≤ ν B (x), or µ A (x) ≥ µ B (x) and ν A (x) ≤ ν B (x) if µ B (x) ≥ ν B (x)

, and (O5)

A is less fuzzy than B if

µ A (x) ≤ µ B (x) and ν A (x) ≤ ν B (x), ∀x ∈ X,
O-monotonicity. It is easy to see that (O4) and (O5) yield two versions of O-monotonicity. These definitions of monotonicity produce particular forms of E IF S :

• E IF S,S n
satisfies the monotonicity based on (O4) as it is stated in [START_REF] Szmidt | New measures of entropy for intuitionistic fuzzy sets[END_REF];

• E IF S,G n
satisfies the monotonicity based on (O4), as it is stated in [START_REF] Guo | On the entropy for Atanassov's intuitionistic fuzzy sets: an interpretation from the perspective of amount of knowledge[END_REF];

• E IF S,B
n satisfies the monotonicity based on (O5) as it could be found in [START_REF] Burillo | Entropy on intuitionistic fuzzy sets and on intervalvalued fuzzy sets[END_REF] where this entropy is given as example.

R-monotonicity.

The measure E IF S,B n satisfies the R-monotonicity if we have

E IF S,B n      x 1 , . . . , x n p 1 , . . . , p n (µ 1 , ν 1 ), . . . , (µ n , ν n )      ≥ E IF S,B n-1      x 1 ∪ x 2 , x 3 , . . . , x n p 1 + p 2 , p 3 , . . . , p n (max(µ 1 , µ 2 ), min(ν 1 , ν 2 )), (µ 3 , ν 3 ), . . . , (µ n , ν n )     
.

considering the union of AIFS as defined in the introduction of this section.

Hereafter, for the sake of simplicity, we note these two measures E IF S,B n and

E IF S,B n-1
respectively.

We have

E IF S,B n -E IF S,B n-1 = 1 -µ 1 -ν 1 + 1 -µ 2 -ν 2 -1 + max(µ 1 , µ 2 ) + min(ν 1 , ν 2 ) = 1 + (max(µ 1 , µ 2 ) -µ 1 -µ 2 ) + (min(ν 1 , ν 2 ) -ν 1 -ν 2 )
and thus

E IF S,B n -E IF SB n-1 = 1 -min(µ 1 , µ 2 ) -max(ν 1 , ν 2 )
This corresponds to the intuitionistic index of the intersection of AIFS, and thus, as a consequence, we have

E IF S,B n -E IF S,B n-1
≥ 0 and E IF S,B n satisfies the R-monotonicity.

Entropies in Artificial Intelligence

Entropies are very commonly used in Artificial Intelligence. Indeed, their monotonicity properties could be one of the main reasons for this success. In the following, some applications of entropies in Artificial Intelligence are presented and the kind of monotonicity involved is highlighted (that could explain why such an entropy is chosen in such applications). Our aim is not to propose a complete review of such applications but to show that an entropy can be considered as a universal tool.

Beyond the search of the maximum entropy common when probabilistic measures are used, monotonicity or maximisation of entropy has also been the core of methods based on non-probabilistic entropies.

Maximum Entropy Principle

A derived utilisation of the concept of monotonicity is the very commonly used maximum entropy principle (Maxent). It was first proposed by Jaynes [START_REF] Jaynes | Information theory and statistical mechanics[END_REF],

in the simple case where weights are not involved, as a way to choose the most appropriate probability distribution to cope with the uncertainty, as the one being "maximally noncommittal with regard to missing information". Jaynes introduced this principle in the case of the Shannon entropy, because of the easiness to solve the optimisation problem of maximising the entropy under

Entropies in Machine Learning

As previously said in the introduction, one of the well-known uses of entropies is in machine learning where the Shannon entropy is very popular.

In supervised learning, a common use of such a measure is dedicated to the learning of a decision tree from a training data set [START_REF] Quinlan | Induction of decision trees[END_REF], [START_REF] Breiman | Classification and Regression Trees[END_REF]. Shannon entropy is not the unique measure to be used in such a process [START_REF] Marsala | Ranking attributes to build fuzzy decision trees: a comparative study of measures[END_REF] but it is one of the most efficient. We can for instance also cite the Gini index of diversity which is also popular [START_REF] Breiman | Classification and Regression Trees[END_REF].

Here, for this kind of tree-like splitting processes, the R-monotonicity is the main property that is sought for. Indeed, measures following this kind of monotonicity property enable the better choice of description attributes when splitting the training set to reduce the uncertainty related to the prediction of the class.

In the building of fuzzy decision trees, other measures have been introduced.

For instance, Renyi's entropy was used in presence of unbalanced datasets in [START_REF] Gajowniczek | Comparison of decision trees with Rényi and Tsallis entropy applied for imbalanced churn dataset[END_REF]. In a fuzzy setting, De Luca and Termini's non-probabilistic entropy [START_REF] De Luca | A definition of a nonprobabilistic entropy in the setting of fuzzy sets theory[END_REF] was used to construct fuzzy decision trees [START_REF] Cios | Continuous ID3 algorithm with fuzzy entropy measures[END_REF] and also for feature selection in classification [START_REF] Luukka | Feature selection using fuzzy entropy measures with similarity classifier[END_REF]. In both approaches, this entropy is used to select attributes bringing the maximum information.

Another example is the use by Yuan's and Shaw's of the measure of ambiguity issued from Higashi and Klir's measure of ambiguity [START_REF] Klir | Fuzzy Sets, Uncertainty, and Information[END_REF] to build fuzzy decision trees [START_REF] Yuan | Induction of fuzzy decision trees[END_REF].

In these cases, with such measures able to handle fuzzy sets, the R-monotonicity is not the main property needed, but the O-monotonicity is more important when used to compare membership functions.

Shannon entropy is also very frequently used as a regularisation term. It is for instance the case in semi-supervised learning [START_REF] Grandvalet | Semi-supervised learning by entropy minimization[END_REF]. In this setting, considering a set of variables involved in the optimisation of a given function, entropy regularisation introduces an entropy term in the function to optimise (either maximise or minimise) in order to lead to a sparse distribution on the value (minimisation of the entropy) or a homogeneous distribution (maximisation of the entropy). Here again, it is the R-monotonicity of the entropy that is called in: the involved process is similar to the one described for the maximum entropy principle.

In unsupervised learning, a popular example of the use of the Shannon entropy lies in a regularisation term for the cost function in fuzzy clustering [START_REF] Kruse | Fundamentals of fuzzy clustering[END_REF].

This regularisation term should be built on the basis of probabilities, weights, fuzzy memberships,... In this case, there are two possibilities, either the entropic regularisation term must be maximised and leads to a uniform distribution of the related values, or it should be minimised to leads to a distribution with null values in order to introduce sparsity in the trained model.

Entropies in Other Applications Domains

One notable use of entropies could be found in biology, in the study of ecological systems. The Shannon entropy, as well as the Gini-Simpson index, could be used to evaluate the diversity of the species in an ecosystem [START_REF] Guiaşu | Entropy in Ecology and Ethology[END_REF], [START_REF] Guiaşu | The weighted Gini-Simpson index: Revitalizing an old index of biodiversity[END_REF]. In this framework, each specie is associated with a probability of occurrence that enables the definition of a distribution of probabilities regarding all the species present in the ecosystem. Here, the Shannon entropy applied to this distribution is a suitable tool to evaluate the diversity of the species and enables the comparison among ecosystems or to model species geographic distributions [START_REF] Austin | Spatial prediction of species distribution: an interface between ecological theory and statistical modelling[END_REF], [START_REF] Phillips | Maximum entropy modeling of species geographic distributions[END_REF].

Similarly to the case of the Maxent principle, it can be highlighted that the R-monotonicity is the property needed in this process.

There are several domains in which the Maxent principle is used. We can cite for instance, non-monotonic reasoning where Shannon entropy is used to build a probability distribution during the decision process in order to select plausible conclusions [START_REF] Goldszmidt | A maximum entropy approach to nonmonotonic reasoning[END_REF]. Another use can be found in description logics [START_REF] Wilhelm | Counting strategies for the probabilistic description logic ALC me under the principle of maximum entropy[END_REF] where this principle is also used to make a choice between models.

Conclusion

Entropy and measures of information have been extensively studied for 70 years. The original quantities dealing with probabilities of events have been extended to take into account fuzzy sets, intuitionistic fuzzy sets and other representation models of uncertainty and imprecision. Most of the proposed measures are only based on a formal analogy between the introduced quantities and classic entropies, in spite of the fact that their purpose is different, entropies measuring the decrease of uncertainty resulting from the occurrence of an event, while fuzzy set related measures evaluate the imprecision of events and the fuzziness or non-specificity of the studied observations.

All these quantities have in common a few or many fundamental properties, depending on the case. Various works have listed such properties, for instance [START_REF] Aczél | On Measures of Information and their Characterizations[END_REF], [START_REF] Klir | Uncertainty-Based Information[END_REF], [START_REF] Klir | Fuzzy Sets, Uncertainty, and Information[END_REF] and shown which quantities satisfy or do not satisfy them. Attempts have also been done to exhibit classes of quantities with a similar behaviour with regard to sets of properties [START_REF] Bouchon | Entropic models[END_REF].

In this paper, we highlight the common property of monotonicity of entropy measures with regard to a refinement of information, showing that the main differences between these quantities come from the diversity of orders defining such a refinement. This paper is not intended to provide a review of all entropy measures existing in the literature, but to clarify the concept of refinement of information and the underlying monotonicity, and to illustrate this paradigm by classic examples in a sample of knowledge representation environments, namely the classic probabilistic one, the fuzzy one, the similarity-based one and the intuitionistic fuzzy framework. A focus is put on the importance of monotonicity when entropies or measures of information are used in Artificial Intelligence.

In the future, we will point out new forms of monotonicity useful in Artificial

Intelligence and we will provide some hints to choose one or the other measure of information in a given context.

  x 2 , . . . , x n p 1 , p 2 , . . . , p

	with S i = p The similarity reflects a point of view on the n events, making explicit to xj ∈U
	which extent they are similar with regard to a given criterion. If we consider
	two different points of view, symbolised by two similarity relations S and S ,
	we can show that this entropy measure is O-monotonous with respect to the
	order (O3):
		x 1 , . . . , x n				x 1 , . . . , x n	
	   	p 1 , . . . , p n S 1 , . . . , S n	   	≤	   	p 1 , . . . , p n S 1 , . . . , S n	   
	if and only if similarities S and S satisfy:
						
						   	= -	xi∈U	p i log S i
				13

n S 1 , S 2 , . . . , S n j S(x i , x j ) for all i = 1, . . . n.

In the following, for the sake of simplicity, wx i is denoted w i when the meaning of i is clear.

the condition of probabilities having a sum equal to 1. His aim was to benefit of a form of monotonicity of the entropy with respect to missing information tolerance. In his work, Jaynes argued [START_REF] Jaynes | Information theory and statistical mechanics[END_REF] that "Mathematically, the maximumentropy distribution has the important property that no possibility is ignored" and he based his proof on the fact that "if all the p i are equal, the quantity

) is a monotonic increasing function of n". In our setting, this Maxent principle can be related to the R-monotonicity of the Shannon entropy: the higher the number of non-null probabilities, the higher the entropy. It is easy to see this fact when considering a set of events on which a distribution of probability should be identified. Let P 1 and P 2 be two distributions of probability on this set of events. Between P 1 and P 2 , Jaynes' principle argues that the distribution that the most covers the set of events should be preferred, which means the one that maximises the number of non-null probabilities associated with the events. In the case where P 1 and P 2 have the same number of non-null probabilities, the one that provides a more homogeneous distribution is preferred (the one that maximises the Shannon entropy). This property is, in fact, a side effect of the Jaynes' maximum entropy principle.

Since the R-monotonicity could be considered as the heart of this principle, it highlights the fact that choosing between two distributions by using Maxent, does not provide any specific information on the relative position of two distributions of probability but only a general information about their relative spread.

In the same spirit, and again in the case where the only available information is provided by probability distributions, Kullback [START_REF] Kullback | Information theory and statistics[END_REF] introduced the concept of minimum discrimination information, corresponding to the minimum value of the Kullback and Leibler's divergence. These two principles have been extended to other entropies or divergences and widely used in Artificial Intelligence.


Ramzi Othman

Amine Ammar

Khalid H Almitani

Reduced modeling Computation of layered soil's harmonic Green functions

Keywords: Ground-borne vibration, Layered viscoelastic soil, Proper Generalized Decomposition (PGD), Model reduction technique, Parametric analysis, Finite-element, Boundary-element, Railway traffic, Earthquakes

Ground-borne vibrations are disturbing to human beings. In order to model and reduce these vibrations, the calculation of the harmonic Green's-functions of the soil is highly required since the most effective numerical solution to predict ground-borne vibration is to couple the finiteelement and the boundary-element methods. In this work, we elaborate a direct space-frequency formulation that is especially suitable of material exhibiting visco-elastic behavior. Moreover, an original formulation of the Proper Generalized Decomposition (PGD) reduction technique was undertaken in the frequency-space domain. This new approach highly shortens time computing costs. It is mainly valuable when a parametric study is concerned. Moreover, the PGD resolution allows providing an off-line solution, thus, the solution is calculated only once.

Introduction

Railway traffic causes vibrations which are mainly due to the uneven nature of railway tracks, misalignment of wheels of the trains or discontinuities of the wheels and tracks [START_REF] Degrande | A numerical model for ground-borne vibrations from underground railway traffic based on a periodic finite elementboundary element formulation[END_REF]. For similar reasons, vibrations are also generated by other terrestrial vehicles traffics. Earthquakes can also initiate vibrations. These vibrations propagate through the soil to the foundations of neighboring buildings and structures, causing disturbances, noise pollution and sometimes damage [START_REF] Chebli | 3D periodic BE-FE model for various transportation structures interacting with soil[END_REF][START_REF] Brookes | The dynamic interaction effects of railway tunnels: Crossrail and the Grand Central Recording Studios[END_REF]. Moreover, they can also have unwelcome effects on the health of occupants and they can disturb the work of vibration-sensitive equipment. To prevent these undesirable effects, ground-borne vibrations have to be reduced or isolated as much as possible. The formulation of a suitable model is difficult because of the number of parameters involved in describing the underground environment. For modelling of vibrations due to the underground train traffic, the tunnel structure and the surrounding soil are assumed much stiffer than the track [START_REF] Fryba | Vibrations of solids and structures under moving loads[END_REF]. Due to this assumption, the railway induced vibrations problem is approached in two steps. Firstly, the wheel-rail-(track) interaction problem is separately solved to derive the dynamic axle forces. These forces are considered as an input external force for the second problem that deals with the track-(tunnel)-soil interaction problem. These problems aim at predicting the vibration propagation in the underground environment. Various models are designed to predict railway traffic-induced vibration [START_REF] Germonpré | Contributions of longitudinal track unevenness and track stiffness variation to railway induced vibration[END_REF]. The numerical solutions obtained earlier were time-consuming because of the slow processing speeds of computers. Analytical approaches were preferred. However, with the help of increasing processing speeds of computers, the calculation time can be reduced admirably. Now, numerical models are highly preferred in the modelling of underground railway-induced vibrations, mainly by using finite/infinite-element models [START_REF] Yang | Soil vibrations caused by underground moving trains[END_REF] and coupled finite element-boundary element (FE-BE) models [START_REF] Gupta | A comparison of two numerical models for the prediction of vibrations from underground railway traffic[END_REF][START_REF] Sheng | Modelling ground vibration from railways using wavenumber finite-and boundary-element methods[END_REF]. Even though two-dimensional numerical models offer computation times that are much shorter than three-dimensional models [START_REF] Andersen | Coupled boundary and finite element analysis of vibration from railway tunnels-a comparison of two-and three-dimensional models[END_REF], they are unable to account for wave propagation along the track. In addition, they cannot accurately simulate the radiation damping of the soil [START_REF] Gupta | A comparison of two numerical models for the prediction of vibrations from underground railway traffic[END_REF]. On the other hand, 3-D numerical models are highly expensive in terms of computation requirements. In addition to theoretical approaches [START_REF] Zhang | A semi-analytical model of the trainfloating slab track-tunnel-soil system considering the non-linear wheel/rail contact[END_REF][START_REF] He | Analytical method for calculation of ground vibration from a tunnel embedded in a multi-layered half-space[END_REF][START_REF] He | Analytical model for vibration prediction of two parallel tunnels in a full-space[END_REF][START_REF] He | A three-dimensional semi-analytical method for calculating vibrations from a moving load on a periodic jointed tunnel[END_REF][START_REF] Yuan | The influence of pore-fluid in the soil on ground vibrations from a tunnel embedded in a layered half-space[END_REF], numerical methods that exploit the invariance in structure (track and tunnel) longitudinal axis direction (2.5D finite/infinite element or FE-BE models) are attracting a particular attention from the scientific community [START_REF] Jin | A 2.5D finite element and boundary element model for the ground vibration from trains in tunnels and validation using measurement data[END_REF][START_REF] Yang | 2.5D vibration of railway-side buildings mitigated by open or infilled trenches considering rail irregularity[END_REF][START_REF] Feng | Preliminary analysis of the interaction between tunnel and soil by the 2.5d fem-bem method[END_REF][START_REF] Yaseri | Javady. 2.5D coupled FEM-SBFEM analysis of ground vibrations induced by train movement[END_REF][START_REF] Godinho | Numerical simulation of soil-structure elastodynamic interaction using iterative-adaptive BEM-FEM coupled strategies[END_REF][START_REF] He | A 2.5-D coupled FE-BE model for the dynamic interaction between saturated soil and longitudinally invariant structures[END_REF]. These techniques are effective when considering very long structures partially or totally embedded in soil [START_REF] Andersen | Coupled boundary and finite element analysis of vibration from railway tunnels-a comparison of two-and three-dimensional models[END_REF][START_REF] Aubry | Periodic BEM and FEM-BEM coupling: applications to seismic behaviour of very long structures[END_REF]. They indeed take into account the soil-structure coupling which is of great importance [START_REF] Chebli | Response of periodic structures due to moving loads[END_REF][START_REF] Chebli | Dynamic response of high-speed ballasted railway tracks: 3D periodic model and in situ measurements[END_REF]. The structure is modeled by techniques of structural dynamics that can be deterministic or probabilistic [START_REF] Arnst | A nonparametric probabilistic model for ground-borne vibrations in buildings[END_REF]. The soil is regarded as a layered half-space, which is made of horizontal elastic homogenous layers that rest on a homogenous elastic half-space. The Green's-function [START_REF] Strukelj | Evaluation of Green's-function for vertical pointload excitation applied to the surface of a layered semi-infinite elastic medium[END_REF] or the fundamental solutions [START_REF] Miklowitz | The Theory of Elastic Waves and Waveguides[END_REF] can be used to describe the dynamic characteristics of a layered half-space. The fundamental solutions are based on the computation of integrals that require the introduction of a fictitious boundary, across which the radiation conditions have to be satisfied [START_REF] Premrov | Solving exterior problems of wave propagation based on an iterative variation of local DtN operators[END_REF]. The Green's functions are approached in two different ways: approximate and analytical solutions. An example of approximate methods is the thin-layer method [START_REF] Kausel | An explicit solution for the Green's-function for dynamic loads in layered media[END_REF]. The analytical methods give solutions that are expressed in terms of integrals that have infinite or semi-infinite integration paths [START_REF] Vostroukhov | Steady-state response of a stratified half-space subjected to a horizontal arbitrary buried uniform load applied at a circular area[END_REF][START_REF] Pak | Three-dimensional Green's-functions for a multilayered half-space in displacement potentials[END_REF][START_REF] Wang | Displacements and stresses due to a vertical point load in an inhomogeneous transversely isotropic halfspace[END_REF][START_REF] Kobayashi | Evaluation of Green's-function on semi-infinite elastic medium[END_REF]. Franssens [START_REF] Franssens | Calculation of the elastodynamic Green's function in layered media by means of a modified propagator matrix method[END_REF] used the propagator matrix method to calculate the inverse Fourier integrands. Cairo et al. [START_REF] Cairo | Analysis of pile groups under vertical harmonic vibration[END_REF] and Kausel [START_REF] Kausel | Generalized stiffness matrix method for layered soils[END_REF] used the Stiffness Matrix Method. Zhang et al. [START_REF] Zhang | An Efficient Numerical Method for Computing Synthetic Seismograms for a Layered Half-space with Sources and Receivers at Close or Same Depths[END_REF] used wavenumber integration method. Shaukath et al. [START_REF] Shaukath | Computation of harmonic green'sfunctions of a homogeneous soil using an axisymmetric finite element method[END_REF] computed Green functions using an axisymmetric finite element approach and a time-explicit integration. In this work, we make a direct space-frequency formulation that is especially suitable of material exhibiting visco-elastic behavior. The advantage of this formulation is that it allows finding the real and complex part of the displacement field and thus permitting to find the frequency response without making any explicit or implicit time integration. Moreover, the Proper-Generalized-Decomposition (PGD) reduction method [START_REF] Ammar | On the space-time separated representation of integral linear viscoelastic models[END_REF][START_REF] Metoui | The proper generalized decomposition for the simulation of delamination using cohesive zone model[END_REF][START_REF] Dumon | Proper Generalized Decomposition method for incompressible Navier-Stokes equations with a spectral discretization[END_REF][START_REF] Boucinha | Space-time proper generalized decompositions for the resolution of transient elastodynamic models[END_REF][START_REF] Ammar | Non-incremental proper generalized decomposition solution of parametric uncoupled models defined in evolving domains[END_REF] is used to drastically reduce the computation time. It is a very powerful technique to deal with large-sized numerical models in a wide range of engineering fields [START_REF] Ibàñez | Multiscale proper generalized decomposition based on the partition of unity[END_REF][START_REF] Vitse | Dealing with a nonlinear material behavior and its variability through PGD models: Application to reinforced concrete structures[END_REF][START_REF] Louf | Fast validation of stochastic structural models using a PGD reduction scheme[END_REF]. However, the PGD method has rarely been used in the space-frequency domain. In the papers [START_REF] Barbarulo | On a PGD model order reduction technique for mid-frequency acoustic[END_REF] and [START_REF] Barbarulo | Proper Generalized Decomposition applied to linear acoustic: A new tool for broad band calculation[END_REF] the dynamic problem is solved with technique which combines the Variational Theory of Complex Rays (VTCR) with Proper Generalized Decomposition (PGD) and does not require the resolution of acoustic problems at many frequencies. The PGD is used to find a representation of the approximate solution which is separated between two variables, the wave propagation direction and the frequency. However this works uses the pressure as a primary variable and does not uses a tensoriel separation of the global operator (which can come back a bit expensive in terms of allocated memory). A formulation of the PGD method in frequency domain was also presented in [START_REF] Haris | Advanced parametric space-frequency separated representations in structural dynamics: A harmonic-modal hybrid approach[END_REF], [START_REF] Aguado | Real-time monitoring of thermal processes by reduced-order modeling[END_REF] and [START_REF] Germoso | Efficient PGD-based dynamic calculation of non-linear soil behavior[END_REF]. The PGD method is then used to compute space-frequency separated representations by considering the frequency as an extra coordinate [START_REF] Haris | Advanced parametric space-frequency separated representations in structural dynamics: A harmonic-modal hybrid approach[END_REF][START_REF] Germoso | Efficient PGD-based dynamic calculation of non-linear soil behavior[END_REF]. Our work is distinguished by: -The use of the viscoelastic constitutive behaviour to introduce naturally the damping. -Consequently the kinematic field (displacement) is naturally formulated with the complex components involving a system which is twice larger than the usally solved in harmonic resolution.

-The operators of our PGD decomposition (with parametrization or using frequency as an extra coordinate) are written in terms of tensor decomposition alleviating the cost of the reformulation of the PGD for each case of decomposition.

Problem statement

In this work, we aim at computing the Green functions of a layered soil, while scaling down the computation time using a model reduction technique. These functions are the displacement response of the layered soil caused by a concentrated force of unit magnitude (Figure 1(a)). The geometry of the problem is

D 5 D 6 D á T U V . 2 Q ë Q ì Q í 2 (a) (b)
Figure 1: Layered soil: (a) unbounded geometry, and (b) bounded geometry a general 3D. We denote by x = (x, y, z) the position vector. Without loss of generality, the soil domain can be defined as x ∈ R 3 ; z ≤ 0 . Consequently, the soil domain is infinite in the x and y directions and semi-infinite in the z direction. However, the unbounded domain is unsuitable for finite element analysis. Thus, it is imperative to introduce virtual boundaries on the three directions. In this work, the bounded soil is defined by (Figure 1(b)):

Ω =    x ∈ R 3 ; 0 ≤ r ≤ L, r = (x 2 + y 2 ), - n j=1 h j ≤ z ≤ 0    , (1)
where L denotes the horizontal distance that separates the concentrated load to the virtual vertical boundary, n is the considered number of layers, and h j is the height of j th (1 ≤ j ≤ n) layer which is defined by:

Ω j = x ∈ Ω; - j k=1 h k ≤ z ≤ - j-1 k=1 h k . (2)
L defines a range of investigation. Indeed, the displacement induced by the concentrated load drastically vanishes as the distance from the load increases. The range of influence of the load is limited to some region around it. However, the virtual boundaries can induce some numerical discrepancies as they create some non-physical and purely numerical wave reflections. On one hand, L should be increased to delay the effects of the boundary-induced wave reflections. On the hand, L should be reduced to decrease the size of the numerical model and the computation time. Using a reduction model technique is a main advantage here as it permits to enlarge L and the size of the numerical model, thus, solving the tradeoff between reducing computation time and reducing the effects of wave reflections occurring at virtual boundaries. The discussion presented above for L is also true for n and h n . Indeed, to make the vertical dimension of the layered soil bound, it is necessary to limit the number of layers to a bounded number n and/or to bound the height of the deepest layer. We are interested in calculating the displacement field u (x, t) = u (x, y, z, t), over the domain Ω, due to a unit concentrated harmonic force. Without any loss of generality, this load can be applied in x 0 = (x 0 = 0, y 0 = 0, z 0 ). This makes that the horizonal dimension is centered around the load position. In Figure 1, the load is depicted along the vertical direction, but it can also take any of the two horizontal directions. Considering symmetry, quarter of the bounded geometry is modeled (Figure 2) in the case of a vertical force. However, in the case of horizontal force, half of the bounded geometry is modeled (Figure 3). The interface between the j th and (j + 1)

th layers is defined as: Γ z j = {x ∈ Ω; z = -j k=1 h k . We also define the following surfaces:

• Γ z 0 = {x ∈ Ω; z = 0} : the upper surface denoted by [START_REF] Germonpré | Contributions of longitudinal track unevenness and track stiffness variation to railway induced vibration[END_REF] in Figure 2,

• Γ z n = {x ∈ Ω; z = - n k=1
h k } : the bottom surface denoted by (4) in Figure 2, • Γ x = {x ∈ Ω; x = 0} : the surface (3) on Figure 2, • Γ y = {x ∈ Ω; y = 0} : the surface (1) on Figure 2, and

• Γ r = x ∈ Ω; r = L, r = (x 2 + y 2 ) : the surface (2) on Figure 2.
In the case of a model involving half of the geometry, i.e., in the case of horizontal forces, we are not dealing with surface (3). The problem is governed by the set of the following equations

ρ (x) ∂ 2 u ∂t 2 -∇ • σ -f = 0, (3)
σ = λtr (ε) 1 + 2µε, and (4)
σ = ∇u + ∇ T u 2 , (5)
defining the momentum equilibrium, the constitutive law and the small strain definition, respectively. In the above equations ρ is the density, σ is the stress tensor, ε is the strain tensor, λ, µ are the Lamé coefficients, f represents the body forces and 1 is the identity tensor.

The boundary conditions are given by:

• A concentrated harmonic force along the vertical direction F z (x 0 ) = F 0 exp(iωt)

• Symmetric conditions imposed to the displacement field on the boundaries Γ x and Γ y

• Damping conditions on ∂Ω damp = Γ z n Γ r that writes σ • n = η u where n is the outward unit normal vector to the surface. This is mainly introduced to reduce wave reflections from the virtual boundaries Γ z n and Γ r that are introduced here to limit the unbounded geometry of the layered soil.

• Between the layers we assume continuity of the displacement field and of the normal stress vector.

Problem formulation

The equation to be solved is:

ρ (x) ∂ 2 u ∂t 2 -∇ • σ -f = 0. ( 6
)
In the absence of volume forces, the variational formulation on Ω writes

Ω ρ (x) V * ∂ 2 u ∂t 2 -V * (∇ • σ) dΩ = 0. (7)
An integration by part is undertaken. With this integration the domain boundary is decomposed onto two parts: (i) the node on which the concentrated force is applied and (ii) the boundary containing damping conditions. Then we obtain

Ω ρ (x) V * ∂ 2 u ∂t 2 dΩ + Ω ε * : σdΩ- ∂Ω damp V * • (σ • n) d (∂Ω) =V * f F (t) , ( 8
)
where n is the outward normal vector to the damping surface, V * f corresponds to the virtual degree of freedom of the node on which the concentrated force is applied and ∂Ω damp = Γ z n Γ r is the part of the boundary on which we can write that

σ • n=η ∂u ∂t . (9)
Thus, the surface integral of Eq. ( 8) becomes

∂Ω damp V * • (σ • n) d (∂Ω) = ∂Ω damp ηV * • ∂u ∂t d (∂Ω) . (10)
For the sake of clarity while constructing the operators, and without loss of generality, we are going to illustrate the 2D case in the coming equations. The 3D extension will not be fundamentally different. The notations of the kinematic field derivatives will be ∂u ∂t = u and ∂ 2 u ∂t 2 = ü. Moreover, in vectoral notation the stress and strain tensors are given by:

σ=   σ xx σ yy σ xy   (11)
and

ε=   ε xx ε yy 2ε xy   . (12)
Eq. ( 8) is then rewritten as follows:

Ω ρ (x) V * üdΩ + Ω ε * T • σdΩ- ∂Ω damp ηV * • ud (∂Ω)=V * f F (t). (13)
Now we are going to write the constitutive equation. Each soil layer is assumed visco-elastic. Consequently, the stress-strain relation is written:

σ = λ (x, ω) tr (ε) 1 + 2µ (x, ω) ε, (14)
where the Lamé constants λ (x, ω) and µ (x, ω) are depending on the vector position x to account for the different layers of the soil; they are also depending on the solicitation frequency to account for the visco-elastic behavior of each layer. Using a fourth order tensor C(x, ω), the constitutive law of Eq. ( 14) writes:

  σ xx σ yy σ xy   = C(x, ω)   ε xx ε yy 2ε xy   . (15)
C(x, ω) accounts for the visco-elastic behavior of each soil layer. As illustration, this fourth order tensor is expressed in the case of plane strain as:

C(x, ω) =   λ(x, ω) + 2µ(x, ω) λ(x, ω) 0 λ(x, ω) λ(x, ω) + 2µ(x, ω) 0 0 0 µ(x, ω)   . (16)
In order to establish the frequency-space discrete formulation, we need to separate the contribution of the two Lamé coefficients

  σ xx σ yy σ xy   = λ (x, ω)   1 1 0 1 1 0 0 0 0     ε xx ε yy 2ε xy   +µ (x, ω)   2 0 0 0 2 0 0 0 1     ε xx ε yy 2ε xy   . (17) From the displacement field u = u x u y , (18)
we can then express the strain and the virtual strain corresponding to the primary field and the test field as

ε =   ∂ ∂x 0 0 ∂ ∂y ∂ ∂y ∂ ∂x   u x u y , (19)
and

ε * = u * v * ∂ ∂x 0 ∂ ∂y 0 ∂ ∂y ∂ ∂x , (20)
respectively. Let N (x) denotes the space interpolation from the discrete displacement field. The components of the continuum displacement field are given by

u x = N T (x) • u x , (21)
and

u y = N T (x) • u y , (22)
where u x and u y contains the nodal components of the kinematic field which depend here on frequency as we have formulated the elastodynamic problem in the frequency-space domain. Consequently, the strain writes:

ε =   ∂ ∂x 0 0 ∂ ∂y ∂ ∂y ∂ ∂x   . N T (x) • u x N T (x) • u y =     ∂N T (x) ∂x 0 0 ∂N T (x) ∂y ∂N T (x) ∂y ∂N T (x) ∂x     . u x u y .
(23) Moreover, the test field can be written as:

u * v * = u * T v * T N (x) N (x) . ( 24
)
Thus, the virtual strain writes:

ε * = u * T v * T ∂N(x) ∂x 0 ∂N(x) ∂y 0 ∂N(x) ∂y ∂N(x) ∂x . ( 25
)
For the sake of clarity, the notation N (x) will simply be denoted N. Accordingly, the visco-elastic deformation energy within the bounded soil domain is evaluated as follows:

Ω ε * T • σdΩ= Ω u * T v * T ∂N ∂x 0 ∂N ∂y 0 ∂N ∂y ∂N ∂x   λ (x, ω)   1 1 0 1 1 0 0 0 0   + µ (x, ω)   2 0 0 0 2 0 0 0 1        ∂N T ∂x 0 0 ∂N T ∂y ∂N T ∂y ∂N T ∂x    u x u y dΩ.
(26) After some simplifications:

Ω ε * T • σdΩ= u * T v * T Ω λ (x, ω) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩ + Ω µ (x, ω) 2 ∂N ∂x ∂N T ∂x + ∂N ∂y ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂x ∂N T ∂y 2 ∂N ∂y ∂N T ∂y + ∂N ∂x ∂N T ∂x dΩ u x u y . (27)
At this stage, it is very important to define the expressions of Lamé coefficients. As we are going to study the case of harmonic excitation:

F (t) = F 0 exp (iωt) , (28)
the Lamé coefficients are defined using complex number representation in order to take into account the visco-elastic properties of the soil. Namely,

λ (x, ω) = j p ′ j (x) q ′ j (ω) + i j p ′′ j (x) q ′′ j (ω) (29) µ (x, ω) = j r ′ j (x) s ′ j (ω) + i j r ′′ j (x) s ′′ j (ω) (30)
j denotes the layer index. This type of definition allows: (i) to take into account different visco-elastic behaviors for the different layers, and (ii) to separate the space and frequency dependency, which is very important later for PGD formulation. Typically if the behavior is homogeneous within each layer the function p ′ j (x) will depend only on the z direction and will be defined as a rectangular function that is equal to 1 on the layer that behavior is given by q ′ j (ω) and vanish (value = 0) on the other layers. Namely,

p ′ j (x) =    1 if - j k=1 h k < z ≤ - j-1 k=1 h k 0 otherwise (31)
This general definition also allows to take into account the multimodal behavior of each layer. In that case, it is possible to write the function q ′ j (ω) as a sum of different function depending on the frequency ω. The same advantages hold for the other terms p ′′ j (x) q ′′ j (ω), r ′ j (x) s ′ j (ω), and r ′′ j (x) s ′′ j (ω). Accordingly, it is possible to rewrite the terms of the last form of Eq. ( 27) of the variational formulation as:

Ω λ (x, ω) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩ = Ω j p ′ j (x) q ′ j (ω) + i j p ′′ j (x) q ′′ j (ω) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩ (32)
This integral also writes:

Ω λ (x, ω) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩ = j Ω p ′ j (x) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩq ′ j (ω) +i j Ω p ′′ j (x) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩq ′′ j (ω) (33)
After some simplifications [START_REF] Cairo | Analysis of pile groups under vertical harmonic vibration[END_REF] where

Ω λ (x, ω) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩ = j P ′ j q ′ j (ω) + i j P ′′ j q ′′ j (ω)
P ( ′ )( ′′ ) j = Ω p ( ′ )( ′′ ) j (x) ∂N ∂x ∂N T ∂x ∂N ∂x ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂y ∂N T ∂y dΩ (35)
In the same way we can write [START_REF] Zhang | An Efficient Numerical Method for Computing Synthetic Seismograms for a Layered Half-space with Sources and Receivers at Close or Same Depths[END_REF] where

Ω µ 2 ∂N ∂x ∂N T ∂x + ∂N ∂y ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂x ∂N T ∂y 2 ∂N ∂y ∂N T ∂y + ∂N ∂x ∂N T ∂x dΩ = j R ′ j s ′ j (ω)+i j R ′′ j s ′′ j (ω)
R ( ′ )( ′′ ) j = Ω r ( ′ )( ′′ ) j (x) 2 ∂N ∂x ∂N T ∂x + ∂N ∂y ∂N T ∂y ∂N ∂y ∂N T ∂x ∂N ∂x ∂N T ∂y 2 ∂N ∂y ∂N T ∂y + ∂N ∂x ∂N T ∂x dΩ (37)
Finally, we obtain one of the most difficult terms of our variational formulation

Ω ε * T .σ dΩ= u * T v * T j P ′ j q ′ j (ω) + i j P ′′ j q ′′ j (ω) + j R ′ j s ′ j (ω) + i j R ′′ j s ′′ j (ω) u x u y (38
) Now we can move to the other terms that will be easier to develop. The inertial term is given by

Ω ρ (x) V * üdΩ= u * T v * T Ω ρ(x) NN T 0 0 NN T dΩ üx üy ( 39
)
where M denotes the mass matrix and is defined as

M = Ω ρ(x) NN T 0 0 NN T dΩ (40)
This matrix takes into account the variation of the densities through the layers by an appropriate definition of ρ(x). Namely,

ρ (x) = ρ j for - j k=1 h k < z ≤ - j-1 k=1 h k (41)
As this material property does not depend on the frequency the resulting mass matrix is given just by one term (contrarily to the previous stiffness matrix

P ( ′ )( ′′ ) j and R ( ′ )( ′′ ) j
). Finally, we need to define a boundary interpolation functions vector L for the damping boundary conditions. L allows to write:

∂Ω damp ηV * • u d (∂Ω) = u * T v * T ∂Ω damp η LL T 0 0 LL T d(∂Ω damp )
ux uy (42) Therefore, the damping matrix is given by

L = ∂Ω damp η LL T 0 0 LL T d(∂Ω damp ). ( 43
)
Remark: It is worth pointing out that, in what will follow, we are going to use an extended matrix L, the size of which is equal to the total number of degrees of freedom. Mainly, lines and columns with zero values will be added to this matrix corresponding to the degrees of freedom that does not belong to the boundary ∂Ω damp . This allows that the definition of the vector u x u y remains constant in the various terms of the variationnal form.

Let now F 0 be the vector that contains 1 on the excited DoF and zero elsewhere. This allows to write:

V * f F (t) = u * T v * T F 0 exp(iωt) (44)
Finally, if we denote by U = u x u y the total discrete degree-of-freedom vector, then we reach the following system after the simplification of the test fields from the variational formulation

Mü+   j P ′ j q ′ j (ω) + i j P ′′ j q ′′ j (ω) + j R ′ j s ′ j (ω) + i j R ′′ j s ′′ j (ω)   U-L U = F 0 exp(iωt) ( 45
)
The displacement is supposed to follow the complex form

U = U ′ + iU ′′ exp(iωt) (46)
That means that the displacement follows the excitation frequency but exhibits a delay depending on the value of ω. Accordingly, the first and second time derivatives of this field write:

U = -ωU ′′ + iωU ′ exp(iωt) (47)
Ü = -ω 2 U ′ -i ω 2 U ′′ exp(iωt) (48)
We can now inject these expressions in the discrete form and simplify by exp(iωt) to obtain the following complex equation:

M -ω 2 U ′ -iω 2 U ′′ + j P ′ j q ′ j (ω) + i j P ′′ j q ′′ j (ω) + j R ′ j s ′ j (ω) + i j R ′′ j s ′′ j (ω) U ′ + iU ′′ - L -ωU ′′ + iωU ′ = F 0 (49
) Finally, separating the real and imaginary parts of the above complex equation and then gathering the obtained two equations into one matricial system, yields:

-ω 2 M + j P ′ j q ′ j (ω) + j R ′ j s ′ j (ω) -j P ′′ j q ′′ j (ω) -j R ′′ j s ′′ j (ω) + ωL j P ′′ j q ′′ j (ω) + j R ′′ j s ′′ j (ω) -ωL -ω 2 M + j P ′ j q ′ j (ω) + j R ′ j s ′ j (ω) U ′ U ′′ = F 0 0 ( 50
) For a given value of ω the obtained the system is linear. Solving this system directly gives U ′ U ′′ . It is also possible to make the matrix of the above system symmetric by multiplying the second line by (-1).

Remark

In order to post-process the results we are going to extract from the solution U ′ U ′′ some significant values (this will be also the case for post processing each mode of the PGD results of the following sections)

Let denotes by

U ′ x , U ′ y , U ′ z (resp. U ′′ x , U ′′ y , U ′′ z
) the three components of U ′ (resp. of U ′′ ). We can then define an amplitude gain along each axes by

A x = u ′ x 2 + u ′′ x 2 (51)
A y = u ′ y 2 + u ′′ y 2 (52)
A z = u ′ z 2 + u ′′ z 2 (53)
We also can define an angular phase according to

φ x = atan(u ′′ x / u ′ x ) (54)
φ y = atan(u ′′ y / u ′ y ) (55)
φ z = atan(u ′′ z / u ′ z ) (56)
Finally it is possible to make a regeneration of the time dependent solution according to

U (t) = U ′ cos (ωt) -U ′′ sin (ωt) (57)
4 Model reduction using PGD decomposition

Without material parameters uncertainty

As stated before, the size of the finite-element problem is large. Sometimes, it is important to bounded soil domain. It is also important to undertake a parametric study. However, the size of the finite-element problem is prohibitive for such study. In order reduce the size of the numerical model and make computation times reasonable, the numerical model will be here solved within the framework of the PGD (Proper Generalized Decomposition) method. The PGD is a numerical technique that allows circumventing the curse of dimensionality in the resolution of a high dimensional problem [START_REF] Ammar | Proper generalized decomposition of time-multiscale models[END_REF][START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | Recent Advances in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF][START_REF] Ammar | On the Convergence of a Greedy Rank-One Update Algorithm for a Class of Linear Systems[END_REF][START_REF] Ammar | An Error Estimator for Separated Representations of Highly Multidimensional Models[END_REF][START_REF] Chinesta | Proper Generalized Decomposition of Multiscale Models[END_REF][START_REF] Ammar | The Proper Generalized Decomposition: A Powerful Tool for Model Reduction[END_REF]. Especially, the case of parametric study could be found in [START_REF] Aghighi | Parametric solution of the Rayleigh-Benard convection model by using the PGD Application to nanofluids[END_REF][START_REF] Ammar | Solving parametric complex fluids models in rheometric flows[END_REF][START_REF] Ammar | Parametric solutions involving geometry: A step towards efficient shape optimization[END_REF]. Following Malick et al. [START_REF] Haris | Advanced parametric space-frequency separated representations in structural dynamics: A harmonic-modal hybrid approach[END_REF] and Germoso et al. [START_REF] Germoso | Efficient PGD-based dynamic calculation of non-linear soil behavior[END_REF] the elastodynamic problem is formulated in the frequency-space domain. The PGD method is then used to solve the elastodynamic problem considering the frequency as an extra coordinate.

In our case we are going to apply such an idea for solving the problem assuming that the value of ω is unknown. This is possible only because the frequency and space variables are independent, which is a main advantage of the frequencyspace formulation. Indeed, it will be less evident to use the PGD method in a time-space formulation as the time variable and space variable are linked in wave propagation problems through characteristic equations.

In order to prepare a kind of general solution in which the displacement field will depends on the position as well as on the value of ω we must write the problem as:

A A AU U U = B B B, (58)
where

A A A = nA j=1 A j x ⊗ A j ω ⊗ A j c ( 59
)
and

B B B = n b j=1 B j x ⊗ B j ω ⊗ B j c . (60)
Here A j x (resp. B j x ) are operators (resp. vectors) of the three coordinate physical space x, A j ω (resp. B j ω ) are operators (resp. vectors) of the frequential space, A j c (resp. B j c ) are operators (resp. vectors) of a 2 DoF complex space, the first coordinate is the real axis and the second one is the imaginary one. ⊗ stands for the tensor product operator. The solutions are searched in the following form

U U U = nu j=1 f j x ⊗ f j ω ⊗ f j c , (61)
where f j c = c ′ j + ic ′′ j . In order to apply the tensor product resolution technique to find U U U we have just to give the tensor form of operators A A A and B B B. The greedy algorithm for building the solution U U U is quickly and is explained in the appendix. Using the following notations:

W 2 = diag -ω 2 , (62)
W 1 = diag (ω) , (63)
Q ( ′ )( ′′ ) j = diag q ( ′ )( ′′ ) j (ω) , (64)
and

S ( ′ )( ′′ ) j = diag s ( ′ )( ′′ ) j (ω) , ( 65
)
where diag is the operator that build a diagonal matrice from all the discrete values of the included function applied individually on each discrete value of ω, the system of Eq. ( 49) becomes defined by

A A A = M ⊗ W 2 ⊗ 1 0 0 1 + j P ′ j ⊗Q ′ j ⊗ 1 0 0 1 + j R ′ j ⊗S ′ j ⊗ 1 0 0 1 + j P ′′ j ⊗Q ′′ j ⊗ 0 -1 1 0 + j R ′′ j ⊗S ′′ j ⊗ 0 -1 1 0 + L⊗W 1 ⊗ 0 1 -1 0 (66)
which is already under the desired form of Eq. ( 59) and also by

B B B = F 0 ⊗    1 . . . 1    ⊗ 1 0 ( 67
)
which is already under the desired form of Eq. ( 60).

The PGD solver provide the solution as detailed in the Appendix:

U U U = A A A\B B B = nu j=1 f j x ⊗ f j ω ⊗ f j c . (68)
And by post processing we obtain:

U U U ′ = nu j=1 f j x ⊗ f j ω ⊗ c ′ j , (69)
and

U U U ′′ = nu j=1 f j x ⊗ f j ω ⊗ c ′′ j . ( 70
)
We also can rewrite again the solution as

U U U = nu j=1 U j ⊗ f j ω , (71)
where U j are the modes obtained by

U j = f j x ⊗ c ′ j + f j x ⊗ c ′′ j . ( 72
)

With material parameters uncertainty

The soil parameters are not well defined. It is important to analyze the response of the layered soil considering that the parameters may vary within a certain range. The computational cost of this type of study is quite prohibitive unless model reduction techniques are used as the elastodynamic problem should be solved multiple times while varying the material parameter within the desired range. Within the framework of the PGD decomposition, it is possible to consider that each material parameter is a variable playing similar role in the decomposition than the space variable x and the frequency ω.

For sake of clarity, we are going to formulate the parametric problem for a soil with two layers. The procedure remains the same for a different parameterization or higher number of layer. As for Section 4.1, the problem is written as follows:

A A AU U U = B B B (73)
where

A A A = nA j=1 A j x ⊗ A j c ⊗ A j E1 ⊗ A j E2 (74) B B B = n b j=1 B j x ⊗ B j c ⊗ B j E1 ⊗ B j E2 (75)
A j x (resp B j x ) are operators (resp vectors) related to three coordinate physical space x, and A j c (resp B j c ) are operators (resp vectors) related to the 2DoF complex space, the first coordinate is the real axis and the second one is the imaginary one. Also, A j E1 (resp B j E1 ) are operators (resp vectors) related to the parametric coordinate E 1 , and A j E1 (resp B j E1 ) are operators (resp vectors) related to the parametric coordinate E 2 We will look to the solution as

U U U = nu j=1 f j x ⊗ f j c ⊗ f j E1 ⊗ f j E2 ( 76
)
f j c = c ′ j + ic ′′ j (77)
In order to apply the tensor product resolution technique to find U U U we just have to give the tensor form of operators A A A and B B B. The greedy algorithm for building the solution U U U is explained in the appendix. We use the following notations (where diag is the operator that build a diagonal matrice from all the discrete values of E 1,2 )

E 1 = diag(E 1 ) (78)
E 2 = diag(E 2 ) (79)
Moreover, 1 E1 the identity matrix of size equal to the number of discrete values of E 1 and 1 E2 the identity matrix of size equal to the number of discrete values of E 2 .

The operator A A A writes:

A A A = M ⊗ -1 0 0 1 ⊗ 1 E1 ⊗ 1 E1 + L ⊗ 0 1 1 0 ⊗ 1 E1 ⊗ 1 E1 + P ′ 1 q ′ 1 + R ′ 1 s ′ 1 ⊗ 1 0 0 -1 ⊗ E 1 ⊗ 1 E2 + P ′ 2 q ′ 2 + R ′ 2 s ′ 2 ⊗ 1 0 0 -1 ⊗ 1 E1 ⊗ E 2 + P ′′ 1 q ′′ 1 + R ′′ 1 s ′′ 1 ⊗ 0 -1 -1 0 ⊗ E 1 ⊗ 1 E2 + P ′′ 2 q ′′ 2 + R ′′ 2 s ′′ 2 ⊗ 0 -1 -1 0 ⊗ 1 E1 ⊗ E 2 , (80)
which is already under the form of Eq. ( 59); whereas, the right-hand operator B B B writes:

B B B = F 0 ⊗ 1 0 ⊗    1 . . . 1    ⊗    1 . . . 1    , (81)
which is already under the form of Eq. ( 60).

Our PGD solver provide the solution (see the appendix)

U U U = A A A\B B B = nu j=1 f j x ⊗ f j c ⊗ f j E1 ⊗ f j E2 . (82)
We also can (if we want) rewrite again the solution as

U U U = nu j=1 U j ⊗ E j (83)
Where U j are the modes obtained by

U j = f j x ⊗ c ′ j + f j x ⊗ c ′′ j (84)
and E j are the combinations of the two 1D parametric modes into the one 2D parametric space

E j = f j E1 ⊗ f j E2 (85)

Results and discussion

We consider in this result section a two-layer soil domain represented in Figure 3.

The unit system used for forces, lengths, densities, and stresses is the SI system. The first and second layer are 5-m and 20-m thick, respectively. The radius of the bounded is considered L = 50 m. The two layers are assumed to have the same density ρ 1,2 = 1980 kg/m 3 and the same Poisson's ratio ν 1,2 = 0.4903. The Young's modulus of the first and second layer are 446 MPa and 250 MPa, respectively. The damping ratio for the first and second layer are 0.042 and 0.039, respectively. To prevent wave reflections at the virtual boundaries, a dynamic viscosity of η = 10 8 is considered. The harmonic force is defined by F (t) = F 0 exp (iωt) where F 0 = 1N. The investigated frequency range is ω = 2.π.50 Hz. Thus, we can define the functions q

′ 1,2 , q ′′ 1,2 , s ′ 1,2 , s ′′ 1,2 using: λ 1,2 = E 1,2 ν 1,2 (1 -2ν 1,2 ) (1 + ν 1,2 ) , µ 1,2 = E 1,2 2 (1 + ν 1,2 ) ( 86
)
And we consider constant values of the frequential functions

q ′ j (ω) = λ j (87) q ′′ j (ω) = 0.042 λ j (88) s ′ j (ω) = µ j (89) s ′′ j (ω) = 0.039 µ j ( 90
)
Even if for the sake of simplicity these functions are taken constant in this application nothing prevents to make them frequency dependent. The reader can for example looks for such function as the real and complex modulus of the Maxwell visco-elastic model given by

q ′ j (ω) = G ′ ω 2 θ 2 j 1 + ω 2 θ 2 j (91) q ′′ j (ω) = G ′′ ωθ j 1 + ω 2 θ 2 j ( 92
)
Similar expression could be used for s ′ j (ω) and s ′′ j (ω). But this requires an experimental identification of the model parameters G ′ , G ′′ and the relaxation time θ j .

Example of non parametric results emphasizing the interest of the frequential formulation of Eq. (50)

The harmonic force is defined on various depths and could be vertical (along z direction) or horizontal (along x direction). The only symmetry plane of the problem is the plane defined by y = 0. Eight cases are investigated and are detailed in Table 1. Figure 3, shows the position of some points on which we have plot the time dependent displacement. Table 1: The eight studied cases 3, 5 and 7, respectively, show the dominance of the horizontal displacements of points 1, 7, 8 and 9, which are the closest to the respective points where the force is applied. Likewise, Figures 5(b), 5(d), 5(f) and 5(h) corresponding to the cases 2, 4, 6 and 8, respectively, show the dominance of the vertical displacements of points 1, 7, 8 and 9, which are the closest to the respective points where the force is applied. However, Figures

and 4(h) corresponding to the cases 2, 4, 6 and 8, respectively, show that the vertical displacements vanish as the corresponding forces are applied in the perpendicular horizontal direction. However, Figures 5(a), 4(c), 4(e) and 4(g) corresponding to the cases 1, 3, 5 and 7, respectively, show that the horizontal displacements vanish as the corresponding forces are applied in the perpendicular vertical direction.

Illustration of the PGD resolution with material parameters uncertainty

An example of result of the tensor form established in equations 24-25 is given in Figures 6789. The harmonic force is defined by an horizontal force F (t) = F 0 exp (iωt) with ω = 2.π.50. We consider that for each layer we have different possible values for young modulus

E 1 ∈ [350 • 10 6 , 450 • 10 6 ], E 2 = [200 • 10 6 , 300 • 10 6 ] (93)
The figure 6 represents the decrease of the weight of the modes. This information is very useful for the comparison of the PGD solution to the full grid solution of the same resolution (using the full tensor product of the FE operators by the material operators). In fact the decrease of the relative value of weight of the 'n th ' enrichment related to the first one implicitly represents the relative decrease of the error of the PGD in relation to the full grid solution. This is the indicator used to estimate the convergence of the PGD solution during the successive enrichments. The functions U j and their associated E j are represented in terms of gain for the first three modes in figures 7-9. The interest of such a representation is mainly to obtain an offline solution. When the user gives the value of E 1 and E 2 then the results comes defined by a linear combination of modes according to

U = nu j=1 α j U j , α j = E j (E 1 , E 2 ) (94)
To emphasize on the time saving obtained by using the PGD method in this case we consider the example where the material parameters are considered as parametric dimensions. Let's take the case where we are interested on problem involving n E1 points in the parametric axis of E 1 and n E2 points in the parametric axis of E 2 . In classical resolution to cover all this space we need n E1 . n E2 resolutions of the linear system (Eq. ( 50)) that involves 6N degrees of freedom where N is the number of nodes. When the PGD is used we have a set of linear system to solve at each iteration of the fixed-point strategy defined by the equation A.16 related to each dimension. The operator of this equation is defined in A.13 and has a size equal to 1. 3N : for the first dimension (which is the spatial dimension)

2. 2 for the second dimension 3. n E1 for the third dimension 4. n E2 for the fourth dimension

In that case the system that will determine the CPU time will be the one related to the spatial dimension. All the other resolutions can be neglected in relation to this one. In general, the fixed point converges in about 5 iterations and the enrichment strategy reaches the global convergence criteria (A.17) within about 10 enrichments. Thus the global cost of the PGD can be estimated to 50 resolutions of 3N dof linear system. This cost is generally much lower than n E1 . n E2 resolutions of 6N dof linear system. But this is not the only advantage of the PGD decomposition. In that case, the PGD resolution allows providing an off-line solution, i.e. a solution that is calculated only once. When the user precise the parameters values the particularization of the solution can be done in real time. This benefit opens up a lot of prospects for touch app development on smartphones for example.

Illustration of the PGD resolution where the frequency is considered as parameter

With the same parameters previously used we work now with PGD decomposition where material parameters are fixed but where frequency is considered as a parameter (see Eq. [START_REF] Aubry | Periodic BEM and FEM-BEM coupling: applications to seismic behaviour of very long structures[END_REF][START_REF] Chebli | Response of periodic structures due to moving loads[END_REF]. An example of result is given in figure 10 and 11. The functions U j are represented in terms of gain for the first four modes in figure 10. In figure 11 we represent f j=1..6 ω . The interest of such a representation is mainly to obtain an offline solution. When the user precise the particular value of ω = ω p then the results comes defined by a linear combination of modes according to

U = nu j=1 α j U j , α j = f j ω (ω p ) (95)
As for the previous case the right plot in figure 11 represents the decrease of the weight of the modes, useful information for estimating the convergence of the PGD successive enrichments. Remark: The different PGD functions U j (which are sometimes called modes) do not actually represent the modes relating to a given frequency. But the solution sought for each frequency is a linear combination of these different functions using the factors f j ω (ω p ).

Conclusion

In this work, we propose a general 3D framework able to give the harmonic response of exited visco-elastic soil within a space-frequency formulation. This formulation allows getting the harmonic solution without any explicit or implicit time integration scheme. 3D results show that we are able to characterize the response for a given frequency in terms of gain, angular phase and reconstituted time dependent displacement field. Displacement of points of interest in our structure (such are points on which constraints on displacement may be imposed) can be accurately represented during time in order to quantify their associated amplitude. The PGD method was also successfully formulated and validated in the frequency-space domain. It can be used for the parametric study where the frequency is a priori unknown. Finally, this new approach allows for a parametric study with a reasonable time computation cost.

Acknowledgement

A General framework of tensor product resolution

Let Ω be a multidimensional domain involving N coordinates (x 1 , x 2 , . . . , x N ) (where each coordinate x j is not necessarily one dimensional). Let consider a weak form of a linear problem given by : a (Ψ (x 1 , x 2 , . . . , x N ) , Ψ * (x 1 , x 2 , . . . , x N )) = b (Ψ * (x 1 , x 2 , . . . , x N )) , (96)

where we are looking to an approximated solution that writes in the continuous form as:

Ψ (x 1 , x 2 , . . . , x N ) = nF j=1 α j F j 1 (x 1 ) F j 2 (x 2 ) . . . F j N (x N ) . ( 97
)
The separated representation is built-up from a projection-enrichment iteration scheme. To this goal, we need to write the continuous expression in a discrete form using the nodal values of each function. This discrete form is given by:

Ψ = nF j=1 α j F j 1 ⊗ F j 2 ⊗ . . . ⊗F j N . (98)
Moreover, it is important to transform the weak form into a discrete form written as:

Ψ * AΨ R = nA k=1        R * T j A k j R j N h = 1 h = j R T h A k h R h        = R * T j KR j (108)
Ψ * AΨ F = nF i=1 nA k=1        α i R * T j A k j F i j N h = 1 h = j R T h A k h F i h        = R * T j V, (109)
and

Ψ * B = nB k=1        R * T j B k j N h = 1 h = j R T h B k h        = R * T j V ′ (110)
In the framework of a fixed-point strategy we are looking, at each iteration, for the solution of the linear system:

KR j + V = V ′ . ( 111
)
We assume that the global convergence is attained when the error

ε = A Ψ -B 2 . (112)
becomes small enough. In terms of minimization, the previous variational formulation is equivalent to the minimization of the next functional

J (T ) = T -A -1 B T A T -A -1 B = T -A -1 B A (113)

Figures 4 and 5

 5 depict the horizontal and vertical displacements, respectively, in the selected points. The displacement is important in the closest point to the applied force and in the same direction as it. The displacement almost vanishes elsewhere. Hence, Figures4(a), 4(c), 4(e) and 4(g) corresponding to the cases 1,

Figure 2 :

 2 Figure 2: Geometry of the studied problem

Figure 4 :

 4 Figure 4: Horizontal displacement: (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (e) Case 5, (f) Case 6, (g) Case 7, and (h) Case 8

Figure 5 :

 5 Figure 5: Vertical displacement: (a) Case 1, (b) Case 2, (c) Case 3, (d) Case 4, (e) Case 5, (f) Case 6, (g) Case 7, and (h) Case 8

Figure 6 :

 6 Figure 6: Decrease of the modes weight during the PGD enrichment.

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 7: Mode number 1 -Space mode represented in terms of gains and parametric mode in the space (E1, E2)

5 Figure 10 :Figure 11 :

 51011 Figure 10: PGD representation of the frequency parametric solution. From top to down: the U j functions associated to the fifth first modes (j = 1..5). For each mode we show the three gains along X, Y and Z (on each column)

  This work was supported by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under the grant No. (D1439-40-135). The authors, therefore, gratefully acknowledge the DSR technical and financial support.

where

and

The projection stage consists on finding the best set of α j coefficients. The associated test function is given by

Substituting Eqs. ( 98), (100), ( 101) and (102) into Eq. (99) yields:

where

and

The enrichment stage includes new candidates for enriching the reduced separated approximation basis such as:

Within a fixed-point alternating direction algorithm, we look at each iteration for the computation of a single discrete function R j assuming all the others known. Thus, when we are looking for R j the test function writes:

The different terms of the discrete weak form read: