Pierre Baraduc
email: pbaraduc@snv.jussieu.fr

Emmanuel Guigon
email: guigon@ccr.jussieu.fr

Where does the population vector of motor cortical cells point during reaching movements?

Visually-guided arm reaching movements are produced by distributed neural networks within parietal and frontal regions of the cerebral cortex. Experimental data indicate that (I) single neurons in these regions are broadly tuned to parameters of movement; (2) appropriate commands are elaborated by populations of neurons; (3) the coordinated action of neurons can be visualized using a neuronal population vector (NPV). However, the NPV provides only a rough estimate of movement parameters (direction, velocity) and may even fail to reflect the parameters of movement when arm posture is changed. We designed a model of the cortical motor command to investigate the relation between the desired direction of the movement, the actual direction of movement and the direction of the NPV in motor cortex. The model is a two-layer self-organizing neural network which combines broadly-tuned (muscular) proprioceptive and (cartesian) visual information to calculate (angular) motor commands for the initial part of the movement of a two-link arm. The network was trained by motor babbling in 5 positions. Simulations showed that (1) the network produced appropriate movement direction over a large part of the workspace; (2) small deviations of the actual trajectory from the desired trajectory existed at the extremities of the workspace; (3) these deviations were accompanied by large deviations of the NPV from both trajectories. These results suggest the NPV does not give a faithful image of cortical processing during arm reaching movements.

INTRODUCTION

When reaching to an object, our brain transforms a visual stimulus on the retina into a finely coordinated motor act. This complex process is subserved in part by distributed neuronal populations within parietal and frontal regions of the cerebral cortex [START_REF] Kalaska | Cerebral cortical mechanisms of reaching movements[END_REF]. Neurons in these areas contribute to coordinate transformations by encoding target position and kinematic parameters of reaching movements in multiple frames of reference and to the elaboration of motor commands by sending directional and positional signals to the spinal cord (Georgopoulos 1996). An ubiquitous feature of cortical populations is that most neurons are broadly tuned to a preferred attribute (e.g. direction) and that tuning curves are uniformly (or regularly) distributed in the attribute space [START_REF] Georgopoulos | On the translation of directional motor cortical commands to activation of muscles via spinal interneuronal systems[END_REF]. Accordingly, a powerful tool to analyse cortical populations is the NPV which describes the behavior of a whole population by a single vector [START_REF] Georgopoulos | On the translation of directional motor cortical commands to activation of muscles via spinal interneuronal systems[END_REF]. [START_REF] Georgopoulos | Neuronal population coding of movement direction[END_REF] have shown that the NPV calculated on a set of directionally tuned neurons in motor cortex points approximately (error f"o.J 15°) in the direction of movement. However, the NPV may fail to indicate the correct direction of movement when the arm is in a particular posture [START_REF] Scott | Changes in motor cortex activity during reaching movements with similar hand paths but different arm postures[END_REF]. These data raise two important questions: (1) how populations of broadly tuned neurons learn to compute a correct sensorimotor transformation? Previous models (Burnod et al. 1992; Bullock et al. 1993; Salinas and Abbott 1995) provided partial solutions to this problem but we still lack a model which closely matches physiological and psychophysical data on reaching movements; (2) Are cortical processes involved in the visual guidance of arm movements readable with the NPV tool? This article provides answers to these questions through a physiologically inspired model of sensorimotor transformations.

MODEL OF THE VISUAL-TO-MOTOR TRANSFORMATION

ARM GEOMETRY

The arm model has voluntarily been chosen simple. It is a planar, two-link arm, with limited (160 degrees) joint excursion at shoulder and elbow. An agonist/antagonist pair is attached at each joint.

INPUT AND OUTPUT eODINGS

No cell is finely tuned to a specific input or output value to mimic the broad tunings or monotonic firing characteristics found in cortical visuomotor areas.

Arm position

By analogy with the role of muscle spindles, proprioceptive sensors are assumed to code muscle length. Arm position is thus represented by the population activity of NT = 20 neurons coding for the length of each agonist or antagonist. The activity of a sensor neuron k is defined by:

Tk = adLn(k))
where LIl ( k) is the length of muscle number n(k) , and ak a piecewise linear sigmoid:

L ~ Ak Ak < L < Ak L? Ak
Sensibility thresholds Ak are uniformly distributed in [Lmin , L max], and the dynamic range is Ak -Ak is taken constant, equal to Lmax -L min .

Desired direction

The direction V of the desired movement in visual space is coded by a population of N x = 50 neurons with cosine tuning in cartesian space. Each visual neuron j thus fires as:

Xj = V• Vj
Vj being the preferred direction of the cell. These 50 preferred directions are chosen uniformly distributed in 2-D space.

Motor Command

In attempt to model the existence of muscular synergies [START_REF] Lemon | The output map of the primate motor cortex[END_REF], we identified motor command with joint movement rather than with muscle contraction . A motor neuron i among Nt = 50 contributes to the effective movement M by its action on a synergy (direction in joint space) Mi. This collective effect is formally expressed by:

M .= LtiMi
where ti is the activity of motor neuron i. The 50 directions of action Mi are supposed uniformly distributed in joint space.

NETWORK STRUCTURE AND LEARNING

STRUCTURE OF THE NETWORK

Information concerning the position of the arm and the desired direction in cartesian space 0-

~, cY ~~~~----~ Cf t;
. . . . . Equation 1 is self-referent; so calculation is done in two steps. The feed-forward input first arrives at time zero when there is no activity in the layer; iterated action of the lateral connections comes into play when this feed-forward input vanishes.

The activity in the somatic layer is then combined with the visual directional information by the output sigma-pi neurons as follows: ti = L Xj Sij j

WEIGHTS AND LEARNING

The only adjustable weights are the Wijk linking the proprioceptive layer to the somatic layer. Connectivity is random and not complete: only 15% of the somatic neurons receive information on arm position. The visuomotor mapping is learnt by modifying the internal representation of the arm.

Motor commands issued by the network are correlated with the visual effect of the movement ("motor babbling"). More precisely, the learning algorithm is a repetition of the following cycle:

1. choice of an arm position among 5 positions (stars on Fig. 2) 2. random emission of a motor command (ti) 3. corresponding visual reafference (Xj) 4. weight modification according to a variant of the delta rule:

c:'Wijk oc (tiXj -Sij) Tk
The random commands are gaussian distributions of activity over the output layer. 5000 learning epochs are sufficient to obtain a stabilized performance. It must be noted that the error between the ideal response of the network and the actual performance never decreases completely to zero, as the constraints of the visuomotor transformation vary over the workspace.

RESULTS

NETWORK PERFORMANCE

Correct learning of the mapping was tested in 21 positions in the workspace in a pointing task toward 16 uniformly distributed directions in cartesian space. Movement directions generated by the network are shown in Fig. 2 (desired direction 0 degree is shown bold). Norm of movement vectors depends on the global activity in the network which varies with arm position and movement direction.

Performance of the network is maximal near the learning positions. However, a good generalization is obtained (directional error 0.3°, SD 12.1°); a bias toward the shoulder can be observed in extreme right or left positions. A similar effect was observed in psychophysical experiments (Ghilardi et a1. 1995). Preferred directions (PO) of output units were computed using a multilinear regression; a perfect cosine tuning was found, which is a consequence of the exact multiplication in sigma-pi neurons. Then, the population vector, the effective movement vector, and the desired movement were compared (Fig. 3) for two different arm configurations A and B marked on Fig. 2. The movement generated by the network (dashed arrow) is close to the desired one (dotted rays) for both arm configurations. However, the population vector (solid arrow) is not always aligned with the movement. The discrepancy between movement and population vector depends both on the direction and the position of the arm: it is maximal for positions near the borders of the workspace as position B. Fig. 3 (position B) shows that the deviations of the population vector are due to the anisotropic distribution of the PDs in cartesian space for given positions.

Difference between direction of action and preferred direction

Marked anisotropy in the distribution of PDs is compatible with accurate performance. To see why, let us call "direction of action" (DA) the motor cell's contribution to the movement. The distribution of DAs presents an anisotropy due to the geometry of the arm. This anisotropy is canceled by the distribution of PDs. Mathematically, if U is a N x 2 matrix of uniformly distributed 2D vectors, the PD matrix is UJ-1 whereas the DA matrix is UJ T , J being the jacobian of the angular-to-cartesian mapping. Difference between DA and PD has been plotted with concentric arcs for four representative neurons at 21 arm positions in Fig. 4. Sign and magnitude of the difference vary continuously over the workspace and often exceed 45 degrees. It can also be noted that preferred directions rotate with the arm as was experimentally noted by (Caminiti et a1. 1991).

DISCUSSION

We first asked how a network of broadly tuned neurons could produce visually guided arm movements. The model proposed here produces a correct behavior over the entire workspace. Biases were observed at the extreme right and left which closely resemble experimental data in humans (Ghilardi et a1. 1995). Single cells in the output layer behave as motor cortical cells do and the NPV of these cells correctly indicated the direction of movement for hand positions in the central region of the workspace (see [START_REF] Caminiti | Making arm movements within different parts of space: the premotor and motor cortical representation of a coordinate system for reaching to visual targets[END_REF]. Models of sensorimotor transformations have already been proposed. However they either considered motor synergies in cartesian coordinates (Burnod et a1. 1992), or used sharply tuned units [START_REF] Bullock | A self-organizing neural model of motor equivalent reaching and tool use by a multijoint arm[END_REF], or motor effects independent of arm position [START_REF] Salinas | Transfer of coded information from sensory to motor networks[END_REF]. Next, the use of the NPV to describe cortical activity was questioned. A fundamental assumption in the calculation of the NPV is that the PD of a neuron is the direction in which the arm would move if the neuron were stimulated. The model shows that the two directions DA and PD do not necessarily coincide, which is probably the case in motor cortex [START_REF] Scott | Changes in motor cortex activity during reaching movements with similar hand paths but different arm postures[END_REF]. It follows that the NPV often points neither in the actual movement direction nor in the desired movement direction (target direction), especially for unusual arm configurations. A maximum-likelihood estimator does not have these flaws; it would however accurately predict the desired movement out of the output unit activities, even for a wrong actual movement. In conclusion: (l) the NPV does not provide a faithful image of cortical visuomotor processes; (2) a correct NPV should be based on the DAs, which cannot easily be determined experimentally;

(3) planning of trajectories in space cannot be realized by the successive recruitment of motor neurons whose PDs sequentially describe the movement.

  Figure 1: Network Architecture

  Figure 2: Performance in a pointing task

Figure 3 :

 3 Figure 3: Actual movement and population vector in two arm positions

Figure 4 :

 4 Figure 4: Difference between direction of action and preferred direction for four units.


Dan Goreac
email: dan.goreac@u-pem.fr

Nonequivalence of Controllability Properties for Piecewise Linear Markov Switch Processes

In this paper we study the exact null-controllability property for a class of controlled PDMP of switch type with switchdependent, piecewise linear dynamics and multiplicative jumps. First, we show that exact null-controllability induces a controllability metric. This metric is linked to a class of backward stochastic Riccati equations. Using arguments similar to the euclidian-valued BSDE in [4], the equation is shown to be equivalent to an iterative family of deterministic Riccati equations that are solvable. Second, we give an example showing that, for switch-dependent coefficients, exact null-controllability is strictly stronger than approximate null-controllability. Finally, we show by convenient examples that no hierarchy holds between approximate (full) controllability and exact null-controllability. The paper is intended as a complement to [15] and [14].

Résumé

Nous étudions la propriété de zéro-contrôlabilité exacte pour une classe de processus de type Markovien à switch ayant des dynamiques linéaires par morceaux à coefficients switchés et bruit multiplicatif. Premièrement, nous montrons que la zérocontrôlabilité exacte induit une métrique de contrôlabilité. Celle-ci est liée à une classe d'équations stochastiques rétrogrades de type Riccati. En employant des arguments similaires aus EDSR classiques ([4]), l'équation de Riccati se réduit à une famille d'équations itératives déterministes de type Riccati qui admettent une solution unique. Deuxièmement, nous présentons un exemple montrant que, pour des systèmes à coefficients switchés, la propriété de zéro-contrôlabilité exacte est strictement plus forte que celle de zéro-contrôlabilité approchée. Finalement, nous montrons, à l'aide d'exemples particuliers, l'impossibilité à établir une hiérarchie entre les propriétés de contrôlabilité approchée (vers une cible arbitraire) et celle de zéro-contrôlabilité exacte.

Introduction

We study the exact null-controllability property for a class of piecewise deterministic Markov processes of switch type. More precisely, our model belongs to Markovian systems consisting of a couple mode/ trajectory (Γ, X) . The mode Γ is a pure jump uncontrolled Markov process corresponding to spikes inducing regime switching. The second component X obeys a controlled linear stochastic differential equation (SDE) with respect to the compensated random measure associated to Γ. The linear coefficients governing the dynamics depend on the current mode.

The exact null-controllability problem concerns criteria allowing one to drive the X T component to zero. This property is of particular importance in the study of regulatory networks (e.g. [START_REF] Cook | Modelling stochastic gene expression: Implications for haploinsufficiency[END_REF], [START_REF] Hasty | Noise-based switches and amplifiers for gene expression[END_REF], [START_REF] Krishna | Stochastic simulations of the origins and implications of long-tailed distributions in gene expression[END_REF], [START_REF] Crudu | Hybrid stochastic simplifications for multiscale gene networks[END_REF], etc.) to distinguish, for example, lytic pathways (e.g. [START_REF] Hasty | Noise-based switches and amplifiers for gene expression[END_REF]).

An extensive literature on controllability is available in different frameworks: finite-dimensional deterministic setting (Kalman's condition, Hautus test [START_REF] Hautus | Controllability and observability conditions of linear autonomous systems[END_REF]), infinite dimensional settings (via invariance criteria in [START_REF] Schmidt | Invariance theory for infinite dimensional linear control systems[END_REF], [START_REF] Curtain | Invariance concepts in infinite dimensions[END_REF], [START_REF] Russell | A general necessary condition for exact observability[END_REF], [START_REF] Jacob | Exact observability of diagonal systems with a finite-dimensional output operator[END_REF], [START_REF] Jacob | On controllability of diagonal systems with one-dimensional input space[END_REF], etc.), Brownian-driven control systems (exact terminal-controllability in [START_REF] Peng | Backward stochastic differential equation and exact controllability of stochastic control systems[END_REF], approximate controllability in [START_REF] Buckdahn | A characterization of approximately controllable linear stochastic differential equations[END_REF], [START_REF] Goreac | A Kalman-type condition for stochastic approximate controllability[END_REF], mean-field Brownian-driven systems in [START_REF] Goreac | Controllability properties of linear mean-field stochastic systems[END_REF], infinite-dimensional setting in [START_REF] Fernández-Cara | On the approximate controllability of a stochastic parabolic equation with a multiplicative noise[END_REF], [START_REF] Sirbu | Null controllability of an infinite dimensional SDE with state-and controldependent noise[END_REF], [START_REF] Barbu | Carleman estimates and controllability of linear stochastic heat equations[END_REF], [START_REF] Goreac | Approximate controllability for linear stochastic differential equations in infinite dimensions[END_REF], etc.), jump systems ( [START_REF] Goreac | A note on the controllability of jump diffusions with linear coefficients[END_REF], [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF], etc.). We refer to [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF] for more details on the literature as well as applications one can address using switch models.

The recent papers [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF] and [START_REF] Goreac | Approximate and approximate null-controllability of a class of piecewise linear markov switch systems[END_REF] consider different characterizations of approximate and approximate nullcontrollability properties for the same class of systems. However, they do not address the question of exact null-controllability (which, in non-stochastic time-homogeneous framework, is identified with approximate null-controllability).

The aim of the present paper is to offer a complement to the research topics in [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF] and [START_REF] Goreac | Approximate and approximate null-controllability of a class of piecewise linear markov switch systems[END_REF]. In [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF], a Riccati-type argument is used to characterize approximate controllability for systems with constant coefficients. Our first aim (in Section 3.1) is to give an answer to the problem left open in [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF]Remark 4] where a family of backward stochastic Riccati equations are presented and absence of results on the solvability is mentioned. As a by-product, the result provides a metric-type characterization of exact null-controllability in Section 3.2. In Section 3.3, we give an example of approximate null-controllable system which fails to be exactly null-controllable. Finally, in Section 3.4, we show by convenient examples that no hierarchy can be established between approximate (full) controllability and exact null-controllability. In other words, we present an example of approximate controllable yet non exactly null-controllable system and an example of null-controllable system which fails to be approximately controllable in certain directions.

We begin with presenting the model and the standing assumptions in Section 2.1. The technical constructions allowing to prove the theoretical results are gathered in Section 2.2. The controllability notions (exact null, approximate null, approximate) are given in Section 2.3. We gather in the same section some useful results on approximate and approximate null-controllability given in [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF] and [START_REF] Goreac | Approximate and approximate null-controllability of a class of piecewise linear markov switch systems[END_REF]. The main results on Riccati BSDEs are given in Section 3.1. The method allowing to deal with this stochastic system is based on the recent ideas in [START_REF] Confortola | Backward stochastic differential equations driven by a marked point process: an elementary approach, with an application to optimal control[END_REF]. As a by-product, the result on Riccati BSDE provides a metric-type characterization of exact null-controllability in Section 3.2. Hierachy (or absence of) between exact null-controllability and approximate null-controllability (resp. approximate full controllability) make the object of Section 3.3 (resp. Section 3.4).

Model and Preliminaries

The Model

We briefly recall the construction of a particular class of pure jump, non explosive processes on a space Ω and taking their values in a metric space (E, B (E)) . Here, B (E) denotes the Borel σ-field of E. The elements of the space E are referred to as modes. These elements can be found in [START_REF] Davis | Markov models and optimization[END_REF] in the particular case of piecewise deterministic Markov processes (see also [START_REF] Brémaud | Point processes and queues : martingale dynamics[END_REF]). To simplify the arguments, we assume that E is finite and we let p ≥ 1 be its cardinal. The process is completely described by a couple (λ, Q) , where λ : E -→ R + and the measure Q : E -→ P (E), where P (E) stands for the set of probability measures on (E, B (E)) such that Q (γ, {γ}) = 0. Given an initial mode γ 0 ∈ E, the first jump time satisfies P 0,γ0 (T 1 ≥ t) = exp (-tλ (γ 0 )) . The process Γ γ0 t := γ 0 , on t < T 1 . The post-jump location γ 1 has Q (γ 0 , •) as conditional distribution. Next, we select the inter-jump time

T 2 -T 1 such that P 0,γ0 T 2 -T 1 ≥ t / T 1 , γ 1 = exp -tλ γ 1 and set Γ γ0 t := γ 1 , if t ∈ [T 1 , T 2 ) . The post-jump location γ 2 satisfies P 0,γ0 γ 2 ∈ A / T 2 , T 1 , γ 1 = Q γ 1 ,
A , for all Borel set A ⊂ E. And so on. To simplify arguments on the equivalent ordinary differential system, following [4, Assumption (2.17)], we will assume that the system stops after a non-random, fixed number M > 0 of jumps i.e. P 0,γ0 (T M+1 = ∞) = 1.

We look at the process Γ γ0 under P 0,γ0 and denote by F 0 the filtration F [0,t] := σ {Γ γ0 r : r ∈ [0, t]} t≥0 . The predictable σ-algebra will be denoted by P 0 and the progressive σ-algebra by P rog 0 . As usual, we introduce the random measure q on Ω × (0, ∞) × E by setting q (ω, A)

= k≥1 1 T k (ω),Γ γ 0 T k (ω) (ω) ∈A , for all ω ∈ Ω, A ∈ B (0, ∞) × B (E) .
The compensated martingale measure is denoted by q. (Further details on the compensator are given in Section 2.2.)

We consider a switch system given by a process (X(t), Γ γ0 (t)) on the state space R N × E, for some N ≥ 1 and the family of modes E. The control state space is assumed to be some Euclidian space R d , d ≥ 1. The component X(t) follows a controlled differential system depending on the hidden variable γ. We will deal with the following model.

(1)

dX x,u s = [A (Γ γ0 s ) X x,u s + B (Γ γ0 s ) u s ] ds + E C Γ γ0 s-, θ X x,u s-q (ds, dθ) , s ≥ 0, X x,u 0 = x.
The operators A (γ) ∈ R N ×N , B (γ) ∈ R N ×d and C (γ, θ) ∈ R N ×N , for all γ, θ ∈ E. For linear operators, we denote by ker their kernel and by Im the image (or range) spaces. Moreover, the control process u : Ω × R + -→ R d is an R d -valued, F 0 -progressively measurable, locally square integrable process. The space of all such processes will be denoted by U ad and referred to as the family of admissible control processes. The explicit structure of such processes can be found in [21, Proposition 4.2.1], for instance. Since the control process does not (directly) intervene in the noise term, the solution of the above system can be explicitly computed with U ad processes instead of the (more usual) predictable processes.

Technical Preliminaries

Before giving the reduction of our backward Riccati stochastic equation to a system of ordinary Riccati differential equations, we need to introduce some notations making clear the stochastic structure of several concepts : final data, predictable and càdlàg adapted processes and compensator of the initial random measure. The notations in this subsection follow the ordinary differential approach from [START_REF] Confortola | Backward stochastic differential equations driven by a marked point process: an elementary approach, with an application to optimal control[END_REF]. Since we are only interested in what happens on [0, T ] , we introduce a cemetery state (∞, γ) which will incorporate all the information after T ∧ T M . It is clear that the conditional law of T n+1 given T n , Γ γ0

Tn is now composed by an exponential part on [T n ∧ T, T ] and an atom at ∞. Similarly, the conditional law of Γ γ0 Tn+1 given T n+1 , T n , Γ γ0

Tn is the Dirac mass at γ if T n+1 = ∞ and given by Q otherwise. Finally, under the assumption P 0,γ0 (T M+1 = ∞) = 1, after T M , the marked point process is concentrated at the cemetery state.

We set

E T : = ([0, T ] × E) ∪ {(∞, γ)}. For every n ≥ 1, we let E T,n ⊂ E T n+1
be the set of all marks of type e = ((t 0 , γ 0 ) , ..., (t n , γ n )) , where

(2) t 0 = 0, (t i ) 0≤i≤n are non-decreasing;

t i < t i+1 , if t i ≤ T ; (t i , γ i ) = (∞, γ) , if t i ≤ T, ∀0 ≤ i ≤ n -1,
and endow it with the family of all Borel sets B n . For these sequences, the maximal time is denoted by |e| := t n . Moreover, by abuse of notation, we set

γ |e| := γ n . Whenever T ≥ t > |e| , we set (3) e ⊕ (t, γ) := ((t 0 , γ 0 ) , ..., (t n , γ n ) , (t, γ)) ∈ E T,n+1 .
By defining (4)

e n := (0, γ 0 ) , T 1 , Γ γ0 T1 , ..., T n , Γ γ0 Tn ,
we get an E T,n -valued random variable, corresponding to our mode trajectories.

A càdlàg process Y continuous except, maybe, at switching times T n and taking its values in a topological vector space S is given by the existence of a family of B n ⊗ B ([0, T ]) /B (S)-measurable functions y n such that, for all e ∈ E T,n , y n (e, •) is continuous on [0, T ] and constant [0, T ∧ |e|] and

(5) If |e| = ∞, then y n (e, •) = 0. Otherwise, on T n (ω) ≤ t < T n+1 (ω) , y t (ω) = y n (e n (ω) , t) , t ≤ T . Similar, an S-valued F-predictable process Z defined on Ω × [0, T ] × E is given by the existence of a family of B n ⊗ B ([0, T ]) ⊗ B (E) /B (S) -measurable functions z n satisfying (6) If |e| = ∞, then z n (e, •, •) = 0. On T n (ω) < t ≤ T n+1 (ω) , z t (ω, γ) = z n (e n (ω) , t, γ) , for t ≤ T , γ ∈ E.
To deduce the form of the compensator, one simply writes

q (ω, dt, dγ) := n≥0 q n en(ω) (dt, dγ) 1 Tn(ω)<t≤Tn+1(ω)∧T such that (7) If n ≥ M, then q n e (dt, dγ) = δ γ (dγ) δ ∞ (dt) . If n ≤ M -1, q n e (dt, dγ) = λ(γ |e| )Q(γ |e| , dγ)1 |e|<∞,t∈[|e|,T ] Leb (dt) + δ γ (dγ) δ ∞ (dt) 1 (|e|<∞,t>T )∪|e|=∞ .
The coefficient function A (Γ γ0 t ) is adapted and can be seen as follows: if |e| = ∞, then A = 0; otherwise, one works with A γ |e| . Similar constructions hold true for C. In fact, the results of the present paper can be generalized to more general path-dependence of the coefficients.

Approximate Controllability, Exact and Approximate Null-Controllability

We will be dealing with the following notions of controllability.

Definition 1 (i) Given the finite time horizon T > 0, the system (1) is said to be approximately controllable (with initial mode γ 0 ∈ E) if, for every final data ξ ∈ L 2 Ω, F [0,T ] , P 0,γ0 ; R N (i.e. F [0,T ] -measurable, square integrable), every initial condition x ∈ R N and every ε > 0, there exists some admissible control process

u ∈ U ad such that E 0,γ0 |X x,u T -ξ| 2 ≤ ε.
(ii) The system (1) is said to be approximately null-controllable if the previous condition holds for ξ = 0.

(iii) The system (1) is said to be (exactly) null-controllable (with initial mode γ 0 ∈ E) if, for every initial condition x ∈ R N there exists some admissible control process u ∈ U ad such that X x,u T = 0, P 0,γ0 -a.s..

The approach of [14, Theorem 1] relies on the duality between the concepts of controllability and observability. For these reasons, one introduces the backward stochastic differential equation.

(8) dY T,ξ t = E Z T,ξ t (θ) q (dt, dθ) -A * (Γ γ0 t ) Y T,ξ t dt -E C * (Γ γ0 t , θ) Z T,ξ t (θ) q (dt, dθ) , Y T,ξ T = ξ ∈ L 2 Ω, F [0,T ] , P 0,γ0 ; R N .
The following characterization follows from standard considerations on the controllability linear operator(s) (cf. [14, Theorem 1]).

Theorem 2 ([14, Theorem 1]) The necessary and sufficient condition for approximate null-controllability (resp. approximate controllability) of ( 1) with initial mode γ 0 ∈ E is that any solution Y T,ξ t , Z T,ξ t (•) of the dual system [START_REF] Davis | Markov models and optimization[END_REF] for which Y T,ξ t ∈ ker B * (Γ γ0 t ) , P 0,γ0 ⊗Leb almost everywhere on Ω × [0, T ] should equally satisfy Y T,ξ 0 = 0, P 0,γ0 -almost surely (resp. Y T,ξ t = 0, P 0,γ0 ⊗Leb -a.s.).

Equivalent assertions are easily obtained by interpreting the system (8) as a controlled, forward one :

(9) dY y,v t = E v t (θ) q (dt, dθ) -A * (Γ γ0 t ) Y y,v t dt - E C * (Γ γ0 t , θ) v t (θ) q (dt, dθ) , Y y,v 0 = y ∈ R N .
The family of admissible control processes is given by v ∈ L 2 q; R N i.e. the space of all

P 0 ⊗ B (E) - measurable, R N -valued functions v s (ω, θ) on Ω × R + × E such that E 0,γ0 T 0 E |v s (θ)| 2 q (ds, dθ) < ∞, for all T < ∞.
Similar duality arguments yield the following characterization of (exact) null-controllability.

Proposition 3 The necessary and sufficient condition for exact null-controllability at time T > 0 of (1) with initial mode γ 0 ∈ E is the existence of a positive constant C T > 0 such that for every initial data y ∈ R N and every v ∈ L 2 q; R N , one has |y|

2 ≤ C T E 0,γ0 T 0 |B * (Γ γ0 t ) Y y,v t | 2 dt .
The proof is quasi-identical to the duality arguments in [14, Theorem 1] by invoking [24, Appendix B, Proposition B.1].

In the remaining of the section, unless stated otherwise, we assume the control matrix B to be modeindependent (constant). Using the explicit construction of BSDE with respect to marked-point processes, an invariance (algebraic) necessary and sufficient criterion for approximate null-controllability has been given in [START_REF] Goreac | Approximate and approximate null-controllability of a class of piecewise linear markov switch systems[END_REF]Theorem 6]. We recall the following invariance concepts (cf. [START_REF] Curtain | Invariance concepts in infinite dimensions[END_REF], [START_REF] Schmidt | Invariance theory for infinite dimensional linear control systems[END_REF]).

Definition 4 Given a linear operator

A ∈R N ×N and a family C = (C i ) 1≤i≤k ⊂ R N ×N , a set V ⊂ R N is said to be (A; C)-invariant if AV ⊂ V + k i=1 Im C i .
We construct a mode-indexed family of linear subspaces of R N denoted by V M,n γ 0≤n≤M, γ∈E by setting (10)

A * (γ) := A * (γ) - E (C * (γ, θ) + I) λ(γ)Q(γ, dθ) and V M,M γ = ker B * ,
for all γ ∈ E, and computing, for every 0

≤ n ≤ M -1, (11) V M,n γ the largest A * (γ) ; (C * (γ, θ) + I) Π V M,n+1 θ : θ ∈ E, Q (γ, θ) > 0 -invariant subspace of ker B * .
Here, Π V denotes the orthogonal projection operator onto the linear space V ⊂ R N . The explicit criterion is the following Theorem 5 ([15, Theorem 6]) The switch system ( 1) is approximately null-controllable (in time T > 0) with γ 0 as initial mode, if and only if the generated set V M,0 γ0 reduces to {0} .

In the same paper [START_REF] Goreac | Approximate and approximate null-controllability of a class of piecewise linear markov switch systems[END_REF], the property of approximate null-controllability for general systems is shown (using convenient examples) to be strictly weaker than approximate controllability. The following sufficient criterion is proven to guarantee the approximate controllability. )-invariant subspace of ker B * is reduced to {0}, for every γ ∈ E. Then, for every T > 0 and every γ 0 ∈ E, the system (1) is approximately controllable in time T > 0.

A Backward Stochastic Riccati Equation Approach to Exact

Null-Controllability

A Riccati Equation

A simple look at [14, Remark 4] shows that a key argument in the analysis of controllability properties resides in a family of backward stochastic Riccati equations. The authors of [START_REF] Goreac | Algebraic invariance conditions in the study of approximate (null-) controllability of Markov switch processes[END_REF]Remark 4] argue that their analysis is limited by solvability of the general BSDE of the form ( 12)

             dK ε,B t = K ε,B t A * (Γ γ0 t ) + A (Γ γ0 t ) K ε,B t -B (Γ γ0 t ) dt + E H ε,B t (θ) q (dt, dθ) + E f ε,B t (θ) * g ε,B t (θ) f ε,B t (θ) -H ε,B t (θ) q (dt, dθ) , where f ε,B t (θ) := C (Γ γ0 t , θ) K ε,B t -H ε,B t (θ) and g ε,B t (θ) := εI + K ε,B t + H ε,B t (θ) -1 , t ∈ [0, T ] . K ε,B t = 0, εI + K ε,B t + H ε,B t (θ) > 0, for almost all t ∈ [0, T ] .
Here, B (Γ γ0 t ) are positive semi-definite matrix. However, by using the structure of the jumps and inspired by [START_REF] Confortola | Backward stochastic differential equations driven by a marked point process: an elementary approach, with an application to optimal control[END_REF], existence of the solution of the previous BSDE will be reduced to a family of itterated (classical) Riccati equations.

The first result gives existence and uniqueness for the solution of the previous equation. Before stating and proving this result, let us concentrate on the specific form of the jump contribution H. We consider a càdlàg process K ε,B continuous except, maybe, at switching times T n . Then, as explained before, this can be identified with a family k n,ε,B . We construct, for every n ≥ 0, [START_REF] Goreac | Controllability properties of linear mean-field stochastic systems[END_REF] k n+1,ε,B (e, t, γ) := k n+1,ε,B (e ⊕ (t, γ) , t) 1 |e|<t and K ε,B t can be obtained by simple integration of the previous quantity with respect to the conditional law of T n+1 , Γ γ0 Tn+1 knowing F Tn . Then, H is simply given by h n,ε,B (e, t, γ) := k n+1,ε,B (e, t, γ) -k n,ε,B (e, t) . The main theoretical contribution of the subsection is the following.

Theorem 7

We assume that P 0,γ0 (T M+1 = ∞) = 1, for some M ≥ 1. For every ε > 0 and every T > 0, the Riccati BSDE (12) admits a unique solution K ε,B

• , H ε,B

• (•) consisting of an S N + -valued (i.e. positive semi-definite) càdlàg process K ε,B • continuous everywhere except, maybe, at jump times and an S N -valued (i.e. symmetric matrix-valued) F-predictable process H ε .

Proof. For notation purposes, we will consider B to be fixed and drop the dependency on B. The proof consists of two steps.

Step 1. Using the previous structure of the candidate to the solution of the Riccati equation, one gets an equivalent system of ordinary Riccati-type differential equations

                 k M,ε (•) = 0, k n,ε (e n , T ) = 0, for every 0 ≤ n ≤ M -1, dk n,ε (e n , t) = k n,ε (e n , t) A * γ |en| + A γ |en| k n,ε (e n , t) -B γ |en| dt
+ E (f n (e n , t, θ)) * g n (e n , t, θ) f n (e n , t, θ) q (dt, dθ) -E k n+1,ε (e n , t, θ) -k n,ε (e n , t) q (dt, dθ) , for t ∈ [0, T ] , where f n (e n , t, θ) = C γ |en| , θ + I k n,ε (e n , t) -k n+1,ε (e n , t, θ) and g n (e n , t, θ) := εI + k n+1,ε (e n ⊕ (t, γ) , t)

-1 1 |e|<t ,
Under the condition that k n,ε (e n , t) ≥ 0, for almost all t ∈ [0, T ] .

The fact that the two systems are indeed equivalent follow from the same arguments as those in [4, Theorem 2]. Step 2. Thus, solvability of the Riccati backward stochastic equation reduces to the solvability of the previous system or, again, to the solvability (in S N + ) of the following equation

• p (t) = p (t) a + a * p (t) -Π + p (t) E b (θ) r -1 t (θ) b * (θ) ν (dθ) p (t) , for t ∈ [t 0 , T ] , p (T ) = εI.
by setting, for a fixed e n (and t > t 0 := |e n |)

a := A * γ |en| -λ γ |en| 1 2 I + E C * γ |en| , θ Q γ |en| , dθ -ε E C * γ |en| , θ + I εI + k n+1,ε (e n ⊕ (t, θ) , t) -1 Q γ |en| , dθ ,
Π := B γ |en| + ε E εI + k n+1,ε (e n ⊕ (t, θ) , t) -1 k n+1,ε (e n ⊕ (t, θ) , t) Q γ |en| , dθ b (θ) := C * γ |en| , θ + I, r t (θ) := εI + k n+1,ε (e n ⊕ (t, θ) , t) and ν (dθ) = λ γ |en| Q γ |en| , dθ
Existence and uniqueness for this equation is standard. Indeed, one notes that Π ≥ 0 and r ≫ 0 (provided that k n+1,ε ≥ 0). If E reduces to a singletone, then this is the classical equation for deterministic control problems (see [START_REF] Yong | Stochastic Controls. Hamiltonian Systems and HJB Equations[END_REF]Chapter 6, Equation 2.34]). The existence and uniqueness is guaranteed by [28, Chapter 6, Corollary 2.10]. For the general case, one assumes that E is given by the standard basis of R p and works with b = ν (e 1 )b e 1 , ..., ν (e p )b (e p ) and

r t :=     r t e 1 0 ... 0 0 r t e 2 ... 0 ... ... ... ... 0 0 ... r t (e p )     ≥ εI ≫ 0.
The proof is complete by descending recurrence over n ≤ M.

First Application: Null-Controllability Metric(s)

Proposition 8 A necessary and sufficient condition for exact null-controllability of (1) with initial mode γ 0 ∈ E at time T > 0 is that the pseudonorm

(14) R N ∋ y -→ p (y)
, where p 2 (y) := inf

v∈L 2 (q;R N ) E 0,γ0 T 0 Π (ker B * (Γ γ 0 t )) ⊥ (Y y,v t ) 2 dt be a norm on R N .
Proof. It is clear that the application p has non-negative values. Homogeneity is a consequence of the equality Y ay,av

• = aY y,v
• , for all y ∈ R N , all a ∈ R and all v ∈ L 2 q; R N (due to the linearity of ( 9)).

To prove the subadditivity, one simply notes

Y y1+y2,v 1 +v 2 • = Y y1,v 1 • + Y y2,v 2 • , for all y 1 , y 2 ∈ R N and all v 1 , v 2 ∈ L 2 q; R N . It follows that p (y 1 + y 2 ) ≤ E 0,γ0 T 0 Π (ker B * (Γ γ 0 t )) ⊥ Y y1+y2,v 1 +v 2 t 2 dt 1 2 ≤ E 0,γ0 T 0 Π (ker B * (Γ γ 0 t )) ⊥ Y y1,v 1 t 2 dt 1 2 + E 0,γ0 T 0 Π (ker B * (Γ γ 0 t )) ⊥ Y y2,v 2 t 2 dt 1 2
, for all v 1 , v 2 ∈ L 2 q; R N . The conclusion follows by taking infimum over such control processes. It follows that p is a pseudonorm (independently of the fact that the system is approximately null-controllable). Necessity follows from Proposition 3 and sufficiency from the equivalence of norms on R N by applying Proposition 3.

Using the form of the Riccati BSDE [START_REF] Goreac | A note on the controllability of jump diffusions with linear coefficients[END_REF], one infers the following explicit condition.

Corollary 9 A necessary and sufficient condition for exact null-controllability of (1) with initial mode γ 0 ∈ E at time T > 0 is that the positive-semidefinite matrix k 0 := inf ε>0 K ε 0 , where, for every ε > 0, K ε is the unique solution of the Riccati equation ( 12) for B := BB * be positive definite. In this case, the metric p given in ( 14) is induced by k 0 i.e. p (y) = k 0 y, y , for all y ∈ R N .

Proof. This result is quite classical (see, e.g. [START_REF] Sirbu | Null controllability of an infinite dimensional SDE with state-and controldependent noise[END_REF] for the Brownian-noise case). For our readers' sake, we sketch the proof. Let us fix, for the time being ε > 0. Then, according to Theorem 7, the Riccati equation ( 12) admits a unique solution. A simple application of Itô's formula (cf. [18, Chapter II, Section 5, Theorem

5.1]) to K ε t Y y,v t , Y y,v t on [0, T ] yields E 0,γ0 T 0 Π [Ker(B * (Γ γ 0 t ))] ⊥ Y y,v t 2 dt = K ε 0 y, y -εE 0,γ0 T 0 |u t | 2 dt + E 0,γ0   T 0 (εI + K ε t + H ε t (θ)) -1 2 f t (θ) Y y,v t- -(εI + K ε t + H ε t (θ)) 1 2 v t (θ) 2 dt   .
One easily notes that inf

v∈L 2 (q;R N ) E 0,γ0 T 0 Π [Ker(B * (Γ γ 0 t ))] ⊥ Y y,v t 2 dt = lim inf ε→0
K ε 0 y, y and the conclusion follows.

Non-equivalence Between Exact and Approximate Null-Controllability

The following example presents a switching system which is approximately null-controllable without being exactly null-controllable.

Example 10 We consider a two-dimensional state space and a one-dimensional control space. Moreover, we consider the mode to switch randomly between three states (for simplicity, E = e 1 , e 2 , e 3 is taken to be the standard basis of R 3 ). The transition measure is given by

Q :=   0 1 0 0 0 1 0 1 0   .
The coefficients are given by

A e 1 = A := 0 0 1 0 , A (γ) := 0 2×2 , if γ = e 1 , B := 1 0 , C (γ) := 0 2×2 .

a) Approximate null-controllability

With the definition [START_REF] Goreac | Approximate controllability for linear stochastic differential equations in infinite dimensions[END_REF], one easily establishes

V M,n γ = ker B * , if γ = e 1 , {0} , if γ = e 1 .
, for every M ≥ 1 and every 0 ≤ n ≤ M. It follows that the system is null-controllable if and only if γ 0 = e 1 . b) Limit of the Riccati equations in the approximate null-controllable case (initial mode γ 0 = e 1 ). Quid est for the controllability metric ? In this case, we recall that the limit of the solutions of the Riccati equations is given by

p 2 (y) := inf v∈L 2 (q;R N ) E 0,γ0 T 0 |B * Y y,v t | 2 dt , where dY y,v t = E v t (θ) q (dt, dθ)-A * (Γ γ0 t ) Y y,v t dt, Y y,v 0 = y ∈ R 2 .
Starting from y := 0 1 , with the feedback control process v ε t := -

1 0 0 0 Y t + 0 e 2ε -1 -e ε +1+ε 1 t∈[ε,2ε] 1 t≤T1 , one gets Y y,v t , 1 0 = 1 -e t 1 t≤ε + 1 -e t + t -ε + 1 -e t-ε e 2ε -1 -e ε + 1 + ε 1 ε≤t≤2ε 1 t≤T1 .
One easily notes that 0 ≥ Y y,v t , 1 0 ≥ c ε := -2 e 2ε -1 . Then, by taking infimum over ε > 0, it follows that p 2 0 1 = 0 and it cannot induce a norm. As consequence, by invoking Corollary 9, the system fails to be (exactly) null-controllable.

Remark 11 Of course, a direct proof of null-controllability can also be given based on the eigenvector 0 1 .

Absence of null-controllability is obvious for γ 0 = e 1 (the system is not even approximately null-controllable).

In the case γ 0 = e 1 , we reason by contradiction. Let us assume that, for some admissible control process u, the system is exactly controllable at time T starting from x 0 = 0 1 . Then, prior to the first jump time,

X x0,u t = 1 0 t 1 0 1 + t 0 1 0 t -s 1 u s 0 ds = t 0 u s ds 1 + t 0 (t -s) u s ds .
Since, on [0, T 1 ) , u is deterministic and square integrable, there exists T > t 0 > 0 such that 1+ t 0 (t -s) u s ds > 1 2 , for every t ≤ t 0 (consequence of the absolute continuity). Hence, on T 1 ≤ t 0 , one gets X x0,u T1-,

0 1 > 1 2 . Second, one notes that X x0,u T1-= X x0,u T1 and d X t , 0 1 = 0 for t ≥ T 1 to conclude that X x0,u T , 0 1 > 1 2 , on T 1 ≤ t 0 . Since P (T 1 ≤ t 0 ) = 1 -e -t0
> 0, it follows that u cannot lead to 0 with full probability. Therefore, although it is approximately null-controllable for some initial modes, the switch system is never exactly null-controllable.

Exact Null-Controllability vs. Approximate (Full) Controllability

It has been shown in [START_REF] Goreac | Approximate and approximate null-controllability of a class of piecewise linear markov switch systems[END_REF] that, in general, approximate controllability is strictly stronger than approximate null-controllability. In the light of the previous example, it is then natural to ask oneself whether the condition on p (given by ( 14)) being a metric implies approximate controllability of the initial system. The answer is negative. We begin with an example of a system governed by an off/on mode which is approximately null-controllable iff the initial mode is set off and is never approximately controllable. We show that, for this system, the Riccati equations give a controllability metric (iff the initial mode is set off).

Example 12

We consider a two-dimensional state space and a one-dimensional control space. Moreover, we consider the mode to switch randomly between inactive 0 and active 1 (i.e. E = {0, 1}). The coefficients are given by

A (γ) = A := 0 0 1 0 , B := 1 0 , C (γ, 1 -γ) := -1 0 γ -1 .
a) Approximate null-controllability One easily establishes V M,n γ = span {γe 2 } , for every M ≥ 1 and every 0 ≤ n ≤ M. It follows that the system is null-controllable if and only if γ 0 = 0. b) Approximate controllability (initial mode γ 0 = 0)

One easily checks that Y

t := 0 Γ 0 t , Z t (•) = 0 (-1) Γ 0 t- , for 0 ≤ t ≤ T satisfies the BSDE (8) with final data Y T = 0 Γ 0 T
. Since this solution stays in ker B * and it is not trivially zero, it follows that the system is (never) approximately controllable. c) Riccati equations in the approximate null-controllable case (initial mode γ 0 = 0) One easily notes that the Riccati equations lead to

dk 0,ε (0, t) = k 0,ε (0, t) 0 1 0 0 + 0 0 1 0 k 0,ε (0, t) - 1 0 0 0 + k 0,ε (0, t) dt + k 1,ε (0, 0, t, 1) εI + k 1,ε (0, 0, t, 1) -1 k 1,ε (0, 0, t, 1) -k 1,ε (0, 0, t, 1) dt, t ∈ [0, T ] .
Since k n+1,ε (0, 0, t, 1) ≥ 0, this solution is at least equal to the one given by K

t := a (t) b (t) b (t) c (t)
, where a (t) := 1 -e t-T , b (t) := (T + 1 -t) e t-T -1 and c (t) := 2 -1 + (T + 1 -t) 2 e t-T . Hence, K 0 is positive definite (for every T > 0) and so is lim inf ε→0+ K 0,ε 0 (to prove this, one simply studies the sign of the function

T → 1 -e -T 2 -2 + T 2 + 2T e -T -(T + 1) e -T -1 2 on R * + ).

Remark 13

The reader may want to note that in the case presented in the previous example, the system is approximately null-controllable if and only if it is (exactly) null-controllable. To give an explicit control leading from x 0 to 0 for the initial setting γ 0 = 0, one proceeds as follows. Prior to the first jump, the system is a deterministic one and given by

dΦ x0,u 1 t = (A -C (0, 1)) Φ x0,u 1 t + Bu 1 t dt, Φ x0,u 1 0 = x 0 ∈ R 2 .
Since Kalman's condition is satisfied for this dterministic system, it is exactly null-controllable at time T > 0.

An explicit control is obtained by considering the (deterministic) controllability Gramian G T , respectively the induced (open-loop) control process u 1 given by G T := T 0 e (A-C(0,1))(T -s) BB * e (A-C(0,1)) * (T -s) ds, u 1 (t) := -B * e (A-C(0,1)) * (T -t) G -1 T e (A-C(0,1)) * T x 0 1 t≤T .

We obtain a stochastic control by setting u p (t) = 0, i.e. we take null-control after the first jumping time. Then, it is obvious that, on T 1 ≥ T, X x0,u T = Φ x0,u 1 T = 0. On T 1 < T, X x0,u T1 = 0 and the conclusion follows.

Finally, one is entitled to ask whether approximate null-controllability implies exact null-controllability. The answer is, again, negative proving that approximate controllability and exact null-controllability are, in general, completely different properties. To illustrate this, let us take, once again, a glance at the first example.

Example 14 We consider a two-dimensional state space and a one-dimensional control space. Moreover, we consider the mode to switch randomly between three states (for simplicity, E = {e 1 , e 2 , e 3 } is taken to be the standard basis of R 3 ). The transition measure is given by Q := As we have seen before (in Example 10), this system is never exactly null-controllable. We assume the system to only jump once. We consider a solution of

1 0 = Φ v T1 , 1 0 -Φ v 1 T1 , 0 1 = -Φ v 1 T1 , 0 1 .
One deduces that Φ v = 0, P 0,e1 -almost surely. Since the process is no longer allowed to jump after T 1 , it follows that the equality actually holds on [0, T ] and the initial system (1) is approximately controllable (cf. Theorem 2).

Proposition 6 (

 6 [START_REF] Goreac | Approximate and approximate null-controllability of a class of piecewise linear markov switch systems[END_REF] Condition 10]) Let us assume that the largest(A * (γ) ; [(C * (γ, θ) + I) Π ker B * : Q (γ, θ) > 0]

  γ) := 0 2×2 , if γ = e 1 , B := 1 0 , C (γ) := 0 2×2

1 ( 1 t= Φ v 1 t , 0 1 .

 111 θ) q (dt, dθ) -A * (Γ γ0 t ) Y y,v t dt, Y y,v 0 = y ∈ R Nthat belongs to ker B * , P 0,e1 -a.s. Due to the approximate null-controllability, it follows that y = 0. Prior to the first jump, v is given by a deterministic functionv 1 = v 1,1 v 1,2and Y y,v t coincides with the deterministicsolution of dΦ v 1 t = -v 1 t -A * Φ v 1 t dt, Φ v 1 0 = 0. Since Φ v 1 t ∈ ker B * = span 0 for all t ∈ [0, T ]) one has, for Leb-almost all t ∈ [0, T ] , -v 1,Hence, at the first jumping time, one has0 = Y y,v T1 ,

and v 1, 1 t= 0 ,

 10 Leb-almost surely on [0, T ] . Then the derivative of Leb-almost surely on [0, T ]. As a consequence, for all (due to right continuity) t ∈ [0, T 1 ] , Y y,v t

This work has been partially supported by he French National Research Agency project PIECE, number ANR-12-JS01-0006.


Rémi Imbach

A Subdivision Solver for Systems of Large Dense Polynomials

Keywords: Interval Arithmetic, Subdivision, Adaptive Multi-Precision, Real Solutions, Large Dense Polynomials Arithmétique par intervalles, subdivision, précision multiple adaptative, solutions réelles, polynômes denses de haut degré et grande bit-size

We describe here the package subdivision solver for the mathematical software SageMath. It provides a solver on real numbers for square systems of large dense polynomials. By large polynomials we mean multivariate polynomials with large degrees, which coefficients have large bit-size. While staying robust, symbolic approaches to solve systems of polynomials see their performances dramatically affected by high degree and bit-size of input polynomials. Available numeric approaches suffer from the cost of the evaluation of large polynomials and their derivatives. Our solver is based on interval analysis and bisections of an initial compact domain of R n where solutions are sought. Evaluations on intervals with Horner scheme is performed by the package fast polynomial for SageMath. The non-existence of a solution within a box is certified by an evaluation scheme that uses a Taylor expansion at order 2, and existence and uniqueness of a solution within a box is certified with krawczyk operator. The precision of the working arithmetic is adapted on the fly during the subdivision process and we present a new heuristic criterion to decide if the arithmetic precision has to be increased.

Introduction

subdivision solver is a solver for square systems of polynomial equations using exhaustive search in an initial bounded real domain given as a box (i.e. a vector of intervals). It is specifically designed to handle systems of large dense polynomials and uses adaptive multiprecision arithmetic to stay robust to hard cases. subdivision solver is proposed as a package for the mathematical software SageMath1 . It is based on the package fast polynomial (see [START_REF] Moroz | Fast polynomial evaluation and composition[END_REF]) for SageMath that provides fast evaluation of polynomials on intervals and uses Horner forms.

The strategies that underlies the design of subdivision solver are:

• using as much as possible double precision interval arithmetic: in this case interval computations are supported by the boost2 interval library that is efficient. When the arithmetic precision is increased, interval computations are supported by the Multi-Precision Floatingpoint Interval library mpfi (see [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the mpfi library[END_REF]) that allows arbitrary arithmetic precision.

• using partial derivatives at second order to obtain sharp interval evaluations; partial derivatives are symbolically computed at the initialization of the solving.

• minimizing evaluations of polynomials and their derivatives that are costly.

The rest of this report is organized as follows. Sec. 2 gives some basics on interval analysis and describes the main tools the solver uses. Sec. 3 presents three criteria that are used to decide to increase arithmetic precision. The two first come from the state of the art, the third is new. Sec. [START_REF] Ramon E Moore | Introduction to Interval Analysis[END_REF] describes the algorithms that we implemented. Sec. 5 proposes some numerical results. The rest of the introduction overviews installation and usage of the package.

A short user manual subdivision solver is a package for SageMath, and works with versions newer than 7.0. Its sole dependency is the package fast polynomial 3 for SageMath that provides fast evaluation of polynomials on intervals and uses Horner forms. After it has been download, move to the directory where it is and install it with sage fast polynomial-0.9.4.spkg. Then download subdivision solver 4 and install it with the command sage subdivision solver-0.0.1.spkg. Thats it. We give here an example of use of subdivision solver on randomly generated polynomials. Above instructions can be either copy-pasted in a sage console to be executed, or copy pasted in a file myfile.sage and executed with the command sage myfile.sage.

The initializer of the class subdivision solver takes as arguments a list of polynomial and the list of variables of the polynomial ring.

The method solve takes as first argument the domain where the real solutions are sought, written as a column vector of RealIntervalField elements. Notice that the arithmetic precision of the first element of the domain determines the initial precision. The arithmetic precision is at least 53. When it is 53, bounds and intervals are double precision floating points, and operations on intervals are directly transposed into operations on double with boost.

To solve the system given as input, the initial domain will be subdivided and/or contracted into sub-domains; the second argument of the method solve is the minimum width of sub-domains to be explored during the resolution. It can be zero, and in this case the process is not guaranteed to terminate, in particular if the system admits a root with multiplicity.

When solving the system given in input, it can arise that arithmetic precision is not sufficient when, for instance, subdividing a box does not help to decide where is a solution due to accumulation of errors in interval computations. Such cases are detected, and the arithmetic precision is increased to face it. The third argument of solve is an upper bound of the arithmetic precision used to compute with intervals. Precision can be arbitrary large. mpfi is used when the precision is larger than 53.

The last argument of the method solve is a string determining whether statistics of the process are printed or not. Its value can be in { 'silent', 'stats', 'stats color' }.

The output of solve is an integer with the following meaning:

0 Each solution lying in the interior of the initial domain has bean isolated in a box; isolating boxes of solutions are obtained with getSolutions.

1 Only a subset of solutions have been isolated; boxes that may contain solutions are obtained with getUndetermined but exploring these boxes will require a larger arithmetic precision than initially allowed.

2 Only a subset of solutions have been isolated; boxes of width smaller than allowed and/or requiring larger arithmetic precision where encountered.

Finally, giving a real number r as argument to the method getSolutions makes width of boxes isolating solutions that are returned smaller than r if maximum precision allows it. Otherwise isolating boxes will have width as small as possible.

Inria 2 Some tools of Interval Analysis

Notations

Let IR be the set of non-empty closed intervals of R. We will denote an element of IR with a lowercase bold letter as x. If x ∈ IR, we will denote by l(x) (resp. u(x)) its lower (resp. upper) bound, by w(x) its width defined as u(x)l(x) and by m(x) its center (that is l(x) + w(x)

2 ). The interior ]l(x), u(x)[ of x will be denoted by int(x) and its boundary x \ int(x) by ∂(x).

Let m ∈ N + * and consider IR m . We call box an element of IR m and we will denote it by an uppercase bold letter as X. Let X = (x 1 , . . . , x m ) be an element of IR m , we will denote by m(X) its center (m(x 1 ), . . . , m(x m )) and by w(X) its width defined as w(X) = max 1≤i≤m w(x i ). We note (int(x 1 ), . . . , int(x m )) the interior of X and ∂(X) its boundary X \ int(X).

Usual arithmetic operators can be extended to intervals and boxes.

Evaluating Functions on Boxes

Consider the polynomial function f : R m → R and let X ∈ IR m . We will denote by f (X) the set {f (X)|X ∈ X}. The following definitions are picked up from [START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF].

We will call interval extension of f a function f :

IR m → IR s.t. X ∈ IR m ⇒ f (X) ⊆ f (X) (1)
and

X ∈ R m ⇒ f (X) = f (X). (2)
We say that an interval extension f is inclusion monotonic if

X ⊆ Y ⇒ f (X) ⊆ f (Y). (3)
For a polynomial f , a natural interval extension of f that is inclusion monotonic is obtained by combining interval operators that intervene in the Horner form of f (see [START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF]). We will denote it 0 f , and 0 f (X) will be the value of the evaluation of 0 f on X ∈ IR m . An other inclusion monotonic interval extension (see [START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF]) that is classically used is the mean value form 1 f defined as

1 f (X) = f (P) + J f (X)(X -P) ( 4
)
where P is the center of X represented as an interval and J f an inclusion monotonic interval extension of the Jacobian matrix J f of f (i.e. the m components of J f are intervals extensions of the m partial derivatives of f ). 0 J f , which components are

0 ∂f x1 , . . . , 0 ∂f xm , is classically used for J f .
Here we will use the interval extension obtained by considering Taylor expansion at order 2 around the center P of X, defined as:

2 f (X) = f (P) + J f (P)(X -P) + 1 2 (X -P) t ( 0 H f (X))(X -P) (5)
where 0 H f is the natural interval extensions of the Hessian matrix H f of f . 2 f is shown in [START_REF] Louis B Rall | Mean value and taylor forms in interval analysis[END_REF] to be an inclusion monotonic interval extension of f , and we will call it evaluation at order two of f . Above definitions are easily extended to functions F : R m → R m . Notice that if F is an interval extension of F , then we have 0 / ∈ F (X) ⇒ F does not vanish on X.

RT n°476

Algorithm 1 CertifyNoSolution(F ,X) Input: A function F : R m → R m and a box X. Output: true or false; if true then F = 0 has no solution in X.

1: if (0 / ∈ 0 F (X) ∩ 2 F (X)) or ( 2 K F (X) ∩ X = ∅) then 2:
return true 3: return false Algorithm 2 CertifyOneSolution(F ,X) Input: A function F : R m → R m and a box X. Output: true or false; if true then X contains a unique solution of F = 0.

1: if 2 K F (X) ⊂ int(X) then 2:
return true 3: return false

The Krawczyk Operator

We consider now a function F : R m → R m defined as F (X) = (f 1 (X), . . . , f m (X)) where f i : R m → R for 1 ≤ i ≤ m and we recall the definition of the Krawczyk operator K F that is a classical tool in interval analysis (see for instance [START_REF] Kearfott | Rigorous global search: continuous problems. Nonconvex optimization and its applications[END_REF]Def. 1.16] or [4, Theo. 8.2] or [9, Sec. 7] ).

K F (X) = P -(J F (P)) -1 F (P) + (I -(J F (P)) -1 J F (X))(X -P) (7)
where P is the center of X represented as an interval and J F is the interval extensions of the Jacobian matrix J F of F . K f has the following properties:

(K 1 ) K F (X) ⊂ int(X) ⇒ F (X) = 0 has one and only one solution in X that is in K F (X) (K 2 ) K F (X) ∩ X = ∅ ⇒ F (X) = 0 has no solution in X (K 3 ) if X and Y are s.t. K F (X) ⊂ int(X) and K F (Y) ⊂ int(Y) then X ∩ Y = ∅ ⇒ X and Y contain the same solution of F (X) = 0. Proof of (K 3 ): because K F is inclusion monotonic (i.e. X ⊆ Y ⇒ K F (X) ⊆ K F (Y)), see [3], one has K F (X ∩ Y) ⊂ int(X ∩ Y) and X ∩ Y contains a unique solution.
When using a mean value extension of J F around the center P of X to compute J F (X) one can rewrite Eq. 7 as

2 K F (X) = P -(J F (P)) -1 (F (P) + H) (8)
where H is a box in IR m which i-th component h i is:

h i = (X -P) t ( 0 H fi (X))(X -P) (9)
where 0 H fi is the natural interval extension of the Hessian matrix H fi of f i . We will call 2 K F (X) the Krawczyk operator at order 2 since it uses second order partial derivatives of F ; we assume here that it satisfies properties (K 1 ), (K 2 ) and (K 3 ). One can use (K 2 ) and Eq. 6 to justify Algo. 1 that certifies that a system F = 0 with F : R m → R m has no solutions on a box X. (K 1 ) justifies Algo. 2 that certifies that F = 0 has a unique solution on X. Algo. 3 decides if the solution of

F = 0 contained in a box X s.t. K F (X) ⊂ int(X) is already contained in a box of a list L containing boxes Y s.t. K F (Y) ⊂ int(Y). It is correct from (K 3 ). Inria Algorithm 3 IsSolInList(F ,X,L) Input: A function F : R m → R m , a box X and a list of boxes L s.t. K F (X) ⊂ int(X) and ∀Y ∈ L, K F (Y) ⊂ int(Y). Output: true or false; if true then the solution of F = 0 contained in X is in a box of L. 1: for Y ∈ L do 2: if X ∩ Y = ∅ then 3:
return true 4: return false

Adapting Arithmetic Precision

Intervals are usually represented by their two bounds represented by floating points, and interval operators by operations on floating points with appropriated rounding policy. In what follows, we will call precision of the arithmetic the number of bits of the mantissa of floating points.

The branch and bound method we did implement uses recursive bisections of boxes until it is possible to certify either the absence or the existence and uniqueness of a solution in each box. Lack of accuracy can intervene during this process either when a box has a component that has roughly the smallest width allowed by the arithmetic precision (i.e. the bounds are consecutive floating points), or when accumulated rounding affects the inclusion monotonicity of an interval extension. These cases have to be detected to avoid infinite computations. When they occur a simple strategy is to double the arithmetic precision and to continue the computations.

We recall in 3.1 two criteria pickud up from [START_REF] Revol | Interval newton iteration in multiple precision for the univariate case[END_REF] that decide if the arithmetic precision is not sufficient, and we present in 3.2 a criterion that is, as far as we know, new and is used as an heuristic in our solver.

Criteria of [7]

[7] presents a simple branch and bound algorithm using adaptive multi-precision to find zeros of an univariate real function f with interval extension f , using interval newton method. At each step, an interval x is bisected in x 1 , x 2 s.t. x = x 1 ∪ x 2 , and the arithmetic precision is increased if one of the following conditions holds:

• w(x 1 ) ≥ w(x) or w(x 2 ) ≥ w(x), • w( f (x 1 )) ≥ w( f (x)) or w( f (x 2 )) ≥ w( f (x)).
The first criterion is satisfied when x has no point of its interior that can be represented as a floating point with actual precision. When the second criterion is satisfied, Eq. 3 is not true and f is not an interval extension of f due to numeric inaccuracy. One can alternatively test if f (x) ⊆ ( f (x 1 ) ∪ f (x 2 )). We adapt here these criteria to a multi-variate context with F : R m → R m , F : IR m → IR m an interval extension of F , X 1 , X 2 and X elements of IR m such that X = X 1 ∪ X 2 . In the branch and bound method presented in Sec. 4, the arithmetic precision is increased if one of the following conditions holds:

(C 1 ) w(X 1 ) ≥ w(X) or w(X 2 ) ≥ w(X), (C 2 ) F (X) ⊆ ( F (X 1 ) ∪ F (X 2 )).
Algo. 4 verifies these conditions and returns true if (C 1 ) or (C 2 ) is satisfied.

RT n°476

Algorithm 4 CheckCond12(F or F ,X,X 1 ,X 2 ) Input: A function F : R m → R m or a list F of functions f : R m → R, three boxes X, X 1 , X 2 . Output: if (C 1 ) or (C 2 ) is satisfied then true; false otherwise. Algorithm 5 CheckCond3(F ,X) Input: A function F : R m → R m , a box X. Output: if (C 3 ) is satisfied then true; false otherwise.

A new criterion

Consider Eqs. ( 7) and [START_REF] Revol | Motivations for an arbitrary precision interval arithmetic and the mpfi library[END_REF]. Whatever its order (we state it here for order two), one can rewrite Krawczyk operator as 2 K F (X) = -( 0 J F (P)) -1 ( 0 F (P)) + P + something [START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF] to highlight that F (P) and J F (P) are interval evaluations. Hence one has

w( 2 K F (X)) = w(( 0 J F (P)) -1 ( 0 F (P))) + w(something) (11)
and a necessary condition for 2 K F (X) to certify the existence of a unique solution in X is w(( 0 J F (P)) -1 ( 0 F (P))) < w(X): otherwise one can not have 2 K F (X) ⊂ int(X). Recall that P is a point (i.e. w(P) = 0) hence one should have, w( 0 F (P)) = 0 and w(( 0 J F (P)) -1 ( 0 F (P))) = 0.

In the branch and bound method presented in Sec. 4, the arithmetic precision is increased if the following condition holds: (C 3 ) w(( 0 J F (P)) -1 ( 0 F (P))) ≥ w(X) and (P -( 0 J F (P)) -1 ( 0 F (P)))) ∩ X = ∅.

The condition w(( 0 J F (P)) -1 ( 0 F (P))) ≥ w(X) can also be due to low values of partial derivatives of F in P. In that case, even if 0 F (P) has small width, its projection P -( 0 J F (P)) -1 ( 0 F (P))) can have large width and can be located "far away" from X. The second part of the condition (C 3 ) is used to filter the latter cases. Algo. 5 verifies (C 3 ) and returns true if it is satisfied.

This condition is used as an heuristic and is checked only if a maximal arithmetic precision specified by the user has not been reached.

Algorithm

We describe in Algo. 6 the adaptive multi-precision subdivision solver we implemented. Its input and output are fully specified in Sec. 4.1. For a given initial box X 0 where solutions of F (X) = 0 are sought, a minimal width ω ∈ R + of boxes that can be explored, an initial arithmetic precision p and a maximal precision p max , it proceeds as follows.

The simple branch and bound algorithm described in Algo. 7 is first applied to X 0 using p as arithmetic precision. Its output are a list of boxes containing a unique solution of F (X) = 0, a list of boxes with width smaller than ω that can contain solutions and a list L prec of boxes such that one of the conditions (C 1 ), (C 2 ) or (C 3 ) holds.

The precision p is then doubled and Algo. 7 is applied to boxes of L prec . This process is iterated until p reaches p max . When Algo. 7 is applied for the last time with p = p max , the condition (C 3 ) is not checked.

Sec. 4.2 describes Algo. 7.

Inria

Algorithm 6 Subdivision Solving With Adaptive Arithmetic Precision

Input: A function F : R m → R m , a box X 0 , a real number ω, an initial precision p, a maximal precision p max as described in Sec 4.1. Output: Two lists L sols , L comp and an integer r as described in Sec 4.1.

1: Let L sols , L ′ sols , L comp , L ′ comp and L prec be empty lists. 2: Push back (X 0 , true) in L prec . 3: while L prec is not empty and p ≤ p max do 4:

(L ′ sols , L ′ comp , L prec ) ← SolveWithFixedPrec(F ,L prec ,ω,p,p max ) 5:
Append boxes of L ′ sols to L sols 6:

Append boxes of L ′ comp to L comp 7:

if 2p < p max or p == p max then 8:

Set p to 2p Set p to p max 11: Determine r as described in Sec 4.1.

12: return L sols , L comp , r

Input and Output of Algo. 6

Our multi-precision branch and bound method accepts mainly as input:

• a polynomial function F (X) = (f 1 (X), . . . , f m (X))
, where X = (x 1 , . . . , x m );

• a box X 0 where solutions of F (X) = 0 are sought;

• a real number ω ∈ R + that is the minimum width of boxes to be explored during the subdivision process;

• an initial precision p and a maximal precision p max of the floating arithmetic used to represent bounds of intervals, given as the number of bits of the mantissa.

The minimal width ω allows to avoid infinite computations arising in particular when F (X) = 0 has non-regular solutions in X 0 . If it is null, the termination of the process is not ensured. The output of our multi-precision branch and bound method consists in:

• a list of boxes L sols containing solutions;

• a list of boxes L comp where other solutions could lie;

• an integer r specifying the status of the solving process.

Consider now the following properties:

(P 1 ) if X s ∈ L sols then it exists a unique X s ∈ X s ∩ int(X 0 ) s.t. F (X s ) = 0; (P 2 ) if X s ∈ int(X 0 ) is s.t. F (X s ) = 0 then it exists at most one X s ∈ L sols s.t. X s ∈ X s ; (P 3 ) if L comp is empty then for each X s ∈ int(X 0 ) s.t. F (X s ) = 0 it exists X s ∈ L sols s.t.
X s ∈ X s ; furthermore there is no solution of F = 0 on ∂(X 0 ).

Properties (P 1 ) and (P 2 ) are always satisfied when the process terminates. When, in addition, (P 3 ) is satisfied, the value of r is 0. When (P 3 ) is not satisfied because of a lack of arithmetic precision, the value of r is 1. Otherwise the value of r is 2.

RT n°476

Branch and Bound Method with Fixed Arithmetic Precision

The branch and bound algorithm described in Algo. 7 can be seen as a Depth First Search in a tree which root corresponds to an initial box X 0 , and children of a node corresponding to a box X correspond either to boxes obtained by bisecting X, i.e. cutting it in two boxes with respect to one of its components, or to a box obtained by contracting X around a solution.

The goal of such an algorithm is to find sub-boxes X of X 0 that contains a unique solution X * of F = 0 s.t. X * ∈ int(X 0 ). Notice that the krawczyk test used in Algo. 2 can only certify existence and uniqueness of a solution in int(X), and solutions lying exactly on borders of subboxes of X 0 could be missed. ǫ-inflation introduced for instance in [START_REF] Stahl | Interval Methods for Bounding the Range of Polynomials and Solving Systems of Nonlinear Equations[END_REF]Sec. 5.9] can among other circumvent this pitfall. It consists in slightly enlarging a box X obtaining X ǫ s.t. X ⊂ X ǫ . Hence Algo. 2 is applied to X ǫ to certify existence and uniqueness of a solution in X ǫ . This mechanism can however lead to find the same solutions in several neighbor boxes. X * ∈ int(X 0 ) is then verified when X ǫ ⊆ int(X 0 ). In order to optimize the number of evaluations of F and its derivatives, the absence of solutions with Algo. 1 is also tested on X ǫ .

A leaf of the tree corresponds (n 1 ) either to a box X s.t. w(X) ≤ ω that can possibly contain solutions but that are smallest than the minimum width ω given in input, (n 2 ) or to a box X that is certified to contain no solution of F = 0 with Algo. 1, (n 3 ) or to a box X s.t. X ǫ contains a unique solution X * of F = 0 and X ǫ ⊆ int(X 0 ).

An inner node of the tree corresponds (n 4 ) either to a box X s.t. existence and uniqueness of a solution is certified in X ǫ , and s.t. X ǫ ∩ ∂(X 0 ) = ∅; this case is identified in Algo. 7 with the boolean variable inflateAndBisect;

(n 5 ) either to a box X where neither the absence nor the existence of a solution can be certified; such a node has two children corresponding to boxes obtained by bisecting X.

When a node corresponding to a box X is visited, it is first checked with Algo. 8 that the actual arithmetic precision p is sufficient. If it is not, X is pushed in the list L prec ; the subtree which root is the actual node will be explored with higher precision providing that p < p max . Notice that the condition (C 3 ) presented in Sec. 3.2 is checked in Algo. 8 only if p < p max .

If the arithmetic precision p is sufficient, the type of the actual node ((n 1 ) or (n 2 ) or . . . or (n 5 )) is determined. The actual node is of type (n 3 ) if CertifyOneSolution(F ,X ǫ ) returns true (see Algo. 2) and X ǫ ⊆ int(X 0 ). If the latter condition is not fulfilled the actual node is of type (n 4 ) and the value of inflateAndBisect is set to false. The sole child of the actual node corresponds to the box 2 K F (X ǫ ) ⊂ X. If X ǫ ⊆ int(X 0 ), it is tested with Algo. 3 that the solution in X ǫ has not already been found. Finally, when the actual node is of type (n 5 ), X is bisected.

Literature proposes different strategies to choose a direction to cut a box X, and some of them are surveyed in [START_REF] Just | Subdivision strategies for boxes in branch-and-bound nonlinear solvers and verification[END_REF]. We did choose here to implement the maximum smear-diameter strategy (see [START_REF] Just | Subdivision strategies for boxes in branch-and-bound nonlinear solvers and verification[END_REF]Sec. 3.1.1,MaxSmearDiam]). In Algo. 7 the bisection is performed by the function Bisect(F ,X,ω), where the direction i to cut X = (x 1 , . . . , x i , . . . , x m ) is chosen so that w(x i ) > ω. The latter function returns a couple (X 1 , X 2 ) of boxes s.t. X 1 = (x 1 , . . . , [l(x i ), m(x i )], . . . , x m ) and X 2 = (x 1 , . . . , [m(x i ), u(x i )], . . . , x m ).

Inria

Algorithm 7 SolveWithFixedPrec(F ,L work ,ω,p,p max ) Input: A function F : R m → R m , a list L work of boxes, a real number ω, a precision p and a maximal precision p max . Output: Three lists L sols , L comp and L prec of boxes.

1: while L work is not empty do 2:

Pop the front of L work and store it in (X, inflateAndBisect).

3:

Algorithm 8 CheckPrec(F ,X,ω,p,p max ,inflateAndBisect) Input: F, X, ω, p, p max as described in Sec. 4.1, and a boolean variable inflateAndBisect.

Output: true if the arithmetic precision has to be increased to treat X; false otherwise. 1: if inflateAndBisect then 2:

(X 1 , X 2 ) = Bisect(F ,X,ω).

3:

Let X ǫ , X 1 ǫ , X 2 ǫ be obtained by ǫ-inflation of X, X 1 , X 2 .

4:

return CheckCond12(F ,X ǫ ,X 1 ǫ ,X 2 ǫ ) or ( CheckCond3(F ,X ǫ ) and p < p max ) 5: else 6:

return (CheckCond3(F ,X) and p < p max ) Table 1: Comparaison of strategies (1),( 2),( 3) and ( 4). m is the number of equations and variables of the system, d is the degree of its polynomial equations. n is the number of boxes that have been explored, t the sequential time in seconds spent for the process on a Intel(R) Core i7-5600U CPU @ 2.60GHz . We note that using evaluation at order 2 (i.e. strategies (1) and ( 2) versus ( 3) and ( 4)) brought an important gain both in terms of number of boxes explored and time, and higher m and d are, higher is this gain. The gain allowed by Krawczyk operator at order 2 (strategy (1) versus strategy (2)) is more perceptible for small values of m and d. Notice that applying strategy (1) instead of strategy (2) does not induce additional evaluations since values of the Hessian matrix are computed for the evaluation at order 2.

http://www.sagemath.org

http://www.boost.org/doc/libs/1_60_0/libs/numeric/interval/doc/interval.htm

http://www.loria.fr/ ~moroz/software.html

http://www.loria.fr/ ~rimbach/ RT n°476

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-0803

if CheckPrec(F ,X,ω,p,p max ,inflateAndBisect) then

4:

Push back X in L prec . if w(X) ≤ ω then 7:

Push back X in L comp . if inflateAndBisect then 10:

Let X ǫ be obtained by ǫ-inflation of X. Let X ǫ be a copy of X.

13:

if CertifyNoSolution(F ,X ǫ ) then 14:

22:

Push (X 1 , true) and (X 2 , true) in the front of L work .

23: return L sols , L comp , L prec

Benchmarks

We propose here some experimental results to justify the chosen strategy for certifying nonexistence or existence and uniqueness of a solution in a box. We compare four possible strategies for functions CertifyNoSolution and CertifyOneSolution :

(1) Evaluation at order 2 and Krawczyk operator at order 2, as described in previous algorithms;

(2) Evaluation at order 2 and Krawczyk operator as defined in Eq. 7: obtained by replacing 2 K F (X) by K F (X) in previous algorithms;

(3) Evaluation at order 1 and Krawczyk operator as defined in Eq. 7: obtained by replacing 2 F by 1 F and 2 K F (X) by K F (X) in previous algorithms;

(4) Natural interval extension and Krawczyk operator as defined in Eq. 7: obtained by replacing 2 F by 0 F and 2 K F (X) by K F (X) in previous algorithms.

We did use these four strategies to solve dense, randomly generated, systems with m polynomial equations of degree d with integer coefficients of bit-size 8, for (m, d) ∈ {(2, 64), (2, 128), [START_REF] Krawczyk | Interval iterations for including a set of solutions[END_REF]16), [START_REF] Krawczyk | Interval iterations for including a set of solutions[END_REF]32), (4, 8), (5, 4)}. Table 1 gives for each strategy and each couple (m, d) the number n of boxes that have been explored and the sequential time t in seconds spent to solve the system.

RT n°476


Oualid Benyahia
email: oualid.benyahia@univ-st-etienne.fr

Christine Largeron
email: christine.largeron@univ-st-etienne.fr

Baptiste Jeudy
email: baptiste.jeudy@univ-st-etienne.fr

Osmar R Zaïane
email: zaiane@cs.ualberta.ca

DANCer: Dynamic Attributed Network with Community Structure Generator

Keywords: Social network mining, Attributed graph, Synthetic data generator

We propose a new generator for dynamic attributed networks with community structure which follow the known properties of real-world networks such as preferential attachment, small world and homophily. After the generation, the different graphs forming the dynamic network as well as its evolution can be displayed in the interface. Several measures are also computed to evaluate the properties verified by each graph. Finally, the generated dynamic network, the parameters and the measures can be saved as a collection of files.

Introduction

The proliferation of complex information networks in diverse fields of application has led to the proposal of a panoply of methods to analyze and discover relevant patterns in these networks. However, evaluating these methods and the comparison of the different approaches are not very easy due to the lack of large real networks with ground truth freely accessible to researchers. The alternative consists in using synthetic data provided by generators. There is a large bibliography regarding generation for static graphs, including the classic Erdős-Rényi (ER) model which generates random graphs or the Barabási-Albert (BA) model that generates random scale-free networks, but very few generators allow the construction of evolving graphs, exhibiting or not a community structure and, none of them takes into account the attribute values of the vertices. The interest of community detection, link prediction and more generally pattern discovery in dynamic networks where vertices are associated with attributes led us to develop the generator DANCer for attributed dynamic graphs with embedded community structure. This generator is an extended version of a previous generator, ANC dedicated to static graphs [START_REF] Largeron | Generating Attributed Networks with Communities[END_REF].

Model

An attributed dynamic network generated by DANCer is represented by (1) a sequence of T attributed graphs G i = (V i , E i ), i ∈ {1, . . . , T }, where V i is a set of vertices, E i a set of undirected edges and where for each vertex v ∈ V i and each real attribute A ∈ A, v A denotes the attribute value of A assigned to vertex v and (2) a sequence of T partitions P i of V i , i ∈ {1, . . . , T } which gives a community for each vertex in the corresponding graph G i , i ∈ {1, . . . , T }. Each partition allows to define a community structure on a graph (i.e., the network at a single timestamp) in such a way that the nodes are grouped into sets densely connected and relatively homogeneous with regard to the attributes, while they are less connected to vertices belonging to other groups, and less similar with regard to their attributes.

The generation of the network is carried out in two phases. In phase one, an initial graph

G 1 = (V 1 , E 1
) is built while respecting the well-known network properties such as preferential attachment, small world or homophily and, in the second phase, this initial graph is modified through two kinds of operations. The first set of operations, called micro operations, consist in removing or adding vertices and edges or updating their attributes whereas the second kind of operations is applied on the communities, i.e., at a macro level. They consist in (1) migrating members of a community to either a new community or an existing one, (2) splitting a community into two new sub-communities and (3) merging two existing communities into a single one.

Software Overview

The user interface has three panels as shown in Fig. 1. In the parameter panel, the user selects the dynamic generator parameters presented in Table 1 3 . Note that a seed is used for the random number generator. It can be saved to reproduce exactly the same network.

The visualization panel allows to display the generated network and its dynamic evolution. Each graph in the sequence can be selected with a timestamp scrollbar and viewed separately (Fig. 1). This panel can also display the size and the evolution of the different communities in the sequence of graphs according to the macro dynamic operations (split, merge and migrate) (see Fig. 2).

The sequence of attributed graphs is built while preserving properties of real networks and several measures, like modularity, clustering coefficient, diameter, expected and observed homophily or within inertia rate are computed on each graph of the dynamic network to describe its properties4 . The changes in these different measures on the sequence of graphs are presented at the bottom of the interface in the measure panel (Fig. 3). The bottom of the user interface includes also a panel displaying the distribution of vertex degrees on each graph of the sequence as shown in Fig. 4.

The generated dynamic network can be saved as a collection of files. For each graph of the sequence, a file indicates the composition of the graph (vertices and edges) and a "parameters" file enumerates all the parameters used by the generator. The graph measures and community dynamics can also be saved in separated files.

Conclusion

The software DANCer and a detailed user manual5 are available under the terms of the GNU Public Licence. Note that our generator can trivially be extended to produce multiplex networks, where all nodes are omnipresent in all levels and intra-level edges connect the representations of a node from one level to the other. This conversion is possible by simply converting each timestamp graph into a layer of the multiplex network and adding the necessary intra-level edges.

Fig. 1 .

 1 Fig. 1. User interface of the generator DANCer.

Fig. 2 .

 2 Fig. 2. Community dynamics display in visualization panel.

Fig. 3 .

 3 Fig. 3. Measures panel of the generator.

Fig. 4 .

 4 Fig. 4. Degree distribution panel.

Table 1 .

 1 Description of the dynamic network generator parameters

	Parameter	Domain Description
	Micro operations
	P roba M icro	[0, 1]	A threshold to select if the micro dynamic updates are
			performed or not
	Add V ertex	[0, 1]	Ratio defining the number of vertices inserted
	Remove V ertex	[0, 1]	Ratio defining the number of vertices removed
	U pdate Attr.	[0, 1]	Ratio defining the number of attributes updated
	Add Btw. Edges	[0, 1]	Ratio defining the number of between edges inserted
	Remove Btw. Edges	[0, 1]	Ratio defining the number of between edges removed
	Add W th. Edges	[0, 1]	Ratio defining the number of within edges inserted
	Remove W th. Edges	[0, 1]	Ratio defining the number of within edges removed
	Macro operations
	P removeEdgeSplit	[0, 1]	Proba. to remove an edge between two vertices in the
			previously same community when splitting a community
	T imestamps	N +	Number of graphs generated
	P roba M erge	[0, 1]	Probability to perform the merge operation
	P roba Split	[0, 1]	Probability to perform the split operation
	P roba M igrate	[0, 1]	Probability to perform the migrate vertices operation

The reader is referred to[START_REF] Largeron | Generating Attributed Networks with Communities[END_REF] for the static network generator parameters http:// journals.plos.org/plosone/article?id=10.1371/journal.pone.0122777.

see[START_REF] Largeron | Generating Attributed Networks with Communities[END_REF] for a more detailed presentation of the properties and corresponding measures.

http://perso.univ-st-etienne.fr/largeron/DANCer_Generator/