Assessing the Resilience of Socio-Ecosystems: Coupling Viability Theory and Active Learning with kd-Trees. Application to Bilingual Societies

1 Introduction

Assessing the resilience of an ecological or social system is becoming a challenge in the context of sustainable development (as stated in [Perrings, 2006]). Traditionally resilience measures the ability of a system to recover after a perturbation (see [Martin, 2004] for an analysis of operational definitions of resilience). In this context the mathematical theory of viability is an interesting framework since it studies the compatibility of dynamical systems and constraints [Aubin et al., 2011]. This framework is used in sustainability studies in order to find control policies that keep the system in a given constraint set, such as the concept of tolerable windows [Bruckner et al., 2003] for climate change studies. When the system evolves outside the desirable constraint set, viability theory can assess whether and how the system can be driven back to desirable states. In particular it can provide the returning time [Doyen and Saint-Pierre, 1997]. In this paper we quantify the resilience as the inverse value of this returning time, as it is done in [Martin, 2004] for ecological systems.

∗The research leading to these results has received funding from the E.C.'s FP7/2009-2013 under grant agreement DREAM n. 222654-2.

Viability analysis is based on the computation of the viability kernel and its capture basin. The viability kernel gathers all the states from which there exists a control function that allows the evolution to stay in the constraint set. Its capture basin gathers the states from which it is possible to reach the viability kernel in finite time. Unfortunately, the complexity of the computation task is exponential with space or time when using a discrete approximation on a grid [SaintPierre, 1994]. Therefore this method is limited to very low dimension (at most 3) or to linear models, so its applicability is very much impaired. Moreover, the complexity of the model in real application can limit the practical use of the method because of computation problems when running the model. Dimension 4 was exceptionally reached in a real food processing application with the help of grid computing (see [Sicard et al., 2012], [Reuillon et al., 2010]).

In the particular case of finite time horizon, [Bonneuil, 2006] proposes a method based on simulated annealing to generate trajectories near the boundary. This method works in higher dimensions but it only produces a set of trajectories so it cannot be used to compute the viability kernel nor resilience level sets. Following the same idea of focusing on the boundary in the general case, [Deffuant et al., 2007] introduces classification functions in Saint-Pierre's viability algorithm in order to reduce the number of calls to the model and to represent the viability set in a more compact way than the traditional regular grid. But the Support Vector Machine (SVM) functions used as classification function in [Deffuant et al., 2007], although very attractive, don't fulfill the conditions of the convergence theorem. So the output sets produced

F-63172 Aubiere Cedex, France

Sophie Martin (3) (3) Irstea, LISC sophie.martin@irstea.fr with the SVM method are not reliable. To overcome these limitations, we also consider the computation of the boundary at each step as an active learning problem, but we use kdtrees as it is done in [Rouquier et al., ] to limit the number of calls to the model. Using kd-trees as classification functions as well as oracle functions to store the successive approximations of the viability sets, it is possible to converge towards the true viability kernel or capture basin.

The paper is organized as follows: Section 2 presents the language competition problem and its translation in the viability framework. We define the resilience and show how it is related to the capture basin defined by the theory. Section 3 is devoted to the computation of the capture basin. It explains how it can be seen as an active learning problem, and it describes the kd-tree algorithm that is used in that purpose. It also gives the proof of convergence of the algorithm towards the output sets of the viability theory. Section 4 presents the result of the computation for the language competition problem and how the computed resilience value can be used to propose action policies that guarantee the language coexistence.

2 Resilience of bilingual societies

Bilingual societies exist, but history tells us that many languages have disappeared. Presently, many languages seem endangered. Some researchers have proposed mathematical models to study language competition (see for instance [Abrams and Strogatz, 2003]). These models are consistent with historical data from past or present endangered languages (such as Welsh, Scottish, Gaelic, Quechua). Models without bilingual population generally conclude that one language eventually becomes extinct. When a bilingual population exists, the survival of the language diversity depends on the relative prestige of languages [Bernard and Martin, 2012].

But the question of resilience of these societies is still pending. Real societies are subjected to migration: both immigration and emigration (due to economic, environmental, educational or safety reasons for instance) disturb the dynamics of the different language speakers and can jeopardize the diversity of a bilingual society. This is the reason why it is interesting to study the ability of a bilingual society to recover from such a perturbation. Here we propose to study this resilience in the viability framework.

2.1 Model of bilingual society

The model from [Bernard and Martin, 2012] describes a bilingual society considering three groups in the population: the monolingual speakers of language A, the monolingual speakers of language B, and the bilingual speakers AB. The model is described by σA the proportion of speakers of A, σB the proportion of speakers of B, and the proportion of bilinguals σAB with σAB = 1 - σA - σB. The rate at which speakers of one language switch to become speakers of the second language depends on the attractiveness of this second language, but transitions between groups concern A or B → AB and conversely. Since it is unlikely for one speaker to switch language from scratch, we have PA→B = PB→A = 0. We then have:

dσA dt = σABPAB→A -σAPA→AB, (with analogous equations for dσB dt and dσAB dt ).

Concerning transitions, monolinguals are sensitive to the size of the monolingual group of the other language when bilinguals are sensitive to the whole group of speakers of the other language. The proportionality factor is the prestige of one language compared to the other (we note s the prestige of language A, so 1 - s is the prestige of language B). A parameter a models how the attractiveness of a language scales with the proportion of speakers. So we have:

PAB→A = (1 - σB)a s and PA→AB = σaB (1 - s). The value of the prestige, s, can evolve with public action: education, advertising, incentive policy, arranged employment, positive discrimination, and so on. Since these kinds of policy take time to give results, ds dt is considered here as a control on the dynamics, bounded in some interval U = [-¯u; ¯u] with ¯u > 0. In this framework the model of a bilingual society is defined by the following equations:

If s is constant in ]0; 1[, the dynamics has three equilibria: (0, 1), (1, 0) which are stable, and an unstable one. Consequently, one language is doomed to become extinct. Therefore it is necessary to apply control policy on s to insure the coexistence.

2.2 Set of desirable states and viability domain

In order to ensure the sustainability of the language diversity, it is necessary to define what is considered as desirable. For instance, an easy way to define language diversity would be to consider that the proportion of monolingual speakers of each language has to be above a threshold σ. The set of desirable set in E = [0, 1]3 is then K = {(σA, σB, s)} such that:

We suppose now that this threshold is set to 0.2. We consider that the boundary of the control variation ¯u is set to 0.1. Parameter a set to 1.31 according to the literature. (It is calibrated in [Abrams and Strogatz, 2003] from historical data).

In this framework the viability theory applies. It can be shown that there exists a subset of K which is a viability domain D: from each state x = (σA, σB, s) in D there exists a control function u(t) that allows the evolution of the dynamics 1 to stay inside D. Figure 1 shows the shape of D.

2.3 Capture basin and resilience value

Outside the viability domain D, the bilingual society may not be viable any longer. In particular, the proportion of one of the languages may fall under the threshold σ = 0.2, which is not considered as a satisfying future. 1

1 Since there exists a viability domain in K, the viability theory states that there must be a viability kernel including D, but actually it

We consider that without additional information, the viability domain D gathers the secure states, since from any starting point in D, we are sure that the evolution can remain in the set of desirable states K (that is, each language is spoken by at least 20% of the population) with the appropriate action policy.

In case of perturbation, the state of the bilingual society can jump from D to a state outside the viability domain. Here we consider that perturbations can affect each variable: Migrations change the value of the proportion of speakers of each language. But the prestige of one language may also be disturbed, for example by criminal cases or impressive achievement concerning one group. In a sustainability viewpoint, the issue is to know if it is possible to reach a secure position in the viability domain, starting from the disturbed state. The set of states from which there exists a control policy that drives the system back to D is called the capture basin of D in the viability theory.

At each time horizon T, it is possible to consider the corresponding capture basin Capt(D, T) = {x | ∃u(.), x(T) ∈ D}. Its boundary ∂Capt(D, T) is the level set of the returning time T to D. In particular we have by definition: Capt(D, 0) = D.

Definition 1 (Resilience of a state) The resilience of a state x outside the viability set D is r(x) = 0 if x /∈ Capt(D). Otherwise, it is r(x) = 1 T with T = min{t| x ∈ Capt(D, t)}.

Definition 2 The resilience of a bilingual society to a given perturbation set P is: min{r(x + y)| x ∈ D, y ∈ P}.

When the capture basin is not greater than the viability domain, then the resilience of the bilingual society is zero: In

could be D itself. In this latter case, outside D, one of the languages is doomed to be spoken by less than a proportion σ = 0.2 of the population.

case of a perturbation, the proportion of speakers of one language will never recover. If the capture basin is greater than the viability domain D, that means that if the state of the society is not as expected yet but still in the capture basin, it can return to the set of desirable states with an appropriate control policy. Computing the capture basin of D will allow to measure resilience and to find these control policies.

3 Computation of the resilience

3.1 Computation of the capture basin: active learning sampling with kd-tree

The computation of the capture basin of the viability domain D can be seen as a particular problem of viability [Aubin, 2001]. Actually, the capture basin can be defined as a projection of the viability kernel of an auxiliary viability problem with one additional dimension, the time dimension. So the same computational limits apply. The classical approach consists in sampling the state space along a grid, and store the coordinates of all viable points, but this is costly in terms of time and memory.

We consider the dynamical model as a black box F defined from E = [0, 1]3 to E = [0, 1]3 . F is the discrete dynamical system corresponding to the model 1. Based on this black box function F , an exact function f(S)(x) from E = [0, 1]3 to {0, 1} is built, which can tell as an oracle whether there exists or not a control u ∈ U such that the trajectory starting at point x with control u belongs to a given set S at the following time step. If this set S is Capt(D, T) then x ∈ Capt(D, T + dt) when f(x) = 1, and x /∈ Capt(D, T + dt) when f(x) = 0.

As in [Rouquier et al., ], we use kd-trees to discover the boundary of Capt(D, T + dt). The kd-tree subdivides the search space E. It focuses on the boundary by construction, so it reduces the number of calls to the oracle f, compared to the extensive search on the grid.

We first build a kd-tree gdt to approximate f(D). So gdt is an approximation of the characteristic function of Capt(D, dt). We then recursively build a kd-tree gT+dt that approximates fg(T) . By construction g(T+dt) is an approximation of the characteristic function of Capt(D, T + dt).

3.2 Algorithm

The kd-tree g = gT is built by refining its root node which contains the whole space E. To this root node is associated a point x0 of D such that an evolution starting at x0 is still in D at time step T . Since D is a viability domain, any point in D is suitable. The root node is labeled to 1, since its associated point x0 is in Capt(D, T).

Nodes that are refined are always divided in two halves along their largest dimension, unless they meet the stopping criterion: a node is not divided any more if the length of its largest side is already h. These are called terminal nodes. With this splitting method, nodes are represented by a vector of half-open intervals (except for border nodes). Since intervals are split in the middle, they can't overlap, they are either disjoint, adjacent or included in one another. Nodes that have to be refined are:

• Border nodes (nodes that are adjacent to the boundary of the search space E) with label 1.

• Adjacent nodes with different labels (they are called critical pairs).

Two nodes are adjacent if and only if on each coordinate axis, one interval is included in (or equal to) the other, except for one coordinate where both intervals are adjacent (their closures have exactly one common point). When a node is refined, a point p is drawn at random (other method can be used) in each child node which is labeled with f(p). .

The sketch of the algorithm is the following:

Algorithm 1 BuildG(root, h)

0. while count ̸= 0 do { 1. count ← 0 2. A←NodeToDivide(root) 3. foreach nodei ∈ A do { 4. if nodei is not a terminal node then do { 5. count ← count + 1 6. Refine(nodei, h) 7. A←A\nodei (a node is refined only once) }}

1. return g = root

Refine(node, h) creates two children at node node, by splitting it along its largest dimension if its greater than h. If node is represented by the vector of intervals ([ai, bi[)1,..,n and is divided along j, its children will be have the same vector values except for the j th coordinate which is [aj, (bj-aj) 2 [ for one child and [ (bj-aj) 2 , bj[ (if bj ̸= 1) for the other ([ (1-aj ) 2 , 1] with bj = 1). A call to the oracle is made for the child that need it in order to attribute a label to the node. (Other methods may require 2 calls).

Algorithm 2 NodeToDivide(node)

1. if node is a leaf then
2. return result }

1. if f(node) = 1 AND node is a border then return {node} 2. else return(∅) 3. else { (node1, node2) ← node.children 4. result ← ∅ 5. result ← result ∪ NodeToDivide(node1) 6. result ← result ∪ NodeToDivide(node2) 7. result ← result ∪ PairsInNodes(node1, node2)

Algorithm 3 PairsInNodes(node1, node2)

node1 and node2 are necessarily adjacent

0. result ← ∅

1. if both node1 and node2 are leaves then
2. if f(node1) ̸= f(node2) then
3. result ← {node1, node2}
4. else if neither node1 nor node2 are leaves then
5. foreach nodei child of node1 do {
6. foreach nodej child of node2 do {
7. if nodei and nodej are adjacent then
8. result ← result ∪ PairsInNodes(nodei, nodej ) }} 9. else do {
9. leaf ← the leaf (either node1 or node2)
10. node← the other node
11. border ←axis and point of adjacency
12. label ← 1 - f (leaf )

14. foreach nodei ∈ BoundaryNodes(node, border, label) do {

15. if leaf and nodei are adjacent then

16. return result

BoundaryNodes(node, border, label) returns all the leaves l within node such that f(l) = label and that share with node the same boundary border.

The characteristic function is built by calling NodeToDivide on the search space E.

3.3 Convergence

The algorithm refines a subtree (generally the root) by splitting leaves in two, so it reaches the stopping criterion in finite time.

We show here that the series of sets defined by the kd-tree algorithm converges towards the viability kernel of the associated viability problem. The method we use verifies the conditions of the convergence theorem 1 proposed in [Deffuant et al., 2007]. This theorem applies when a discrete viability kernel is approximated with a classification procedure.

A general viability problem is defined by a controlled dynamics in E ⊂ Rp modeled by a set of equations such as (1), and a compact set of constraints K ⊂ E, such as (2) in which we want the dynamics to evolve:

The viability kernel is the subset of states from which there exists a control function u(t) that allows the evolution of the dynamics to stay inside K: Viab(K) = {x ∈ K | ∃u(.) | x(t) ∈ K ∀t ∈ [0,+∞[}. Considering a time interval dt, the discrete dynamical system F associated to 3 is a correspondence E ↦→ E which associates to a state x its set of successors. We assume that F is µ-Lipschitz, which means that the images of two vectors x and y can't diverge from more than µd(x, y):

As for the language competition model, we consider that an oracle f is available, to tell whether at the following time step a state evolves towards a given set S. In order to guarantee the convergence toward the viability kernel it is necessary to consider an augmented set; The oracle is assumed to be defined from E to {0, 1} in the following way: ˆ f(S)(x) = 1 if and only if there exists u ∈ U(x) such that d(F(x), S) ≤ µβ(h), otherwise ˆ f(x) = 0 (where β(h) depends on the grid and tends to 0 as h tends to 0). This means that the oracle ˆ f is based on the characteristic function of S ∪ (∪x∈SB(x, µ)).

We want to prove that algorithm 1 converges to the characteristic function of the viability kernel viabF (K) of the discrete dynamic (4) when its stopping criterion h tends to 0 and when the viability kernel has sufficiently smooth properties.

Theorem [Deffuant et al., 2007] assures that the series of sets built with a learning procedure verifying conditions (5)

and (6) converges towards the viability set. In order to apply this theorem, we consider a discrete grid Kh induced by the stopping criterion h. For example, the points of the grid can be centered on each terminal leaf of a totally refined tree. We obviously have: ∀x ∈ K, ∃xh ∈ Kh, such that d(x, xh) ≤ β(h), with β(h) → 0 when h → 0. (Actually β(h) = h2 √p). In the following we omit the stopping criterion index h when it is not necessary. At each time step we define a discrete set Kn+1 ⊂ Kn ⊂ Kh, and we build a kd-tree g on the sets Kn+1 , Kh\Kn+1 with algorithm 1. Kn+1 gathers the points of Kn such that the dynamics (4) send them near g(Kn ) (at a distance less than µβ(h)). We note gn the kd-tree built with the subsets of Kn and Kh\Kn that correspond to points xh for which the oracle ˆ f has been called when refining the tree. We note Gn = {x ∈ K | gn (x) = 1}. We have K0 = Kh and G0 = K. Then:

Since Kn+1 ⊂ Kn ⊂ Kh it converges to a limit Lh. We note gh the kd-tree built on the subsets Lh and Kh\Lh, and Gh = {x ∈ K | gh(x) = 1}.

Theorem 1 (fact) If the learning procedure g verifies the following conditions (5) and (6), then Gh converges to Viab(K) when h → 0.

We consider that the viability kernel viabF (K) (and its complementary set in K) verify the following conditions: We note Vϵ(S) = {x ∈ S | B(x, ϵ) ⊂ S} (V is an open space in S)

The path-connectivity of the viability kernel and its complementary are insured when K is simply connected and sufficiently regular. Conditions 8 insure that when the stopping criterion h is small enough, then there are no more thin tentacle that can't be seen through the grid.

Theorem 2 The series of characteristic functions defined by the kd-trees build with algorithm 1 converges to the characteristic function of the viability kernel when the stopping criterion tends to 0.

Proof: The learning procedure g based on algorithm 1 verifies the conditions of Theorem 1.

We first prove condition (6). Let x be a state in K such that gn+1 (x) = 0. This means that x belongs to a leaf L whith labeling point xh ∈ K\Kn+1 . If L is a terminal node, then it is divided to the stopping criterion, and we have: d(x, xh) ≤ β(h). If the leaf L of x is not a terminal node, let us consider the terminal node N to which x would belong if the node L was divided into a subtree to the stopping criterion. Let yh be the labeling point of N: we have d(x, yh) ≤ β(h). We also have N ⊂ L. Since L is not a terminal node, the closest labeling point y of a terminal node with label 1 is at least separated from L by another terminal node with label 0. So necessarily yh ̸= y and so yh ∈ Kn+1 and condition (6) is fulfilled. Since the kd-tree algorithm 1 proceeds in the same way whatever the label of the nodes, the condition (5) is also fulfilled. ·

The oracle ˆ f is defined with the following procedure: Let g be a kd-tree build with algorithm 1. g is the characteristic function of G. We build a kd-tree ˆ g by changing the label of critical leaves to 0, then refining the tree again with algorithm 1. This adds a level of new critical leaves with label 0 (refined from larger nodes with label 0), that are adjacent to the leaves whose label have been modified. If we note G = {x ∈ K, g(x) = 1} and H = {x ∈ K, ˆ g(x) = 0}, we have: h√p ≥ d(G,H) ≥ h. This operation is repeated k times until kh ≥ µβ(h). So the use of ˆ g(k) is a bounded approximation of the oracle ˆ f. This guarantees the convergence of the process to the viability kernel.

4 Preservation of language diversity

4.1 Full resilience of bilingual societies

We used algorithm 1 in order to compute an approximation of the capture basin of the viability domain D of the language competition problem. 2

From the viability domain D, the algorithm computes the successive approximations of the capture basin. Figure 2 shows the capture basin of D at different time horizons. The boundaries are the level sets of the resilience.

2 For this experiment we use directly f as the oracle function and not ˆ f, and the points associated to the leaves are drawn at random. It reduces the computation time (but the convergence is no longer strictly decreasing).

Figure 3 shows the evolution of the volume of the different approximation of the capture basin in language competition problem when increasing the time horizon. Since we necessarily have σA + σB ≤ 1, the volume of the state space is 0.5. We can see that the volume of the capture basin actually converges to 0.5. This means that bilingual societies are fully resilient: the capture basin of the viability domain covers the entire space. So whatever the perturbations, there

always exists a control policy that allows to return in the set of desirable states in less than 25 time steps, so the global resilience is r = 0.04.

With this information it is possible to accept temporary undesirable situations from the language diversity viewpoint, such as receiving refugees, which can shift the present state in D to a state outside D. For instance if the initial state is (σA = 0.28, σB = 0.48, s = 0.40) it can be decided to accept a disturbance to the state (σA = 0.22, σB = 0.54, s = 0.40) since its returning time step is only 4.

4.2 Choosing a policy after a disturbance

If a disturbance is inescapable (for instance aging of population, etc.), the computation of the capture basin provides control functions in order to return to the set of desirable states D. For example, we consider a bilingual society which state is supposed to be (σA = 0.35, σB = 0.35, s = 0.2). This means that the proportion of bilingual speakers is 30%. This state is inside D, which means that an appropriate control of s can maintain the bilinguism property. We suppose now that a perturbation occurs, like a large immigration of speakers of language A. It changes the proportion of speakers of language B from 0.35 to 0.41 and the proportion of speakers of language A to 0.30. The new state P = (0.30, 0.41, 0.2) is not in D. P is very close to D but even with the maximal control available ds dt = 0.1, the dynamics drives the system outside the constraint set K (at time step 3), as it can be seen on figure 4. Since the system is resilient, we already know that it is possible to return to the set of desirable states K. The resilience of P is r(P) = 0.14, which means that it can return in 7 time steps to point Q = (0.2, 0.61, 0.90).

5 Conclusion

We have presented in this article a new method to compute the level sets of the resilience defined in the framework of the viability theory. This method uses kd-trees with a stopping criterion h to store the characteristic function of the capture basin of the viability kernel (which provides the level set of the resilience). It uses a slightly modified kd-tree to approximate the oracle function that allow to compute these sets. The number of calls to the oracle function is limited to the number of nodes of the kd-trees. The building function concentrates these calls at points near the boundary of the viability set. We have proved that this algorithm converges towards the viability kernel when the stopping criterion h tends to 0. Moreover, the algorithm that codes the dynamics is a black box for the algorithm that build the viability kernel (or the capture basin). This modular aspect guarantees the reusability of our approach to other dynamics (the code is available as free software). We have shown how this method can be used in a language competition problem. We have shown how a resilience value can be defined. We have computed the level sets of the resilience and we have shown how it can be used to propose action policy in order to guarantee the language coexistence. This method focuses on the boundary of the viability sets rather than on the sets themselves, so it allows to consider state space with one additional dimension But it still suffers the curse of dimensionality. Further work is in progress in order to parallelize the refining part of the algorithm. Another limit is that in order to provide efficient and realistic action policies it is necessary to know the distance of a state to the boundary of the viability kernel. The algorithm can easily provide an approximation of this distance, so further work also concentrates on this point.